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

正文

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


完成本教程后,您将了解:

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

自回归模型

自回归模型,即线性回归,基于历史输入值的线性组合对输出值进行建模。
模型如下:
y t = c + ϕ 1 ∗ y t − 1 + ϕ 2 ∗ y t − 2 + . . . + ϕ p ∗ y t − p + ϵ t y_t = c + \phi_1*y_{t-1}+\phi_2*y_{t-2}+...+\phi_p*y_{t-p}+\epsilon_t yt=c+ϕ1yt1+ϕ2yt2+...+ϕpytp+ϵt
式中: y t y_t yt是时间序列 Y \textbf{Y} Yt时刻的观测值, ϕ t \phi_t ϕt是通过对训练数据优化模型(例如最小二乘)得到的系数, ϵ t \epsilon_t ϵt是t时刻的残差, c c c为模型的常数项。

移动平均模型

移动平均模型是基于移动平均过程,是一种常见的模拟时间序列过程。移动平均模型的输出变量是随机项的当前值和各种过去值线性组合。模型如下:
y t = c + ϵ t + θ 1 ∗ ϵ t − 1 + θ 2 ∗ ϵ t − 2 + . . . + θ q ∗ ϵ t − q y_t = c +\epsilon_t+ \theta_1*\epsilon_{t-1}+\theta_2*\epsilon_{t-2}+...+\theta_q*\epsilon_{t-q} yt=c+ϵt+θ1ϵt1+θ2ϵt2+...+θqϵtq
其中 ϵ t \epsilon_t ϵt理解为均值为零的不相关正态分布的随机变量(实质是一个 innovation process)。

自回归移动平均模型

自回归移动平均模型就是上述两种模型的组合,模型如下:
y t = c + ϕ 1 ∗ y t − 1 + ϕ 2 ∗ y t − 2 + . . . + ϕ p ∗ y t − p + ϵ t + θ 1 ∗ ϵ t − 1 + θ 2 ∗ ϵ t − 2 + . . . + θ q ∗ ϵ t − q y_t = c + \phi_1*y_{t-1}+\phi_2*y_{t-2}+...+\phi_p*y_{t-p}+\epsilon_t+ \theta_1*\epsilon_{t-1}+\theta_2*\epsilon_{t-2}+...+\theta_q*\epsilon_{t-q} yt=c+ϕ1yt1+ϕ2yt2+...+ϕpytp+ϵt+θ1ϵt1+θ2ϵt2+...+θqϵtq

部署模型

本文仅开发了简单的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 ∗ y t − 1 + 0.077105 ∗ y t − 2 + ϵ t y_t = 2.3181 + 0.71548*y_{t-1}+0.077105*y_{t-2}+\epsilon_t yt=2.3181+0.71548yt1+0.077105yt2+ϵt
即:
y t − ϵ t = y ^ t = 2.3181 + 0.71548 ∗ y t − 1 + 0.077105 ∗ y t − 2 y_t-\epsilon_t=\hat{y}_t = 2.3181 + 0.71548*y_{t-1}+0.077105*y_{t-2} ytϵt=y^t=2.3181+0.71548yt1+0.077105yt2
实施预测可以直接使用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 ∗ ϵ t − 1 + 0.3554 ∗ ϵ t − 2 + ϵ t y_t = 11.184 + 0.75973*\epsilon_{t-1}+0.3554*\epsilon_{t-2}+\epsilon_t yt=11.184+0.75973ϵt1+0.3554ϵt2+ϵt
即:
y t − ϵ t = y ^ t = 11.184 + 0.75973 ∗ ϵ t − 1 + 0.3554 ∗ ϵ t − 2 y_t-\epsilon_t=\hat{y}_t = 11.184 + 0.75973*\epsilon_{t-1}+0.3554*\epsilon_{t-2} ytϵt=y^t=11.184+0.75973ϵt1+0.3554ϵt2
实施预测可以直接使用forecast()函数:

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

预报结果为:
在这里插入图片描述
使用forecat()函数自动进行预测,结果如下:
在这里插入图片描述

使用手动获取参数按照上述公式进行预测:
这里我注意到, ϵ t − 1 \epsilon_{t-1} ϵt1 ϵ t − 2 \epsilon_{t-2} ϵt2都是没法观测的,未知的。拿 ϵ t − 1 \epsilon_{t-1} ϵt1为例,我们想要求出 ϵ t − 1 \epsilon_{t-1} ϵt1就需要知道 y ^ t − 1 \hat{y}_{t-1} y^t1,因为 ϵ t − 1 = y t − 1 − y ^ t − 1 \epsilon_{t-1} =y_{t-1}-\hat{y}_{t-1} ϵt1=yt1y^t1,但是 y ^ t − 1 = 11.184 + 0.75973 ∗ ϵ t − 2 + 0.3554 ∗ ϵ t − 3 \hat{y}_{t-1}=11.184 + 0.75973*\epsilon_{t-2}+0.3554*\epsilon_{t-3} y^t1=11.184+0.75973ϵt2+0.3554ϵt3,从这里可以看出这是个递归的过程,需要设置初值,迭代进行确定,我这里直接将时段初的 ϵ \epsilon ϵ设置为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);

模型估计结果:
在这里插入图片描述
公式:
y ^ t = 0.073296 + 1.233 ∗ y t − 1 − 0.23973 ∗ y t − 2 − 0.64267 ∗ ϵ t − 1 − 0.23219 ∗ ϵ t − 2 \hat{y}_t = 0.073296 + 1.233*y_{t-1}-0.23973*y_{t-2}-0.64267*\epsilon_{t-1}-0.23219*\epsilon_{t-2} y^t=0.073296+1.233yt10.23973yt20.64267ϵt10.23219ϵt2
手动预报和自动预报结果对比:
在这里插入图片描述
有较小的差距。

参考链接

  • 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
评论 17
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值