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).


1 Answer 1


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[0]
    s2 = parameters[1]
    eta = parameters[2]
    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[1]) / dt, m, sqrt(2 * ols[1] * v / m)]

    CIRparams, fevals = fmin(CIR_LOGLIK, parameters, args=(dt, data, p[0]), 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[0]
y_ = par_CIR[1]
eta = par_CIR[2]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.