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);
