function [tau, N, std_err] = mlmc(X0, a, b, B, T, M, epsilon, options)
if nargin < 8
options = '';
end
 if any(options == 'A')
f = @(dW, X, h) X + a * h + b * dW;
g = @(X) b;
else
g = @(X) b * X;
if any(options == 'E')
f = @(dW, X, h) X .* (1 + a * h + b * dW);
else
f = @(dW, X, h) X .* (1 + a * h + b * dW + 0.5 * b^2 * (dW.^2 ? h));
end
end
L = 0;
N_initial = 100;
converged = false;
A = zeros(3,1);...,
while ¬ converged
dif = mlmc_sim(X0, f, g, B, T, M, L, N_initial, options);
A(1, L+1) = N_initial;
A(2, L+1) = dif(1);
A(3, L+1) = dif(2);
vars = A(3, :) ./ A(1, :) ? (A(2, :) ./ A(1, :)).^2;
N = max(N_initial, ceil(2 * sqrt(vars ./ (M.^(0:L))) * sum(sqrt(vars .* ...
(M.^(0:L)))) / epsilon^2));
for l = 0:L
dN = N(l+1) ? A(1, l+1);
if dN > 0
dif = mlmc_sim(X0, f, g, B, T, M, l, dN, options);
A(1, l+1) = A(1, l+1) + dN;
A(2, l+1) = A(2, l+1) + dif(1);
A(3, l+1) = A(3, l+1) + dif(2);
end
end
if L > 1 && M^L ? 16
if M^L > 1024
converged = true;
else
converged = max( abs(A(2, L) / A(1, L)) / M, ...
                 abs(A(2, L+1) / A(1, L+1)) ...
                ) < (M?1)* epsilon / sqrt(2);
end
end
L = L + 1;
end
vars = (A(3, :) ? A(2, :).^2 ./ A(1, :)) ./ (A(1, :) ? 1);
86
std_err = sqrt(sum(vars ./ A(1, :)));
tau = sum(A(2, :) ./ A(1, :));
end
function [dif, fine] = mlmc_sim(X0, f, g, B, T, M, l, N, options)
K_f = M^l;
K_c = K_f / M;
h_f = T / K_f;
h_c = T / K_c;
sqrt_h_f = sqrt(h_f);
sum_f = 0;
sum_dif = 0;
square_sum_f = 0;
square_sum_dif = 0;
if any(options == 's')
simple = true;
minimum = false;
probability = false;
elseif any(options == 'm')
simple = false;
minimum = true;
probability = false;
else
simple = false;
minimum = false;
probability = true;
end
122
123
for n_ = 0:10000:N
n = min(10000, N ? n_);
if n == 0
break
end
X_f = X0 * ones(1, n);
X_c = X_f;
tau_f = T * ones(1, n);
tau_c = T * ones(1, n);
if probability
tau_f = zeros(1, n);
tau_c = zeros(1, n);
P_f = ones(1, n);
P_c = ones(1, n);
end
if l == 0
dW = sqrt_h_f * randn(1, n);
X_old = X_f;
X_f = f(dW, X_f, h_f);
if simple
tau_f(X_f < B) = T/2;
elseif minimum
U = rand(1, n);
mins = 0.5 * (X_f + X_old ? sqrt((X_f ? X_old).^2 ? 2 * h_f * g(X_old).^2 .* ...
log(U)));
tau_f(mins < B) = T/2;
else
p = exp(?2* (X_old ? B) .* (X_f ? B) ./ (g(X_old).^2 * h_f));
p(X_f < B) = 1;
tau_f = p * T/2 + (1 ? p) * T;
end
else
for i = 1:K_c
dW_f = sqrt_h_f * randn(M, n);
if minimum
U = rand(M, n);
end
for m = 1:M
X_old_f = X_f;
X_f = f(dW_f(m, :), X_f, h_f);
if simple
tau_f(X_f < B) = min(((i ? 1) * M + m ? 0.5) * h_f, tau_f(X_f < B));
elseif minimum
mins_f = 0.5 * (X_f + X_old_f ? sqrt((X_f ? X_old_f).^2 ? 2 * h_f * ...
g(X_old_f).^2 .* log(U(m, :))));
tau_f(mins_f < B) = min(((i ? 1) * M + m ? 0.5) * h_f, tau_f(mins_f < B));
else
p_f = exp(?2* (X_old_f ? B) .* (X_f ? B) ./ (g(X_old_f).^2 * h_f));
p_f(X_f < B) = 1;
tau_f = tau_f + P_f .* p_f * ((i ? 1) * M + m ? 0.5) * h_f;
P_f = P_f .* (1 ? p_f);
end
end
X_old_c = X_c;
dW_c = sum(dW_f);
X_c = f(dW_c, X_c, h_c);
if simple
tau_c(X_c < B) = min((i ? 0.5) * h_c, tau_c(X_c < B));
elseif minimum || probability
b = g(X_old_c);
if M == 2
X_half_c = 0.5 * (X_old_c + X_c ? b .* (dW_f(2, :) ? dW_f(1, :)));
elseif M == 4
X_half_c = 0.5 * (X_old_c + X_c ? b .* (dW_f(3, :) + dW_f(4, :) ? ...
dW_f(1, :) ? dW_f(2, :)));
X_quarter_1_c = 0.5 * (X_old_c + X_half_c ? b .* (dW_f(2, :) ? dW_f(1, :)));
X_quarter_3_c = 0.5 * (X_half_c + X_c ? b .* (dW_f(4, :) ? dW_f(3, :)));
end
if minimum
if M == 2
mins_1 = 0.5 * (X_half_c + X_old_c ? sqrt((X_half_c ? X_old_c).^2 ? ...
h_c * b.^2 .* log(U(1, :))));
mins_2 = 0.5 * (X_c + X_half_c ? sqrt((X_c ? X_half_c).^2 ? h_c * b.^2 ...
.* log(U(2, :))));
mins_c = min(mins_1, mins_2);
elseif M == 4
mins_1 = 0.5 * (X_quarter_1_c + X_old_c ? sqrt((X_quarter_1_c ? ...
X_old_c).^2 ? 0.5 * h_c * b.^2 .* log(U(1, :))));
mins_2 = 0.5 * (X_half_c + X_quarter_1_c ? sqrt((X_half_c ? ...
X_quarter_1_c).^2 ? 0.5 * h_c * b.^2 .* log(U(2, :))));
mins_3 = 0.5 * (X_quarter_3_c + X_half_c ? sqrt((X_quarter_3_c ? ...
X_half_c).^2 ? 0.5 * h_c * b.^2 .* log(U(3, :))));
mins_4 = 0.5 * (X_c + X_quarter_3_c ? sqrt((X_c ? X_quarter_3_c).^2 ? ...
0.5 * h_c * b.^2 .* log(U(4, :))));
mins_c = min(min(mins_1, mins_2), min(mins_3, mins_4));
else
U_c = rand(1, n);
mins_c = 0.5 * (X_c + X_old_c ? sqrt((X_c ? X_old_c).^2 ? 2 * h_c * ...
b.^2 .* log(U_c)));
end
tau_c(mins_c < B) = min((i ? 0.5) * h_c, tau_c(mins_c < B));
else
if M == 2
p1_c = exp(?4 * (X_old_c ? B) .* (X_half_c ? B) ./ (g(X_old_c).^2 * h_c));
p2_c = exp(?4 * (X_half_c ? B) .* (X_c ? B) ./ (g(X_old_c).^2 * h_c));
p_c = p1_c .* (1 ? p2_c) + p2_c;
p_c(X_half_c < B | X_c < B) = 1;
elseif M == 4
p1_c = exp(?8 * (X_old_c ? B) .* (X_quarter_1_c ? B) ./ (g(X_old_c).^2 ...
* h_c));
p2_c = exp(?8 * (X_quarter_1_c ? B) .* (X_half_c ? B) ./ ...
(g(X_old_c).^2 * h_c));
p3_c = exp(?8 * (X_half_c ? B) .* (X_quarter_3_c ? B) ./ ...
(g(X_old_c).^2 * h_c));
p4_c = exp(?8 * (X_quarter_3_c ? B) .* (X_c ? B) ./ (g(X_old_c).^2 * ...
h_c));
p_c = 1 ? (1 ? p1_c) .* (1 ? p2_c) .* (1 ? p3_c) .* (1 ? p4_c);
p_c(X_quarter_1_c < B | X_half_c < B | X_quarter_3_c < B | X_c < B) = 1;
else
p_c = exp(?2* (X_old_c ? B) .* (X_c ? B) ./ (g(X_old_c).^2 * h_c));
p_c(X_c < B) = 1;
end
tau_c = tau_c + P_c .* p_c * (i ? 0.5) * h_c;
P_c = P_c .* (1 ? p_c);
end
end
end
if probability
tau_f = tau_f + P_f * T;
tau_c = tau_c + P_c * T;
end
end
sum_dif = sum_dif + sum(tau_f ? tau_c);
sum_f = sum_f + sum(tau_f);
square_sum_dif = square_sum_dif + sum((tau_f ? tau_c).^2);
square_sum_f = square_sum_f + sum(tau_f.^2);
end
dif = [sum_dif, square_sum_dif];
fine = [sum_f, square_sum_f];
if l == 0
dif = fine;
end
end