Price a digital option by integration

A structured-products desk is selling a note that pays $1 million in three months if an index then stands above 4,400. The index is at 4,250 today. The payoff is all or nothing at the strike: a cash-or-nothing call, also called a digital.

The desk's pricing library has no formula for it yet. Price it by integrating against the model's density for the index at expiry, then check the price with a simulation.

In the pricing model the index level at expiry $T$ is

$$S_T = S_0\,e^{(r - \sigma^2/2)T + \sigma\sqrt{T}\,Z}$$

where $S_0$ is today's level, $r$ the continuous rate, $\sigma$ the volatility and $Z$ a standard normal number.

Your task

1. Write lognormal_density(x, S0, T, r, sigma): the density of $S_T$ at x, which is 0 for $x \le 0$. Raise ValueError unless S0, T and sigma are all positive.
2. Write digital_by_integral(S0, K, T, r, sigma, payout): the price today, payout times $e^{-rT}$ times $P(S_T > K)$. Find the probability by integrating the density with quad.
3. Write digital_by_simulation(S0, K, T, r, sigma, payout, n, seed): the pair (price, standard error) from n simulated levels. Draw the normals with np.random.default_rng(seed).standard_normal(n). Raise ValueError if n is below 1.

digital_by_integral must stay accurate for a weekly digital deep in the money, and for a 10-year one at high volatility. With $T = 0$ or $\sigma = 0$ there is no randomness: the price is the discounted payout if $S_0e^{rT} > K$, and 0 otherwise. Raise ValueError if S0 or K is not positive, or if T or sigma is negative.

The data

SPOT, STRIKE, EXPIRY, RATE, VOL and PAYOUT hold the note's terms, with rates and volatility as annual decimals and the expiry in years.

Hint: The density

Case study 5.4.6 shows that $\ln S_T$ is normal, with mean $m = \ln S_0 + (r - \frac{1}{2}\sigma^2)T$ and standard deviation $s = \sigma\sqrt{T}$. Lesson 5.4.4 gives the lognormal density as $\frac{\varphi(z)}{s\,x}$, with $z = \frac{\ln x - m}{s}$.

Hint: Keeping quad on target

Over $[K, \infty)$, quad can miss a narrow peak far from $K$ and return 0. Substitute $x = e^{m + sz}$ (lesson 2.6.1): then $dx = s\,x\,dz$ and the integrand becomes $\varphi(z)$. Then cut the infinite end, as case study 2.6.5 does, keeping $z$ within about 10 of the centre at both ends.

Hint: The simulation

Case study 5.4.6 simulates lognormal prices from standard normals. Case study 4.7.5 gives the standard error of a simulated share $\hat{p}$ as $\sqrt{\hat{p}(1 - \hat{p})/n}$; scale it by the discounted payout, as the price is. The simulated price should land within 2 or 3 standard errors of the integral.