% Solution of Lab session on Time Series Prediction
clear; close all;

N=100;
t=[1:N]';
C=[1 0.5 0];
A=[1 -1.8 0.8075];
sigmae2=1;

% 1. Covariance
myARMA=idpoly(A,[],C);
g=impulse(myARMA,1000);
nlag=20;
Ry=sigmae2*xcorr(g,nlag);

% 2. Simulation
e=iddata([],sqrt(sigmae2)*randn(size(t)));
y=sim(myARMA,e);

% 3. Sample covariance
yt=y.y-mean(y.y);
Rys=xcorr(yt,nlag,'unbiased');
taus=[-nlag:1:nlag];
figure, plot(taus,Ry,'*',taus,Rys,'+r')
xlabel('\tau'); ylabel('R_y'); title('Covariance of y(t): true (blue), sample (red)')
% try with higher N!

% 4. Wiener predictor
step=1;
ypred=predict(myARMA,y,step);
MSE1=sum((ypred.y-y.y).^2)/(N)
q1=deconv([C zeros(1,step-1)],A);
trueMSE1=sigmae2*sum(q1.^2)
figure, plot(t,y.y,'b',t,ypred.y,'g')
xlabel('time'); ylabel('y(t)'); title('1-step predictor, y(t): true (blue), predicted (green)')

step=3;
ypred3=predict(myARMA,y,step);
MSE3=sum((ypred3.y-y.y).^2)/(N)
q3=deconv([C zeros(1,step-1)],A);
trueMSE3=sigmae2*sum(q3.^2)
figure, plot(t,y.y,'b',t,ypred3.y,'g')
xlabel('time'); ylabel('y(t)'); title('3-step predictor, y(t): true (blue), predicted (green)')

% 5. Example with unstable zero
Cnew=[1 2 0];
A=[1 -1.8 0.8075];
myARMAnew=idpoly(A,[],Cnew);
ynew=sim(myARMAnew,e);

G=tf(Cnew,A,1);
PhiY=G*G'*sigmae2;
[Hsf,sigmae2sf]=spectralfact(PhiY);
[Csf,Asf]=tfdata(Hsf,'v');
myARMAnew_sf=idpoly(Asf,[],Csf);

step=1;
yprednew=predict(myARMAnew_sf,ynew,step);
MSE1=sum((yprednew.y-ynew.y).^2)/(N)
q1=deconv([Csf zeros(1,step-1)],A);
trueMSE1=sigmae2sf*sum(q1.^2)
figure, plot(t,ynew.y,'b',t,yprednew.y,'g')
xlabel('time'); ylabel('y(t)'); title('1-step predictor, y(t): true (blue), predicted (green)')

figure, plot(t,ypred.y-y.y,'b',t,yprednew.y-ynew.y,'r')
xlabel('time'); ylabel('error'); title('Comparison of 1-step prediction error in items 4 (blue) and 5 (red)')

