非参数贝叶斯波动率估计

自然结构θk |θk-1.~ IG（α，θk-1/α). 另一张是S.GUGUSHVILI、F.H.VAN DER MEULEN、M.SCHAUER和P.SPREIJ100 200 300 400 10-410-21001021041061081010α = 40.0, αζ= 20.0100 200 300 40010-410-21001021041061081010α = 30.0, αζ= 30.0100 200 300 40010-410-21001021041061061081010α=20.0，αζ=40.0图1。马尔可夫链{θk}的实现，α=40，αζ=20（左面板）和α=30，αζ=30（中面板）和α=20，αζ=40（右面板）。在所有情况下，θ固定为500。另一方面，可以通过设置θk |θk实现θk之间的正相关-1.~ IG（α，（αθk-1)-1） ，但这会导致难以处理的后验计算。在本研究中，使用潜在变量ζk对IGMC的定义优先考虑了这两个要点。更多讨论请参见（Cemgil和Dikmen，2007）。备注3。将漂移系数Bt设置为零有效地导致假设过程X具有独立的（高斯）增量。实际上，由于实际应用中的偏差通常为非零，因此过程的增量是相关的，因此所有观测值都包含一些关于每个时间点t的波动率值的间接信息∈ [0，T]。另一方面，假设在syields上的IGMC优先于系数{θk}的后验依赖性，这将有助于用较小的样本量n进行推断。有关说明，请参见第4节。2.3. 吉布斯取样器。可通过采用（5）的直接计算验证θk和ζk的完整条件分布为反伽马、θk |ζk、ζk+1~ IG公司α+αζ，αζk+αζk+1, k=2，N- 1,(6)θ|ζ~ IG公司α+ αζ, α+αζζ,（7） θN |ζN~ IG公司α、 αζN,（8） ζk |θk，θk-1.~ IG公司αζ+α，αζθk-1+αθk, k=2，N、 （9）详见第6节。

接下来，马尔可夫链{θk}的有效转移核由（Cemgil和Dikmen，2007）中的公式（4）给出，是逆伽马分布的比例混合；然而，它的准确表达对我们的目的没有直接关系。如（Cemgil和Dikmen，2007）第700页所述，根据参数值α，αζ，链{θk}可能呈现增加或减少的趋势。我们通过在图1中绘制{θk}实现α和αζ的不同值来说明这一点。在波动性估计的背景下，如果可以获得关于波动性趋势的先验信息，则该特征具有吸引力。（Cemgil和Dikmen，2007）中的推断是使用平均场变分Bayes方法进行的，参见（Blei等人，2017）。在此，我们描述了一种基于吉布斯抽样的全非参数贝叶斯波动率估计方法（参见，例如，（Gelfand Smith，1990）和（Geman and Geman，1984）），参见（Cemgil等人，2007）。算法在值ζ、…、处初始化，ζN，例如，由优先规范（5）生成。为了推导θk的全部条件的更新公式，definezk=kmXi=（k-1） m+1Yi，n，k=1，N- 1，ZN=nXi=（N-1） m+1Yi，n。使用此符号，来自（3）满足度（10）Ln（θ）的可能性∝ θ-（m+r）/2Nexp-nZN2TθNN-1Yk=1θ-米/平方公里-nZk2Tθk.使用这个公式和方程（6），并回顾反伽玛密度（2）的形式，可以看到θk，k=2，…，的更新分布，N- 1，isIGα+αζ+m，αζk+αζζk+1+nZk2T,而对于θ和θNareIGα+αζ+m，α+αζζ+nZ2T, IG公司α+m+r，αζN+nZN2T,分别地接下来，将使用公式（9）更新潜变量ζk。

除了通过θk’s的先前值之外，ζk’s的更新步骤不直接涉及数据Xn。最后，对θk’s和ζk’s重复这两个吉布斯步骤很多次（直到链可以被评估为合理收敛），这会给出θk’s的后验样本。使用后验推断可以按照通常的方式进行，例如，通过计算θk的样本后验平均值以及样本分位数，分别通过平方波动率s的边际贝叶斯可信带提供点估计和不确定性量化。后验样本平方根的类似计算可用于获得波动率函数s本身的点估计和可信带。2.4. 超参数选择。我们首先假设箱子的数量N以某种方式被控制，我们只需要处理超参数α、αζ和α，它们控制马尔可夫链的先验属性。在（Cemgil和Dikmen，2007）中，引入了IGMC先验知识，但未讨论超参数选择指南。在（Cemgil和Dikmen，2008）中，超参数是通过在那里研究的特定问题（音频去噪和单声道音频源分离）中手动调整的。另一个可行的解决方案是尝试超参数α、αζ和α的几种不同组合，如果只是为了验证不同结论对超参数变化的敏感性。Dikmen和8 S.GUGUSHVILI、F.H.VAN DER MEULEN、M.SCHAUER和P.SPREIJCemgil（2010年3月）讨论了一些进一步的超参数优化方法。在（Cemgil和Dikmen，2008）中，建议通过最大似然法优化超参数；实际实现依赖于EM算法（见（Dempster et al.，1977）），一些其他细节见（Dikmen和Cemgil，2008）。

换句话说，proposalin（Dikmen和Cemgil，2008）相当于使用经验Bayes方法（参见，例如，（Efron，2010），（Robbins，1956）和（Robbins，1964））。后者的使用非常广泛，通常会产生良好的实际效果，但理论上对该方法的理解仍然不够充分，除了白噪声模型之类的玩具模型（但是，参见（Donnet et al.，2018）和（Petrone et al.，2014）了解其他上下文中的一些结果）。在实际方面，在我们的例子中，给定序列{ζk}和{θk}的维数相当高，即2N- 1如果N较大，且闭合形式下的边际可能性不可用，则该方法预计将是计算密集型的。因此，先验地，没有理由不尝试通过为超参数配备先验来代替完全贝叶斯方法，而这实际上是我们目前工作中的默认方法。然而，相应的全条件分布结果是非标准的，这就需要在吉布斯采样器中使用Metropolis-Hastings步骤（参见，例如，（Hastings，1970），（Metropolis et al.，1953）和（Tierney，1994））。我们在第6节中提供了必要的详细信息。最后，我们简要讨论了超参数N的选择。如（Gugushvili et al.，2017jun）所述，在实践中，建议使用（Gugushvili et al.，2017jun）中的理论结果（建议 nλ/（2λ+1），如果真实波动率函数为λ-H¨older平滑），则同时尝试n的几个值。不同的N均提供未知波动率的信息，但分辨率不同；更多讨论见（Gugushvili等人，2017年6月）第5节。

正如我们将在第3节的模拟示例中看到的那样，对于N的选择，使用IGMC先验的推断结论非常可靠。这是因为通过超参数α和αζ，IGMC先验有一个额外的层来控制应用的平滑度；当α和αζ配备先验知识（如上所述）时，实际上可以从数据中学习它们。3、本节中的合成数据示例计算是用编程语言Julia完成的，请参见（Bezanson et al.，2017）。为了测试我们的估计方法的实际性能，我们使用了一个具有挑战性的例子，使用了来自（Donoho和Johnstone，1995）的blocks函数。作为第二个例子，我们考虑了Cox-Ross-Ingersoll模型的情况。以下各小节给出了准确的详细信息。我们在区间[0，1]上有800 001个等距点的网格上使用Euler格式，以获得（1）不同漂移系数和色散系数组合的解。然后对这些样本进行二次抽样，以获得每个示例中的4000个观察值。超参数α设置为0.1，而对于其他两个超参数，我们假设α=αζ，并在之前使用了一个diff-use-IG（0.3，0.3），以下特别提到的情况除外。使用第2节中的吉布斯采样器进行推断，使用Metropolis-Hastings步骤更新超参数α。后者使用了一个独立的高斯随机游走方案，并进行了缩放，以确保约50%的接受率；见第6节。Gibbs取样器运行了200000次迭代，我们使用了1000个样本的老化。在每个示例中，weNONPARAMETRIC贝叶斯波动率估计90.0 0.2 0.4 0.6 0.8 1.0-5051015X0.0 0.2 0.4 0.6 0.8 1.001020030S3图2。

过程X的样本路径来自（11）（左面板）和相应的波动率函数s（右面板）。绘制从系数ξk的中心后间期获得的95%边缘可信带=√θk.3.1。块功能。作为我们的第一个例子，我们考虑了波动率函数由来自（Donoho和Johnstone，1995）的blocks函数给出的情况。对于正性的垂直偏移，定义如下：（11）s（t）=10+3.655606×Xj=1hjK（t- tj），t∈ [0，1]，其中K（t）=（1+sgn（t））/2，{tj}=（0.1，0.13，0.15，0.23，0.25，0.4，0.44，0.65，0.76，0.78，0.81），{hj}=（4，-5, 3, -4, 5, -4.2, 2.1, 4.3, -3.1, 2.1, -4.2).该函数是非参数回归中一个具有挑战性的基准示例：它通常非常平滑，但在空间上不均匀，并且具有突变的特征（参见（Wasserman，2006）中的第9章）。与非参数回归不同，我们设置中的噪声（维纳过程）应被视为乘法，与s成比例，而不是加法，这与s取较大值的事实相结合，使推理问题更加复杂。我们这里的主要目标是比较基于IGMC先验的方法与基于IIG先验的方法的性能（Gugushvili等人，2017年6月）。为了完成设计规范，我们的漂移系数被选为相当强的线性漂移B（x）=-10倍+20倍。在图2中，我们绘制了块函数（11）和此模拟运行中使用的进程X的相应实现。图3的左侧和右侧面板对比了在N=160和α=αζ=20与N=320和α=αζ=40的情况下使用GMC之前获得的结果。这些曲线图表明，增加N会影响未平滑的先前实现，这可以通过增加αζ、α的值来平衡，而αζ、α具有相反的平滑效果。

因此，事实上，两个情节看起来非常相似。图4的左上和右上面板给出了使用IIG基于先验的方法获得的估计结果（Gugushvili等人，2017年6月）。箱数再次为N=160和N=320，在这两种情况下，我们对（平方）波动率函数的系数使用了不同的独立IG（0.1，0.1）先验值10 S.GUGUSHVILI、F.H.VAN DER MEULEN、M.SCHAUER和P.SPREIJ0.0.2 0.4 0.6 0.8 1.0020406080S3边际95%cred。频带0.0 0.2 0.4 0.6 0.8020406080S3边际95%cred。带状图3。使用N=160，α=αζ=20（左图）和N=320，α=αζ=40（右图），来自（11）的波动率函数s，具有IGMC先验的95%边缘可信带；在这两种情况下，α=0.1。（详见（Gugushvili等人，2017年6月）。这些结果必须与之前使用IGMC获得的结果进行对比，绘制在图4的左下和右下面板中，其中我们假设α=0.1和α=αζ~ IG（0.3，0.3）。

图4得出了以下结论：o尽管IGMC和IIG方法都在全球范围内恢复了波动率函数的形状，但IIG方法会导致推断结论中更大的不确定性，这反映在更广泛的边际置信区间中。在N=320的情况下，这种影响尤其明显，其中IIG先前的带宽几乎无法进行推断基于IGMC prior的频带看起来比IIG prior的频带更“规则”。o将结果与图3进行比较，我们看到了将超参数α、αζ与先验值相匹配的好处：图3中的可信带不能充分捕捉函数s在点t=0.2之前和之后的两个倾角，因为s完全落在可信带之外。因此，图3中使用了不正确的平滑量基于IIG先验的方法对箱子编号的选择很敏感：使用N=160个箱子比较图4的左上面板和使用N=320个箱子的右上面板，其中可信范围更宽。另一方面，基于IGMC先验的方法会自动重新平衡其使用的平滑量和不同数量的箱子N，这要感谢hyperprior在参数α、αζ上的作用；事实上，图4中底部的两个图看起来很相似。3.2. CIR模型。如前几节所述，我们的核心估计程序假设波动率函数是确定性的。然而，在本小节中，为了测试我们的方法的适用性限制和未来扩展的可能性，我们将其应用于波动率函数是随机的情况。

年（Mykland，2012年）的研究为这种方法提供了支持，但我们将重点放在实际方面，并将相应的理论研究推迟到另一个场合。非参数贝叶斯波动率估计110.0 0.2 0.4 0.6 0.8 1.0020406080S3边际95%可信区间。频带0.0 0.2 0.4 0.6 0.8020406080S3边际95%cred。频带0.0 0.2 0.4 0.6 0.8 1.0020406080S3边际95%cred。频带0.0 0.2 0.4 0.6 0.8020406080S3边际95%cred。带状图4。使用N=160（左上面板）和N=320个料仓（右上面板），具有IIG之前IG（0.1，0.1）的叠加95%边缘可信带的波动率函数s。使用N=160（左下面板）和N=320个仓位（右下面板），来自（11）的波动率函数与GMC之前的95%边际可信区间叠加；在这两种情况下，α=0.1和α=αζ~IG（0.3，0.3）。具体而言，我们考虑了Cox-Ross-Ingersoll（CIR）模型或平方根过程，（12）dXt=（η- ηXt）dt+ηpXtdWt，X=X>0，t∈ [0，T]。这里η、η、η>0是模型的参数。这一分化过程于（Feller，1951a）和（Feller，1951b）引入，并作为短期利率模型在金融领域广受欢迎，见（Cox et al.，1985）。条件2η>η确保了X的严格正性和遍历性。来自（1）的波动率函数现在对应于随机过程t7的实现→ η√Xt，其中X解CIR方程（12）。我们取任意参数值（13）η=6，η=3，η=2，x=1。图5的左面板绘制了X的样本路径，而相应的波动率在同一图的右面板中给出。在图6中，我们使用N=160和N=320bin以及超参数规格α=0.1和α=αζ显示了使用IGMC先验获得的估计结果~ IG（0.3，0.3）。

从该图中得出的结论是，我们的贝叶斯方法以一种相当令人满意的方式捕获了已实现波动率的总体形状。12 S.GUGUSHVILI、F.H.VAN DER MEULEN、M.SCHAUER和P.SPREIJ0.0.2 0.4 0.6 0.8 1.00123456X0.0 0.2 0.4 0.6 0.8 1.00123456s（ω）图5。过程X的样本路径来自（12）（leftpanel）和相应的实现波动率函数s（ω）（rightpanel）。参数值在（13）中给出。0.0.2 0.4 0.6 0.8 1.00123456s边际95%可信区间。带0.0 0.2 0.4 0.6 0.8 1.00123456s边缘95%cred。带状图6。使用N=160（左面板）和N=320个料仓（右面板），来自（12）的波动率函数s与IGMC先验的95%边际可信带叠加；在这两种情况下，α=0.1和α=αζ~ IG（0.3，0.3）。4、道琼斯工业平均指数在本节中，我们对时间序列变化点检测中的经典数据集进行了重新分析；例如，参见（Chen和Gupta，2012），（Diaz，1982），（Hsu，1977），（Hsu，1979）和（Iacus，2008）。具体而言，我们考虑1971年7月2日至1974年8月2日期间道琼斯工业平均指数的周收盘价。总的来说，有162个观测值可用，它们构成了一个相对较小的样本，因此推理问题相当重要。数据可作为sde包（见（Iacus，2016））中的数据集DWJ在R中访问（见（R Core Team，2017））。请参见图7的左面板以获取可视化效果。在（Iacus，2008）的weeklydata Xti中，i=1，n、 转换为返回Yti=（Xti-Xti公司-1） /Xti-1，并执行了（De Gregorio和Iacus，2008）的最小二乘变化点估计程序。在R中复制相应的计算机代码导致1973年3月16日的变化点估计。