Euler-Maruyama 方法数值算例

pdf原文件

将随机微分方程写成积分形式有
$X(t)=X_0+\int_0 ^t f(X(s))ds+\int_0 ^t g(X(s))dW(s), 0\le t\le T. (1)$

其中 $f, g$ 都是标量函数并且初始条件 $X_0$ 是随机变量。这里的积分采用的是 Ito 积分,方程的解 $X(t)$ 是随 $t$ 变化的随机变量。现在定义一种求解上述方程的数值方法,并且将数值方法在时间步长趋于 0 时所得的随机变量 $X(t)$ 作为方程的解。将方程写为微分形式有

$dX(t) = f (X(t))dt + g(X(t))dW (t), X(0) = X_0 , 0\le t \le T. (2)$

在区间 $[0, T ]$ 采用一个算法来求解上面微分方程,首先进行离散处理。使 $\Delta t = T /L$ ,其中 $L$ 为某一正整数,并且有 $\tau _j= j \Delta t$ 。数值近似 $X(\tau _j)$ 表示成 $X_j$ 。Euler-Maruyama 具有如下形式

$X_j = X_{j-1} + f (X_{j-1})\Delta t + g(X_{j-1})(W(\tau _j)- W (\tau _{j-1} )), j = 1, 2, ..., L. (3)$

方程 (3) 可以由积分形式来理解,积分形式为

$X(\tau _j ) = X(\tau _{j-1} ) +\int_{\tau _{j-1}} ^{\tau _j} f (X(s))ds +\int _{\tau _{j-1}} ^{\tau _j} g(X(s))dW (s). (4)$

方程 (3) 中每一项都和方程 (4) 相对应,当 $g\equiv 0$ 时,方程 (3) 退化成 Euler 方法。

首先我们将离散 Brownian 路径来生成方程 (3) 所需要的增量 $W (\tau _j)- W (\tau_{j-1})$ ,为了方便起见,我们将时间步长选为 Brownian 路径增量 $\delta t$ 的整数倍 $R\delta t$ ,其中 $R\ge 1$ 为整数。这样可以使得采用 Euler-Maruyama 计算的时间点序列 $\{ \tau _j \}$ 是包含在 Brownian 路径的时间序列 $\{ \tau _j \}$ 中的。列举一个将 Euler-Maruyama 方法应用到线性 SDE 的例子

$dX(t) = \lambda X(t)dt + \mu X(t)dW (t), X(0) = X_0. (5)$

其中 $\lambda$ 和 $\mu$ 都是实常数;因此在方程 (2) 中 $f (X) = \lambda X$ 和 $g(X) = \mu X$ 。这个是金融数学中的Black-Scholes 偏微分方程可以由方程 (5) 得到 [2]。方程 (5) 的解析解为 [3]

$X(t) = X(0) \exp((\lambda -\mu ^2)t + \mu W (t)). (6)$

取 $\lambda = 2, \mu = 1, X_0 = 1$ 。令 Brownian 路径在区间 $[0, 1]$ 之间和 $\delta t = 2^{-8}$ ,并且按照方程 (6) 得到方程的解析解 $Xture$ 。然后采用 Euler-Maruyama 计算时采用时间步 $\Delta t = R\delta t$ ,令 $R = 4$ 。采用Euler-Maruyama 的一般步长所对应的增量 $W (\tau _j ) -W (\tau _{j-1} )$ 的数值为

$W (\tau _j ) - W (\tau_{j-1}) = W (jR\delta t)-W ((j-1)R\delta t) =\sum _{k=jR-R+1} ^{jR} dW_k.$

由此得到 Euler-Maruyama 方法的数值解。得到结果如图 1

此时在 $t=T$ 的点数值解和解析解的误差为 $emerr = 0.6907$ ,当取 $R = 1$ 时,如图 1,在 $t = T$ 的点数值解和解析解的误差为 $emerr = 0.0821$ 。

Matlab程序

randn('state',100)
lambda=2;mu=1;Xzero=1;
T=1;N=5000;dt=1/N;
dW=sqrt(dt)*randn(1,N);
W=cumsum(dW);Xtrue=Xzero*exp((lambda-0.5*mu^2)*([dt:dt:T])+mu*W);
plot([0:dt:T],[Xzero,Xtrue],'k-'),hold on;R=1;Dt=R*dt;L=N/R;
Xem=zeros(1,L);
Xtemp_em=Xzero;
for j=1:Lt=j*Dt;Winc=sum(dW(R*(j-1)+1:R*j));Xtemp_em=Xtemp_em+Dt*lambda*Xtemp_em+mu*Xtemp_em*Winc; Xem(j)=Xtemp_em;
endplot([0:Dt:T],[Xzero,Xem],'k--*'),hold off;
xlabel('{$t$}','Interpreter','latex','FontSize',12)
ylabel('{$X_t$}','Interpreter','latex','FontSize',16,'Rotation',0,'HorizontalAlignment','right')emerr=abs(Xem(end)-Xtrue(end));

参考文献

[1] Desmond J. Higham, An Algorithmic Introduction to Numerical Simulation of Stochatic Differential Equations, SIAM REVIEW, 43(3), 2001, 525-546.

[2] J. C. Hull, Options, Futures, and Other Derivatives, 6th ed., Pearson PH, 2009, 307-309.

[3] X. Mao, Stochatic Differential Equations and Applications, Horwood, Chichester, 1997, 105.

Euler-Maruyama 方法数值算例相关推荐

matlab设计降维状态观测器,降维状态观测器的两种设计方法及算例
第 17卷第 3期 2005年 6月军械工程学院学报 Journal of Ordnance Engineering College Vo1．17 No．3 Jun．,2005 文章编 ...
综合能源系统通用建模及规划方法研究-算例-笔记
综合能源系统通用建模及规划方法研究--算例分析区域综合能源系统规划对象包括能量枢纽和区域能源配网,前者的规划需综合考虑系统结构和设备型号,为区域综合能源系统规划的核心内容,为了更直接说明本文提出 ...
【微分方程数值解】有限差分法（二）两点边值问题数值算例（附python代码）
1.两点边值问题形式一般的两点边值问题形式为: −uxx=f(x)-u_{xx}=f(x)−uxx=f(x) u(a)=α,u(b)=βu(a)=\alpha,u(b)=\betau(a)=α,u ...
【微分方程数值解】常微分方程（一）欧拉方法和改进欧拉方法（附python算例，封装类）
欧拉方法与改进欧拉方法一.算法原理对给定微分方程 {y′=f(x,y)y(x0)=y0(1)\begin{cases} y' = f(x,y)\\ y(x_0) = y_0 \end{cases} ...
有限元方法入门：有限元方法简单的二维算例（三角形剖分）
有限元方法简单的二维算例(三角形剖分) 文章目录有限元方法简单的二维算例(三角形剖分) 算例描述变分问题有限元离散问题转化有限元三要素参考单元与一般单元一般单元上的形函数一般单元上的积 ...
有限元方法入门：有限元方法简单的二维算例（矩形剖分）
#有限元方法简单的二维算例(矩形剖分) 算例描述我们对下述椭圆边值问题 \label{eq1}{−Δu=fu∣∂Ω=0\left \{ \begin{aligned} & -\Delta u ...
万字字符长文带你了解遗传算法（有几个算例源码）
一.遗传算法的基本概念简单而言,遗传算法使用群体搜索技术,将种群代表一组问题解, 通过对当前种群施加选择.交叉和变异等一系列遗传操作来产生新-一代的种群,并逐步使种群进化到包含近似最优解的状态.由于 ...
动态规划原理介绍（附7个算例，有代码讲解）
动态规划思想动态规划(Dynamicprogramming)是一种通过把原问题分解为相对简单的子问题的方式求解复杂问题的方法. 动态规划常常适用于有重叠子问题和最优子结构性质的问题,动态规划方法所耗 ...
qr分解求线性方程组_梯度下降求解线性方程组算例设计
凸二次优化问题 Theory. 设是实对称正定矩阵, ,则求解凸二次优化问题等价于求解线性方程组 . Proof. 二次型二阶可导,极小值点处梯度为零,现对优化的目标函数求梯度.二次型本质上具有: ...

Euler-Maruyama 方法数值算例

Euler-Maruyama 方法数值算例相关推荐

最新文章

热门文章