function M = mcmmpe_armax(D, si)
    n_alfa = 60;
    n_beta = 60;
    na = si(1);
    nb = si(2);
    nc = si(3);
    nk = si(4);

    N=250;
    model_arx = arx(D, [n_alfa n_beta 1]);

    epsilon_tilda = pe(model_arx, D);
    epsilon_tilda = epsilon_tilda.OutputData;
    y = D.OutputData;
    u = D.InputData;

    phi = zeros(N,na+nb+nc);
    for n=1:N
        for i=1:na
            if (n-i)<= 0
                phi(n,i) = 0;
            else
                phi(n,i) = -y(n-i);
            end
        end

        for i = 1:nb
            if (n-i)<=0
                phi(n,i+na)=0;
            else
                phi(n,i+na) = u(n-i);
            end
        end

        for i = 1: nc
            if (n-i)<=0
                phi(n,i+na+nb)= 0;
            else
                phi(n,i+na+nb)= epsilon_tilda(n-i);
            end
        end
    end

    a = (1/N)*(phi'*phi);
    b = (1/N)*(phi'*y);

    theta_hat = a\b

    A_estimat = [1 theta_hat(1) theta_hat(2)];
    B_estimat = [0 theta_hat(3) theta_hat(4)];
    C_estimat = [1 theta_hat(5) theta_hat(6)];

    M = idpoly(A_estimat,B_estimat,C_estimat);
end



























clc

N = 250;
A = [1 -1.5 0.7];
B = [0 1 0.5];
C = [1 -1 0.2];

u_id = sign(randn(N,1));
e_id = randn(N,1);
y_id = filter(B,A,u_id) + filter(C,A,e_id);

u_va = sign(randn(N,1));
e_va = randn(N,1);
y_va = filter(B,A,u_va) + filter(C,A,e_va);


Did = iddata(y_id, u_id);
Dva = iddata(y_va, u_va);

si = [2 2 2 1];
M = mcmmpe_armax(Did, si)

eps_id = pe(M,Did);
eps_id = eps_id.OutputData;
disp_id = var(eps_id);

eps_va = pe(M,Dva);
eps_va = eps_va.OutputData;
disp_va = var(eps_va);

E_id = 100 * (1 - sqrt( sum(eps_id.^2) / sum( (y_id - sum(y_id)/N).^2)));
E_va = 100 * (1 - sqrt( sum(eps_va.^2) / sum( (y_va - sum(y_va)/N).^2)));


[y_sim_id,~, ~] = compare(Did, M);
y_sim_id = y_sim_id.OutputData;

[y_sim_va,~, ~] = compare(Dva, M);
y_sim_va = y_sim_va.OutputData;

figure;
hold on;
plot(y_id, 'b', 'DisplayName', 'Iesire masurata');
plot(y_sim_id, 'r', 'DisplayName', 'Iesire simulata');
hold off;

title('Set identificare: Iesire masurata si simulata');
xlabel('Timp');
ylabel('Amplitudine');
legend show;
grid on;
text(100,13, sprintf('Fitness: %.2f%%', E_id), 'FontSize', 10, 'Color', 'k', 'BackgroundColor', 'w');

figure;
hold on;
plot(y_va, 'b', 'DisplayName', 'Iesire masurata');
plot(y_sim_va, 'r', 'DisplayName', 'Iesire simulata');
hold off;

title('Set validare: Iesire masurata si simulata');
xlabel('Timp');
ylabel('Amplitudine');
legend show;
grid on;
text(100,13, sprintf('Fitness: %.2f%%', E_va), 'FontSize', 10, 'Color', 'k', 'BackgroundColor', 'w');


figure;
plot(eps_id, 'b');
title('Set identificare: Eroare de predictie');
xlabel('Timp');
ylabel('Amplitudine');
grid on;
text(200,2, sprintf('Dispersia: %.2f%%', disp_id), 'FontSize', 10, 'Color', 'k', 'BackgroundColor', 'w');

figure;
plot(eps_va, 'b');
title('Set validare: Eroare de predictie');
xlabel('Timp');
ylabel('Amplitudine');
grid on;
text(200,2, sprintf('Dispersia: %.2f%%', disp_va), 'FontSize', 10, 'Color', 'k', 'BackgroundColor', 'w');