一個尋找 SDE 解析解的經驗方法——借助 SymPy 計算機代數系統
目錄
-
- Ito 公式與轉換
- 猜測 \(f\) 的形式
- 若幹案例
- 案例一：幾何布朗運動
- 案例二
- 案例三
- 案例四
- 案例五：随機 Gompertzian 模型
- 案例六
- 案例七：Log Mean-Reverting 模型
- 案例八：特定參數的 Cox Ingersoll Ross 模型
- 參考文獻
- 附錄：SymPy 代碼
摘要：本文記述了一種找到 SDE 解析解的經驗方法，并附帶了輔助符号計算的 SymPy 代碼。
一維 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\) 的函數：
&\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
此時可以稱 \(Y_t\) 是“簡單”SDE。
至此，尋找解析解的過程轉換成了尋找一個非平凡函數 \(f(t,x)\)，滿足上述兩個偏微分方程。
從最簡單的形式入手，猜測 \(f(t,x)\) 符合乘法形式，即
\[f(t, x) = F(t)G(x)
那麼，偏微分方程組簡化為：
& \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
從直覺上看，突破口在第二個等式上，從第二個等式先解出 \(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}
代入到方程組中得到
&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)\) 的一個非平凡解是
&f(t, x) = \ln x\\
&F(t)=1, G(x) = \ln x
那麼
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}
對于 \(d X_t = \frac{3}{4}t^2X_t^2 dt + tX_t^{3/2} dB_t\) 來說（文獻【1】），
\mu(t,x) = t x^{3/2}\\
\nu(t, x) = \frac{3}{4}t^2x^2
&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
&f(t, x) = x^{-1/2}\\
&F(t)=1, G(x) = x^{-1/2}
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}
對于 \(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】），
\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})
&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
&f(t, x) = x^{-r+1}\\
&F(t)=1, G(x) = x^{-r+1}
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}}
對于 \(d X_t = X^3 dt + X^2 dB_t\) 來說（文獻【1】），
\mu(t,x) = x^2\\
\nu(t, x) = x^3
&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
&f(t, x) = x^{-1}\\
&F(t)=1, G(x) = x^{-1}
dY_t &= d X_t^{-1}\\
&= 0 dt - 1 d B_t\\
Y_t &= - B_t +C\\
X_t &= \frac{1}{- B_t +C}
對于 \(d X_t = \left(-b X_t \ln X_t \right) dt + cX_t dB_t\) 來說（文獻【2】），
\mu(t,x) = cx\\
\nu(t, x) = -bx\ln x
&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
&f(t, x) = e^{bt}\ln x\\
&F(t)=e^{bt}, G(x) = \ln(x)
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)
對于 \(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】），
\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}}
&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
&f(t, x) = x^{\frac{1}{4}}\\
&F(t)=1, G(x) = x^{\frac{1}{4}}
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
對于 \(d X_t = \eta X_t(\theta(t) - \ln X_t) dt + \rho X_t dB_t\) 來說（文獻【3】），
\mu(t,x) = \rho x\\
\nu(t, x) = \eta x(\theta(t) - \ln x)
&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
&f(t, x) = e^{\eta t} \ln x\\
&F(t)=e^{\eta t}, G(x) = \ln x
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)
對于 \(d X_t = \alpha (\beta - X_t) dt + \sigma X_t^{\frac{1}{2}} dB_t\) 來說（文獻【3】），
\mu(t,x) = \sigma x^{\frac{1}{2}} \\
\nu(t, x) = \alpha (\beta - x)
&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
\(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}\)，此時
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
因為這個特定參數的 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
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))
★ 持續學習 ★ 堅持創作 ★
![多屬性決策模型詳解（matlab）[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)