Vasicek Model Simulation with Python

In this article we will outline the Vasicek Model for interest rate derivatives pricing, describe its mathematical formulation, implement and carry out a Monte Carlo simulation using Python and discuss a few real world applications of the model in quantitative finance.

Recently on QuantStart we wrote a tutorial article that discussed the mean-reverting Ornstein-Uhlenbeck process, outlining some of its applications as well as providing some Python snippets to generate sample paths. In this article we are going to introduce the Vasicek Model, which is example of a one-factor short rate model used to model interest rate behaviour for interest rate derivatives pricing. We will outline its mathematical properties and then use Python to discretise and simulate a number of sample interest rate paths using the model.

Since the Vasicek Model is a direct application of the Ornstein-Uhlenbeck (OU) process, it is worth familiarising yourself with the prior article in order to gain some intuition as to the process dynamics.

What is the Vasicek Model?

The Vasicek Model is one of the most fundamental models in quantitative finance for modeling interest rate dynamics. It provides a framework for describing how interest rates evolve over time in a stochastic manner.

The model belongs to a family of models known as 'mean-reverting' processes, which suggests that interest rates tend to move toward a long-term average level, even though they exhibit short-term fluctuations. The Vasicek Model has gained popularity due to its relatively simple structure and the ease with which it can be used to simulate sample interest rate paths.

The Vasicek Model was introduced by Oldrich Vasicek in 1977. It is formulated as a type of Stochastic Differential Equation (SDE), which captures the random fluctuations in interest rates as well as their tendency to revert to a long-term mean.

In mathematical terms, the Vasicek Model can be expressed as:

\begin{eqnarray} dr_t = \theta (\mu - r_t) dt + \sigma dW_t \end{eqnarray}

Where:

  • $r_t$: The instantaneous short-term interest rate at time $t$. $r_0$ represents the initial rate at the start of the observation period.
  • $\theta$: The speed of mean reversion, indicating how quickly the rate reverts to the mean. A higher value of $\theta$ means that the interest rate will return to the average faster, while a lower value indicates a slower adjustment.
  • $\mu$: The long-term mean level of the interest rate. This is the central value around which the interest rate fluctuates. It represents the average level to which the interest rate is expected to revert over time.
  • $\sigma$: The volatility or standard deviation of interest rate changes. This parameter controls the level of randomness or uncertainty in the movement of the interest rate. Higher volatility means more unpredictable and larger fluctuations in the interest rate.
  • $W_t$: A Wiener process (standard Brownian motion) that introduces randomness

Note: We have replicated the same notation from the previous article on the Ornstein-Uhlenbeck process and it can be seen that the equations are identical.

This equation essentially means that the change in the interest rate, $dr_t$, consists of two parts: a deterministic part, $\theta (\mu - r_t) dt$, which pulls the rate back towards the mean, and a stochastic part, $\sigma dW_t$, which introduces random fluctuations.

Applications in Quantitative Finance

The Vasicek Model has various practical applications in quantitative finance, particularly in the field of fixed income and interest rate modeling. The model can be used to price zero-coupon bonds and coupon bonds by modeling the evolution of interest rates over time. It also provides a framework to understand the term structure of interest rates. The Vasicek Model is also often used to price interest rate derivatives such as interest rate swaps, swaptions, and caps/floors.

It is also a key tool in financial risk management, since financial institutions use the Vasicek Model for assessing the risk of interest rate movements on portfolios of fixed-income assets and for determining the probability of default on various assets.

Advantages and Limitations of the Model

The Vasicek model has several advantages that make it a popular choice in quantitative finance, but it also comes with certain limitations. They advantages are listed below:

  • Mean reversion: The model captures the real-world behavior of interest rates, which often tend to revert to a long-term average. This feature makes it more realistic compared to models that assume constant rates.
  • Simplicity: Despite being mathematically grounded, the Vasicek model is relatively simple and easy to implement, making it suitable for practical applications.
  • Closed-form solution: The Vasicek model has a known analytical solution, which allows for easier analysis and computation.

The Vasicek Model also possesses some limitations:

  • Negative interest rates: One significant limitation of the Vasicek model is that it allows for the possibility of negative interest rates, which is not always realistic in practice. In reality, interest rates are generally non-negative, although some recent periods have seen negative rates in certain economies.
  • Constant volatility: The model assumes a constant level of volatility over time (i.e. $\sigma$ is a fixed value), which is not an accurate representation of real financial markets, where volatility can change in response to various economic events.
  • Lack of flexibility: The Vasicek model may be too simplistic for capturing complex interest rate dynamics, particularly in markets where rates exhibit more intricate behaviors. More advanced models, such as the Cox-Ingersoll-Ross (CIR) model, address some of these limitations.

Mathematical Overview

In this section we will delve into the mathematical aspects of the Vasicek Model, exploring its structure and deriving an analytical solution.

The Vasicek model is a type of stochastic differential equation (SDE) that describes the evolution of the short-term interest rate $r_t$ over time. Let us recap the SDE from above:

\begin{eqnarray} dr_t = \theta (\mu - r_t) dt + \sigma dW_t \end{eqnarray}

The model essentially has two components. The first is the deterministic term $\theta (\mu - r_t) dt$, which drives the mean-reverting behavior. The second is the stochastic term $\sigma dW(t)$, which introduces randomness into the system. The Vasicek model suggests that the short-term interest rate will be pulled towards the mean interest rate $\mu$ at a rate governed by $\theta$, while still being subject to random fluctuations determined by $\sigma$.

The mean-reverting behaviour is captured by the term $\theta (\mu - r_t)$, where:

  • If $r_t > \mu$, the term $(\mu - r_t)$ is negative, causing $r_t$ to decrease.
  • If $r_t < \mu$, the term $(\mu - r_t)$ is positive, causing $r_t$ to increase.

Hence, the parameter $\theta$ determines how strongly the interest rate is pulled back towards $\mu$. Larger values of $\theta$ indicate a faster reversion to the mean, while smaller values suggest a slower adjustment.

The stochastic term $\sigma dW_t$ represents the random fluctuations in the interest rate, where $W_t$ is a Wiener process with properties:

\begin{eqnarray} E[dW_t] = 0 \quad \text{and} \quad Var(dW_t) = dt \end{eqnarray}

This component ensures that the Vasicek model can capture the unpredictable movements of interest rates over time, making it realistic for financial modeling.

Analytical Solution and Properties of the Model

The Vasicek model is one of the few stochastic processes that can be solved analytically. By applying Itô's lemma we can derive an explicit solution for $r_t$. The solution to the SDE is given by:

\begin{eqnarray} r_t = r_0 e^{-\theta t} + \mu(1 - e^{-\theta t}) + \sigma \int_0^t e^{-\theta (t - s)} dW_s \end{eqnarray}

Where:

  • $r_0 e^{-\theta t}$: The influence of the initial rate $r_0$ diminishes over time.
  • $\mu (1 - e^{-\theta t})$: Represents the long-term mean component, which gradually dominates as $t \to \infty$.
  • $\sigma \int_0^t e^{-\theta (t - s)} dW_s$: The accumulated impact of the stochastic shocks over time.

Note: This is exactly the same solution as for the previous article on Ornstein-Uhlenbeck simulation.

As $t \to \infty$ the variance approaches:

\begin{eqnarray} Var(r_\infty)) = \frac{\sigma^2}{2\theta} \end{eqnarray}

This demonstrates that the Vasicek Model has stationary distributions, meaning that it does not grow to infinity.

Python Implementation and Simulation

In this section, we will demonstrate how to implement and simulate the Vasicek model using Python. We will first set up the necessary environment, then implement the model from scratch, simulate interest rate paths using Monte Carlo methods and then finally visualise the results.

Setting Up the Python Environment

Firstly, if not already present in your Python environment, it will be necessary to install the following Python libraries using pip:

pip install numpy scipy pandas matplotlib

Once installed, you can import them in your Python script:

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.stats import norm

Implementing the Vasicek Model

Next, we will implement the Vasicek model equation. Recall from above that the Vasicek model is represented by the following stochastic differential equation (SDE):

\begin{eqnarray} dr_t = \theta (\mu - r_t) dt + \sigma dW_t \end{eqnarray}

The Euler-Maruyama method, a numerical technique for approximating solutions to SDEs, will be used to discretise the Vasicek model over a set of time steps. This discretisation approximates the continuous process by considering small, discrete time steps $\Delta t$. The discretisation is as follows:

\begin{eqnarray} r_{t+\Delta t} = r_t + \theta (\mu - r_t) \Delta t + \sigma \sqrt{\Delta t} \epsilon_t \end{eqnarray}

Here, $\epsilon_t$ is a random variable drawn from a standard normal distribution (i.e., $\epsilon_t \sim \mathcal{N}(0, 1)$). This discretisation allows us to iteratively compute the value of $r_t$ over time, simulating the behavior of the Vasicek Model process.

Note: This is an identical discretisation to that of the Ornstein-Uhlenbeck process covered in the previous article.

Below is a Python function to implement this:

def vasicek_model(theta, mu, sigma, r0, T, dt):
    """
    Simulate the Vasicek model using the Euler-Maruyama method.

    Parameters
    ----------
    theta : `float`
        The speed of mean reversion
    mu : `float`
        The long-term mean level
    sigma : `float`
        The volatility of the random aspect
    r0 : `float`
        The initial interest rate
    T : `float`
        The total time horizon
    dt : `float`
        The time step size

    Returns
    -------
    `np.ndarray`
        The array of simulated interest rates over time
    """
    N = int(T / dt)  # Number of time steps
    rates = np.zeros(N)  # Pre-allocate an array large enough for the sample path
    rates[0] = r0  # Set the initial rate

    for t in range(1, N):  # Skip the initial rate (start at 1, not 0)
        dr = theta * (mu - rates[t - 1]) * dt + sigma * np.sqrt(dt) * np.random.normal()
        rates[t] = rates[t - 1] + dr
    return rates

This function simulates the interest rate $r_t$ from time $t = 0$ to $T$ with a specified step size $dt$. As mentioned above, the parameters $\theta$, $\mu$, and $\sigma$ represent the mean reversion speed, long-term mean, and volatility, respectively. The initial interest rate is $r_0$.

Simulating Interest Rate Paths Using Monte Carlo Methods

Monte Carlo simulation involves generating multiple interest rate paths to understand the behavior of the Vasicek model under different scenarios. This technique provides insights into the variability and uncertainty of the interest rate evolution. Below is an example of how to simulate multiple paths:

def simulate_vasicek_paths(theta, mu, sigma, r0, T, dt, num_simulations):
    """
    Simulate multiple interest rate paths using the Vasicek model.

    Parameters
    ----------
    theta : `float`
        The speed of mean reversion
    mu : `float`
        The long-term mean level
    sigma : `float`
        The volatility of the random aspect
    r0 : `float`
        The initial interest rate
    T : `float`
        The total time horizon
    dt : `float`
        The time step size
    num_simulations : `int`
        The number of Monte Carlo simulations to run

    Returns
    -------
    `pd.DataFrame`
        The Pandas DataFrame containing simulated interest rate paths
    """
    N = int(T / dt)  # Number of time steps
    all_simulations = np.zeros((N, num_simulations))  # Pre-allocate the two-dimensional array of multiple paths
    
    for i in range(num_simulations):
        all_simulations[:, i] = vasicek_model(theta, mu, sigma, r0, T, dt)
    
    return pd.DataFrame(
        all_simulations,
        columns=[f'Simulation {i+1}' for i in range(num_simulations)]
    )

This function simulates multiple paths using the previously defined vasicek_model function. Each column in the resulting DataFrame represents an independent simulation of the interest rate paths.

Visualizing the Results and Analyzing Model Behavior

The following function utilises Matplotlib to plot the Monte Carlo simulation of multiple paths:

def plot_vasicek_paths(df, T, dt):
    """
    Visualize simulated interest rate paths.

    Parameters
    ----------
    df : `pd.DataFrame`
        The Pandas DataFrame containing simulated paths
    T : `float`
        The total time horizon
    dt : `float`
        The time step size

    Returns
    -------
    None
    """
    plt.figure(figsize=(12, 6))
    time_points = np.linspace(0, T, int(T / dt))

    for column in df.columns:
        plt.plot(time_points, df[column], lw=1.0, alpha=0.6)

    plt.title('Simulated Vasicek Interest Rate Paths')
    plt.xlabel('Time')
    plt.xlim(0.0, 1.0)
    plt.ylabel('Interest Rate')
    plt.show()

We can then run the entire code with the following example parameter values, by setting $\theta = 2.0$, $\mu = 0.05$, $\sigma = 0.02$, and $r_0 = 0.03$ to simulate 10 different interest rate paths over a one-year period with a step size $dt = 0.001$:

if __name__ == "__main__":
    theta = 2.0  # Speed of mean reversion
    mu = 0.05  # Long-term mean interest rate (5%)
    sigma = 0.02  # Volatility of the random component (2%)
    r0 = 0.03  # Initial value of the interest rate (3%)
    T = 1.0  # Time horizon (1 year)
    dt = 0.001  # Time step size
    num_simulations = 10  # Number of sample paths to generate

    # Simulate the paths
    simulated_paths = simulate_vasicek_paths(theta, mu, sigma, r0, T, dt, num_simulations)
    
    # Plot the results
    plot_vasicek_paths(simulated_paths, T, dt)

The output of this code is visualised in the following figure:

Monte Carlo simulation of multiple sample paths from the Vasicek Model plotted with Python

We have deliberately set the initial interest rate at 3% and the longer-term mean interest rate at 5%, in order to highlight the mean-reverting behaviour of the Vasicek Model. It can be seen that all of the paths generated initially trend upward away from the 3% initial value, but soon revert towards the longer term rate at 5% near $T=1$. The mean reversion speed $\theta=2.0$ has been chosen specifically to highlight the rapid reversion to the mean once the path deviates from the mean value.

Summary and Next Steps

In this article we described the Vasicek Model, as an application of an Ornstein-Uhlenbeck process. We described its mathematical formulation and derived its analytical solution. We provided a basic Python implementation of a Monte Carlo simulation of multiple realisations of discretised sample paths from this model. In subsequent articles we will extend this model to the Cox-Ingersoll-Ross (CIR) and Hull-White models that add more realistic interest rate behaviours.

Full Code

# vasicek-model.py

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.stats import norm


def vasicek_model(theta, mu, sigma, r0, T, dt):
    """
    Simulate the Vasicek model using the Euler-Maruyama method.

    Parameters
    ----------
    theta : `float`
        The speed of mean reversion
    mu : `float`
        The long-term mean level
    sigma : `float`
        The volatility of the random aspect
    r0 : `float`
        The initial interest rate
    T : `float`
        The total time horizon
    dt : `float`
        The time step size

    Returns
    -------
    `np.ndarray`
        The array of simulated interest rates over time
    """
    N = int(T / dt)  # Number of time steps
    rates = np.zeros(N)  # Pre-allocate an array large enough for the sample path
    rates[0] = r0  # Set the initial rate

    for t in range(1, N):  # Skip the initial rate (start at 1, not 0)
        dr = theta * (mu - rates[t - 1]) * dt + sigma * np.sqrt(dt) * np.random.normal()
        rates[t] = rates[t - 1] + dr
    return rates


def simulate_vasicek_paths(theta, mu, sigma, r0, T, dt, num_simulations):
    """
    Simulate multiple interest rate paths using the Vasicek model.

    Parameters
    ----------
    theta : `float`
        The speed of mean reversion
    mu : `float`
        The long-term mean level
    sigma : `float`
        The volatility of the random aspect
    r0 : `float`
        The initial interest rate
    T : `float`
        The total time horizon
    dt : `float`
        The time step size
    num_simulations : `int`
        The number of Monte Carlo simulations to run

    Returns
    -------
    `pd.DataFrame`
        The Pandas DataFrame containing simulated interest rate paths
    """
    N = int(T / dt)  # Number of time steps
    all_simulations = np.zeros((N, num_simulations))  # Pre-allocate the two-dimensional array of multiple paths
    
    for i in range(num_simulations):
        all_simulations[:, i] = vasicek_model(theta, mu, sigma, r0, T, dt)
    return pd.DataFrame(all_simulations, columns=[f'Simulation {i+1}' for i in range(num_simulations)])


def plot_vasicek_paths(df, T, dt):
    """
    Visualize simulated interest rate paths.

    Parameters
    ----------
    df : `pd.DataFrame`
        The Pandas DataFrame containing simulated paths
    T : `float`
        The total time horizon
    dt : `float`
        The time step size

    Returns
    -------
    None
    """
    plt.figure(figsize=(12, 6))
    time_points = np.linspace(0, T, int(T / dt))

    for column in df.columns:
        plt.plot(time_points, df[column], lw=1.0, alpha=0.6)

    plt.title('Simulated Vasicek Interest Rate Paths')
    plt.xlabel('Time')
    plt.xlim(0.0, 1.0)
    plt.ylabel('Interest Rate')
    plt.show()


if __name__ == "__main__":
    theta = 2.0  # Speed of mean reversion
    mu = 0.05  # Long-term mean interest rate (5%)
    sigma = 0.02  # Volatility of the random component (2%)
    r0 = 0.03  # Initial value of the interest rate (3%)
    T = 1.0  # Time horizon (1 year)
    dt = 0.001  # Time step size
    num_simulations = 10  # Number of sample paths to generate

    # Simulate the paths
    simulated_paths = simulate_vasicek_paths(theta, mu, sigma, r0, T, dt, num_simulations)
    
    # Plot the results
    plot_vasicek_paths(simulated_paths, T, dt)