# Parametric inference for multidimensional hypoelliptic diffusion with full observations

Multidimensional hypoelliptic diffusions arise naturally as models of neuronal activity. Estimation in those models is complex because of the degenerate structure of the diffusion coefficient. We build a consistent estimator of the drift and variance parameters with the help of a discretized log-likelihood of the continuous process in the case of fully observed data. We discuss the difficulties generated by the hypoellipticity and provide a proof of the consistency of the estimator. We test our approach numerically on the hypoelliptic FitzHugh-Nagumo model, which describes the firing mechanism of a neuron.

02/24/2020

### Adaptive and non-adaptive estimation for degenerate diffusion processes

We discuss parametric estimation of a degenerate diffusion system from t...
05/28/2021

### Parameter estimation in CKLS model by continuous observations

We consider a stochastic differential equation of the form dr_t = (a - b...
04/29/2020

### Adaptive tests for parameter changes in ergodic diffusion processes from discrete observations

We consider the adaptive test for the parameter change in discretely obs...
07/24/2018

### Contrast function estimation for the drift parameter of ergodic jump diffusion process

In this paper we consider an ergodic diffusion process with jumps whose ...
04/17/2019

### Nonparametric drift estimation for diffusions with jumps driven by a Hawkes process

We consider a 1-dimensional diffusion process X with jumps. The particul...
09/06/2021

### Unbiased Estimation of the Hessian for Partially Observed Diffusions

In this article we consider the development of unbiased estimators of th...
02/12/2018

### Estimating Diffusion With Compound Poisson Jumps Based On Self-normalized Residuals

This paper considers parametric estimation problem of the continuous par...

## 1 Introduction

Hypoelliptic diffusions naturally occur in various applications, most notably in neuroscience, molecular physics and in mathematical finance. In particular, multidimensional models of a neuron population, such as stochastic approximation of a Hawkes process (Ditlevsen and Löcherbach, 2017), or exotic models of option pricing (Malliavin and Thalmaier, 2006) are described by a hypoelliptic diffusion.

The main difference to the classical (or elliptic) setting is that in the hypoelliptic case the dimensionality of the noise is lower than the dimensionality of the system of stochastic differential equations (SDE), which describes the process. Hypoellipticity can be intuitively explained in the following way: though the covariance matrix of noise is singular due to a degenerate diffusion coefficient, smooth transition density with respect to the Lebesgue measure still exists. That is the case when the noise is propagated to all the coordinates through the drift term.

Properties of hypoelliptic diffusions significantly differ from those of elliptic ones, when all coordinates are driven by a Gaussian noise. Thus they are more difficult to study. The first problem is that each coordinate has a variance of different order. It is the main cause why classical numerical approximation methods do not work well with hypoelliptic diffusions. In particular, it is proven that for hypoelliptic systems the classical Euler-Maruyama scheme does not preserve ergodic properties of the true process (Mattingly et al., 2002). The second problem is the degenerate diffusion coefficient. As the explicit form of the transition density is often unknown, parametric inference is usually based on its discrete approximation with the piece-wise Gaussian processes (see, for example Kessler (1997)). But in the hypoelliptic case this approach cannot be applied directly because the covariance matrix of the approximated transition density is not invertible.

Now let us be more specific. Consider a two-dimensional system of stochastic differential equations of the form:

 {dXt=a1(Xt,Yt;θ)dtdYt=a2(Xt,Yt;θ)dt+b(Xt,Yt;σ)dWt, (1)

where , is the drift term, is the diffusion coefficient,

is a standard Brownian motion defined on some probability space

,

is the vector of the unknown parameters, taken from some compact set

.

The goal of this paper is to estimate the parameters of (1) from discrete observations of both coordinates and . It is achieved in two steps: first, we consider a discretization scheme in order to approximate the transition density of the continuous process preserving the ergodic property, and then we propose an estimation technique which maximizes the likelihood function of the discrete approximate model. Let us discuss the solutions proposed by other authors for hypoelliptic systems of different types.

Several works treat the parametric inference problem for a particular case of system (1). It is natural to introduce first the class of stochastic Damping Hamiltonian systems, also known as Langevin equations (Gardiner and Collett, 1985). These hypoelliptic models arise as a stochastic expansion of 2-dimensional deterministic dynamical systems — for example, the Van der Pol oscillator (Van der Pol, 1920) perturbed by noise. They are defined as the solution of the following SDE:

 {dXt=YtdtdYt=a2(Xt,Yt;θ)dt+b(Xt,Yt;σ)dWt. (2)

The particular case of Hamiltonian systems with and is considered in Ozaki (1989), where the link between the continuous-time solution of (2) and the corresponding discrete model is obtained with the so-called local linearization scheme. The idea of this scheme is the following: for a system of SDE with a non-constant drift and a constant variance, its solution can be interval-wise approximated by a system with a linear drift, and the original covariance matrix being expanded by adding higher-order terms. It allows to construct a quasi Maximum Likelihood Estimator. Pokern et al. (2007) attempt to solve the problem of the non-invertibility of the covariance matrix for the particular case of system (2) with a constant variance with the help of Itô-Taylor expansion of the transition density. The parameters are then estimated with a Gibbs sampler based on the discretized model with the noise propagated into the first coordinate with order . This approach allows to estimate the variance coefficient, but it is not efficient for estimating the parameters of the drift term. In Samson and Thieullen (2012) it is shown that a consistent estimator for fully and partially observed data can be constructed using only the discrete approximation of the second equation of the system (2). This method works reasonably good in practice even for more general models when it is possible to convert a system (1) to a simpler form (2). However, the transformation of the observations sampled from the continuous model (1) requires the prior knowledge of the parameters involved in the first equation. The other particular case of (1), when and the drift term is linear and thus the transition density is known explicitly, is treated in Le-Breton and Musiela (1985). A consistent maximum likelihood estimator is then constructed in two steps — first, a covariance matrix of the process is estimated from available continuous-time observations, and then it is used for computing the parameters of the drift term. Few other works are also devoted to a non-parametric estimation of the drift and the variance terms (Cattiaux et al., 2014, 2016). To the best of our knowledge, for systems (1) the only reference is Ditlevsen and Samson (2017). They construct a consistent estimator using a discretization scheme based on a Itô-Taylor expansion. But the estimation is conducted separately for each coordinate, so it requires partial knowledge of the parameters of the system.

In this paper we propose a new estimation method, adjusting the local linearization scheme described in Ozaki (1989) developed for the models of type (2) to the more general class of SDEs (1). Under the hypoellipticity assumption this scheme propagates the noise to both coordinates of the system and allows to obtain an invertible covariance matrix. We start with describing the discretization scheme, approximating the transition density and proposing a contrast estimator based on the discretized log-likelihood. While we attempt to estimate the parameters included in the drift and diffusion coefficient simultaneuously, we also explain in which cases and how the estimation can be splitted. Then we study the convergence of the scheme and prove the consistency of the proposed estimator based on the 2-dimensional contrast. We also detail the different speeds of convergence of the diffusion parameters, which are implied by the hypoellipticity. To the best of our knowledge, the proof of this consistency is the first in the literature. We finish with some numerical experiments, testing the proposed approach on the hypoelliptic FitzHugh-Nagumo model and compare it to other estimators.

This paper is organized as follows: Section 2 presents the model and assumptions. Discrete model is introduced in Section 3. The estimators are studied in Section 4 and illustrated numerically in Section 5. We close with Section 6, devoted to conclusions and discussions. Formal proofs are gathered in Appendix.

## 2 Models and assumptions

### 2.1 Notations

We consider system (1). We assume that both variables are discretely observed at equally spaced periods of time on some finite time interval . The vector of observations at time is denoted by , where is the value of the process at the time . We further assume that it is possible to draw a sufficiently large and accurate sample of data, i.e that may be arbitrary large, and the partition size — arbitrary small. Let us also introduce the vector notations:

 dZt=A(Zt;θ)dt+B(Zt;σ)d~Wt,Z0=ω0,t∈[0,T] (3)

where , is a standard two-dimensional Brownian motion defined on the filtered probability space, is a -measurable dimensional random vector. Matrices and represent, respectively, the drift and the diffusion coefficient, that is and

 B(Zi;σ)=(000b(Zt;σ)). (4)

Throughout the paper we use the following abbreviations: and . We suppress the dependency on the parameter , when its value is clear from context, otherwise additional indices are introduced. True values of the parameters are denoted by . We also adopt the notations from Pokern et al. (2007) and refer to the variable which is directly driven by Gaussian noise as ”rough”, and to as ”smooth”.

### 2.2 Assumptions

Further, we are working under the following assumptions:

• Functions and have bounded partial derivatives of every order, uniformly in . Furthermore .

• Global Lipschitz and linear growth conditions. :

 ∥A(Zt;θ)−A(Zs;θ)∥+∥B(Zt;σ)−B(Zs;σ)∥≤K∥Zt−Zs∥ ∥A(Zt;θ)∥2+∥B(Zt;σ)∥2≤K2(1+∥Zt∥2),

where is the standard Euclidean norm. Further, denote by the initial value of the process , then .

• Process is ergodic and there exists a unique invariant probability measure

with finite moments of any order.

• Both functions and are identifiable, that is .

Assumption (A1) ensures that the weak Hörmander condition is satisfied, thus the system is hypoelliptic in the sense of stochastic calculus of variations (Nualart, 2006, Malliavin and Thalmaier, 2006). In order to prove it we first write the coefficients of system (3) as two vector fields, converting (3) from Itô to the Stratonovich form:

Then their Lie bracket is equal to

 [A0,A1]=(∂ya1∂xa2(x,y;θ)−12∂xb(x,y;σ)∂2xyb(x,y;σ)).

By (A1) the first element of this vector is not equal to , thus we conclude that and generate . That means that the weak Hörmander condition is satisfied and as a result the transition density for system (3) exists, though not necessarily has an explicit form.

(A2) implies the existence and uniqueness in law of the weak solution of the system (3) Karatzas and Shreve (1987) and (A4) is a standard condition which is needed to prove the consistency of the estimator.

(A3) ensures that we can apply the weak ergodic theorem, that is, for any continuous function with polynomial growth at infinity:

 1T∫T0f(Zs)ds⟶T→∞ν0(f)a.s.

We do not investigate the conditions under which process is ergodic, as it is not the main focus of this work. Ergodicity of the stochastic damping Hamiltonian system (2) is studied in Wu (2001). Conditions for a wider class of hypoelliptic SDEs can be found in Mattingly et al. (2002). It is also important to know that if the process is ergodic, then its sampling is also ergodic (Genon-Catalot et al., 2000).

## 3 Discrete model

### 3.1 Derivation

Now we introduce a discretization scheme which approximates the transition density of system (3). It is a generalized version of the local linearization scheme presented in Ozaki (1989) for systems (2). Recall that the process is observed at equally spaced periods of time of size . On each interval we consider a new process described by the homogeneous linear system with the constant diffusion coefficient. That is, we use the approximation on each interval of length , where is the Jacobian of the drift coefficient , computed at the beginning of such an interval. In other words, instead of working with (3), we study systems of the following type:

 dZs=JτZsds+B(Zτ;σ)d~Ws,Z0=Zτ,s∈(τ,τ+Δ], (5)

where the initial value is given by the observation of the true process at time . Note that the value of the diffusion matrix is fixed at time . Solution of (5) has an explicit form:

 Zs=ZτeJτs+∫sτeJτ(s−v)B(Zτ;σ)d~Wv,∀s∈(τ,τ+Δ].

Then the first moment and the covariance matrix of the process on each interval are given, respectively, by:

 E[Zs]=ZτeJτs, (6) Σ(Zs;θ,σ2)=E[(∫sτeJτ(s−v)B(Zτ;σ)d~Wv)(∫sτeJτ(s−v)B(Zτ;σ)d~Wv)T]. (7)

Continuous representation of the drift and the variance terms (6) - (7) is not convenient for the numerical implementation, so it has to be discretized. For the drift term the discretization is straightforward and follows directly from the definition of the matrix exponent. For the covariance matrix we use the following proposition, proof of which is postponed to appendix:

###### Proposition 1.

The second-order Taylor approximation of matrix defined in (7) has the following form:

 b2(Zτ;σ)⎛⎜⎝(∂ya1)2Δ33(∂ya1)Δ22+(∂ya1)(∂ya2)Δ33(∂ya1)Δ22+(∂ya1)(∂ya2)Δ33Δ+(∂ya2)Δ22+(∂ya2)2Δ33⎞⎟⎠+O(Δ4), (8)

where the derivatives are computed at time .

Note that the noise in the first coordinate appears only through . Thus, under assumption (A1) matrix (7) has rank 2, while the original covariance matrix is of rank 1. This fact allows us to use techniques developed for non-degenerate Gaussian diffusions. However, note that this matrix is still highly ill-conditioned, as its determinant is of order .

At this point we give up the continuous time setting and move to the discrete process. Let us denote the approximated form (8) of matrix (7) by . Then the approximation for the solution of (5) is given by:

 Zi+1=¯A(Zi;θ)+¯B(Zi;θ,σ)Ξi, (9)

where is a standard Gaussian 2-dimensional random vector, is any matrix such that , is an approximation of the conditional expectation of the drift (6), given by:

 ¯A(Zi;θ)=Zi+ΔA(Zi;θ)+Δ22JiA(Zi;θ)+⋯+Δkk!JkiA(Zi;θ)+O(Δk+1) (10)

Thus, on each small interval the discretized process can be approximated by the Gaussian process such that .

### 3.2 Properties of the discrete model

The difference between the true process and its approximation cannot be computed explicitly since the transition density of the process is in general unknown. But the moments of the process can be approximated by a moment generator function. For sufficiently smooth and integrable function

we know that:

 E(f(Zt+Δ)|Zt=z)=j∑i=0Δii!Lif(z)+O(Δj+1), (11)

where is the times iterated generator of model (3), given by

 Lf(z)=(∂zf(z))A(z)+122▽Bf(z),

where is a weighted Laplace type operator. Thus, the value of the process is approximated by:

 E(Zt+Δ|Zt=z)=z+ΔA(z;θ)+Δ22LA(z;θ)+O(Δ3), (12)

which coincides with (10) up to the terms of order . Adding more terms to (10) does not improve the accuracy of the scheme unless the model is linear. Further we assume that . Now let us denote by the first element of the vector , and by the second. We have the following proposition:

###### Proposition 2 (Weak convergence of the local linearization scheme).

For the following holds:

 E(Xi+1−¯A1(Zi;θ0)) =O(Δ3) E(Yi+1−¯A2(Zi;θ0)) =O(Δ3) E(Xi+1−¯A1(Zi;θ0))2 =(∂ya1)2θ0Δ33b2(Zi;σ0)+O(Δ4) E(Yi+1−¯A2(Zi;θ0))2 =Δb2(Zi;σ0)+O(Δ2) =(∂ya1)θ0Δ22b2(Zi;σ0)+O(Δ3)
###### Proof.

It follows directly from the comparison between the scheme (8) - (10) and the moment generator function (11). More detailed expansion of (11) can be found in Kloeden et al. (2003). ∎

For systems (3) it is suitable to approximate the drift term up to the order . It is not sufficient to use a first-order approximation, as in that case the variance of the first coordinate would be suppressed by the error in the drift term. However, for elliptic systems it is often enough to approximate a drift only up to .

## 4 Parameter estimation

Let us introduce a contrast function for system (3). In the elliptic case this function is defined as times the log-likelihood of the discretized model (Florens-Zmirou (1989), Kessler (1997)). In hypoelliptic case, however, we must modify this criterion taking into account the specific structure of the covariance matrix. Most notably, the contrast is obtained by dividing the first part by :

 L(θ,σ2;Z0:N)=12N−1∑i=0(Zi+1−¯A(Zi;θ))TΣ−1Δ(Zi;θ,σ2)(Zi+1−¯A(Zi;θ))+N−1∑i=0logdet(ΣΔ(Zi;θ,σ2)). (13)

Then the estimator is defined as:

 (^θ,^σ2)=argminθ,σ2L(θ,σ2;Z0:N) (14)

For system (2) this correction is not proposed nor in Ozaki (1989), nor in Pokern et al. (2007). We justify this bias theoretically in Lemma 1. However, the 2-dimensional criterion (14) is tricky to analyze because of the different orders of variance for the first and the second coordinate. As a result, under the scaling which allows to estimate the parameters included in the rough coordinate, we cannot say anything about the parameters from the smooth coordinate and vice versa. It does not introduce any additional restrictions on the numerical implementation, but must be taken into account during the theoretical study. When both equations are driven by the same parameters, the task is simplified as the parameters can be then estimated from a one-dimensional criteria, which involves only one of the two equations.

Let us denote the parameters in the first coordinate by , and in the second by . We start with considering the variance term. In order to prove the convergence of the estimator for the , we fix the value of the parameter in the first coordinate in order to focus on the parameters which directly regulate the diffusion term of the original process. It results in the following Lemma:

###### Lemma 1.

Under assumptions (A1)-(A4), and while the following holds:

 1NLN,ΔN(θ,σ2;Z0:N)Pθ⟶∫(b2(z;σ0)b2(z;σ)+logb2(z;σ))ν0(dz) (15)

Then we proceed with the drift term. Note that the first statement of the Lemma 2 suffices to obtain the consistency of the estimator, but only in the case when both equations are driven by the same parameters. However, in combination with the second term it allows us to establish the main result of the paper:

###### Lemma 2.

Under assumptions (A1)-(A4), and the following holds:

 limN→∞,ΔN→0ΔNN[LN,ΔN(θ,σ20;Z0:N)−LN,ΔN(θ0,σ20;Z0:N)]Pθ⟶6∫(a1(z;θ0)−a1(z;θ))2b2(z;σ20)(∂ya1)2θν0(dz)

Note that consistency for each term is obtained under different scalings. Scaling is standard for the diffusion coefficient, while — for the drift parameter (Kessler (1997)). But is very unusual. That means that each parameter converges to its real value with different speed. It is the main theoretical difficulty encountered while constructing a two-dimensional contrast, compared to 2 different one-dimensional contrasts proposed in Ditlevsen and Samson (2017). Finally, we establish the following theorem:

###### Theorem 1.

Under assumptions (A1)-(A4) and and the following holds:

 (^θ,^σ2)Pθ⟶(θ0,σ20)
###### Proof.

Let us start with the variance term. The result follows from Lemma 1. Denote the integral on the right side of (15) by . We can choose some subsequence such that converges to some . By the definition of the estimator we know that . But we also know that and thus . It proves the consistency of .

Consistency for the drift term essentially follows from Lemma 2 and from the identifiability of the drift functions. We can state that there exists a subsequence which converges to . At the same time, as is identifiable, thus . That means that the estimator is consistent with respect to the parameters of the first coordinate. Proof for the is analogous. ∎

Another way to unify the scaling and simplify the task of parameter estimation is to consider a quadratic variance of the drift term. This approach does not allow to estimate the parameters of the diffusion term, but is effective for the parameters of the drift. Its main advantage consists in the fact that the computation of the contrast function does not require finding an inverse of matrix (8), which is ill-conditioned and can cause numerical problems for small values of . Consider

 ^θQV=argminθ1N−1N−1∑i=0∥∥Zi+1−¯A(Zi;θ)∥∥2, (16)

with being defined in (10).

###### Proposition 3.

Under assumptions (A1)-(A4) the following holds:

 ^θQVPθ⟶θ0

What about the diffusion term, in the case when , parameter can be computed explicitly with the help of the sample covariance matrix. Good properties of this approach for the elliptic case are proven in Florens-Zmirou (1989). When it comes to the hypoelliptic systems, this approach must be modified, as the discretization step of order does not allow to compute the terms of order , which represent the propagated noise. However, the value of can still be inferred from the observations of the rough coordinate by computing

 ¯σ2=1NΔN−1∑i=0(Yi+1−Yi)2f(Xi,Yi) (17)

## 5 Simulation study

### 5.1 The model

The two estimators and are evaluated on the simulation study with a hypoelliptic stochastic neuronal model called FitzHugh-Nagumo model (Fitzhugh, 1961). It is a simplified version of the Hodgkin-Huxley model (Hodgkin and Huxley, 1952), which describes in a detailed manner activation and deactivation dynamics of a spiking neuron. First it was studied in the deterministic case, then it was improved by adding two sources of noise to the both coordinates, what results in an elliptic SDE. However, it is often argued that only ion channels are perturbed by noise, while the membrane potential depends on them in a deterministic way. This idea leads to a 2-dimensional hypoelliptic diffusion. In this paper we consider a hypoelliptic version with noise only in the second coordinate as studied in Leon and Samson (2017). More precisely, the behaviour of the neuron is defined through the solution of the system

 {dXt=1ε(Xt−X3t−Yt−s)dtdYt=(γXt−Yt+β)dt+σdWt, (18)

where the variable represents the membrane potential of the neuron at time , and is a recovery variable, which could represent channel kinetic. Parameter is the magnitude of the stimulus current and is often known in experiments, is a time scale parameter and is typically significantly smaller than , since moves ”faster” than . Parameters to be estimated are .

Hypoellipticity and ergodicity of (18) are proven in Leon and Samson (2017). In Jensen et al. (2012) it is proven that is unidentifiable when only one coordinate is observed. Parametric inference for elliptic FitzHugh-Nagumo model both in fully- and partially observed case is investigated in Jensen (2014). The same problem, but for the hypoelliptic setting is studied in Ditlevsen and Samson (2017).

### 5.2 Experimental design

In order to make our experiments more representative, we consider two different settings: an excitatory and an oscillatory behaviour. For the first regime, the drift parameters are set to and the diffusion coefficient , and for the second and . The diffusion coefficient does not change the behaviour pattern, only the ”noisiness” of the observations.

We organize the trials as follows: first, we generate 100 trajectories using formula (9) for each set of parameters with and . Then we downsample the sequence and work only with each 10-th value of the process, so that and . We estimate the parameters corresponding to the contrast given by (13). We refer to this method as linearized contrast. For the ”quadratic variance” estimation (QV) we do the following: we estimate the parameter explicitly from the observations of the second variable by (17), and then compute the parameters of the drift by minimizing (16). In addition, we compare both methods to the 1.5 strong order scheme (Ditlevsen and Samson, 2017), based on two separate estimators for each coordinate.

The minimization of the criterions is conducted with the optim function in R with the Conjugate Gradient method. In Table 1

we present the mean value of the estimated parameters and the standard deviation (in brackets). Figure

LABEL:FitzHugh illustrates the estimation densities. Linearized estimator is depicted in blue, Quadratic-Variance in red, 1.5 scheme in green.

The first thing is that the estimation of the diffusion coefficient with the 2-dimensional linearized estimator is biased in both sets of data, even though the contrast is corrected with respect to the hypoellipticity of the system. This bias does not appear in the one-dimensional criteria and when the value is directly computed from the observations. Thus its origin may be explained by the dimensionality of the system. Parameters of the second coordinate and are estimated efficiently with all three methods, though the 1.5 scheme provides a more accurate estimation. It is expected, since one of the parameters is fixed to its real value. However, in the case of , 1-dimensional criteria does not score better than the linearized and QV estimators. This parameter seems to be underestimated in the case of 1.5 scheme, and a bit overestimated with the linearization scheme, as well as the diffusion coefficient. During the simulation study it is also observed that is the most sensitive to the initial value with which the optim function is initialized, since it directly regulates the amount of noise which is propagated to the first coordinate.

## 6 Conclusions

The proposed estimator successfully generalizes parametric inference methods developed for models of type (2) for more general class (1). Numerical study shows that it can be used with no prior knowledge of the parameters. It is the most prominent advantage of our method over the analogous works.

From the theoretical point of view, our estimators also reveal good properties. Both linearized and the quadratic variance contrasts are consistent. We did not study the question of the asymptotic normality of the estimator, since in multidimensional case this question is much more harder to treat in comparison to the elliptic case because of the different orders of variance for each coordinate.

The most important direction of prospective work is the adaptation of the estimation method to the case when only the observations of the first coordinate are available. Under proper conditions it must be possible to couple the contrast minimization with one of the existing filtering methods and estimate the parameters of the system (at least, partially). It would allow to face real experimental data.

Another point is the generalization of the contrast to systems of higher dimension. In practice we deal with arbitrary number of rough and smooth variables, and the general rule which describes the behaviour of the contrast in that case is not yet deriven.

## 7 Acknowledgments

Author’s work was financially supported by LabEx PERSYVAL-Lab and LABEX MME-DII. Sincere gratitude is expressed to Adeline Leclercq-Samson and Eva Löcherbach for numerous discussions, helpful remarks and multiple rereadings of the paper draft.

## References

• Cattiaux et al. (2014) Cattiaux, P., León, J. R., and Prieur, C. (2014). Estimation for stochastic damping hamiltonian systems under partial observation. ii drift term. ALEA, 11(1):p–359.
• Cattiaux et al. (2016) Cattiaux, P., León, J. R., Prieur, C., et al. (2016). Estimation for stochastic damping hamiltonian systems under partial observation. iii. diffusion term. The Annals of Applied Probability, 26(3):1581–1619.
• Ditlevsen and Löcherbach (2017) Ditlevsen, S. and Löcherbach, E. (2017). Multi-class oscillating systems of interacting neurons. SPA, 127:1840–1869.
• Ditlevsen and Samson (2017) Ditlevsen, S. and Samson, A. (2017). Hypoelliptic diffusions: discretization, filtering and inference from complete and partial observations. submitted.
• Fitzhugh (1961) Fitzhugh, R. (1961). Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal, 1(6):445–466.
• Florens-Zmirou (1989) Florens-Zmirou, D. (1989). Approximate discrete-time schemes for statistics of diffusion processes. Statistics: A Journal of Theoretical and Applied Statistics, 20(4):547–557.
• Gardiner and Collett (1985) Gardiner, C. W. and Collett, M. J. (1985). Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation. Phys. Rev. A, 31:3761–3774.
• Genon-Catalot and Jacod (1993) Genon-Catalot, V. and Jacod, J. (1993). On the estimation of the diffusion coefficient for multi-dimensional diffusion processes. In Annales de l’IHP Probabilités et statistiques, volume 29, pages 119–151.
• Genon-Catalot et al. (2000) Genon-Catalot, V., Jeantheau, T., and Larédo, C. (2000).

Stochastic volatility models as hidden markov models and statistical applications. bernoulli 6 1051–1079.

Mathematical Reviews (MathSciNet), 10:3318471.
• Hodgkin and Huxley (1952) Hodgkin, A. L. and Huxley, A. F. (1952). A quantitative description of membrane currents and its application to conduction and excitation in nerve. Journal of Physiology-London, 117(4):500–544.
• Jensen (2014) Jensen, A. C. (2014). Statistical Inference for Partially Observed Diffusion Processes. Phd thesis, University of Copenhagen.
• Jensen et al. (2012) Jensen, A. C., Ditlevsen, S., Kessler, M., and Papaspiliopoulos, O. (2012). Markov chain monte carlo approach to parameter estimation in the fitzhugh-nagumo model. Physical Review E, 86(4):041114.
• Karatzas and Shreve (1987) Karatzas, I. and Shreve, S. E. (1987). Brownian Motion and Stochastic Calculus. Graduate Texts in Mathematics. Springer, 1 edition.
• Kessler (1997) Kessler, M. (1997). Estimation of an ergodic diffusion from discrete observations. Scandinavian Journal of Statistics, 24(2):211–229.
• Kloeden et al. (2003) Kloeden, P. E., Platen, E., and Schurz, H. (2003). Numerical solution of SDE through computer experiments. Universitext. Springer.
• Le-Breton and Musiela (1985) Le-Breton, A. and Musiela, M. (1985). Some parameter estimation problems for hypoelliptic homogeneous gaussian diffusions. Banach Center Publications, 16(1):337–356.
• Leon and Samson (2017) Leon, J. R. and Samson, A. (2017). Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate. Annals of Applied Probability. to appear.
• Malliavin and Thalmaier (2006) Malliavin, P. and Thalmaier, A. (2006). Stochastic calculus of variations in mathematical finance. Springer finance. Springer, 1 edition.
• Mattingly et al. (2002) Mattingly, J. C., Stuart, A. M., and Higham, D. J. (2002). Ergodicity for sdes and approximations: locally lipschitz vector fields and degenerate noise. Stochastic processes and their applications, 101(2):185–232.
• Nualart (2006) Nualart, D. (2006). Malliavin Calculus and Related Topics. Springer. New York.
• Ozaki (1989) Ozaki, T. (1989). Statistical identification of nonlinear random vibration systems. Journal of Applied Mechanics, 56:186–191.
• Pokern et al. (2007) Pokern, Y., Stuart, A. M., and Wiberg, P. (2007). Parameter estimation for partially observed hypoelliptic diffusions. J. Roy. Stat. Soc., 71(1):49–73.
• Samson and Thieullen (2012) Samson, A. and Thieullen, M. (2012). Contrast estimator for completely or partially observed hypoelliptic diffusion. Stochastic Processes and their Applications, 122:2521–2552.
• Van der Pol (1920) Van der Pol, B. (1920). A theory of the amplitude of free and forced triode vibrations. Radio Review, 1(1920):701–710.
• Wu (2001) Wu, L. (2001). Large and moderate deviations and exponential convergence for stochastic damping hamiltonian systems. Stochastic processes and their applications, 91(2):205–238.

## 8 Appendix

### 8.1 Properties of the scheme

###### Proof of the Proposition 1.

Let us consider each integral of (8) separately. Denote:

 Wt+Δ=∫t+ΔteJ(t+Δ−s)B(Zt;σ)d~Ws,

where we suppress the dependency of the Jacobian of the starting point on the interval in order to keep notations simple. Recalling the Jacobian of system (3) and the definition of the matrix exponent, we have:

 Wt+Δ=∫t+Δt(I+J(t+Δ−s)+O(Δ2))B(Zt;σ)d~Ws==∫t+Δt[(1+∂xa1(t+Δ−s)∂ya1(t+Δ−s)∂xa2(t+Δ−s)1+∂ya2(t+Δ−s))+O(Δ2)](01)b(Zt;σ)dWs=b(Zt;σ)⎡⎣∂ya1∫t+Δt(t+Δ−s)dWs+O(Δ2)∫t+ΔtdWs+∂ya2∫t+Δt(t+Δ−s)dWs+O(Δ2)⎤⎦

Then we can calculate :

 E[Wt+ΔW′t+Δ]=b2(Zt;σ)E⎛⎝Σ(1)ΔΣ(12)ΔΣ(12)ΔΣ(2)Δ⎞⎠+O(Δ4),

where entries are given by:

 Σ(1)Δ =(∂ya1)2[∫t+Δt(t+Δ−s)dWs]2 Σ(12)Δ =(∂ya1∫t+Δt(t+Δ−s)dWs)(∫t+ΔtdWs+∂ya2∫t+Δt(t+Δ−s)dWs) Σ(2)Δ =(∫t+ΔtdWs+∂ya2∫t+Δt(t+Δ−s)dWs)2

The first entry can be easily calculated by the Itô isometry:

 E[Σ(1)Δ]=(∂ya1)2E[∫t+Δt(t+Δ−s)dWs]2=(∂ya1)2∫t+Δt(t+Δ−s)ds=(∂ya1)2Δ33

Now consider the product of two stochastic integrals in the remaining terms. Assume for simplicity that . From the properties of the stochastic integrals (Karatzas and Shreve (1987)), it is straightforward to see that:

 E⎡⎣limn→∞∑ti,ti−1∈[0,Δ](Δ−s)(Wti−Wti−1)∑ti,ti−1∈[0,Δ](Wti−Wti−1)⎤⎦==limn→∞∑ti,ti−1∈[0,Δ](Δ−s)E[(Wti−Wti−1)2]=∫Δ0(Δ−s)ds=Δ22

That gives the proposition. ∎

### 8.2 Auxiliary results

In this section we introduce an index in the notation for the time step in order to highlight that it depends on the experimental design. Whenever this dependency is not important, the old notations are used. We start with an important Lemma which links the sampling and the probabilistic law of the continuous process:

###### Lemma 3 (Kessler (1997)).

Let and , let be such that is differentiable with respect to and , with derivatives of polynomial growth in uniformly in . Then:

 1NN∑i=1f(Zi;θ)Pθ0⟶∫f(z;θ)ν0(dz) as N→∞ % uniformly in θ.

Lemma is proven in Kessler (1997) for the one-dimensional case. However, as its proof is based only on ergodicity of the process and the assumptions analogous to ours, and not on the discretization scheme or dimensionality, we take it for granted without giving a formal generalization for a multi-dimensional case. Then proposition 2 in combination with the continuous ergodic theorem and Lemma 3 allow us to establish the following important result:

###### Lemma 4.

Let be a function with the derivatives of polynomial growth in , uniformly in . Then:

1. Assume .

 1NΔ3NN−1∑i=0f(Zi;θ)(∂ya1(Zi;θ0))2(Xi+1−¯A1(Zi;θ0))2PΘ⟶13∫f(z;θ)b2(z;σ0)ν0(dz)
2. Assume .

 1NΔNN−1∑i=0f(Zi;θ)(Yi+1−¯A2(Zi;θ0))2PΘ⟶∫f(z;θ)b2(z;σ0)ν0(dz)
3. Assume .

 1NΔ2Nn−1∑i=0f(Zi;θ)(∂ya1(Zi;θ0))(Xi+1−¯A1(Zi;θ0))(Yi+1−¯A2(Zi;θ0))PΘ⟶12∫f(z;θ)b2(z;σ0)ν0(dz)
###### Proof.

Let us denote:

 ζ(1)i =1NΔ3Nf(Zi;θ)(∂ya1(Zi;θ0))2(Xi+1−¯A1(Zi;θ0))2 ζ(2)i =1NΔNf(Zi;θ)(Yi+1−¯A2(Zi;θ0))2 ζ(1,2)i =1NΔ2Nf(Zi;θ)∂ya1(Zi;θ0)(Xi+1−¯A1(Zi;θ0))(Yi+1−¯A2(Zi;θ0))

Thanks to Proposition 2 we know that:

 Eθ0[ζ(1)i|Fi]=13NN−1∑i=0f(Zi;θ)b2(Zi;σ0)+O(ΔN).

Then from Lemma 3 it follows that for uniformly in :

 Eθ0[ζ(1)i|Fi]PΘ⟶13∫f(z;θ)b2(z;σ0)ν0(dz).

The same applies for terms and . ∎

Let us introduce an auxiliary lemma which repeats Lemma 3 in Ditlevsen and Samson (2017). Its proof is based on Lemma 9 from Genon-Catalot and Jacod (1993) and Lemma 3:

###### Lemma 5.

Let be a function with derivatives of polynomial growth in , uniformly in .

1. Assume and . Then

 1NΔ2NN−1∑i=0f(Zi;θ)(Xi+1−¯A1(Zi;θ0))Pθ⟶0

uniformly in .

2. Assume and . Then

 1NΔNN−1∑i=0f(Zi;θ)(Yi+1−¯A2(Zi;θ0))Pθ⟶0

uniformly in .

3. and . Then

 1NN−1∑i=0f(Zi;θ)(Yi+1−¯A2(Zi;θ0))Pθ⟶0

uniformly in .

###### Proof.

Expectation of the first sum tends to zero for and due to Proposition 2. Convergence for is due Lemma 9 in Genon-Catalot and Jacod (1993) and uniformity in follows the proof of Lemma 10 in Kessler (1997). The second and the third assertions are proven in the same way, with the proper scaling. ∎

### 8.3 Consistency of the estimator

###### Proof of Lemma 1.

We can split the contrast in the following sum:

 limN→∞,ΔN→012NLN,ΔN(θ,σ2;Z0:N)=limN→∞,ΔN→0[3T1−3T2+T3+T4]

where terms are given by follows:

 T1 =1NN−1∑i=0(Xi+1−¯A1(Zi;θ))2Δ3Nb2(Zi;σ)(∂ya1)2θ T2 =1NN−1∑i=0(Xi+1−¯A1(Zi;θ))(Yi+1−¯A2(Zi;θ))Δ2Nb2(Zi;σ)(∂ya1)θ T3 =1NN−1∑i=0(Yi+1−¯A2(Zi;θ))2ΔNb2(Zi;σ) T4 =1N