Michael Brenndoerfer | Data & AI, Private Equity, Technology

Michael Brenndoerfer

Itô's Lemma and Stochastic Calculus › Code Implementation

Michael BrenndoerferNovember 23, 202511 min read

Quantitative FinanceAll chapters

Part 5 of 15Stochastic Calculus and Pricing Theory

Chapter 27Itô's Lemma and Stochastic Calculus › Code Implementation

Reading section 6 of 11 · Code Implementation · 11 minSection 6 of 11 · 11 min

Section 6 of 11 · 11 min read

Itô's Lemma and Stochastic Calculus › Code ImplementationOne section · Sections in this chapter6 / 11 ↓

Sections in this chapter11

Reading tools Definitions & inline links

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

Let's implement simulations to verify our analytical results from Itô's Lemma.

In[6]:

Code

import numpy as np

np.random.seed(42)


def simulate_gbm(S0, mu, sigma, T, n_steps, n_paths):
    """
    Simulate geometric Brownian motion paths using Euler-Maruyama discretization.

    Parameters:
    -----------
    S0 : float - Initial price
    mu : float - Drift rate
    sigma : float - Volatility
    T : float - Time horizon
    n_steps : int - Number of time steps
    n_paths : int - Number of simulation paths

    Returns:
    --------
    t : array - Time points
    S : array - Simulated price paths (n_steps+1 x n_paths)
    W : array - Brownian motion paths (n_steps+1 x n_paths)
    """
    dt = T / n_steps
    t = np.linspace(0, T, n_steps + 1)

    # Generate random increments
    dW = np.random.normal(0, np.sqrt(dt), (n_steps, n_paths))

    # Build Brownian motion paths
    W = np.zeros((n_steps + 1, n_paths))
    W[1:, :] = np.cumsum(dW, axis=0)

    # Simulate price paths using Euler-Maruyama
    S = np.zeros((n_steps + 1, n_paths))
    S[0, :] = S0

    for i in range(n_steps):
        S[i + 1, :] = S[i, :] * (1 + mu * dt + sigma * dW[i, :])

    return t, S, W

Verifying the GBM SolutionLink Copied

We derived that $S_T = S_0 \exp((\mu - \frac{1}{2}\sigma^2)T + \sigma W_T)$ . Let's compare the Euler-Maruyama simulation with the exact solution.

In[7]:

Code

import numpy as np

# Parameters
S0 = 100
mu = 0.10  # 10% annual drift
sigma = 0.20  # 20% annual volatility
T = 1.0  # 1 year
n_steps = 252  # Daily steps
n_paths = 10000

# Simulate
t, S_euler, W = simulate_gbm(S0, mu, sigma, T, n_steps, n_paths)

# Compute exact solution at terminal time
S_exact = S0 * np.exp((mu - 0.5 * sigma**2) * T + sigma * W[-1, :])

# Calculate statistics for comparison
euler_mean = np.mean(S_euler[-1, :])
exact_mean = np.mean(S_exact)
euler_std = np.std(S_euler[-1, :])
exact_std = np.std(S_exact)
euler_median = np.median(S_euler[-1, :])
exact_median = np.median(S_exact)
euler_p5 = np.percentile(S_euler[-1, :], 5)
exact_p5 = np.percentile(S_exact, 5)
euler_p95 = np.percentile(S_euler[-1, :], 95)
exact_p95 = np.percentile(S_exact, 95)

Out[8]:

Console

Comparison of Euler-Maruyama vs Exact GBM Solution at T=1
============================================================
Statistic                Euler-Maruyama              Exact
------------------------------------------------------------
Mean................           110.3381           110.3435
Std Dev.............            22.1806            22.1897
Median..............           108.2752           108.2616
5th Percentile......            77.9855            77.9685
95th Percentile.....           150.0485           150.0308

The Euler-Maruyama approximation closely matches the exact solution, with small differences due to the discretization error.

Out[9]:

Visualization

Fan of 50 muted GBM price paths expanding over one year, with an emphasized mean path and dashed expected-value curve. — Geometric Brownian motion sample paths simulated via Euler-Maruyama. The muted ensemble expands over time, illustrating increasing uncertainty while remaining strictly positive. The emphasized mean path tracks the dashed theoretical expected value, confirming the log-normal growth dynamics.

Out[10]:

Visualization

Filled Euler-Maruyama histogram and outlined exact-solution histogram closely overlapping across terminal GBM prices. — Terminal price distribution at T=1 comparing Euler-Maruyama simulation with the exact analytical solution. The histogram estimates overlap significantly, validating that the discretized simulation correctly captures the log-normal distribution predicted by the exact solution derived via Itô's Lemma.

Verifying the Log-Return DistributionLink Copied

Itô's Lemma tells us that $\ln(S_T/S_0) \sim \mathcal{N}((\mu - \frac{1}{2}\sigma^2)T, \sigma^2 T)$ .

In[11]:

Code

import numpy as np
from scipy import stats

# Compute log returns from exact solution
log_returns = np.log(S_exact / S0)

# Theoretical parameters
theoretical_mean = (mu - 0.5 * sigma**2) * T
theoretical_std = sigma * np.sqrt(T)

# Simulated statistics
sim_mean = np.mean(log_returns)
sim_std = np.std(log_returns)

# Normality test
_, p_value = stats.shapiro(np.random.choice(log_returns, 5000, replace=False))

Out[12]:

Console

Log-Return Distribution Analysis
==================================================

Theoretical (from Itô's Lemma):
  Mean:     0.080000
  Std Dev:  0.200000

Simulated:
  Mean:     0.078533
  Std Dev:  0.199612

Shapiro-Wilk normality test p-value: 0.7918
  (High p-value supports normality)

The simulated log-returns match the theoretical distribution predicted by Itô's Lemma, and the Shapiro-Wilk test confirms normality.

Out[13]:

Visualization

Histogram of simulated log-returns with a theoretical normal density curve overlaid, showing close agreement around the theoretical mean. — Distribution of log-returns for geometric Brownian motion. The histogram matches the theoretical normal density curve, confirming that log-returns are normally distributed as predicted by the transformation from Itô's Lemma.

Out[14]:

Visualization

Q-Q plot with simulated quantiles aligned closely along the diagonal reference line, confirming normality of GBM log-returns. — Q-Q plot of log-returns illustrating normality. The close alignment with the diagonal reference line validates that the simulated returns follow the theoretical normal distribution derived from Itô's Lemma.

Visualizing the Itô CorrectionLink Copied

The Itô correction $-\frac{1}{2}\sigma^2$ has a significant impact, especially for high volatility. Let's visualize this.

In[15]:

Code

import numpy as np

sigmas = np.linspace(0, 0.8, 100)
mu_fixed = 0.10

# Expected simple return (approximately mu*T for small dt)
expected_simple_return = mu_fixed

# Expected log return
expected_log_return = mu_fixed - 0.5 * sigmas**2

Out[16]:

Visualization

Line chart showing expected log-return decreasing quadratically as volatility increases while expected return stays constant. — Impact of the Itô correction on expected returns. As volatility increases, the gap between the expected simple return (blue) and the expected log-return (red) widens due to the convexity adjustment $-\frac{1}{2}\sigma^2$. This demonstrates why highly volatile assets can have positive expected simple returns but negative expected compound growth rates.

This graph illustrates why high-volatility assets can have positive expected returns but negative expected log-returns. For example, with $\mu = 10\%$ and $\sigma = 40\%$ , the expected log-return is only $10\% - \frac{1}{2}(40\%)^2 = 2\%$ . At $\sigma \approx 45\%$ , the expected log-return becomes negative despite the positive drift! This phenomenon has large implications for long-term investing and helps explain why highly volatile assets can be poor long-term investments even when their expected single-period returns are positive.

Verifying the Exponential MartingaleLink Copied

We showed that $Y_t = e^{\sigma W_t - \frac{1}{2}\sigma^2 t}$ is a martingale. Let's verify that $\mathbb{E}[Y_T] = Y_0 = 1$ .

In[17]:

Code

import numpy as np

# Simulate the exponential martingale
sigma_m = 0.3
T_m = 2.0
n_paths_m = 50000

# Generate terminal Brownian motion values
W_T = np.random.normal(0, np.sqrt(T_m), n_paths_m)

# Compute Y_T
Y_T = np.exp(sigma_m * W_T - 0.5 * sigma_m**2 * T_m)

# Also compute without the correction (not a martingale)
Y_T_wrong = np.exp(sigma_m * W_T)

# Compute statistics
mean_Y_T = np.mean(Y_T)
mean_Y_T_wrong = np.mean(Y_T_wrong)
theoretical_wrong = np.exp(0.5 * sigma_m**2 * T_m)

Out[18]:

Console

Exponential Martingale Verification
==================================================

With Itô correction (true martingale):
  E[Y_T] = 0.998446  (should be 1.0)

Without correction (not a martingale):
  E[exp(σW_T)] = 1.092474
  Theoretical: exp(½σ²T) = 1.094174

The exponential martingale with the correction has expected value 1, while the uncorrected version has expected value $e^{\frac{1}{2}\sigma^2 T}$ . This difference is precisely the Itô correction at work. The correction term in the exponent exactly compensates for the convexity of the exponential function, so the process has no drift.

Out[19]:

Visualization

Line chart of 30 exponential martingale paths fluctuating around the red dashed line at 1, with a black mean path confirming zero drift. — Exponential martingale $e^{\sigma W_t - \frac{1}{2}\sigma^2 t}$ paths over 200 time steps. The process fluctuates around 1 (red dashed line), with the mean path confirming zero drift — the Itô correction exactly cancels the convexity effect of the exponential function.

Line chart of 30 uncorrected exponential paths trending upward, with a black mean path and red dashed theoretical expectation curve rising above 1. — Uncorrected exponential process $e^{\sigma W_t}$ paths over 200 time steps. Without the Itô correction, the process drifts upward and the mean path follows the theoretical expectation $e^{\frac{1}{2}\sigma^2 t}$ (red dashed line).