Skip to content

Commit 4348139

Browse files
Update rng usage in monte_carlo.md
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
1 parent 5ee005d commit 4348139

1 file changed

Lines changed: 20 additions & 19 deletions

File tree

‎lectures/monte_carlo.md‎

Lines changed: 20 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,8 @@ We will use the following imports.
4545
```{code-cell} ipython3
4646
import numpy as np
4747
import matplotlib.pyplot as plt
48-
from numpy.random import randn
48+
49+
rng = np.random.default_rng()
4950
```
5051

5152

@@ -168,9 +169,9 @@ $$
168169
169170
S = 0.0
170171
for i in range(n):
171-
X_1 = np.exp(μ_1 + σ_1 * randn())
172-
X_2 = np.exp(μ_2 + σ_2 * randn())
173-
X_3 = np.exp(μ_3 + σ_3 * randn())
172+
X_1 = np.exp(μ_1 + σ_1 * rng.standard_normal())
173+
X_2 = np.exp(μ_2 + σ_2 * rng.standard_normal())
174+
X_3 = np.exp(μ_3 + σ_3 * rng.standard_normal())
174175
S += (X_1 + X_2 + X_3)**p
175176
S / n
176177
```
@@ -183,9 +184,9 @@ We can also construct a function that contains these operations:
183184
def compute_mean(n=1_000_000):
184185
S = 0.0
185186
for i in range(n):
186-
X_1 = np.exp(μ_1 + σ_1 * randn())
187-
X_2 = np.exp(μ_2 + σ_2 * randn())
188-
X_3 = np.exp(μ_3 + σ_3 * randn())
187+
X_1 = np.exp(μ_1 + σ_1 * rng.standard_normal())
188+
X_2 = np.exp(μ_2 + σ_2 * rng.standard_normal())
189+
X_3 = np.exp(μ_3 + σ_3 * rng.standard_normal())
189190
S += (X_1 + X_2 + X_3)**p
190191
return (S / n)
191192
```
@@ -210,9 +211,9 @@ To make it faster, let's implement a vectorized routine using NumPy.
210211

211212
```{code-cell} ipython3
212213
def compute_mean_vectorized(n=1_000_000):
213-
X_1 = np.exp(μ_1 + σ_1 * randn(n))
214-
X_2 = np.exp(μ_2 + σ_2 * randn(n))
215-
X_3 = np.exp(μ_3 + σ_3 * randn(n))
214+
X_1 = np.exp(μ_1 + σ_1 * rng.standard_normal(n))
215+
X_2 = np.exp(μ_2 + σ_2 * rng.standard_normal(n))
216+
X_3 = np.exp(μ_3 + σ_3 * rng.standard_normal(n))
216217
S = (X_1 + X_2 + X_3)**p
217218
return S.mean()
218219
```
@@ -399,7 +400,7 @@ M = 10_000_000
399400
Here is our code
400401

401402
```{code-cell} ipython3
402-
S = np.exp(μ + σ * np.random.randn(M))
403+
S = np.exp(μ + σ * rng.standard_normal(M))
403404
return_draws = np.maximum(S - K, 0)
404405
P = β**n * np.mean(return_draws)
405406
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
520521
521522
h = h0
522523
for t in range(n):
523-
s[t+1] = s[t] + μ + np.exp(h) * randn()
524-
h = ρ * h + ν * randn()
524+
s[t+1] = s[t] + μ + np.exp(h) * rng.standard_normal()
525+
h = ρ * h + ν * rng.standard_normal()
525526
526527
return np.exp(s)
527528
```
@@ -583,8 +584,8 @@ def compute_call_price(β=default_β,
583584
h = h0
584585
# Simulate forward in time
585586
for t in range(n):
586-
s = s + μ + np.exp(h) * randn()
587-
h = ρ * h + ν * randn()
587+
s = s + μ + np.exp(h) * rng.standard_normal()
588+
h = ρ * h + ν * rng.standard_normal()
588589
# And add the value max{S_n - K, 0} to current_sum
589590
current_sum += np.maximum(np.exp(s) - K, 0)
590591
@@ -629,7 +630,7 @@ def compute_call_price_vector(β=default_β,
629630
s = np.full(M, np.log(S0))
630631
h = np.full(M, h0)
631632
for t in range(n):
632-
Z = np.random.randn(2, M)
633+
Z = rng.standard_normal((2, M))
633634
s = s + μ + np.exp(h) * Z[0, :]
634635
h = ρ * h + ν * Z[1, :]
635636
expectation = np.mean(np.maximum(np.exp(s) - K, 0))
@@ -706,8 +707,8 @@ def compute_call_price_with_barrier(β=default_β,
706707
option_is_null = False
707708
# Simulate forward in time
708709
for t in range(n):
709-
s = s + μ + np.exp(h) * randn()
710-
h = ρ * h + ν * randn()
710+
s = s + μ + np.exp(h) * rng.standard_normal()
711+
h = ρ * h + ν * rng.standard_normal()
711712
if np.exp(s) > bp:
712713
payoff = 0
713714
option_is_null = True
@@ -744,7 +745,7 @@ def compute_call_price_with_barrier_vector(β=default_β,
744745
h = np.full(M, h0)
745746
option_is_null = np.full(M, False)
746747
for t in range(n):
747-
Z = np.random.randn(2, M)
748+
Z = rng.standard_normal((2, M))
748749
s = s + μ + np.exp(h) * Z[0, :]
749750
h = ρ * h + ν * Z[1, :]
750751
# Mark all the options null where S_n > barrier price

0 commit comments

Comments
 (0)