Eigenfunction martingale estimating functions and filtered data for drift estimation of discretely observed multiscale diffusions.

We propose a novel method for drift estimation of multiscale diffusion processes when a sequence of discrete observations is given. For the Langevin dynamics in a two-scale potential, our approach relies on the eigenvalues and the eigenfunctions of the homogenized dynamics. Our first estimator is derived from a martingale estimating function of the generator of the homogenized diffusion process. However, the unbiasedness of the estimator depends on the rate with which the observations are sampled. We therefore introduce a second estimator which relies also on filtering the data, and we prove that it is asymptotically unbiased independently of the sampling rate. A series of numerical experiments illustrate the reliability and efficiency of our different estimators.

Introduction

Learning models from data is a problem of fundamental importance in modern applied mathematics. The abundance of data in many application areas such as molecular dynamics, atmosphere/ocean science makes it possible to develop physics-informed data driven methodologies for deriving models from data (Raissi et al. 2019; Yang et al. 2021; Zhang et al. 2019). Naturally, most problems of interest are characterised by a very high-dimensional state space and by the presence of many characteristic length and time scales. When it is possible to decompose the state space into the resolved and unresolved degrees of freedom, then one is usually interested in the derivation of a model for the resolved degrees of freedom, while treating the unresolved scales as noise. Clearly, these reduced models are stochastic, often described by stochastic differential equations (SDEs). The goal of this paper is to derive rigorous and systematic methodologies for learning coarse-grained models that accurately describe the dynamics at macroscopic length and time scales from noisy observations of the full, unresolved dynamics. We apply the proposed methodologies to simple models of fast/slow SDEs for which the theory of homogenization exists, that enables us to study the inference problem in a rigorous and systematic manner.

In many applications, the available data are noisy, not equidistant and certainly not compatible with the coarse-grained model. The presence of observation noise and of the model-data mismatch renders the problem of learning macroscopic models from microscopic data highly ill-posed. Several examples from econometrics (market microstructure noise) (Aït-Sahalia et al. 2006) and molecular dynamics show that standard algorithms, e.g. maximum likelihood or quadratic variation for the diffusion coefficient, are asymptotically biased and they fail to estimate correctly the parameters in the coarse-grained model. In a series of earlier works, this problem was studied using maximum likelihood techniques with subsampled data (Pavliotis and Stuart 2007; Papavasiliou et al. 2009), methodologies based on the method of moments (Krumscheid et al. 2013; Kalliadasis et al. 2015; Krumscheid et al. 2015), quadratic programming approaches (Crommelin and Vanden-Eijnden 2006a) as well as Bayesian approaches (Abdulle and Di Blasio 2020; Abdulle et al. 2020). We also mention the pioneering work on estimating the integrated stochastic volatility in the presence of market microstructure noise (Aït-Sahalia et al. 2006; Zhang et al. 2005). In particular, in Aït-Sahalia and Jacod (2014) the authors analyse the correct interplay between the intensity of the microstructure noise and the optimal rates of convergence.

The main observation in Pavliotis and Stuart (2007); Papavasiliou et al. (2009) is that when the maximum likelihood estimator (MLE) of the fast/slow system is evaluated at the full data, then the MLE becomes asymptotically biased; in fact, the original data are not compatible with the homogenized equation, and therefore data need to be preprocessed, for instance under the form of subsampling. On the other hand, when the MLE is evaluated at appropriately subsampled data, then it becomes asymptotically unbiased. Although this is an interesting theoretical observation (see also later developments in Spiliopoulos and Chronopoulou (2013)), it does not lead to an efficient algorithm. The reason for this is that the performance of the estimator depends very sensitively on the choice of the sampling rate. In addition, the optimal sampling rate is not known and is strongly dependent on the problem under investigation. Furthermore, subsampling naturally leads to an increase in the variance, unless appropriate variance reduction methodologies are used.

In a recent work (Abdulle et al. 2021), we addressed the problem of lack of robustness of the MLE with subsampling algorithm by introducing an appropriate filtering methodology that leads to a stable and robust algorithm. In particular, rather than subsampling the original trajectory, we smoothed the data by applying an appropriate linear time-invariant filter from the exponential family and we modified the MLE by inserting the new filtered data. This new estimator was thus independent of the subsampling rate and also asymptotically unbiased and robust with respect to the parameters of the filter.

However, the assumption that the full path of the solution is observed is not realistic in most applications. In fact, in all real problems one can only obtain discrete measurements of the diffusion process. Hence, in this paper we focus on the problem of learning the coarse-grained homogenized model assuming that we are given discrete observations from the microscopic model. In this paper, we use the martingale estimating functions that were introduced in Bibby and Sø rensen (1995), where the authors study drift estimation for discrete observations of one-scale processes and show that estimators based on the discretized continuous-version likelihood function can be strongly biased. They therefore propose martingale estimating functions obtained by adjusting the discretized continuous-version score function by its compensator which leads to unbiased estimators. Moreover, in Kessler and Sørensen (1999) a different type of martingale estimating function, which is dependent on the eigenvalues and eigenfunctions of the generator of the stochastic process, is introduced and asymptotic unbiasedness and normality are proved. Furthermore, another inference methodology that uses spectral information is proposed in Crommelin and Vanden-Eijnden (2006b). Their approach consists of inferring the drift and diffusion functions of a diffusion process by minimizing an objective function which measures how close the generator is to having a reference spectrum which is obtained from the time series through the construction of a discrete-time Markov chain. This idea has been further expanded in several directions in Crommelin and Vanden-Eijnden (2011).

In this paper, we propose a new estimator for learning homogenised SDEs from noisy discrete data that is based on the martingale estimators that were introduced in Kessler and Sørensen (1999). The main idea is to consider the eigenvalues and eigenfunctions of the generator of the homogenized process. This new estimator is asymptotically unbiased only if the distance between two consecutive observations is not too small compared with the multiscale parameter describing the fastest scale, i.e. if data are compatible with the homogenized model. Therefore, in order to obtain unbiased approximations independently of the sampling rate with which the observations are obtained, we propose a second estimator which, in addition to the original observations, relies also on filtered data obtained following the filtering methodology presented in Abdulle et al. (2021). We observe that smoothing the original data makes observations compatible with the homogenized process independently of the rate with which they are sampled and hence this second estimator gives a black-box tool for parameter estimation.

Our main contributions

The main goal of this paper is to propose new algorithms based on martingale estimating functions and filtered data for which we can prove rigorously that they are asymptotically unbiased and not sensitive with respect to, e.g. the sampling rate and the observation error. In particular, we combine two main ideas:We prove theoretically and observe numerically that the estimator without filtered data is asymptotically unbiased if:Moreover, we show that the estimator with filtered data corrects the bias caused by a sampling rate smaller than the multiscale parameter, and therefore, it is asymptotically unbiased independently of the sampling rate.

the use of martingale estimating functions for discretely observed diffusion processes based on the eigenvalues and the eigenfunctions of the generator of the homogenized process, which was originally presented for one-scale problems in Kessler and Sørensen (1999);

the filtering methodology for smoothing the data in order to make them compatible with the homogenized model, which was introduced in Abdulle et al. (2021).

the observations are taken at the homogenized regime, i.e. the sampling rate is independent of the parameter measuring scale separation;

the observations are taken at the multiscale regime, i.e. the sampling rate is dependent on the fastest scale, and the sampling rate is bigger than the multiscale parameter.

Outline. The rest of the paper is organized as follows. In Sect. 2, we present the Langevin dynamics and its corresponding homogenized equation and we introduce the two proposed estimators based on eigenvalues and eigenfunctions of the generator with and without filtered data. In Sect. 3, we present the main results of this work, i.e. the asymptotic unbiasedness of the two estimators, and in Sect. 4 we perform numerical experiments which validate the efficacy of our methods. Section 4.1 is devoted to the proof of the main results which are presented in Sect. 3. Finally, in Appendix, we show some technical results which are employed in the analysis and we explain some details about the implementation of the proposed methodology.

Problem setting

In this work, we study the following class of multiscale diffusion processes. Consider the following two-scale SDE, observed over the time interval [0, T]where describes the fast scale, and are, respectively, the drift and diffusion coefficients and is a standard one-dimensional Brownian motion. The functions and are the slow-scale and fast-scale parts of the potential, and they are assumed to be smooth. Moreover, we also assume p to be periodic with period . We remark that our setting can be considered as a semi-parametric framework similar to the one of Krumscheid et al. (2013). The components of the potential function V, in fact, can be viewed as basis functions for a truncated expansion (e.g. Taylor series or Fourier expansion) of the unknown slow-scale potential , where the components of the unknown drift term contain the generalized Fourier coefficients, i.e.We also mention that assuming a parametric form for the potential V is a technique usually employed in the statistics literature in order to regularize the likelihood function and obtain a parametric approximation of the actual MLE of V, which does not exist in general (Pokern et al. 2009).

Remark 2.1

For clarity of the presentation, we focus our analysis on scalar multiscale diffusions with a finite number of parameters in the drift that have to be learned from data. Nevertheless, we remark that all the following theory can be generalized to the case of multidimensional diffusion processes in , for which we provide further details in Appendix C and an example in Sect. 4.5. However, the problem becomes more complex and computationally expensive from a numerical viewpoint and it can be prohibitive if the dimension d is too large, since the methodology proposed in this paper requires the solution of the eigenvalue problem for the generator of a d-dimensional diffusion process.

The theory of homogenization (see, for example, Bensoussan et al. 2011, Chapter 3 or Pavliotis and Stuart 2008, Chapter 18) guarantees the existence of the following homogenized SDE whose solution is the limit in law of the solutions of (2.1) as random variables in where , . The coefficient has the explicit formulawithand where the function is the unique solution with zero-mean with respect to the measure of the differential equationendowed with periodic boundary conditions. In particular, for one-dimensional diffusion processes, we havewhich impliesOur goal is to derive estimators for the homogenized drift coefficient A based on multiscale data originating from (2.1). In this work, we consider the same setting as Abdulle et al. (2021), which is summarized by the following assumption.

Assumption 2.2

The potentials p and V satisfy

and is L-periodic for some ,

and each component is polynomially bounded from above and bounded from below, and there exist such that

is Lipschitz continuous, i.e. there exists a constant such that

Let us remark that, under Assumption 2.2, it has been proved in Pavliotis and Stuart (2007) that both processes (2.1) and (2.2) are geometrically ergodic and their invariant measure has a density with respect to the Lebesgue measure. In particular, let us denote by and the densities of the invariant measures of and , respectively, defined byand

Remark 2.3

The value of the initial condition in the SDE (2.1) is important neither for the numerical experiments nor for the following analysis and can be chosen arbitrarily. In fact, the process is geometrically ergodic and therefore it converges to its invariant distribution with density exponentially fast for any initial condition.

Drift estimation problem. Consider uniformly distributed observation times , set and let be a realization of the solution of (2.1). We then assume to know a sample of the realization where , and we aim to estimate the drift coefficient A of the homogenized equation (2.2). First, since we deal with discrete observations of stochastic processes, we employ martingale estimating functions based on eigenfunctions, which have already been studied for problems without a martingale structure in Kessler and Sørensen (1999). Second, by observing that if the time-step is too small with respect to the multiscale parameter , then the data could be compatible with the full dynamics rather than with the coarse-grained model, we also adopt the filtering methodology presented in Abdulle et al. (2021), which has been proved to be beneficial for correcting the behaviour of the maximum likelihood estimator (MLE) in the setting of continuous observations.

Martingale estimating functions based on eigenfunctions

We first remark that a general theory for martingale estimating functions exists and is thoroughly outlined in Bibby and Sø rensen (1995). They appear to be appropriate for multiscale problems due to their robustness properties. In this paper, we develop martingale estimating functions based on the eigenfunctions of the generator of the process, since the theory of the eigenvalue problem for elliptic differential operators and the multiscale analysis of this eigenvalue problem are well developed. Let be the set of admissible drift coefficients for which Assumption 2.2(ii) is satisfied. To describe our methodology, we consider the solution of the homogenized process (2.2) with a generic parameter instead of the exact drift coefficient A:which, according to (2.6), has invariant measureThe generator of (2.7) is defined for all as:where the subscript denotes the dependence of the generator on the unknown drift coefficient a. From the well-known spectral theory of diffusion processes and under our assumptions on the potential V, we deduce that has a countable set of eigenvalues (see, for example, Hansen et al. 1998). In particular, let be the sequence of eigenvalue-eigenfunction couples of the generator which solve the eigenvalue problemwhich, due to (2.9), is equivalent toand where the eigenvalues satisfy and the eigenfunctions form an orthonormal basis for the weighted space . We mention in passing that, by making a unitary transformation, the eigenvalue problem for the generator of the Langevin dynamics can be transformed to the standard Sturm–Liouville problem for Schrödinger operators (Pavliotis 2014, Chapter 4). We now state a formula, which has been proved in Kessler and Sørensen (1999) and will be fundamental in the rest of the paperwhere is the constant distance between two consecutive observations. We now discuss how this eigenvalue problem can be used for parameter estimation. Let J be a positive integer and let be J arbitrary functions possibly dependent on the parameter a, which satisfy Assumption 2.5(i)(ii) stated below, and define for the martingale estimating functionThen, given a set of observations , we consider the score function defined byThis function can be seen as an approximation in terms of eigenfunctions of the true score function, i.e. the gradient of the log-likelihood function with respect to the unknown parameter. The full derivation of a martingale estimating function as an approximation of the true score function is given in detail in Bibby and Rensen (1995, Sect. 2). The first step is a discretization of the gradient of the continuous-time log-likelihood, which yields a biased estimating function. Hence, the next step is adjusting this function by adding its compensator in order to obtain a zero-mean martingale. Moreover, by using the eigenfunctions of the generator, it is shown in Kessler and Sørensen (1999) that this approach is suitable for scalar diffusion processes with no multiscale structure, i.e. processes with a single characteristic length/time scale. In fact, by a classical result for ergodic diffusion processes (Pavliotis 2014, Sect. 4.7), any function in the space weighted by the invariant measure can be written as an infinite linear combination of the eigenfunctions of the generator of the diffusion process.

Remark 2.4

In the construction of the martingale estimating function , we omitted the first index because, for ergodic diffusion processes, the first eigenvalue is zero, , and its corresponding eigenfunction is constant, , and hence, they would give independently of the function . Therefore, it would not provide us with any information about the unknown parameters in the drift.

The estimator . The first estimator we propose for the homogenized drift coefficient A is given by the solution of the M-dimensional nonlinear systemAn intuition on why is a good score function is given by the following result. Let be the score function where the observations of the slow variable of the multiscale process are replaced by the homogenized ones, then due to equation (2.12)which means that the zero of the expectation of the score function with homogenized observations is exactly the drift coefficient of the effective equation. In Algorithm 1, we summarize the main steps for computing the estimator and further details about the implementation can be found in Appendix B. We finally introduce the following technical assumption which will be employed in the analysis.

Assumption 2.5

The following hold for all and for all : where the dot denotes either the Jacobian matrix or the gradient with respect to a.

is continuously differentiable with respect to a for all ;

all components of , , , are polynomially bounded;

the slow-scale potential V is such that , , , and all components of , , are polynomially bounded;

Remark 2.6

In Kessler and Sørensen (1999) the authors propose a method to choose the functions in order to obtain optimality in the sense of Godambe and Heyde (1987): This optimal set of functions can be seen as the projection of the score function onto the set of martingale estimating functions obtained by varying the function . For the class of diffusion processes for which the eigenfunctions are polynomials, the optimal estimating functions can be computed analytically. In fact, they are related to the moments of the transition density, which can be computed explicitly. Moreover, another procedure is to choose functions which depend only on the unknown parameter and which minimize the asymptotic variance. This approach is strongly related to the asymptotic optimality criterion considered by Heyde and Gay (1989). For further details on how to choose these functions we refer to Kessler and Sørensen (1999), and we remark that their calculation requires additional computational cost. Nevertheless, the theory we develop is valid for all functions which satisfy Assumptions 2.5(i) and 2.5(ii) and we observed in practice that choosing simple functions independent of the unknown parameter, e.g. monomials of the form with , is sufficient to obtain satisfactory estimations. We also remark that in one dimension we can characterize completely all diffusion processes whose generator has orthogonal polynomials as eigenfunctions (Bakry et al. 2014, Sect. 2.7). Partial results in this directions also exist in higher dimensions.

The filtering approach

We now go back to our multiscale SDE (2.1) and, inspired by Abdulle et al. (2021), we propose a second estimator for the homogenized drift coefficient by filtering the data. In particular, we modify by filtering the observations and inserting the new data into the score function in order to take into account the case when the step size is too small with respect to the multiscale parameter . Let us consider the exponential kernel defined asfor which a rigorous theory has been developed in Abdulle et al. (2021). We remark that this exponential kernel is a low-pass filter, which cuts the high frequencies and highlights the slowest components. We then define the filtered observations choosing and computing the weighted average for all where the fast-scale component of the original multiscale trajectory is eliminated, and we define the new score function as a modification of (2.14), i.e.

Remark 2.7

Notice that the filtered data only partially replace the original data in the definition of the score function. This idea is inspired by Abdulle et al. (2021) where the same approach is used with the maximum likelihood estimator. The importance of keeping also the original observations becomes apparent in the proofs of the main results. However, a simple intuition is provided by equation (2.12). This equation is essential in order to obtain the unbiasedness of the estimators when the sampling rate is independent of the multiscale parameter , but it is not valid for the filtered process.

The estimator . The second estimator is given by the solution of the M-dimensional nonlinear systemThe main steps to compute the estimator are highlighted in Algorithm 2 and additional details about the implementation can be found in Appendix B. Note that (2.16) can be rewritten asWe introduce its continuous version which will be employed in the analysisWe remark that the joint process satisfies the system of multiscale SDEsand, using the theory of homogenization, when goes to zero it converges in law as a random variable in to the two-dimensional process , which solvesMoreover, it has been proved in Abdulle et al. (2021) that the two-dimensional processes and are geometrically ergodic and their respective invariant measures have densities with respect to the Lebesgue measure denoted respectively by and . Let us finally remark that given discrete observations we can only compute , but the theory, which has to be employed for proving the convergence results, has been studied for the continuous-time process .

Remark 2.8

The only difference in the construction of the estimators and is the fact that the latter requires filtered data, which are obtained from discrete observations, and thus, it is computationally more expensive. Therefore, when it is possible to use the estimator without filtered data, it is preferable to employ it.

Main results

In this section, we present the main results of this work, i.e. the asymptotic unbiasedness of the proposed estimators. We first need to introduce the following technical assumption, which is a nondegeneracy hypothesis related to the use of the implicit function theorem for the functions (2.14) and (2.17) in the limit as .

Assumption 3.1

Let A be the homogenized drift coefficient of equation (2.2). Then, the following hold where is the invariant measure of the couple , whose existence is guaranteed by Lemma A.2, and is the gradient of the stochastic process in (2.7) with respect to the drift coefficient a.

,

,

,

,

Remark 3.2

The nondegeneracy Assumption 3.1, which is analogous to Condition 4.2(a) in Kessler and Sørensen (1999), holds true in all nonpathological examples and does not constitute an essential limitation on the range of validity of the results proved in this paper. Further details about the necessity of this assumption for the analysis of the proposed estimator will be given in Sect. 5.2.

The proofs of the following two main theorems are the focus of Sect. 5.

Theorem 3.3

Let J be a positive integer. Under Assumptions 2.2, 2.5, 3.1 and if is independent of or with , there exists such that for all , an estimator which solves the system exists with probability tending to one as . Moreover,where A is the homogenized drift coefficient of equation (2.2).

Theorem 3.4

Let J be a positive integer. Under Assumptions 2.2, 2.5, 3.1 and if is independent of or with and , , there exists such that for all an estimator which solves the system exists with probability tending to one as . Moreover,where A is the homogenized drift coefficient of equation (2.2).

Remark 3.5

Notice that in both Theorem 3.3 and Theorem 3.4 the order of the limits is important and they cannot be interchanged. In fact, we first consider the large data limit, i.e. the number of observations N tends to infinity, and then we let the multiscale parameter vanish. Moreover, in Theorem 3.4 the values and are not allowed because of technicalities in the proof, but we observe numerically that the estimator works well also in these two particular cases.

These two theorems show that both estimators based on the multiscale data from (2.1) converge to the homogenized drift coefficient A of (2.2). Since the analysis is similar for the two cases, we will mainly focus on the second score function with filtered observations and at the end of each step we will state the differences with respect to the estimator without pre-processed data.

Remark 3.6

Since the main goal of this work is the estimation of the effective drift coefficient A, in the numerical experiments and in the following analysis we will always assume the effective diffusion coefficient to be known. Nevertheless, we remark that our methodology can be slightly modified in order to take into account the estimation of the effective diffusion coefficient too. In fact, the parameter a can be replaced by the parameter where a stands for the drift and s stands for the diffusion, yielding nonlinear systems of dimension corresponding to (2.15) and (2.18). The proofs of the asymptotic unbiasedness of the new estimators and can be adjusted analogously. For completeness, we provide a more detailed explanation and a numerical experiment illustrating this approach in Sect. 4.6.

A particular case

Before analysing the general framework, let us consider the simple case of the Ornstein–Uhlenbeck process, i.e. let the dimension of the parameter and let . Then, the multiscale SDE (2.1) becomesand its homogenized version isLetting , then the eigenfunctions and the eigenvalues satisfyThe solution of the eigenvalue problem can be computed explicitly (see Pavliotis 2014, Sect. 4.4); we haveand satisfies the recurrence relationwith and . It is also possible to prove by induction thatLet us consider the simplest case with only one eigenfunction, i.e. , and , which impliesThen, the score functions (2.14) and (2.17) becomeThe solutions of the equations and can be computed analytically and are given byandComparing these estimators with the discrete MLE defined in Pavliotis and Stuart (2007) without filtered data asand the discrete MLE with filtered datawe notice that they coincide in the limit as vanishes. We remark that we are comparing our estimator with the discrete MLE instead of the analytical formula for the MLE in continuous time since we assume that we are observing our process at discrete times. Therefore, the continuous time MLE has to be approximated using the available discrete data (Pavliotis 2014, Sect. 5.3). In the following theorems we show the asymptotic limit of the estimators. We do not provide a proof for these results since Theorem 3.7 and Theorem 3.9 are particular cases of Theorem 3.3 and Theorem 3.4 respectively, and Theorem 3.8 follows from the proof of Theorem 3.3 as highlighted in Remark 5.11.

Theorem 3.7

Let be independent of or with . Then, under Assumption 2.2, the estimator (3.1) satisfieswhere A is the drift coefficient of the homogenized equation (2.2).

Theorem 3.8

Let be independent of or with . Then, under Assumption 2.2, the estimator (3.1) satisfieswhere is the drift coefficient of the homogenized equation (2.1).

Theorem 3.9

Let be independent of or with , . Then, under Assumption 2.2, the estimator (3.2) satisfieswhere A is the drift coefficient of the homogenized equation (2.2).

Remark 3.10

Notice that it is possible to write different proofs for Theorems 3.7, 3.8 and 3.9, which take into account the specific form of the estimators, and thus show stronger results. In fact, if the distance between two consecutive observations is independent of the multiscale parameter , then the convergences in the statements do not only hold in probability, but also almost surely. We expect that almost sure convergence can be proved for a larger class of equations, but are neither aware of related literature showing such a stronger result, nor have been able to prove it.

Numerical experiments

In this section, we present numerical experiments which confirm our theoretical results and show the power of the martingale estimating functions based on eigenfunctions and filtered data to correct the unbiasedness caused by discretization and the fact that we are using multiscale data to fit homogenized models. Moreover, we present a sensitivity analysis with respect to the number N of observations and the number J of eigenvalues and eigenfunctions taken into account. In the experiments that we present data are generated employing the Euler–Maruyama method with a fine time step h, in particular we set . Letting , we generate data for and we select a sequence of observations , where and with . In view of Remark 2.3, we do not require stationarity of the multiscale dynamics; hence, we always set the initial condition to be . Notice that the time step h is only used to generate numerically the original data and has to be chosen sufficiently small in order to have a reliable approximation of the continuous path. However, the distance between two consecutive observations is the rate at which we sample the data, which we assume to know, from the original trajectory. In order to compute the filtered data , we employ equation (2.19). We repeat this procedure for different realizations of Brownian motion and we plot the average of the drift coefficients computed by the estimators. We finally remark that in order to compute our estimators we need the value of the diffusion coefficient of the homogenized equation. In all the numerical experiments, we compute it exactly using the formula for the coefficient K given by the theory of homogenization, but we also remark that its value could be estimated employing the subsampling technique presented in Pavliotis and Stuart (2007) or modifying the estimating function as explained in Remark 3.6.

Sensitivity analysis with respect to the number of observations

We consider the multiscale Ornstein–Uhlenbeck process, i.e. equation (2.1) with , and we take , the multiscale parameter , the drift coefficient and the diffusion coefficient . Notice that for this choice of the slow-scale potential the technical assumptions required in the main Theorems 3.3, 3.4 can be easily checked. We plot the results computed by the estimator with and and we then divide the analysis in two cases: “small” and “big”.

Let us first consider “small”, i.e. with , and take . In Fig. 1, we plot the results of the estimator as a function of the number of observations N. We remark that in this case the number of observations needed to reach convergence is strongly dependent and inversely proportional to the distance between two consecutive observations. This means that in order to reach convergence we need the final time T to be sufficiently large independently of . In fact, when the distance is small, the discrete observations are a good approximation of the continuous trajectory and therefore what matters most is the length T of the original path rather than the number N of observations.

Fig. 1

Sensitivity analysis with respect to the number N of observations for different values of , for the estimator with

In order to study the case “big”, i.e. , we set with , and take . Figure 2 shows that in this case the number of observations needed to reach convergence is an increasing function of . Therefore, in order to have a reliable approximation of the drift coefficient of the homogenized equation, the final time T has to be chosen depending on . This is justified by the fact that, differently from the previous case, the discrete data are less correlated and therefore they do not well approximate the continuous trajectory. In particular, when the distance between two consecutive observations is very large, then in practice we need a huge amount of data because a good approximation of the unknown coefficient is obtained only if the final time T is very large.

Fig. 2

Sensitivity analysis with respect to the number N of observations for different values of , for the estimator with

Sensitivity analysis with respect to the number of eigenvalues and eigenfunctions

Let us now consider equation (2.1) with four different slow-scale potentialsThe other functions and parameters of the SDE are chosen as in the previous subsection, i.e. , , and . Moreover, we set and and we vary . The functions appearing in the estimating function are given by for all .

In Fig. 3, where we plot the values computed by and , we observe that the number J of eigenvalues and eigenfunctions slightly improve the results, in particular for the fourth potential, but the estimation stabilizes when the number of eigenvalues J is still small, e.g. . Therefore, in order to reduce the computational cost, it seems to be preferable not to take large values of J. This is related to how quickly the eigenvalues grow and, therefore, how quickly the corresponding exponential terms decay. The rigorous study of the accuracy of the spectral estimators as a function of the number of eigenvalues and eigenfunctions that we take into account will be investigated elsewhere.

Fig. 3

Sensitivity analysis with respect to the number J of eigenvalues and eigenfunctions for different slow-scale potentials, for the estimators and

Verification of the theoretical results

We consider the same setting as in the previous subsection, i.e. equation (2.1) with slow-scale potentials given by (4.1) and , , and . Moreover, we set , and , and we choose the distance between two successive observations to be with .

In Fig. 4, we compare our martingale estimator without filtered data with the discrete maximum likelihood estimator denoted . The MLE does not provide good results for two reasons:Nevertheless, as observed in these numerical experiments and investigated in greater detail in Pavliotis and Stuart (2007), there exists an optimal value of such that works well, but this value is not known a priori and is strongly dependent on the problem, hence this technique is not robust. Figure 4 shows that the second issue, i.e. when is relatively big, can be solved employing , an estimator for discrete observations, and that filtering the data is not needed as proved in Theorem 3.3.

Fig. 4

Comparison between the discrete maximum likelihood estimator presented in Pavliotis and Stuart (2007) and our estimator with without filtered data as a function of the distance between two successive observations for different slow-scale potentials

if is small, more precisely if with , sampling the data does not completely eliminate the fast-scale components of the original trajectory; therefore, since we are employing data generated by the multiscale model, the estimator is trying to approximate the drift coefficient of the multiscale equation, rather than the one of the homogenized equation;

if is relatively big, in particular if with , then we are taking into account only the slow-scale components of the original trajectory, but a bias is still introduced because we are discretizing an estimator which is usually used for continuous data.

Then, in order to solve also the first problem, in Fig. 5 we compare with our martingale estimator with filtered data. We observe that inserting filtered data in the estimator allows us to disregard the fast-scale components of the original trajectory and to obtain good approximations of the drift coefficient A of the homogenized equation independently of , as already shown in Theorem 3.4. In particular, we notice that the results still improve even for big values of if we employ the estimator based on filtered data. Finally, as highlighted in Remark 5.11, we observe that the limiting value of the estimator as the number of observations N goes to infinity and the multiscale parameter vanishes is strongly dependent on the problem and cannot be computed theoretically. However, if we consider the slow-scale potential , i.e. the multiscale Ornstein-Uhlenbeck process, then the limit, as proved in Theorem 3.8, is the drift coefficient of the multiscale equation.

Fig. 5

Comparison between our two estimators without filtered data and with filtered data with as a function of the distance between two successive observations for different slow-scale potentials

Multidimensional drift coefficient

In this experiment, we consider a multidimensional drift coefficient, in particular we set . We then consider the bistable potential, i.e.and the fast-scale potential . We choose the exact drift coefficient of the multiscale equation (2.1) to be and the diffusion coefficient to be . We also set the number of eigenfunctions , the function , the distance between two consecutive observations and the final time . We then compute the estimator after observations and in Fig. 6 we plot the result of the experiment for the cases and . Since we are analysing the case independent of , filtering the data is not necessary and therefore we consider the estimator which is computationally less expensive to compute.

Fig. 6

Evolution in time of the estimator with for a two-dimensional drift coefficient

We observe that the estimation is approaching the exact value A of the drift coefficient of the homogenized equation as the number of observations increases, until it starts oscillating around the true value . Moreover, we notice that the time needed to reach a neighbourhood of A is smaller when the multiscale parameter is closer to its vanishing limit. In Table 1, we report the absolute error defined aswhere denotes the Euclidean norm, varying the number of observations N for the two values of the multiscale parameter.

Table 1

Absolute error defined in (4.2) between the homogenized drift coefficient A and the estimator with for a two-dimensional drift coefficient

N	100	200	300	400	500	600	700	800	900	1000
ε=0.1	0.742	0.395	0.215	0.201	0.093	0.036	0.011	0.027	0.034	0.028
ε=0.05	0.086	0.031	0.019	0.031	0.018	0.049	0.081	0.085	0.055	0.053

Multidimensional stochastic process: interacting particles

In this section, we consider a system of d interacting particles in a two-scale potential, a problem with a wide range of applications which has been studied in Gomes and Pavliotis (2018). For and for all , consider the system of SDEsIn this paper we fix the number of particles and study the performance of our estimators as vanishes. The very interesting problem of inference for mean field SDEs, obtained in the limit as , will be investigated elsewhere. It can be shown (see, for example, Gomes and Pavliotis 2018, Sect. 2.1 and Duncan and Pavliotis (2016); Delgadino et al. (2021)) that converges in law as goes to zero to the solution of the homogenized systemwhere and K is defined in (2.3). Moreover, the first eigenvalue and eigenfunction of the generator of the homogenized system can be computed explicitly and they are given, respectively, byHence, letting independent of , given a sequence of observations , we can express the estimators analyticallyLet us now set , , and . We then simulate system (4.3) for different final times and approximate the drift coefficient A of the homogenized system (4.4) for and . In Figs. 7 and 8, we plot the results, respectively, of the estimators with and with for two different values of . As expected, we observe that our estimator provides a better approximation of the unknown coefficient A when the time T increases and that this value stabilizes after approximately .

Fig. 7

Evolution in time of the estimator with for a d-dimensional system of interacting particles with sampling rate

Fig. 8

Evolution in time of the estimator with for a d-dimensional system of interacting particles with sampling rate

Simultaneous inference of drift and diffusion coefficients

As highlighted by Remark 3.6, a small modification of our methodology allows us to estimate the diffusion coefficient, in addition to drift coefficients. Define the parameter , whose exact value is given by , where A and are the drift and diffusion coefficients of the homogenized equation, respectively. Then, the eigenvalue problem reads for all where the eigenvalues and eigenfunctions are now dependent on the new parameter . Accordingly, also the functions can be chosen dependent on both the drift and diffusion coefficients and, moreover, they have to take values in , i.e. . Therefore, the new score functions and are defined from , which is the set of admissible parameters , to and thus give nonlinear systems of dimension . Finally, the solutions and of the systems are the estimators of both the drift and diffusion coefficients of the homogenized equation. In fact, small modifications in the proofs of the main results, in particular in the notation, yield the asymptotic unbiasedness of the estimators under the same conditions, i.e.Consider now the same setting of Sect. 4.1, i.e. the multiscale Ornstein–Uhlenbeck potential with , , , and let us assume that both the drift and diffusion coefficients are unknown. We remark that in this case we have . Then, set the final time , the sampling rate and the number of eigenfunctions and eigenvalues . Moreover, we choose the functions . Since the distance between two consecutive observations is independent of the multiscale parameter , we consider the estimator without filtered data. In Fig. 9, we plot the evolution of our estimator varying the number of observations N for two different values of , in particular and . We observe that if the multiscale parameter is smaller, then the number of observations needed to obtain a reliable approximation of the unknown parameters is lower.

Fig. 9

Simultaneous inference of drift and diffusion coefficient for the estimator with

Asymptotic unbiasedness

In this section, we prove our main results. The plan of the proof is the following:We first define the Jacobian matrix of the function introduced in (2.13) with respect to a: which will be employed in the following and where denotes the outer product in and the dot denotes either the Jacobian matrix or the gradient with respect to a, e.g. . Then note that, under Assumption 2.2, due to ergodicity and stationarity and by Bibby and Rensen (1995, Lemma 3.1) we haveand where and denotes, respectively, that and are distributed according to their invariant distribution. We remark that the invariant distribution exists due to Lemma A.2. By equation (5.1) the Jacobian matrices of and with respect to a are given byand

we first study the limiting behaviour of the score functions and defined in (2.14) and (2.17) as the number of observations N goes to infinity, i.e. as the final time T tends to infinity;

we then show the continuity of the limit of the score functions obtained in the previous step and we compute their limits as the multiscale parameter vanishes (Sect. 5.1);

we finally prove our main results, i.e. the asymptotic unbiasedness of the drift estimators (Sect. 5.2).

Continuity of the limit of the score function

In this section, we first prove the continuity of the functions and . We then study the limit of these functions for . As the proof for the filtered and the non-filtered are similar, we will concentrate on the filtered case and comment on the non-filtered case. Before entering into the proof, we give two preliminary technical lemmas which will be used repeatedly and whose proof can be found, respectively, in Appendix A.1 and Appendix A.3.

Lemma 5.1

Let be defined in (2.16) and distributed according to the invariant measure of the process . Then, for any there exists a constant uniform in such that

Lemma 5.2

Let be a function which is polynomially bounded along with all its derivatives. Then,where satisfies for all and for a constant independent of and

We start here with a continuity result for the score function and its Jacobian matrix with respect to the unknown parameter.

Proposition 5.3

Under Assumption 2.5, the functions and defined in (5.2) and (5.3), where can be either independent of or with , are continuous.

Proof

We only prove the statement for , then the argument is similar for . Letting and , we want to show thatBy the triangle inequality, we havethen we divide the proof in two steps and we show that the two terms vanish.

Step 1: as .

Since and are continuously differentiable with respect to a for all , respectively, due to Assumption 2.5 and Lemma A.4, then also is continuously differentiable with respect to a. Therefore, by the mean value theorem for vector-valued functions, we haveThen, letting be a constant independent of , since and are polynomially bounded still by Assumption 2.5 and , and have bounded moments of any order by Pavliotis and Stuart (2007, Corollary 5.4) and Lemma 5.1, we obtainwhich implies that vanishes as goes to both if is independent of and if .

Step 2: as .

If is independent of , then we haveand the right-hand side vanishes due to the continuity of for all and the continuity of the solution of a stochastic differential equation with respect to a parameter (see Krylov 2009, Theorem 2.8.1). Let us now consider the case with and let us assume, without loss of generality, that . Denoting and applying Itô’s lemma we have for all then we can writewhere is given byLet be independent of and notice that since is bounded, are polynomially bounded and and have bounded moments of any order by Pavliotis and Stuart (2007, Corollary 5.4) and Lemma 5.1, applying Hölder’s inequality we obtain Therefore, by the continuity of the solution of a stochastic differential equation with respect to a parameter (see Mishura et al. 2010) and due to the bound (5.4), we deduce thatwhich implies that vanishes as goes to and concludes the proof.

Remark 5.4

Notice that the proof of Proposition 5.3 can be repeated analogously for the functions and without filtered data in order to prove their continuity.

Next we study the limit as vanishes and we divide the analysis in two cases. In particular, we consider independent of and with . In the first case (Proposition 5.5), data are sampled at the homogenized regime ignoring the fact that the they are generated by a multiscale model, while in the second case (Proposition 5.7) the distance between two consecutive observations is proportional to the multiscale parameter and thus, data are sampled at the multiscale regime preserving the multiscale structure of the full path.

Proposition 5.5

Let the functions and be defined in (5.2) and (5.3) and let be independent of . Under Assumption 2.5 and for any , we have

We only prove the statement for , then the argument is similar for . By the triangle inequality, we havewherewhich vanishes due to the first step of the proof of Proposition 5.3 andLet us remark that the convergence in law of the joint process to the joint process by Lemma A.2 implies the convergence in law of the triple to the triple since , and , . Therefore, we havewhich implies the desired result.

Remark 5.6

Similar results to Proposition 5.3 and Proposition 5.5 can be shown for the estimator without filtered data. In particular we have that and are continuous in andSince the proof is analogous, we do not report here the details.

Proposition 5.7

Let the functions and be defined in (5.2) and (5.3) and let with and , . Under Assumption 2.5 and for any , we have where the generator is defined in (2.9).

, where

, where

We only prove the statement for , then the argument is similar for . By the triangle inequality, we havethen we need to show that the two terms vanish and we distinguish two cases.

Case 1: .

Applying Lemma 5.2 to the functions for all and noting thatsinceis a martingale with , we have where satisfies for a constant independent of and and for all We now study the three terms separately. First, by Cauchy–Schwarz inequality, since is polynomially bounded, has bounded moments of any order by Lemma 5.1 and due to (5.5) we obtain Let us now focus on for which we havewhere is distributed according to the invariant measure of the continuous process andBy the mean value theorem for vector-valued functions, we haveand since are polynomially bounded, , , have bounded moments of any order, respectively, by Pavliotis and Stuart (2007, Corollary 5.4), Abdulle et al. (2021, Lemma C.1) and Lemma 5.1 and applying Hölder’s inequality and Corollary A.3 we obtainMoreover, notice that by homogenization theory (see Abdulle et al. 2021, Sect. 3.2) the joint process converges in law to the joint process and thereforewhich together with (5.7) and (5.8) yieldsWe now consider and we follow an argument similar to . We first havewhere the first term in the right-hand side converges due to homogenization theory and the second one is bounded byTherefore, we obtainwhich together with (5.6) and (5.9) implieswhich shows that vanishes as goes to zero. Let us now consider . Following the first step of the proof of Proposition 5.3, we havewhere assumes values in the line connecting a and , and repeating the same computation as above we obtainwhich together with (5.10) gives the desired result.

Case 2: .

Let be distributed according to the invariant measure of the continuous process and defineThen, we haveand we first bound the remainder . Applying Itô’s lemma to the process with the functions for each we haveand observing thatsinceis a martingale with , we obtain By the mean value theorem for vector-valued functions, we haveand since are polynomially bounded, , , have bounded moments of any order, respectively, by Pavliotis and Stuart (2007, Corollary 5.4), Abdulle (2021, Lemma C.1) and Lemma 5.1 and applying Hölder’s inequality, we obtainfor a constant independent of and . We repeat a similar argument for and to getwhich together with (5.14) yieldMoreover, applying Lemma 5.2 and proceeding similarly to the first part of the first case of the proof, we havewhich together with (5.15) and Corollary A.3 impliesLet us now consider . Replacing equation (5.12) into the definition of in (5.11) and observing that similarly to (5.13), it holdswe obtainWe rewrite inside the integrals employing equation (2.21) and Itô’s lemmahence due to stationarity we havewhereandSince and are polynomially bounded, is bounded and and have bounded moments of any order, respectively, by Pavliotis and Stuart (2007, Corollary 5.4) and Abdulle et al. (2021, Lemma C.1), is bounded byLet us now move to and let us define the functionswhere and are, respectively, the densities with respect to the Lebesgue measure of the invariant distributions of the joint processes and and and are their marginals with respect to the first component. Integrating by parts we havewhich impliesEmploying the last equation in the proof of Lemma 3.5 in Abdulle et al. (2021) with and , we haveand thus we obtainLetting go to zero and due to homogenization theory, it followsthen applying formula (5.19) for the homogenized equation, i.e. with and and replaced by A and , and integrating by parts we have Therefore, we obtainwhich together with (5.11), (5.17) and bounds (5.16) and (5.18) implies that vanishes as goes to zero. Finally, analogously to the first case we can show that also vanishes, concluding the proof.

Remark 5.8

A similar result to Proposition 5.7 can be shown for the estimator without filtered data only if , i.e. the first case in the proof. In particular, we have where the generator is defined in (2.9). Since the proof is analogous, we do not report here the details. On the other hand, if , we can show that The proof is omitted since it is similar to the second case of the proof of Proposition 5.7.

, where

, where

, where

, where

Proof of the main results

Let us remark that we aim to prove the asymptotic unbiasedness of the proposed estimators, i.e. their convergence to the homogenized drift coefficient A as the number of observations N tends to infinity and the multiscale parameter vanishes. Therefore, we study the limit of the score functions and their Jacobian matrices as and evaluated in the desired limit point A.

We first analyse the case independent of and we consider the limit of Proposition 5.5 and Remark 5.6 evaluated in . Then, due to equation (2.12) we getand similarly we obtainOn the other hand, if is a power of , we study the limit of Proposition 5.7 and Remark 5.8 evaluated in and by (2.10) we haveMoreover, differentiating equation (2.12) with respect to a, we getwhere the process satisfiesTherefore, due to (2.12) and (5.23), we haveandThen, due to Lemma A.4, we can differentiate the eigenvalue problem (2.11) with respect to a and deduce thatwhere the dot denotes the gradient with respect to a and which together with (2.11) impliesandBefore showing the main results, we need two auxiliary lemmas, which in turn rely on the technical Assumption 3.1, which can now be rewritten as: Since the proofs of the two lemmas are similar, we only write the details of the first one.

,

,

,

.

Lemma 5.9

Under Assumption 2.5 and Assumption 3.1, there exists such that for all there exists such that if is independent of or with and , Moreover

Let us first extend the functions and by continuity in with their limit given by Proposition 5.5 and Proposition 5.7 depending on and note that due to (5.21) if is independent of and (5.22) otherwise, we haveMoreover, by (5.24), (5.25) and Assumption 3.1, we know thatTherefore, since the functions and are continuous by Proposition 5.3, the implicit function theorem (see Hurwicz and Richter 2003, Theorem 2) gives the desired result.

Lemma 5.10

Under Assumption 2.5 and Assumption 3.1, there exists such that for all there exists such that if is independent of or with Moreover,

We are now ready to prove the asymptotic unbiasedness of the estimators, i.e. Theorems 3.3 and 3.4. We only prove Theorem 3.4 for the estimator with filtered data. The proof of Theorem 3.3 for the estimator without filtered data is analogous and is omitted here.

Proof of Theorem 3.4

We need to show for a fixed : We first note that by Lemma 5.9 we haveWe then follow the steps of the proof of Bibby and Rensen (1995, Theorem 3.2). Due to Barndorff–Nielsen and Sorensen (1994, Theorem A.1), claims (i) and (ii) hold true if we verify thatand as where is a positive definite covariance matrix andfor small enough such that . Result (5.27) is a consequence of Florens–Zmirou (1989, Theorem 1). We then have where the right-hand side vanishes by Bibby and Rensen (1995, Lemma 3.3) and the continuity of (Proposition 5.3), implying result (5.26). Hence, we proved (i) and (ii), which conclude the proof of the theorem.

the existence of the solution of the system with probability tending to one as ;

in probability with .

Remark 5.11

Notice that if with and we do not employ the filter, in view of (5.20) and following the same proof of Theorem 3.4, we could compute the asymptotic limit of as N goes to infinity and vanishes if we knew such thatThe value of cannot be found analytically since it is, in general, different from the drift coefficients and A of the multiscale and homogenized equations (2.1) and (2.2). Nevertheless, we observe that in the simple scale of the multiscale Ornstein–Uhlenbeck process we have .

Conclusion

In this work, we presented new estimators for learning the effective drift coefficient of the homogenized Langevin dynamics when we are given discrete observations from the original multiscale diffusion process. Our approach relies on a martingale estimating function based on the eigenvalues and eigenfunctions of the generator of the coarse-grained model and on a linear time-invariant filter from the exponential family, which is employed to smooth the original data. We studied theoretically the convergence properties of our estimators when the sample size goes to infinity and the multiscale parameter describing the fastest scale vanishes. In Theorem 3.3 and Theorem 3.4, we proved, respectively, the asymptotic unbiasedness of the estimators with and without filtered data. We remark that the former is not robust with respect to the sampling rate at finite multiscale parameter, while the estimator with filtered data is robust independently of the sampling rate. We analysed numerically the dependence of our estimators on the number of observations and the number of eigenfunctions employed in the estimating function noticing that the first eigenvalues in magnitude are sufficient to approximate the drift coefficient. Moreover, we performed several numerical experiments, which highlighted the effectiveness of our approach and confirmed our theoretical results. We believe that eigenfunction estimators can be very useful in applications, for example to multiparticle systems and their mean field limit (Gomes and Pavliotis 2018), since the eigenvalue problem for the generator of a reversible Markov process is a very well-studied problem. This means, in particular, that it is possible to study rigorously the proposed estimators and to prove asymptotic unbiasedness and asymptotic normality. Furthermore, in order to be able to assess the accuracy of the estimators, we could analyse its rate of convergence with respect to both the number of observations and the fastest scale. This is a highly nontrivial problem since it first requires the development of a fully quantitative periodic homogenization theory and we will return to this problem in future work. Finally, we think that it would be interesting to extend our estimators to the nonparametric framework and consider more general multiscale models.