View difference between Paste ID: Khp3T952 and EasayNG2
SHOW: | | - or go back to the newest paste.
1-
%EM Euler-Maruyama method on linear SDE
1+
%EM Euler-Maruyama method on linear SDE
2-
%
2+
%
3-
% SDE is dX = lambda*X dt + mu*X dW, X(0) = Xzero,
3+
% SDE is dX = lambda*X dt + mu*X dW, X(0) = Xzero,
4-
% where lambda = 2, mu = 1 and Xzero = 1.
4+
% where lambda = 2, mu = 1 and Xzero = 1.
5-
%
5+
%
6-
% Discretized Brownian path over [0,1] has dt = 2^(-8).
6+
% Discretized Brownian path over [0,1] has dt = 2^(-8).
7-
% Euler-Maruyama uses timestep R*dt.
7+
% Euler-Maruyama uses timestep R*dt.
8-
randn(’state’,100)
8+
randn(’state’,100)
9-
lambda = 2; mu = 1; Xzero = 1; % problem parameters
9+
lambda = 2; mu = 1; Xzero = 1; % problem parameters
10-
T = 1; N = 2^8; dt = 1/N;
10+
T = 1; N = 2^8; dt = 1/N;
11-
dW = sqrt(dt)*randn(1,N); % Brownian increments
11+
dW = sqrt(dt)*randn(1,N); % Brownian increments
12-
W = cumsum(dW); % discretized Brownian path
12+
W = cumsum(dW); % discretized Brownian path
13-
Xtrue = Xzero*exp((lambda-0.5*mu^2)*([dt:dt:T])+mu*W);
13+
Xtrue = Xzero*exp((lambda-0.5*mu^2)*([dt:dt:T])+mu*W);
14-
plot([0:dt:T],[Xzero,Xtrue],’m-’), hold on
14+
plot([0:dt:T],[Xzero,Xtrue],’m-’), hold on
15-
R = 4; Dt = R*dt; L = N/R; % L EM steps of size Dt = R*dt
15+
R = 4; Dt = R*dt; L = N/R; % L EM steps of size Dt = R*dt
16-
Xem = zeros(1,L); % preallocate for efficiency
16+
Xem = zeros(1,L); % preallocate for efficiency
17-
Xtemp = Xzero;
17+
Xtemp = Xzero;
18-
for j = 1:L
18+
for j = 1:L
19-
Winc = sum(dW(R*(j-1)+1:R*j));
19+
Winc = sum(dW(R*(j-1)+1:R*j));
20-
Xtemp = Xtemp + Dt*lambda*Xtemp + mu*Xtemp*Winc;
20+
Xtemp = Xtemp + Dt*lambda*Xtemp + mu*Xtemp*Winc;
21-
Xem(j) = Xtemp;
21+
Xem(j) = Xtemp;
22-
end
22+
end
23-
plot([0:Dt:T],[Xzero,Xem],’r--*’), hold off
23+
plot([0:Dt:T],[Xzero,Xem],’r--*’), hold off
24-
xlabel(’t’,’FontSize’,12)
24+
xlabel(’t’,’FontSize’,12)
25-
ylabel(’X’,’FontSize’,16,’Rotation’,0,’HorizontalAlignment’,’right’)
25+
ylabel(’X’,’FontSize’,16,’Rotation’,0,’HorizontalAlignment’,’right’)
26
emerr = abs(Xem(end)-Xtrue(end))