% 
% 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/' ); 

Y_all = load_monthly_car_sales(); n_all = length(Y_all); 

fh = figure; p_data = plot( Y_all, '-bx' ); axis tight; hold on; 
xlabel( 'months since Jan 1960' ); ylabel( 'monthly sales of cars in Quebec' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/example_4_1_plt', 'epsc' ); 

s = 12; 

% Part (a) globally constant linear trend model with seasonal indicators: 
%
X_all = gen_global_linear_seasonal_indicator_X(s,n_all); 

% get a SUBSET of the data: 
n = 96; Y = Y_all(1:n); X = X_all(1:n,:);

[beta,sC,s_hat,t_stats,SSE,SSR,SSTO,R2,MSR,MSE,F_ratio,MI] = wwx_regression( X, Y, 0 );

% plot the predictions:
y_hat = X * beta; 
figure(fh); p_glri = plot( y_hat, '-r.' ); 
legend( [p_data,p_glri], {'raw data', 'GLR indicator'}, 'location', 'northwest' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/example_4_1_plt', 'epsc' ); 

% make a table of results:
[ beta, sC, t_stats ] 

% plot the model extrapolated: 
% 
T_extrap = (n+1):n_all;  
Y_extrap = X_all((n+1):end,:) * beta; 
figure(fh); plot( T_extrap, Y_extrap, 'k-', 'markersize', 10 ); line( [n,n], [5,27] ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/example_4_1_plt', 'epsc' ); 

% compute the autocorrelation function representing the residuals ... these are the same as 
% the autocorrelation function representing the ONE-STEP-AHEAD forecast errors when we have a 
% global model 
%
et  = Y - X * beta; 
et2 = (et - mean(et))/std(et); et2=et2(:); n_r_k = round(0.2*length(et));
r_k = aasampleunbiasedautoc(et2.',n_r_k);
sf=figure; sh_glr_ind=stem( 0:(n_r_k-1), r_k, 'rx' ); hold on; 
xlabel( 'lag (k)' ); ylabel( 'sample based autocorrelation r_k' ); axis tight; 
legend( [sh_glr_ind], {'GLR indicator'} ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/example_4_1_r_k', 'epsc' ); 

% plot the value of the standard errors of the sample autocorrelations: 
se_r_k = 1.96 / sqrt(n);

figure(sf); hold on; plot( 0:(n_r_k-1), -2*se_r_k*ones(n_r_k,1), 'r-' ); 
figure(sf); hold on; plot( 0:(n_r_k-1), +2*se_r_k*ones(n_r_k,1), 'r-' ); axis tight; 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/example_4_1_r_k', 'epsc' ); 

return;