Vasicek Model for Interest Rates

Reference: https://en.wikipedia.org/wiki/Vasicek_model (1977)

Multidimensional Stochastic Processes

import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
import pandas as pd
import math

Do also:

• Hull–White model
• CIR model
• Merton (Ho-Lee) model, Black-Karasinski (Black-Derman-Toy)

It is a mathematical model describing the evolution of interest rates. It is a one-factor short rate model, interest rate movements are driven by only one source of market risk. The model can be used in the valuation of interest rate derivatives.

It is basically an AR(1) process, GBM is modified to introduce mean reversion. Interest rates exist in a limited range.

SDE:

$dr_t = a(b - r_t)dt + \sigma\sqrt{dt} z_t $

Terminal value asymptotic distribution is $ N(b, \sigma^2/(2a))$.

b = long-term mean

a = speed of reversion

$\sigma^2/(2a) $ = long-term volatility

We do not need to impose that $r_t$ be positive values.

There is an application to European short-term interest rates below.

#Simulate I paths: 
T = 500
I = 100
dt = 1/50

#What are plausible parameter values?
sigma = 0.05
r0 = 0.1
a = 2.0
b= 0.05
theta = [a,b,sigma]

def Vasicek(theta):
    r = np.zeros((T+1, I))
    r[0, :] = r0

    for t in range(1, T+1):
        z=npr.standard_normal(I)
        r[t, :]=r[t-1,:] +theta[0]*(theta[1]-r[t-1,:])*dt+theta[2]*np.sqrt(dt)*z
    return r                                                                                                           
                                                                                                               
r = Vasicek(theta)  

LTStd = sigma/math.sqrt(2*a); LTStd

0.025

df = pd.DataFrame(r)
df.describe()
# Google total re

	0	1	2	3	4	5	6	7	8	9	...	90	91	92	93	94	95	96	97	98	99
count	501.000000	501.000000	501.000000	501.000000	501.000000	501.000000	501.000000	501.000000	501.000000	501.000000	...	501.000000	501.000000	501.000000	501.000000	501.000000	501.000000	501.000000	501.000000	501.000000	501.000000
mean	0.059204	0.054405	0.051587	0.047544	0.050570	0.054908	0.037890	0.055462	0.055035	0.056309	...	0.046472	0.044013	0.056133	0.047333	0.065450	0.058638	0.054307	0.057164	0.056380	0.067504
std	0.019852	0.024084	0.022124	0.022873	0.027052	0.024540	0.025908	0.023603	0.027678	0.026316	...	0.032923	0.022456	0.020498	0.024077	0.020605	0.032983	0.023667	0.017914	0.023256	0.025123
min	0.000888	-0.004094	-0.015962	-0.005261	-0.019419	0.002830	-0.025742	0.002684	-0.027665	-0.001042	...	-0.022755	-0.005974	-0.000160	-0.029049	0.021355	-0.023941	0.001382	0.003884	-0.002807	0.003058
25%	0.045800	0.038191	0.034475	0.030246	0.031330	0.036914	0.018555	0.040508	0.036093	0.036223	...	0.020800	0.027683	0.043369	0.035306	0.050765	0.035458	0.036575	0.044896	0.037637	0.051428
50%	0.057479	0.053558	0.054008	0.044126	0.053603	0.051413	0.033922	0.052219	0.054166	0.052734	...	0.045622	0.042219	0.057408	0.049285	0.064577	0.061600	0.055739	0.056636	0.058870	0.069294
75%	0.072240	0.068289	0.067567	0.064948	0.067051	0.068173	0.056600	0.070473	0.076379	0.074553	...	0.070734	0.058411	0.070588	0.060916	0.077397	0.085859	0.070212	0.069149	0.073330	0.083729
max	0.120733	0.120920	0.106779	0.100000	0.105701	0.131934	0.101044	0.120201	0.118075	0.131398	...	0.122284	0.114923	0.104723	0.109325	0.126098	0.121032	0.115072	0.113119	0.118506	0.133691

8 rows × 100 columns

rT = r[-1]

fig, ax1 = plt.subplots()#
ax1.hist(rT, bins=50); ax1.set_title('Vasicek process - Normal Distribution')

Text(0.5, 1.0, 'Vasicek process - Normal Distribution')

df = pd.DataFrame(rT)
df.describe()

	0
count	100.000000
mean	0.051208
std	0.025854
min	-0.005500
25%	0.037485
50%	0.047219
75%	0.067146
max	0.131666

plt.plot(r[:, :10], lw = 1.5)
plt.xlabel('# periods'); plt.ylabel('Index Level'); plt.title("Vasicek Process with mean 0.05")
plt.grid(True)

#Autocorrelation is very high among state variables and dies out very slowly.
np.corrcoef(r[:-1,0], r[1:,0])

array([[1.        , 0.93207337],
       [0.93207337, 1.        ]])

#Correlation should be negative between state variables and subsequent return, and proportional with a and rt.
dr = r[1:,0]/r[:-1,0]-1
dr.shape
np.corrcoef(r[:-1,0], dr) #?

array([[ 1.        , -0.18279462],
       [-0.18279462,  1.        ]])

Application to Short-Term Interest Rates

Though LIBOR, Euribor, and the federal funds rate are concerned with the same action, i.e. interbank loans, there is an important distinction: the federal funds rate is a target interest rate that is set by the FOMC for implementing U.S. monetary policies.

For this reason I shall use Eonia. The Eonia rate is the 1-day interbank interest rate for the Euro zone. In other words, it is the rate at which banks provide loans to each other with a duration of 1 day. Therefore Eonia can be considered as the 1 day Euribor rate.

As of 1 October 2019 EONIA is calculated with a reformed methodology tracking the euro short-term rate (EONIA: Euro Interbank Offered Rate). Following a recommendation made by the working group on euro risk-free rates on 14 March 2019, as of 2 October for the trade date 1 October 2019 the European Money Market Institute (EMMI) changed the way it calculates the EONIA. Rate for the overnight maturity is calculated as the euro short-term rate plus a spread of 8.5 basis points. Eonia will be discontinued on Jan. 3, 2022.

ESTR Data http://webstat.banque-france.fr/en/quickview.do?SERIES_KEY=285.ESTR.B.EU000A2X2A25.WT

#import csv data
lines = open('./ESTR.csv','r').readlines()
lines = [line.replace("", "") for line in lines]
lines[:4]

['Date,ESTR\n', '01.10.19,-0.549\n', '02.10.19,-0.551\n', '03.10.19,-0.555\n']

#import data from the csv file with read_csv
estr = pd.read_csv('./ESTR.csv', index_col=0, parse_dates=True, sep=",", dayfirst=True,  false_values='0',)
np.round(estr.head())

	ESTR
Date
2019-10-01	-1.0
2019-10-02	-1.0
2019-10-03	-1.0
2019-10-04	-1.0
2019-10-07	-1.0

nulls = estr[ estr['ESTR'] == '-' ].index
# Delete these row indexes from dataFrame
estr.drop(nulls , inplace=True)

estr.sort_index(axis = 0) 
estr['ESTR'] = estr['ESTR'].astype(float) #for some reason the data appeared as strings

plt.plot(estr['ESTR'])

[<matplotlib.lines.Line2D at 0x116264978>]

estr.describe()

	ESTR
count	231.000000
mean	-0.542087
std	0.007155
min	-0.555000
25%	-0.548000
50%	-0.541000
75%	-0.538000
max	-0.511000

es = np.array(estr['ESTR'])
des = es[1:]/es[:-1]-1 #% change in short term rates  

np.corrcoef(des,es[:-1])
#First-order autocorrelation is positive!

array([[1.        , 0.30008802],
       [0.30008802, 1.        ]])

#Bird's Eye View Calibration:
#Long-term mean:
b = np.mean(es)
#Long-term standard deviation is a function of a and sigma:
#LTstd = np.std(es)
#We can use the autocorrelation to find a:
#beta = -a = cov(rt, drt)/var(rt)
VCV = np.cov(es[:-1], des)
a = -VCV[1,0]/np.var(es)
sigma = np.std(es)*np.sqrt(2*a)

print(round(a,5), round(b,5), round(sigma,5))

#Oops! First-order autocorrelation is positive, so it is not a Vasicek process.

-0.33843 -0.54209 nan

/Users/purnur/miniconda3/envs/website/lib/python3.6/site-packages/ipykernel_launcher.py:10: RuntimeWarning: invalid value encountered in sqrt
  # Remove the CWD from sys.path while we load stuff.

import statsmodels

from statsmodels import tsa

#ACF = tsa.stattools.acf(es, nlags=10, fft=False)
#pd.DataFrame(ACF, columns=['ACF'])

#PACF = tsa.stattools.pacf(es, nlags=10)
#pd.DataFrame(PACF, columns=['PACF'])

# statsmodels.tsa.stattools.pacf(es, nlags=10)
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# Create figure
fig, (ax1, ax2) = plt.subplots(2,1, figsize=(8,8))
# Make ACF plot
plot_acf(es, lags=10, zero=False, ax=ax1)
# Make PACF plot
plot_pacf(es, lags=10, zero=False, ax=ax2)
plt.show()

Indeed, it seems the ESTR process is an AR(2) process, therefore it cannot be modeled by the Vasicek model, which fits an AR(1) process.

#One example of a Vasicek model with random parameters
T = len(es)-1
I = 1
dt = 1/50
theta = [0, 0.04, 0.2]

r0 = es[0]

r = Vasicek(theta)

plt.plot(r, 'r')

[<matplotlib.lines.Line2D at 0x1a1917d7f0>]

plt.plot(es, r,'b.')

[<matplotlib.lines.Line2D at 0x1a19239710>]

#Just for the sake of the exercise, I shall try to fit The Vasicek model parameters using fmin.

ordd = 2
N=10

f = lambda theta: Vasicek(theta) #I minimize the norm for each value of the short-term rate

#g = lambda theta: np.linalg.norm(es-f(theta), ord=ordd)**ordd
def g(theta):
    g=0
    for n in range(0,N):
        g += np.linalg.norm(es-f(theta), ord=ordd)**ordd
    return g

import scipy.optimize
thetaopt = scipy.optimize.fmin(func=g, x0=[0.,0.,0.],ftol=1e-7,xtol=1e-7)
print('Parameters a,b,sigma for Euro short-term rate', thetaopt)

Warning: Maximum number of function evaluations has been exceeded.
Parameters a,b,sigma for Euro short-term rate [ 0.0051964  -0.00347781 -0.0012848 ]

ropt = Vasicek(thetaopt)
plt.plot(es, ropt,'b.')

[<matplotlib.lines.Line2D at 0x1a19313828>]

plt.plot(ropt)

[<matplotlib.lines.Line2D at 0x1a1938feb8>]

The results do not appear implausible, in spite of using the wrong model and obtaining weird parameter values like negative speed of reversion to the mean and negative instantaneous standard deviation.

#Normality tests

import scipy as sp
from scipy import stats

#Shapiro-Wilk normality test

from scipy.stats import shapiro
stat, p = shapiro(es)
print('Statistics=%.3f, p=%.3f' % (stat, p))
# interpret
alpha = 0.05
if p > alpha:
	print('Sample looks Gaussian (fail to reject H0)')
else:
	print('Sample does not look Gaussian (reject H0)')

Statistics=0.965, p=0.000
Sample does not look Gaussian (reject H0)

#Jarque-Bera normality test

JB, p = stats.jarque_bera(es)
print('Statistics=%.3f, p=%.3f' % (JB, p))
# interpret
alpha = 0.05
if p > alpha:
	print('Sample looks Gaussian (fail to reject H0)')
else:
	print('Sample does not look Gaussian (reject H0)')

Statistics=16.441, p=0.000
Sample does not look Gaussian (reject H0)