Hindawi Publishing Corporation

Journal of Probability and Statistics

Volume 2016, Article ID 5191583, 8 pages
http://dx.doi.org/10.1155/2016/5191583

Research Article
Parameter Estimation in Mean Reversion Processes with
Deterministic Long-Term Trend

Freddy H. Marín Sánchez and Verónica M. Gallego

Department of Mathematical Sciences, Universidad Eafit, Carrera 49 No. 7 Sur 50, Medellin, Colombia

Correspondence should be addressed to Freddy H. Marı́n Sánchez; [email protected]

Received 12 April 2016; Accepted 23 July 2016

Academic Editor: Chin-Shang Li

Copyright © 2016 F. H. M. Sánchez and V. M. Gallego. This is an open access article distributed under the Creative Commons
Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is
properly cited.

This paper describes a procedure based on maximum likelihood technique in two phases for estimating the parameters in mean
reversion processes when the long-term trend is defined by a continued deterministic function. Closed formulas for the estimators
that depend on observations of discrete paths and an estimation of the expected value of the process are obtained in the first phase.
In the second phase, a reestimation scheme is proposed when a priori knowledge exists of the long-term trend. Some experimental
results using simulated data sets are graphically illustrated.

1. Introduction data series that help to detect unobservable information from

the process.
In the literature it is common to find studies and applications One of the first difficulties in the application of parameter
of mean reversion models, at both practical and theoretical estimation techniques is that it is necessary to have informa-
levels, especially in fields such as energy markets, commodi- tion with the same time resolution that is specified in the
ties, and interest rates. In these models there is a long-term model. Since models of stochastic differential equations are
trend which acts as an attractor making the process to oscil- formulated in continuous time and that the observed paths
late around it and there is a random component that adds of the process can be obtained only in discrete time, it is nec-
volatility to the movement. The specific characteristics of the essary to discretize the EDE model. An approach frequently
models will depend on the structure of the trend and the used for this purpose is given by Euler-Maruyama. With
volatility of the process. this methodology a new discrete-time process is obtained,
Some common models of mean reversion are those in with which it is possible to make inferences that are valid
which all the parameters in the process are constant, such for continuous time model, as in the case of estimating the
as Ornstein-Uhlenbeck model, the CIR model proposed by parameters.
Cox et al. [1], and in general the models obtained from the There are different procedures that can be implemented
CKLS structure proposed by Chan et al. [2]. In addition to such as methods based on distributional moments [5], Kernel
these models there are other methods not so common in methods [6, 7], least squares, and maximum likelihood
the literature in which the reversion is determined not by a [8–10]. Parameter estimation becomes fundamental too for
parameter but a deterministic function or another stochastic modeling control systems. In the literature it is possible to
process [3, 4]. In general, the fit of the models to specific find different methodologies for those issues such as those
situations is done by assigning values to the parameters and developed for Meng and Ding [11], where the parameter
functions, either from a priori knowledge of the problem or estimation in an equation-error autoregressive (EEAR) sys-
from parameter estimation techniques based on historical tem is done through a transformation of the model to an
2 Journal of Probability and Statistics

equivalent one removing the autoregressive term and to 2. Mean Reversion Processes with
obtain the parameters of the original model the principle of Deterministic Long-Term Trend
equivalence is used. Another methodology is proposed for
Cao and Liu [12] where the parameter estimation in a power Mean reversion processes of one factor with constant param-
system is carried out from a recursive procedure to estimate eters and mean given by a deterministic function can be
simultaneously all the parameters from techniques based on written as
the hierarchical identification principle, using gradient algo- 𝛾
rithms and least squared algorithms. Complementarily, Ding 𝑑𝑋𝑡 = 𝛼 (𝜇 (𝑡) − 𝑋𝑡 ) 𝑑𝑡 + 𝜎𝑋𝑡 𝑑𝐵𝑡 ; 𝑡 ∈ [0, 𝑇] (1)
et al. [13] used methods based on the gradient algorithms and
with initial condition 𝑋0 = 𝑥, where 𝛼 > 0, 𝜎 > 0, and
least squares iterative algorithms for system identification in
𝛾 ∈ [0, 3/2] are constants, 𝜇(𝑡) is a deterministic function,
output error (OE) and output error moving average (OEMA)
and {𝐵𝑡 }𝑡≥0 is a Unidimensional Standard Brownian Motion
systems, where the parameters that depend on unknown
defined on a complete probability space (Ω, F, P).
variables are computed using estimates of these unknown
𝛼 parameter is called reversion rate, 𝜇(𝑡) is the mean
variables from previously estimated parameters.
reversion level or trend of long-run equilibrium, 𝜎 is the
In the area of multirate system identification, Ding et parameter associated with the volatility, and 𝛾 determines the
al. [14, 15] used the polynomial transformation technique to sensitivity of the variance to the level of 𝑋𝑡 .
deal with the identification problem for dual-rate stochastic Model (1) is a generalization of the models CKLS, Chan
systems, while Sahebsara et al. [16] discussed the parameter et al. [2], where the mean reversion level is a deterministic
estimation of multirate multi-input multioutput systems. function that captures the trend of the process. In a general
Also, Ding and Chen [17] presented the combined parameter way, 𝜇(𝑡) plays the role of an attractor at each point 𝑡 in the
and state estimation algorithms of the lifted state-space sense that, when 𝑋𝑡 > 𝜇(𝑡) the trend term 𝛼(𝜇(𝑡) − 𝑋𝑡 ) <
models for general dual-rate systems based on the hierar- 0 and therefore 𝑋𝑡 decreases and when 𝑋𝑡 < 𝜇(𝑡) a similar
chical identification method. Shi et al. [18] gave a crosstalk argument establishes that 𝑋𝑡 grows.
identification algorithm for multirate xDSL FIR systems. To establish the relation between the mean reversion
In general, the auxiliary model identification idea has function 𝜇(𝑡) and the expected value of the process 𝐸[𝑋𝑡 ], we
been used to solve the identification problem of dual-rate can consider (1) written in integral form as
single-input single-output systems as shown in [19].
𝑡 𝑡
For more detailed information about multi-innovation 𝑋𝑡 − 𝑋0 = 𝛼 ∫ (𝜇 (𝑠) − 𝑋𝑠 ) 𝑑𝑠 + 𝜎 ∫ 𝑋𝑠𝛾 𝑑𝐵𝑠 . (2)
stochastic gradient algorithm for multi-input multioutput 0 0
systems and for multirate multi-input systems as well as 𝑡
auxiliary models based multi-innovation extended stochastic Taking the expected value, 𝐸[𝑋𝑡 ] − 𝐸[𝑋0 ] = 𝛼 ∫0 (𝜇(𝑠) −
identification theory (see [20] and references therein). We 𝐸[𝑋𝑠 ])𝑑𝑠 and thus obtain the ordinary differential equation:
will concentrate on one-factor stochastic model in which the
parameters are estimated using maximum likelihood based 𝑚̇ (𝑡) = 𝛼 (𝜇 (𝑡) − 𝑚 (𝑡)) ; 𝑚 (𝑡) = 𝐸 [𝑋𝑡 ] . (3)
on the discretized model when the long-term trend is given
by a deterministic function. In this case we must estimate The solution of this equation is 𝑚(𝑡) = 𝑚(0)𝑒−𝛼𝑡 +
𝑡
both the trend function and the parameters. To estimate the 𝑒 ∫0 𝑒𝛼𝑠 𝜇(𝑠)𝑑𝑠.
−𝛼𝑡

long-term trend function some convolution techniques and To illustrate graphically the relation between 𝜇(𝑡) and
numerical differentiation are used, whereas the normality 𝑚(𝑡), we use 𝜇(𝑡) = 𝑎 + 𝑏 sin(2𝜋𝑡/𝑃) as an example. In this
properties of residuals resulting from discretization are used case the solution to (3) is given by
for estimating parameters by maximum likelihood and in this
way it is possible to obtain closed formulas for estimators 𝛼𝛽𝑏
𝑚 (𝑡) = 𝑎 + [𝑚0 − 𝑎 + ] 𝑒−𝛼𝑡
based on the observations of a path of the process and the 𝛼2 + 𝛽2
estimation of the long-term trend. (4)
𝛼𝑏
In the estimation process there may be some bias in the + 2 [𝛼 sin (𝛽𝑡) − 𝛽 cos (𝛽𝑡)] ,
estimates with respect to the actual values despite the fact 𝛼 + 𝛽2
that they are obtained by the method of maximum likeli- where 𝛽 = 2𝜋/𝑃.
hood. When this situation occurs it is necessary to develop Figure 1 shows the dynamic behavior of 𝜇(𝑡), 𝑚(𝑡) and
alternative methodologies to have a greater adjustment to the one path of 𝑋𝑡 for some specific values in the parameters.
estimates based on the initial estimates [21]. From the Figure it can be seen that the expected value of
This paper is organized as follows. In Section 2 a descrip- the process mimics the long-term mean: that is just the way
tion of mean reversion processes is presented for cases when mean reversion works.
the long-term trend is given by a continued deterministic
function. The first phase of the estimation technique is shown 3. Phase 1: Parameter Estimation
in Section 3. Section 4 presents some examples and prelim-
inary results for the estimation technique and in Section 5 Similar to the procedure established in Marı́n Sánchez and
(second phase), a procedure for reestimating the parameters Palacio [10] and Huzurbazar [22] the main objective is to
from additional a priori knowledge is described. obtain closed formulas for the parameter estimators from
Journal of Probability and Statistics 3

3 Note that 𝑌𝑡 , 𝑖 = 1, 2, 3, . . . , 𝑇, are independent random

2.5 variables distributed 𝑁(0, 𝜎2 Δ).
2 The variable 𝑌𝑡 depends on the observations 𝑋 = (𝑋0 ,
1.5 𝑋1 , . . . , 𝑋𝑇 )󸀠 of the sample path of the process and on the
1 unobservable quantities 𝑚(𝑡) and its derivative 𝑚(𝑡). ̇ Taking
0.5 into account that the sample path resumes all the information
0 that is known from the process, it is necessary to estimate the
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 expected value 𝑚(𝑡) from those observations. As in Tifenbach
Time (t) [4], p. 31-32, we can assume that the expected value of the
process can be approximated using a convolution1 over the
Xt 𝜇(t)
sample path, so the variable 𝑚(𝑡) will depend only on the
m(t)
observations 𝑋 = (𝑋0 , 𝑋1 , . . . , 𝑋𝑇 )󸀠 .
Figure 1: Dynamic behavior for 𝑚(𝑡), 𝜇(𝑡), and 𝑋𝑡 , where 𝑋0 = 2, Thus, we assume that the expected value of the process 𝑋𝑡
𝛼 = 37, 𝜎 = 0.7, 𝛾 = 1, 𝑎 = 1.5, 𝑏 = 1, 𝑃 = 1, and 𝛽 = 2𝜋 with can be approximated by taking a convolution of the sample
Δ𝑡 = 1/250 for a total of 250 observations. path 𝑋 = (𝑋0 , 𝑋1 , . . . , 𝑋𝑇 ) for some particular convolution
function 𝑐(𝑡) given by

discrete observations of a sample path of the process. This 𝑁

can be achieved from the discretized maximum likelihood 𝑚𝑖 = ∑ 𝑐𝑘 𝑋𝑖−𝑘 . (9)
function, which is obtained from the process that results from 𝑘=−𝑁
the Euler-Maruyama discretization scheme applied to (1) and
using the normal and independent increments properties The most common convolution is the moving average, where
of the Brownian Motion over the residuals of the discrete 𝑐𝑘 constants are given by 𝑐𝑘 = 1/(2𝑁 + 1) for all 𝑁, which
process. is denoted by (𝑐𝑡1 ). Other convolutions can be defined varying
Consider the differential equation (1) with the initial value the relative weight of each observation: for example, they
𝛾
𝑋(0) = 𝑥. 𝛼(𝜇(𝑡) − 𝑋t ) and 𝜎𝑋𝑡 are, respectively, the trend can be defined from the coefficients of odd rows of Pascal’s
and diffusion functions that depend on 𝑡, 𝑋𝑡 , and a vector triangle which is denoted by (𝑐𝑡2 ). In this case the central
of unknown parameters 𝜃 = [𝛼, 𝜎, 𝛾]. Now suppose that the observations have a greater weight than more distant obser-
conditions that guarantee the existence and uniqueness of the vations and for this reason the convolution is close to the path.
solution Mao [23] are satisfied, so the expected value of 𝑋𝑡 Those specific cases of convolutions are presented in Figure 2.
exists and 𝜇(𝑡) can be defined using 𝑚(𝑡) and its derivative The convolution is used for the purpose of smoothing the
̇ (see (3)). In this way,
𝑚(𝑡) sample path, eliminating the noise and capturing the trend.
This objective can also be achieved using other techniques
𝑚̇ (𝑡) like filtering. In Figure 2 the Hodrick-Prescott [24] filter
𝜇 (𝑡) = 𝑚 (𝑡) + . (5) hp
𝛼 denoted by (𝑐𝑡 ) was used to obtain an estimation of 𝑚(𝑡) that
Replacing 𝜇(𝑡) in (1) results in are smoother than the convolutions (𝑐𝑡1 ) and (𝑐𝑡2 ).
At this point we have that the expected value 𝑚(𝑡) can be
𝑚̇ (𝑡) 𝛾 approximated from the sample path 𝑋 = (𝑋0 , 𝑋1 , . . . , 𝑋𝑇 ).
𝑑𝑋𝑡 = 𝛼 (𝑚 (𝑡) + − 𝑋𝑡 ) 𝑑𝑡 + 𝜎𝑋𝑡 𝑑𝐵𝑡 ; Now, to get the derivative of 𝑚(𝑡) we can apply numerical
𝛼 (6)
derivatives techniques like Taylor’s theorem: for example, the
𝑡 ∈ [0, 𝑇] . derivative at the observation 𝑖 can be defined using a three-
point rule as
Now we can assume that we have the discrete observa-
tions 𝑋𝑡 = 𝑋(𝜏𝑡 ) for the times {𝜏0 , 𝜏1 , . . . , 𝜏𝑇 }, where Δ 𝑡 = 2𝑚𝑖+1 − 3𝑚𝑖 + 𝑚𝑖−1
𝜏𝑡 − 𝜏𝑡−1 ≥ 0, 𝑡 = 1, 2, . . . , 𝑇. Using the numerical scheme of 𝑚̇ 𝑖 = , (10)
Euler-Maruyama on (6) it follows that Δ

𝑋 (𝜏𝑡 ) = 𝑋 (𝜏𝑡−1 ) for 𝑖 = 1, 2, . . . , 𝑇 − 1. For 𝑖 = 0 the derivate is given by 𝑚̇ 0 =

(𝑚1 − 𝑚0 )/Δ and for 𝑖 = 𝑇 by 𝑚̇ 𝑇 = (𝑚𝑇 − 𝑚𝑇−1 )/Δ. This
+ [𝛼 (𝑚 (𝜏𝑡−1 ) − 𝑋 (𝜏𝑡−1 )) + 𝑚̇ (𝜏𝑡−1 )] Δ 𝑡 (7) procedure can be done using a different number of points to
calculate the derivate at each point.
+ 𝑋𝛾 (𝜏𝑡−1 ) 𝜀𝜏𝑡 , With the estimation of 𝑚(𝑡) and 𝑚(𝑡) ̇ we have one
realization of 𝑌𝑡 and now we can define its normal density
where 𝜀𝜏𝑡 ≈ 𝑁(0, 𝜎2 Δ 𝑡 ). function:
Define the new variable 𝑌𝑡 :
1/2
1
𝑓 (𝑌𝑖 : 𝜃) = ( )
𝑋 − 𝑋𝑡−1 − [𝛼 (𝑚𝑡−1 − 𝑋𝑡−1 ) + 𝑚̇ 𝑡−1 ] Δ 2𝜋𝜎2 Δ
𝑌𝑡 = 𝑡 𝛾 (11)
𝑋𝑡−1 (8) −1 𝑋 − 𝑋𝑖−1 − [𝛼 (𝑚𝑖−1 − 𝑋𝑖−1 ) + 𝑚̇ 𝑖−1 ] Δ
2
⋅ exp [ 2
( 𝑖 𝛾 ) ].
2𝜎 Δ 𝑋𝑖−1
with Δ constant.
4 Journal of Probability and Statistics

3 3
2 2
1 1
0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time (t) Time (t)

Xt E[Xt ] Xt E[Xt ]
ct1 ct2
(a) (b)
3
2
1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time (t)
Xt E[Xt ]
hp
ct
(c)

Figure 2: Dynamic behavior for process 𝑋𝑡 , the expected value 𝑚(𝑡), and three types of convolution. (a) Moving average (𝑐𝑡1 ) with 𝑁 = 20.
hp
(b) Pascal’s weights (𝑐𝑡2 ) with 𝑁 = 20. (c) Hodrick-Prescott filter (𝑐𝑡 ) with 𝜆 hp = 40000. For all the cases Δ𝑡 = 1/250 for a total of 250
observations and the parameters of 𝑋𝑡 are the same as in Figure 1.

Since each 𝑌𝑖 is independent the joint density function ̂ is determined by the equation:
The estimator 𝜎
can be expressed as the product of their marginal densities:
̂
𝜎
𝑇
𝑓 (𝑌1 , 𝑌2 , . . . , 𝑌𝑇 ) = ∏𝑓 (𝑌𝑘 : 𝜃) . (12) 2 (17)
1 𝑇 𝑋𝑖 − 𝑋𝑖−1 − [̂
𝛼 (𝑚𝑖−1 − 𝑋𝑖−1 ) + 𝑚̇ 𝑖−1 ] Δ
𝑘=1 =√ ∑( 𝛾 ).
𝑇Δ 𝑖=1 𝑋𝑖−1
Consequently, the likelihood function is given by
𝑇/2 ̂ we can now approximate
Finally, with the estimation of 𝛼
1 −1
𝐿 (𝜃 | {𝑌𝑡 }) = ( 2
) exp [ 2 ̂ (𝑡) with (5).
𝜇
2𝜋𝜎 Δ 2𝜎 Δ
(13)
𝑇
𝑋 − 𝑋𝑖−1 − [𝛼 (𝑚𝑖−1 − 𝑋𝑖−1 ) + 𝑚̇ 𝑖−1 ] Δ
2 4. Preliminary Results
⋅ ∑( 𝑖 𝛾 ) ].
𝑖=1 𝑋𝑖−1 The aim of this section is to show the performance of the
estimators defined in the previous section. To make this
Taking the natural logarithm on both sides it follows that we will simulate some paths of a process satisfying (1) with
−𝑇 a set of parameters 𝛼 and 𝜎 and giving a deterministic
log 𝐿 = log (2𝜋𝜎2 Δ) form to the 𝜇(𝑡) function and then find the estimators
2
2 (14) with the established procedure in order to determine—using
𝑇
1 𝑋 − 𝑋𝑖−1 − [𝛼 (𝑚𝑖−1 − 𝑋𝑖−1 ) + 𝑚̇ 𝑖−1 ] Δ statistical properties—if the estimates are approximate to
− 2 [∑ ( 𝑖 𝛾 ) ].
2𝜎 Δ 𝑖=1 𝑋𝑖−1 actual parameters.
We analyze the performance of the estimators for sinu-
Now consider the problem of maximizing the likelihood soidal and parabolic 𝜇(𝑡) trend functions and for a set of
function given by values of 𝛼 and 𝜎 when 𝛾 = 1. Each generated path consists
𝜕 log (𝐿) of 250 observations with Δ𝑡 = 1/250. The statistics are
= 0,⃗ ie, 𝜃̂ = arg max (log (𝐿)) . (15) calculated for 1000 different simulations. Considering that for
𝜕𝜃 𝜃 the parameters estimation it is necessary to find the expected
If we assumed that the value of 𝛾 is known a priori, the value of the process from each path, we compare two different
estimation process leads to the estimated parameter 𝛼 given convolution functions and the Hodrick-Prescott filter. For
by numerical derivation a five-point derivation rule is used.
The values of the parameters and 𝜇(𝑡) functions used in
2𝛾
∑𝑇𝑖=1 ((𝑋𝑖 − 𝑋𝑖−1 − 𝑚̇ 𝑖−1 Δ) (𝑚𝑖−1 − 𝑋𝑖−1 ) /𝑋𝑖−1 ) the simulations are presented in Table 1.
̂=
𝛼 𝛾 2
. (16) In Table 2 we present some basic statistics (mean 𝜃 and
∑𝑇𝑖=1 [(𝑚𝑖−1 − 𝑋𝑖−1 ) /𝑋𝑖−1 ] Δ standard deviation 𝑠𝜃 ) for the estimators obtained from
Journal of Probability and Statistics 5

Table 1: Simulation parameters for sinusoidal and parabolic trend (2) Define a functional structure 𝑓(Θ, 𝑡) from a priori
functions. knowledge of the process.
Case 1 Case 2 (3) Estimate the vector parameters Θ minimizing the
𝛼 30 30 error function defined by
𝜎 0.4 0.13
𝑇
𝑋0 2 1.5 Error (Θ) = ∑ [𝑓 (Θ, 𝑡) − 𝜇
2
̂ (𝑡)] . (18)
𝜇(𝑡) 1.5 + sin(2𝜋𝑡) −0.5𝑡2 + 𝑡 + 1.5 𝑡=1

As an example, suppose the functional form of the long-term

the simulations. The convolution functions used in the mean is defined as
simulations are the ones defined for the weights (𝑐𝑘1 ) and
(𝑐𝑘3 ). 𝑐𝑘1 weights are the same defined in Section 3 for the 𝑓 (Θ, 𝑡) = 𝑎𝑡2 + 𝑏𝑡 + 𝑐, (19)
moving average and 𝑐𝑘3 weights are given by 𝑐𝑘,𝑁 3
= (2𝑁 −
|𝑘|)/𝑁(3𝑁 + 1) for 𝑘 from −𝑁 to 𝑁. In this convolution where Θ = [𝑎, 𝑏, 𝑐]󸀠 . In this case the error function would be
the central observation has a relative weight twice the most given by
distant observation.
𝑇
From Table 2 we can conclude that the estimation tech- 2
𝜇 (𝑡) − (𝑎𝑡2 + 𝑏𝑡 + 𝑐)] .
Error (Θ) = ∑ [̂ (20)
nique is working well, especially in the first example (sinu- 𝑡=0
soidal trend). The relative error of the estimators is less than
2% in the case of the parameter 𝜎 but there is a bias in the To estimate Θ it is necessary to solve the system of equations
parameter 𝛼 that depends on the convolution function and at defined by ∇ Error(Θ) = 0, which is equivalent to the system
the same time on the parameter 𝑁 used in such convolution. of linear equations 𝑀Θ = 𝑑, where
In Figure 3 we can see the performance of the estimation
of the deterministic function 𝜇(𝑡). In that figure an estimation ∑𝑡4 ∑𝑡3 ∑𝑡2
from a random path is shown instead of the average. The [ ]
[ ]
estimates are close to the theoretical function in both cases [∑𝑡3 ∑𝑡2 ∑𝑡 ]
𝑀=[ ]
but having some differences at the beginning and the final of [ ]
[ 2 ]
the observations. One form to measure the fit between the ∑𝑡 ∑𝑡 (𝑇 + 1)
theoretical function and its estimation is using the Root Mean [ ]
(21)
Square error (RMS). Table 3 presents that measure for the
proposed estimation in cases 1 and 2. ∑𝑡2 𝜇
̂ (𝑡)
[ ]
[ ]
𝜇 (𝑡) ]
[ ∑𝑡̂
𝑑=[ ].
5. Phase 2: Re-Estimation [ ]
[ ]
From the simulations we can see that the estimate of 𝜇(𝑡) 𝜇 (𝑡)
∑̂
[ ]
is not very accurate and that differences come from 𝑚 ̂ (𝑡)
̂̇
and 𝑚(𝑡), and in turn this affects the estimation of other A similar procedure can be implemented under the assump-
parameters. tion that 𝑓(Θ, 𝑡) = 𝑎 + 𝑏 sin((2𝜋/𝑃)𝑡), to find Θ = [𝑎, 𝑏]󸀠 if 𝑃
If it were possible to have all the paths of the process it is known. A comparison between 𝜇(𝑡), 𝜇est (𝑡), and 𝜇rest (𝑡) =
would not be necessary for numerical approximation of the ̂ 𝑡) is shown in Figure 4 for both cases and for a random
𝑓(Θ,
expected value or its derivative to capture information from path of the process defined with the same set of parameters
𝜇(𝑡), and it could be enough to perform an average of all of Table 1.
paths at each time 𝑡𝑖 and that could generate pretty accurate Having a best estimate of 𝜇(𝑡) allows a reestimation of the
estimations. As this condition is not the usual but only has parameters 𝛼 and 𝜎. To achieve this we must perform a similar
a path, the procedure detailed above is necessary and this procedure to that proposed in Section 4, defining
involves the errors mentioned in the estimates of 𝑚(𝑡) and
̇
𝑚(𝑡). 𝑋𝑡 − 𝑋𝑡−1 − 𝛼 (𝜇𝑡−1 − 𝑋𝑡−1 ) Δ
𝑌𝑡 = 𝛾 Δ constant. (22)
To correct this deficiency it is proposed to carry out 𝑋𝑡−1
a reestimate procedure trying to adjust a function to the
initial estimate from a priori knowledge of the process. This So,
recalculation allows a correction in the other estimators while
providing a functional form that allows simulations of future ̂̂ − 𝑋 ) Δ/𝑋2𝛾 )
∑𝑇𝑖=1 ((𝑋𝑖 − 𝑋𝑖−1 ) (𝜇 𝑖−1 𝑖−1 𝑖−1
̂̂ =
𝛼 ,
paths of the process. ̂̂ − 𝑋 ) /𝑋𝛾 ]2 Δ
The following procedure is proposed to perform the ∑𝑇𝑖=1 [(𝜇 𝑖−1 𝑖−1 𝑖−1
reestimate of 𝜇(𝑡): (23)
2
𝑇 ̂̂ (𝜇
𝑋𝑖 − 𝑋𝑖−1 − 𝛼 ̂̂ − 𝑋 ) Δ
(1) Get an estimate of 𝜇̂ (𝑡) according to the procedure ̂̂ = √ 1 ∑ (
𝜎
𝑖−1 𝑖−1
).
𝛾
described in Section 4. 𝑇Δ 𝑖=1 𝑋𝑖−1
6 Journal of Probability and Statistics

Table 2: Results obtained with 1000 simulations with 250 observations each. The real values of the parameters are shown in Table 1. For case
hp
1 in the convolutions we use 𝑁 = 22 and for case 2 𝑁 = 60. For both cases we set 𝜆 = 62500 in 𝑐𝑘 .

Case 1 Case 2
hp hp
𝜇1 (𝑡) 𝑐𝑘1 𝑐𝑘3 𝑐𝑘 𝜇2 (𝑡) 𝑐𝑘1 𝑐𝑘3 𝑐𝑘
𝛼 24.8540 31.3347 32.8771 𝛼 31.1160 35.8249 51.1093
𝛼 𝛼
𝑠𝛼 3.2218 3.9947 5.1406 𝑠𝛼 5.8254 6.6984 10.1414
𝜎 0.4002 0.3942 0.4023 𝜎 0.1293 0.1281 0.1267
𝜎 𝜎
𝑠𝜎 0.0177 0.0174 0.0174 𝑠𝜎 0.0060 0.0059 0.0060

3 2

2 1.8

1 1.6

0 1.4
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

𝜇 𝜇
𝜇est 𝜇est
(A) (A)

3 2

2 1.8

1 1.6

0 1.4
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

𝜇 𝜇
𝜇est 𝜇est

(B) (B)
3 2

2 1.8

1 1.6

0 1.4
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

𝜇 𝜇
𝜇est 𝜇est

(C) (C)
(a) (b)

Figure 3: Estimated 𝜇1 (𝑡) and 𝜇2 (𝑡). The left side correspond to Case 1 (sinusoidal) and the right side to Case 2 (parabolic). The (A) plots are
hp
obtained with the convolution (𝑐𝑘1 ), (B) plots with (𝑐𝑘3 ) and (C) plots with (𝑐𝑘 ). The parameters involved in this figure are presented in Tables
1 and 2.
Journal of Probability and Statistics 7

3 2

2 1.8

1 1.6

0 1.4
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Time (t) Time (t)

𝜇 𝜇rest 𝜇 𝜇rest
𝜇est 𝜇est
(a) (b)

Figure 4: Comparison of dynamic behavior of 𝜇(𝑡), 𝜇est (𝑡), and 𝜇rest (𝑡). In case (a) 𝜇1 (𝑡) = 1.5+sin(2𝜋𝑡) and in case (b) 𝜇2 (𝑡) = −0.5𝑡2 +𝑡+1.5.

Table 3: Results of RMS error for each estimation 𝜇1 (𝑡) and 𝜇2 (𝑡). Table 4: Estimated and reestimated parameters. The estimations are
hp
calculated with the convolution (𝑐𝑘 ). Results obtained with 1000
RMS Case 1 Case 2 simulations with 250 observations each for 𝜇1 and 𝜇2 as defined in
𝑐𝑘1 0.0782 0.0266 Table 1: case 1 𝛼 = 30 and 𝜎 = 0.4, case 2 𝛼 = 30 and 𝜎 = 0.13.
𝑐𝑘3 0.0578 0.0238 ̂̂ ̂̂
hp 𝜇1 (𝑡) ̂
𝛼 𝛼 𝜇2 (𝑡) ̂
𝛼 𝛼
𝑐𝑘 0.0767 0.0140
𝛼 24.8540 29.6620 𝛼 31.1160 30.4616
𝛼 𝛼
𝑠𝛼 3.2218 3.4518 𝑠𝛼 5.8254 5.6828
𝜇1 (𝑡) ̂
𝜎 𝜎̂̂ 𝜇2 (𝑡) ̂
𝜎 𝜎̂̂
Table 4 summarizes the results obtained from the first
estimation and the reestimation procedure for the cases 𝜎 0.4002 0.4021 𝜎 0.1293 0.1300
𝜎 𝜎
where 𝜇(𝑡) has a sinusoidal form and a parabolic form and 𝑠𝜎 0.0177 0.0187 𝑠𝜎 0.0060 0.0060
the other parameters are defined in Table 1.
From the results of the above example we can observe Table 5: Results of RMS error for the estimation and reestimation
how the reestimation procedure reduces the bias in the initial of 𝜇(𝑡): case 1 sinusoidal trend and case 2 parabolic trend.
estimation of 𝛼 for both sinusoidal and parabolic cases. The
final bias in 𝛼 is approximately 1.5% of the real value of RMS Case 1 Case 2
the parameter, reducing almost 90% of the error in the first 𝜇est (𝑡) 0.0961 0.0137
estimation in the case of sinusoidal function and near 60% 𝜇rest (𝑡) 0.0261 0.0058
for the parabolic function. The reestimations of 𝜎 do not
significantly change from the original estimate. For the case
of 𝜇(𝑡) the RMS error is present in Table 5 for a random path because it allows us to adjust the initial estimates of the
of the process. In both cases the error in the reestimation is parameters making more reliable the results obtained with
less than the 50% of the error present in the initial estimation. the adjusted models.
On the other hand, in the development of the proposed
6. Conclusions and Comments methodology it can be seen that it allows us to obtain closed
formulas for the maximum likelihood estimators in mean
This paper shows how starting from a path of a stochastic reversion process when the long-term trend is defined by a
process is possible to find an approximation of the parameters deterministic function, decreasing computational effort and
that define the underlying dynamics It is important to allowing greater understanding of procedures.
highlight that, for modeling the dynamic behavior of prices
of some commodities or interest rates, only one realization Competing Interests
of the process in discrete time is available and it represents
the only source of information to capture the nature of the The authors declare that they have no competing interests.
process. Therefore, it becomes very important to obtain and
interpret the invisible and visible information found on the Endnotes
path to make projections more adjusted to the future process
behavior being studied. 1. A convolution function 𝑐(𝑡) is an even, continuous, and
It is essential that a model that seeks to reproduce the real valued function satisfying the fact that 𝑐(𝑡) ⩾ 0 and
∞
behavior of a process that is generated in the real world ∫−∞ 𝑐(𝑡)𝑑𝑡 = 1. A convolution of a curve 𝑠(𝑡) is another
incorporates a priori information provided by an expert in the function 𝑙(𝑡) defined by
topic and that this additional information allows the model
to better fit the process. This prior information is key in ∞

the development of the methodology proposed in this paper, 𝑙 (𝑡) = ∫ 𝑐 (𝑤) 𝑠 (𝑡 − 𝑤) 𝑑𝑤, (∗)
−∞
8 Journal of Probability and Statistics

where 𝑐(𝑡) is a convolution function. [17] F. Ding and T. Chen, “Hierarchical identification of lifted state-
In the discrete case the convolution of two sequences can space models for general dual-rate systems,” IEEE Transactions
be written as on Circuits and Systems. I. Regular Papers, vol. 52, no. 6, pp. 1179–
1187, 2005.
∞
[18] Y. Shi, F. Ding, and T. Chen, “Multirate crosstalk identification
𝑙𝑖 = ∑ 𝑐𝑘 𝑠𝑖−𝑘 . (∗∗) in xDSL systems,” IEEE Transactions on Communications, vol.
𝑘=−∞ 54, no. 10, pp. 1878–1886, 2006.
[19] F. Ding and T. Chen, “Combined parameter and output estima-
References tion of dual-rate systems using an auxiliary model,” Automatica,
vol. 40, no. 10, pp. 1739–1748, 2004.
[1] J. C. Cox, J. Ingersoll, and S. Ross, “A theory of the term structure [20] L. Han, J. Sheng, F. Ding, and Y. Shi, “Auxiliary model identifi-
of interest rates,” Econometrica, vol. 53, no. 2, pp. 385–407, 1985. cation method for multirate multi-input systems based on least
[2] K. C. Chan, G. A. Karolyi, F. A. Longstaff, and A. B. Sanders, “An squares,” Mathematical and Computer Modelling, vol. 50, no. 7-
empirical comparison of alternative models of the short-term 8, pp. 1100–1106, 2009.
interest rate,” The Journal of Finance, vol. 47, no. 3, pp. 1209–1227, [21] J. Yu, “Bias in the estimation of the mean reversion parameter in
1992. continuous time models,” Journal of Econometrics, vol. 169, no.
[3] D. Pilipovic, Energy Risk: Valuing and Managing Energy Deriva- 1, pp. 114–122, 2012.
tives, McGraw-Hill, New York, NY, USA, 2007. [22] V. S. Huzurbazar, “The likelihood equation, consistency and the
[4] B. Tifenbach, Numerical Methods for Modeling Energy Spot maxima of the likelihood function,” Annals of Eugenics, vol. 14,
Prices, University of Calgary, 2000. pp. 185–200, 1948.
[5] T. Koulis and A. Thavaneswaran, “Inference for interest rate [23] X. Mao, Stochastic Differential Equation and Applications, Hor-
models using Milstein’s approximation,” Journal of Mathemat- wood Publishing Limited, England, UK, 1997.
ical Finance, vol. 3, no. 1, pp. 110–118, 2013. [24] R. J. Hodrick and E. C. Prescott, “Postwar U.S business cycles: an
[6] D. Florens-Zmirou, “On estimating the diffusion coefficient empirical investigation,” Journal of Money, Credit and Banking,
from discrete observations,” Journal of Applied Probability, vol. vol. 29, no. 1, pp. 1–16, 1997.
30, no. 4, pp. 790–804, 1993.
[7] G. J. Jiang and J. L. Knight, “A nonparametric approach to the
estimation of diffusion processes, with an application to a short-
term interest rate model,” Econometric Theory, vol. 13, no. 5, pp.
615–645, 1997.
[8] K. B. Nowman, “Gaussian estimation of single-factor contin-
uous time models of the term structure of interest rates,” The
Journal of Finance, vol. 52, no. 4, pp. 1695–1706, 1997.
[9] J. Yu and P. Phillips, “Corrigendum to ‘a Gaussian approach
for continuous time models of the shortterm interest rate’,”
Econometrics Journal, vol. 14, pp. 126–129, 2011.
[10] F. H. Marı́n Sánchez and J. S. Palacio, “Gaussian estimation of
one-factor mean reversion processes,” Journal of Probability and
Statistics, vol. 2013, Article ID 239384, 10 pages, 2013.
[11] D. Meng and F. Ding, “Model equivalence-based identification
algorithm for equation-error systems with colored noise,” Algo-
rithms, vol. 8, no. 2, pp. 280–291, 2015.
[12] Y. Cao and Z. Liu, “Signal frequency and parameter estimation
for power systems using the hierarchical identification princi-
ple,” Mathematical and Computer Modelling, vol. 52, no. 5-6, pp.
854–861, 2010.
[13] F. Ding, P. X. Liu, and G. Liu, “Gradient based and least-
squares based iterative identification methods for OE and
OEMA systems,” Digital Signal Processing, vol. 20, no. 3, pp.
664–677, 2010.
[14] F. Ding, P. X. Liu, and Y. Shi, “Convergence analysis of
estimation algorithms for dual-rate stochastic systems,” Applied
Mathematics and Computation, vol. 176, no. 1, pp. 245–261, 2006.
[15] F. Ding, P. X. Liu, and H. Z. Yang, “Parameter identification
and intersample output estimation for dual-rate systems,” IEEE
Transactions on Systems, Man, and Cybernetics, Part A: Systems
and Humans, vol. 38, no. 4, pp. 966–975, 2008.
[16] M. Sahebsara, T. Chen, and S. L. Shah, “Frequency-domain
parameter estimation of general multi-rate systems,” Computers
and Chemical Engineering, vol. 30, no. 5, pp. 838–849, 2006.
Advances in Advances in Journal of Journal of
Operations Research
Hindawi Publishing Corporation
Decision Sciences
Hindawi Publishing Corporation
Applied Mathematics
Hindawi Publishing Corporation
Algebra
Hindawi Publishing Corporation
Probability and Statistics
Hindawi Publishing Corporation
http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014

The Scientific International Journal of

World Journal
Hindawi Publishing Corporation
Differential Equations
Hindawi Publishing Corporation
http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014

Submit your manuscripts at

http://www.hindawi.com

International Journal of Advances in

Combinatorics
Hindawi Publishing Corporation
Mathematical Physics
Hindawi Publishing Corporation
http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014

Journal of Journal of Mathematical Problems Abstract and Discrete Dynamics in

Complex Analysis
Hindawi Publishing Corporation
Mathematics
Hindawi Publishing Corporation
in Engineering
Hindawi Publishing Corporation
Applied Analysis
Hindawi Publishing Corporation
Nature and Society
Hindawi Publishing Corporation
http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014

International
Journal of Journal of
Mathematics and
Mathematical
Discrete Mathematics
Sciences

Journal of International Journal of Journal of

Hindawi Publishing Corporation Hindawi Publishing Corporation Volume 2014

Function Spaces
Hindawi Publishing Corporation
Stochastic Analysis
Hindawi Publishing Corporation
Optimization
Hindawi Publishing Corporation
http://www.hindawi.com Volume 2014 http://www.hindawi.com http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014