采用时间序列预测股价变化

2,540 阅读20分钟
原文链接: www.biaodianfu.com

时间序列简介

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

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

序列的平稳性

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

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

  • 严平稳:如果对所有的时刻t,任意正整数k和任意k个正整数(t_1,t_2,...,t_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_tr_{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,相关性越高。只不过这里计算的时候对向量做了去均值处理,即中心化操作。而不是直接用向量X, Y计算。

自相关函数 (Autocorrelation Function, ACF)

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

r_tr_{t-l}的相关系数称为r_t的间隔为l的自相关系数,通常记为\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,则间隔为l的样本自相关系数的估计为:

    \[\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)渐进服从自由度为m\chi^2分布

决策规则:

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

即,Q(m)的值大于自由度为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,随机变量X_tX_s的协方差为零,则称其为纯随机过程

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

线性时间序列

时间序列{r_t},如果能写成:

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

\mur_t的均值, psi_0=1,\{a_t\}为白噪声序列,则我们称{r_t} 为线性序列。其中a_t称为在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^2a_t的方差。因为Var(r_t)一定小于正无穷,因此\{\psi_i^2\}必须是收敛序列,因此满足 i \to \infty  时, \psi_i^2 \to 0。即,随着i的增大,远处的扰动a_{t-i}r_t的影响会逐渐消失。

自回归(AR)模型

在上节中,我们计算了部分数据端的ACF,看表可知间隔为1时自相关系数是显著的。这说明在t-1时刻的数据r_{t-1},在预测t时刻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模型t时刻的扰动或者说新息,则可以发现,该模型,使用了过去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}为白噪声序列,pq都是非负整数。AR和MA模型都是ARMA(p,q)的特殊形式

利用向后推移算子B(后移算子 B,即上一时刻),上述模型可写成:

    \[(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流程类似,这里就不再详述了。

打赏作者 微信支付标点符 wechat qrcode 支付宝标点符 alipay qrcode