Calibration of the Vasicek Model: An Step by Step Guide

Victor Bernal A.

April 12, 2016

victor.bernal@mathmods.eu

Abstract
In this report we present 3 methods for calibrating the OrnsteinUhlenbeck process to a data set. The
model is described and the sensitivity analysis with respect to changes in the parameters is performed.
In particular the Least Squares Method, the Maximum Likelihood Method and the Long Term Quantile
Method are presented in detail.

Introduction
The OrnsteinUhlenbeck process [3] (named after Leonard Ornstein and George Eugene Uhlenbeck), is a
stochastic process that, over time, tends to drift towards its long-term mean: such a process is called mean-
reverting. It can also be considered as the continuous-time analogue of the discrete-time AR(1) process where
there is a tendency of the walk to move back towards a central location, with a greater attraction when the
process is further away from the center.
Vasicek assumed that the instantaneous spot Interest Rate under the real world measure evolves as an Orstein-
Uhlenbeck process with constant coecients [5]. The most important feature which this model exhibits is the
mean reversion,which means that if the interest rate is bigger than the long run mean µ, then the coecient
λ makes the drift become negative so that the rate will be pulled down in the direction of µ. Similarly, if the
interest rate is smaller than the long run mean. Therefore, the coecient λ is the speed of adjustment of the
interest rate towards its long run level. This model is of particular interest in nance because there are also
compelling economic arguments in favor of mean reversion. When the rates are high, the economy tends to
slow down and borrowers require less funds. Furthermore, the rates pull back to its equilibrium value and
the rates decline. On the contrary when the rates are low, there tends to be high demand for funds on the
part of the borrowers and rates tend to increase. One unfortunate consequence of a normally distributed
interest rate is that it is possible for the interest rate to become negative.
In this article we start from the Euler Maruyana discretization scheme for the Vasicek process and the
sensitivity analysis, then we present in detail 3 well known methods (Least squares Method, Maximum
Likelihood Method and the Long Term Quantile Method) for calibrating the model's parameter to a data
set. We also refer the reader to [4] where some of these techniques are applied but using a dierent scheme.

1 Euler Scheme and Sensitivity Analysis

The stochastic dierential equation (SDE) for the Ornstein-Uhlenbeck process is given by

drt = λ (µ − rt ) dt + σdWt (1.1)

with λ the mean reversion rate, µ the mean, and σ the volatility. The solution of the model is

ˆt
rt = r0 exp (−λt) + µ (1 − exp (−λt)) + σ exp (−λt) dWt (1.2)
0

Here the interest rates are normally distributed and the expectation and variance are given by

1
 Eo [rt ] = r0 exp (−λt) + µ (1 − exp (−λt)) (1.3)
and

σ2
V ar [rt ] = (1 − exp (−2λt)) (1.4)
2λ
σ2
as t → ∞, the limit of expected rate and variance, will converge to µ and 2λ respectively. The Euler
Maruyana Scheme for this models is
√
rt+δt = rt + λ (µ − rt ) δt + σ δtN (0, 1) (1.5)
the process can go negative with probability
 √ 
P (rt+δt ≤ 0) = P rt + λ (µ − rt ) δt + σ δtN (0, 1) ≤ 0 (1.6)
  
rt + λ (µ − rt ) δt
= P N (0, 1) ≤ − √ (1.7)
σ δt
 
rt + λ (µ − rt ) δt
=Φ − √ (1.8)
σ δt
Some of the parameters play a big role in the pricing of a nancial derivatives or in the forecasting of a process,
while some of them do not aect them so much. Therefore, depending on what we want to price or forecast,
it is important to check the sensitivity of the models with respect to dierent parameters. The corresponding
sensitivity analysis is performed as presented in [6]. Lets consider a two outcomes of a process which dier
exclusively in a perturbation of one of the parameters (but they have the same stochastic realization N (0, 1)),
then for λ we have that
√
r̃t+δt = rt + [λ + ∆λ] (µ − rt ) δt + σ δtN (0, 1) (1.9)
so
r̃t+δt − rt+δt = ∆λ (µ − rt ) δt (1.10)
when λ is increased the variance (1.4) decreases. So, the change in the reversion coecient will not aect
the short rate (1.3) in long term, just eect the time which is necessary for the interest rate to come back to
the long term mean. Therefore, λ is important in the pricing of the nancial instruments which are aected
by the volatility (1.4), but are not dependent on the long term expected value (1.3) of the simulated interest
rate. For µ is performed as
√
r̃t+δt = rt + λ ([µ + ∆µ] − rt ) δt + σ δtN (0, 1) (1.11)
so

r̃t+δt − rt+δt = λ (∆µ − rt ) δt (1.12)

for σ is performed as
√
r̃t+δt = rt + λ (µ − rt ) δt + [σ + ∆σ] δtN (0, 1) (1.13)
so √
r̃t+δt − rt+δt = ∆σ δtN (0, 1) (1.14)
Because of the standard Brownian motion, in the long term, the eect of the change in σ does not aect the
expected value of the interest rate (1.3), but it increases the variance (1.4).
%% Euler Scheme Vasicek
clc
clear a l l
close a l l

%% Set t h e seed a t 123

2
 Vasicek Process
0.9

0.8

0.7

0.6

Interest rate
0.5

0.4

0.3

0.2

0.1
0 20 40 60 80 100
time

Figure 1.1: Vasicek Process

rng ( 1 2 3 )

%% Parameters o f t h e Model
lambda = 0 . 3 ;% R e v e r t i o n c o e f f i c i e n t
N=3; % Number o f s i m u l a t i o n s
mu= 0 . 7 ; % Long term Mean
sigma = 0 . 0 5 ; % V o l a t i l i t y
d e l t a _ t =1; % Time s t e p
T=100; % Time l e n g h t
c o l o r =[ ' b ' , ' r ' , 'm' ] ; % c o l o r
n=T. / d e l t a _ t ; % Number o f time s t e p s
j =1;

%% S i m u l a t i n g Vasicek Euler Scheme

S0 = 0 . 1 ; % S t a r t i n g p o i n t
i =2;
S(1)= S0 ; % S t a r t i n g p o i n t
while i <n+1
S ( i )= S ( i − 1) + lambda . ∗ ( mu−S ( i − 1)) ∗ d e l t a _ t + sigma . ∗ ( sqrt ( d e l t a _ t ) . ∗ randn ( 1 ) ) ;
i=i +1;
end

2 Least Squares Calibration

The idea of least squares is that we choose parameter estimates that minimize the average squared dierence
between observed and predicted values. That is, we maximize the t of the model to the data by choosing
the model that is closest, on average, to the data. Rewriting (1.5) we have
√
rt+δt = rt (1 − λδt) + λµδt + σ δtN (0, 1) (2.1)

The relationship between consecutive observations rt+δt and rt is linear with a iid normal random term 

rt+δt = art + b +  (2.2)

or  
a
(2.3)
 
[rt+δt ] = rt 1 +
b
where
a = 1 − λδt (2.4)

3
 Least Squares Fitting
1

0.9

0.8

0.7

0.6

ri+1
0.5

0.4

0.3

0.2

0.1
0.4 0.5 0.6 0.7 0.8 0.9 1
ri

Figure 2.1: Least Squares Fitting

b = λµδt (2.5)
√
 = σ δtN (0, 1) (2.6)
where using a least squares tting as described in the Appendix.
n
 n n

1
P P P
(ri+1 ri ) − n ri ri+1
â = i=1 n  ni=1 ni=1  (2.7)
P 2 1 P P
ri − n ri ri
i=1 i=1 i=1
and
n n
!
1 X X
b̂ = ri+1 − ari (2.8)
n i=1 i=1
so we can estimate the model's parameters as

1−a
λ= (2.9)
δt

b
µ= (2.10)
1−a

V ar ()
σ2 = (2.11)
δt
we refer the reader to the Appendix (4) where the derivation and the implemented code is presented.

3 Maximum Likelihood Calibration

In maximum likelihood estimation, we search over all possible sets of parameter values for a specied model
to nd the set of values for which the observed sample was most likely. That is, we nd the set of parameter
values that, given a model, were most likely to have given us the data that we have in hand. The distribution
of rt+δt in the Euler scheme is given by
2
1 [rt+δt − (rt + λ (µ − rt ) δt)]
f (rt+δt | rt , µ, λ, σ) = √ exp − (3.1)
2πσ 2 δt 2σ 2 δt
so the log-likelihood
n
Y Xn
L = ln f (ri | ri−1 , µ, λ, σ) = ln f (ri | ri−1 , µ, λ, σ) (3.2)
i=1 i=1

4
 Random Noise
0.1

0.08

0.06

0.04

0.02

Noise
0

−0.02

−0.04

−0.06

−0.08

1
Time

Figure 3.1: Noise Box Plot

n
" #
2
X 1 [ri − (ri−1 + λ (µ − ri−1 ) δt)]
= ln √ exp − (3.3)
i=1 2πσ 2 δt 2σ 2 δt
n   2
X 1 [ri − (ri−1 + λ (µ − ri−1 ) δt)]
= ln √ − (3.4)
i=1 2πσ 2 δt 2σ 2 δt
h√ i−n X n 2
[ri − (ri−1 + λ (µ − ri−1 ) δt)]
= ln 2πσ 2 δt − 2
(3.5)
i=1
2σ δt
obtaining

n 2
n   X [ri − (ri−1 + λ (µ − ri−1 ) δt)]
L=− ln 2πσ 2 δt − (3.6)
2 i=1
2σ 2 δt

using the procedure described in the Appendix (4) the parameters are estimated giving

n  
1 X ri − ri−1 (1 − λδt)
µ̂ = (3.7)
n i=1 λδt

n
P
(ri − ri−1 ) (µ − ri−1 )
1 i=1
λ= n (3.8)
δt P 2
(µ − ri−1 )
i=1

n 2
1 X [ri − (ri−1 + λ (µ − ri−1 ) δt)]
σ2 = (3.9)
n i=1 δt

we refer the reader to the Appendix (4) where the derivation and the implemented code is presented.
% Y−B( 1 ) .X−B(2)= Noise
de_trended=Y−B ( 1 ) . ∗ X−B( 2 )

figure ( 5 )
h i s t f i t ( de_trended , 1 3 )
figure ( 6 )
p r o b p l o t ( ' normal ' , de_trended )
grid on

5
 Probability plot for Normal distribution

0.995
0.99

0.95
0.9

0.75

Probability
0.5

0.25

0.1
0.05

0.01
0.005

−0.1 −0.05 0 0.05 0.1

Data

Figure 3.2: Normal Plot

4 Long Term Quantile Method

The major assumption in this model is that the quantiles from the historical data are representative for
quantiles in the future. Therefore, a 95% condence interval is taken from the historical data and the
parameters in the short interest rate model are chosen such that in 95% of the cases the generated interest
rates will fall within the condence interval taken from the historical data.

σ2
 
rt ∼ N r0 exp (−λt) + µ (1 − exp (−λt)) , (1 − exp (−2λt)) (4.1)
2λ
The dierence between the long term standard deviation and the deviation at t is
r r
σ2 σ2
∆σ = (1 − exp (−2λt)) − (4.2)
2λ 2λ
using Taylor expansion for exp (−2λt) ≪ 1
r
σ2
  
1
∆σ = 1+ exp (−2λt) − 1 (4.3)
2λ 2
so
∆σ 1
= exp (−2λt)  1 (4.4)
σ 2
the dierence between the long term mean and the mean at t is

∆µ = (r0 − µ) exp (−λt) (4.5)

so  
∆µ r0
= − 1 exp (−λt) (4.6)
µ µ
The condence interval is the one such that
 
σ σ
P µ − 1.96 √ ≤ lim rt ≤ µ + 1.96 √ = 0.95 (4.7)
2λ t→∞ 2λ
calling
σ
q̃0.95 = µ + 1.96 √ (4.8)
2λ
and
σ
q̃0.05 = µ − 1.96 √ (4.9)
2λ

6
we can obtain the parameters as

q̃0.95 + q̃0.05
µ= (4.10)
2
and

2
σ 2 (1.96)
λ=2 2 (4.11)
(q̃0.95 − q̃0.05 )

%% Montecarlo S i m u l a t i o n
j =1;

while j <10000
S(1)= S0 ; % S t a r t i n g p o i n t

while i <n+1
S ( i )= S ( i − 1) + lambda . ∗ ( mu−S ( i − 1)) ∗ d e l t a _ t + sigma . ∗ ( sqrt ( d e l t a _ t ) . ∗ randn ( 1 ) ) ;
i=i +1;
end

MC( j )=S (T ) ;
j=j +1;
end

figure ( 9 )
h= h i s t f i t (MC, 1 2 )
t i t l e ( ' D i s t r i b u t i o n o f t h e V a s i c e k p r o c e s s by Montecarlo S i m u l a t i o n ' ) grid on
the error associated with λ is
2 2
ln λ = ln 2 (1.96) + ln σ 2 − ln (x − y) (4.12)
using partial dierentiation we have that the maximum percentual error (for a detailed discussion on errors
analysis see [1, 2]) is given by
   
4λ  4σ 
 + 2  (4q̃0.95 + 4q̃0.05 )  (4.13)
 
= 2 
λ σ   (q̃0.95 − q̃0.05 ) 

z=q u a n t i l e (MC, 0 . 9 5 ) y=q u a n t i l e (MC, 0 . 0 5 )

% Mu
( z+y ) . / 2
% Lambda
2 ∗ ( 1 . 9 6 . ∗ sigma / ( z−y ) ) . ^ 2

7
 Distribution of the Vasicek process by Montecarlo Simulation
2500

2000

1500

1000

500

0
0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9

Figure 4.1: Montecarlo Histogram

Appendix
Least Squares Fitting
The residuals for the model are given by

Ri = ri+1 − (ari + b) (4.14)

This method minimizes the sum of squared residuals, which is given by
n
X Xn n
X n
X
2 2
S= Ri2 ; = (ri+1 ) + (ari + b) − 2 (ri+1 (ari + b)) (4.15)
i=1 i=1 i=1 i=1

The least squares estimators for the parameters can then be found by dierentiating S with respect to these
parameters and setting these derivatives equal to zero. For b we have that
n n
∂S X X
=2 (ari + b) − 2 ri+1 (4.16)
∂b i=1 i=1

n
X n
X
ari + nb = ri+1 (4.17)
i=1 i=1

n n
!
1 X X
b= ri+1 − ari (4.18)
n i=1 i=1

For a we have that

n n
∂S X X
=2 (ari + b) ri − 2 (ri+1 ri ) (4.19)
∂a i=1 i=1

isolating
n
X n
X
(ari + b) ri = (ri+1 ri ) (4.20)
i=1 i=1

then
n
X n
X Xn
ari2 + bri = (ri+1 ri ) (4.21)
i=1 i=1 i=1

using 4.18 we have !

n n n n n
X 1 X X X X
ari2 + ri ri+1 − ari = (ri+1 ri ) (4.22)
i=1
n i=1 i=1 i=1 i=1

8
is equal to ! !
n n n n n n
X 1 X X 1 X X X
ari2 − a ri ri + ri ri+1 = (ri+1 ri ) (4.23)
i=1
n i=1 i=1
n i=1 i=1 i=1

grouping " !# !
n n n n n n
X 1 X X X 1 X X
a ri2 − ri ri = (ri+1 ri ) − ri ri+1 (4.24)
i=1
n i=1 i=1 i=1
n i=1 i=1

nally
n
 n n

1
P P P
(ri+1 ri ) − n ri ri+1
i=1 i=1 i=1
a= n
 n n
 (4.25)
1
ri2
P P P
− n ri ri
i=1 i=1 i=1

%% C a l i b r a t i o n u s i n g Least Squares r e g r e s s i o n

% P l o t S_i vs S_i −1
figure ( 3 )
Y= S ( 2 : end ) ; % removing f i r s t p o i n t
X=S ( 1 : end − 1);% removing t h e l a s t p o i n t
plot (Y, X, ' . ' )
l s l i n e % l e a s t squares l i n e
ylabel ( ' r_i+1 ' ) xlabel ( ' r_i ' )
t i t l e ( ' Least Squares F i t t i n g ' )
grid on

%Rewrite t h e o f f s e t term
o f f s e t=o n e s ( s i z e ( S , 2 ) − 1 , 1 ) ;
new_X=[X' , o f f s e t ] ;
B = new_X\Y' % S o l v e s in t h e Least Squares s e n s e

est_lambda=1−B ( 1 ) . / d e l t a _ t
est_mu= B( 2 ) . / ( 1 − B( 1 ) )

Maximum Likelihood Fitting

An Estimator for µ
n
∂L X 2 [ri − (ri−1 + λ (µ − ri−1 ) δt)] λµδt
= − (4.26)
∂µ i=1 2σ 2 δt
n
X [ri − (ri−1 (1 − λδt) + λµδt)] λµ
= − =0 (4.27)
i=1
σ2
n
X
[ri − (ri−1 (1 − λδt) + λµδt)] = 0 (4.28)
i=1
n
X n
X
λµδt = [ri − ri−1 (1 − λδt)] (4.29)
i=1 i=1
n  
1 X ri − ri−1 (1 − λδt)
µ̂ = (4.30)
n i=1 λδt
An Estimator for

9
 n
∂L X [ri − (ri−1 + λ (µ − ri−1 ) δt)] (µ − ri−1 ) δt
= − (4.31)
∂λ i=1 σ 2 δt
n
X [ri − (ri−1 + λ (µ − ri−1 ) δt)] (µ − ri−1 )
= − (4.32)
i=1
σ2
h i
2
Xn (ri − ri−1 ) (µ − ri−1 ) − λ (µ − ri−1 ) δt
= − =0 (4.33)
i=1
σ2
n
X Xn
2
(ri − ri−1 ) (µ − ri−1 ) = λ (µ − ri−1 ) δt (4.34)
i=1 i=1
n
P
(ri − ri−1 ) (µ − ri−1 )
1 i=1
λ= n (4.35)
δt P 2
(µ − ri−1 )
i=1

An estimator for σ
n 2
∂L n X [ri − (ri−1 + λ (µ − ri−1 ) δt)]
=− 2
(4πσδt) − −2 (4.36)
∂σ 4πσ δt i=1
2σ 3 δt
n 2
n X [ri − (ri−1 + λ (µ − ri−1 ) δt)]
=− + =0 (4.37)
σ i=1 σ 3 δt
n 2
n X [ri − (ri−1 + λ (µ − ri−1 ) δt)]
= (4.38)
σ i=1 σ 3 δt
n 2
1 X [ri − (ri−1 + λ (µ − ri−1 ) δt)]
σ2 = (4.39)
n i=1 δt

%% C a l i b r a t i o n u s i n g Maximum L i k e l i h o o d E s t i m a t o r s
n = length ( S) − 1;
Sx = sum ( S ( 1 : end − 1) ) ;
Sy = sum ( S ( 2 : end ) ) ;
Sxx = sum ( S ( 1 : end − 1).^2 ) ;
Sxy = sum ( S ( 1 : end − 1). ∗ S ( 2 : end ) ) ;
Syy = sum ( S ( 2 : end ) . ^ 2 ) ;

mu = ( Sy ∗ Sxx − Sx ∗ Sxy ) / ( n ∗ ( Sxx − Sxy ) − ( Sx^2 − Sx ∗ Sy ) ) ;

lambda = ( ( Sxy − mu∗ Sx − mu∗ Sy + n ∗mu^2) / ( Sxx −2∗mu∗ Sx + n ∗mu^2) ) / d e l t a _ t ;
a = 1− lambda ∗ d e l t a _ t ;
sigmah2 = ( Syy − 2 ∗ a ∗ Sxy + a^2 ∗ Sxx − 2 ∗mu∗ (1 − a ) ∗ ( Sy − a ∗ Sx ) + n ∗mu^2 ∗ (1 − a ) ^ 2 ) / n ;
sigma = sqrt ( sigmah2 ∗ 2 ∗ lambda/(1 − a ^ 2 ) ) ;
end

Acknowledgements
This work has been done under the collaboration of Aylin Chakaroglu, MSc. (Societé General RISQ/-
MAR/RIM team) who we wish to thank for the help in the revision of this report.

10
References
[1] Philip R Bevington and D Keith Robinson. Data reduction and error analysis. McGraw-Hill, 2003.

[2] William Lichten. Data and error analysis in the introductory physics laboratory. Allyn & Bacon, 1988.

[3] George E Uhlenbeck and Leonard S Ornstein. On the theory of the brownian motion. Physical review,
36(5):823, 1930.

[4] Emile Van Elen. Term structure forecasting, 2010.

[5] Oldrich Vasicek. An equilibrium characterization of the term structure. Journal of nancial economics,
5(2):177188, 1977.

[6] S Zeytun and A Gupta. A Comparative Study of the Vasicek and the CIR Model of the Short Rate. ITWM
Kaiserslautern, Germany, 2007.

11