前言

上次介绍了时间序列的一些基础小知识和AR模型。部分小伙伴指出公式还是多了些，所以本篇开始将删繁趋简，重点介绍如何用python实现这些模型~ 关于公式的推导和理解，有兴趣的可以参阅参考文献

参考文献

1.《金融时间序列分析》 第2版 Ruey S.Tsay著 王辉、潘家柱 译

2.Time Series analysis tsa

三、滑动平均(MA)模型

这里我们直接给出MA(q)模型的形式：

c0为一个常数项。这里的at,是AR模型t时刻的扰动或者说新息，则可以发现，该模型，使用了过去q个时期的随机干扰或预测误差来线性表达当前的预测值。

3.1 MA模型的性质3.1.1 平稳性

MA模型总是弱平稳的，因为他们是白噪声序列（残差序列）的有限线性组合。因此，根据弱平稳的性质可以得出两个结论：

3.1.2 自相关函数

对q阶的MA模型，其自相关函数ACF总是q步截尾的。因此MA(q)序列只与其前q个延迟值线性相关，从而它是一个“有限记忆”的模型。

这一点可以用来确定模型的阶次，后面会介绍。

3.1.3 可逆性

当满足可逆条件的时候，MA(q)模型可以改写为AR(p)模型。这里不进行推导，给出1阶和2阶MA的可逆性条件。

1阶：

2阶：

3.2 MA的阶次判定

我们通常利用上面介绍的第二条性质：MA(q)模型的ACF函数q步截尾来判断模型阶次。示例如下：

使用上证指数的日涨跌数据（2013年1月至2014年8月）来进行分析，先取数据：

from scipy import  statsimport statsmodels.api as sm  # 统计相关的库import numpy as npimport pandas as pdimport matplotlib.pyplot as pltIndexData = DataAPI.MktIdxdGet(indexID=u"",ticker=u"000001",beginDate=u"20130101",endDate=u"20140801",field=u"tradeDate,closeIndex,CHGPct",pandas="1") IndexData = IndexData.set_index(IndexData['tradeDate'])data = np.array(IndexData['CHGPct']) # 上证指数日涨跌IndexData['CHGPct'].plot(figsize=(15,5))

可以看出，序列看上去是弱平稳的。下面我们画出序列的ACF：

fig = plt.figure(figsize=(20,5))ax1=fig.add_subplot(111)fig = sm.graphics.tsa.plot_acf(data,ax=ax1)

我们发现ACF函数在43处截尾，之后的acf函数均在置信区间内，我们判定该序列MA模型阶次为43阶。

3.3 建模和预测

由于sm.tsa中没有单独的MA模块，我们利用ARMA模块，只要将其中AR的阶p设为0即可。

函数sm.tsa.ARMA中的输入参数中的order(p,q)，代表了AR和MA的阶次。模型阶次增高，计算量急剧增长，因此这里就建立10阶的模型作为示例，如果按上一节的判断阶次来建模，计算时间过长。

我们用最后10个数据作为out-sample的样本，用来对比预测值。

order = (0,10)train = data[:-10]test = data[-10:]tempModel = sm.tsa.ARMA(train,order).fit()

我们先来看看拟合效果,计算

delta = tempModel.fittedvalues - trainscore = 1 - delta.var()/train.var()print score

0.0278740087833

可以看出，score远小于1，拟合效果不好。

然后我们用建立的模型进行预测最后10个数据

predicts = tempModel.predict(371, 380, dynamic=True)print len(predicts)comp = pd.DataFrame()comp['original'] = testcomp['predict'] = predictscomp.plot()

10

可以看出，建立的模型效果很差，预测值明显小了1到2个数量级！就算只看涨跌方向，正确率也不足50%。该模型不适用于原数据。

没关系，如果真的这么简单能预测指数日涨跌才奇怪了~

关于MA的内容只做了简单介绍，下面主要介绍ARMA模型

四、ARMA模型

在有些应用中，我们需要高阶的AR或MA模型才能充分地描述数据的动态结构，这样问题会变得很繁琐。为了克服这个困难，提出了自回归滑动平均(ARMA)模型。

基本思想是把AR和MA模型结合在一起，使所使用的参数个数保持很小。

模型的形式为：

其中，{at}为白噪声序列，p和q都是非负整数。AR和MA模型都是ARMA(p,q)的特殊形式。

利用向后推移算子B，上述模型可写成:

（后移算子B，即上一时刻）

这时候我们求rt的期望，得到：

和上期我们的AR模型一毛一样。因此有着相同的特征方程：

该方程所有解的倒数称为该模型的特征根，如果所有的特征根的模都小于1，则该ARMA模型是平稳的。

有一点很关键：ARMA模型的应用对象应该为平稳序列！ 我们下面的步骤都是建立在假设原序列平稳的条件下的~

4.1 识别ARMA模型阶次4.1.1 PACF、ACF 判断模型阶次

我们通过观察PACF和ACF截尾，分别判断p、q的值。(限定滞后阶数50)

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)

可以看出，模型的阶次应该为(27,27)。然而，这么高的阶次建模的计算量是巨大的。

为什么不再限制滞后阶数小一些？如果lags设置为25，20或者更小时，阶数为(0,0)，显然不是我们想要的结果。

综合来看，由于计算量太大，在这里就不使用(27,27)建模了。采用另外一种方法确定阶数。

4.2 信息准则定阶

关于信息准则，上一期其实有过一些介绍：

目前选择模型常用如下准则：（其中L为似然函数，k为参数数量，n为观察数）

我们常用的是AIC准则，AIC鼓励数据拟合的优良性但是尽量避免出现过度拟合(Overfitting)的情况。所以优先考虑的模型应是AIC值最小的那一个模型。

下面，我们分别应用以上3种法则，为我们的模型定阶，数据仍然是上证指数日涨跌幅序列：

为了控制计算量，我们限制AR最大阶不超过6，MA最大阶不超过4。 但是这样带来的坏处是可能为局部最优。

sm.tsa.arma_order_select_ic(data,max_ar=6,max_ma=4,ic='aic')['aic_min_order']  # AIC

(3, 3)

sm.tsa.arma_order_select_ic(data,max_ar=6,max_ma=4,ic='bic')['bic_min_order']  # BIC

(0, 0)

result = sm.tsa.arma_order_select_ic(data,max_ar=6,max_ma=4,ic='hqic')['hqic_min_order'] # HQIC

可以看出，AIC求解的模型阶次为(3,3)。我们这里就以AIC准则为准~ 至于到底哪种准则更好，小伙伴们可以分别建模进行对比~

4.2 模型的建立及预测

我们使用上一节AIC准则求解的模型阶次(3,3)来建立ARMA模型，源数据为上证指数日涨跌幅数据，最后10个数据用于预测。

order = (3,3)train = data[:-10]test = data[-10:]tempModel = sm.tsa.ARMA(train,order).fit()

同样的，先来看看拟合效果：

delta = tempModel.fittedvalues - trainscore = 1 - delta.var()/train.var()print score

0.0492668871573

如果对比之前建立的AR、MA模型，可以发现拟合精度上有所提升，但仍然是不够看的级别~

接着来看预测效果：

predicts = tempModel.predict(371, 380, dynamic=True)print len(predicts)comp = pd.DataFrame()comp['original'] = testcomp['predict'] = predictscomp.plot()

10

可以看出，虽然还是准头很低，不过相比之前的MA模型，只看涨跌的话，胜率为55.6%，效果还是好了不少~

五、 ARIMA模型

到目前为止，我们研究的序列都集中在平稳序列。即ARMA模型研究的对象为平稳序列。如果序列是非平稳的，就可以考虑使用ARIMA模型。

ARIMA比ARMA仅多了个"I"，代表着其比ARMA多一层内涵：也就是差分

一个非平稳序列经过d次差分后，可以转化为平稳时间序列。d具体的取值，我们得分被对差分1次后的序列进行平稳性检验，若果是非平稳的，则继续差分。直到d次后检验为平稳序列。

5.1 单位根检验

ADF是一种常用的单位根检验方法，它的原假设为序列具有单位根，即非平稳，对于一个平稳的时序数据，就需要在给定的置信水平上显著，拒绝原假设。

下面给出示例：我们先看上证综指的日指数序列：

data2 = IndexData['closeIndex'] # 上证指数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.1704489，大于显著性水平。原假设：序列具有单位根即非平稳。不能被拒绝。因此上证指数日指数序列为非平稳的。

我们将序列进行1次差分后再次检验！

data2Diff = data2.diff()  # 差分data2Diff.plot(figsize=(15,5))

从图形来看，序列近似平稳序列，我们来进行ADF检验：

temp = np.array(data2Diff)[1:] # 差分后第一个值为NaN,舍去t = sm.tsa.stattools.adfuller(temp)  # ADF检验print "p-value:   ",t[1]

p-value: 2.31245750144e-30

可以看出，p-value非常接近于0，拒绝原假设，因此，该序列为平稳的。

可见，经过1次差分后的序列是平稳的，对于原序列，d的取值为1即可。

5.2 ARIMA(p,d,q)模型阶次确定

上面一小节我们确定了差分次数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)

可以看出，模型的阶次为(27,27)，还是太高了。建模计算了太大。我们再看看AIC准则

sm.tsa.arma_order_select_ic(temp,max_ar=6,max_ma=4,ic='aic')['aic_min_order']  # AIC

(2, 2)

根据AIC准则，差分后的序列的ARMA模型阶次为(2,2)。因此，我们要建立的ARIMA模型阶次(p,d,q) = (2,1,2)

5.3 ARIMA模型建立及预测

根据上一结确定的模型阶次，我们对差分后序列建立ARMA(2,2)模型：

order = (2,2)data = np.array(data2Diff)[1:] # 差分后，第一个值为NaNrawdata = np.array(data2)train = data[:-10]test = data[-10:]model = sm.tsa.ARMA(train,order).fit()

我们先看差分序列的ARMA拟合值。

plt.figure(figsize=(15,5))plt.plot(model.fittedvalues,label='fitted value')plt.plot(train[1:],label='real value')plt.legend(loc=0)

delta = model.fittedvalues - trainscore = 1 - delta.var()/train[1:].var()print score

0.0397490005618

再看对差分序列的预测情况：

predicts = model.predict(10,381, dynamic=True)[-10:]print len(predicts)comp = pd.DataFrame()comp['original'] = testcomp['predict'] = predictscomp.plot(figsize=(8,5))

10

可以看出，差分序列ARMA模型的拟合效果和预测结果并不好，预测值非常小，这代表什么？这代表对于新的值，这里认为它很接近上一时刻的值。

这个影响可能来自模型阶次，看来模型阶次还是得尝试更高阶的。这里就不建模了，（计算时间太长），大家有兴趣可以试试高阶的模型。

然后我们将预测值还原（即在上一时刻指数值的基础上加上差分差值的预估）：

rec = [rawdata[-11]]pre = model.predict(371, 380, dynamic=True) # 差分序列的预测for i in range(10):    rec.append(rec[i]+pre[i])

plt.figure(figsize=(10,5))plt.plot(rec[-10:],'r',label='predict value')plt.plot(rawdata[-10:],'blue',label='real value')    plt.legend(loc=0)

我们发现，由于对差分序列的预测很差，还原到原序列后，预测值几乎在预测前一个值上小幅波动。

小结一下：

目前为止，我们接触的模型其实核心就是ARMA模型，其他模型都是其特殊形式，目前为止，对指数日涨跌序列或者指数日序列的预测都没有特别好的效果。这类模型的阶次确立以及建立高阶模型是主要难点。