Chapter 4: Stochastic Differential Equations and Wiener Process

Stochastic Differential Equations and Wiener Process

📖 Reading Time: 20-25 minutes 📊 Difficulty: Beginner 💻 Code Examples: 0 📝 Exercises: 0

🌐 EN | 🇯🇵 JP | Last sync: 2025-11-16

Fundamentals of Mathematics Dojo > Probability Theory and Stochastic Processes > Chapter 4

4.1 Wiener Process (Brownian Motion)

📐 Definition: Wiener Process (Standard Brownian Motion)
A stochastic process \(\{W(t)\}_{t \geq 0}\) is called a Wiener process if:

  1. \(W(0) = 0\) (starts from zero)
  2. Independent increments: \(W(t_1), W(t_2) - W(t_1), W(t_3) - W(t_2), \ldots\) are independent
  3. Stationary increments: \(W(t+s) - W(s) \sim N(0, t)\)
  4. Continuous paths: \(W(t)\) is a continuous function in \(t\)

Properties: \[E[W(t)] = 0, \quad Var(W(t)) = t, \quad Cov(W(s), W(t)) = \min(s, t)\]

💻 Code Example 1: Simulation of Wiener Process (Brownian Motion)

import numpy as np import matplotlib.pyplot as plt from scipy import stats # Wiener process simulation np.random.seed(42) def simulate_wiener_process(T, dt): """Wiener process simulation""" n_steps = int(T / dt) t = np.linspace(0, T, n_steps + 1) # Increments dW ~ N(0, dt) dW = np.random.normal(0, np.sqrt(dt), n_steps) # Build Wiener process by cumulative sum W = np.concatenate([[0], np.cumsum(dW)]) return t, W # Parameters T = 10 # End time dt = 0.01 # Time step # Simulate multiple paths n_paths = 10 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # (1) Sample paths for _ in range(n_paths): t, W = simulate_wiener_process(T, dt) axes[0, 0].plot(t, W, alpha=0.6, linewidth=1.5) axes[0, 0].axhline(y=0, color=‘red’, linestyle=’—’, linewidth=2, label=‘E[W(t)] = 0’) axes[0, 0].set_xlabel(‘Time t’, fontsize=11) axes[0, 0].set_ylabel(‘W(t)’, fontsize=11) axes[0, 0].set_title(‘Sample Paths of Wiener Process’, fontsize=12, fontweight=‘bold’) axes[0, 0].legend() axes[0, 0].grid(alpha=0.3) # (2) Variance evolution in time n_simulations = 1000 t_checkpoints = [1, 2, 5, 10] for t_check in t_checkpoints: W_values = [] for _ in range(n_simulations): t, W = simulate_wiener_process(t_check, dt) W_values.append(W[-1]) W_values = np.array(W_values) # Histogram axes[0, 1].hist(W_values, bins=50, alpha=0.4, density=True, label=f’t={t_check}’) # Theoretical distribution N(0, t) x = np.linspace(-6, 6, 1000) for t_check in t_checkpoints: theoretical_pdf = stats.norm.pdf(x, 0, np.sqrt(t_check)) axes[0, 1].plot(x, theoretical_pdf, linewidth=2, linestyle=’—’) axes[0, 1].set_xlabel(‘W(t)’, fontsize=11) axes[0, 1].set_ylabel(‘Density’, fontsize=11) axes[0, 1].set_title(‘Distribution of W(t) at Different Times’, fontsize=12, fontweight=‘bold’) axes[0, 1].legend() axes[0, 1].grid(alpha=0.3) # (3) Evolution of mean and variance over time t, W_sample = simulate_wiener_process(T, dt) n_paths_analysis = 500 all_paths = [] for _ in range(n_paths_analysis): t, W = simulate_wiener_process(T, dt) all_paths.append(W) all_paths = np.array(all_paths) mean_W = all_paths.mean(axis=0) std_W = all_paths.std(axis=0) axes[1, 0].plot(t, mean_W, color=‘#667eea’, linewidth=2.5, label=‘Observed Mean’) axes[1, 0].fill_between(t, mean_W - std_W, mean_W + std_W, alpha=0.3, color=‘#667eea’, label=‘Observed ±1σ’) axes[1, 0].plot(t, np.zeros_like(t), ‘r—’, linewidth=2, label=‘Theoretical Mean=0’) axes[1, 0].plot(t, np.sqrt(t), ‘g—’, linewidth=2, label=‘Theoretical σ=√t’) axes[1, 0].plot(t, -np.sqrt(t), ‘g—’, linewidth=2) axes[1, 0].set_xlabel(‘Time t’, fontsize=11) axes[1, 0].set_ylabel(‘W(t)’, fontsize=11) axes[1, 0].set_title(‘Evolution of Mean and Variance over Time’, fontsize=12, fontweight=‘bold’) axes[1, 0].legend() axes[1, 0].grid(alpha=0.3) # (4) Covariance structure Cov(W(s), W(t)) = min(s, t) s_values = np.linspace(0, T, 100) t_values = np.linspace(0, T, 100) S, T_grid = np.meshgrid(s_values, t_values) # Theoretical covariance Cov_theory = np.minimum(S, T_grid) im = axes[1, 1].imshow(Cov_theory, extent=[0, T, 0, T], origin=‘lower’, cmap=‘viridis’, aspect=‘auto’) axes[1, 1].set_xlabel(‘Time s’, fontsize=11) axes[1, 1].set_ylabel(‘Time t’, fontsize=11) axes[1, 1].set_title(‘Cov(W(s), W(t)) = min(s, t)’, fontsize=12, fontweight=‘bold’) plt.colorbar(im, ax=axes[1, 1]) plt.tight_layout() plt.show() print(“Verification of Wiener Process Properties:”) print(”=“*60) for t_check in t_checkpoints: W_values = [] for _ in range(n_simulations): t, W = simulate_wiener_process(t_check, dt) W_values.append(W[-1]) W_values = np.array(W_values) print(f”t={t_check}:”) print(f” Theoretical: E[W(t)]={0:.4f}, Var(W(t))={t_check:.4f}”) print(f” Observed: E[W(t)]={W_values.mean():.4f}, Var(W(t))={W_values.var():.4f}“)

4.2 Drift and Volatility

📐 Definition: Wiener Process with Drift
A Wiener process with drift coefficient \(\mu\) and volatility coefficient \(\sigma\): \[X(t) = X_0 + \mu t + \sigma W(t)\] Expectation and Variance: \[E[X(t)] = X_0 + \mu t, \quad Var(X(t)) = \sigma^2 t\] Drift \(\mu\) represents a deterministic trend, while volatility \(\sigma\) represents the magnitude of random fluctuations.

💻 Code Example 2: Effects of Drift and Volatility

Simulation of Wiener process with drift np.random.seed(42) def simulate_drift_wiener(X0, mu, sigma, T, dt): """Simulation of Wiener process with drift""" t, W = simulate_wiener_process(T, dt) X = X0 + mu * t + sigma * W return t, X # Parameter settings X0 = 100 T = 5 dt = 0.01 # Different combinations of drift and volatility configs = [ {‘mu’: 0, ‘sigma’: 1, ‘label’: ‘μ=0, σ=1 (Standard Wiener)’}, {‘mu’: 2, ‘sigma’: 1, ‘label’: ‘μ=2, σ=1 (Positive Drift)’}, {‘mu’: -1, ‘sigma’: 1, ‘label’: ‘μ=-1, σ=1 (Negative Drift)’}, {‘mu’: 0, ‘sigma’: 3, ‘label’: ‘μ=0, σ=3 (High Volatility)’}, ] fig, axes = plt.subplots(2, 2, figsize=(14, 10)) axes = axes.flatten() for idx, config in enumerate(configs): mu = config[‘mu’] sigma = config[‘sigma’] # Simulate multiple paths for _ in range(20): t, X = simulate_drift_wiener(X0, mu, sigma, T, dt) axes[idx].plot(t, X, alpha=0.4, linewidth=1.5, color=‘#667eea’) # Expectation t_theory = np.linspace(0, T, 1000) E_X = X0 + mu * t_theory axes[idx].plot(t_theory, E_X, ‘r—’, linewidth=2.5, label=f’E[X(t)]={X0}+{mu}t’) # ±1σ confidence band std_X = sigma * np.sqrt(t_theory) axes[idx].fill_between(t_theory, E_X - std_X, E_X + std_X, alpha=0.2, color=‘red’, label=‘E[X]±σ’) axes[idx].set_xlabel(‘Time t’, fontsize=11) axes[idx].set_ylabel(‘X(t)’, fontsize=11) axes[idx].set_title(config[‘label’], fontsize=12, fontweight=‘bold’) axes[idx].legend() axes[idx].grid(alpha=0.3) plt.tight_layout() plt.show() # Statistical verification print(“\nStatistical Verification of Wiener Process with Drift:”) print(”=“60) for config in configs: mu = config[‘mu’] sigma = config[‘sigma’] # Final values of multiple paths n_sims = 1000 final_values = [] for _ in range(n_sims): t, X = simulate_drift_wiener(X0, mu, sigma, T, dt) final_values.append(X[-1]) final_values = np.array(final_values) print(f”\n{config[‘label’]}:”) print(f” Theoretical: E[X({T})]={X0 + muT:.2f}, Var(X({T}))={sigma**2*T:.2f}”) print(f” Observed: E[X({T})]={final_values.mean():.2f}, Var(X({T}))={final_values.var():.2f}“)

4.3 Stochastic Differential Equations (SDE) and Euler-Maruyama Method

📐 Definition: Stochastic Differential Equation (SDE)
General form of a stochastic differential equation: \[dX(t) = \mu(X, t) , dt + \sigma(X, t) , dW(t)\] Where:

Euler-Maruyama Method (numerical solution): \[X_{n+1} = X_n + \mu(X_n, t_n) \Delta t + \sigma(X_n, t_n) \sqrt{\Delta t} , Z_n\] where \(Z_n \sim N(0, 1)\) are independent standard normal random variables.

💻 Code Example 3: SDE Numerical Solution by Euler-Maruyama Method

Implementation of Euler-Maruyama method def euler_maruyama(mu_func, sigma_func, X0, T, dt): """ Numerical solution of SDE by Euler-Maruyama method Parameters: ----------- mu_func : function Drift term μ(X, t) sigma_func : function Diffusion term σ(X, t) X0 : float Initial value T : float End time dt : float Time step """ n_steps = int(T / dt) t = np.linspace(0, T, n_steps + 1) X = np.zeros(n_steps + 1) X[0] = X0 for i in range(n_steps): dW = np.random.normal(0, np.sqrt(dt)) X[i+1] = X[i] + mu_func(X[i], t[i]) * dt + sigma_func(X[i], t[i]) * dW return t, X # Example 1: Linear SDE dX = -θX dt + σ dW (special case of Ornstein-Uhlenbeck process) theta = 0.5 sigma = 1.0 mu_func = lambda X, t: -theta * X sigma_func = lambda X, t: sigma np.random.seed(42) X0 = 5 T = 10 dt = 0.01 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # (1) Simulation of multiple paths for _ in range(20): t, X = euler_maruyama(mu_func, sigma_func, X0, T, dt) axes[0, 0].plot(t, X, alpha=0.4, linewidth=1.5, color=‘#667eea’) axes[0, 0].axhline(y=0, color=‘red’, linestyle=’—’, linewidth=2, label=‘Equilibrium’) axes[0, 0].set_xlabel(‘Time t’, fontsize=11) axes[0, 0].set_ylabel(‘X(t)’, fontsize=11) axes[0, 0].set_title(‘Linear SDE: dX = -θX dt + σ dW’, fontsize=12, fontweight=‘bold’) axes[0, 0].legend() axes[0, 0].grid(alpha=0.3) # (2) Comparison of different time steps dts = [0.1, 0.01, 0.001] np.random.seed(42) for dt_test in dts: t, X = euler_maruyama(mu_func, sigma_func, X0, T, dt_test) axes[0, 1].plot(t, X, linewidth=2, alpha=0.7, label=f’Δt={dt_test}’) axes[0, 1].set_xlabel(‘Time t’, fontsize=11) axes[0, 1].set_ylabel(‘X(t)’, fontsize=11) axes[0, 1].set_title(‘Effect of Time Step’, fontsize=12, fontweight=‘bold’) axes[0, 1].legend() axes[0, 1].grid(alpha=0.3) # (3) Convergence to mean (mean reversion) n_paths = 500 all_paths = [] for _ in range(n_paths): t, X = euler_maruyama(mu_func, sigma_func, X0, T, dt) all_paths.append(X) all_paths = np.array(all_paths) mean_X = all_paths.mean(axis=0) std_X = all_paths.std(axis=0) axes[1, 0].plot(t, mean_X, color=‘#667eea’, linewidth=2.5, label=‘Observed Mean’) axes[1, 0].fill_between(t, mean_X - std_X, mean_X + std_X, alpha=0.3, color=‘#667eea’, label=‘±1σ’) # Theoretical mean (analytical solution): E[X(t)] = X0 * exp(-θt) E_X_theory = X0 * np.exp(-theta * t) axes[1, 0].plot(t, E_X_theory, ‘r—’, linewidth=2, label=‘Theoretical Mean’) axes[1, 0].set_xlabel(‘Time t’, fontsize=11) axes[1, 0].set_ylabel(‘E[X(t)]’, fontsize=11) axes[1, 0].set_title(‘Convergence to Mean (Mean Reversion)’, fontsize=12, fontweight=‘bold’) axes[1, 0].legend() axes[1, 0].grid(alpha=0.3) # (4) Stationary distribution # Distribution after sufficient time (t→∞ gives N(0, σ²/(2θ))) final_values = all_paths[:, -1] axes[1, 1].hist(final_values, bins=50, density=True, alpha=0.6, color=‘#667eea’, edgecolor=‘black’, label=‘Observed’) # Theoretical stationary distribution x = np.linspace(final_values.min(), final_values.max(), 1000) stationary_var = sigma**2 / (2 * theta) stationary_pdf = stats.norm.pdf(x, 0, np.sqrt(stationary_var)) axes[1, 1].plot(x, stationary_pdf, ‘r-’, linewidth=2.5, label=f’N(0, {stationary_var:.2f})’) axes[1, 1].set_xlabel(‘X(∞)’, fontsize=11) axes[1, 1].set_ylabel(‘Density’, fontsize=11) axes[1, 1].set_title(‘Stationary Distribution’, fontsize=12, fontweight=‘bold’) axes[1, 1].legend() axes[1, 1].grid(alpha=0.3) plt.tight_layout() plt.show() print(“Verification of Linear SDE:”) print(f”Theoretical stationary variance: σ²/(2θ) = {stationary_var:.4f}”) print(f”Observed stationary variance: {final_values.var():.4f}“)

4.4 Geometric Brownian Motion (Stock Price Model)

📊 Theorem: Geometric Brownian Motion
Geometric Brownian motion is defined by the following SDE: \[dS(t) = \mu S(t) , dt + \sigma S(t) , dW(t)\] Analytical Solution (by Itô’s formula): \[S(t) = S_0 \exp\left[\left(\mu - \frac{\sigma^2}{2}\right)t + \sigma W(t)\right]\] Expectation: \(E[S(t)] = S_0 e^{\mu t}\)
Variance: \(Var(S(t)) = S_0^2 e^{2\mu t}(e^{\sigma^2 t} - 1)\) Applications: Black-Scholes model (stock prices, exchange rates)

💻 Code Example 4: Geometric Brownian Motion (Stock Price Model)

Simulation of geometric Brownian motion np.random.seed(42) def geometric_brownian_motion(S0, mu, sigma, T, dt): """Geometric Brownian motion simulation (using analytical solution)""" n_steps = int(T / dt) t = np.linspace(0, T, n_steps + 1) # Generate Wiener process dW = np.random.normal(0, np.sqrt(dt), n_steps) W = np.concatenate([[0], np.cumsum(dW)]) # Analytical solution S = S0 * np.exp((mu - 0.5 * sigma2) * t + sigma * W) return t, S # Parameters (stock price model) S0 = 100 # Initial stock price mu = 0.1 # Drift (annual return rate) sigma = 0.2 # Volatility (annual) T = 1 # 1 year dt = 1/252 # 1 trading day fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # (1) Sample paths n_paths = 50 for _ in range(n_paths): t, S = geometric_brownian_motion(S0, mu, sigma, T, dt) axes[0, 0].plot(t, S, alpha=0.4, linewidth=1.5, color=‘#667eea’) # Expectation t_theory = np.linspace(0, T, 1000) E_S = S0 * np.exp(mu * t_theory) axes[0, 0].plot(t_theory, E_S, ‘r—’, linewidth=2.5, label=f’E[S(t)]={S0}e^{{μt}}’) axes[0, 0].set_xlabel(‘Time (years)’, fontsize=11) axes[0, 0].set_ylabel(‘Stock Price S(t)’, fontsize=11) axes[0, 0].set_title(‘Geometric Brownian Motion (Stock Price Model)’, fontsize=12, fontweight=‘bold’) axes[0, 0].legend() axes[0, 0].grid(alpha=0.3) axes[0, 0].set_ylim([0, 200]) # (2) Distribution of log returns # log(S(t)/S0) ~ N((μ - σ²/2)t, σ²t) n_sims = 5000 final_prices = [] for _ in range(n_sims): t, S = geometric_brownian_motion(S0, mu, sigma, T, dt) final_prices.append(S[-1]) final_prices = np.array(final_prices) log_returns = np.log(final_prices / S0) axes[0, 1].hist(log_returns, bins=50, density=True, alpha=0.6, color=‘#667eea’, edgecolor=‘black’, label=‘Simulation’) # Theoretical distribution x = np.linspace(log_returns.min(), log_returns.max(), 1000) theoretical_mean = (mu - 0.5 * sigma2) * T theoretical_std = sigma * np.sqrt(T) theoretical_pdf = stats.norm.pdf(x, theoretical_mean, theoretical_std) axes[0, 1].plot(x, theoretical_pdf, ‘r-’, linewidth=2.5, label=‘Theoretical’) axes[0, 1].set_xlabel(‘Log Return log(S(T)/S₀)’, fontsize=11) axes[0, 1].set_ylabel(‘Density’, fontsize=11) axes[0, 1].set_title(‘Distribution of Log Returns’, fontsize=12, fontweight=‘bold’) axes[0, 1].legend() axes[0, 1].grid(alpha=0.3) # (3) Comparison of different volatilities sigmas = [0.1, 0.2, 0.3, 0.4] colors = [‘blue’, ‘green’, ‘orange’, ‘red’] for sig, color in zip(sigmas, colors): for _ in range(10): t, S = geometric_brownian_motion(S0, mu, sig, T, dt) axes[1, 0].plot(t, S, alpha=0.3, linewidth=1.5, color=color) # Expectation for each volatility axes[1, 0].plot(t_theory, S0 * np.exp(mu * t_theory), ’—’, linewidth=2, color=color, label=f’σ={sig}’) axes[1, 0].set_xlabel(‘Time (years)’, fontsize=11) axes[1, 0].set_ylabel(‘Stock Price S(t)’, fontsize=11) axes[1, 0].set_title(‘Comparison with Different Volatilities’, fontsize=12, fontweight=‘bold’) axes[1, 0].legend() axes[1, 0].grid(alpha=0.3) # (4) Monte Carlo pricing (European call option) K = 100 # Strike price r = 0.05 # Risk-free rate # Option payoff max(S(T) - K, 0) payoffs = np.maximum(final_prices - K, 0) option_price_mc = np.exp(-r * T) * payoffs.mean() # Black-Scholes formula (analytical solution) from scipy.stats import norm as norm_dist d1 = (np.log(S0/K) + (r + 0.5sigma**2)T) / (sigmanp.sqrt(T)) d2 = d1 - sigmanp.sqrt(T) option_price_bs = S0norm_dist.cdf(d1) - Knp.exp(-r*T)*norm_dist.cdf(d2) axes[1, 1].hist(payoffs, bins=50, alpha=0.6, color=‘#667eea’, edgecolor=‘black’) axes[1, 1].axvline(payoffs.mean(), color=‘red’, linestyle=’—’, linewidth=2, label=f’Mean Payoff={payoffs.mean():.2f}’) axes[1, 1].set_xlabel(‘Payoff max(S(T)-K, 0)’, fontsize=11) axes[1, 1].set_ylabel(‘Frequency’, fontsize=11) axes[1, 1].set_title(‘Distribution of Option Payoff’, fontsize=12, fontweight=‘bold’) axes[1, 1].legend() axes[1, 1].grid(alpha=0.3) plt.tight_layout() plt.show() print(“Verification of Geometric Brownian Motion:”) print(f”Theoretical: E[S(T)] = {S0 * np.exp(mu * T):.2f}”) print(f”Observed: E[S(T)] = {final_prices.mean():.2f}”) print(f”\nOption Price:”) print(f”Monte Carlo: {option_price_mc:.4f}”) print(f”Black-Scholes: {option_price_bs:.4f}“)

4.5 Ornstein-Uhlenbeck Process (Mean Reversion)

📐 Definition: Ornstein-Uhlenbeck Process
The Ornstein-Uhlenbeck process is an SDE with mean-reverting properties: \[dX(t) = \theta(\mu - X(t)) , dt + \sigma , dW(t)\] Where:

Analytical Solution: \[X(t) = X_0 e^{-\theta t} + \mu(1 - e^{-\theta t}) + \sigma \int_0^t e^{-\theta(t-s)} dW(s)\] Stationary distribution: \(X(\infty) \sim N\left(\mu, \frac{\sigma^2}{2\theta}\right)\)

💻 Code Example 5: Ornstein-Uhlenbeck Process (Mean Reversion)

Simulation of Ornstein-Uhlenbeck process np.random.seed(42) def ornstein_uhlenbeck(X0, theta, mu, sigma, T, dt): """Ornstein-Uhlenbeck process simulation""" mu_func = lambda X, t: theta * (mu - X) sigma_func = lambda X, t: sigma return euler_maruyama(mu_func, sigma_func, X0, T, dt) # Parameters X0 = 10 theta = 2.0 # Speed of reversion mu = 5 # Long-term mean sigma = 1.5 # Volatility T = 10 dt = 0.01 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # (1) Mean reversion from different initial values initial_values = [0, 3, 5, 7, 10] colors = plt.cm.viridis(np.linspace(0, 1, len(initial_values))) for X0_test, color in zip(initial_values, colors): for _ in range(5): t, X = ornstein_uhlenbeck(X0_test, theta, mu, sigma, T, dt) axes[0, 0].plot(t, X, alpha=0.4, linewidth=1.5, color=color) axes[0, 0].axhline(y=mu, color=‘red’, linestyle=’—’, linewidth=2.5, label=f’Long-term Mean μ={mu}’) axes[0, 0].set_xlabel(‘Time t’, fontsize=11) axes[0, 0].set_ylabel(‘X(t)’, fontsize=11) axes[0, 0].set_title(‘Visualization of Mean Reversion’, fontsize=12, fontweight=‘bold’) axes[0, 0].legend() axes[0, 0].grid(alpha=0.3) # (2) Comparison of different reversion speeds thetas = [0.5, 1.0, 2.0, 5.0] for theta_test in thetas: for _ in range(10): t, X = ornstein_uhlenbeck(X0, theta_test, mu, sigma, T, dt) axes[0, 1].plot(t, X, alpha=0.3, linewidth=1.5, label=f’θ={theta_test}’ if _ == 0 else None) axes[0, 1].axhline(y=mu, color=‘red’, linestyle=’—’, linewidth=2) axes[0, 1].set_xlabel(‘Time t’, fontsize=11) axes[0, 1].set_ylabel(‘X(t)’, fontsize=11) axes[0, 1].set_title(‘Effect of Reversion Speed θ’, fontsize=12, fontweight=‘bold’) axes[0, 1].legend() axes[0, 1].grid(alpha=0.3) # (3) Stationary distribution n_paths = 5000 final_values = [] for _ in range(n_paths): t, X = ornstein_uhlenbeck(X0, theta, mu, sigma, T, dt) final_values.append(X[-1]) final_values = np.array(final_values) axes[1, 0].hist(final_values, bins=50, density=True, alpha=0.6, color=‘#667eea’, edgecolor=‘black’, label=‘Simulation’) # Theoretical stationary distribution N(μ, σ²/(2θ)) x = np.linspace(final_values.min(), final_values.max(), 1000) stationary_var = sigma2 / (2 * theta) stationary_pdf = stats.norm.pdf(x, mu, np.sqrt(stationary_var)) axes[1, 0].plot(x, stationary_pdf, ‘r-’, linewidth=2.5, label=f’N({mu}, {stationary_var:.2f})’) axes[1, 0].set_xlabel(‘X(∞)’, fontsize=11) axes[1, 0].set_ylabel(‘Density’, fontsize=11) axes[1, 0].set_title(‘Stationary Distribution’, fontsize=12, fontweight=‘bold’) axes[1, 0].legend() axes[1, 0].grid(alpha=0.3) # (4) Time evolution of mean and variance all_paths = [] for _ in range(1000): t, X = ornstein_uhlenbeck(X0, theta, mu, sigma, T, dt) all_paths.append(X) all_paths = np.array(all_paths) mean_X = all_paths.mean(axis=0) var_X = all_paths.var(axis=0) # Theoretical mean and variance E_X_theory = X0 * np.exp(-theta * t) + mu * (1 - np.exp(-theta * t)) Var_X_theory = (sigma2 / (2theta)) * (1 - np.exp(-2theta*t)) axes[1, 1].plot(t, mean_X, color=‘#667eea’, linewidth=2.5, label=‘Observed Mean’) axes[1, 1].plot(t, E_X_theory, ‘r—’, linewidth=2, label=‘Theoretical Mean’) axes[1, 1].fill_between(t, mean_X - np.sqrt(var_X), mean_X + np.sqrt(var_X), alpha=0.2, color=‘#667eea’, label=‘Observed ±1σ’) axes[1, 1].fill_between(t, E_X_theory - np.sqrt(Var_X_theory), E_X_theory + np.sqrt(Var_X_theory), alpha=0.2, color=‘red’) axes[1, 1].set_xlabel(‘Time t’, fontsize=11) axes[1, 1].set_ylabel(‘E[X(t)]’, fontsize=11) axes[1, 1].set_title(‘Time Evolution of Mean and Variance’, fontsize=12, fontweight=‘bold’) axes[1, 1].legend() axes[1, 1].grid(alpha=0.3) plt.tight_layout() plt.show() print(“Verification of Ornstein-Uhlenbeck Process:”) print(f”Theoretical stationary distribution: N({mu}, {stationary_var:.4f})”) print(f”Observed stationary distribution: N({final_values.mean():.4f}, {final_values.var():.4f})“)

4.6 Simulation of Stochastic Integrals

💻 Code Example 6: Simulation of Stochastic Integrals

Simulation of Itô integral: I(t) = ∫₀ᵗ f(s) dW(s) np.random.seed(42) def stochastic_integral(f_func, T, dt): """Simulation of stochastic integral ∫₀ᵗ f(s) dW(s)""" n_steps = int(T / dt) t = np.linspace(0, T, n_steps + 1) # Wiener process increments dW = np.random.normal(0, np.sqrt(dt), n_steps) # Calculation of stochastic integral (Itô integral) I = np.zeros(n_steps + 1) for i in range(n_steps): I[i+1] = I[i] + f_func(t[i]) * dW[i] return t, I T = 5 dt = 0.01 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # (1) Case f(t) = 1: ∫₀ᵗ dW(s) = W(t) f_const = lambda t: 1 n_paths = 20 for _ in range(n_paths): t, I = stochastic_integral(f_const, T, dt) axes[0, 0].plot(t, I, alpha=0.4, linewidth=1.5, color=‘#667eea’) axes[0, 0].set_xlabel(‘Time t’, fontsize=11) axes[0, 0].set_ylabel(‘∫₀ᵗ dW(s)’, fontsize=11) axes[0, 0].set_title(‘Itô Integral: ∫₀ᵗ dW(s) = W(t)’, fontsize=12, fontweight=‘bold’) axes[0, 0].grid(alpha=0.3) # (2) Case f(t) = t: ∫₀ᵗ s dW(s) f_linear = lambda t: t for _ in range(n_paths): t, I = stochastic_integral(f_linear, T, dt) axes[0, 1].plot(t, I, alpha=0.4, linewidth=1.5, color=‘#764ba2’) axes[0, 1].set_xlabel(‘Time t’, fontsize=11) axes[0, 1].set_ylabel(‘∫₀ᵗ s dW(s)’, fontsize=11) axes[0, 1].set_title(‘Itô Integral: ∫₀ᵗ s dW(s)’, fontsize=12, fontweight=‘bold’) axes[0, 1].grid(alpha=0.3) # (3) Property of Itô integral: E[I(t)] = 0 n_sims = 1000 final_values_const = [] final_values_linear = [] for _ in range(n_sims): t, I_const = stochastic_integral(f_const, T, dt) t, I_linear = stochastic_integral(f_linear, T, dt) final_values_const.append(I_const[-1]) final_values_linear.append(I_linear[-1]) axes[1, 0].hist(final_values_const, bins=50, alpha=0.6, density=True, color=‘#667eea’, edgecolor=‘black’, label=‘∫₀ᵗ dW(s)’) axes[1, 0].hist(final_values_linear, bins=50, alpha=0.6, density=True, color=‘#764ba2’, edgecolor=‘black’, label=‘∫₀ᵗ s dW(s)’) axes[1, 0].axvline(0, color=‘red’, linestyle=’—’, linewidth=2, label=‘E[I]=0’) axes[1, 0].set_xlabel(‘I(T)’, fontsize=11) axes[1, 0].set_ylabel(‘Density’, fontsize=11) axes[1, 0].set_title(‘Distribution of Itô Integral’, fontsize=12, fontweight=‘bold’) axes[1, 0].legend() axes[1, 0].grid(alpha=0.3) # (4) Variance of Itô integral: Var(∫₀ᵗ f(s) dW(s)) = ∫₀ᵗ f(s)² ds # Case f(t) = t: Var = ∫₀ᵗ s² ds = t³/3 theoretical_var = T**3 / 3 print(“Properties of Itô Integral:”) print(f”∫₀ᵗ dW(s): E={np.mean(final_values_const):.6f}, Var={np.var(final_values_const):.4f} (Theoretical Var={T:.4f})”) print(f”∫₀ᵗ s dW(s): E={np.mean(final_values_linear):.6f}, Var={np.var(final_values_linear):.4f} (Theoretical Var={theoretical_var:.4f})”) # Quadratic variation of Itô integral # (∫₀ᵗ f(s) dW(s))² - ∫₀ᵗ f(s)² ds is a martingale squared_integrals = [] quadratic_variations = [] for _ in range(n_sims): t, I = stochastic_integral(f_linear, T, dt) squared_integrals.append(I[-1]**2) # ∫₀ᵗ s² ds = t³/3 quadratic_variations.append(I[-1]2 - T3/3) axes[1, 1].hist(quadratic_variations, bins=50, alpha=0.6, color=‘#667eea’, edgecolor=‘black’) axes[1, 1].axvline(np.mean(quadratic_variations), color=‘red’, linestyle=’—’, linewidth=2, label=f’Mean={np.mean(quadratic_variations):.4f}’) axes[1, 1].set_xlabel(‘I² - [I,I]’, fontsize=11) axes[1, 1].set_ylabel(‘Frequency’, fontsize=11) axes[1, 1].set_title(‘Quadratic Variation (Martingale Property)’, fontsize=12, fontweight=‘bold’) axes[1, 1].legend() axes[1, 1].grid(alpha=0.3) plt.tight_layout() plt.show()

4.7 Langevin Equation

💻 Code Example 7: Langevin Equation

Langevin equation (overdamped): dv = -γv dt + σ dW # Application to materials science: motion of Brownian particles np.random.seed(42) # Parameters gamma = 1.0 # Friction coefficient sigma = np.sqrt(2 * gamma) # Volatility (by fluctuation-dissipation theorem) T_sim = 20 dt = 0.01 # Velocity simulation (Ornstein-Uhlenbeck process) def simulate_velocity(v0, gamma, sigma, T, dt): """Velocity simulation by Langevin equation""" mu_func = lambda v, t: -gamma * v sigma_func = lambda v, t: sigma return euler_maruyama(mu_func, sigma_func, v0, T, dt) fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # (1) Time evolution of velocity v0 = 5 n_paths = 50 for _ in range(n_paths): t, v = simulate_velocity(v0, gamma, sigma, T_sim, dt) axes[0, 0].plot(t, v, alpha=0.4, linewidth=1.5, color=‘#667eea’) axes[0, 0].axhline(y=0, color=‘red’, linestyle=’—’, linewidth=2, label=‘Equilibrium Velocity=0’) axes[0, 0].set_xlabel(‘Time’, fontsize=11) axes[0, 0].set_ylabel(‘Velocity v(t)’, fontsize=11) axes[0, 0].set_title(‘Time Evolution of Velocity by Langevin Equation’, fontsize=12, fontweight=‘bold’) axes[0, 0].legend() axes[0, 0].grid(alpha=0.3) # (2) Stationary distribution of velocity (Maxwell-Boltzmann distribution) n_sims = 5000 final_velocities = [] for _ in range(n_sims): t, v = simulate_velocity(v0, gamma, sigma, T_sim, dt) final_velocities.append(v[-1]) final_velocities = np.array(final_velocities) axes[0, 1].hist(final_velocities, bins=50, density=True, alpha=0.6, color=‘#667eea’, edgecolor=‘black’, label=‘Simulation’) # Theoretical distribution (by fluctuation-dissipation theorem N(0, σ²/(2γ))) x = np.linspace(final_velocities.min(), final_velocities.max(), 1000) stationary_var = sigma2 / (2 * gamma) stationary_pdf = stats.norm.pdf(x, 0, np.sqrt(stationary_var)) axes[0, 1].plot(x, stationary_pdf, ‘r-’, linewidth=2.5, label=f’N(0, {stationary_var:.2f})’) axes[0, 1].set_xlabel(‘Velocity v’, fontsize=11) axes[0, 1].set_ylabel(‘Density’, fontsize=11) axes[0, 1].set_title(‘Stationary Distribution of Velocity (Maxwell-Boltzmann)’, fontsize=12, fontweight=‘bold’) axes[0, 1].legend() axes[0, 1].grid(alpha=0.3) # (3) Time evolution of position (integrating velocity) def simulate_position(v0, gamma, sigma, T, dt): """Position simulation (time integration of velocity)""" t, v = simulate_velocity(v0, gamma, sigma, T, dt) x = np.cumsum(v) * dt return t, x, v for _ in range(20): t, x, v = simulate_position(v0, gamma, sigma, T_sim, dt) axes[1, 0].plot(t, x, alpha=0.4, linewidth=1.5, color=‘#764ba2’) axes[1, 0].set_xlabel(‘Time’, fontsize=11) axes[1, 0].set_ylabel(‘Position x(t)’, fontsize=11) axes[1, 0].set_title(‘Position of Brownian Particle (Diffusion)’, fontsize=12, fontweight=‘bold’) axes[1, 0].grid(alpha=0.3) # (4) Mean Square Displacement (MSD) # Measurement of diffusion coefficient n_paths_msd = 500 all_positions = [] for _ in range(n_paths_msd): t, x, v = simulate_position(0, gamma, sigma, T_sim, dt) all_positions.append(x) all_positions = np.array(all_positions) msd = np.mean(all_positions2, axis=0) # Theoretical MSD (long-time limit): = 2Dt, D = σ²/(2γ²) = 1/γ D_theory = 1 / gamma msd_theory = 2 * D_theory * t axes[1, 1].plot(t, msd, color=‘#667eea’, linewidth=2.5, label=‘Observed MSD’) axes[1, 1].plot(t, msd_theory, ‘r—’, linewidth=2, label=f’Theoretical MSD=2Dt (D={D_theory:.2f})’) axes[1, 1].set_xlabel(‘Time’, fontsize=11) axes[1, 1].set_ylabel(”, fontsize=11) axes[1, 1].set_title(‘Mean Square Displacement (Measurement of Diffusion Coefficient)’, fontsize=12, fontweight=‘bold’) axes[1, 1].legend() axes[1, 1].grid(alpha=0.3) plt.tight_layout() plt.show() # Measurement of diffusion coefficient # Estimate diffusion coefficient from slope of MSD = 2Dt from scipy.stats import linregress slope, intercept, r_value, p_value, std_err = linregress(t[100:], msd[100:]) D_measured = slope / 2 print(“Diffusion by Langevin Equation:”) print(f”Theoretical diffusion coefficient: D = σ²/(2γ²) = {D_theory:.4f}”) print(f”Measured diffusion coefficient: D = {D_measured:.4f}”) print(f”Einstein relation: D = kT/γ (fluctuation-dissipation theorem)”)

💡 Note: The Langevin equation is used to model various phenomena in materials science, including motion of Brownian particles, diffusion of colloids, and dynamics of polymers. The fluctuation-dissipation theorem relates friction and random force intensity in thermal equilibrium states.

Exercises

📝 Exercise 1: Properties of Wiener Process
For a Wiener process W(t):

  1. Verify by simulation that W(2) - W(1) and W(4) - W(3) are independent
  2. Verify by simulation that Cov(W(2), W(5)) = min(2, 5) = 2
  3. Confirm by numerical differentiation that the path of W(t) is nowhere differentiable

📝 Exercise 2: Application of Geometric Brownian Motion
Consider a geometric Brownian motion with initial stock price S₀=100, μ=0.08, σ=0.25:

  1. Simulate the stock price distribution after 1 year and visualize with a histogram
  2. Calculate the price of a European call option with strike price K=110 using Monte Carlo method
  3. Compare with the Black-Scholes formula and evaluate the error

📝 Exercise 3: Analysis of Ornstein-Uhlenbeck Process
For the OU process dX = 3(2-X)dt + 1.5 dW:

  1. Run a simulation starting from X₀=5 and visualize convergence to the mean
  2. Compare the stationary distribution with theoretical values
  3. Investigate how convergence time changes when varying the reversion speed θ

← Chapter 3 Chapter 5 →

Disclaimer

← Previous Chapter Back to Series Index Next Chapter →