Vasicek One Factor Model

import numpy as np

def vasicek_model(r0, K, theta, sigma, T=1., N=10, seed=28):    
    np.random.seed(seed)
    dt = T/float(N)    
    rates = [r0]
    for i in range(N):
        dr = K*(theta-rates[-1])*dt + sigma*np.random.normal()
        rates.append(rates[-1] + dr)
    return range(N+1), rates

if __name__ == "__main__":
    x, y = vasicek_model(0.01875, 0.20, 0.01, 0.012, 10., 200)

    import matplotlib.pyplot as plt
    plt.plot(x,y)
    plt.show()

import numpy as np
import pandas as pd 
from typing import Any
import matplotlib.pyplot as plt

def Brownian_process(T = 1, dt = 0.001, rho = None):
  N = int(T / dt)
  if not rho:
    W = np.zeros(N) # preallocate the output array holding the sample paths with the inital point

    for i in range(1, N):
      W[i] = W[i-1] + np.random.normal(scale = np.sqrt(dt)) # Increment a random normal at each step

    return W
  
  if rho is not None:  # if rho is defined, that means that the output will be a 2-dimensional Brownian motion
    W_1 = np.zeros(N)
    W_2 = np.zeros(N)

    for i in range(1, N): # Simulate two Brownian motions and construct a joint bivariate Brownian motion with parameter rho
      Z_1 = np.random.normal(scale = np.sqrt(dt))
      Z_2 = np.random.normal(scale = np.sqrt(dt))
      Z_3 = rho * Z_1 + np.sqrt(1- rho**2) * Z_2
            
      W_1[i] = W_1[i-1] + Z_1 # First Brownian motion
      W_2[i] = W_2[i-1] + Z_3 # Second Brownian motion
            
    return [W_1, W_2]

# Example 
# The economist wants to generate discrete sample paths of two Brownian motions with a correlation coefficient of 0.6. 
# The Brownian motions needs to start at 0 at time 0 and ends at time T = 3 with an increment of 0.5.
x = range(0,6)
y = Brownian_process(3, 0.5, 0.6)
plt.plot(x, y[0], label = 'Path 1')
plt.plot(x, y[1], label = 'Path 2')
plt.xlabel('Time step')
plt.title('Two Brownian Motion Paths with Rho equal to 0.6')
plt.legend()
plt.show()

def Vasicek_One_Factor(r0: float = 0.1, a: float = 1.0, b: float = 0.1, sigma: float = 0.2, T: int = 52, dt = 0.1) -> pd.DataFrame:
  """
  :param r0: float, starting interest rate of the Vasicek process
  :param a: float, speed of reversion which is a parameter that characterizes the velocity at which such trajectories will regroup around b in time
  :param b: float, long term mean level. All future trajectories of  r will evolve around mean level b in the long run
  :param sigma: float, instantaneous volatility measures instant by instant the amplitude of randomness entering the system
  :param T: integer, end modelling time. From 0 to T the time series runs.
  :param dt: float, increment of time that the process runs on. Ex. dt = 0.1 then the time series is 0, 0.1, 0.2,...
  :return: N x 2 pd.DataFrame where index is modeling time and values are a realisation of the underlying's interest rates
  """
  N = int(T / dt)  # number of subintervals of length 1/dt between 0 and max modeling time T

  time, delta_t = np.linspace(0, T, num=N, retstep=True)

  brownian_process = Brownian_process(T, dt)

  r = np.ones(N) * r0

  for t in range(1, N):
    r[t] = r[t - 1] + a * (b - r[t - 1]) * dt + sigma * (brownian_process[t] - brownian_process[t - 1])

  dict = {'Time': time, 'Interest Rate': r}
  interest_rate_simulation = pd.DataFrame.from_dict(data=dict)
  interest_rate_simulation.set_index('Time', inplace=True)

  return interest_rate_simulation

Vasicek_One_Factor(0.1, 0.8, 0.08, 0.2, 10, 0.5)