% 
% Generates seasonal ARIMA models for the gas usage data set.
% 
% Written by:
% -- 
% John L. Weatherwax                2009-04-21
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; drawnow; clear functions; clear; rehash; clc; 

if( exist('/home/weatherw/Trash/Poularikas/AdaptiveFilteringCode') ) addpath('/home/weatherw/Trash/Poularikas/AdaptiveFilteringCode'); end
addpath( '../Chapter2/' ); 
addpath( '../Chapter4/' ); 

[Y_all] = load_gas_usage(); n_all = length(Y_all);

n = n_all - 10; Y = Y_all(1:n); % hold out some simple data for testing (if desired)

fh = figure; p_data = plot( 1:n, Y, '-bx' ); axis tight; hold on; 
xlabel( 'Time' ); ylabel( 'Monthly Residential Gas Usage in Iowa City' );
saveas( gcf, '../../WriteUp/Graphics/Chapter6/month_gas_data_plt', 'epsc' ); 

% plot the sample autocorrelation function of the timeseries directly 
plot_SACF( Y, 1 ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter6/month_gas_data_direct_sacf', 'epsc' ); 

% observe that this series is non-stationary (has a changing mean) ... plot the first seasonal (s=12) difference 

Yd = sdiff(Y,12,1); nd = length(Yd); 

fhd = figure; p_data = plot( 1:nd, Yd, '-bx' ); axis tight; hold on; 
xlabel( 'Time' ); ylabel( 'The first seasonal difference of the gas usage data' );
saveas( gcf, '../../WriteUp/Graphics/Chapter6/month_gas_diff_data_plt', 'epsc' ); 

% plot the sample autocorrelation function of the timeseries directly 
r_k = plot_SACF( Yd, 1 ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter6/month_gas_data_sdiff_sacf', 'epsc' ); 

% assume a particular SARIMA model based on the above plots ... (0,0,1)(0,1,1)_{12}

% estimate the mean value if the first seasonal 
[muhat,sigmahat,muci,sigmaci] = normfit(Yd,0.05);
fprintf('mean estimate = %10.6f; standard error = %10.6f\n',muhat,diff(muci));

%
% See the file model_gas_section_6_7_2.R for code to : 
%

% estimate the coefficients of the model (in R) ... 

% compute the one-step-ahead predicitions, z_n(1), and their residuals using the estimated 
% model coefficients

% compute the Portmanteau statistic Q ... skipped for time 

% compute three forecasts (and their variances)

return;