将随机微分方程写成积分形式有
其中都是标量函数并且初始条件是随机变量。这里的积分采用的是 Ito 积分,方程的解是随变化的随机变量。现在定义一种求解上述方程的数值方法,并且将数值方法在时间步长趋于 0 时所得的随机变量作为方程的解。将方程写为微分形式有
在区间采用一个算法来求解上面微分方程,首先进行离散处理。使 ,其中为某一正整数,并且有。数值近似 表示成。Euler-Maruyama 具有如下形式
方程 (3) 可以由积分形式来理解,积分形式为
方程 (3) 中每一项都和方程 (4) 相对应,当时,方程 (3) 退化成 Euler 方法。
首先我们将离散 Brownian 路径来生成方程 (3) 所需要的增量,为了方便起见,我们将时间步长选为 Brownian 路径增量的整数倍,其中为整数。这样可以使得采用 Euler-Maruyama 计算的时间点序列是包含在 Brownian 路径的时间序列中的。列举一个将 Euler-Maruyama 方法应用到线性 SDE 的例子
其中和都是实常数;因此在方程 (2) 中和 。这个是金融数学中的Black-Scholes 偏微分方程可以由方程 (5) 得到 [2]。方程 (5) 的解析解为 [3]
取 。令 Brownian 路径在区间之间和,并且按照方程 (6) 得到方程的解析解 。然后采用 Euler-Maruyama 计算时采用时间步,令。采用Euler-Maruyama 的一般步长所对应的增量的数值为
由此得到 Euler-Maruyama 方法的数值解。得到结果如图 1
此时在的点数值解和解析解的误差为,当取时,如图 1,在的点数值解和解析解的误差为 。
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:L
t=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;
end
plot([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.