CQF笔记M1L5仿真和操作随机微分方程

Module 1 Building Blocks of Quant Finance
- Lecture 5 Simulating and Manipulating Stochastic Differential Equations
- - 操作随机微分方程
  - 转移概率密度函数
  - - 前向方程
    - 对数正态随机游走
    - 稳态分布
    - 反向方程
  - 仿真对数正态随机游走
  - 生成相关的随机数

Module 1 Building Blocks of Quant Finance

Lecture 5 Simulating and Manipulating Stochastic Differential Equations

操作随机微分方程

dG=a(G,t)dt+b(G,t)dXdG = a(G,t) \ dt + b(G,t) \ dXdG=a(G,t) dt+b(G,t) dX
称为GGG的随机微分方程或者dGdGdG的随机游走，由两部分组成：

a(G,t)dta(G,t)dta(G,t)dt是确定性部分，dtdtdt的系数称为drift或者growth
b(G,t)dXb(G,t)dXb(G,t)dX是随机部分, dXdXdX的系数称为diffusion或者volatility

资产价格和几何布朗运动

dS=μSdt+σSdX(1)dS=\mu S dt + \sigma S dX \tag{1}dS=μSdt+σSdX(1)

drift a(S,t)=μSa(S,t)=\mu Sa(S,t)=μS
diffusion b(S,t)=σSb(S,t)=\sigma Sb(S,t)=σS

这个随机过程叫几何布朗运动Geometric Brownian Motion(GMB)或指数布朗运动Exponential Brownian Motion(EMB)，广泛用于资产定价

期权价格

df=f(X+dX)−f(x)=dfdXdX+12d2fdX2dt(2)df = f(X+dX)-f(x) = \frac{df}{dX} dX + \frac{1}{2} \frac{d^2f}{dX^2} dt \tag{2}df=f(X+dX)−f(x)=dXdfdX+21dX2d2fdt(2)
其中用到limdt→0dX2=dtlim_{dt \to 0} dX^2 = dtlimdt→0dX2=dt

重写公式(1)(1)(1)为 dSS=μdt+σdX\begin{aligned} \frac{dS}{S} = \mu dt + \sigma dX \end{aligned}SdS=μdt+σdX

dSdSdS反映很短时间间隔dtdtdt内的价格变化, dSS\begin{aligned} \frac{dS}{S} \end{aligned}SdS表示资产收益率
μ\muμ 是资产收益率的平均值
σ\sigmaσ 反映资产收益率的波动性（标准差）
dXdXdX 是布朗运动的增量，服从正态分布dX∼N(0,dt)dX \sim N(0, dt)dX∼N(0,dt)

期权价格是随机过程的函数V(S)V(S)V(S), 进行泰勒展开得到 dV=dVdSdS+12d2VdS2dS2\begin{aligned} dV=\frac{dV}{dS} dS + \frac{1}{2} \frac{d^2V}{dS^2} dS^2 \end{aligned}dV=dSdVdS+21dS2d2VdS2

忽略dt32,dt2dt^\frac{3}{2}, dt^2dt23,dt2等高阶无穷小, 得到dS2=σ2S2dtdS^2=\sigma ^2 S^2 dtdS2=σ2S2dt，带入得到
dV=(μSdVdS+12σ2S2d2VdS2)dt+(σSdVdS)dX(3)\begin{aligned} dV = (\mu S \frac{dV}{dS} + \frac{1}{2} \sigma^2 S^2 \frac{d^2V}{dS^2})dt + (\sigma S \frac{dV}{dS})dX \end{aligned} \tag{3}dV=(μSdSdV+21σ2S2dS2d2V)dt+(σSdSdV)dX(3)
这是一个新的随机微分方程，同样包括确定性部分和随机部分

实例1：对数
由于几何布朗运动的右侧有S，无法直接通过积分求解，定义中间函数V(S)=logSV(S)=logSV(S)=logS，有

dVdS=1Sd2VdS2=−1S2\begin{aligned} & \frac{dV}{dS} = \frac{1}{S} \\ & \frac{d^2V}{dS^2} = - \frac{1}{S^2} \\ \end{aligned}dSdV=S1dS2d2V=−S21

由公式(3)(3)(3)得到
d(logS)=(μ−12σ2)dt+σdXd(logS) = (\mu - \frac{1}{2} \sigma ^2)dt + \sigma dXd(logS)=(μ−21σ2)dt+σdX
两边积分，得到
∫0td(logS)=∫0t(μ−12σ2)dτ+∫0tσdX=(μ−12σ2)t+σ(X(t)−X(0))\begin{aligned} \int_0^t d(logS) &= \int_0^t(\mu - \frac{1}{2} \sigma ^2)d\tau + \int_0^t \sigma dX \\ &= (\mu - \frac{1}{2} \sigma ^2)t + \sigma (X(t) - X(0)) \end{aligned}∫0td(logS)=∫0t(μ−21σ2)dτ+∫0tσdX=(μ−21σ2)t+σ(X(t)−X(0))
假设X(0)=0,S(0)=S0X(0) = 0, S(0) = S_0X(0)=0,S(0)=S0,
S(t)=S0e(μ−12σ2)t+σX(t)=S0e(μ−12σ2)t+σϕt\begin{aligned} S(t) &= S_0 e^{(\mu - \frac{1}{2} \sigma ^2)t + \sigma X(t)} &=S_0 e^{(\mu - \frac{1}{2} \sigma ^2)t + \sigma \phi \sqrt{t}} \end{aligned}S(t)=S0e(μ−21σ2)t+σX(t)=S0e(μ−21σ2)t+σϕt

实例2：均值复归

考察 dr=γ(rˉ−r)dt+σdXdr = \gamma(\bar{r} - r)dt + \sigma dXdr=γ(rˉ−r)dt+σdX

γ\gammaγ 是回复速度
rˉ\bar{r}rˉ 是均值（长期均值）

这个微分方程同样无法直接通过积分求解rrr，令u=r−rˉu = r - \bar{r}u=r−rˉ, 两边同时乘以eγte^{\gamma t}eγt, 使用乘法法则，得到d(ueγt)=σeγtdXd(ue^{\gamma t}) = \sigma e^{\gamma t} dXd(ueγt)=σeγtdX, 两边积分得到
u(t)=u(0)e−γt+σ∫0teγ(s−t)dXs=u(0)e−γt+σ(X(t)−γ∫0tX(s)eγ(s−t)ds)\begin{aligned} u(t) &= u(0)e^{-\gamma t} + \sigma \int_0^t e^{\gamma (s-t)} dX_s \\ &= u(0)e^{-\gamma t} + \sigma (X(t) - \gamma \int_0^t X(s)e^{\gamma (s-t)} ds) \end{aligned}u(t)=u(0)e−γt+σ∫0teγ(s−t)dXs=u(0)e−γt+σ(X(t)−γ∫0tX(s)eγ(s−t)ds)

转移概率密度函数

dy=A(y,t)dt+B(y,t)dXdy = A(y, t)dt + B(y, t)dXdy=A(y,t)dt+B(y,t)dX

注意概率密度函数p(y,t;y′,t′)p(y,t;y',t')p(y,t;y′,t′)的定义
Prob(a<y′<battimet′∣yattimet)=∫abp(y,t;y′,t′)dy′\begin{aligned} Prob(a<y'<b \ at \ time \ t' | y\ at\ time\ t) = \int_a^b p(y,t;y',t')dy' \end{aligned}Prob(a<y′<b at time t′∣y at time t)=∫abp(y,t;y′,t′)dy′

前向方程

已知：当前时刻ttt的值为yyy
求解：未来时刻t′t't′的值为y′y'y′的概率(分布)

通过三叉树思想可以得到前向方程和反向方程, 参考泰勒级数和转移概率密度函数

∂p∂t′=12∂2∂y′2(B(y′,t′)2p)−∂∂y′(A(y′,t′)p)\begin{aligned} \frac{\partial p}{\partial t'} = \frac{1}{2} \frac{\partial ^2}{\partial {y'}^2}(B(y',t')^2 p) - \frac{\partial}{\partial y'}(A(y',t')p) \end{aligned}∂t′∂p=21∂y′2∂2(B(y′,t′)2p)−∂y′∂(A(y′,t′)p)

这个方程称为Fokker–Planck或者forward Kolmogorov方程

对数正态随机游走

dS=μSdt+σSdXdS=\mu S dt + \sigma S dXdS=μSdt+σSdX

A(S,t)=μSA(S,t)=\mu SA(S,t)=μS
B(S,t)=σSB(S,t)=\sigma SB(S,t)=σS

前向方程变为
∂p∂t′=12∂2∂y′2(σ2S′2p)−∂∂y′(μS′p)\begin{aligned} \frac{\partial p}{\partial t'} = \frac{1}{2} \frac{\partial ^2}{\partial {y'}^2}(\sigma ^ 2 {S'}^2 p) - \frac{\partial}{\partial y'}(\mu S' p) \end{aligned}∂t′∂p=21∂y′2∂2(σ2S′2p)−∂y′∂(μS′p)

求解得到
p(S,t;S′,t′)=1σS′2π(t′−t)e−(log(S/S′)+(μ−12σ2)(t′−t))2/2σ2(t′−t)\begin{aligned} p(S, t; S', t') = \frac{1}{\sigma S' \sqrt{2 \pi (t' - t)}} e^{-(log(S/S')+(\mu - \frac{1}{2}\sigma^2)(t'-t))^2/2\sigma ^ 2(t'-t)} \end{aligned}p(S,t;S′,t′)=σS′2π(t′−t)1e−(log(S/S′)+(μ−21σ2)(t′−t))2/2σ2(t′−t)

对数随机游走的概率密度分布

下图是上图表达的是同一个意思

稳态分布

当t′→∞t' \to \inftyt′→∞时, p(y,t;y′,t′)p(y, t; y', t')p(y,t;y′,t′)与初始时刻ttt的状态yyy无关。

利率、通胀和波动率的SDE模型呈现稳态分布
对数正态分布随机游走会变为无界或者退化为0（非稳态）

稳态分布p∞(y′)p_{\infty}(y')p∞(y′)满足微分方程

12d2dy′2(B2p∞)−ddy′(Ap∞)=0\begin{aligned} \frac{1}{2} \frac {d^2}{{dy'}^2}(B^2 p_{\infty}) - \frac{d}{dy'}(Ap_{\infty}) = 0 \end{aligned}21dy′2d2(B2p∞)−dy′d(Ap∞)=0

Vasicek model

dr=γ(rˉ−r)dt+σdXdr = \gamma(\bar{r} - r)dt + \sigma dXdr=γ(rˉ−r)dt+σdX

带入稳态分布的微分方程
12σ2d2p∞dy′2−γddr′((rˉ−r′)p∞)=0\begin{aligned} \frac{1}{2} \sigma ^2 \frac {d^2 p_{\infty}}{{dy'}^2} - \gamma \frac{d}{dr'}((\bar{r}-r')p_{\infty}) = 0 \end{aligned}21σ2dy′2d2p∞−γdr′d((rˉ−r′)p∞)=0

解方程得到
p∞=1σγπe−γ(rˉ−r′)2σ2\begin{aligned} p_{\infty} = \frac{1}{\sigma} \sqrt{\frac{\gamma}{\pi}} e^{-\frac{\gamma(\bar{r} - r')^2}{\sigma ^2}} \end{aligned}p∞=σ1πγe−σ2γ(rˉ−r′)2

说明利率
r∼N(rˉ,σ22γ)r \sim N(\bar{r}, \frac{\sigma^2}{2 \gamma})r∼N(rˉ,2γσ2)

TODO: 疑问: Vasicek模型本身满足稳态分布，还是假设他满足稳态分布

反向方程

backward Kolmogorov equation

∂p∂t+12B(y,t)2∂p∂y2+A(y,t)∂p∂y=0\begin{aligned} \frac{\partial p}{\partial t} + \frac{1}{2} B(y,t)^2 \frac {\partial ^ p} {\partial y ^ 2} + A(y, t) \frac {\partial p} {\partial y} = 0 \end{aligned}∂t∂p+21B(y,t)2∂y2∂p+A(y,t)∂y∂p=0

仿真对数正态随机游走

dS=μSdt+σSdXdS=\mu S dt + \sigma S dXdS=μSdt+σSdX

变换为离散时间
Si+1=Si(1+μδt+σϕδ1/2)\begin{aligned} S_{i+1} = S_i (1+\mu \delta t + \sigma \phi \delta ^ {1/2}) \end{aligned}Si+1=Si(1+μδt+σϕδ1/2)

Euler method 用电子表格软件仿真这个模型

输入参数

初始值 S0S_0S0
时间步长 δt\delta tδt
drift rate μ\muμ
波动率 σ\sigmaσ
步数

标准正态分布发生器 ϕ\phiϕ
在Excel中，

精确但是慢的方法：用函数NORMSINV(RAND())产生标准正态分布的随机数
不精确但是快的方法：(∑i=112\begin{aligned} \sum_{i=1}^{12} \end{aligned}i=1∑12 RAND()) - 6
这里是用均匀分布之和呈现正态分布的特性，只取了12次均匀分布的和，所以不太精确

仿真其他随机游走

基本思路都是先离散化，再用正态分布随机数发生器产生数据

生成相关的随机数

生成不相关的随机数

生成不相关的随机数ϵ1\epsilon_1ϵ1和ϵ2\epsilon_2ϵ2
组合随机数
ϕ1=ϵ1\phi_1 = \epsilon_1ϕ1=ϵ1
ϕ2=ρϵ1+1−ρ2ϵ2\phi_2 = \rho \epsilon_1 + \sqrt{1-\rho^2} \epsilon_2ϕ2=ρϵ1+1−ρ2ϵ2$
ϕ1\phi_1ϕ1和ϕ2\phi_2ϕ2的相关性为ρ\rhoρ

独立正态分布Xi∼N(μi,σi2)Xi \sim N(\mu_i, \sigma_i^2)Xi∼N(μi,σi2)的加权和满足正态分布
∑i=1nωiXi∼N(∑i=1nωiμi,∑i=1nωi2σi2)\begin{aligned} \sum_{i=1}^{n} \omega_i X_i \sim N(\sum_{i=1}^{n} \omega_i \mu_i, \sum_{i=1}^{n} \omega_i ^2 \sigma_i ^ 2) \end{aligned}i=1∑nωiXi∼N(i=1∑nωiμi,i=1∑nωi2σi2)

以上几个结论都可以通过独立分布的相关性/协方差为0证明