%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example: New Plant and Equipment Expenditure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%% Import  the data for Example New Plant and Equipment Expenditures  to matlab

load New_Plant;

%%%%  Using first 44 data to estimate, the other 8 data are used for
%%%%  computer forecast error and check model

n=44;

%%%% plot a figure for original data 

subplot(1,2,1);

plot(New_Plant(1:n));

title({'Quarterly new plant and equipment expenditure in U.S. industries';'(in billions of dollars),first quarter 1964 to fourth quarter1976'});

xlabel('Year');

ylabel('New plant equipment expenditure (billons $)');

%%%% Using first n-12 data to modeling data, the other 12 data used to
%%%% check forecast error

%%%%  change data to z_t, get the ordinary least square estimate for beta

z=log(New_Plant(1:n))';

%%%%% plot the new data after logarithm transformation figures from
%%%%% original data 

subplot(1,2,2);

plot(z);

title({'Logarithm of quarterly new plant and equipment expenditure in U.S. industries';'(in billions of dollars, first quarter 1964 to fourth quarter1976'});

xlabel('Year');

ylabel('Logarithm of New plant equipment expenditure (billons $)');



%%%%  Construct matrix X using indication function 

X=[ones(n,1) (1:n)'];

for i=1:3
    
  X=[X (mod((1:n)',4)==i) ];
 
end

%%% compute the number of parameters

p=3+2;

%%%  get the least square estimate of beta

beta=inv(X'*X)*X'*z;

%%% computing residuals for model and the estimate of sigma
 
Error=z-X*beta;

sigma=sqrt(sum(Error.^2)/(n-p));


%%%  Compute standard deviation of the estimate of beta

beta_error=sigma*sqrt(diag(inv(X'*X)));

%%% compute the value of t-statistcs for beta

t_ratio=beta./beta_error;

%%% Show the estimate of beta, its standard deivation and t-statistics on
%%% the screen 

[beta beta_error t_ratio]


%%%%  Compute the auto-correlation of residuals from lag=1 to 12

for i=1:8
    r(i)=auto_corr(Error,i);
end

%%%  compute SSTO, SSR, SSE, MSR, MSE to construct ANOVA table

STTO=z'*z-n*(mean(z))^2;

SSR=beta'*X'*z-n*(mean(z))^2;

SSE=Error'*Error;

MSR=SSR/(p-1);

MSE=SSE/(n-p);

%%% Show ANOVA table on the screen

[SSR 4 MSR; SSE 39 MSE; STTO 43 0]

%%% Show the table for autocorrelation of residuals on the screen

[1:8; r(1:8)]

%%% Plot a figue for residuals

figure

stem(Error)

title({'Time series plot of the residuals from the model'})

xlabel('Time');

ylabel('Residulas e_t');




%%%% Compute SSE from alpha=0.01 to 0.2

alpha=0.1:0.1:1.5;

m=length(alpha);

sse=zeros(m,1);

r_error=zeros(m,8);

error_mean=zeros(m,1);

error_std=zeros(m,1);


% sse: SSE(alpha), 
% t_error: auto_correlation of one step forecast error,
% error_mean: mean of forecast error. 
% error_std: standard deviation of the estimate of the mean of forecast error

for i=1:m
   
    if alpha(i)~=1
      
        [sse(i), r_error(i,:), error_mean(i), error_std(i)]=SSE3(z,alpha(i),New_Plant);
    
    else
        
        [sse(i), r_error(i,:), error_mean(i), error_std(i)]=SSE4(z,alpha(i),New_Plant);
     
    end
end

%%% Plot a figure to show SSE with alpha, and select alpha which minimize
%%% SSE. 

%%% Show values of alpha and SSE(alpha) on the screen by table and graph

[alpha' sse]

figure

plot(alpha,sse);

title({'Plot of SSE(\alpha) from locally constant seasonal moldes using indicator functions';'for Example New Plant and Equipment Expenditure'});

xlabel('\alpha');

ylabel('SSE(\alpha)');