% Compute predict error by Canadian data using an AR(2) model X(t)=a+b*X(t-1)+c*X(t-2)+error
% Parameters, a, b,c, are taken from the AR(2) model fit to the Canadian
% lynx dataset from t=1 to t=82. Then compute the predict error using
% Canadian data from t=83 to t=114. 



load lynx.dat % import lynx data set
log_lynx=log10(lynx); % make log10 transformation for dataset
subplot(3,1,1); % decompose the figure to two sub-figure, and active the first sub-figure
plot(log_lynx(83:114),'.-');
title('Number of trap lynx'); % make a figure for log-lynx data

% Using AR(2) model to fit lynx data


% Contruct linear model form, Y=X(t), X=[1 X(t-1) X(t-2)], then
% Y=X*beta+error, beta=[a b c]'  Using dataset from X(1) to X(82). Other
% data from X(83) to X(114) will be used for check prediction error of
% AR(2)

Y=log_lynx;
X0=ones(80,1);
X1=Y(2:81);
X2=Y(1:80);
Y=Y(3:82);
X=[X0 X1 X2];

% Estiamte beta by formula  beta=(X'X)^{-1}X'Y 

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

% Estimate  error's standard deviation using formula
% sigma^2= ||Y-X*beta||^2/(n-p), here n=112, p=3.

sigma=sqrt((Y-X*beta)'*(Y-X*beta)/(112-3));

% predicte data from t=83 to t=114 using the estimate from lynx data (t=1 to t=82) fitting, the
% first two data use same as the original lynx data X(81),X(82). 
pred_data=zeros(114,1);

pred1_data=zeros(114,1);

pred_data(81:82)=log_lynx(81:82);

pred1_data(81:82)=log_lynx(81:82);

% repeated to produce predition X(t) from t=83 to t=114.

% method 1, using the observation data to X(t-1) to predict X(t)
N=114;
for i=83:114
    pred_data(i)=beta(1)+beta(2)*log_lynx(i-1)+beta(3)*log_lynx(i-2);
end

% method 2, only using the observation data X(81) and X(82) predict the
% future data to t=114. 
for i=83:114
    pred1_data(i)=beta(1)+beta(2)*pred1_data(i-1)+beta(3)*pred1_data(i-2);
end

% check mean and median of prediction error. 

MMSE_mean=mean((pred_data(83:114)-log_lynx(83:114)).^2);
MMSE_median=median((pred_data(83:114)-log_lynx(83:114)).^2);

MMSE1_mean=mean((pred1_data(83:114)-log_lynx(83:114)).^2);
MMSE1_median=median((pred1_data(83:114)-log_lynx(83:114)).^2);



subplot(3,1,2); % active sub-figure 2
plot(pred_data(83:114),'.-'); % plot the simulated dataset
title('A predtion from an AR(2) model by method 1'); % added a title 

subplot(3,1,3); % active sub-figure 2
plot(pred1_data(83:114),'.-'); % plot the prediction dataset of method 2
title('A predtion from an AR(2) model by method 2'); % added a title 

