How to simulate a multivariable autoregressive model forecast

16 views (last 30 days)
John on 6 Oct 2016
Answered: L L on 10 Jul 2018
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');

John on 12 Oct 2016
Edited: John on 12 Oct 2016
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);