Tidy Finance - CIR Model Calibration using Python (2024)

CIR Model - Overview

The CIR model aims to capture the dynamics of interest rates, offering a powerful alternative to the Vasicek model.

The Vasicek model can be mathematically defined by an Ornstein–Uhlenbeck process with an added drift term. This is a stochastic process particularly suitable for interest rate dynamics, owing to its inherent properties of stationarity and mean-reversion. These attributes align with the behavior commonly observed in interest rates.

The SDE (stochastic differential equation) for the Vasicek model is the following (where \(dW_t\) denotes the Wiener process): \[dr_t=k(\theta-r_t)dt+\sigma dW_t\]

Moreover, the model yields an intuitive interpretation of its parameters:

• \(\theta\) is the long-run average

• \(k\) is the intensity with which the process returns to its long-run average

• \(\sigma\) is the instantaneous volatility

In contrast, the CIR model differs by introducing a novel component into the diffusion part of the process, precisely \(\sqrt{r_t}\), as depicted below: \[dr_t=k(\theta-r_t)dt+\sigma\sqrt{r_t}dW_t\]

This seemingly minor adjustment has profound implications, as it ensures that interest rates remain strictly positive. However, this modification also engenders notable distinctions between the two models.

Without delving into mathematical intricacies, it is imperative to recognize that while the Vasicek model adheres to a Gaussian process, the CIR model deviates from this, leading to a substantially more intricate conditional distribution. However, it is not necessary to grasp the intricacies of the distribution to develop Ordinary Least Squares (OLS) estimates, as elucidated in this article.

Model Discretization

The initial step necessitates discretization to manage the CIR model effectively and derive OLS estimates. This task is efficiently achieved through the Euler-Maruyama Scheme, a numerical method enabling the transformation of the stochastic differential equation (SDE) into a discrete-time equation. Essentially, this scheme involves approximating each differential term as a finite difference, with the accuracy of the approximation improving as the time step decreases.

This leads to the following equation: \[r_{t + \delta t} - r_t = k(\theta - r_t)\delta t + \sigma \sqrt{r_t}N(0, \delta t)\]

Following some straightforward manipulations, it can be rewritten as: \[\frac{{r_{t + \delta t} - r_t}}{{\sqrt{r_t}}} = \frac{k\theta \delta t}{\sqrt{r_t}} - k\sqrt{r_t} \delta t + \sigma \sqrt{\delta t} N(0,1)\]

This can be interpreted as: \[y_i = \beta_1 z_{1,i} + \beta_2 z_{2,i} + \epsilon_i\]

Where: \[\begin{aligned} y_i &= \frac{r_{t+\delta t}-r_t}{\sqrt{r_t}} \\ \beta_1 &= k\theta \\ \beta_2 &= -k \\ z_{1,i} &= \frac{\delta t}{\sqrt{r_t}} \\ z_{2,i} &= \sqrt{r_t}\delta t \\ \epsilon_i &= \sigma \sqrt{\delta t} N(0,1)\end{aligned}\]

Estimating the parameters of interest then becomes straightforward: \[\begin{aligned}\hat{k} &= -\hat{\beta_2} \\\hat{\theta} &= \frac{\hat{\beta_1}}{\hat{k}} \\\hat{\sigma^2} &= \frac{\hat{Var(\epsilon)}}{\delta_t}\end{aligned}\]

Python Implementation

Upon reviewing the mathematical details, the next step is to translate them into Python code. I will define two functions:

• simulate_cir(): This function simulates a path based on theoretical parameters. It is important to note that in the CIR model, the term \(\sqrt{r_t}\) requires that \(r_t\) remains non-negative. However, since we are working with a discrete form of the model, simulation errors may arise due to negative rates stemming from the random component of the model. These errors need to be addressed. This is facilitated by the following piece of code: np.sqrt(max(0, r)), which guarantees that no negative interest rate is passed under the square root.
• ols_cir(): This function performs the estimates based on a given array of values.

# import librariesimport numpy as npfrom sklearn.linear_model import LinearRegressionimport matplotlib.pyplot as pltimport seaborn as sns# function to simulate pathsdef simulate_cir(k, theta, sigma, r0, T, N): # populate an empty array dt = T / N interest_rate_paths = np.zeros(N+1) interest_rate_paths[0] = r0  for t in range(1, N+1): Z = np.random.randn() r = interest_rate_paths[t-1] interest_rate_paths[t] = r + k * (theta-r) * dt + sigma * np.sqrt(dt) * np.sqrt(max(0, r)) * Z  return interest_rate_paths# function to estimate parametersdef ols_cir(data, dt): # define variables Nsteps = len(data) rs = data[:Nsteps - 1]  rt = data[1:Nsteps]  # model initialization model = LinearRegression() # feature engineering to fit the theoretical model y = (rt - rs) / np.sqrt(rs) z1 = dt / np.sqrt(rs) z2 = dt * np.sqrt(rs) X = np.column_stack((z1, z2)) # fit the model model = LinearRegression(fit_intercept=False) model.fit(X, y) # calculate the predicted values (y_hat), residuals and the parameters y_hat = model.predict(X) residuals = y - y_hat beta1 = model.coef_[0]  beta2 = model.coef_[1] # get the parameter of interest for CIR k0 = -beta2 theta0 = beta1/k0 sigma0 = np.std(residuals)/np.sqrt(dt)  return k0, theta0, sigma0

The subsequent step involves applying these functions to a specific case by fixing theoretical parameters for the CIR model. In this instance, I will opt for relatively standard parameter values. However, a more formal analysis would require exploring a range of values for each parameter to evaluate the robustness of the estimators under different conditions.

k_true = 5 # True mean reversion speedtheta_true = 0.05 # True long-run meansigma_true = 0.03 # True volatility of interest ratesr0_true = 0.3 # True initial interest rateT = 1 # Time horizonN = 100 # Number of time stepsdt = T/N# simulated pathnp.random.seed(123)series = simulate_cir(k_true, theta_true, sigma_true, r0_true, T, N)estimates = ols_cir(series, dt)# print resultsprint(f"The theoretical parameters are: k={k_true}, theta={theta_true}, sigma={sigma_true}")print(f"The estimates are: k={round(estimates[0], 3)}, theta={round(estimates[1], 3)}, sigma={round(estimates[2], 3)}")

The theoretical parameters are: k=5, theta=0.05, sigma=0.03The estimates are: k=5.078, theta=0.051, sigma=0.034

Monte-Carlo Results

The final step entails conducting a Monte Carlo Simulation to thoroughly scrutinize how the estimators perform. The model comprises two components: a deterministic one and a stochastic one. In the preceding step, I calibrated the model using a single vector of Gaussian samples, which generated a particular path. Now, I will repeat this process with several different vectors of Gaussian samples, generating multiple paths while maintaining the same real parameters.

If the estimators are reasonably accurate, we should observe a convergence of the Monte Carlo estimates means towards the real parameters, indicating the unbiasedness of the estimators in principle.

It is important to note that due to discretization errors, there is a possibility of generating paths with negative interest rates. To address this issue, I will adopt the approach suggested in Orlando, Mininni, and Bufalo (2020)3, which involves shifting all series by the 99th percentile before the calibration. This allows, in most cases, to work with a positive series of data. It is crucial to highlight that if we intend to use the calibrated model for predictions, we must subsequently subtract the same quantity from the predictions. This adjustment is implemented in the code via an if statement within the while loop.

# define variablesn_scenarios = 10000paths = np.zeros(shape=(n_scenarios, N+1))OLS_estimates = np.zeros(shape=(n_scenarios,3))# perform simulationsnp.random.seed(123)count = 0while count < n_scenarios: # fill the array with the simulations series = simulate_cir(k_true, theta_true, sigma_true, r0_true, T, N)  # Check if there's at least one negative number in the series if np.any(series < 0): # Calculate the value of the 99th percentile percentile_99 = np.percentile(series, 99) # Add the value of the 99th percentile to each element of the series series += percentile_99 paths[count, :] = series  # OLS estimation estimates = ols_cir(series, dt) OLS_estimates[count,:] = estimates  # update counter count += 1# plot of the estimates distributionscolours = ["blue", "red", "green"] parameters = ["k", "Theta", "Sigma"] real_values = [k_true, theta_true, sigma_true]# Set the seaborn stylesns.set_style("whitegrid")for i in range(3): # Using range(3) for iterating over indices # Set up the plot plt.figure(figsize=(7, 5)) sns.histplot(OLS_estimates[:, i], bins=500, color=colours[i], alpha=0.7, kde=True) # Plot mean as black dot line mean_estimate = np.mean(OLS_estimates[:, i]) plt.axvline(mean_estimate, color="black", linestyle="--", label="Monte-Carlo Mean") # Plot specified value as purple dot line plt.axvline(real_values[i], color="purple", linestyle="-.", label="Theoretical Value") plt.xlabel("Estimates", fontsize="medium") plt.ylabel("Frequency", fontsize="medium") plt.title(parameters[i], fontsize="large") plt.legend(title_fontsize="small", fontsize="x-small") plt.show()

Conclusions and Real-World Application

To summarise, this article has covered:

1. Introduction to the CIR model

2. Implementation of OLS estimators for model calibration

3. Evaluation of estimator performance via Monte Carlo simulations

To conclude, let’s include a final snippet code to apply the routine to real-world data. I will fetch the 13-week Treasury Bill Rate using the yfinance package and then apply the ols_cir function to calibrate the model. Here are a couple of insights on the code:

• Data covers all of 2023
• Following convention, dt is set to 1/252 since there are approximately 252 trading days in a year, and interest rates are expressed in annual terms

import yfinance as yf# Ticker symbol for US interest rateticker_symbol = "^IRX" # 13 Week Treasury Bill Rate# Fetch datainterest_rate_data = yf.download(ticker_symbol, start="2023-01-01", end="2024-01-01")interest_rate_data = interest_rate_data["Close"]/100 # divide by 100 as the data are expressed in percentagedt = 1/252# model calibrationestimates = ols_cir(interest_rate_data.values, dt)print(f"The estimates are: k={round(estimates[0], 3)}, theta={round(estimates[1], 3)}, sigma={round(estimates[2], 3)}")

[*********************100%%**********************] 1 of 1 completedThe estimates are: k=6.442, theta=0.052, sigma=0.03

Figure1: Results of Monte Carlo Simulation for parameter k.Figure2: Results of Monte Carlo Simulation for parameter Theta.Figure3: Results of Monte Carlo Simulation for parameter Sigma.