Stochastic Calculus for Quantitative Finance

1import numpy as np
2import matplotlib.pyplot as plt
3
4def simulate_brownian_motion(T: float, n_steps: int, n_paths: int = 1) -> np.ndarray:
5    """
6    Simulate standard Brownian motion paths.
7    
8    Args:
9        T: Time horizon
10        n_steps: Number of time steps
11        n_paths: Number of paths to simulate
12        
13    Returns:
14        Array of shape (n_paths, n_steps + 1) containing paths
15    """
16    dt = T / n_steps
17    # Generate increments: dW ~ N(0, sqrt(dt))
18    dW = np.random.normal(0, np.sqrt(dt), size=(n_paths, n_steps))
19    # Cumulative sum to get paths, prepend zeros for W_0 = 0
20    W = np.concatenate([np.zeros((n_paths, 1)), np.cumsum(dW, axis=1)], axis=1)
21    return W
22
23# Simulate and plot
24T, n_steps = 1.0, 1000
25paths = simulate_brownian_motion(T, n_steps, n_paths=5)
26t = np.linspace(0, T, n_steps + 1)
27
28plt.figure(figsize=(12, 6))
29for i in range(5):
30    plt.plot(t, paths[i], alpha=0.7, label=f'Path {i+1}')
31plt.xlabel('Time')
32plt.ylabel('W(t)')
33plt.title('Standard Brownian Motion Paths')
34plt.legend()
35plt.grid(True, alpha=0.3)
36plt.show()
37

1def euler_maruyama(
2    x0: float,
3    mu: callable,
4    sigma: callable,
5    T: float,
6    n_steps: int,
7    n_paths: int = 1
8) -> np.ndarray:
9    """
10    Euler-Maruyama discretization for general SDEs.
11    
12    Args:
13        x0: Initial value
14        mu: Drift function mu(x, t)
15        sigma: Diffusion function sigma(x, t)
16        T: Time horizon
17        n_steps: Number of time steps
18        n_paths: Number of paths
19        
20    Returns:
21        Array of shape (n_paths, n_steps + 1)
22    """
23    dt = T / n_steps
24    sqrt_dt = np.sqrt(dt)
25    
26    X = np.zeros((n_paths, n_steps + 1))
27    X[:, 0] = x0
28    
29    t = np.linspace(0, T, n_steps + 1)
30    
31    for i in range(n_steps):
32        dW = np.random.normal(0, 1, n_paths) * sqrt_dt
33        X[:, i+1] = X[:, i] + mu(X[:, i], t[i]) * dt + sigma(X[:, i], t[i]) * dW
34    
35    return X
36
37# Example: Ornstein-Uhlenbeck process (mean-reverting)
38# dX_t = theta * (mu - X_t) dt + sigma dW_t
39theta, mu_ou, sigma_ou = 0.5, 1.0, 0.3
40
41def drift_ou(x, t):
42    return theta * (mu_ou - x)
43
44def diffusion_ou(x, t):
45    return sigma_ou * np.ones_like(x)
46
47paths_ou = euler_maruyama(0.5, drift_ou, diffusion_ou, T=5.0, n_steps=1000, n_paths=100)
48
49# Plot mean and confidence bands
50t = np.linspace(0, 5.0, 1001)
51mean_path = np.mean(paths_ou, axis=0)
52std_path = np.std(paths_ou, axis=0)
53
54plt.figure(figsize=(12, 6))
55plt.plot(t, mean_path, 'b-', linewidth=2, label='Mean')
56plt.fill_between(t, mean_path - 2*std_path, mean_path + 2*std_path, 
57                 alpha=0.3, label='95% CI')
58plt.axhline(y=mu_ou, color='r', linestyle='--', label='Long-term mean')
59plt.xlabel('Time')
60plt.ylabel('X(t)')
61plt.title('Ornstein-Uhlenbeck Process (Mean Reversion)')
62plt.legend()
63plt.grid(True, alpha=0.3)
64plt.show()
65

1def verify_ito_lemma():
2    """
3    Verify Itô's lemma numerically for f(X_t) = X_t^2 where dX_t = mu dt + sigma dW_t
4    """
5    mu, sigma = 0.1, 0.3
6    X0 = 1.0
7    T, n_steps = 1.0, 10000
8    dt = T / n_steps
9    sqrt_dt = np.sqrt(dt)
10    
11    # Simulate X_t
12    X = np.zeros(n_steps + 1)
13    X[0] = X0
14    
15    for i in range(n_steps):
16        dW = np.random.normal(0, sqrt_dt)
17        X[i+1] = X[i] + mu * dt + sigma * dW
18    
19    # Direct computation: f(X_T) = X_T^2
20    f_direct = X[-1]**2
21    
22    # Itô's lemma: df = (2*mu*X + sigma^2) dt + 2*sigma*X dW
23    # Integrated: f(X_T) = f(X_0) + integral of drift + integral of diffusion
24    f_ito = X0**2
25    for i in range(n_steps):
26        dW = (X[i+1] - X[i] - mu * dt) / sigma  # Recover dW
27        df_drift = (2 * mu * X[i] + sigma**2) * dt
28        df_diffusion = 2 * sigma * X[i] * dW
29        f_ito += df_drift + df_diffusion
30    
31    print(f"Direct computation: {f_direct:.6f}")
32    print(f"Itô's lemma: {f_ito:.6f}")
33    print(f"Difference: {abs(f_direct - f_ito):.6e}")
34    
35verify_ito_lemma()
36# Output: Difference < 1e-10 (numerical precision)
37

1def european_call_monte_carlo(
2    S0: float,
3    K: float,
4    r: float,
5    sigma: float,
6    T: float,
7    n_sims: int = 100000
8) -> tuple[float, float]:
9    """
10    Price European call option using Monte Carlo (Feynman-Kac).
11    
12    Returns:
13        (price, standard_error)
14    """
15    # Simulate terminal stock prices under risk-neutral measure
16    # S_T = S_0 * exp((r - 0.5*sigma^2)*T + sigma*sqrt(T)*Z)
17    Z = np.random.normal(0, 1, n_sims)
18    ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
19    
20    # Payoff: max(S_T - K, 0)
21    payoffs = np.maximum(ST - K, 0)
22    
23    # Discount to present value
24    price = np.exp(-r * T) * np.mean(payoffs)
25    std_error = np.exp(-r * T) * np.std(payoffs) / np.sqrt(n_sims)
26    
27    return price, std_error
28
29# Example
30S0, K, r, sigma, T = 100, 100, 0.05, 0.2, 1.0
31price_mc, se = european_call_monte_carlo(S0, K, r, sigma, T)
32
33# Compare with Black-Scholes analytical formula
34from scipy.stats import norm
35
36def black_scholes_call(S0, K, r, sigma, T):
37    d1 = (np.log(S0/K) + (r + 0.5*sigma**2)*T) / (sigma*np.sqrt(T))
38    d2 = d1 - sigma*np.sqrt(T)
39    return S0 * norm.cdf(d1) - K * np.exp(-r*T) * norm.cdf(d2)
40
41price_bs = black_scholes_call(S0, K, r, sigma, T)
42
43print(f"Monte Carlo: {price_mc:.4f} ± {se:.4f}")
44print(f"Black-Scholes: {price_bs:.4f}")
45print(f"Difference: {abs(price_mc - price_bs):.4f}")
46

1def simulate_measure_change(
2    S0: float,
3    mu: float,
4    r: float,
5    sigma: float,
6    T: float,
7    n_steps: int,
8    n_paths: int = 10000
9) -> tuple[np.ndarray, np.ndarray]:
10    """
11    Simulate stock paths under both physical and risk-neutral measures.
12    
13    Returns:
14        (paths_P, paths_Q): Paths under P and Q measures
15    """
16    dt = T / n_steps
17    sqrt_dt = np.sqrt(dt)
18    
19    # Simulate under physical measure P
20    paths_P = np.zeros((n_paths, n_steps + 1))
21    paths_P[:, 0] = S0
22    
23    for i in range(n_steps):
24        dW_P = np.random.normal(0, sqrt_dt, n_paths)
25        paths_P[:, i+1] = paths_P[:, i] * (1 + mu * dt + sigma * dW_P)
26    
27    # Transform to risk-neutral measure Q
28    # W^Q = W^P + theta * t, where theta = (mu - r) / sigma
29    theta = (mu - r) / sigma
30    paths_Q = np.zeros((n_paths, n_steps + 1))
31    paths_Q[:, 0] = S0
32    
33    for i in range(n_steps):
34        # Under Q: dS = r*S*dt + sigma*S*dW^Q
35        dW_Q = np.random.normal(0, sqrt_dt, n_paths)
36        paths_Q[:, i+1] = paths_Q[:, i] * (1 + r * dt + sigma * dW_Q)
37    
38    return paths_P, paths_Q
39
40# Simulate
41S0, mu, r, sigma, T = 100, 0.15, 0.05, 0.2, 1.0
42paths_P, paths_Q = simulate_measure_change(S0, mu, r, sigma, T, 252, 10000)
43
44# Compare terminal distributions
45print(f"Under P: E[S_T] = {np.mean(paths_P[:, -1]):.2f} (theoretical: {S0 * np.exp(mu * T):.2f})")
46print(f"Under Q: E[S_T] = {np.mean(paths_Q[:, -1]):.2f} (theoretical: {S0 * np.exp(r * T):.2f})")
47

1def asian_call_monte_carlo(
2    S0: float,
3    K: float,
4    r: float,
5    sigma: float,
6    T: float,
7    n_steps: int,
8    n_sims: int = 100000
9) -> float:
10    """Price arithmetic Asian call option."""
11    dt = T / n_steps
12    sqrt_dt = np.sqrt(dt)
13    
14    payoffs = np.zeros(n_sims)
15    
16    for sim in range(n_sims):
17        S = np.zeros(n_steps + 1)
18        S[0] = S0
19        
20        for i in range(n_steps):
21            Z = np.random.normal()
22            S[i+1] = S[i] * np.exp((r - 0.5*sigma**2)*dt + sigma*sqrt_dt*Z)
23        
24        avg_price = np.mean(S)
25        payoffs[sim] = max(avg_price - K, 0)
26    
27    return np.exp(-r * T) * np.mean(payoffs)
28
29price_asian = asian_call_monte_carlo(100, 100, 0.05, 0.2, 1.0, 252)
30print(f"Asian call price: {price_asian:.4f}")
31

1def barrier_call_monte_carlo(
2    S0: float,
3    K: float,
4    B: float,
5    r: float,
6    sigma: float,
7    T: float,
8    n_steps: int,
9    n_sims: int = 100000
10) -> float:
11    """Price up-and-out barrier call option."""
12    dt = T / n_steps
13    sqrt_dt = np.sqrt(dt)
14    
15    payoffs = np.zeros(n_sims)
16    
17    for sim in range(n_sims):
18        S = S0
19        knocked_out = False
20        
21        for _ in range(n_steps):
22            Z = np.random.normal()
23            S = S * np.exp((r - 0.5*sigma**2)*dt + sigma*sqrt_dt*Z)
24            
25            if S >= B:
26                knocked_out = True
27                break
28        
29        if not knocked_out:
30            payoffs[sim] = max(S - K, 0)
31    
32    return np.exp(-r * T) * np.mean(payoffs)
33
34price_barrier = barrier_call_monte_carlo(100, 100, 120, 0.05, 0.2, 1.0, 252)
35print(f"Barrier call price: {price_barrier:.4f}")
36

1def black_scholes_pde_solver(
2    S0: float,
3    K: float,
4    r: float,
5    sigma: float,
6    T: float,
7    S_max: float = 200,
8    n_S: int = 100,
9    n_t: int = 100
10) -> float:
11    """
12    Solve Black-Scholes PDE using explicit finite difference.
13    
14    PDE: dV/dt + 0.5*sigma^2*S^2*d^2V/dS^2 + r*S*dV/dS - r*V = 0
15    """
16    # Grid setup
17    dS = S_max / n_S
18    dt = T / n_t
19    
20    S = np.linspace(0, S_max, n_S + 1)
21    V = np.zeros((n_t + 1, n_S + 1))
22    
23    # Terminal condition: V(S, T) = max(S - K, 0)
24    V[-1, :] = np.maximum(S - K, 0)
25    
26    # Boundary conditions
27    # V(0, t) = 0 (call worthless if S=0)
28    # V(S_max, t) = S_max - K*exp(-r*(T-t)) (deep in the money)
29    
30    # Coefficients for explicit scheme
31    for j in range(n_t - 1, -1, -1):
32        for i in range(1, n_S):
33            # Finite difference approximations
34            dV_dS = (V[j+1, i+1] - V[j+1, i-1]) / (2 * dS)
35            d2V_dS2 = (V[j+1, i+1] - 2*V[j+1, i] + V[j+1, i-1]) / (dS**2)
36            
37            # Explicit update
38            V[j, i] = V[j+1, i] - dt * (
39                0.5 * sigma**2 * S[i]**2 * d2V_dS2 +
40                r * S[i] * dV_dS -
41                r * V[j+1, i]
42            )
43        
44        # Boundary conditions
45        V[j, 0] = 0
46        V[j, -1] = S_max - K * np.exp(-r * (T - j*dt))
47    
48    # Interpolate to get price at S0
49    idx = np.searchsorted(S, S0)
50    if idx == 0 or idx == len(S):
51        return V[0, idx]
52    
53    # Linear interpolation
54    w = (S0 - S[idx-1]) / (S[idx] - S[idx-1])
55    return (1 - w) * V[0, idx-1] + w * V[0, idx]
56
57price_pde = black_scholes_pde_solver(100, 100, 0.05, 0.2, 1.0)
58print(f"PDE solver price: {price_pde:.4f}")
59print(f"Black-Scholes analytical: {black_scholes_call(100, 100, 0.05, 0.2, 1.0):.4f}")
60

1#include <vector>
2#include <cmath>
3#include <random>
4#include <algorithm>
5
6class SDESolver {
7private:
8    std::mt19937 rng_;
9    std::normal_distribution<double> normal_;
10    
11public:
12    SDESolver(unsigned seed = 42) : rng_(seed), normal_(0.0, 1.0) {}
13    
14    // Simulate GBM path
15    std::vector<double> simulateGBM(
16        double S0,
17        double mu,
18        double sigma,
19        double T,
20        int n_steps
21    ) {
22        std::vector<double> path(n_steps + 1);
23        path[0] = S0;
24        
25        double dt = T / n_steps;
26        double sqrt_dt = std::sqrt(dt);
27        
28        for (int i = 0; i < n_steps; ++i) {
29            double dW = normal_(rng_) * sqrt_dt;
30            path[i+1] = path[i] * std::exp(
31                (mu - 0.5 * sigma * sigma) * dt + sigma * dW
32            );
33        }
34        
35        return path;
36    }
37    
38    // European call option pricing via Monte Carlo
39    std::pair<double, double> europeanCallMC(
40        double S0,
41        double K,
42        double r,
43        double sigma,
44        double T,
45        int n_sims
46    ) {
47        double sum = 0.0;
48        double sum_sq = 0.0;
49        
50        double sqrt_T = std::sqrt(T);
51        double drift = (r - 0.5 * sigma * sigma) * T;
52        
53        for (int i = 0; i < n_sims; ++i) {
54            double Z = normal_(rng_);
55            double ST = S0 * std::exp(drift + sigma * sqrt_T * Z);
56            double payoff = std::max(ST - K, 0.0);
57            
58            sum += payoff;
59            sum_sq += payoff * payoff;
60        }
61        
62        double mean_payoff = sum / n_sims;
63        double variance = (sum_sq / n_sims) - (mean_payoff * mean_payoff);
64        double std_error = std::sqrt(variance / n_sims);
65        
66        double price = std::exp(-r * T) * mean_payoff;
67        double se = std::exp(-r * T) * std_error;
68        
69        return {price, se};
70    }
71};
72
73// Usage
74int main() {
75    SDESolver solver;
76    auto [price, se] = solver.europeanCallMC(100.0, 100.0, 0.05, 0.2, 1.0, 1000000);
77    
78    std::cout << "Price: " << price << " ± " << se << std::endl;
79    return 0;
80}
81

1def european_call_antithetic(S0, K, r, sigma, T, n_sims):
2    """Monte Carlo with antithetic variates."""
3    Z = np.random.normal(0, 1, n_sims // 2)
4    
5    # Positive and negative samples
6    ST_pos = S0 * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*Z)
7    ST_neg = S0 * np.exp((r - 0.5*sigma**2)*T - sigma*np.sqrt(T)*Z)
8    
9    payoffs = (np.maximum(ST_pos - K, 0) + np.maximum(ST_neg - K, 0)) / 2
10    
11    return np.exp(-r * T) * np.mean(payoffs)
12

1def simulate_correlated_gbm(
2    S0: np.ndarray,
3    mu: np.ndarray,
4    sigma: np.ndarray,
5    corr_matrix: np.ndarray,
6    T: float,
7    n_steps: int
8) -> np.ndarray:
9    """
10    Simulate correlated GBM processes.
11    
12    Args:
13        S0: Initial prices (n_assets,)
14        mu: Drifts (n_assets,)
15        sigma: Volatilities (n_assets,)
16        corr_matrix: Correlation matrix (n_assets, n_assets)
17        
18    Returns:
19        Paths of shape (n_steps + 1, n_assets)
20    """
21    n_assets = len(S0)
22    dt = T / n_steps
23    sqrt_dt = np.sqrt(dt)
24    
25    # Cholesky decomposition for correlated normals
26    L = np.linalg.cholesky(corr_matrix)
27    
28    S = np.zeros((n_steps + 1, n_assets))
29    S[0] = S0
30    
31    for i in range(n_steps):
32        # Generate correlated normal samples
33        Z = np.random.normal(0, 1, n_assets)
34        Z_corr = L @ Z
35        
36        for j in range(n_assets):
37            S[i+1, j] = S[i, j] * np.exp(
38                (mu[j] - 0.5*sigma[j]**2)*dt + sigma[j]*sqrt_dt*Z_corr[j]
39            )
40    
41    return S
42