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

lambda = 1;

if(nargin < 4)
    nr = 100;
end

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

teta = zeros(4, nr);

for i=1:nr

    [ry, k] = xcov(y(:, i),'biased');
    ry = ry(k>=0);

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

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

    [ruy, k] = xcov(u(:,i), y(:, i), 'biased');
    ruy = ruy(k>=-2);

    R = [ry(1) ry(2) -ryu(2) -ryu(3);
        ry(2) ry(1) -ryu(1) -ryu(2);
        -ruy(3) -ruy(4) ru(1) ru(2);
        -ruy(2) -ruy(3) ru(2) ru(1)];

    r = [-ry(2);
         -ry(3);
          ruy(2);
          ruy(1)];

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

    eroare = y(:,i) - [[0; -1*y(1:(N-1),i)] [0; 0; -1*y(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 ARX[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