Added Option Pricing Exercises

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<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}
+```

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 <functions>`, we discussed the concept of {ref}`recursive function calls <recursive_functions>`.
 
-Try to write a recursive implementation of homemade bisection function {ref}`described above <bisect_func>`.
+Try to write a recursive implementation of the homemade bisection function {ref}`described above <bisect_func>`.
 
 Test it on the function {eq}`root_f`.
 ```