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

phi_hat   = 0.4;
sigma_hat = 0.18; 
z_49      = 33.4; z_nm1 = z_49; % n minus 1 
z_50      = 33.9; z_n   = z_50; % n 

% Part (a):
%

z_l = []
for ll=1:5,
  z_np1 = (1+phi_hat)*z_n - phi_hat*z_nm1; % n plus 1 
  z_l(ll) = z_np1; 
  fprintf('l=%5d; z_np1= %10.6f\n',ll,z_np1); 
  z_nm1 = z_n;
  z_n   = z_np1; 
end

percent    = 0.95;
alpha      = 1 - percent; 
u_lambda_2 = norminv( 1-alpha/2, 0, 1 );   % <- the lambda/2 "percent" point of the standard normal 

% compute the confidence intervals around the above points: 
%
phi_l = cumsum( phi_hat.^( 0:4 ) );        % <- the values of $\psi_l$ for $l=0,1,2,3,4$

V_l = ( sigma_hat^2 ) * cumsum( phi_l.^2 ) % <- the values of the variances of the l-step-ahead look aheads $l=1,2,\cdots,5$

ci_l = [ 1:5; z_l - u_lambda_2 * sqrt( V_l ) ; z_l + u_lambda_2 * sqrt( V_l ) ] ; 

% Part (b):
%

z_51 = 34.2; % <- the new measurement 

% the new predictions are given by:
%
z_lp1 = z_l(2:end) + phi_l(2:end) * ( z_51 - z_l(1) ); 

phi_l = phi_l(1:end-1); % <- we now only predict for 1 <= l <= 4

V_l = ( sigma_hat^2 ) * cumsum( phi_l.^2 ) % <- the values of the variances of the l-step-ahead look aheads $l=1,2,\cdots,5$

ci_lp1 = [ 1:4; z_lp1 - u_lambda_2 * sqrt( V_l ) ; z_lp1 + u_lambda_2 * sqrt( V_l ) ] ; 

figure; 

plot( 2:5, ci_l(2,2:end), 'x-r' ); hold on; plot( 2:5, ci_l(3,2:end), 'x-r' ); 
plot( 2:5, ci_lp1(2,:), 'x-g' ); plot( 2:5, ci_lp1(3,:), 'x-g' ); 
xlabel( 'l (look ahead index)' ); ylabel( 'confidence interval for z_{51}(l)' ); 
saveas(gcf, '../../WriteUp/Graphics/Chapter5/before_and_after_ci', 'epsc' ); 

return;