如何MATLAB实现用ARIMA模型输出参数实施预测

正文

自回归（AR）模型、移动平均（MA）模型、自回归移动平均（ARMA）和自回归差分移动平均（ARIMA）模型是时间序列模型，它们主要是使用历史时间步的观测值作为回归方程的输入，以预测下一时间步的值。这是一个非常简单的想法，可以导致对一系列时间序列问题的准确预测。在本教程中，您将了解如何使用MATLAB实现时间序列预测模型。

---

完成本教程后，您将了解：

• 如何部署一个时间序列模型并进行预测。
• 如何获取已经估计的时间序列模型的参数实施直接进行预测。
• 结合公式更加深入了解自回归移动平均等时间序列模型。

自回归模型

自回归模型，即线性回归，基于历史输入值的线性组合对输出值进行建模。
模型如下：
$$y_t=c+{\varphi }_1\ast y_{t-1}+{\varphi }_2\ast y_{t-2}+...+{\varphi }_p\ast y_{t-p}+{ϵ}_t$$
式中：$y_t$是时间序列$\text{Y}$t时刻的观测值，${\varphi }_t$是通过对训练数据优化模型（例如最小二乘）得到的系数，${ϵ}_t$是t时刻的残差，$c$为模型的常数项。

移动平均模型

移动平均模型是基于移动平均过程，是一种常见的模拟时间序列过程。移动平均模型的输出变量是随机项的当前值和各种过去值线性组合。模型如下：
$$y_t=c+{ϵ}_t+{\theta }_1\ast {ϵ}_{t-1}+{\theta }_2\ast {ϵ}_{t-2}+...+{\theta }_q\ast {ϵ}_{t-q}$$
其中${ϵ}_t$理解为均值为零的不相关正态分布的随机变量（实质是一个 innovation process）。

自回归移动平均模型

自回归移动平均模型就是上述两种模型的组合，模型如下：
$$y_t=c+{\varphi }_1\ast y_{t-1}+{\varphi }_2\ast y_{t-2}+...+{\varphi }_p\ast y_{t-p}+{ϵ}_t+{\theta }_1\ast {ϵ}_{t-1}+{\theta }_2\ast {ϵ}_{t-2}+...+{\theta }_q\ast {ϵ}_{t-q}$$

部署模型

本文仅开发了简单的AR、MA、ARMA和ARIMA模型，参数没有进行优化，用于演示，参数只要稍加调整，就可获得更好的预测效果。

数据

在示例中使用的最低日温度数据。

• 下载最低温度数据集
  下载数据集到当前工作目录中，命名为“daily_minimum_temperatures.csv”将。
  下面的代码将数据集作为一个数组加载。

clc;
clear;
% load data
file = fopen("daily-minimum-temperatures.csv");
fmt = '"%u-%u-%u" %f'
if file>0    
     series = textscan(file,fmt,'Delimiter',',','HeaderLines',1);
     % close the file
     fclose(file);         
end
y = series{:,4}; % 仅取数值使用
plot(y);

然后创建数据集的折线图:

AR模型测试

下面演示，建立一个AR(2)模型对未来7天的值进行预测，同时写出方程，并通过取系数和滞后值的点积来计算手动输出值。
给出参数数量p = 2（即模型形式），建立模型，然后估计参数

% model
AR_Order = 2;
MA_Order = 0;
AR2 = arima(AR_Order, 0, MA_Order);
EstMdl = estimate(AR2,y);

估计模型的结果直接会在窗口输出：

所以估计得到的方程为：
$$y_t=2.3181+0.71548\ast y_{t-1}+0.077105\ast y_{t-2}+{ϵ}_t$$
即：
$$y_t-{ϵ}_t={\hat{y}}_t=2.3181+0.71548\ast y_{t-1}+0.077105\ast y_{t-2}$$
实施预测可以直接使用forecast()函数：

step = 7;
auto_fore = forecast(EstMdl,step,'Y0',y);
auto_fore'

输出结果：

使用手动获取参数按照上述公式进行预测：

mannual_fore = size(1:step); %预分配内存
history = y;
for i=1:step
    lags = history(end-AR_Order+1:end); % 获取滞后项目
    lags = rot90(lags,2); % 翻转一下顺序和系数要对应
    yhat = cell2mat(EstMdl.AR)*lags + EstMdl.Constant; % 可以理解为上述公式的矩阵形式
    history = [history; yhat]; %将预测值加入到历史数据中，因为下一时段的滚动预测需要用到上一个时段的预测值
end
mannual_fore = history(end-step+1:end);
mannual_fore’

输出结果：

可以将上述结果进行相减，看是否是完全一致的：

完全一致。

MA模型测试

与AR模型类似，先建模估计参数，下面以MA(2)模型为例

% model
AR_Order = 0;
MA_Order = 2;
MA2 = arima(AR_Order, 0, MA_Order);
EstMdl = estimate(MA2,y);

估计结果：

所以估计得到的方程为：
$$y_t=11.184+0.75973\ast {ϵ}_{t-1}+0.3554\ast {ϵ}_{t-2}+{ϵ}_t$$
即：
$$y_t-{ϵ}_t={\hat{y}}_t=11.184+0.75973\ast {ϵ}_{t-1}+0.3554\ast {ϵ}_{t-2}$$
实施预测可以直接使用forecast()函数：

step = 7;
auto_fore = forecast(EstMdl,step,'Y0',y);

预报结果为：

使用forecat()函数自动进行预测，结果如下：

使用手动获取参数按照上述公式进行预测：
这里我注意到，${ϵ}_{t-1}$和${ϵ}_{t-2}$都是没法观测的，未知的。拿${ϵ}_{t-1}$为例，我们想要求出${ϵ}_{t-1}$就需要知道${\hat{y}}_{t-1}$，因为${ϵ}_{t-1}=y_{t-1}-{\hat{y}}_{t-1}$，但是${\hat{y}}_{t-1}=11.184+0.75973\ast {ϵ}_{t-2}+0.3554\ast {ϵ}_{t-3}$，从这里可以看出这是个递归的过程，需要设置初值，迭代进行确定，我这里直接将时段初的$ϵ$设置为0，然后进行模拟求出历史数据的残差序列，再实施预测，MATLAB源码可能不是这么干的，我这块也不太懂，暂时这么处理：

mannual_fore = size(1:step);
history = y;
[len,~] = size(history);
residuals = size(1:len);
residuals(1:MA_Order) = 0; % 按照残差的阶数将初始值设置为0

使用这个初值进行模拟求解，找出历史数据的其他初值，基本思路就是先计算模拟值，然后使用观测值和模拟计算出残差：

for i=MA_Order+1:len
    resids = residuals(i-2:i-1);
    resids = rot90(resids,2);
    resid = history(i) - (cell2mat(EstMdl.MA)*resids' + EstMdl.Constant);
    residuals(i) = resid;
end

计算出残差序列之后，根据上述公式和残差序列计算样本外的预测值：

for i=1:step
    resids = residuals(end-MA_Order+1:end);
    resids = rot90(resids,2);
    yhat = cell2mat(EstMdl.MA)*resids' + EstMdl.Constant;
    residuals = [residuals, 0];
    history = [history; yhat];
    mannual_fore = history(end-step+1:end);
    mannual_fore'
end

手动预测的输出结果为：

结合上面手动和自动的预测结果可以看出，两者相同，然后将两者相减

可以看出来结果不是完全一样，第二个结果有点差距，总体来看差距很小。

ARMA模型测试

ARMA模型的测试类比上述AR和MA组合即可，这里不再赘述，仅给出代码、公式和运行结果，使用的示例是ARMA(2,0,2)。
代码：

clc;
clear;
% load data
file = fopen("daily-minimum-temperatures.csv");
fmt = '"%u-%u-%u" %f'
if file>0    
     series = textscan(file,fmt,'Delimiter',',','HeaderLines',1);
     % close the file
     fclose(file);         
end
y = series{:,4};
%plot(y);
% model
AR_Order = 2;
MA_Order = 2;
MA1 = arima(AR_Order, 0, MA_Order);
EstMdl = estimate(MA1,y);
step = 10;
auto_fore = forecast(EstMdl,step,'Y0',y);
mannual_fore = size(1:step);
history = y;
[len,~] = size(history);
residuals = size(1:len);
max_order = max(MA_Order,AR_Order);
residuals(1:max_order) = 0;

for i=max_order+1:len
    lags = history(i-AR_Order:i-1);
    lags = rot90(lags,2);
    resids = residuals(i-MA_Order:i-1);
    resids = rot90(resids,2);
    resid = history(i) - (cell2mat(EstMdl.AR)*lags + cell2mat(EstMdl.MA)*resids' + EstMdl.Constant);
    residuals(i) = resid;
end

for i=1:step
    lags = history(end-AR_Order+1:end);
    lags = rot90(lags,2);
    resids = residuals(end-MA_Order+1:end);
    resids = rot90(resids,2);
    yhat = cell2mat(EstMdl.AR)*lags + cell2mat(EstMdl.MA)*resids' + EstMdl.Constant;
    residuals = [residuals, 0];
    history = [history; yhat];
end
mannual_fore = history(end-step+1:end);

模型估计结果：

公式：
$${\hat{y}}_t=0.073296+1.233\ast y_{t-1}-0.23973\ast y_{t-2}-0.64267\ast {ϵ}_{t-1}-0.23219\ast {ϵ}_{t-2}$$
手动预报和自动预报结果对比：

有较小的差距。

参考链接

• https://en.wikipedia.org/wiki/Moving-average_model
• https://www.mathworks.com/help/econ/arma-model.html?searchHighlight=arma&s_tid=srchtitle
• https://machinelearningmastery.com/make-manual-predictions-arima-models-python/
• https://encyclopediaofmath.org/wiki/Stochastic_process,_renewable