使用R语言随机波动模型SV处理时间序列中的随机波动率

采样函数svsample期望其输入数据y是数字矢量,而没有任何缺失值(NA),如果提供其他任何内容,则会引发错误。

由Kaizong Ye,Liao Bao撰写

在y包含零的情况下,发出警告,并在进行辅助混合采样之前,将大小为sd(y)/ 10000的小偏移常数添加到平方收益上。

但是,我们通常建议完全避免零回报,例如通过预先降低零回报。下面是如何使用样本数据集exrates1准备数据的说明

×

MCMC概述

从名字我们可以看出,MCMC由两个MC组成,即蒙特卡罗方法(Monte Carlo Simulation,简称MC)和马尔科夫链(Markov Chain ,也简称MC)。要弄懂MCMC的原理我们首先得搞清楚蒙特卡罗方法和马尔科夫链的原理。我们将用三篇来完整学习MCMC。在本篇,我们关注于蒙特卡罗方法。

02

蒙特卡罗方法引入

蒙特卡罗原来是一个赌场的名称,用它作为名字大概是因为蒙特卡罗方法是一种随机模拟的方法,这很像赌博场里面的扔骰子的过程。最早的蒙特卡罗方法都是为了求解一些不太好求解的求和或者积分问题。比如积分:

如果我们很难求解出f(x)的原函数,那么这个积分比较难求解。当然我们可以通过蒙特卡罗方法来模拟求解近似值。如何模拟呢?假设我们函数图像如下图:

则一个简单的近似求解方法是在[a,b]之间随机的采样一个点。比如x0,然后用f(x0)代表在[a,b]区间上所有的f(x)的值。那么上面的定积分的近似求解为:

当然,用一个值代表[a,b]区间上所有的f(x)的值,这个假设太粗糙。那么我们可以采样[a,b]区间的n个值:x0,x1,…xn−1,用它们的均值来代表[a,b]区间上所有的f(x)的值。这样我们上面的定积分的近似求解为:

虽然上面的方法可以一定程度上求解出近似的解,但是它隐含了一个假定,即x在[a,b]之间是均匀分布的,而绝大部分情况,x在[a,b]之间不是均匀分布的。如果我们用上面的方法,则模拟求出的结果很可能和真实值相差甚远。

怎么解决这个问题呢? 如果我们可以得到x在[a,b]的概率分布函数p(x),那么我们的定积分求和可以这样进行:

上式最右边的这个形式就是蒙特卡罗方法的一般形式。当然这里是连续函数形式的蒙特卡罗方法,但是在离散时一样成立。

可以看出,最上面我们假设x在[a,b]之间是均匀分布的时候,p(xi)=1/(b−a),带入我们有概率分布的蒙特卡罗积分的上式,可以得到:

也就是说,我们最上面的均匀分布也可以作为一般概率分布函数p(x)在均匀分布时候的特例。那么我们现在的问题转到了如何求出x的分布p(x)对应的若干个样本上来。

03

条概率分布采样

上一节我们讲到蒙特卡罗方法的关键是得到x的概率分布。如果求出了x的概率分布,我们可以基于概率分布去采样基于这个概率分布的n个x的样本集,带入蒙特卡罗求和的式子即可求解。但是还有一个关键的问题需要解决,即如何基于概率分布去采样基于这个概率分布的n个x的样本集。 

对于常见的均匀分布uniform(0,1)是非常容易采样样本的,一般通过线性同余发生器可以很方便的生成(0,1)之间的伪随机数样本。而其他常见的概率分布,无论是离散的分布还是连续的分布,它们的样本都可以通过uniform(0,1)的样本转换而得。比如二维正态分布的样本(Z1,Z2)可以通过通过独立采样得到的uniform(0,1)样本对(X1,X2)通过如下的式子转换而得:

其他一些常见的连续分布,比如t分布,F分布,Beta分布,Gamma分布等,都可以通过类似的方式从uniform(0,1)得到的采样样本转化得到。在python的numpy,scikit-learn等类库中,都有生成这些常用分布样本的函数可以使用。

不过很多时候,我们的x的概率分布不是常见的分布,这意味着我们没法方便的得到这些非常见的概率分布的样本集。那这个问题怎么解决呢?

04

接受—拒绝采样

对于概率分布不是常见的分布,一个可行的办法是采用接受-拒绝采样来得到该分布的样本。既然 p(x) 太复杂在程序中没法直接采样,那么我设定一个程序可采样的分布 q(x) 比如高斯分布,然后按照一定的方法拒绝某些样本,以达到接近 p(x) 分布的目的,其中q(x)叫做 proposal distribution。

具体采用过程如下,设定一个方便采样的常用概率分布函数 q(x),以及一个常量 k,使得 p(x) 总在 kq(x) 的下方。如上图。

首先,采样得到q(x)的一个样本z0,采样方法如第三节。然后,从均匀分布(0,kq(z0))中采样得到一个值u。如果u落在了上图中的灰色区域,则拒绝这次抽样,否则接受这个样本z0。重复以上过程得到n个接受的样本z0,z1,…zn−1,则最后的蒙特卡罗方法求解结果为:

整个过程中,我们通过一系列的接受拒绝决策来达到用q(x)模拟p(x)概率分布的目的。

05

蒙特卡罗方法小结

使用接受-拒绝采样,我们可以解决一些概率分布不是常见的分布的时候,得到其采样集并用蒙特卡罗方法求和的目的。但是接受-拒绝采样也只能部分满足我们的需求,在很多时候我们还是很难得到我们的概率分布的样本集。比如:

1)对于一些二维分布p(x,y),有时候我们只能得到条件分布p(x|y)和p(y|x)和,却很难得到二维分布p(x,y)一般形式,这时我们无法用接受-拒绝采样得到其样本集。

2)对于一些高维的复杂非常见分布p(x1,x2,…,xn),我们要找到一个合适的q(x)和k非常困难。

从上面可以看出,要想将蒙特卡罗方法作为一个通用的采样模拟求和的方法,必须解决如何方便得到各种复杂概率分布的对应的采样样本集的问题。


准备数据

采样函数svsample期望其输入数据y是数字矢量,而没有任何缺失值(NA),如果提供其他任何内容,则会引发错误。在y包含零的情况下,发出警告,并在进行辅助混合采样之前,将大小为sd(y)/ 10000的小偏移常数添加到平方收益上。

但是,我们通常建议完全避免零回报,例如通过预先降低零回报。下面是如何使用样本数据集exrates1准备数据的说明。

图1提供了该数据集中时间序列的可视化。


热门课程

R语言数据分析挖掘必知必会

面对扑面而来的数据浪潮,包含Google、Facebook等国际企业,都已采用R语言进行数据分析

探索课程


除了现实世界的数据外,还可以使用内置的数据生成器svsim。此函数仅产生SV流程的实现,并返回svsim类的对象,该对象具有自己的print,summary和plot方法。

可以看出,该函数调用主MCMC采样器并将其输出转换为与coda兼容的对象。后者的完成主要是出于兼容性的考虑,并且可以直接访问在那里实施的收敛诊断检查。

svsample的返回值是svdraws类型的对象,该对象是具有八个元素的命名列表,其中包含(1)参数在para中绘制,(2)潜在的对数波动率,(3)初始潜在的对数波动率绘制latent0,(4)y中提供的数据,(5)运行时中的采样运行时,(6)先验中的先验超参数,(7)细化中的细化值,以及(8)这些图的汇总统计信息,以及一些常见的转换。


评估输出并显示结果

按照常规做法,可使用svdraws对象的print和summary方法。每个参数都有两个可选参数showpara和showlatent,用于指定应显示的输出。如果showpara为TRUE(默认设置),则会显示参数绘制的值/摘要。如果showlatent为TRUE(默认值),则显示潜在变量绘制的值/摘要。在下面的示例中,仅显示参数绘制的摘要。

(1)volplot:绘制潜在波动率的后分位数,以百分比表示,即随时间变化的100 exp(ht = 2)后验分布的经验分位数。除了必需的svsample-object本身之外,此函数还带有几个可选参数。
volplot或help(volplot)访问的相应帮助文档。常用的可选可选参数包括n步前波动率预测的预测,x轴上标签的日期以及一些图形参数。下面的代码片段显示了一个典型示例,图3显示了其输出。


(2)paratraceplot:显示θ中包含的参数的轨迹图。图5显示了一个示例。


 (3)paradensplot:显示θ中包含的参数的核密度估计。为了更快地绘制较大的后验样本,应将此参数设置为FALSE。如果参数showprior为TRUE(默认值),则先验分布通过虚线灰色线指示。图6显示了从汇率提取数据集中获得的EUR-USD汇率的示例输出。

svdraws对象的通用绘图方法将上述所有绘图合并为一个绘图。可以使用上述所有参数。请参见图7。

 为了提取标准化残差,可以在给定的svdraws对象上使用残差/残差方法。使用可选的参数类型,可以指定摘要统计的类型。当前,类型允许为“平均值”或“中位数”,其中前者对应于默认值。此方法返回svresid类的实向量,其中包含每个时间点所请求的标准化残差的摘要统计量。还有一种绘图方法,当通过参数origdata给定时,提供了将标准化残差与原始数据进行比较的选项。请参见下面的代码,对于相应的输出,请参见图8。


可下载资源

关于作者

Kaizong Ye拓端研究室(TRL)的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

​非常感谢您阅读本文,如需帮助请联系我们!

 
QQ在线咨询
售前咨询热线
15121130882
售后咨询热线
0571-63341498