clear; close all;

%% 3. Simulation
N=1000;
t=[1:N]';
C=[1 -3/4 0];
A=[1 -3/2 3/4];
sigmae2=8;
myARMA=idpoly(A,[],C);
e=iddata([],sqrt(sigmae2)*randn(size(t)));
y=sim(myARMA,e);

%% 4. Wiener predictor
step=2;
ypred=predict(myARMA,y,step);
sampleMSE=sum((ypred.y-y.y).^2)/(N)    % sample MSE error
q=deconv([C zeros(1,step-1)],A);
trueMSE=sigmae2*sum(q.^2)              % true MSE error
figure, plot(t,y.y,'b',t,ypred.y,'g')

%% 5. Numerical estimate of MSE
Nreal=20;
for h=1:Nreal,
    e=iddata([],sqrt(sigmae2)*randn(size(t)));
    y=sim(myARMA,e);
    ypred=predict(myARMA,y,step);
    sampleMSE(h)=sum((ypred.y-y.y).^2)/(N);
end
estMSE=mean(sampleMSE)

