MAXIMUM LIKELIHOOD ESTIMATION OF THE COX-INGERSOLL-ROSS PROCESS : THE MATLAB IMPLEMENTATION

The square root diffusion process is widely used for modeling interest rates behaviour. It is an underlying process of the well-known Cox-Ingersoll-Ross term structure model (1985). We investigate maximum likelihood estimation of the square root process (CIR process) for interest rate time series. The MATLAB implementation of the estimation routine is provided and tested on the PRIBOR 3M time series. 1 CIR Process for Interest Rate Modeling A continuous-time model in finance typically rest on one or more stationary diffusion processes {Xt, t ≥ 0}, with dynamics represented by stochastic differential equations: dXt = μ(Xt)dt + σ(Xt)dWt, (1) where {Wt, t ≥ 0} is a standard Brownian motion. The functions μ(·) and σ2(·) are, respectively, the drift and the diffusion functions of the process. The fundamental process in interest rate modeling is the square root process given by the following stochastic differential equation: drt = α(μ− rt)dt + √ rtσdWt, (2) where rt is the interest rate and θ ≡ (α, μ, σ) are model parameters. The drift function μ(rt, θ) = α(μ − rt) is linear and possess a mean reverting property, i.e. interest rate rt moves in the direction of its mean μ at speed α. The diffusion function σ(rt, θ) = rtσ is proportional to the interest rate rt and ensures that the process stays on a positive domain. The square root process (2) is the basis for the Cox, Ingersoll, and Ross short-term interest rate model [1] and therefore often denoted as the CIR process in the financial literature. 1.1 CIR process densities If α, μ, σ are all positive and 2αμ ≥ σ2 holds, the CIR process is well-defined and has a steady state (marginal) distribution. The marginal density is gamma distributed. For maximum likelihood estimation of the parameter vector θ ≡ (α, μ, σ) transition densities are required. The CIR process is one of few cases, among the diffusion processes, where the transition density has a closed form expression. We follow the notation given in [1] on page 391. Given rt at time t the density of rt+∆t at time t + ∆t is p(rt+∆t|rt; θ, ∆t) = ce−u−v( v u ) q 2 Iq(2 √ uv), (3) This work was supported by grant of IGA VSE nr. IG410046 where c = 2α σ2(1− e−α∆t) , u = crte, v = crt+∆t, q = 2αμ σ2 − 1, and Iq(2 √ uv) is modified Bessel function of the first kind and of order q. The transition density (3) has been originally derived in [2]. Sometimes, it is particularly useful to work with a transformation st+∆t = 2crt+∆t. We can easily derive that the transition density of st+∆t is g(st+∆t|st; θ, ∆t) = g(2crt+∆t|2crt; θ, ∆t) = 1 2c p(rt+∆t|rt; θ, ∆t), (4) which is the the noncentral χ2 distribution with 2q + 2 degrees of freedom and non-centrality parameter 2u. 2 Maximum Likelihood Implementation in MATLAB Parameter estimation is carried out on interest rate time series with N observations {rti , i = 1 . . . N}. We consider equally spaced observations with ∆t time step. 2.1 Likelihood function The likelihood function for interest rate time series with N observations is