动态线性回归模型

简介

高斯状态空间模型——通常被称为结构时间序列或不可观测成分模型——提供了一种将时间序列分解为多个不同成分的方法。如果误差项联合服从高斯分布,则可以使用卡尔曼滤波器以闭合形式提取这些成分,并通过预测误差分解和最大似然估计来估计参数。

我们可以在状态空间框架中支持动态线性回归

\[y_{t} = \boldsymbol{x}_{t}^{'}\boldsymbol{\beta}_{t} + \epsilon_{t}\]
\[\boldsymbol{\beta}_{t} = \boldsymbol{\beta}_{t-1} + \boldsymbol{\eta}_{t}\]
\[\epsilon_{t} \sim N\left(0,\sigma_{\epsilon}^{2}\right)\]
\[\boldsymbol{\eta}_{t} \sim N\left(\boldsymbol{0},\Sigma_{\eta}\right)\]

示例

在金融投资组合构建中,我们经常关注股票的\(\beta\)系数,该系数可用于构建收益的系统性组成部分。但这个系数可能并非静态值。对于正态分布的收益(!),我们可以使用基于卡尔曼滤波和平滑算法的动态线性回归模型来追踪其演变。首先让我们获取一些超额收益数据。我们将以亚马逊股票(AMZN)为例,并使用标普500指数代表"市场"。

from pandas_datareader import DataReader
from datetime import datetime

a = DataReader('AMZN',  'yahoo', datetime(2012,1,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 = ["Amazon Returns"]

spy = DataReader('SPY',  'yahoo', datetime(2012,1,1), datetime(2016,6,1))
spy_returns = pd.DataFrame(np.diff(np.log(spy['Adj Close'].values)))
spy_returns.index = spy.index.values[1:spy.index.values.shape[0]]
spy_returns.columns = ['S&P500 Returns']

one_mon = DataReader('DGS1MO', 'fred',datetime(2012,1,1), datetime(2016,6,1))
one_day = np.log(1+one_mon)/365

returns = pd.concat([one_day,a_returns,spy_returns],axis=1).dropna()
excess_m = returns["Amazon Returns"].values - returns['DGS1MO'].values
excess_spy = returns["S&P500 Returns"].values - returns['DGS1MO'].values
final_returns = pd.DataFrame(np.transpose([excess_m,excess_spy, returns['DGS1MO'].values]))
final_returns.columns=["Amazon","SP500","Risk-free rate"]
final_returns.index = returns.index

plt.figure(figsize=(15,5))
plt.title("Excess Returns")
x = plt.plot(final_returns);
plt.legend(iter(x), final_returns.columns);
http://www.pyflux.com/notebooks/GaussianStateSpace/output_43_1.png

我们在此定义一个动态线性回归如下:

model = pf.DynReg('Amazon ~ SP500', data=final_returns)

我们也可以使用更高级别的封装器,它允许我们指定分布族,不过如果我们选择一个非高斯分布族,模型将以不同的方式进行估计(不是通过卡尔曼滤波器):

model = pf.DynamicGLM('Amazon ~ SP500', data=USgrowth, family=pf.Normal())

接下来我们估计潜在变量。对于这个例子,我们将使用最大似然点质量估计 \(z^{MLE}\)

x = model.fit()
x.summary()

Dynamic Linear Regression
====================================== =================================================
Dependent Variable: Amazon             Method: MLE
Start Date: 2012-01-04 00:00:00        Log Likelihood: 2871.5419
End Date: 2016-06-01 00:00:00          AIC: -5737.0838
Number of observations: 1101           BIC: -5722.0719
========================================================================================
Latent Variable         Estimate   Std Error  z        P>|z|    95% C.I.
======================= ========== ========== ======== ======== ========================
Sigma^2 irregular       0.0003
Sigma^2 1               0.0
Sigma^2 SP500           0.0024
========================================================================================

我们可以使用plot_fit()来绘制样本内拟合图:

model.plot_fit(figsize=(15,15))
http://www.pyflux.com/notebooks/GaussianStateSpace/output_47_0.png

第三张图展示了\(\beta_{AMZN}\)。在初始预热期过后,2013年\(\beta\)值在1上方徘徊,尽管2014年它与市场表现高度相关。最近它又稳定下来,在1上方浮动。第四张图显示了剩余的收益残差成分(不包括\(\alpha\))。

类描述

class DynLin(formula, data)

动态线性回归模型

参数 类型 描述
formula string 使用Patsy符号指定回归模型
data pd.DataFrame 包含单变量时间序列

属性

latent_variables

一个包含模型潜在变量信息的pf.LatentVariables()对象,包括先验设置、任何拟合值、初始值和其他潜在变量信息。当模型被拟合时,这里就是潜在变量被更新/存储的地方。有关此对象内属性的信息以及访问潜在变量信息的方法,请参阅潜在变量文档。

方法

adjust_prior(index, prior)

调整模型潜在变量的先验分布。潜在变量及其索引可以通过打印附加到模型实例的latent_variables属性来查看。

参数 类型 描述
index int 要更改的潜变量索引
prior pf.Family instance Prior distribution, e.g. pf.Normal()

返回: void - 修改模型的 latent_variables 属性

fit(method, **kwargs)

估计模型的潜在变量。用户选择一个推断选项,该方法会返回一个结果对象,同时更新模型的latent_variables属性。

参数 类型 描述
method str 推断选项:例如 'M-H' 或 'MLE'

请参阅文档中的贝叶斯推断和经典推断部分,了解完整的推断选项列表。可以输入与所选特定推断模式相关的可选参数。

返回: 包含估计潜在变量信息的pf.Results实例

plot_fit(**kwargs)

绘制模型对数据的拟合情况。可选参数包括figsize,即绘图图形的尺寸。

返回 : void - 显示一个matplotlib绘图

plot_ppc(T, nsims)

绘制后验预测检查的直方图,使用用户选择的差异度量。此方法仅在使用贝叶斯推断进行拟合时有效。

参数 类型 描述
T function Discrepancy, e.g. np.mean or np.max
nsims int PPC需要进行多少次模拟

返回值: void - 显示一个matplotlib绘图

plot_predict(h, oos_data, past_values, intervals, **kwargs)

绘制模型的预测结果,并附带置信区间。

参数 类型 描述
h int 预测向前多少步
oos_data pd.DataFrame 包含h步外生变量的数据框
past_values int 要绘制的历史数据点数量
intervals boolean 是否绘制区间

需要明确的是,oos_data参数应该是一个与初始化模型实例时使用的初始数据框格式相同的DataFrame。原因在于,要预测未来值,您需要指定关于未来外生变量的假设。例如,如果您预测h步向前,该方法将从oos_data中获取前h行数据,并提取您在patsy公式中指定的外生变量值。

可选参数包括figsize - 图表绘制的尺寸。请注意 如果您使用最大似然估计或变分推断,显示的区间将不会 反映潜在变量的不确定性。只有Metropolis-Hastings方法能提供完全贝叶斯 预测区间。由于平均场推断的局限性(无法考虑后验相关性), 变分推断的贝叶斯区间不予显示。

返回 : void - 显示一个matplotlib绘图

plot_predict_is(h, fit_once, fit_method, **kwargs)

绘制模型在样本内的滚动预测。这意味着用户假装数据的最后一部分是样本外的,并在每个时间段后进行预测并评估其表现。用户可以选择是在开始时一次性拟合参数,还是在每个时间步都进行拟合。

参数 类型 描述
h int 使用多少个先前的时间步
fit_once boolean 是否只拟合一次,还是每个时间步都拟合
fit_method str 选择哪种推断方法,例如'MLE'

可选参数包括figsize - 要绘制的图形尺寸。h是一个整数,表示要模拟性能的前几步。

返回 : void - 显示一个matplotlib绘图

plot_sample(nsims, plot_data=True)

绘制模型后验预测密度的样本。此方法仅适用于通过贝叶斯推断拟合的模型。

参数 类型 描述
nsims int 需要抽取的样本数量
plot_data boolean 是否同时绘制真实数据

返回 : void - 显示一个matplotlib绘图

plot_z(indices, figsize)

返回潜在变量及其相关不确定性的绘图。

参数 类型 描述
indices int or list 要绘制的潜变量索引
figsize tuple matplotlib图形的大小

返回 : void - 显示一个matplotlib绘图

ppc(T, nsims)

返回后验预测检验的p值。此方法仅在您使用贝叶斯推断进行拟合时才有效。

参数 类型 描述
T function Discrepancy, e.g. np.mean or np.max
nsims int PPC需要进行多少次模拟

返回值: int - 差异检验的p值

predict(h, oos_data, intervals=False)

返回模型预测结果的DataFrame。

参数 类型 描述
h int 预测向前多少步
oos_data pd.DataFrame 包含h步外生变量的数据框
intervals boolean 是否返回预测区间

需要明确的是,oos_data参数应该是一个与初始化模型实例时使用的原始数据框格式相同的DataFrame。原因在于,要预测未来值,您需要指定关于未来外生变量的假设。例如,如果要预测未来h步,该方法将从oos_data中获取前5行数据,并采用您在patsy公式中指定为外生变量的值。

请注意,如果您使用最大似然估计或变分推断,显示的区间将不会反映潜在变量的不确定性。只有Metropolis-Hastings方法能提供完全贝叶斯预测区间。由于平均场推断的局限性(未考虑后验相关性),变分推断的贝叶斯区间未予显示。

返回 : pd.DataFrame - 模型预测结果

predict_is(h, fit_once, fit_method)

返回模型样本内滚动预测的DataFrame。

参数 类型 描述
h int 使用多少个先前的时间步
fit_once boolean 是否只拟合一次,还是每个时间步都拟合
fit_method str 选择哪种推断方法,例如'MLE'

返回 : pd.DataFrame - 模型预测结果

sample(nsims)

返回从后验预测密度中抽取的数据的np.ndarray数组。此方法仅在您使用贝叶斯推断拟合模型时才有效。

参数 类型 描述
nsims int 进行多少次后验抽样

返回值 : np.ndarray - 来自后验预测密度的样本。

simulation_smoother(beta)

返回从Durbin和Koopman(2002)模拟平滑器中抽取的数据的np.ndarray数组。

参数 类型 描述
beta np.array 潜在变量的np.array数组

如果已经拟合了模型,建议直接使用model.latent_variables.get_z_values()作为beta输入。

返回 : np.ndarray - 来自模拟平滑器的样本

参考文献

Durbin, J. 和 Koopman, S. J. (2002). 一种简单高效的状态空间时间序列分析模拟平滑器。Biometrika, 89(3):603–615.

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