【MATLAB第60期】【更新中】基于MATLAB的ARMAX具有外生回归因子的移动平均自回归模型

【MATLAB第60期】【更新中】基于MATLAB的ARMAX具有外生回归因子的移动平均自回归模型

版本更新：

2023/7/29版本：
1.增加自定义参数，方便直接套数据运行。

pre_num=3;%预采样数据个数
learn_pr=0.85;  %训练数据比例（不包括预采样数据）
mm=pre_num;%输入响应数据个数    
nn=pre_num;%输出响应数据个数   

2.增加ARIMAX模型参数自动选择功能
（1）可手动选择ADF或者KPSS平稳性检验方式
存在不足：目前无法解决D>0的问题，所以优先选择D=0的检验方法，本文选用KPSS检验
[p, d, q ] = fit_model( Y,learn_num, test_num );
（2）检验p、d、q是否满足后续正常运行条件

if pre_num>=p+d
    disp('--------------------------------------------');
     fprintf('ARIMAX(%d, %d, %d)满足运行要求', p, d, q);
else
    disp('--------------------------------------------');
     fprintf('ARIMAX(%d, %d, %d)不满足运行要求,', p, d, q);
     disp('请增大pre_num预采样数值');
end

3.增加参数评估结果（命令行窗口展示）

最优模型参数为：ARIMAX (2, 0, 0)
ARIMAX(2,0,0) Model (Gaussian Distribution)
 
    Effective Sample Size: 53
    Number of Estimated Parameters: 7
    LogLikelihood: -164.451
    AIC: 342.901
    BIC: 356.693
 
                Value     StandardError    TStatistic    PValue
                ______    _____________    __________    ______

    Constant    -25.21        10.11          -2.49        0.01 
    AR{1}         0.30         0.08           3.63        0.00 
    AR{2}        -0.06         0.06          -1.01        0.31 
    Beta(1)       3.24         0.32          10.27        0.00 
    Beta(2)       0.00         0.00           3.21        0.00 
    Beta(3)       1.92         0.58           3.31        0.00 
    Variance     29.01         7.63           3.81        0.00 

4.增加置信区间绘图及评价（默认95%，可自行调整）

picp2 = PICP (ci, Y(end-test_num+1:end));
pimw2 = PIMWP(ci, Y(end-test_num+1:end));
disp(['测试集的区间覆盖率为:', num2str(picp2), '。区间平均宽度百分比为:', num2str(pimw2)])
测试集的区间覆盖率为:1。区间平均宽度百分比为:3.212

版本不足：

1.未解决D>0的问题，即只能处理通过检验方法判定的平稳的时间序列数据。
2.p、q筛选区间最大值目前只能为2，且缺少d结果的选择（d通过平稳性检验提前确定）。

一、简要介绍

1. ARMAX模型相比ARMA考虑了影响因素 ，即可以实现基于时间序列数据的回归预测。
2. 目前，ARMAX预测未来功能存在困难，本篇文章不予介绍。大致思路需要通过时间滞后构造数据，使前时间段的X预测后时间段的Y，即多步预测。
3. 此示例展示如何将时间序列中的时间划分为预采样期T0、训练期Ty和预测期Tf，并显示了如何提供适当数量的观测值来初始化用于估计和预测的动态模型。
4. 通过定义ARMA模型中的参数，可实现ARIMAX和SARIMAX模型。本文介绍最基础的ARMAX模型。
   

二、导入数据

1. 本篇文章案例数据采用3输入1输出，62个样本（1962-2023年）。
2. 本文用table格式打开，方便对时间进行处理。

clear all
data=readtable('数据集.xlsx');
DataTable = table2timetable(data);%将DataTable转换为时间表。
varnames = ["Y" "X1" "X2" "X3" ];
Tbl = rmmissing(DataTable(:,varnames));%通过应用列表删除从数据中删除所有前导NaN。
T = size(Tbl,1) %总样本量
Y = Tbl.Y; %因变量
X = Tbl{:,varnames(2:end)};%变量

三、建立模型

为了训练和预测模型，估计必须有足够的预采样数据来初始化自回归项，同样，要从训练模型中预测，预测必须有足够的预采样样本。
此外，预测期中的回归分量需要预测历史数据或未来的预测数据Y，那么需要有与之对应的X，不然无法预测未来。

1. 本文考虑一个ARMAX（1,2）模型，该模型以X1、X2、X3为外生变量，预测Y。
2. 将样本的时间线划分为预采样、训练和预测时段。将模型拟合到训练样本，并使用预采样数据来初始化自回归项。然后，根据训练模型对Y进行预测。
3. 指定预采样数据，以初始化自回归项。一般预采样数据个数为Mdl.P，因为p在之前就设置好了，所以手动设置为1。 -
4. 指定训练数据，选择2-56作为训练数据。输入数据XEst则为X2-X56，输出为Y2-Y56
5. 指定预测数据，57-62共6个数据进行测试。输入数据则为X57-X62，输出为Y57-Y62

idxpresample = 1;%预采样数据y0es个数，1
idxestimate = 2:56;%训练数据yest 个数,55
idxforecast = 57:T;%预测数据个数 ，6

1. 建立ARMAX（1,2）模型

Mdl = arima(1,0,2); % P D Q
%ARIMAX(1,0,2) Model (Gaussian Distribution)

Effective Sample Size: 55
Number of Estimated Parameters: 8
LogLikelihood: -162.152
AIC: 340.303
BIC: 356.362

            Value     StandardError    TStatistic    PValue
            ______    _____________    __________    ______

Constant    -28.86        12.92          -2.23        0.03 
AR{1}         0.20         0.05           4.02        0.00 
MA{1}         0.65         0.16           4.19        0.00 
MA{2}         0.05         0.18           0.31        0.76 
Beta(1)       3.42         0.28          12.37        0.00 
Beta(2)       0.00         0.00           3.19        0.00 
Beta(3)       1.96         0.76           2.57        0.01 
Variance     21.30         5.16           4.13        0.00 

1. 在训练样本结束时指定必要的观测值作为样本前数据进行预测，需指定训练期的数据，且数据个数至少为1，本文取两个，即训练输入的最后2个值X55-X56和训练输出的最后2个值Y55-Y56。
2. 预测数据，假设预测的数量为M，则M必须小于等于XF的个数，不然无法运行。本文指定M=6，预测期的输入变量XF为X57-X62。

 [yf,ymse] = forecast(Mdl,M);

1. 置信区间预测

ci = yf + 1.96*[-sqrt(ymse) sqrt(ymse)];

1. 绘图。因年份较多，故只展示后面一半的数据。

yrs = year(Tbl.Time(round(T/2):end));%绘制后半部分的响应数据和预测。

figure;
plot(yrs,Tbl.Y(round(T/2):end),"b","LineWidth",2);
hold on
plot(yrs(end-size(idxforecast,2)+1:end),yf,"r--","LineWidth",2);
h = gca;
px = yrs([end - size(idxforecast,2)+1 end end end - size(idxforecast,2)+1]);
py = h.YLim([1 1 2 2]);
hp = patch(px,py,[0.9 0.9 0.9]);
uistack(hp,"bottom");
axis tight
title("ARMAX模型");
legend(["预测区段" "实际值" "预测值"])

四、效果展示

五、代码获取

后台私信回复“60期”即可获取下载方式。