# Calibration of Heston version of CIR

I'd like to calibrate a variant of Heston model for interest rates which is describe by this couple of SDE

\begin{aligned}dr_t&=a(b-r_t)+\sqrt{r_t}\sigma_t dW_t^1 \\ d\sigma_t&=k(\theta-\sigma_t)+\sqrt{\sigma_t}\zeta dW_t^2 \end{aligned}

with real markets data.

First of all I discretize the equations getting

\begin{aligned}r_t-r_{t-1}&=a(b-r_t)+\sqrt{r_{t-1}}\sigma_t \varepsilon_t \\ \sigma_t-\sigma_{t-1}&=k(\theta-\sigma_t)+\sqrt{\sigma_{t-1}}\zeta \epsilon_t \end{aligned}

And now I don't exactly how to proceed. I have some ideas:

• Try to apply a VAR (Vector AutoRegressive) model: rearraing the discretized equations I should get

\begin{aligned} \frac{r_t}{\sqrt{r_{t-1}}}&=ab+(a-1)\sqrt{r_{t-1}}+\sigma_t \varepsilon_t \\ \frac{\sigma_t}{\sqrt{\sigma_{t-1}}}&=k\theta+(k-1)\sqrt{\sigma_{t-1}}+\zeta \epsilon_t \end{aligned} or \begin{aligned} \frac{r_t-r_{t-1}}{\sqrt{r_{t-1}}}&=\frac{a(b-r_t)}{\sqrt{r_{t-1}}}+\sigma_t \varepsilon_t \\ \frac{\sigma_t-\sigma_{t-1}}{\sqrt{\sigma_{t-1}}}&=\frac{k(\theta-\sigma_t)}{\sqrt{\sigma_{t-1}}}+\zeta\epsilon_t \end{aligned} but I don't know whether this choice is suitable and which of these two proposals is the best one because of the square root and because the left sides of the equations are a trasformation of the right ones

• Apply OLS to the $\sigma_t$ equation, plug the results into the $r_t$ equation and the run again OLS or ML. The first part should not be a problem as the volatility is constant, $\zeta$, and OLS should be good rearraning in this way

$$\frac{\sigma_t-\sigma_{t-1}}{\sqrt{\sigma_t}}=\frac{k(\theta-\sigma_t)}{\sqrt{\sigma_t}}+\zeta\varepsilon_t$$ In the second equation ($r_t$ with $\sigma_t$ plugged in) I don't think OLS are a good choice because the model is heteroskedastic and $\varepsilon_t$ and $\epsilon_t$ may be correlated a priori.

Do you have any suggestion for the choice of the calibration technique? Moreover, how can I practical implement this calibration? I can use MatLab, Stata and Python (which would be my optimal choice if possible).

I have attached below a Python implementation of the CIR calibration step by MLL. After this has been run you can continue to estimate $\mu$ (drift) and $\rho$ (correlation between Brownian motions).

# CIR Calibration Script
import numpy as np
from numpy import abs, r_, var, ones, sqrt, exp,log, real, mean
from numpy.linalg import solve
from scipy.optimize import fmin
from scipy.special import ive

# Parameters
t_ = len(data) # Length time series
p = (1 / t_) * ones((1, t_))  # Probabilities
t_obs = 252 * [...]  # Years
delta_t = [...]  # Time-step

def CIR_LOGLIK(parameters, dt, data, p):
kappa = parameters
s2 = parameters
eta = parameters
c = (2 * kappa) / ((eta ** 2) * (1 - exp(-kappa * dt)))
q = ((2 * kappa * s2) / (eta ** 2)) - 1
u = c * exp(-kappa * dt) * data[:-1]
v = c * data[1:]

maxll = np.sum(- p[1:] * log(c) + p[1:] * (u + v) - p[1:] * log(v / u) * q / 2 - p[1:] * log(
ive(q, 2 * sqrt(u * v)))
- p[1:] * abs(real(2 * sqrt(u * v))))
return maxll

def FitCIR(data, dt, parameters=None, p=None):
if parameters is None:
x = r_[ones((1, t_ - 1)), data[np.newaxis, :t_ - 1]]
ols = solve(x @ x.T, x @ data[1:t_].T)
m = mean(data)
v = var(data, ddof=1)
parameters = r_[-log(ols) / dt, m, sqrt(2 * ols * v / m)]

CIRparams, fevals = fmin(CIR_LOGLIK, parameters, args=(dt, data, p), maxfun=1000, xtol=1e-8, ftol=1e-8,
disp=False, retall=True)
return CIRparams

par_CIR = FitCIR(y[0, -t_obs:], delta_t, None, p)
kappa = par_CIR
y_ = par_CIR
eta = par_CIR