在金融世界中，绝大多数的金融数据，从股票的价格到GDP数值，都是一个随时间变化的变量，其在不同时点的不同值构成的序列即所谓时间序列。同时，金融时间序列又与其他普通的时间序列有较大区别：金融理论以及金融时间序列都包含不确定因素，例如，资产波动率有各种不同的定义，对于一个股票收益率序列，波动率是不能直接观察到的。
正因为带有不确定性，统计的理论和方法在金融时间序列分析中起到重要作用。因此，时间序列分析是研究金融数据的重要工具，本系列将由浅入深地介绍金融时间序列分析的常用方法以及Python工具，并展示效果。
本篇的主要内容有：
简而言之：
对某一个或者一组变量 x(t) 进行观察测量，将在一系列时刻 $t_1,t_2,⋯,t_n$ 所得到的离散数字组成的序列集合，称之为时间序列。
例如: 某股票A从2015年6月1日到2016年6月1日之间各个交易日的收盘价，可以构成一个时间序列；某地每天的最高气温可以构成一个时间序列。
时间序列可能的一些特征:
趋势：是时间序列在长时期内呈现出来的持续向上或持续向下的变动。
季节变动：是时间序列在一年内重复出现的周期性波动。它是诸如气候条件、生产条件、节假日或人们的风俗习惯等各种因素影响的结果。
循环波动：是时间序列呈现出得非固定长度的周期性变动。循环波动的周期可能会持续一段时间，但与趋势不同，它不是朝着单一方向的持续变动，而是涨落相同的交替波动。
不规则波动：是时间序列中除去趋势、季节变动和周期波动之后的随机波动。不规则波动通常总是夹杂在时间序列中，致使时间序列产生一种波浪形或震荡式的变动。只含有随机波动的序列也称为平稳序列。
# 获取平安银行收盘价数据
raw_data = D.history_data(['000001.SZA'],'2015-01-01','2018-01-01',['close'])[['date','close']].set_index('date')
raw_data['close_diff_1'] = raw_data['close'].diff(1) # 1阶差分处理
raw_data['close_diff_2'] = raw_data['close_diff_1'].diff(1) # 2阶差分处理
T.plot(raw_data[['close', 'close_diff_1', 'close_diff_2']],
# high、low显示在第一栏，高度40%，open、close显示在第二栏，其他的在最后一栏
panes=[['close', '40%'], ['close_diff_1', '40%']],
# height=500，设置高度为500
options={'chart':{'height': 500},
'title': {'text': '价格走势及一阶差分、二阶差分'}
})
上图中第一张图为中国银行部分年份的收盘指数，是一个非平稳时间序列；而下面两张为平稳时间序列(当然这里没有检验,只是为了让大家看出差异,关于检验序列的平稳性后续会讨论)
细心的读者已经发现，下面两张图，实际上是对第一个序列做了差分处理，方差和均值基本平稳，成为了平稳时间序列，后面我们会谈到这种处理。
下面可以给出平稳性的定义了：
严平稳：
如果对所有的时刻 t，任意正整数 k 和任意 k 个正整数 $$ \large (t_1,t_2......t_k), \\ \large (r_{t_1},r_{t_2},......r_{t_k})$$ 的联合分布与 $$ \large (r_{t_1 + t},r_{t_2+t},......r_{t_k + t}) $$ 的联合分布相同，我们称时间序列 {$r_t$} 是严平稳的。
也就是，$$ \large (r_{t_1},r_{t_2},......r_{t_k})$$ 的联合分布在时间的平移变换下保持不变，这是个很强的条件，难以用经验方法验证。因此经常假定的是平稳性的一个较弱的方式
弱平稳：
若时间序列 {$r_t$} 满足下面两个条件:$$ \large E(r_t) = \mu, \mu 是常数\\ \large Cov(r_t,r_{t-l}) = \gamma_l, \gamma_l 只依赖于l $$ 则时间序列 {$r_t$} 为弱平稳的。即该序列的均值，$r_t$与$r_{t-l}$ 的协方差不随时间而改变，l为任意整数。
在金融数据中，通常我们所说的平稳序列，是指弱平稳的序列。
可能会有读者疑问为何要“强行”引入弱稳定性的概念，其实时间序列的弱平稳和协方差的绝对可加性，都是为了保证大数定理在时间序列上成立，这是所有统计分析的基础。
差分
回头我们再谈之前说的差分操作：
差分（这里为前向），就是求时间序列{$r_t$}在 t时刻的值 $r_t$ 与 t-1 时刻的值$r_{t-1}$ 的差不妨记做 dt,则我们得到了一个新序列{dt}，为一阶差分，对新序列{dt}再做同样的操作，则为二阶差分。
通常非平稳序列可以经过d次差分，处理成弱平稳或者近似弱平稳时间序列。回头看上图，我们发现二阶差分得到的序列比一阶差分效果更好。
对于两个向量，我们希望定义它们是不是相关。一个很自然的想法，用向量与向量的夹角来作为距离的定义，夹角小，就距离小，夹角大，就距离大。
早在中学数学中，我们就经常使用余弦公式来计算角度：$$ \large cos<\vec a , \vec b> = \frac {\vec a \cdot \vec b}{\lvert \vec a \rvert \lvert \vec b \rvert} $$
而对于，$$\large \vec a \cdot \vec b$$ 我们叫做内积，例如$$ \large (x_1,y_1) \cdot (x_2,y_2) = x_1x_2 + y_1y_2 $$
我们再来看相关系数的定义公式,X和Y的相关系数为： $$ \large \rho_{xy} = \frac{Cov(X,Y)}{\sqrt{Var(X)Var(Y)}}$$
而我们的根据样本的估计计算公式为; $$ \large \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计算。
对于减去均值操作，并不影响角度计算，是一种“平移”效果，如下图所示：
# 导入相关模块
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
a = pd.Series([9,8,7,5,4,2])
b = a - a.mean() # 去均值
plt.figure(figsize=(10,4))
a.plot(label='a')
b.plot(label='mean removed a')
plt.legend()
相关系数度量了两个向量的线性相关性，而在平稳时间序列 {$r_t$} 中，我们有时候很想知道，$r_t$ 与它的过去值 $r_{t-i}$ 的线性相关性。 这时候我们把相关系数的概念推广到自相关系数。
$r_t$ 与 $r_{t-l}$ 的相关系数称为 $r_t$ 的间隔为 l 的自相关系数，通常记为 $ρ_l$ 。具体的： $$\large \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)} $$ 这里用到了弱平稳序列的性质: $$ \large Var(r_t)=Var(r_{t-l})$$
对一个平稳时间序列的样本 {$r_t$}， 1 <= t <= T,则间隔为 l 的样本自相关系数的估计为: $$ \large \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 $$
则函数 $$ \large \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 $$
混成检验统计量：$$\large Q(m) = T(T+2) \sum_{l=1}^{m} \frac{\hat{\rho_l}^{2}}{T-l}$$ Q(m)渐进服从自由度为m的 $$ \chi^2 $$分布
决策规则： $$ \large Q(m) > \chi_\alpha^{2} ,拒绝H0$$ 即,Q(m)的值大于自由度为m的卡方分布100(1-α)分位点时，我们拒绝H0。
大部分软件会给出Q(m)的 p-value，则当 p-value 小于等于显著性水平 α时拒绝H0。
下面给出示例：
from scipy import stats
import statsmodels.api as sm # 统计相关的库
data = raw_data['close'] # 平安银行收盘价
m = 10 # 我们检验10个自相关系数
acf,q,p = sm.tsa.acf(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')
output
我们取显著性水平为0.05，可以看出，所有的p-value都小于0.05；则我们拒绝原假设H0。
因此，我们认为该序列，即平安银行收盘价序列，是序列相关的
我们再来看看同期平安银行的日收益率序列：
data2 =raw_data['close'].pct_change().dropna() # 平安银行收益率
m = 10 # 我们检验10个自相关系数
acf,q,p = sm.tsa.acf(data2,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')
output
可以看出，p-value均大于显著性水平0.05。我们选择假设H0，即，平安银行日收益率序列没有显著的相关性
随机变量$X(t)$（t=1，2，3……），如果是由一个不相关的随机变量的序列构成的，即对于所有S不等于T，随机变量$X_t和X_s$的协方差为零，则称其为纯随机过程。
对于一个纯随机过程来说，若其期望和方差均为常数，则称之为白噪声过程。白噪声过程的样本实称成为白噪声序列，简称白噪声。之所以称为白噪声，是因为他和白光的特性类似，白光的光谱在各个频率上有相同的强度，白噪声的谱密度在各个频率上的值相同。特别地，若$X_t$还服从均值为0、方差为$\sigma^2$的正态分布，则称这个序列为高斯白噪声。
时间序列{$r_t$}，如果能写成： $$ \large r_t = \mu + \sum_{i=0}^{\infty}\psi_ia_{t-i} \\ \large \mu为r_t 的均值， \psi_0=1,\{a_t\}为白噪声序列$$ 则我们称{$r_t$} 为线性序列。其中$a_t$称为在t时刻的新息(innovation)或扰动(shock)
很多时间序列具有线性性，即是线性时间序列，相应的有很多线性时间序列模型，例如接下来要介绍的AR、MA、ARMA，都是线性模型，但并不是所有的金融时间序列都是线性的
对于弱平稳序列，我们利用白噪声的性质很容易得到$r_t$的均值和方差： $$ \large E(r_t) = \mu , Var(r_t) = \sigma_a^2 \sum_{i=0}^{\infty} \psi_i^{2} \\ \large \sigma_a^2为a_t的方差。$$
因为 $Var(r_t)$一定小于正无穷，因此 $$ \{\psi_i^2\}$$ 必须是收敛序列，因此满足 $$ i \to \infty 时， \psi_i^2 \to 0$$
即，随着i的增大，远处的扰动 $a_{t-i}$ 对 $r_t$ 的影响会逐渐消失。
到目前为止介绍了一些基本知识和概念，如平稳性、相关性、白噪声、线性序列，介绍的过程中并没有太深入，目前来说“够用”了，一些细节会在后面章节补充。下面开始介绍一些线性模型。
在第一部分中，我们计算了平安银行部分数据段的ACF，看表可知间隔为1时自相关系数是显著的。这说明在t-1时刻的数据$r_{t-1}$ ，在预测t时刻时的 $r_t$ 时可能是有用的
根据这点我们可以建立下面的模型：$$ \large r_t = \phi_0 + \phi_1r_{t-1} + a_t$$
其中{$a_t$}是白噪声序列，这个模型与简单线性回归模型有相同的形式，这个模型也叫做一阶自回归(AR)模型，简称AR(1)模型
从AR(1)很容易推广到AR(p)模型：$$ \large r_t = \phi_0 + \phi_1r_{t-1} + \cdots + \phi_pr_{t-p}+ a_t \qquad (2.0)$$
我们先假定序列是弱平稳的，则有; $$\large E(r_t) = \mu \qquad \large Var(r_t) = \gamma_0 \qquad \large Cov(r_t,r_{t-j})=\gamma_j , 其中\mu,\gamma_0 是常数 $$
因为{$a_t$}是白噪声序列，因此有：$$ \large E(a_t) = 0, Var(a_t) = \sigma_a^2$$
所以有：$$ \large 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})=...= u$，从而： $$ \large \mu = \phi_0 + \phi_1\mu + \phi_2\mu+ \cdots +\phi_p\mu\\ \large E(r_t)=\mu=\frac{\phi_0}{1-\phi_1 - \phi_2 - \cdots -\phi_p} \qquad (2.1)$$
对式2.1，假定分母不为0， 我们将下面的方程称为特征方程: $$ \large 1 - \phi_1x - \phi_2x^2 - \cdots -\phi_px^p = 0$$
下面我们就用该方法，检验平安银行日收益率序列的平稳性
temp = np.array(data2) # 载入收益率序列
model = sm.tsa.AR(temp) # 建立AR模型
results_AR = model.fit() # 训练模型一
T.plot(pd.DataFrame({'收益率序列':temp[20:],'AR模型拟合值':results_AR.fittedvalues}))
我们可以看看模型有多少阶
print(len(results_AR.roots))
可以看出，自动生成的AR模型是20阶的。关于阶次的讨论在下节内容，我们画出模型的特征根，来检验平稳性
pi,sin,cos = np.pi,np.sin,np.cos
r1 = 1
theta = np.linspace(0,2*pi,360)
x1 = r1*cos(theta)
y1 = r1*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()
可以看出，所有特征根都在单位圆内，则序列为平稳的！
一般有两种方法来决定p：
第一种：利用偏相关函数(Partial Auto Correlation Function,PACF)
第二种：利用信息准则函数
对于偏相关函数的介绍，这里不详细展开，重点介绍一个性质：
AR(p)序列的样本偏相关函数是 p 步截尾的。
所谓截尾，就是快速收敛应该是快速的降到几乎为0或者在置信区间以内。
具体我们看下面的例子，还是以之前平安银行日收益率序列
fig = plt.figure(figsize=(20,5))
ax1=fig.add_subplot(111)
fig = sm.graphics.tsa.plot_pacf(temp,lags=100,ax=ax1)
我们看出，按照截尾来看，和之前调用的自动生成AR模型，阶数为20比较一致。 当然，我们很少会用这么高的阶次。。
现在有以上这么多可供选择的模型，我们通常采用AIC法则。我们知道：增加自由参数的数目提高了拟合的优良性，AIC鼓励数据拟合的优良性但是尽量避免出现过度拟合(Overfitting)的情况。所以优先考虑的模型应是AIC值最小的那一个。赤池信息准则的方法是寻找可以最好地解释数据但包含最少自由参数的模型。不仅仅包括AIC准则，目前选择模型常用如下准则：
AIC = -2 ln(L) + 2 k 中文名字：赤池信息量 akaike information criterion
BIC = -2 ln(L) + ln(n)*k 中文名字：贝叶斯信息量 bayesian information criterion
HQ = -2 ln(L) + ln(ln(n))*k hannan-quinn criterion
下面我们来测试下3种准则下确定的p，仍然用平安银行日收益率序列。为了减少计算量，我们只计算间隔前10个看看效果。
aicList = []
bicList = []
hqicList = []
for i in range(1,11): #从1阶开始算
order = (i,0) # 这里使用了ARMA模型，order 代表了模型的(p,q)值，我们令q始终为0，就只考虑了AR情况，此时arma(i,0) = ar(i)。
tempModel = sm.tsa.ARMA(temp,order).fit()
aicList.append(tempModel.aic)
bicList.append(tempModel.bic)
hqicList.append(tempModel.hqic)
T.plot(pd.DataFrame({'aicList':aicList,'bicList':bicList,'hqicList':hqicList}),title='阶数p的选择依据')
可以看出，3个准则在第5点均取到最小值，也就是说，p的最佳取值应该在5，我们只计算了前10个，结果未必正确。
whatever，我们的目的是了解方法。
当然，利用上面的方法逐个计算是很耗时间的，实际上，有函数可以直接按照准则计算出合适的 order，这个是针对 ARMA模型的，我们后续再讨论。
根据式2.0，如果模型是充分的，其残差序列应该是白噪声，根据我们第一章里介绍的混成检验，可以用来检验残差与白噪声的接近程度。
我们先求出残差序列：
delta[:10].tolist()
delta = results_AR.fittedvalues - temp[20:] # 残差
T.plot(pd.Series(delta.tolist()), title='residual error')
然后我们检查它是不是接近白噪声序列
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', "AC", "Q", "P-value"])
output = output.set_index('lag')
output
观察p-value可知，该序列可以认为没有相关性，近似得可以认为残差序列接近白噪声。
我们使用下面的统计量来衡量拟合优度：$$ \large R^2 = 1 - \frac{残差的平方和}{总的平方和}$$
但是，对于一个给定的数据集，R2是用参数个数的非降函数，为了克服该缺点，推荐使用调整后的R2: $$ \large Adj R^2 = 1- \frac{残差的平方}{r_t的方差} $$
它的值在0-1之间，越接近1，拟合效果越好。
下面我们计算之前对平安银行日收益率的AR模型的拟合优度。
score = 1 - delta.var()/temp[20:].var()
print('拟合优度为：',score)
可以看出，模型的拟合程度并不好，当然，这并不重要，也许是这个序列并不适合用AR模型拟合。
我们首先得把原来的样本分为训练集和测试集，再来看预测效果，还是以之前的数据为例：
train = temp[:-10]
test = temp[-10:]
output = sm.tsa.AR(train).fit()
predicts = output.predict(len(train), len(train)+9, dynamic=True)
print(len(predicts))
comp = pd.DataFrame()
comp['original'] = temp[-10:] # 测试集
comp['predict'] = predicts
comp
这里我们直接给出MA(q)模型的形式：
$$ \large r_t = c_0 + a_t - \theta_1a_{t-1} - \cdots - \theta_qa_{t-q}$$
$c_0$为一个常数项。这里的$a_t$,是AR模型t时刻的扰动或者说信息，则可以发现，该模型，使用了过去q个时期的随机干扰或预测误差来线性表达当前的预测值。
对q阶的MA模型，其自相关函数ACF总是q步截尾的。 因此MA(q)序列只与其前q个延迟值线性相关，从而它是一个“有限记忆”的模型。
这一点可以用来确定模型的阶次，后面会介绍。
当满足可逆条件的时候，MA(q)模型可以改写为AR(p)模型。这里不进行推导，给出1阶和2阶MA的可逆性条件。
1阶：$$ \large |\theta_1|< 1$$
2阶： $$ \large |\theta_2| < 1, \theta_1 + \theta_2 < 1$$
我们通常利用上面介绍的第二条性质：MA(q)模型的ACF函数q步截尾来判断模型阶次。示例如下：
使用平安银行的日涨跌数据（2015年1月至2018年1月）来进行分析，先取数据：
from scipy import stats
import statsmodels.api as sm # 统计相关的库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
raw_data = D.history_data(['000001.SZA'],'2015-01-01','2018-01-01',['close'])[['date','close']].set_index('date')
raw_data['return'] = raw_data.pct_change()
raw_data = raw_data.dropna()
T.plot(raw_data['return'],title='平安银行涨跌幅')
可以看出，序列看上去是弱平稳的。下面我们画出序列的ACF：
data = np.array(raw_data['return']) # 平安银行涨跌幅
fig = plt.figure(figsize=(20,5))
ax1=fig.add_subplot(111)
fig = sm.graphics.tsa.plot_acf(data,ax=ax1)
我们发现ACF函数在20处截尾，之后的acf函数均在置信区间内，我们判定该序列MA模型阶次为20阶。
由于sm.tsa中没有单独的MA模块，我们利用ARMA模块，只要将其中AR的阶p设为0即可。
函数sm.tsa.ARMA中的输入参数中的order(p,q)，代表了AR和MA的阶次。 模型阶次增高，计算量急剧增长
我们用最后10个数据作为out-sample的样本，用来对比预测值。
order = (0,20)
train = data[:-10] # 训练集
test = data[-10:] # 测试集
tempModel = sm.tsa.ARMA(train,order).fit() # 拟合模型
我们先来看看拟合效果,计算$$ \large Adj R^2 = 1- \frac{残差的平方}{r_t的方差} $$
delta = tempModel.fittedvalues - train
score = 1 - delta.var()/train.var()
print('拟合优度：',score)
可以看出，score远小于1，拟合效果不好。
然后我们用建立的模型进行预测最后10个数据
predicts = tempModel.predict(len(train), len(train)+9, dynamic=True)
print(len(predicts))
comp = pd.DataFrame()
T.plot(pd.DataFrame({'original':test,'predict':predicts}),title='测试集上的预测效果')
可以看出，建立的模型效果很差，预测值明显小了1到2个数量级！就算只看涨跌方向，正确率也不足50%。该模型不适用于原数据。
没关系，如果真的这么简单能股票的涨跌才奇怪了~
1.《金融时间序列分析》 第2版 Ruey S.Tsay著 王辉、潘家柱 译
2.Time Series analysis tsa http://nipy.bic.berkeley.edu/nightly/statsmodels/doc/html/tsa.html