Solve the heat equation on a grid

An options desk wants a grid pricer it can trust before it builds one for products that have no formula. The quant lead asks for the simplest one first. It should solve the heat equation with the explicit scheme, know its own stability limit exactly, and match Black–Scholes on plain calls and puts.

The step size must never be guessed. Too long a step and the scheme blows up; too short and the pricer is slow for nothing.

Your task

1. Write explicit_heat(u0, r, steps, left=None, right=None): the explicit scheme for $u_t = u_{xx}$ from the profile u0, whose first and last values are the ends, with $r = \frac{\Delta t}{\Delta x^2}$. left and right are each None, which holds that end at its starting value, or a function f(n) giving its value after step n. Return the profile after steps steps, and leave u0 unchanged.
2. Write stability_limit(cells): the largest $r$ at which no wave grows on a rod of $M$ = cells cells with fixed ends. One step multiplies the grid wave $\sin(j\pi x)$ by $1 - 4r\sin^2\!\left(\frac{j\pi}{2M}\right)$, for $j = 1, \ldots, M - 1$.
3. Write heat_price(S, K, T, rate, sigma, kind="call", cells=200, width=6.0, ratio=0.4), the price of a European "call" or "put". rate is the interest rate, called $r$ in case study 6.4.4.
   • Change variables as in case study 6.4.4, so that $V = Ke^{ax + b\tau}u$ with $u_s = u_{xx}$ and $s = \frac{1}{2}\sigma^2\tau$.
   • Use cells + 1 points in $x$, centred on today's $x_0 = \ln(S/K)$ and reaching width times $\sigma\sqrt{T}$ on each side. Take the fewest equal steps from $s = 0$ to $\frac{1}{2}\sigma^2T$ that keep $\frac{\Delta s}{\Delta x^2}$ at or below ratio.
   • At the two ends, set $V$ after each step to what the option is nearly worth there. Far below the strike a call is worth 0 and a put $Ke^{-r\tau} - S$. Far above it a call is worth $S - Ke^{-r\tau}$ and a put 0.
   • Read the price at the middle point. At $T = 0$, return the payoff.

explicit_heat raises ValueError if r is not positive, or if steps is not a whole number of at least 0. It raises it too when u0 has fewer than three points. stability_limit raises it unless cells is a whole number of at least 2.

heat_price raises it if S, K or sigma is not positive, T is negative, or kind is neither "call" nor "put". It also raises it if cells is not an even whole number, or if ratio is not positive or is above stability_limit(cells).

The data

There is no market data. The code provides bs_price(S, K, T, rate, sigma, kind), the Black–Scholes formula, to check your prices against. The tests use the six-month call and put of case study 6.4.4 and five other options.

Hint: Where to start

Lesson 6.5.1 writes the explicit update as one line of slices, u[1:-1] + r * (u[2:] - 2 * u[1:-1] + u[:-2]). It also finds the factor $1 - 4r\sin^2(\pi\Delta x/2)$ for the first sine wave. Ask which wave has the factor furthest below 0.

Hint: From the option to the rod

Case study 6.4.4 gives $a = \frac{1}{2} - \frac{r}{\sigma^2}$ and $b = -r - \frac{1}{2}\sigma^2a^2$, with $r$ the interest rate, and the call's starting profile $\max(e^{(1 - a)x} - e^{-ax}, 0)$. Lesson 6.4.2 shows that $e^{cx}$ grows into $e^{cx + c^2s}$, which checks your end values in $u$.

Hint: Counting the steps

Round the number of steps up, so that the actual ratio never exceeds ratio. Lesson 6.5.1's tent table shows what a ratio above the limit does.