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 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 with initial condition , where , , , and are constants and 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 with .

In general, mean reversion processes are stationary. For example, note that (1) can be written in integral form as

Taking the expected value, and thus we obtain the differential equation

The solution of this equation is , so .

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 , 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 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 , it is said that the process has additive noise so that we get the linear stochastic differential equation

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 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 it is said that the process has proportional noise and gives the nonhomogeneous linear stochastic differential equation given by

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 , we obtain the nonlinear stochastic differential equation

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 and , equation (1) takes the form 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).

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, assuming that volatility stays constant on the intervals ; .

In this way, we obtain an equation of the form where the conditional distribution .

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 where , () satisfies the conditions ; () and This states that the logarithm of the Gaussian likelihood function may be written as where is a vector whose elements can be calculated.

The optimization of is carried out using the methods developed by Bergstrom [2630] 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

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

They introduce a sequence of positive numbers which deliver the required time changes. For any fixed constant , let

Then a succession of points is constructed using the iterations with to obtain where .

Note that according to (16), (12) takes the form

Since the data are observed in discrete time, the time change formula cannot be applied directly and therefore must be approximated by

This should be selected so that is the volatility of the error term unconditional in the equation

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 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 , . 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

In the context of the likelihood, the estimation of is based on the likelihood function where are the Markov transition densities, or and is the joint density function.

Although the differential equation (1) can be solved analytically on , the functions and are not available in closed form, and, therefore, the estimate can not be accomplished, except for the case where .

Note that the exact discretized solution of (1) takes the form which can be approximated by where .

Thus, where , corresponding to the numerical approximation of Euler.

Now define the new variable

Note that , are independent random variables distributed normal, .

In terms of the density function can be written:

Since each is independent can be expressed the joint density function as their marginal densities product, that is,

Consequently, the likelihood function is given by

Taking the natural logarithm on both sides, it follows that

Now consider the problem of maximizing the likelihood function given by

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 where Thus, the estimated parameters are given by

Additionally, the parameter is determined by an equation independent of the system of (31) and is completely determined by the estimated parameters:

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

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.

To perform the first comparison parameters are set as specified in Table 2 and obtained a series simulated. Then you set the value to 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.

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 , , , . The results are the average of 512 estimations of monthly series with 1000 observations.

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 , , , and .)

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.