%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))