Hull–White Simulation with Euler–Maruyama Discretization

# Is this the correct discretisation of the Hull-White SDE for building a python model?


# Is this the correct discretisation of the Hull-White SDE for building a python model?












I've tried to build a basic one-factor Hull-White model using python, which I've done by trying to discretise the characteristic SDE.

According to my notes, the Hull-White SDE is

$$ dr_t = \alpha (\mu(t) - r_t) dt + \sigma dW_t $$

where $\alpha$ is a parameter determining the rate of mean-reversion, $\mu(t)$ is a function whose value is the mean-reversion level at time $t$, $\sigma$ is a parameter determining the volatility of the rates produced by the model, and $W_t$ is a standard Brownian motion process.

I've derived the function $\mu(t)$ to match the current observed yield curve by setting it equal to $$ \mu(t) = f_0(0, t) + \alpha^{-1} \frac{\partial}{\partial t}f_0(0, t) + \frac{\sigma^2}{2 \alpha^2} (1 - e^{-2 \alpha t}) $$ where $$ f_0(0, t) := - \frac{\partial}{\partial t} \ln B_0(0, t) $$ and $B_0(0, t)$ is an interpolation function (polynomial of degree 4) which closely replicates the observed prices of ZCBs.

To get my $\mu(t)$ I've not discretised anything. Instead, I've just fitted my polynomial, differentiated this polynomial manually, and plugged it into the above formula.

My yield curve looks like this:

and this gives me a $\mu(t)$ that looks like this:

I've arbitrarily set $\alpha = 0.1$ and $\sigma = 0.1$ because I don't know what these parameters would normally look like (please let me know if these are stupid choices).

The model I've tried to build uses 1 year time steps so I've tried to discretise the SDE as follows: $$ \Delta r_t = \alpha (\mu(t) - r_t) + \sigma X $$ for $t \in [0, 39]$ where $X$ is a random variate simulated from a $N(0, 1)$ distribution and $r_0$ is the current instantaneous spot rate.

However, when running this model I get the following result:

This is clearly very wrong!

I suspect that the issue is the discretisation of the HW model. Am I correct to be suspicious of this? If so, what is the correct discretisation of the HW SDE?

For reference, my code looks as follows:

```
import numpy as np
from math import exp, sqrt
from matplotlib import pyplot as plt

def hull_white_one_fac(r0, alpha, sigma, mu):
    dw = np.random.normal(size=40)
    dt = np.arange(0, 40)
    r = [r0]
    drt = [alpha * (mu[t] - r[-1]) * dt[t] + sigma * dw[t] for t in range(0, 40)]
    return r0 + np.cumsum(drt)

def hull_white_mu(yc, alpha, sigma):
    zcb = [(1+ir)**-t for ir, t in zip(yc, range(0, 40))]
    p = np.polyfit(range(0, 40), zcb, 4)
    #int_zcb = [p[0] * x**4 + p[1] * x**3 + p[2] * x**2 + p[3] * x + p[4] for x in range(0, 39)]
    int_fwd = [- (x * (x * (4 * p[0] * x + 3 * p[1]) + 2 * p[2]) + p[3])/(x * (x * (x * (p[0] * x + p[1]) + p[2]) + p[3]) + p[4]) for x in range(0, 40)]
    diff_int_fwd = [(-x * (4 * p[0] * x + 3 * p[1]) - x * (8 * p[0] * x + 3 * p[1]) - 2 * p[2])/(x * (x * (x * (p[0] * x + p[1]) + p[2]) + p[3]) + p[4]) - ((x * (x * (p[0] * x + p[1]) + p[2]) + x * (x * (p[0] * x + p[1]) + x * (2 * p[0] * x + p[1]) + p[2]) + p[3]) * (-x * (x * (4 * p[0] * x + 3 * p[1]) + 2 * p[2]) - p[3]))/(x * (x * (x * (p[0] * x + p[1]) + p[2]) + p[3]) + p[4])**2 for x in range(0, 40)]
    mu = [int_fwd_t + 1 / alpha * diff_int_fwd_t + sigma**2 / (2 * alpha**2) * (1 - exp(-2 * alpha * t)) for int_fwd_t, diff_int_fwd_t, t in zip(int_fwd, diff_int_fwd, range(0, 40))]
    return mu

if __name__ == '__main__':
    yc = [0.0119570319370659, 0.013932701452346, 0.0155402499730131, 0.0170571519360345, 0.0184980810473157,
          0.0198172195995803, 0.0209961254997071, 0.0220294317532332, 0.0229150185187919, 0.0236552422871417,
          0.0242647899419137, 0.0247564999995886, 0.0251425189780907, 0.0254435362824561, 0.0256775380876704,
          0.0258580239023412, 0.0259953222775315, 0.0260974691351983, 0.0261708094460174, 0.026220418849926,
          0.0262504051003476, 0.0262641274856628, 0.0262643591390461, 0.0262534088604249, 0.0262332156314007,
          0.0262054249416492, 0.0261714378113969, 0.0261324526863991, 0.0260895002981996, 0.0260434715696152,
          0.0259951401334477, 0.0259451806407383, 0.0258941837503939, 0.0258426684823296, 0.0257910924605755,
          0.0257398604560517, 0.0256892937533311, 0.0256395310703303, 0.0255906600981819, 0.0255427589299506,
          0.0254958972440893]
    a = 0.1
    s = 0.1

    mu = hull_white_mu(yc, a, s)

    hw = hull_white_one_fac(0.0119570319370659, a, s, mu)

    plt.plot(mu)
    plt.title("mu(t)")
    plt.show()

    plt.plot(yc)
    plt.title("Yield Curve")
    plt.show()

    plt.plot(hw)
    plt.title("Hull-White Realisation")
    plt.show()
```

## Answer by Nap D. Lover (score 3, accepted)

https://quant.stackexchange.com/a/74112

I am led to believe that your discretization is fine but uses too large of a step-size. This answer will tackle discretizing SDEs in general and then apply it to the Hull-White model.

The simplest and most well-known method for discretizing SDEs is the Euler Maruyama scheme. Given the one-dimensional SDE $$dX_t = \mu(t, X_t) dt + \sigma(t, X_t) dB_t$$ we can generate sample paths over the interval $[0,T]$ by discretizing this interval into $n+1$ points: $$0=t_0 < t_1<\dotsc < t_n = T.$$ Next let $Z_1,Z_2, \dotsc$ be an IID collection of standard Normal RVs $\mathcal{N}(0,1)$. Then iterate $$X_{t_i} = X_{t_{i-1}} +\mu(t_{i-1}, X_{t_{i-1}}) h+ \sqrt{h} \sigma(t_{i-1}, X_{t_{i-1}})Z_i,$$ where $h=(t_i-t_{i-1})/n$ is the time step-size between successive nodes, which we assume is constant. See the link for discussions of order of convergence and numerical error, etc.

Here is python code to do so for arbitrary $\mu, \sigma$ that are functions of time and the state:

```
import numpy as np
def euler_maruyama(x0, T, mu, sigma, n):
    """ Assumes x0, mu and sigma are all scalars.
    """
    x = np.zeros(n + 1)
    x[0] = x0
    h = T / n
    for i in range(n):
        x[i + 1] = x[i] + mu(i * h, x[i]) * h + sigma(i * h, x[i]) * np.sqrt(h) * np.random.normal()
    return x
```

For Hull-White, we have the specific SDE $$dr_t = \alpha(\mu(t)-r_t) dt + \sigma dW_t,$$ where $\sigma, \alpha$ are constants. Applying the Euler scheme gives, and using the shorthand $r_i = r_{t_i}$ $$r_i = r_{i-1}+\alpha(\mu_{i-1} -r_{i-1})h +\sqrt{h} \sigma Z_i$$

You can just as well implement this from scratch in python. But let's rather implement the drift and volatility and pass them to our Euler Maruyama function:

```
import numpy as np
import matplotlib.pyplot as plt
r0 = 0.9
alpha = 0.1
sigma0 = 0.1
T = 39
n = 100000

def mu(t, x):
    return alpha * (0.5 - r0)

def sigma(t, x):
    return sigma0

x = euler_maruyama(r0, T, mu, sigma, n)
t = np.linspace(0, T, n + 1)
fig = plt.figure()
plt.plot(t, x)
plt.show()
```

Here we assume $\mu(t)=\mu$, the long-term equilibrium rate, is constant. This produces a sample path like:

Notice that negative rates are allowed--this is baked into the design of the SDE. Other models, like Cox-Ingersoll-Ross avoid this by using a diffusion term $\sigma \sqrt{r_t} dW_t.$ We show one more simulation with a specific function $\mu(t,r)=\alpha(m(t)-r)$ where $m$ is given by $$m(t)=(t^2-(t-50)^3)0.0001:$$

All of these plots were produced with $T=39$ and $n=100000$ time steps, so the time step-size is $$h=T/n=0.00039.$$ Now the accuracy of the Euler-Maruyama scheme depends on having $h$ be small. But if you are curious of seeing a simulation with $h=1$, so $n=39$, then we obtain (using the same long term rate function as before):

Please comment for questions, clarifications or corrections.