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

Y = load_us_lumber(); Y = Y/1e3; n = length(Y); 

fh = figure; os = plot( Y, '-kx' ); grid on; axis tight; 
xlabel('time'); ylabel('annual lumber production (K)'); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_4_lp', 'epsc' ); 

% pick an initial value smooth S_0^{[1]}:
S_0 = Y(1); 
S_0 = mean(Y); 
S_0 = median(Y); 

% pick value of the smoothing factor w:
w = 0.9; 
[ Ys, et, MSE ] = simple_exp_smoothing( Y, w, S_0 ); 

figure(fh); hold on; wp9 = plot( Ys, '-ob' ); 
legend( [os,wp9], {'original data', 'discount \omega=0.9'} ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_4_lp', 'epsc' ); 

% find the optimal relaxation factor: 
%
omega_grid = linspace( 0.7, 0.99, 100 ); 
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

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/prob_3_4_mse_plot', 'epsc' ); 

[min_value,min_indx] = min( mse_omega ); w_opt = omega_grid(min_indx);

[ Ys, et, MSE ] = simple_exp_smoothing( Y, w_opt, S_0 ); 

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

% look at the error residuals:
mean(et)

figure; stem(et, 'x', 'markersize', 10 ); title( 'error residuals' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_4_err_residuals', 'epsc' ); 

addpath('/home/weatherw/Trash/Poularikas/AdaptiveFilteringCode');
n_r_k = round(0.2*length(et));
et2 = (et-mean(et))/std(et);
r_k = aasampleunbiasedautoc( et2, n_r_k );

figure; stem( 0:(length(r_k)-1), r_k, 'x', 'markersize', 10 ); 
two_sde = 2*sqrt(1/n);
hold on; 
plot( 0:(length(r_k)-1),  two_sde * ones(1,length(r_k)), '-r' ); 
plot( 0:(length(r_k)-1), -two_sde * ones(1,length(r_k)), '-r' ); 
title( 'sample autocorrelations' );
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_4_samp_autocorrelatio', 'epsc' );