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

时间序列简介

在数学上，随机过程被定义为一族时间随机变量，即{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流程类似，这里就不再详述了。