Ergodic properties of the solution to a fractional stochastic heat equation, with an application to diffusion parameter estimation

The paper deals with a stochastic heat equation driven by an additive fractional Brownian space-only noise. We prove that a solution to this equation is a stationary and ergodic Gaussian process. These results enable us to construct a strongly consistent estimator of the diffusion parameter.


Introduction
In the last four decades, there have been enormous advances in the study of random field solutions to stochastic partial differential equations (SPDEs) driven by general Brownian noises. The starting point of this theory was the seminal work by Walsh [38]. Most of the results published till now are mainly focused on the analysis of heat and wave equations perturbed by Gaussian white noises in time with a fairly general spatial correlation [3,8,16,17,31,33]. At the same time some publications are devoted to SPDEs driven by fractional-type noises [20,22,28,30,[34][35][36].
The interest in equations of such type is mainly due to the development of analysis and the theory of random processes on one hand, and various applications to some real phenomena and real situations on the other hand. For instance, various types of SPDEs provide suitable models in the study of population growth [18], some climate and oceanographic phenomena [1,23], mathematical finance [7].
Consequently, statistical methods are required to calibrate SPDE models from given observations. In particular, a problem of parameter estimation based on discrete observations of a solution to an SPDE attracted considerable interest very recently. It was first investigated in [29]. Later, applying similar methods, parabolic SPDEs including the stochastic heat equation had been studied in [6,10,14]. In these articles estimators based on power variations of time-increments of the solution were constructed and central limit theorems were proved. In [24] an adaptive maximum likelihood type estimator of the coefficient parameter was proposed for a parabolic linear second order SPDE. However, there are still basic questions which are not settled in the statistical literature on SPDEs, see the recent review [13] for details. This paper deals with the stochastic heat equation, which is a typical example of SPDE. This equation and its numerical approximation has been extensively studied in [11,19,21,25,37]. Also, in [5,28] estimators for the drift parameter based on power variation of the solution to the fractional stochastic heat equation was constructed.
Specifically, this paper is devoted to the following stochastic heat equation with a fractional Brownian noise B H x : We prove the stationarity and ergodicity of its solution u(t, x) as a function of the spatial variable x by analyzing the behavior of the covariance function. Then these properties are applied for construction of a strongly consistent estimator for the unknown diffusion parameter σ . We extend the results of [2], where a similar problem for SPDE with white noise was studied. We would like to emphasize that in the present paper we study a stochastic heat equation driven by a fractional space-only noise. Unlike SPDEs with space-time noise, equations with purely spatial Gaussian noise are studied not so extensively in the statistical literature. We can mention only the papers [2,12,15], where the parameter estimation for equations with space-only white noise was investigated. At the same time, such type of noise is an important type of stationary perturbations, see discussion and examples in [27].
Concerning the sampling scheme, we assume that the solution u(t, x) is observed at equidistant spatial points for a fixed time. On one hand, in many practical applications the solution indeed is observed only at some discrete space points; e. g. temperature of a heated body, velocity of a turbulent flow, instantaneous forward rates where the space variable corresponds to time until maturity [14], see also [29]. On the other hand, it turns out that in order to estimate the diffusion coefficient, it is enough to observe the solution at one time instant. Such situation is quite common for statistical inference for SPDEs, see [13,Section 3]. However, it is possible to incorporate the additional information of observing the solution discretely in time by taking the (weighted) average of the estimators similarly to [6] or [14].
The paper is organized as follows. In Section 2 we recall the definition of mild solution to SPDE (1) and the required theory of stochastic integration with respect to fractional Brownian motion. In Section 3 we calculate the covariance and variance of the solution and establish its stationarity and ergodicity. Based on these results, in Section 4 we construct a strongly consistent estimator of the parameter σ by discrete observations of a trajectory of the solution. Numerical results are presented in Section 5. All proofs are collected in Appendix A.

Preliminaries
Let ( , F, P) be a complete probability space. Let B H = B H x , x ∈ R be a twosided fractional Brownian motion with the Hurst index H ∈ (0, 1), that is, a centered Gaussian process, starting at 0, with the covariance function Following [38], we define a solution to SPDE (1) by where G is Green's function of the heat equation: Due to the Hölder properties of a fractional Brownian motion and Green's function, the integral with respect to the fractional Brownian motion in (3) exists as the pathwise Riemann-Stieltjes integral. Let us briefly recall some basic facts concerning the integration of Hölder continuous functions.
Let a < b. Denote by C λ ([a, b]) the space of λ-Hölder continuous functions,  [40, Thms. 4.2.1 and 3.1]). It is well known that in the case of Hölder continuous functions the Riemann-Stieltjes integral coincides with Young's integral [39] and with the generalized Lebesgue-Stieltjes integral [40].
This theory can be easily applied to the pathwise stochastic integration with respect to a fractional Brownian motion, since the sample paths of B H belong to C β ([a, b]) a. s., for any β ∈ (0, H ) (see, e. g., [32,Sec. 1.16]). On the other hand, it is not hard to see that for any fixed t > 0, s > 0 and x ∈ R, the function where G 2 denotes the partial derivative of G with respect to the spatial variable: It is known (see, e. g., the proof of Prop. A.1 in [9]) that for all γ > H, and the solution (3) can be written in the form This form of the solution does not contain a stochastic integral with respect to a fractional Brownian motion, hence, it is more convenient for further calculations.

Properties of solution
In this section we will investigate properties of the solution u related to its covariance structure. In particular, we will establish stationarity and ergodicity.

Covariance and variance. Exact formulae
We start with deriving explicit expressions for the variance and covariance of u(t, ·). The asymptotic behavior of the covariance function will be investigated in the next subsection.

Proposition 1. For fixed t ∈ [0, T ], u(t, ·) is a stationary Gaussian process with covariance function
The next result gives a simpler expression for the covariance function of the solution. This expression contains a single integral over R instead of double integral over R 2 .

Proposition 2. The covariance function R(t, x) can be represented in the form
dw ds dr, (7) where w α := |w| α sign w.

Remark 1.
In the case H > 1 2 , it is possible to integrate by parts once more and rewrite the formula for R(t, x) in the form dw ds dr.
where denotes the gamma function.

Upper bound for covariance function and its asymptotic behavior
Our next goal is to establish ergodicity of u(t, ·). Since u(t, ·) is a stationary Gaussian process (by Proposition 1), it suffices to show that R(t, x) → 0, as x → ∞. The following proposition plays a crucial role in the proof of ergodicity.

Proposition 4.
For t > 0 and x > 0, the covariance function R(t, x) admits the following upper bound: where C H is a positive constant depending only on H .
Proposition 4 implies that the covariance function R(t, x) of the solution u(t, x) vanishes as x → +∞. Since u(t, ·) is a stationary Gaussian process, this yields the following result. Corollary 1. For a fixed t > 0, the random process {u(t, x), x ∈ R} is ergodic.

Diffusion parameter estimation
In this section we apply previous results to the following statistical problem. Assume that for a fixed time t > 0 and a fixed step δ > 0, the random field u given by (3) is observed at the points x k = kδ, k = 1, . . . , N. Our aim is to construct a strongly consistent estimator for σ based on these observations. By the results of the previous section, the field u(t, x), x ∈ R is strictly stationary and ergodic. Therefore, for any Borel function g : R → R such that E|g(u(t, 0))| < ∞, thanks to ergodic theorem, it holds that u(t, 0)) , a. s., as N → ∞.
This gives the idea to consider the following estimator for σ 2 : see Proposition 3. Taking into account (9), we have the following theorem.
Theorem 1.σ 2 N is a strongly consistent estimator for the parameter σ 2 as N → ∞. Remark 3. The Hurst parameter H ∈ (0, 1) is assumed to be known. It can be estimated independently of σ with the help of quadratic variations, see [4,26].

Simulation
In this section we illustrate the quality of the estimator with the help of simulation experiments. We generate the trajectories of the solution u(t, x) to the equation (1) by the discretization of the formula (3). Since the parameter σ enters into the righthand side of (3) multiplicatively, it suffices to consider the case σ = 1 only. Also we choose t = 1 as the observation time. For each value of the Hurst index H , we simulate 100 sample paths of the solution u (1, x), using the partition with the step x = 2 −7 . First, we integrate Green's function numerically with respect to s, then   Tables 1-3. We see that the estimates converge to the true value of σ 2 , and their standard deviations tend to zero. Hence these simulations confirm the consistency. However the rate of convergence for H = 0.9 is not very high. Also we see that the best results are obtained for δ = 1. Probably, the horizon of observations T = Nδ is more important for the quality of the estimators than the step δ.

A Proofs
Proof of Proposition 1. Using the formula (5) for the solution, we can write the covariance in the form Inserting the explicit expression (2) for the covariance function of a fractional Brownian motion, we obtain cov u(t, z), u(t, z + x) The term I 1 can be represented as a product of two integrals: It is not hard to see that Similarly, one can prove that I 2 = 0. Thus, we have Finally, by the change of variables y = y − z, v = v − z, we obtain cov u(t, z), u(t, z + x) Since the expression does not depend on z, the process u(t, ·) is stationary. Moreover, we see that the formula (6) holds.

Proof of Proposition 2.
In order to obtain (7), we need to transform the inner integral in (6), namely By the substitution v − y = w, we have Using the explicit expression for G 2 , see (4), we get Let us transform the sum under the exponent: .
Hence, (11) is rewritten in the form The inner integral can be computed by applying the following equalities for the moments of a normal random variable: if ξ ∼ N (μ, ς 2 ), then Combining this with (12), we obtain Let us transform each integral separately, applying integration by parts. First, we rewrite I − (s, r) as follows: dw.
Now let us transform the second integral in the same way: Thus, Taking into account that, by (6) and (10), R(t, x) = − σ 2 2 t 0 t 0 I (s, r) ds dr, we conclude the proof.
Proof of Proposition 3. By Propositions 1 and 2, we have that dw ds dr. (13) Let us consider the inner integral. Taking into account that the integrand is an even function, and using the substitution Inserting this expression into (13), we arrive at It remains to compute the integral Combining this with (14), we get the announced equality for E u(t, x) 2 .
In order to prove Proposition 4, we start with an auxiliary result.

Lemma 1.
For any α > 0 and β > 0, there exists a constant C > 0 such that for all t > 0 and x > 0, Proof. By the substitution t−s Note that the function f (x) = x −α exp − β x is continuous on (0, +∞), and . Consequently, Proof of Proposition 4. In this proof C will denote a generic positive constant that may depend on H ; its value is not important and may change from one line to another. According to Proposition 2, dw ds dr dw ds dr (16) Let us consider each term separately. In order to bound J 1 , we write because exp 2wx 2(2t−s−r) ≤ 1 when w ≤ 0 and x ≥ 0. Further, we transform the integral in the right-hand side of (18), using the substitution z = w 2 2(2t−s−r) , w = √ 2(2t − s − r)z, dw = √ 2(2t − s − r) 1 2 √ z dz. We obtain Hence, by (18), Inserting this bound into (15), we arrive at ds dr.
The integrals in the right-hand side can be bounded with the help of Lemma 1. We obtain |J 1 | ≤ Cσ 2 t 2 x 2H −2 .