在有些应用中，我们需要高阶的AR或MA模型才能充分地描述数据的动态结构，这样问题会变得很繁琐。为了克服这个困难，提出了自回归滑动平均(ARMA)模型。
基本思想是把AR和MA模型结合在一起，使所使用的参数个数保持很小。
模型的形式为：$$ \large r_t = \phi_0 + \sum_{i=1}^p\phi_ir_{t-i} + a_t + \sum_{i=1}^q\theta_ia_{t-i} $$
其中，{$a_t$}为白噪声序列，p和q都是非负整数。AR和MA模型都是ARMA(p,q)的特殊形式。
利用向后推移算子$B$，上述模型可写成: $$ \large (1-\phi_1B - \cdots -\phi_pB^p)r_t = \phi_0 + (1-\theta_1B-\cdots - \theta_qB^q)a_t$$
（后移算子$B$，即上一时刻）
这时候我们求$r_t$的期望，得到：$$ \large E(r_t) = \frac{\phi_0}{1-\phi_1-\cdots - \phi_p}$$
和上期我们的AR模型一模一样。因此有着相同的特征方程：$$\large 1 - \phi_1x - \phi_2x^2 - \cdots -\phi_px^p = 0$$
有一点很关键：ARMA模型的应用对象应该为平稳序列 我们下面的步骤都是建立在假设原序列平稳的条件下的
from scipy import stats
import statsmodels.api as sm # 统计相关的库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
raw_data = D.history_data(['000001.SZA'],'2015-01-01','2018-01-01',['close'])[['date','close']].set_index('date')
data = np.array(raw_data.pct_change().dropna()) # 平安银行日涨跌
fig = plt.figure(figsize=(20,10))
ax1=fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(data,lags=30,ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(data,lags=30,ax=ax2)
可以看出，模型的阶次应该为(28,28)。然而，这么高的阶次建模的计算量是巨大的。
为什么不再限制滞后阶数小一些？如果lags设置为25，20或者更小时，阶数为(0,0)，显然不是我们想要的结果。
综合来看，由于计算量太大，在这里就不使用(28,28)建模了。采用另外一种方法确定阶数。
关于信息准则，上一期其实有过一些介绍：
目前选择模型常用如下准则： （其中L为似然函数，k为参数数量，n为观察数）
AIC = -2 ln(L) + 2 k 中文名字：赤池信息量 akaike information criterion
BIC = -2 ln(L) + ln(n)*k 中文名字：贝叶斯信息量 bayesian information criterion
HQ = -2 ln(L) + ln(ln(n))*k hannan-quinn criterion
我们常用的是AIC准则，AIC鼓励数据拟合的优良性但是尽量避免出现过度拟合(Overfitting)的情况。所以优先考虑的模型应是AIC值最小的那一个模型。
下面，我们分别应用以上3种法则，为我们的模型定阶，数据仍然是平安银行日涨跌幅序列：
为了控制计算量，我们限制AR最大阶不超过6，MA最大阶不超过4。 但是这样带来的坏处是可能为局部最优。
print("AIC", sm.tsa.arma_order_select_ic(data,max_ar=6,max_ma=4,ic='aic')['aic_min_order']) # AIC
print("BIC", sm.tsa.arma_order_select_ic(data,max_ar=6,max_ma=4,ic='bic')['bic_min_order']) # BIC
print("HQIC", sm.tsa.arma_order_select_ic(data,max_ar=6,max_ma=4,ic='hqic')['hqic_min_order']) # HQIC
可以看出，三种方法下求解的模型阶次为(3,2)。
我们使用上一节AIC准则求解的模型阶次(3,2)来建立ARMA模型，源数据为平安银行日涨跌幅数据，最后10个数据用于预测。
order = (3,2)
train = data[:-10]
test = data[-10:]
tempModel = sm.tsa.ARMA(train,order).fit()
同样的，先来看看拟合效果：
delta = tempModel.fittedvalues - train
score = 1 - delta.var()/train.var()
print('拟合优度：',score)
如果对比之前建立的AR、MA模型，可以发现拟合精度上有所提升，但仍然是不够看的级别~
接着来看预测效果：
predicts = tempModel.predict(len(train), len(train)+9, dynamic=True)
print(len(predicts))
comp = pd.DataFrame({'original':[i[0] for i in np.reshape(test,(10,1)).tolist()],'predict':predicts.tolist()})
T.plot(comp,title='模型预测与真实值')
可以看出，虽然还是准头很低，不过相比之前的MA模型，只看涨跌的话，胜率为50%。
到目前为止，我们研究的序列都集中在平稳序列。，即ARMA模型研究的对象为平稳序列。如果序列是非平稳的，就可以考虑使用ARIMA模型。
ARIMA比ARMA仅多了个"I"，代表着其比ARMA多一层内涵：也就是差分
一个非平稳序列经过d次差分后，可以转化为平稳时间序列。d 具体的取值，我们得分被对差分1次后的序列进行平稳性检验，若果是非平稳的，则继续差分。直到d次后检验为平稳序列。
ADF是一种常用的单位根检验方法，它的原假设为序列具有单位根，即非平稳，对于一个平稳的时序数据，就需要在给定的置信水平上显著，拒绝原假设。
下面给出示例：我们先看平安银行的日指数序列：
data2 = raw_data['close'] # 平安银行
data2.plot(figsize=(15,5))
看图形，这里显然是非平稳的。接着我们进行ADF单位根检验。
temp = np.array(data2)
t = sm.tsa.stattools.adfuller(temp) # ADF检验
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%']
output
可以看出，p-value为0.553824，大于显著性水平。原假设：序列具有单位根即非平稳。不能被拒绝。 因此平安银行日指数序列为非平稳的。
我们将序列进行1次差分后再次检验！
data2Diff = data2.diff(1) # 差分
T.plot(data2Diff, title='一阶差分')
从图形来看，序列近似平稳序列，我们来进行ADF检验：
temp = np.array(data2Diff)[1:] # 差分后第一个值为NaN,舍去
t = sm.tsa.stattools.adfuller(temp) # ADF检验
print("p-value: ", t[1])
可以看出，p-value非常接近于0，拒绝原假设，因此，该序列为平稳的。
可见，经过1次差分后的序列是平稳的，对于原序列，d的取值为1即可。
上面一小节我们确定了差分次数d，接下来，我们就可以将差分后的序列建立ARMA模型。
首先，我们还是尝试PACF和ACF来判断p、q。
temp = np.array(data2Diff)[1:] # 差分后第一个值为NaN,舍去
fig = plt.figure(figsize=(20,10))
ax1=fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(temp,lags=30,ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(temp,lags=30,ax=ax2)
可以看出，模型的阶次为(20,20)，还是太高了。建模计算了太大。我们再看看AIC准则
sm.tsa.arma_order_select_ic(temp,max_ar=6,max_ma=4,ic='aic')['aic_min_order'] # AIC
根据AIC准则，差分后的序列的ARMA模型阶次为(3,2)。因此，我们要建立的ARIMA模型阶次(p,d,q) = (3,1,2)
根据上一结确定的模型阶次，我们对差分后序列建立ARMA(3,2)模型：
order = (3,2)
data = np.array(data2Diff)[1:] # 差分后，第一个值为NaN
train = data[:-10]
test = data[-10:]
model = sm.tsa.ARMA(train,order).fit()
我们先看差分序列的ARMA拟合值。
T.plot(pd.DataFrame({'real value':train,'fitted value':model.fittedvalues}),title='训练误差')
delta = model.fittedvalues - train
score = 1 - delta.var()/train[1:].var()
print('拟合优度： ',score)
再看对差分序列的预测情况：
predicts = model.predict(len(train),len(train)+9, dynamic=True)
print(len(predicts))
T.plot(pd.DataFrame({'real value':test,'fitted value':predicts}),title='测试误差')
可以看出，差分序列ARMA模型的拟合效果和预测结果并不好，预测值非常小，这代表什么？这代表对于新的值，这里认为它很接近上一时刻的值。
这个影响可能来自模型阶次，看来模型阶次还是得尝试更高阶的。这里就不建模了，（计算时间太长），大家有兴趣可以试试高阶的模型。
1.《金融时间序列分析》 第2版 Ruey S.Tsay著 王辉、潘家柱 译
2.Time Series analysis tsa http://nipy.bic.berkeley.edu/nightly/statsmodels/doc/html/tsa.html