高阶格式的递归边际量子化

然而，量化器在与每个码字相关联的Voronoi区域上累积概率质量，并且由于这些区域在尾部大于分布体中，因此在尾部中的码字的累积概率质量成比例地大于在体中的码字的累积概率质量。2.2.2 N个单中心卡方分布通常，非中心卡方分布必须使用Besselffunction来指定，但当自由度等于1时，情况并非如此。特别是，考虑随机变量X=（Z+u），其中Z~ N（0，1）。然后是X~ χ′2（1，λ）是一个非中心的i平方分布，具有一个自由度，n关于中心参数λ=u。此外，在x上∈ R+，我们有Fx（x）=√x个φ（x+）+φ（x-)FX（x）=Φ（x+）- Φ（x-)MX（x）=（1+λ）Φ（x+）- Φ（x-)+ φ（x+）x-- φ（x-)x+，其中φ（·）和Φ（·）分别是标准正常PDF和CDF，x±=±√x个-√λ。这意味着我们可以使用标准的正态PDF和CDF表示具有一个自由度的非中心卡方分布，从而允许高效地计算量化方案。这对于我们在本文后面实现高阶RMQ方案时的计算效率非常重要。在实现VQ时，在左、右极限下评估这些函数时必须小心。为确保收敛，我们将fX（0）=fX（0）=MX（0）=0，fX(∞) = 0，外汇(∞) = 1和MX(∞) = 1+λ。当然，当x为负时，这三个函数都为零。γn给出了Γ（0）的良好初始猜测=（3）+√λ） nN型对于√λ<2.55nN+1- 2.5+√λ对于√λ≥ 2.5，对于1≤ n≤ N、 图4显示了基数N=50的三个量化器示例，用于一系列非中心参数范围内具有一个自由度的非中心chisquared分布。请注意，对于非中心性参数的某些值（例如。

中心点的分布是双峰的，而对于较大的值，它类似于高斯分布。3递归边际量化考虑随机微分方程dxt=a（Xt）dt+b（Xt）dWt，X=X规定的连续时间差∈ R、 （3）定义过滤概率空间(Ohm, F，（Ft）t∈[0，T]，P），具有a和b充分光滑的正则函数，以确保弱解的存在。有趣的问题是：如何最佳地近似Xtk：Ohm → R、 一段时间内离散点tk∈ [0，T]，当Xtkis的分布未知时？通常，这是通过使用SDE的离散时间近似方案执行蒙特卡罗实验来实现的，最简单的方案是Euler Maruyama【Maruyama，1955年】updateeXk+1=eXk+a（eXk）t+b（eXk）√tZk+1=：U（eXk，Zk+1），0 5 10 15 20 25 3000.51密度λ=100.10.2量化器0 5 10 15 20 25 30 35 4000.050.10.15密度λ=500.020.040.06量化器0 20 40 60 80 100 12000.02密度λ=5000.020.04量化器图4：具有不同非中心参数一个自由度的非中心卡方分布的三个示例。对于0≤ k<k，其中t=t/K，独立Zk+1~ N（0，1），初始值ex=x。Pag\'es和Sagna【2015】的创新之处在于表明，基于量化这些更新的递归过程是可能的。由于是高斯分布，因此可以使用矢量量化来获得Γ，这是上述方案第一步的最佳量化网格。然而，必须找到一种方法来量化EXK+1的递增（边际）分布。已知exk的分布，量化器Γk+1的失真可以写为Γk+1= 呃eXk+1-bXk+1i=EhEheXk+1-bXk+1eXkii=EhEhU（eXk，Zk+1）-bXk+1eXkii=ZREhU（x，Zk+1）-bXk+1idP（eXk≤ x） 。（4） 不幸的是，我们不知道k＞1时exk的确切分布。

Pag\'es和Sagna【2015】的主要结果表明，如果使用先前量化的Bxk分布，而不是连续的Bxk分布，则结果过程会发生变化。此外，与此过程相关的误差由一个常数限定，该常数取决于使用的参数。结果，可以将（4）中的积分重写为码字sin量化器Γ及其相关概率上的s um。然后，可以计算畸变的近似值Γk+1≈ DΓk+1:=NkXi=1EhU（γik，Zk+1）-bXk+1iPbXk=γik.这里，NK是量化器Γkat时间步长k的基数，允许变化。根据D的定义Γk+1我们现在可以指定牛顿-拉斐逊迭代，以便在时间步长k+1计算量化器，从而使失真最小化。给定时间tk处的量化器，表示为列向量Γk，以及相关概率，PbXk=γik对于1≤ 我≤ Nk，量化器Γk+1在时间tk+1的牛顿-拉斐逊迭代由Γ（n+1）k+1=Γ（n）k+1给出-h类DΓ（n）k+1我-1.DΓ（n）k+1, （5） 对于0≤ n<nmax。在进一步发展数学之前，我们暂停一下，以直观地解释RMQ算法是如何进行的。图5描述了发生的过程。顶部面板在时间步k处移动量化器。在每个码字上，传播高斯Neuler更新（第二个面板）。在面板三中，这些更新通过相关原始码字的概率进行加权，并求和以产生时间步k+1处的边缘密度，如最终面板所示。与该边缘密度相关的分布是在时间步长k+1处被量化以产生量化器的分布。重复此过程，直到最终生成量化器。现在我们开始推导牛顿-拉斐逊迭代所需的量。

为了总结符号，我们以asU（γik，Zk+1）=：Uik+1=mikZik+1+cik，（6）量化器（前一时间步）条件密度加权条件密度边际密度的形式编写更新图5：RMQ算法的说明。其中mik：=b（γik）√t和cik：=γik+a（γik）t、 带Zik+1~ N（0，1）同分布于Zk+1。这里，我们为随机变量Zik+1引入了一个新的索引，表示它可能依赖于γik。这在Euler更新的情况下是多余的，因为Zk+1是一个标准的正态随机变量，与起点无关。当我们分析更一般的情况时，这种更一般的表示法是必要的（见第4节）。我们分别用fZik+1和fZik+1表示相应的密度和分布函数。有了这个符号，exk+1isFeXk+1（x）=ZRP（U（y，Zk+1）的近似边缘分布≤ x） dPeXk公司≤ y≈NkXi=1P（Uik+1≤ x） PbXk=γik=NkXi=1H类(-mik）+sgn（mik）FZik+1x个- 锡克米克PbXk=γik, （7） 其中H（·）是Heaviside阶跃函数，sgn（·）是signum函数。最后一步是因为倒数第二行的左手概率可能是writenASP（Uik+1≤ x）=PZik+1≤x个-锡克米克对于mik≥ 01- PZik+1≤x个-锡克米克对于mik<0。需要注意的是，对于Euler更新，MIK是SDE波动率的代表，通常保证其正性，在这种情况下，可以简化（7）。然而，我们坚持这种公式，因为在一般情况下，MIK可能无法保证为正。畸变梯度的要素D（Γk+1）可以写成D（Γk+1）γjk+1=2NkXi=1EhI{Uik+1∈Rj（Γk+1）}γjk+1- Uik+1iPbXk=γik= 2NkXi=1ZUik+1∈Rj（Γk+1）γjk+1- Uik+1dP（Zik+1≤ x） PbXk=γik, （8） 其中1≤ j≤ Nk+1是跟踪tk+1量化器元素的索引，与tk量化器关联的索引是i。

（8）中的积分边界现在必须用积分变量表示。此处，如第2节，Uik+1∈ Rj（Γk+1）等于不等式Rj-k+1<Uik+1≤ rj+k+1，rj±k+1=（γj±1k+1+γjk+1），（9）和r1-k+1=-∞ 和rNk+1+k+1=∞ 根据定义。确定条件归一化区域边界，ri，j±k+1=rj±k+1- cikmik，允许用随机变量Zik+1asUik+1写出不等式∈ Rj（Γk+1）=ri，j-k+1<Zik+1≤ ri，j+k+1用于mik≥ 0ri，j-k+1>Zik+1≥ ri，j+k+1对于mik<0，现在可以将其用作进行积分的范围。直接计算（8），tk+1时畸变梯度的每个元素由D（Γk+1）γjk+1=2NkXi=1h（γjk+1-cik）sgn（mik）FZik+1（ri，j+k+1）- FZik+1（ri，j-k+1）-|迈克|MZik+1（ri，j+k+1）- MZik+1（ri，j-k+1）iPbXk=γik.Fur thermore，三对角黑森犬的对角线，D（Γk+1），由下式得出D（Γk+1）γjk+1=NkXi=12 sgn（mik）FZik+1（ri，j+k+1）- FZik+1（ri，j-k+1）+2 | mik | fZik+1（ri，j+k+1）（γjk+1- γj+1k+1）+2 | mik | fZik+1（ri，j-k+1）（γj-1k+1- γjk+1）PbXk=γik,超对角线和次对角线元素由D（Γk+1）γjk+1γj+1k+1=NkXi=12 | mik | fZik+1（ri，j+k+1）（γjk+1- γj+1k+1）PbXk=γik和D（Γk+1）γjk+1γj-1k+1=NkXi=12 | mik | fZik+1（ri，j-k+1）（γj-1k+1- γjk+1）PbXk=γik,分别地虽然这些表达式可能看起来很复杂，但它们只是密度f函数、累积分布函数和随机变量Zik+1的第一个较低部分期望的总和，因此当这些函数已知时，很容易计算。除了初始猜测之外，现在已经提供了实现牛顿迭代（5）所需的所有细节。

在本文考虑的所有应用中，我们假设Nk=N表示1≤ k≤ K、 并使用前一时间步中的量化器作为初始猜测，即Γ（0）K+1=ΓK.3.1有效实施，如第2.1节所述，其中我们为VQ所需的牛顿-拉斐逊迭代提供了有效的矩阵公式，RMQ也符合矩阵规范。这有助于简化和计算效率的实现。除了猜测Γk+1之外，我们还需要以下时间索引列向量[mk]i=mik，[ck]i=cik，1≤ 我≤ Nkand公司[Γk+1]i=γi+1k+1- γik+1，1≤ 我≤ Nk+1- 1、保留概率pk=[P（bXk=γk），…，P（bXk=γNkk）]的行向量，长度为d的行向量用jd表示。

除了Γk+1，必须在每次Newton-Raphson迭代之前重新计算，其他向量每个时间步计算一次。在每次牛顿-拉斐逊迭代之前，必须根据Γk+1的新估计计算三个矩阵：转移概率的Nk×Nk+1矩阵[Pk+1]i，j=P（bXk+1=γjk+1 | bXk=γik）=sgn（mik）hFZik+1（ri，j+k+1）- FZik+1（ri，j-k+1）i，另一个尺寸相同的矩阵，具有较低的偏矩值[Mk+1]i，j=MZik+1（ri，j+k+1）- MZik+1（ri，j-k+1）（10）和Nk×Nk+1- 1正区域边界处的密度值矩阵[fk+1]i，j=fZik+1（ri，j+k+1）。然后，可以根据这些向量和矩阵将时间步k+1处的失真函数的梯度写为D（Γk+1）= 2件装（千+千）- ckjNk+1oPk+1- （| mk | jNk+1）o Mk+1, （11） 在哪里o 是阿达玛（或元素）产品。（三对角）Hessian矩阵的超对角和次对角元素，D（Γk+1），由向量of=-主键|mk公司|o-1j（Nk+1-（1）ofk+1o (Γk+1 jnk）, （12） 主对角线由hmain=2pkPk+1给出+ho OFF | 0+0 | ho OFF, （13） 在哪里o - 指数中的1表示元素方向的逆。方程式（11）、（12）和（13）提供了实现（5）中牛顿-拉斐逊迭代所需的必要组件。经过必要的迭代次数后，用Pk+1=pkPk+1计算与最终量化器Γk+1相关的概率，其中Pk+1必须根据最终的Γk+1重新计算。

因此，此处给出的矩阵公式允许将R MQ解释为在齐次离散时间马尔可夫链中的传播，其中Γkre表示时间步k处的马尔可夫状态，处于这些状态的概率为pk，并且相关的转移概率矩阵为pk+1。有时在文献中，时间步长k和k+1之间的转移概率矩阵表示为Pk，k+1，我们选择省略第一个指数。3.2示例作为RMQ算法的第一个示例，图6显示了正则（风险中性）几何布朗运动过程的量化器在时间上的演化DST=rStdt+σStdWt，（14）50 100 150 200 250 300 350量化器00.0020.0040.0060.0080.01概率GBM Euler RMQ00Time3000.50.005概率量化器200GBM Euler RMQ（3D）10010.01图6：GBM过程量化器的时间演变。使用参数S=100、r=5%和σ=30%。使用的RMQ参数为T=1，K=12，对于所有K，t=t/K和Nk=200，对于VQ算法，nmax=50，对于RMQ算法，nmax=5。除非另有说明，这些参数用于本文其他计算中使用几何布朗运动的地方。在这些图中，量化器的颜色表示相关的时间步长，蓝色的线接近初始时间，绿色的线接近最终时间。本公约始终保持不变。4高阶RMQ方案鉴于上一节中R MQ算法的一般公式，我们现在探索高阶扩展：Milstein方案和简化的弱阶2.0方案。

只要可以为所有i和k计算CDF、PDF和随机变量Zik+1的较低部分期望值，就可以将SDE的任何数值格式写成有效形式（6），并与RMQ算法一起使用。4.1 Milstein方案Milstein【1975】方案f或（3）中的SDE由eXk+1=eXk+a（eXk）给出t+b（eXk）√tZk+1+b（eXk）b′（eXk）t型Zk+1- t型,对于0≤ k<k，其中t=t/K和Zk+1~ N（0，1），初始值ex=x。通过完成平方，可以将此m写入xk+1=eXk+a（eXk）-b（eXk）b′（eXk）t型-b（eXk）b′（eXk）-1+b（eXk）b′（eXk）t型Zk+1+√tb′（eXk）-1..因此，更新（6）可以以asUik+1=mikZik+1+cik所需的形式写入，其中MIK=b（γik）b′（γik）tandcik=γik+a（γik）-b（γik）b′（γik）t型-b（γik）b′（γik）-1、随机变量Zik+1为非中心卡方分布，具有一个自由度和非中心参数λik+1=√tb′（γik）-重要的是要注意，与Euler Maruyama情况不同，随机变量Zik+1的分布~ χ′2（1，λik+1）现在取决于co-dewordγik。虽然Milstein格式的强收敛阶为1，但与Euler格式相比，两种格式的强收敛阶均为1。因此，虽然Milstein格式在某种程度上比Euler格式更精确，但我们需要对更高的弱阶收敛进行不同的更新，这是在近似金融收益预期时所需要的。我们现在探讨这种格式。4.2弱阶2.0 Taylor格式虽然不可能以所需的格式编写弱阶2.0 Taylor格式，但Kloeden和Platen[1999]简化的弱阶2.0格式是可行的。

此模式由eXk+1=eXk+a（eXk）给出t+b（eXk）√tZk+1+b（eXk）b′（eXk）t（Zk+1- （1）+a′（eXk）b（eXk）+a（eXk）b′（eXk）+b′（eXk）b（eXk）(t） Zk+1+a（eXk）a′（eXk）+a′（eXk）b（eXk）(t） ，对于0≤ k<k，此处为w t=t/K和Zk+1~ N（0，1），初始值ex=x。同样，使用平方完成以所需的形式写入此更新，Uik+1=mikZik+1+cik，其中mik=b（γik）b′（γik）t-6-5-4-3-2日志2(t） -22-20-18-16-14-12-10-8-6-4log2（|力矩误差|）GBM第一动量仪，β=0.998Milstein，β=0.997阶2.0，β=1.856-6-5-4-3-2log2(t） -22-20-18-16-14-12-10-8-6-4log2（|力矩误差|）CEV第一力矩计，β=0.997Milstein，β=0.997阶2.0，β=1.885图7：GBM和CEV的第一力矩的收敛。andcik=γik+a（γik）-b（γik）b′（γik）t+a（γik）a′（γik）+a′（γik）b（γik）(t）-b（γik）+a′（γik）b（γik）+a（γik）b′（γik）+b′（γik）b（γik）t型2b（γik）b′（γik）。这里，Zik+1再次是具有一个自由度的非中心卡方分布，非中心参数由λik+1=b（γik）给出+a′（γik）b（γik）+a（γik）b′（γik）+b′（γik）b（γik）tb（γik）b′（γik）p(t） ！，或者更简洁地说，Zik+1~ χ′2（1，λik+1）。4.3示例为了说明上述方案的准确性，RMQ算法应用于几何布朗运动（如（14）所述）和恒定方差弹性（CEV）过程。CEV工艺的SDE isdSt=rStdt+σSαtdWt，在以下示例中，选择的工艺特定参数为S=100，r=5%，α=0.7，σ=σLNS1-α、 σ表示瞬时对数正态波动率σLN=30%。对于GBM，使用第3.2节中的参数。对于GBM和CEV，也使用了该节中提到的RMQ特定参数，除非另有说明。在图7中，我们给出了弱阶收敛的数值证据。