原文链接：点这里
在了解了AR和MA模型后，我们将进一步学习结合了这两者的ARIMA模型，ARIMA在时间序列数据分析中有着非常重要的地位。但在这之前，让我们先来看ARIMA的简化版ARMA模型（Autoregressive moving average model），ARMA同样是结合了AR和MA模型，公式如下：
yt = c + φ1*yt-1 + φ2*yt-2 + ... +φp*yt-p + et + θ1*et-1 + θ2*et-2 + ... +θq*et-q
可以看出，ARMA模型就是AR和MA的简单结合，同时包含了历史数值项和错误项。由于AR模型对时间序列有平稳性要求，ARMA模型也存在这个限制，因此我们将其拓展到ARIMA模型，引入的差分概念是一种获得时间序列的方法。最常使用的一种差分方法是计算当前项和前项的差值，获得一组新的时间序列。由于ARIMA强大的功能，这使得在统计领域有着非常广泛的应用，尤其适用于机器学习和深度学习难以应用的小样本数据集上，同时也应当警惕过拟合的情况发生。
当我们拿到一个陌生的数据集时，可能会思考应当使用AR，MA，还是结合AR和MA的模型，一个简单的准则是画出ACF和PACF曲线，根据其特征来判断。
| 函数 | AR(p) | MA(q) | ARMA or ARIMA |
| ACF曲线 | 缓缓下降 | 在间隔=q后突然下降 | 没有明显截止点 |
| PACF曲线 | 在间隔=p后突然下降 | 缓缓下降 | 没有明显截止点 |
在实际分析过程中，我们也可以一律使用ARIMA模型，因为AR,MA,ARMA都是它的一种特殊情况。
ARIMA有三个参数p、d、q，写作ARIMA（p,d,q），其中p代表AR(p)，自回归阶数，d代表Integrated (d)，差分阶数，q代表MA(q)，移动平均阶数。我们应当确保这三个参数尽可能的小避免过拟合，一个可供参数的准则是，不要让d超过2，p和q超过5，并且p和q尽量保证一个是模型主导项，另一个相对较小。
ARIMA有一些经典的特殊例子：
- ARIMA(0,0,0) 是白噪声模型
- ARIMA(0,1,0) 是随机行走模型
- ARIMA(0,1,1) 是指数平滑模型
- ARIMA(0,2,2) 是Holt线性趋势模型
实战
如何确定ARIMA模型的参数通常有两种方法，手动拟合法和自动拟合法。
手动拟合法是一种经验方法，最著名的方法被称为Box-Jenkins method，这是一个多步迭代过程，分为以下几个步骤：
- 结合对数据集的可视化和领域知识选择一组合适的模型
- 拟合模型
- 根据模型表现微调参数，以解决模型表现不佳的地方
通过一个简单的例子来学习手动box-jenkins方法的使用过程。
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
from pmdarima.arima import auto_arima
from pmdarima.arima import ADFTest
import statsmodels.api as sm
# 模拟ARMA数据，在完成最终拟合前，请假装不知道样本时间序列的参数
ar = np.array([1, -0.8, 0.4])
ma = np.array([1, -0.7])
arma_process = sm.tsa.ArmaProcess(ar, ma)
y = arma_process.generate_sample(1000)
plt.plot(y)
plt.show()
# # step1 画出acf,pacf曲线，根据经验规则，没有明显的cutoff，模型应当包含AR和MA项
# plot_acf(y)
# plot_pacf(y)
# plt.show()
# # step2 用最简单的arima(1,0,1)尝试拟合
# mod = ARIMA(y, order=(1, 0, 1))
# res = mod.fit()
#
# # step3 画出拟合残差的acf和pacf
# plot_acf(res.resid)
# plot_pacf(res.resid)
# plt.show()
# # step4 从上面对residual绘制acf和pacf可以看到，还存在有显著的自相关性，说明说明最开始尝试的ARIMA(1,0,1)模型并没有充分表达原数据
# # 尝试增加模型复杂度
# a2m1 = ARIMA(y, order=(2, 0, 1)).fit()
# plot_acf(a2m1.resid)
# plot_pacf(a2m1.resid)
# plt.show()
# step5 通过将AR项从一阶改成二阶，这个更复杂的模型取得了不错的效果，残差的acf和pacf几乎看不到其他相关性存在
# 为了验证是否还能获得更好的模型，继续修改参数进行测试，一个简单快速的验证方式是计算真实值和预测值的相关系数
def model_performance(p, d, q):
model = ARIMA(y, order=(p, d, q)).fit()
df = pd.DataFrame({'predict': model.predict(), 'true': y})
corr = df.corr().loc['predict', 'true']
print("ARIMA({0},{1},{2}) performance: {3} ".format(p, d, q, corr))
model_performance(1, 0, 1)
model_performance(2, 0, 1)
model_performance(1, 0, 2)
model_performance(2, 0, 2)
model_performance(2, 0, 3)
# ARIMA(1,0,1) performance: 0.2979967322503139
# ARIMA(2,0,1) performance: 0.45403917281587897
# ARIMA(1,0,2) performance: 0.4399028235919454
# ARIMA(2,0,2) performance: 0.4560342835471262
# ARIMA(2,0,3) performance: 0.46020925616191444
print(ARIMA(y, order=(2, 0, 1)).fit().params)
# [-0.0094566 0.80238963 -0.41757641 -0.68993587 0.96598217]
# 和我们最开始创建的已经很接近了可以看到ARIMA(2,0,1)的效果有了明显提升，而当再增加模型复杂度时，效果没有明显提升。当然这只是一个极为粗略的手动拟合的例子，真实场景则需要更为精细的优化。
从上面的例子中可以看出手动拟合法较为复杂，对经验要求很高，一些统计学家也长期致力于这类方法的优化思路。但同时也受到了一些人的批评，认为这种纯手工的方法效率太低，而且很依赖经验，不是最优的方法。
下面就来学习如何进行自动拟合。
自动拟合法
# 自动拟合法
# step1 准备数据
sales_data = pd.read_csv('../../srouce/data/retail_sales.csv')
sales_data['date'] = pd.to_datetime(sales_data['date'])
sales_data.set_index('date', inplace=True)
sales_data.plot()
plt.show()
# step2 检查平稳性
adf_test = ADFTest()
print(adf_test.should_diff(sales_data))
# (0.01, False)
# 结果表明不平稳，提示我们需要引入差分项
# step3 划分训练集和测试集
train = sales_data[:60]
test = sales_data[60:]
# step4 拟合模型
arima_model = auto_arima(train, start_p=0, d=1, start_q=0, max_p=5, max_d=5, max_q=5,
start_P=0, D=1, start_Q=0, max_P=5, max_D=5, max_Q=5, m=12,
seasonal=True, trace=True, n_fits=50)
# 查看模型结果
print(arima_model.summary())
# step5 预测时间序列，并和真实值比较
pred = pd.DataFrame(arima_model.predict(n_periods=12), columns=['predicted'], index=test.index)
plt.plot(train)
plt.plot(test, label='true')
plt.plot(pred, label='predict')
plt.legend()
plt.show()
