【MATLAB第60期】源码分享 | 基于MATLAB的ARMAX具有外生回归因子的移动平均自回归模型
一、简要介绍
- ARMAX模型相比ARMA考虑了影响因素 ，即可以实现基于时间序列数据的回归预测。
- 目前，ARMAX预测未来功能存在困难，本篇文章不予介绍。大致思路需要通过时间滞后构造数据，使前时间段的X预测后时间段的Y，即多步预测。
- 此示例展示如何将时间序列中的时间划分为预采样期T0、训练期Ty和预测期Tf，并显示了如何提供适当数量的观测值来初始化用于估计和预测的动态模型。
- 通过定义ARMA模型中的参数，可实现ARIMAX和SARIMAX模型。本文介绍最基础的ARMAX模型。
二、导入数据
- 本篇文章案例数据采用3输入1输出，62个样本（1962-2023年）。
- 本文用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，不然无法预测未来。
- 本文考虑一个ARMAX（1,2）模型，该模型以X1、X2、X3为外生变量，预测Y。
- 将样本的时间线划分为预采样、训练和预测时段。将模型拟合到训练样本，并使用预采样数据来初始化自回归项。然后，根据训练模型对Y进行预测。
- 指定预采样数据，以初始化自回归项。一般预采样数据个数为Mdl.P，因为p在之前就设置好了，所以手动设置为1。 -
- 指定训练数据，选择2-56作为训练数据。输入数据XEst则为X2-X56，输出为Y2-Y56
- 指定预测数据，57-62共6个数据进行测试。输入数据则为X57-X62，输出为Y57-Y62
idxpresample = 1;%预采样数据y0es个数，1
idxestimate = 2:56;%训练数据yest 个数,55
idxforecast = 57:T;%预测数据个数 ，6
- 建立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，本文取两个，即训练输入的最后2个值X55-X56和训练输出的最后2个值Y55-Y56。
- 预测数据，假设预测的数量为M，则M必须小于等于XF的个数，不然无法运行。本文指定M=6，预测期的输入变量XF为X57-X62。
[yf,ymse] = forecast(Mdl,M);
- 置信区间预测
ci = yf + 1.96*[-sqrt(ymse) sqrt(ymse)];
- 绘图。因年份较多，故只展示后面一半的数据。
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期”即可获取下载方式。
