% 
% 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_housing_starts(); n_all = length(Y_all); 

fh = figure; p_data = plot( Y_all, 'bx' ); axis tight; hold on; 
xlabel( 'months since Jan 1966' ); ylabel( 'monthly housing starts' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/prob_4_6_global_plt', 'epsc' ); 

s = 12; 

%
% GLOBALLY constant models: 
%

% 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 = 120; 
n = 108; 
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/prob_4_6_global_plt', 'epsc' ); 

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

% 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/prob_4_6_global_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/Chapter3/prob_4_6_global_r_k', 'epsc' ); 

% 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, 'r-' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/prob_4_6_global_plt', 'epsc' ); 

MSE_glb_ind = mean( ( Y_all( (n+1):end ) - Y_extrap ).^2 ); 

% Part (b) globally constant linear trend model with trigonometric terms:
%
% m=12 is way to many harmonics: 
%
m = 4; 
X_all = gen_global_linear_seasonal_trigonometric_X(s,m,n_all); 

% get a SUBSET of the data: 
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 insample predictions:
y_hat = X * beta; 
figure(fh); p_glrt = plot( y_hat, '-c' ); 
legend( [p_data,p_glri,p_glrt], {'raw data', 'GLR indicator', 'GLR trigonometric'}, 'location', 'northwest' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/prob_4_6_global_plt', 'epsc' ); 

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

% 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);
figure(sf); sh_glr_trig=stem( 0:(n_r_k-1), r_k, 'cx' ); hold on; 
xlabel( 'lag (k)' ); ylabel( 'sample based autocorrelation r_k' ); axis tight; 
legend( [sh_glr_ind,sh_glr_trig], {'GLR indicator', 'GLR trig'} ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/prob_4_6_global_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/prob_4_6_global_r_k', 'epsc' ); 

% 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, 'c-' ); line( [n,n],[30,140] ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/prob_4_6_global_plt', 'epsc' ); 

MSE_glb_trig = mean( (Y_all( (n+1):end ) - Y_extrap).^2 ); 

%
% LOCALLY constant models: 
%

fh = figure; p_data = plot( Y_all, 'bx', 'markersize', 3 ); axis tight; hold on; 
xlabel( 'months since Jan 1955' ); ylabel( 'monthly sales of U.S. cars' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/prob_4_6_local_plt', 'epsc' ); 

% Part (c) locally constant linear trend model with optimal alpha and seasonal indicators: 
%
omega_grid = linspace(0.1,0.9,250); 
[w_opt,z_hat_1,beta_hats,et,MSE,omega_grid,MSE_omega] = locally_constant_indicator_model_optimum(Y,s,omega_grid); 

SSE_omega = n * MSE_omega; 
figure; plot( omega_grid, SSE_omega, '-x' ); grid on; 
xlabel('\alpha'); ylabel('SSE(\alpha)'); title( 'SSE as a function of \alpha' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/prob_4_6_indicator_omega_vs_SSE', 'epsc' ); 

% plot this optimal relaxation factor (in-sample): 
figure(fh); ph_opt_o_ind = plot( 1:n, z_hat_1, 'r-' ); 
legend( [p_data,ph_opt_o_ind], {'raw data', 'Opt. Local Indicator'}, 'location', 'northwest' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/prob_4_6_local_plt', 'epsc' ); 

% compute the autocorrelation function representing the residuals ... these are the same as 
% the autocorrelation function representing the ONE-STEP-AHEAD forecast errors of our given model
%
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_loc_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_loc_ind], {'Local indicator'} ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/prob_4_6_local_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/prob_4_6_local_r_k', 'epsc' ); 

% lets compute this model for the second half of the timeseries (the out of sample data): 
%
beta_0 = beta_hats(:,end); 
[z_hat_2,beta_hats_2,et_2,MSE_2] = locally_constant_indicator_model(Y_all(n+1:end),s,w_opt,beta_0);
T_extrap = (n+1):n_all;  
figure(fh); plot( T_extrap, z_hat_2, 'r-' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/prob_4_6_local_plt', 'epsc' ); 

MSE_loc_ind = MSE_2; 


% Part (d) locally constant linear trend model with optimal alpha and trigonometric functions:
%
omega_grid = linspace(0.1,0.9,250); 
[w_opt,z_hat_1,beta_hats,et,MSE,omega_grid,MSE_omega] = locally_constant_trigonometric_model_optimum(Y,s,m,omega_grid); 

SSE_omega = n * MSE_omega; 
figure; plot( omega_grid, SSE_omega, '-x' ); grid on; 
xlabel('\alpha'); ylabel('SSE(\alpha)'); title( 'SSE as a function of \alpha' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/prob_4_6_trig_omega_vs_SSE', 'epsc' ); 

% plot this optimal relaxation factor (in-sample): 
figure(fh); ph_opt_o_trig = plot( 1:n, z_hat_1, 'c-' ); 
legend( [p_data,ph_opt_o_ind,ph_opt_o_trig], {'raw data', 'Opt. Local Indicator', 'Opt. Local Trig'}, 'location', 'northwest' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/prob_4_6_local_plt', 'epsc' ); 

% compute the autocorrelation function representing the residuals ... these are the same as 
% the autocorrelation function representing the ONE-STEP-AHEAD forecast errors of our given model
%
et2 = (et - mean(et))/std(et); et2=et2(:); n_r_k = round(0.2*length(et));
r_k = aasampleunbiasedautoc(et2.',n_r_k);
figure(sf); sh_loc_trig=stem( 0:(n_r_k-1), r_k, 'cx' ); hold on; 
xlabel( 'lag (k)' ); ylabel( 'sample based autocorrelation r_k' ); axis tight; 
legend( [sh_loc_ind,sh_loc_trig], {'Local indicator', 'Local Trig'} ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/prob_4_6_local_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/prob_4_6_local_r_k', 'epsc' ); 

% lets compute this model for the second half of the timeseries (the out of sample data): 
%
beta_0 = beta_hats(:,end); 
[z_hat_2,beta_hats_2,et_2,MSE_2] = locally_constant_trigonometric_model(Y_all(n+1:end),s,m,w_opt,beta_0);
T_extrap = (n+1):n_all;  
figure(fh); plot( T_extrap, z_hat_2, 'c-' ); line( [n,n],[30,140] ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/prob_4_6_local_plt', 'epsc' ); 

MSE_loc_trig = MSE_2; 

% display all of the MSE for each of the methods: 
%
fprintf('cummulative MSE results for each method'); 
[ MSE_glb_ind, MSE_glb_trig, MSE_loc_ind, MSE_loc_trig ] 

return;