时间序列入门指南

简介

时间序列分析是统计学和计量经济学的一个子领域。时间序列数据\(y_{t}\)按时间\(t\)索引并按顺序排列。这带来了独特的挑战，包括数据内部的自相关性、数据点的不可交换性以及数据和参数的非平稳性。由于数据的顺序特性，时间序列分析具有特定的目标。我们可以将这些目标总结为：对时间序列进行描述（涉及潜在成分或感兴趣的特征），以及预测（旨在产生合理的未来预测）（Harvey, 1990）。

从开始到结束，我们可以将时间序列建模置于一个遵循Box循环(Blei, D.M. 2014)精神的框架中。具体来说，我们：

1. 为时间序列数据构建模型
2. 对模型进行推理
3. 检查模型拟合度，执行评估与批评
4. 修改模型，重复直到满意
5. 使用模型进行回顾和预测

下面我们以摩根大通股票数据为例，概述一个模型构建过程，其中索引为日频数据。考虑以下时间序列数据：

import pandas as pd
import numpy as np
from pandas_datareader.data import DataReader
from datetime import datetime

a = DataReader('JPM',  'yahoo', datetime(2006,6,1), datetime(2016,6,1))
a_returns = pd.DataFrame(np.diff(np.log(a['Adj Close'].values)))
a_returns.index = a.index.values[1:a.index.values.shape[0]]
a_returns.columns = ["JPM Returns"]

a_returns.head()

	JPM收益率
2006-06-02	0.005264
2006-06-05	-0.019360
2006-06-06	-0.013826
2006-06-07	-0.003072
2006-06-08	0.002364

数据的索引对此数据具有意义；我们不能简单地“洗牌”，否则可能会失去诸如季节性、趋势、周期和其他组成部分的有意义的依赖关系。

第一步：数据可视化

由于时间序列具有顺序性，绘制数据图表可以帮助我们了解其特性。我们还可以绘制数据的自相关图（以及数据的变换）来理解序列中是否存在自相关性。最后，在这一阶段，我们可以思考可能解释序列变化的潜在特征。

对于我们的股市数据，我们可以先绘制数据：

plt.figure(figsize=(15, 5))
plt.ylabel("Returns")
plt.plot(a_returns)
plt.show()

可以看出，该时间序列的波动性随时间变化，并在市场动荡时期（如2008年金融危机期间）呈现聚集性特征。通过绘制收益率及其平方值的自相关函数图，我们可以获得更深入的洞察：

import pyflux as pf
import matplotlib.pyplot as plt
pf.acf_plot(a_returns.values.T[0])
pf.acf_plot(np.square(a_returns.values.T[0]))

平方收益显示出强烈的自相关证据。自相关持续存在并在多个滞后期间衰减的事实表明波动率中存在自回归效应。对于收益而言，自相关的证据较弱，尽管第一个滞后项是显著的。

第二步：提出模型

我们推理一个能够解释数据变异的模型，并指定我们对模型参数的任何先验信念。我们观察到了波动聚集的证据。建模这种效应的一种方法是通过GARCH模型来描述波动性（Bollerslev, T. 1986）。

\[y_{t} \sim N\left(\mu,\sigma_{t}\right)\]

\[\sigma_{t}^{2} = \omega + \alpha\epsilon_{t}^{2} + \beta{\sigma_{t-1}^{2}}\]

我们将对这个模型进行贝叶斯推断，因此需要指定一些先验分布。我们将通过对数变换确保\(\omega > 0\)，并对\(\alpha, \beta\)使用截断正态先验：

my_model = pf.GARCH(p=1, q=1, data=a_returns)
print(my_model.latent_variables)

  Index    Latent Variable     Prior           Prior Hyperparameters   V.I. Dist  Transform
  ======== =================== =============== ======================= ========== ==========
  0        Vol Constant        Normal          mu0: 0, sigma0: 3       Normal     exp
  1        q(1)                Normal          mu0: 0, sigma0: 0.5     Normal     logit
  2        p(1)                Normal          mu0: 0, sigma0: 0.5     Normal     logit
  3        Returns Constant    Normal          mu0: 0, sigma0: 3       Normal     None

my_model.adjust_prior(1, pf.TruncatedNormal(0.01, 0.5, lower=0.0, upper=1.0))
my_model.adjust_prior(2, pf.TruncatedNormal(0.97, 0.5, lower=0.0, upper=1.0))

第三步：执行推理

第三步我们需要决定如何进行模型的推断。下面我们使用Metropolis-Hastings方法对我们的GARCH模型进行近似推断。我们还绘制了潜在变量\(\alpha\)和\(\beta\)：

result = my_model.fit('M-H', nsims=20000)

Tuning complete! Now sampling.
Acceptance rate of Metropolis-Hastings is 0.33865

my_model.plot_z([1,2])

第四步：评估模型拟合度

接下来我们评估模型的拟合效果，并确定是否可以进一步改进模型。对于时间序列，可视化拟合的最简单方法是将序列与其预测值进行对比绘图；我们还可以检查样本外表现。如果我们寻求进一步的模型改进，我们将回到第二步并继续。一旦满意，我们就进入第五步。

下面我们绘制GARCH模型的拟合结果，可以观察到它捕捉到了序列中的波动聚集现象：

my_model.plot_fit(figsize=(15,5))

我们还可以绘制后验预测密度的样本：

my_model.plot_sample(nsims=10, figsize=(15,7))

我们可以看到样本（彩色部分）似乎捕捉到了数据中的变化（方形数据点）。

我们还可以对生成序列的特征进行后验预测检查(PPC)，例如峰度：

from scipy.stats import kurtosis
my_model.plot_ppc(T=kurtosis)

看起来我们生成的数据低估了序列的峰度。这并不奇怪，因为我们假设了正态分布的收益率，因此可能需要考虑其他波动率模型。

第五步：分析与预测

当我们对模型感到满意后，就可以用它来分析历史时间序列并进行预测。对于我们的GARCH模型，从前面的拟合图中可以看出，主要捕捉到的波动期是在2007-2008年金融危机期间，以及2011年末欧债危机期间。我们还可以用该模型获得未来预测：

my_model.plot_predict(h=30, figsize=(15,5))

参考文献

Blei, D. M. (2014). 构建、计算、评估、迭代：基于隐变量模型的数据分析方法. 《统计学年评及其应用》, 1, 203–232.

Bollerslev, T. (1986). 广义自回归条件异方差性. 计量经济学杂志. 四月, 31:3, 第307-27页.

哈维 A. C. (1990). 《预测、结构时间序列模型与卡尔曼滤波》. 剑桥大学出版社.