非平稳反射布朗运动的精确模拟

因此，问题归结为生成（η）的精确样本, sup0≤T≤Z（t）给定() = x、 下一个定理提供了过程定律ZBr，以及astandard Brownian bridge，它被用作验收/拒收方案中的验收概率。定理6.1。假设γ（s）是一个连续可微函数，设τa=inf{t≥0:| Z（t）- Z（0）|≥ a} 。然后，对于一些常数c，我们得到（Z（t））0≤T≤∧τa∈ dωZ() = 十、Q（Wt）0≤T≤∧τa∈ dωW= 十、= cP（I=1 |ω）×eψ（ω），（6.1），其中ψ（ω）=如果τa为0≥ 2(-τa）（十）- §a）- （十）- ~a-Rτaγ（u）du）如果τa<,§a=Z（τa）=±a，P（I=1 |ω）在（5.3）中定义。证据在附录中。观察ψ（ω）≤ ~m/2（x+a+m）),式中，m=sup0≤T≤γ（t）。因此，很容易想象构造一个成功概率与exp（ψ（ω））成正比的伯努利随机变量。通过修改子程序1中的步骤2、3和4，可以构造指示器I。主要区别在于以W为条件的采样（τa，Wκ，Wκ，…，Wκb）= x、 κ，κ，κb<τa∧  是速率为2ma的泊松过程的跳跃时间。该程序包括：样品（τa∧ , Wτa∧) 鉴于W= x、 样本κ，κ，κb<τa∧ , 这是速率为2mA的泊松过程的跳跃时间。考虑两种不同的情况：-如果τa≤ , 样本Wκ，Wκ，Wκb条件on（τa，Wτa）。-如果τa≤ , 样本Wκ，Wκ，Wκb在W上的条件= x和τa≤ .生成（η）样本的程序的其他步骤, m()) 与子程序1中的步骤5和6完全相同。第2步是可行的，通过使用定理6.2，我们计算布朗桥和标准布朗运动的第一次击中时间分布之间的似然比。定理6.2。设τa=inf{t≥ 0:| Wt |≥ a} 成为标准布朗运动的第一次击球时间。

假设r=|x | a<1，那么我们有q（τa）≥ |W= 十）=1.- E-2(1-r）· p（0，1，, 1.- r） ，其中p（s，x，t，y）在（C.1）中定义。此外，对于某些常数c，我们有q（τa）∈ dt，Wτa=~a | W= x） Q（τa）∈ dt，Wτa=~a）=cg十、- ~a√ - T≤C√2πe·| a- |x| |（6.3）每t, 式中，g（·）是高斯函数，a=±a。函数p（s，x，t，y）的定义和定理的证明见附录。这个定理有助于上述步骤中的第2步。如果| x |>a，很明显τa<. 假设| x |<a.对指示剂I（τa）进行采样≥ ) 可以通过生成两个独立的均匀随机变量U和V来实现。IfU<p（0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,，, 1.- r） 和V<1.- E-2(1-r）< 1为真，我们设置τa≥ , 确定U是否<p，, 1.- Chen和Huang（2013）中的命题4.1在有限时间内实现了2r）。定理6.2的第二部分有助于生成（τa，Wτa）的样本，该样本以τa< 而W= x、 我们可以使用接受/拒绝方法。按照Burqand-Jones（2008）提出的方法，可以生成标准布朗运动的（τa，Wτa）样本。根据方程式（6.3），该样本可能被接受，概率与g成正比十、-~a√-T作为（τa，Wτa）的样本，假设τa< 而W= x、 上述程序的第2步很清楚。此外，如果τa≤ , 生成Wκ，Wκ，Wκb，Wτa通过Chen和Huang（2013）中的程序（17）可以对τais进行补充。现在，如果条件（5.4）成立，我们可以指定I=1。在本部分的其余部分中，我们将详细介绍如何在一系列实例κ、κ、。

，κb<τa∧  以W为条件= x和τ在这种情况下 < τa.LetfWt=1.-aWτa-taif Wτa=a1+aWτa-taif Wτa=-a、 由于布朗运动的自相似性和时间反转特性，很容易验证fw=（fWt:0）≤ T≤τaa）是一个布朗过程，假设≤fWt≤ 2和FWTB+1=1-xa。因此，我们感兴趣的是在instancest=τaa的序列上生成样本OFW≥ T≥ . . . , ≥ 肺结核≥ tb+1，其中ti=τa-κIa对于i=1，b、 设V=（Vt:tb+1）≤ T≤ t） 假设Vtb+1=1，则为布朗曲线-xa和Vt=1。Devroye（2010年，第6节）讨论了（Vt，…，Vtb）的精确采样。我们可以生成随机向量（Vt，…，Vtb）的样本，并将其接受为（fWt，…，fWtb）的样本，概率按命题6.1计算。验收决定可通过Chen和Huang（2013年，第4.4节）中提出的方法确定。提议6.1。无论如何≥ T≥ . . . , ≥ 肺结核≥ tb+1和0≤ Yyb≤ （fWt，…，fWtb）的联合条件分布关于（Vt，…，Vtb）有以下似然比：QfWt∈ 嗯，fWtb∈ 戴布0≤fWt≤ 2，fWtb+1=yb+1，fWt=1Q及物动词∈ 嗯，Vtb∈ 戴布0≤ Vt，Vtb+1=yb+1，Vt=1= ~cbYi=0p（ti+1，yi+1；ti，yi），其中yb+1=1-xa，y=1，且c>0是一个常数。这个命题的证明类似于Chen和Huang（2013）中的定理4.2。给定骨架（fWκ，…，fWκb，W), 采样sup0≤T≤WT和maximumtime的位置类似于子程序1的步骤6。

联合采样（ui，tui）是有效的，其中ui=supκi≤T≤κi+1Wt和tui是[κi，κi+1]上最大值的位置，条件是τaas为布朗曲流的最大值；参见Devroye（2010）中的maxmeander算法。7算法2的分析在本节中，我们分析了算法2，并表明算法以多项式时间终止。该算法生成（vt，M（t））的精确样本的运行时间最多为O（1/？）γ，与t无关。因此，生成三元组（vt，M（t），Y（t））的精确样本的算法的运行时间最多为O（1/？）γ+O（log（t）），其中O（log（t））是读取输入和计算Y（t）的运行时间，条件是（vt，M（t））。这里的关键思想是构造一个恒定的漂移布朗运动，控制Z（s）。引理7.1。设mk为直到时间βk的过程Z（t）的最大值，其中βk等于Z（βk）<ζk-光盘然后，Z（t）<mk表示所有t<αk+1。具体来说，如果αk+1=∞, 那么mk是所有t的Z（t）的全局最大值≥ 0.证明。根据假设（3.2），我们有Z（t）=Z（βk）+B（t）- B（βk）+Ztβkγ（u）du≤ Z（βk）+B（t）- B（βk）+d- （t）- βk）·γ=Uβkt。观察Uβkβk=Z（βk）+d≤ mk- （c）- 1） d.回想一下αk+1=inf{t≥ βk:Ut≥ mk}。因此，对于所有的t<αk+1，我们有Z（t）<mk。通过子程序1可以生成βk+1和mk+1的精确样本。注意，过程uβkis是一个具有恒定漂移的布朗运动。因此，我们可以很容易地根据以下引理对命中时间αk+1进行采样。引理7.2。设x=mk- Uβkβk>0。那么，P（αk+1=∞|Uβk）=1- 经验(-2′γx）和pαk+1≤ t+βk | Uβkβk，αk+1<∞= P搞笑x |γ|，x≤ T, （7.1）式中，IG（·，·）表示具有均值x |γ|和形状参数x的逆高斯分布。备注7.3。从逆高斯分布中提取的数据可以非常高效地生成；见Devroye（1986年，第四章）中描述的算法。证据

根据假设（3.2），\'γ>0，我们从Karatzas和Shreve（1991，第297页）得到αk+1≥ T+βk | Uβk= Φx+’γT√T- 经验(-2′γx）Φ-x+’γT√T.因此，P（αk+1=∞|Uβk）=1- 经验(-2′γx）和pαk+1≤ T+βk | Uβk，αk+1<∞= 经验(-2′γx）Φ-x+’γT√T+ Φ-十、- γT√T,这是逆高斯分布的分布。由于这两个引理，我们可以很容易地观察到算法终止于一个有限的迭代。算法在迭代k中终止的概率（即αk+1=∞)至少是1- 经验(-2′γ（c）- 1） d）。此外，如果αk+1=∞, 那么ζkis是所有t的Z（t）的最大值≥ 因此，在每个步骤k，程序以至少恒定的概率终止。因此，算法几乎肯定会在有限的时间内终止。我们在下面的定理中总结了这个结果。定理7.1。对于有限时间内的每2次，γ终止。此外，预期的迭代次数最多为（1）-经验(-2（c）- 1） d′γ）-生成（vt，M（t））精确样本的预期运行时间为每t的O（1/？）γ∈ [0, ∞].证据很明显，生成（vt，M（t））的精确样本的运行时间对于每个t是有界的≥ 0以生成（v）精确样本的运行时间为界∞, M(∞)).因此，我们展示了最近案例的结果。观察uβkβk=Zβk+d≤ mk- （c）- 1） 因此，根据引理7.2，算法在下一步终止的概率αk+1=∞|Uβk= 1.- 实验（2）（mk）- Uβk）’γ）≥ 1.- exp（2（c）- 1） d′γ）。因此，由K表示的终止前的迭代次数由一个几何变量控制。所以我们有[K]≤ (1 - 经验(-2（c）- 1） d′γ）-1.假设αK=∞, 我们有≤ Ut<Mk，适用于所有t>βK。

由于mk是时间βk之前的过程的最大值，我们得出结论≥0Z（t）=mK。根据定理5.2，生成ζk=supαk样本的预期运行时间≤子程序1得出的t<βkZ（t）每1个为O（1/|γ|）≤ K≤ K.由于E[K]=O（1/|γ），算法2的预期运行时间是O1/γ.最后，值得一提的是，可以通过改变βk:βk=inf{t的更新规则来提高算法的性能≥  + αk:Z（t）≤ MK- cd}，（7.2）其中mk=支持≤+αk-1Z（t）。在这里 是一个参数，用于确保每个步骤不会太短。仿真实验表明，选择合理 提高了运行时间。8.数值实验我们通过数值实验说明了精确采样方法的有效性和相对性能。我们将精确算法2应用于漂移系数γ（u）=cos（2πu）的RBM- 0.5.换句话说，dXt=（cos（2πt）- 0.5）dt+dB（t）+dLt。（8.1）漂移系数γ（u）=cos（2πu）- 0.5是周期为1的周期函数，且¨γ=-Zcos（2πu）- 0.5du=0.5>0。因此，假设3.2成立。回想一下thatX（n）=> M(∞) = 监督≥0Zt（cos（2πu）- .5） du+B（t）. （8.2）在本实验中，我们比较了离散化方法和生成M样本的精确算法(∞). 传统的离散化技术只能逼近M的样本(∞);精确算法返回M的精确样本(∞).8.1离散化方法我们将我们的算法的精确绘图与简单离散化方案的近似绘图进行比较。离散化（8.1）的一种简单方法是：Xti+1=max（Xti+γ（ti）（ti+1- ti+B（ti+1）- B（ti），0），其中i=1，2。步长δ>0。然而，Asmussen等人（1995年）表明，这种离散化方案有很大的偏差，偏差至少为δ1/2级。

作为替代方案，我们对依赖时间的RBM采用类似于L’epingle（1995）的离散化方案。在Proposition A.1中，我们证明了该方案的偏差为δ阶。在这里，我们快速地概述了隔离方案。对于时间步长δ>0，定义分段常数漂移布朗运动^Zδ={^Zδ（t）}t≥0.考虑离散网格点ti=δi，i=0，1，2。。设^Zδ（t）=0，且forti≤ T≤ ti+1^Zδ（t），^Zδ（ti）+γi（t- ti）+B（t）- B（ti），（8.3），其中γi=δRti+1tiγ（u）du。联合分发^Zδ（ti+1），supti≤T≤ti+1^Zδ（t）给定^Yδ（ti）已知。因此，过程的精确采样及其在网格点上的最大值是可能的。Asmussen和Glynn（2007年，第302页）给出了该算法的详细信息。通过选择一个大的T=T和足够小的δ，随机变量sup0≤T≤tM^Zδ（t）近似(∞).在命题A.1中，我们讨论了如何在时间步长δ和试验次数之间分配离散化方法的计算预算。我们表明，对于离散化的Firstorder方法，将时间步长的数量增加到复制次数的第四根是渐近最优的。然而，最优的比例常数是未知的。8.2结果我们生成了M(∞) 使用精确算法2和离散化方案在（8.2）中定义。对于离散化方案，我们使用不同的增量δ和固定时间范围T=35。表1显示了从M获得200000张图纸所需的时间(∞) 对于不同增量δ的精确算法和离散化方案。此外，还包括Kolmogorov–Smirnov检验的p值，该检验将不同增量δ的Euler方案的近似样本与精确样本进行比较。

p值表示样本是否来自同一分布。表8.2报告了E[M]估计值的比较(∞)] 用这两种方法。根据定理A.1，步长设置为δ=N-, 其中N是模拟试验的数量。表中第四列和第五列显示了E[M]的估计值(∞] 以及不同数量试验的90%置信区间。标准误差（SE）估计为模拟输出的样本标准偏差除以试验次数的平方根。偏差由估计值的期望值与E[M]的真实值之间的差异给出(∞)].精确方法产生的估计偏差为零。因此，使用精确方法生成的200万次试验估计了真实值。离散化方案的偏差通过使用800000条轨迹进行估计。表的第8列报告了PSE+偏差计算的均方根误差（RMSE）。最后一列显示了生成不同数量试验所需的计算时间。在一台配备英特尔酷睿Duo 3.16GHz处理器和4GB内存的服务器上进行了模拟。代码是用MATLAB版本（R2011b）编写的。值得注意的是，精确算法比离散化方法有效得多。偏差为零，误差较小，甚至算法也提高了运行时间。图3比较了精确和离散化方法的收敛速度。该方法具有最优的收敛速度；估计量的RMSE以O（1）的速率减小/√t） ，其中t是计算预算。离散化方案的收敛速度为O（t）-2/5），确认定理A.1。

由于离散化方案的收敛速度较慢，我们可以得出结论，离散化偏差是显著的。表1：通过Kolmogorov–Smirnov检验比较离散化方案和精确方法生成的样本。δp值时间（秒）离散化-1E-152 21-2E-20 29-四十点零零七六二-六十点零零二二二一-一百点一四五二三零五-120.231 11200精确方法——1305 E[M]的模拟结果(∞] 根据模型（8.1）。Kolmogorov–Smirnov检验的p值，带有零假设，即精确和相应的近似绘制来自相同的分布。参考Asmussen，S.（1992年）。重型运输中的排队模拟。运筹学数学，17（1）：84–111。Asmussen，S.，Glynn，P.，和Pitman，J.（1995年）。一维反射布朗运动模拟中的离散化误差。《应用概率年鉴》，第875-896页。Asmussen，S.和Glynn，P.W.（2007）。随机模拟：算法与分析。斯普林格，纽约。Asmussen，S.和Thrisson，H.（1987年）。周期性队列的马尔可夫链方法。应用概率杂志，24（1）：第215-225页。Bambos，N.和Walrand，J.（1989年）。在具有周期性输入的队列上。应用概率杂志，26（2）：381-389页。巴斯，R.F.（2011）。随机过程，第33卷。剑桥大学出版社。A.贝斯科、O.帕帕斯皮利奥普洛斯和G.O.罗伯茨（2006）。扩散样本路径的回顾性精确模拟及其应用。伯努利，12（6）：1077-1098。Beskos，A.和Roberts，G.O.（2005）。精确模拟差异。《应用可能性年鉴》，15（4）：2422-2444。Burq，Z.A.和Jones，O.D.（2008）。第一次通过时布朗运动的模拟。《模拟中的数学与计算机》，77（1）：64-71。陈北和黄，Z.（2013）。布朗运动驱动的随机微分方程的局部化和精确模拟。运筹学数学，38（3）：591-616。德夫罗耶，L。

(1986). 非均匀随机变量生成。斯普林格，纽约。Devroye，L.（2010年）。关于与布朗运动和布朗弯曲有关的一些分布的精确模拟算法。《应用概率与统计学的最新发展》，第1-35页。斯普林格。吉塞克，K.和斯梅洛夫，D.（2013）。跳跃差异的精确采样。运筹学，61（4）：894-907。哈里森，J.M.（1985）。布朗运动和随机流动系统。威利，纽约。哈里森，J.M.和莱蒙，A.J.（1977）。周期队列的极限定理。《应用可能性杂志》，14（3）：第566-576页。海曼，D.P.和惠特，W.（1984）。具有时变到达率的队列的渐近行为。应用概率杂志，21（1）：第143-156页。伊格哈特，D.L.和惠特，W.（1970）。重传输中的多信道队列提高了不适用概率，2（1）：150–177。I.Karatzas和Shreve，S.（1991年）。布朗运动与随机微积分，第113卷。斯普林伯格，纽约。L’epingle，D.（1995年）。反射随机微分方程的Euler格式。《模拟中的数学与计算机》，38（1）：119-126。梅西，W.（1981）。非固定队列。斯坦福大学博士论文。惠特·W.（1989）。规划排队模拟。管理科学，35（11）：1341-1366。离散化方案的最佳收敛速度本附录给出了通过第8小节讨论的离散化方案近似生成样本时，在选择试验次数N和时间步长δ之间的最佳权衡。1.假设我们想要计算α=E[M（T）]，其中M（T）=sup0≤T≤TB（t）+Ztγ（u）duby-montecarlo模拟。然而，我们生成了^Mδ=sup0的样本≤T≤T^Zδ（T），其中^Zδ（T）由（8.3）使用长度δ的时间步长定义。作为α的估计量，我们计算^α（δ，N）=NNXj=1^Mδj，其中^Mδ，^Mδ是^Mδ的i.i.d.拷贝序列。