ARIMA模型——本质上是error和t-?时刻数据差分的线性模型！！！如果数据序列是非平稳的，并存在一定的增长或下降趋势，则需要对数据进行差分处理!ARIMA（p，d，q）称为差分自回归移动平均模型，AR是自回归， p为自回归项； MA为移动平均，q为移动平均项数，d为时间序列成为平稳时所做的差分次数
https://www.cnblogs.com/bradleon/p/6827109.html 文章里写得非常好，需详细看。尤其是arima的举例！
可以看到：ARIMA本质上是error和t-?时刻数据差分的线性模型！！！
ARIMA模型全称为自回归积分滑动平均模型(Autoregressive Integrated Moving Average Model,简记ARIMA)，是由博克思(Box)和詹金斯(Jenkins)于70年代初提出一著名时间序列（Time-series Approach）预测方法 [1] ，所以又称为Box-Jenkins模型、博克思-詹金斯法。其中ARIMA（p，d，q）称为差分自回归移动平均模型，AR是自回归， p为自回归项； MA为移动平均，q为移动平均项数，d为时间序列成为平稳时所做的差分次数。所谓ARIMA模型，是指将非平稳时间序列转化为平稳时间序列，然后将因变量仅对它的滞后值以及随机误差项的现值和滞后值进行回归所建立的模型。ARIMA模型根据原序列是否平稳以及回归中所含部分的不同，包括移动平均过程（MA）、自回归过程（AR）、自回归移动平均过程（ARMA）以及ARIMA过程。
优点： 模型十分简单，只需要内生变量而不需要借助其他外生变量。
缺点：
1.要求时序数据是稳定的（stationary），或者是通过差分化(differencing)后是稳定的。
2.本质上只能捕捉线性关系，而不能捕捉非线性关系。
注意，采用ARIMA模型预测时序数据，必须是稳定的，如果不稳定的数据，是无法捕捉到规律的。比如股票数据用ARIMA无法预测的原因就是股票数据是非稳定的，常常受政策和新闻的影响而波动。
严谨的定义： 一个时间序列的随机变量是稳定的，当且仅当它的所有统计特征都是独立于时间的（是关于时间的常量）。
判断的方法：
- 稳定的数据是没有趋势(trend)，没有周期性(seasonality)的; 即它的均值，在时间轴上拥有常量的振幅，并且它的方差，在时间轴上是趋于同一个稳定的值的。
- 可以使用Dickey-Fuller Test进行假设检验。
ARIMA模型有三个参数:p,d,q。
- p--代表预测模型中采用的时序数据本身的滞后数(lags) ,也叫做AR/Auto-Regressive项
- d--代表时序数据需要进行几阶差分化，才是稳定的，也叫Integrated项。
- q--代表预测模型中采用的预测误差的滞后数(lags)，也叫做MA/Moving Average项
ARIMA建模基本步骤
- 获取被观测系统时间序列数据；
- 对数据绘图，观测是否为平稳时间序列；对于非平稳时间序列要先进行d阶差分运算，化为平稳时间序列；
- 经过第二步处理，已经得到平稳时间序列。要对平稳时间序列分别求得其自相关系数ACF 和偏自相关系数PACF，通过对自相关图和偏自相关图的分析，得到最佳的阶层 p 和阶数 q
- 由以上得到的d、q、p，得到ARIMA模型。然后开始对得到的模型进行模型检验。
具体例子会在另一篇文章中给出。
from:https://www.cnblogs.com/bradleon/p/6827109.html
预测程序
使用ARIMA模型对裙子长度预测
from：https://www.cnblogs.com/ECJTUACM-873284962/p/7379717.html
1、加载数据
skirts <- scan("http://robjhyndman.com/tsdldata/roberts/skirts.dat", skip=5)
str(skirts)
head(skirts)
boxplot(skirts)
length(skirts)
2、把数据转化为是时间序列
skirts_ts <- ts(skirts, start=c(1886), frequency=1)
1)查看时间序列对应的时间
skirts_ts
2)画出时间序列图
plot.ts(skirts_ts)
从图可知：女人裙子边缘的直径做成的时间序列数据，从 1866 年到 1911 年在平均值上是不平稳的
3、做差分得到平稳序列
1)做时间序列的一阶差分
skirts_diff <- diff(skirts_ts, differences = 1)
plot.ts(skirts_diff)
从一阶差分的图中可以看出，数据仍是不平稳的，继续差分
2)做时间序列的二阶差分
skirts_diff2 <- diff(skirts_ts, differences = 2)
plot.ts(skirts_diff2)
二次差分后的时间序列在均值和方差上看起来是平稳了
4、找到合适的ARIMA模型
寻找 ARIMA(p,d,q)中合适的 p 值和 q
1)自相关图ACF
acf(skirts_diff2, lag.max = 20)
acf(skirts_diff2, lag.max = 20, plot = F)
自相关图显示滞后1阶自相关值基本没有超过边界值，虽然5阶自相关值超出边界，那么很可能属于偶然出现的，而自相关值在其他上都没有超出显著边界， 而且我们可以期望 1 到 20 之间的会偶尔超出 95%的置信边界。 自相关图5阶后结尾
2)偏相关图PACF
pacf(skirts_diff2, lag.max = 20)
pacf(skirts_diff2, lag.max = 20, plot = F)
偏自相关值选1阶后结尾
故我们的ARMIA模型为armia（1,2,5
3)使用auto.arima()函数，自动获取最佳的ARIMA模型
library(forecast)
auto.arima(skirts_ts, ic=c("aicc", "aic", "bic"), trace = T)
Best model: ARIMA(1,2,0)
5、建立ARIMA模型：并对比arima(1, 2, 0)与arima(1, 2, 5)模型
1）arima(1, 2, 0)模型
(skirts_arima <- arima(skirts_ts, order = c(1, 2, 0)))
aic = 391.33
2）arima(1, 2, 5)模型
(skirts_arima <- arima(skirts_ts, order = c(1, 2, 5)))
aic = 381.6
AIC是赤池消息准则SC是施瓦茨准则，当两个数值最小时，则是最优滞后分布的长度。我们进行模型选择时，AIC值越小越好。所以arima(1, 2, 5)模型较好
6、预测：预测5年后裙子的边缘直径
(skirts_forecast <- forecast.Arima(skirts_arima, h=5, level = c(99.5)))
plot.forecast(skirts_forecast)
ARIMA 模型实践
from：https://www.jianshu.com/p/4130bac8ebec
模型具体的理论知识就不再做过多说明了，来个实际的例子吧。
ARIMA 模型对湖北省 GDP 的实证分析及预测
这里的例子是采用了一篇论文的数据，【ARIMA模型在湖北省GDP预测中的应用】，可以去中国知网搜索篇名进行下载。
| 年份 | GDP |
|---|---|
| 1978 | 151.00 |
| 1979 | 188.46 |
| 1980 | 199.38 |
| ... | ... |
| 2013 | 24668.49 |
数据的平稳性处理及检验
这里我们用 Python 对数据进行分析处理建模。
画出原始数据的时间路径图
#-*- coding:utf-8 -*-
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.graphics.tsaplots import acf,pacf,plot_acf,plot_pacf
from statsmodels.tsa.arima_model import ARMA
time_series = pd.Series([151.0, 188.46, 199.38, 219.75, 241.55, 262.58, 328.22, 396.26, 442.04, 517.77, 626.52, 717.08, 824.38, 913.38, 1088.39, 1325.83, 1700.92, 2109.38, 2499.77, 2856.47, 3114.02, 3229.29, 3545.39, 3880.53, 4212.82, 4757.45, 5633.24, 6590.19, 7617.47, 9333.4, 11328.92, 12961.1, 15967.61])
time_series.index = pd.Index(sm.tsa.datetools.dates_from_range('1978','2010'))
time_series.plot(figsize=(12,8))
plt.show()
由上图我们可以看出，这个时间序列是呈指数形式的，波动性比较大，不是稳定的时间序列，一般对于这种指数形式的数据，可以对其取对数，将其转化为线性趋势。
time_series = np.log(time_series)
time_series.plot(figsize=(8,6))
plt.show()
由上图可以看出，去了对数之后的时间路径图明显具有线性趋势，为了确定其稳定性，对取对数后的数据进行 adf 检验
t=sm.tsa.stattools.adfuller(time_series, )
output=pd.DataFrame(index=['Test Statistic Value', "p-value", "Lags Used", "Number of Observations Used","Critical Value(1%)","Critical Value(5%)","Critical Value(10%)"],columns=['value'])
output['value']['Test Statistic Value'] = t[0]
output['value']['p-value'] = t[1]
output['value']['Lags Used'] = t[2]
output['value']['Number of Observations Used'] = t[3]
output['value']['Critical Value(1%)'] = t[4]['1%']
output['value']['Critical Value(5%)'] = t[4]['5%']
output['value']['Critical Value(10%)'] = t[4]['10%']
print(output)
检验结果如下
| 检验项 | 检验结果 |
|---|---|
| Test Statistic Value | 0.807369 |
| p-value | 0.991754 |
| Lags Used | 1 |
| Number of Observations Used | 31 |
| Critical Value(1%) | -3.66143 |
| Critical Value(5%) | -2.96053 |
| Critical Value(10%) | -2.61932 |
由上表可知，t 统计量要大于任何置信度的临界值，因此认为该序列是非平稳的，所以再对序列进行差分处理，发现差分之后的序列基本达到稳定，如下图所示，并且通过了 ADF 检验，检验结果见下表。
time_series = time_series.diff(1)
time_series = time_series.dropna(how=any)
time_series.plot(figsize=(8,6))
plt.show()
t=sm.tsa.stattools.adfuller(time_series)
output=pd.DataFrame(index=['Test Statistic Value', "p-value", "Lags Used", "Number of Observations Used","Critical Value(1%)","Critical Value(5%)","Critical Value(10%)"],columns=['value'])
output['value']['Test Statistic Value'] = t[0]
output['value']['p-value'] = t[1]
output['value']['Lags Used'] = t[2]
output['value']['Number of Observations Used'] = t[3]
output['value']['Critical Value(1%)'] = t[4]['1%']
output['value']['Critical Value(5%)'] = t[4]['5%']
output['value']['Critical Value(10%)'] = t[4]['10%']
print(output)
| 检验项 | 检验结果 |
|---|---|
| Test Statistic Value | -3.52276 |
| p-value | 0.00742139 |
| Lags Used | 0 |
| Number of Observations Used | 31 |
| Critical Value(1%) | -3.66143 |
| Critical Value(5%) | -2.96053 |
| Critical Value(10%) | -2.61932 |
确定自相关系数和平均移动系数（p,q）
根据时间序列的识别规则，采用 ACF 图、PAC 图，AIC 准则（赤道信息量准则）和 BIC 准则（贝叶斯准则）相结合的方式来确定 ARMA 模型的阶数, 应当选取 AIC 和 BIC 值达到最小的那一组为理想阶数。
plot_acf(time_series)
plot_pacf(time_series)
plt.show()
r,rac,Q = sm.tsa.acf(time_series, qstat=True)
prac = pacf(time_series,method='ywmle')
table_data = np.c_[range(1,len(r)), r[1:],rac,prac[1:len(rac)+1],Q]
table = pd.DataFrame(table_data, columns=['lag', "AC","Q", "PAC", "Prob(>Q)"])
print(table）
根据上面的几个图，我们可以先取 p=1, q=2。进行模型估计,结果见下图。
p,d,q = (1,1,2)
arma_mod = ARMA(time_series,(p,d,q)).fit(disp=-1,method='mle')
summary = (arma_mod.summary2(alpha=.05, float_format="%.8f"))
print(summary)
这里的 p和q 参数可以调整，然后找出最佳的（AIC最小，BIC最小），经过比较， p=0,q=1 为理想阶数。
这里有一个自动取 p和q 的函数，如果要自动定阶的话，可以采用
(p, q) =(sm.tsa.arma_order_select_ic(dta,max_ar=3,max_ma=3,ic='aic')['aic_min_order'])
#这里需要设定自动取阶的 p和q 的最大值，即函数里面的max_ar,和max_ma。ic 参数表示选用的选取标准，这里设置的为aic,当然也可以用bic。然后函数会算出每个 p和q 组合(这里是(0,0)~(3,3)的AIC的值，取其中最小的,这里的结果是(p=0,q=1)。
残差和白噪声检验
个人感觉这个就是对模型 ARIMA(0,1,1) 的残差序列 arma_mod.resid 进行 ADF 检验。
arma_mod = ARMA(time_series,(0,1,1)).fit(disp=-1,method='mle')
resid = arma_mod.resid
t=sm.tsa.stattools.adfuller(resid)
output=pd.DataFrame(index=['Test Statistic Value', "p-value", "Lags Used", "Number of Observations Used","Critical Value(1%)","Critical Value(5%)","Critical Value(10%)"],columns=['value'])
output['value']['Test Statistic Value'] = t[0]
output['value']['p-value'] = t[1]
output['value']['Lags Used'] = t[2]
output['value']['Number of Observations Used'] = t[3]
output['value']['Critical Value(1%)'] = t[4]['1%']
output['value']['Critical Value(5%)'] = t[4]['5%']
output['value']['Critical Value(10%)'] = t[4]['10%']
print(output)
#结果如下
Test Statistic Value -3.114
p-value 0.025534
Lags Used 1
Number of Observations Used 30
Critical Value(1%) -3.66992
Critical Value(5%) -2.96407
Critical Value(10%) -2.62117
当然这里也可以画出 acf 图和 pacf 图。
模型预测
arma_model = sm.tsa.ARMA(time_series,(0,1)).fit(disp=-1,maxiter=100)
predict_data = arma_model.predict(start=str(1979), end=str(2010+3), dynamic = False)
预测结果还原
对预测出来的数据，进行逆差分操作（由原始数据取对数后的数据加上预测出来的数据），然后再取指数即可还原。
| 年份 | 2011年 | 2012年 | 2013年 |
|---|---|---|---|
| 实际值 | 19632.26 | 22250.45 | 24668.49 |
| 预测值 | 19314.03 | 22415.10 | 26014.08 |
上图最后3个为预测值，然后查询2011年到2013年湖北GDP的实际值，可以进行对照
| 年份 | 2011年 | 2012年 | 2013年 |
|---|---|---|---|
| 实际值 | 19632.26 | 22250.45 | 24668.49 |
| 预测值 | 19314.03 | 22415.10 | 26014.08 |
小结
从预测对结果看，2011年到2013年的预测结果和实际的差别不大。这个模型在短期预测结果比较好。模型处理主要还是应用了Python 第三方库 statsmodels 中的模型算法，其中还有很多细节，可以查阅相关文档，这里只是简单的应用了一下，由于代码都是一小段一小段写的，很乱，只提供了一些片段供参考。
参考资料
python时间序列分析
时间序列实战（一）
用ARIMA模型做需求预测
python 时间序列分析之ARIMA
ARIMA模型文档
作者：熙淺
链接：https://www.jianshu.com/p/4130bac8ebec
來源：简书
简书著作权归作者所有，任何形式的转载都请联系作者获得授权并注明出处。
