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通过平稳性检验提前确定）。
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，不然无法预测未来。
idxpresample = 1;%预采样数据y0es个数，1
idxestimate = 2:56;%训练数据yest 个数,55
idxforecast = 57:T;%预测数据个数 ，6
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
[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期”即可获取下载方式。