How to simulate a multivariable autoregressive model forecast

Question

John 2016-10-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/306142-how-to-simulate-a-multivariable-autoregressive-model-forecast

回答： L L 2018-7-10

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

John 2016-10-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/306142-how-to-simulate-a-multivariable-autoregressive-model-forecast#answer_238663

编辑：John 2016-10-12

在 MATLAB Online 中打开

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