%   
% 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

[T,Y] = load_student_enrollment(); n = length(Y); 

% Test simple exponential smoothing: 
% 
fh = figure; os = plot( T, Y/1.e3, '-gx' ); grid on; axis tight; 
xlabel('time'); ylabel('student enrollment (K)'); title( 'University of Iowa student enrollment, 1951/52 to 1979/1980' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/example_3_4_ts_plot', 'epsc' ); 

% pick an initial value of the smoothing factor w to try:
w = 0.9; 
[ S_n, dum1, dum2, et, MSE ] = double_exp_smoothing( Y, w );

figure(fh); hold on; wp9 = plot( T, S_n/1.e3, '-or' ); 
legend( [os,wp9], {'original data', 'double exponential smooth (\omega=0.9)'}, 'location', 'northwest' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/example_3_4_ts_plot', 'epsc' ); 

% find the optimal relaxation factor for simple exponential smoothing: 
%
omega_grid = linspace( 0.01, 0.9, 500 ); 
mse_omega  = zeros(1,length(omega_grid));
for ii=1:length(omega_grid),
  [ Ys, dum1, dum2, et, mse_omega(ii) ] = double_exp_smoothing( Y, omega_grid(ii) );
end

figure; 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; 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/example_3_4_mse_plot', 'epsc' ); 

[min_value,min_indx] = min( mse_omega ); w_opt = omega_grid(min_indx);
fprintf('optimal omega = %10.6f; alpha = %10.6f\n', w_opt, 1 - w_opt );

%w_opt = 1 - 0.87; 
[ Ys, dum1, dum2, et, MSE ] = double_exp_smoothing( Y, w_opt );

figure(fh); hold on; wp_opt = plot( T, Ys/1.e3, '-b' ); 
legend( [os,wp9,wp_opt], {'original data', 'discount \omega=0.9', 'optimal discount'}, 'location', 'northwest' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/example_3_4_ts_plot', 'epsc' ); 

% plot 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=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/example_3_4_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/example_3_4_r_k', 'epsc' );