时间序列简介
在数学上,随机过程被定义为一族时间随机变量,即{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流程类似,这里就不再详述了。
打赏作者 微信支付 支付宝