% 
% 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/' ); 
addpath( '../Chapter4/' ); % plot_SACF.m

Y_all = load_yield_data(); n_all = length(Y_all); 

fh = figure; p_data = plot( Y_all, '-bx' ); axis tight; hold on; 
xlabel( 'months since 1961' ); ylabel( 'difference between mortgate yield and goverment load yields' );
%saveas( gcf, '../../WriteUp/Graphics/Chapter5/yield_data_plt', 'epsc' ); 

% note this data is clearly non-stationary ... 

r_k = plot_SACF( Y_all ); 
axis( [0.25,20,-0.5,0.9] ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter5/yield_data_sacf', 'epsc' ); 

% attempt to find the optimal difference to take ... not the minimum at d=1 => first order differences 
Y=Y_all; figure; plot( [ std(Y), std(diff(Y)), std(diff(Y,2)), std(diff(Y,3)), std(diff(Y,4)), std(diff(Y,5)) ], '-x' ); 

figure; plot( diff(Y), 'x' ); 

% override with the book's values: 
r_k = [ 0.841     0.683     0.584     0.515     0.457     0.427     0.405     0.386     0.361     0.321   
        0.329     0.338     0.337     0.294     0.231     0.166     0.126     0.062     0.047     0.042  ].';
r_k = r_k(:); 

% now lets compute the sample partial autocorrelation function using this data:
phi_kk = spacf( r_k ); n_k = length(phi_kk); 

sf=figure; stem( 1:n_k, phi_kk,  'bx', 'linewidth', 2 ); 
xlabel( 'lag (k)' ); ylabel( 'sample partial autocorrelations phi_kk' ); axis tight; 

% plot the value of the standard errors of the sample autocorrelations: 
n = n_all; se_r_k = sqrt(1/n);

figure(sf); hold on; plot( 1:(n_k), -2.0 * se_r_k*ones(1,n_k), 'r-' ); 
figure(sf); hold on; plot( 1:(n_k), +2.0 * se_r_k*ones(1,n_k), 'r-' ); axis( [0.25, 20, -0.4, 1.0 ] );
saveas( gcf, '../../WriteUp/Graphics/Chapter5/yield_data_spacf', 'epsc' );