Gaussian Estimation of One-Factor Mean Reversion Processes

Abstract

We propose a new alternative method to estimate the parameters in one-factor mean reversion processes based on the maximum likelihood technique. This approach makes use of Euler-Maruyama scheme to approximate the continuous-time model and build a new process discretized. The closed formulas for the estimators are obtained. Using simulated data series, we compare the results obtained with the results published by other authors.

1. Introduction

Different applications of stochastic differential equations have attracted the attention of many researchers in relation to the theoretical advances needed for modeling and simulation of the dynamic behavior of the main variables affecting the evolution of some phenomena. Applications in fields of economics, finance, physics, insurance, and others have been drivers of these developments. Particularly, through diffusion processes have been carried out modeling interest rates and commodity prices, where they used in their modeling processes mean reversion. Some trends of these studies are based on some a priori knowledge to determine the structure of the model to be used, while other alternatives are part of a general structure, and statistical procedures through the particular structure are determined to represent the phenomenon. Whichever approach is used, in general, in the specification of the model, it is necessary to estimate parameters to achieve a proper fit of the model.

Although some parameter estimation methods have been extensively studied in diffusion processes, this has a serious problem, which is that the specification of the models is done in continuous time but the data that is available for adjusting the parameters are discretely sampled in time. A common solution to this problem is to discretize the continuous-time model based on Euler-Maruyama scheme and perform procedures on the theoretical formulation resulting in discrete time. From this scheme alternatives can be proposed for parameter estimation approach among which stands out the maximum likelihood estimation.

Following the classification proposed by Yu and Phillips [1], different parameter estimation procedures can be classified globally into three categories according to how to proceed from the model formulation: discretization of continuous process and estimation of the parameters with the resulting discretized process, estimation of the transition probability function for the continuous-time model, and independent estimations of the trend and volatility components by kernel methods.

Thus, among the methods that are based on the discretization of the continuous-time model is Nowman [2], who considers a stochastic differential equation with two unknown variables in the drift term and two in the volatility, elasticity allowing the variance. Similarly, Yu and Phillips [1] present an improvement to the method proposed by Nowman [2], which is based on a nonuniform sampling of the observations, and the estimation process makes use of the equivalence of the methods of maximum likelihood and least squares under the conditions laid down in the model. Alternatively, Yu [3] proposes a new scheme that is based on the idea of transforming a stochastic differential equation (EDE) in its discrete-time equivalent, an AR(1) to determine the formulas that approximate the bias error present in the estimation of one of the parameters.

Some of the works that could be classified as methods in which it seeks to find the transition probability function and use Martingale properties to approximate the likelihood function are Durham and Gallant [4], Valdivieso et al. [5], and Elerian et al. [6]. Durham and Gallant [4] propose a scheme based on simulation. Establishing different model modifications, Pedersen [7] improves computational performance, and characteristics of the model are studied using simulated synthetic series. Koulis and Thavaneswaran [8] study combined martingale estimation functions for discretely observed interest rates models. They derive closed-form expressions for the first four conditional moments of the CIR model via Ito’s formula and for extended versions of the CIR model using the Milstein’s approximation and construct the optimal estimating functions. They also show that combined estimating functions have more information when the conditional mean and variance of the observed process depend on the parameter of interest. Valdivieso et al. [5] proposed a method for estimating the transfer function through the inversion of the characteristic function of the process and the use of Fourier transform. This methodology is developed for some particular cases of D-OU processes (Ornstein-Uhlenbeck processes with broken car distribution D). Elerian et al. [6] used simulation-based techniques for estimating the transition probability depending on the sample to incorporate auxiliary dummy observations. The parameters are obtained by simulation and optimization using numerical algorithms.

Methods based on kernel techniques approximate density functions of the components and spreading tendency of the process. The works of Jiang and Knight [9] and Florens-Zmirou [10] are within this framework.

In this paper, we propose an alternative scheme for parameter estimation in one-factor mean reversion models, using the Euler-Maruyama discretization. Using the properties of discretized error normality, the maximum likelihood technique to obtain closed formulas for the estimators applies. We present some performance measures of the maximum likelihood estimators obtained with the proposed method in simulated processes, and we compare them with the results published by Nowman [2], Yu and Phillips [1], and Durham and Gallant [4].

This paper is organized as follows. In Section 2, a brief description of mean reversion processes of one factor with constant parameters. Some Gaussian estimation techniques are presented in Section 3. Section 4 develops the proposed method, and Section 5 presents some experimental numerical results with different sampling rates.

2. Mean Reversion Processes

Mean reversion processes have been considered as a natural choice to describe the dynamic behavior of interest rates in Vasicek [11], Brennan and Schwartz [12], Cox et al. [13], and Chan et al. [14]; daily prices of some commodities such as oil and natural gas in Lari-Lavassani et al. [15] have a fundamental role in modeling spot prices when they exhibited long-term trends, such as electrical energy (Pilipovic [16]).

Mean reversion processes of one factor with constant parameters can be written as 𝑑𝑋_𝑡=𝛼𝜇−𝑋_𝑡𝑑𝑡+𝜎𝑋^𝛾_𝑡𝑑𝐵_𝑡;𝑡∈[0,𝑇] (1) with initial condition 𝑋₀=𝑥, where 𝛼>0, 𝜎>0, 𝜇∈ℝ, and 𝛾∈[0,3/2] are constants and {𝐵_𝑡}_𝑡≥0 is a unidimensional standard Brownian motion defined on a complete probability space (Ω,ℱ,ℙ).

The 𝛼 parameter is called reversion rate, 𝜇 is the mean reversion level or trend of long-run equilibrium, 𝜎 is the parameter associated with the volatility, and 𝛾 determines the sensitivity of the variance to the level of 𝑋_𝑡.

The model (1) is known as CKLS proposed by Chan et al. [14] and is a generalization of the models of interest rate in the short term with constant parameters in which 𝛾 = 0, 1/2, and 1 or 𝛾=3/2 with 𝛼=0.

In general, mean reversion processes are stationary. For example, note that (1) can be written in integral form as 𝑋_𝑡−𝑋₀=𝛼^𝑡₀𝜇−𝑋_𝑠𝑑𝑠+𝜎^𝑡₀𝑋^𝛾_𝑠𝑑𝐵_𝑠. (2)

Taking the expected value, 𝐸[𝑋_𝑡]−𝐸[𝑋₀]=𝛼∫^𝑡₀(𝜇−𝐸[𝑋_𝑠])𝑑𝑠 and thus we obtain the differential equation 𝑚̇(𝑡)=𝛼𝜇−𝑚(𝑡),𝑚(𝑡)=𝐸𝑋_𝑡. (3)

The solution of this equation is 𝑚(𝑡)=𝜇+(𝑋₀−𝜇)𝑒^−𝛼𝑡, so lim_𝑡→∞𝐸[𝑋_𝑡]=lim_𝑡→∞𝑚(𝑡)=𝜇.

It can be concluded then that 𝑚(𝑡)=𝜇 is the equilibrium level for the linear differential equation. In other words, the expected value of the process 𝑚(𝑡) returns in the long term to 𝜇. In a finite-time horizon, 𝜇 plays the role of attractor in the sense that when 𝑋_𝑡>𝜇 the trend term 𝛼(𝜇−𝑋_𝑡)<0, and therefore, 𝑋_𝑡 decreases and when 𝑋_𝑡<𝜇 a similar argument establishes that 𝑋_𝑡 grows.

The stationarity condition is the key to perform the Gaussian estimation of parameters 𝛼,𝜎, and 𝜇, assuming that the elasticity (𝛾=0,1/2,1,3/2) is given a priori.

Although the general model (1) derives various interest rate models studied in the literature, our interest is focused primarily on those models that are strictly stationary.

For example, when 𝛾=0, it is said that the process has additive noise so that we get the linear stochastic differential equation 𝑑𝑋_𝑡=𝛼𝜇−𝑋_𝑡𝑑𝑡+𝜎𝑑𝐵_𝑡;𝑋₀=𝑥,𝑡∈[0,𝑇]. (4)

This model is also known as Ornstein-Uhlenbeck mean reversion process type used by Vasicek [11] in modeling the equilibrium price of the coupon. This Gaussian process has been widely used in the valuation of options on bonds, futures, options on futures, and other financial derivative instruments, while in Lari-Lavassani et al. [15] it has been used to study the dynamic behavior of the prices of oil and natural gas.

Although the model of Vasicek [11] is analytically tractable, it has some drawbacks. Since the interest rate in the short term is normally distributed, for each 𝑡∈[0,𝑇] exists a probability that 𝑋_𝑡 is negative, and this is not reasonable from an economic standpoint. Another drawback of this model is that conditional volatility assumes changes in the interest rate that are constants, independent of the level of 𝑋_𝑡.

If 𝛾=1 it is said that the process has proportional noise and gives the nonhomogeneous linear stochastic differential equation given by 𝑑𝑋_𝑡=𝛼𝜇−𝑋_𝑡𝑑𝑡+𝜎𝑋_𝑡𝑑𝐵_𝑡;𝑋₀=𝑥,𝑡∈[0,𝑇]. (5)

This model is used by Brennan and Schwartz [12] to obtain a numerical model of convertible bond prices. This process is also used by Courtadon [17] in developing a model for option pricing on discount bonds.

When 𝛾=1/2, we obtain the nonlinear stochastic differential equation 𝑑𝑋_𝑡=𝛼𝜇−𝑋_𝑡𝑑𝑡+𝜎𝑋_𝑡^1/2𝑑𝐵_𝑡,2𝛼𝜇≥𝜎²,𝑋₀=𝑥,𝑡∈[0,𝑇]. (6)

This process is naturally linked to the 𝜒 centered-square distribution and is the (square root) version of the model of Brennan and Schwartz [12], proposed by Cox et al. [18], which produces types of instantaneous short-term interest strictly positive.

This model has also been used extensively in the development of pricing models for derivative instruments sensitive to interest rates. Examples include pricing models for mortgage-backed in Dunn and McConnell [19], models of options for discount bonds in Cox et al. [18], future and option pricing models for futures in Ramaswamy and Sundaresan [20], the swap valuation models in Sundaresan [21], and option pricing models for the convenience yield in Longstaff [22].

Finally, for the special case in which 𝛾=3/2 and 𝛼=0, equation (1) takes the form 𝑑𝑋_𝑡=𝜎𝑋_𝑡^3/2𝑑𝐵_𝑡,𝑋₀=𝑥,𝑡∈[0,𝑇]. (7) This model is introduced by Cox et al. [13] in their study of variable-rate (VR) securities. A similar model is also used by Constantinides and Ingersoll [23] for value bonds in the presence of taxes.

Table 1 shows the summary of the models associated with (1) that are of interest in this work (Vasicek, Brennan-Schwartz, CIR-SR,CIR-VR, and CKLS).

Table 1 

Summary of the models.

3. Gaussian Estimation of Diffusion Models

Many theoretical models in economics and finance are formulated in continuous time as either a diffusion or diffusion systems. For this reason, identification, estimation, and, therefore, the investigation of the asymptotic properties of samples of the estimators of these processes turn out to be a difficult task. These problems arise mainly due to the lack of availability of continuous sampling observations. In fact, the estimation of diffusion models has been considered in the literature for many years, and in most scientific journals it only refers to the estimate considering continuous sampling observations, although the data that these models described can only be sampled at discrete points in time (See, e.g., Jiang and Knight [9] and references therein.)

There are various techniques and methods for estimation of models of mean reversion of a single factor; with different structural forms in the trend and distribution functions, most simulation techniques-based on statistical and quite rigorous arguments. The work of Jiang and Knight [9] proposes a parametric identification and estimation procedure for diffusion processes based on discrete sampling observations. The nonparametric kernel estimator for diffusion functions developed in this work deals with general diffusion processes and avoids any functional form specification for either the trend function or distribution function. They show that under certain regularity conditions, the nonparametric estimator function is consistent, and timely dissemination also asymptotically follows a mixed normal distribution. In Jiang and Knight [24], there are some numerical experimental results obtained by a Monte Carlo study using the results of Jiang and Knight [9], which include parameter estimation for Ornstein-Uhlenbeck process. In this sense, the lack of analytical solutions has hampered empirical evidence and validation of alternative hypotheses about these continuous-time models.

Although the maximum likelihood method is a widely used parametric method, the analytical expressions for the likelihood functions exist only for processes of mean reversion with additive noise. In fact, in most of the work on estimation of diffusion models, processes become more complex processes of Ornstein-Uhlenbeck type. For example, Giet and Lubrano [25] consider a diffusion process in which the elasticity constant of the nonlinear term of diffusion can be freely chosen, achieving after reducing it a transformation to an Ornstein-Uhlenbeck process. The major drawback is that there is no degree of freedom in the form of the trend function which is completely determined by the shape of the distribution function and the choice of the underlying linear model. Valdivieso et al. [5] proposed a methodology to estimate maximum likelihood parameters of a one-dimensional stationary process of Ornstein-Uhlenbeck type, which is driven by general Levy process and which is constructed through a distribution D self-decomposed. This approach is based on the inversion of the characteristic function and using the discrete fractional fast Fourier transform. The maximization of the likelihood function is performed using numerical methods.

The framework proposed in our work has some features in common with Nowman [2], Yu and Phillips [1], and Elerian et al. [6]. The similarity of our approach lies in the use of Euler numerical scheme, using discretely sampled data and maximum likelihood estimation, but the main difference is that in our proposal does not need any numerical optimization because it is possible to obtain a closed formula for the estimators.

For example, in Elerian et al. [6] developed a technique based on Bayesian estimation of nonlinear EDEs approximating the Markov transition densities of broadcasts, including Ornstein-Uhlenbeck processes and elasticity models for interest rates, when observations are discretely sampled. The frame of the estimate is based on the introduction of auxiliary data to complete the diffusion latent missing between each pair of measurements. Combine methods of Monte Carlo Markov Chain algorithm based on the Metropolis-Hastings, in conjunction with the numerical scheme of Euler-Maruyama, to sample the posterior distribution of latent data and model parameters.

Another alternative for Gaussian estimation was proposed by Nowman [2]. In this paper, an estimate using the exact discretization of the model is proposed, 𝑑𝑋_𝑡=𝛼𝜇−𝑋_𝑡𝑑𝑡+𝜎𝑋^𝛾_()𝑠Δ𝑑𝐵_𝑡,𝑠Δ≤𝑡<(𝑠+1)Δ, (8) assuming that volatility stays constant on the intervals [𝑠Δ,(𝑠+1)Δ); 𝑠=0,1,….

In this way, we obtain an equation of the form 𝑋_𝑡=𝑒^−𝛼Δ𝑋_(𝑡−Δ)+𝜇1−𝑒^−𝛼Δ+𝜂_𝑡, (9) where the conditional distribution 𝜂_𝑡∣𝐹^𝐵_𝑡−1≈𝑁(0,(𝜎²/2𝛼)(1−𝑒^−2𝛼Δ){𝑋_𝑡−1}^2𝛾).

However, despite the obtainable approximate conditional distribution, this approach does not correspond to the Euler discretization applied to (1). This estimation technique requires the discrete model corresponding to (2) which is given by 𝑋_𝑡=𝑒^−𝛼𝑋_(𝑡−1)+𝜇1−𝑒^−𝛼+𝜂_𝑡,(𝑡=0,1,…,𝑇), (10) where 𝜂_𝑡, (𝑡=0,1,…,𝑇) satisfies the conditions 𝐸[𝜂_𝑠𝜂_𝑡]=0; (𝑠≠𝑡) and 𝐸𝜂²_𝑡=^𝑡_𝑡−1𝑒^{−2(𝑡−𝜏)𝛼}𝜎²𝑋_(𝑡−1)^2𝛾=𝜎𝑑𝜏₂2𝛼1−𝑒^−2𝛼𝑋_(𝑡−1)^2𝛾=𝑚²_𝑡𝑡. (11) This states that the logarithm of the Gaussian likelihood function may be written as 𝐿(𝜃)=^𝑇_𝑡=12log𝑚_𝑡𝑡+𝑋_𝑡−𝑒^−𝛼𝑋_(𝑡−1)+𝜇1−𝑒^−𝛼𝑚_𝑡𝑡²=^𝑇_𝑡=12log𝑚_𝑡𝑡+𝜀²_𝑡, (12) where 𝜀=[𝜀₁,…,𝜀_𝑇] is a vector whose elements can be calculated.

The optimization of 𝐿(𝜃) is carried out using the methods developed by Bergstrom [26–30] and Nowman [2].

Yu and Phillips [1] propose a Gaussian estimation using random time changes which is based on the idea that any continuous time martingale can be written as a Brownian motion after a suitable time change. It considers the explicit solution (1) to obtain the expression 𝑋_(𝑡+ℎ)=𝜇1−𝑒^−𝛼ℎ+𝑒^−𝛼ℎ𝑋_(𝑡)+^ℎ₀𝜎𝑒^{−𝛼(ℎ−𝜏)}𝑋_(𝑡+𝜏)^𝛾𝑑𝐵_𝜏forℎ>0. (13)

Let 𝑀(ℎ)=𝜎∫^ℎ₀𝑒^{−𝛼(ℎ−𝜏)}{𝑋_(𝑡+𝜏)}^𝛾𝑑𝐵_𝜏.  𝑀(ℎ) is a continuous martingale (In Yu and Phillips [1],   𝑀(ℎ) was assumed as a continuous martingale, but this assumption is generally not warranted, and [𝑀]_ℎ does not induce a Brownian motion as Yu and Phillips showed in [31]. In addition, they establish that a simple decomposition splits the error process into a trend component and a continuous martingale process.) with quadratic variation [𝑀]_ℎ=𝜎^2ℎ₀𝑒^{−2(ℎ−𝜏)𝛼}𝑋_(𝑡+𝜏)^2𝛾𝑑𝜏. (14)

They introduce a sequence of positive numbers {ℎ_𝑗} which deliver the required time changes. For any fixed constant 𝑎>0, let ℎ_𝑗+1=inf𝑠∣𝑀_𝑗𝑠≥𝑎=inf𝑠∣𝜎^2𝑠₀𝑒^{−2(𝑠−𝜏)𝛼}𝑋_{(𝑡𝑗+𝜏)}^2𝛾𝑑𝜏≥𝑎. (15)

Then a succession of points {𝑡_𝑗} is constructed using the iterations 𝑡_𝑗+1=𝑡_𝑗+ℎ_𝑗+1 with 𝑡₁=0 to obtain 𝑋_(𝑡𝑗+1)=𝜇1−𝑒^{−𝛼ℎ𝑗+1}+𝑒^{−𝛼ℎ𝑗+1}𝑋_(𝑡𝑗)+𝑀_ℎ𝑗+1, (16) where 𝑀_ℎ𝑗+1=𝐵(𝑎)≈𝑁(0,𝑎).

Note that according to (16), (12) takes the form 𝐿(𝜃)=1𝑁_𝑗+𝑋2log𝑎_(𝑡𝑗+1)−𝜇1−𝑒^{−𝛼ℎ𝑗+1}−𝑒^{−𝛼ℎ𝑗+1}𝑋_(𝑡𝑗)𝑎². (17)

Since the data are observed in discrete time, the time change formula cannot be applied directly and therefore must be approximated by ℎ_𝑗+1=Δmin𝑠∣^𝑠_𝑖=1𝜎²𝑒^{−2𝛼(𝑠−𝑖)Δ}𝑋_{(𝑡𝑗+𝑖Δ)}^2𝛾≥𝑎. (18)

This should be selected so that 𝑎 is the volatility of the error term unconditional in the equation 𝑋_(𝑡+Δ)=𝜇1−𝑒^−𝛼Δ+𝑒^−𝛼Δ𝑋_(𝑡)+𝜀. (19)

The implementation of this method is based on the estimate given by Nowman [2] to find the parameters estimates 𝜎 and 𝛾. Furthermore, the estimation of the remaining parameters 𝜇 and 𝛼 is carried out using the numerical optimization algorithm by Powell [32].

4. Gaussian Estimation of One-Factor Mean Reversion Processes with Constant Parameters

The development of this section will focus on the description of the proposed technique for the Gaussian estimation of mean reversion models of a single factor, with linear and nonlinear functions of diffusion constant parameters and linear trend functions with constant parameters.

As the density function for most diffusion processes is not known analytically, the proposed idea is not to find or approximate the transition density function or the marginal density function of the process. Otherwise, we use the numerical approximation of Euler on the EDE to build a new discrete process, Huzurbazar [33], and use the fact that in this approximation the normalized residuals resulting from Brownian process are independent standard normal random variables. As pointed out by Pedersen [7], Kohatsu-Higa and Ogawa [34], and Bally and Talay [35], the density function of the discrete process obtained with Euler scheme converges to the density function of the continuous process, and as discussed by Zellner [36], these residues can be interpreted as a vector of unknown parameters with the advantage of knowing their distribution. So we build a new variable for the models defined in (1) so that it can adapt to the observed data and later can be evaluated by examining the normalized residuals that will be useful in the estimation process. By knowing the density function of this new process we can find discretized likelihood function and use the conventional method of maximum likelihood estimation to find a closed formula for the parameters.

Consider the stochastic differential equation (1), with the initial condition 𝑋₀=𝑥. The terms 𝛼(𝜇−𝑋_𝑡) and 𝜎𝑋^𝛾_𝑡 being the trend and diffusion functions, respectively, that depend on 𝑋_𝑡 and a vector of unknown parameters 𝜃=[𝛼,𝜇,𝜎,𝛾], where {𝐵_𝑡}_𝑡≥0 is a unidimensional standard Brownian motion defined on a complete probability space (Ω,ℱ,ℙ). The trend function describes the movement of the process due to changes in time, while the diffusion function measures the magnitude of the random fluctuations around the trend function. Suppose that the conditions under which we guarantee the existence and uniqueness of the solution of Mao [37] are satisfied. Assume further the process observations in discrete time 𝑋_𝑡=𝑋(𝜏_𝑡) at times {𝜏₀,𝜏₁,…,𝜏_𝑇} where Δ_𝑡=𝜏_𝑡−𝜏_𝑡−1≥0, 𝑡=1,2,…,𝑇. The objective is to estimate the unknown parameter vector 𝜃 given the historical information of the sample 𝑋=(𝑋₀,𝑋₁,…,𝑋_𝑇)^′ with equally spaced data. Note that the explicit solution of (1) is given by 𝑋_𝑡=𝜇+𝑋₀−𝜇𝑒^−𝛼𝑡+^𝑡₀𝜎𝑒^{−𝛼(𝑡−𝑠)}𝑋^𝛾_𝑠𝑑𝐵_𝑠. (20)

In the context of the likelihood, the estimation of 𝜃 is based on the likelihood function log𝐿𝑋₁,𝑋₂,…,𝑋_𝑇∣𝑋₀,𝜃=^𝑇_𝑡=1log𝑔𝑋_𝑡+1∣𝑋_𝑡,𝜃, (21) where 𝑔(𝑋_𝑡+1∣𝑋_𝑡,𝜃) are the Markov transition densities, or 𝐿(𝜃∣𝑋₁,𝑋₂,…,𝑋_𝑇)=𝑓(𝑋₁,𝑋₂,…,𝑋_𝑇,𝜃) and 𝑓(𝑋₁,𝑋₂,…,𝑋_𝑇,𝜃) is the joint density function.

Although the differential equation (1) can be solved analytically on [0,𝑇], the functions 𝑔(𝑋_𝑡+1∣𝑋_𝑡,𝜃) and 𝑓(𝑋₁,𝑋₂,…,𝑋_𝑇,𝜃) are not available in closed form, and, therefore, the estimate can not be accomplished, except for the case where 𝛾=0.

Note that the exact discretized solution of (1) takes the form 𝑋𝜏_𝑡=𝜇+𝑋𝜏_𝑡−1−𝜇𝑒^{−𝛼(𝜏𝑡−𝜏𝑡−1)}+^𝜏_𝜏𝑡−1𝜎𝑒^{−𝛼(𝜏𝑡−𝑠)}𝑋^𝛾_(𝑠)𝑑𝐵_𝑠 (22) which can be approximated by 𝑋𝜏_𝑡=𝜇+𝑋𝜏_𝑡−1−𝜇𝑒^{−𝛼Δ𝑡}+𝜎𝑒^{−𝛼Δ𝑡}𝑋^𝛾𝜏_𝑡−1Δ𝐵_𝑡, (23) where Δ𝐵_𝑡=𝐵_𝜏𝑡−𝐵_𝜏𝑡−1.

Thus, 𝑋𝜏_𝑡=𝑋𝜏_𝑡−1+𝛼𝜇−𝑋𝜏_𝑡−1Δ_𝑡+𝑋^𝛾𝜏_𝑡−1𝜀_𝜏𝑡, (24) where 𝜀_𝜏𝑡∼𝑁(0,𝜎²Δ_𝑡), corresponding to the numerical approximation of Euler.

Now define the new variable 𝑌_𝑖=𝑋_𝑖−𝑋_𝑖−1−𝛼𝜇−𝑋_𝑖−1Δ𝑋^𝛾_𝑖−1=𝜀_𝑖,Δconstant. (25)

Note that 𝑌_𝑖, 𝑖=1,2,3,…,𝑇 are independent random variables distributed normal, 𝑁(0,𝜎²Δ).

In terms of the density function can be written: 𝑓𝑌₁=1∶𝜃2𝜋𝜎²Δ^1/2×exp−12𝜎²Δ𝑋₁−𝑋₀−𝛼𝜇−𝑋₀Δ𝑋^𝛾₀²𝑓𝑌₂=1∶𝜃2𝜋𝜎²Δ^1/2×exp−12𝜎²Δ𝑋₂−𝑋₁−𝛼𝜇−𝑋₁Δ𝑋^𝛾₁²⋮𝑓𝑌_𝑇=1∶𝜃2𝜋𝜎²Δ^1/2×exp−12𝜎²Δ𝑋_𝑇−𝑋_𝑇−1−𝛼𝜇−𝑋_𝑇−1Δ𝑋^𝛾_𝑇−1². (26)

Since each 𝑌_𝑖 is independent can be expressed the joint density function as their marginal densities product, that is, 𝑓𝑌₁,𝑌₂,…,𝑌_𝑇=𝑓𝑌₁∶𝜃𝑓𝑌₂∶𝜃⋯𝑓𝑌_𝑇∶𝜃. (27)

Consequently, the likelihood function is given by 𝐿𝜃∣𝑌_𝑇=12𝜋𝜎²Δ^𝑇/2×exp−12𝜎²Δ^𝑇_𝑖=1𝑋_𝑖−𝑋_𝑖−1−𝛼𝜇−𝑋_𝑖−1Δ𝑋^𝛾_𝑖−1². (28)

Taking the natural logarithm on both sides, it follows that log𝐿=−𝑇2log2𝜋𝜎²Δ−12𝜎²Δ^𝑇_𝑖=1𝑋_𝑖−𝑋_𝑖−1−𝛼𝜇−𝑋_𝑖−1Δ𝑋^𝛾_𝑖−1². (29)

Now consider the problem of maximizing the likelihood function given by 𝜕log(𝐿)=𝜕𝜃0,i.e.,𝜃=Arg_𝜃maxlog(𝐿). (30)

The estimation process leads to the estimated parameters 𝛼 and 𝜇, are the solutions of which the next equation system, where it is assumed that the value of 𝛾 is known. Consider ()()𝐴−𝐸1−𝛼Δ−𝛼𝜇Δ𝐵=0𝐶−𝐵1−𝛼Δ−𝛼𝜇Δ𝐷=0, (31) where 𝐴=^𝑇_𝑖=1𝑋_𝑖𝑋_𝑖−1𝑋^2𝛾_𝑖−1,𝐵=^𝑇_𝑖=1𝑋_𝑖−1𝑋^2𝛾_𝑖−1,𝐶=^𝑇_𝑖=1𝑋_𝑖𝑋^2𝛾_𝑖−1,𝐷=^𝑇_𝑖=11𝑋^2𝛾_𝑖−1,𝐸=^𝑇_𝑖=1𝑋_𝑖−1𝑋^𝛾_𝑖−1². (32) Thus, the estimated parameters are given by 𝛼=𝐸𝐷−𝐵²−𝐴𝐷+𝐵𝐶𝐸𝐷−𝐵²Δ,()𝜇=𝐴−𝐸1−𝛼Δ.𝛼𝐵Δ (33)

Additionally, the parameter 𝜎 is determined by an equation independent of the system of (31) and is completely determined by the estimated parameters: 𝜎=1𝑇Δ^𝑇_𝑖=1𝑋_𝑖−𝑋_𝑖−1−𝛼𝜇−𝑋_𝑖−1Δ𝑋^𝛾_𝑖−1². (34)

This estimation scheme has several important qualities. First, you come to a closed formula for Gaussian estimation without numerical optimization methods. Second, the computational implementation is simple and does not require sophisticated routines to perform calculations very extensive.

5. Simulations

This section will present comparisons of the results with those published by other authors using the techniques described above. Comparisons are made with the estimate based on simulated data sets to analyze statistical properties of the proposed estimators.

Yu and Phillips [1] present the results of a simulation study with the Monte Carlo technique in comparing their method with Nowman’s [2]. To make this comparison are suggested three scenarios: to resolve data monthly, weekly, and daily. In all cases the basic pattern of square root is defined as 𝑑𝑋_𝑡=𝜆+𝛽𝑋_𝑡𝑑𝑡+𝜎𝑋^𝛾_𝑡𝑑𝐵_𝑡;with𝜆>0,𝛽<0,𝜎>0,𝛾=0.5. (35)

Note that with appropriate substitution of this model parameters are equivalent to those described in (1) (The equivalence of models (1) and (35) is satisfied by having 𝜆=𝛼𝜇 and 𝛽=−𝛼. Originally in Yu and Phillips [1], the parameter 𝜆 in (35) is defined as 𝛼 but was changed to avoid confusion with the model (1).)

The specific parameters used in each stage are shown in Table 2.

Table 2 

Monte Carlo study parameters.

To perform the first comparison parameters are set as specified in Table 2 and obtained a series simulated. Then you set the value to 𝛾=0.5 and proceed to the estimation by the method of Section 4 to obtain 𝜆, 𝛽, and 𝜎. This procedure was repeated 1000 times (as in Yu and Phillips [1]). Some statistical estimators are presented in Tables 3, 4, and 5.

Table 3 

Comparison  of monthly data series with Nowman, Yu-Phillips, and proposed methods.

Table 4 

Comparison of weekly data series with Nowman, Yu-Phillips, and proposed methods.

Table 5 

Comparison  of daily data series with Nowman, Yu-Phillips, and proposed methods.

Tables 3, 4, and 5 show that the method proposed by Nowman [2] gives good estimates for the parameters 𝜎 and 𝛾, as manifested by Yu and Phillips [1] it becomes clear that a starting point for the method is defined by the latter. Furthermore, through the resulting estimate of the parameter 𝜎, our method becomes much more accurate than Nowman [2], introducing a bias less than 0.2% for the monthly data series and a small value for cases of weekly and daily series.

In the case of 𝜆 and 𝛽 are parameters observed inefficient performance (as discussed Yu and Phillips [1, page 220]) in the method of Nowman, 𝜆 yielding bias of 87%, 47% and 64% for data monthly, weekly and daily, respectively. Also, for 𝛽 the bias present in the estimate is 94%, 46%, and 54% for the respective monthly, weekly, and daily series. The method proposed by Yu and Phillips has a better performance than Nowman, measured from the bias in the parameters, reducing it by 15%, 8%, and 16% for the case of 𝜆 and 5%, 6%, and 6% in the case of 𝛽 for the respective monthly, weekly, and daily series. Our method has a performance similar to that of Yu and Phillips for daily data series, reducing the bias regarding Nowman method in 15% for 𝜆 and 4% for 𝛽. For weekly data series for reducing 𝜆 bias, regarding the method of Nowman, is 7% whereas for 𝛽 is 6%, as reduction Yu and Phillips method.

An important point to note in the results obtained by our method is the condition that the parameter 𝛾 is known a priori. This is necessary in order to find closed formulas for the estimators.

In the approach to the model proposed by Nowman and Yu-Phillips, 𝑑𝑋_𝑡=(𝜆+𝛽𝑋_𝑡)𝑑𝑡+𝜎𝑋^𝛾_𝑡𝑑𝐵_𝑡 is implied to have particular difficulty in estimating the parameter of speed of reversion. To make this condition clear, Table 6 shows the results obtained by the method of Durham and Gallant [4] and our method to estimate the parameters of the model 𝑑𝑋_𝑡=𝛼(𝜇+𝑋_𝑡)𝑑𝑡+𝜎𝑋^𝛾_𝑡𝑑𝐵_𝑡 with 𝛼>0, 𝜇>0, 𝜎>0, 𝛾=0.5. The results are the average of 512 estimations of monthly series with 1000 observations.

Table 6 

Comparison of Durham and Gallant, and proposed methods.

The results shown in Table 6 show that the errors in estimation of 𝜇 and 𝜎 are small (less than 1%) for both methods. On the other hand, the estimated speed of reversion has an error of about 10% resulting in Durham and Gallant, while our method has an error close to 7%. From this analysis, it is clear that the parameter 𝛼 has a higher error in estimating than the other parameter estimates.

In this regard, the estimation error of 𝛼 in Yu [3] provides an analysis of the bias in the parameter estimation inherent reversion total observation time which is sampled in the continuous process, for the special case, where the long-term average 𝜇 is zero. Reference is made to the fact that the bias in parameter depends asymptotically on total observation time which of the process and not the number of observations made in this period.

Figures 1, 2, and 3 show the results of the estimation of the parameters for different observation horizons with different sampling frequencies. (The results are obtained as the average of the individual estimates on 1000 simulated series in 4 sampling scenarios: (a) bimonthly series with 500 observations, (b) monthly series with 1000 observations; (c) weekly series with 4333 observations and (d) daily series with 21,000 observations. The reference model for the simulations was 𝑑𝑋_𝑡=𝛼(𝜇−𝑋_𝑡)𝑑𝑡+𝜎𝑋^𝛾_𝑡𝑑𝐵_𝑡, with 𝛼=1.5, 𝜇=0.5, 𝛾=0.5, and 𝜎=0.4.)

Figure 1 

𝛼 estimation results. 𝑑𝑡 = 1/6, 1/12, 1/52, and 1/252, 𝑛 = 500, 1000, 4333, and 21000, respectively.

Figure 2 

𝜇 estimation results. 𝑑𝑡 = 1/6, 1/12, 1/52, and 1/252, 𝑛 = 500, 1000, 4333, and 21000, respectively.

Figure 3 

𝜎 estimation results. 𝑑𝑡 = 1/6, 1/12, 1/52, and 1/252, 𝑛 = 500, 1000, 4333, and 21000, respectively.

Figure 1 shows the evolution of the estimate of 𝛼 for observations from 10 to 83 years. Behavior of the estimates in the four scenarios can be established empirically by the dependence of the bias on the estimate and the horizon of observation as they have more years of observation bias decreases, with no evidence of a significant difference between the estimates different frequencies.

From Figure 2 it is observed that after few observations, virtually independent of the year of observation, estimations for 𝜇 are obtained with a bias of less than 1%, showing the efficiency of the estimation of the long-term average reversion processes average. In the case of estimation of 𝜎, Figure 3 shows that the bias in the parameter estimates is more dependent on the number of observations which is the total time of observation of the process. This is reflected in the fact that the estimates are in series with a higher sampling rate reaches a lower bias for the same observation horizon than lower sample rate and therefore have fewer observations. On the other hand, as for the case of the estimate of 𝜇, the bias present in the estimates of 𝜎 is practically negligible (less than 1%) in most scenarios.

6. Conclusions and Comments

This paper presents a new alternative for Gaussian estimation models of mean reversion of one factor under the Euler numerical approximation that includes, among others, some models of interest rates in the short term.

Monte Carlo simulations show a good approximation of the diffusion term parameter but less approximation in a trend term parameter. The estimation results have almost the same order as that in Nowman’s [2] and in some cases better than that in Durham and Gallant’s [4]. It should be noted that the proposed methodology outperforms the estimate of 𝜎 compared to the other alternatives presented. These results become more relevant since in practice the estimation of volatility is extremely useful in the evaluation of financial derivatives, in which the trend term loses significance when it is absorbed by the processes of change measurement and valuation risk neutral.

The estimate based on this new method, though does not eliminate the bias in estimation of the reversion rate parameter for a small sample, is in some cases comparable with the results in Yu and Phillips [1], if samples have longer time horizons.

This estimation scheme has several important qualities. First, we come to a closed formula for Gaussian estimation without numerical optimization methods. Second, the computational implementation is simpler and does not require sophisticated routines to perform extensive calculations.