贵州大学数学与统计学院，贵州 贵阳

收稿日期：2023年4月10日；录用日期：2023年6月8日；发布日期：2023年6月15日

摘要

本文使用多水平蒙特卡洛梯度投影优化方法求解一类随机最优控制问题。蒙特卡洛方法，是一种广泛使用的方法，是求解最优控制问题的一种常用方法，但它有自身的局限性，收敛速度慢。我们选取收敛速度较快的多水平蒙特卡洛方法。对这种方法的理论进行分析，之后将该方法用于具体实例，求解金融问题，并通过数值实验验证方法的有效性。

关键词

随机微分方程，最优控制，多水平蒙特卡洛方法，梯度投影，期权定价

Application of Optimal Control Based on Multilevel Monte Carlo Method in Finance

Peng Lei

School of Mathematics and Statistics, Guizhou University, Guiyang Guizhou

Received: Apr. 10^th, 2023; accepted: Jun. 8^th, 2023; published: Jun. 15^th, 2023

ABSTRACT

This paper uses a multi-level Monte Carlo gradient projection optimization method to solve a class of stochastic optimal control problems. Monte Carlo method is a widely used method for solving optimal control problems, but it has its own limitations and slow convergence speed. We selected a multi-level Monte Carlo method with fast convergence speed to analyze the theory of this method, and then applied it to specific examples to solve financial problems. The effectiveness of the method was verified through numerical experiments.

Keywords:Stochastic Differential Equation, Optimal Control, Multilevel Monte Carlo Method, Gradient Projection, Option Pricing

Copyright © 2023 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

随机最优控制问题在工程、金融和金融领域得到了广泛的应用，许多学者研究了不同情况下的随机最优控制问题，并将其使用在了求解金融问题。对于随机微分方程，有时候很难找到显式解，因此数值逼近是一个很好的选择。近年来有了一些方法，比如说：有限维随机规划 [1] [2] ，随机最大值原理 [3] ，Bellman动态规划原理 [4] 和鞅方法 [5] ，近年来广泛地应用于数值求解。

金融安全是国家安全的重要体现，当初亚洲金融危机，世界都受其影响，因此对于金融问题的研究很重要，也是国家安全的保障。因此用随机微分方程这一个工具，去研究金融问题，用数值模拟去求解随机微分方程，是本文研究的目的。

蒙特卡洛方法是一种常见的求解随机问题的方法，但收敛速度慢。于是对它加以改进，改进的方法有很多，如多水平蒙特卡洛方法 [6] 和拟蒙特卡洛(QMC)方法 [7] 。对于一个随机微分状态方程，当控制是确定之时，求解随机最优控制问题应当考虑它的不确定性，此时梯度经常包含求期望，多水平蒙特卡洛方法是求解这类问题的一种好方法。

本文将多水平蒙卡洛方法与梯度投影方法相结合，求解随机最优控制问题。为了减少梯度计算带来的统计误差和离散误差的影响，我们使用MLMC方法来估计每次迭代中的梯度，其中的均方误差(MSE)被动态地更新。对这种方法给出理论依据，将金融问题转化为随机最优控制问题，进行数值模拟。

2. 随机最优控制问题

给定一个完全的概率空间 $\left(\Omega ,F,{\left\{F_{\text{t}}\right\}}_{\text{t}\ge \text{0}},P\right)$ ， $L_F^{\text{2}}\left(\left[\text{0},T\right];R\right)$ 是一个实值平方可积的 $F_t$ -适应过程空间，于是有 ${\Vert y\Vert }_{L^2\left(\Omega ,L^2\left[0,T\right]\right)}<\infty$ 。 ${\left\{F_t\right\}}_{t\ge 0}$ 是一维标准布朗运动 ${\left\{W_t\right\}}_{t\ge 0}$ 生成的自然过滤。定义：

${\Vert y\Vert }_{L^2\left(\Omega ,L^2\left[0,T\right]\right)}={\left({\displaystyle {\int }_{\Omega }{\displaystyle {\int }_0^T{\left|y_t\right|}^2\text{d}t\text{d}P\left(\omega \right)}}\right)}^{\frac{1}{2}}$ (2.1)

我们给定如下的一个目标函数：

$\min J\left(y,u\right)={\displaystyle {\int }_0^{\text{T}}E\left[h\left(y\right)\right]\text{d}t}+{\displaystyle {\int }_0^{\text{T}}j\left(u\right)\text{d}t}$ (2.2)

其中h和j都是具有一阶连续导数的光滑函数， $u\in U_{ad}$ 是一个确定性控制， $U_{ad}$ 是一个在控制空间 $L^2\left(0,T\right)$ 下的闭凸控制集。 $y\left(u\right)\in L_F^2\left(\left[0,T\right];R\right)$ 是由如下随机微分方程生成的随机过程：

$\text{d}y=f\left(t,y,u\right)\text{d}t+g\left(t,y\right)\text{d}{W}_t,y\left(0\right)=y_0$ (2.3)

设f是一个关于 $t,y,u$ 连续可微的函数，g是一个关于 $t,y$ 连续可微的函数。可得方程(2.2)存在唯一解 $y(\cdot )\in L_F^2\left(\left[0,T\right];R\right)$ ， $\left(y_0,u(\cdot )\right)\in R\times U_{ad}$ 。于是函数 $y(\cdot )$ 是关于 $u(\cdot )$ 相对应的， $u^{\ast }$ 是方程(2.2)的最优解。这个问题称之为随机最优控制问题。

3. 梯度投影方法

函数 $J\left(u\right)=J\left(y\left(u\right),u\right)$ ， $y\left(u\right)$ 是(2.2)的解。假设 $J\left(u\right)$ 是一个凸函数，U是一个希尔伯特空间，K是U上的一个闭凸子集。 $b\left(\cdot ,\cdot \right)$ 是一个对称且正定的双线性形式，定义B： $U\to U$ 满足 $\left(Bu,v\right)=b\left(u,v\right)$ 。根据投影运算符 $P_k:U\to K$ 以及一阶优化的必要条件，可以得到

$u^{\ast }=P_K\left(u^{\ast }-\rho B^{-1}{J}^{\prime }\left(u^{\ast }\right)\right)$ (3.1)

其中 $J^{\prime }\left(u\right)$ 是由 $J\left(u\right)$ 生成的Gateaux导数， $\rho$ 一个正常数 [8] 。

设 $0=t_0^N<t_1^N<\cdots <t_N^N=T$ ， $t_{n+1}^N-t_n^N=T/N=:h$ ， $I_n^N=\left[t_{n-1}^N,t_n^N\right]$ ，则分段常数空间 $U_N$ 可以如下定义：

$U_N=\left\{u\in U:u={\displaystyle \underset{n=1}{\overset{N}{\sum }}{\alpha }_n{X}_{I_n^N}a.e.,{\alpha }_n\in R}\right\}$ (3.2)

其中 $X_{I_n^N}$ 是 $I_n^N$ 的示性函数。于是公式(3.1)可以做如下的近似计算：

$u^{\ast ,N}=P_{K_N}\left(u^{\ast ,N}-\rho B^{-1}{J}^{\prime }\left(u^{\ast ,N}\right)\right)$ (3.3)

其中 $K_N=K\cap U_N$ 。

为了对(2.1)和(2.2)数值近似求解，我们可以采用下面的迭代公式：

$\{\begin{array}{l}b\left(u_{i+0.5}^N,v\right)=b\left(u_i^N,v\right)-{\rho }_i\left({J^{\prime }}_N\left(u_i^N\right),v\right),v\in K_N,\\ u_{i+1}^N=P_{K_N}^b\left(u_{i+0.5}^N\right),\end{array}$ (3.4)

${\rho }_i$ 是迭代步长， ${J^{\prime }}_N$ 是 $J^{\prime }(\cdot )$ 的数值近似。用下面式子表示 ${J^{\prime }}_N$ 与 $J^{\prime }(\cdot )$ 之间的误差

${\epsilon }_N=\underset{i}{\sup }\Vert J^{\prime }\left(u_i^N\right)-{J^{\prime }}_N\left(u_i^N\right)\Vert$ (3.5)

于是有如下定理3.1。

定理3.1：若 $J^{\prime }(\cdot )$ 满足利普西茨条件，并且关于 $u^{\ast }$ 和 $u^{\ast ,N}$ 是一直单调的，于是存在正常数c和C满足

$\Vert J^{\prime }\left(u^{\ast }\right)-J^{\prime }\left(v\right)\Vert \le C\Vert u^{\ast }-v\Vert \forall v\in K,$ (3.6)

$\left(J^{\prime }\left(u^{\ast }\right)-J^{\prime }\left(v\right),u^{\ast }-v\right)\ge c{\Vert u^{\ast }-v\Vert }^2\forall v\in K,$ (3.7)

$\Vert J^{\prime }\left(u^{\ast ,N}\right)-J^{\prime }\left(v\right)\Vert \le C\Vert u^{\ast ,N}-v\Vert \forall v\in K_N,$ (3.8)

$\left(J^{\prime }\left(u^{\ast ,N}\right)-J^{\prime }\left(v\right),u^{\ast ,N}-v\right)\ge c{\Vert u^{\ast ,N}-v\Vert }^2\forall v\in K_N,$ (3.9)

假设

$N\to \infty$ 时， ${\epsilon }_N=\underset{i}{\sup }\Vert J^{\prime }\left(u_i^N\right)-{J^{\prime }}_N\left(u_i^N\right)\Vert \to 0$ ， (3.10)

对于一些常数 $0<\delta <1$ ，我们可选择 ${\rho }_i$ ，使得满足

$0<1-2c{\rho }_i+\left(1+2C\right){\rho }_i^2\le {\delta }^2$ (3.11)

则(3.4)是收敛的。更准确的，我们有

当 $i\to \infty ,N\to \infty$ 时， $\Vert u^{\ast }-u_i^N\Vert \to 0$ 。

根据随机最大值原理，可以用下面的修正公式方便地计算目标函数的梯度：

$-d_p=\left[h^{\prime }\left(y\right)+p{f^{\prime }}_y\left(t,y,u\right)-p{g^{\prime }}_y{\left(t,y\right)}^2\right]\text{d}t+p{g^{\prime }}_y\left(t,y\right)\text{d}{W}_t,p\left(T\right)=0,$ (3.12)

$J\left(u\right)$ 的导数可以由下面式子表示：

$J^{\prime }\left(u\right)v={\displaystyle {\int }_0^{\text{T}}\left(E\left[p{f^{\prime }}_u\left(t,y,u\right)\right]+j^{\prime }\left(u\right)\right)v\text{d}t},\forall v\in U$

根据Riesz表示定理，可得

$J^{\prime }\left(u\right)=E\left[p{f^{\prime }}_u\left(t,y,u\right)\right]+j^{\prime }\left(u\right).$ (3.13)

上述就是梯度投影方法。由于梯度包含期望，因此期望估计的准确性直接影响梯度投影算法的收敛性。通常，它需要大量样本。特别地，当扩散系数 $g\left(t,y\right)$ 大时样本数量需要更多。因此，我们使用MLMC方法来估计期望。

4. 基于梯度投影的蒙特卡洛方法

我们考虑期望 $E\left[A\right]=E\left[f\left(y\right)g\left(p\right)\right]$ ，其中 $y,p$ 时状态方程和(3.12)的解。我们假设 $f,g$ 有连续的导数，通过 $A_{h_L}=f\left(y_{h_L}\right)g\left(p_{h_L}\right)$ 作数值近似，在理论分析之前，我们做如下的假设：

H1：y的误差估计满足 ${\Vert I_L{y}_{h_L}-y\Vert }_{L^2\left(\Omega ,L^2\left[0,T\right]\right)}\underset{˜}{<}{h}_L^{{\beta }_y}$ 。

H2：修正的状态p的误差估计满足 ${\Vert I_L{p}_{h_L}-p\Vert }_{L^2\left(\Omega ,L^2\left[0,T\right]\right)}\underset{˜}{<}{h}_L^{{\beta }_p}$ 。

H3：对样本估计的计算成本 $C_L$ 满足 $C_L\underset{˜}{<}{h}_L^{-{\gamma }_1}+h_L^{-{\gamma }_2}$ 。

在这假设下我们给出如下的定理。

定理4.1：假设H1和H2成立， $A\in L^2\left(\Omega ,L^2\left[0,T\right]\right)$ ，f和g满足利普西茨连续，则多水平蒙特卡洛法(MLMC)对A的期望误差满足：

${\Vert E\left[A\right]-Y\Vert }_{L^2\left(\Omega ,L^2\left[0,T\right]\right)}\underset{˜}{<}{h}_L^{{\beta }_y}+h_L^{{\beta }_p}+{\displaystyle \underset{l=0}{\overset{L}{\sum }}{N}_l^{-\frac{1}{2}}\left(h_l^{{\beta }_y}+h_l^{{\beta }_p}\right)}.$ (4.1)

定理4.2：假设H1和H2和H3都成立， $A\in L^2\left(\Omega ,L^2\left[0,T\right]\right)$ ，f和g满足利普西茨连续，则通过下面样本数 ${\left\{N_l\right\}}_{l=0}^L$ 的选择，可以得到MLMC估计 [7] 。

$N_l=\{\begin{array}{l}O\left({\left(h_l^{-{\gamma }_1}+h_l^{-{\gamma }_2}\right)}^{-\frac{2}{3}}{\left(h_l^{{\beta }_y}+h_l^{{\beta }_p}\right)}^{\frac{2}{3}}{\left(h_L^{{\beta }_y}+h_L^{{\beta }_p}\right)}^{-2}\right),\tau >0,\\ O\left({\left(h_l^{-{\gamma }_1}+h_l^{-{\gamma }_2}\right)}^{-\frac{2}{3}}{\left(h_l^{{\beta }_y}+h_l^{{\beta }_p}\right)}^{\frac{2}{3}}{\left(L+1\right)}^2{\left(h_L^{{\beta }_y}+h_L^{{\beta }_p}\right)}^{-2}\right),\tau =0,\\ O\left({\left(h_l^{-{\gamma }_1}+h_l^{-{\gamma }_2}\right)}^{-\frac{2}{3}}{\left(h_l^{{\beta }_y}+h_l^{{\beta }_p}\right)}^{\frac{2}{3}}{h}_L^{\frac{2\tau }{3}}{\left(h_L^{{\beta }_y}+h_L^{{\beta }_p}\right)}^{-2}\right),\tau <0.\end{array}$ (4.2)

其中

$\tau =\min \left\{2{\beta }_y-{\gamma }_1,2{\beta }_p-{\gamma }_1,{\beta }_p+{\beta }_y-{\gamma }_1,2{\beta }_y-{\gamma }_2,2{\beta }_p-{\gamma }_2,{\beta }_p+{\beta }_y-{\gamma }_2\right\}.$ (4.3)

有 ${\Vert E\left[A\right]-Y\Vert }_{L^2\left(\Omega ,L^2\left[0,T\right]\right)}\underset{˜}{<}{h}_L^{{\beta }_y}+h_L^{{\beta }_p}$ 。

当 $L\to \infty$ 时，MLMC的总的计算成本 $C_{mlmc}$ 是渐进有界的，如下：

$C_{mlmc}\underset{˜}{<}\{\begin{array}{l}{\left(h_L^{{\beta }_y}+h_L^{{\beta }_p}\right)}^{-2},\tau >0,\\ {\left(L+1\right)}^3{\left(h_L^{{\beta }_y}+h_L^{{\beta }_p}\right)}^{-2},\tau =0,\\ h_L^{\tau }{\left(h_L^{{\beta }_y}+h_L^{{\beta }_p}\right)}^{-2},\tau <0.\end{array}$ (4.4)

基于梯度投影的MLMC方法的最优控制，它的算法如下：

1) 输入 $\tau ,{\epsilon }^{\left(1\right)},i_{\max },u^0,\eta$ ，给出迭代误差、初始MSE、最大迭代步长和控制函数。

2) 用MLMC和MSE ${\epsilon }^{\left(i\right)}$ 去估算梯度 ${J^{\prime }}_N\left(u^{\left(i-1\right)}\right)$ 。

3) 用梯度投影(3.4)去更新 $u^{\left(i\right)}$ 。

4) 如果 $\left|\Vert u^{\left(i\right)}-u^{\left(i-1\right)}\Vert \right|\le \tau$ ，返回 $u^{\left(i\right)}$ ，开始程序下一步；否则结束。

5) 如果有 ${\epsilon }^{\left(i\right)}>\eta \left|\Vert u^{\left(i\right)}-u^{\left(i-1\right)}\Vert \right|$ 或者 ${\epsilon }^{\left(i\right)}<{\eta }^2\left|\Vert u^{\left(i\right)}-u^{\left(i-1\right)}\Vert \right|$ 成立，决定下一步的MSE ${\epsilon }^{\left(i\text{+1}\right)}=\eta \left|\Vert u^{\left(i\right)}-u^{\left(i-1\right)}\Vert \right|$ 或者 ${\epsilon }^{\left(i\text{+1}\right)}=\text{max}\left(\tau ,\eta \left|\Vert u^{\left(i\right)}-u^{\left(i-1\right)}\Vert \right|\right)$ 。如果条件不成立，则 ${\epsilon }^{\left(i\text{+1}\right)}={\epsilon }^{\left(i\right)}$ 。

其中 $u,v$ 是2个服从某种确定分布的随机向量，h是离散步长。 $\left|\Vert \cdot \Vert \right|$ 是通过这样定义的内积产生的， ${\left(u,v\right)}_{\Omega ,\left[0,T\right]}={\displaystyle {\int }_{\Omega }h{u}^{\text{T}}}v\text{d}P\left(\Omega \right)$ 。

如果通过上述算法产生的梯度 ${J^{\prime }}_N\left(u^{\left(i\right)}\right)$ ，定理3.1的条件都成立，当 $0<\delta <1$ 时， ${\rho }_i$ 由这 $0<1-2c{\rho }_i+\left(1+2C\right){\rho }_i^2\le {\delta }^2$ 决定。我们有(3.4)的基于多水平蒙特卡洛估计是收敛的，于是有

$\left|\Vert u^{\ast }-u^{\left(i\right)}\Vert \right|\to 0\left(i\to \infty \right)$ (4.5)

证明：根据(3.8) (3.9)和定理3.1，有

$\left|\Vert J^{\prime }\left(u\right)-J^{\prime }\left(v\right)\Vert \right|\le C\left|\Vert u-v\Vert \right|$ (4.6)

${\left(J^{\prime }\left(u\right)-J^{\prime }\left(v\right),u-v\right)}_{\Omega ,\left[0,T\right]}\le c\left|\Vert u-v\Vert \right|\forall u,v\in R^N$ (4.7)

可以得到

${\left|\Vert u^{\left(i+1\right)}-u^{\ast }\Vert \right|}^2\le \delta {\left|\Vert u^{\left(i\right)}-u^{\ast }\Vert \right|}^2+{\left|\Vert J^{\prime }\left(u^{\left(i\right)}\right)-{J^{\prime }}_N\left(u^{\left(i\right)}\right)\Vert \right|}^2,$ (4.8)

其中 $0<\delta <1$ ， $u^{\ast }$ 是(2.2)和(2.3)的精确解。根据上述算法有

$\left|\Vert J^{\prime }\left(u^{\left(i\right)}\right)-{J^{\prime }}_N\left(u^{\left(i\right)}\right)\Vert \right|\le \left|\Vert u^{\left(i\right)}-u^{\left(i-1\right)}\Vert \right|$ (4.9)

我们根据定理有 $i\to \infty$ 时， $\left|\Vert u^{\left(i\right)}-u^{\left(i-1\right)}\Vert \right|\to 0$ 。

因此可得

$\left|\Vert u^{\ast }-u^{\left(i\right)}\Vert \right|\to 0\left(i\to \infty \right)$ .

5. 金融上的应用与数值模拟

随机最优控制问题在工程、金融和金融领域得到了广泛的应用，梯度投影算法下的多水平蒙特卡洛方法可以用在金融问题上，求解金融问题。接下来根据具体问题，实例分析投影梯度下的蒙特卡洛方法的应用，并进行数值模拟。

例子：这个例子来自于金融中的期权定价问题，根据B-S模型 [9] [10] 可得

$\text{d}{y}_t=u\left(t\right)y\left(t\right)\text{d}t+\sigma y\left(t\right)\text{d}{W}_t,y\left(0\right)=y_0$ (5.1)

可将期权定价转化为下面的优化控制问题：

$\underset{u\in L^2\left(0,T\right)}{\min }J\left(u\right)=\frac{1}{2}{\displaystyle {\int }_0^{\text{T}}E\left[{\left(y-y_d\right)}^2\right]\text{d}t}+\frac{1}{2}{\displaystyle {\int }_0^{\text{T}}{u}^2\text{d}t,}$ (5.2)

其中 $\sigma$ 为常数。对这个问题求优化问题可以转化为求解下面的优化控制问题：

$\{\begin{array}{l}\text{d}y=uy\text{d}t+\sigma y\text{d}{W}_t,y\left(0\right)=y_0,\\ -\text{d}p=\left(y-y_d+up-p{\sigma }^2\right)\text{d}t+\sigma p\text{d}{W}_t,p\left(T\right)=0,\\ u=-E\left[py\right]\end{array}$ (5.3)

最优控制问题的精确解如下：

$u^{\ast }=\frac{T-t}{\frac{1}{y_0}-Tt+\frac{1}{2}{t}^2},y_d=\frac{{\text{e}}^{{\sigma }^{2}t}-{\left(T-t\right)}^2}{\frac{1}{y_0}-Tt+\frac{1}{2}{t}^2}+1$ (5.4)

为了用投影梯度算法下的多水平蒙特卡洛法求解优化控制问题，我们设置以下参数进行数值模拟实验：

$u^{\left(0\right)}=1,T=1,y_0=1,\sigma =0.05,0.1,0.2,\tau ={10}^{-4},{\epsilon }^{\left(0\right)}=5e-2,\eta =0.5$ 。

当MLMC方法被用于估算梯度时，参数 ${\alpha }_y,{\alpha }_p,{\beta }_y,{\beta }_p,{\gamma }_1,{\gamma }_2$ 需要被模拟，由于欧拉算法用于数值模拟伴随方程和状态方程，相关参数已经理论上确定了。因此

${\alpha }_y={\alpha }_p=1,{\beta }_y={\beta }_p=0.5,{\gamma }_1={\gamma }_2=1$ .

数值实验结果见图1。当 $\sigma =0.2$ 时，实验计算过程的部分信息，见表。

6. 结论

从图1可以知道：当 $\sigma =0.05,0.1,0.2$ 时，这种算法时收敛的，最终的误差 $\Vert u^{\left(i\right)}-u^{\ast }\Vert$ 分别是 $6.4\text{e}-4,1.2\text{e}-3,9.9\text{e}-3$ 。当 $\sigma =0.2$ 时迭代的总时间为 $2.4\text{e}+4$ 秒。因为我们选择负梯度方向，当 $\tau$ 逐步增加时，图中最后一排体现出误差缓慢的增加 $\left(\sigma =0.1,\sigma =0.2\right)$ 。

从表1可以知道：此算法的精度高，用时也比较少，梯度投影下的MLMC能高效准确的数值模拟最优控制问题，为金融中的期权定价问题提供一种好的数值计算方案。

表1. 实验计算过程的部分信息(σ = 0.2)

图1. 数值实验结果(h = 1/64)

本文讨论了梯度投影下的MLMC的随机最优问题，并将之用于金融问题的实例分析中，数值实验体现了这种方法的高效性和准确性，为分析金融问题提供了一种好的手段。这种方法为求解高维随机最优控制问题提供一种可能，也可以用在随机倒向微分方程上面，也可以求解高维下的期权定价问题。

文章引用

雷 鹏. 一种基于多水平蒙特卡罗方法的最优控制在金融上的应用
Application of Optimal Control Based on Multilevel Monte Carlo Method in Finance[J]. 运筹与模糊学, 2023, 13(03): 1610-1617. https://doi.org/10.12677/ORF.2023.133162

参考文献