[量化金融] 使用阻尼BFGS更新构建Metropolis Hastings提案 [推广有奖]

Samp。pMH0-0.12-89±62 0.62 55pMH1-0.37-113±53 1.14 130dBFGS-0.76-24±7 2.16 53iBFGS FLIP 0.52 0.95 42±17 1.83 77iBFGS reg 0.45 0.95 34±15 1.77 61iBFGS hyb 0.60 0.95 40±17 1.84 76eBFGS hyb 0.62 1.00 34±19 1.89 64pMH0-0.04-745±177 0.3 211pMH1-0.25-275±0.3 1031 DBFGS-0.32-39±29 0.3 14iBFGS hyb 0.33 0.93 79±43 0.3 27表1：性能统计由于平均值超过25，蒙特卡罗法适用于MH的不同提案。5.1具有合成数据的LGSS模型我们首先考虑线性高斯状态空间（LGSS）模型，因为可以使用卡尔曼平滑器精确解决状态推断问题。这使我们能够使用精确递归计算对数后验概率及其梯度，这将指示各种MH方案的最佳性能。模型由xt+1 | xt给出~ Nxt+1；u+φ（xt- u），σv, （17a）yt | xt~ Nyt；xt，0.5, （17b）θ={u，φ，σv}和u∈ R、 φ∈ (-1，1）和σv∈ R+。使用参数{0.2、0.5、1.0}，从模型中模拟了一个由T=500个观测值组成的合成数据集。我们使用25个蒙特卡罗模拟，使用相同的数据计算IF（12），并估计不同方案的计算时间。表1总结了这些模拟的中值结果，其中包括接受率、修正的Hessian估计分数、IF的最大值和每次迭代所需的时间以及从后验获得一个有效样本所需的时间。对于IF，我们还提供IQR（25%和75%分位数之间的距离）。BFGS类型建议书的基准是MH0和MH1的预处理版本，分别表示为pMH0和pMH1。在pMH1中，我们设置H-1（θ）=P，其中P表示使用导频运行计算的后验协方差的估计值。

在pMH0中，我们使用与pmh1相同的方法，但也设置G（θ）=0。使用三种不同的BFGS方案：（d）amped，（i）忽略曲率条件和（e）加强φ后验估计值0.40 0.50 0.60 0.700 2 4 6 8 10 120 50 100 150 2000.40 0.50 0.60 0.70迭代φ0 100 300 5000.0 0 0.2 0.4 0.6 0.8 1.0迭代φ后验估计值0.40 0.50 0.70 0 2 4 8 10 120 50 100 150 2000.40 0.50 0.60 0.70迭代φ0 100 300 5000.0 0 0 0 0 0 0 0.2 0.4 0.6 0.8 1.0迭代φ后验估计值0.40 0.50 0.60 0.7002 4 6 8 10 120 50 100 150 2000.40 0.50 0.60 0.70φ0 100 300 5000.0 0.2 0.4 0.6 0.8 1.0φ的迭代图1：使用pMH0（顶部）、pMH1（中部）和dBFGS（底部）对φ进行后验估计（左侧）、迹线图（中部）和ACF估计（右侧）。左侧和中间图中的虚线表示估计的后验平均值。右图中的虚线表示95%的置信区间。如果左图中的灰色线表示先验分布。曲率条件。对于后两种情况，Hessian估计值通常为负值，因此需要进行一些修正。我们采用附录B中概述的三种不同方法：（FLIP）Ping负特征值，（reg）标定估计值和（hyb）rid方法。每个有效样本的时间（以毫秒为单位）显示在表中最右边的列中。我们注意到，所提出的基于阻尼BFGS的曲率信息局部逼近方法性能最好。从后部产生一个样本需要53毫秒。这比pMH0小，pMH0将是此设置中的标准方法。

此外，我们注意到，其他BFGS类型方案需要对Hessian进行大量修正，这是不可取的，这可能会引入数值不稳定性。图1显示了模拟研究中特定情况下的后验估计、马尔可夫链跟踪及其相应的ACF。我们注意到，与pMH0（上部）和pMH1（中部）相比，BFGS基方案（下部）的混合效果更好。此外，由于马尔可夫链中的M步依赖性，BFGS方案的ACF表现出与其他两个方案相当不同的行为。通过比较后验估计，我们得出结论，所提出的方法能够以正确的参数和合理的方差为中心生成良好的估计。5.2重新审视LGSS模型我们重复相同的实验，如Dahlin等人【2015a】所述，将卡尔曼滤波器和平滑器替换为粒子滤波器和固定滞后粒子平滑器。同样，表1总结了所有计时的结果，现在以秒表示。在这种情况下，pMH0和pMH1的表现更差，因为对数后验值及其梯度的估计值不一致。然而，所提出的方法表现良好，并在14秒后生成一个有效样本。与其他方法相比，这是一个大幅度的减少，尤其是与Zhangand Sutton【2011】和Dahlin等人【2015b】提出的拟牛顿方法相比，至少增加了两倍。因此，所提出的方法优于设计方案的标准方法和当前最先进的拟牛顿方案。5.3 SV模型和比特币价格数据为了证明所提方法的实际应用，我们考虑了估计2015年11月7日至2017年11月7日期间比特币价格波动性的问题。

图2的上半部分显示了对数回报（两天之间比特币价格的百分比变化）为绿点。请注意，日志返回的平均值为零，但方差随时间而变化。We-15-10-5 0 5 10 15时间观测11月15日5月16日11月16日5月17日11-4-2 0 2 4 6时间日志波动性Nov 15 5月16日11月16日5月17日17μ后验估计0.5 1.5 2.50.0 0.5 1.0 1.5 1.5 2.0φ后验估计0.85 0.950 5 10 15 20 25σ后验估计0.3 0.5 0.70 2 4 6 8ρ后验估计-0.3 0.0.30 1 2 4 5 6 7图2：顶部：比特币的日志收益率（绿色）及其估计值95%使用模型和对数波动率估计的置信区间（橙色）。底部是dBFGS获得的u（pruple）、φ（magents）、σv（绿色）和ρ（黄色）的后验估计值。虚线和灰线分别表示估计的后验平均值和先验分布。目的是通过模型捕捉波动性的这种变化（因为这种变化是已知的xt+1年xt公司~ Nxt+1年;u+φ（xt- u),σvρρexp（xt）,这就是所谓的带有杠杆作用的随机波动率（SV）模型。此处，未知潜在参数与LGSS模型中的相关ρ相同∈ (-1，1），即θ={u，φ，σv，ρ}。AIM是在给定数据的情况下估计对数波动率，这可以通过对参数的后验估计进行边缘化来实现，参见Andrieu等人【2010】或Dahlin和Sch¨on【2017】。图2的中间图给出了对数volaility（潜在状态）的结果估计。请注意，比特币价格随时间而变化，当比特币价格波动且每日变化较大时，其变化幅度较大。此外，对数波动率过程的平均值相当大，这意味着对数波动率通常很大（与股票价格相比）。

相关性ρ很可能非常接近于零。这与-0.8，这是股票的典型值，这意味着股票价格的大幅下跌会增加波动性（因为投资者出售资产）。这些信息在许多金融应用中非常有用，例如比特币期货定价，以及计算银行和金融机构向监管机构提交的各种风险度量。6结论数值例子表明，所提出的方法优于许多现有方法，可以为MH创建良好的建议。此外，我们想再次强调，提议的方法基本上不需要试运行，这是所有预处理方法所必需的。此外，阻尼BFGS方法始终提供了Hessian的正有限估计，因此不需要潜在的数值不稳定Hessian校正。此外，梯度信息对于大维参数空间中的估计至关重要，这在SSM和传递函数模型中很常见。所有这些好处都有可能使MH被广泛用于识别动力系统。在本文的范围内，有许多有趣的未来工作途径。SR1更新【Nocedal和Wright，2006年，第6.2章】是BFGS的一种替代方案，已知BFGS在许多情况下可提供更准确的黑森估计值。此外，来自优化的信赖域方法可能有助于MH解决数值稳定性问题。最后，对于参数空间较大的模型，需要进行更广泛的数值计算。

在这种情况下，需要替代或更好的粒子平滑算法来获得合理的梯度估计。本文使用的源代码和数据可从GitHub获得https://github.com/compops/qnmh-sysid2018/并通过Docker（参见README.md）。参考C。Andrieu和J.Thoms。关于自适应MCMC的教程。《统计与计算》，18（4）：343–3732008。C、 Andrieu、A.Doucet和R.Holenstein。粒子马尔可夫链蒙特卡罗方法。皇家统计学会杂志：B辑（统计方法），72（3）：269–3422010。O、 Capp\'e、e.Moulines和T.Ryd\'en。隐马尔可夫模型中的推理。Springer Verlag，2005年。J、 Dahlin和T.B Sch¨on。开始使用粒子Metropolis Hastings进行非线性动力学模型中的推理。《统计软件杂志》，2017年。正在印刷。J、 Dahlin、F.Lindsten和T.B.Schon。粒子Metropolis Hastings使用梯度和Hessianinformation。《统计与计算》，25（1）：81–922015a。J、 Dahlin、F.Lindsten和T.B.Schon。准牛顿粒子Metropolis Hastings。第17届IFAC系统识别（SYSID）研讨会论文集，第981-986页，中国北京，2015年10月b。A、 Doucet和A.Johansen。粒子过滤与平滑教程：十五年后。InD.Crisan和B.Rozovsky，《牛津非线性滤波手册》编辑。牛津大学出版社，2011年。A、 Doucet、M.K.Pitt、G.Deligiannidis和R.Kohn。使用无偏似然估计时，马尔可夫链蒙特卡罗的有效实现。Biometrika，102（2）：295–3132015。J、 Durbin和S.J.Koopman。利用状态空间方法进行时间序列分析。牛津大学出版社，第2版，2012年。M、 Girolami和B.Calderhead。Riemann流形Langevin和Hamilton蒙特卡罗方法。皇家统计学会杂志：B辑（统计方法），73（2）：1–372011年。五十、 永。

系统识别：用户理论。普伦蒂斯·霍尔，1999年。C、 Nemeth、C.Sherlock和P.Fearnhead。粒子Metropolis调整的Langevin算法。Biometrika，103（3）：701–7172016。J、 Nocedal和S.Wright。数值优化。Springer Verlag，2006年第2版。五、 彼得卡。贝叶斯系统识别。Automatica，17（1）：41–531981年。G、 Poyiadjis、A.Doucet和S.S.Singh。状态空间模型中分数和观测信息矩阵的粒子近似及其在参数估计中的应用。Biometrika，98（1）：65–802011。C、 罗伯特。贝叶斯选择。Springer Verlag，2007年。C、 罗伯特和卡塞拉。蒙特卡罗统计方法。Springer Verlag，第2版，2004年。G、 罗伯茨和罗森塔尔。朗之万离散近似的最佳缩放。皇家统计学会杂志：B辑（统计方法），60（1）：255–2681998。T、 B.Schon、A.Wills和B.Ninness。非线性状态空间模型的系统识别。Automatica，47（1）：39–492011年。Y、 张和C.A.萨顿。马尔可夫链蒙特卡罗的拟牛顿方法。2011年12月，西班牙格拉纳达，2011年神经信息处理系统会议记录。A实施细节在第5.1节和第5.2节中，我们使用带有RTS平滑器的标准卡尔曼滤波器来计算对数后验概率及其梯度。此外，使用N=1000个粒子的自举粒子滤波器和具有滞后10的固定滞后平滑器，有关算法描述，请参见Dahlin等人[2015a]。使用Doucet等人【2015年】的结果选择粒子数。为了简单起见，我们在真实参数中初始化所有MH算法，并将其运行K=10000次迭代，然后丢弃前3000次asburn in。MH0和MH1的预处理矩阵是通过多次试运行来计算的。

步长是使用现有的拇指规则选择的【Roberts和Rosenthal，1998年】作为= 1.37和= 使用Kalman方法和= 1.48和= 0.47当使用粒子方法时【Nemethet等人，2016年】。对于qMH方法，我们使用内存长度M=20，并使用随机游走建议，用于具有步长的第一个M次迭代 = 所有三个参数均为0.01。步长= 0.5用于初始M次迭代后的所有qMH算法。对LGSS模型进行重新参数化，以使马尔可夫链中的所有参数不受约束（能够假定任何实值），由φ=tanh（(R)φ），σv=exp（(R)σv）给出，其中θ={u，(R)φ，(R)σv}是马尔可夫链的新状态。变量的这种变化将aJacobian术语引入了接受概率，参见Dahlin和Sch¨on【2017年，第6.3.2节】。最后，我们使用以下先验密度u~ N（0，1），φ~ T不适用(-1,1）（0.5,1），σv~ G（2，2），其中T N（a，b）（·）表示[a，b]上的截断高斯分布，G（a，b）表示平均a/b的Gammadistribution。在第5.3节中，我们使用了与LGSS模型相同的设置，但将N增加到1500。此外，我们将先验值（p（u）保持不变）略微更改为φ~ T不适用(-1,1）（0.95，0.05），σv~ G（2，10），ρ~ T不适用(-1,1）（0,1），更好地反映了现实数据中通常存在的参数值。相关性ρ以与φ相同的方式重新参数化。比特币数据计算为yt=100[对数（st）-日志（st-1） ，其中St表示从https://www.quandl.com/BITSTAMP/USD.BHessian校正第一种方法通过∑？QN公司ψk，M= ∑QNψk，M- 2λminIp，λmindenoting为最小（负）特征值。这使所有特征值都变为正。

第二种方法计算特征值∑？QN公司ψk，M= Q∧| Q-其中Q和∧|分别表示特征向量矩阵和∑QN（ψk，M）特征值绝对值的对角矩阵。第三种方法是Dahlin等人[2015a]的混合方法。其中，估计值被后验协方差∑的全局近似值代替？QN公司ψk，M= Σ-1emp，其中∑emp表示使用老化阶段后半部分计算的后验协方差的样本估计。