%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example: Car_Sales
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%% Import Car_Sales data to matlab

load Car_Sales;

%%%%  Compute the sample size of data

n=length(Car_Sales);


%%% Plot a figure for original data

plot(1960:1/12:1969-1/12, Car_Sales);

title('Monthly car sales in Quebec; January 1960 to December 1968');

xlabel('Year');

ylabel('Car Sales (1000s)');

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

n=n-12;

z=Car_Sales(1:n)';


%%%%  Construct matrix X using indication function 

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

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

%%% compute the number of parameters

p=11+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-14));


%%%  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 estimate of beta, its standard deviation 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:12
    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 12 MSR; SSE 83 MSE; STTO 95 0]

%%% Show the value of autocorrelations of residuals on the screen

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

%%%% compute forecast error for last 12 data

zl=Car_Sales(n+1:n+12);

z_forecast=zeros(1,12);

for i=1:12
    xi=[1 i+n zeros(1,11)];
    
    if i~=12
       xi(2+i)=1;
    end
    
    z_forecast(i)=xi*beta;
end

forecast_error=zl-z_forecast;

%%%% Show last 12 data, its forecast and the forecast error on the screen

[zl; z_forecast; forecast_error]
   



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Another look at Car Sales using Trignonmetric function to modeling
%%% Seasonal component

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%  Construct matrix X using indication function 

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

t=(1:n)';

for i=1:5
    
  Xt=[Xt sin(2*pi*i/12*t) cos(2*pi*i/12*t)];
 
end

Xt=[Xt cos(pi*t)];

beta_t=inv(Xt'*Xt)*Xt'*z;

%%% computing residuals for model and the estimate of sigma
 
Error_t=z-Xt*beta_t;

%%%% cmputing sum of squares for every regression component

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

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

SSR_trend=(n-1)*var(Xt(:,2)*beta_t(2));

SSR_1=(n-1)*var(Xt(:,3:4)*beta_t(3:4));

SSR_2=(n-1)*var(Xt(:,5:6)*beta_t(5:6));

SSR_3=(n-1)*var(Xt(:,7:8)*beta_t(7:8));

SSR_4=(n-1)*var(Xt(:,9:10)*beta_t(9:10));

SSR_5=(n-1)*var(Xt(:,11:12)*beta_t(11:12));

SSR_6=(n-1)*var(Xt(:,13)*beta_t(13));
 
SSE_t=Error'*Error;

%%%% Show every sum square of every regression component and its degree of freedom on the screen 

[SSR_t 12; SSR_trend 1; SSR_1 2; SSR_2 2; SSR_3 2; SSR_4 2; SSR_5 2; SSR_6 1; SSE_t 83; STTO_t 95]


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Modeling Data using first two hamonic 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%  Construct matrix X using indication function 

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

t=(1:n)';

for i=1:2
    
  Xt=[Xt sin(2*pi*i/12*t) cos(2*pi*i/12*t)];
 
end

beta_t=inv(Xt'*Xt)*Xt'*z;

%%% computing residuals for model and the estimate of sigma
 
Error_t=z-Xt*beta_t;

sigma_t=sqrt(sum(Error_t.^2)/(n-6));


%%%  Compute standard deviation of the estimate of beta

beta_error_t=sigma_t*sqrt(diag(inv(Xt'*Xt)));


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

t_ratio_t=beta_t./beta_error_t;

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

[beta_t beta_error_t t_ratio_t]



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

for i=1:12
    r_t(i)=auto_corr(Error_t,i);
end

%%% show the value of autocorrelation of residuals for the model  using
%%% first two hamonic functions

[1:12; r_t(1:12)]


%%%% compute forecast error for last 12 data

zl=Car_Sales(n+1:n+12);

z_forecast_t=zeros(1,12);

for t=1:12
    xi=[1 t+n sin(2*pi*1/12*(n+t)) cos(2*pi*1/12*(n+t)) sin(2*pi*2/12*(n+t)) cos(2*pi*2/12*(n+t))];    
    z_forecast_t(t)=xi*beta_t;
end

forecast_error_t=zl-z_forecast_t;

%%% show last 12 data, their forecast value and forecast error on screen 

[zl; z_forecast_t; forecast_error_t]




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

alpha=0.01:0.01:0.20;

m=length(alpha);

sse=zeros(m,1);

r_error=zeros(m,6);

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
   
    [sse(i), r_error(i,:), error_mean(i), error_std(i)]=SSE2(z,alpha(i));
   
end


%%% Show the table for SSE and alpha 

[alpha' sse]

%%% Plot a figure for SSE with alpha

figure

plot(alpha,sse);

title({'Plot of SSE(\alpha) from locally constant seasonal moldes';'using trigonometric functions for Example Car Sales'});

xlabel('\alpha');

ylabel('SSE(\alpha)');