时间序列预测在金融领域的运用

时间序列简介

 

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

长期趋势(T)
是在较长时期内受某种根本性因素作用而形成的总的变动趋势

季节变动(S)
是在一年内随着季节的变化而发生的有规律的周期性变动。它是诸如气候条件、生产条件、节假日或人们的风俗习惯等各种因素影响的结果。

循环变动(C)
是时间序列呈现出得非固定长度的周期性变动。循环波动的周期可能会持续一段时间,但与趋势不同,它不是朝着单一方向的持续变动,而是涨落相同的交替波动。

不规则变动(I)
是时间序列中除去趋势、季节变动和周期波动之后的随机波动。不规则波动通常总是夹杂在时间序列中,致使时间序列产生一种波浪形或震荡式的变动。只含有随机波动的序列也称为平稳序列。

序列的平稳性

 

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

 

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

严平稳
:如果对所有的时刻t,任意正整数k和任意k个正整数
,
的联合分布与
的联合分布相同,我们称时间序列是严平稳的。也就是,
的联合分布在时间的平移变换下保持不变,这是个很强的条件。

宽平稳
:若时间序列
满足下面两个条件:
,是常数
,只依赖于l则时间序列
为弱平稳的。即该序列的均值,与的协方差不随时间而改变,l为任意整数。

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

 

差分,就是求时间序列
在t时刻的值与t−1时刻的值的差,则我们得到了一个新序列
,为一阶差分,对新序列
再做同样的操作,则为二阶差分。通常非平稳序列可以经过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))

 

 

相关系数和自相关函数

 

相关系数

 

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

 

 

对于
我们叫做内积,例如:

 

 

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

 

 

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

 

 

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

 

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

 

自相关函数 (Autocorrelation Function, ACF)

 

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

 

与的相关系数称为的间隔为的自相关系数,通常记为:

 

 

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

 

 

对一个平稳时间序列的样本
,,则间隔为的样本自相关系数的估计为:

 

 

则函数
称为的 样本自相关函数(ACF)

 

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

 

混成检验

 

原假设:

 

 

备择假设:

 

 

混成检验统计量:

 

 


渐进服从自由度为的
分布

 

决策规则:

 

 

即,
的值大于自由度为的卡方分布
分位点时,我们拒绝。

 

大部分软件会给出
的p-value,则当p-value小于等于显着性水平时拒绝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;则我们拒绝原假设。因此,我们认为该序列是序列相关的。我们再来看看同期的日收益率序列:

 

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。我们选择假设,即收益率序列没有显着的相关性。

 

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

 

白噪声序列

 

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

 

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

 

线性时间序列

 

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

 

 

为的均值,
为白噪声序列,则我们称{} 为线性序列。其中称为在时刻的 新息(innovation)
或 扰动(shock)。

 

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

 

 


为的方差。因为
一定小于正无穷,因此
必须是 收敛序列
,因此满足时,
。即,随着的增大,远处的扰动对的影响会逐渐消失。

 

自回归(AR)模型

 

在上节中,我们计算了部分数据端的ACF,看表可知间隔为1时自相关系数是显着的。这说明在时刻的数据,在预测时刻时可能是有用的!

 

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

 

 

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

 

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

 

 

AR(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模型:

 

 

 

 

 

 

 

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

 

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

 

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

 

 

查看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中有相应的函数,所以公式仅作了解:

 

 

 

 

 

 

我们的目标是,选择阶数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²值在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)模型的形式:

 

 

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

 

MA模型的性质

 

平稳性

 

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

 

 

 

自相关函数

 

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

 

可逆性

 

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

 

1阶:

 

2阶:

 

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模型结合在一起,使所使用的参数个数保持很小。模型的形式为:

 

 

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

 

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

 

 

这时候我们求的期望,得到:

 

 

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

 

 

该方程 所有解的倒数称为该模型的特征根
,如果所有的 特征根的模都小于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流程类似,这里就不再详述了。

发表评论

电子邮件地址不会被公开。 必填项已用*标注