Stochastic Progress

Brownian Motion 模拟

布朗运动又被称为Wienner过程。

标准布朗运动的数学性质

设连续时间随机过程\(B_t:0\leq t\leq T\)是\([0,T)\)上的标准布朗运动，

• \(B_0=0\)

• 独立增量性：对于有限个时刻\(0\leq t_1< t_2<\dots<t_n\leq T,\)随机变量 \[ B_{t_2}-B_{t_1},B_{t_3}-B_{t_2},\dots,B_{t_n}-B_{t_{n-1}} \] 是独立的。

• 正态性：对任意的\(0\leq s<t<T,B_t-B_s\sim N(0,t-s)\).

设\(B_0=0.\)对\(j=1,2,...,N,\)做以下几步

• \(t_j=t_{j-1}+\Delta t\)
• 产生\(Z_j\sim N(0,1)\)
• \(\Delta B_j=Z_j\sqrt{\Delta t}\)
• \(B_j=B_{j-1}+\Delta B_j\)

set.seed(123)
N <- 100 # 时间[0,T]分割的份数
T <- 1
Delta <- T/N # 时间间隔
B <- numeric(N+1) #初始化0向量
t <- seq(0,T,length=N+1)
for(i in 2:(N+1))
  B[i] <- B[i-1]+rnorm(1)*sqrt(Delta)
plot(t,B,type="l", main="Brownian Motion",ylim = c(-1,1))

假设\(\{X_t,0\leq t\leq T\}\)是关于布朗运动生成的事件流适应的随机过程，满足\(\int_{0}^T\mathbb{E}(X_s^2)ds<+\infty.\)则\(X\)的\(It\hat{o}\)积分定义为： \[ I_t(X)=\int_{0}^TX_sdB_s=\lim_{||\Pi_n||\rightarrow 0}\sum_{i=0}^{n-1} X_{t_i}(B_{t_{i+1}}-B_{t_i}) \] 例如： \[ \begin{align} \int_{t_0}^{t_1}dB_s & =\lim_{n\rightarrow \infty}\sum_{i=0}^{n-1}(B_{t_{i+1}}-B_{t_i})\\ & = B_t-B_{t_0} \end{align} \]

随机过程\(\{X_t,0\leq t\leq T\}\)可以写成如下形式： \[ X_t=X_0+\int_0^t u(s)ds+\int_{0}^tv(s)dB_s \] 其中\(u(t,\omega),v(s,\omega)\)是两个适应过程，且满足 \[ \mathbb{P}\{\int_{0}^T|u(t,\omega)|dt<\infty\}=1\quad and\quad \mathbb{P}\{\int_{0}^Tv^2(t,\omega)dt<\infty\}=1 \] 则\(\{X_t,0\leq t\leq T\}\)被称为\(It\hat{o}\)过程。

随机微分方程（以下简称SDE）

\[ X_t=X_0+\int_{0}^t \mu(s,X_s)ds+\int_{0}^t\sigma(s,X_s)dB_s \] 通常写成微分形式 \[ dX_t = \mu(t,X_t)dt+\sigma(t,X_t)dB_t \] 其中\(\mu(\cdot)\)被称作漂移项，\(\sigma(\cdot)\)被称作扩散项。

下面介绍几个SDE

布朗运动： \[ dX_t=\mu dt+\sigma dB_t \] 几何布朗运动： $$

dX_t = X_t dt+X_tdB_t $$ ## SDE的解 如果SDE的漂移项和扩散项满足下面的条件

• 全局Lipschitz条件：对于所有的\(x,y\in \mathbf{R},t\in[0,T],\)存在常数\(K<\infty\)使得 \[ |\mu(t,x)-\mu(t,y)|+|\sigma(t,x)-\sigma(t,y)|<K|x-y| \]
• 线性增长条件：对于所有的\(x\in\mathbf{R},t\in [0,T],\)存在常数\(C<\infty\)使得 \[ |\mu(t,x)|+|\sigma(t,x)|<C(1+|x|) \]

则SDE存在唯一的连续的强解使得 \[ \mathbb{E}\{\int_0^T |X_t|^2dt\}<\infty. \]

布朗运动

\[ dX_t = \mu dt+\sigma dB_t \] 在给定初值\(X_{t_0}\)的条件下，可以求出方程的解为 \[ X_t=X_0+\mu(t-t_0)+\sigma(B_t-B_{t_0}) \] ### 几何布朗运动

$$

dX_t = X_t dt+X_tdB_t \[ 在给定初值$X_{t_0}$的条件下，可以求出方程的解: \] \[\begin{align} X_t &= X_{t_0}\exp\{(\mu-\frac{1}{2}\sigma^2)(t-t_0)+\sigma(B_t-B_{t_0}) \} \end{align}\] $$

\(It\hat{o}\)引理

假设\(X_t\)满足SDE: \[ dX_t=\mu(t,X_t)dt+\sigma(t,X_t)dB_t. \] \(f(X)\)是\(X\)的函数，则 \[ df(X)=f_X(X)dX+\frac{1}{2}f_{XX}(dX)^2 \] \((dX)^2\)展开时，有如下法则：

\[ dt\times dt=dt\times dB_t=dB_t\times dt=0;dB_t\times dB_t=dt. \] ## \(It\hat{o}\)引理的应用

假设\(X_t\)满足SDE:

\[ dX_t=\mu X_tdt+\sigma X_tdB_t \] 如何解\(X_t?\)

设\(Y_t = \ln X_t,\frac{\partial Y_t}{\partial X_t} = \frac{1}{X_t},\frac{\partial^2 Y_t}{\partial (X_t)^2}=-\frac{1}{X_t^2},\)由\(|t\hat{o}\)引理知： \[ \begin{align} dY_t&=\frac{\partial Y_t}{\partial X_t}dX_t +\frac{1}{2}\frac{\partial^2 Y_t}{\partial X_t^2}(dX_t)^2\\ &=\frac{1}{X_t}(\mu X_tdt+\sigma X_tdB_t)+\frac{1}{2}(-\frac{1}{X_t^2})\sigma^2X_t^2dt\\ &=\mu dt+\sigma dB_t-\frac{1}{2}\sigma^2dt\\ &=(\mu-\frac{1}{2}\sigma^2)dt+\sigma dB_t \end{align} \] 于是可以看出\(Y_t\)是布朗运动，这是因为它服从布朗运动的随机微分方程。

我们选择的\(Y_t\)看似是比较有技巧的,选取的原则是由于\(dX_t\)表达式中有\(X_t\)的一次项，所以我们选择的对\(Y_t\)求导可以消掉一次项。

\(Y_t\)服从布朗运动，于是有 \[ \begin{align} Y_t&=Y_{t_0}+\Big(\mu-\frac{1}{2}\sigma^2\Big)(t-t_0)+\sigma(B_{t}-B_{t_0})\\ X_t&=\exp(Y_t)\\ &=X_{t_0}\exp\Big[\Big(\mu-\frac{1}{2}\sigma^2\Big)(t-t_0)+\sigma(B_{t}-B_{t_0})\Big] \end{align} \]

几何布朗运动的模拟

设定初始值为\(X_0\), 对\(j=1,2,...,N,\)做以下几步

• \(t_j=t_{j-1}+\Delta t\)
• 产生\(Z_j\sim N(0,1)\)
• \(\Delta B_j=Z_j\sqrt{\Delta t}\)
• \(X_j=X_{j-1}\exp\Big((\mu-\frac{1}{2}\sigma^2)\Delta t+\sigma \Delta B_j\Big)\)

set.seed(111)
mu <- 1
sigma <- 0.5
x <- 10
N <- 100
Delta <- T/N #时间间隔
B <- numeric(N+1) #初始化向量  这里的t_0=0,B_{t_0}=0
t <- seq(0,T,length = N+1)
for(i in 2:(N+1))
  B[i] <- B[i-1]+rnorm(1)*sqrt(Delta)

S <- x*exp((mu-sigma^2/2)*t+sigma*B)

plot(t,S,type="l",main="几何布朗运动")

并不是所有的SDE都能解出显式解，更多的SDE只能通过迭代求出数值解。求SDE数值解实际上也就是模拟出解的路径。

Euler格式

\[ Y_{i+1}=Y_i+\mu(t_i,Y_i)(t_{i+1}-t_i)+\sigma(t_i,Y_i)(B_{i+1}-B_i) \] ### Milstein格式

\[ \begin{align} Y_{i+1}&=Y_i+\mu(t_i,Y_i)(t_{i+1}-t_i)+\sigma(t_i,Y_i)(B_{i+1}-B_i)\\ &+\frac{1}{2}\sigma(t_i,Y_i)\sigma_x(t_i,Y_i)\{(B_{i+1}-B_i)^2-(t_{i+1}-t_i)\} \end{align} \] 其中\(B_{i+1}-B_i=\sqrt{t_{i+1}-t_i}Z_{i+1},\quad i=0,...,n-1.\)

SDE的参数估计

考虑下面的SDE: \[ dX_t=\mu(X_t;\theta)dt+\sigma(X_t,;\theta)dB_t \] 其中\(\theta\in \Theta\subset \mathbf{R}^p\)是\(p\)维参数，\(\theta_0\)是参数的真实值。函数\(\mu:\mathbf{R}\times\Theta\rightarrow \mathbf{R}\)和\(\sigma:\mathbf{R}\times\Theta\rightarrow(0,+\infty)\)是已知的，并且使得SDE的解存在。如果我们观测到该SDE解的一条轨迹的离散采样\(\{X_n,n\in \mathbb{N}\},\)用这些数据来估计参数\(\theta,\)得到的参数的估计为\(\hat{\theta}\).

极大似然估计

假设\(x_0,...,x_N\)是\(X_t\)在均匀离散时刻\(t_i=\Delta t\)的观测，其中\(i=0,1,...,N,\Delta t=T/N.\)

令\(p(t_k,x_k|t_{k-1},x_{k-1};\theta)\)是从\((t_{k-1},x_{k-1})\)到\((t_k,x_k)\)的传递概率密度（条件密度？），假设初始状态的概率密度为\(p_0(x_0|\theta),\)则似然函数为：

\[ f(\theta)=p_0(x_0|\theta)\prod_{k=1}^Np(t_k,x_k|t_{k-1},x_{k-1};\theta). \] 对数似然函数为：

\[ L(\theta)=\sum_{k=1}^N\ln p(t_k,x_k|t_{k-1},x_{k-1};\theta)+\ln p_0(x_0|\theta). \] 对SDE的Euler近似，有

\[ X_{k}=X_{k-1}+\mu(t_{k-1},x_{k-1};\theta)\Delta t+\sigma(t_{k-1},x_{k-1};\theta)\sqrt{\Delta t}\eta_k. \]

其中\(\eta_k\sim N(0,1).\)由上式可得，

\[ p(t_k,x_k|t_{k-1},x_{k-1};\theta)=\frac{1}{\sqrt{2\pi\sigma_k}}\exp\Big(-\frac{(x_k-\mu_k)^2}{2\sigma_k^2}\Big). \] 其中\(\mu_k=x_{k-1}+\mu(t_{k-1},x_{k-1};\theta),\sigma_k=\sigma(t_{k-1},x_{k-1};\theta)\sqrt{\Delta t}\).

考虑美洲鹤1939–1985年的种群数据，假设\(t\)时刻种群大小\(X_t\)满足SDE:

\[ dX_t=\theta_1X_tdt+\theta_2\sqrt{X_t}dB_t,\quad X_0=18. \]