一个寻找 SDE 解析解的经验方法——借助 SymPy 计算机代数系统

Posted 2020-11-26 量化与固收札记

目录

• 一个寻找 SDE 解析解的经验方法——借助 SymPy 计算机代数系统
  • Ito 公式与转换
  • 猜测 \(f\) 的形式
  • 若干案例
    • 案例一：几何布朗运动
    • 案例二
    • 案例三
    • 案例四
    • 案例五：随机 Gompertzian 模型
    • 案例六
    • 案例七：Log Mean-Reverting 模型
    • 案例八：特定参数的 Cox Ingersoll Ross 模型
  • 参考文献
  • 附录：SymPy 代码

一个寻找 SDE 解析解的经验方法——借助 SymPy 计算机代数系统

摘要：本文记述了一种找到 SDE 解析解的经验方法，并附带了辅助符号计算的 SymPy 代码。

Ito 公式与转换

一维 SDE 的形式如下：

\[ d X_t = \nu(t, X_t) dt + \mu(t,X_t)d B_t \]

经验解法的核心是找到一个非平凡函数 \(f(t,x)\)，使得 \(Y_t = f(t, X_t)\) 的解析解可以轻松获得，然后用 \(f\) 的逆变换得到 \(X_t\) 的解析解。

应用 Ito 公式，得到 \(Y_t\) 的微分形式：

\[ \begin{aligned} dY_t &= d f(t,X_t)\ &= \left(\frac{\partial f}{\partial t} + \frac{\partial f}{\partial x}\nu + \frac{1}{2}\frac{\partial^ 2 f}{\partial x^2}\mu^2 \right)dt + \frac{\partial f}{\partial x}\mu d B_t\ &= P_1 dt + P_2 d B_t \end{aligned} \]

若要 \(Y_t\) 的解析解可以轻松得到，可以要求 \(\frac{\partial P_1}{\partial x} = 0\) 并且 \(\frac{\partial P_2}{\partial x} = 0\)，即要求 \(P_1\) 和 \(P_2\) 只是 \(t\) 的函数：

\[ \begin{aligned} &\frac{\partial^ 2 f}{\partial t \partial x} + \frac{\partial^ 2 f}{\partial x^2}\nu + \frac{\partial f}{\partial x}\frac{\partial \nu}{\partial x} + \frac{1}{2}\left(\frac{\partial^ 3 f}{\partial x^3}\mu^2 + 2\frac{\partial^ 2 f}{\partial x^2}\frac{\partial \mu}{\partial x}\mu \right) = 0 \&\frac{\partial^ 2 f}{\partial x^2}\mu + \frac{\partial f}{\partial x} \frac{\partial \mu}{\partial x}= 0 \end{aligned} \]

此时可以称 \(Y_t\) 是“简单”SDE。

至此，寻找解析解的过程转换成了寻找一个非平凡函数 \(f(t,x)\)，满足上述两个偏微分方程。

猜测 \(f\) 的形式

从最简单的形式入手，猜测 \(f(t,x)\) 符合乘法形式，即

\[ f(t, x) = F(t)G(x) \]

那么，偏微分方程组简化为：

\[ \begin{aligned} & \frac{d F}{d t}\frac{d G}{d x} + F\left(\frac{d^2 G}{d x^2}\nu + \frac{d G}{d x}\frac{\partial \nu}{\partial x} + \frac{1}{2}\frac{d^3 G}{d x^3}\mu^2 + \frac{d^2 G}{d x^2}\frac{\partial \mu}{\partial x}\mu \right) = 0 \& F\left(\frac{d^2 G}{d x^2}\mu + \frac{d G}{d x} \frac{\partial \mu}{\partial x}\right)= 0 \end{aligned} \]

从直觉上看，突破口在第二个等式上，从第二个等式先解出 \(G\)，进而解出 \(F\)。

若干案例

案例一：几何布朗运动

对于几何布朗运动 \(d X_t = r(t) X_t dt + \sigma(t) X_t dB_t\) 而言，

\[ \begin{cases} \mu(t,x) = \sigma(t) x\\nu(t, x) = r(t) x \end{cases} \]

代入到方程组中得到

\[ \begin{aligned} &F \left(r x \frac{d^{2}G}{d x^{2}} + r \frac{dG}{d x} + \frac{\sigma^{2}}{2} x^{2} \frac{d^{3}G}{d x^{3}} + \sigma^{2} x \frac{d^{2}G}{d x^{2}}\right) + \frac{dF}{d t} \frac{dG}{d x}=0\&F \left(\sigma x \frac{d^{2}G}{d x^{2}} + \sigma \frac{dG}{d x}\right)=0 \end{aligned} \]

\(f(t,x)\) 的一个非平凡解是

\[ \begin{aligned} &f(t, x) = \ln x\&F(t)=1, G(x) = \ln x \end{aligned} \]
那么

\[ \begin{aligned} dY_t &= d \ln(X_t)\ &= (r - \frac{1}{2}\sigma^2 )dt + \sigma d B_t\Y_t &= \int_0^t r(s) - \frac{1}{2}\sigma^2(s) ds + \int_0^t \sigma(s) dB_s + C\X_t &= e^{\int_0^t r(s) - \frac{1}{2}\sigma^2(s) ds + \int_0^t \sigma(s) dB_s + C} \end{aligned} \]

案例二

对于 \(d X_t = \frac{3}{4}t^2X_t^2 dt + tX_t^{3/2} dB_t\) 来说（文献【1】），

\[ \begin{cases} \mu(t,x) = t x^{3/2}\\nu(t, x) = \frac{3}{4}t^2x^2 \end{cases} \]

代入到方程组中得到

\[ \begin{aligned} &F \left(\frac{t^{2}}{2} x^{3} \frac{d^{3}G}{d x^{3}} + \frac{9}{4} t^{2} x^{2} \frac{d^{2}G}{d x^{2}} + \frac{3}{2} t^{2} x \frac{dG}{d x}\right) + \frac{dF}{d t} \frac{dG}{d x}=0\&F \left(t x^{\frac{3}{2}} \frac{d^{2}G}{d x^{2}} + \frac{3}{2} t \sqrt{x} \frac{dG}{d x}\right)=0 \end{aligned} \]

\(f(t,x)\) 的一个非平凡解是

\[ \begin{aligned} &f(t, x) = x^{-1/2}\&F(t)=1, G(x) = x^{-1/2} \end{aligned} \]
那么

\[ \begin{aligned} dY_t &= d X_t^{-1/2}\ &= - \frac{1}{2} t d B_t\Y_t &= - \frac{1}{2} \int_0^t sdB_s + C\X_t &= \frac{1}{\left(-\frac{1}{2}\int_0^t sdB_s + C\right)^2} \end{aligned} \]

案例三

对于 \(d X_t = \frac{1}{2}(c^2(t)rX^{2r-1} - c^2(t)X^{r})dt + c^2(t)X^{r} dB_t, (r\ne1)\) 来说（文献【1】），

\[ \begin{cases} \mu(t,x) = c^2(t)x^{r}\\nu(t, x) = \frac{1}{2}(c^2(t)rx^{2r-1} - c^2(t)x^{r}) \end{cases} \]

代入到方程组中得到

\[ \begin{aligned} &F \left(c^{2} r x^{2 r - 1} \frac{d^{2}G}{d x^{2}} + \frac{c^{2} r}{2 x} \left(- x^{r} + x^{2 r - 1} \left(2 r - 1\right)\right) \frac{dG}{d x} + \frac{c^{2}}{2} x^{2 r} \frac{d^{3}G}{d x^{3}} + \frac{c^{2}}{2} \left(r x^{2 r - 1} - x^{r}\right) \frac{d^{2}G}{d x^{2}}\right) + \frac{dF}{d t} \frac{dG}{d x}=0\&F \left(c r x^{r - 1} \frac{dG}{d x} + c x^{r} \frac{d^{2}G}{d x^{2}}\right)=0 \end{aligned} \]

\(f(t,x)\) 的一个非平凡解是

\[ \begin{aligned} &f(t, x) = x^{-r+1}\&F(t)=1, G(x) = x^{-r+1} \end{aligned} \]
那么

\[ \begin{aligned} dY_t &= d X_t^{1-r}\ &= \frac{c^{2}}{2} \left(r - 1\right)dt+ c \left(1 - r\right) d B_t\Y_t &= \int_0^t \frac{c^{2}(s)}{2} \left(r - 1\right) ds + \int_0^t c(s) \left(1 - r\right) d B_s +C\X_t &= \left( \int_0^t \frac{c^{2}(s)}{2} \left(r - 1\right) ds + \int_0^t c(s) \left(1 - r\right) d B_s +C\right)^{\frac{1}{1-r}} \end{aligned} \]

案例四

对于 \(d X_t = X^3 dt + X^2 dB_t\) 来说（文献【1】），

\[ \begin{cases} \mu(t,x) = x^2\\nu(t, x) = x^3 \end{cases} \]

代入到方程组中得到

\[ \begin{aligned} &F \left(\frac{x^{4}}{2} \frac{d^{3}G}{d x^{3}} + 3 x^{3} \frac{d^{2}G}{d x^{2}} + 3 x^{2} \frac{dG}{d x}\right) + \frac{dF}{d t} \frac{dG}{d x}=0\&F \left(x^{2} \frac{d^{2}G}{d x^{2}} + 2 x \frac{dG}{d x}\right)=0 \end{aligned} \]

\(f(t,x)\) 的一个非平凡解是

\[ \begin{aligned} &f(t, x) = x^{-1}\&F(t)=1, G(x) = x^{-1} \end{aligned} \]
那么

\[ \begin{aligned} dY_t &= d X_t^{-1}\ &= 0 dt - 1 d B_t\Y_t &= - B_t +C\X_t &= \frac{1}{- B_t +C} \end{aligned} \]

案例五：随机 Gompertzian 模型

对于 \(d X_t = \left(-b X_t \ln X_t \right) dt + cX_t dB_t\) 来说（文献【2】），

\[ \begin{cases} \mu(t,x) = cx\\nu(t, x) = -bx\ln x \end{cases} \]

代入到方程组中得到

\[ \begin{aligned} &F \left(- b x \ln{\left(x \right)} \frac{d^{2}G}{d x^{2}} - b \left(\ln{\left(x \right)} + 1\right) \frac{dG}{d x} + \frac{c^{2}}{2} x^{2} \frac{d^{3}G}{d x^{3}} + c^{2} x \frac{d^{2}G}{d x^{2}}\right) + \frac{dF}{d t} \frac{dG}{d x}=0\&F \left(c x \frac{d^{2}G}{d x^{2}} + c \frac{dG}{d x}\right)=0 \end{aligned} \]

\(f(t,x)\) 的一个非平凡解是

\[ \begin{aligned} &f(t, x) = e^{bt}\ln x\&F(t)=e^{bt}, G(x) = \ln(x) \end{aligned} \]
那么

\[ \begin{aligned} dY_t &= d (e^{bt}\ln X_t)\ &= -\frac{c^{2}}{2} e^{b t} dt - c e^{b t} d B_t\Y_t &= -\frac{c^{2}}{2b}e^{bt} - c\int_0^t e^{bs} dB_s +C\X_t &= \exp\left(-\frac{c^{2}}{2b} - ce^{-bt}\int_0^t e^{bs} dB_s +Ce^{-bt}\right) \end{aligned} \]

案例六

对于 \(d X_t = \left(\alpha(t)X_t^{\frac{3}{4}} + \frac{3}{8} \beta^2 X_t^{\frac{1}{2}} \right) dt + \beta X_t^{\frac{3}{4}} dB_t\) 来说（文献【3】），

\[ \begin{cases} \mu(t,x) = \beta x^{\frac{3}{4}}\\nu(t, x) = \alpha(t)x^{\frac{3}{4}} + \frac{3}{8} \beta^2 x^{\frac{1}{2}} \end{cases} \]

代入到方程组中得到

\[ \begin{aligned} &F \left(\frac{\beta^{2}}{2} x^{\frac{3}{2}} \frac{d^{3}G}{d x^{3}} + \frac{3}{4} \beta^{2} \sqrt{x} \frac{d^{2}G}{d x^{2}} + \left(\frac{3 \alpha}{4 \sqrt[4]{x}} + \frac{3 \beta^{2}}{16 \sqrt{x}}\right) \frac{dG}{d x} + \left(\alpha x^{\frac{3}{4}} + \frac{3}{8} \beta^{2} \sqrt{x}\right) \frac{d^{2}G}{d x^{2}}\right) + \frac{dF}{d t} \frac{dG}{d x}=0\&F \left(\beta x^{\frac{3}{4}} \frac{d^{2}G}{d x^{2}} + \frac{3 \beta}{4 \sqrt[4]{x}} \frac{dG}{d x}\right)=0 \end{aligned} \]

\(f(t,x)\) 的一个非平凡解是

\[ \begin{aligned} &f(t, x) = x^{\frac{1}{4}}\&F(t)=1, G(x) = x^{\frac{1}{4}} \end{aligned} \]
那么

\[ \begin{aligned} dY_t &= d (X_t^{\frac{1}{4}})\ &= \frac{\alpha}{4} dt + \frac{\beta}{4} d B_t\Y_t &= \int_0^t \frac{1}{4}\alpha(s) ds+ \frac{\beta}{4} B_t +C\X_t &= \left(\int_0^t \frac{1}{4}\alpha(s) ds+ \frac{\beta}{4} B_t +C \right)^4 \end{aligned} \]

案例七：Log Mean-Reverting 模型

对于 \(d X_t = \eta X_t(\theta(t) - \ln X_t) dt + \rho X_t dB_t\) 来说（文献【3】），

\[ \begin{cases} \mu(t,x) = \rho x\\nu(t, x) = \eta x(\theta(t) - \ln x) \end{cases} \]

代入到方程组中得到

\[ \begin{aligned} &F \left(\eta x \left(\theta - \ln{\left(x \right)}\right) \frac{d^{2}G}{d x^{2}} + \eta \left(\theta - \ln{\left(x \right)} - 1\right) \frac{dG}{d x} + \frac{\rho^{2}}{2} x^{2} \frac{d^{3}G}{d x^{3}} + \rho^{2} x \frac{d^{2}G}{d x^{2}}\right) + \frac{dF}{d t} \frac{dG}{d x}=0\&F \left(\rho x \frac{d^{2}G}{d x^{2}} + \rho \frac{dG}{d x}\right)=0 \end{aligned} \]

\(f(t,x)\) 的一个非平凡解是

\[ \begin{aligned} &f(t, x) = e^{\eta t} \ln x\&F(t)=e^{\eta t}, G(x) = \ln x \end{aligned} \]
那么

\[ \begin{aligned} dY_t &= d (e^{\eta t} \ln X_t)\ &= \left(\eta \theta - \frac{\rho^{2}}{2}\right) e^{\eta t} dt + \rho e^{\eta t} d B_t\Y_t &= \int_0^t \left(\eta \theta(s) - \frac{\rho^{2}}{2}\right) e^{\eta s} ds + \int_0^t \rho e^{\eta s} B_s +C\X_t &= \exp \left( e^{-\eta t}\int_0^t \left(\eta \theta(s) - \frac{\rho^{2}}{2}\right) e^{\eta s} ds + e^{-\eta t}\int_0^t \rho e^{\eta s} B_s +Ce^{-\eta t} \right) \end{aligned} \]

案例八：特定参数的 Cox Ingersoll Ross 模型

对于 \(d X_t = \alpha (\beta - X_t) dt + \sigma X_t^{\frac{1}{2}} dB_t\) 来说（文献【3】），

\[ \begin{cases} \mu(t,x) = \sigma x^{\frac{1}{2}} \\nu(t, x) = \alpha (\beta - x) \end{cases} \]

代入到方程组中得到

\[ \begin{aligned} &F \left(\alpha \left(\beta - x\right) \frac{d^{2}G}{d x^{2}} - \alpha \frac{dG}{d x} + \frac{x}{2} \sigma^{2} \frac{d^{3}G}{d x^{3}} + \frac{\sigma^{2}}{2} \frac{d^{2}G}{d x^{2}}\right) + \frac{dF}{d t} \frac{dG}{d x}=0\&F \left(\sigma \sqrt{x} \frac{d^{2}G}{d x^{2}} + \frac{\sigma}{2 \sqrt{x}} \frac{dG}{d x}\right) \end{aligned} \]

\(G(x)\) 的一个非平凡解是 \(\sqrt x\)，把 \(G\) 代入到第一个等式得到：

\[ 8 x \frac{dF}{d t} - \left(4 \alpha x + 4 \alpha \beta - \sigma^{2}\right) F=0 \]
如果 \(4\alpha \beta = \sigma^2\)，那么 \(F(t)\) 的一个非平凡解是 \(e^{\frac{\alpha}{2} t}\)，此时
\[ \begin{aligned} dY_t &= d (e^{\frac{\alpha}{2} t} \sqrt X_t)\ &= 0 dt + \frac{\sigma}{2} e^{\frac{\alpha}{2} t} d B_t\Y_t &= \int_0^t \frac{\sigma}{2} e^{\frac{\alpha}{2} s} B_s +C\X_t &= e^{-\alpha t}\left( \int_0^t \frac{\sigma}{2} e^{\frac{\alpha}{2} s} B_s +C \right)^2 \end{aligned} \]

因为这个特定参数的 CIR 模型存在解析解，它也许会成为金融工程计算中一个不错的控制变量。

参考文献

1. Analytical solutions for stochastic differential equations via Martingale process
2. Exact Solutions of Stochastic Differential Equations
3. Exact Solvability of Stochastic Differential Equations Driven Finite Activit Levy Processes

附录：SymPy 代码

import sympy as sp
from sympy.abc import alpha, beta, eta, theta, rho, sigma, b, c, F, m, x, r, t, G

one = sp.Integer(1)
two = sp.Integer(2)
three = sp.Integer(3)
four = sp.Integer(4)
eight = sp.Integer(8)

dG = sp.Derivative(G, x)
dG2 = sp.Derivative(G, x, 2)
dG3 = sp.Derivative(G, x, 3)
dG4 = sp.Derivative(G, x, 4)
dF = sp.Derivative(F, t)
dF2 = sp.Derivative(F, t, 2)

# case 1
# mu = sigma * x
# nu = r * x

# case 2
# mu = t * x ** (three / two)
# nu = three / four * t ** 2 * x ** 2

# case 3
# mu = c * x ** r
# nu = one / two * (c ** 2 * r * x ** (2 * r - 1) - c ** 2 * x ** r)

# case 4
# mu = x ** 2
# nu = x ** 3

# case 5
# mu = c * x
# nu = -b * x * sp.ln(x)

# case 6
# mu = beta * x ** (three / four)
# nu = alpha * x ** (three / four) + three / eight * beta ** 2 * x ** (one / two)

# case 7
# mu = rho * x
# nu = eta * x * (theta - sp.ln(x))

# case 8
mu = sigma * x ** (one / two)
nu = alpha * (beta - x)

# 方程组

dMu = mu.diff(x)
dMu2 = mu.diff(x, 2)
dNu = nu.diff(x)

eq1 = dF * dG + F * (dG2 * nu + dG * dNu + one / two * dG3 * mu ** 2 + dG2 * dMu * mu)
eq2 = F * (dG2 * mu + dG * dMu)

print(sp.latex(
    sp.powsimp(eq1),
    long_frac_ratio=1,
    ln_notation=True))
print(sp.latex(
    sp.powsimp(eq2),
    long_frac_ratio=1,
    ln_notation=True))

# 求解 G

Gf = sp.symbols('G', cls=sp.Function)
de = sp.Eq(Gf(x).diff(x) * dMu + Gf(x).diff(x, 2) * mu, 0)

solveG = sp.dsolve(de)

print(sp.latex(
    solveG,
    long_frac_ratio=1,
    ln_notation=True))

# 化简关于 F 的微分方程

f = sp.symbols('F', cls=sp.Function)
g = sp.sqrt(x)
dg = g.diff(x)
dg2 = g.diff(x, 2)
dg3 = g.diff(x, 3)

sim_eq1 = f(t).diff(t) * dg + f(t) * (dg2 * nu + dg * dNu + one / two * dg3 * mu ** 2 + dg2 * dMu * mu)

print(sp.latex(
    sp.simplify(sim_eq1),
    long_frac_ratio=1,
    ln_notation=True))

# 计算 P1 和 P2

Ff = sp.sqrt(x)

p1 = Ff.diff(t) + Ff.diff(x) * nu + one / two * Ff.diff(x, 2) * mu ** 2
p2 = Ff.diff(x) * mu

print(sp.latex(
    sp.simplify(p1),
    long_frac_ratio=1))
print(sp.latex(
    sp.simplify(p2),
    long_frac_ratio=1))