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

lambda = 1;

% Verificarea corectitudinii parametrilor de apel ai rutinei
if(nargin < 4)
    nr = 100;
end

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

teta = zeros(4, nr);

% Aplicare metoda MCMMP
for i=1:nr

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

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

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

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

    % Cacul matrice de covarianta a datelor (R)
    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)];

    % Calcul vector de covarianta a datelor (r)
    r = [-ry(2);
         -ry(3);
          ruy(2);
          ruy(1)];

    % Estimare parametrii pentru realizarea curenta
    t = R \ r;
    teta(:, i) = t;

    % Estimare zgomot alb si dispersie
    % La eroare avem matrice Toeplitz
    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

% Calcul medie parametrii estimati
a1_m = mean(teta(1,:));
a2_m = mean(teta(2,:));
b1_m = mean(teta(3,:));
b2_m = mean(teta(4,:));

% Afisare grafic
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');
text(2, 1.2, ['Parametri adevarati a2 a3 b2 b3 : ', sprintf('%9.4f',[A(2) A(3) B(2) B(3)])]);
text(2, 0.8, ['Parametri estimati  a1m a2m b1m b2m: ', sprintf('%9.4f',[a1_m a2_m b1_m b2_m])]);

end