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

[T,Y] = load_series_1(); n = size(Y,1); alpha = 0.05; 

% make the requested plots: 
% 
fh=figure; plot( T, Y, 'og' ); title( 'Iowa nonfarm income' ); grid on; %axis tight; 
saveas( gcf, '../../WriteUp/Graphics/Chapter2/prob_2_19_y', 'epsc' ); 

lfh=figure; plot( T, log(Y), 'og' ); title( 'Log(Iowa nonfarm income)' ); grid on; axis tight; 
saveas( gcf, '../../WriteUp/Graphics/Chapter2/prob_2_19_ly', 'epsc' ); 

% fit the regression model \log(y_t) = \beta_0 + \beta_1 t :
% 
X = [ ones(n,1), T ];  % <- create the augmented matrix we will regress on 
[beta,sC,s_hat,t_stats,SSE,SSR,SSTO,R2,MSR,MSE,F_ratio,MI] = wwx_regression( X, log(Y), 1 ); p=1; 

ly_hat = X * beta; 

figure(fh); hold on; 
plot( T, exp(ly_hat), 'rx' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter2/prob_2_19_y', 'epsc' ); 

figure(lfh); hold on; 
plot( T, ly_hat, 'rx' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter2/prob_2_19_ly', 'epsc' ); 

% now consider the residuals (in log space): 
e = log(Y) - ly_hat; 
figure; plot( T, e, 'x' ); title( 'residuals' ); 
xlabel( 'time' ); ylabel( 'Ln(y) - hat(Ln(y))' ); axis tight; 
saveas( gcf, '../../WriteUp/Graphics/Chapter2/prob_2_19_res', 'epsc' ); 

% calculate the correlation of these residuals: 
n_autocorr = round(0.2*length(e));
if( 1 ) 
  if( exist('/home/weatherw/Trash/Poularikas/AdaptiveFilteringCode') ) addpath('/home/weatherw/Trash/Poularikas/AdaptiveFilteringCode'); end
  %aasamplebiasedautoc(e.',10)
  e2 = e - mean(e); e2 = e2/std(e2); 
  r_k = aasampleunbiasedautoc(e2.',n_autocorr);
else
  r_k = autocor(e,n_autocorr);
end

sf=figure; sh=stem( 0:(n_autocorr-1), r_k, 'x' ); 
xlabel( 'lag (k)' ); ylabel( 'sample based autocorrelation r_k' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter2/prob_2_19_r_k', 'epsc' ); 

% plot the value of the standard errors of the sample autocorrelations: 
se_r_k = 1.96 / sqrt(n);
fprintf('two sigma confidence bounds on the magintude of sample autocorrelations = %10.6f\n', se_r_k);

figure(sf); hold on; plot( 0:(n_autocorr-1), se_r_k*ones(n_autocorr,1), 'r-' ); axis tight; 
saveas( gcf, '../../WriteUp/Graphics/Chapter2/prob_2_19_r_k', 'epsc' ); 

% predict the next for quarters and their confidence intervals: 
ly_predicted = []; y_predicted = [];
ly_ci        = []; y_ci        = []; 
for qi=1:4, 
  t_extrap  = T(end) + qi; 
  x         = [ 1; t_extrap ];
  ly_extrap = (x.') * beta; 
  y_extrap  = exp(ly_extrap);  
  v         = tinv( 1 - 0.5*alpha, n-p-1 ) * s_hat * sqrt( 1 + (x.') * MI * x ); 
  
  ly_predicted = [ ly_predicted; ly_extrap ]; y_predicted  = [ y_predicted;   y_extrap ];     % <- save this information
  ly_ci        = [ ly_ci ; ly_extrap-v, ly_extrap+v ]; y_ci = [ y_ci ; exp(ly_extrap-v), exp(ly_extrap+v) ]; 
  fprintf('quarter = %5d; log ci = [%6.2f, %6.2f, %6.2f]; direct ci = [%10.2f, %10.2f, %10.2f]\n',...
          qi, ly_extrap-v,ly_extrap,ly_extrap+v, exp(ly_extrap-v),y_extrap,exp(ly_extrap+v));
end

figure(fh); plot( T(end)+[1:4], y_predicted, 'k.', 'markersize', 8 ); 
plot( T(end)+[1:4], y_ci(:,1), 'k.', 'markersize', 8 ); 
plot( T(end)+[1:4], y_ci(:,2), 'k.', 'markersize', 8 ); axis tight; 
saveas( gcf, '../../WriteUp/Graphics/Chapter2/prob_2_19_y', 'epsc' ); 

figure(lfh); plot( T(end)+[1:4], ly_predicted, 'k.', 'markersize', 8 ); 
plot( T(end)+[1:4], ly_ci(:,1), 'k.', 'markersize', 8 ); 
plot( T(end)+[1:4], ly_ci(:,2), 'k.', 'markersize', 8 ); axis tight; 
saveas( gcf, '../../WriteUp/Graphics/Chapter2/prob_2_19_ly', 'epsc' );