%
clear all; clc;

u = sqrt(2)*randn(2000,1);
e = sqrt(0.02)*randn(2000,1);

yy =zeros(2000,1);

bigloop = 30;

for bigi = 1:bigloop,
for ii = 1:2000
    if (ii == 1)
        yy(ii) = e(ii);
    elseif (ii == 2)
        yy(ii) = 1.991*yy(ii-1) + e(ii);
    elseif ii == 3
        yy(ii) = 1.991*yy(ii-1) - 2.203*yy(ii-2) + e(ii);
    elseif ii == 4
        yy(ii) = 1.991*yy(ii-1) - 2.203*yy(ii-2) + 1.841*yy(ii-3)...
            + 0.103*u(ii-3) + e(ii);
    elseif ii >= 5
        yy(ii) = 1.991*yy(ii-1) - 2.203*yy(ii-2) + 1.841*yy(ii-3)...
            - 0.894*yy(ii-4) + 0.103*u(ii-3) + 0.181*u(ii-4) + e(ii);
    end    
end
    phix = zeros(1,6);
    RR = zeros(6,6);
    FF = zeros(6,1);
for ii = 5:2000
    phix = [-yy(ii-1); -yy(ii-2); -yy(ii-3); -yy(ii-4); u(ii-3); u(ii-4)];
    RR = RR + phix*phix';
    FF = FF + phix*yy(ii);   
end
    Rp = RR/(2000);
    Fp = FF/(2000);
    theta(bigi,:) = Rp\Fp;
end; 

trueTheta = [ -1.991 2.203 -1.841 0.894 0.103 0.181]';
meanTheta = sum(theta,1)/bigloop;
Bias = norm(trueTheta-meanTheta');
disp('Bias = '); disp(Bias) ;
Varx = zeros(1,6);
for ii = 1:bigloop,
    Varx = Varx + theta(ii,:)-meanTheta;
end
Varm = norm(Varx.^2/bigloop);
disp('Variance = '); disp(Varm);

