# Hybrid Monte Carlo

## Markov Models

In this notebook we analyse the evolution of HJM / Markov models. These models arise e.g. from Quasi Gaussian rates models or Markov Futures models.

We assume the model is driven by $d$-dimensional state variable $x$ and $d\times d$-dimensional auxilliary variable $y$.

The state variable $x$ follows a mean-reverting dynamic
$$
  dx(t) = \left[ \theta(y,t) - \chi(t) x(t) \right] dt + \sigma(t)^\top dW(t).
$$
Model parameters are $d\times d$-matrix of volatilities $\sigma(t)$ and diagonal $d\times d$-matrix of mean reversion speed parameters $\chi(t)$. The $d$-dimensional drift vector $\theta(y,t)$ couples state variables $x$ and auxilliary variables $y$.

The auxilliary variables $y$ follow the convection equation
$$
  dy(t) = \left[ \sigma(t)^\top \sigma(t) - \chi(t) y(t) - y(t) \chi(t) \right] dt.
$$

### Quasi-Gaussian Model

For Quasi-Gaussian rates models the drift term becomes
$$
  \theta(y,t) = y(t) \mathbb{1}.
$$

### Markov Futures Model

For Markov Futures models the drift term becomes
$$
  \theta(y,t) = \frac{1}{2} \left[ y(t)\chi(t) - \sigma(t)^\top \sigma(t) \right] \mathbb{1}.
$$

### Notation

We use the following functions to structure our notation:
$$
  H(s,t)   = \text{diag}\left\{ h_i(s,t)\right\},
  \quad \text{with} \quad
  h_i(s,t) = \exp\left\{ -\int_s^t \chi_i(u)du \right\}
$$
and
$$
  G(s,t)  = \int_s^t H(s,u)du = \int_s^t H(u,t)du.
$$
Functions $G$ and $H$ are related via $G'(s,t)=\partial G(s,t)/\partial t = H(s,t)$.

Moreover, we denote the instantanous variance as
$$
  V(t) = \sigma(t)^\top \sigma(t).
$$


## Evolution of State and Auxilliary Variables

In this section we derive a representation of the state and auxilliary variables. This representation is the basis for Monte-Carlo simulation.

We consider a time intervall from $s$ to $t$ ($s<t$). To simplify notation and derivation we assume model parameters are constant on that intervall. That is $\sigma(u)=\sigma$ and $\chi(u)=\chi$ for $u\in (s, t]$.

Note that even though model parameters are assumed (piece-wise) constant the drift term $\theta(u,t)$ is modelled exact. This is required because this quantity grows monotonically in time.




### Auxilliary Variables

The dynamics for $y$ are
$$
  dy(t) = \left[ V(t) - \chi(t) y(t) - y(t) \chi(t) \right] dt.
$$
The solution to this linear systems of ODEs is given by
$$
 y(t) = H(s,t)y(s)H(s,t) +
        \underbrace{\int_s^t H(u,t) V(u) H(u,t)du}_{M(s,t)}.
$$

For constant parameters the integral $M(s,t) = \left[ M_{i,j}(s,t) \right]$ is given as
$$
  M_{i,j}(s,t) = V_{i,j} \frac{1 - e^{-(\chi_i + \chi_j)(t-s)}}{\chi_i + \chi_j}.
$$
The elements of the first term $H(s,t)y(s)H(s,t)$ are $y_{i,j}(s)e^{-(\chi_i + \chi_j)(t-s)}$. This yields a representation of the elements $y_{i,j}(t)$ which we will use later on.

We get the representation of auxilliary variable elements as
\begin{align}
  y_{i,j}(t) &= y_{i,j}(s)e^{-(\chi_i + \chi_j)(t-s)} + 
                V_{i,j} \frac{1 - e^{-(\chi_i + \chi_j)(t-s)}}{\chi_i + \chi_j} \\
             &= \left[ y_{i,j}(s) - \frac{V_{i,j}}{\chi_i + \chi_j} \right] 
                 e^{-(\chi_i + \chi_j)(t-s)} +
                 \frac{V_{i,j}}{\chi_i + \chi_j}. \\
             &= a_{i,j} e^{-(\chi_i + \chi_j)(t-s)} + b_{i,j}
\end{align}
with $b_{i,j} = V_{i,j} / (\chi_i + \chi_j)$ and $a_{i,j} = V_{i,j} - b_{i,j}$.

Substituting $B = \left[ b_{i,j} \right]_{i,j}$ we get the matrix representation $A = y(s)-B$ and
$$
  y(t) = H(s,t) A H(s,t) + B.
$$

### State Variable

We have
$$
  dx(t) = \left[ \theta(y,t) - \chi(t) x(t) \right] dt + \sigma(t)^\top dW(t).
$$

The solution is given as
$$
 x(t) = H(s,t) \left[ x(s) + \int_s^t H(s,u)^{-1}
        \left[ \theta(y,u) du + \sigma(u)^\top dW(u) \right] \right].
$$

Covariance conditional on time $s$ becomes
$$
  \text{Cov}_s \left[x(t)\right] = \int_s^t H(u,t) \sigma(u)^\top \sigma(u) H(u,t) du.
$$
It turns out that
$$
  \text{Cov}_s \left[x(t)\right] = \int_s^t H(u,t) V(u) H(u,t) du = M(s,t).
$$

Conditional expectation is given as
$$
  \mathbb{E}_s\left[x(t)\right] = H(s,t) x(s) + \int_s^t H(u,t) \theta(y,u) du.
$$

For general drift functions $\theta(y,t)$ and with the representation $y(t) = H(s,t) A H(s,t) + B$ we can formulate $f(u) =  H(u,t) \theta\left(y(u),u\right)$. Then we can apply standard quadrature methods to calculate $\int_s^t f(u) du$. This typically works well since functions involved are smooth.

### Integrated Auxilliary Variable

The conditional expectation of the state variable requires calculation of $\int_s^t H(u,t) \theta(y,u) du$. In relevant applications the drift function $\theta(y,u)$ is an affine function of $y(u)$. In such situations we can calculate the integral.

As an intermediate step we calculate $\int_s^t H(u,t) y(u) du$. Using $y(u)=H(s,u) A H(s,u) + B$ yields
\begin{align}
  \int_s^t H(u,t) y(u) du
  &= \int_s^t H(u,t) \left[ H(s,u) A H(s,u) + B \right] du \\
  &= \int_s^t H(s,t) A H(s,u) du + \int_s^t H(u,t)B du \\
  &= H(s,t) A \int_s^t H(s,u) du + \int_s^t H(u,t) du B \\
  &= H(s,t) A G(s,t) + G(s,t) B.
\end{align}

We can also substitute back $A = y(s)-B$. This yields
\begin{align}
  \int_s^t H(u,t) y(u) du
  &= H(s,t) \left[ y(s) - B \right] G(s,t) + G(s,t) B \\
  &= H(s,t) y(s) G(s,t) - H(s,t) B G(s,t) + G(s,t) B \\
  &= H(s,t) y(s) G(s,t) + \left[G(s,t) B - H(s,t) B G(s,t) \right].
\end{align}

We have a closer look at the term $C = \left[G(s,t) B - H(s,t) B G(s,t) \right]$. The elements $C_{i,j}$ are
\begin{align*}
  C_{i,j} &= \frac{1-e^{-\chi_i(t-s)}}{\chi_i} \frac{V_{i,j}}{\chi_i + \chi_j} - 
             e^{-\chi_i(t-s)} \frac{V_{i,j}}{\chi_i + \chi_j} \frac{1-e^{-\chi_j(t-s)}}{\chi_j} \\
          &= \frac{V_{i,j}}{\chi_i + \chi_j} \left[ 
               \frac{1-e^{-\chi_i(t-s)}}{\chi_i} - e^{-\chi_i(t-s)} \frac{1-e^{-\chi_j(t-s)}}{\chi_j}
               \right].
\end{align*}
For the diagonal elements $C_{i,i}$ we further get
\begin{align*}
  C_{i,i} &= \frac{V_{i,j}}{2 \chi_i} \left[ 
               \frac{1-e^{-\chi_i(t-s)}}{\chi_i} - e^{-\chi_i(t-s)} \frac{1-e^{-\chi_i(t-s)}}{\chi_i} \right] \\
          &= \frac{V_{i,j}}{2 \chi_i} \frac{1-e^{-\chi_i(t-s)}}{\chi_i} \left[ 
               1 - e^{-\chi_i(t-s)} \right] \\
          &= \frac{V_{i,j}}{2} \left[ \frac{1-e^{-\chi_i(t-s)}}{\chi_i} \right]^2.
\end{align*}


### Quasi Gaussian Model Drift and Expectation

The drift term in Quasi Gaussian rates models is just $\theta(y,t) = y(t) \mathbb{1}$. This yields the conditional expectation
$$
  \mathbb{E}_s\left[x(t)\right] = H(s,t) \left[ x(s) + A G(s,t) \mathbb{1} \right] + G(s,t) B  \mathbb{1}
$$
or equivalently
$$
  \mathbb{E}_s\left[x(t)\right] = H(s,t) \left[ x(s) + y(s) G(s,t) \mathbb{1} \right] + \left[G(s,t) B - H(s,t)  B G(s,t) \right] \mathbb{1}.
$$


### Futures Model Drift and Expectation

For Futures Model the drift is
$\theta(y,t) = \frac{1}{2} \left[ y(t)\chi(t) - \sigma(t)^\top \sigma(t) \right] \mathbb{1}$.
We can use linearity of the integral operator and get
\begin{align}
  \int_s^t H(u,t) \theta(y,u) du
  &= \int_s^t H(u,t) \frac{1}{2} \left[ y(t)\chi(t) - \sigma(t)^\top \sigma(t) \right] \mathbb{1} du \\
  &= \frac{1}{2} \int_s^t H(u,t) \left[ y(t)\chi - V \right] \mathbb{1} du \\
  &= \frac{1}{2} \left[ \int_s^t H(u,t) y(t) du  \chi\mathbb{1} -  \int_s^t H(u,t) du V  \mathbb{1} \right] \\
  &= \frac{1}{2} \left\{ \left[ H(s,t) A G(s,t) + G(s,t) B \right] \chi\mathbb{1} -  G(s,t) V  \mathbb{1} \right\}.
\end{align}

This gives the conditional expectation of the state variable as
$$
  \mathbb{E}_s\left[x(t)\right] = H(s,t) x(s) +
  \frac{1}{2} \left\{ \left[ H(s,t) A G(s,t) + G(s,t) B\right] \chi -  G(s,t) V \right\}\mathbb{1}.
$$



## Futures Reconstruction

Using the Markov Model representation the futures price $F(t,T)$ can be reconstructed from state variables $x(t)$ and auxilliary variables $y(t)$ via
$$
  F(t,T) = F(0,T) * e^{X(t,T)}
$$
with
$$
  X(t,T) = \mathbb{1}^\top H(t,T) \left[ x(t) +
      \frac{1}{2} y(t) \left(I - H(t,T)\right)\mathbb{1}
      \right].
$$

## Model Test


```python
import sys
sys.path.append('../')  # make python find our modules
import numpy as np

from hybmc.models.MarkovFutureModel import MarkovFutureModel
```

We set up a simple 2-factor model.


```python
d = 2
times = np.array([0.0])
sigmaT = np.zeros([1,2,2])
sigmaT[0,:,:] = np.array([[ 0.10, 0.15 ],
                          [ 0.20, 0.25 ]])
chi = np.array([0.10, 0.20])
#
model = MarkovFutureModel(None,d,times,sigmaT,chi) # future curve is not yet implemented
```

Using the model we simulate the state variables.


```python
from hybmc.simulations.McSimulation import McSimulation
times = np.linspace(0.0, 5.0, 2)
nPaths = 2**13
seed = 314159265359
sim = McSimulation(model,times,nPaths,seed,False,showProgress=True)
```

Now we can calculate futures prices and analyse the terminal distrivution.


```python
idx = -1
t = times[idx]
dT = 1.0
T = t + dT
F = np.array([ model.futurePrice(t,T,X,None) for X in sim.X[:,idx,:] ])
Fav = np.mean(F)
sigma = np.std(F) / np.sqrt(t)
print('t: %6.2f, T: %6.2f, F: %6.4f, sigma: %6.4f' % (t,T,Fav,sigma) )
```

## Andersen Parametrisation

We test another parametrisation taken from Andersen 2008, sec. 7.6.


```python
def sigmaT(sigma_0, sigma_infty, rho_infty):
    h1 = -sigma_infty + rho_infty * sigma_0
    h2 = sigma_0 * np.sqrt(1.0 - rho_infty**2)
    hi = sigma_infty
    return np.array([ [h1, h2], [hi, 0.0] ])

def chi(kappa):
    return np.array([ kappa, 0.0 ])
```


```python
# Andersen example
kappa       = 1.35
sigma_0     = 0.50
sigma_infty = 0.17
rho_infty   = 0.50

sigmaT_ = sigmaT(sigma_0,sigma_infty,rho_infty)
chi_    = chi(kappa)
print(sigmaT_)
print(chi_)
```


```python
d = 2
times = np.array([0.0])
model = MarkovFutureModel(None,d,times,np.array([ sigmaT_ ]),chi_)
```


```python
times = np.linspace(0.0, 5.0, 2)
nPaths = 2**13
seed = 314159265359
sim = McSimulation(model,times,nPaths,seed,False,showProgress=True)
```


```python
idx = -1
t = times[idx]
dT = 1.0
T = t + dT
F = np.array([ model.futurePrice(t,T,X,None) for X in sim.X[:,idx,:] ])
Fav = np.mean(F)
sigma = np.std(F) / np.sqrt(t)
print('t: %6.2f, T: %6.2f, F: %6.4f, sigma: %6.4f' % (t,T,Fav,sigma) )
```


```python

```

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