%MILSTRONG Test strong convergence of Milstein: vectorized % % Solves dX = r*X*(K-X) dt + beta*X dW, X(0) = Xzero, % where r = 2, K= 1, beta = 1 and Xzero = 0.5. % % Discretized Brownian path over [0,1] has dt = 2^(-11). % Milstein uses timesteps 128*dt, 64*dt, 32*dt, 16*dt (also dt for reference). % % Examines strong convergence at T=1: E | X_L - X(T) |. % Code is vectorized: all paths computed simultaneously. randn('state',100) r = 2; K = 1; beta = 0.25; Xzero = 0.5; % problem parameters T = 1; N = 2^(9); % % Simulation time T and number of time steps N at finest % resolution M = 500; % number of paths sampled Nres=5; % number of resolution levels dt = T/N; % time step at finest resolution level dW = sqrt(dt)*randn(M,N); % Brownian increments Xmil = zeros(M,Nres); % preallocate array for p = 1:Nres R=2^(p-1); Dt = R*dt; L = N/R; % L timesteps of size Dt = R dt Xtemp = Xzero*ones(M,1); % start all trajectories at Xzero for j = 1:L Winc = sum(dW(:,R*(j-1)+1:R*j),2); % Brownian increments over time intervals of size Dt Xtemp = Xtemp + Dt*r*Xtemp.*(K-Xtemp) + beta*Xtemp.*Winc ... + 0.5*beta^2*Xtemp.*(Winc.^2 - Dt); % Milstein scheme end Xmil(:,p) = Xtemp; % store Milstein solution at t =1 end Xref = Xmil(:,1); % Reference solution is taken to be the one at % finest time step Xerr = abs(Xmil(:,2:Nres) - repmat(Xref,1,Nres-1)); % Error in each path mean(Xerr); % Calculate Mean pathwise (strong) errors at each % resolution level Dtvals = dt*2.^(1:Nres-1); % Milstein timesteps used subplot(224) % lower RH picture loglog(Dtvals,mean(Xerr),'b*-'), hold on loglog(Dtvals,Dtvals,'r--'), hold off % reference slope of 1 axis([1e-3 1e-1 1e-4 1]) xlabel('\Delta t') ylabel('Sample average of | X(T) - X_L |') title('milstrong.m','FontSize',10) %%%% Least squares fit of error = C * Dt^q %%%% A = [ones(Nres-1,1), log(Dtvals)']; rhs = log(mean(Xerr)'); sol = A\rhs; q = sol(2); resid = norm(A*sol - rhs); figure(1) disp(sprintf('The least squares fit for the slope of the strong Milstein error in the log-log plot is %f, with residual %f.',q,resid))