Numerical Analysis of Yield Curves Implied by Two-Factor Interest Rate Models

Numerical Analysis of Yield Curves Implied by Two-Factor Interest Rate Models

VERONIKA CHRONHOLM

Department of Mathematical Sciences

C

HALMERS

U

NIVERSITY OF

T

ECHNOLOGY

Thesis for the degree of Master of Science

Numerical Analysis of Yield Curves Implied by Two-Factor Interest Rate

Models

VERONIKA CHRONHOLM

Department of Mathematical Sciences Chalmers University of Technology

University of Gothenburg

Gothenburg, Sweden 2021

© Veronika Chronholm, 2021.

Supervisor: Simone Calogero, Department of Mathematical Sciences Examiner: Moritz Schauer, Department of Mathematical Sciences

Department of Mathematical Sciences

Chalmers University of Technology and University of Gothenburg SE-412 96 Gothenburg

Telephone +46 (0)31- 772 10 00

Typeset in L

^A

TEX

Gothenburg, Sweden 2021

Numerical Analysis of Yield Curves Implied by Two-Factor Interest Rate Models Veronika Chronholm

Department of Mathematical Sciences

Chalmers University of Technology and University of Gothenburg

Abstract

We investigate the yield curves implied by coupon bonds in models where the market short

rate is given by a two-factor stochastic model. Specifically, we investigate generalisations

of the two-factor Vasicek, Cox-Ingersoll-Ross, and mixed models where the two Brownian

motions that feature in each model are allowed to have nonzero constant correlation. We

also study the two-factor Rendlemann-Bartter model with nonzero constant correlation. In

all these models, we manage to recreate the four yield curve shapes commonly discussed in

the literature; normal, steep, inverted, and flat. We also investigate how some of the interest

rate model parameters affect the qualitative properties of the yield curves produced, and

compare the yield curves obtained in the different models.

Acknowledgements

I would like to thank my supervisor Simone Calogero for his support and guidance through- out this project.

Veronika Chronholm, Gothenburg, May 2021

Contents

1 Introduction 1

1.1 Outline . . . . 2

1.2 Method and Numerics . . . . 2

2 Background 3 2.1 Zero Coupon Bonds and Coupon Bonds . . . . 3

2.2 Risk Neutral Pricing . . . . 3

2.3 Yield to Maturity and Yield Curves . . . . 4

3 Interest Rate Models 7 3.1 Theory . . . . 7

3.1.1 A Brief Note on the Existence and Uniqueness of Solutions to Multi- dimensional SDE . . . . 7

3.1.2 Introduction of Interest Rate Models . . . . 7

3.1.3 Mean Reversion of the Interest Rate . . . . 11

3.1.4 Meaning of Parameter Values . . . . 14

3.2 Numerical Methods . . . . 14

3.2.1 Euler-Maruyama Finite Difference Schemes . . . . 14

3.2.2 Modelling SDE with Correlated Brownian Motions . . . . 17

3.2.3 Numerical correlation of the factors X1 and X2 . . . . 19

4 Zero Coupon Bond Pricing 21 4.1 Theory . . . . 21

4.1.1 The Pricing PDE . . . . 21

4.1.2 Affine yield models . . . . 23

4.2 Numerical Methods . . . . 27

4.2.1 Monte-Carlo Methods . . . . 28

4.2.2 Numerical Comparison of Zero Coupon Bond Prices . . . . 29

4.2.3 Comparison of ZCB Prices with Different Correlation . . . . 33

4.2.4 Comparison of ZCB Prices in Different Interest Rate Models . . . . . 34

5 Yield curves 35 5.1 Theory . . . . 35

5.1.1 Interpretation of yield curve shapes . . . . 35

5.1.2 Calculating the yield to maturity . . . . 36

5.1.3 Choice of Coupon Values . . . . 37

5.2 Numerical Results . . . . 37

5.2.1 Recreation of Yield Curve Shapes in Different Interest Rate Models . 38 5.2.2 Comparison of Yield Curves In Different Models . . . . 38

5.2.3 Comparison of Yield Curves with Different Correlation . . . . 39

5.2.4 Recreating Normal Yield Curves with Negative Yield for Short Matu-

rities . . . . 40

5.2.5 Mixed Model Comparison . . . . 41

5.2.6 Rendleman-Bartter Model Without Mean Reversion . . . . 42

6 Summary and Discussion 45 Bibliography 47 A Sample Code 49 A.1 Vasicek Multi-Level Monte-Carlo Code . . . . 49

A.1.1 Functions . . . . 50

A.2 CIR Model Code . . . . 54

A.2.1 Functions . . . . 56

B Convergence Plots 59

1

Introduction

The yield curve associated with a class of bonds tells us about the cost of borrowing and lending money by investing in those bonds. Historically, changes in the shape of the yield curve have happened before a significant economic change [6] [17], and it is thus of interest to be able to recreate the different shapes that the yield curve can take in a given numerical model. In the literature, the four different shapes are referred to as normal, inverted, steep, and flat or humped.

A yield curve plots the yield to maturity of bonds with varying maturities versus their maturity. The yield to maturity of a bond is the constant continuously compounded interest rate at which all future payments of the bond should be discounted to obtain its present value.

Institutions such as the Swedish National Debt Office and the United States Treasury issue coupon bonds as a way to finance public spending. It is thus of interest to specifically investigate coupon bond yield curves. A coupon bond is a financial contract which promises to pay some amount of money – the face value – at maturity, and some other amounts of money – the coupons – at other times between the present date and the time of maturity.

The aim of this thesis is to numerically investigate the coupon bond yield curves implied by two-factor interest rate models for the short rate. The short rate is defined as the continuously compounded, annualized interest rate at which money can be borrowed instan- taneously and for an infinitesimally short time. The interest rate models considered are thus models in which the short rate is an affine function of two stochastic processes, which solve a system of two coupled stochastic differential equations. A benefit of two-factor models is that they allow for two different sources of uncertainty in the model, in the form of two Brownian motions, which can be either independent or correlated.

When modelling the interest rate, it might be of interest whether the chosen model can recreate the different qualitative behaviours of the yield curve. We thus attempt to recreate the main shapes that the yield curve can take in the different interest rate models that we study. We also investigate how allowing for correlation between the Brownian motions in the interest rate model affects the bond prices as well as the behaviour of the yield curve in the models.

The focus lies on the class of models known as two-factor affine interest rate models, namely

the two-factor Vasicek, Cox-Ingersoll-Ross, and mixed models, where we allow for nonzero

constant correlation between the two driving Brownian motions in each model. We also

consider one non-affine model, the two-factor Rendlemann-Bartter model with constant

correlation between the driving Brownian motions.

1.1 Outline

In chapter 2 we present some relevant background on financial mathematics. The concepts of zero coupon bonds and coupon bonds and the relationship between their prices are covered.

We also discuss risk neutral pricing and the risk neutral pricing formula, and how to calculate the yield to maturity of a coupon bond.

Chapters 3, 4, and 5 are each divided into two parts; one theoretical and one that covers nu- merical methods. In the theoretical part of each chapter, theoretical concepts and properties are introduced, whereas relevant parts of the numerical methods and results are presented in the numerical methods part.

In chapter 3, we introduce the interest rate models that we study and discuss some of their theoretical properties, such as mean reversion, and existence and uniqueness of solu- tions for the Cox-Ingersoll-Ross and mixed models, where at least one factor in the model must remain posititive for all times t. We also discuss numerical methods for solving the stochastic differential equations of the models and thus generating paths of the interest rate.

Specifically, we discuss the Euler-Maruyama scheme and some modifications of it for the Cox-Ingersoll-Ross model and the mixed model.

In chapter 4 we discuss the pricing of zero coupon bonds using the risk neutral pricing formula. We derive the partial differential equation for the price of a zero coupon bond, and discuss the case of affine yield models where the PDE can be separated into a system of ordinary differential equations. We also discuss how to use Monte-Carlo methods to numerically estimate the zero coupon bond price given by the risk neutral pricing formula, and compare the prices obtained using this method with the prices obtained by numerically solving the system of ODE that results from the pricing PDE, where possible.

In chapter 5, we discuss in more detail how to calculate the yield to maturity of a coupon bond, and how we choose the values of the coupons. In the numerical part of the chapter, our main numerical results are presented, in the form of numerically generated yield curves from the different interest rate models.

In chapter 6, we summarise the numerical results and present a brief discussion.

1.2 Method and Numerics

The main numerical methods employed in this thesis are Monte-Carlo and multi-level Monte-

Carlo methods, in combination with Euler-Maruyama schemes to generate paths of the

solutions to the stochastic differential equations that make up the interest rate models. To

solve the system of ODE that the zero coupon bond pricing PDE separates into, the MatLab

solver ode45 has been used. All code is written in MatLab version 9.7 for Windows. Some

sample code displaying a few of the different algorithms used can be found in the appendix,

as can a few numerical convergence plots.

2

Background

2.1 Zero Coupon Bonds and Coupon Bonds

A zero coupon bond (ZCB) with face value K and maturity T is a contract that promises to pay the amount K at time T . Since holding a portfolio containing one share of a zero coupon bond with face value K and maturity T , is equivalent holding a portfolio with K shares of a zero coupon bond with face value 1 and maturity T , we can without loss of generality limit our discussion to bonds with face value 1.

A coupon bond with face value K, maturity T , and coupons {c

₁

, ..., c

_N

} paid at times {t

1

, ..., t

_N

}, where 0 < t

1

< ... < t

_N

= T , is a contract which promises to pay the amount c

_k

at time t

_k

for k = 1,2,...,(N − 1) and the amount (K + c

_N

) at maturity. Bonds with maturity 1 year or shorter do not pay any coupons, so a coupon bond with maturity 1 year simply pays the amount K at maturity, thus being equivalent to a zero coupon bond with maturity T and face value K.

Holding a coupon bond at time t gives the holder no right to coupons paid out before time t. As in the case of zero coupon bonds, it is enough to consider coupon bonds with face value 1.

For a time t ∈ [0,T ], holding a coupon bond which pays coupons as described above is equivalent to holding a portfolio consisting of shares of zero coupon bonds with face value 1 and with maturities {t

k

, k : t

k

> t}, where the number of shares of the coupon bond with maturity t

N

= T is (1 + c

N

) and the number of shares of the zero coupon bond with maturity t

k

is c

k

for all other k. In other words, if we let k(t) be the smallest index such that t

k

> t, then the value B

_c

(t,T ) of the coupon bond at time t is

B

c

(t, T ) =

N −1

X

i=k(t)

c

i

B(t,t

i

) + (1 + c

N

)B(t,T ) (2.1)

where B(t,t

k

) is the price of the zero coupon bond with face value 1 and maturity t

k

at time t. Thus, to calculate the price of a coupon bond at time t in a given market model, it is enough to be able to calculate the values at time t of zero coupon bonds with different maturities.

2.2 Risk Neutral Pricing

A sufficient condition for a market model to be arbitrage-free is the existence of a risk-neutral probability measure under which the discounted value of any portfolio in the market is a martingale.[16]

The risk-neutral pricing formula can be derived by considering a portfolio consisting of

one share of a given asset. The condition that the discounted value of the portfolio is a martingale, is then equivalent to requiring that the discounted asset price be a martingale in the risk-neutral probability measure. The discount factor D(t) solves the differential equation dD(t) = −R(t)D(t)dt, where R(t) is the interest rate, or short rate. The solution to this differential equation is

D(t) = exp 

− Z

t

0

R(s)ds 

. (2.2)

By requiring that the discounted price of a zero coupon bond at time t, D(t)B(t,T ) be a martingale in the risk-neutral probability measure we get

D(t)B(t,T ) = ˜ E[D(T )B(T ,T ) | F (t)], (2.3) where F (t) is the filtration generated by the Brownian motion. Since we have restricted the discussion to zero coupon bonds with face value 1, by definition B(T,T ) = 1. Further, D(t) is F (t)-measurable, meaning that we can divide both sides of the equation by D(t) and bring it inside the expectation on the right hand side. Thus, we obtain

B(t,T ) = ˜ E

 exp 

− Z

T

t

R(s)ds 

| F (t)



. (2.4)

In particular, the price at time zero of the zero coupon bond is

B(0,T ) = ˜ E

 exp 

− Z

T

0

R(s)ds  

. (2.5)

In the above, the expectation is taken in the risk-neutral probability measure.

To go from the physical probability measure to the risk-neutral probability measure, one needs to determine the market price of risk. To do this, there needs to be a risky asset – such as a stock – in the market from which one can derive the so called market price of risk equations. Since we are working with a model where the dynamics of the short rate are specified instead, we don’t have a risky asset given in the model from which we could get the market price of risk equations. To get around this problem, any interest rate model discussed in this thesis will be formulated directly in the risk-neutral probability measure.

2.3 Yield to Maturity and Yield Curves

In the case of a zero coupon bond with face value 1 and maturity T , the yield to maturity at a time 0 ≤ t ≤ T is the quantity Y (t,T ) for which it holds that

B(t,T ) = e

−Y (t,T )(T −t)

, (2.6) where B(t,T ) is the price of the zero coupon bond at time t.

The yield to maturity at time t of a coupon bond with maturity T is the quantity Y

c

(t,T )

for which

2. Background

B

_c

(t,T ) =

N −1

X

i=k(t)

c

_i

e

^{−Yc(t,T )(ti−t)}

+ (1 + c

_N

)e

^−Yc(t,T )(T −t)

. (2.7)

Thus, if we know the price of the coupon bond at time t, we have an equation where the only unknown is Y

c

(t,T ). Since B

c

(t,T ) can be calculated from the prices of zero coupon bonds through equation (2.1), to compute the yield to maturity of a coupon bond in a given model, we simply need to be able to price zero coupon bonds and then solve equation (2.7) for the yield.

The yield to maturity at time zero, Y

_c

(0,T ) is given by

B

c

(0,T ) =

N −1

X

i=1

c

i

e

^{−Yc(0,T )ti}

+ (1 + c

N

)e

^{−Yc(0,T )T}

. (2.8)

If we measure time in years, and let both the maturity and all times at which a coupon is paid be integer numbers of years, equation (5.1) is simply a polynomial equation for x = e

^{−Yc(0,T )}

, namely

N −1

X

k=1

c

k

x

^k

+ (1 + c

N

)x

^N

− B

c

(0,T ) = 0. (2.9)

This polynomial equation will be discussed in more detail in chapter 5.

3

Interest Rate Models

3.1 Theory

3.1.1 A Brief Note on the Existence and Uniqueness of Solutions to Multidimensional SDE

Consider a two-dimensional stochastic differential equation

dX(t) = b(X(t), t)dt + σ(X(t),t)dW(t), (3.1) where σ is a 2 × 2 matrix, b is a two-dimensional vector, and W(t) is a two-dimensional Brownian motion. The SDE then has a unique strong global solution if the Lipschitz condi- tion

|b(x,t) − b(y,t)| + |σ(x,t) − σ(y,t)| < K

T

|x − y| (3.2)

and the linear growth condition

|b(x,t)| + |σ(x,t)| ≤ K

T

(1 + |x|) (3.3) both hold for all t such that 0 ≤ t ≤ T , and for all T, and where the constant K

T

only depends on T [13]. In the above, the norm of the vector b is defined as |b| = pb

²₁

+ b

²₂

, and the norm of the matrix σ is defined by |σ|

²

= Tr(σσ

^>

). We note that the conditions (3.2) and (3.3) are sufficient but not all necessary for the existence of a unique solution.

3.1.2 Introduction of Interest Rate Models

The studied interest rate models are two-factor interest rate models, where the market short rate is modelled by a stochastic process R(t) that is a function of two factors X

₁

(t) and X

₂

(t) which solve some stochastic differential equations. Specifically, we study models of the form

R(t) = δ

₀

+ δ

1

X

₁

(t) + δ

2

X

₂

(t) (3.4a)

dX

1

(t) = (µ

1

− λ

11

X

1

(t) − λ

12

X

2

(t))dt + σ

1

X

1

(t)

^γ1

dW

1

(t) (3.4b)

dX

₂

(t) = (µ

₂

− λ

₂₁

X

₁

(t) − λ

₂₂

X

₂

(t))dt + σ

₂

X

₂

(t)

^γ2

dW

₂

(t). (3.4c)

W

₁

(t) and W

₂

(t) are standard Brownian motions, i.e. continuous-time stochastic processes which have the following properties[13].

• W

i

(0) = 0

• Independent increments: For t > s, W

_i

(t) − W

_i

(s) is independent of the values taken up to time s, or of W

i

(u) for 0 ≤ u < s.

• Normal increments: W

i

(t)−W

i

(s) is normally distributed with mean zero and variance t − s.

• Continuity of paths: The stochastic process W

i

(t) for t ≥ 0 has almost surely contin- uous paths.

In this thesis, we consider the case where the Brownian motions W

₁

(t) and W

₂

(t) have constant correlation ρ ∈ [−1,1], and thus dW

₁

dW

₂

= ρdt. The parameters d

₀

, d

₁

, d

₂

, µ

₁

, µ

₂

, λ

₁₁

, λ

₁₂

, λ

₂₁

, λ

₂₂

, σ

₁

, and σ

₂

are real constants, and γ

₁

, γ

₂

∈ [0,1]. Here, we study the two-factor Vasicek, Cox-Ingersoll-Ross (CIR), and mixed models, as well as the two-factor Rendleman-Bartter model.

The two factor Vasicek model

In the two-factor Vasicek model, γ

1

= γ

2

= 0 and the model is thus given by

R(t) = δ

0

+ δ

1

X

1

(t) + δ

2

X

2

(t) (3.5a) dX

₁

(t) = (µ

₁

− λ

₁₁

X

₁

(t) − λ

₁₂

X

₂

(t))dt + σ

₁

dW

₁

(t) (3.5b) dX

2

(t) = (µ

2

− λ

21

X

1

(t) − λ

22

X

2

(t))dt + σ

2

dW

2

(t). (3.5c) In the Vasicek model, both the factors X

1

(t) and X

2

(t) and the interest rate R(t) are allowed to take both positive and negative values. The conditions (3.2) and (3.3) hold for all x, y – and we note that the drift and diffusion are time imdependent – so the stochastic differential equations (3.5b)-(3.5c) have a unique strong solution.

The two factor CIR model

In the two-factor CIR model, γ

₁

= γ

₂

= 0.5 and the model is given by

R(t) = δ

₀

+ δ

₁

X

₁

(t) + δ

₂

X

₂

(t) (3.6a) dX

1

(t) = (µ

1

− λ

11

X

1

(t) − λ

12

X

2

(t))dt + σ

1

p X

1

dW

1

(t) (3.6b) dX

2

(t) = (µ

2

− λ

21

X

1

(t) − λ

22

X

2

(t))dt + σ

2

p X

2

dW

2

(t). (3.6c)

Since both factors {X

₁

(t)} and {X

₂

(t)} feature under a square root in the case of the CIR model, both factors must remain positive for the interest rate to remain real.

Furthermore, the Lipschitz condition (3.2) holds everywhere except at zero, where the deriva- tive of the square root function is infinite. Thus, if we under some condition on the coeffi- cients of the SDE can be sure that the process does not hit zero, we can say that the SDE has a unique strong solution.

In the case of the one-dimensional CIR model,

dX(t) = a(b − X(t))dt + σ p

X(t)dW (t), (3.7)

3. Interest Rate Models

X(t) > 0 for all t almost surely when the so-called Feller condition ab >

^σ₂²

holds, as long as the initial condition X(0) > 0, and where the constants a and b are assumed to be positive.

Thus, under this condition, the one-dimensional CIR SDE has a unique strong solution.

In the case of the two-factor CIR model, we first consider the case when σ

1

= σ

2

= 0, when the system of SDE reduces to a system of ODE to get a necessary condition for positivity on the drift coefficients. We then get the system of ODE

dx

1

(t)

dt = µ

1

− λ

11

x

1

(t) − λ

12

x

2

(t) (3.8a) dx

₂

(t)

dt = µ

2

− λ

21

x

1

(t) − λ

22

x

2

(t). (3.8b) For the function x

1

(t) to not become negative, we need its derivative to be nonnegative when x

1

(t) = 0, and similarly the derivative of x

2

(t) must be nonnegative when x

2

(t) = 0, meaning

µ

₁

− λ

12

x

₂

≥ 0 (3.9)

µ

2

− λ

21

x

1

≥ 0 (3.10)

This holds when µ

₁

≥ 0 and λ

₁₂

≤ 0, as long as x

₂

(t) ≥ 0. Similarly, for x

₂

(t) not to become negative, we require that µ

₂

≥ 0 and λ

₂₁

≤ 0.

Now, in the case where σ

1

, σ

2

6= 0, these conditions on the drift coefficients are not sufficient for the stochastic processes X

1

(t) and X

2

(t) to remain (strictly) positive. Similarly to in the one-factor CIR model, we also require a condition which involves the relationship between the drift parameters and the volatility parameters σ

1

and σ

2

. A Feller condition for two- dimensional CIR processes is not widely studied in the literature, but a condition that can be applied to the two-factor CIR model with non-correlated Brownian motions is presented and proved in [7]. The condition is as follows. Consider a system of SDE of the form

dX(t) = (µ − ΛX(t))dt + Σ pν

₁

(X(t)) 0 0 pν

2

(X(t))



dW(t), (3.11a)

where W(t) is a standard two-dimensional Brownian motion with zero correlation between the components, Σ and Λ are 2×2 constant matrices and µ is a constant vector. Furthermore,

ν

i

(x) = a

i

+ b

i

· x. (3.12)

The conditions for maintaining strict positivity of the two-dimensional process X(t) are then formulated as follows; for i = 1,2

For all x such that ν

i

(x) = 0, b

^>_i

(µ − Λx) > 1

2 b

^>_i

ΣΣ

^>

b

i

(3.13) For all j, if (b

^>_i

Σ)

_j

6= 0, then ν

_i

= ν

_j

(3.14)

In the case of the two-factor CIR model (3.6), µ = µ

₁

µ

₂



, Λ = λ

₁₁

λ

₁₂

λ

₂₁

λ

₂₂



, Σ = σ

₁

0 0 σ

₂

 , and ν

i

(x) = x

i

, i.e. a

1

= a

2

= 0, b

1

= 1

0



, and b

2

= 0 1



. The condition (3.13) then

becomes

∀x : x

1

= 0, 1 0

"

µ

1

µ

2



− λ

11

λ

12

λ

21

λ

22

  0 x

2

 #

> 1

2 1 0
σ

₁²

0 0 σ

²₂

 1 0



(3.15)

∀x : x

2

= 0, 0 1

"

µ

₁

µ

2



− λ

₁₁

λ

12

λ

21

λ

22

 x

₁

0

 #

> 1

2 0 1
σ

₁²

0 0 σ

²₂

 0 1



(3.16)

or equivalently

µ

₁

− λ

₁₂

x

₂

> 1

2 σ

²₁

(3.17a)

µ

2

− λ

21

x

1

> 1

2 σ

²₂

. (3.17b)

We note that these conditions are a stricter version of the conditions (3.9). Furthermore, if we as previously impose λ

₁₂

, λ

₂₁

≤ 0, we can simply require that µ

₁

>

¹₂

σ

₁²

and µ

₂

>

¹₂

σ

²₂

. This is enough as long as x

₁

,x

₂

≥ 0.

The condition (3.14) becomes trivial in the case of the two-factor CIR model (3.6), as we shall now see. For i = 1,

1 0
σ

1

0 0 σ

2



= σ

1

0
. (3.18)

Only the first component is nonzero, meaning that it must hold that ν

1

= ν

1

, which is clearly true. The same argument can be made for i = 2.

Thus a sufficient condition for the processes X

1

(t) and X

2

(t) to remain strictly positive for all t almost surely, is to require that the conditions (3.17) hold or alternatively, as we have chosen to do here, that λ

12

, λ

21

≤ 0 and µ

1

>

¹₂

σ

²₁

and µ

2

>

¹₂

σ

₂²

. Under these conditions, the two-dimensional SDE (3.6) has a unique strong solution[7]. This is the case since the processes X

1

(t) and X

2

(t) never hit zero, which is the only point where the Lipschitz condition does not hold.

Finally, in the two-factor CIR model (3.6), the parameters δ

0

, δ

1

, and δ

2

are restricted to δ

0

≥ 0, and δ

1

,δ

2

> 0, so that the interest rate remains positive.

The two factor mixed model

In the two-factor mixed model, one factor follows a CIR model and one factor follows a Vasicek model, so for example γ

₁

= 0.5 and γ

₂

= 0 so that the model is given by

R(t) = δ

0

+ δ

1

X

1

(t) + δ

2

X

2

(t) (3.19a) dX

1

(t) = (µ

1

− λ

11

X

1

(t) − λ

12

X

2

(t))dt + σ

1

p X

1

(t)dW

1

(t) (3.19b)

dX

2

(t) = (µ

2

− λ

21

X

1

(t) − λ

22

X

2

(t))dt + σ

2

dW

2

(t). (3.19c)

In the mixed model, one factor – here the factor X

₁

(t) – must remain positive, whereas the

other factor – here X

₂

(t) – is allowed to become negative. To ensure that the factor X

₁

(t)

remains positive in theory, we for simplicity set the parameter λ

₁₂

= 0, and get the model

3. Interest Rate Models

R(t) = δ

₀

+ δ

₁

X

₁

(t) + δ

₂

X

₂

(t) (3.20a) dX

₁

(t) = (µ

₁

− λ

₁₁

X

₁

(t))dt + σ

₁

p

X

₁

(t)dW

₁

(t) (3.20b) dX

2

(t) = (µ

2

− λ

21

X

1

(t) − λ

22

X

2

(t))dt + σ

2

dW

2

(t). (3.20c) We can then apply the Feller condition for the one factor CIR model to equation (3.20b) and note that X

₁

(t) > 0 for all t almost surely as long as µ

₁

>

^σ₂²¹

, and X(0) > 0. Under this condition, the SDE (3.20b) has a unique strong solution. Consequently, we can say that the system of SDE (3.20) has a unique strong solution.

The two factor Rendleman-Bartter model

Finally, we also consider the two-factor Rendleman-Bartter model, where γ

1

= γ

2

= 1, and the model is given by

R(t) = δ

0

+ δ

1

X

1

(t) + δ

2

X

2

(t) (3.21a) dX

₁

(t) = (µ

₁

− λ

₁₁

X

₁

(t) − λ

₁₂

X

₂

(t))dt + σ

₁

X

₁

(t)dW

₁

(t) (3.21b) dX

2

(t) = (µ

2

− λ

21

X

1

(t) − λ

22

X

2

(t))dt + σ

2

X

2

(t)dW

2

(t). (3.21c) The drift and volatility fulfil the conditions (3.2) and (3.3), and thus a unique strong solution exists.

What is commonly referred to as the Rendleman-Bartter interest rate model, is the one- factor model

dR(t) = θR(t)dt + γR(t)dW (t), (3.22) where the interest rate follows a geometric Brownian motion[11]. In this model the interest rate follows the same dynamics as a stock does in what is commonly called the Black-Scholes model. One problem with using the Black-Scholes model for interest rate dynamics, is that the process R(t) is then not mean-reverting.

Here, we will mainly consider a modified version of the two-factor Rendleman-Bartter model which has the mean-reverting property, but we will also briefly consider the non-mean- reverting two-factor version of the model for comparison.

In the next section, we find the conditions on the coefficients in the general interest rate model (3.4) under which the interest rate is mean-reverting. In particular, under these conditions, the two-factor Rendleman-Bartter model discussed here is naturally also mean- reverting.

3.1.3 Mean Reversion of the Interest Rate

An important property of any interest rate model is that it replicates the mean-reverting property observed in interest rates in real markets. This means that in the limit as t → ∞, the mean of the interest rate converges to a constant, finite value. We will now derive the conditions under which the two factor Vasicek model, CIR model, mixed model, and Rendleman-Bartter model are mean-reverting.

For the factors X

1

(t) and X

2

(t) to be mean reverting, the matrix of coefficients

Λ = λ

₁₁

λ

₂₁

λ

₁₂

λ

₂₂



(3.23)

must have strictly positive eigenvalues[16]. We will now see why this is the case, and find the value that the expectation converges to. The stochastic differential equations (3.4b)-(3.4c) written in integral form are

X

1

(t) = X

1

(0) + Z

t

0

(µ

1

− λ

11

X

1

(s) − λ

12

X

2

(s))ds + Z

t

0

σ

1

X

1

(s)

^γ1

dW

1

(s) (3.24a)

X

2

(t) = X

2

(0) + Z

t

0

(µ

2

− λ

21

X

1

(s) − λ

22

X

2

(s))ds + Z

t

0

σ

2

X

2

(s)

^γ2

dW

2

(s). (3.24b)

Taking the expectation of both equations and using the linearity of expectation, we get

E[X

1

(t)] = X

₁

(0) + E

 Z

t 0

(µ

₁

− λ

₁₁

X

₁

(s) − λ

₁₂

X

₂

(s))ds

 + E

 Z

t 0

σ

₁

X

₁

(s)

^γ1

dW

₁

(s)

 (3.25a) E[X

2

(t)] = X

2

(0) + E

 Z

t 0

(µ

2

− λ

21

X

1

(s) − λ

22

X

2

(s))ds

 + E

 Z

t 0

σ

2

X

2

(s)

^γ2

dW

2

(s)

 . (3.25b) We note that the expectation of the Itô integrals is zero, by the martingale property of the Itô integral, and then exchange the order of integration and expectation in the remaining terms to get

E[X

1

(t)] = X

1

(0) + Z

t

0

(µ

1

− λ

11

E[X

1

(s)] − λ

12

E[X

2

(s)])ds (3.26a)

E[X

2

(t)] = X

₂

(0) + Z

t

0

(µ

₂

− λ

₂₁

E[X

1

(s)] − λ

₂₂

E[X

2

(s)])ds, (3.26b)

having also used the linearity of expectation again. Noting that the expectations of X

1

(t) and X

2

(t) will both be deterministic functions of t, which we can label E[X

¹

(t)] = v

1

(t) and E[X

2

(t)] = v

2

(t) we get the coupled ordinary differential equations

d

dt v

1

(t) = µ

1

− λ

11

v

1

(t) − λ

12

v

2

(t) (3.27a) d

dt v

2

(t) = µ

2

− λ

21

v

1

(t) − λ

22

v

2

(t) (3.27b) or in vector form

d

dt v(t) + Λv(t) = µ (3.28)

where v(t) = v

₁

(t) v

2

(t)



, µ = µ

₁

µ

2



, and Λ is the matrix of coefficients defined in equation

(3.23). The solution to (3.28) is given by the general solution to the homogeneous problem,

3. Interest Rate Models

plus a particular solution. Since the right hand side is a constant vector, a particular solution that is also a constant vector can easily be found using the method of undetermined coefficients. Using the ansatz v

^p

(t) = a

₁

a

₂



, where a

₁

and a

₂

are real constants, and noting that clearly the time derivative of a constant vector is zero, we find

Λv

^p

= µ =⇒ v

^p

= Λ

⁻¹

µ, (3.29)

where Λ

⁻¹

denotes the matrix inverse of Λ.

The general homogeneous solution is found by finding the eigenvalues and eigenvectors of the matrix Λ. Denoting the eigenvalues by r

1

and r

2

, and the eigenvectors by ξ

1

and ξ

2

, the general homogeneous solution is given by

v

^h

(t) = c

₁

ξ

1

e

^−r1t

+ c

₂

ξ

2

e

^−r2t

, (3.30) where c

1

and c

2

are real constants determined by the initial condition. Thus, the general solution to (3.28) is given by

v(t) = Λ

⁻¹

µ + c

1

ξ

1

e

^−r1t

+ c

2

ξ

2

e

^−r2t

. (3.31) Considering the limit as t → ∞, we see that as long as r

1

and r

2

are real, the limit exists if and only if r

1

and r

2

are strictly positive, and in that case

t→∞

lim v(t) = Λ

⁻¹

µ. (3.32)

We recall that r

₁

and r

₂

are the eigenvalues of Λ, and thus we have shown that the general interest rate model (3.4) is mean reverting if and only if the matrix Λ has strictly positive eigenvalues. In that case,

t→∞

lim E[X(t)] = Λ

⁻¹

µ, (3.33)

where X(t) = X

1

(t) X

2

(t)



. We note also that the inverse of the matrix Λ always exists if the matrix has positive eigenvalues. Since the relationship between the short rate R(t) and the two factors X

1

(t) and X

2

(t) is given by (3.4a), we have

t→∞

lim E[R(t)] = δ

0

+ δ

^>

Λ

⁻¹

µ, (3.34)

where δ

^>

= δ

₁

δ

₂

. It is thus important that the parameters of the interest rate model are chosen such that the quantity on the right hand side of (3.34) is of the order of magnitude of an interest rate.

We note that in the case of the mixed model, where we have restricted to the case where

λ

12

= 0, the matrix Λ is triangular, and its eigenvalues are thus given by its diagonal entries.

Thus the requirement that Λ have positive eigenvalues is in the case of the mixed model equivalent to requiring that λ

₁₁

, λ

₂₂

> 0.

In the case of the Rendleman-Bartter model, where we want to investigate both the mean- reverting and non-mean-reverting versions of the model, for the latter the matrix of coeffi- cients Λ should be chosen to have negative eigenvalues.

3.1.4 Meaning of Parameter Values

As we have just seen, the quantity Λ

⁻¹

µ features in the limit as t → ∞ of the expectation of the interest rate. If we let the parameters δ

₁

and δ

₂

be such that the interest rate is the weighted average of the two factors X

₁

(t) and X

₂

(t), this means that for the mean interest rate to converge to a sensible value, the quantity Λ

⁻¹

µ should be of the order of magnitude of an interest rate, as should the parameter δ

₀

. By the order of magnitude of an interest rate, we here mean roughly between −0.1 and 0.1, or between −10 percent and +10 percent.

Naturally, the initial values of the factors X

1

and X

2

, X

1

(0), and X

2

(0) should also be of the order of magnitude of an interest rate.

Finally, we should discuss the sensible range of values for the parameters σ

1

and σ

2

, which determine the scale of the volatility of the interest rate. We choose σ

1

and σ

2

to be positive, and between 1 and 10 percent.

3.2 Numerical Methods

3.2.1 Euler-Maruyama Finite Difference Schemes

To simulate paths of the interest rate R(t), we numerically solve the SDE (3.4). This is easily done using an Euler-Maruyama scheme, which for the general interest rate model (3.4) is given by

X

1

(t

n+1

) = X

1

(t

n

) + (µ

1

− λ

11

X

1

(t

n

) − λ

12

X

2

(t

n

)(t

n+1

− t

n

)

+ σ

1

X

1

(t

n

)

^γ1

(W

1

(t

n+1

) − W

1

(t

n

)) (3.35a) X

2

(t

n+1

) = X

2

(t

n

) + (µ

1

− λ

21

X

1

(t

n

) − λ

22

X

2

(t

n

))(t

n+1

− t

n

)

+ σ

₂

X

₂

(t

_n

)

^γ2

(W

₂

(t

_n+1

) − W

₂

(t

_n

)). (3.35b) In the case of the two factor Vasicek model and the two factor Rendleman-Bartter model, we simply use the Euler-Maruyama scheme (3.35) as it is. However, for the two factor CIR and mixed models some modifications need to be made.

The Brownian increments W

i

(t

n+1

) − W

i

(t

n

) are normally distributed with mean zero and variance t

_n+1

− t

n

, and are thus simulated by generating normally distributed random num- bers.

We have established conditions on the model parameters under which the factors that appear under square roots remain positive in the CIR and mixed models. However, even under these conditions the discretised versions of the processes appearing in the Euler-Maruyama scheme have a nonzero probability of becoming negative. If this happens in a numerical simulation, the next step in the scheme cannot be computed, as it would involve taking the square root of a negative number, yielding a complex number. This frequently happens in numerical simulations of the two-factor CIR model using the non-modified Euler-Maruyama scheme, especially for larger time steps.

A simple way to avoid this problem would be to keep generating each step in the scheme

with new random numbers until it yields a positive result, but doing this gives an algorithm

3. Interest Rate Models

that takes a random amount of time to run. Instead, we make modifications to the Euler- Maruyama scheme to ensure the positivity of the factors. Many different ways to do this exist, but here we will focus on two different modifications of the scheme.

The first is the so-called symmetrised Euler-Maruyama scheme, where one simply takes the absolute value of the regular scheme to obtain

X

1

(t

n+1

) = |X

1

(t

n

) + (µ

1

− λ

11

X

1

(t

n

) − λ

12

X

2

(t

n

)(t

n+1

− t

n

) + σ

1

p X

1

(t

n

)(W

1

(t

n+1

) − W

1

(t

n

))| (3.36a) X

2

(t

n+1

) = |X

2

(t

n

) + (µ

1

− λ

21

X

1

(t

n

) − λ

22

X

2

(t

n

))(t

n+1

− t

n

)

+ σ

2

p X

2

(t

n

)(W

2

(t

n+1

) − W

2

(t

n

))| (3.36b) for the CIR model. For the mixed model, we will use a half-symmetrised scheme of the form

X

₁

(t

n+1

) = |X

₁

(t

n

) + (µ

₁

− λ

₁₁

X

₁

(t

n

))(t

n+1

− t

n

) + σ

₁

p

X

₁

(t

_n

)(W

₁

(t

_n+1

) − W

₁

(t

_n

))| (3.37a) X

₂

(t

n+1

) = X

₂

(t

n

) + (µ

₁

− λ

₂₁

X

₁

(t

n

) − λ

₂₂

X

₂

(t

n

))(t

n+1

− t

n

)

+ σ

2

(W

2

(t

n+1

) − W

2

(t

n

)). (3.37b) Note that in the scheme for the mixed model above, we have set the parameter value λ

12

to zero, as discussed previously.

Strong convergence results for the symmetrised Euler-Maruyama scheme in the one dimen- sional case can be found in [2], for values of the exponent γ between 0.5 and 1.

A benefit of the symmetrised scheme is that the approximation processes preserve the pos- itivity of the factors, which not all schemes do.

The second scheme that we will consider is the scheme suggested in [14]. This scheme also preserves positivity in the approximation process. The main idea of the scheme is to sample random numbers from a different distribution, rather than from a normal distribution, thus creating a Markov chain approximation scheme which is guaranteed to remain positive.

Weak convergence of the whole path of the approximation process to that of the weak solution to the SDE is shown in [14] in the general multidimensional case.

In the general case, the scheme is formulated as follows. Consider the d-dimensional SDE

dX(t) = b(t, X(t))dt + σ(t, X(t))dW(t) (3.38) where the drift b : R

⁺

× R

^d

→ R

^d

and the diffusion σ : R

⁺

× R

^d

→ R

^d

⊗ R

^d

are both continuous functions, and R

^d

⊗ R

^d

is the space of d × d matrices. We also require that the SDE has a unique weak solution for each initial condition X(0) ∈ E := R

^m+

× R

^d−m

, such that X(t) ∈ E for all t ≥ 0. The approximating scheme with constant time step

_n¹

is then given by

X(k + 1) = X(k) + 1 n b  k

n , X(k)  + 1

√ n ˜ σ  k

n , X(k) 

(ε

k

− α) (3.39)

where the continuous function ˜ σ : R

⁺

×E → R

^d

⊗R

^d

is defined by the relation σσ

^>

= ˜ σΣ˜ σ

^>

,

and where Σ is a symmetric semi-definite positive d × d matrix. Further, ε

k

for k = 0,1,2...

are independent identically distributed random vectors with mean vector α and covariance matrix Σ, such that P(˜ σ(t,x)ε

_k

∈ E ∀(t, x) ∈ R

+

× E) = 1.

Furthermore, the following condition on the choice of mean vector and smallest allowed value of n, where

¹_n

is the time step must hold;

inf

(t,x)∈R+×E

 x + 1

n b(t,x) − 1

√ n ˜ σ(t, x)α



∈ E ∀n ≥ n

₀

. (3.40)

We note that both in the two factor CIR model and the mixed model, drift and diffusion functions are continuous. We have also previously noted that the SDE in both models have unique strong solutions, which implies that they also have unique weak solutions[13].

For the CIR model, we have

˜

σ = σ

₁

pX

₁

(t) 0 0 σ

2

pX

2

(t)



, (3.41)

b = µ

₁

− λ

11

X

1

− λ

12

X

2

µ

2

− λ

21

X

1

− λ

22

X

2



, (3.42)

and

Σ = 1 ρ ρ 1



, (3.43)

since the correlation between the two Brownian motions in the model is ρ.

Thus, the approximation scheme for the two factor CIR model is given by

X

1

(k + 1) = X

1

(k) + 1

n (µ

1

− λ

11

X

1

(k) − λ

12

X

2

(k)) + 1

√ n σ

1

p X

1

(k)(ε

k,1

− α

1

) (3.44a)

X

₂

(k + 1) = X

₂

(k) + 1

n (µ

₂

− λ

21

X

₁

(k) − λ

₂₂

X

₂

(k)) + 1

√ n σ

₂

p

X

₂

(k)(ε

_k,2

− α

2

) (3.44b)

The two factor CIR model is discussed specifically in [14], and it is found that the condition (3.40) holds when n

0

> max(λ

11

, λ

22

), and 0 < α

i

< 2 q

µ

i

(1 −

^λ_nⁱⁱ

0

) for i = 1,2, and the random vector ε has non-negative components. Further, it is possible to construct a non-negative random vector ε with mean vector α and covariance matrix Σ if and only if

−ρ ≤ α

1

α

2

. Thus it is possible to choose n

0

and α such that the condition (3.40) holds as long as −ρ < 4 √

µ

1

µ

2

.

Moving on to the mixed model, we have

σ = ˜ σ

1

pX

1

(t) 0

0 σ

2



, (3.45)

b =

 µ

1

− λ

11

X

1

µ

2

− λ

21

X

1

− λ

22

X

2



, (3.46)

3. Interest Rate Models

and the same correlation matrix Σ as in the CIR model.

Thus, the approximation scheme for the two factor mixed model is

X

₁

(k + 1) = X

₁

(k) + 1

n (µ

₁

− λ

₁₁

X

₁

(k)) + 1

√ n σ

₁

p

X

₁

(k)(ε

_k,1

− α

₁

) (3.47a) X

2

(k + 1) = X

2

(k) + 1

n (µ

2

− λ

21

X

1

(k) − λ

22

X

2

(k)) + 1

√ n σ

2

(ε

k,2

− α

2

). (3.47b)

The condition (3.40) becomes

inf

(x1,x2)∈R+×R

 x + 1

n b(t,x) − 1

√ n ˜ σ(t, x)α



∈ R

+

× R ∀n ≥ n

0

, (3.48)

and we note that it is enough to consider the condition on the first component, since clearly the second component remains real if the first component remains positive. Thus we have the condition

inf

x₁∈R+

 x

₁

+ 1

n (µ

₁

− λ

₁₁

x

₁

) − 1

√ n σ

₁

√ x

₁

α

₁



, (3.49)

which holds when we choose n

0

> λ

₁₁

and 0 < α

1

≤

_σ²

1

q µ

₁

(

^λ_n¹¹

0

), as discussed in [14].

Clearly we need to make a specific choice for the distribution of the random vector ε to implement the schemes (3.44) and (3.44). In [14], the suggestion of using a Bernoulli- type random vector is put forward, and it is noted that this is a good choice since it is computationally cheap to generate Bernoulli distributed random numbers. We shall make the same choice here, and outline the details below. Specifically, the distribution of ε is chosen so that for i = 1,2

α

i

ε

i

α

²_i

+ 1 ∼ Bern

 α

²_i

1 + α

²_i



, (3.50)

where ε

i

are the components of the random vector ε, and α

i

are the components of its mean.

Equivalently,

P(ε

i

= 0) = 1

1 + α

²_i

and P



ε

i

= α

_i²

+ 1 α

i



= α

²_i

1 + α

²_i

. (3.51) We note that it must hold that the quantity

_1+α^α2i2

i

∈ (0,1), as it is a probability.

To generate the random vector ε, we generate each of its components by first generating a uniformly distributed random number r

i

in the interval (0,1), and then setting ε

i

=

^α2i_α⁺¹

i

if r

i

<

_1+α^α2i2

i

and ε

i

= 0 otherwise.

3.2.2 Modelling SDE with Correlated Brownian Motions

In all the interest rate models that we consider, the correlation ρ between the two driving

Brownian motions is potentially nonzero. To numerically simulate the paths of the interest

rate factors X

₁

and X

₂

, as well as the interest rate itself, we must thus be able to numerically generate random numbers with constant correlation ρ ∈ (−1,1).

Given two independent Brownian motions W

1

(t) and W

2

(t), the process W

^ρ

(t) given by

W

^ρ

(t) = ρW

1

(t) + p

1 − ρ

²

W

₂

(t) (3.52)

is also a Brownian motion, and the correlation between W

₁

(t) and W

^ρ

(t) is equal to ρ[5]. Or equivalently, given a two-dimensional Brownian motion B(t) with independent components, the two-dimensional Brownian motion W(t) given by

W(t) = 1 0

ρ p

1 − ρ

²



B(t) (3.53)

has the correlation matrix

1 ρ ρ 1



. (3.54)

In the regular Euler-Maruyama schemes for the Vasicek model and the Rendleman-Bartter model, and in the symmetrised schemes for the CIR and mixed models, we can simply multiply each independent two dimensional Brownian motion increment by the matrix

1 0

ρ p

1 − ρ

²



(3.55)

to generate the correlated Brownian motion increments, since we are using normally dis- tributed random numbers with mean zero to model the Brownian increments.

In the weak scheme taken from [14], we instead use Bernoulli-type random numbers with nonzero mean, and then subtract the mean value. In this case, some more care is required.

Given a random vector ε with independent components ε

1

and ε

2

that follow the distribution (3.51), with α

1

= α

2

= a, the random vector ˜ ε defined by

˜

ε = 1 0

ρ p

1 − ρ

²



ε (3.56)

has correlation matrix

1 ρ ρ 1



(3.57) as required, and mean vector

˜ α =

 a

a(ρ + p 1 − ρ

²

)



. (3.58)

Thus, in the schemes (3.44) and (3.47) the mean vector ˜ α should be subtracted, rather than

the mean vector α of the original random vector ε with independent components. We note

that when ρ = 0 the two coincide.

3. Interest Rate Models

3.2.3 Numerical correlation of the factors X1 and X2

When the matrix of coefficients Λ is diagonal, the only source of correlation between two factors X

1

(t) and X

2

(t) in the general interest rate model should be the correlation ρ between the two Brownian motions. We now compare the numerical correlation between the factors X

1

(t) and X

2

(t) to the theoretical correlation ρ between the driving Brownian motions.

We make comparisons in all the four studied interest rate models, and using the different discretisation schemes that have been discussed.

For the numerical correlation, we generate 10,000 paths each of X

1

and X

2

, and then look at the numerical correlation between X

1

(T ) and X

2

(T ), where T is the endpoint of the simulated paths. It can be seen in table 3.1 that the numerical correlation agrees well with the correlation ρ between the Brownian motions in all cases.

Vasicek Euler-Maruyama Scheme

ρ -1 -0.9 -0.5 0 0.5 0.9 1

Num. Correlation -0.9978 -0.8977 -0.5047 -0.0048 0.5029 0.8989 1.0 CIR Weak Scheme

ρ -1 -0.9 -0.5 0 0.5 0.9 1

Num. Correlation -0,9743 -0.8786 -0.4888 0.0080 0.4960 0.8988 1.0 CIR Symmetrised Scheme

ρ -1 -0.9 -0.5 0 0.5 0.9 1

Num. Correlation -0,8793 -0.7995 -0.4643 0.0124 0.4801 0.8943 1.0 Mixed Model Weak Scheme

ρ -1 -0.9 -0.5 0 0.5 0.9 1

Num. Correlation -0,9933 -0.8933 -0.5022 4 · 10

⁻⁴

0.4994 0.8960 0.9934 Mixed Model Symmetrised Scheme

ρ -1 -0.9 -0.5 0 0.5 0.9 1

Num. Correlation -0,9727 -0.8757 -0.4818 0.0074 0.4854 0.8742 0.9731 Rendleman-Bartter Euler-Maruyama Scheme

ρ -1 -0.9 -0.5 0 0.5 0.9 1

Num. Correlation -0,9941 -0.8924 -0.5028 0.0129 0.4854 0.8998 1.00

Table 3.1: Table showing comparisons between the numerical correlation between the

factors X

₁

(t) and X

₂

(t) and the correlation ρ between the two Brownian motions in each

model. The correlation agrees well for all interest rate models, and for all different numerical

schemes studied.

4

Zero Coupon Bond Pricing

4.1 Theory

4.1.1 The Pricing PDE

Here, we derive the pricing PDE for a zero coupon bond in a two-factor interest rate model of the form studied in this thesis.

Recall the general interest rate model (3.4) that was introduced in chapter 3,

R(t) = δ

0

+ δ

1

X

1

(t) + δ

2

X

2

(t) (4.1a) dX

₁

(t) = (µ

1

− λ

11

X

₁

(t) − λ

12

X

₂

(t))dt + σ

1

(X

1

(t))

^γ1

dW

₁

(t) (4.1b) dX

2

(t) = (µ

2

− λ

21

X

1

(t) − λ

22

X

2

(t))dt + σ

2

(X

2

(t))

^γ2

dW

2

(t), (4.1c) and the risk-neutral pricing formula for the zero coupon bond with maturity T

B(t,T ) = ˜ E

 exp 

− Z

T

t

R(s)ds 

| F (t)



. (4.2)

Recall further that the risk-neutral pricing formula was derived by requiring that the dis- counted ZCB price be a martingale in the risk-neutral probability measure, and that the interest rate model is formulated in the risk neutral probability measure.

Since the short rate R(t) is a function of the factors X

1

(t) and X

2

(t), which solve the stochastic differential equations (4.1b)-(4.1c) and are thus Markov processes, there must exist some function f (t,x

1

,x

2

) such that

B(t,T ) = f (t,X

₁

,X

₂

). (4.3)

Since the discounted bond price is a martingale in the risk neutral probability measure, the differential d(D(t)B(t,T )) must have dt-term zero. Here, as previously, D(t) is the discounting factor, which is given by

D(t) = exp 

− Z

t

0

R(s)ds 

, (4.4)

and satisfies the differential equation dD(t) = −R(t)D(t)dt.