R中的单变量和多变量随机波动率建模

更准确地说，让m表示横截面维度，相应的协方差矩阵∑t包含m（m+1）/2自由度，m中的一个二次项。表1说明了m的各种值的“维度诅咒”。打破这种诅咒的一种方法是使用潜在因子，从而实现∑t.3.1的稀疏表示。因子SV模型潜在因子模型体现了这样一种理念，即即使是高维系统也只受少数随机性来源的驱动。这几个随机性来源控制着几个因素，而这些因素又解释了观测值之间的相互作用。此外，潜在因子模型为动态协方差矩阵估计提供了有效工具。它们允许减少未知量的数量。带有r因子的传统潜在因子模型暗示了明确的运行时讨论，请参见Kastner和Frühwirth Schnatter（2014）以及Hosszejni和Kastner（2019）。有关在单个MCMC链中使用多核计算的可能性以及这样做时潜在的速度增益，请参阅Kastner（2019）。Darjus Hosszejni，Gregor Kastner 7m∑tper数据点的∑tfree元素1 110 55 5.5100 5050 50.51000 500 500 500.5表1：时变协方差矩阵∑tfor组件系列不同数量的∑tper数据点的∑tfree元素的绝对和相对数量m∑tper数据点的∑tfree元素10 44 4.4100 494.941000 4994 4.994表2：自由元素的绝对和相对数量时变协方差矩阵的元素∑tin a因子模型，用于不同数量的组分序列m和数量的因子r=4。分解∑t=_∑t+(R)∑t，（6），其中秩（_∑t）=r<m，且∑是包含异向误差方差的对角矩阵。

对称_∑tgives上的秩假设上升到因子分解_∑t=ψψ>，其中ψ∈ Rm×R包含mr- r（r- 1） /2个自由元素（参见Higham 1990中的pivotedCholesky算法）。因此，m（r+1）- r（r- 1） /2自由元素保留在∑t中，现在在m中仅呈线性。表2说明了m和r=4的各种值的“破坏维数诅咒”。下面，我们将介绍factorstochvol包中使用的FactorSV模型。我们对观测值yt=（yt1，…，ytm）>建模如下。yt |β，∧，ft，'∑t~ Nm（β+λft，(R)∑t），ft |∧∑t~ Nr（0，∑t），（7）式中，ft=（ft1，…，ftr）>是因子向量，β=（β，…，βm）>是观测规范平均值，且∧∈ Rm×ris一个包含因子载荷的高矩阵。协方差矩阵∑tand∑tare都是对角的，表示独立的vanilla-SV过程。“∑t=诊断（exp（”“ht1），exp（Phtm）），<<∑t=诊断（exp（Pht1，…，exp（Phtr）），<<hti~ N（(R)ui+(R)ui（(R)ht）-1，我- ui），？σi），i=1，m、 htj~ N（|uj+|uj（|ht）-1，j- uj），￠σj），j=1，r、 （8）为了从贝叶斯的角度对因子SV进行更为理论化的处理，读者可以参考Pitt和Shephard（1999）、Aguilar和West（2000）、Chib、Nardari和Shephard（2006）以及Han（2006）。基于方程7，我们可以将方程6重新表述为∑t=∧∧∑t∧>+(R)∑t，（9）8建模RFM中的随机波动性，其中几个识别问题是显而易见的：因子的顺序、符号和尺度是不确定的。更具体地说，对于大小为r×r的任何广义置换矩阵XP，我们发现另一个有效分解∑t=∧∧∧∑t（∧）>+\'∑t，其中∧=P-1和▄∑t=P▄∑tP>。我们通过将其对数方差水平调整为零，即j=1时的μj=0，…，来解决因子尺度的模糊性，r、 符号和顺序识别可以通过限制因子载荷矩阵∧来实现。

在限制∧的系数TochVol中有几个选项可用，详情见第5.2.3.2节。需要为平均值、潜在对数方差过程和因子载荷矩阵∧指定先验分布先验值。我们选择βj~ N（bβ，bβ），独立于j=1，m、 对于Bβ的较小值，这会使βj向Bβ收缩；对于Bβ的大值，先验知识是相当缺乏信息的。对数方差过程与第2.2节中的单变量情况具有相同的先验规范。对于∧，目前在factorstochvol中实现了三种类型的先验。所有三个都可以用∧ij的形式书写~ N（0，τij）独立于每个适用i∈ {1，…m}和j∈ {1，…，r}。首先，我们可以先验地确定所有τij——不一定是相同的值。这将为载荷矩阵的每个元素生成一个正常的先验。第二种类型是一种等级先验，其已被开发用于诱导更灵活和潜在更强的收缩，即∧ij |τij~ N（0，τij），τij |λi~ G（a，aλi/2）。（10） 该分布被Griffn和Brown（2010）称为正态伽马先验分布，意味着传统方差V（λij |λi）为2/λi，无条件的过剩峰度为3/a。a的值被视为一个结构参数，由用户确定，其中选择最小值（.1）将强制向零收缩，而选择较大值（&1）将强制收缩较小。案例a=1是一种特殊案例，称为贝叶斯套索先验（Park and Casella，2008）。根据λi的数据估计参数λiis~ G（c，d）。第三种类型是对第二种类型的轻微修改。

由于因子荷载矩阵∧每行中的方差可被视为来自同一基础分布的“随机效应”，因此方程10中的先验知识会导致行收缩，并进行元素自适应。类似地，也可以考虑采用单元自适应的柱状收缩，即∧ij |τij~ N（0，τij），τij |λj~ G（a，aλj/2），（11）与相应的先验λj~ G（c，d）。3.3. 估计因子SV模型中的贝叶斯估计建立在stochvol中的单变量vanilla SV实现基础上，具有多个效率提升级别。为了缓解高维潜在的缓慢收敛问题，通过利用ASIS的几种变体的采样器进行。在factorstochvol中实现的采样细节在广义置换矩阵中进行了描述，该矩阵具有置换矩阵的零-非零模式，但允许其具有任何非零值，而不仅仅是一个值。因此，广义置换矩阵不仅可以置换，还可以缩放和切换其乘法器的符号。Darjus Hosszejni，Gregor Kastner 9Kastner et al.（2017年，因子负荷使用高斯先验）以及Kastner（2019年，因子负荷使用分层收缩先验）。与stochvol类似，为了使计算时间即使在高维情况下也能承受，factorstochvol的主采样器是用C++编写的。它使用R包rcpps来简化R和C++之间的通信。单变量SV部件借用自Stochvol，并通过其C/C++-级更新函数update\\u fast\\u SV（）进行接口。这样，就避免了在每次MCMC迭代时在解释的R代码和编译的C++代码之间移动。stochvol软件包stochvol软件包分别通过其采样例程svsample（）、svtsample（）、svlsample（）和svtlsample（）提供了设置单变量SV、SVt、SVl和SVtl模型的方法。

在下文中，我们介绍了建议使用stochvol的工作流。首先，我们讨论使用默认设置的估计、可视化和预测。然后，我们展示了如何调整先验超参数的值以及如何配置采样机制。4.1. 准备数据并运行MCMC采样器我们估计了三个模型，它们举例说明了stochvol的功能和用户界面。利用一揽子计划中的汇率数据，我们以三种不同的方式对2008年3月1日至2012年3月1日期间的欧元兑瑞士法郎汇率（1欧元兑瑞士法郎的价格）进行建模（1028个数据点）。具有SV残差的AR（1）模型第一个示例是具有SV残差的AR（1）模型，即方程式1变为yt | yt-1，β，β，ht~ N（β+βyt-1，膨胀（ht）），ht+1，ht~ N（u+Д（ht- u), σ).利用该模型，我们检验了汇率是否遵循SV的随机游走。在这种情况下，我们预计β和β的后验值将分别集中在0和1左右。为了估计AR（1）-SV模型，我们需要将输入y准备为长度n的数字序列，并将其作为第一个输入参数传递给svsample（），如下所示：R>set。种子（1）R>库（“stochvol”）R>数据（“exrates”）R>索引<-其中（exrates$日期>=截止日期（“2008-03-01”）&+exrates$日期<=截止日期。日期（“2012-03-01”））R>CHF\\U价格<-exrates$CHF[ind]R>res\\U sv<-svsample（CHF\\U价格，designmatrix=“ar1”）我们设置designmatrix=“ar1”以使用AR（1）规范。

更一般地，对于常数平均模型，DesignMatrix可以采用“ar0”形式的字符值，对于AR（1）、AR（2）等，可以采用“ar1”、“ar2”等形式的字符值。10用SVt残差对RConstant均值模型中的随机波动性进行建模第二个例子是一个带有SVt残差的常数均值模型，即方程2 becomesyt |β，ht，ν~ tν（β，exp（ht/2）），ht+1，ht~ N（u+Д（ht- u), σ).如果返回是重尾的，则ν的大部分后部质量集中在低值上，例如小于20。否则，几乎没有证据表明存在高峰度。我们通过在之前计算的CHF\\U价格上应用logret（）来计算对数回报。然后，为了估计具有重尾SV残差的常数平均模型，我们将对数返回向量传递给svtsample（），并将designmatrix设置为“ar0”。R> 设置。seed（2）R>CHF\\u logret<-100*logret（CHF\\u price）R>res\\u svt<-svtsample（CHF\\u logret，designmatrix=“ar0”）SVl residuals多元回归第三个示例是一个带有截距、两个回归器和SVlresiduals的多元回归模型；也就是说，等式3变成了ytht+1！ht，ζ，xt1xt2！，β,ββ!~ Nβ+βxt1+βxt2u+Д（ht- u)!, Σρ!,∑ρ=exp（ht）ρσexp（ht/2）ρσexp（ht/2）σ！。为了举例说明，我们将欧元兑瑞士法郎的对数收益回归到同期的对数收益onEURUSD和EURJPY，分别为每美元1欧元和日元1欧元。为了使用stochvol估计多元回归模型，我们需要准备一个维数为n×K的数值矩阵，其中行对应于时间点，列对应于协变量。

我们创建一个截距作为X的第一列，并将第二列和第三列分别设置为欧元兑美元对数收益和欧元兑日元对数收益；最后，我们在多元回归中使用X列作为协变量。R> 设置。seed（3）R>X<-cbind（常数=1100*logret（汇率$美元[ind]），+100*logret（汇率$日元[ind]）R>res\\U svl<-svlsample（CHF\\U logret，designmatrix=X）4.2。将结果可视化通常，模型参数和潜在量的联合后验分布标志着贝叶斯分析的目标。为了检验它，我们可以查看汇总统计数据和各种类型的边际后验分布可视化。此外，建议检查Darjus Hosszejni，Gregor Kastner 11估计的波动率百分比（5%/50%/95%后分位数）50 1002008-12-11 2009-09-25 2010-07-08 2011-04-14 2012-01-245000 15000-6.-2 mu（thin=1）5000 150000.985 0.995 phi（thin=1）5000 150000.06 0.10 0.14 sigma（thin=1）5000 15000跟踪-0.4-0.1 0.2 rho的轨迹（薄=1）-6.-2 0 2 4 60.0 0.2 0.4 mu0.985 0.990 0.995 1.0000 100 250 phi0.06 0.10 0.140 10 20 30 sigma密度-0.5-0.3-0.1 0.10 2 Rho的密度图1：估计模型的默认图。顶行显示每日波动率（百分比）100 exp（h/2）的后验值（通过其中位数（黑色）以及5%和95%分位数（灰色）的汇总。其余面板总结了参数u、Д、σ和ρ的马尔可夫链。特别是，中间一行显示迹线图，底部一行显示先前（灰色，虚线）和之后（黑色，实心）密度。马尔可夫链用于解决可能的收敛问题——这通常通过研究后验量的跟踪图来实现。因此，受coda包的启发，stochvol提供了R泛型函数plot（）和summary（）的自身实例。

为了介绍stochvol提供的用于分析MCMC输出的工具，我们在本节剩余部分简要检查了第三个示例（带SVl误差的多元回归）的结果。首先，我们绘制估计的输出。R> 绘图（res\\u svl，showobs=FALSE，dates=exrates$date[ind[-1]]）结果如图1所示。我们在第一行中看到了波动率后验密度的总结。除了中位数外，我们还通过每个时间点的5%和95%分位数获得不确定性的量化。在第二行中，我们可以跟踪SV参数的马尔可夫链的演化。在本例中，它们是u、Д、σ和ρ。最后，我们分别在第三行ingray和black中看到参数的先验密度图和后验密度图。它们表现出高度的持久性和显著的杠杆作用。接下来，我们观察AR系数。R> 对于（i in seq\\u len（3））{+尾码：：traceplot（svbeta（res\\u svl）[，i]）+尾码：：densplot（svbeta（res\\u svl）[，i]，show.obs=FALSE）12 R5000 10000 15000 20000中的随机波动率建模-0.06 0.00迭代-0.06-0.04-0.02 0.00 0.020 10 30N=20000带宽=0.0015365000 10000 15000 20000-0.15-0.05迭代次数-0.20-0.15-0.10-0.05 0.000 5 10N=20000带宽=0.003135000 10000 15000 200000.15 0.25迭代0.15 0.20 0.25 0.300 5 15N=20000带宽=0.002825图2：p（β| y）后验提取的跟踪图和估计核密度+}结果如图2所示。在左侧，我们在跟踪图中没有发现任何收敛或混合问题的迹象。

在右侧，我们看到β、β和β的后验概率均不集中在0左右，因此协变量似乎对因变量有影响。作为最后一步，我们将打印估算结果的数字摘要。R> 摘要（res\\u svl，show潜伏=FALSE）#########svdraws对象摘要######先验分布：##mu ~正态（均值=0，sd=100）###（phi+1）/2 ~ Beta（a=5，b=1.5）#######σ^ 2 ~伽马（形状=0.5，速率=0.5）##nu ~无穷大##rho ~ Beta（a=4，b=4）##Beta ~多变量极大（…）######2000年磨合后存储20000张MCMC图纸无稀释。#####SV参数的后验图（细化=1）：##平均sd 5%50%95%ESSDarjus Hosszejni，Gregor Kastner 13###mu-1.2554 1.4048-2.9497-1.5371 1.727 113##phi 0.9973 0.0023 0.9927 0.9979 1.000 376#sigma 0.0855 0.0123 0.0668 0.0846 0.108 344#rho-0.1290 0 0.0877-0.2731-0.1307 0.014 226##exp（mu/2）0.7458 0.9817 0.2288 0.4637 2.372 113##sigma^2 0.0075 0.0022 0.0045 0.0072 0.012344#####回归系数的后验图（细化=1）：##平均sd 5%50%95%ESS##beta\\u 0-0.025 0.011-0.042-0.025-0.0078 9506##beta_1-0.095 0.021-0.130-0.095-0.0599 3703##beta_2 0.227 0.019 0.195 0.227 0.2589 2452为了简洁起见，我们将show潜伏=FALSE设置为不打印所有1027个潜伏状态。输出显示了老化的长度和提取的数量，参数的先前规格，以及参数u、Д、σ和ρ的边际后验分布的简明摘要，此外还显示了波动率exp（u/2）和σ的水平，以及回归系数β的向量。此后验汇总表由后验平均值和标准偏差、5%、50%和95%分位数的列组成。

用户可以通过参数分位数将0到1之间的值序列传递给svsample（）、svtsample（）、svlsample（）或svtlsample（），从而影响显示的分位数。表中最后一列描述了所谓的有效样本量（ESS），这是一种衡量聚合MCMC链质量的指标。形式上，马尔可夫链C的ESS通过/（1+2P）定义∞s=1ρe fff（s）），其中M是C的长度，ρe fff（s）表示C元素之间滞后s的自相关函数。原则上，ESS是具有与我们（边际）链相同蒙特卡罗误差的连续不相关链的样本量。直观峰值，这意味着ESS是获得的独立且分布相同的绘图的数量，并给出了我们的链对后空间探索的程度。ESS值越高，混合效果越好。4.3. Stochvol预测我们使用我们的估计模型来预测数据集中剩余天数的日志收益。为此，我们首先准备接下来24天的协变量，并通过通用predict（）函数的argumentnewdata和估计输出传递它们。请注意，我们需要25天的价格数据才能获得24次回报。R> 设置。种子（4）R>pred\\u ind<-seq（尾部（ind，1），长度。out=25）R>pred\\u X<-cbind（常数=1100*logret（汇率$美元[pred\\u ind]），+100*logret（汇率$日元[pred\\u ind]）R>pred\\u svl<-predict（res\\u svl，24，newdata=pred\\u X）由于我们可以访问未来日志回报的整个分布，我们可以通过分位数量化预测的不确定性。在下面的代码片段中，我们可视化了k=1，…，的k步超前预测分布，24，以及真正观察到的值。