% 
% Generates ARIMA models for the demand for repair parts 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/' ); 

Y_all = load_demand_repair_data(); n_all = length(Y_all); 

n = n_all - 4; Y = Y_all(1:n); % hold out some data

fh = figure; p_data = plot( 1:n, Y, '-bx' ); axis tight; hold on; 
xlabel( 'Time' ); ylabel( 'montly demand for repair parts' );
saveas( gcf, '../../WriteUp/Graphics/Chapter5/demand_data_plt', 'epsc' ); 

% plot the log of this data set: 
Y = log(Y); 
fhl = figure; p_data = plot( 1:n, Y, '-bx' ); axis tight; hold on; 
xlabel( 'Time' ); ylabel( 'montly demand for repair parts' );
saveas( gcf, '../../WriteUp/Graphics/Chapter5/log_demand_data_plt', 'epsc' ); 

% plot the sample autocorrelation function of the non-statonary series 
plot_SACF( Y ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter5/log_demand_sacf', 'epsc' ); 

% observe that this series is non-stationary (has a changing mean) ... plot the first difference 
% -- compute the sample variance of the first few differences
fprintf('std(Y)= %10.6f, std(diff(Y))= %10.6f, std(diff(Y,2))= %10.6f, std(diff(Y,3))= %10.6f\n',std(Y),std(diff(Y)),std(diff(Y,2)),std(diff(Y,3)))

nd = n-1; Yd = diff(Y); 

fhd = figure; p_data = plot( 1:nd, Yd, '-bx' ); axis tight; hold on; 
xlabel( 'Time' ); ylabel( 'The first difference in the log of the demand for repair parts' );
saveas( gcf, '../../WriteUp/Graphics/Chapter5/diff_log_demand_plt', 'epsc' ); 

% plot the sample partial autocorrelation function
r_k = plot_SACF( Yd ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter5/diff_log_demand_plt_sacf', 'epsc' ); 

% plot the sample partial autocorrelation function
plot_SPACF(r_k(2:end),nd);
saveas( gcf, '../../WriteUp/Graphics/Chapter5/diff_log_demand_plt_spacf', 'epsc' ); 

% estimate the coefficients of the model (in R) ... see the file model_demand_section_5_8_3.R

% compute the one-step-ahead predicitions, z_n(1), and their residuals using the estimated 
% model coefficients
% 
theta = 0.55;
sigma = sqrt(0.024);

a = zeros(1,nd+1);
% now 
% a(1) is a_{1} is set to zero 
% a(2) is a_{2} is the computed residual 
% ...
% a(123) is a_{123} (the last residual)
%

znof1 = zeros(1,nd); 
% now
% znof1(1) is z_{1}(1) 
% znof1(2) is z_{2}(1) 
% ....
% znof1(121) is z_{121}(1)
% znof1(122) is z_{122}(1) ... last element filled ... predicts the last value z_{123}
%

znof1(1) = Y(1);              % <- a_1 (the first residual) is assumed zero 
a(2)     = Y(2) - znof1(1);   % <- a_2 (the second residual) 
for ii=2:nd,
 znof1(ii) = Y(ii) - theta * a(ii);
 a(ii+1)   = Y(ii+1) - znof1(ii); 
end

% Note: 
% exp(znof1(1)) ... should be compared with Y(2)
% exp(znof1(2)) ... should be compared with Y(3) ... etc
% 
figure(fh); plot( 2:n, exp(znof1(1:end)), '-rx' ); 
title( 'Raw Demand Data (in blue) and MA(1) predictions (in red)' );
saveas( gcf, '../../WriteUp/Graphics/Chapter5/log_demand_data_plt', 'epsc' ); 

% plot the sample autocorrelation function of the residuals 
plot_SACF( a );
saveas( gcf, '../../WriteUp/Graphics/Chapter5/log_demand_data_res_SACF', 'epsc' ); 

% compute the Portmanteau statistic Q ... skipped for time 

% compute three forecasts (and their variances) ... skipped for time