lectures/monte_carlo.md
            
                      Original file line number
                      Diff line number
                      Diff line change
                  
    @@ -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

    ```

    @@ -180,12 +181,14 @@ 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 * 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)

    ```

    @@ -195,7 +198,7 @@ def compute_mean(n=1_000_000):
  
    Now let's call it.

    ```{code-cell} ipython3

    compute_mean()

    compute_mean(rng=rng)

    ```

    @@ -209,18 +212,20 @@ 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):

        X_1 = np.exp(μ_1 + σ_1 * randn(n))

        X_2 = np.exp(μ_2 + σ_2 * randn(n))

        X_3 = np.exp(μ_3 + σ_3 * randn(n))

    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))

        S = (X_1 + X_2 + X_3)**p

        return S.mean()

    ```

    ```{code-cell} ipython3

    %%time

    compute_mean_vectorized()

    compute_mean_vectorized(rng=rng)

    ```

    @@ -232,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)

    ```

    @@ -399,7 +404,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}")

    @@ -514,14 +519,16 @@ $$ 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)

        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)

    ```

    @@ -537,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)

    @@ -575,16 +582,19 @@ 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):

            s = np.log(S0)

            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)

    @@ -593,7 +603,7 @@ def compute_call_price(β=default_β,
  
    ```{code-cell} ipython3

    %%time

    compute_call_price()

    compute_call_price(rng=rng)

    ```

    @@ -624,12 +634,14 @@ 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):

            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))

    @@ -639,7 +651,7 @@ def compute_call_price_vector(β=default_β,
  
    ```{code-cell} ipython3

    %%time

    compute_call_price_vector()

    compute_call_price_vector(rng=rng)

    ```

    @@ -650,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)

    ```

    @@ -696,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):

    @@ -706,8 +721,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

    @@ -722,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)

    ```

    @@ -739,12 +754,15 @@ 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)

        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

    @@ -757,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}
[monte_carlo.md] Update np.random → Generator API #741

New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Open
Chihiro2000GitHub wants to merge 3 commits into main from update-rng-monte-carlo
+54 −36
Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

[monte_carlo.md] Update np.random → Generator API #741

Diff view

Diff view

There are no files selected for viewing

Uh oh!

Uh oh!

[monte_carlo.md] Update np.random → Generator API #741

Are you sure you want to change the base?

[monte_carlo.md] Update np.random → Generator API #741

Uh oh!

Uh oh!

Diff view

Diff view

There are no files selected for viewing

Uh oh!
