Estimation of all parameters in the reflected Orntein-Uhlenbeck process from discrete observations
aa r X i v : . [ m a t h . S T ] S e p Estimation of all parameters in the reflectedOrntein-Uhlenbeck process from discrete observations
Yaozhong Hu a , Yuejuan Xi b, ∗ a Department of Mathematical and Statistical Sciences, University of Alberta at EdmontonEdmonton, Alberta Canada, T6G 2G1 b School of Mathematical Sciences, Nankai University, Tianjin, PR China, 300071
Abstract
Assuming that a reflected Ornstein-Uhlenbeck state process is observed at dis-crete time instants, we propose generalized moment estimators to estimate alldrift and diffusion parameters via the celebrated ergodic theorem. With thesampling time step h > n tends to infinity.This provides a complete solution to an open problem left in Hu et al. [5]. Keywords:
Reflected Ornstein-Uhlenbeck process; Ergodic theorem; Spectralrepresentation of transition density; Strong consistency; Asymptotic normality.
1. Introduction
On a filtered probability space (Ω , P , F , {F t } t ≥ ) let W = { W ( t ) } t ≥ bea one-dimensional standard Brownian motion. All the processes mentioned inthis paper will be adapted to {F t } t ≥ . We consider the (reflected) Ornstein-Uhlenbeck (ROU) process, reflected at zero, which is defined by the following ∗ Corresponding author
Email addresses: [email protected] (Yaozhong Hu), [email protected] (Yuejuan Xi) ne-dimensional stochastic differential equation (SDE): dX t = κ ( θ − X t ) dt + σdW t + L t , t ∈ R + = { x, x ≥ } ,X = x ∈ R + , (1.1)where κ, θ, σ ∈ (0 , ∞ ) are constants and L t is the minimal continuous increasingprocess which ensures that X t ≥ t ≥
0. The ROU process is a usefulstochastic model in finance and queue theory (cf. Linetsky [7], Ward and Glynn[8] and the references therein).This paper will concern with the statistical estimation problem for the pa-rameters κ, θ, σ from the observations. In most practical situations the obser-vations of the process { X t , t ≥ } can be made only at discrete time instants t k = kh , k = 1 , , · · · , and usually the time interval h between consecutiveobservations cannot be made arbitrarily small. To deal with this situation anergodic type of estimator to estimate κ and θ is proposed in a previous work Huet al. [5] and the strong consistency and asymptotic normality of the estimatorsare also obtained there. However, as pointed out in Hu et al. [5] they wereunable to estimate σ (or σ ) by using the ergodic type estimator and insteadthey proposed to use ˆ σ c,n := nh P nk =1 ( X ( k +1) h − X kh ) as the estimator of σ . Let us also mention a work on the estimation of the parameters κ, θ, σ forthis ROU when continuous observation is available (Bo et al. [2]). This wouldrequire that h → σ c,n → σ ). This paper will fill this gap. We shall introduce an ergodic typeestimator to estimate σ (and hence we can estimate all the parameters κ, θ, σ simultaneously) and prove the strong consistency and asymptotic normality forall estimators (including the estimator ˆ σ n for σ ) regardless the (fixed) valueof h . This work is motivated by a recent work of Cheng et al. [4], where theOrnstein-Uhlenbeck process has no reflection, but the Brownian motion wasreplaced by a stable process.Now let us describe our ergodic estimators for all parameters. It is well-known that there is a unique invariant probability density function π ( x ) of X t f we havelim n →∞ n n X k =1 f ( X t k ) = Z R + f ( x ) π ( x ) dx, (1.2)and the invariant probability density function π ( x ) has the following explicitexpression (see Hu et al. [5]): π ( x ) = √ κσ φ ( √ κ ( x − θ ) σ )1 − Φ( − √ κθσ ) , (1.3) φ ( x ) = e − | x | / √ π is the standard normal probability density function, andΦ( x ) = R x −∞ φ ( u ) du is the standard normal distribution function. As observedin Hu et al. [5] the invariant measure π ( x ) remains the same function if thequantities θ and κσ remain unchanged. Thus, we cannot expect to use (1.2) toestimate κ and σ simultaneously. To this end and motivated by Cheng et al.[4] we shall use the ergodic theorem for X kh , X ( k +1) h , which states that for anyintegrable function f : R → R ,lim n →∞ n n X k =1 f ( X kh , X ( k +1) h ) = E f ( ˜ X , ˜ X h ) = Z R f ( x, y ) π ( x ) p h ( x, y ) dxdy, (1.4)where ˜ X is a random variable independent of the Brownian motion W andhaving the invariant probability density π ( x ), ˜ X is the solution to (1.1) withinitial random variable ˜ X , and p h ( x, y ) is the transition density of X .With some specific choices of f in (1.2) and (1.4) we can obtain our ergodicestimators, whose detailed construction is given in the next section, where thestrong consistency and asymptotic normality are also obtained.Section 3 will provide a numerical example which demonstrates the conver-gence results of our estimators and which also demonstrates that ˆ σ c,n does notconverge.
2. Strong consistency and asymptotic normality
In this section, we aim to construct the estimators for all the parameters κ, θ, σ of the ROU process { X t , t ≥ } given by (1.1) based on discrete obser-vations { X t , , · · · , X t n } , where t k = kh with the observation time interval h Lemma 2.1.
The h -skeleton sampled chain { X kh : k ≥ } is ergodic. Namely,for any initial value x ∈ R + and f ∈ L ( R + ) and g ∈ L ( R ) we have lim N →∞ n n X k =1 f ( X kh ) = E [ f ( X ∞ )] = Z R + f ( x ) π ( x ) dx, a.s., lim n →∞ n n X k =1 g ( X kh , X ( k +1) h ) = E g ( ˜ X , ˜ X h )= Z R g ( x, y ) π ( x ) p h ( x, y ) dxdy, a.s., (2.1)(2.2) where ˜ X is a random variable independent of the Brownian motion W andhaving the invariant probability density π , ˜ X is the solution to (1.1) with initialrandom variable ˜ X , and p h ( x, y ) is the transition density of X . As illustrated in Hu et al. [5] it is impossible to use (2.1) alone to estimate allthe parameters κ, θ, σ . So we take f ( x ) = x and f ( x ) = x in (2.1) and we take g ( x, y ) = xy in (2.2) to obtain a system of three equations to determine the pa-rameters κ, θ, σ . Some elementary computations yield the following expressionsfor the stationary moments of the invariant measure. E [ X ∞ ] = θ + σ √ κ φ ( √ κθσ )1 − Φ( − √ κθσ ) , E [ X ∞ ] = σ κ + θ + θ σ √ κ φ ( √ κθσ )1 − Φ( − √ κθσ ) , E ( ˜ X ˜ X h ) = Z R xyπ ( x ) p h ( x, y ) dxdy . (2.3)However, to our best knowledge, there is no compact explicit form for the tran-sition probability density p t ( x, y ). We shall use the following spectral represen-tation for the transition density p t ( x, y ) derived in Linetsky [7] p t ( x, y ) = π ( y ) + m ( y ) ∞ X i =1 e − λ i t ϕ i ( x ) ϕ i ( y ) , t > , where the notations are described as follows.41) π ( x ) is the stationary density given by (1.3) and m ( x ) is the speed measuredefined by m ( x ) = 2 σ e − κ ( θ − x ) /σ . (2) The eigenvalues 0 < λ < λ < · · · < λ i < · · · are roots of H λ/κ − ( −√ κθ/σ ) = 0 , where H is the Hermite function (see Lebedev [6]).(3) The normalized eigenfunctions ϕ i ( x ) are given by ϕ i ( x ) = ± κ / σ / e κθ / (2 σ ) H λ i /κ ( √ κ ( x − θ ) /σ ) q λ i △ i H λ i /κ ( −√ κθ/σ ) , i = 1 , , · · · where △ i = ∂H ν − ( −√ κθ/σ ) ∂ν | ν = λ i /κ .Now we replace E ( X ∞ ) , E ( X ∞ ) , E ( ˜ X ˜ X h ) in (2.3) by their sample approxima-tions to yield n n X k =1 X kh = θ + σ √ κ φ ( √ κθσ )1 − Φ( − √ κθσ ) , n n X k =1 X kh = σ κ + θ + θ σ √ κ φ ( √ κθσ )1 − Φ( − √ κθσ ) , n n X k =1 X kh X ( k +1) h = Z R xyπ ( x ) p h ( x, y ) dxdy. (2.4)This is a system of three equations for the three unknown parameters. Weexpect that it would give a unique solution ˆ κ n , ˆ θ n , ˆ σ n , which we call the ergodicestimators of the parameters. The system is still complicated to analyze andto be solved. We will further simplify it. To this end we denote u = θ and v = √ κθσ . Then the first two equations in (2.4) depends only on u and v andthey give a unique solution ˆ u n and ˆ v n . We then write p h ( x, y ) = p h ( x, y ; κ, θ, σ )as a kernel p h ( x, y ; u, v, σ ) depending on parameters u, v, σ . Finally, we replacethe parameters u and v in kernel p h ( x, y ; u, v, σ ) by the obtained values ˆ u n andˆ v n , then the third equation in (2.4) becomes one equation for one unknown σ .This greatly simplifies the computations.5o summarize the above discussion, we have transformed the system (2.4)into the following system of equations. n n X k =1 X kh = u + uv φ ( v )1 − Φ( − v ) , n n X k =1 X kh = u v + u + u v φ ( v )1 − Φ( − v ) , n n X k =1 X kh X ( k +1) h = Z ∞ Z ∞ xyp h ( x, y ) π ( x ) dxdy . (2.5)The right-hand side of the above third equation depends on u , v and σ bysubstituting θ and κ by θ = u and κ = v σ u into the expression of π and p h ( x, y ). Define ˜ λ i = λ i /σ . For sake of the numerical computation we writethe dependence explicitly as follows: p h ( x, y ) = π ( y ) + m ( y ) ∞ X i =1 e − ˜ λ i σ h ϕ i ( x ) ϕ i ( y ) , (2.6)where(1) π and m are given by π ( x ) = vu φ (cid:16) vu x − v (cid:17) / [1 − Φ( − v )] , (2.7)and m ( x ) = 2 σ e − v / v x/u − v x / (2 u ) , (2.8)(2) The eigenvalues 0 < ˜ λ < ˜ λ < · · · < ˜ λ i < · · · are roots of H u ˜ λ/v − ( − v/ √
2) = 0 . (2.9)(3) The eigenfunctions are given by ϕ i ( x ) = ± σ ( v/ ( √ u )) / e v / H u ˜ λ i /v (( vu x − v ) / √ q λ i △ i H u ˜ λ i /v ( − v/ √ , i = 1 , , · · · (2.10)with △ i = ∂H ν − ( − v/ √ ∂ν | ν =2 u ˜ λ i /v .Now we summarize our discussion as follows. Construction of the ergodic estimators for all parameters κ, θ, σ :6i) Solve the first two equations in the system (2.5) to obtain ˆ u n and ˆ v n .(ii) Substitute the obtained ˆ u n and ˆ v n into the transition probability kernel p h ( x, y ) according to (2.6)-(2.10) to obtain the third equation in the system(2.5), which now contains only one unknown σ and solve it to obtain ˆ σ n .(iii) Solve ˆ u n = ˆ θ n and ˆ v n = √ κ n ˆ θ n ˆ σ n to obtainˆ θ n = ˆ u n and ˆ κ n = ˆ v n ˆ σ n θ n . (2.11) Remark 2.1.
In numerical computation, we shall need to take finite terms inthe spectral representation of the transition probability function (in our numer-ical simulation we take about twelve terms and the results are satisfactory).The Hermite functions and the roots of the Hermite functions can be handledby the standard mathematical software package. The system (2.5) of algebraicequations does not give an explicit solution. There are many standard methodsto solve it, such as the Newton-Raphson iteration method.To study the strong consistency and the asymptotic normality, we denote theright-hand sides of the equation in the system (2.5) by g ( u, v ), g ( u, v ), and g ( u, v, σ ), respectively. Denote M ,n = 1 n n X k =1 X kh , M ,n = 1 n n X k =1 ( X kh ) , M ,n = 1 n n X k =1 X kh X ( k +1) h . Then the equation (2.5) can be rewritten as g ( u, v ) = M ,n , g ( u, v ) = M ,n , g ( u, v, σ ) = M ,n . (2.12)( g also depends on h which is fixed) Or we write g ( u, v, σ ) = M n , (2.13)where g = ( g , g , g ) T and M n = ( M ,n , M ,n , M ,n ) T . Denote by J ( u, v, σ ) the determinant of the Jacobian of g . Then J ( u, v, σ ) = J ( g , g ) ∂∂σ g ( u, v, σ ) , (2.14)7here J ( g , g ) is the determinant of the Jacobian of g ad g . Hu et al. [5]proved that J ( g , g ) is never 0. If ∂∂σ g ( u, v, σ ) is not singular in some domain D ⊆ R , then by the inverse function theorem, for any ( u, v, σ ) ∈ D , g =( g , g , g ) has a unique inverse in a neighbourhood ( u, v, σ ). If ( u, v, σ ) (orequivalently, ( κ, θ, σ )) are the true parameters, then by Lemma 2.1, we seewhen n is sufficiently large ( M ,n , M ,n , M ,n ) will be in the neighbourhood of g ( u, v, σ ). This means when n is sufficiently large the equation (2.12) has asolution.Thus, the critical question now is to find a domain D such that ∂∂σ g ( u, v, σ )is not singular on D . This is an elementary analysis problem. The explicitexpression of the derivative of g ( u, v, σ ) with respect to σ can be obtained (seeRemark 2.1 below). However, this expression is complicated and it is hard toobtain the domain of ( u, v, σ ) so that inside this domain this derivative is notsingular. We shall proceed as follows to reduce the ∂∂σ g ( u, v, σ ) from a functionof three variables u, v, σ to a function of one variable σ .Since the first two equations in (2.5) is independent of σ , as indicated abovewe can solve them without considering the third equation in (2.5). Hu et al. [5]proved that there exist continuous inverse mapping ( h , h ) of ( g , g ) : R → R such that the ergodic estimators defined byˆ u n := h ( M ,n , M ,n ) , ˆ v n := h ( M ,n , M ,n ) (2.15)converge almost surely to the true parameters u = h ( g ( u, v ) , g ( u, v )) = θ, v = h ( g ( u, v ) , g ( u, v )) = √ κθσ . After the estimators ˆ u n and ˆ v n have been obtained, we can substitute theminto the g . Thus g ( σ ) = g (ˆ u n , ˆ v n , σ ) and g ′ ( σ ) = ∂∂σ g (ˆ u n , ˆ v n , σ ) will befunctions of single variable σ . We can plot the derivative function g ′ ( σ ) inan interval D σ that is as large as we believe it contains the true parameter σ (we shall plot g ′ ( σ ) for some value of u and v in next section). If g ′ isnever equal to 0 on D σ , then the solution to g (ˆ u n , ˆ v n , σ ) = M ,n is uniqueon D σ (if not then there are two different points σ < σ in D σ such that8 (ˆ u n , ˆ v n , σ ) = g (ˆ u n , ˆ v n , σ ) = M ,n . By the mean value theorem there is a σ ∈ [ σ , σ ] ⊆ D σ such that ∂∂σ g (ˆ u n , ˆ v n , σ ) = 0).If g ′ is not singular on D σ , then the third equation (2.12) has a uniquesolution on D σ , which gives the ergodic estimator ˆ σ n = h (ˆ u n , ˆ v n , M ,n ) of σ ,where h (ˆ u n , ˆ v n , · ) is the continuous inverse of g (ˆ u n , ˆ v n , · ). By Lemma 2.1 it iseasy to see that ˆ σ n → σ a.s.Now we summarize the above discussion as the following theorem. Theorem 2.1. (i) The first two equations of the system (2.5) have a uniquesolution pair (ˆ u n , ˆ v n ) = ( h ( M ,n , M ,n ) , h ( M ,n , M ,n )) .(ii) If ∂∂σ g (ˆ u n , ˆ v n , σ ) ( g is defined by the right-hand side of the third equationin (2.5) ) is not singular on some interval σ ∈ D σ which contains the trueparameter σ , then when n is sufficiently large the third equation of (2.5) ,namely, g (ˆ u n , ˆ v n , σ ) = M ,n = n X k =1 X kh X ( k +1) h (2.16) has a unique solution ˆ σ n .(iii) (ˆ κ n , ˆ θ n , ˆ σ n ) T → ( κ, θ, σ ) T almost surely as n → ∞ , where ˆ θ n and ˆ κ n aregiven by (2.11) . Next, we study the joint asymptotic behavior of the all estimators (ˆ θ n , ˆ κ n , ˆ σ n ). Theorem 2.2.
Let ∂∂σ g (ˆ u n , ˆ v n , σ ) be nonsingular on some interval σ ∈ D σ which contains the true parameter σ . Then, the estimators (ˆ θ n , ˆ κ n , ˆ σ n ) satisfythe following asymptotic normality property: √ n ((ˆ θ n , ˆ κ n , ˆ σ n ) T − ( θ, κ, σ ) T ) ⇒ N (0 , Σ) , where Σ is a covariance matrix defined in (2.17) below.Proof. For any nice function f and g , denote σ fg :=Cov( f ( ˜ X , ˜ X h ) , g ( ˜ X , ˜ X h )) + ∞ X k =1 [Cov( f ( ˜ X , ˜ X h ) , g ( ˜ X kh , ˜ X ( K +1) h ))+ Cov( g ( ˜ X , ˜ X h ) , f ( ˜ X kh , ˜ X ( k +1) h ))] . f ( x, y ) = x , f ( x, y ) = y , f ( x, y ) = xy and denote˜Σ := ( σ f k f l ) ≤ k,l ≤ . Then an application of the multivariate Markov chain central limit theorem (e.g.Brooks et al. [3, Section 1.8.1]) yields √ n (( M ,n , M ,n , ˜ M ,n ) T − ( E X ∞ , E X ∞ , E ˜ X ˜ X h ) T ) d → N (0 , ˜Σ ) . To simplify notations, introduce the following two mappings: h : ( x , x , x ) ( h ( x , x ) , h ( x , x ) , h ( x , x , x )) , and η : ( x , x , x ) ( x , x x x , x ) , where η is the inverse transform of (2.11). Then from delta method, we have √ n ( h ( M ,n , M ,n , ˜ M ,n ) T − h ( κ, θ, σ )) T d → N (0 , ¯Σ)where ¯Σ = ∇ h (Θ) ˜Σ ∇ h (Θ) T . Finally, applying the delta method again, wearrive at the asymptotic behavior of the ergodic estimators: √ n ( η ( h ( M ,n , M ,n , ˜ M ,n )) T − η ( h (Θ)) T ) ⇒ N (0 , Σ)where Σ = ∇ η ( h (Θ)) ˜Σ ∇ η ( h (Θ)) T (2.17)completing the proof of the theorem. (cid:3) Remark 2.1.
We mention that we do not know the monotonicity of g ′ ( σ ) intheory. But we can observe it numerically. Since σ >
0, to investigate the signof g ′ ( σ ) is equivalent to discuss the sign of ∂∂σ g ( u, v, σ ). ∂∂σ g ( u, v, σ ) = − h Z ∞ Z ∞ xy [ m ( y ) ∞ X i =1 ˜ λ i e − ˜ λ i σ h ϕ i ( x ) ϕ i ( y )] dxdy. For fixed u and v an example of the values of h ∂∂σ g ( u, v, σ ) is plotted inFigure 1. It shows that the partial derivatives are always less than zero on theconcerned interval. 10 .2 0.4 0.6 0.8 1.0 σ - - Figure 1: The graph of h ∂∂σ g ( u, v, σ ) for u = 1, v = √ / . h = 0 .
3. Numerical experiments
In this section, we present a numerical experiment to illustrate our method.In Table 1 we set the following true parameters σ = 0 . κ = 1, θ = 1. The timestep is fixed by h = 0 .
5. In the experiments, we use the truncation p N,h ( x, y ) = π ( x ) + N X i =1 e − ˜ λ i σ h ϕ i ( x ) ϕ i ( y )in (2.6). Here we take N = 12. It can be seen that our estimators for allthe parameters, including ˆ σ n are strongly consistent. On the other hand wealso include variation estimator ˆ σ c,n = qP nk =1 ( X ( k +1) h − X kh ) / ( nh ), whichis observed not consistent. Acknowledgement
Y. Hu is supported by an NSERC discovery grant and a startup fund of Uni-versity of Alberta. Y. Xi is supported by the National Natural Science Foun-dation of China (Grant No. 11631004, 71532001) and the China ScholarshipCouncil. 11 able 1: The estimators (ˆ κ, ˆ θ, ˆ σ, ˆ σ c ) for different values of n . n( × )2 3 4 5 6 8ˆ κ θ σ σ c References [1] Billingsley, P., 1961. Statistical inference for Markov processes. Universityof Chicago Press.[2] Bo, L., Wang, Y., Yang, X., Zhang, G., 2011. Maximum likelihood estima-tion for reflected Ornstein-Uhlenbeck processes. J. Statist. Plann. Inference141, 588–596. doi: .[3] Brooks, S., Gelman, A., Jones, G.L., Meng, X.L. (Eds.), 2011. Handbook ofMarkov chain Monte Carlo. Chapman & Hall/CRC Handbooks of ModernStatistical Methods, CRC Press, Boca Raton, FL. doi: .[4] Cheng, Y., Hu, Y., Long, H., 2020. Generalized moment estimators for α -stable Ornstein-Uhlenbeck motions from discrete observations. Stat. In-ference Stoch. Process. 23, 53–81. doi: .[5] Hu, Y., Lee, C., Lee, M.H., Song, J., 2015. Parameter estimation for reflectedOrnstein-Uhlenbeck processes with discrete observations. Stat. InferenceStoch. Process. 18, 279–291. doi: .[6] Lebedev, N.N., 1965. Special functions and their applications. RevisedEnglish edition. Translated and edited by Richard A. Silverman, Prentice-Hall, Inc., Englewood Cliffs, N.J. 127] Linetsky, V., 2005. On the transition densities for reflected diffusions. Adv.in Appl. Probab. 37, 435–460. doi: .[8] Ward, A.R., Glynn, P.W., 2003. A diffusion approximation fora Markovian queue with reneging. Queueing Syst. 43, 103–128.doi:10.1023/A:1021804515162