鞅最优运输问题的计算方法

dDk：=d∏kr=1n（r），（8）等价于以下LP问题：maxp∈RD+，δ∈RD+。。。，δN-1.∈RDN公司-1+Xi，。。。，iNpi，。。。，股份有限公司xixNiN公司s、 t.Xi，。。。，ik-1，ik+1，。。。，iNpi，。。。，iN=αkik，用于ik∈ Ikand k=1，N- δki，。。。，ik，j≤Xik+1，。。。，iNpi，。。。，在里面xk+1ik+1，j- xkik，j≤ δki，。。。，ik，j，用于ik∈ Ik，j∈ J和k=1，N、 Xi，。。。，ik，jδki，。。。，ik，j≤ εn，对于k=1，N- 1，我们回忆xkik的地方=xkik，1，xkik，d.在获得了一般的收敛结果之后，我们接下来讨论Pεn（un）的收敛速度问题。我们给出了实线上单步模型收敛速度的估计。据我们所知，下面的误差界是文献中的第一个。定理2.5。设N=2，d=1，或等效地，u=（u，ν）和c:R→ R、 除了定理2.2的条件外，我们假设sup（x，y）∈Ryyc（x，y）< ∞ ν有一个有限的秒矩。然后存在C>0，这样Pεn（un，νn）-P（u，ν）≤ C inf>0λn（R），对于所有n≥ 1，imsart aap版本。2014年10月16日文件：数字MOT\\U修订版。tex日期：2019年4月8日MOT问题的计算方法7其中λn：（0，∞) → R由λn（R）给出：=（R+1）εn+Z(-∞,-R）∪（R，∞)|y |-Rν（dy）。特别是，收敛速度与εnif supp（ν）有界成正比。我们将定理2.2和2.5的证明推迟到第4节，并在本节结束时讨论如何将定理2.2应用于解决金融中出现的其他约束OT问题。备注2.6。一般来说，分布u，un当d≥ 2、对于k=1，N、 让Sk：=S（1）k，S（d）k, 式中，S（i）k代表时间k时ithstock的价格。那么，实际上，只有看涨期权的价格S（i）k- K+, orput选项K-S（i）k+, 对于一组有限的罢工，K在市场上流动可用。

即使对所有可能的罢工K引用看涨期权，也只能得到分布uK，iofS（i）K。因此，这导致了一个改进的优化问题。表示~uk：=（uk，1，…，uk，d）和~u：=（~uk）1≤k≤N、 设Mε（~u）为ε的集合-拟鞅测度psating PoS（i）k-1=uk，i，对于k=1，N和i=1，d、 然后，我们用pε（~u）：=supP定义优化问题∈Mε（~u）EPc（S，…，SN）.（9） ε=0的问题（9）首先由Lim提出，并在[36]中被称为多重鞅最优运输。虽然本文重点讨论P（u）的数值计算，但我们强调定理2.2允许立即扩展到近似P（~u）。P（u）的数值格式：概率离散化。受定理2.2和推论2.4的启发，我们接下来开发了一种基于边缘分布的可测离散化计算P（u）的数值格式。关键是选择合适的序列（un）n≥1当k=1时，N，（a）uNK在有限集上受支持xkik：ik∈ Ik,（b） 重量unk{xkik}可以显式计算，也可以用先验已知的精度近似计算，（c）W的上界unk，uk很容易获得。如上所述，该问题与概率度量的最佳量化密切相关，其目标是最佳逼近给定的概率度量u∈ P通过具有给定数量支撑点的离散度量。对于给定u，其nd-量化un与（xi）1相关≤我≤nd公司 Rdand（Ei）1≤我≤ndis由un（dx）定义：=Pndi=1u[Ei]δxi（dx），其中（Ei）1≤我≤ndis au-分区，即u工程安装∩Ej= 0表示所有i 6=j和u∪1.≤我≤ndEi公司= 因此，nd-u的最佳量化是inf的解决方案ndXi=1ZEix个- xiu（dx）,（10） inf接管所有（xi）1≤我≤ndandu-隔板（Ei）1≤我≤nd。

我们陈述了Graf和Luschgy【25】的收敛结果，另见【26，20】。imsart aap版本。2014年10月16日文件：数字MOT\\U修订版。tex日期：1998年4月8日MOT问题计算方法第3.1条（Graf和Luschgy）。对于每个n≥ 1，（10）中的inf可以通过nd实现-最佳量化器（x*i） 1个≤我≤ndand（E*i） 1个≤我≤nd。Letun*成为相应的优化器，然后limn→∞nW（un*, u）存在且不确定。以下为u：=limn→∞nW（un*, u）存在nu≥ 1使W（un*, u) ≤ （Au+1）/n，对于所有n≥ nu。尽管它们在理论上很有吸引力，但在实践中，使用最佳量化器是有问题的，因为上面的关键量，例如Au和nu，通常是未知的。类似地，通常情况下，数量un*{x*i}很难精确计算或以规定的精度近似计算。为了克服这些困难，我们采用了两种不同的离散化方法，这两种方法都可以在实践中实现。我们的第一种方法，我们称之为确定性离散化，适用于给定边缘u，u在能够计算出针对它们的积分的意义上。情况就是这样，例如，当u，un已知密度函数。第二种方法称为随机离散化，适用于我们能够从边缘取样的情况。在整个第3节中，我们需要以下可积条件，它对应于电力期权的市场价格是有限的。假设3.2（θth-力矩）。存在θ>1和Mθ<∞ 因此Zrd | x |θuN（dx）≤ Mθ。注意，根据Jensen不等式，上述条件表示RRD | x |θuk（dx）≤ Mθfork=1，N、 此外，每当我们考虑以下通用度量u时，我们还将假设满足假设3.2.3.1。确定性离散化。我们设计了一种简单的离散化方法，当supp（u）有界时，该方法具有与最优量化相同的渐近效率。

我们假设u是已知的，因为概率u【E】对于所有E都是已知的∈B（Rd）。我们从这个理想的设置开始，然后考虑已知密度的情况，这允许以一定的精度计算uE。步骤1：截断。对于R>0，让BR rD表示由br定义的方框：=nx=（x，…，xd）：| xi |≤ R、 对于i=1，做然后一个有x个∈ Rd：| x |≤ R BR公司x个∈ Rd：| x |≤ 博士. 取R，使u【BR】>0，并将u截断为概率度量uR（dx）：=1BR（x）u（dx）+u【BcR】δ（dx），其中BcR：=Rd\\BR。显然，BR支持uRis。考虑从u中提取的随机变量X，观察1BR（X）X是根据uR分布的。通过定义巴斯斯坦距离，我们得到了W（uR，u）≤ EBR（X）X- 十、=ZBcR | x |u（dx）≤ Mθ/Rθ-1，（11）产生limR→∞W（uR，u）=0。imsart aap版本。2014年10月16日文件：数字MOT\\U修订版。tex日期：2019年4月8日MOT问题的计算方法9步骤2：离散化。表示为Ohmn Rd由所有q=（q，…，qd）的元素SQ/n组成的可数子空间∈ Zd。对于每个q∈ Zd，我们用V（q/n）表示 Rd x=（x，…，xd）的子集，使得bnxc=q，即对于i=1，…，bnxic=qi，d、 a的位置∈ R、 美国银行∈ Z是小于或等于a的最大整数。我们构造了一个概率度量u（n），其支持度包含在Ohmnbyu（n）{q/n}:= uV（q/n）. 然后u（n）∈ 所有f的P满意度∈ ∧，Zfdu（n）=Xq∈Zdf（q/n）u（n）{q/n}=Zf（n）du，（12），其中f（n）：Rd→ R由f（n）（x）定义：=fbnxc/n. 这意味着鉴于（4）Wu（n），u= supf公司∈ΛZfdu（n）-Zfdu≤ supf公司∈∧Zf（n）- fdu≤ d/n，其中第二个不等式为（12）。请注意，如果supp（u）是有界的，那么supp也是有界的u（n）,距离Wu（n），u为1/n级，与W（un）相同*, u).步骤3：选择参数。

将u替换为uRin步骤2，一个hasWu（n）R，u≤ Wu（n）R，uR+ W（uR，u）≤ d/n+Mθ/Rθ-1、根据杨氏不等式得出| a |γ+| b |θ≥ γ1/γθ1/θ| ab |对于所有a、b∈ R、 其中γ>1是θ的共轭数，即1/θ+1/γ=1。分别设置a=（d/n）1/γ和b=Mθ/Rθ-1.1/θ，它认为d/n+Mθ/Rθ-1.≥ （γd）1/γ（θMθ）1/θ/（Rn）1/γ，Rθ可以达到质量-1=θMθn/γd.自supp基数u（n）R是比例（Rn）D，它决定了相应LP问题中的变量数量，设置r=Rn：=θMθn/γd1/(θ-1） 得出固定计算复杂度的最佳上界，即Wu（n）Rn，u≤ γd/n。分别用ukF替换u，k=1，N，我们获得uN=（unk）1≤k≤n按照上述步骤，unk：=u（n）k，Rn。那么定理2.2 yieldslimn→∞PNγd/n（un）=P（u）。备注3.3。通常，un可能不再属于P, 即使u∈ P. 当d=1时，第4.2节给出了一个保持递增凸阶的显式离散化。在最近的平行研究中，Alfonsi等人【1】研究了构建un的方法∈ P.通过上述分析，我们可以构建近似的度量值u，并将这些值uV（q/n）已知所有q/n∈ Ohmn、 这可能是可能的，例如，当u是原子时，但总的来说，我们需要讨论如何很好地近似这些值。我们这样做是为了测量密度函数，即u（dx）=ρ（x）dx。在这种情况下，对于某些xq，一个简单的点估计ρ（xq）/nd∈ V（q/n），提供接近u（n）的自然候选值{q/n}. 然而，使用定理2.2，我们需要以显式和非渐进的方式限制结果度量和u（n）之间的Wasserstein距离。如前所述，为了简单起见，我们将u截断，将R设置为整数m。

在上支持u（n）和▄u（n）mbeOhmn∩ b定义如下：如果0 6=q/n∈ Bm，然后是u（n）m{q/n}:=ZV（q/n）ρ（x）dx和￠u（n）m{q/n}:= ρ（xq）/nd，imsart aap版本。2014年10月16日文件：数字MOT\\U修订版。tex日期：201910年4月8日MOT问题的计算方法∈ V（q/n）任意选择，u（n）m{0}:= 1.-Xq/n6=0u（n）m{q/n}和￠u（n）m{0}:= 1.-Xq/n6=0u（n）m{q/n},其中，上述总和确实为u（n）m{q/n}= u（n）m{q/n}= 0表示q/n/∈ Bm。A以上，Wu（n）m，u≤ d/n+Mθ/Mθ-1以下给出了W的上限u（n）m，u（n）m.提案3.4。假设3.2成立。假设ρ是连续的，或者说，对于每个R>0，存在κR:[0，∞) → 非递减的R，使得κR（0）=0和ρ（x）-ρ（y）≤ κR（| x- y |），对于所有x，y∈ BR。（i） 如果u具有有界支撑，即supp（u） br对于某些R>0，则为新u（n）m，u≤ εn：=d/n+2dd（R+1）d+1κR+1（d/n），对于所有m≥ bRc+1。（13） （ii）如果ρ是一致连续的，即存在一个一致的κ=κrf，对于所有R>0，则u（n）m，u≤ εm，n：=d/n+mθ/mθ-1+2ddmd+1κ（d/n）。（14） （iii）如果{xq}06=q/n∈b满足度ρ（xq）≤ ρ（x）表示所有x∈ V（q/n），然后新u（n）m，u≤ τm，n：=d/n+mθ/mθ-1+inf1≤j≤mnddjd+1κj（d/n）+4Mθ/jθ-1o。（15） 证明。在条件（i）下，一个具有u=umas m≥ bRc+1，因此Wu（n）m，u=Wu（n）m，um≤ d/n.注释也支持u（n）m, 支持u（n）m BbRc+1定义。对于任何f∈ ∧，它适用于m≥ bRc+1Zfd？u（n）m-Zfdu（n）m=X06=q/n∈BbRc+1f（q/n）ZV（q/n）ρ（x）-ρ（xq）dx公司≤X06=q/n∈BbRc+1 | q/n | ZV（q/n）κR+1（d/n）dx≤ 2dd（R+1）d+1κR+1（d/n），产生Wu（n）m，u（n）m≤ 2dd（R+1）d+1κR+1（d/n），进一步（13）。至于（ii），我们推断u（n）m，u（n）m≤ 2ddmd+1κ（d/n），使用相同的Bm参数，从而获得（14）。或者，假设第三个条件成立。

对于每个整数1≤ j≤ m、 一个有Zfd？u（n）m-Zfdu（n）m≤X06=q/n∈Bj | q/n | ZV（q/n）ρ（x）-ρ（xq）dx公司+Xq/n∈Bm\\Bj | q/n | ZV（q/n）ρ（x）-ρ（xq）dx公司≤ 2ddjd+1κj（d/n）+2ZBcjd+| x|ρ（x）dx≤ 2ddjd+1κj（d/n）+4Mθ/jθ-1、因此Wu（n）m，u（n）m≤ inf1≤j≤m级ddjd+1κj（d/n）+4Mθ/jθ-1.导出了和（15）。imsart aap版本。2014年10月16日文件：数字MOT\\U修订版。tex日期：2019年4月8日MOT问题的计算方法11通过简单的计算，一个具有limm→∞画→∞εm，n=limm→∞画→∞τm，n=0。因此，存在合适的序列（mn）n≥1和（nm）m≥1，这样Limn→∞εmn，n=limn→∞τmn，n=0和limm→∞εm，nm=limm→∞τm，nm=0。请注意，前面的选项mn：=bRnc可能不会产生limn→∞εbRnc，n=limn→∞τbRnc，n=0，这些序列必须根据ρ计算。然而，如果ρ为L-Lipschitz，则有εm，n=d/n+mθ/mθ-1+2dmdmd+1L/n，我们推断取（mn）n是必要的≥1或（nm）m≥1如此限制→∞md+1n/n=0或limm→∞md+1/nm=0。将所有内容放在一起，分别取ukin代替u，表示k=1，N、 上述程序产生一个测量向量|u（N）m=u（n）k，m1.≤k≤N、 在命题3.4的条件下，我们有limn→∞PNεmn，nu（n）mn= 画→∞PNτmn，nu（n）mn= P（u）或limm→∞PNεm，nmu（nm）m= limm公司→∞PNτm，nmu（nm）m= P（u）。3.2. 随机离散化。我们现在考虑一种不同的离散化程序，该程序适用于有一个黑箱根据u生成独立随机变量的情况。提供i.i.d.u序列-分布随机变量（Xn）n≥1，定义经验测量值μnbyμn（dx）：=nXi=1nδXi（dx）。定义^u是一个随机度量，并遵循Glivenko-Cantelli定理，参见Fournier和Guillin[22]，limn→∞W（un，u）=0几乎可以肯定，limn→∞EW（un，u）=0.通过将u替换为ukF，构建随机度量值^unk，对于k=1，N并设置^uN：=（^unk）1≤k≤N

与定理2.2相比，我们现在得到了一个随机收敛结果。提案3.5。让定理2.2的条件成立。给定序列（εm）m≥1.(0, ∞) 收敛到零，一个有limm→∞画→∞Pεm（un）=P（u）几乎可以肯定。此外，对于任何子序列（^nm）m≥1这样的下午≥1E级W⊕（μm，u）/εm<∞, 几乎可以肯定→∞Pεm（uunm）=P（u）。证据对于每个固定的εm>0，一个具有以下推论4.3的不等式Pεm（μun）-P（u）≤ Lip（c）εm+P2εm（u）-P（u）每当W⊕（un，u）≤ εm，（16）其中，我们记得Lip（c）是c的Lipschitz常数。应用Glivenko-Cantelli定理，可以得到limn→∞Pεm（μun）-P（u）≤ Lip（c）εm+P2εm（u）-P（u）=：δ最大值。利用命题4.4，我们得到了limm→∞P2εm（u）=P（u），第一个断言的ConvergencerResult如下。对于任何δ>0，存在mδ，使得εm≤ 所有m的δ≥ mδ。在语句中使用^nma，我们有xm≥mδPhPεm（^u^nm）-P（u）> δi≤Xm公司≥mδPhW⊕^^nm，u> εmi≤Xm公司≥mδEW⊕（μm，u）/εm，根据Borel-Cantelli引理，意味着limm→∞Pεm（^unm）=P（u）几乎可以肯定。imsart aap版本。2014年10月16日文件：数字MOT\\U修订版。tex日期：201912年4月8日。MOT问题的计算方法。确切地说，如果εm>0，则必须关注LP问题Pεm（μun），以便⊕（un，u）≤ εmoccurs具有高概率。为此，为了选择合适的序列（εm）m≥1和（^nm）m≥1，我们需要量化EW⊕（un，u）. 幸运的是，文献[22]中的定理1在假设3.2下提供了这样的估计，尽管省略了一些情况，例如d=1、2和θ=2。为了完整性，我们将所有情况都考虑在内，将此结果表述为引理3.6。引理3.6（[22]）。假设3.2成立。

存在C（θ，d）>0，使得EW⊕un，u≤χ（n）表示所有n≥ 1，式中χ（n）：=NC（θ，d）n1/θ-1如果d=1且1<θ<2，（1+对数n）n-1/2如果d=1，θ=2，n-1/2如果d=1且θ>2，n1/θ-1如果d=2且1<θ<2，1+（对数n）n-1/2如果d=2，θ=2，（1+对数n）n-1/2如果d=2且θ>2，n1/θ-1如果d≥ 3和1<θ<d/（d- 1） ，（1+对数n）n-1/dif d≥ 3和θ=d/（d- 1） ，n-1/dif d≥ 3和θ>d/（d-1).因此，一个人有PW⊕un，u> εm≤ χ（n）/εm。我们注意到，常数C（θ，d）在[22]中没有明确规定。为了实现我们的方案，我们需要确定^nmand，为此，我们必须显式地计算C（θ，d）。这是可能的，主要遵循了[22]中的论点，但冗长乏味，推迟到附录A。3.3. 数值示例。我们现在讨论一些具体的MOT问题，以说明定理2.2以及我们的离散化方案是如何应用的。我们的大多数示例都提供了一个封闭形式的优化器或其分析特性，这允许我们验证数值结果。我们记得，对于d=1，我们只写x=x和Sk=Sk。示例3.7。Beiglb¨ock和Juillet在[5]一个特殊问题中研究了N=2，d=1，c（x，y）=h（x- y） ，其中h:R→ R有一个严格的凸导数。Let（u，ν）∈ P. [5]的定理1.7表明，如果u有密度ρ，则存在两个可测函数ξ±：R→ 使唯一优化器P*∈ ξ±上支持P（u，ν）的M（u，ν），即P*（dx，dy）=u（dx）x个-ξ-（x） ξ+（x）-ξ-（x） Δξ+（x）（dy）+ξ+（x）-xξ+（x）-ξ-（x） Δξ-（x） （dy）,式中ξ-（十）≤ x个≤ ξ+（x），ξ+（x）<ξ+（x）和ξ-（十）/∈ξ-（x） ，ξ+（x）对于所有x，x∈ R，x<x。我们想用数值说明上述结果。设ρ是由ρ（x）定义的截断伽玛函数：=1[0,1]（x）x3/2e-x/C，其中C：=Zx3/2e-xdx>1/5。imsart aap版本。

2014年10月16日文件：数字MOT\\U修订版。tex日期：2019年4月8日MOT问题的计算方法13下一步构造ν（dx）=σ（y）dy byσ（y）：=ρ（y/2）/6+4ρ（2y）/3=1[0,2]（y）（y/2）3/2e-y/2/C+1[0,1/2]（y）（2y）3/2e-2年/年。根据结构，一个具有（u，ν）∈ P, supp（u）=[0，1]和supp（ν）=[0，2]。此外，可以验证ρ和σ是L-L=7的supp（u）和supp（ν）上的Lipschitz。将第3.1节的分解应用于u和ν，我们得到了|u（n）和|ν（n），支持识别号：0≤ i<n和电话：0≤ j<2n, 定义为u（n）识别号（i/n）:=(1 -Pn编号-1k=1ρ（xk）/n如果i=0，ρ（xi）/n如果1≤ i<n，￠ν（n）日本/日本:=(1 -P2n-1k=1σ（yk）/n如果j=0，σ（yj）/n如果1≤ j<2n，其中xi∈识别号，（识别号+1）/识别号和yj∈j/n，（j+1）/n对于i=1，n-1和j=1，2n个-1、根据命题3.4，W⊕（|u（n），|ν（n）），（u，ν）≤ （3L+2）/n=：εn。那么相应的LP问题如下：max（pi，j）∈R2n+n-1Xi=02n-1Xj=0pi，jh（一）-j） /不适用s、 t.2n-1Xj=0pi，j=αni，对于i=0，n- 1，n-1Xi=0pi，j=βnj，对于j=0，2n个-1，n-1Xi=02n个-1Xj=0pi，jj/n-αnii/n≤ εn，式中αni：=u（n）识别号（i/n）和βnj：=￠ν（n）日本/日本. 以h（x）：=ex为例，我们使用Gurobi solver解决了LP问题，结果如图1所示。左侧窗格显示值图1。计算，例如3.7。第一个窗格显示10的值Pεn（|u（n），|ν（n））≤ n≤ 第二个窗格绘制n=100时优化器的热图。Pεn（|u（n），|ν（n））10≤ n≤ 200，这在数值上显示了Pεn（|u（n），|ν（n））的收敛性。右侧窗格显示优化器的热图（p*i、 j）对于n=100。我们看到theimsart aap版本。