function [] = mcmmp_oe2(A, B, N, nr)

lambda = 1;

if(nargin < 4)
    nr = 100;
end

u = sign(randn(N, nr));
e = randn(N, nr) * lambda^2;
yi = filter (B,A, u);
v = e;
y = yi + v;
y2 = e - y;

teta = zeros(4, nr);

for i=1:nr

    [ry2, k] = xcov(y2(:, i),'biased');
    ry2 = ry2(k>=0);

    [ru, k] = xcov(u(:, i),'biased');
    ru = ru(k>=0);

    [ryu, k] = xcov(y(:, i), u(:,i), 'biased');
    ryu = ryu(k>=1);

    [ry2u, k] = xcov(y2(:, i), u(:,i), 'biased');
    ry2u = ry2u(k>=0);

    [ruy2, k] = xcov(u(:, i), y2(:,i), 'biased');
    ruy2 = ruy2(k>=0);

    [ryy2, k] = xcov(y(:, i), y2(:,i), 'biased');
    ryy2 = ryy2(k>=-1);

    R = [ry2(1) ry2(2) ry2u(1) ry2u(2);
        ry2(2) ry2(1) ruy2(2) ry2u(1);
        ruy2(1) ruy2(2) ru(1) ru(2);
        ry2u(2) ruy2(1) ru(2) ru(1)];

    r = [ryy2(3);
         ryy2(4);
         ryu(1);
         ryu(2)];

    t = R \ r;
    teta(:, i) = t;

    eroare = y(:,i) - [[0; y2(1:(N-1),i)] [0; 0; y2(1:(N-2),i)] [0; u(1:(N-1),i)] [0; 0; u(1:(N-2), i)]]*t;
    deviatiastandard(i, 1) = norm(eroare)/sqrt(N-size(teta,1));
    dispersia(i,1) = deviatiastandard(i)^2;
end

a1_m = mean(teta(1,:));
a2_m = mean(teta(2,:));
b1_m = mean(teta(3,:));
b2_m = mean(teta(4,:));

figure
plot(dispersia)
hold on
plot(zeros(nr, 1) + lambda^2);
hold off
title('Model OE[2, 2]');
xlabel('Numarul realizarii');
ylabel('Varianta');
legend('zgomot estimat', 'zgomot real');

x_limits = xlim;
y_limits = ylim;

x_text = mean(x_limits);
y_text1 = y_limits(2) - 0.05 * (y_limits(2) - y_limits(1));
y_text2 = y_limits(2) - 0.1 * (y_limits(2) - y_limits(1));

text(x_text, y_text1, ['Parametri adevarati: ', sprintf('%9.4f',[A(2) A(3) B(2) B(3)])], ...
    'HorizontalAlignment', 'center', 'VerticalAlignment', 'top', 'FontWeight', 'bold');
text(x_text, y_text2, ['Parametri estimati: ', sprintf('%9.4f',[a1_m a2_m b1_m b2_m])], ...
    'HorizontalAlignment', 'center', 'VerticalAlignment', 'top', 'FontWeight', 'bold');

end