开始使用粒子Metropolis Hastings进行推理

2010）算法，其中还包括吉布斯采样的粒子版本（Andrieu等人2010；Lindsten等人2014）。PMCMC算法是当前研究的一个主题，已经给出了许多改进其性能的措施，有关这些措施的一些示例，请参见第6节。在本教程中，我们使用了第2.5节中讨论的u的独立提案。这本质上意味着所有粒子过滤器都是独立的。然而，直观地说，通过关联粒子滤波器可能会有所收获，因为对于θ的微小变化，状态估计通常非常相似。结果表明，关联u会导致对数似然估计值的正相关，从而降低接受概率估计值的方差。在实践中，这意味着N不需要像u在迭代之间不相关时那样快速地用T进行缩放。这在T较大时尤其有用，因为它降低了PMH的计算成本。有关更多信息和源代码，请参见Dahlin等人（2015a）、Choppala等人（2016）和Deligiannidis等人（2018）。5.应用示例：估算OMXS30中的波动率我们继续具体应用PMH算法来推断随机波动率（SV；Hull and White 1987）模型的参数。这是一个非线性SSM，从PMH开始，在非线性动态Mo-delsGaussian噪声中进行推理，这种模型中的推理是一个重要问题，因为对数波动率（该模型中的潜在状态）有助于风险管理和各种金融合同的定价。更多信息请参见Tsay（2005）和Hull（2009）。SV模型的具体参数由byx给出~ N十、u，σv1- φ, （26a）xt+1 | xt~ Nxt+1；u+φ（xt）- u），σv, （26b）yt | xt~ Nyt；0，exp（xt）, （26c）其中参数由θ={u，φ，σv}表示。

给你，先生∈ R、 φ∈ [-1,1]和σv∈R+分别表示状态过程的平均值、持续性和标准偏差。请注意，该模型与LGSS模型非常相似，但这里的状态变量会缩放观测噪声的方差。因此，我们有零均值的高斯观测和exp（xt/2）给出的依赖于状态的标准偏差。在计量经济学中，波动率是标准差的另一个词，因此xtis被称为对数波动率。这个模型中的测量值是所谓的对数返回，yt=100 logstst-1.= 100日志（st）- 原木（圣-1),其中，St表示某些金融资产（如指数、股票或商品）在时间t的价格。这里，{st}Tt=1是纳斯达克OMXS30指数的每日收盘价，即斯德哥尔摩证券交易所交易量最大的30只股票的加权平均值。我们从Quandl提取2012年1月2日至2014年1月2日期间的数据。图5的顶部显示了由此产生的日志返回。请注意日志返回中的时变持续波动性，即小变化和大变化周期。这被称为挥发性聚类效应，是SV模型旨在捕捉的真实世界数据的特征之一。从（26）来看，当|φ|接近1，且σvis较小时，可以实现这一点。因为这些选择会导致一阶自回归过程，具有很大的自相关性。本应用程序的目标是从观测数据y1:T估计参数θ和对数波动率x0:T。我们可以使用PMH作为两个样本，从参数的后部估计这两个量，并且可以在算法的每次迭代中获得状态。

为了完成SV模型，我们假设了一些基于参数通常范围的领域知识的参数先验，即p（u）=N（u；0，1），p（φ）=tn[-1,1]（φ；0.95,0.05），p（σv）=G（σv；2,10）。这里，G（a，b）表示形状为a且标度为b的伽马分布，即预期值a/b.5.1。修改实施粒子过滤器和PMH的实施需要适应这种新模型。我们通过替换两种算法框架中的部分代码，概述了必要的修改。结果实现和源代码可在函数中找到。数据可从以下网站下载：https://www.quandl.com/data/NASDAQOMX/OMXS30.JohanDahlin，Thomas B.Sch"on分别提出了23个微粒过滤器模型和微粒过滤器模型。在文章过滤器中，我们需要修改除重采样之外的所有步骤。在初始化过程中，我们需要通过替换：<initializeStates>=AncestorIndex[，1]<-1:NoParticle SnormizedWeights[，1]=1/NoParticle Sparticles[，1]<-rnorm（noParticles，mu，sigmav/sqrt（1-phi^2））来模拟初始粒子系统。此外，我们需要为粒子过滤器实现选择（粒子）建议分布和权重函数。对于SV模型，我们不能像LGSS模型那样计算最优选择的闭式表达式。

因此，采用了一个bPF，它对应于利用状态动力学作为byx（i）t | xa（i）tt给出的方案（21）-1.~ fθxtxa（i）tt-1.= Nxt；u + φxa（i）tt-1.- u, σv,观测模型作为权函数byW（i）t=gθytx（i）t= Nyt；0，expx（i）t.这两种选择导致估计器bxnt=πt[x]更改为xnt=NXi=1w（i）tx（i）t。对粒子过滤器的这三种更改是通过替换来实现的：<传播粒子>=part1<-mu+phi*（粒子[new祖先，t-1]-mu）粒子[，t]<-part1+rnorm（noParticles，0，sigmav）<weightParticles>=yhattrance<-exp（粒子[，t]/2）权重[，t]-dnorm（y[t-1]，yhatMean，yhatVariance，log=TRUE）xHatFiltered[t]<-sum（归一化权重[，t]*粒子[，t]）我们还需要推广PMH代码，使其具有多个参数。这很简单，我们建议读者参考源代码进行必要的更改。主要的变化是用矩阵θ和θ代替向量φ和φ。也就是说，马尔可夫链的状态现在是三维的，对应于θ中的三个元素。此外，参数建议被选为以先前参数θ（k）为中心的多变量高斯分布-1） 协方差矩阵∑=diag(), 哪里 表示包含三个元素的向量。这是通过以下方式实现的：<proposeParameters>=thetaProposed[k，]<-rmvnorm（1，mean=theta[k-1，]，sigma=stepSize）24开始使用PMH进行非线性动力学模型中的推理接受概率的计算也被修改，以包括新的先验分布。对于该模型，要求|φ|<1且σv>0，以确保其稳定且标准偏差为正。

这是通过以下方式实现的：<computeAcceptProbability>=if（（abs（thetapposed[k，2]）<1.0）和（（thetapposed[k，3]>0.0））{res<-particleFilterSVmodel（y，thetapposed[k，2]）logLikelihoodProposed[k]<-res$loglikelikelihoodxhatFilteredProposed[k，]<-res$xHatFiltered}priorMu priorMu（theaposed[k，1]，0，1，log=TRUE）priorMu-priorMu-priorMu-priorMu-loglikelikelikelikelikelikelikelihoodproposed[k，1]，TRUE]priorPhi<-dnorm（thetaProposed[k，2]，0.95，0.05，log=TRUE）priorPhi<-priorPhi-dnorm（theta[k-1，2]，0.95，0.05，log=TRUE）priorSigmaV<-dgamma（thetaProposed[k，3]，2，10，log=TRUE）priorSigmaV<-priorSigmaV<-priorMu[k，3]，2，10，log=TRUE）先验<-priorSigmaV[k]概率（Previor+likelihoodDifference）acceptProbability<-acceptProbability*（abs（thetaProposed[k，2]）<1.0）acceptProbability<-acceptProbability*（thetaProposed[k，3]>0.0）此外，<acceptRejectStep>需要进行修改，以考虑到theta是avector，有关详细信息，请参阅源代码。5.2. 估计对数波动率考虑到θ的不确定性，可以使用PMH的伪边缘视图来计算对数波动率的估计值。粒子滤波器是给定u的一种确定性算法。因此，随机变量u相当于滤波器分布的估计值。结果是，通过将状态估计包含在PMH生成的马尔可夫链中来计算状态估计非常简单。这是通过修改粒子过滤器的代码来实现的。在算法运行一次之后，我们在时间T以w（1:N）T给出的概率对一个粒子进行采样，然后沿着祖先谱系返回到T=0，并在状态空间中提取相应的路径。

这使我们能够在利用过去和未来信息时更好地估计对数波动率。因此，这通常会减少估计值的方差。这是通过使用存储的重采样祖先指数和采样状态轨迹来完成的，替换为：<returnEstimates>=ancestorIndex<-sample（noParticles，1，prob=normalizedWeights[，T]）sampledTrajectory<-cbind（ancestorIndex[ancestorIndex]，1:（T+1））xHatFiltered<-particles[sampledTrajectory]列表（xHatFiltered=xHatFiltered，loglikelibility=loglikelike）Johan Dahlin，Thomas B.Sch"on 250 100 200 300 400 500-4-20 0 2 4timelog-returns0 100 200 300 400 500-2-1 0 1 2 TimeLog波动率估计μ后验估计-1.0-0.5 0.0 0 0.5 0 0.5 1.0 1.0 52500 3000 3500 4000-1.0-0.0 0 0 0 0 0 0.5 1.0迭代μ后验估计0 20 40 60 80 100-0.2 0.2 0.6 1.0迭代。88 0.92 0.96 1.000 5 10 15 20 252500 3000 3500 40000.88 0.92 0.96 1.00迭代φ0 20 40 60 80 100-0.2 0.2 0.6 1.0迭代φσV后验估计0。0.1 0.2 0.3 0.40 2 4 6 8 102500 3000 3500 40000.0 0.1 0.2 0.3 0.4迭代σv0 20 40 60 80 100-0.2 0.2 0.6 1.0迭代σvFigure 5：上图：2012年1月2日至2014年1月2日期间纳斯达克OMXS30指数95%置信区间的每日对数收益率（深绿色）和估计对数波动率（橙色）。底部：从PMH获得的u（紫色）、φ（品红）和σv（绿色）的后验估计值（左）、马尔科夫链轨迹（中）和相应的ACF（右）。左图和中图中的虚线和实线分别表示参数后验平均值和参数前验值。26从PMH开始，在函数particleFilterSVmodel的非线性动力学模型中进行推理。

然后，以与PMH中候选参数相同的方式处理xHatFiltered中的采样状态轨迹，有关详细信息，请参阅源代码。因此我们可以通过以下方法计算参数和对数波动率的后验平均值及其相应的标准偏差：R>thestationary<-pmhOutput$theta[noburniterations:noIterations，]R>thetaHatMean<-colMeans（thestationary）R>thetahats标准偏差<-apply（thestationary，2，sd）R>xhatstational<-pmhOutput$xhatsfilteredmean[2500:7500，]R>xhatsfilteredmean<-colMeans（xhatstational）R> xhatFilteredStandardDeviation<-应用（xhatStational，2，sd），其中pmhOutput表示particleMetropolishingsSVModel的输出变量。5.3. 数值说明现在，我们使用本节开头介绍的OMXS30数据，将粒子过滤模型和粒子过滤模型应用于SV模型。完整的实现在函数/代码示例3_sv中提供。由此产生的后验估计值、轨迹和ACF如图7所示。我们看到，马尔可夫链和后验分布明显集中在后验平均估计Bθ={-0.23,0.97,0.15}，标准偏差{0.37,0.02,0.06}。这证实了我们的观点，即对数波动率是一个缓慢变化的过程，因为φ接近1，σ相对较小。图7的第二行还计算并给出了该参数下的对数波动率估计值。我们注意到，在t=100和t=370附近，波动性更大，这与对数收益率中的情况很好地对应。6.选定的高级PMH主题在本节中，我们概述了第5节中实施的一些可能的改进和最佳实践。我们将讨论初始化、收敛诊断以及如何改进所谓的马尔可夫链混合。

对于后者，我们考虑三种不同的方法：调整随机行走方案、重新参数化模型和调整粒子数。本节通过对更高级的提案分发和后处理的小型调查得出结论。6.1. 初始化在后验概率较大的区域初始化马尔可夫链以获得有效的算法非常重要。理论上，我们可以在参数空间的任何地方初始化链，并且仍然可以得到收敛的马尔可夫链。然而，在实践中，当参数后验值假设较小或对θ的变化相对不敏感时，可能会出现数值困难。因此，建议尝试获得参数的粗略估计，以初始化更接近后验模式的链。在LGSS模型中，我们可以使用标准优化算法结合卡尔曼滤波器来估计后验信号的模式。然而，这在一般SSM中是不可能的，因为只有后验噪声估计可以从粒子约翰·达林、托马斯·B·舍恩27滤波器中获得。因此，标准的优化方法可能陷入局部极小或无法完全收敛。缓解这一问题的一种方法是同步扰动和随机逼近（SPSA；Spall 1987）算法，该算法可与粒子滤波器结合使用，以估计后验模态。6.2. 诊断收敛一般很难证明马尔可夫链已达到其平稳状态，并且所获得的样本实际上是πθ（θ）的样本。一个简单的解决方案是在参数空间的不同点初始化算法，并比较产生的后验概率。

如果得出的估计值相似，这可能表明链已达到其稳定状态。另一种方法是利用Kolmogorov-Smirnov（KS；Massey 1951）测试来确定后验估计值在磨合后不会改变。这是通过将样品{θ（k）}Kk=1分为三个部分来实现的：老化和来自固定相的两组相等数量的样品。由于KS测试需要不相关的样本，可以使用细化来降低自相关。马尔可夫链通过保留每个样本来细化，其中选择l，使得两个保留的连续样本之间的自相关性可以忽略不计。然后进行标准KS测试，以得出两个减薄隔板中的样本是否来自相同（固定）分布。Robert和Casella（2009）和Gelman等人（2013）讨论了诊断收敛性的其他方法。改进的混合标准蒙特卡罗方法假设可以从感兴趣的分布或某些仪器分布中获得独立样本。然而，从PMH获得的样本是相关的，因为它们是由马尔可夫链生成的实现。直观地说，与独立样本相比，这些相关样本包含的目标分布信息更少。也就是说，如果马尔可夫链中的自相关较大，大多数样本都是相似的，并且目标没有完全被探索。这种自相关也被称为混合，这是MCMC文献中一个非常重要的概念。这是比较不同MCMC算法效率的关键数量。结果表明，混合通过缩放估计的渐近方差影响MCMC算法的性能。

这从（9）中控制bπKθ[~n]的CLT中可以明显看出，在某些假设下，它的形式是√Kπθ[φ] -bπKθ[~n]D-→ N0，VARπ[~n]·IACTθ（1:K）, K→ ∞, （27）当VARπ[~n]=πθ[（~n）时-πθ[φ])] < ∞. 给你，艾卡特θ（1:K）表示马尔可夫链的积分自相关时间（IACT），计算为ACF下的面积。IACT的一种解释是，它估计了MCMC算法在获得两个不相关样本之间的迭代次数。因此，如果这些假设包括马尔可夫链应该足够快地到达其静止状态，那么IACT就是一个。这一要求被称为一致遍历性，这意味着马尔可夫链的分布与其平稳分布之间的总变分范数应以几何速率减少，而不考虑马尔可夫链的初始化，更多细节参见例如Tierney（1994）。28非线性动力学模型中推理的PMH入门是独立的，例如重要性抽样。因此，最小化IACT与最大化MCMC算法中的混合相同。在图5的右栏中，我们展示了SVmodel中三个参数的ACF。该信息可用于计算相应的iact byIACT（θ（1:K））=1+2∞Xτ=1ρτθ（1:K）,式中ρτ=E[（θ（k）- πθ[θ]（θ（k+τ）- πθ[θ])]/πθ[(θ - πθ[θ]）表示θ（1:K）在滞后τ处的自相关系数。实际上，我们无法准确计算IACT，因为我们不知道所有可能滞后的自相关系数。此外，当τ较大且K相当小时，很难估计ρτ。这个问题的一个可能解决方案是将AFF估计值减少到一定的滞后，只考虑小于该极限L的滞后。