UA MATH575B 数值分析下VI 统计物理的随机模拟方法1

模拟SDE的解

以最简单的随机微分方程为例，考虑随机游走
$$d{x}_t={\xi }_{t}dt$$
其中${\xi }_t$是白噪声，假设$⟨{\xi }_{t_1},{\xi }_{t_2}⟩=2D\delta (t_1-t_2)$，其中$D>0$，$\delta$是Dirac函数。假设随机过程$x_t$的概率密度函数为$\rho (t,x)$，它满足Fokker-Planck方程：
$$\frac{\partial \rho }{\partial t}=D\frac{{\partial }^2\rho }{\partial x^2}$$
这个是一维的扩散方程，$D$叫做扩散系数。这个方程的一个解是
$$\rho (t,x)=\frac{1}{\sqrt{4\pi Dt}}\exp \left(-\frac{x^2}{4Dt}\right)$$
即$x_t$是高斯过程，分布可以记为$N(0,2Dt)$。

现在对时间做离散化，记time step为$\tau$，简单回顾一下解数值ODE的Runge-Kutta方法 (RK)：对于动力系统$\dot{x}=f(t,x)\in {\mathbb{R}}^m$，定义RK的系数为$A\in {\mathbb{R}}^{s\times s},b\in {\mathbb{R}}^s,c\in {\mathbb{R}}^s$
$$\begin{gathered} k_i=hf(t+c_{i}h,x(t)+\sum _{j=1}^s{A}_{ij}{k}_j) \\ x(t+h)=x(t)+\sum _{j=1}^s{b}_j{k}_j \end{gathered}$$
其中$s$是RK方法的阶数，$h$是步长，根据上述递推关系，可以计算出$x(t)$的序列。

对于随机微分方程$d{x}_t={\xi }_{t}dt$，考虑RK2：
$$A=\left[\begin{array}{cc}0& 0\\ \frac{1}{2}& 0\end{array}\right],b=[0,1{]}^T,c=[0,\frac{1}{2}{]}^T,h=\tau$$
则
$$\begin{gathered} k_1=\tau {\xi }_i,k_2=\tau {\xi }_{i+\frac{1}{2}} \\ {x}_{i+1}=x_i+\tau {\xi }_{i+\frac{1}{2}} \end{gathered}$$
$$x_{i+1}=x_i+\tau {\xi }_{i+\frac{1}{2}}$$
其中$⟨{\xi }_i,{\xi }_j⟩=2D{\delta }_{ij}/\tau$，其中${\delta }_{ij}$为Kronecker符号，所有${\xi }_i$都是独立同分布的。因此可以将上式写成：
$$\begin{gathered} x_{i+1}=x_i+\sqrt{2D\tau }\zeta \\ \zeta \sim N(0,1) \end{gathered}$$

要解这个随机过程可以用随机模拟的方法：

1. 第一步，从标准正态分布中抽一个样本$\zeta$；
2. 第二步，根据$x_{i+1}=x_i+\sqrt{2D\tau }\zeta$更新随机过程的值；

重复这两步，直到达到设置的时间上限。

由此我们可以总结出模拟一个随机微分方程的一般性做法：

1. 用数值ODE的方法导出SDE解的递推关系
2. 生成正态分布的随机数，根据递推关系模拟这个SDE

例1下面的SDE表示几何随机游走
$$d{x}_t=x_t{\xi }_{t}dt$$
其中${\xi }_t$是白噪声，假设$⟨{\xi }_{t_1},{\xi }_{t_2}⟩=2D\delta (t_1-t_2)$，$D>0$。用这个例子介绍两种离散化的规则：

1. Ito规则：$x_{i+1}=x_i+\sqrt{2D\tau }{x}_i\zeta$
2. Stratonovich规则：$x_{i+1}=x_i+\sqrt{2D\tau }\frac{x_i+x_{i+1}}{2}\zeta$

在这个例子中，Ito规则导出的结果是$E{x}_{i+1}=E{x}_i=x_0$，这个比较明显。但Stratonovich规则导出的是
$$\begin{gathered} x_{i+1}=x_i\frac{1+\sqrt{2D\tau }\frac{\zeta }{2}}{1-\sqrt{2D\tau }\frac{\zeta }{2}}=x_i(1+\sqrt{2D\tau }\frac{\zeta }{2})(1+\sqrt{2D\tau }\frac{\zeta }{2}+2D\tau \frac{{\zeta }^2}{4}+o({\zeta }^2)) \\ =x_i(1+\sqrt{2D\tau }\zeta +D\tau {\zeta }^2+o({\zeta }^2))\to x_i(1+D\tau +\sqrt{2D\tau }\zeta ) \end{gathered}$$
因此
$$E{x}_{i+1}=E{x}_i(1+D\tau )$$
因此这两个规则模拟出来的结果会有明显差别。

用蒙特卡罗方法估计SDE解的自相关函数

随着时间流逝，$x_t$信息减少的规律可以用自相关函数来表示，
$$K_{XX}(\tau )=⟨(x_{t+\tau }-x_t{)}^2⟩={\int }_t^{t+\tau }d{t}_1{\int }_t^{t+\tau }d{t}_2⟨{\xi }_{t_1}{\xi }_{t_2}⟩=2D\tau$$
随机游走的自相关函数显然只与两个时点的差有关，这种随机过程叫平稳随机过程。如果SDE的解比较复杂，甚至根本写不出解析解，很难直接计算自相关函数，我们可以用Monte Carlo方法估计自相关函数：
$$K_{XX}(\tau )=E{x}_t{x}_{t+\tau }\approx \frac{1}{T}\sum _{i=1}^T{x}_i{x}_{i+\tau }$$
其中$x_i,x_{i+\tau }$可以用随机模拟得到。

例2 一个简单的单位根过程，假设扩散系数$D=\frac{1}{2}$
$$d{x}_t=(-x_t+{\xi }_t)dt$$
参考SDE的第一讲，这个方程定义了一个Ornstein-Uhlenbeck过程，它的解是
$$x_t={\int }_{-\infty }^t{e}^{-(t-s)}{\xi }_{s}ds$$
根据这个解计算自相关函数，第一行到第二行用到的性质是Ito Isometry：
$$\begin{gathered} K_{XX}(\tau )=E{x}_t{x}_{t+\tau }=E{\int }_{-\infty }^t{e}^{-(t-s)}{\xi }_{s}ds{\int }_{-\infty }^{t+\tau }{e}^{-(t+\tau -r)}{\xi }_{r}dr \\ ={\int }_{-\infty }^t{\int }_{-\infty }^{t+\tau }{e}^{-(t-s)}{e}^{-(t+\tau -r)}E{\xi }_s{\xi }_{r}dsdr \\ ={\int }_{-\infty }^t{\int }_{-\infty }^{t+\tau }{e}^{-(t-s)}{e}^{-(t+\tau -r)}\delta (r-s)dsdr=\frac{1}{2}\exp (-\tau ),\tau \ge 0 \end{gathered}$$

如果我们不知道Ito Isometry，要计算自相关函数的解析解就会比较麻烦，这时就可以用Monte Carlo方法来估计自相关函数。先用RK2得到SDE的数值算法：
$$x_{i+1}=e^{-\tau }{x}_i+\sqrt{\tau }\zeta \approx (1-\tau )x_i+\sqrt{\tau }\zeta$$
根据这个递推关系模拟SDE的解，然后根据
$$K_{XX}(\tau )=E{x}_t{x}_{t+\tau }\approx \frac{1}{T}\sum _{i=1}^T{x}_i{x}_{i+\tau }$$
估计自相关函数。