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



load lynx.dat % import lynx data set
log_lynx=log10(lynx); % make log10 transformation for dataset
subplot(2,1,1); % decompose the figure to two sub-figure, and active the first sub-figure
plot(log_lynx,'.-');
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]'

Y=log_lynx;
X0=ones(112,1);
X1=Y(2:113);
X2=Y(1:112);
Y=Y(3:114);
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));

% Simulate a new time series using the estimate from lynx data fitting, the
% first two data use same as the original lynx data 

sim_data(1:2)=log_lynx(1:2);

% repeated to produce X(t) from t=3 to t=114
N=114;
for i=3:114
    sim_data(i)=beta(1)+beta(2)*sim_data(i-1)+beta(3)*sim_data(i-2)+ normrnd(0,sigma,[1 1]);
end


subplot(2,1,2); % active sub-figure 2
plot(sim_data,'.-'); % plot the simulated dataset
title('A realized from an AR(2) model'); % added a title 

