The paper is devoted to a stochastic heat equation with a mixed fractional Brownian noise. We investigate the covariance structure, stationarity, upper bounds and asymptotic behavior of the solution. Based on its discrete-time observations, we construct a strongly consistent estimator for the Hurst index H and prove the asymptotic normality for H<3/4. Then assuming the parameter H to be known, we deal with joint estimation of the coefficients at the Wiener process and at the fractional Brownian motion. The quality of estimators is illustrated by simulation experiments.
Stochastic partial differential equationmixed fractional Brownian motionHurst index estimationstrong consistencyasymptotic normality60G2260H1562F1062F12KR was supported by the Sydney Mathematical Research Institute under Ukrainian Visitors Program.Introduction
The paper is devoted to parameter estimation in a stochastic heat equation of the following form: ∂u∂t−12·∂2u∂x2(t,x)=σB˙xH+κW˙x,t>0,x∈R,0,\hspace{0.2778em}x\in \mathbb{R},\]]]>u(0,x)=0,x∈R. The right-hand side of (1) is a mixed fractional noise. It consists of two independent stochastic processes, namely, a fractional Brownian motion BH={BxH,x∈R} with Hurst parameter H∈(0,1) and a Wiener process W={Wx,x∈R}; σ and κ are some positive coefficients. The solution to the problem (1)–(2) is understood in mild sense; the precise definition will be given in Section 2.
The mixed fractional Brownian motion was first introduced by P. Cheridito [8] in order to model financial markets that are simultaneously arbitrage-free and have a memory. The properties of this process were studied in [36]. We refer to the book [28] for a detailed presentation of the modern theory in this area. More involved mixed fractional models described by stochastic differential equations are the subject of numerous publications [13, 35, 33, 17] appeared during last decades. Such equations can be used to model processes on financial markets, where two principal random noises influence the prices. The first source of randomness is the stock exchange itself with thousands of agents. The noise coming from this source can be assumed white and is best modeled by a Wiener process. The second one has the financial and economical background and can be modeled by a fractional Brownian motion BH. Stochastic partial differential equations with such noises can be used, in particular, for the modeling of instantaneous forward rates, where the space variable corresponds to time until maturity [9, 22]. Such equations arise also in geophysics, especially in physical oceanography [30] and in geostatistics [34]. For example, in models for sea surface temperature, noise terms can represent various heat fluxes and ocean processes [29].
Existence, uniqueness and properties of solutions for stochastic differential equations with mixed noise were studied in various papers [12, 16, 17, 25–27]. A stochastic heat equation with white and fractional noises was investigated in [24]. Several approaches to parameter identification in simple linear mixed fractional models for various observation schemes were proposed in [6, 14, 23, 20, 21]. The problem of drift parameter estimation in a mixed stochastic differential equation of a general form was studied in [19]. The statistical problems for the mixed fractional Vasicek model were investigated in the recent papers [7, 31].
Similarly to our previous papers the solution u(t,x) is observed at equidistant spatial points for a several fixed time instants. On one hand, there are many practical cases where the solution is observed at some discrete space points such as temperature of a heated body or velocity of a turbulent flow. In many cases the measurements with a high space resolution are available, but the time series are short. For example, this is the case for satellite observations of sea surface temperature, see [30]. For this reason, it is suitable to assume that u(t,x) is observed at a large number of space points xk and only few different time instants ti. On the other hand, observing the solution at three time points is enough to construct estimators for unknown parameters σ, κ and H. Nevertheless, the additional information of observing the solution discretely in time can be consolidated by taking the (weighted) average of the estimators similarly to [5] or [9].
The present paper is devoted to the problem of estimating unknown parameters H, σ, κ in the equation (1), by discrete observations of its solution u(t,x). The results of this paper are an extension of our previous works [2, 3], where the problems of estimating H and σ were studied for the equation (1) with the fractional Brownian motion only (that is, for κ=0). The diffusion parameter estimator for SPDE with white noise and its properties was considered in [4]. Similarly to the mentioned articles, in the present paper we start with proving stationarity and ergodicity of the solution u(t,x) as a function of the spatial variable x by analyzing the behavior of the covariance function. Based on these results we construct a strongly consistent estimator of H (assuming that the parameters σ and κ are unknown). The asymptotic normality of this estimator is proved for any H∈(0,12)∪(12,34). Then we consider the problem of estimating the couple of parameters (σ,κ) when the value of H is known. We prove strong consistency of the estimator and investigate its asymptotic normality.
The paper is organized as follows. In Section 2 we introduce the definition of a mild solution to SPDE (1) and present its properties. Furthermore, we prove the limit theorem for it, needed for establishing properties of statistical estimators. The statistical problems are investigated in Section 3. In Subsection 3.1 we construct an estimator of the Hurst index H and prove its strong consistency and asymptotic normality. Subsection 3.2 is devoted to the estimators of the parameters σ and κ and their asymptotic properties. Numerical results are presented and discussed in Section 4.
Preliminaries
Assume that BH=BxH,x∈R is a two-sided fractional Brownian motion with the Hurst index H∈(0,1), while W=Wx,x∈R is a Wiener process, independent of BH. Let G be Green’s function of the heat equation, that is
G(t,x)=12πtexp−x22t,ift>0,δ0(x),ift=0.0\text{,}\\ {} {\delta _{0}}(x),\hspace{1em}& \hspace{2.5pt}\text{if}\hspace{2.5pt}t=0.\end{array}\right.\]]]>
Similarly to [2–4] (see also [10] and the references cited therein), we define a solution to SPDE (1) in a mild sense as follows.
The random field u(t,x),t≥0,x∈R defined by
u(t,x)=σ∫0t∫RG(t−s,x−y)dByHds+κ∫0t∫RG(t−s,x−y)dWyds
is called a solution to SPDE (1)–(2).
As shown in [2], both stochastic integrals in (3) exist as pathwise Riemann–Stieltjes integrals. This fact follows from the Hölder regularity of the integrands and integrators. Namely, Green’s function is obviously Lipshitz continuous, while sample paths of the fractional Brownian motion are Hölder continuous up to order H. Such regularity guarantees the existence of the first integral in (3). The second integral is also well defined, since the integrand is square integrable, see [4, Theorem 2.1].
The next proposition summarizes basic properties of the solution u(t,x). These properties, especially stationarity and ergodicity, enable us to construct and investigate statistical estimators for H, κ and σ.
Letu=u(t,x),t∈[0,T],x∈Rbe a solution to (1) defined by (3). Then the following properties hold.
For allt,s∈[0,T]andx,z∈R,cov(u(t,z),u(s,x+z))=cov(u(t,0),u(s,x))=12π∫0t∫0s(q+r)−32∫Rσ2Hy2H−1+κ22×(signy)(y−x)exp−(y−x)22(q+r)dydqdr.Consequently, for fixed distinct pointst1,…,tn∈[0,T], the stochastic processu(t1,x)⋮u(tn,x),x∈R, is a multivariate stationary Gaussian process.
The variance ofu(t,x)is given byE[u(t,x)2]=σ2vt(H)+κ2vt12,t>0,x∈R,0,\hspace{0.2778em}x\in \mathbb{R},\]]]>wherevt(H)=cHtH+1,cH=2H+1(2H−1)Γ(H+12)π(H+1),
Γ denotes the gamma function.
For allt,s∈[0,T]andx>00$]]>, the covariance function admits the following upper bound:cov(u(t,0),u(s,x))≤CHtsσ2x2H−2+κ2x−1,whereCHis a positive constant depending only on H.
For allt,s∈[0,T]andx∈R,cov(u(t,x),u(s,x))=σ22HΓ(H+12)(t+s)H+1−tH+1−sH+1π(H+1)+κ2232(t+s)32−t32−s323π.
For a fixedt>00$]]>, the random processu(t,x),x∈Ris ergodic.
The proposition follows from the corresponding results for the equation with pure fractional noise studied in [2, 3]. Indeed, all the statements are based on the properties of the covariance function of the solution. However, since BH and W are independent, we see that this covariance function can be represented as
cov(u(t,x),u(s,z))=σ2cov(ub(t,x),ub(s,z))+κ2cov(uw(t,x),uw(s,z)),
where
ub(t,x)=∫0t∫RG(t−s,x−y)dByHds,uw(t,x)=∫0t∫RG(t−s,x−y)dWyds.
Then combining the equality (9) with statements of [3, Prop. 2.2], we immediately obtain formulas (4), (7) and (8). The equality (5) follows from (9) and [2, Prop. 3]. Finally, the last statement of the proposition holds, because the solution u(t,x),x∈R is a stationary Gaussian process, whose covariance function vanishes as x→∞, according to (8). Hence, the process u(t,·) is ergodic. □
Let us fix some δ>00$]]> and consider the following sequence:
VN(t)=1N∑k=1Nu(t,kδ)2,t>0,N∈N.0,\hspace{0.2778em}N\in \mathbb{N}.\]]]>
The sequence (10) will serve as a statistic for parameter estimation problems in Section 3. We introduce the following notation in addition to (6):
μ(t):=EVN(t)=E[u(t,0)2]=σ2vt(H)+κ2vt12,ρtsH(k)=cov(u(t,kδ),u(s,0)),rts(H)=2∑k=−∞∞ρtsH(k)2,t,s>0.0.\end{array}\]]]>
The next theorem describes an asymptotic behavior of the stochastic process VN. It gives a law of large numbers and a central limit theorem for its finite-dimensional distributions (VN(t1),…,VN(tn)) as N→∞. This result is crucial for construction of the estimators and for establishing their asymptotic properties.
LetH∈(0,1).
For anyt>00$]]>,VN(t)→μ(t)a. s., asN→∞.
If additionallyH∈(0,34), then for anyn≥1and any distinct positivet1,…,tn,NVN(t1)−μ(t1)⋮VN(tn)−μ(tn)→dN(0,R)asN→∞,whereN(0,R)is the multivariate normal distribution with zero mean and the covariance matrixR=rtitj(H)i,j=1n.
1. The ergodic theorem implies that for any t>00$]]>1N∑k=1Nu(t,kδ)2→E[u(t,0)2]a. s. asN→∞,
which is equivalent to (12).
2. In order to prove the convergence (13), we shall apply the Cramér–Wold theorem. In other words, we need to show that for all α1,…,αn∈R, the convergence
N∑i=1nαiVN(ti)−μ(ti)→dN0,s2
holds with the asymptotic variance
s2=∑i=1nαi2rtiti(H)+2∑1≤i<j≤nαiαjrtitj(H).
Representing VN as the sum (10) and using (11), we rewrite (14) in the form
1N∑k=1N∑i=1nαiu(ti,kδ)2−Eu(ti,kδ)2→dN(0,s2).
Further, since u(t1,kδ)⋮u(tn,kδ),x∈R, is a multivariate stationary Gaussian sequence, the convergence (16) can be established by application of the multivariate Breuer–Major theorem [1, Theorem 4]. In order to verify the conditions of this theorem, note that the function F(x1,…,xn)=∑i=1nαixi2 has Hermite rank 2, therefore we need to check the convergence of the series:
∑k=−∞∞ρtitjH(k)2<∞,i,j=1,…,n.
It follows immediately from the upper bound (7) that these power series converge if and only if H∈(0,34). Thus, the assumptions of [1, Theorem 4] are satisfied, whence the convergence (16) holds with the following asymptotic variance:
s2=Var(F(u(t1,0),…,u(tn,0)))+2∑k=1∞cov(F(u(t1,0),…,u(tn,0)),F(u(t1,kδ),…,u(tn,kδ))).
Now we must only check that the asymptotic variance s2 satisfies (15). By the definition of the function F, we have
s2=Var∑i=1nαiu(ti,0)2+2∑k=1∞cov∑i=1nαiu(ti,0)2,∑j=1nαju(tj,kδ)2=∑i=1nαi2Var(u(ti,0)2)+2∑1≤i<j≤nαiαjcovu(ti,0)2,u(tj,0)2+2∑k=1∞∑i=1nαi2covu(ti,0)2,u(ti,kδ)2+2∑k=1∞∑1≤i<j≤nαiαj(covu(ti,0)2,u(tj,kδ)2+covu(tj,0)2,u(ti,kδ)2).
Now we can use the following well-known fact: if ξ and η are zero-mean Gaussian random variables, then cov(ξ2,η2)=2cov(ξ,η)2, in particular, Var(ξ2)=2Var(ξ)2 (this is a corollary of the Isserlis theorem [18]). Then we get
s2=2∑i=1nαi2ρtitiH(0)2+4∑1≤i<j≤nαiαjρtitjH(0)2+4∑k=1∞∑i=1nαi2ρtitiH(k)2+4∑k=1∞∑1≤i<j≤nαiαjρtitjH(k)2+ρtjtiH(k)2.
Taking into account the equality ρtsH(k)=ρstH(−k), we may rewrite this expression in the more compact form:
s2=2∑i=1nαi2∑k=−∞+∞ρtitiH(k)2+4∑1≤i<j≤nαiαj∑k=−∞+∞ρtitjH(k)2.
Thus the equality (15) is verified. This completes the proof of Theorem 1. □
Parameter estimation
Let us consider the following statistical problem. It is supposed that for fixed t1,…,tn and fixed step δ>00$]]>, the random field u given by (3) is observed at the points xk=kδ, k=1,…,N. So the observations have the following form:
u(ti,kδ),i=1,…,n,k=1,…,N.
Our aim is to estimate the unknown parameters H, σ and κ. We shall do this it two steps. We start with construction of a strongly consistent estimator of H that does not depend on κ and σ. Also, we shall establish its asymptotic normality. Then assuming that H is known, we shall estimate the parameters σ and κ.
In what follows we assume that H≠12, because otherwise the model is non-identifiable. The parameters σ and κ are assumed to be positive.
Estimation of H
In order to estimate H without knowledge of σ and κ, it suffices to observe u(ti,xk) only at three time instants t1<t2<t3. If we write (12) at these points and replace convergences with equalities, we shall get the following system of equations
VN(t1)=σ2cHt1H+1+κ2c12t13/2,VN(t2)=σ2cHt2H+1+κ2c12t23/2,VN(t3)=σ2cHt3H+1+κ2c12t33/2.
Excluding the unknown parameter κ from the system, we obtain
t2−3/2VN(t2)−t1−3/2VN(t1)=σ2cHt2H−12−t1H−12,t3−3/2VN(t3)−t2−3/2VN(t2)=σ2cHt3H−12−t2H−12.
Then excluding σ we arrive at the following estimating equation for H:
t3H−12−t2H−12t2H−12−t1H−12=t3−3/2VN(t3)−t2−3/2VN(t2)t2−3/2VN(t2)−t1−3/2VN(t1).
The solution of (19) (if exists), can be viewed as an estimator of H.
Note that the left-hand side of (19) is indeterminate for H=1/2. However, it is easy to see by l’Hôpital’s rule that there exists the limit
limH→12t3H−12−t2H−12t2H−12−t1H−12=limH→12t3H−12logt3−t2H−12logt2t2H−12logt2−t1H−12logt1=logt3−logt2logt2−logt1.
Therefore, one may define by continuity
f(H):=t3H−1/2−t2H−1/2t2H−1/2−t1H−1/2,ifH≠12,logt3−logt2logt2−logt1ifH=12.
Then the estimator of H is defined as
HˆN=f(−1)t3−3/2VN(t3)−t2−3/2VN(t2)t2−3/2VN(t2)−t1−3/2VN(t1),
where f(−1) denotes the inverse function of f. In order to prove its existence, we need to establish that f:R→(0,∞) is a one-to-one function. This is true, since f is always a strictly increasing function (see Fig. 1) as shown in the following lemma.
The function f(H) for t1=1, t2=2, t3=3
For any0<t1<t2<t3, the functionf:R→(0,∞)defined by (20) is strictly increasing with respect to H.
We prove the statement for the case H∈(12,∞). The interval (−∞,12) is considered similarly. The derivative of f with respect to H is equal to
f′(H)=t2H−12−t1H−12t3H−12logt3−t2H−12logt2t2H−12−t1H−122−t3H−12−t2H−12t2H−12logt2−t1H−12logt1t2H−12−t1H−122.
Therefore, it suffices to prove the inequality
t2H−12−t1H−12t3H−12logt3−t2H−12logt2>t3H−12−t2H−12t2H−12logt2−t1H−12logt1.\left({t_{3}^{H-\frac{1}{2}}}-{t_{2}^{H-\frac{1}{2}}}\right)\left({t_{2}^{H-\frac{1}{2}}}\log {t_{2}}-{t_{1}^{H-\frac{1}{2}}}\log {t_{1}}\right).\end{aligned}\]]]>
In order to establish (23), observe that the function h(x)=xlogx, x>00$]]>, is strictly convex (indeed, its second derivative h″(x)=1/x>00$]]>). This means that for any α∈(0,1), x>00$]]> and y>00$]]>,
αh(x)+(1−α)h(y)>h(αx+(1−α)y).h\big(\alpha x+(1-\alpha )y\big).\]]]>
Let us take
x=t3H−12>0,y=t1H−12>0,andα=t2H−12−t1H−12t3H−12−t1H−12∈(0,1).0,\hspace{1em}y={t_{1}^{H-\frac{1}{2}}}>0,\hspace{1em}\text{and}\hspace{1em}\alpha =\frac{{t_{2}^{H-\frac{1}{2}}}-{t_{1}^{H-\frac{1}{2}}}}{{t_{3}^{H-\frac{1}{2}}}-{t_{1}^{H-\frac{1}{2}}}}\in (0,1).\]]]>
Then
1−α=t3H−12−t2H−12t3H−12−t1H−12,αx+(1−α)y=t2H−12,
and (24) becomes
t2H−12−t1H−12t3H−12−t1H−12(H−12)t3H−12logt3+t3H−12−t2H−12t3H−12−t1H−12(H−12)t1H−12logt1>(H−12)t2H−12logt2,(H-\frac{1}{2}){t_{2}^{H-\frac{1}{2}}}\log {t_{2}},\end{aligned}\]]]>
which is equivalent to (23). □
The above lemma yields that the estimator HˆN is well defined at least for sufficiently large N (when the right-hand side of estimating equation (19) becomes positive). The asymptotic properties of HˆN are summarized in the following theorem, which is the first main result of the paper.
1. For anyH∈(0,12)∪(12,1),HˆNis a strongly consistent estimator of the parameter H asN→∞.
2. ForH∈(0,12)∪(12,34), the estimatorHˆNis asymptotically normal:NHˆN−H→dN(0,ς2)asN→∞,whereς2=1D2σ4cH2∑i,j=13rtitj(H)LiLj,L1=t3H−12−t2H−12t13/2,L2=t1H−12−t3H−12t23/2,L3=t2H−12−t1H−12t33/2,D=t2H−12−t1H−12t3H−12logt3−t2H−12logt2−t3H−12−t2H−12t2H−12logt2−t1H−12logt1.
The inequality (23) from the proof of Lemma 1 implies that D>00$]]> for all H≠1/2.
1. The strong consistency follows from the construction of the estimator. Indeed, (12) implies that
t3−3/2VN(t3)−t2−3/2VN(t2)t2−3/2VN(t2)−t1−3/2VN(t1)→f(H)a. s., asN→∞.
Then the convergence HˆN→H a. s. as N→∞ follows from (21) and (25) due to the continuity of f(−1).
2. By taking expectations in the equalities (18), we get the following relations
t2−3/2μ(t2)−t1−3/2μ(t1)=σ2cHt2H−12−t1H−12,t3−3/2μ(t3)−t2−3/2μ(t2)=σ2cHt3H−12−t2H−12,
whence
t3−3/2μ(t3)−t2−3/2μ(t2)t2−3/2μ(t2)−t1−3/2μ(t1)=t3H−12−t2H−12t2H−12−t1H−12=f(H),
or
H=f(−1)t3−3/2μ(t3)−t2−3/2μ(t2)t2−3/2μ(t2)−t1−3/2μ(t1).
Therefore,
NHˆN−H=N(g(VN(t1),VN(t2),VN(t3))−g(μ(t1),μ(t2),μ(t3))),
where
g(x1,x2,x3)=f(−1)t3−3/2x3−t2−3/2x2t2−3/2x2−t1−3/2x1.
In order to derive the desired asymptotic normality from the convergence (13), we shall apply the delta method. Denoting
h(x)=ddxf(−1)(x)=1f′(f(−1)(x)),
we see that the partial derivatives of g equal
g1′(x1,x2,x3)=ht3−3/2x3−t2−3/2x2t2−3/2x2−t1−3/2x1·t1−3/2t3−3/2x3−t2−3/2x2t2−3/2x2−t1−3/2x12,g2′(x1,x2,x3)=−ht3−3/2x3−t2−3/2x2t2−3/2x2−t1−3/2x1·t2−3/2t3−3/2x3−t1−3/2x1t2−3/2x2−t1−3/2x12,g3′(x1,x2,x3)=ht3−3/2x3−t2−3/2x2t2−3/2x2−t1−3/2x1·t3−3/2t2−3/2x2−t1−3/2x1.
By the delta method, we derive from (13) the convergence (2) with
ς2=∑i,j=13rtitj(H)gi′gj′(μ(t1),μ(t2),μ(t3)).
It remains to prove that gi′(μ(t1),μ(t2),μ(t3))=Li/(Dσ2cH), i=1,2,3. Let us consider in detail the proof of this equality for i=1, the cases i=2 and i=3 are considered similarly. Using successively (27), (28) and (22), we get
ht3−3/2x3−t2−3/2x2t2−3/2x2−t1−3/2x1|xi=μ(ti)=h(f(H))=1f′(H)=t2H−12−t1H−122D.
After inserting this expression into (29) and taking into account the relations (26), we obtain
g1′(μ(t1),μ(t2),μ(t3))=t2H−12−t1H−122D·t1−3/2t3H−12−t2H−12σ2cHt2H−12−t1H−122=t1−3/2t3H−12−t2H−12Dσ2cH=L1Dσ2cH.
Note also that the above representation yields g1′(x1,x2,x3)≠0 in the neighborhood of the point (μ(t1),μ(t2),μ(t3)), which is necessary for applying the delta method. The derivatives g2′ and g3′ are considered similarly. □
The estimator of H was obtained as a solution to some exponential equation. However, it would be more convenient for applications and modeling to have the explicit form of the estimator. It turns out that in some particular cases it is possible to solve the estimating equation explicitly. Let us consider such example in more detail.
Assume that t1=h>00$]]>, t2=2h, t3=4h. Substituting these values in the definition of f, we get
f(H)=t3H−12−t2H−12t2H−12−t1H−12=4H−12hH−12−2H−12hH−122H−12hH−12−hH−12=2H−12.
Therefore, f(−1)(x)=12+log2x, x>00$]]>, consequently (21) becomes1
We use log2+ rather than log2 so that the right-hand side of (30) is always well defined. The function log2+x is defined as log2x if x>00$]]> and 0 if x≤0. Note that the expression under log2+ in (30) eventually becomes positive, since in converges to 2H−1/2.
HˆN=12+log2+t3−3/2VN(t3)−t2−3/2VN(t2)t2−3/2VN(t2)−t1−3/2VN(t1).
In this case
L1=4H−12−2H−12h2−H,L2=−4H−12−1232h2−H,L3=2H−12−18h2−H,
and
D=h2H−12H−12−14H−12log(4h)−2H−12log(2h)−h2H−14H−12−2H−122H−12log(2h)−logh=h2H−14H−12−2H−122H−12log(4h)−log(2h)−2H−12log(2h)+logh=h2H−12H−122H−12−12H−12log2−log2=h2H−12H−122H−12−12log2.
Denoting ℓi=D−1Lih1+Hlog2, we arrive at the following result.
Lett1=h>00$]]>,t2=2h,t3=4h. Then the estimatorHˆNcan be written in the explicit form (30). In this case Theorem2holds withς2=1σ4cH2h2+2H(log2)2∑i,j=13rtitj(H)ℓiℓj,whereℓ1=12H−12−1,ℓ2=−2H−12+12H+12H−12−1,ℓ3=12H+522H−12−1.
Evidently, the explicit form of the estimator can be obtained also in a slightly more general case, when t1=h, t2=ah, t3=a2h with some a>00$]]>. This leads to the estimator of the form (30) with loga+ instead of log2+.
Estimation of σ and κ
Now we assume that the Hurst index H is known and investigate the estimation of the coefficients σ and κ. From the first two equations of (17), one can derive the following estimators:
σˆN2=t1−3/2VN(t1)−t2−3/2VN(t2)cHt1H−1/2−t2H−1/2,κˆN2=t1−1−HVN(t1)−t2−1−HVN(t2)c12t11/2−H−t21/2−H.
Now we are ready to formulate and prove the second main result of the paper.
1. For anyH∈(0,12)∪(12,1),(σˆN2,κˆN2)is a strongly consistent estimator of the parameter(σ2,κ2)asN→∞.
2. ForH∈(0,12)∪(12,34), the estimator(σˆN2,κˆN2)is asymptotically normal:NσˆN2−σ2κˆN2−κ2→dN(0,Σ)asN→∞,where the asymptotic covariance matrix Σ consists of the following elements:Σ11=t1−3(rt1t1(H)+rt1t2(H))+t2−3(rt1t2(H)+rt2t2(H))cH2t12H−1−2(t1t2)H−12+t22H−1,Σ12=Σ21=t1−52−H(rt1t1(H)+rt1t2(H))+t2−52−H(rt1t2(H)+rt2t2(H))cHc122−t1H−12t212−H−t112−Ht2H−12,Σ22=t1−2−H(rt1t1(H)+rt1t2(H))+t2−2−H(rt1t2(H)+rt2t2(H))c122t11−2H−2(t1t2)12−H+t21−2H.
Using the definition (31) of the estimator σˆN2, we rewrite the error σˆN2−σ in the following form
σˆN2−σ2=t1−3/2VˆN(t1)−t2−3/2VˆN(t2)−σ2cHt1H−12+σ2cHt2H−12cHt1H−12−t2H−12=1cHt1H−12−t2H−12(t1−3/2VˆN(t1)−σ2cHt1H+1−κ2c12t13/2−t2−3/2VˆN(t2)−σ2cHt2H+1−κ2c12t23/2)=t1−3/2VˆN(t1)−μ(t1)−t2−3/2VˆN(t2)−μ(t2)cHt1H−12−t2H−12,
where the last equality follows from (11). Similarly, one can represent κˆN2−κ2 as follows:
κˆN2−κ2=t1−1−HVˆN(t1)−μ(t1)−t2−1−HVˆN(t2)−μ(t2)c12t112−H−t212−H.
Hence, we see that the random vector in the left-hand side of (32) is a linear transformation of the left-hand side of (13) (for n=2), namely
NσˆN2−σ2κˆN2−κ2=t1−3/2cHt1H−12−t2H−12−t2−3/2cHt1H−12−t2H−12t1−1−Hc12t112−H−t212−H−t2−1−Hc12t112−H−t212−HNVˆN(t1)−μ(t1)NVˆN(t2)−μ(t2).
Therefore, taking into account the convergence (13), we conclude that (33) weakly converges in distribution to a bivariate normal distribution with the following covariance matrix:
Σ=t1−3/2cHt1H−12−t2H−12−t2−3/2cHt1H−12−t2H−12t1−1−Hc12t112−H−t212−H−t2−1−Hc12t112−H−t212−Hrt1t1(H)rt1t2(H)rt1t2(H)rt2t2(H)×t1−3/2cHt1H−12−t2H−12t1−1−Hc12t112−H−t212−H−t2−3/2cHt1H−12−t2H−12−t2−1−Hc12t112−H−t212−H=Σ11Σ12Σ12Σ22.
This completes the proof. □
Maximum likelihood estimation and Fisher information
In this subsection we analyze the efficiency of the estimator (σˆN2,κˆN2) by comparing it to the maximum likelihood estimator (MLE). Note that MLE is hard to compute, however, it is possible to identify the corresponding Fischer information matrix. We use for construction of MLE the same observations as in the previous subsection, namely let the observation vector be
XN=(u(t1,δ),u(t1,2δ),…,u(t1,Nδ),u(t2,δ),u(t2,2δ),…,u(t2,Nδ))⊤.
Obviously, XN has 2N-dimensional Gaussian distribution with probability density
f(XN,θ)=(2π)−NdetΓN−12exp−12XN⊤Γ−1XN,
where ΓN is the covariance matrix of XN, that is,
ΓN=ρt1t1H(0)⋯ρt1t1H(N−1)ρt2t1H(0)⋯ρt2t1H(N−1)⋮⋱⋮⋮⋱⋮ρt1t1H(N−1)⋯ρt1t1H(0)ρt2t1H(N−1)⋯ρt2t1H(0)ρt2t1H(0)⋯ρt2t1H(N−1)ρt2t2H(0)⋯ρt2t2H(N−1)⋮⋱⋮⋮⋱⋮ρt2t1H(N−1)⋯ρt2t1H(0)ρt2t2H(N−1)⋯ρt2t2H(0).
Due to (9), this matrix can be decomposed as ΓN=σ2ΓNb+κ2ΓNw, where ΓNb and ΓNw are the covariance matrices of
(ub(t1,δ),ub(t1,2δ),…,ub(t1,Nδ),ub(t2,δ),ub(t2,2δ),…,ub(t2,Nδ))⊤
and
(uw(t1,δ),uw(t1,2δ),…,uw(t1,Nδ),uw(t2,δ),uw(t2,2δ),…,uw(t2,Nδ))⊤,
respectively. The log-likelihood function is
ℓ(XN,θ)=−Nlog(2π)−12log(detΓN)−12XN⊤ΓN−1XN.
Then, MLE of θ=(σ2,κ2) is obtained as the solution to the following system of equations: ∂ℓ(XN,θ)∂σ2=−12tr(ΓN−1ΓNb)+12XN⊤ΓN−1ΓNbΓN−1XN=0,∂ℓ(XN,θ)∂κ2=−12tr(ΓN−1ΓNw)+12XN⊤ΓN−1ΓNwΓN−1XN=0 (here and after we use the differentiation formulas of a matrix with respect to given parameter,2
For square matrices X and Y, ∂(XY)=(∂X)Y+X(∂Y), ∂(X−1)=−X−1(∂X)X−1, ∂(ln(detX))=tr(X−1∂X), ∂(tr(X))=tr(∂X).
see, e. g., [32] for more details).
The maximum likelihood estimator θˆNmle of θ can hardly be written in the explicit form, since the estimating equations involve the inverse matrix ΓN−1, which depends nonlinearly on σ2 and κ2. However, using the general theory of maximum likelihood estimation for dependent observations [11], it is possible to establish the asymptotic normality in the form
TN(θ)12θˆNmle−θ→dN(0,I2)asn→∞,
where TN(θ) is the Fisher information matrix and I2 is the 2×2 identity matrix. The rigorous proof of (36) as well as a careful analysis of the asymptotic behavior of TN(θ) requires an additional investigation. To the best of our knowledge, even for much simpler model of the mixed fractional Brownian motion, this problem has not been completely solved yet, see the recent paper [15] for details. Therefore, here we restrict ourselves to the identification of the matrix TN(θ).
The Fisher information matrixTN(θ)has the formTN(θ)=12trΓN−1ΓNb212trΓN−1ΓNbΓN−1ΓNw12trΓN−1ΓNbΓN−1ΓNw12trΓN−1ΓNw2.
In order to identify TN(θ), let us calculate the second derivatives. Note that
∂∂σ2ΓN−1ΓNb=∂∂σ2ΓN−1ΓNb=−ΓN−1ΓNbΓN−1ΓNb=−ΓN−1ΓNb2,∂∂σ2ΓN−1ΓNbΓN−1=∂∂σ2ΓN−1ΓNbΓN−1+ΓN−1ΓNb∂∂σ2ΓN−1=−ΓN−1ΓNb2ΓN−1−ΓN−1ΓNbΓN−1ΓNbΓN−1=−2ΓN−1ΓNb2ΓN−1.
Hence,
∂2ℓ(XN,θ)∂(σ2)2=12trΓN−1ΓNb2−XN⊤ΓN−1ΓNb2ΓN−1XN.
Taking expectations, we obtain that the corresponding element of the Fisher information matrix equals
−E∂2ℓ(XN,θ)∂(σ2)2=−12trΓN−1ΓNb2+EXN⊤ΓN−1ΓNb2ΓN−1XN=12trΓN−1ΓNb2,
since for any matrix A=(aij)i,j=1,…2N, we have the equality EXN⊤AXN=∑i,jaijEXiXj=trAΓN. Arguing as above, one can write the derivatives
∂2ℓ(XN,θ)∂(κ2)2=12trΓN−1ΓNw2−XN⊤ΓN−1ΓNw2ΓN−1XN∂2ℓ(XN,θ)∂σ2∂κ2=12trΓN−1ΓNwΓN−1ΓNb−12XN⊤ΓN−1ΓNbΓN−1ΓNw+ΓNwΓN−1ΓNbΓN−1XN,
and calculate their expectations, identifying other elements of TN(θ). □
1. Similarly to the previous subsection, in the case H=12 it is impossible to estimate both parameters, σ2 and κ2, simultaneously. Only estimation of the sum σ2+κ2 is possible. In this case ΓNb=ΓNw, therefore the estimation equations (34) and (35) coincide.
2. The results of this subsection are valid for any other observations vector of the form X=(u(ti,xk),i=1…,M,k=1,…N) and its covariance matrix Γ (with decomposition Γ=σ2Γb+κ2Γw).
3. Similar approach can be applied to the case, when H is unknown, that is, to the problem of estimation of all three parameters σ2, κ2 and H.
Simulations
Let us illustrate the theoretical properties of the estimator by some numerical results. We consider the model with the coefficients σ=κ=1 for various values of H. For each value of the Hurst index H, we simulate 50 sample paths of the solution u(t,x) to the equation (1). The trajectories of a solution are generated by the discretization of the formula (3).
We choose t1=0.25, t2=0.5, t3=1 as the moments of observations, so that the conditions of Corollary 1 are satisfied and the estimator of H can be computed by the explicit formula (30). For each ti we observe u(ti,kδ), k=1,…,N, with the step δ=1.
Means and standard deviations of the estimator HˆN
N
28
29
210
211
212
H=0.1
Mean
−0.0121
−0.0299
0.0353
0.0527
0.0682
Std. dev.
0.5237
0.5552
0.2061
0.1242
0.0752
H=0.2
Mean
0.1616
0.1415
0.1910
0.1961
0.1799
Std. dev.
0.3419
0.2243
0.1540
0.0897
0.0692
H=0.3
Mean
0.1543
0.2845
0.2685
0.2955
0.2999
Std. dev.
0.4451
0.4177
0.1930
0.1314
0.0781
H=0.4
Mean
0.2254
0.2725
0.3129
0.2854
0.3313
Std. dev.
0.8089
0.8803
0.4661
0.1608
0.1299
H=0.6
Mean
0.4563
0.3495
0.5266
0.5618
0.5775
Std. dev.
0.7108
0.9010
0.2926
0.1992
0.1108
H=0.7
Mean
0.6384
0.7024
0.7160
0.7151
0.6980
Std. dev.
0.3391
0.1512
0.1065
0.0848
0.0511
H=0.8
Mean
0.8160
0.8042
0.8073
0.8022
0.8074
Std. dev.
0.1929
0.0922
0.0644
0.0467
0.0333
H=0.9
Mean
0.8583
0.8722
0.8815
0.8939
0.8958
Std. dev.
0.1216
0.0772
0.0671
0.0474
0.0334
Means and standard deviations of the estimator σˆN2 for σ=1, κ=1
N
28
29
210
211
212
H=0.1
Mean
1.0084
1.0511
1.0571
1.0405
1.0440
Std. dev.
0.3971
0.2933
0.1829
0.1355
0.0909
H=0.2
Mean
1.0963
1.0478
1.0407
1.0183
1.0170
Std. dev.
0.3742
0.2370
0.1888
0.1261
0.1014
H=0.3
Mean
1.0107
1.0437
0.9616
0.9939
1.0032
Std. dev.
0.4530
0.3618
0.2454
0.1673
0.1216
H=0.4
Mean
0.8117
1.0042
1.0566
1.0964
1.0660
Std. dev.
0.9423
0.7003
0.4789
0.3048
0.2167
H=0.6
Mean
1.0905
1.1579
1.0793
1.0667
1.0673
Std. dev.
0.6755
0.5750
0.3875
0.2967
0.2154
H=0.7
Mean
0.9536
1.0545
1.0653
1.0106
0.9889
Std. dev.
0.4650
0.3498
0.2526
0.1649
0.1202
H=0.8
Mean
1.0171
1.0216
1.0287
1.0252
1.0095
Std. dev.
0.4619
0.2805
0.2465
0.1720
0.1368
H=0.9
Mean
1.1814
1.1236
1.0738
1.0391
1.0308
Std. dev.
0.8814
0.7882
0.6851
0.4922
0.3681
Means and standard deviations of the estimator κˆN2 for σ=1, κ=1
N
28
29
210
211
212
H=0.1
Mean
1.0472
1.0133
1.0005
0.9953
0.9886
Std. dev.
0.1664
0.1483
0.0963
0.0606
0.0442
H=0.2
Mean
0.9608
0.9720
0.9732
0.9876
0.9859
Std. dev.
0.2430
0.1712
0.1338
0.0834
0.0720
H=0.3
Mean
0.9899
0.9638
1.0355
1.0075
1.0032
Std. dev.
0.4197
0.3056
0.2059
0.1472
0.1017
H=0.4
Mean
1.2052
1.0076
0.9417
0.9092
0.9365
Std. dev.
0.8747
0.6671
0.4480
0.2988
0.2188
H=0.6
Mean
0.8934
0.8285
0.9070
0.9258
0.9279
Std. dev.
0.6869
0.5600
0.3733
0.2869
0.2137
H=0.7
Mean
1.0483
0.9526
0.9519
0.9863
1.0040
Std. dev.
0.4351
0.3250
0.2228
0.1448
0.1078
H=0.8
Mean
0.9996
0.9938
0.9905
0.9848
0.9974
Std. dev.
0.3113
0.1941
0.1429
0.1031
0.0795
H=0.9
Mean
1.0017
0.9877
0.9767
0.9870
0.9945
Std. dev.
0.2985
0.2311
0.2072
0.1480
0.1084
Table 1 reports the means and standard deviations of HˆN for various H and N. We see that the estimates converge to the true value of the parameter H. However, the convergence is much slower compared to the estimation of H in the pure fractional case (i. e. κ=0), which was considered in [3]. The results become poorer when H approaches 1/2, or when H is close to zero. It’s worth mentioning that the best performance of HˆN is observed for large values of H (0.8 and 0.9), for which the asymptotic normality does not hold.
The means and standard deviations of the estimates σˆN2 and κˆN2 are reported in Tables 2–3. Here we also clearly see that the estimators converge to the true values of the parameters, however the results for both estimators become worse when H is close to 1/2. We observe that unlike HˆN, the estimator σˆ2 converges slowly for H=0.9, demonstrating better results for small H. Note that the similar situation for the coefficient at the fractional Brownian motion takes place in the pure fractional case, see [2].
ReferencesArcones, M.A.: Limit theorems for nonlinear functionals of a stationary Gaussian sequence of vectors. 22(4), 2242–2274 (1994). https://doi.org/10.1214/aop/1176988503Avetisian, D., Ralchenko, K.: Ergodic properties of the solution to a fractional stochastic heat equation, with an application to diffusion parameter estimation. 7(3), 339–356 (2020). https://doi.org/10.15559/20-VMSTA162Avetisian, D.A., Ralchenko, K.V.: Estimation of the Hurst and diffusion parameters in fractional stochastic heat equation. 104, 61–76 (2021). https://doi.org/10.1090/tpms/1145Avetisian, D.A., Shevchenko, G.M.: Estimation of diffusion parameter for stochastic heat equation with white noise. 3, 9–16 (2018). https://doi.org/10.17721/1812-5409.2018/3.1Bibinger, M., Trabs, M.: Volatility estimation for stochastic PDEs using high-frequency observations. 130(5), 3005–3052 (2020). MR4080737. https://doi.org/10.1016/j.spa.2019.09.002Cai, C., Chigansky, P., Kleptsyna, M.: Mixed Gaussian processes: A filtering approach. 44(4), 3032–3075 (2016). MR3531685. https://doi.org/10.1214/15-AOP1041Cai, C.H., Huang, Y.Z., Sun, L., Xiao, W.L.: Maximum likelihood estimation for mixed fractional Vasicek processes. 6(1), 44 (2022). https://doi.org/10.3390/fractalfract6010044Cheridito, P.: Mixed fractional Brownian motion. 7(6), 913–934 (2001). MR1873835. https://doi.org/10.2307/3318626Cialenco, I., Huang, Y.: A note on parameter estimation for discretely sampled SPDEs. 20(3), 2050016–28 (2020). MR4101083. https://doi.org/10.1142/S0219493720500161Cialenco, I., Kim, H.-J.: Parameter estimation for discretely sampled stochastic heat equation driven by space-only noise. 143, 1–30 (2022). MR4332773. https://doi.org/10.1016/j.spa.2021.09.012Crowder, M.J.: Maximum likelihood estimation for dependent observations. 38(1), 45–53 (1976). MR0403035Da Silva, J.L., Erraoui, M.: Mixed stochastic differential equations: Existence and uniqueness result. 31(2), 1119–1141 (2018). MR3803926. https://doi.org/10.1007/s10959-016-0738-9Ditlevsen, S., De Gaetano, A.: Mixed effects in stochastic differential equation models. 3(2), 137–153 (2005). MR2259358Dozzi, M., Mishura, Y., Shevchenko, G.: Asymptotic behavior of mixed power variations and statistical estimation in mixed models. 18(2), 151–175 (2015). MR3348583. https://doi.org/10.1007/s11203-014-9106-5Dufitinema, J., Pynnönen, S., Sottinen, T.: Maximum likelihood estimators from discrete data modeled by mixed fractional Brownian motion with application to the Nordic stock markets. 51(9), 5264–5287 (2022). MR4491681. https://doi.org/10.1080/03610918.2020.1764581Guerra, J., Nualart, D.: Stochastic differential equations driven by fractional Brownian motion and standard Brownian motion. 26(5), 1053–1075 (2008). MR2440915. https://doi.org/10.1080/07362990802286483Han, M., Xu, Y., Pei, B.: Mixed stochastic differential equations: averaging principle result. 112, 106705–7 (2021). MR4145479. https://doi.org/10.1016/j.aml.2020.106705.Isserlis, L.: On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables. 12(1/2), 134–139 (1918). https://doi.org/10.1093/biomet/12.1-2.134Kozachenko, Y., Melnikov, A., Mishura, Y.: On drift parameter estimation in models with fractional Brownian motion. 49(1), 35–62 (2015). MR3304366. https://doi.org/10.1080/02331888.2014.907294.Kubilius, K., Mishura, Y., Ralchenko, K.: . Bocconi & Springer Series, vol. 8, p. 390. Bocconi University Press, MilanSpringer, Cham (2017). MR3752152. https://doi.org/10.1007/978-3-319-71030-3Kukush, A., Lohvinenko, S., Mishura, Y., Ralchenko, K.: Two approaches to consistent estimation of parameters of mixed fractional Brownian motion with trend. 25(1), 159–187 (2022). MR4419677. https://doi.org/10.1007/s11203-021-09252-6Markussen, B.: Likelihood inference for a discretely observed stochastic partial differential equation. 9(5), 745–762 (2003). MR2047684. https://doi.org/10.3150/bj/1066418876Mishura, Y., Ralchenko, K., Shklyar, S.: Maximum likelihood drift estimation for gaussian process with stationary increments. 46(3–4), 67–78 (2017). https://doi.org/10.17713/ajs.v46i3-4.672Mishura, Y., Ralchenko, K., Shevchenko, G.: Existence and uniqueness of mild solution to stochastic heat equation with white and fractional noises. 98, 149–170 (2019). https://doi.org/10.1090/tpms/1068Mishura, Y.S., Shevchenko, G.M.: Rate of convergence of Euler approximations of solution to mixed stochastic differential equation involving Brownian motion and fractional Brownian motion. 19(4), 387–406 (2011). MR2871847. https://doi.org/10.1515/ROSE.2011.021Mishura, Y.S., Shevchenko, G.M.: Existence and uniqueness of the solution of stochastic differential equation involving Wiener process and fractional Brownian motion with Hurst index H>1/21/2$]]>. 40(19–20), 3492–3508 (2011). MR2860753. https://doi.org/10.1080/03610926.2011.581174Mishura, Y., Shevchenko, G.: Mixed stochastic differential equations with long-range dependence: Existence, uniqueness and convergence of solutions. 64(10), 3217–3227 (2012). MR2989350. https://doi.org/10.1016/j.camwa.2012.03.061Mishura, Y., Zili, M.: p. 194. ISTE Press, LondonElsevier Ltd, Oxford (2018). MR3793191Patrizio, C.R., Thompson, D.W.: Understanding the role of ocean dynamics in midlatitude sea surface temperature variability using a simple stochastic climate model. 35(11), 3313–3333 (2022). https://doi.org/10.1175/JCLI-D-21-0184.1Piterbarg, L., Rozovskii, B.: Maximum likelihood estimators in the equations of physical oceanography. In: , vol. 39, pp. 397–421. Birkhäuser Boston, Boston, MA (1996)Prakasa Rao, B.L.S.: Maximum likelihood estimation in the mixed fractional Vasicek model. 22 (2021). https://doi.org/10.1007/s41096-020-00094-8Schott, J.R.: , 2nd edn. Wiley Series in Probability and Statistics. Wiley-Interscience. John Wiley & Sons, Hoboken, NJ (2005)Shevchenko, G.: Mixed stochastic delay differential equations. 89, 181–195 (2014). MR3235184. https://doi.org/10.1090/S0094-9000-2015-00944-3Vergara, R.C.: Development of geostatistical models using stochastic partial differential equations. PhD thesis, MINES, Paris Tech (2018). http://cg.ensmp.fr/bibliotheque/public/CARRIZO_These_02513.pdfWiqvist, S., Golightly, A., McLean, A.T., Picchini, U.: Efficient inference for stochastic differential equation mixed-effects models using correlated particle pseudo-marginal algorithms. 157, 107151–26 (2021). MR4192029. https://doi.org/10.1016/j.csda.2020.107151Zili, M.: On the mixed fractional Brownian motion. , 32435–9 (2006). MR2253522. https://doi.org/10.1155/JAMSA/2006/32435