状态空间中快速近似推理的贝叶斯优化

我们利用第6.2节中的程序估计每种油的αsv模型参数。假设过滤后的残差（17）是独立的，如果已知x1:t，则其分布为A（bαi，1）乘以（16）。阶段1独立于每个资产进行，因此直接以并行方式实施。1998年2000年2002年2004年2006年2010-20年10月10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日10日。0.2 0.4 0.6 0.8 1.00.0 0.2 0.4 0.6 0.8 1.0布伦特杜拜1998 2000 2002 2004 2006 2008 2010-20-10 0 10 20日期日志返回（迪拜）0 1 2 3 4 5 6 7平滑日志波动性0。0.2 0.4 0.6 0.8 1.00.0 0.2 0.4 0.6 0.8 1.0 Brentmaya1998 2000 2002 2004 2006 2008 2010-20-10 0 10 20日期日志返回（maya）0 1 2 3 4 5 6 7平滑日志波动性0。0.2 0.4 0.6 0.8 1.00.0 0.2 0.4 0.6 0.8 1.0Dubaimaya图5：左图：1997年1月10日至2010年6月4日期间布伦特（绿色）、迪拜（橙色）和玛雅（紫色）石油的每周日志回报。阴影区域表示根据αsv模型，对数收益率的95%置信区间。深灰色曲线表示对数波动率的估计值。右图：相应的转换过滤残差，但i.-20-10 0 10 30 40 50timelog-Returns 1997 2000 2003 2006 2009图6:VaR0的估计值。99（et），用于使用gsv（洋红）和带有Student-t copula的αsv（绿色）模型的三种石油期货的等权投资组合。虚线表示估算和验证数据的划分。与pmh相比，使用该算法将该阶段的计算时间从每个资产的小时或天减少到大约半小时。我们遵循麦克尼尔等人[2010年，第231-231页]对残差中的相关性进行建模。这相当于通过残差的经验cdf将概率变换应用到三维超立方体中，见图5的右部分。

转换后的残差{bu1:T，1，…，bu1:T，d}然后通过Student的st-copula函数组合，以找到联合分布的模型。学生t-copula中的自由度和相关矩阵分别使用映射和通过肯德尔τ的矩匹配来估计。最后，我们利用copula模型及其边际来估计每种石油的var。信用风险值α∈ （0，1）由Var′α（et）=infn定义- et∈ R:G(-et）≥ αo，即最小的损失-损失（负对数回报）超过-etisno大于（1- α). 我们采用蒙特卡罗方法来估计VaRα（et），方法是：（i）从copula模拟，（ii）通过基于经验cdf的分位数变换获得模拟的过滤残差，（iii）通过估计的波动率和（17）的倒数计算得出的对数收益。我们对每种资产重复10万次，然后计算对应于VaR0的经验分位数。99（et），用于同等权重的投资组合。然而，更高级的投资组合权重方案可以很容易地实现。由此产生的var估计如图6所示，其中我们将αsv模型与使用gs以类似方式估计的gsvmodel（15）进行了比较。我们注意到，两个模型的估计值相差很大，尤其是在数据系列的末尾。对验证数据的回溯测试给出了两个模型的0个违规，预计违规数量为2.3。gsa算法需要-20 0 40 TimeLog-returns 1926 1936 1946 1956 1976 1986 1996 2006图7：日志返回（灰点）和VaR0的估计值。99（et），用于使用gsv（品红）和αsv（绿色）模型以及Student t-copula的等权重股票组合。

灰线表示估算和验证数据的划分。30岁左右- 60分钟来推断三种资产中每种资产的模型，这种推断可以直接在并行架构上运行。pmhalgorithm所需的相应计算时间约为18- 24小时。6.4通过一个最终的数值例子计算股票组合的VAR，以说明所提出的方法也可以应用于大型投资组合。我们考虑的数据是Kenneth French提供的30个行业组合。投资组合包括1926年9月至2015年12月期间30个不同工业部门的月度日志报告。同样，我们将数据划分为一个估计集和一个验证集，分别有716个和358个数据点。我们采用与第6.3节相同的程序，结果如图7所示。结论是相似的，两种var估计都非常相似。对于αsv和gsv模型，回溯测试分别给出0和1个违规。预计违规次数为3.5次。与pmh相比，gsa的计算速度与第6.3节中的相似，因为观测次数相似。7结论我们提出了gsa算法，这是一种近似具有难治性似然性的ssms后验分布的算法。第6节中提供的插图表明，所提出的算法相当准确，并且显示出计算成本的大幅降低。我们以一到两个数量级的速度获得了类似的后验估计，这将计算时间从几十小时减少到几十分钟。此外，gsa似乎对对数后验估计中的噪声非常鲁棒，这通常会导致pmh算法有时陷入停滞，而spsa算法收敛缓慢。

总的来说，这表明所提出的算法是一种高效的推理方法，使得实践者可以在应用工作中使用具有难以解决的可能性的模型，例如具有α稳定裕度的copula模型。未来的工作包括：（i）采用gp的稀疏表示；（ii）开发新的量化函数；（iii）利用更好的定制协方差函数；（iv）引入自适应方法来选择公差水平. Huber[2014]和Bijl等人[2015]讨论了GPS序列和稀疏表示的一些想法。设计新的采集函数，以便在map估计的同时获得对数后验矩或高阶矩的Hessian估计值，这也很有趣。这可以改进参数后验的拉普拉斯近似，并为其他后验近似打开大门，例如使用倾斜的Student t分布。此外，定制gp先验的一种有趣方法是利用非平稳协方差函数[Paciorek和Schervish，2004]。这将捕捉到这样一个事实，即对数后验值在参数空间的某些部分通常会迅速下降，而在其他部分几乎会下降。最后，加入一些适应性的方法来进行选择是有益的 在smc abc算法中，减少用户的选择数量。如前所述，我们最初实施了Del Moral等人[2012]和Calvet and Czellar[2015]提出的自适应算法。然而，当使用中等N值来近似对数后验值时，我们遇到了一些问题。估计值表现出很大的偏差，导致使用pmh对θ进行有偏差的估计。因此，我们选择不去适应 能够保持N小得多。

因此，探索使用中等N值对对数后验值进行良好估计的其他自适应模式是很有趣的。确认这项工作得到了以下支持：动力系统概率建模（合同编号：621-20135524）和林奈中心cadics，两者均由瑞典研究委员会资助。模拟是在瑞典国家计算基础设施（snic）atLink–oping大学提供的资源上进行的。J.Dahlin感谢Fredrik Lindsten、Neil Lawrence和CarlHenrik Ek进行了富有成效的讨论。还要感谢Joerg M.Gablonsky、Abraham Lee、Per A.Brodtkorband和GPy团队提供了他们的Python实现。实现细节GPO算法：我们利用GPy软件包[GPy作者，2014]计算GPy预测后验概率并估计gp先验超参数。在本文中，我们假设一个零先验平均函数m（θ）=0，以及偏差和Mat’ern 5/2协方差函数的组合。这种选择对应于对数后验的先验，具有一些非零均值和两个连续导数。这些都是合理的假设，因为这种平滑度是在拉普拉斯近似中假设的。我们使用经验贝叶斯（eb）估计超参数，即通过优化关于λ的边缘似然。在sga中也可以使用使用切片采样[Murray et al.，2010]或smc[Svensson et al.，2015]的超参数边缘化的更先进方案。对于（14）中的采集规则，我们使用ζ=0.01和∑=0.01Ip。我们使用L=50个样本初始化gpo算法，这些样本是使用拉丁超立方体采样获得的，AbrahamLee编写的实现可在https://pypi.python.org/pypi/pyDOE.

使用Joerg M.Gablonsky编写的直接实现解决了算法2中第8行和第11行的优化问题，可从https://pypi.python.org/pypi/DIRECT/.最后，对于第12行，我们使用了Per A.Brodtkorb提供的Pythonimplementation，网址为https://pypi.python.org/pypi/Numdifftools.Section6.1：我们在gs中使用N=2000个粒子，并且N=2000个粒子具有高斯密度和公差水平 = 0.20和ψ（x）=gsa中的x，得出图2中的结果。在初始化之后，我们对K=450次迭代运行GPOALGolmethods，并在每25次迭代之前重新估计GPOALGolmethods的超参数。gpo算法ΘGPOis的搜索空间由u给出∈ (0, 1), φ ∈ （0,1）和σv∈ (0.01, 1). 我们使用以下先前的密度p（u）~ N（u；0,0.2），p（φ）~ T N(-1,1）（φ；0.9,0.05），p（σv）~ G（σv；2，20），其中tn（a，b）（·）表示[a，b]上的截断高斯分布，G（a，b）表示平均a/b的伽马分布。对于pmh，我们使用N=2000个粒子的smc abc算法来估计对数后验分布。我们在θ={0.10,0.95,0.12}中初始化pmh，并运行M=15000次迭代的算法（将前5000次作为老化丢弃）。参数建议选择为q（θ|θk）=N（θ；θk，∑q），∑q=2.562·10-4.diag（137,7,38），由渐近经验法则得出，参见Dahlin和Sch¨on[2015]，并通过试验运行对后协方差进行了估计。对于spsa，我们使用N=2000个粒子，并遵循Spall[1998]的规定，通过试验运行选择超参数a=0.001、c=0.30、a=35、α=0.602和γ=0.101。第6.2节：实际数据计算为yt=100[log（st）- 原木（圣-1） ]中，其中St表示未来合同的价格。我们遵循Yildirim等人[2014]的观点，应用ψ（x）=arctan（x）来稳定似然方差（和梯度估计）。

两步方法适用于以下数据：https://www.quandl.com/CHRIS/ICE_KC2.sample从零均值对称α-稳定分布A（α，γ）。首先，我们采样v（1）t~ Exp（1）和v（2）t~ U(-π/2, π/2). 然后，我们通过应用转换ˇyt=γsin获得一个样本（当α6=1时）αv（2）tcos（v（2）t）1/α“cos(α - 1） v（2）tv（1）t#1-αα.有关生成α稳定随机数的更多信息，请参见Nolan[2003]。我们使用N=2000个粒子，其高斯密度具有容差水平 = smc abcto中的0.10估计对数后验概率。我们使用与之前相同的设置运行gpo算法，但添加了α∈ （1.2,2）到搜索空间和p（α）~ B（α/2；20，2）表示先验分布，其中B（a，B）表示β分布。我们在θ={0.22,0.93,0.25,1.55}中初始化pmh，并使用试验运行asq（θ|θk）=N（θ；θk，∑q），∑q=2.562·10选择参数建议-3.诊断（26,1,9,11）。第6.3节和第6.4节：石油和股票投资组合示例的大多数设置相同。我们在smc和smc abc算法中使用N=5000个粒子，并保留其余设置Assets 6.1和6.2。对于股票投资组合示例，我们将u的搜索空间更改为（0,4），因为此数据集中的每周日志回报可能比其他数据集中的每日日志回报大得多。copula中自由度的映射估计是通过使用auniform先验的拟牛顿解算器获得的。参考资料c。安德烈、A·杜塞特和R·霍伦斯坦。粒子马尔可夫链蒙特卡罗方法。theRoyal统计学会杂志：B辑（统计方法），72（3）：269-34220010。H.Bijl、J-W.van Wingerden、T.B.Sch¨on和M.Verhaegen。使用FITC和PITC近似的在线稀疏高斯过程回归。《第17届IFAC系统识别研讨会论文集》（SYSID），第703-708页，中国北京，2015年10月。E.布罗库、V.M.科拉和N。

德弗雷塔斯。关于昂贵代价函数的贝叶斯优化的教程，应用于主动用户建模和分层强化学习。印前，2010年。arXiv:1012.2599v1。A.D.公牛。高效全局优化算法的收敛速度。机器学习研究杂志，12:2879–29042011。可从以下网址获取数据：http://freakonometrics.free.fr/oil.xls.Data网址：http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html.L.E.卡尔维特和V.泽拉。近似贝叶斯计算过滤的精确方法。《金融计量经济学杂志》，13（4）：798-838，2015年。答：查彭蒂埃。公共财政政策。技术报告，2015年5月。统一资源定位地址https://hal.archives-ouvertes.fr/hal-01151233.J.Dahlin和F.Lindsten。用于参数推断的基于粒子滤波的高斯过程优化。2014年8月在南非开普敦举行的第19届IFAC世界大会上发表。J·达林和T·B·施农。开始使用粒子Metropolis Hastings在非线性模型中进行推理。印前：，2015年。arXiv:1511.01707v4。T.A.迪恩和S.S.辛格。近似贝叶斯估计的渐近行为。印前：，2011年。arXiv:1105.3655v1。P.德尔·莫勒尔、A.杜塞特和A.贾斯拉。一种用于近似贝叶斯计算的自适应序贯蒙特卡罗方法。《统计与计算》，22（5）：1009-102012。A.杜塞特和A.约翰森。粒子过滤与平滑教程：十五年后。InD.Crisan和B.Rozovsky主编，《牛津非线性滤波手册》。牛津大学出版社，2011年。J·德宾和S·J·库普曼。用状态空间方法进行时间序列分析。牛津大学出版社，第2版，2012年。E.埃利希、A.贾斯拉和N.坎塔斯。具有难处理似然性的隐马尔可夫模型的无梯度参数估计。《应用概率中的方法和计算》，第17（2）：315–3492015页。麻省理工大学。

古特曼和J·科兰德。基于模拟器的统计模型无似然推理的贝叶斯优化。机器学习研究杂志，17（125）：1-472016。M.F.胡伯。递归高斯过程：在线回归和学习。《模式识别快报》，2014年第45期，第85-91页。A.贾斯拉、S.S.辛格、J.S.马丁和E.麦考伊。通过近似贝叶斯计算进行过滤。《统计与计算》，22（6）：1223-12372012。H.乔。基于copula模型的两阶段估计方法的渐近有效性。多变量分析杂志，94（2）：401–4192005。D.R.琼斯、C.D.佩特纳和B.E.斯图克曼。没有Lipschitz常数的Lipschitz优化。优化理论与应用杂志，79（1）：157–181，1993年。D·J·利佐特。实用的贝叶斯优化。艾伯塔大学博士论文，2008年。L·容格。系统识别：面向用户的理论。普伦蒂斯·霍尔，1999年。J-M.马林、P.普洛、C.P.罗伯特和R.J.赖德。近似贝叶斯计算方法。《统计与计算》，22（6）：1167-11801212。A.J.麦克尼尔、R.弗雷和P.恩布雷切斯。定量风险管理：概念、技术和工具。普林斯顿大学出版社，2010年。E.米兹和M.韦林。GPS-ABC：高斯过程替代近似贝叶斯计算。2014年7月，加拿大魁北克市，第30届艺术情报不确定性会议记录。I.默里、R.亚当斯和D.麦凯。椭圆切片取样。2010年5月，意大利撒丁岛，第13届国际艺术情报与统计会议（AISTATS）论文集，第541-548页。J.诺兰。稳定分布：重尾数据的模型。Birkhauser，2003年。C·J·帕西奥里克和M·J·舍尔维什。

高斯过程回归的非平稳协方差函数。《2004年神经信息处理系统（NIPS）会议记录》，第273-280页，加拿大温哥华，2004年12月。帕诺夫和斯波科尼。半参数问题的有限样本Bernstein von Mises定理。贝叶斯分析，10（3）：665-7102015。M.K.皮特、R.S.席尔瓦、P.乔达尼和R.科恩。基于粒子滤波的马尔可夫链蒙特卡罗模拟方法的一些性质。计量经济学杂志，171（2）：134-1512012。C·E·拉斯穆森和C·K·I·威廉姆斯。机器学习的高斯过程。麻省理工学院出版社，2006年。J·C·斯帕尔。随机优化的同步摄动算法的实现。IEEE航空航天和电子系统学报，34（3）：817–823，1998年。S.V.斯托扬诺夫、B.Racheva Iotova、S.T.拉切夫和F.J.法博齐。波动市场风险估计的随机模型：一项调查。运筹学年鉴，176（1）：293-3092010。A.斯文森、J.达林和T.B.肖恩。用序贯蒙特卡罗法边缘化高斯过程超参数。2015年12月，在墨西哥坎昆举行的第六届IEEE国际多传感器自适应处理（CAMSAP）计算进展研讨会上发表。GPy的作者们。GPy：Python中的高斯过程框架。http://github.com/SheffieldML/GPy, 2014.E.巴斯克斯和J.贝特。具有固定均值和协方差函数的预期改进算法的收敛性。《统计规划与推理杂志》，140（11）：3088–30952010。S.N.伍德。噪声非线性生态动力系统的统计推断。《自然快报》，466:1102–11042010。S.Yildirim、S.S.Singh、T.Dean和A.Jasra。用序贯蒙特卡罗法在具有难处理可能性的隐马尔可夫模型中进行参数估计。《计算与图形统计杂志》，24（3）：846–8652014。