随机微分方程与 Ito Lemma 的关系

一、随机微分方程的由来
二、随机微分方程的求解
三、随机积分的仿真
四、小结

一、随机微分方程的由来

随机微分方程是我们研究随机现象的重要数学工具，我们可以通过一个简单的例子来进行说明。

例1（The Vasicek interest rate model）

VIR 模型是借贷利率的标准模型，它描述的是瞬时利率（instantaneous interest rate，rtr_trt）在均值( μ\muμ )附近上下波动的现象，该模型的数学描述是：
drt=c[μ−rt]dt+σdBt,t∈[0,T](1)dr_t=c[\mu-r_t]dt+\sigma dB_t,\quad t\in [0,T]\tag{1} drt=c[μ−rt]dt+σdBt,t∈[0,T](1)
其中 [μ−rt][\mu-r_t][μ−rt] 表示 rtr_trt 偏离 μ\muμ 的程度，c 是反映回归速度的一个参数，而 BtB_tBt 表示标准布朗运动，σ\sigmaσ 反映的是随机波动的波动性，该数学模型很好地反映了：rtr_trt reverts to the mean μ\muμ. （1）也可以表示为积分形式：
rt=r0+c∫0t[μ−rs]ds+σ∫0tdBs(2)r_t=r_0 + c\int_0^t[\mu-r_s]ds+\sigma\int_0^t dB_s\tag{2} rt=r0+c∫0t[μ−rs]ds+σ∫0tdBs(2)
这是一个典型的线性随机微分方程的求解过程，其目的是要找到满足（1）或（2）的随机过程。

二、随机微分方程的求解

随机微分方程的求解，总体上有两种方法：

1、公式法

2、数值求解方法

本文主要叙述的是“公式法”，我们通过该方法可以得到微分方程的显示解，即一个描述微分方程解的函数表达式。在求解的过程中，我们主要要利用的工具就是 Ito Lemma，它可以将随机微分方程转化为非随机微分方程，让求解过程绕开 Ito 积分难题。那什么是 Ito 积分呢？如下就是它的表达式：
I(t)=∫0tCsdBs若 Cs 是布朗过程 Bs，则有：Sn=∑i=1nBti−1ΔiBI(t)=lim⁡n→∞Sn=12(Bt2−t)(3)I(t)=\int_0^t C_sdB_s \tag{3}\\ \text{若 Cs 是布朗过程 Bs，则有：}\\ S_n = \sum_{i=1}^n B_{t_{i-1}}\Delta_i B \\ I(t) = \lim_{n\rightarrow\infin} S_n = \frac 12(B_t^2-t) I(t)=∫0tCsdBs若 Cs 是布朗过程 Bs，则有：Sn=i=1∑nBti−1ΔiBI(t)=n→∞limSn=21(Bt2−t)(3)
Ito 积分的结果实质上是一个随机过程，该积分有四个性质：

1、The Ito stochastic integral has expectation zero.

2、The Ito stochastic integral satisfies the isometry property:
E(∫0tCsdBs)2=∫0tECs2ds,t∈[0,T](4)E\left( \int_0^t C_s dB_s \right)^2 = \int_0^t EC_s^2 ds,\quad t\in [0,T]\tag{4} E(∫0tCsdBs)2=∫0tECs2ds,t∈[0,T](4)
3、The Ito stochastic integral is linear, and it is linear on adjacent intervals

4、The process I(C)I(C)I(C) has continuous sample paths.

若积分的对象 CsC_sCs 不是布朗过程，则该积分并不容易求出来，要解决随机微分方程求Ito积分的困难，就需要求助于 Ito Lemma，即 Ito 引理，该引理有4个形式：

1、简单版本
f(Bt)−f(Bs)=∫stf′(Bx)dBx+12∫stf′′(Bx)dBx(5)f(B_t)-f(B_s)=\int_s^t f'(B_x)dB_x + \frac 12\int_s^tf''(B_x)dB_x\tag{5} f(Bt)−f(Bs)=∫stf′(Bx)dBx+21∫stf′′(Bx)dBx(5)
2、扩展1
f(t,Bt)−f(s,Bs)=∫st[f1(x,Bx)+12f22(x,Bx)]dx+∫stf2(x,Bx)dBx(6)f(t,B_t)-f(s,B_s)=\int_s^t \left[f_1(x,B_x)+\frac 12 f_{22}(x,B_x)\right] dx+\int_s^t f_2(x,B_x)dB_x\tag{6} f(t,Bt)−f(s,Bs)=∫st[f1(x,Bx)+21f22(x,Bx)]dx+∫stf2(x,Bx)dBx(6)
3、扩展2

Ito 过程：
Xt=X0+∫0tAs(1)ds+∫0tAs(2)dBs(7)X_t = X_0+\int_0^t A_s^{(1)}ds+\int_0^t A_s^{(2)}dB_s\tag{7} Xt=X0+∫0tAs(1)ds+∫0tAs(2)dBs(7)
Ito 扩展2：
f(t,Xt)−f(s,Xs)=∫st[f1+Ay(1)f2+12[Ay(2)]2f22]dy+∫stA2(2)f2dBy(8)f(t,X_t)-f(s,X_s)=\int_s^t\left[f_1+A_y^{(1)}f_2+\frac 12[A_y^{(2)}]^2f_{22}\right]dy+\int_s^t A_2^{(2)}f_2 dB_y\tag{8} f(t,Xt)−f(s,Xs)=∫st[f1+Ay(1)f2+21[Ay(2)]2f22]dy+∫stA2(2)f2dBy(8)
扩展2是最常用的 Lemma 形式。

4、扩展3：
f(t,Xt(1),Xt(2))−f(s,Xs(1),Xs(2))f(t,X_t^{(1)},X_t^{(2)})-f(s,X_s^{(1)},X_s^{(2)}) f(t,Xt(1),Xt(2))−f(s,Xs(1),Xs(2))
我在这里就不列其完整方程了，在【1】中有详细说明。按照这种方式，还可以扩展到更多维的随机过程上。

从上面的扩展2形式来看，Lemma 需要一个二次可导的函数 fff ，此函数往往是我们自己创造的，为了消掉积分中那些难点。现在，我们就用 Lemma 求解例1方程，步骤如下：

1、待解方程
rt=r0+c∫0t[μ−rs]ds+∫0tσdBsSupposert=r0+∫0t[c1(s)rs+c2(s)]ds+∫0t[σ1(s)rs+σ2(s)]dBsWe getc1(s)=−c,c2(s)=cμand σ1(s)=0,σ2(s)=σr_t = r_0+c\int_0^t[\mu - r_s]ds+\int_0^t\sigma dB_s \\ \text{Suppose}\\ r_t = r_0 + \int_0^t[c_1(s)r_s+c_2(s)]ds + \int_0^t[\sigma_1(s)r_s+\sigma_2(s)]dB_s\\ \text{We get}\\ c_1(s)=-c,\quad c_2(s)=c\mu\quad \text{and} \ \sigma_1(s)=0,\quad \sigma_2(s)=\sigma rt=r0+c∫0t[μ−rs]ds+∫0tσdBsSupposert=r0+∫0t[c1(s)rs+c2(s)]ds+∫0t[σ1(s)rs+σ2(s)]dBsWe getc1(s)=−c,c2(s)=cμand σ1(s)=0,σ2(s)=σ
2、造一个方程
Yt=f(t,Xt)=exp{−∫0tc1(s)ds}Xt=y(t)XtSuppose: rs=Xs,y(t)=exp{−∫0tc1(s)ds}y′(t)=−c1(t)y(t)Y_t=f(t,X_t)=exp\left\{ -\int_0^tc_1(s)ds \right\}X_t = y(t)X_t\\ \text{Suppose: } \ r_s = X_s ,\ y(t)=exp\left\{ -\int_0^tc_1(s)ds \right\} \\ y'(t)=-c_1(t)y(t) Yt=f(t,Xt)=exp{−∫0tc1(s)ds}Xt=y(t)XtSuppose: rs=Xs, y(t)=exp{−∫0tc1(s)ds}y′(t)=−c1(t)y(t)
3、应用 Lemma 的扩展2
f1=y′(t)Xt,f2=y(t),f22=0A(1)=c[μ−rs],A(2)=σYt−Y0=f(t,Xt)−f(0,X0)=∫st[f1+Ay(1)f2+12[Ay(2)]2f22]dy+∫stA2(2)f2dBy=∫0t[y′(s)Xs+c[μ−Xs]y(s)+12σ20]ds+∫0tσy(s)dBs=∫0t[−cy(s)Xs+c[μ−Xs]y(s)]ds+∫0tσy(s)dBs=∫0tcμy(s)ds+∫0tσy(s)dBs=cu∫0ty(s)ds+σ∫0ty(s)dBsf_1=y'(t)X_t,\quad f_2=y(t),\quad f_{22} = 0\\ A^{(1)}=c[\mu - r_s],\quad A^{(2)}=\sigma\\ Y_t-Y_0=f(t,X_t)-f(0,X_0)=\int_s^t\left[f_1+A_y^{(1)}f_2+\frac 12[A_y^{(2)}]^2f_{22}\right]dy+\int_s^t A_2^{(2)}f_2 dB_y\\ =\int_0^t\left[y'(s)X_s+ c[\mu - X_s]y(s)+\frac 12\sigma^20\right]ds+\int_0^t \sigma y(s)dB_s\\ =\int_0^t\left[-cy(s)X_s+ c[\mu - X_s]y(s)\right]ds+\int_0^t \sigma y(s)dB_s\\ =\int_0^tc\mu y(s)ds+\int_0^t \sigma y(s)dB_s = cu\int_0^ty(s)ds+\sigma\int_0^ty(s)dBs f1=y′(t)Xt,f2=y(t),f22=0A(1)=c[μ−rs],A(2)=σYt−Y0=f(t,Xt)−f(0,X0)=∫st[f1+Ay(1)f2+21[Ay(2)]2f22]dy+∫stA2(2)f2dBy=∫0t[y′(s)Xs+c[μ−Xs]y(s)+21σ20]ds+∫0tσy(s)dBs=∫0t[−cy(s)Xs+c[μ−Xs]y(s)]ds+∫0tσy(s)dBs=∫0tcμy(s)ds+∫0tσy(s)dBs=cu∫0ty(s)ds+σ∫0ty(s)dBs
这里的消元是关键，通过所构造的函数特殊形式，在代入Lemma时，将积分号中的 XsX_sXs 消除了。代入：c1(t)=−cc_1(t)=-cc1(t)=−c, 有：
y(t)=exp{−∫0tc1(s)ds}=ect∫0tecsds=1c(ect−1)y(t)=exp\left\{ -\int_0^tc_1(s)ds \right \}=e^{ct}\\ \int_0^t e^{cs}ds=\frac 1c(e^{ct}-1) y(t)=exp{−∫0tc1(s)ds}=ect∫0tecsds=c1(ect−1)
代入上式，有：
Yt=Y0+μ(ect−1)+σ∫0tecsdBs=Xtectrt=Xt=r0e−ct+μ(1−e−ct)+σe−ct∫0tecsdBsY_t=Y_0+\mu(e^{ct}-1)+\sigma\int_0^t e^{cs}dB_s=X_te^{ct}\\ r_t=X_t=r_0e^{-ct}+\mu(1-e^{-ct})+\sigma e^{-ct}\int_0^te^{cs}dB_s Yt=Y0+μ(ect−1)+σ∫0tecsdBs=Xtectrt=Xt=r0e−ct+μ(1−e−ct)+σe−ct∫0tecsdBs
这就是例1微分方程的解。

通过 Lemma 解微分方程的求解的过程其实就是通过设计一个函数 f(t,Xt)f(t,X_t)f(t,Xt) ，其代入扩展2 式后，能够消除积分号中的 XsX_sXs ，从而得到解。消元就是用 Lemma 的目的。

再看一个例子：

例2：The Homogeneous Linear Equation:
Xt=X0+∫0tc1(s)Xsds+∫0tσ1(s)XsdBsX_t = X_0+\int_0^t c_1(s)X_sds+\int_0^t\sigma_1(s)X_sdB_s Xt=X0+∫0tc1(s)Xsds+∫0tσ1(s)XsdBs
解：

令：Yt=lnXt=f(t,Xt)Y_t=\mathbb {ln} X_t=f(t,X_t)Yt=lnXt=f(t,Xt)，于是有：
f(t,x)=lnx,f1(t,x)=0,f2(t,x)=x−1,f22(t,x)=−x−2f(t,x)=\mathbb{ln} x,\quad f_1(t,x)=0, \quad f_2(t,x)=x^{-1}, \quad f_{22}(t,x)=-x^{-2} f(t,x)=lnx,f1(t,x)=0,f2(t,x)=x−1,f22(t,x)=−x−2
代入扩展2，有：
Yt−Y0=f(t,Xt)−f(0,X0)=∫0t[f1+Ay(1)f2+12[Ay(2)]2f22]dy+∫0tA2(2)f2dBydYt=[c1(t)−12σ1(t)x2x−2]dt+x−1xσ1(t)dBt=[c1(t)−0.5∗σ1(t)]dt+σ1(t)dBtY_t-Y_0=f(t,X_t)-f(0,X_0)=\int_0^t\left[f_1+A_y^{(1)}f_2+\frac 12[A_y^{(2)}]^2f_{22}\right]dy+\int_0^t A_2^{(2)}f_2 dB_y\\ dY_t = [c_1(t)-\frac 12\sigma_1(t)\color{#F00}{x^2x^{-2}}]dt + \color{#f00}{x^{-1}x}\sigma_1(t) dB_t=[c_1(t)-0.5*\sigma_1(t)]dt + \sigma_1(t) dB_t Yt−Y0=f(t,Xt)−f(0,X0)=∫0t[f1+Ay(1)f2+21[Ay(2)]2f22]dy+∫0tA2(2)f2dBydYt=[c1(t)−21σ1(t)x2x−2]dt+x−1xσ1(t)dBt=[c1(t)−0.5∗σ1(t)]dt+σ1(t)dBt
这里有明显的消元过程。最后有：
Xt=X0exp{∫0t[c1(s)−0.5∗σ1(s)]ds+∫0tσ1(s)dBs}X_t=X_0\mathbb {exp}\left\{ \int_0^t [c_1(s)-0.5*\sigma_1(s)]ds + \int_0^t\sigma_1(s) dB_s\right\} Xt=X0exp{∫0t[c1(s)−0.5∗σ1(s)]ds+∫0tσ1(s)dBs}
上式右边并无包含 XtX_tXt ，皆为已知函数的积分，这便是解。

三、随机积分的仿真

例1 的解为：
rt=Xt=r0e−ct+μ(1−e−ct)+σe−ct∫0tecsdBs(9)r_t=X_t=r_0e^{-ct}+\mu(1-e^{-ct})+\sigma e^{-ct}\int_0^te^{cs}dB_s \tag{9} rt=Xt=r0e−ct+μ(1−e−ct)+σe−ct∫0tecsdBs(9)
令：σ=0.5,r0=5.05,μ=5,c=0.2\sigma=0.5,r_0=5.05,\mu=5,c=0.2σ=0.5,r0=5.05,μ=5,c=0.2

1、生成 Brownian 运动

N = 1000
M = np.linspace(2,400,200).astype('int')
t = np.linspace(0,2*np.pi, 1000)Z = np.random.normal(0,1,400)fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(221)for m in M[:20]:Bt = Z[0]/np.sqrt(2*np.pi) * tfor n in range(1,m):Bt = Bt + Z[n]/np.sqrt(np.pi) * 2 *  np.sin(n*t/2) / nax.plot(t, Bt)ax = fig.add_subplot(222)
Bt = Z[0]/np.sqrt(2*np.pi) * t
for n in range(1,39):Bt = Bt + Z[n]/np.sqrt(np.pi) * 2 *  np.sin(n*t/2) / n
ax.plot(t, Bt)    ax = fig.add_subplot(223)
Bt = Z[0]/np.sqrt(2*np.pi) * t
for n in range(1,99):Bt = Bt + Z[n]/np.sqrt(np.pi) * 2 *  np.sin(n*t/2) / n
ax.plot(t, Bt) ax = fig.add_subplot(224)
Bt = Z[0]/np.sqrt(2*np.pi) * t
for n in range(1,399):Bt = Bt + Z[n]/np.sqrt(np.pi) * 2 *  np.sin(n*t/2) / n
ax.plot(t, Bt)
plt.show()

sigma = 0.01
r0 = 5.05
mu = 5def riemann_stieltjes_sum(f, bs):rs_sum = np.zeros(f.shape[0])for i in range(1,rs_sum.shape[0]):rs_sum[i] = rs_sum[i-1] + (bs[i]-bs[i-1])*f[i]return rs_sumt = np.asarray([i/10.0 for i in range(0,1000)])
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)c = 1
rt = r0*np.exp(-1*c*t) + mu*(1-np.exp(-1*c*t)) + sigma*np.exp(-1*c*t)*riemann_stieltjes_sum(np.exp(c*t), Bt)
ax.plot(t,rt, label='c=1')c = 0.1
rt = r0*np.exp(-1*c*t) + mu*(1-np.exp(-1*c*t)) + sigma*np.exp(-1*c*t)*riemann_stieltjes_sum(np.exp(c*t), Bt)
ax.plot(t,rt,label='c=0.1')
plt.legend()
plt.show()

从图中可以看到的是，瞬时利率( rtr_trt )围绕均值 μ\muμ 上下波动，反映回归速度的参数 c 控制着回到均值的时间。

四、小结

随机微分方程的表达形式有两种：积分形式和微分形式。其实，微分形式是不存在的，因为布朗过程的微分是不存在的，它只是积分形式的简化书写形式，其简化形式，比较利于我们去理解模型，因而，我们不用随机积分方程来命名，而喜欢用随机微分方程来对此概念命名。

随机微分方程是我们对现实世界随机过程进行建模的利器，一般的随机微分方程是没有显示解的，只要特殊的、简单的线性微分方程有显式解，但解决任何问题都遵循从简单到复杂，解决了简单的线性随机微分方程，可以为我们研究更复杂的模型提供解决思路和基础。

最后，【1】是一本非常好的关于随机微分方程的参考书，概念清楚，内容逻辑铺垫充分，而且很薄。

参考：

[1] “Elementary Stochastic Calculus with Finance in View” (《随机分析基础》)，Thomas Mikosch, World Scientific, ISBN 978-7-5100-0524-4/O.740, 1998