This post explains how to use moment matching to reduce variance in Monte Carlo simulation of the Hull-White term structure model.

In an earlier blog post on how the Hull-White Monte Carlo simulations are notorious for not coverging with some of the expected moments. In this post, I would like to touch upon a variance reduction technique called *moment matching* that can be employed to fix this issue of convergence.

The idea behind moment matching is rather simple. Lets consider the specific example of short rate model. For the short rate model, it is known that the average of stochastic discount factors generated from each path has to agree with the model (or the give yield curve) discount factors. The idea of moment matching is to correct the short rates generated by the term structure model such that the average of stochastic discount factors from the simulation matches the model discount factors.

In [1]:

```
import QuantLib as ql
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.integrate import cumtrapz
ql.__version__
```

Out[1]:

'1.11'

In [2]:

```
sigma = 0.01
a = 0.001
timestep = 360
length = 30 # in years
forward_rate = 0.05
day_count = ql.Thirty360()
todays_date = ql.Date(15, 1, 2015)
```

In [3]:

```
ql.Settings.instance().evaluationDate = todays_date
yield_curve = ql.FlatForward(
todays_date,
ql.QuoteHandle(ql.SimpleQuote(forward_rate)),
day_count)
spot_curve_handle = ql.YieldTermStructureHandle(yield_curve)
```

`generate_paths`

function below is the time grid and a matrix of short rates generated by the model. This is discussed in detaio in the Hull-White simulation post.

In [4]:

```
hw_process = ql.HullWhiteProcess(spot_curve_handle, a, sigma)
rng = ql.GaussianRandomSequenceGenerator(
ql.UniformRandomSequenceGenerator(
timestep,
ql.UniformRandomGenerator(125)))
seq = ql.GaussianPathGenerator(hw_process, length, timestep, rng, False)
```

In [5]:

```
def generate_paths(num_paths, timestep):
arr = np.zeros((num_paths, timestep+1))
for i in range(num_paths):
sample_path = seq.next()
path = sample_path.value()
time = [path.time(j) for j in range(len(path))]
value = [path[j] for j in range(len(path))]
arr[i, :] = np.array(value)
return np.array(time), arr
```

Here is a plot of the generated short rates.

In [6]:

```
num_paths = 128
time, paths = generate_paths(num_paths, timestep)
for i in range(num_paths):
plt.plot(time, paths[i, :], lw=0.8, alpha=0.6)
plt.title("Hull-White Short Rate Simulation")
plt.show()
```

The model zero coupon bond price $B(0, T)$ is given as: $$B(0, T) = E\left[\exp\left(-\int_0^T r(t)dt \right) \right]$$ where $r(t)$ is the short rate generated by the model. The expectation of the stochastic discount factor at time $T$ is the price of the zero coupon bond at that time. In a simulation the paths are generated in a time grid and the discretization introduces some error. The empirical estimation of the zero coupon bond price from a Monte Carlo simulation $\hat{B}(0, t_m)$ maturing at time $t_m$ is given as:

$$\hat{B}(0, t_m) = \frac{1}{N}\sum_{i=1}^{N} \exp\left(-\sum_{j=0}^{m-1} \hat{r}_i(t_j)[t_{j+1}-t_j] \right)$$

where $\hat{r}_i(t_j)$ is the short rate for the path $i$ at time $t_j$ on the time grid. The expression for the moment matched short rates is given as [1]:

$$ r^c_i(t_j) = \hat{r}_i(t_j) + \frac{\log \hat{B}(0, t_{j+1}) - \log \hat{B}(0, t_{j})}{t_{j+1} - t_j} - \frac{\log B(0, t_{j+1}) - \log B(0, t_{j})}{t_{j+1} - t_j}$$

In [7]:

```
def stoch_df(paths, time):
return np.mean(
np.exp(-cumtrapz(paths, time, initial=0.)),axis=0
)
B_emp = stoch_df(paths, time)
logB_emp = np.log(B_emp)
B_yc = np.array([yield_curve.discount(t) for t in time])
logB_yc = np.log(B_yc)
```

In [8]:

```
deltaT = time[1:] - time[:-1]
deltaB_emp = logB_emp[1:]-logB_emp[:-1]
deltaB_yc = logB_yc[1:] - logB_yc[:-1]
```

In [9]:

```
new_paths = paths.copy()
new_paths[:,1:] += (deltaB_emp/deltaT - deltaB_yc/deltaT)
```

In [10]:

```
plt.plot(time,
stoch_df(paths, time),"r-.",
label="Original", lw=2)
plt.plot(time,
stoch_df(new_paths, time),"g:",
label="Corrected", lw=2)
plt.plot(time,B_yc, "k--",label="Market", lw=1)
plt.title("Zero Coupon Bond Price")
plt.legend()
```

Out[10]:

<matplotlib.legend.Legend at 0x7f78aa2f1990>

In [11]:

```
def alpha(forward, sigma, a, t):
return forward + 0.5* np.power(sigma/a*(1.0 - np.exp(-a*t)), 2)
avg = [np.mean(paths[:, i]) for i in range(timestep+1)]
new_avg = [np.mean(new_paths[:, i]) for i in range(timestep+1)]
plt.plot(time, avg, "r-.", lw=3, alpha=0.6, label="Original")
plt.plot(time, new_avg, "g:", lw=3, alpha=0.6, label="Corrected")
plt.plot(time,alpha(forward_rate, sigma, a, time), "k--", lw=2, alpha=0.6, label="Model")
plt.title("Mean of Short Rates")
plt.legend(loc=0)
```

Out[11]:

<matplotlib.legend.Legend at 0x7f78a8a3d2d0>

[1] Paul Glasserman, *Monte Carlo Methods in Financial Engineering*, Chapter 4

quantlib python finance

- Valuing Interest Rate Caps and Floors Using QuantLib Python
- Valuing Bonds with Credit Spreads in QuantLib Python
- Introduction to QuantLib Python
- American Option Pricing with QuantLib and Python
- QuantLib Python Notebooks On Docker