% 
% Generates a locally constant seasonal model using seasonal indicators on the 
% new plant equipment 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/' ); 

Y_all = load_quarterly_new_plant_equipment(); n_all = length(Y_all); Y_all = log( Y_all ); 

fh = figure; p_data = plot( Y_all, '-bx' ); axis tight; hold on; 
xlabel( 'quarters since 1964' ); ylabel( 'quartly new plant equipment orders' );
saveas( gcf, '../../WriteUp/Graphics/Chapter4/loc_const_indicator_plt', 'epsc' ); 

s = 4; 

% downsample to a smaller number of elements
n = 44; Y = Y_all(1:n); n_all = length(Y_all); 

alpha = 0.2; omega = 1 - alpha; 
[z_hat_1,beta_hats] = locally_constant_indicator_model(Y,s,omega); 

figure(fh); ph_trial_omega = plot( z_hat_1, 'r.-' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/loc_const_indicator_plt', 'epsc' ); 

% find the optimal relaxation factor: 
%
alpha_grid = [ 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.2, 1.3, 1.5, 1.6 ]; omega_grid = 1 - alpha_grid; 
[w_opt,z_hat_1,beta_hats,et,MSE,omega_grid,MSE_omega] = locally_constant_indicator_model_optimum(Y,s,omega_grid); 

% display [ alpha, SSE(alpha) ]
SSE_omega = n * MSE_omega; 
[ 1-omega_grid(:), n*MSE_omega(:) ] 

figure; plot( alpha_grid, SSE_omega, '-x' ); grid on; 
axis( [0.2, 1.6, 0.03, 0.1] ); 
xlabel('\alpha'); ylabel('SSE(\alpha)'); title( 'SSE as a function of \alpha' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/loc_const_indicator_sse_plt', 'epsc' ); 

% plot this optimal relaxation factor (in-sample): 
figure(fh); ph_opt_omega = plot( 1:n, z_hat_1, 'g.-' ); 
legend( [p_data,ph_trial_omega,ph_opt_omega], {'original data', '\alpha=0.2', '\alpha=1.2=optimal'}, 'location', 'northwest' )
saveas( gcf, '../../WriteUp/Graphics/Chapter4/loc_const_indicator_plt', '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);

figure(fh); ph_opt_outsample = plot( (n+1):n_all, z_hat_2, '-gd', 'markersize', 10 );
legend( [p_data,ph_trial_omega,ph_opt_omega,ph_opt_outsample], {'original data', '\alpha=0.2', '\alpha=1.2=optimal (insample)', 'optimal \alpha (outsample)'}, 'location', 'northwest' )
title('locally linear model fit to the new plant data set'); line( [n,n], [2.3,3.7] ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/loc_const_indicator_plt', 'epsc' ); 

% plot or display the sample autocorrelations
% 
plot_SACF( et );
title( 'SACF: only the lag k=8 maybe significant' );