# Optimal control for parameter estimation in partially observed hypoelliptic
stochastic differential equations

Quentin Clairon
University of Bordeaux, Inria Bordeaux Sud-Ouest,
Inserm, Bordeaux Population Health Research Center, SISTM Team, UMR1219,
F-33000 Bordeaux, France
Vaccine Research Institute, F-94000 Créteil, France
and
Adeline Samson
Univ. Grenoble Alpes, CNRS, Inria, Grenoble Institute of Engineering, LJK,
38000 Grenoble, France  The author gratefully acknowledges the LabEx PERSYVAL-
Lab (ANR-11-LABX-0025-01) and MIAI@Grenoble Alpes, (ANR-19-P3IA-0003).

###### Abstract

We deal with the problem of parameter estimation in stochastic differential
equations (SDEs) in a partially observed framework. We aim to design a method
working for both elliptic and hypoelliptic SDEs, the latters being
characterized by degenerate diffusion coefficients. This feature often causes
the failure of contrast estimators based on Euler Maruyama discretization
scheme and dramatically impairs classic stochastic filtering methods used to
reconstruct the unobserved states. All of theses issues make the estimation
problem in hypoelliptic SDEs difficult to solve. To overcome this, we
construct a well-defined cost function no matter the elliptic nature of the
SDEs. We also bypass the filtering step by considering a control theory
perspective. The unobserved states are estimated by solving deterministic
optimal control problems using numerical methods which do not need strong
assumptions on the diffusion coefficient conditioning. Numerical simulations
made on different partially observed hypoelliptic SDEs reveal our method
produces accurate estimate while dramatically reducing the computational price
comparing to other methods.

Keywords: Stochastic differential equations, Parameter estimation,
Hypoellipticity, Optimal control theory

##  1 Introduction

Stochastic differential equations (SDEs) are widely used in several fields to
model stochastic temporal dynamics. We focus on hypoelliptic SDE, that is when
the diffusion matrix is not of full rank but its solution has a smooth
density. More precisely, we consider a process
$Z_{t}=(V_{t},U_{t})\in\mathbb{R}^{d}$, of dimension $d$ on a time interval
$[0,T]$, defined as the solution of the following hypoelliptic SDE:

$\begin{array}[]{l}dV_{t}=g_{\theta}(V_{t},U_{t},t)dt\\\
dU_{t}=h_{\theta}(V_{t},U_{t},t)dt+B_{\sigma}(Z_{t},t)dW_{t}\\\
Z(0)=Z_{0}\end{array}$ (1)

where $V_{t}$ is the $d_{V}$-dimensional component corresponding to the smooth
part of the process, that is no directly affected by stochastic perturbations,
and $U_{t}$ is the $d_{U}$-dimensional rough part, $W_{t}$ is a $d_{U}$
dimensional Brownian motion acting on the system through the $d_{U}$ squared
diffusion matrix $B_{\sigma}(Z_{t},t)$. In some applications, only the smooth
component can be measured.

Let us give some examples of such models. Stochastic damping Hamiltonian
systems or Langevin equations have been developed to describe a particle
moving in an environment defined by a potential Wu (2001). The system is
usually two dimensional, the first equation describes the position of the
particle and the other its velocity. The noise is degenerate because it acts
only on the velocity and not on the position. Only the position is directly
measured. We can cite applications of these models in molecular dynamics
Leimkuhler and Matthews (2015), paleoclimate research Ditlevsen et al. (2002),
neural field models Coombes and Byrne (2019); Ditlevsen and Löcherbach (2017).
Another application is neuronal models of membrane potential dynamics.
Examples are the hypoelliptic FitzHugh-Nagumo (FHN) model Gerstner and Kistler
(2002); DeVille et al. (2005); Leon and Samson (2018), the hypoelliptic
Hodgkin-Huxley model Goldwyn and Shea-Brown (2011); Tuckwell and Ditlevsen
(2016), or synaptic-conductance based models with stochastic channel dynamics
Paninski et al. (2012); Ditlevsen and Greenwood (2013). Only the first
equation, corresponding to the membrane potential, can be measured by intra-
cellular recordings. It is therefore important to develop estimation methods
for this class of models.

In the rest of paper, we assume the drift functions $g_{\theta}$ and
$h_{\theta}$ depend on a $d_{\theta}$ dimensional parameter $\theta$ and the
matrix $B_{\sigma}$ on the $d_{\sigma}$ dimensional parameter $\sigma$ known
as volatility. Partial observations of $(Z_{t})$ are defined as:

$Y_{t}=CZ_{t}$

where $C$ is the $d_{C}\times d$ observation matrix. We assume $(Y_{t})$ is
discretely observed on the interval $\left[0,\,T\right]$ at times
$0=t_{0}<\ldots<t_{n}=T$ and without measurement error. Our aim is to estimate
the unknown parameter $\psi=(\theta,\sigma)$ of model (1) using the discrete
and partial observations $(Y_{0},\ldots,Y_{n})$ and possibly without knowing
the initial condition $Z_{0}.$

Parametric estimation of hypoelliptic models faces several difficulties. When
the first equation reduces to $g_{\theta}(V_{t},U_{t})=U_{t}$, different
contrasts estimators haven been proposed Genon-Catalot et al. (2000); Gloter
(2000); Ditlevsen and Sørensen (2004); Gloter (2006); Samson and Thieullen
(2012); Lu et al. (2016); Leon et al. (2019). This specific case allows to
deal with partial observations, the second coordinate $U_{t}$ being replaced
by increments of $V_{t}$. When the drift $g_{\theta}$ is more complex, we deal
with several entangled issues. First, solutions of the SDE are generally non
explicit. Numerical schemes are used to approximate the discretized process
and derive an estimation criteria. However, compared to elliptic system, the
degeneracy of the noise complicates the statistical estimation. Thus
estimation cost function that would be derived directly from the Euler-
Maruyama fails because the covariance matrix
$\varGamma_{\sigma}^{T}\varGamma_{\sigma}$ where
$\varGamma_{\sigma}(Z_{t},t)=\left(0_{d_{V},d_{U}},B_{\sigma}(Z_{t},t)\right)^{t}$,
is not full rank. Gloter and Yoshida (2021) propose a local Gaussian
approximation of the transition density in the case of totally observed
coordinates. Several papers propose to use a likelihood based on higher order
schemes of approximation, the 1.5 order scheme Ditlevsen and Samson (2019),
the Local Linearization scheme Melnykova (2020). In the case of partial
observations, the unobserved coordinates $U_{t}$ has to be filtered or
imputed, either in a Bayesian settings with a Markov Chain Monte Carlo
algorithm Pokern et al. (2009); Graham et al. (2019) or in a frequentist
settings with a particle filter or sequential Monte Carlo coupled with a
Stochastic Approximation Expectation Maximization (SAEM) algorithm Ditlevsen
and Samson (2014, 2019). A recent simulation filter has been proposed by
Bierkens et al. (2020) for hypoelliptic diffusions. Still the theoretical
assessment and practical efficiency of such numerical methods critically rely
on $\varGamma_{\sigma}$ being well conditioned. These methods are also time
consuming.

In this work, we propose an estimation method based on optimal control theory
which applies to nonlinear hypoelliptic systems and partial observations. The
idea of using optimal control theory is to treat the coordinate filtering
problem as a deterministic tracking problem. Then, numerical devices coming
from control theory allow us to estimate the optimal control bringing the
observed system states the closest possible to the observations. This was
proposed for ODE system Clairon and Brunel (2017, 2019); Clairon (2021); Iolov
et al. (2017). For SDE, the only reference is Clairon and Samson (2020) which
uses linear-quadratic theory, a particular case of the Pontryagin maximum
principle, to solve the tracking problem. But this limits the application of
the method to linear systems. Also, this method separatly estimates $\sigma$
and $\theta$ via a nested procedure because of identifiability issues.
Volatility $\sigma$ was estimated via the minimization of an outer criteria
demanding repeated estimation of $\theta$ based on the optimization of an
inner criteria. This dramatically increases the computation time of the method
when $\sigma$ is unknown. Additionaly, the criteria used for $\sigma$
inference was properly defined only for linear SDEs.

In this work, we construct an estimation procedure which simultaneously
estimate $(\theta,\sigma)$ by relying on a contrast function bypassing
identifiability issues for $\sigma$ estimation encountered in Clairon and
Samson (2020). This contrast function is well defined for linear and nonlinear
SDEs, both elliptics and hypoelliptics. Also, the control problem solutions
required by our approach are easily obtained for non-linear SDEs via an
adaptation of the method proposed by Cimen and Banks (2004a).

Our statistical criteria exploits the hypoelliptic nature of (1) which assumes
”enough” interactions between $V_{t}$ and $U_{t}$ such that the Wiener process
$W_{t}$ influences all the coordinates of the SDE, potentially with some
delay. More precisely, we consider a lagged contrast function linking $W_{t}$
and the observed part of the system state $CZ_{t+m_{B}}$ after a delay
$m_{B}>0$ such that the dependence $CZ_{t+m_{B}}$with respect to $W_{t}$ is
regular enough to ensure a well-conditioned variance for $CZ_{t+m_{B}}$. This
requires to know explicitly the dependence of $CZ_{t+m_{B}}$ to $W_{t}$. For
this we rely on the so-called pseudo-linear representation of the discretized
SDE (1) Cimen (2008). The scheme is expressed at each time $t$ as a function
of the $m_{B}$ previous lagged states. The discretized process is no more
Markovian, this is a direct consequence of hypoellipticity. From this lagged
formulation of the discretized process, we deduce a lagged contrast which is
the pseudo-likelihood of the representation. This statistical contrast depends
on the whole state variable, even the unobserved ones. We thus need to filter
or predict the unobserved coordinate $U_{t}$. We propose a predictor which is
a balance between fidelity to the observations and to the solution of SDE (1).
This predictor is defined via an optimal control problem that can be solved
easily without assumption on $\varGamma_{\sigma}$ rank. The optimal control
problem is solved with the linear quadratic theory for linear SDEs and for
nonlinear models with an adaptation of the procedure developed by Cimen and
Banks (2004b) based on the pseudo-linear representation. Our criteria is
somewhat similar to the Generalized profiling introduced in Ramsay et al.
(2007) for Ordinary differential equations but here dedicated to SDEs with
potentially degenerate diffusions. We also consider a criteria making a
tradeoff between data and model fidelity to regularize the inverse problem of
parameter estimation. This method provides estimators for the parameters and
for the unobserved coordinate $U_{t}$. The simulation study illustrates the
good performances of our method.

Section 2 studies the hypoellipticity of the continuous process, provides a
discrete approximation scheme and the corresponding concept of
hypoellipticity. In section 3, we derive the expression of our statistical
criteria which defines our estimator. The construction of the state predictor
and the optimal control numerical methods are presented in section 4. Then, in
section 5 we conduct a simulation study on three hypoelliptic models to
evaluate the practical accuracy of our method as well as its computational
efficiency.

##  2 Model and its discrete approximation

This Section describes the properties of hypoelliptic systems and the
connexity property that allows to propagate the noise to the smooth
coordinates. Then we introduce the approximate Euler-Maruyama scheme and the
lagged pseudo-linear formulation, the key of our estimation method.

###  2.1 Models and its assumptions

Let us denote:

$f_{\theta}(Z_{t},t)=\left(\begin{array}[]{c}g_{\theta}(V_{t},U_{t},t)\\\
h_{\theta}(V_{t},U_{t},t)\end{array}\right),\,\varGamma_{\sigma}(Z_{t},t)=\left(\begin{array}[]{c}0_{d_{V},d_{U}}\\\
B_{\sigma}(Z_{t},t)\end{array}\right)$

the drift and diffusion coefficient of the diffusion process $Z_{t}$
respectively, then

$dZ_{t}=f_{\theta}(Z_{t},t)dt+\varGamma_{\sigma}(Z_{t},t)dW_{t}.$

We assume $B_{\sigma}$ is full rank and we consider two classes of models
called elliptic and hypoelliptic.

Elliptic models are defined by $d_{V}=0$. In that case, the Brownian process
directly acts on the whole system, the diffusion coefficient is not
degenerated. Girsanov formula applies and the process has a continuous smooth
density.

Hypoelliptic models are defined with $d_{V}>0$ and verify the Hörmander
condition detailed below. Let us first defined the Lie brackets. For $f$ a
function with values in $\mathbb{R}^{d}$, $f_{l}$ stands for its $l$-th
component. For a matrix $\varGamma$, $\varGamma_{j}$ denotes its $j-$th
column.

######  Definition 1.

The Lie bracket of two functions $f,g:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is
defined as

$\left[f,\,g\right]_{l}=\frac{\partial g_{l}(z)}{\partial
z}^{T}f(z)-\frac{\partial f_{l}(z)}{\partial z}^{T}g(z).$

Hypoellipticity is then defined as follows:

######  Definition 2.

(Hörmander condition) Let us consider the set $\mathcal{L}$ defined
iteratively as:

  * •

Initialization step: $\mathcal{L}$ is composed of vectors
$L^{0}=\varGamma_{j}$, for $j=1,\ldots d_{U}$.

  * •

Generalization step at iteration $k$, if $L^{k-1}\in\mathcal{L}$ then vectors
$L^{k}$ defined by $L^{k}=\left[f(Z_{t},t),\,L^{k-1}\right]$ and by
$L^{k}=\left[\varGamma_{j}(Z_{t},t),\,L^{k-1}\right]$ for $j=1,\ldots d_{U}$
belong to $\mathcal{L}$.

If at some iteration, $\mathcal{L}$ spans $\mathbb{R}^{d}$, the weak Hörmander
condition is fulfilled and the system is said hypoelliptic.

###  2.2 Hypoellipticity properties

We derive a necessary condition on the drift $g_{\theta}$ respected by
hypoelliptic systems (the proof is given in Appendix 1 in Supplementary
Materials):

######  Proposition 1.

(Connexity property) If the SDE (1) is hypoelliptic then for each
$l\in\left\llbracket 1,d_{V}\right\rrbracket$, it exists $j\in\left\llbracket
1,d_{U}\right\rrbracket$ and a sequence $\left(q_{1},\ldots
q_{m_{l}}\right)\in\left\llbracket 1,\,d_{V}\right\rrbracket^{m_{l}}$ with
$m_{l}\leq d_{V}$ such that

$\frac{\partial g_{\theta,l}(Z_{t},t)}{\partial V_{q_{1}}}\frac{\partial
g_{\theta,q_{1}}(Z_{t},t)}{\partial V_{q_{2}}}\cdots\frac{\partial
g_{\theta,q_{m_{l}-1}}(Z_{t},t)}{\partial V_{q_{m_{l}}}}\frac{\partial
g_{\theta,q_{m_{l}}}(Z_{t},t)}{\partial U_{j}}\neq 0.$ (2)

This can be seen as a connexity property between the set of smooth variables
and rough ones. It means that each variable $V_{l,t}$ is influenced by at
least one component of $U_{t}$ at least undirectly through a path of auxiliary
smooth variables ($V_{q_{1}},\ldots,V_{q_{m_{l}}}$). Of course, from a given
sequence $(q_{1},\ldots,q_{m_{l}})$, it is possible to construct bigger ones.
To remove ambiguity $m_{l}$ denotes the length of the shortest possible path
such that the connexity property (2) holds.

Although the diffusion coefficient $\varGamma_{\sigma}(Z_{t},t)$ is rank
deficient, the solution of (1) has a smooth density with respect to the
Lebesgue measure.

If the drift $f(Z_{t})$ satisfies a dissipativity condition

$<f(z),z>\leq\alpha-\delta\|z\|^{2},\forall z\in\mathbb{R}^{d}$

where $\alpha,\delta>0$, the function $L(z)=1+\|z\|^{2}$ is a Lyapounov
function. It ensures the geometric ergodicity of the solution of the SDE. It
means that $(Z(t))_{t\in[0,T]}$ converges exponentially fast to a unique
invariant distribution.

However, usually, the solution of SDE (1) is not explicit and a numerical
scheme is needed to approximate the solution.

###  2.3 Euler-Maruyama discretization

Let us introduce the Euler-Maruyama discretized process ($\tilde{Z}_{i}$, for
$i=1,\ldots,n$) of SDE (1). To simplify the notations, we assume that the time
points are equidistant, that is $\triangle_{i}=t_{i+1}-t_{i}=\triangle,\forall
i=1,\ldots,n$. Then the discretized approximate process is defined by:

$\begin{array}[]{l}\tilde{Z}_{i+1}=\tilde{Z}_{i}+\triangle
f_{\theta}(\tilde{Z}_{i},t_{i})+\sqrt{\triangle}\varGamma_{\sigma}(\tilde{Z}_{i},t_{i})u_{i}\\\
\tilde{Z}_{0}=Z_{0}\end{array}$ (3)

where $u_{i}\sim N(0,I_{m})$ is the normalized increment of the Brownian
motion at $t_{i}$, i.e.
$u_{i}=\frac{1}{\sqrt{\triangle}}\left(W_{t_{i+1}}-W_{t_{i}}\right)$ and
$u=\left(u_{0},\ldots,u_{n-1}\right)$.

The connexity property can be illustrated on this Euler-Maruyama
discretisation. Proposition 2 proves that $m_{l}+1$ iterations of (3) are
required to have $\tilde{V}_{l,i+m_{l}+1}$ affected by some entry
$\tilde{U}_{j,i}$ of the rough component.

######  Proposition 2.

Let us assume that for each $l\in\left\llbracket 1,d_{{}_{V}}\right\rrbracket$
it exists a unique minimal length sequence $\left(q_{1},\ldots
q_{m_{l}}\right)$ such that the connexity property holds. Then the Euler-
Maruyama approximation ($\tilde{Z}_{i}$, for $i=1,\ldots,n$) is such that for
each $l\in\left\llbracket 1,d_{{}_{V}}\right\rrbracket$ it exists
$j\in\left\llbracket 1,d_{U}\right\rrbracket$ such that:

$\frac{\partial\tilde{V}_{l,i+k}}{\partial\tilde{U}_{j,i}}=0,\quad\mbox{for
}k=1,\ldots,m_{l}\quad\mbox{and}\quad\frac{\partial\tilde{V}_{l,i+m_{l}+1}}{\partial\tilde{U}_{j,i}}\neq
0$

So, in model (3), each coordinate of $\tilde{V}_{i+m_{B}+1}$ with
$m_{B}:=\max_{l}{m_{l}}$ is influenced by $\tilde{U}_{i}$. It also means that
$m_{B}+1$ is the minimal delay to propagate the noise $u_{i}$ from the rough
components to all the smooth ones in the discretized model (3).

Let us introduce a new observability assumption that reinforces
hypoellipticity for a discrete process. It requires that the observed
variables give enough information about the diffusion process, or equivalently
that the $Y_{t}$ are influenced by all the rough components.

######  Assumption 1.

(H1): The $d_{o}\times d_{U}$ matrix

$C\left(\prod_{l=1}^{m_{B}}\mathbf{A}_{\theta,i+l}\right)\varGamma_{\sigma}(z_{i},t_{i})$

is full rank for all $1\leq i\leq n-m_{B}$ and for all possible values
$\left\\{z_{i}\right\\}_{0\leq l\leq m_{B}}\in\mathbb{R}^{d\times m_{B}}.$

Let us detail this assumption for elliptic and some hypoelliptic models.

  * •

Elliptic models verify $m_{B}=0$. Assumption (H1) imposes that
$C\varGamma_{\sigma}(\tilde{Z}_{i},t_{i})$ is full rank. It ensures a smooth
transition density for the discrete process ($Y_{i}$). This is similar to the
continuous process.

  * •

We focus on hypoelliptic SDE with $\frac{\partial g_{\theta}}{\partial
V}(V_{t},U_{t},t)$ full rank and $Y_{t}=V_{t}$. This corresponds to the case
$m_{B}=1$, that is all the components of $V_{t}$ have a direct interaction
with at least one component of $U_{t}$ and stochasticity only acts on the
unobserved variables. This corresponds to the 1-step hypoellipticity defined
by Buckwar et al. (2021).

###  2.4 Lagged pseudo-linear formulation of the discretized process

Proposition 2 proves the existence of a minimal delay $m_{B}$. For estimation
purpose, we need the explicit link between $\tilde{Z}_{i+m_{B}+1}$ and
$u_{i}$. For this, we rely on what we call later the lagged pseudo-linear
formulation of (3).

Let us introduce the matrix $A_{\theta}(\tilde{Z}_{i},t_{i})$ and the vector
$r_{\theta}(t_{i})$ such that
$A_{\theta}(\tilde{Z}_{i},t_{i})\tilde{Z}_{i}+r_{\theta}(t_{i})=f_{\theta}(\tilde{Z}_{i},t_{i})$.
This allows to reformulate (3) with a pseudo-linear expression:

$\left\\{\begin{array}[]{lll}\tilde{Z}_{i+m_{B}+1}&=&\mathbf{A}_{\theta,i+m_{B}}\,\tilde{Z}_{i+m_{B}}+\triangle
r_{\theta}(t_{i+m_{B}})\\\
&+&\sqrt{\triangle}\,\varGamma_{\sigma}(\tilde{Z}_{i+m_{B}},t_{i+m_{B}})\,u_{i+m_{B}}\\\
\tilde{Z}_{0}&=&Z_{0}\end{array}\right.$ (4)

with $\mathbf{A}_{\theta,i}:=I_{d}+\triangle A_{\theta}(\tilde{Z}_{i},t_{i})$.
We can replace $\tilde{Z}_{i+m_{B}}$ by its state-space expression (3) and
iterate to obtain the lagged pseudo-linear formulation.

######  Proposition 3.

The lagged pseudo-linear formulation of the discretized process (3) is defined
as follows. For $i=0,\ldots,n-m_{B}-1$:

$\begin{array}[]{lll}\tilde{Z}_{i+m_{B}+1}&=&\prod_{r=0}^{m_{B}}\mathbf{A}_{\theta,i+r}\tilde{Z}_{i}\\\
&+&\sum_{r=0}^{m_{B}}\left(\prod_{l=r+1}^{m_{B}}\mathbf{A}_{\theta,i+l}\right)\triangle
r_{\theta}(t_{i+r})\\\
&+&\sum_{r=0}^{m_{B}}\left(\prod_{l=r+1}^{m_{B}}\mathbf{A}_{\theta,i+l}\right)\sqrt{\triangle}\varGamma_{\sigma}(\tilde{Z}_{i+r},t_{i+r})u_{i+r}.\end{array}$
(5)

In the last expression the multiplication symbol $\prod_{l=0}^{m_{B}-1}$
denotes the sequential left multiplication i.e.
$\prod_{j=a}^{b}\mathbf{A}_{\theta,j}=\mathbf{A}_{\theta,b}\times\mathbf{A}_{\theta,b-1}\times\ldots\times\mathbf{A}_{\theta,a}$
with the convention $\prod_{j=b+1}^{b}\mathbf{A}_{\theta,j}=I_{d}.$

The pseudo-linear formulation gives the explicit link between
$\tilde{Z}_{i+m_{B}+1}$ and $u_{i}$. The Markovian property is no more
fulfilled for the discretised process: the smooth density is obtained only if
we express $\tilde{Z}_{i+m_{B}+1}$ as a function of the $m_{B}$ previous
states. This is a direct consequence of hypoellipticity.

##  3 Cost function construction

Our aim is to estimate the unknown parameter $\psi=(\theta,\sigma)$ of model
(1) using the discrete and partial observations $(Y_{0},\ldots,Y_{n})$ and
possibly without knowing the initial condition $Z_{0}.$

We first assume that the process $Z$ is observed and construct an estimation
contrast based on the lagged pseudo-linear representation. Then, we propose to
filter $Z$ by optimal control.

###  3.1 Statistical criteria assuming $Z$ observed

Let us consider in this Section that the variables $Z$ are observed. We can
use formula (5) to build a statistical criteria. Let us denote
$t_{i+m_{B}}^{i}=\left(t_{i},\ldots,t_{i+m_{B}}\right)$,
$\tilde{Z}_{i:i+m_{B}}=\left(\tilde{Z}_{i}\ldots,\tilde{Z}_{i+m_{B}}\right)$
and define

$\begin{array}[]{l}X_{i:i+m_{B}+1}=Y_{i+m_{B}+1}-C\left(\prod_{r=0}^{m_{B}}\mathbf{A}_{\theta,i+r}\right)\tilde{Z}_{i}-\sum_{r=0}^{m_{B}}F_{i+r:i+m_{B}}\\\
F_{i+r:i+m_{B}}=\triangle
C\left(\prod_{l=r+1}^{m_{B}}\mathbf{A}_{\theta,i+l}\right)r_{\theta}(t_{i+r})\\\
G_{i+r:i+m_{B}}=\sqrt{\triangle}C\left(\prod_{l=r+1}^{m_{B}}\mathbf{A}_{\theta,i+l}\right)\varGamma_{\sigma}(Z_{i+r}^{d},t_{i+r}).\end{array}$

Then we have

$\begin{array}[]{l}X_{i:i+m_{B}}=\sum_{r=0}^{m_{B}}G_{i+r:i+m_{B}}u_{i+r}.\end{array}$
(6)

Under assumption (H1), the matrix

$\varSigma_{i:i+m_{B}}=\sum_{r=0}^{m_{B}}G_{i+r:i+m_{B}}G_{i+r:i+m_{B}}^{T}$

is nonsingular. Then, conditionnally on $\tilde{Z}_{i:i+m_{B}}$, we have

$X_{i:i+m_{B}+1}\mid\tilde{Z}_{i:i+m_{B}}\sim
N(0_{d_{o},1},\,\varSigma_{i:i+m_{B}})$

and we can define the following cost function based on the log likelihood:

$H^{m_{B}}(\psi\mid
Y,Z)=-log\left(\mathbb{P}_{\psi}\left\\{X_{i:i+m_{B}+1}\mid\tilde{Z}_{i:i+m_{B}}\right\\}_{i\in\left\llbracket
1,\,n-m_{B}-1\right\rrbracket}\right)$ (7)

and the associated estimator:

$\begin{array}[]{lll}\widehat{\psi}&=&\arg\min_{\psi}H^{m_{B}}(\psi\mid
Y,Z).\end{array}$ (8)

From equation (6), we see that $X_{i:i+m_{B}+1}$ depends at most of
$\left(u_{i},\ldots,\,u_{i+m_{B}}\right)$. Thus the sequence
$(X_{i:i+m_{B}+1})$ for varying $i$ can be composed of dependant terms which
makes $H^{m_{B}}$ not easy to calculate in this general setting. We detail in
the next subsection some cases where $H^{m_{B}}$ has an explicit form.

###  3.2 Simplification of $H^{m_{B}}$ for some models

Note first that when the sequence of $(X_{i:i+m_{B}+1})$ is composed of
independent terms, the expression of $H^{m_{B}}$ simplifies to

$H^{m_{B}}(\psi\mid
Y,Z)=\sum_{i=1}^{n-m_{B}-1}\left(X_{i:i+m_{B}+1}^{T}\varSigma_{i:i+m_{B}}^{-1}X_{i:i+m_{B}+1}+\log\left(\left|\varSigma_{i:i+m_{B}}\right|\right)\right).$
(9)

This is the case for elliptic and 1-step hypoelliptic systems. Let us give
more detailed expressions for $\varSigma_{i:i+m_{B}}$ for these two cases.

####  3.2.1 Elliptic system

Here, we have $m_{B}=0$:

$\begin{array}[]{l}X_{i:i+1}=Y_{i+1}-C\left(\mathbf{A}_{\theta,i}\tilde{Z}_{i}+\triangle
r_{\theta}(t_{i})\right)\end{array}$

and $X_{i:i+1}\mid\tilde{Z}_{i}\sim N(0_{d_{o},1},\,\varSigma_{i:i})\textrm{
with }\varSigma_{i:i}=\triangle
C\varGamma_{\sigma}(\tilde{Z}_{i},t_{i})\varGamma_{\sigma}(\tilde{Z}_{i},t_{i})^{T}C^{T}.$
Thus, we fall back to the classic contrast estimator of elliptic SDE.

####  3.2.2 1-step hypoelliptic system

In the case where $m_{B}=1$ and $V_{t}=CZ_{t}$, we end up with
$C\varGamma_{\sigma}(\tilde{Z}_{i},t_{i})=0$. So condition (H1) is respected
if
$CA_{\theta}(\tilde{Z}_{i+1},t_{i+1})\varGamma_{\sigma}(\tilde{Z}_{i},t_{i})$
is of full rank. Then $X_{i:i+m_{B}+1}$ is given by:

$\begin{array}[]{l}X_{i:i+2}=Y_{i+2}-C\mathbf{A}_{\theta,i+1}\mathbf{A}_{\theta,i}\end{array}\tilde{Z}_{i}-C\left(\triangle\mathbf{A}_{\theta,i+1}r_{\theta}(t_{i})+\triangle
r_{\theta}(t_{i+1})\right)$

and follows the Gaussian law $X_{i:i+2}\mid\tilde{Z}_{i:i+1}\sim
N(0_{d_{o},1},\,\varSigma_{i:i+1})$ with:

$\varSigma_{i:i+1}=\triangle
CA_{\theta}(\tilde{Z}_{i+1},t_{i+1})\varGamma_{\sigma}(\tilde{Z}_{i},t_{i})\varGamma_{\sigma}(\tilde{Z}_{i},t_{i})^{T}A_{\theta}(\tilde{Z}_{i+1},t_{i+1})C^{T}.$

##  4 State predictor $\overline{Z}_{\theta,\sigma}$

Now, we propose a state predictor of the hidden process $Z$ using optimal
control theory and then we plug it in the statistical criteria.

###  4.1 Formalization of the related optimal control problem

In this section, we denote $(\tilde{Z}(Z_{0},u))$ the solution of (4) for a
given sequence $u=\left(u_{0},\ldots,u_{n-1}\right)$ and a given initial
condition $Z_{0}$. Similarly to trajectory fitting estimators Kutoyants
(1991); Dietz (2001), a first predictor could be the solution $\overline{Z}$
of (4) which is the closest to the observations. This solution corresponds to
the sequence $u_{Z_{0}}^{M}$ such that:

$\displaystyle u_{Z_{0}}^{M}$ $\displaystyle=$
$\displaystyle\arg\min_{u}\left\\{\sum_{i=0}^{n}\left\|C\tilde{Z}_{i}(Z_{0},u)-Y_{i}\right\|_{2}^{2}\right\\}$
(10)

when the initial condition $Z_{0}$ is fixed. However, this optimization
problem is ill-posed. The solution is not necessarily unique and would not
always be continuous as a function of the parameters and observations.

We thus propose to introduce a penalization term into the optimization problem
which will lead to a Tikhonov regularized version of it. This is known to
remove sources of ill-posedness Engl et al. (2009), in particular it re-
establishes uniqueness and the continuity of the solution for linear SDEs. The
penalized problem is the following:

$\displaystyle\overline{u}_{Z_{0}}$ $\displaystyle:=$
$\displaystyle\arg\min_{u}\left\\{\sum_{i=0}^{n}\left\|C\tilde{Z}_{i}(Z_{0},u)-Y_{i}\right\|_{2}^{2}-\frac{2}{w}\log
P(u)\right\\}$

where $P(u)$ is the density of the increment $u$. The term $-\frac{2}{w}\log
P(u)$ penalizes the sequences which are unlikely to be a realization of the
Brownian motion and $w>0$ makes the balance between the model and data
fidelity. We have

$\log
P(u)=\log\prod_{i=0}^{n-1}P(u_{i})\varpropto\log\prod_{i=0}^{n-1}e^{-\frac{1}{2}u_{i}^{T}u_{i}}=-\sum_{i=0}^{n-1}\frac{1}{2}u_{i}^{T}u_{i}$

thus

$\overline{u}_{Z_{0}}=\arg\min_{u}\left\\{\sum_{j=1}^{n-1}\left(\left\|C\tilde{Z}_{i}(Z_{0},u)-Y_{i}\right\|_{2}^{2}+\frac{1}{w}u_{i}^{T}u_{i}\right)+\left\|C\tilde{Z}_{n}(Z_{0},u)-Y_{n}\right\|_{2}^{2}\right\\}.$

The predictor corresponds to the solution $\overline{u}_{Z_{0}}$ of the
following deterministic discrete optimal control problem:

$\begin{array}[]{ll}\textrm{Minimize in
$u$:}&C_{w}(u|Y;Z_{0})=\sum_{j=1}^{n-1}\left(\left\|C\tilde{Z}_{i}(Z_{0},u)-Y_{i}\right\|_{2}^{2}+\frac{1}{w}u_{i}^{T}u_{i}\right)+\left\|C\tilde{Z}_{n}(Z_{0},u)-Y_{n}\right\|_{2}^{2}\\\
\textrm{Subject to:
}&\left\\{\begin{array}[]{l}\tilde{Z}_{i+1}(Z_{0},u)=\mathbf{A}_{\theta,i}\tilde{Z}_{i}(Z_{0},u)+\Delta_{i}r_{\theta}(t_{i})+\sqrt{\triangle}\varGamma_{\sigma}(\tilde{Z}_{i}(Z_{0},u),t_{i})u_{i}\\\
\tilde{Z}_{0}=Z_{0}.\end{array}\right.\end{array}$ (11)

Thus, if we are able to solve problem (11) to derive $\overline{u}_{Z_{0}}$,
we can solve (4) to obtain the related state predictor
$\overline{Z}_{Z_{0}}:=\tilde{Z}(Z_{0},\overline{u}_{Z_{0}})$. However, our
optimal control still depends on $Z_{0}$ which is potentially unknown. If
required, its estimation is bypassed by profiling the cost $C_{w}$ and
defining as initial condition estimator:

$\widehat{Z_{0}}=\arg\min_{Z_{0}}\left\\{\min_{u}C_{w}(u|Y;Z_{0})\right\\}.$

Let us now denote $\overline{u}:=\overline{u}_{\widehat{Z_{0}}}$ the
corresponding optimal control and
$\overline{Z}:=\overline{Z}_{\widehat{Z_{0}}}$ the related state-space
predictor.

We propose to solve (11) with numerical control theory methods. For linear
models, (11) is a Linear-Quadratic problem. The control theory then ensures:

  * •

existence and uniqueness of the solutions $\overline{u}$ and $\overline{Z}$,
which are continuous functions of $\left(\theta,\sigma\right)$ and $Y$, under
mild regularity hypothesis on $A_{\theta}$ and $\varGamma_{\sigma}$ (in
particular, $\varGamma_{\sigma}$ is not required to be of full rank),

  * •

computation of $\overline{u}$, $\overline{Z}$ and, if required,
$\widehat{Z_{0}}$ by solving a finite difference equation.

For nonlinear models, we use an adaptation of Cimen and Banks (2004b) which
substitutes the original problem with a finite sequence of Linear-Quadratic
ones to benefit of the points listed above. The linear and nonlinear cases are
detailed in the next subsection.

###  4.2 Numerical methods for solving (11)

We have defined our state predictor as solution of a deterministic tracking
problem. Several algorithms or tools have been developed to solve this kind of
control problem. In this section, $\theta$ is fixed at a given value. Then to
simplify the notations, we omit it.

####  4.2.1 Linear models

We consider the case where
$\mathbf{A}_{i}:=\mathbf{A}_{\theta,i}=I_{d}+\triangle
A_{\theta}(\tilde{Z}_{i},t_{i}):=I_{d}+\triangle A_{\theta}(t_{i})$ and
$\varGamma\left(t_{i}\right):=\varGamma_{\sigma}\left(\tilde{Z}_{i},t_{i}\right)$
for $i=1,\ldots,n$. In this framework $\overline{u}_{Z_{0}}$ is derived by
solving a finite difference equation known as Riccati equation defined by

$\begin{array}[]{lll}E_{n}&=&C^{T}C\\\ h_{n}&=&-C^{T}Y_{n}\end{array}$

and for $i=n-1,\ldots,1$:

$\begin{array}[]{lll}E_{i}&=&\mathbf{A}_{i}^{T}E_{i+1}\mathbf{A}_{i}+C^{T}C-\triangle\mathbf{A}_{i}^{T}E_{i+1}\varGamma(t_{i})G(E_{i+1})\varGamma(t_{i})^{T}E_{i+1}\mathbf{A}_{i}\\\
h_{i}&=&\Delta_{i}\mathbf{A}_{i}^{T}E_{i+1}r(t_{i})+\mathbf{A}_{i}^{T}h_{i+1}-C^{T}Y_{i}\\\
&-&\triangle\mathbf{A}_{i}^{T}E_{i+1}\varGamma(t_{i})G(E_{i+1})\varGamma(t_{i})^{T}\left(h_{i+1}+\Delta_{i}E_{i+1}r(t_{i})\right)\\\
\end{array}$ (12)

with $G(E_{i+1})=\left[\frac{1}{w}\times
I_{d}+\triangle\varGamma(t_{i})^{T}E_{i+1}\varGamma(t_{i})\right]^{-1}$. Then
$\overline{u}_{Z_{0},i}$ is given by:

$\overline{u}_{Z_{0},i}=-\sqrt{\triangle}G(R_{k+1})\varGamma(t_{i})^{T}\left(E_{i+1}\left(\mathbf{A}_{i}\overline{Z}_{i}+\Delta_{i}r(t_{i})\right)+h_{i+1}\right)$
(13)

and $\overline{u}_{Z_{0},i}$ is used in (1) to obtain $\overline{Z}_{Z_{0}}.$
Moreover, the initial condition estimator $\widehat{Z_{0}}$ is given by
$\widehat{Z_{0}}=-\left(E_{0}\right)^{-1}h_{0}.$ Computational details are
given in Appendix 2.

####  4.2.2 Nonlinear models

For non-linear model, we propose to apply the previous algorithm thanks to the
pseudo-linear representation. We replace the original problem by a sequence of
Linear-Quadratic control problems solved iteratively until a convergence
criteria is verified. This method is an adaptation of Cimen and Banks (2004b)
to the case of discrete models. We propose the following algorithm to solve
(11).

  1. 1.

Initialisation $\forall i\in\left\llbracket
0,\,n\right\rrbracket,\,\overline{Z_{i}^{0}}=Z_{0}^{r}$ where $Z_{0}^{r}$ is
an arbitrary starting point if $Z_{0}$ is unknown or
$\overline{Z_{0}^{0}}=Z_{0}$ otherwise.

  2. 2.

At iteration $l$, compute $\left(\overline{Z^{l}},\,\overline{u^{l}}\right)$
by solving the Linear-Quadratic optimal control problem derived from the
pseudo-linear representation of the SDE:

$\begin{array}[]{ll}\textrm{Minimize in
$u$:}&C_{w}^{l}(u|Y;Z_{0})=\sum_{j=1}^{n-1}\left(\left\|C\tilde{Z}_{i}(Z_{0},u)-Y_{i}\right\|_{2}^{2}+\frac{1}{w}u_{i}^{T}u_{i}\right)+\left\|C\tilde{Z}_{n}(Z_{0},u)-Y_{n}\right\|_{2}^{2}\\\
\textrm{Subject to:
}&\left\\{\begin{array}[]{l}\tilde{Z}_{i+1}(Z_{0},u)=\mathbf{A}_{i}^{l}\tilde{Z}_{i}(Z_{0},u)+\Delta_{i}r(t_{i})+\sqrt{\triangle}\varGamma^{l}(t_{i})u_{i}\\\
\tilde{Z}_{0}=Z_{0}\end{array}\right.\end{array}$ (14)

with $\mathbf{A}_{i}^{l}=I_{d}+\triangle A(\overline{Z_{i}^{l-1}},t_{i})$,
$\varGamma^{l}\left(t_{i}\right):=\varGamma\left(\overline{Z_{i}^{l-1}},t_{i}\right)$
and $\overline{Z^{l-1}}$ the state variable corresponding to the optimal
control $\overline{u^{l-1}}$, the minimizer of the cost
$C_{w}^{l-1}(u|Y;Z_{0})$.

  3. 3.

If
$\sum_{i=0}^{n}\left\|\overline{Z_{i}^{l}}-\overline{Z_{i}^{l-1}}\right\|_{2}^{2}<\varepsilon$
then stop otherwise go back to step 2.

  4. 4.

Use $\overline{Z}\simeq\overline{Z^{l}}$ as state variable predictor.

The interest of this algorithm is that at each iteration $l$, problem (14) can
be solved using Linear-Quadratic theory which ensures:

  * •

the existence and uniqueness of the solution problem $\overline{u^{l}}$,

  * •

that $\overline{u^{l}}$ and the corresponding state predictor
$\overline{Z^{l}}$ can be computed by solving the Riccati equation (12) with
$\mathbf{A}_{i}^{l}:=I_{d}+\triangle A_{\theta}(\overline{Z_{i}^{l-1}},t_{i})$
and
$\varGamma^{{}^{l}}\left(t_{i}\right):=\varGamma\left(\overline{Z_{i}^{l-1}},t_{i}\right)$.

###  4.3 Selection of weights $w$

A data driven selection method for $w$ needs to be specified. Let us denote
$\overline{u}_{w}$ the solution of the control problem (11) obtained for a
given weight value $w$ and the corresponding estimator
$\widehat{\psi}=\left(\widehat{\theta},\widehat{\sigma}\right)$. The sequence
$\overline{u}_{w}=(\overline{u}_{w,0},\ldots\overline{u}_{w,n-1})$ is supposed
to mimic increments of a Brownian motion. So ideally
$\left\|\overline{u}_{w,i}\right\|_{2}^{2}\sim\chi^{2}(d_{U})$ with
$\chi^{2}(d_{U})$ the $\chi^{2}$ distribution with $d_{U}$ degrees of freedom.
Thus the i.i.d sequence
$(\left\|\overline{u}_{w,0}\right\|_{2}^{2},\ldots,\left\|\overline{u}_{w,n-1}\right\|_{2}^{2})$
has ideally a density proportional to
$\prod_{i}\left\|\overline{u}_{w,i}\right\|_{2}^{2\left(\frac{d_{U}}{2}-1\right)}e^{-\left\|\overline{u}_{w,i}\right\|_{2}^{2}/2}.$
Based on that, we choose the optimal weight $\widehat{w}$ among a set of
values $W$ which maximizes the external criteria:

$K(w)=\prod_{i}\left\|\overline{u}_{w,i}\right\|_{2}^{2}{}^{\left(\frac{d_{U}}{2}-1\right)}e^{-\left\|\overline{u}_{w,i}\right\|_{2}^{2}/2}.$
(15)

There are some similarities between our method and Generalized Profiling
Ramsay et al. (2007). In both cases, estimation is based on a nested
optimization procedure. 1/ The hyperparameter $w$ balancing the model and data
fidelity is chosen via the minimization of an outer criteria $K$ given by
(15). 2/ For a given $w$, the structural parameter $\psi$ is estimated by
minimizing the middle criteria $H^{m_{B}}$ given by (8). 3/For a given set
$(w,\psi)$, a state variable predictor is computed by optimizing the inner
criteria $C_{w}$ given by (11).

###  4.4 Summary of the estimation procedure

We summarize the method with a pseudo-algorithm formalism:

  * •

Outer criteria: estimation of optimal weight $\widehat{w}$ defined by:

$\widehat{w}:=\arg\min_{w\in W}K(w).$

  * •

Middle criteria: estimate $\widehat{\psi}$ of $\psi$ defined by:

$\widehat{\psi}:=\arg\min_{\psi}H^{m_{B}}(\psi\mid Y,\overline{Z}_{\psi}).$

  * •

Inner criteria: compute state predictor $\overline{Z}_{\psi}$ via algorithm
presented in section 4.2.1 for linear SDEs or section 4.2.2 for nonlinear
ones.

##  5 Simulation study

###  5.1 Experimental design for the simulations

A simulation study is conducted to analyze the practical accuracy of our
method and its computational cost on three partially observed hypoelliptic
systems, one linear and two non-linear ones.

Each system is observed on an interval $[0,T]$ and is sampled at $n$ times
uniformly at every $\frac{T}{(n-1)}$ time points with Euler-Maruyama scheme.
Three different values for the set $(T,n)$ are tested to quantify the effects
of the interval length $T$ and the sample size $n$ on estimation accuracy.
Estimation results are given in terms of empirical bias and variance computed
after Monte- Carlo simulations based on $N_{MC}=1000$ trials. The sample size
has an important impact on the computational efficiency of the method. Thus,
we also give the mean computational time for a given $w$ .

###  5.2 Examples

####  5.2.1 Monotone cyclic feedback system

We consider a neuronal monotone cyclic feedback system proposed in Ditlevsen
and Löcherbach (2017). This model describes the oscillatory behavior of a
system of three populations of neurons in interaction:

$\begin{array}[]{l}dX_{1}(t)=\left(-\nu X_{1}(t)+X_{2}(t)\right)dt\\\
dX_{2}(t)=\left(-\nu X_{2}(t)+X_{3}(t)\right)dt\\\ dX_{3}(t)=-\nu
X_{3}(t)dt+cdW_{t}\end{array}$ (16)

where $X_{1}(t),X_{2}(t)$ and $X_{3}(t)$ are the limit dynamics of each
population.

We consider partial observations with $C=\left(1\,0\,0\right)$. The model is
defined by $g_{\theta}(x_{1},x_{2},x_{3})=\left(\begin{array}[]{l}-\nu
x_{1}+x_{2}\\\ -\nu x_{2}+x_{3}\end{array}\right)$ and
$h_{\theta}(x_{1},x_{2},x_{3})=-\nu x_{3}$,
$A_{\theta}(t)=\left(\begin{array}[]{ccc}-\nu&1&0\\\ 0&-\nu&1\\\
0&0&-\nu\end{array}\right),$ $r_{\theta}(t)=\left(0\,0\,0\right)^{T}$ and
$\varGamma_{\sigma}(t)=\left(0\,0\,c\right)^{T}.$ The system is hypoelliptic.
Applying the generation step of the weak Hörmander condition gives
$L_{0}=\left(0\,0,\,c\right)^{T}$,
$L_{1}=\left[f_{\theta}(Z),L_{0}\right]^{T}=\left(0\,-c\>\nu c\right)^{T}$ and
$L_{2}=\left[f_{\theta}(Z),L_{1}\right]^{T}=\left(c\,0\,\nu c^{2}\right)^{T}$.
Matrix $(L_{0}\,L_{1}\,L_{2})$ spans $\mathbb{R}^{3}.$

Interestingly the need for a second iteration to ensure hypoellipticity is
mimicked by our connexity condition (2) which needs the auxiliary variable
$x_{2}$ to verify $\frac{\partial g_{\theta,1}}{\partial x_{2}}\frac{\partial
g_{\theta,2}}{\partial x_{3}}=1\neq 0$.

Let us now detail assumption (H1). For $m_{B}=2$, we have $G_{i+2:i+m_{B}}=0$,
$G_{i+1:i+m_{B}}=0$ and $G_{i:i+m_{B}}=\sqrt{\triangle}c$. Thus $m_{B}=2$
implies hypothesis (H1) and $X_{i:i+m_{B}}$ only depends on $u_{i}$. So the
estimator is defined as the minimizer of $H^{m_{B}}$ given by equation (8).

A thousand simulations are performed with initial conditions
$\left(X_{1}(0),X_{2}(0),X_{3}(0)\right)=\left(0,0,0\right)$ and true
parameter values set to $\nu=0.2$ and $c=0.15$. Figure 1 illustrates a
simulation on the observation interval
$\left[0,\,T\right]=\left[0,\,20\right].$

Figure 1: Simulated trajectory of the neuronal monotone cyclic feedback system
(16) with parameters $\nu=0.2$ and $c=0.15$. $X_{1}:$ Solid line, $X_{2}:$
dashed line, $X_{3}:$ dotted line.

As said before, the estimation of $\left(\nu,c\right)$ is made from the
observation of $X_{1}$ only. The initial condition are considered known. We
choose $w$ among $W=\left\\{10^{15},10^{20},\,10^{25},\,10^{30}\right\\}.$
Results are given in Table 1 for different values of $T$ and $n$:
$(T,n)=(10,10^{3})$, $(100,10^{3})$, $(10,10^{4})$. The estimators have a good
precision which increases with $T$ and mesh refinement. Also, we point out the
computation time is expressed in terms of second for the case $n=10^{3}$ where
it was expressed in terms of hours in Clairon and Samson (2020) for an
equivalently complex linear model (Harmonic Oscillator model, section 6.2.1).

Table 1: Cyclic feedback model. Estimation of parameters from 1000 simulated trajectories (mean and variance) and mean computational time for a given $w$. | $\nu$ | $c$ | Comp Time
---|---|---|---
True value | 0.2 | 0.15 |
$T=10,n=10^{3}$ | 0.23 (5e-3) | 0.14 (1e-5) | 9s
$T=100,n=10^{3}$ | 0.21 (6e-4) | 0.15 (1e-6) | 19s
$T=10,n=10^{4}$ | 0.22 (2e-3) | 0.15 (1e-6) | 2min50s

####  5.2.2 Hypoelliptic FitzHugh-Nagumo model

We consider now the hypoelliptic neuronal model, which is a minimal
representation of a spiking neuron model such as the Hodgkin-Huxley model Leon
and Samson (2018). It is defined as:

$\begin{array}[]{l}dV_{t}=\frac{1}{\varepsilon}(V_{t}-V_{t}^{3}-U_{t}+s)dt\\\
dU_{t}=(\gamma V_{t}-U_{t}+\beta)dt+\sigma dW_{t}\end{array}$ (17)

where the variable $V_{t}$ represents the membrane potential of a neuron at
time $t$, and $U_{t}$ is a recovery variable, which could represent channel
kinetics. Parameter $s$ is the magnitude of the stimulus current. Often $s$
represents injected current and is thus controlled in a given experiment. It
is therefore assume known and set to $s=0$.

We consider partially observations with $C=\left(1\,0\right)$. The model is
non-linear. Several choices for $A_{\theta}(Z_{t},t)$ are possible, we take
$A_{\theta}(Z_{t},t)=\left(\begin{array}[]{cc}(1-V_{t}^{2})/\varepsilon&-1/\varepsilon\\\
\gamma&-1\end{array}\right)$ and $r_{\theta}(t)=\left(\begin{array}[]{c}0\\\
\beta\end{array}\right)$. Since
$\varGamma_{\sigma}(t)=\left(0\,\sigma\right)^{T}$, we got
$C\varGamma_{\sigma}=0$, $\frac{\partial}{\partial
U}(\frac{1}{\varepsilon}(V_{t}-V_{t}^{3}-U_{t}))=-\frac{1}{\epsilon}\neq 0$
and
$CA_{\theta}(\tilde{Z}_{i+1},t_{i+1})\varGamma_{\sigma}(\tilde{Z}_{i},t_{i})=-\sigma/\varepsilon\neq
0$. The model belongs to the case described in Section (3.2.2) and we can use
the corresponding simplified expression of $H^{m_{B}}$ for the parameter
estimation.

A thousand simulations are performed with initial conditions set to
$\left(V_{0},\,U_{0}\right)=\left(0,0\right)$ and true parameter values
$\left(\varepsilon,\gamma,\beta\right)=\left(0.1,1.5,0.8\right)$ and
$\sigma=0.3$. Figure 2 illustrates a simulation on the observation interval
$\left[0,\,T\right]=\left[0,\,20\right].$

Figure 2:  Simulated trajectory of the neuronal hypoelliptic FitzHugh-Nagumo
model (17) with parameters
$\left(\varepsilon,\gamma,\beta\right)=\left(0.1,1.5,0.8\right)$ and
$\sigma=0.3$. $V:$ solid line, $U:$ dashed line.

The estimation is made from the observation of $V$ only and the initial
conditions are considered unknown and need to be estimated. We choose $w$
among $W=\left\\{10^{16},\,10^{18},\,10^{20},\,10^{25}\right\\}.$

Results are given in Table 2 for $(T,n)=(1,10^{3}),(10,10^{3}),(1,10^{4})$. As
in the previous example, we observe the bias and variance decreasing with $T$
and $n.$ Regarding the computation time and comparative accuracy with other
methods, we recall in Table 3, the results obtained for the case
$T=10,n=10^{3}$ by Clairon and Samson (2020) and Ditlevsen and Samson (2019)
(they are originally presented in Clairon and Samson (2020) table 4, section
6.2.2). We obtain estimation with equivalent or higher accuracy and always
with significantly reduced computational cost.

Table 2: Hypoelliptic FitzHugh-Nagumo model. Estimation of parameters from 1000 simulated trajectories (mean and variance) and mean computational time for a given $w$.  | $\varepsilon$ | $\gamma$ | $\beta$ | $\sigma$ | Comp Time
---|---|---|---|---|---
True value | 0.1 | 1.5 | 0.8 | 0.3 |
$T=1,n=10^{3}$ | 0.11 (2e-3) | 2.13 (1e-1) | 1.18 (3e-1) | 0.32 (2e-2) | 1min10s
$T=10,n=10^{3}$ | 0.09 (2e-5) | 1.58 (6e-2) | 0.87 (5e-2) | 0.29 (2e-4) | 3min15s
$T=1,n=10^{4}$ | 0.10 (4e-4) | 1.90 (8e-2) | 1.10 (2e-1) | 0.31 (3e-3) | 8min30s
Table 3: Hypoelliptic FitzHugh-Nagumo model. Estimation results for the case $T=10,n=10^{3}$ with our method (first row), Clairon and Samson (2020) (second row), Ditlevsen and Samson (2019) third row.  | $\varepsilon$ | $\gamma$ | $\beta$ | $\sigma$ | Comp Time
---|---|---|---|---|---
True value | 0.1 | 1.5 | 0.8 | 0.3 |
$\widehat{\psi}$ | 0.09 (2e-5) | 1.58 (6e-2) | 0.87 (5e-2) | 0.29 (2e-4) | 3min15s
Clairon and Samson (2020) | 0.09 (1e-4) | 1.27 (2e-2) | 0.65 (1e-2) | 0.31 (4e-3) | 45min
Ditlevsen and Samson (2019) | 0.10 (1e-4) | 1.59 (2e-2) | 0.87 (2e-2) | 0.31 (4e-4) | 6h00min

####  5.2.3 Synaptic-conductance model

We consider the conductance-based model with diffusion synaptic input defined
in Ditlevsen and Samson (2019). It describes the voltage dynamics accross the
membrane of a neuron:

$\begin{array}[]{l}C_{c}dV_{t}=(-G_{L}(V_{t}-V_{L})-G_{E,t}(V_{t}-V_{E})-G_{I,t}(V_{t}-V_{I})+I_{inj})dt\\\
dG_{E,t}=-\frac{1}{\tau_{E}}(G_{E,t}-g_{E})+\sigma_{E}\sqrt{G_{E,t}}dW_{E,t}\\\
dG_{I,t}=-\frac{1}{\tau_{I}}(G_{I,t}-g_{I})+\sigma_{I}\sqrt{G_{I,t}}dW_{I,t}\end{array}$
(18)

where $C_{c}$ is the total capacitance, $G_{L}$, $G_{E}$ and $G_{I}$ are the
leak, excitation and inhibition conductances, $V_{L}$, $V_{E}$ and $V_{I}$ are
their respective reversal potentials, and $I_{inj}$ is the injected current.
The conductances $G_{E,t}$ and $G_{I,t}$ are assumed to be stochastic
functions of time, where $W_{E,t}$ and $W_{I,t}$ are two independent Brownian
motions. The square root in the diffusion coefficient ensures that the
conductances stay positive. Parameters $\tau_{E}$, $\tau_{I}$ are time
constants, $g_{E}$, $g_{I}$ are the mean conductances and
$\sigma_{E},\sigma_{I}$ the diffusion coefficients. Here
$U_{t}=(G_{E,t},G_{I,t})$. We assume to know the capacitance, the reversal
potentials, and the injected current.

We consider partial observations with only component $V_{t}$ observed which
corresponds to the observation matrix $C=\left(1\,0\,0\right)$. For the
pseudo-linear representation, we choose
$A_{\theta}(Z_{t},t)=\left(\begin{array}[]{ccc}-G_{L}/C_{c}&-(V_{t}-V_{E})/C_{c}&-(V_{t}-V_{I})/C_{c}\\\
0&-1/\tau_{E}&0\\\ 0&0&-1/\tau_{I}\end{array}\right)$,
$r_{\theta}(t)=\left(\begin{array}[]{c}\left(G_{L}V_{L}+I_{inj}\right)/C_{c}\\\
g_{E}/\tau_{E}\\\ g_{I}/\tau_{I}\end{array}\right).$ Since
$g_{\theta}(V,G_{E},G_{I})=(-G_{L}(V-V_{L})-G_{E}(V-V_{E})-G_{I}(V-V_{I})+I_{inj})/C_{c}$,
$C=\left(1\,0\,0\right)$ and
$\varGamma_{\sigma}(Z_{t},t)=\left(\begin{array}[]{ccc}0&\sigma_{E}\sqrt{G_{E,t}}&0\\\
0&0&\sigma_{I}\sqrt{G_{I,t}}\end{array}\right)^{T},$ we got
$C\varGamma_{\sigma}=\left(0\,0\right)$, $\frac{\partial
g_{\theta}}{\partial\left(G_{E},G_{I}\right)}(V,G_{E},G_{I})$ is of full rank
and
$CA_{\theta}(\tilde{Z}_{i+1},t_{i+1})\varGamma_{\sigma}(\tilde{Z}_{i},t_{i})\neq
0$. It enters the framework of section (3.2.2).

A thousand simulations are performed with initial conditions set to
$\left(V_{0},\,G_{E,0},\,G_{I,0}\right)=\left(-60,\,10,\,1\right)$ and true
parameters values set to
$\left(G_{L},V_{L},V_{E},V_{I},I_{inj},g_{E}\right)=\left(50,-70,0,-80,-60,17.8\right)$,
$\left(\tau_{E},\,\tau_{I},\,g_{I}\right)=\left(0.5,1,9.4\right)$ and
$\left(\sigma_{E},\sigma_{I}\right)=\left(0.1,\,0.1\right)$. The estimation of
$\left(\tau_{E},\,\tau_{I},\,g_{I},\sigma_{E},\sigma_{I}\right)$ is made from
the observation of $V$ only. For the sake of identifiability, the initial
conditions are assumed known. We choose $w$ among $W=\left\\{10^{8},5\times
10^{8},\,10^{9},\,5\times 10^{9}\right\\}.$

Results are given in Table 4 for $(T,n)=(20,10^{3}),(200,10^{3}),(20,10^{4})$.
We observe a difference in terms of accuracy between the estimation of
$\theta=(\tau_{E},\tau_{I},g_{I})$ and
$\sigma=\left(\sigma_{E},\sigma_{I}\right)$, the latter suffering from more
bias than the former. This is understandable and already noticed in Ditlevsen
and Samson (2019). Indeed we have to estimate the diffusion of two sources of
stochastic disturbance from only one resulting signal. Because of the number
of involved state variables and parameters to estimate, the ratio of
observed/unobserved states and diffusion matrix structure, This estimation
problem is more complex than the previous ones as the number of parameters to
estimate is higher and the ratio of observed/unobserved states is also higher.
This explains the higher computational time for the case $T=20,n=10^{4}$.

Table 4: Synaptic-conductance model. Estimation of parameters from 1000 simulated trajectories (mean and variance) and mean computational time for a given $w$.  | $\tau_{E}$ | $\tau_{I}$ | $g_{I}$ | $\sigma_{E}$ | $\sigma_{I}$ | Comp Time
---|---|---|---|---|---|---
True value | 0.5 | 1 | 9.4 | 0.1 | 0.1 |
$T=20,n=10^{3}$ | 0.51 (2e-3) | 1.11 (0.04) | 9.45 (0.02) | 0.05 (3e-4) | 0.15 (1e-3) | 1min40s
$T=200,n=10^{3}$ | 0.47 (4e-4) | 1.11 (0.01) | 9.41 (2e-3) | 0.09 (7e-6) | 0.14 (8e-5) | 2min00s
$T=20,n=10^{4}$ | 0.52 (2e-3) | 1.08 (0.04) | 9.40 (0.02) | 0.11 (4e-5) | 0.07 (1e-3) | 1h40min

##  6 Discussion

In this work, we propose a new estimation method which gives an unifying
framework for elliptic and hypoelliptic systems, partially or fully observed.
For this, we rely on a lagged discretization of the original SDE which let
enough time to the stochastic perturbations to affect all the state variables.
By doing so, we have constructed a statistical criteria based on a well-
defined density even for partially observed hypoelliptic SDEs. This criteria
requires a state variable predictor obtained by solving a control problem
balancing data and model fidelity. The numerical procedures used to solve it
do not require $B_{\sigma}(Z_{t},t)$ to be of full rank and so are well
adapted to hypoelliptic systems. Because of this, the deterministic control
perspective constitutes a relevant alternative to MCMC approaches and explain
the reasonable computational cost of our method. It only provides the
pointwise estimator $\overline{Z}$ we need for $\psi$ estimation and does not
aim to reconstruct its whole distribution. Now, we conclude this work by
presenting refinments and extensions of the presented method which will be
investigated in the future.

Our method requires to select an hyperparameter $w$. To bypass this, we aim to
define in future works our predictor $\overline{u}$ as a maximum a posteriori
estimator (MAP) in a functional space. For this, we will rely on the work of
Dashti et al. (2013) in which the MAP is defined as the solution of a new
deterministic optimal control problem where our regularization term
$\frac{2}{w}\ln P(u)$ is replaced by $\left\|u\right\|_{E}.$ Here, $E$ denotes
the so-called Cameron-Martin space which is totally determined by the
covariance operator of the Brownian motion $u$ and does not involve a nuisance
parameter.

Also, we think our method is well suited for generalization to SDEs driven by
processes $W_{t}$ different from the Wiener ones. All it requires would be to
know their densities to modify $H^{m_{B}}$ accordingly and the penalization
term appearing in the control problem. Interestingly in this general setting,
to a given type of SDE will be associated a given type of deterministic
optimal control problem.

We end this work by going one step further about the connection between our
estimation problem and control theory by pointing out its similarity with the
issue of structural assessment of a system controllability as exposed in
Daoutidis and Kravaris (1992). Their problem is the following: given an ODE
$\dot{x}=f(x)+Bu$, can we control the behavior of some system outputs
$y_{i}=h_{i}(x)$ knowing that $u=\left\\{u_{j}\right\\}_{j\in\left\llbracket
1,d_{u}\right\rrbracket}$ only affect directly a subset of $x$? This leads to
define in a similar way as our $m_{l}$’s the integers $r_{ij}$, named relative
orders, quantifying the sluggishness of $y_{i}$ response to $u_{j}$ variation.
Then, conditions for controllability are formulated via the non-singularity of
a so-called characteristic matrix $C(x)$ constructed from the $r_{ij}$’s and
Lie derivative in a way mirroring the matrix $M$ appearing in the proof of
proposition 1. From this, we hope to derive $m_{B}$ in a less exploratory
manner directly from a given model structure and make explicit its link with
the required number of iteration of generalization step 2 to fulfill Hörmander
condition.

## Code/Software

Our estimation method is implemented in R and a code reproducing the examples
of Section 5 is available on a GitHub repository located here.

## References

  * Bierkens et al. (2020) Bierkens, J., F. van der Meulen, and M. Schauer (2020).  Simulation of elliptic and hypo-elliptic conditional diffusions.  Adv Applied Proba 52.
  * Buckwar et al. (2021) Buckwar, E., A. Samson, M. Tamborrino, and I. Tubikanec (2021).  Splitting methods for sdes with locally lipschitz drift. an illustration on the fitzhugh-nagumo model.  arXiv preprint arXiv:2101.01027.
  * Cimen (2008) Cimen, T. (2008).  State-dependent riccati equation (sdre) control: A survey.  IFAC Proceedings 41, 3761–3775.
  * Cimen and Banks (2004a) Cimen, T. and S. Banks (2004a).  Global optimal feedback control for general nonlinear systems with nonquadratic performance criteria.  Systems and Control Letters 53, 327–346.
  * Cimen and Banks (2004b) Cimen, T. and S. Banks (2004b).  Nonlinear optimal tracking control with application to super-tankers for autopilot design.  Automatica 40, 1845–1863.
  * Clairon (2021) Clairon, Q. (2021).  A regularization method for the parameter estimation problem in ordinary differential equations via discrete optimal control theory.  Journal of Statistical Planning and Inference 210, 1–19.
  * Clairon and Brunel (2017) Clairon, Q. and N. Brunel (2017).  Optimal control and additive perturbations help in estimating ill-posed and uncertain dynamical systems.  JASA.
  * Clairon and Brunel (2019) Clairon, Q. and N. J.-B. Brunel (2019).  Tracking for parameter and state estimation in possibly misspecified partially observed linear ordinary differential equations.  Journal of Statistical Planning and Inference 199, 188–206.
  * Clairon and Samson (2020) Clairon, Q. and A. Samson (2020).  Optimal control for estimation in partially observed elliptic and hypoelliptic linear stochastic differential equations.  Statistical Inference for Stochastic Processes 23(1), 105–127.
  * Coombes and Byrne (2019) Coombes, S. and A. Byrne (2019).  Lecture Notes in Nonlinear Dynamics in Computational Neuroscience: from Physics and Biology to ICT, Chapter Next generation neural mass models.  PoliTO Springer Series. Springer.
  * Daoutidis and Kravaris (1992) Daoutidis, P. and C. Kravaris (1992).  Structural evaluation of control configurations for multivariable nonlinear processes.  Chemical engineering science 47(5), 1091–1107.
  * Dashti et al. (2013) Dashti, M., K. J. Law, A. M. Stuart, and J. Voss (2013).  Map estimators and their consistency in bayesian nonparametric inverse problems.  Inverse Problems 29(9), 095017.
  * DeVille et al. (2005) DeVille, R., E. Vanden-Eijnden, and C. Muratov (2005, SEP).  Two distinct mechanisms of coherence in randomly perturbed dynamical systems.  Physical Review E 72(3, 1).
  * Dietz (2001) Dietz, H. (2001).  Asymptotic behaviour of trajectory fitting estimators for certain non-ergodic sde.  Statistical Inference for Stochastic Processes 4, 249–258.
  * Ditlevsen et al. (2002) Ditlevsen, P., S. Ditlevsen, and K. Andersen (2002).  The fast climate fluctuations during the stadial and interstadial climate states.  Annals of Glaciology 35, 457–462.
  * Ditlevsen and Greenwood (2013) Ditlevsen, S. and P. Greenwood (2013, AUG).  The Morris-Lecar neuron model embeds a leaky integrate-and-fire model.  Journal of Mathematical Biology 67(2), 239–259.
  * Ditlevsen and Löcherbach (2017) Ditlevsen, S. and E. Löcherbach (2017).  Multi-class oscillating systems of interacting neurons.  Stochastic Processes and Their Application 127, 1840–1869.
  * Ditlevsen and Löcherbach (2017) Ditlevsen, S. and E. Löcherbach (2017).  Multi-class oscillating systems of interacting neurons.  Stochastic Processes and Their Applications 127, 1840–1869.
  * Ditlevsen and Samson (2014) Ditlevsen, S. and A. Samson (2014).  Estimation in the partially observed stochastic morris-lecar neuronal model with particle filter and stochastic approximation methods.  Annals of Applied Statistics 2, 674–702.
  * Ditlevsen and Samson (2019) Ditlevsen, S. and A. Samson (2019).  Hypoelliptic diffusions: discretization, filtering and inference from complete and partial observations.  J Royal Statistical Society B 81, 361–384.
  * Ditlevsen and Sørensen (2004) Ditlevsen, S. and M. Sørensen (2004).  Inference for observations of integrated diffusion processes.  Scand. J. Statist. 31(3), 417–429.
  * Engl et al. (2009) Engl, H. W., C. Flamm, P. Kügler, J. Lu, S. Müller, and P. Schuster (2009).  Inverse problems in systems biology.  Inverse Problems 25(12).
  * Genon-Catalot et al. (2000) Genon-Catalot, V., T. Jeantheau, and C. Larédo (2000, 12).  Stochastic volatility models as hidden markov models and statistical applications.  Bernoulli 6(6), 1051–1079.
  * Gerstner and Kistler (2002) Gerstner, W. and W. Kistler (2002).  Spiking Neuron Models.  Cambridge University Press.
  * Gloter (2000) Gloter, A. (2000).  Discrete sampling of an integrated diffusion process and parameter estimation of the diffusion coefficient.  ESAIM Probab. Statist. 4, 205–227.
  * Gloter (2006) Gloter, A. (2006).  Parameter estimation for a discretely observed integrated diffusion process.  Scand. J. Statist. 33(1), 83–104.
  * Gloter and Yoshida (2021) Gloter, A. and N. Yoshida (2021).  Adaptive estimation for degenerate diffusion process.  Electronic J Stat 15, 1424–1472.
  * Goldwyn and Shea-Brown (2011) Goldwyn, J. H. and E. Shea-Brown (2011, NOV).  The What and Where of Adding Channel Noise to the Hodgkin-Huxley Equations.  PLOS Computational Biology 7(11).
  * Graham et al. (2019) Graham, M. M., A. H. Thiery, and A. Beskos (2019).  Manifold markov chain monte carlo methods for bayesian inference in a wide class of diffusion models.  arXiv: Computation.
  * Iolov et al. (2017) Iolov, A., S. Ditlevsen, and A. Longtin (2017).  Optimal design for estimation in diffusion processes from first hitting times.  SIAM J. Uncertainty Quantification 5, 88–110.
  * Kutoyants (1991) Kutoyants, Y. (1991).  Minimum distance parameter estimation for diffusion type observation.  Comptes rendus de l Academie des sciences 312, 637.
  * Leimkuhler and Matthews (2015) Leimkuhler, B. and C. Matthews (2015).  Molecular Dynamics with deterministic and stochastic numerical methods, Volume 39 of Interdisciplinary Applied Mathematics.  Springer International Publishing Switzerland.
  * Leon et al. (2019) Leon, J., L. Rodriguez, and R. Ruggiero (2019).  Consistency of a likelihood estimator for stochastic damping hamiltonian systems. totally observed data.  ESAIM PS 23, 1–36.
  * Leon and Samson (2018) Leon, J. and A. Samson (2018).  Hypoelliptic stochastic fitzhugh-nagumo neuronal model: mixing, up-crossing and estimation of the spike rate.  Annals of Applied Probability 28(4), 2243–2274.
  * Lu et al. (2016) Lu, F., K. Lin, and A. Chorin (2016).  comparison of continuous and discrete-time data-based modeling for hypoelliptic systems.  Comm App Math Comp Sci 11, 187–216.
  * Melnykova (2020) Melnykova, A. (2020).  Parametric inference for hypoelliptic ergodic diffusions with full observations.  Statistical Inference for Stochastic Processes 23, 595–635.
  * Paninski et al. (2012) Paninski, L., M. Vidne, B. DePasquale, and D. Fereira (2012).  Inferring synaptic inputs given a noisy voltage trace via sequential Monte Carlo methods.  Journal of Computational Neuroscience 33(1), 1–19.
  * Pokern et al. (2009) Pokern, Y., A. Stuart, and P. Wiberg (2009).  Parameter estimation for partially observed hypoelliptic diffusions.  J. Roy. Stat. Soc. B 71(1), 49–73.
  * Ramsay et al. (2007) Ramsay, J., G. Hooker, J. Cao, and D. Campbell (2007).  Parameter estimation for differential equations: A generalized smoothing approach.  Journal of the Royal Statistical Society (B) 69, 741–796.
  * Samson and Thieullen (2012) Samson, A. and M. Thieullen (2012).  Contrast estimator for completely or partially observed hypoelliptic diffusion.  Stochastic Processes and Their Applications 122, 2521–2552.
  * Tuckwell and Ditlevsen (2016) Tuckwell, H. C. and S. Ditlevsen (2016, OCT).  The Space-Clamped Hodgkin-Huxley System with Random Synaptic Input: Inhibition of Spiking by Weak Noise and Analysis with Moment Equations.  Neural Computation 28(10), 2129–2161.
  * Wu (2001) Wu, L. (2001).  Large and moderate deviations and exponential convergence for stochastic damping hamiltonian systems.  Stochastic Process. Appl. 91, 205–238.