如何科学地做一名巫师——使用ARMA做时间序列预测全流程

你有没有遇到过这样的问题：我有一段数据，它是随着时间等间隔采样的，现在想用某种方法预测出后续一段时间的趋势。这就是所谓的时间序列的预测问题。

时间序列预测的应用主要是在经济领域，比如预测股市行情，预测GDP走势；也可以用来预测机械性能退化趋势或者其他诸如某个量随时间变化这样的场景当中。

AR/MA/ARMA模型是分析时间序列的重要方法，在本篇文章中将重点介绍如何使用ARMA(p,q)模型对时间序列信号进行一次完整的预测。而对于ARMA(p,q)模型的相关理论知识则不进行重点介绍。没有涉及到的理论部分，建议阅读相关专业书籍。

本文的目的是，对于不了解时间序列与ARMA(p,q)的同学，可以套用文中的方法得到想要的预测结果；对于了解理论知识的同学，通过代码加深理解。

所以，这不是一篇ARMA模型入门的文章，而是一篇应用ARMA进行预测的傻瓜式示范教程，就算你对里边的理论完全不懂，相信也照样可以依样画葫芦做出预测（虽然不推荐这样！）。

1.导入实验数据

close all
clear all
load Data_EquityIdx   %纳斯达克综合指数
len = 120;
Y = DataTable.NASDAQ(1:len);
plot(Y)

2.平稳性检验

使用ARMA模型要求时间序列必须是平稳的，所以第一步是对原始数据进行平稳性检验。检验方法有很多种，包括ADF、KPSS、P-P等。这里用ADF检验和KPSS检验。

y_h_adf = adftest(Y)
y_h_kpss = kpsstest(Y)

得到y_h_adf =0，y_h_kpss =1，没有通过检验。

此时对Y取差分，再进行检验：

% 一阶差分，结果平稳。如果依旧不平稳的话，再次求差分，直至通过检验
Yd1 = diff(Y);
yd1_h_adf = adftest(Yd1)
yd1_h_kpss = kpsstest(Yd1)

得到yd1_h_adf =1，yd1_h_kpss =0，通过检验了。

如果此时并没有通过检验的话，可以取二阶差分，再不行就三阶，直到信号变得平稳，此时差分的阶数就是ARIMA(p,i,q)中的i的值。上述例子中一阶差分就平稳了，所以i=1。

3.确定ARMA模型阶数

模型阶数的确定方法主要包括两种：（1）ACF和PACF法。（2）基于准则的确定法。（都是我自己起的名）

（1）ACF和PACF法

所谓ACF即自相关函数，PACF为偏自相关函数。

% ACF和PACF法，确定阶数
figure
autocorr(aimY)
figure
parcorr(aimY)

做出的结果如下两图：

自相关函数ACF

偏自相关函数PACF

怎样根据这两站图确定p和q呢，有这样一张图：

通过ACF和PACF一般特征

对于我们关注的ARMA(p,q)，通俗地说，PACF最后一个在蓝线外（即阈值外）的Lag值就是p值；ACF最后一个在蓝线外（即阈值外）的Lag值就是q值。

如果按照这个标准，p似乎要定到13，q也要定到12，这个值有些偏大了。

（2）基于准则的确定法

如果方法（1）不理想，可以用方法（2），这个方法不用纠结，可以直接给出推荐结果。代码如下：

% 通过AIC，BIC等准则暴力选定阶数
max_ar = 3;
max_ma = 3;
[AR_Order,MA_Order] = ARMA_Order_Select(Y,max_ar,max_ma,1);     

ARMA_Order_Select是我自己写的一个函数，只需要输入信号、最大p值、最大q值和差分阶数就可以得到推荐p和q。

运行结果为：

AR_Order =    3
MA_Order =    3

即p=1,q=3。

4.残差检验

为了确保确定的阶数合适，还需要进行残差检验。残差即原始信号减掉模型拟合出的信号后的残余信号。如果残差是随机正态分布的、不自相关的，这说明残差是一段白噪声信号，也就说明有用的信号已经都被提取到ARMA模型中了。

Mdl = arima(AR_Order, 1, MA_Order);  %第二个变量值为1，即一阶差分
EstMdl = estimate(Mdl,data);
[res,~,logL] = infer(EstMdl,data);   %res即残差

stdr = res/sqrt(EstMdl.Variance);
figure('Name','残差检验')
subplot(2,3,1)
plot(stdr)
title('Standardized Residuals')
subplot(2,3,2)
histogram(stdr,10)
title('Standardized Residuals')
subplot(2,3,3)
autocorr(stdr)
subplot(2,3,4)
parcorr(stdr)
subplot(2,3,5)
qqplot(stdr)

残差检验

上图为残差检验的结果图。Standardized Residuals是查看残差是否接近正态分布，理想的残差要接近正态分布；ACF和PACF检验残差的自相关和偏自相关，理想的结果应该在图中不存在超出蓝线的点；最后一张QQ图是检验残差是否接近正太分布的，理想的结果中蓝点应该靠近红线。

除了上述图像检验方法，还可以通过Durbin-Watson对相关性进行检验：

% Durbin-Watson 统计是计量经济学分析中最常用的自相关度量
diffRes0 = diff(res);  
SSE0 = res'*res;
DW0 = (diffRes0'*diffRes0)/SSE0 % Durbin-Watson statistic，该值接近2，则可以认为序列不存在一阶相关性。

运算结果为1.9997，这个值越接近2越说明残差不存在一阶相关性。

上述检验可以证明，残差接近正态分布，且相互独立，可以认为ARMA建模符合要求。

5.预测

万事俱备，终于来到最后一步。

%% 5.预测
step = 300; %预测步数为300
[forData,YMSE] = forecast(EstMdl,step,'Y0',data);   %matlab2018及以下版本写为Predict_Y(i+1) = forecast(EstMdl,1,'Y0',Y(1:i)); 
lower = forData - 1.96*sqrt(YMSE); %95置信区间下限
upper = forData + 1.96*sqrt(YMSE); %95置信区间上限

figure()
plot(data,'Color',[.7,.7,.7]);
hold on
h1 = plot(length(data):length(data)+step,[data(end);lower],'r:','LineWidth',2);
plot(length(data):length(data)+step,[data(end);upper],'r:','LineWidth',2)
h2 = plot(length(data):length(data)+step,[data(end);forData],'k','LineWidth',2);
legend([h1 h2],'95% 置信区间','预测值',...
	     'Location','NorthWest')
title('Forecast')
hold off

多步预测结果

上图中灰线为用来训练的1200个数据，黑线为未来值的预测，红线为95%置信区间上下限。也就是说未来真实值有95%的概率落在这个范围内。

可以看到使用ARIMA方法进行长期预测的结果是趋势性的。

6.封装

上述整个过程可以封装成一个函数，如下：

%% 进行使用ARIMA进行预测的函数
function [forData,lower,upper] = Fun_ARIMA_Forecast(data,step,max_ar,max_ma,figflag)
%  使用ARIMA进行预测的函数（使用n阶差分、不使用对数），可以直接调用，差分阶数自动确定。
% 输入：
% data为待预测数据，一维数据，最小11个数据。但是数据长度处于11~15时依旧可能出现报错的情况。
% step为拟预测步数
% max_ar 为最大p值
% max_ma 为最大q值
% figflag 为画图标志位，'on'为画图，'off'为不画
% 输出：
% forData为预测结果，其长度等于step
% lower为预测结果的95%置信下限值
% upper为预测结果的95%置信上限值

直接调用就可以出结果。

使用上边封装好的函数试一下循环迭代的单步预测效果：其思路就是比如从30开始，用前30个数据预测第31个；再用前31个预测第32个...用前100个预测第101个。然后再把所有单步预测结果画出来和原数据做对比。

%% 使用Fun_ARIMA_Forecast函数实现预测
% 单步预测
close all
clear all
addpath ../funs %将funs文件夹添加进路径
load Data_EquityIdx   %纳斯达克综合指数
len = 100;
data = DataTable.NASDAQ(1:len); %如果要替换数据，将此处data替换即可。
forData1 = zeros(1,len); %全部初始化为0
for i = 30:len
    forData1(i+1) = Fun_ARIMA_Forecast(data(1:i),1,2,2,'off');
end
figure()
plot(31:len,data(31:len))
hold on
plot(31:len,forData1(31:len))
legend('原始数据','单步预测数据')

画出的结果如下（对比从30~100的数据段）：

如果要获取本文的完整代码，可以在下述链接中免费获取：

ARIMA代码说明文档（公开版） | 工具箱文档

ARIMA时间序列预测图形界面版软件（试用版） | 工具箱文档

为了不再开坑，一篇讲完了

参考：

MATLAB 文档中心 - MathWorks 中国

【时间序列分析】ARMA预测GDP的python实现