How to simulate a multivariable autoregressive model forecast

Below is the example for arx command, followed by an estimate of the response. How to modify the example and forecast to represent a system with two outputs with no inputs (two coupled time series) for armax?

%%MO Example 
clear; close all; clc;
A0 = {[1  -1.5  0.7 ], [0 -.5]
    [0 -.04 .01], [1  -1.4  .6 .001]};
A=A0;
m0 = idpoly(A);
e = iddata([],randn(300,size(A,1)),'ts',1);
y = sim(m0,e);
plot(y);
data = y;
na = [2 1; 2 3];
nb = 0*na;
nk = 0*na+1;
useARX = false;
if (useARX)
    sys = arx(data,na);
else
    nc = [2; 3];
    sys = armax(data,[na nc]);
end
t = get(data,'SamplingInstants');
%%Plot the predicted results
nToPredict = 5;
nToEval = 40;
[yc,fit,x0] = compare(sys,data(1:nToEval),nToPredict);

figure;
plot(t(1:nToEval),y.y(1:nToEval,1),'k-', ...
    t(1:nToEval),y.y(1:nToEval,2),'b-', ...
    t(1:nToEval),yc.y(:,1),'k:x', ...
    t(1:nToEval),yc.y(:,2),'b:x');
hold all;
yp = predict(sys,data(1:40),nToPredict);
tp = get(yp,'SamplingInstants');
plot(tp,yp.y(:,1),'ko', ...
    tp,yp.y(:,2),'bo')

%%Estimate the forecast
% sys = Discrete-time AR model:                                      
%   Model for output "y1": A(z)y_1(t) = - A_i(z)y_i(t) + e_1(t)
%     A(z) = 1 - 1.475 z^-1 + 0.6946 z^-2                      
%                                                              
%     A_2(z) = -0.5373 z^-1                                    
%                                                              
%   Model for output "y2": A(z)y_2(t) = - A_i(z)y_i(t) + e_2(t)
%     A(z) = 1 - 1.428 z^-1 + 0.6401 z^-2 + 0.06096 z^-3       
%                                                              
%     A_1(z) = -0.108 z^-1 + 0.05012 z^-2                      
% Sample time: 1 seconds
% ARMAX
% sys =Discrete-time ARMA model:                                        
%   Model for output "y1": A(z)y_1(t) = - A_i(z)y_i(t) + C(z)e_1(t)
%     A(z) = 1 - 1.487 z^-1 + 0.6851 z^-2                          
%                                                                  
%     A_2(z) = -0.4975 z^-1                                        
%                                                                  
%     C(z) = 1 + 0.07709 z^-1 - 0.05804 z^-2                       
%                                                                  
%   Model for output "y2": A(z)y_2(t) = - A_i(z)y_i(t) + C(z)e_2(t)
%     A(z) = 1 - 2.417 z^-1 + 2.053 z^-2 - 0.6112 z^-3             
%                                                                  
%     A_1(z) = -0.041 z^-1 + 0.03667 z^-2                          
%                                                                  
%     C(z) = 1 - 1.027 z^-1 - 0.01797 z^-2 + 0.1476 z^-3           
%                                                                  
% Sample time: 1 seconds

[~,ny]=size(data);
est_y = nan*yp.y;
max_na = max(na(:));
for it=1:(length(tp)-nToPredict+1)
    if (it > max_na )
        %%Save previous y
        prev_y = data.y(it-(1:max_na),:);
        est_y_prev = nan(1,ny);
        for iPred = 1:nToPredict
            for i_out = 1:ny
                azy = 0;
                for j_in = 1:ny
                    a_coeff = sys.A{i_out,j_in};
                    n_a_coeff = length(a_coeff);
                    azy = azy - a_coeff(2:end) * prev_y(1:(n_a_coeff-1),j_in);
                end
                est_y_prev(1,i_out) = 1/sys.A{i_out,i_out}(1)*azy ;
            end
            % Save y for next step
            prev_y(2:max_na,:) = prev_y(1:(max_na-1),:);
            prev_y(1,:) = est_y_prev;
        end
        est_y(it+(nToPredict-1),:) = est_y_prev;
    end
end

% Add the estimate to the plot
hold all;
hp=plot(tp,est_y,'gs');
set(hp,'MarkerSize',10);
legend('data',sprintf('compare %.1f%%',min(abs(fit))),'predict','estimated','location','best');

Answer 1

[John]

With assistance from Matlab support, the following code was created to compare the original data with forecast values from the functions compare and predict, and also manually calculated forecast values, for arx and armax. This illustrates how to manually calculated forecast values for ARX and ARMAX for multiple output time series data.

nToPredict = 3;
rng(101)
A = {[1  -1.5  0.7 ], [0 -.5]
    [0 -.04 .01], [1  -1.4  .6 .001]};
m0 = idpoly(A);
e = iddata([],randn(300,size(A,1)),'ts',1);
y = sim(m0,e);
ny = size(m0,1);

plot(y);

na = [2 1; 2 3];
useARX = true;
if (useARX)
    sys = arx(y, na);
else
    nc = [2; 3];
    sys = armax(y,[na nc]);
end
t = get(y,'SamplingInstants');
%%Plot the predicted results
nToEval = 40;
[yc,fit,x0] = compare(sys,y(1:nToEval),nToPredict);

%%Estimate the forecast
% sys = Discrete-time AR model:
%   Model for output "y1": A(z)y_1(t) = - A_i(z)y_i(t) + e_1(t)
%     A(z) = 1 - 1.475 z^-1 + 0.6946 z^-2
%
%     A_2(z) = -0.5373 z^-1
%
%   Model for output "y2": A(z)y_2(t) = - A_i(z)y_i(t) + e_2(t)
%     A(z) = 1 - 1.428 z^-1 + 0.6401 z^-2 + 0.06096 z^-3
%
%     A_1(z) = -0.108 z^-1 + 0.05012 z^-2
% Sample time: 1 seconds
% ARMAX
% sys =Discrete-time ARMA model:
%   Model for output "y1": A(z)y_1(t) = - A_i(z)y_i(t) + C(z)e_1(t)
%     A(z) = 1 - 1.487 z^-1 + 0.6851 z^-2
%
%     A_2(z) = -0.4975 z^-1
%
%     C(z) = 1 + 0.07709 z^-1 - 0.05804 z^-2
%
%   Model for output "y2": A(z)y_2(t) = - A_i(z)y_i(t) + C(z)e_2(t)
%     A(z) = 1 - 2.417 z^-1 + 2.053 z^-2 - 0.6112 z^-3
%
%     A_1(z) = -0.041 z^-1 + 0.03667 z^-2
%
%     C(z) = 1 - 1.027 z^-1 - 0.01797 z^-2 + 0.1476 z^-3
%
% Sample time: 1 seconds

if (useARX)
    nmax = max(na(:));
else
    nmax = max(max([na, nc]));
end
ypred = predict(sys,y(1:40),nToPredict);
tp = get(ypred,'SamplingInstants');
ymeas = y.y;
est_y = nan*ypred.y;
A = sys.A; C = sys.C;
yp = zeros(length(t),ny);
y0 = ymeas(1:nmax,:);
est_y(1:nmax,:) = y0;
err = zeros(size(yp));
for ct=(nmax+1):(length(tp)-nToPredict+1)
    %%Extract previous y
    prev_y = ymeas(ct-(1:nmax),:);
    err_ct = err(ct-(1:nmax),:);
    for iPred = 1:nToPredict
        est_y_prev = nan(1,ny);
        for i_out = 1:ny
            azy = 0;
            for j_in = 1:ny
                Ap = A{i_out,j_in}(2:end);
                azy = azy -  Ap * prev_y(1:na(i_out,j_in),j_in);
            end
            if (~useARX)
                azy = azy + C{i_out}(2:end)*err_ct(1:nc(i_out),i_out);
            end
            est_y_prev(1,i_out) = 1/A{i_out,i_out}(1)*azy ;
        end
        % Save y for next step
        prev_y = [est_y_prev; prev_y(1:end-1,:)];
        err_ct = [zeros(1, ny); err_ct(1:end-1,:)];
        if iPred==1, yp1 = est_y_prev; end
    end
    est_y(ct,:) = est_y_prev;
    err(ct,:) = ymeas(ct,:) - yp1;
end
figure;
hp_raw = plot(t(1:nToEval),y.y(1:nToEval,1),'b-d', ...
    t(1:nToEval),y.y(1:nToEval,2),'k-d', ...
    t(1:nToEval),yc.y(:,1),'b:x', ...
    t(1:nToEval),yc.y(:,2),'k:x');
%
hold all;
hp_predict=plot(tp,ypred.y(:,1),'bo', ...
    tp,ypred.y(:,2),'ko');
set(hp_predict,'MarkerSize',6);

% Add the estimate to the plot
hp_est=plot(t(((nmax+2):size(est_y,1))+nToPredict-1),est_y((nmax+2):end,1),'bs', ...
    t(((nmax+2):size(est_y,1))+nToPredict-1),est_y((nmax+2):end,2),'ks');
set(hp_est,'MarkerSize',10);
legend([hp_raw(1); hp_raw(3); hp_predict(1); hp_est(1)], ...
    'raw','compare','predict','estimated','location','best');
if (useARX)
    titleStr = 'ARX';
else
    titleStr = 'ARMAX';
end
titleStr = sprintf('Model %s, n to predict %d', titleStr, nToPredict);

title(titleStr);

Answer 2

[L L]

If there are two input and there is only one input, what are the shape of na, nb, nc, nk. Thank you.