%   
% 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; clc; clear; clear functions; 

% calculate the correlation of these residuals:
if( exist('/home/weatherw/Trash/Poularikas/AdaptiveFilteringCode') ) addpath('/home/weatherw/Trash/Poularikas/AdaptiveFilteringCode'); end

Y_save = load_computer_software_sales(); n = length(Y_save); 

% downsample our data: 
n = 60; Y = Y_save(1:n); Y_test = Y_save( (n+1):end );

%--
% Test SIMPLE EXPONENTIAL SMOOTHING: 
%--
fh = figure; os = plot( Y, '-gx' ); grid on; axis tight; 
xlabel('time'); ylabel('monthly sales of software packages' ); title( 'exponential smoothings' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_plt', 'epsc' ); 

% pick an initial value of for the smooth
S_0 = Y(1);

% pick an initial value of the smoothing factor w to try:
if( 0 ) 
  w = 0.9; 
  [ S_n, et, MSE ] = simple_exp_smoothing( Y, w, S_0 );

  figure(fh); hold on; wp9 = plot( S_n, '-or' ); 
  legend( [os,wp9], {'original data', 'simple smooth discount \omega=0.9'}, 'location', 'northwest' ); 
  saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_plt', 'epsc' ); 
end

% find the optimal relaxation factor for simple exponential smoothing: 
%
omega_grid = linspace( 0.1, 0.99, 500 ); 
mse_omega  = zeros(1,length(omega_grid));
for ii=1:length(omega_grid),
  [ Ys, et, mse_omega(ii) ] = simple_exp_smoothing( Y, omega_grid(ii), S_0 ); 
end

fo=figure; plot( omega_grid, mse_omega, 'bx-' ); 
xlabel( 'relaxation factor (omega)' ); ylabel( 'Mean Square Error' ); title( 'MSE as a function of relaxation \omega' );
grid on; axis tight; 
axis( [ omega_grid(1), omega_grid(end), 1700, 2400 ] ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_mse_plot', 'epsc' ); 

[min_value,min_indx] = min( mse_omega ); w_opt = omega_grid(min_indx);
fprintf('Simple Exponential Smoothing Optimal Omega = %10.6f\n', w_opt); 
[ Ys, et, MSE ] = simple_exp_smoothing( Y, w_opt, S_0 ); 

figure(fh); hold on; wp_ses_opt = plot( Ys, '-ob' ); 
legend( [os,wp_ses_opt], {'original data', 'SES (optimal discount)'}, 'location', 'northwest' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_plt', 'epsc' ); 

% plot the autocorrelation of the residuals: 
%
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_ses=stem( 0:(n_r_k-1), r_k, 'x' ); 
xlabel( 'lag (k)' ); ylabel( 'sample based autocorrelation r_k' ); %axis tight; 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_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_3_14_r_k', 'epsc' ); 

% predict the out-of-sample points: 
%
S_0 = Ys(end); 
[ ses_Ys_test, ses_et_test, ses_MSE_test ] = simple_exp_smoothing( Y_test, w_opt, S_0 ); 

%--
% Test DOUBLE EXPONENTIAL SMOOTHING:   
%--

% find the optimal relaxation factor for double exponential smoothing: 
%
omega_grid = linspace( 0.1, 0.99, 500 ); 
mse_omega  = zeros(1,length(omega_grid));
for ii=1:length(omega_grid),
  [ zHats, S_1, S_2, et, mse_omega(ii) ] = double_exp_smoothing( Y, omega_grid(ii) ); 
end


figure(fo); hold on; plot( omega_grid, mse_omega, 'rx-' ); 
xlabel( 'relaxation factor (omega)' ); ylabel( 'Mean Square Error' ); title( 'MSE as a function of relaxation \omega' );
grid on; axis tight; 
axis( [ omega_grid(1), omega_grid(end), 1700, 2400 ] ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_mse_plot', 'epsc' ); 

[min_value,min_indx] = min( mse_omega ); w_opt = omega_grid(min_indx);
fprintf('Double Exponential Smoothing Optimal Omega = %10.6f\n', w_opt); 
[ zHats, S_1, S_2, et, MSE ] = double_exp_smoothing( Y, w_opt );

figure(fh); hold on; wp_des_opt = plot( zHats, '-or' ); 
legend( [os,wp_ses_opt,wp_des_opt], {'original data', 'SES', 'DES'}, 'location', 'northwest' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_plt', 'epsc' ); 

% plot the residuals: 
%
et2 = (et - mean(et))/std(et);  n_r_k = round(0.2*length(et));
r_k = aasampleunbiasedautoc(et2.',n_r_k); 
figure(sf); hold on; sh_des=stem( 0:(n_r_k-1), r_k, 'rx' ); 
xlabel( 'lag (k)' ); ylabel( 'sample based autocorrelation r_k' ); %axis tight; 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_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-' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_r_k', 'epsc' ); 

% predict the out-of-sample points: 
%
S_0_1 = S_1(end); 
S_0_2 = S_2(end);
[ des_Ys_test, S_1_test, S_2_test, des_et_test, des_MSE_test ] = double_exp_smoothing( Y_test, w_opt, S_0_1, S_0_2 ); 

%--
% TRIPLE EXPONENTIAL SMOOTHING
%--

% find the optimal relaxation factor for double exponential smoothing: 
%
omega_grid = linspace( 0.1, 0.99, 500 ); 
mse_omega  = zeros(1,length(omega_grid));
for ii=1:length(omega_grid),
  [ zHats, S_1, S_2, S_3, et, mse_omega(ii) ] = triple_exp_smoothing( Y, omega_grid(ii) ); 
end

figure(fo); hold on; plot( omega_grid, mse_omega, 'kx-' ); 
xlabel( 'relaxation factor (omega)' ); ylabel( 'Mean Square Error' ); title( 'MSE as a function of relaxation \omega' );
grid on; 
%axis( [ omega_grid(1), omega_grid(end), 1700, 2400 ] ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_mse_plot', 'epsc' ); 

[min_value,min_indx] = min( mse_omega ); w_opt = omega_grid(min_indx);
fprintf('Triple Exponential Smoothing Optimal Omega = %10.6f\n', w_opt); 
[ zHats, S_1, S_2, S_3, et, MSE ] = triple_exp_smoothing( Y, w_opt );

figure(fh); hold on; wp_tes_opt = plot( zHats, '-ok' ); 
legend( [os,wp_ses_opt,wp_des_opt,wp_tes_opt], {'original data', 'Simple ES', 'Double ES', 'Triple ES'}, 'location', 'northwest' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_plt', 'epsc' ); 

% plot the residuals: 
%
et2 = (et - mean(et))/std(et);  n_r_k = round(0.2*length(et));
r_k = aasampleunbiasedautoc(et2.',n_r_k); 
figure(sf); hold on; sh_tes=stem( 0:(n_r_k-1), r_k, 'kx' ); 
xlabel( 'lag (k)' ); ylabel( 'sample based autocorrelation r_k' ); %axis tight; 
legend( [sh_ses,sh_des,sh_tes],  {'Simple ES', 'Double ES', 'Triple ES'}, 'location', 'northwest' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_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-' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_r_k', 'epsc' ); 

% predict the out-of-sample points: 
%
S_0_1 = S_1(end); 
S_0_2 = S_2(end);
S_0_3 = S_3(end);
[ tes_Ys_test, S_1_test, S_2_test, S_3_test, tes_et_test, tes_MSE_test ] = triple_exp_smoothing( Y_test, w_opt, S_0_1, S_0_2, S_0_3 ); 

% compare the three methods: 

fprintf('the values themselves\n');
[ Y_test, ses_Ys_test, des_Ys_test, tes_Ys_test ] 

fprintf('the residuals\n');
[ ses_et_test, des_et_test, tes_et_test ] 

fprintf('the MSE\n'); 
[ ses_MSE_test, des_MSE_test, tes_MSE_test ] 

% lets plot the predictions for each of the three models: 
figure(fh); plot( (n+1):(n+length(Y_test)), Y_test, '-og', 'markersize', 10 ); hold on; 
plot( (n+1):(n+length(Y_test)), ses_Ys_test, '-ob'); 
plot( (n+1):(n+length(Y_test)), des_Ys_test, '-or'); 
plot( (n+1):(n+length(Y_test)), tes_Ys_test, '-ok'); axis tight; 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_plt', 'epsc' );