Python金融分析（十）：贝叶斯线性回归

发表于 2019-05-08 分类于他山攻错阅读次数： Disqus：

引言

之前关于贝叶斯学派的东西接触得比较少，主要是原来上学的时候觉得先验概率这个东西不靠谱，不过既然碰到了还是了解一下。

贝叶斯统计

了解统计学发展史的可能知道，统计学里有两大著名流派——频率学派和贝叶斯学派。平时我们接触到的可能是频率学派的内容更多一些，但是贝叶斯学派近几十年来也越发得到重视。贝叶斯方法与频率学派（古典学派、经典方法）有着有着基本的区别。

在经典方法中，参数 $\theta$ 被认为是一个未知、但固定的量。从以 $\theta$ 为参数的总体中抽取一组随机样本，基于样本的观测值来获得关于 $\theta$ 的知识，也就是通过样本推断总体。在贝叶斯方法中，参数 $\theta$ 则是一个随机变量，其变化可以用一个概率分布，即先验分布（prior distribution）来描述，实际上这个分布反映了试验者的主观认知。然后从以 $\theta$ 为参数的总体中抽取样本，通过贝叶斯法则利用这些样本去修正试验者的主观认知（先验分布），得到后验分布。（来源：《统计推断》）

贝叶斯定理

贝叶斯定理描述了后验概率与先验概率之间的关系，其基本的公式为：

P(A|B)=\frac{P(B|A)P(A)}{P(B)} \tag{1}

如果把A看成我们对某个事物客观规律（概率）的认识，把B看作在某种程度上与A关联的实际事件，那么：

$P(A)$ 就是A的先验概率，也就是我们对客观规律的初始认知；
$P(B)$ 就是B的先验概率，例如是现实中某类事件发生的历史统计，也叫做标准化常量（normalized constant）；
$P(B|A)$ 就是已知A发生B的条件概率，即B的后验概率，也就是似然性；
$P(A|B)$ 就是A的后验概率，也就是我们通过事件B的发生情况及A、B间的似然性，对初始认知的修正。

我们知道 $\theta$ 的函数:

L(\theta|x)=f(x|\theta)

为 $\theta$ 的似然函数，通过令 $L(\theta|x)$ 最大得到 $\theta$ 估计值的方法，就是极大似然估计。所以上面的似然性，其实就反映了A和B的关联程度，但是在这里还有一个被称为标准化常量的 $P(B)$ ，如果我们 $P(B|A)/P(B)$ 就可以得到标准似然度，相当于将关联程度进行了标准化。为什么这么说呢？我们令 $(1)$ 两边除以 $P(A)$ ：

\frac{P(A|B)}{P(A)}=\frac{P(B|A)}{P(B)}

也就是

\frac{P(A,B)}{P(A)P(B)}=\frac{P(B,A)}{P(A)P(B)} \tag{2}

(2)两边均为标准似然度，其含义也就一目了然了。

将贝叶斯定理用在连续概率密度上，就可以得到：

\begin{aligned} f(x|y)&=\frac{f(x,y}{f(y)}=\frac{f(y|x)f(x)}{f(y)}\\ &=\frac{f(y|x)f(x)}{\int_{ - \infty }^{ + \infty }{f(y|x)f(x)dx} } \end{aligned}

贝叶斯线性回归

说实话，我对贝叶斯回归也不是非常熟悉，贝叶斯线性回归是用来估计线性回归参数的一种方法。

我们经常用到的最小二乘法和极大似然估计，实际上都是基于频率学派的思想，即模型的参数是固定的，因变量是由输入、模型（参数）和误差项共同决定的，当然两者的区别是优化目标不同（残差的表示不同）。

而贝叶斯线性回归，则是基于贝叶斯学派的思想，将模型的参数看成是服从某种概率分布的随机变量，在初始分布即先验分布的基础上，利用收集的样本，通过贝叶斯法则不断迭代更新后验分布的参数（这里面用到了共轭先验），得到最终的后验分布作为参数服从的概率分布。在这个地方，其实跟极大似然估计还是有点关系的：

极大似然估计的时候，是参数的似然函数
贝叶斯线性回归的时候，是样本的似然函数

具体的公式有点复杂，有兴趣的可以查阅Bishop的《Pattern Recognition and Machine Learning》。

共轭先验

如前所述，这里，我们首先要假设参数服从某种分布，那么如何确定先验概率呢？理论上来讲，如果我们没有任何先验知识，这个先验分布实际上是可以任意选择的。但是，为了更方便地简化计算，最理想的情况就是让先验分布和似然分布有相同的形式，也就是两者为共轭分布，此时先验分布叫做似然分布的共轭先验。共轭先验的好处主要在于代数上的方便性，可以直接给出后验分布的封闭形式，否则的话只能数值计算。共轭先验也有助于获得关于似然函数如何更新先验分布的直观印象。

所有的指数分布族都具有共轭先验，而其中的高斯分布族在高斯似然函数下与其自身共轭 (自共轭)，所以我们假设参数的先验分布都服从高斯分布。

MCMC过程

有了先验分布以后，接下来的问题就是：如何估计后验分布呢？与频率学派的目标不同，贝叶斯线性回归要估计的不是某一个数，而是一个概率分布，所以不能够用一般的最优化方法来解决。解决这个问题当然有相应的方法，我们只介绍其中一种用来估计概率分布的方法，那就是马尔科夫链蒙特卡洛（MCMC）方法。

我们常常遇到这样的问题：模型构建好之后，有一个概率分布P(X)（称为目标分布），但是不能显式的给出其表达，只能生成一系列符合这个分布的样本。这种问题称为“采样”。

对于一个给定的概率分布P(X)，若是要得到其样本，通过构造一条马尔可夫链，使得该马尔可夫链的平稳分布为P(X)，这样，无论其初始状态为何值，假设记为 $x_0$ ，那么随着马尔科夫过程的转移，得到了一系列的状态值： $x_0,x_1,\ldots,x_n,x_{n+1},\ldots$ ，如果这个马尔可夫过程在第n步时已经收敛，那么分布P(X)的样本即为 $x_n,x_{n+1},\ldots$ 。

贝叶斯方法中，关注的是后验概率 $P(\theta|X)$ 。在给定观测X的情况下，需要估计系统参数 $\theta$ 。如果后验概率没有明确表达，或者由于多重积分难以计算，则无法直接求解 $\theta$ ，但是可以生成一系列符合此概率的分布。在已知样本的情况下，由贝叶斯公式可知,

P(\theta|X) = P(\theta)P(X|\theta)/P(X) \propto P(\theta)P(X|\theta)

因此，对后验概率的MCMC采样实际上是通过对 $P(\theta)P(X|\theta)$ 进行采样完成的。

更进一步地，如果似然函数本身也很复杂呢？这好像涉及到了Approximate Bayesian computation（ABC），有兴趣的可以参阅：
https://en.wikipedia.org/wiki/Approximate_Bayesian_computation

频率学派做法

无论是最小二乘法还是极大似然估计、矩估计的做法，对于噪声服从高斯分布的线性回归，其参数估计的结果应该是一样的。我们以之前介绍过的OLS为例，可以看到估计的参数还是与实际的参数非常接近的：

x = np.linspace(0, 10, 500)
y = 4 + 2 * x + np.random.standard_normal(len(x)) * 2
reg = np.polyfit(x, y, 1)
print(reg)
plt.figure(figsize=(10, 6))
plt.scatter(x, y, c=y, marker='v', cmap='coolwarm')
plt.plot(x, reg[1] + reg[0] * x, lw=2.0)
plt.colorbar()
plt.xlabel('x')
plt.ylabel('y')
plt.show()

[1.97451806 4.1131656 ]

output_66_1

PyMC3实现贝叶斯回归

通过PyMC3，Python生态系统提供了一个全面的软件包，可以在技术上实现贝叶斯统计和概率编程。接下来，我们利用PyMC3包进行贝叶斯回归。

PyMC3包

在介绍PyMC3之前，首先要提一下概率编程。所谓概率编程，简而言之，就是用计算机编程的方法实现概率分布的自动化推理。PyMC3就是Python的概率编程包，它允许用户使用各种数值方法拟合贝叶斯模型，其中最著名的就是马尔科夫蒙特卡洛（MCMC）和变分推理（variational inference, VI）。

PyMC从2003年开始开发，最开始旨在构建Metropolis-Hastings采样算法，推广MCMC的应用。PyMC3做了重大更新，使用Theano作为后端，能够支持更多基于梯度的MCMC方法（包括Hamiltonian Monte Carlo (HMC), the No U-turn Sampler (NUTS), and Stein Variational Gradient Descent）和变分推理方法（包括自动微分变分推理（ADVI）和算子变分推理（OPVI））。未来的PyMC4将以TensorFlow作为计算后端。

以上源自其官网（https://docs.pymc.io/history.html ）介绍

具体来说，我们会用到PyMC3中的以下几个方法：

find_MAP()通过最优化方法找到局部后验最大的点估计值作为采样算法的起点
NUTS()实现了新一代的MCMC采样算法“efficient No-U-Turn Sampler with dual averaging” (NUTS)
sample()利用find_MAP()确定的起始点和NUTS()确定的采样步长进行采样

具体代码如下，有机会再详细学一下这个的用法：

1	import pymc3 as pm

%%time
with pm.Model() as model:
    # For the regression line y_hat=alpha+beta*x
    # Define the priors
    alpha = pm.Normal('alpha', mu=0, sd=20)
    beta = pm.Normal('beta', mu=0, sd=10)
    sigma = pm.Uniform('sigma', lower=0, upper=10)

    # Specify the regression line
    y_hat = alpha + beta * x

    # Define the likelihood function
    # For the likelihood, assume a normal distribution with a mean of y_hat
    # and a uniformly distributed standard deviation of between 0 and 10, i.e. sigma.
    likelihood = pm.Normal('y', mu=y_hat, sd=sigma, observed=y)

    # Inference
#     start = pm.find_MAP() # Finds the starting value by optimization.
#     step = pm.NUTS()  # Instantiates the MCMC algorithm.
    trace = pm.sample(100, tune=1000, progressbar=True, verbose=False)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sigma, beta, alpha]
Sampling 4 chains: 100%|██████████████████████████████████████████| 404000/404000 [05:36<00:00, 1200.61draws/s]


Wall time: 6min 39s

1	pm.summary(trace)

	mean	sd	mc_error	hpd_2.5	hpd_97.5	n_eff	Rhat
alpha	4.115037	0.176746	0.004062	3.781424	4.467545	1937.477138	1.000214
beta	1.974356	0.030586	0.000737	1.915640	2.034231	1716.633197	1.000227
sigma	1.971559	0.064254	0.001298	1.847498	2.098978	2352.041599	1.001912

1 2	pm.traceplot(trace, lines={'alpha': 4, 'beta': 2, 'sigma': 2}) plt.show()

output_78_0

1	range(0, len(trace), 5)

range(0, 1000, 5)

plt.figure(figsize=(15, 9))
plt.scatter(x, y, c=y, marker='v', cmap='coolwarm')
plt.colorbar()
plt.xlabel('x')
plt.ylabel('y')
for i in range(0, len(trace), 250):
    plt.plot(x, trace['alpha'][i] + trace['beta'][i] * x)

output_80_0

比较最小二乘法的结果、贝叶斯回归的结果和真值，我们发现最小二乘法和贝叶斯回归的结果是非常接近的，但是与真值仍然有一定的差距（尤其是alpha），这主要是由于样本数目有限和噪声的存在造成的。

真实数据示例

接下来，我们使用真实的数据进行分析，我们选择两个ETF：黄金ETF——GLD和黄金股ETF——GDX，显然两者明显具有相关性，这就是我们进行线性回归的基础。

etf_gdx = pddata.DataReader(name='GDX', data_source='yahoo',
                              start='2010-01-01', end='2019-05-05',
                              retry_count=3, pause=0.1,
                              session=None, access_key=None)
etf_gld = pddata.DataReader(name='GLD', data_source='yahoo',
                              start='2010-01-01', end='2019-05-05',
                              retry_count=3, pause=0.1,
                              session=None, access_key=None)
data = {'gdx': etf_gdx['Adj Close'],
        'gld': etf_gld['Adj Close']
       }
data = pd.DataFrame(data).dropna()
data = (data / data.iloc[0] * 100)
data.head()

	gdx	gld
Date
2009-12-31	100.000000	100.000000
2010-01-04	103.246033	102.320385
2010-01-05	104.241496	102.227192
2010-01-06	106.773410	103.913899
2010-01-07	106.254056	103.270899

1 2	data.plot(figsize=(15,9)) plt.show()

output_85_0

由于我们要做的是GDX和GLD之间的回归，所以应该看一下两者的散点图，我们以GLD指数为横轴，GDX指数为纵轴作图，并利用颜色的渐变表示日期：

# Convert thge DatetimeIndex to matplotlib dates
mpl_dates = mpl.dates.date2num(data.index.to_pydatetime())
plt.figure(figsize=(15, 9))
plt.scatter(data['gld'], data['gdx'], c=mpl_dates, marker='o', cmap='coolwarm')
plt.xlabel('GLD')
plt.ylabel('GDX')
plt.colorbar(ticks=mpl.dates.DayLocator(interval=250),
             format=mpl.dates.DateFormatter('%d %b %y'))
plt.show()

output_87_0

从图上来看，最近6年以来，GDX和GLD有一定的正相关关系，接下来我们就实现两者的贝叶斯线性回归。

%%time
with pm.Model() as model:
    # For the regression line y_hat=alpha+beta*x
    # Define the priors
    alpha = pm.Normal('alpha', mu=0, sd=20)
    beta = pm.Normal('beta', mu=0, sd=10)
    sigma = pm.Uniform('sigma', lower=0, upper=30)

    # Specify the regression line
    y_hat = alpha + beta * data['gld']

    # Define the likelihood function
    # For the likelihood, assume a normal distribution with a mean of y_hat
    # and a uniformly distributed standard deviation of between 0 and 10, i.e. sigma.
    likelihood = pm.Normal('y', mu=y_hat, sd=sigma, observed=data['gdx'])

    # Inference
#     start = pm.find_MAP() # Finds the starting value by optimization.
#     step = pm.NUTS()  # Instantiates the MCMC algorithm.
    trace = pm.sample(100, tune=1000, progressbar=True, verbose=False)

Only 100 samples in chain.
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sigma, beta, alpha]
Sampling 4 chains: 100%|███████████████████████████████████████████████| 4400/4400 [00:12<00:00, 346.99draws/s]
The acceptance probability does not match the target. It is 0.89675341184176, but should be close to 0.8. Try to increase the number of tuning steps.


Wall time: 1min 7s

1	pm.summary(trace)

	mean	sd	mc_error	hpd_2.5	hpd_97.5	n_eff	Rhat
alpha	-90.840964	3.126800	0.215263	-96.486945	-85.137745	154.696026	1.024367
beta	1.352514	0.025732	0.001747	1.305662	1.399505	158.630489	1.024574
sigma	21.670428	0.330073	0.019059	21.067100	22.316189	218.890762	0.996822

1 2	pm.traceplot(trace) plt.show()

output_91_0

plt.figure(figsize=(15, 9))
plt.scatter(data['gld'], data['gdx'], c=mpl_dates,
         marker='o', cmap='coolwarm')
plt.xlabel('GLD')
plt.ylabel('GDX')
for i in range(len(trace)):
    plt.plot(data['gld'],
             trace['alpha'][i] + trace['beta'][i] * data['gld'])
plt.colorbar(ticks=mpl.dates.DayLocator(interval=250),
          format=mpl.dates.DateFormatter('%d %b %y'))
plt.show()

output_92_0

时变参数与随机游走

从上面的图我们可以发现，虽然GDX和GLD存在一定的正相关，但是两者的关系是随时间变化的，也就意味着参数可能是时变的。如何描述这种参数漂移的现象呢？一种做法是采用随机游走模型来描述参数的变化：

\alpha_t \sim N(\alpha_{t-1}, \sigma_\alpha^2)\\ \beta_t \sim N(\beta_{t-1}, \sigma_\beta^2)

与时不变的模型相比，我们又多出来 $\sigma_\alpha^2$ $\sigma_\beta^2$ 两个参数，那么对 $\alpha$ 和 $\beta$ 也需要采样：

subsample = 1
counts = int(len(data) / subsample)
total = subsample * counts

model_randomwalk = pm.Model()
with model_randomwalk:
    # std of random walk
    sigma_alpha = pm.Exponential('sig_alpha', 50., testval=.1)
    sigma_beta = pm.Exponential('sig_beta', 50., testval=.1)

    alpha = pm.GaussianRandomWalk('alpha', sd=sigma_alpha,
                                  shape=counts)
    beta = pm.GaussianRandomWalk('beta', sd=sigma_beta,
                                 shape=counts)

然后，再得到数据的似然函数：

with model_randomwalk:
    # Bring the parameter vectors to interval length.
    alpha_r = np.repeat(alpha, subsample)
    beta_r = np.repeat(beta, subsample)
    # Define the regression model.
    regression = alpha_r + beta_r * data['gld'].values[:total]
    sd = pm.Uniform('sd', 0, 20)   # The prior for the standard deviation.
    likelihood = pm.Normal('GDX', mu=regression, sd=sd,
                           observed=data['gdx'].values[:total])

%%time
import scipy.optimize as sco
with model_randomwalk:
    trace_rw = pm.sample(100, tune=2000, progressbar=True)

Only 100 samples in chain.
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sd, beta, alpha, sig_beta, sig_alpha]
Sampling 4 chains: 100%|████████████████████████████████████████████████| 8400/8400 [27:52<00:00,  2.91draws/s]
The acceptance probability does not match the target. It is 0.8833541435609861, but should be close to 0.8. Try to increase the number of tuning steps.
The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.
The acceptance probability does not match the target. It is 0.5117834416794408, but should be close to 0.8. Try to increase the number of tuning steps.
The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.
The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.
The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.
The gelman-rubin statistic is larger than 1.4 for some parameters. The sampler did not converge.
The number of effective samples is smaller than 10% for some parameters.


Wall time: 28min 52s

1 2	sh = np.shape(trace_rw['alpha']) sh

(400, 2350)

我们可以看一下，如果认为参数服从随机游走的话，参数以及拟合的回归方程随时间变化的情况：

import datetime as dt
part_dates = np.linspace(min(mpl_dates),
                         max(mpl_dates), sh[1])
index = [dt.datetime.fromordinal(int(date)) for
         date in part_dates]
alpha = {'alpha_%i' % i: v for i, v in
         enumerate(trace_rw['alpha']) if i < 100}
beta = {'beta_%i' % i: v for i, v in
        enumerate(trace_rw['beta']) if i < 100}
df_alpha = pd.DataFrame(alpha, index=index)
df_beta = pd.DataFrame(beta, index=index)
df_alpha.plot(color='b', alpha=0.05, legend=False,figsize=(15, 9), title='alpha')
df_beta.plot(color='r', alpha=0.05, legend=False,figsize=(15, 9), title='beta')
plt.show()

output_103_0

output_103_1

plt.figure(figsize=(15, 9))
plt.scatter(data['gld'], data['gdx'], c=mpl_dates, marker='o', cmap='coolwarm')
plt.colorbar(ticks=mpl.dates.DayLocator(interval=250),
          format=mpl.dates.DateFormatter('%d %b %y'))
plt.xlabel('gld')
plt.ylabel('gdx')
x = np.linspace(min(data['gld']), max(data['gld']))
for i in range(0,sh[1],50):  
    alpha_rw = np.mean(trace_rw['alpha'].T[i])
    beta_rw = np.mean(trace_rw['beta'].T[i])
    plt.plot(x, alpha_rw + beta_rw * x, '--', lw=0.7, alpha=.7, color=plt.cm.coolwarm(i / sh[1]))

output_104_0

虽然我们是用时变参数的方法完成了估计，但是这样实际上有可能会产生过拟合的问题：两者之间毕竟没有特别强的线性关系，许多其他的因素并没有考虑到其中。这些参数的变化完全有可能是因为其他因素的干扰。

风格轮动

而两者的关系可能是多模态的：10年的数据与12年以后的数据明显服从两个不同的正相关关系，中间的少数无规律的散点可以认为是两个模态切换的过渡过程。这也让我想到了股市里面的风格轮动，之前找工作的时候曾经用R写过用马尔科夫跳变系统描述风格轮动，有机会的话可以整理一下。