% 
% Generates a locally constant seasonal model using trigonometric functions 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_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/loc_const_trig_plt', 'epsc' ); 

s = 12; m = 2; 

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

alpha = 0.1; omega = 1 - alpha; 
[z_hat_1,beta_hats] = locally_constant_trigonometric_model(Y,s,m,omega); 

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

% find the optimal relaxation factor: 
%
alpha_grid = [ 0.03:0.01:0.1, 0.15, 0.2 ]; omega_grid = 1 - alpha_grid; 
[w_opt,z_hat_1,beta_hats,et,MSE,omega_grid,MSE_omega] = locally_constant_trigonometric_model_optimum(Y,s,m,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_trig_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.1', '\alpha=optimal'}, 'location', 'northwest' )
saveas( gcf, '../../WriteUp/Graphics/Chapter4/loc_const_trig_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_trigonometric_model(Y_all(n+1:end),s,m,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.1', '\alpha=optimal (insample)', 'optimal \alpha (outsample)'}, 'location', 'northwest' )
title('locally linear model fit to the new plant data set'); line( [n,n], [5,27] ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter4/loc_const_trig_plt', 'epsc' ); 

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