diff --git a/lectures/parallelization.md b/lectures/parallelization.md index dd7cde0a..3e28e4e6 100644 --- a/lectures/parallelization.md +++ b/lectures/parallelization.md @@ -486,3 +486,119 @@ 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 + +n, β, K = 20, 0.99, 100 +μ, ρ, ν, S0, h0 = 0.0001, 0.1, 0.001, 10, 0 + +@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..9d99472f 100644 --- a/lectures/scipy.md +++ b/lectures/scipy.md @@ -440,12 +440,144 @@ We leave you to investigate the [set of available routines](http://docs.scipy.or ## Exercises +The first few exercises 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-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 + +$$ +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} +: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 +:class: dropdown +``` + +Here's one possible solution + +```{code-cell} ipython3 +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(μ)) + +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 +``` + +```{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 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`. ```