LESSON 3.6 Numerical Methods for SDEs – Euler-Maruyama and Milstein Schemes

1. LEARNING OBJECTIVES

By the end of this lesson, you will be able to:

Understand the need for numerical methods in solving SDEs.
Implement the Euler-Maruyama (EM) scheme for SDEs.
Analyze the strong and weak convergence properties of the EM scheme.
Implement the Milstein scheme for improved convergence.
Understand the difference between strong and weak convergence.
Apply the Euler-Maruyama scheme to Geometric Brownian Motion.
Implement variance reduction techniques (antithetic variates, control variates).
Apply numerical methods to price options using Monte Carlo simulation.

2. WHY NUMERICAL METHODS ARE NEEDED

Many SDEs in finance do not have closed-form solutions:

The CIR process has no closed-form solution.
The Heston model has no closed-form solution.
Jump-diffusion models have no closed-form solution.
Path-dependent options require Monte Carlo simulation.

Numerical methods allow us to approximate the solution to SDEs.

3. THE EULER-MARUYAMA (EM) SCHEME

3.1 The Scheme

Given the SDE:

dX_t = μ(t, X_t) dt + σ(t, X_t) dW_t

The Euler-Maruyama scheme discretizes time into N steps of size Δt = T / N:

X_{i+1} = X_i + μ(t_i, X_i) Δt + σ(t_i, X_i) ΔW_i

Where:

ΔW_i = W_{i+1} – W_i ~ N(0, Δt)
X_i is the approximation of X at time t_i = i Δt

3.2 Implementation Algorithm

Choose N (number of steps) and Δt = T / N.
Set X_0 = x_0.
For i = 0, 1, …, N-1:
a. Generate Z_i ~ N(0, 1).
b. Set ΔW_i = Z_i * √Δt.
c. Update X_{i+1} = X_i + μ(t_i, X_i) Δt + σ(t_i, X_i) ΔW_i.
Return X_N as the approximation of X_T.

3.3 Convergence Properties

Strong Convergence (Pathwise): EM has strong order 0.5:

E[|X_T – X_N|] ≤ C * Δt^(0.5)
Weak Convergence (Distributional): EM has weak order 1.0:

|E[g(X_T)] – E[g(X_N)]| ≤ C * Δt

3.4 Financial Application – Monte Carlo Option Pricing

To price a European option:

Generate M independent paths using the EM scheme.
Compute the payoff for each path: H_i = g(X_T^{(i)}).
Average the payoffs: H_avg = (1/M) Σ H_i.
Discount: V_0 = e^{-rT} * H_avg.

4. THE MILSTEIN SCHEME

4.1 The Scheme

The Milstein scheme improves the strong convergence order to 1.0:

X_{i+1} = X_i + μ(t_i, X_i) Δt + σ(t_i, X_i) ΔW_i + (1/2) σ(t_i, X_i) σ'(t_i, X_i) ((ΔW_i)² – Δt)

Where σ'(t, x) = ∂σ/∂x.

4.2 Derivation of the Milstein Scheme

The Milstein scheme is derived from the Ito-Taylor expansion of the SDE. It includes the second-order term in the diffusion coefficient.

4.3 Convergence Properties

Strong Convergence: Milstein has strong order 1.0:

E[|X_T – X_N|] ≤ C * Δt
Weak Convergence: Milstein has weak order 1.0 (same as EM for weak convergence).

4.4 When to Use Milstein

Use the Milstein scheme when:

Path-dependent options require strong convergence (e.g., Asian options, barrier options).
The volatility function depends on the state variable (σ'(x) ≠ 0).
You need higher accuracy for the same number of time steps.

5. STRONG VS. WEAK CONVERGENCE

5.1 Strong Convergence (Pathwise)

Measures the expected error between the true solution and the numerical approximation at each point along the path:

E[|X_T – X_N|] ≤ C * Δt^γ

Where γ is the strong order.

5.2 Weak Convergence (Distributional)

Measures the error in expectations:

|E[g(X_T)] – E[g(X_N)]| ≤ C * Δt^β

Where β is the weak order.

5.3 Comparison

Scheme	Strong Order	Weak Order
Euler-Maruyama	0.5	1.0
Milstein	1.0	1.0

5.4 Choosing a Scheme

For pricing European options (weak convergence is sufficient): EM is adequate.
For path-dependent options or hedging (strong convergence is required): Use Milstein.

6. VARIANCE REDUCTION TECHNIQUES

6.1 Antithetic Variates

For each path generated using Z_i, generate a second path using -Z_i.

Advantages:

Reduces variance by exploiting symmetry.
Simple to implement.
No extra random number generation.

Variance Reduction:

Var(V_antithetic) ≤ (1/2) Var(V_standard)

6.2 Control Variates

Use a known analytical solution (or a highly accurate approximation) as a control variate.

Let V be the option price we want to estimate. Let C be the price of a similar option with a known closed-form solution (e.g., a European option with the same parameters).

The control variate estimator is:

V_CV = V + β * (C_known – C_estimated)

Where:

V is the Monte Carlo estimate.
C_estimated is the Monte Carlo estimate of the control option.
C_known is the closed-form price of the control option.
β is the optimal coefficient:

β = – Cov(V, C) / Var(C)

6.3 Importance Sampling

Change the probability measure to sample more efficiently in regions that contribute more to the expectation.

This is used for pricing out-of-the-money options and rare event simulation.

6.4 Stratified Sampling

Partition the sample space into strata. Sample each stratum proportionally to its probability. This reduces variance by ensuring that each region is sampled adequately.

7. APPLICATION TO GEOMETRIC BROWNIAN MOTION

7.1 EM Scheme for GBM

For GBM:

dS_t = μ S_t dt + σ S_t dW_t

The EM scheme is:

S_{i+1} = S_i + μ S_i Δt + σ S_i ΔW_i

7.2 Milstein Scheme for GBM

For GBM, σ(t, S) = σ S, so σ'(t, S) = σ.

The Milstein scheme is:

S_{i+1} = S_i + μ S_i Δt + σ S_i ΔW_i + (1/2) σ² S_i ((ΔW_i)² – Δt)

7.3 Comparison of EM and Milstein for GBM

For GBM, we have a closed-form solution:

S_T = S_0 * exp( (μ – (1/2) σ²) T + σ W_T )

We can compare the EM and Milstein approximations to the true solution.

8. MONTE CARLO SIMULATION FOR OPTION PRICING

8.1 Algorithm for European Call Option

Set parameters: S_0, K, T, r, σ, M, N.
Compute Δt = T / N.
For m = 1, …, M:
a. Generate N random variables Z_i ~ N(0, 1).
b. Initialize S_0^{(m)} = S_0.
c. For i = 0, …, N-1:
- ΔW_i = Z_i * √Δt
- Update S_{i+1}^{(m)} = S_i^{(m)} + r S_i^{(m)} Δt + σ S_i^{(m)} ΔW_i
  d. Compute payoff: H_m = max(S_N^{(m)} – K, 0).
Compute average payoff: H_avg = (1/M) Σ H_m.
Discount: V_0 = e^{-rT} * H_avg.

8.2 Algorithm for Asian Option

For an Asian option, the payoff depends on the average price:

H = max( (1/N) Σ_{i=1}^N S_i – K, 0 )

The algorithm is the same, but the payoff is calculated using the average of all simulated prices.

8.3 Algorithm for Barrier Option

For a barrier option, the payoff is triggered only if the barrier is not hit:

H = max(S_T – K, 0) * 1_{max_{t≤T} S_t < B}

We must track the maximum price along the path to determine if the barrier was hit.

9. PRACTICAL IMPLEMENTATION

A. Euler-Maruyama for GBM:

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

def euler_maruyama_gbm(S0, mu, sigma, T, n_steps, n_paths):
    dt = T / n_steps
    S = np.zeros((n_paths, n_steps + 1))
    S[:, 0] = S0
    for i in range(n_steps):
        dW = np.random.normal(0, np.sqrt(dt), n_paths)
        S[:, i+1] = S[:, i] + mu * S[:, i] * dt + sigma * S[:, i] * dW
    return S

def milstein_gbm(S0, mu, sigma, T, n_steps, n_paths):
    dt = T / n_steps
    S = np.zeros((n_paths, n_steps + 1))
    S[:, 0] = S0
    for i in range(n_steps):
        dW = np.random.normal(0, np.sqrt(dt), n_paths)
        S[:, i+1] = S[:, i] + mu * S[:, i] * dt + sigma * S[:, i] * dW + 0.5 * sigma**2 * S[:, i] * (dW**2 - dt)
    return S

def exact_gbm(S0, mu, sigma, T, n_steps, n_paths):
    dt = T / n_steps
    t = np.linspace(0, T, n_steps + 1)
    W = np.zeros((n_paths, n_steps + 1))
    W[:, 1:] = np.cumsum(np.random.normal(0, np.sqrt(dt), (n_paths, n_steps)), axis=1)
    S = S0 * np.exp((mu - 0.5 * sigma**2) * t + sigma * W)
    return S

# Parameters
S0 = 100
mu = 0.08
sigma = 0.20
T = 1.0
n_steps = 100
n_paths = 1000

# Generate paths
S_EM = euler_maruyama_gbm(S0, mu, sigma, T, n_steps, n_paths)
S_Mil = milstein_gbm(S0, mu, sigma, T, n_steps, n_paths)
S_Exact = exact_gbm(S0, mu, sigma, T, n_steps, n_paths)

# Compare terminal distributions
print(f"Exact: Mean={np.mean(S_Exact[:, -1]):.2f}, Var={np.var(S_Exact[:, -1]):.2f}")
print(f"EM:     Mean={np.mean(S_EM[:, -1]):.2f}, Var={np.var(S_EM[:, -1]):.2f}")
print(f"Mil:    Mean={np.mean(S_Mil[:, -1]):.2f}, Var={np.var(S_Mil[:, -1]):.2f}")

# Error analysis
error_EM = np.mean((S_Exact[:, -1] - S_EM[:, -1])**2)
error_Mil = np.mean((S_Exact[:, -1] - S_Mil[:, -1])**2)
print(f"EM MSE: {error_EM:.4f}")
print(f"Mil MSE: {error_Mil:.4f}")

B. Monte Carlo Option Pricing with EM:

def euro_call_price_em(S0, K, T, r, sigma, n_steps, n_paths):
    dt = T / n_steps
    S = np.zeros((n_paths, n_steps + 1))
    S[:, 0] = S0
    for i in range(n_steps):
        dW = np.random.normal(0, np.sqrt(dt), n_paths)
        S[:, i+1] = S[:, i] + r * S[:, i] * dt + sigma * S[:, i] * dW
    payoffs = np.maximum(S[:, -1] - K, 0)
    price = np.exp(-r * T) * np.mean(payoffs)
    std_error = np.exp(-r * T) * np.std(payoffs) / np.sqrt(n_paths)
    return price, std_error

# Parameters
S0 = 100
K = 100
T = 1.0
r = 0.03
sigma = 0.20
n_steps = 50
n_paths = 10000

price, std_error = euro_call_price_em(S0, K, T, r, sigma, n_steps, n_paths)
print(f"Option Price: {price:.4f}")
print(f"Standard Error: {std_error:.4f}")

# Compare with Black-Scholes
from scipy.stats import norm
d1 = (np.log(S0 / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
bs_price = S0 * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
print(f"Black-Scholes Price: {bs_price:.4f}")
print(f"Difference: {price - bs_price:.4f}")

C. Variance Reduction with Antithetic Variates:

def euro_call_price_antithetic(S0, K, T, r, sigma, n_steps, n_paths):
    dt = T / n_steps
    # Generate half the paths, then use antithetic
    half_paths = n_paths // 2
    S1 = np.zeros((half_paths, n_steps + 1))
    S2 = np.zeros((half_paths, n_steps + 1))
    S1[:, 0] = S0
    S2[:, 0] = S0

    for i in range(n_steps):
        Z = np.random.normal(0, 1, half_paths)
        dW1 = Z * np.sqrt(dt)
        dW2 = -Z * np.sqrt(dt)
        S1[:, i+1] = S1[:, i] + r * S1[:, i] * dt + sigma * S1[:, i] * dW1
        S2[:, i+1] = S2[:, i] + r * S2[:, i] * dt + sigma * S2[:, i] * dW2

    payoffs1 = np.maximum(S1[:, -1] - K, 0)
    payoffs2 = np.maximum(S2[:, -1] - K, 0)
    all_payoffs = np.concatenate([payoffs1, payoffs2])
    price = np.exp(-r * T) * np.mean(all_payoffs)
    std_error = np.exp(-r * T) * np.std(all_payoffs) / np.sqrt(n_paths)
    return price, std_error

# Compare standard vs antithetic
price_std, se_std = euro_call_price_em(S0, K, T, r, sigma, n_steps, n_paths)
price_anti, se_anti = euro_call_price_antithetic(S0, K, T, r, sigma, n_steps, n_paths)

print(f"Standard: Price={price_std:.4f}, SE={se_std:.4f}")
print(f"Antithetic: Price={price_anti:.4f}, SE={se_anti:.4f}")
print(f"Variance Reduction: {1 - se_anti**2 / se_std**2:.2%}")

D. Convergence Analysis:

def convergence_analysis(S0, mu, sigma, T, n_paths, n_steps_list):
    errors_EM = []
    errors_Mil = []
    dt_list = []

    for n_steps in n_steps_list:
        dt = T / n_steps
        dt_list.append(dt)
        S_Exact = exact_gbm(S0, mu, sigma, T, n_steps, n_paths)
        S_EM = euler_maruyama_gbm(S0, mu, sigma, T, n_steps, n_paths)
        S_Mil = milstein_gbm(S0, mu, sigma, T, n_steps, n_paths)

        error_EM = np.mean((S_Exact[:, -1] - S_EM[:, -1])**2)
        error_Mil = np.mean((S_Exact[:, -1] - S_Mil[:, -1])**2)
        errors_EM.append(error_EM)
        errors_Mil.append(error_Mil)

    return dt_list, errors_EM, errors_Mil

n_steps_list = [10, 20, 50, 100, 200, 500]
dt_list, errors_EM, errors_Mil = convergence_analysis(S0, mu, sigma, T, 1000, n_steps_list)

# Plot convergence
plt.figure(figsize=(10, 6))
plt.loglog(dt_list, errors_EM, 'o-', label='Euler-Maruyama')
plt.loglog(dt_list, errors_Mil, 's-', label='Milstein')
plt.loglog(dt_list, [dt**0.5 for dt in dt_list], '--', label='O(Δt^0.5)')
plt.loglog(dt_list, [dt for dt in dt_list], '--', label='O(Δt)')
plt.xlabel('Δt')
plt.ylabel('MSE')
plt.title('Convergence Analysis')
plt.legend()
plt.show()

10. SUMMARY FOR THE FINANCE PRACTITIONER

Euler-Maruyama is the simplest numerical scheme for SDEs. It has strong order 0.5 and weak order 1.0.
Milstein improves strong convergence to order 1.0 by adding the (ΔW)² – Δt term. Use Milstein for path-dependent options.
Strong convergence is pathwise; weak convergence is distributional. For European option pricing, weak convergence is sufficient.
Variance reduction techniques (antithetic variates, control variates) reduce Monte Carlo variance without increasing the number of paths.
Convergence analysis confirms that EM converges at rate O(Δt^0.5) and Milstein at O(Δt).

Diploma in Financial Technology LESSON 3.6 Numerical Methods for SDEs – Euler-Maruyama and Milstein Schemes

Diploma in Financial Technology