一个寻找 SDE 解析解的经验方法——借助 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)=0 \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 模型存在解析解，它也许会成为金融工程计算中一个不错的控制变量。
参考文献
- Analytical solutions for stochastic differential equations via Martingale process
- Exact Solutions of Stochastic Differential Equations
- Exact Solvability of Stochastic Differential Equations Driven Finite Activit Levy Processes
附录：SymPy 代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
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))
本文由作者按照 CC BY 4.0 进行授权