一个寻找 SDE 解析解的经验方法——借助 SymPy 计算机代数系统
目录
摘要：本文记述了一种找到 SDE 解析解的经验方法，并附带了辅助符号计算的 SymPy 代码。
Ito 公式与转换
一维 SDE 的形式如下：
[d X_t =
u(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) 的微分形式：
[egin{aligned}
dY_t &= d f(t,X_t)\
&= left(frac{partial f}{partial t} + frac{partial f}{partial x}
u + frac{1}{2}frac{partial^ 2 f}{partial x^2}mu^2
ight)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) 的函数：
[egin{aligned}
&frac{partial^ 2 f}{partial t partial x} + frac{partial^ 2 f}{partial x^2}
u + frac{partial f}{partial x}frac{partial
u}{partial x} +
frac{1}{2}left(frac{partial^ 3 f}{partial x^3}mu^2 + 2frac{partial^ 2 f}{partial x^2}frac{partial mu}{partial x}mu
ight) = 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)
]
那么，偏微分方程组简化为：
[egin{aligned}
& frac{d F}{d t}frac{d G}{d x} + Fleft(frac{d^2 G}{d x^2}
u + frac{d G}{d x}frac{partial
u}{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
ight) = 0 \
& Fleft(frac{d^2 G}{d x^2}mu + frac{d G}{d x} frac{partial mu}{partial x}
ight)= 0
end{aligned}
]
从直觉上看，突破口在第二个等式上，从第二个等式先解出 (G)，进而解出 (F)。
若干案例
案例一：几何布朗运动
对于几何布朗运动 (d X_t = r(t) X_t dt + sigma(t) X_t dB_t) 而言，
[egin{cases}
mu(t,x) = sigma(t) x\
u(t, x) = r(t) x
end{cases}
]
代入到方程组中得到
[egin{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}}
ight) + 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}
ight)=0
end{aligned}
]
(f(t,x)) 的一个非平凡解是
[egin{aligned}
&f(t, x) = ln x\
&F(t)=1, G(x) = ln x
end{aligned}
]
那么
[egin{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】），
[egin{cases}
mu(t,x) = t x^{3/2}\
u(t, x) = frac{3}{4}t^2x^2
end{cases}
]
代入到方程组中得到
[egin{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}
ight) + 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}
ight)=0
end{aligned}
]
(f(t,x)) 的一个非平凡解是
[egin{aligned}
&f(t, x) = x^{-1/2}\
&F(t)=1, G(x) = x^{-1/2}
end{aligned}
]
那么
[egin{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
ight)^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 e1)) 来说（文献【1】），
[egin{cases}
mu(t,x) = c^2(t)x^{r}\
u(t, x) = frac{1}{2}(c^2(t)rx^{2r-1} - c^2(t)x^{r})
end{cases}
]
代入到方程组中得到
[egin{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
ight)
ight) 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}
ight) frac{d^{2}G}{d x^{2}}
ight) + 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}}
ight)=0
end{aligned}
]
(f(t,x)) 的一个非平凡解是
[egin{aligned}
&f(t, x) = x^{-r+1}\
&F(t)=1, G(x) = x^{-r+1}
end{aligned}
]
那么
[egin{aligned}
dY_t &= d X_t^{1-r}\
&= frac{c^{2}}{2} left(r - 1
ight)dt+ c left(1 - r
ight) d B_t\
Y_t &= int_0^t frac{c^{2}(s)}{2} left(r - 1
ight) ds + int_0^t c(s) left(1 - r
ight) d B_s +C\
X_t &= left( int_0^t frac{c^{2}(s)}{2} left(r - 1
ight) ds + int_0^t c(s) left(1 - r
ight) d B_s +C
ight)^{frac{1}{1-r}}
end{aligned}
]
案例四
对于 (d X_t = X^3 dt + X^2 dB_t) 来说（文献【1】），
[egin{cases}
mu(t,x) = x^2\
u(t, x) = x^3
end{cases}
]
代入到方程组中得到
[egin{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}
ight) + 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}
ight)=0
end{aligned}
]
(f(t,x)) 的一个非平凡解是
[egin{aligned}
&f(t, x) = x^{-1}\
&F(t)=1, G(x) = x^{-1}
end{aligned}
]
那么
[egin{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 ight) dt + cX_t dB_t) 来说（文献【2】），
[egin{cases}
mu(t,x) = cx\
u(t, x) = -bxln x
end{cases}
]
代入到方程组中得到
[egin{aligned}
&F left(- b x ln{left(x
ight)} frac{d^{2}G}{d x^{2}} - b left(ln{left(x
ight)} + 1
ight) 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}}
ight) + 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}
ight)=0
end{aligned}
]
(f(t,x)) 的一个非平凡解是
[egin{aligned}
&f(t, x) = e^{bt}ln x\
&F(t)=e^{bt}, G(x) = ln(x)
end{aligned}
]
那么
[egin{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} - cint_0^t e^{bs} dB_s +C\
X_t &= expleft(-frac{c^{2}}{2b} - ce^{-bt}int_0^t e^{bs} dB_s +Ce^{-bt}
ight)
end{aligned}
]
案例六
对于 (d X_t = left(alpha(t)X_t^{frac{3}{4}} + frac{3}{8} eta^2 X_t^{frac{1}{2}} ight) dt + eta X_t^{frac{3}{4}} dB_t) 来说（文献【3】），
[egin{cases}
mu(t,x) = eta x^{frac{3}{4}}\
u(t, x) = alpha(t)x^{frac{3}{4}} + frac{3}{8} eta^2 x^{frac{1}{2}}
end{cases}
]
代入到方程组中得到
[egin{aligned}
&F left(frac{eta^{2}}{2} x^{frac{3}{2}} frac{d^{3}G}{d x^{3}} + frac{3}{4} eta^{2} sqrt{x} frac{d^{2}G}{d x^{2}} + left(frac{3 alpha}{4 sqrt[4]{x}} + frac{3 eta^{2}}{16 sqrt{x}}
ight) frac{dG}{d x} + left(alpha x^{frac{3}{4}} + frac{3}{8} eta^{2} sqrt{x}
ight) frac{d^{2}G}{d x^{2}}
ight) + frac{dF}{d t} frac{dG}{d x}=0\
&F left(eta x^{frac{3}{4}} frac{d^{2}G}{d x^{2}} + frac{3 eta}{4 sqrt[4]{x}} frac{dG}{d x}
ight)=0
end{aligned}
]
(f(t,x)) 的一个非平凡解是
[egin{aligned}
&f(t, x) = x^{frac{1}{4}}\
&F(t)=1, G(x) = x^{frac{1}{4}}
end{aligned}
]
那么
[egin{aligned}
dY_t &= d (X_t^{frac{1}{4}})\
&= frac{alpha}{4} dt + frac{eta}{4} d B_t\
Y_t &= int_0^t frac{1}{4}alpha(s) ds+ frac{eta}{4} B_t +C\
X_t &= left(int_0^t frac{1}{4}alpha(s) ds+ frac{eta}{4} B_t +C
ight)^4
end{aligned}
]
案例七：Log Mean-Reverting 模型
对于 (d X_t = eta X_t( heta(t) - ln X_t) dt + ho X_t dB_t) 来说（文献【3】），
[egin{cases}
mu(t,x) =
ho x\
u(t, x) = eta x( heta(t) - ln x)
end{cases}
]
代入到方程组中得到
[egin{aligned}
&F left(eta x left( heta - ln{left(x
ight)}
ight) frac{d^{2}G}{d x^{2}} + eta left( heta - ln{left(x
ight)} - 1
ight) frac{dG}{d x} + frac{
ho^{2}}{2} x^{2} frac{d^{3}G}{d x^{3}} +
ho^{2} x frac{d^{2}G}{d x^{2}}
ight) + frac{dF}{d t} frac{dG}{d x}=0\
&F left(
ho x frac{d^{2}G}{d x^{2}} +
ho frac{dG}{d x}
ight)=0
end{aligned}
]
(f(t,x)) 的一个非平凡解是
[egin{aligned}
&f(t, x) = e^{eta t} ln x\
&F(t)=e^{eta t}, G(x) = ln x
end{aligned}
]
那么
[egin{aligned}
dY_t &= d (e^{eta t} ln X_t)\
&= left(eta heta - frac{
ho^{2}}{2}
ight) e^{eta t} dt +
ho e^{eta t} d B_t\
Y_t &= int_0^t left(eta heta(s) - frac{
ho^{2}}{2}
ight) e^{eta s} ds + int_0^t
ho e^{eta s} B_s +C\
X_t &= exp left(
e^{-eta t}int_0^t left(eta heta(s) - frac{
ho^{2}}{2}
ight) e^{eta s} ds + e^{-eta t}int_0^t
ho e^{eta s} B_s +Ce^{-eta t}
ight)
end{aligned}
]
案例八：特定参数的 Cox Ingersoll Ross 模型
对于 (d X_t = alpha (eta - X_t) dt + sigma X_t^{frac{1}{2}} dB_t) 来说（文献【3】），
[egin{cases}
mu(t,x) = sigma x^{frac{1}{2}} \
u(t, x) = alpha (eta - x)
end{cases}
]
代入到方程组中得到
[egin{aligned}
&F left(alpha left(eta - x
ight) 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}}
ight) + 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}
ight)=0
end{aligned}
]
(G(x)) 的一个非平凡解是 (sqrt x)，把 (G) 代入到第一个等式得到：
[8 x frac{dF}{d t} - left(4 alpha x + 4 alpha eta - sigma^{2}
ight) F=0
]
如果 (4alpha eta = sigma^2)，那么 (F(t)) 的一个非平凡解是 (e^{frac{alpha}{2} t})，此时
[egin{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
ight)^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 代码
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))