时间序列简介

在数学上，随机过程被定义为一族时间随机变量，即{x(t),t∈T},其中T表示时间t的变动范围。当T={0,±1,±2,…}时，此类随机过程x(t)是离散时间t的随机函数，称为时间序列。时间序列的构成要素有：长期趋势，季节变动，循环变动，不规则变动：

长期趋势(T)是在较长时期内受某种根本性因素作用而形成的总的变动趋势
季节变动(S)是在一年内随着季节的变化而发生的有规律的周期性变动。它是诸如气候条件、生产条件、节假日或人们的风俗习惯等各种因素影响的结果。
循环变动(C)是时间序列呈现出得非固定长度的周期性变动。循环波动的周期可能会持续一段时间，但与趋势不同，它不是朝着单一方向的持续变动，而是涨落相同的交替波动。
不规则变动(I)是时间序列中除去趋势、季节变动和周期波动之后的随机波动。不规则波动通常总是夹杂在时间序列中，致使时间序列产生一种波浪形或震荡式的变动。只含有随机波动的序列也称为平稳序列。

序列的平稳性

在百度词条中是这样粗略的讲的：平稳时间序列粗略地讲，一个时间序列，如果均值没有系统的变化（无趋势）、方差没有系统变化，且严格消除了周期性变化，就称之是平稳的。

时间序列的平稳性，和其它随机过程一样，分为严平稳和宽平稳。

严平稳：如果对所有的时刻t，任意正整数k和任意k个正整数, $(r_{t1},r_{t2},...,r_{tk})$ 的联合分布与 $(r_{t1+t},r_{t2+t},...,r_{tk+t})$ 的联合分布相同，我们称时间序列 ${r_t}$ 是严平稳的。也就是， $(r_{t1},r_{t2},...,r_{tk})$ 的联合分布在时间的平移变换下保持不变，这是个很强的条件。
宽平稳：若时间序列 $\{r_t\}$ 满足下面两个条件： $E(r_t)=\mu$ ， $\mu$ 是常数 $Cov(r_t,r_{t-l}) = \gamma _l$ ， $\gamma _l$ 只依赖于l则时间序列 $\{r_t\}$ 为弱平稳的。即该序列的均值，与 $r_{t−l}$ 的协方差不随时间而改变，l为任意整数。

由于严平稳的实现较为困难，所以一般来说平稳时间序列指的是宽平稳。如果是非平稳时间序列，可以通过多次差分的方法处理为宽平稳序列。

差分，就是求时间序列 $\{r_t\}$ 在t时刻的值 r_t 与t−1时刻的值 $r_{t−1}$ 的差 d_t ，则我们得到了一个新序列 $\{d_t\}$ ，为一阶差分，对新序列 $\{d_t\}$ 再做同样的操作，则为二阶差分。通常非平稳序列可以经过d次差分，处理成弱平稳或者近似弱平稳时间序列。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
 
import tushare as ts
pro = ts.pro_api()
df = pro.daily(ts_code='600519.SH',start_date='20180101') #茅台股价
 
df = df[['trade_date','close']]
df['trade_date'] = pd.to_datetime(df['trade_date'])
data = df.set_index('trade_date')
 
data['change_pct'] =''
for i in range(0,len(data['close'])-1):
    data.iloc[i, 1] = 100*(data['close'][i+1] - data['close'][i])/data['close'][i] # 计算每日收益率（涨跌幅）
data['change_pct'] = data['change_pct'][0:-1] 
data['close_diff_1'] = data['close'].diff(1)
data = data[0:-1] # 由于changepct值最后一行缺失，因此去除最后一行
data['change_pct'] = data['change_pct'].astype('float64') #使changepct数据类型和其他三列保持一致
 
data.plot(subplots=True,figsize=(12,12))

相关系数和自相关函数

相关系数

对于两个向量，我们希望定义它们是不是相关。一个很自然的想法，用向量与向量的夹角作为距离的定义，夹角小，就距离小，夹角大，就距离大。我们经常使用余弦公式来计算角度：

$\cos<\vec{a},\vec{b}>=\frac{\vec{a}\cdot \vec{b}}{\lvert \vec{a} \rvert \lvert \vec{b} \rvert}$

对于 $\vec{a}\cdot \vec{b}$ 我们叫做内积，例如：

$(x_1,y_1)\cdot (x_2,y_2)={x_1}{x_2}+{y_1}{y_2}$

相关系数的定义公式,X和Y的相关系数为：

$\rho_{xy} = \frac{Cov(X,Y)}{\sqrt{Var(X)Var(Y)}}$

根据样本的估计计算公式为：

$\rho_{xy} = \frac{\sum_{t=1}^{T}(x_t- \bar x)(y_t-\bar y)}{\sqrt{ \sum_{t=1}^{T}(x_t- \bar x)^{2}\sum_{t=1}^{T}(y_t- \bar y)^{2}}} = \frac{\overrightarrow{(X- \bar x)} \cdot \overrightarrow{(Y- \bar y)}}{\lvert \overrightarrow{(X- \bar x)} \rvert \lvert \overrightarrow{(Y- \bar y)} \rvert}$

相关系数实际上就是计算了向量空间中两个向量的夹角，协方差是去均值后两个向量的内积。

如果两个向量平行，相关系数等于1或者-1，同向的时候是1，反向的时候就是-1。如果两个向量垂直，则夹角的余弦就等于0，说明二者不相关。两个向量夹角越小，相关系数绝对值越接近1，相关性越高。只不过这里计算的时候对向量做了去均值处理，即中心化操作。而不是直接用向量, 计算。

自相关函数（Autocorrelation Function, ACF）

相关系数度量了两个向量的线性相关性，而在平稳时间序列 $\{r_t\}$ 中，我们有时候很想知道， r_t 与它的过去值 $r_{t-i}$ 的线性相关性。这时候我们把相关系数的概念推广到自相关系数。

r_t 与 $r_{t-l}$ 的相关系数称为 r_t 的间隔为的自相关系数，通常记为 $\rho_l$ ：

$\rho_l = \frac{Cov(r_t,r_{t-l})} {\sqrt{Var(r_t)Var(r_{t-l})}} = \frac{Cov(r_t,r_{t-l})} {Var(r_t)}$

这里用到了弱平稳序列的性质：

$Var(r_t) = Var(r_{t-l})$

对一个平稳时间序列的样本 $\{r_t\}$ ， $1\le t\le T$ ,则间隔为的样本自相关系数的估计为：

$\hat \rho_l = \frac{\sum_{t=l+1}^{T}(r_t- \bar r)(r_{t-l}-\bar r)}{ \sum_{t=1}^{T}(r_t- \bar r)^{2}}, 0 \leqslant l \leqslant T-1$

则函数 $\hat \rho_1,\hat \rho_2 ,\hat \rho_3...$ 称为 r_t 的样本自相关函数(ACF)

当自相关函数中所有的值都为0时，我们认为该序列是完全不相关的；因此，我们经常需要检验多个自相关系数是否为0。

混成检验

原假设：

$H0: \rho_1 = ... = \rho_m=0$

备择假设：

$H1: \exists i \in \{1,...,m\}, \rho_i \ne 0$

混成检验统计量：

$Q(m) = T(T+2) \sum_{l=1}^{m} \frac{\hat{\rho_l}^{2}}{T-l}$

Q(m) 渐进服从自由度为的 $\chi^2$ 分布

决策规则：

$Q(m) > \chi_\alpha^{2} ,拒绝H_0$

即, Q(m) 的值大于自由度为的卡方分布 $100(1-\alpha)$ 分位点时，我们拒绝 H_0 。

大部分软件会给出 Q(m) 的p-value，则当p-value小于等于显著性水平 $\alpha$ 时拒绝H0。

示例：

import statsmodels.api as sm
close_data = data['close'] # 上证指数每日收盘价
m = 10 # 检验10个自相关系数
 
acf,q,p = sm.tsa.acf(close_data, nlags=m, qstat = True)
 
# np.r_是按列连接两个矩阵，就是把两矩阵上下相加，要求列数相等，类似于pandas中的concat()。
# np.c_是按行连接两个矩阵，就是把两矩阵左右相加，要求行数相等，类似于pandas中的merge()。
out = np.c_[range(1,11), acf[1:], q, p]
 
output = pd.DataFrame(out, columns=['lag', 'ACF', 'Q', 'P-value'])
output = output.set_index('lag')
print(output)

lag        ACF            Q        P-value
1.0   0.963283   212.504727   3.903196e-48
2.0   0.928599   410.863584   6.054892e-90
3.0   0.893812   595.463129  9.703299e-129
4.0   0.858595   766.569712  1.337073e-164
5.0   0.820359   923.482644  2.201222e-197
6.0   0.778788  1065.538402  5.955998e-227
7.0   0.744518  1195.960070  5.278705e-254
8.0   0.712521  1315.960270  8.343090e-279
9.0   0.686764  1427.955165  7.044640e-302
10.0  0.658661  1531.448754   0.000000e+00

我们取显著性水平为0.05，可以看出，所有的p-value都小于0.05；则我们拒绝原假设 H_0 。因此，我们认为该序列是序列相关的。我们再来看看同期的日收益率序列：

change_data = data['change_pct']
m = 10 # 检验10个自相关系数
 
acf,q,p = sm.tsa.acf(change_data,nlags=m,qstat=True)  ## 计算自相关系数 及p-value
out = np.c_[range(1,11), acf[1:], q, p]
output=pd.DataFrame(out, columns=['lag', "AC", "Q", "P-value"])
output = output.set_index('lag')
print(output)

lag         AC         Q   P-value
1.0   0.010467  0.025092  0.874137
2.0  -0.018733  0.105814  0.948468
3.0  -0.015576  0.161875  0.983496
4.0  -0.023547  0.290575  0.990414
5.0   0.034672  0.570871  0.989298
6.0  -0.133740  4.760213  0.574915
7.0  -0.047866  5.299284  0.623491
8.0  -0.085952  7.045520  0.531729
9.0   0.064843  8.043923  0.529726
10.0 -0.032424  8.294713  0.600074

可以看出，p-value均大于显著性水平0.05。我们选择假设 H_0 ，即收益率序列没有显著的相关性。

白噪声序列和线性时间序列

白噪声序列

随机变量 X(t) (t=1,2,3…)，如果是由一个不相关的随机变量的序列构成的，即对于所有S不等于T，随机变量和的协方差为零，则称其为纯随机过程。

对于一个纯随机过程来说，若其期望和方差均为常数，则称之为白噪声过程。白噪声过程的样本实称成为白噪声序列，简称白噪声。之所以称为白噪声，是因为他和白光的特性类似，白光的光谱在各个频率上有相同的强度，白噪声的谱密度在各个频率上的值相同。

线性时间序列

时间序列{ r_t }，如果能写成：

$r_t = \mu + \sum_{i=0}^{\infty}\psi_ia_{t-i}$

$\mu$ 为 r_t 的均值， $psi_0=1,\{a_t\}$ 为白噪声序列，则我们称{ r_t } 为线性序列。其中 a_t 称为在时刻的新息(innovation)或扰动(shock)。

很多时间序列具有线性性，即是线性时间序列，相应的有很多线性时间序列模型，例如接下来要介绍的AR、MA、ARMA，都是线性模型，但并不是所有的金融时间序列都是线性的。对于弱平稳序列，我们利用白噪声的性质很容易得到 r_t 的均值和方差：

$E(r_t) = \mu , Var(r_t) = \sigma_a^2 \sum_{i=0}^{\infty} \psi_i^{2}$

$\sigma_a^2$ 为 a_t 的方差。因为 Var(r_t) 一定小于正无穷，因此 $\{\psi_i^2\}$ 必须是收敛序列，因此满足 $i \to \infty$ 时， $\psi_i^2 \to 0$ 。即，随着的增大，远处的扰动 $a_{t-i}$ 对 r_t 的影响会逐渐消失。

自回归(AR)模型

在上节中，我们计算了部分数据端的ACF，看表可知间隔为1时自相关系数是显著的。这说明在 t-1 时刻的数据 $r_{t-1}$ ，在预测时刻 r_t 时可能是有用的！

根据这点我们可以建立下面的模型：

$r_t = \phi_0 + \phi_1r_{t-1} + a_t$

其中{ a_t }是白噪声序列，这个模型与简单线性回归模型有相同的形式，这个模型也叫做一阶自回归(AR)模型，简称AR(1)模型。

从AR(1)很容易推广到AR(p)模型：

$r_t = \phi_0 + \phi_1r_{t-1} + \cdots + \phi_pr_{t-p}+ a_t \qquad (2.0)$

AR(p)模型的特征根及平稳性检验

我们先假定序列是弱平稳的，则有;

$E(r_t) = \mu \qquad Var(r_t) = \gamma_0 \qquad Cov(r_t,r_{t-j})=\gamma_j$

其中 $\mu$ , $\gamma_0$ 是常数。因为{ a_t }是白噪声序列，因此有：

$E(a_t) = 0, Var(a_t) = \sigma_a^2$

所以有：

$E(r_t) = \phi_0 + \phi_1E(r_{t-1}) + \phi_2E(r_{t-2}) + \cdots + \phi_pE(r_{t-p})$

根据平稳性的性质，又有 $E(r_t)=E(r_{t-1})=...=\mu$ ，从而：

$\mu = \phi_0 + \phi_1\mu + \phi_2\mu+ \cdots +\phi_p\mu\\$

$E(r_t)=\mu=\frac{\phi_0}{1-\phi_1 - \phi_2 - \cdots -\phi_p} \qquad$

公式中假定分母不为0，我们将下面的方程称为特征方程:

$1 - \phi_1x - \phi_2x^2 - \cdots -\phi_px^p = 0$

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

下面我们就用该方法，检验日收益率序列的平稳性：

change_data = data['change_pct'] # 载入收益率序列
change_data = np.array(change_data, dtype=np.float)
model = sm.tsa.AR(change_data)
results_AR = model.fit()
plt.figure(figsize=(12,6))
plt.plot(change_data, 'b', label = 'changepct')
plt.plot(results_AR.fittedvalues, 'r', label = 'AR model')
plt.legend()

我们可以看看模型有多少阶：

print(len(results_AR.roots)) # 模型阶数

可以看出，自动生成的AR模型是15阶的。关于阶次的讨论在下节内容，我们画出模型的特征根，来检验平稳性：

r1 = 1
theta = np.linspace(0, 2*np.pi, 360)
x1 = r1 * np.cos(theta)
y1 = r1 * np.sin(theta)
plt.figure(figsize=(6,6))
plt.plot(x1, y1, 'k')  # 画单位圆
roots = 1 / results_AR.roots  # 计算results_AR.roots后取倒数
for i in range(len(roots)):
    plt.plot(roots[i].real, roots[i].imag, '.r', markersize=8)  #画特征根
plt.show()

可以看出，所有特征根都在单位圆内，则序列为平稳的。

确定AR模型的阶数

对于线性回归模型的阶数，一般认为随着阶数的升高，会提高模型的精度。但是过高的阶数会造成模型过于复杂。下面介绍两种确定AR模型阶数的方法：偏相关函数（pACF）和信息准则函数（AIC、BIC、HQIC）。

偏相关函数(pACF)

首先引入阶数依次升高的AR模型：

$r_t = \phi_{0,1} + \phi_{1,1}r_{t-1} + \epsilon_{1t}$

$r_t = \phi_{0,2} + \phi_{1,2}r_{t-1} + \phi_{2,2}r_{t-2} + \epsilon_{2t}$

$r_t = \phi_{0,3} + \phi_{1,3}r_{t-1} + \phi_{2,3}r_{t-2} + \phi_{3,3}r_{t-3} + \epsilon_{3t}$

$r_t = \phi_{0,4} + \phi_{1,4}r_{t-1} + \phi_{2,4}r_{t-2} + \phi_{3,4}r_{t-3} + \phi_{4,4}r_{t-4} + \epsilon_{4t}$

$\dots \dots \dots \dots$

等式右边第一项为常数项，最后一项为误差项， $\hat{\phi}_{i,i}$ 为{rt}的间隔为i的样本偏自相关函数。 $\hat{\phi}_{i,i}$ 表示的是在第i-1行的基础上添加的一个阶数对{rt}的贡献。对于AR模型，当超过某个值p时，偏自相关函数应逐渐趋近于0，这既符合我们的直观认识，也有相应的理论证明可以找到。

对于偏相关函数的介绍，这里不详细展开，重点介绍一个性质：AR(p)序列的样本偏相关函数是p步截尾的。所谓截尾，就是快速收敛应该是快速的降到几乎为0或者在置信区间以内。

下面对2018年茅台股价涨跌幅作出偏自相关函数图，在5%的显著水平下，取置信区间为

$\frac{2}{\sqrt{n}} = \frac{2}{\sqrt{243}}$

查看Pacf图：

fig = plt.figure(figsize=(12,6))
ax1 = fig.add_subplot(111)
fig = sm.graphics.tsa.plot_pacf(change_data, lags=40, ax=ax1)

很遗憾，从上图中并不太容易确定一个明显的阶数p。不过确实可以看出，随着时间的增加，偏相关函数逐渐趋近于0. 这个结果给了我们两点思路：首先偏相关函数趋近于0，说明运用AR模型模拟该序列是有一定的合理性；其次上图很难确定阶数p或阶数p很大，说明模型可能有MA（滑动平均，将在后面介绍）的部分。它至少给了我们一些感性的认识。

信息准则函数

目前比较常用的信息准则函数包括赤池信息量（AIC），贝叶斯信息量（BIC）和汉南-奎因信息量（HQIC）。他们都基于似然函数L和参数的数量k，其中n表示样本大小。由于Python中有相应的函数，所以公式仅作了解：

$AIC = -2 \ln L + 2 k$

$BIC = -2 \ln L + \ln n \cdot k$

$HQIC = -2 \ln L + \ln(\ln n) \cdot k$

我们的目标是，选择阶数k使得信息量的值最小。以上这么多可供选择的模型，我们通常采用AIC法则。我们知道：增加自由参数的数目提高了拟合的优良性，AIC鼓励数据拟合的优良性但是尽量避免出现过度拟合(Overfitting)的情况。所以优先考虑的模型应是AIC值最小的那一个。

下面我们来测试下3种准则下确定的p，仍然用日收益率序列。为了减少计算量，我们只计算间隔前7个看看效果。

aicList = []
bicList = []
hqicList = []
for i in range(1,8):  #从1阶开始算
    order = (i,0)  # 使用ARMA模型，其中MA阶数为0，只考虑了AR阶数。
    chgdataModel = sm.tsa.ARMA(change_data,order).fit()
    aicList.append(chgdataModel.aic)
    bicList.append(chgdataModel.bic)
    hqicList.append(chgdataModel.hqic)
 
plt.figure(figsize=(10,6))
plt.plot(aicList,'ro--',label='aic value')
plt.plot(bicList,'bo--',label='bic value')
plt.plot(hqicList,'ko--',label='hqic value')
plt.legend(loc=0)

计算信息准则函数的时候，我们发现随着阶数的增加，运算速度越来越慢，所以只取了前7阶，可以看出AIC给出的阶数为0，BIC和HQIC给出的阶数为0. 当然，前面提到只用AR模型预测该序列是不准确的，且只计算前7阶的结果也不足以判断最佳阶数。

当然，利用上面的方法逐个计算是很耗时间的，实际上，有函数可以直接按照准则计算出合适的order，这个是针对ARMA模型的，我们后续再讨论。

AR模型的检验

如果模型是充分的，其残差序列应该是白噪声，根据我们第一章里介绍的混成检验，可以用来检验残差与白噪声的接近程度。先求出残差序列：

delta = results_AR.fittedvalues - change_data[15:] # 取得误差项
plt.figure(figsize=(10,6))
# plt.plot(change_data[17:],label='original value')
# plt.plot(results_AR.fittedvalues,label='fitted value')
plt.plot(delta,'r',label=' residual error')
plt.legend(loc=0)

然后我们检查它是不是接近白噪声序列:

acf,q,p = sm.tsa.acf(delta,nlags=10,qstat=True)  # 计算自相关系数 及p-value
out = np.c_[range(1,11), acf[1:], q, p]
output=pd.DataFrame(out, columns=['lag', "ACF", "Q", "P-value"])
output = output.set_index('lag')
print(output)

lag        ACF         Q   P-value
1.0   0.001303  0.000363  0.984792
2.0  -0.003122  0.002459  0.998771
3.0   0.013582  0.042318  0.997714
4.0  -0.014884  0.090414  0.999008
5.0   0.014418  0.135770  0.999656
6.0   0.007551  0.148269  0.999936
7.0  -0.003774  0.151407  0.999990
8.0  -0.014438  0.197556  0.999996
9.0   0.022536  0.310547  0.999996
10.0 -0.014879  0.360051  0.999999

观察p-value可知，该序列可以认为没有相关性，近似得可以认为残差序列接近白噪声。

拟合优度

调整后的拟合优度检验：

$AdjR^2 = 1 - \frac{\hat{\sigma}^2_{\epsilon}}{\hat{\sigma}^2_{r_t}}$

右边分子代表残差的方差，分母代表 ${r_t}$ 的方差。adjR²值在0到1之间，越大说明拟合效果越好。计算该序列的adjR²值，发现效果并不好，这也符合我们的预期。

下面我们计算之前对日收益率的AR模型的拟合优度：

adjR = 1 - delta.var()/change_data[15:].var()
print(adjR)

输出值为：0.05178424488182409，可以看出，模型的拟合程度并不好，当然，这并不重要，也许是这个序列并不适合用AR模型拟合。

AR模型的预测

我们首先得把原来的样本分为训练集和测试集，再来看预测效果，还是以之前的数据为例：

train = change_data[:-10]
test = change_data[-10:]
output = sm.tsa.AR(train).fit()  
output.predict()
 
predicts = output.predict(216, 225, dynamic=True) #共226个数据
compare = pd.DataFrame()
compare['original'] = change_data[-10:]
compare['predict'] = predicts
print(compare)

   original   predict
0  1.608146  0.853237
1  0.388352 -0.585180
2 -1.726237  0.435650
3  1.406797  0.482424
4 -0.406002  0.967699
5 -3.883607 -0.727744
6 -1.830801  0.412699
7 -0.174712 -0.076781
8 -2.877610  0.223773
9 -1.677702 -0.055838

compare.plot()

该模型的预测结果是很差。那么我们是不是可以通过其他模型获得更好的结果呢？我们将在下一部分继续介绍。

滑动平均（moving-average, MA）模型

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

$r_t = c_0 + a_t - \theta_1a_{t-1} - \cdots - \theta_qa_{t-q}$

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

MA模型的性质

平稳性

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

$E(r_t) = c_0$

$Var(r_t) = (1+ \theta_1^2 + \theta_2^2 + \cdots +\theta_q^2)\sigma_a^2$

自相关函数

对q阶的MA模型，其自相关函数ACF总是q步截尾的。因此MA(q)序列只与其前q个延迟值线性相关，从而它是一个“有限记忆”的模型。这一点可以用来确定模型的阶次，后面会介绍。

可逆性

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

1阶： $|\theta_1|< 1$

2阶： $|\theta_2| < 1, \theta_1 + \theta_2 < 1$

MA的阶次判定

之前我们知道，通过计算偏自相关函数（pACF）可以粗略确定AR模型的阶数。而对于MA模型，需要采用自相关函数（ACF）进行定阶。

fig2 = plt.figure(figsize=(12,6))
ax2 = fig2.add_subplot(111)
fig2 = sm.graphics.tsa.plot_acf(change_data,lags=40,ax=ax2)

可以发现还是不太明显，这里以6阶为列，作出MA模型与原日涨跌幅的图像：

results_MA = sm.tsa.ARMA(change_data, (0,6)).fit() # 取ARMA模型中AR的阶数为0，MA的阶数为6
plt.figure(figsize=(12,6))
plt.plot(change_data, 'b', label = 'changepct')
plt.plot(results_MA.fittedvalues, 'r', label = 'MA model')
plt.legend()

MA模型的预测

由于sm.tsa中没有单独的MA模块，我们利用ARMA模块，只要将其中AR的阶p设为0即可。函数sm.tsa.ARMA中的输入参数中的order(p,q)，代表了AR和MA的阶次。模型阶次增高，计算量急剧增长，因此这里就建立6阶的模型作为示例。我们用最后10个数据作为out-sample的样本，用来对比预测值。

train = change_data[:-10]
test = change_data[-10:]
chgdataModel_MA = sm.tsa.ARMA(train, (0,6)).fit() # 取ARMA模型中AR的阶数为0，MA的阶数为6

我们先来看看拟合效果：

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

可以看出，score为 0.028525493755446107远小于1，拟合效果不好。然后我们用建立的模型进行预测最后10个数据：

predicts_MA = chgdataModel_MA.predict(216, 225, dynamic=True)
compare_MA = pd.DataFrame()
compare_MA['original'] = test
compare_MA['predict_MA'] = predicts_MA
print(compare_MA)

可以看出，建立的模型效果很差，就算只看涨跌方向，正确率也不足。该模型不适用于原数据。没关系，如果真的这么简单能预测股价日涨跌才奇怪了~关于MA的内容只做了简单介绍，下面主要介绍ARMA模型。

自回归-滑动平均（ARMA）模型

在有些应用中，我们需要高阶的AR或MA模型才能充分地描述数据的动态结构，这样问题会变得很繁琐。为了克服这个困难，提出了自回归滑动平均(ARMA)模型。基本思想是把AR和MA模型结合在一起，使所使用的参数个数保持很小。模型的形式为：

$r_t = \phi_0 + \sum_{i=1}^p\phi_ir_{t-i} + a_t + \sum_{i=1}^q\theta_ia_{t-i}$

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

利用向后推移算子（后移算子，即上一时刻），上述模型可写成:

$(1-\phi_1B - \cdots -\phi_pB^p)r_t = \phi_0 + (1-\theta_1B-\cdots - \theta_qB^q)a_t$

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

$E(r_t) = \frac{\phi_0}{1-\phi_1-\cdots - \phi_p}$

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

$1 - \phi_1x - \phi_2x^2 - \cdots -\phi_px^p = 0$

该方程所有解的倒数称为该模型的特征根，如果所有的特征根的模都小于1，则该ARMA模型是平稳的。有一点很关键：ARMA模型的应用对象应该为平稳序列！下面的步骤都是建立在假设原序列平稳的条件下的。

识别ARMA模型阶次

PACF、ACF 判断模型阶次

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

fig = plt.figure(figsize=(12,12))
ax1=fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(change_data,lags=30,ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(change_data,lags=30,ax=ax2)

可以看出，模型的阶次应该为(6,24)。然而，这么高的阶次建模的计算量是巨大的。为什么不再限制滞后阶数小一些？如果lags设置为25，20或者更小时，阶数为(0,0)，显然不是我们想要的结果。

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

信息准则定阶

关于信息准则，已有过一些介绍。下面，我们分别应用以上3种法则，为我们的模型定阶，数据仍然是日涨跌幅序列。为了控制计算量，我们限制AR最大阶不超过6，MA最大阶不超过4。但是这样带来的坏处是可能为局部最优。

import warnings
warnings.filterwarnings("ignore")
sm.tsa.arma_order_select_ic(change_data,max_ar=6,max_ma=4,ic='aic')['aic_min_order']  # AIC
#(5, 4)
 
sm.tsa.arma_order_select_ic(change_data,max_ar=6,max_ma=4,ic='bic')['bic_min_order']  # BIC
#(0, 0)
 
sm.tsa.arma_order_select_ic(change_data,max_ar=6,max_ma=4,ic='hqic')['hqic_min_order'] # HQIC
#(0, 0)

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

ARMA模型的预测和检验

预测的流程和MA流程一致，这里不再重复。预测结果的拟合优度得分仍然很低，但是从预测图上来看，好像比之前强了一点，至少在中间部分效果还不错。

可能的原因有：

信息准则定阶时，为减少计算量，我们设定了AR和MA模型的最大阶数，这样得到的阶数只是局部最优解。
ARMA模型不足以预测股价的变化趋势（这也正常，不然人人都可以赚大钱了）

差分自回归滑动平均（ARIMA）模型

到目前为止，我们研究的序列都集中在平稳序列。即ARMA模型研究的对象为平稳序列。如果序列是非平稳的，就可以考虑使用ARIMA模型。ARIMA比ARMA仅多了个”I”，代表着其比ARMA多一层内涵：也就是差分。一个非平稳序列经过d次差分后，可以转化为平稳时间序列。d具体的取值，我们得分被对差分1次后的序列进行平稳性检验，若果是非平稳的，则继续差分。直到d次后检验为平稳序列。

由于整体的流程与ARMA流程类似，这里就不再详述了。

打赏作者微信支付

支付宝