# Hybrid Monte Carlo

## Affine Short Rate Models

In this notebook we analyse yield curve modelling based on affine term structure models. We start with a classical CIR model. Then we analyse initial yield curve calibration via deterministic shift extension. Finally, we also analyse the impact of square root processes on volatility smile.


```python
import sys
sys.path.append('../')  # make python find our modules
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go
import QuantLib as ql
```

## CIR Model Properties

As a first step we wet up a CIR model and analyse modelled yield curves and volatilities.


```python
from hybmc.models.AffineShortRateModel import CoxIngersollRossModel, quadraticExponential, cirMoments

r0         = 0.02
chi_       = 0.07
theta_     = 0.05
sigma_     = 0.07
cirModel = CoxIngersollRossModel(r0,chi_,theta_,sigma_,quadraticExponential(1.5))
```

We have a look at the *intitial* yield curve implied by the model.


```python
dT = 1.0/365.0
f = lambda t, T, rt : np.log(cirModel.zeroBondPrice(t,T,rt) / cirModel.zeroBondPrice(t,T+dT,rt)) / dT
```


```python
T = np.linspace(0.0,20.0,21)

X0 = cirModel.initialValues()
f_ = np.array([ f(0.0,T_,cirModel.r0) for T_ in T ])    
curve = pd.DataFrame([ T, f_ ]).T
curve.columns = ['T', 'f(0,T)']
fig = px.line(curve, x='T', y='f(0,T)')
fig.show()
```

Next we check *future* model-implied curves.


```python
T0 = 5.0
shortRates = np.linspace(0.02, 0.04, 5)
for r in reversed(shortRates):
    f_ = np.array([ f(T0,T0+T_,r) for T_ in T ])
    fig.add_trace(go.Scatter(x=T0+T, y=f_, mode='lines', name='r=%6.4f'%r))
fig.show()
```

The humped shape looks much better compared to Hull-White model.

We are also interested in the volatility of rates.

Zero bonds are given by $P(t,T,r) = \exp(-B_{CIR}(t,T) r(t) + A_{CIR})(t,T)$. Future zero rates from $T_0$ to $T_1$ are defined also

\begin{align}
 F(t;T_0,T_0) &= \frac{1}{T_1-T_0} \log\left( \frac{P\left(t,T_0,r(t)\right)}{P\left(t,T_1,r(t)\right)} \right) \\
              &= \frac{-\left[ B_{CIR}(t,T_0) - B_{CIR}(t,T_1)\right] r(t) + A_{CIR}(t,T_0) - A_{CIR}(t,T_1) }{T_1-T_0}.
\end{align}

In particular, we get for the variance of $F(T_0,T_0,T_1)$

$$
 Var\left[ F(T_0,T_0,T_1) \right] = \left[ \frac{B_{CIR}(T_0,T_1)}{T_1-T_0} \right]^2 \cdot  Var\left[ r(T_0) \right].
$$

This yields the proxy ATM swap rate volatility
$$
  \sigma(T_0,T_1) = \underbrace{\frac{B_{CIR}(T_0,T_1)}{T_1-T_0}}_{\lambda(T_0,T_1)} 
                    \underbrace{\sqrt{ \frac{Var\left[ r(T_0) \right]}{T_0} } }_{\sigma_{CIR}}
$$


```python
lambda_ = lambda T0,T1 : cirModel.ricattiAB(T0,T1,0.0,1.0)[1] / (T1 - T0)
expiryTimes = np.linspace(1.0, 10.0,10)
swapTerms = np.linspace(1.0, 10.0,10)
scalings = pd.DataFrame([ [T0, dT, lambda_(T0,T0+dT)] for T0 in expiryTimes for dT in swapTerms ],columns=['T0', 'dT', 'scaling'])
#fig = go.Figure(data=[go.Surface(x=scalings.T0,y=scalings.dT,z=scalings.scaling)])
fig = px.scatter_3d(scalings, x='T0', y='dT', z='scaling')
fig.show()
```


```python
sigma_CIR = lambda T0 : np.sqrt(cirMoments(cirModel.r0,T0,cirModel.chi(0.0),cirModel.theta(0.0),cirModel.sigma(0.0))[1] / T0)
vols = pd.DataFrame([ [T0, sigma_CIR(T0)] for T0 in expiryTimes], columns=['T0','sigma_CIR'])
fig = px.line(vols,x='T0',y='sigma_CIR')
fig.show()
```

## Yield Curve Fit


```python
import QuantLib as ql
today = ql.Settings.instance().evaluationDate
curveData = pd.DataFrame([[ 0.0,   5.0,   10.0,   20.0   ],
                          [ 0.020, 0.028,  0.033,  0.035 ]]).T
curveData.columns = ['T', 'f']
curveData['Date'] = [ today + int(t*365) for t in curveData['T'] ]
yts = ql.ForwardCurve(curveData['Date'],curveData['f'],ql.Actual365Fixed())
```


```python
fMarket = lambda T : yts.forwardRate(T,T,ql.Continuous).rate()
curve['fM(0,T)'] = [ fMarket(T) for T in curve['T'] ]
fig = go.Figure()
fig.add_trace(go.Scatter(x=curve['T'], y=curve['f(0,T)'], mode='lines', name='f(0,T)'))
fig.add_trace(go.Scatter(x=curve['T'], y=curve['fM(0,T)'], mode='lines', name='fM(0,T)'))
fig.show()
```


```python
zeroRateCir = lambda T : -np.log(cirModel.zeroBondPrice(0.0,T,cirModel.r0))/T
zeroRateYts = lambda T : -np.log(yts.discount(T)) / T
zeros = pd.DataFrame(np.linspace(0.1,20,200),columns=['T'])
zeros['CIR'] = [ zeroRateCir(T) for T in zeros['T'] ]
zeros['Yts'] = [ zeroRateYts(T) for T in zeros['T'] ]
fig = go.Figure()
fig.add_trace(go.Scatter(x=zeros['T'], y=zeros['CIR'], mode='lines', name='CIR'))
fig.add_trace(go.Scatter(x=zeros['T'], y=zeros['Yts'], mode='lines', name='Yts'))
fig.show()
```


```python
from hybmc.models.ShiftedRatesModel import ShiftedRatesModel
shiModel = ShiftedRatesModel(yts,cirModel)
```


```python
from hybmc.simulations.McSimulation import McSimulation

times = np.linspace(0.0,20.0,21)
nPaths = 2**13
seed = 3141
simCir = McSimulation(cirModel,times,nPaths,seed,showProgress=True)
simShi = McSimulation(shiModel,times,nPaths,seed,showProgress=True)
#
dT = 0.0
zcbCir = np.mean(np.array([
        [ cirModel.zeroBond(times[t],times[t]+dT,simCir.X[p,t,:],None) / cirModel.numeraire(times[t],simCir.X[p,t,:]) for t in range(len(times)) ]
        for p in range(nPaths) ]), axis=0)
zcbShi = np.mean(np.array([
        [ shiModel.zeroBond(times[t],times[t]+dT,simShi.X[p,t,:],None) / shiModel.numeraire(times[t],simShi.X[p,t,:]) for t in range(len(times)) ]
        for p in range(nPaths) ]), axis=0)
#
mcZeroCir = [ -np.log(df)/T for df,T in zip(zcbCir,times) ]
mcZeroShi = [ -np.log(df)/T for df,T in zip(zcbShi,times) ]
fig.add_trace(go.Scatter(x=times[1:], y=mcZeroCir[1:], mode='markers', name='CIR'))
fig.add_trace(go.Scatter(x=times[1:], y=mcZeroShi[1:], mode='markers', name='Yts'))
fig.show()
```


```python

```

07a176f32c974dac0befc96ffd0ab389cf56c437
ipynb
Jupyter Notebook
doc/AffineShortRateModels.ipynb
sschlenkrich/HybridMonteCarlo
72f54aa4bcd742430462b27b72d70369c01f9ac4
MIT
2021-08-18T18:34:41.000Z
2021-12-24T07:05:19.000Z
doc/AffineShortRateModels.ipynb
sschlenkrich/HybridMonteCarlo
72f54aa4bcd742430462b27b72d70369c01f9ac4
MIT
doc/AffineShortRateModels.ipynb
sschlenkrich/HybridMonteCarlo
72f54aa4bcd742430462b27b72d70369c01f9ac4
MIT
2021-01-31T11:41:19.000Z
2022-03-25T19:51:20.000Z
Qwen/Qwen-72B
1. YES
2. YES