From 43481398d33683cdbe2ed1087b87a6ed7a5e9d63 Mon Sep 17 00:00:00 2001 From: Chihiro Watanabe Date: Mon, 18 May 2026 16:50:38 +0900 Subject: [PATCH 1/2] Update rng usage in monte_carlo.md Co-Authored-By: Claude Sonnet 4.6 --- lectures/monte_carlo.md | 39 ++++++++++++++++++++------------------- 1 file changed, 20 insertions(+), 19 deletions(-) diff --git a/lectures/monte_carlo.md b/lectures/monte_carlo.md index 9f66c2296..72a331af6 100644 --- a/lectures/monte_carlo.md +++ b/lectures/monte_carlo.md @@ -45,7 +45,8 @@ We will use the following imports. ```{code-cell} ipython3 import numpy as np import matplotlib.pyplot as plt -from numpy.random import randn + +rng = np.random.default_rng() ``` @@ -168,9 +169,9 @@ $$ S = 0.0 for i in range(n): - X_1 = np.exp(μ_1 + σ_1 * randn()) - X_2 = np.exp(μ_2 + σ_2 * randn()) - X_3 = np.exp(μ_3 + σ_3 * randn()) + X_1 = np.exp(μ_1 + σ_1 * rng.standard_normal()) + X_2 = np.exp(μ_2 + σ_2 * rng.standard_normal()) + X_3 = np.exp(μ_3 + σ_3 * rng.standard_normal()) S += (X_1 + X_2 + X_3)**p S / n ``` @@ -183,9 +184,9 @@ We can also construct a function that contains these operations: def compute_mean(n=1_000_000): S = 0.0 for i in range(n): - X_1 = np.exp(μ_1 + σ_1 * randn()) - X_2 = np.exp(μ_2 + σ_2 * randn()) - X_3 = np.exp(μ_3 + σ_3 * randn()) + X_1 = np.exp(μ_1 + σ_1 * rng.standard_normal()) + X_2 = np.exp(μ_2 + σ_2 * rng.standard_normal()) + X_3 = np.exp(μ_3 + σ_3 * rng.standard_normal()) S += (X_1 + X_2 + X_3)**p return (S / n) ``` @@ -210,9 +211,9 @@ To make it faster, let's implement a vectorized routine using NumPy. ```{code-cell} ipython3 def compute_mean_vectorized(n=1_000_000): - X_1 = np.exp(μ_1 + σ_1 * randn(n)) - X_2 = np.exp(μ_2 + σ_2 * randn(n)) - X_3 = np.exp(μ_3 + σ_3 * randn(n)) + X_1 = np.exp(μ_1 + σ_1 * rng.standard_normal(n)) + X_2 = np.exp(μ_2 + σ_2 * rng.standard_normal(n)) + X_3 = np.exp(μ_3 + σ_3 * rng.standard_normal(n)) S = (X_1 + X_2 + X_3)**p return S.mean() ``` @@ -399,7 +400,7 @@ M = 10_000_000 Here is our code ```{code-cell} ipython3 -S = np.exp(μ + σ * np.random.randn(M)) +S = np.exp(μ + σ * rng.standard_normal(M)) return_draws = np.maximum(S - K, 0) P = β**n * np.mean(return_draws) print(f"The Monte Carlo option price is approximately {P:3f}") @@ -520,8 +521,8 @@ def simulate_asset_price_path(μ=default_μ, S0=default_S0, h0=default_h0, n=def h = h0 for t in range(n): - s[t+1] = s[t] + μ + np.exp(h) * randn() - h = ρ * h + ν * randn() + s[t+1] = s[t] + μ + np.exp(h) * rng.standard_normal() + h = ρ * h + ν * rng.standard_normal() return np.exp(s) ``` @@ -583,8 +584,8 @@ def compute_call_price(β=default_β, h = h0 # Simulate forward in time for t in range(n): - s = s + μ + np.exp(h) * randn() - h = ρ * h + ν * randn() + s = s + μ + np.exp(h) * rng.standard_normal() + h = ρ * h + ν * rng.standard_normal() # And add the value max{S_n - K, 0} to current_sum current_sum += np.maximum(np.exp(s) - K, 0) @@ -629,7 +630,7 @@ def compute_call_price_vector(β=default_β, s = np.full(M, np.log(S0)) h = np.full(M, h0) for t in range(n): - Z = np.random.randn(2, M) + Z = rng.standard_normal((2, M)) s = s + μ + np.exp(h) * Z[0, :] h = ρ * h + ν * Z[1, :] expectation = np.mean(np.maximum(np.exp(s) - K, 0)) @@ -706,8 +707,8 @@ def compute_call_price_with_barrier(β=default_β, option_is_null = False # Simulate forward in time for t in range(n): - s = s + μ + np.exp(h) * randn() - h = ρ * h + ν * randn() + s = s + μ + np.exp(h) * rng.standard_normal() + h = ρ * h + ν * rng.standard_normal() if np.exp(s) > bp: payoff = 0 option_is_null = True @@ -744,7 +745,7 @@ def compute_call_price_with_barrier_vector(β=default_β, h = np.full(M, h0) option_is_null = np.full(M, False) for t in range(n): - Z = np.random.randn(2, M) + Z = rng.standard_normal((2, M)) s = s + μ + np.exp(h) * Z[0, :] h = ρ * h + ν * Z[1, :] # Mark all the options null where S_n > barrier price From d594d79d5c2798570bdfbe259ec0ea3797d30f29 Mon Sep 17 00:00:00 2001 From: Chihiro Watanabe Date: Mon, 25 May 2026 11:45:03 +0900 Subject: [PATCH 2/2] Pass rng as explicit argument in monte_carlo.md functions Co-Authored-By: Claude Sonnet 4.6 --- lectures/monte_carlo.md | 51 +++++++++++++++++++++++++++-------------- 1 file changed, 34 insertions(+), 17 deletions(-) diff --git a/lectures/monte_carlo.md b/lectures/monte_carlo.md index 72a331af6..a7d72eb8c 100644 --- a/lectures/monte_carlo.md +++ b/lectures/monte_carlo.md @@ -181,7 +181,9 @@ S / n We can also construct a function that contains these operations: ```{code-cell} ipython3 -def compute_mean(n=1_000_000): +def compute_mean(n=1_000_000, rng=None): + if rng is None: + rng = np.random.default_rng() S = 0.0 for i in range(n): X_1 = np.exp(μ_1 + σ_1 * rng.standard_normal()) @@ -196,7 +198,7 @@ def compute_mean(n=1_000_000): Now let's call it. ```{code-cell} ipython3 -compute_mean() +compute_mean(rng=rng) ``` @@ -210,7 +212,9 @@ But the code above runs quite slowly. To make it faster, let's implement a vectorized routine using NumPy. ```{code-cell} ipython3 -def compute_mean_vectorized(n=1_000_000): +def compute_mean_vectorized(n=1_000_000, rng=None): + if rng is None: + rng = np.random.default_rng() X_1 = np.exp(μ_1 + σ_1 * rng.standard_normal(n)) X_2 = np.exp(μ_2 + σ_2 * rng.standard_normal(n)) X_3 = np.exp(μ_3 + σ_3 * rng.standard_normal(n)) @@ -221,7 +225,7 @@ def compute_mean_vectorized(n=1_000_000): ```{code-cell} ipython3 %%time -compute_mean_vectorized() +compute_mean_vectorized(rng=rng) ``` @@ -233,7 +237,7 @@ We can increase $n$ to get more accuracy and still have reasonable speed: ```{code-cell} ipython3 %%time -compute_mean_vectorized(n=10_000_000) +compute_mean_vectorized(n=10_000_000, rng=rng) ``` @@ -515,7 +519,9 @@ $$ s_{t+1} = s_t + \mu + \exp(h_t) \xi_{t+1} $$ Here is a function to simulate a path using this equation: ```{code-cell} ipython3 -def simulate_asset_price_path(μ=default_μ, S0=default_S0, h0=default_h0, n=default_n, ρ=default_ρ, ν=default_ν): +def simulate_asset_price_path(μ=default_μ, S0=default_S0, h0=default_h0, n=default_n, ρ=default_ρ, ν=default_ν, rng=None): + if rng is None: + rng = np.random.default_rng() s = np.empty(n+1) s[0] = np.log(S0) @@ -538,7 +544,7 @@ titles = 'log paths', 'paths' transforms = np.log, lambda x: x for ax, transform, title in zip(axes, transforms, titles): for i in range(50): - path = simulate_asset_price_path() + path = simulate_asset_price_path(rng=rng) ax.plot(transform(path)) ax.set_title(title) @@ -576,7 +582,10 @@ def compute_call_price(β=default_β, n=default_n, ρ=default_ρ, ν=default_ν, - M=10_000): + M=10_000, + rng=None): + if rng is None: + rng = np.random.default_rng() current_sum = 0.0 # For each sample path for m in range(M): @@ -594,7 +603,7 @@ def compute_call_price(β=default_β, ```{code-cell} ipython3 %%time -compute_call_price() +compute_call_price(rng=rng) ``` @@ -625,8 +634,10 @@ def compute_call_price_vector(β=default_β, n=default_n, ρ=default_ρ, ν=default_ν, - M=10_000): - + M=10_000, + rng=None): + if rng is None: + rng = np.random.default_rng() s = np.full(M, np.log(S0)) h = np.full(M, h0) for t in range(n): @@ -640,7 +651,7 @@ def compute_call_price_vector(β=default_β, ```{code-cell} ipython3 %%time -compute_call_price_vector() +compute_call_price_vector(rng=rng) ``` @@ -651,7 +662,7 @@ Now let's try with larger $M$ to get a more accurate calculation. ```{code-cell} ipython3 %%time -compute_call_price(M=10_000_000) +compute_call_price(M=10_000_000, rng=rng) ``` @@ -697,7 +708,10 @@ def compute_call_price_with_barrier(β=default_β, ρ=default_ρ, ν=default_ν, bp=default_bp, - M=50_000): + M=50_000, + rng=None): + if rng is None: + rng = np.random.default_rng() current_sum = 0.0 # For each sample path for m in range(M): @@ -723,7 +737,7 @@ def compute_call_price_with_barrier(β=default_β, ``` ```{code-cell} ipython3 -%time compute_call_price_with_barrier() +%time compute_call_price_with_barrier(rng=rng) ``` @@ -740,7 +754,10 @@ def compute_call_price_with_barrier_vector(β=default_β, ρ=default_ρ, ν=default_ν, bp=default_bp, - M=50_000): + M=50_000, + rng=None): + if rng is None: + rng = np.random.default_rng() s = np.full(M, np.log(S0)) h = np.full(M, h0) option_is_null = np.full(M, False) @@ -758,7 +775,7 @@ def compute_call_price_with_barrier_vector(β=default_β, ``` ```{code-cell} ipython3 -%time compute_call_price_with_barrier_vector() +%time compute_call_price_with_barrier_vector(rng=rng) ``` ```{solution-end}