1 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–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_{x}^{H}}$:
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.
(1)
\[ \begin{array}{c}\displaystyle \left(\frac{\partial u}{\partial t}-\frac{1}{2}\Delta u\right)(t,x)=\sigma {\dot{B}_{x}^{H}},\hspace{1em}t>0,\hspace{0.2778em}x\in \mathbb{R},\\ {} \displaystyle u(0,x)=0,\hspace{1em}x\in \mathbb{R}.\end{array}\]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.
2 Preliminaries
Let $(\Omega ,\mathcal{F},\mathsf{P})$ be a complete probability space. Let ${B^{H}}=\left\{{B_{x}^{H}},x\in \mathbb{R}\right\}$ be a two-sided fractional Brownian motion with the Hurst index $H\in (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^{\lambda }}([a,b])$ the space of λ-Hölder continuous functions, $\lambda \in (0,1]$. If $f\in {C^{\alpha }}([a,b])$ and $g\in {C^{\beta }}([a,b])$ with $\alpha +\beta >1$, then the integral ${\textstyle\int _{a}^{b}}f(x)\hspace{0.1667em}dg(x)$ exists as the Riemann–Stieltjes integral
\[ {\int _{a}^{b}}f(x)\hspace{0.1667em}dg(x)=\underset{\left|\pi \right|\to 0}{\lim }\sum \limits_{i}f({x_{i}^{\ast }})\big(g({x_{i+1}})-g({x_{i}})\big),\]
where $\pi =\{a={x_{0}}\le {x_{0}^{\ast }}\le {x_{1}}\le \cdots \le {x_{n-1}}\le {x_{n-1}^{\ast }}\le {x_{n}}=b\}$, $\left|\pi \right|={\max _{i}}\left|{x_{i+1}}-{x_{i}}\right|$, see, e. g., [40, Thm. 4.2.1]. Moreover, the integral ${\textstyle\int _{a}^{b}}g(x)\hspace{0.1667em}df(x)$ also exists, and the integration-by-parts formula holds:
\[ {\int _{a}^{b}}f(x)\hspace{0.1667em}dg(x)=f(b)g(b)-f(a)g(a)-{\int _{a}^{b}}g(x)\hspace{0.1667em}df(x)\]
(see [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^{\beta }}([a,b])$ a. s., for any $\beta \in (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\in \mathbb{R}$, the function $G(t-s,x-\cdot )$ is (globally) Lipschitz, i. e., it belongs to ${C^{\alpha }}([a,b])$ with $\alpha =1$. This means that for almost all $\omega \in \Omega $, the integral ${\textstyle\int _{a}^{b}}G(t-s,x-y)\hspace{0.1667em}d{B_{y}^{H}}$ exists in the Riemann–Stieltjes sense, and
It is known (see, e. g., the proof of Prop. A.1 in [9]) that for all $\gamma >H$, $\frac{{B_{x}^{H}}}{{\left|x\right|^{\gamma }}}\to 0$ a. s., as $\left|x\right|\to \infty $. Therefore, $G(t-s,x-y)\hspace{0.1667em}{B_{y}^{H}}\to 0$, a. s., as $y\to \pm \infty $. Hence,
\[\begin{aligned}{}{\int _{a}^{b}}G(t-s,x-y)d{B_{y}^{H}}=G(t-s,x-b){B_{b}^{H}}& -G(t-s,x-a){B_{a}^{H}}\\ {} & +{\int _{a}^{b}}{B_{y}^{H}}{G^{\prime }_{2}}(t-s,x-y)\hspace{0.1667em}dy,\end{aligned}\]
where ${G^{\prime }_{2}}$ denotes the partial derivative of G with respect to the spatial variable:
(4)
\[ {G^{\prime }_{2}}(t,x)=\frac{\partial }{\partial x}G(t,x)=-\frac{x}{\sqrt{2\pi }{t^{3/2}}}\exp \left\{-\frac{{x^{2}}}{2t}\right\},\hspace{1em}t>0,\hspace{0.2778em}x\in \mathbb{R}.\]
\[ {\int _{-\infty }^{+\infty }}G(t-s,x-y)d{B_{y}^{H}}={\int _{-\infty }^{+\infty }}{B_{y}^{H}}{G^{\prime }_{2}}(t-s,x-y)\hspace{0.1667em}dy,\]
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.3 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.
3.1 Covariance and variance. Exact formulae
We start with deriving explicit expressions for the variance and covariance of $u(t,\cdot )$. The asymptotic behavior of the covariance function will be investigated in the next subsection.
Proposition 1.
For fixed $t\in [0,T]$, $u(t,\cdot )$ is a stationary Gaussian process with covariance function
(6)
\[\begin{aligned}{}R(t,x)& :=\operatorname{cov}\big(u(t,0),u(t,x)\big)\\ {} & =-\frac{{\sigma ^{2}}}{2}{\int _{0}^{t}}\hspace{-0.1667em}\hspace{-0.1667em}{\int _{0}^{t}}\hspace{-0.1667em}\hspace{-0.1667em}{\int _{\mathbb{R}}}\hspace{-0.1667em}{\int _{\mathbb{R}}}{G^{\prime }_{2}}(t-s,-y){G^{\prime }_{2}}(t-r,x-v)|y-v{|^{2H}}\hspace{0.1667em}dr\hspace{0.1667em}ds\hspace{0.1667em}dv\hspace{0.1667em}dy.\end{aligned}\]The next result gives a simpler expression for the covariance function of the solution. This expression contains a single integral over $\mathbb{R}$ instead of double integral over ${\mathbb{R}^{2}}$.
Proposition 2.
The covariance function $R(t,x)$ can be represented in the form
where ${w^{\langle \alpha \rangle }}:={\left|w\right|^{\alpha }}\operatorname{sign}w$.
(7)
\[\begin{aligned}{}R(t,x)=\frac{{\sigma ^{2}}H}{\sqrt{2\pi }}{\int _{0}^{t}}\hspace{-0.1667em}\hspace{-0.1667em}& {\int _{0}^{t}}{(2t-s-r)^{-\frac{3}{2}}}\\ {} & \times {\int _{\mathbb{R}}}{w^{\langle 2H-1\rangle }}(w-x)\exp \left\{-\frac{{(w-x)^{2}}}{2(2t-s-r)}\right\}\hspace{0.1667em}dw\hspace{0.1667em}ds\hspace{0.1667em}dr,\end{aligned}\]Remark 1.
In the case $H>\frac{1}{2}$, it is possible to integrate by parts once more and rewrite the formula for $R(t,x)$ in the form
\[\begin{aligned}{}R(t,x)=\frac{{\sigma ^{2}}H(2H-1)}{\sqrt{2\pi }}{\int _{0}^{t}}\hspace{-0.1667em}\hspace{-0.1667em}& {\int _{0}^{t}}{(2t-s-r)^{-\frac{1}{2}}}\\ {} & \times {\int _{\mathbb{R}}}{\left|w\right|^{2H-2}}\exp \left\{-\frac{{(w-x)^{2}}}{2(2t-s-r)}\right\}\hspace{0.1667em}dw\hspace{0.1667em}ds\hspace{0.1667em}dr.\end{aligned}\]
3.2 Upper bound for covariance function and its asymptotic behavior
Our next goal is to establish ergodicity of $u(t,\cdot )$. Since $u(t,\cdot )$ is a stationary Gaussian process (by Proposition 1), it suffices to show that $R(t,x)\to 0$, as $x\to \infty $. The following proposition plays a crucial role in the proof of ergodicity.
Proposition 4 implies that the covariance function $R(t,x)$ of the solution $u(t,x)$ vanishes as $x\to +\infty $. Since $u(t,\cdot )$ is a stationary Gaussian process, this yields the following result.
4 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 $\delta >0$, the random field u given by (3) is observed at the points ${x_{k}}=k\delta $, $k=1,\dots ,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\in \mathbb{R}$ is strictly stationary and ergodic. Therefore, for any Borel function $g:\mathbb{R}\to \mathbb{R}$ such that $\mathsf{E}|g(u(t,0))|<\infty $, thanks to ergodic theorem, it holds that
This gives the idea to consider the following estimator for ${\sigma ^{2}}$:
where
see Proposition 3. Taking into account (9), we have the following theorem.
(9)
\[ \frac{1}{N}{\sum \limits_{k=1}^{N}}g(u(t,{x_{k}}))\to \mathsf{E}\big[g(u(t,0))\big],\hspace{1em}\text{a.s., as}\hspace{2.5pt}N\to \infty .\]5 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 right-hand side of (3) multiplicatively, it suffices to consider the case $\sigma =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 $\Delta x={2^{-7}}$. First, we integrate Green’s function numerically with respect to s, then we calculate the stochastic integral with respect to the fractional Brownian motion ${B_{y}^{H}}$ replacing it by the sum (we integrate over $[x-4,x+4]$ instead of $\mathbb{R}$).
Table 1.
Means and standard deviations of the estimator ${\hat{\sigma }_{N}^{2}}$ for the observations step $\delta =1$
N | ${2^{4}}$ | ${2^{5}}$ | ${2^{6}}$ | ${2^{7}}$ | ${2^{8}}$ | ${2^{9}}$ | ${2^{10}}$ | |
H = 0.1 | Mean | 1.0470 | 1.0120 | 1.0017 | 0.9981 | 1.0073 | 1.0132 | 1.0135 |
Std. dev | 0.3880 | 0.2619 | 0.1812 | 0.1254 | 0.0892 | 0.0658 | 0.0455 | |
H = 0.3 | Mean | 0.9991 | 1.0161 | 0.9994 | 0.9932 | 0.9992 | 0.9964 | 0.9951 |
Std. dev | 0.4326 | 0.2920 | 0.1928 | 0.1467 | 0.1088 | 0.0699 | 0.0502 | |
H = 0.5 | Mean | 1.0051 | 0.9613 | 0.9576 | 0.9883 | 0.9833 | 0.9893 | 0.9915 |
Std. dev | 0.4292 | 0.2601 | 0.2009 | 0.1343 | 0.0858 | 0.0731 | 0.0479 | |
H = 0.7 | Mean | 0.9476 | 0.9532 | 0.9916 | 1.0299 | 1.0187 | 1.0038 | 1.0049 |
Std. dev | 0.5587 | 0.4144 | 0.2941 | 0.2363 | 0.1780 | 0.1161 | 0.0775 | |
H = 0.9 | Mean | 0.9910 | 1.0163 | 0.9722 | 0.9800 | 0.9588 | 0.9543 | 0.9602 |
Std. dev | 0.9989 | 0.8766 | 0.7046 | 0.5660 | 0.4976 | 0.4934 | 0.3416 |
Table 2.
Means and standard deviations of the estimator ${\hat{\sigma }_{N}^{2}}$ for the observations step $\delta =\frac{1}{2}$
N | ${2^{4}}$ | ${2^{5}}$ | ${2^{6}}$ | ${2^{7}}$ | ${2^{8}}$ | ${2^{9}}$ | ${2^{10}}$ | |
H = 0.1 | Mean | 1.1065 | 1.0451 | 1.0291 | 1.0155 | 1.0182 | 1.0325 | 1.0268 |
Std. dev | 0.6812 | 0.5189 | 0.3806 | 0.2639 | 0.1821 | 0.1363 | 0.0948 | |
H = 0.3 | Mean | 1.0415 | 1.0292 | 1.0489 | 1.0109 | 1.0064 | 1.0150 | 1.0078 |
Std. dev | 0.7277 | 0.5332 | 0.3749 | 0.2667 | 0.1944 | 0.1325 | 0.0875 | |
H = 0.5 | Mean | 1.0695 | 0.9728 | 1.0597 | 1.0671 | 1.0327 | 1.0205 | 1.0079 |
Std. dev | 0.9274 | 0.6440 | 0.4363 | 0.2916 | 0.1742 | 0.1215 | 0.1086 | |
H = 0.7 | Mean | 0.9555 | 0.9919 | 0.9982 | 0.9838 | 0.9821 | 1.0414 | 1.0228 |
Std. dev | 0.9722 | 0.7945 | 0.5504 | 0.3966 | 0.2896 | 0.2088 | 0.1417 | |
H = 0.9 | Mean | 0.9763 | 0.9302 | 0.9012 | 0.9731 | 0.9667 | 0.9365 | 0.9323 |
Std. dev | 1.3018 | 1.1639 | 0.8602 | 0.7813 | 0.7033 | 0.5201 | 0.4076 |
Table 3.
Means and standard deviations of the estimator ${\hat{\sigma }_{N}^{2}}$ for the observations step $\delta =\frac{1}{4}$
N | ${2^{4}}$ | ${2^{5}}$ | ${2^{6}}$ | ${2^{7}}$ | ${2^{8}}$ | ${2^{9}}$ | ${2^{10}}$ | |
H = 0.1 | Mean | 0.8753 | 1.0704 | 1.1010 | 1.0429 | 1.0229 | 0.9875 | 0.9893 |
Std. dev | 0.8769 | 0.9028 | 0.6155 | 0.4187 | 0.3364 | 0.2298 | 0.1550 | |
H = 0.3 | Mean | 1.1555 | 1.0970 | 1.0627 | 0.9430 | 0.9469 | 0.9717 | 0.9923 |
Std. dev | 1.0885 | 0.8464 | 0.6896 | 0.4895 | 0.3406 | 0.2733 | 0.1872 | |
H = 0.5 | Mean | 0.9159 | 0.9157 | 0.9598 | 0.9836 | 1.0244 | 1.0295 | 1.0139 |
Std. dev | 1.2274 | 0.8237 | 0.6379 | 0.5033 | 0.4382 | 0.3087 | 0.2228 | |
H = 0.7 | Mean | 0.9758 | 0.9743 | 0.9323 | 0.8939 | 0.9806 | 1.0009 | 1.0069 |
Std. dev | 1.3749 | 1.1514 | 0.9001 | 0.6895 | 0.5734 | 0.4143 | 0.3277 | |
H = 0.9 | Mean | 0.6687 | 0.6877 | 0.7213 | 0.7765 | 0.8427 | 0.8468 | 0.9224 |
Std. dev | 0.8962 | 0.8343 | 0.7729 | 0.7604 | 0.7558 | 0.6637 | 0.6916 |
We study the performance of the estimator ${\hat{\sigma }_{N}^{2}}$ for various values of the Hurst index H and three values of the observation step $\delta =1,0.5,0.25$. The means and standard deviations of the estimates are reported in Tables 1–3. We see that the estimates converge to the true value of ${\sigma ^{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 $\delta =1$. Probably, the horizon of observations $T=N\delta $ is more important for the quality of the estimators than the step δ.