From 01922a26ceac4821ec8c67a7367ed11fcf0cc28a Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Mon, 3 Oct 2022 06:29:00 +1100 Subject: [PATCH 01/11] Added Option Pricing Exercises --- lectures/parallelization.md | 107 ++++++++++++++++++++++++++++++++ lectures/scipy.md | 119 ++++++++++++++++++++++++++++++++++++ 2 files changed, 226 insertions(+) diff --git a/lectures/parallelization.md b/lectures/parallelization.md index dd7cde0a..5e3a29d7 100644 --- a/lectures/parallelization.md +++ b/lectures/parallelization.md @@ -486,3 +486,110 @@ on the number of CPUs on your machine.) ```{solution-end} ``` + + +```{exercise} +:label: parallel_ex2 + +In {doc}`our lecture on SciPy`, we discussed pricing a call option in a +setting where the underlying stock price had a simple and well-known +distribution. + +Here we discuss a more realistic setting. + +We recall that the price of the option obeys + +$$ P = \beta^n \mathbb E \max\{ S_n - K, 0 \} $$ + +where + +1. $\beta$ is a discount factor, +2. $n$ is the expiry date, +2. $K$ is the strike price and +3. $\{S_t\}$ is the price of the underlying asset at each time $t$. + +Suppose that `n, β, K = 20, 0.99, 100`. + +Assume that the stock price obeys + +$$ \ln \frac{S_{t+1}}{S_t} = \mu + \sigma_t \xi_{t+1} $$ + +where + +$$ + \sigma_t = \exp(h_t), + \quad + h_{t+1} = \rho h_t + \nu \eta_{t+1} +$$ + +Here $\{\xi_t\}$ and $\{\eta_t\}$ are IID and standard normal. + +(This is a **stochastic volatility** model, where the volatility $\sigma_t$ +varies over time.) + +Use the defaults `μ, ρ, ν, S0, h0 = 0.0001, 0.1, 0.001, 10, 0`. + +(Here `S0` is $S_0$ and `h0` is $h_0$.) + +By generating $M$ paths $s_0, \ldots, s_n$, compute the Monte Carlo estimate + +$$ + \hat P_M + := \beta^n \mathbb E \max\{ S_n - K, 0 \} + \approx + \frac{1}{M} \sum_{m=1}^M \max \{S_n^m - K, 0 \} + $$ + + +of the price, applying Numba and parallelization. + +``` + + +```{solution-start} parallel_ex2 +:class: dropdown +``` + + +With $s_t := \ln S_t$, the price dynamics become + +$$ s_{t+1} = s_t + \mu + \exp(h_t) \xi_{t+1} $$ + +Using this fact, the solution can be written as follows. + + +```{code-cell} ipython3 +from numpy.random import randn +M = 10_000_000 + +@njit(parallel=True) +def compute_call_price_parallel(β=β, + μ=μ, + S0=S0, + h0=h0, + K=K, + n=n, + ρ=ρ, + ν=ν, + M=M): + current_sum = 0.0 + # For each sample path + for m in prange(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() + # And add the value max{S_n - K, 0} to current_sum + current_sum += np.maximum(np.exp(s) - K, 0) + + return β**n * current_sum / M +``` + +Try swapping between `parallel=True` and `parallel=False` and noting the run time. + +If you are on a machine with many CPUs, the difference should be significant. + +```{solution-end} +``` diff --git a/lectures/scipy.md b/lectures/scipy.md index 6e165ef9..8bd24328 100644 --- a/lectures/scipy.md +++ b/lectures/scipy.md @@ -440,6 +440,125 @@ We leave you to investigate the [set of available routines](http://docs.scipy.or ## Exercises +The first few execises concern pricing a European call option under the +assumption of risk neutrality. The price satisfies + +$$ P = \beta^n \mathbb E \max\{ S_n - K, 0 \} $$ + +where + +1. $\beta$ is a discount factor, +2. $n$ is the expiry date, +2. $K$ is the strike price and +3. $\{S_t\}$ is the price of the underlying asset at each time $t$. + +For example, if the call option is to buy stock in Amazon at strike price $K$, the owner has the right (but not the obligation) to buy 1 share in Amazon at price $K$ after $n$ days. + +The payoff is therefore $\max\{S_n - K, 0\}$ + +The price is the expectation of the payoff, discounted to current value. + + +```{exercise} +:label: sp_ex01 + +Suppose that $S_n$ has the [log-normal](https://en.wikipedia.org/wiki/Log-normal_distribution) distribution with parameters $\mu$ and $\sigma$. Let $f$ denote the density of this distribution. Then + +$$ P = \beta^n \int_0^\infty \max\{x - K, 0\} f(x) dx $$ + +Plot the function + +$$g(x) = \beta^n \max\{x - K, 0\} f(x)$$ + +over the interval $[0, 400]$ when `μ, σ, β, n, K = 4, 0.25, 0.99, 10, 40`. + +Hint: From `scipy.stats` you can import `lognorm` and then use `lognorm(x, σ, scale=np.exp(μ)` to get the density $f$. +``` + +```{solution-start} sp_ex01 +:class: dropdown +``` + +Here's one possible solution + +```{code-cell} ipython3 +from scipy.integrate import quad +from scipy.stats import lognorm + +def g(x): + return β**n * np.maximum(x - K, 0) * lognorm.pdf(x, σ, scale=np.exp(μ)) + +x_grid = np.linspace(0, 400, 1000) +y_grid = g(x_grid) + +fig, ax = plt.subplots() +ax.plot(x_grid, y_grid, label="$g$") +ax.legend() +plt.show() +``` + +```{solution-end} +``` + +```{exercise} +:label: sp_ex02 + +In order to get the option price, compute the integral of this function numerically using `quad` from `scipy.optimize`. + +``` + +```{solution-start} sp_ex02 +:class: dropdown +``` + +The integral and hence the price is + +```{code-cell} ipython3 +P, error = quad(g, 0, 1_000) +print(f"The numerical integration based option price is {P:3f}") +``` + +```{solution-end} +``` + +```{exercise} +:label: sp_ex03 + +Try to get a similar result using Monte Carlo to compute the expectation term in the option price, rather than `quad`. + +In particular, use the fact that if $S_n^1, \ldots, S_n^M$ are independent +draws from the lognormal distribution specified above, then, by the law of +large numbers, + +$$ \mathbb E \max\{ S_n - K, 0 \} + \approx + \frac{1}{M} \sum_{m=1}^M \max \{S_n^m - K, 0 \} + $$ + +Set `M = 10_000_000` + +``` + +```{solution-start} sp_ex03 +:class: dropdown +``` + +Here is one solution: + +```{code-cell} ipython3 +M = 10_000_000 +S = np.exp(μ + σ * np.random.randn(M)) +return_draws = np.maximum(S - K, 0) +P = β**n * np.mean(return_draws) +print(f"The Monte Carlo option price is {P:3f}") +``` + + +```{solution-end} +``` + + + ```{exercise} :label: sp_ex1 From 63d5477599748c04b4800d3eb3ac744119259065 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Mon, 3 Oct 2022 08:22:42 +1100 Subject: [PATCH 02/11] fix default values --- lectures/parallelization.md | 5 ++++- lectures/scipy.md | 2 ++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/lectures/parallelization.md b/lectures/parallelization.md index 5e3a29d7..cb3e2428 100644 --- a/lectures/parallelization.md +++ b/lectures/parallelization.md @@ -538,7 +538,7 @@ $$ := \beta^n \mathbb E \max\{ S_n - K, 0 \} \approx \frac{1}{M} \sum_{m=1}^M \max \{S_n^m - K, 0 \} - $$ +$$ of the price, applying Numba and parallelization. @@ -562,6 +562,9 @@ Using this fact, the solution can be written as follows. from numpy.random import randn M = 10_000_000 +n, β, K = 20, 0.99, 100 +μ, ρ, ν, S0, h0 = 0.0001, 0.1, 0.001, 10, 0 + @njit(parallel=True) def compute_call_price_parallel(β=β, μ=μ, diff --git a/lectures/scipy.md b/lectures/scipy.md index 8bd24328..9a4e8375 100644 --- a/lectures/scipy.md +++ b/lectures/scipy.md @@ -485,6 +485,8 @@ Here's one possible solution from scipy.integrate import quad from scipy.stats import lognorm +μ, σ, β, n, K = 4, 0.25, 0.99, 10, 40 + def g(x): return β**n * np.maximum(x - K, 0) * lognorm.pdf(x, σ, scale=np.exp(μ)) From 7b25619899127e24479f87527298e7f1af68c861 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Mon, 3 Oct 2022 14:11:04 +1100 Subject: [PATCH 03/11] update hint admonitions --- lectures/scipy.md | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/lectures/scipy.md b/lectures/scipy.md index 9a4e8375..85594565 100644 --- a/lectures/scipy.md +++ b/lectures/scipy.md @@ -459,8 +459,9 @@ The payoff is therefore $\max\{S_n - K, 0\}$ The price is the expectation of the payoff, discounted to current value. -```{exercise} +```{exercise-start} :label: sp_ex01 +``` Suppose that $S_n$ has the [log-normal](https://en.wikipedia.org/wiki/Log-normal_distribution) distribution with parameters $\mu$ and $\sigma$. Let $f$ denote the density of this distribution. Then @@ -472,7 +473,13 @@ $$g(x) = \beta^n \max\{x - K, 0\} f(x)$$ over the interval $[0, 400]$ when `μ, σ, β, n, K = 4, 0.25, 0.99, 10, 40`. -Hint: From `scipy.stats` you can import `lognorm` and then use `lognorm(x, σ, scale=np.exp(μ)` to get the density $f$. +```{hint} +:class: dropdown + +From `scipy.stats` you can import `lognorm` and then use `lognorm(x, σ, scale=np.exp(μ)` to get the density $f$. +``` + +```{exercise-end} ``` ```{solution-start} sp_ex01 From 456b0f24d634fccef3d38a03308d883618d589f7 Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Mon, 3 Oct 2022 15:04:10 +1100 Subject: [PATCH 04/11] misc --- lectures/scipy.md | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/lectures/scipy.md b/lectures/scipy.md index 85594565..ddecf714 100644 --- a/lectures/scipy.md +++ b/lectures/scipy.md @@ -520,11 +520,9 @@ In order to get the option price, compute the integral of this function numerica :class: dropdown ``` -The integral and hence the price is - ```{code-cell} ipython3 P, error = quad(g, 0, 1_000) -print(f"The numerical integration based option price is {P:3f}") +print(f"The numerical integration based option price is {P:.3f}") ``` ```{solution-end} @@ -573,7 +571,7 @@ print(f"The Monte Carlo option price is {P:3f}") In {ref}`this lecture `, we discussed the concept of {ref}`recursive function calls `. -Try to write a recursive implementation of homemade bisection function {ref}`described above `. +Try to write a recursive implementation of the homemade bisection function {ref}`described above `. Test it on the function {eq}`root_f`. ``` From 95939e2ceb5c4904370da03edf281754ff9664f8 Mon Sep 17 00:00:00 2001 From: mmcky Date: Mon, 3 Oct 2022 15:29:29 +1100 Subject: [PATCH 05/11] Update lectures/parallelization.md --- lectures/parallelization.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lectures/parallelization.md b/lectures/parallelization.md index cb3e2428..6c8781be 100644 --- a/lectures/parallelization.md +++ b/lectures/parallelization.md @@ -553,7 +553,9 @@ of the price, applying Numba and parallelization. With $s_t := \ln S_t$, the price dynamics become -$$ s_{t+1} = s_t + \mu + \exp(h_t) \xi_{t+1} $$ +$$ +s_{t+1} = s_t + \mu + \exp(h_t) \xi_{t+1} +$$ Using this fact, the solution can be written as follows. From 05de98aee0167a8d89794d01da7fc844306792cf Mon Sep 17 00:00:00 2001 From: mmcky Date: Mon, 3 Oct 2022 15:29:36 +1100 Subject: [PATCH 06/11] Update lectures/scipy.md --- lectures/scipy.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lectures/scipy.md b/lectures/scipy.md index ddecf714..7be3f638 100644 --- a/lectures/scipy.md +++ b/lectures/scipy.md @@ -443,7 +443,9 @@ We leave you to investigate the [set of available routines](http://docs.scipy.or The first few execises concern pricing a European call option under the assumption of risk neutrality. The price satisfies -$$ P = \beta^n \mathbb E \max\{ S_n - K, 0 \} $$ +$$ +P = \beta^n \mathbb E \max\{ S_n - K, 0 \} +$$ where From 970cc2c13a8c020d0090e492a835841cb2253514 Mon Sep 17 00:00:00 2001 From: mmcky Date: Mon, 3 Oct 2022 15:29:44 +1100 Subject: [PATCH 07/11] Update lectures/scipy.md --- lectures/scipy.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lectures/scipy.md b/lectures/scipy.md index 7be3f638..d5a1bf4f 100644 --- a/lectures/scipy.md +++ b/lectures/scipy.md @@ -471,7 +471,9 @@ $$ P = \beta^n \int_0^\infty \max\{x - K, 0\} f(x) dx $$ Plot the function -$$g(x) = \beta^n \max\{x - K, 0\} f(x)$$ +$$ +g(x) = \beta^n \max\{x - K, 0\} f(x) +$$ over the interval $[0, 400]$ when `μ, σ, β, n, K = 4, 0.25, 0.99, 10, 40`. From 9e301b16c2475d02a7506ee9b86c1cb2955d2a93 Mon Sep 17 00:00:00 2001 From: mmcky Date: Mon, 3 Oct 2022 15:29:50 +1100 Subject: [PATCH 08/11] Update lectures/scipy.md --- lectures/scipy.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lectures/scipy.md b/lectures/scipy.md index d5a1bf4f..6c4474a8 100644 --- a/lectures/scipy.md +++ b/lectures/scipy.md @@ -467,7 +467,9 @@ The price is the expectation of the payoff, discounted to current value. Suppose that $S_n$ has the [log-normal](https://en.wikipedia.org/wiki/Log-normal_distribution) distribution with parameters $\mu$ and $\sigma$. Let $f$ denote the density of this distribution. Then -$$ P = \beta^n \int_0^\infty \max\{x - K, 0\} f(x) dx $$ +$$ +P = \beta^n \int_0^\infty \max\{x - K, 0\} f(x) dx +$$ Plot the function From af639872f234c8d062e10223644a4ebfac6f0ff1 Mon Sep 17 00:00:00 2001 From: mmcky Date: Mon, 3 Oct 2022 15:30:57 +1100 Subject: [PATCH 09/11] Update lectures/parallelization.md --- lectures/parallelization.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lectures/parallelization.md b/lectures/parallelization.md index 6c8781be..f842fe2d 100644 --- a/lectures/parallelization.md +++ b/lectures/parallelization.md @@ -512,7 +512,9 @@ Suppose that `n, β, K = 20, 0.99, 100`. Assume that the stock price obeys -$$ \ln \frac{S_{t+1}}{S_t} = \mu + \sigma_t \xi_{t+1} $$ +$$ +\ln \frac{S_{t+1}}{S_t} = \mu + \sigma_t \xi_{t+1} +$$ where From 1774cb339e5227c82367b571f81eb8e0c6dee098 Mon Sep 17 00:00:00 2001 From: mmcky Date: Mon, 3 Oct 2022 15:31:02 +1100 Subject: [PATCH 10/11] Update lectures/parallelization.md --- lectures/parallelization.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lectures/parallelization.md b/lectures/parallelization.md index f842fe2d..3e28e4e6 100644 --- a/lectures/parallelization.md +++ b/lectures/parallelization.md @@ -499,7 +499,9 @@ Here we discuss a more realistic setting. We recall that the price of the option obeys -$$ P = \beta^n \mathbb E \max\{ S_n - K, 0 \} $$ +$$ +P = \beta^n \mathbb E \max\{ S_n - K, 0 \} +$$ where From 2c01f9d541040690a3bcf2dbd301956aa39916d8 Mon Sep 17 00:00:00 2001 From: mmcky Date: Mon, 3 Oct 2022 17:27:51 +1100 Subject: [PATCH 11/11] Update lectures/scipy.md Fixed Typo --- lectures/scipy.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lectures/scipy.md b/lectures/scipy.md index 6c4474a8..9d99472f 100644 --- a/lectures/scipy.md +++ b/lectures/scipy.md @@ -440,7 +440,7 @@ We leave you to investigate the [set of available routines](http://docs.scipy.or ## Exercises -The first few execises concern pricing a European call option under the +The first few exercises concern pricing a European call option under the assumption of risk neutrality. The price satisfies $$