# Parametric Estimation of Ordinary Differential Equations with Orthogonality Conditions

Nicolas J.-B. Brunel, Quentin Clairon, Florence d'Alché-Buc

► To cite this version:

Nicolas J.-B. Brunel, Quentin Clairon, Florence d'Alché-Buc. Parametric Estimation of Ordinary Differential Equations with Orthogonality Conditions. 2013. [hal-00713368v2](https://hal.science/hal-00713368v2)

HAL Id: hal-00713368

https://hal.science/hal-00713368v2

Preprint submitted on 28 Oct 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Parametric Estimation of Ordinary Differential Equations with
Orthogonality Conditions

Nicolas J-B. Brunel<fnref n="1" />, Quentin Clairon<fnref n="2" />, Florence d'Alché-Buc<fnref n="3" />,<fnref n="4" />

4

Abstract

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

Differential equations are commonly used to model dynamical deterministic systems in appli-
cations. When statistical parameter estimation is required to calibrate theoretical models to data,
classical statistical estimators are often confronted to complex and potentially ill-posed optimization
problem. As a consequence, alternative estimators to classical parametric estimators are needed for
obtaining reliable estimates. We propose a gradient matching approach for the estimation of para-
metric Ordinary Differential Equations observed with noise. Starting from a nonparametric proxy of
a true solution of the ODE, we build a parametric estimator based on a variational characterization
of the solution. As a Generalized Moment Estimator, our estimator must satisfy a set of orthogonal
conditions that are solved in the least squares sense. Despite the use of a nonparametric estimator,
we prove the root-*n* consistency and asymptotic normality of the Orthogonal Conditions estimator.
We can derive confidence sets thanks to a closed-form expression for the asymptotic variance, and
we give a practical way to optimize the variance by adaptive reweighting. Finally, we compare
our estimator in several experiments in order to show its versatility and relevance with respect to
classical Gradient Matching and Nonlinear Least Squares estimators.

**Key-words:** Gradient Matching, Non-parametric statistics, Plug-in Property, Variational formulation, Sobolev Space.

<footnote n="1">¹ ENSIIE & Laboratoire Statistique et Génome, Université d'Évry Val d'Essonne, UMR CNRS 8071 - USC INRA - FRANCE</footnote>

<footnote n="2">² Laboratoire Analyse et Probabilités, Université d'Évry Val d'Essonne - FRANCE</footnote>

<footnote n="3">³ Laboratoire IBISC, Université d'Évry Val d'Essonne - FRANCE</footnote>

<footnote n="4">⁴ INRIA-Saclay, LRI, Université Paris Sud, UMR CNRS 8623 - FRANCE</footnote>
# Supplementary Material II: Experiments - Tables and Figures

## Contents

<table><tbody><tr><td><strong>1 Linear ODE: α-pinene</strong></td><td style="text-align: right;"><strong>3</strong></td></tr><tr><td style="padding-left: 2em;">1.1 Known initial condition</td><td style="text-align: right;">3</td></tr><tr><td style="padding-left: 2em;">1.2 Unknown initial condition</td><td style="text-align: right;">4</td></tr><tr><td><strong>2 Nonlinear ODEs</strong></td><td style="text-align: right;"><strong>5</strong></td></tr><tr><td style="padding-left: 2em;">2.1 Ricatti equation</td><td style="text-align: right;">5</td></tr><tr><td style="padding-left: 4em;">2.1.1 Known change-point time T<sub>r</sub></td><td style="text-align: right;">5</td></tr><tr><td style="padding-left: 4em;">2.1.2 Unknown change-point time T<sub>r</sub></td><td style="text-align: right;">8</td></tr><tr><td style="padding-left: 2em;">2.2 FitzHugh-Nagumo equation</td><td style="text-align: right;">9</td></tr><tr><td><strong>3 Real data analysis</strong></td><td style="text-align: right;"><strong>11</strong></td></tr><tr><td style="padding-left: 2em;">3.1 Influenza virus growth and migration model</td><td style="text-align: right;">11</td></tr><tr><td style="padding-left: 2em;">3.2 Blowfly models</td><td style="text-align: right;">13</td></tr></tbody></table>

## Introduction

The supplementary material contains additional figures and tables to the paper "Parametric Estimation of Ordinary Differential Equations with Orthogonality Conditions". It gives more information about the results of the Monte Carlo simulations for the comparison of the TS, NLS and OC estimators. Indeed, for the evaluation of estimator accuracy, we compute the Absolute Relative Error (ARE) defined by

$$ARE = \frac{1}{N_{MC}} \sum_{i=1}^{N_{MC}} \frac{|\theta^* - \hat{\theta}_i|}{|\theta^*|}.$$

Additionally, we compute the weighted Orthogonal Condition estimators obtained by minimizing the criterion $Q_{n,L}^W(\theta) = e_L(\hat{\phi}, \theta)^\top W e_L(\hat{\phi}, \theta)$, where the weight matrix computed by the IRWOC
# Supplementary Material I:

## Proofs of results in “Parametric estimation of Ordinary Differential Equations with Orthogonality Conditions”

## Contents

1 Ordinary Differential Equations 1

2 Assumptions 1

3 Consistency 2

4 Asymptotics 4

5 From optimal Weighting to a practical algorithm 7

### 1 Ordinary Differential Equations

**Cauchy-Lipschitz theory for non-autonomous ODE.** If the vector $(t, x) \mapsto f(t, x, \theta)$ satisfies the following conditions on $[0, 1] \times \mathcal{X}$:

(a) $f$ is $L^2$-Caratheodory, i.e. $x \mapsto f(t, x, \theta)$ is continuous for $t$ a.e in $[0, 1]$, $t \mapsto f(t, x, \theta)$ is measurable for all $x \in \mathcal{X}$ and $\forall c > 0, \exists h_c(\cdot, \theta) \in L^2 \setminus |x| \le c \Rightarrow |f(t, x, \theta)| \le h_c(t, \theta)$,

(b) $f$ is $L^2$-Lipschitz, i.e. $\exists a(\cdot, \theta) \in L^2 / \forall x, x' \in \mathcal{X}, |f(t, x, \theta) - f(t, x', \theta)| \le a(t, \theta) |x - x'|$.

Then, there exists a unique solution to the Initial Value Problem (theorem 3.4 in [4])

$$ \begin{cases} \dot{\phi}(t) = \mathbf{f}(t, \phi(t), \theta) \\ \phi(0) = \phi_0 \end{cases} \qquad (1) $$

### 2 Assumptions

**Condition C1 (a)** $\Theta$ is a compact set of $\mathbb{R}^p$ and $\theta^*$ is an interior point of $\Theta$, $\mathcal{X}$ is an open subset of $\mathbb{R}^2$; (b) $(t, x) \mapsto f(t, x, \theta^*)$ is $L^2$-Lipschitz and $L^2$-Caratheodory.

**Condition C2 (a)** $(Y_i, t_i)$ are i.i.d. with variance $V(Y|T=t) = \Sigma_\epsilon = \sigma^2 I_2$; (b) For every $K$, there is a non-singular constant matrix $B$ such that for $P^K = B_p^K(t)$; (i) the smallest eigenvalue of $E[P^K(T)P^K(T)^\top]$ is bounded away from zero uniformly in $K$ and (ii) there is a sequence of constants $\zeta_0(K)$ satisfying $\sup_t |P^K(t)| \le \zeta_0(K)$ and $K = K(n)$ such that $\zeta_0(K)^2 K/n \to 0$ as $n \to \infty$; (c) There are $\alpha, \mathbf{c}_{1,K}, \mathbf{c}_{2,K}$ such that $\|\phi_j^* - p^K \mathbf{c}_{j,K}\|_\infty = \sup_{[0,1]} |\phi_j^*(t) - p^K(t)^\top \mathbf{c}_{j,K}| = O(K^{-\alpha})$.

**Condition C3** There exists $D > 0$, such that the $D$-neighborhood of the solution range $\mathcal{D} = \{x \in \mathbb{R}^2 | \exists t \in [0, 1], |x - \phi^*(t)| < D\}$ is included in $\mathcal{X}$ and $f$ is $C^2$ in $(x, \theta)$ on $\mathcal{D} \times \Theta$ for $t$ in $[0, 1]$ a.e. Moreover, the derivatives of $f$ w.r.t $x$ and $\theta$ (with obvious notations) $f_x, f_\theta, f_{xx}, f_{x\theta}$ and $f_{\theta\theta}$ are $L^2$ uniformly bounded on $\mathcal{D} \times \Theta$ by $L^2$ functions $\bar{h}_x, \bar{h}_\theta, \bar{h}_{x\theta}, \bar{h}_{xx}$ and $\bar{h}_{\theta\theta}$ (respectively).

**Condition C4** Let $(\varphi_\ell)_{\ell \ge 1}$ be an orthonormal sequence of $C^1$ functions in $H_0^1$.

**Condition C5** $\theta^*$ is the unique global minimizer of $Q_\mathcal{F}^*$ and $\inf_{|\theta - \theta^*| > \epsilon} Q_\mathcal{F}^*(\theta) > 0$.
# Parametric Estimation of Ordinary Differential Equations with Orthogonality Conditions

Nicolas J-B. Brunel¹, Quentin Clairon², Florence d'Alché-Buc³,⁴

May, 04 2013

## Abstract

Differential equations are commonly used to model dynamical deterministic systems in applications. When statistical parameter estimation is required to calibrate theoretical models to data, classical statistical estimators are often confronted to complex and potentially ill-posed optimization problem. As a consequence, alternative estimators to classical parametric estimators are needed for obtaining reliable estimates. We propose a gradient matching approach for the estimation of parametric Ordinary Differential Equations observed with noise. Starting from a nonparametric proxy of a true solution of the ODE, we build a parametric estimator based on a variational characterization of the solution. As a Generalized Moment Estimator, our estimator must satisfy a set of orthogonal conditions that are solved in the least squares sense. Despite the use of a nonparametric estimator, we prove the root-n consistency and asymptotic normality of the Orthogonal Conditions estimator. We can derive confidence sets thanks to a closed-form expression for the asymptotic variance. Finally, the OC estimator is compared to classical estimators in several (simulated and real) experiments and ODE models in order to show its versatility and relevance with respect to classical Gradient Matching and Nonlinear Least Squares estimators. In particular, we show on a real dataset of influenza infection that the approach gives reliable estimates. Moreover, we show that our approach can deal directly with more elaborated models such as Delay Differential Equation (DDE).

**Key-words:** Gradient Matching, Nonparametric statistics, Methods of Moments, Plug-in Property,
# 1 Introduction

Differential Equations are a standard mathematical framework for modeling dynamics in physics, chemistry, biology, engineering sciences, etc and have proved their efficiency in describing the real world. A classical model is the Ordinary Differential Equation defined for $t \in [0, 1]$,

$$ \dot{x}(t) = \mathbf{f}(t, x(t), \theta) \qquad (1.1) $$

where dot indicates derivative with respect to time. $\mathbf{f}$ is a time-dependent vector field from $\mathcal{X} \subset \mathbb{R}^d$ to $\mathbb{R}^d$ which is parametrized by a parameter $\theta \in \Theta \subset \mathbb{R}^p$, $d, p \ge 1$. An important task is then the estimation of the parameter $\theta$ from real data. [21] proposed a significant improvement to this statistical problem, and gave motivations for further statistical studies. We are interested in the definition and in the optimality of a statistical procedure for the estimation of the parameter $\theta$ from noisy observations $y_1, \dots, y_n \in \mathbb{R}^d$ of a solution at times $t_1 < \dots < t_n$.

Most works deal with Initial Value Problems (IVP), i.e. with ODE models having a given (possibly unknown) initial value $x(0) = x_0$. There exists then a unique solution $\phi(\cdot, x_0, \theta)$ of (1.1) defined on the interval $[0, 1]$, that depends smoothly on $x_0$ and $\theta$. The estimation of $\theta$ is a classical problem of nonlinear regression, where we regress $y$ on the time $t$. If $x_0$ is known, the Nonlinear Least Square Estimator $\hat{\theta}^{NLS}$ (NLSE) is obtained by minimizing

$$ Q_n^{LS}(\theta) = \sum_{i=1}^{n} |y_i - \phi(t_i, x_0, \theta)|^2 \qquad (1.2) $$

where $|\cdot|$ is the classical Euclidean norm. The NLSE, Maximum Likelihood Estimator or more general M-estimators [23] are commonly used because of their good statistical properties (root-$n$ consistency, asymptotic efficiency), but they come with important computational difficulties (repeated ODE integrations and presence of multiple local minima) that can decrease their interest. We refer to [21] for a detailed overview of the previous works in this field and of some of the problems encountered in estimating ODE. Global optimization methods are then often used, such as simulated annealing, evolutionary algorithms ([15] for a comparison of such methods). Other classical estimators are obtained by interpreting noisy ODEs as state-space models: filtering and smoothing techniques can be used for parameter
algorithm, as introduced and motivated in section 5 of *Supplementary Material I*. We also compute the coverage probabilities of the estimator (NLS, OC and OC with an optimally weighted weight matrix): we look at the influence of the sample size *n* and of the noise level confidence ellipses derived from the asymptotic (Gaussian) approximation.

The supplementary material contains additional tables and figures for the models introduced in the paper: α-pinene, Riccati, and also for the real data example: Influenza virus growth and Blowfly populations dynamics. Moreover, we have tested our approach with the “FitzHugh-Nagumo model” for neuron dynamics. This model is chosen as it is an example of 2D nonlinear ODE with nonlinear dependence in the parameters (with a periodic solution). This model was introduced by Ramsay et al. as a benchmark for ODE estimation. We recall then briefly the model, and we compare the 3 approaches discussed in the paper. Concerning the Riccati equation, we plot the objective functions minimized by NLS and OC, in order to explain the difference of the NLS and OC estimators when the change-point time is not known.

Comments

From the simulations, we can check that the confidence sets given by NLS are too stringent, and NLS overestimates the coverage probabilities, as the (Monte-Carlo) estimated coverage probabilities are quite smaller than 95% (for all simulated models: α-pinene, Riccati, FitzHugh-Nagumo). On the contrary for OC, we can check the reliability of the derived confidence sets, as the estimated coverage probability is around 95% (or even higher). This is partly due to the estimated variance of OC that is higher than the NLS (as we can see by comparing the sum of estimated parameters variances $Tr(V(\hat{\theta}))$). Hence, this also confirms the quality of the Gaussian approximation derived from the asymptotic analysis of the paper.

The Optimally Weighted OC was introduced in order to ameliorate the variance of the OC esti-
mator, in order to reduce the size of the Confidence sets while preserving the coverage probabilities.
Although this version should improve on the unweighted OC estimator, this is not always true in
practice, and the “Optimal” estimator can be worse. This is mainly due to the fact that we approx-
imate the variance by linearization technics, and we do not control the error propagation of this
approximation in the IRWOC algorithm. Hence, we are not sure that we obtain the optimal matrix.
Moreover the control of the nonsingularity of the variance matrix can force us to reduce the number
of orthogonal conditions, which can induce a loss of inference power. Nevertheless, the optimal
weighting approaches improve upon the simple OC estimator for α-pinene and FitzHugh-Nagumo;
in addition the optimal weighting estimator gives the best estimate for the Influenza Virus Growth
model.
**Condition C6** There exists $L_0$ such that for $L \ge L_0$, $J_{\theta,L}(g, \theta)$ is full rank in a neighborhood of $(\phi^*, \theta^*)$.

**Condition C7 (a)** The times $T_1, \dots, T_n$ have a density $\pi$ w.r.t. Lebesgue measure such $0 < c < \pi < C < \infty$; 
(b) $E[\epsilon^4] < \infty$.

**Condition C8** For $\ell = 1, \dots, L$, $\theta \in \Theta$, there exists $\tilde{\beta}_{K_\ell}$ in $\mathbb{R}^{K_\ell}$ with $\| \frac{f_{x}(\cdot, \phi^{*}, \theta)\phi_{\ell} + \dot{\phi}_{\ell}}{\pi} - \tilde{\beta}_{K_\ell}^{\top} p^{K_\ell} \|_{L^2} \to 0$.

**Condition C9 (a)** The series estimator is a regression spline with a uniform knot sequence $(\tau_{1,K}, \dots, \tau_{N_K,K})$ defining the spline basis $\mathbf{p}^K$ satisfies $\max_i |\tau_{i+1,K} - \tau_{i,K}| \to 0$ as $K \to 0$; 
(b) For all $\theta \in \Theta$, for $\ell = 1 \dots L$, $v_\ell : t \mapsto \frac{f_{x}(t, \phi^{*}(t), \theta)\phi_{\ell}(t) + \dot{\phi}_{\ell}(t)}{\pi(t)}$ is $C^1$.

## 3 Consistency

**Theorem 3.1.** If conditions C1 to C6 are satisfied, then

$$ \hat{\theta}_{n,L} - \theta_{L}^{*} = O_{P}(1) $$

and the bias $B_L = \theta_L^* - \theta^*$ tends to zero as $L \to \infty$.

In particular, if we use the sine basis and if $\mathcal{E}(\phi^*, \theta)$ is in $H^1$ for all $\theta$, then $B_L = o(\frac{1}{L})$.

*Proof.* The classical proof for the consistency of an M-estimator $\hat{\theta}_{n,L} = \arg\min_{\theta \in \Theta} Q_{n,L}(\theta)$ such as the orthogonal conditions estimator relies basically on two stages: the uniform convergence of $Q_{n,L}(\theta)$ towards $Q_L^*(\theta) = \sum_{\ell=1}^L |e_\ell(\phi^*, \theta)|^2$ for $\theta$ in $\Theta$, and the fact that the true parameter $\hat{\theta}_{n,L}$ is a unique isolated global maximum (by theorem 5.7 in [6]). In a first step, we assume that $Q_L^*(\theta)$ has a unique global minimum denoted $\theta_L^*$, and we will show that this is indeed the case, and that $\theta_L^*$ is not very far from the true parameter $\theta^*$. We have to show that $\sup_\theta |Q_{n,L}(\theta) - Q_L^*(\theta)| \to 0$ as $n \to \infty$. From the simple additive expression of $Q_{n,L}(\theta)$, we see that it suffices to show the uniform convergence of $|e_\ell(\hat{\phi}, \theta)|^2$ to $|e_\ell(\phi^*, \theta)|^2$. From the inequality $|a^2 - b^2| \le |a - b|(|a| + |b|)$, we have

$$ \left| |e_{\ell}(\hat{\phi}, \theta)|^2 - |e_{\ell}(\phi^{*}, \theta)|^2 \right| \le |e_{\ell}(\hat{\phi}, \theta) - e_{\ell}(\phi^{*}, \theta)| \left( |e_{\ell}(\hat{\phi}, \theta)| + |e_{\ell}(\phi^{*}, \theta)| \right). $$

As $\hat{\phi}$ is such that for $j = 1, 2$ we have $\|\hat{\phi}_j - \phi_j^*\|_\infty = O_P(\zeta_0(K) (\sqrt{K/n} + K^{-\alpha}))$ this means that $\hat{\phi}(t) \in \mathcal{D} = \{x \in \mathbb{R}^2 | \exists t \in [0, 1], |x - \phi^*(t)| \le D\}$ with a probability tending to 1. This shows that with a probability tending to 1, we have $|f(t, \hat{\phi}(t), \theta)| \le h(t, \theta) \le H(t)$ (because $f$ is uniformly $L^2$-Caratheodory), moreover $\hat{\phi}$ is also bounded by $M + D > 0$ (because $\phi^*$ is bounded as a continuous function on $[0, 1]$) and

$$ \begin{aligned} |e_{\ell}(\hat{\phi}, \theta)| + |e_{\ell}(\phi^{*}, \theta)| &\le |\langle \mathcal{E}(\hat{\phi}, \theta), \varphi_{\ell} \rangle| + |\langle \hat{\phi}, \dot{\varphi}_{\ell} \rangle| + |\langle \mathcal{E}(\phi^{*}, \theta), \varphi_{\ell} \rangle| + |\langle \phi^{*}, \dot{\varphi}_{\ell} \rangle| \\ &\le 2 \|H\|_{L^2} + \|\dot{\varphi}_{\ell}\|_{L^2} (2M + D) \end{aligned} $$

Hence $|e_\ell(\hat{\phi}, \theta)| + |e_\ell(\phi^*, \theta)|$ is uniformly bounded in probability, so that the uniform convergence of the criterion boils down to the uniform convergence of $|e_\ell(\hat{\phi}, \theta) - e_\ell(\phi^*, \theta)|$. We can re-write $e_\ell(\hat{\phi}, \theta) - e_\ell(\phi^*, \theta) = \langle \mathcal{E}(\hat{\phi}, \theta) - \mathcal{E}(\phi^*, \theta), \varphi_\ell \rangle + \langle \hat{\phi} - \phi^*, \dot{\varphi}_\ell \rangle$ and it is clear that the second right-hand side term converges uniformly to zero in probability for all $\ell \ge 1$. Consequently, we just have to check that $\mathcal{E}(\hat{\phi}, \theta) - \mathcal{E}(\phi^*, \theta)$ converges uniformly in $\theta$ to 0 in probability. First of all, we remark that $g \mapsto \langle \mathcal{E}(g, \theta), \varphi_\ell \rangle$ is a continuous function from $B_\infty(\phi^*, D) = \{g \in C([0, 1]) \mid \|g - \phi^*\|_\infty \le D\}$ to $\mathbb{R}$ (w.r.t the sup-norm), because

$$ \begin{aligned} |\langle \mathcal{E}(g, \theta), \varphi_\ell \rangle - \langle \mathcal{E}(g', \theta), \varphi_\ell \rangle| &\le |a(\cdot, \theta)| |g - g'|, \varphi_\ell| \\ &\le |a(\cdot, \theta), \varphi_\ell| \|g - g'\| \end{aligned} $$

with $a(\cdot, \theta) \in L^2$ for all $\theta$, as the vector field is $L^2$-Lipschitz. By the continuous mapping theorem, we get the point-wise convergence of $\langle \mathcal{E}(\hat{\phi}, \theta) - \mathcal{E}(\phi^*, \theta), \varphi_\ell \rangle$, and the hard part consist in the uniform convergence. This
## 1 Introduction

### 1.1 Problem position and motivations

Differential Equations are a standard mathematical framework for modeling dynamics in physics, chemistry, biology, engineering sciences, etc and have proved their efficiency in describing the real world. Such models are defined thanks to a time-dependent vector field $f$, defined on the state-space $\mathcal{X} \subset \mathbb{R}^d$ and that depends on a parameter $\theta \in \Theta \subset \mathbb{R}^p$, $d, p \ge 1$. The vector field is then a function from $[0, 1] \times \mathcal{X} \times \Theta$ to $\mathbb{R}^d$. If $\phi(t)$ is the current state of the system, the time evolution is given by the following Ordinary Differential Equation, defined for $t \in [0, 1]$ by:

$$ \dot{\phi}(t) = f(t, \phi(t), \theta) \qquad (1.1) $$

where dot indicates derivative with respect to time. An important task is then the estimation of the parameter $\theta$ from real data. [30] proposed a significant improvement to this statistical problem, and gave motivations for further statistical studies. We are interested in the definition and in the optimality of a statistical procedure for the estimation of the parameter $\theta$ from noisy observations $y_1, \dots, y_n \in \mathbb{R}^d$ of a solution at times $t_1 < \dots < t_n$.

Most works deal with Initial Value Problems (IVP), i.e. with ODE models having a given (possibly unknown) initial value $\phi(0) = \phi_0$. There exists then a unique solution $\phi(\cdot, \phi_0, \theta)$ to the ODE (1.1) defined on the interval $[0, 1]$, that depends smoothly on $\phi_0$ and $\theta$.

The estimation of $\theta$ is a classical problem of nonlinear regression, where we regress $y$ on the time $t$.

<footnote n="1">1 Variational formulation, Sobolev Space.</footnote>

<footnote n="1">¹ ENSIIE & Laboratoire Statistique et Génome, Université d'Evry Val d'Essonne, UMR CNRS 8071 - USC INRA - FRANCE</footnote>

<footnote n="2">² Laboratoire Analyse et Probabilités, Université d'Evry Val d'Essonne - FRANCE</footnote>

<footnote n="3">³ Laboratoire IBISC, Université d'Evry Val d'Essonne - FRANCE</footnote>

<footnote n="4">⁴ INRIA-Saclay, LRI, Université Paris Sud, UMR CNRS 8623 - FRANCE</footnote>
inference [5], which can provide estimates with reduced computational complexity [20, 11, 10].

1

Nevertheless, the difficulty of the optimization problem is the outward sign of the illposedness of the
inverse problem of ODE parameter estimation, [6]. Hence some improvements on classical estimation
have been proposed by regularizing the statistical inference in an appropriate way. Starting from different
methods used for solving ODEs, different estimators can be developed based on a mixture of nonpara-
metric estimation and collocation approximation. This gives rise to Gradient Matching (or Two-Step)
estimators that consists in approximating the solution $\phi$ with a basis expansion, such as cubic splines.
The rationale is to estimate nonparametrically the solution $\phi$ by $\hat{\phi} = \sum_{k=1}^{L} \hat{c}_k B_k$ so that we can also
estimate the derivative $\hat{\dot{\phi}}$. An estimator of $\theta$ can be obtained by looking for the parameter that makes
$\hat{\phi}$ satisfy the differential equation (1.1) in the best possible manner. Two different methods have been
proposed, based on a $L^2$ distance between $\hat{\dot{\phi}}$ and $f(t, \hat{\phi}, \theta)$: The first one, called the two-step method, was
originally proposed by [25], and has been particularly developed in (bio)chemical engineering [13, 26, 18].
It avoids the numerical integration of the ODE and usually gives rise to simple optimization program
and fast procedures that usually performs well in practice. The statistical properties of this two stage
estimator (and several variants) have been studied in order to understand the influence of nonparametric
technics to estimate a finite dimensional parameter [4, 12, 8, 27]. While keeping the same kind of numer-
ical approximation of the solution, [21] proposed a second method based on the generalized smoothing
approach for determining at the same time the parameter $\theta$ and the nonparametric estimation $\hat{\phi}$. The
essential difference between these two approaches is that the nonparametric estimator in the generalized
smoothing approach is computed adaptively with respect to the parametric model, whereas two-step
estimators are "model-free smoothing".

We introduce here a new estimator that can be seen as an improvement of the previous two-step
estimators. It uses also a nonparametric proxy $\hat{\phi}$, but we modify the criterion used to identify the
ODE parameter (i.e. the second step). The initial motivation is to get a closed-form expression for
the asymptotic variance and confidence sets. The most notable feature of the proposed method is the
use of a variational characterization of the solution of the differential equations instead of the classical
point-wise one. This characterization is general and can cover a greater number of situations. Thanks to
its computational tractability, we can give a precise description of the asymptotics and give the bias and
# 1 Linear ODE: α-pinene

## 1.1 Known initial condition

<table><caption>Table 1: MSE, Asymptotic Variance & ARE for α-pinene model with known Initial Condition</caption><thead><tr><th>×10<sup>-2</sup></th><th colspan="6">MSE</th><th colspan="5">Tr (V(θ̂))</th></tr><tr><th>(n, σ)</th><th>TS</th><th>OC</th><th>OCopt</th><th>OC,0</th><th>OC,0 opt</th><th>NLS</th><th>OC</th><th>OCopt</th><th>OC,0</th><th>OC,0 opt</th><th>NLS</th></tr></thead><tbody><tr><td>(400, 3)</td><td>0.72</td><td>0.05</td><td>0.04</td><td>0.04</td><td>0.03</td><td>0.02</td><td>0.04</td><td>0.04</td><td>0.04</td><td>0.05</td><td>0.02</td></tr><tr><td>(400, 8)</td><td>2.28</td><td>0.22</td><td>0.27</td><td>0.25</td><td>0.22</td><td>0.10</td><td>0.95</td><td>1.91</td><td>1.20</td><td>2.23</td><td>0.12</td></tr><tr><td>(200, 3)</td><td>1.19</td><td>0.27</td><td>0.33</td><td>0.30</td><td>0.32</td><td>0.03</td><td>0.09</td><td>0.20</td><td>0.13</td><td>0.24</td><td>0.03</td></tr><tr><td>(200, 8)</td><td>2.95</td><td>0.44</td><td>0.51</td><td>0.37</td><td>0.47</td><td>0.18</td><td>2.66</td><td>4.96</td><td>2.68</td><td>5.48</td><td>0.27</td></tr><tr><td>(50, 3)</td><td>2.39</td><td>0.27</td><td>0.27</td><td>0.26</td><td>0.26</td><td>0.16</td><td>1.37</td><td>2.75</td><td>1.58</td><td>3.31</td><td>0.16</td></tr><tr><td>(50, 8)</td><td>4.54</td><td>1.03</td><td>0.94</td><td>0.93</td><td>0.89</td><td>0.68</td><td>7.96</td><td>8.20</td><td>7.27</td><td>9.13</td><td>1.68</td></tr></tbody></table><table><thead><tr><th>×10<sup>-2</sup></th><th colspan="6">ARE</th></tr><tr><th>(n, σ)</th><th>TS</th><th>OC</th><th>OCopt</th><th>OC,0</th><th>OC,0 opt</th><th>NLS</th></tr></thead><tbody><tr><td>(400, 3)</td><td>105.85</td><td>20.45</td><td>19.83</td><td>18.99</td><td>17.64</td><td>11.61</td></tr><tr><td>(400, 8)</td><td>213.69</td><td>47.90</td><td>49.48</td><td>50.77</td><td>45.29</td><td>28.65</td></tr><tr><td>(200, 3)</td><td>150.53</td><td>37.16</td><td>40.36</td><td>36.38</td><td>37.40</td><td>16.58</td></tr><tr><td>(200, 8)</td><td>235.62</td><td>57.50</td><td>60.64</td><td>58.54</td><td>59.97</td><td>38.93</td></tr><tr><td>(50, 3)</td><td>220.59</td><td>51.19</td><td>52.41</td><td>51.13</td><td>51.72</td><td>34.49</td></tr><tr><td>(50, 8)</td><td>283.02</td><td>112.52</td><td>110.44</td><td>107.21</td><td>107.77</td><td>76.52</td></tr></tbody></table>

<figure><img src="image_1.png" /><figcaption>Figure 1: Coverage Probabilities for the 95% Confidence Ellipse for α-pinene model with known Initial Condition</figcaption></figure>
can be proven by controlling the oscillations in $\theta$ of the process $\mathcal{E} = (\mathcal{E}(\theta))_{\theta \in \Theta} = (\langle \mathcal{E}(\hat{\phi}, \theta) - \mathcal{E}(\phi^*, \theta), \varphi_\ell \rangle)_{\theta \in \Theta}$ and by using theorem 18.11 in [6], as convergence in probability towards a constant is equivalent to weak convergence. In order to show that the process $\mathcal{E}$ converges weakly to 0 in the space $C(\Theta)$ of continuous functions on $\Theta$ equipped with the supremum norm, we have to check that for all $k$ and all $\theta_1, \dots, \theta_k$ in $\Theta$

$$ (\mathcal{E}(\theta_1), \dots, \mathcal{E}(\theta_k)) \rightsquigarrow (0, \dots, 0) \qquad (2) $$

and that for all $\epsilon, \alpha > 0$ there exists a partition of $\Theta_1, \dots, \Theta_K$ of $\Theta$ such that

$$ \limsup_{n \to \infty} P \left( \sup_k \sup_{\theta, \theta' \in \Theta_k} |\mathcal{E}(\theta) - \mathcal{E}(\theta')| \ge \alpha \right) \le \epsilon. \qquad (3) $$

The first condition (2) is a direct consequence of the point-wise convergence in probability of $\mathcal{E}(\theta)$ in $\Theta$. Concerning the second condition, from we have for $t$ in $[0, 1]$ a.e.

$$ \begin{aligned} \left| f(t, \hat{\phi}(t), \theta) - f(t, \hat{\phi}(t), \theta') \right| & \le \sum_{i=1}^{p} \left| f_{\theta_i}(t, \hat{\phi}(t), \tilde{\theta}_n(t)) \right| \left| \theta_i - \theta_i' \right| \\ & \le \bar{a}'(t) \left| \theta - \theta' \right| \end{aligned} $$

with a probability tending to 1 ($\tilde{\theta}_n(t)$ being a parameter between $\theta$ and $\theta'$). As a consequence the following inequality

$$ |\mathcal{E}(\theta) - \mathcal{E}(\theta')| \le \langle \bar{a}' |\theta - \theta'|, |\varphi_\ell| \rangle \le |\theta - \theta'| $$

is true with a probability tending to 1. Since $\Theta$ is a compact set, it is possible to find a finite partition $\Theta_1, \dots, \Theta_K$ of $\Theta$ such that the diameter of $\Theta_i$ is smaller than an arbitrary $\alpha$ independently of $n$. This ensures that condition (3) is also satisfied and the uniform convergence of $Q_{n,L}(\theta)$ to $Q_L^*(\theta)$ can be claimed. Now, we relate $\theta_L^* - \theta^*$ to the approximation quality of $\mathcal{E}_F(\phi^*, \theta)$ by the basis $(\phi_\ell)_{\ell \ge 1}$. We remark first that we have $|\mathcal{E}_F(\phi^*, \theta)| \le 2\bar{h}(t)$: this implies that for all $\ell$, and all $\theta$ we have $|e_\ell^*(\theta)| \le 2\bar{h}_\ell$ with $\bar{h}_\ell = \langle \bar{h}, \varphi_\ell \rangle$. The global rate of convergence of the series $\sum_\ell e_\ell^*(\theta)$ is controlled uniformly by the rate of $\sum_{\ell \ge 1} H_\ell^2$, denoted $r_L^2 = \sum_{\ell > L} H_\ell^2$. By orthogonality, we can write

$$ \| \mathcal{E}_F(\phi^*, \theta) \|_{L^2}^2 = \| \mathcal{E}_L(\phi^*, \theta) \|_{L^2}^2 + \| \mathcal{R}_L(\phi^*, \theta) \|_{L^2}^2 $$

which means that is $\|\mathcal{E}_L(\phi^*, \theta)\|_{L^2}^2$ is a perturbation of the function $\theta \mapsto \|\mathcal{E}_F(\phi^*, \theta)\|_{L^2}^2$ by the function $\theta \mapsto -\|\mathcal{R}_L(\phi^*, \theta)\|_{L^2}^2$. This perturbation $\|\mathcal{R}_L(\phi^*, \theta)\|_{L^2}^2$ is uniformly dominated by $r_L$, hence it becomes possible to relate the two minima. From assumption C5, we know $\|\mathcal{R}_L(\phi^*, \theta)\|_{L^2}^2$ is differentiable and we compute a series decomposition of its gradient thanks to $\partial_\theta \mathcal{R}_L(\phi^*, \theta) = \sum_{\ell > L} \partial_\theta e_\ell(\phi^*, \theta) \varphi_\ell$:

$$ \partial_\theta \| \mathcal{R}_L(\phi^*, \theta) \|_{L^2}^2 = 2 \mathcal{R}_L(\phi^*, \theta) \partial_\theta \mathcal{R}_L(\phi^*, \theta). $$

and we recall that $\mathcal{R}_L(\phi^*, \theta)$ converges uniformly to 0 and $\partial_\theta \mathcal{R}_L(\phi^*, \theta)$ is uniformly bounded in $\theta$ on $\Theta$ (as a continuous function on the compact set $\Theta$). Starting from this last remark, we use the Implicit Function Theorem to the continuously differentiable function $G(\epsilon, \theta) = \partial_\theta \|\mathcal{E}_F(\phi^*, \theta)\|_{L^2}^2 - 2\epsilon \partial_\theta \mathcal{R}_L(\phi^*, \theta)$, in order to get a Taylor expansion. We denote $\theta_\epsilon$ the solution to $G(\epsilon, \theta_\epsilon)$, and we remark that in particular $G(0, \theta^*) = 0$. Thus, there exists $\epsilon_0, \delta_0 > 0$ and a function $\psi: [-\epsilon_0, \epsilon_0] \to B(\theta^*, \delta_0)$ such that $\psi(\epsilon) = \theta_\epsilon$ (i.e. $G(\epsilon, \psi(\epsilon)) = 0$). We can also compute the first order variation of $\psi: \psi(\epsilon) = \psi(0) + \epsilon\psi'(0) + o(\epsilon)$ where

$$ \psi'(0) = -2 \left( \partial_{\theta\theta} \| \mathcal{E}_F(\phi^*, \theta) \|_{L^2}^2 \right)^{-1} \partial_\theta \mathcal{R}_L(\phi^*, \theta^*). $$

Since $\mathcal{R}_L(\phi^*, \theta)$ converges uniformly to 0, there exists $L_0 > 0$, such that for $L > L_0$, $|\mathcal{R}_L(\phi^*, \theta)| \le \epsilon_0$ so that we can apply the linearization above to $G(\mathcal{R}_L(\phi^*, \theta), \theta) = \partial_\theta \|\mathcal{E}_F(\phi^*, \theta)\|_{L^2}^2 - \partial_\theta \|\mathcal{R}_L(\phi^*, \theta)\|_{L^2}^2$. We obtain that the minima $\theta_L^*$ and $\theta^*$ are such that

$$ \theta_L^* = \theta^* + \mathcal{R}_L(\phi^*, \theta) 2 \left( \partial_{\theta\theta} \| \mathcal{E}_F(\phi^*, \theta) \|_{L^2}^2 \right)^{-1} \partial_\theta \mathcal{R}_L(\phi^*, \theta^*) + o(r_L). $$

This implies that $|\theta_L^* - \theta^*| = O(r_L)$. If we use the sine basis, this means that $\mathbf{P}_F H = \sum_\ell \bar{h}_\ell \sqrt{2} \sin(\pi \ell t)$ but it is also in $H_0^1$ (as it is in $H^1$), then we have $\sum_{\ell \ge 1} \ell^2 \bar{h}_\ell^2$ and $r_L^2 = o(\frac{1}{L^2})$. $\square$
1. If $\phi_0$ is known, the Nonlinear Least Square Estimator $\hat{\theta}^{NLS}$ (NLSE) is obtained by minimizing

$$Q_n^{LS}(\theta) = \sum_{i=1}^{n} |y_i - \phi(t_i, \phi_0, \theta)|^2 \qquad (1.2)$$

2. where $|\cdot|$ is the classical Euclidean norm. The NLSE, Maximum Likelihood Estimator or more general 3. M-estimators [36] are commonly used because of their good statistical properties (root-$n$ consistency, 4. asymptotic efficiency), but they come with important computational difficulties (repeated ODE inte- 5. grations and presence of multiple local minima) that can decrease their interest. We refer to [30] for a 6. detailed overview of the previous works in this field. An adapted NLS estimator (dedicated the specific 7. difficulties of ODEs) is also introduced and studied in [43].

8. Global optimization methods are then often used, such as simulated annealing, evolutionary algo- 9. rithms ([22] for a comparison of such methods). Other classical estimators are obtained by interpreting 10. noisy ODEs as state-space models: filtering and smoothing techniques can be used for parameter infer- 11. ence [9], which can provide estimates with reduced computational complexity [29, 17, 16].

12. Nevertheless, the difficulty of the optimization problem is the outward sign of the ill-posedness of the 13. inverse problem of ODE parameter estimation, [12]. Hence some improvements on classical estimation 14. have been proposed by adding regularization constraints in an appropriate way.

15. Starting from different methods used for solving ODEs, different estimators can be developed based 16. on a mixture of nonparametric estimation and collocation approximation. This gives rise to Gradient 17. Matching (or Two-Step) estimators that consists in approximating the solution $\phi$ with a basis expansion 18. $\{B_1, \dots, B_K\}$, such as cubic splines. The rationale is to estimate nonparametrically the solution $\phi$ by 19. $\hat{\phi} = \sum_{k=1}^{L} \hat{c}_k B_k$ so that we can also estimate the derivative $\hat{\dot{\phi}}$. An estimator of $\theta$ can be obtained by 20. looking for the parameter that makes $\hat{\phi}$ satisfy the differential equation (1.1) in the best possible manner. 21. Two different methods have been proposed, based on a $L^2$ distance between $\hat{\dot{\phi}}$ and $f(t, \hat{\phi}, \theta)$: The first 22. one, called the two-step method, was originally proposed by [38], and has been particularly developed 23. in (bio)chemical engineering [20, 39, 28]. It avoids the numerical integration of the ODE and usually 24. gives rise to simple optimization program and fast procedures that usually performs well in practice. 25. The statistical properties of this two stage estimator (and several variants) have been studied in order 26. to understand the influence of nonparametric techniques to estimate a finite dimensional parameter
1 variance of the estimator. We also give a way to ameliorate adaptively our estimator and to compute
2 asymptotic confidence intervals.

3 First, we give the ODE model and main assumptions, we introduce our estimator, and show its
4 consistency. Then, we provide a detailed description of the asymptotics, by proving its root-n consistency
5 and asymptotic normality. Based on the asymptotic approximation, we give a closed-form expression of
6 the asymptotic variance, and we address the problem of obtaining the best variance through the choice
7 of an appropriate weighting matrix. Finally, we provide some insights into the practical behavior of the
8 estimator through simulations. The objective of the experiments part is to show its interest with respect
9 to the nonlinear least squares and classical gradient matching estimators.

## 2 Differential Equation Model and Gradient Matching

### 2.1 ODE models and Gradient Matching

12 For ease of readability, we focus on a two-dimensional system of ODEs. In our case, as there is no
13 computational and theoretical differences between the situation $d = 2$ and $d > 2$, there is no lack
14 of generality by this assumption. We consider noisy observations $Y_1, \dots, Y_n \in \mathbb{R}^2$ of the function $\phi^*$
15 measured at random times $t_1 < \dots < t_n \in [0, 1]$:

$$Y_i = \phi^*(t_i) + \epsilon_i \qquad (2.1)$$

16 where $\epsilon_1, \dots, \epsilon_n$ are i.i.d with $E(\epsilon_i) = 0$ and $V(\epsilon_i) = \sigma^2 I_2$. We suppose that the regression function
17 $\phi^*$ belongs to the Sobolev space $H^1 = \{u \in L^2([0, 1]) \mid \dot{u} \in L^2([0, 1])\}$, and $\phi^*$ is a (Caratheodory or
18 generalized) solution to the parametrized Ordinary Differential Equation (1.1), i.e. there exists a true
19 parameter $\theta^* \in \Theta \subset \mathbb{R}^p$ such that for $t \in [0, 1]$ almost everywhere (a.e.)

$$\dot{\phi}^*(t) = \mathbf{f}(t, \phi^*(t), \theta^*) \qquad (2.2)$$

20 where $\mathbf{f} = (f_1, f_2)$ is a vector field from $[0, 1] \times \mathcal{X} \times \Theta$ to $\mathbb{R}^2$, where $\mathcal{X} \subset \mathbb{R}^2$. The classical Cauchy-
21 Lipschitz theory of Initial Value Problem (IVP) can be extended to rougher vector field in time $t$ (see
## 1.2 Unknown initial condition

<table><caption>Table 2: MSE, Asymptotic Variance & ARE for α-pinene model with unknown initial conditions</caption><thead><tr><th>×10<sup>-2</sup></th><th colspan="6">MSE</th><th colspan="5">Tr (V(θ))</th></tr><tr><th>(n, σ)</th><th>TS</th><th>OC</th><th>OCopt</th><th>OC,1</th><th>OC,1 opt</th><th>NLS</th><th>OC</th><th>OCopt</th><th>OC,1</th><th>OC,1 opt</th><th>NLS</th></tr></thead><tbody><tr><td>(400, 3)</td><td>0.25</td><td>0.11</td><td>0.25</td><td>0.11</td><td>0.52</td><td>0.07</td><td>0.10</td><td>0.11</td><td>0.10</td><td>0.09</td><td>0.06</td></tr><tr><td>(400, 8)</td><td>1.07</td><td>0.85</td><td>0.85</td><td>0.56</td><td>0.73</td><td>0.50</td><td>1.06</td><td>0.94</td><td>0.82</td><td>0.45</td><td>0.61</td></tr><tr><td>(200, 3)</td><td>0.6</td><td>0.37</td><td>0.50</td><td>0.23</td><td>0.34</td><td>0.14</td><td>0.25</td><td>0.25</td><td>0.20</td><td>0.17</td><td>0.14</td></tr><tr><td>(200, 8)</td><td>1.64</td><td>1.42</td><td>1.18</td><td>0.83</td><td>1.68</td><td>1.34</td><td>2.36</td><td>2.30</td><td>1.64</td><td>0.50</td><td>1.54</td></tr><tr><td>(50, 3)</td><td>1.33</td><td>1.31</td><td>1.18</td><td>0.80</td><td>1.04</td><td>0.69</td><td>1.63</td><td>1.46</td><td>1.02</td><td>0.39</td><td>0.76</td></tr><tr><td>(50, 8)</td><td>3.64</td><td>2.11</td><td>1.91</td><td>1.79</td><td>3.44</td><td>1.96</td><td>5.34</td><td>4.35</td><td>2.20</td><td>0.81</td><td>4.38</td></tr></tbody></table><table><thead><tr><th>×10<sup>-2</sup></th><th colspan="6">ARE</th></tr><tr><th>(n, σ)</th><th>TS</th><th>OC</th><th>OCopt</th><th>OC,1</th><th>OC,1 opt</th><th>NLS</th></tr></thead><tbody><tr><td>(400, 3)</td><td>60.13</td><td>32.66</td><td>39.51</td><td>32.58</td><td>56.48</td><td>24.13</td></tr><tr><td>(400, 8)</td><td>127.82</td><td>86.64</td><td>88.53</td><td>79.02</td><td>112.08</td><td>63.52</td></tr><tr><td>(200, 3)</td><td>88.75</td><td>52.84</td><td>59.02</td><td>47.48</td><td>75.10</td><td>34.70</td></tr><tr><td>(200, 8)</td><td>158.17</td><td>117.44</td><td>113.69</td><td>108.70</td><td>177.77</td><td>98.85</td></tr><tr><td>(50, 3)</td><td>138.31</td><td>99.87</td><td>97.36</td><td>89.22</td><td>121.73</td><td>69.27</td></tr><tr><td>(50, 8)</td><td>247.57</td><td>161.97</td><td>164.31</td><td>165.55</td><td>214.20</td><td>144.83</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>Figure 2: Coverage Probabilities for the 95% Confidence Ellipse for α-Pinene with unknown Initial Conditions</figcaption></figure>
# 4 Asymptotics

The asymptotics of the estimators $\hat{\theta}_{n,L}$ are obtained by two (successive) linearizations of the orthogonal conditions $e_L(\hat{\phi}, \hat{\theta}_{n,L})$ around $(\phi^*, \theta^*)$. We detail the proof of proposition 4.1 and theorem 4.1.

## Section 4.1: Proofs

**Proposition 4.1.** If conditions C1-C6 are satisfied, then

$$ \left[ \mathbf{J}_{\theta,L} (\hat{\phi}, \hat{\theta}_{n,L})^{\top} \tilde{\mathbf{J}}_{\theta,L} \right]^{-1} \mathbf{J}_{\theta,L} (\hat{\phi}, \hat{\theta}_{n,L})^{\top} \xrightarrow{P} \mathbf{M}_L^* = \left[ \mathbf{J}_{\theta,L}^{*\top} \mathbf{J}_{\theta,L}^{*} \right]^{-1} \mathbf{J}_{\theta,L}^{*\top} \quad (4) $$

where the matrix $\tilde{\mathbf{J}}_{\theta,L}$ is the Jacobian $\mathbf{J}_{\theta,L}$ evaluated at a point $\tilde{\theta}$ between $\theta^*$ and $\hat{\theta}_{n,L}$. Moreover, we have

$$ \hat{\theta}_{n,L} - \theta_L^* = -\mathbf{M}_L^* \mathbf{e}_L(\hat{\phi}, \theta_L^*) + o_P(1). \quad (5) $$

*Proof.* If conditions C1-C6 are satisfied, then

$$ \left[ \mathbf{J}_{\theta,L} (\hat{\phi}, \hat{\theta}_{n,L})^{\top} \tilde{\mathbf{J}}_{\theta,L} \right]^{-1} \mathbf{J}_{\theta,L} (\hat{\phi}, \hat{\theta}_{n,L})^{\top} \xrightarrow{P} \mathbf{M}_L^* = \left[ \mathbf{J}_{\theta,L}^{*\top} \mathbf{J}_{\theta,L}^{*} \right]^{-1} \mathbf{J}_{\theta,L}^{*\top} \quad (6) $$

and we have

$$ \hat{\theta}_{n,L} - \theta_L^* = -\mathbf{M}_L^* \mathbf{e}_L(\hat{\phi}, \theta^*) + o_P(1). \quad (7) $$

The first order condition implies that $\hat{\theta}_{n,L}$ satisfies

$$ \mathbf{J}_{\theta,L} (\hat{\phi}, \hat{\theta}_{n,L})^{\top} \mathbf{e}_L (\hat{\phi}, \hat{\theta}_{n,L}) = 0. \quad (8) $$

We develop a Taylor expansion of $f(t, \hat{\phi}(t), \theta)$ around $\theta^*$ of order 1 for $t$ a.e. in $[0, 1]$: the orthogonal condition $e_\ell(\hat{\phi}, \theta)$ can be decomposed as $e_\ell(\hat{\phi}, \theta^*) + \langle f_\theta(\cdot, \hat{\phi}, \tilde{\theta}), \varphi_\ell \rangle$ where $\tilde{\theta}$ is on straight line between $\theta^*$ and $\hat{\theta}_{n,L}$. We can write it in vector form

$$ \mathbf{e}_L(\hat{\phi}, \hat{\theta}_{n,L}) = \mathbf{e}_L(\hat{\phi}, \theta^*) + \tilde{\mathbf{J}}_{\theta,L}(\hat{\theta}_{n,L} - \theta_L^*). \quad (9) $$

$\tilde{\mathbf{J}}_{\theta,L}$ is a matrix $\mathbb{R}^{L \times p}$ with entries $\langle f_{\theta_i}(\cdot, \hat{\phi}, \tilde{\theta}), \varphi_\ell \rangle$. Thus, if we premultiply (9) by $\mathbf{J}_{\theta,L}(\hat{\phi}, \hat{\theta}_{n,L})^\top$, we get the asymptotic expansion

$$ 0 = \mathbf{J}_{\theta,L} (\hat{\phi}, \hat{\theta}_{n,L})^{\top} \mathbf{e}_L(\hat{\phi}, \theta_L^*) + \mathbf{J}_{\theta,L} (\hat{\phi}, \hat{\theta}_{n,L})^{\top} \tilde{\mathbf{J}}_{\theta,L} (\hat{\theta}_{n,L} - \theta_L^*). \quad (10) $$

This shows that the key result for relating the behavior of $(\hat{\theta}_{n,L} - \theta_L^*)$ to the behavior of $\mathbf{e}_L(\hat{\phi}, \theta^*)$ is the convergence in probability of $\mathbf{J}_{\theta,L}(\hat{\phi}, \hat{\theta}_{n,L})^\top$. This is indeed the case because the matrix value function $(g, \theta) \mapsto J_{\theta,L}(g, \theta)$ is continuous w.r.t the uniform norm on $\mathcal{D} \times \Theta$ (because functions $\frac{\partial}{\partial \theta_i} f$ are uniformly Lipschitz) and $(\hat{\phi}, \hat{\theta}_{n,L})$ converges in probability to $(\phi^*, \theta^*)$. Finally, since the matrix $\mathbf{J}_{\theta,L}^{*\top} \mathbf{J}_{\theta,L}^{*}$ is nonsingular, by the continuous mapping theorem, we see that

$$ \left[ \mathbf{J}_{\theta,L} (\hat{\phi}, \hat{\theta}_{n,L})^{\top} \tilde{\mathbf{J}}_{L} \right]^{-1} \mathbf{J}_{n,L} (\hat{\theta}_{n,L})^{\top} $$

converges in probability to $[\mathbf{J}_{\theta,L}^* \mathbf{J}_{\theta,L}^*]^{-1} \mathbf{J}_{\theta,L}^{*\top} = \mathbf{M}_L^*$. Thus, this implies that $\hat{\theta}_{n,L} - \theta_L^* = -\mathbf{M}_L^* \mathbf{e}_{n,L}(\theta^*) + o_P(1)$. $\square$
[8, 19, 14]. While keeping the same kind of numerical approximation of the solution, [30] proposed a second method based on the generalized smoothing approach for determining at the same time the parameter $\theta$ and the nonparametric estimation $\hat{\phi}$. The essential difference between these two approaches is that the nonparametric estimator in the generalized smoothing approach is computed adaptively with respect to the parametric model, whereas two-step estimators are “model-free smoothing”.

We introduce here a new estimator that can be seen as an improvement and a generalization of the previous two-step estimators. It uses also a nonparametric proxy $\hat{\phi}$, but we modify the criterion used to identify the ODE parameter (i.e. the second step). The initial motivations are

• to get a closed-form expression for the asymptotic variance and confidence sets,

• to reduce sensitivity to the estimation of the derivative in Gradient Matching approaches,

• to take into account explicitly time-dependent vector field, with potential discontinuities in time.

The most notable feature of the proposed method is the use of a variational formulation of the differential equations instead of the classical point-wise one, in order to generate conditions to satisfy. This formulation is rather general and can cover a greater number of situations: we come up with a generic class of estimator of Differential Equations (e.g Ordinary, Delay, Partial, Differential-Algebraic), that can incorporate relatively easily prior knowledge about the true solution. In addition to the versatility of the method, the criterion is built in order to offer computational tractability, that implies that we can give a precise description of the asymptotics and give the bias and variance of the estimator. We also give a way to ameliorate adaptively our estimator and to compute asymptotic confidence intervals.

First, we introduce the statistical ODE-based model and main assumptions, we motivate and describe our estimator, and show its consistency. Then, we provide a detailed description of the asymptotics, by proving its root-n consistency and asymptotic normality. Based on the asymptotic approximation, we give a closed-form expression of the asymptotic variance, and we address the problem of obtaining the best variance through the choice of an appropriate weighting matrix. Finally, we provide some insights into the practical behavior of the estimator through simulations and by considering two real-data examples. The objective of the experiments parts is to show the interest of OC with respect to the nonlinear least squares and classical gradient matching estimators.
1 theorem 3.4 in [17]) under the following assumptions

2 IVP (a) $f$ is $L^2$-Caratheodory, i.e. $x \mapsto f(t, x, \theta)$ is continuous for $t$ a.e in $[0, 1]$, $x \mapsto f(t, x, \theta)$ is measurable for all $x \in \mathcal{X}$ and $\forall c > 0, \exists h_c(\cdot, \theta) \in L^2 \setminus |x| \le c \implies |f(t, x, \theta)| \le h_c(t, \theta)$,

4 IVP (b) $f$ is $L^2$-Lipschitz, i.e. $\exists a(\cdot, \theta) \in L^2 / \forall x, x' \in \mathcal{X}, |f(t, x, \theta) - f(t, x', \theta)| \le a(t, \theta) |x - x'|$.

5 In that case, for each $\theta \in \Theta$ and $x_0 \in \mathcal{X}$, we can guarantee that there exists a unique solution in $H^1$ to the IVP

$$ \left\{ \begin{array}{l} \dot{x}(t) = f(t, x(t), \theta) \\ x(0) = x_0 \end{array} \right. \qquad (2.3) $$

7 As a generalized solution is in $H^1$ and satisfies the ODE only almost everywhere, it motivates the use of weak derivative, that give a neat way to express the notion of a.e solution. Moreover, weak derivatives have been proved to be a very convenient way to solve differential equation models as it permits to introduce quite naturally variational formulation, in particular for Boundary Value Problems and Partial Differential Equations (a typical example is Galerkin approximation, [3]). If we denote the inner product of $L^2$ as $\forall \varphi, \psi \in L^2([0, 1])$, $\langle \varphi, \psi \rangle = \int_0^1 \varphi(t)\psi(t)dt$, the weak derivative of the function $g$ in $H^1$ is not defined point-wise but as the function $\dot{g} \in L^2$ satisfying $\langle \dot{g}, \psi \rangle = -\langle g, \dot{\psi} \rangle$, for all function $\varphi$ in $C^1$ with support included in $]0, 1[$ (denoted $C_C^1([0, 1])$). Of course, if $t \mapsto \phi(t, x_0, \theta)$ is a $C^1$ function on $]0, 1[$, the classical derivative $\dot{\phi}$ is also the weak derivative.

16 We introduce then the weak form or variational formulation of the ODE (2.3). A weak solution $g$ to (2.3) is a function in $H^1$ such that $\forall \varphi \in C_C^1([0, 1])$

$$ \int_0^1 \mathbf{f}(t, g(t), \theta)\varphi(t)dt + \int_0^1 g(t)\dot{\varphi}(t)dt = 0 \qquad (2.4) $$

18 In our case, a generalized solution is a weak solution and we know that there exists at least one weak solution to (2.4). For the statistical inference task, we use only the variational formulation as a necessary condition.
# 2 Nonlinear ODEs

## 2.1 Ricatti equation

We recall that the asymptotic criterion for the nonlinear least squares is

$$ 
\begin{aligned}
R_{NLS}^*(\theta, x_0) &= E_{T,Y} \left[ \|Y - \phi(T; \theta, x_0)\|^2 \right] \\
&= \sigma^2 + \int_0^1 \|\phi(t; \theta, x_0) - \phi(t; \theta^*, x_0^*)\|^2 \pi(t) dt
\end{aligned}
$$

whereas the asymptotic criterion for OC is

$$ Q_L^*(\theta) = \sum_{\ell=1}^{L} e_{\ell}(\phi^*, \theta)^2. $$

When the Ricatti equation is simply $\dot{x} = ax^2 + c\sqrt{t}$, with unknown a and c, the criterion $R_{NLS}$ is nonlinear in $\theta$ (no closed-form), whereas the $Q_L^*(\theta)$ is a simple quadratic form, which is simple to optimize. In that case, we can check (by Monte-Carlo simulations) that the optimization of the empirical counterpart of $R_{NLS}(\theta, x_0)$ gives better results (when we start close to $\theta^*$) than $Q_L^*(\theta)$. Nevertheless, when the parametrized ODE is $\dot{x} = ax^2 + c\sqrt{t} - d'c1_{[T_r;14]}$, the statistical inverse problem becomes much harder to solve. This can be seen by plotting the criterion function $R_{NLS}^*(\theta, x_0)$ which is a rough function, that makes it difficult to optimize. At the contrary, the function $Q_L^*(\theta)$ is quite smooth and exhibits (locally) only one minimum at $\theta^*$. In order to visualize the 2 objective functions, we fix $a = a^*$ and $c = c^*$ and we consider $d' \in [0, 3]$ and $T_r \in [10, 14]$ (for NLS, we have also to fix $x_0 = x_0^* = -1$). We recall that we have $d'^* = 2$ and $T_r^* = 11$.

### 2.1.1 Known change-point time $T_r$
## section 4.3: Proofs

For ease of reading, we restate theorem 4.1, and we give a sketch of the proof that relies essentially on the properties of plug-in estimators for series estimators obtained in [3]. First, we give a preliminary result about the asymptotics of the Orthogonal Conditions $e_L(\hat{\phi}, \theta) - e_L(\phi^*, \theta)$ for all $\theta \in \Theta$.

**Proposition 4.2 - Asymptotic Normality of Orthogonal Conditions for general series estimators.**
Under conditions C1-C3, C7 and C8 and if K is taken such that $\sqrt{n}K^{-\alpha} \to 0$ and

1. $f$ is linear (in state $x$),

2. OR, $f$ is nonlinear (in state $x$) and $\frac{\zeta(K)^4K^2}{n} \to 0$,

then for all $\theta \in \Theta$,

$$ \sqrt{n} (e_L(\hat{\phi}, \theta) - e_L(\phi^*, \theta)) \rightsquigarrow N(0, V_{e,L}(\theta)) $$

Moreover, the asymptotic variance can be estimated by $\hat{V}_{e,L}(\theta) = A(\phi^*, \theta)\hat{V}(\hat{c}_1)A(\phi^*, \theta)^\top + B(\phi^*, \theta)\hat{V}(\hat{c}_2)B(\phi^*, \theta)^\top$.

*Proof.* We have to show a kind of “delta method” result. Hence we need to compute the Fréchet derivative of the functional $g \mapsto e_L(g, \theta)$ around $\phi^*$ for each $\theta$. For simplicity, we consider the scalar functional $e_\ell(g, \theta)$. From section 4.3, we know that the Fréchet derivative for the sup norm in $\mathcal{D} = \{x \in \mathbb{R}^2 | \exists t \in [0, 1], |x - \phi^*(t)| < D\}$ is given, for all $h \in L^2$

$$ De_\ell(g, \theta).h = \int_0^1 \{f_x(t, g(t), \theta) \varphi(t) + \dot{\varphi}(t)\} h(t) dt. $$

Hence, the Fréchet derivative can be represented as a scalar product with a gradient vector $v(g)$, w.r.t. the $L^2$inner product with respect to the distribution $\pi$ of times $T$ (bounded away from 0, see condition C7), i.e

$$ \begin{aligned} De_\ell(g, \theta).h &= \int_0^1 \frac{\{f_x(t, g(t), \theta) \varphi(t) + \dot{\varphi}(t)\}}{\pi(t)} h(t) \pi(t) dt \\ &= \langle v(g), h \rangle_{L^2(\pi)} \end{aligned} $$

with $v(g)(t) = \frac{\{f_x(t, g(t), \theta)\varphi(t) + \dot{\varphi}(t)\}}{\pi(t)}$. Moreover, this gradient is in $L^2(\pi)$, because $f_x$ is dominated, $\varphi, \dot{\varphi}$ and $\pi^{-1}$ are bounded above.

It is possible to derive a central limit theorem if the gradient function $v(\phi^*)$ can be itself approximated by a series expansion in the same basis $p_{kK}$, $k = 1, \dots, K$ as the function $\phi^*$. We need to be sure that the approximation is still valid under a random design, which is immediate, as $\pi$ is equivalent to Lebesgue measure $c \le \pi \le C$. Indeed, let $\tilde{v}_K = \beta_K^\top p_K$ being an approximation of $v$ (according to condition C8), then we have

$$ \int_0^T (v(t) - \tilde{v}_K(t))^2 \pi(t) dt \le C \int_0^T (v(t) - \tilde{v}_K(t))^2 dt \to 0 $$

In section 4.3, we have checked that the Fréchet derivative satisfies the following properties (thanks to the regularity of the vector field, condition C3):

1. $\theta \in \Theta$, $|e_\ell(g+h, \theta) - e_\ell(g, \theta) - De_\ell(g, \theta).h| \le C \|h\|_\infty^2$, and

2. $|De_\ell(g, \theta).h - De_\ell(g', \theta).h| \le C \|h\|_\infty \|g - g'\|_\infty$ with $C$, a constant independent of $\theta, \epsilon$ and $g, g'$.

The same conditions are satisfied for the vector of conditions $e_L(\hat{\phi}, \theta)$, so conditions of theorem 3 in [3] are satisfied, and we can claim that even for a nonlinear functional, we have a root-n rate. Let $v_L(\theta)$ be the Jacobian of the vector of orthogonal conditions $e_L(g, \theta)$ evaluated at $\phi^*$, with entries $v_\ell(\theta) = \frac{f_x(\cdot, \phi^*, \theta)\varphi_\ell + \dot{\varphi}_\ell}{\pi}$. Then, we can claim that

$$ \sqrt{n} (e_L(\hat{\phi}, \theta) - e_L(\phi^*, \theta)) \rightsquigarrow N(0, V_{e,L}(\theta)) $$

and the asymptotic variance is given by $V_{e,L}(\theta) = E[v_L(\theta)v_L(\theta)^\top \sigma^2]$. Moreover, the asymptotic variance can be estimated by the plugin estimate as described in section 4.3 i.e $\hat{V}_{e,L}(\theta) = A(\phi^*, \theta)\hat{V}(\hat{c}_1)A(\phi^*, \theta)^\top + B(\phi^*, \theta)\hat{V}(\hat{c}_2)B(\phi^*, \theta)^\top$. $\square$
## 1.2 Examples

We motivate our work in detail by presenting two models that are relatively common and simple but that nevertheless causes difficulties for estimation.

### 1.2.1 Ricatti ODE

The (scalar) Ricatti equation is defined by a quadratic vector field $f(t, x) = a(t)x^2 + b(t)x + c(t)$ where $a(\cdot), b(\cdot), c(\cdot)$ are time-varying functions. This equation arises naturally in control theory for solving linear-quadratic control problem [35], or in mathematical finance, in the analysis of stochastic interest rate models [7]. We consider one of the simplest case where $a$ is constant, $b = 0$ and $c(t) = c\sqrt{t}$. The objective is to estimate parameters $a, c$ from the noisy observations $y_i = \phi(t_i) + \epsilon_i$ for $t_i \in [0, 14]$. Here the true parameters are $a = 0.11, c = 0.09$ and $\phi_0 = -1$, and one can see the solution and simulated observations in figure 1.1. Although the solution $\phi$ is smooth in the parameters, there exists no closed form and simulations are required for implementing NLS and classical approaches. The hard part in this equation is due to the extreme sensitivity of the squared term in the vector field: for small differences in the parameters or initial condition, the solution can explode before reaching the final time $T = 14$.

1.1. Explosions are not due to numerical instability but to the failure of (theoretical) existence of a global solution on the entire interval (e.g the tangent function is solution of $\dot{\phi} = \phi^2 + 1$, $\phi(0) = 0$ and explodes at $t = \frac{\pi}{2}$). The explosions have to be handled in estimation algorithms and this slows down the exploration of the parameter space (which can be difficult for high-dimensional state or parameter spaces). Nevertheless, we show in the experiment part that NLS or Gradient Matching can do well for parameter estimation, but some additional difficulties does appear when the time-dependent function $c(\cdot)$ has abrupt changes. We consider the case where $c(t) = c\sqrt{t} - d'1_{[T_r,T]}$, $T_r$ is a change-point time, with $d' > 0$. This situation is classical (e.g in engineering) where some input variables $t \mapsto u(t)$ modify the evolution of the system $\dot{\phi} = f(t, \phi(t)) + u(t)$ (typically it can be the introduction of a new chemical species in a reactor at time $T_r$), see figure 1.1. The Cauchy-Lipschitz theory for existence and uniqueness of solutions to time-discontinuous ODE is extended straightforwardly with measure theoretic arguments [35]. The (generalized) solution is defined almost everywhere and belongs to a Sobolev space. For sake of completeness, we provide a generalized version of the Cauchy-Lipschitz theorem for IVP in
## 2.2 Definition

We define a new gradient matching estimator based on (2.4): starting from a nonparametric estimator $\hat{\phi}$, computed from the observations $(t_i, y_i)$, $i = 1, \dots, n$, we want to find the parameter $\theta$ that minimizes the discrepancy between the parametric derivative $t \mapsto f(t, \hat{\phi}(t), \theta)$ and a nonparametric estimate of the derivative, e.g. $\dot{\hat{\phi}}$. A classical discrepancy measure is the $L^2$ distance, that gives rise to the two-step estimator $\hat{\theta}^{TS}$ defined as $\hat{\theta}^{TS} = \arg \min_{\theta \in \Theta} R_{n,w}(\theta)$ where

$$R_{n,w}(\theta) = \int_{0}^{1} |\dot{\hat{\phi}}(t) - f(t, \hat{\phi}(t), \theta)|^2 w(t) dt. \qquad (2.5)$$

This estimator is consistent for several usual nonparametric estimators [4, 12, 8], but the use of a positive weight function $w$ vanishing at the boundaries ($w(0) = w(1) = 0$) is needed to get the classical parametric root-$n$ rate. This approach has the computational advantage of decoupling each equation and enables to have simple estimation procedures. For this reason and ease of readability, we consider only the estimation of the parameter $\theta_1$ when $f$ can be written $f(t, x, \theta) = (f_1(t, x, \theta_1), f_2(t, x, \theta_2))^{\top}$ and $\theta = (\theta_1, \theta_2)^{\top}$ ($\theta_i \in \mathbb{R}^{p_i}$ and $p_1 + p_2 = p$). The joint estimation of $\theta = (\theta_1, \theta_2)^{\top}$ can be done by stacking the observations into a single column: there is no consequence on the asymptotics, but the estimator covariance matrix has to be slightly modified in order to take into account the correlations between the two equations $f_1$ and $f_2$. Having said that, we write simply $f = f_1$ and $\theta = \theta_1$ and we consider only one equation $\dot{x}_1 = f(t, x, \theta)$. We use a nonparametric estimator $\hat{\phi} = (\hat{\phi}_1, \hat{\phi}_2)$ of $\phi^* : [0, 1] \to \mathbb{R}^2$.

Starting from (2.4), a reasonable estimator $\hat{\theta}$ should satisfy the weak formulation

$$\forall \varphi \in C_C^1([0, 1]), \quad \int_0^1 f(t, \hat{\phi}(t), \hat{\theta}) \varphi(t) dt + \int_0^1 \hat{\phi}_1(t) \dot{\varphi}(t) dt = 0. \qquad (2.6)$$

The vector space $C_C^1([0, 1])$ is not tractable for variational formulation, and one prefers Hilbert space with a structure related to $L^2$. In our case, we use $H_0^1 = \{h \in H^1|h(0) = h(1) = 0\}$ which has a simple description within $L^2$: an orthonormal basis is given by the sine functions $t \mapsto \sqrt{2}\sin(\ell\pi t)$, $\ell \ge 1$ and we have

$$H_0^1 = \left\{ \sum_{\ell=1}^{\infty} a_{\ell} \sqrt{2} \sin(\ell \pi t) \left| \sum_{\ell=1}^{\infty} \ell^2 a_{\ell}^2 < \infty \right. \right\} \qquad (2.7)$$
<figure><img src="image_3.png" /><figcaption>Figure 3: Contour plots for Asymptotic Statistical Criterions in Ricatti ODE: (Left) NLS criterion $R_{NLS}^*(d', T_r)$; (Right) OC criterion $Q_L^*(d', T_r)$ ($d'^* = 2, T_r^* = 11$).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 4: Contour plots for Asymptotic Statistical Criterions in Ricatti ODE: (Left) NLS criterion $R_{NLS}^*(d', T_r)$; (Right) OC criterion $Q_L^*(d', T_r)$, ($d'^* = 2, T_r^* = 11$)</figcaption></figure>
When splines with a uniform sequence of knots on [0, 1] are used for estimating $\phi^*$, it is possible to replace condition **C8** by an easier condition to check based on classical regularity assumptions. Indeed, splines are good universal approximating space for smooth functions: typically, if the meshsize $\delta_K$ of the knots of the splines basis $\mathbf{p}^K$ is such that $\delta_K \to 0$ and the function $v_\ell : t \mapsto \frac{f_x(t, \phi^*(t), \theta)\phi_\ell(t)+\phi_\ell(t)}{\pi(t)}$ is $C^1$, then one can find a sequence of splines $s_\ell^K$ such that $\|v_\ell - s_\ell^K\|_{L^2} \to 0$ as $K \to \infty$, see theorem 7.3 in [1]. In that case, we can replace condition **C8** by condition **C9**.

**Proposition 4.3** - Asymptotic Normality of Orthogonal Conditions for splines. We suppose that the function $\phi^*$ is $C^s$ with $s \ge 3$ and that conditions C1, C2(a), C3, C7 and C9 are satisfied, then if $\sqrt{n}K^{-s} \to 0$ and

1. $f$ is linear (in state $x$) and $\frac{K^2}{n} \to 0$,

2. OR, $f$ is nonlinear (in state $x$) and $\frac{K^4}{n} \to 0$

then for all $\theta \in \Theta$,

$$ \sqrt{n} (e_L(\hat{\phi}, \theta) - e_L(\phi^*, \theta)) \rightsquigarrow N(0, V_{e,L}(\theta)) $$

and the asymptotic variance can be estimated by $\widehat{V}_{e,L}(\theta) = A(\phi^*, \theta)\widehat{V}(\hat{c}_1)A(\phi^*, \theta)^\top + B(\phi^*, \theta)\widehat{V}(\hat{c}_2)B(\phi^*, \theta)^\top$.

*Proof.* This proposition is an application of proposition 4.2 to splines with uniform knots. It corresponds also to theorem 9 in [3], which is a general theorem for root-n consistency estimation of functional estimated with (uniform) splines. With this series estimators, conditions **C2** is satisfied, as the condition **C8** (with condition **C9**). $\square$

With propositions 4.1, 4.2 and 4.3, the proof of the theorem 4.1 is rather straightforward as it follows the lines of classical proof for Generalized Moments Estimators.

**Theorem 4.1** - Root-n consistency & asymptotic normality of OC Estimators. If either the following conditions are satisfied:

$\mathbf{p}^K$ is a general series estimators Under conditions C1-C8 and if $f$ is a linear vector field or, $f$ is a non-linear vector field and $K$ is chosen such that $\frac{\zeta_0(K)^4K^2}{n} \to 0$

$\mathbf{p}^K$ is a uniform knot splines Under conditions C1-C2(a), C3-C7, C9 and if $f$ is a linear vector field and $\frac{K^2}{n} \to 0$, or $f$ is a nonlinear vector field and $\frac{K^4}{n} \to 0$

Then $\hat{\theta}_{n,L}$ is such that

$$ \sqrt{n} (\hat{\theta}_{n,L} - \theta_L^*) \rightsquigarrow N(0, \mathbf{V}_L^*) \quad (11) $$

with

$$ \mathbf{V}_L^* = \mathbf{M}_L^* \mathbf{V}_{e,L}^* \mathbf{M}_L^{*\top}. \quad (12) $$

where $\mathbf{V}_{e,L}^* = V_{e,L}(\theta_L^*)$. The asymptotic variance can be estimated by $\widehat{\mathbf{M}}_L \widehat{V}_{e,L}(\hat{\theta}_{n,L}) \widehat{\mathbf{M}}_L^\top \xrightarrow{P} \mathbf{V}_L^*$. In particular, if we use regression splines and $t \mapsto f(t, \phi^*(t), \theta)$ is $C^s$ on [0, 1] with $s \ge 3$, then (11) holds with $K$ such that $\sqrt{n}K^{-s} \to 0$ and $n^{-1}K^4 \to 0$. Moreover, if $L = L(n) \to \infty$, $n \to \infty$ is chosen such that the bias $\mathbf{B}_{L(n)} = O(n^{-1/2})$, then we have

$$ \hat{\theta}_{n,L(n)} - \theta^* = O_P(n^{-1/2}) \quad (13) $$

In particular, this is the case when the test functions $\phi_\ell$ are the sine basis, and $L(n) = O(n^\alpha)$ with $\alpha > 1/2$.

*Proof.* We start with the asymptotic representation given in proposition 4.1:

$$ \hat{\theta}_{n,L} - \theta_L^* = -\mathbf{M}_L^* \mathbf{e}_L(\hat{\phi}, \theta_L^*) + o_P(1) $$

As the parameter $\theta_L^*$ is the global minima of $\theta \mapsto Q_L^*(\theta) = \|\mathbf{e}_L(\phi^*, \theta)\|_2^2$, the first order condition implies that $\mathbf{J}_{\theta,L}(\phi^*, \theta_L^*)^\top \mathbf{e}_L(\phi^*, \theta_L^*) = 0$. By definition of $\mathbf{M}_L^*$, we obtain that $\mathbf{M}_L^* \mathbf{e}_L(\phi^*, \theta_L^*) = 0$ and we can write

$$ (\hat{\theta}_{n,L} - \theta_L^*) = -\mathbf{M}_L^* (\mathbf{e}_L(\hat{\phi}, \theta_L^*) - \mathbf{e}_L(\phi^*, \theta_L^*)) + o_P(1) $$
1 *Supplementary Material I.* This abrupt change causes some difficulties in estimating non-parametrically
2 the solution and its derivative, which can make Gradient Matching less precise. We consider then
3 the estimation of the two additional parameters *d*' and $T_r$. Hence, the parameter estimation problem
4 can be seen as a change-point detection problem, where the solution *φ* still depends smoothly in the
5 parameters. Nevertheless, in the case of the joint estimation of *a*, *c*, *d*' and $T_r$, the particular influence
6 of the parameter $T_r$ makes the problem much more difficult to deal with for classical approaches as it
7 is suggested by the objective functions in *Supplementary Material II*. The variational formulation for
8 model estimation gives a seamless approach for estimating models which possess time discontinuities.

<figure><img src="image_2.png" /><figcaption>Figure 1.1: Solutions of Ricatti ODE with noisy observations ($n = 50, \sigma = 0.4$). Left figure is smooth time-dependent ODE. Right figure has a change-point at time $T_r = 5$ ($d' = 1$).</figcaption></figure>

9

10

1.2.2 Dynamics of Blowfly populations

11
12
13
14
15
16
17
The modeling of the dynamics of population is a classical topic in ecology an more generally in biology.
Differential Equations can describe very precisely the mechanics of evolution, with birth, death and
migration effects. The case of single-species models is the easiest case to consider, as interactions with
rest of the world can be limited, and the acquisition of reliable data is easier. In the 50s, Nicholson
measured quite precisely the dynamics of a blowfly population, known as Nicholson’s experiments [26].
The data are relatively hard to model, and it is common to use Delay Differential Equation (DDE)
whose general form is Ṅ(t) = f (N(t), N(t − τ), θ), in order to account for the almost chaotic behavior
1 Hence, it suffices to consider a countable number of orthogonal conditions (2.6) defined, for instance,
2 with the test functions $\varphi_\ell = \sqrt{2} \sin(\ell \pi t)$, $\forall \ell \ge 1$:

$$
\mathcal{C}_{\ell}(\theta) : \int_{0}^{1} f(t, \hat{\phi}(t), \hat{\theta}) \varphi_{\ell}(t) dt + \int_{0}^{1} \hat{\phi}(t) \dot{\varphi}_{\ell}(t) dt = 0. \qquad (2.8)
$$

3 More generally, we consider a family of orthonormal functions $\varphi_\ell \in H_0^1$, with $\ell \ge 1$, and we introduce
4 the vector space $\mathcal{F} = \overline{\operatorname{span}\{\varphi_\ell, \ell \ge 1\}}$. The vector space $\mathcal{F}$ may not be necessarily dense in $H_0^1$, as
5 the functions $\varphi_\ell$ could be chosen for computational tractability or because of a natural interpretation
6 (for instance B-splines, polynomials, wavelets, ad-hoc functions, ...). For this reason, we introduce the
7 orthogonal decomposition of $H_0^1 = \mathcal{F} \oplus \mathcal{F}^\perp$ where we can have $\mathcal{F} \neq H_0^1$. In general, an estimator $\hat{\theta}$
8 satisfying $\mathcal{C}_\ell(\hat{\theta})$ for $\ell \ge 1$ also approximately satisfies (2.6). However in practice, we will use a finite set
9 of orthogonal constraints defined by $L$ test functions ($L > p$).
10 In order to discuss the influence of the choice of $\mathcal{F}$ and of finite dimensional subspace spanned by
11 $\varphi_1, \dots, \varphi_L$ we introduce the nonlinear operator $\mathcal{E}: (g, \theta) \mapsto \mathcal{E}(g, \theta)$, such that $t \mapsto \mathcal{E}(g, \theta)(t) =$
12 $f(t, g(t), \theta)$.

10 In order to discuss the influence of the choice of $\mathcal{F}$ and of finite dimensional subspace spanned by
11 $\varphi_1, \dots, \varphi_L$ we introduce the nonlinear operator $\mathcal{E}: (g, \theta) \mapsto \mathcal{E}(g, \theta)$, such that $t \mapsto \mathcal{E}(g, \theta)(t) =$
12 $f(t, g(t), \theta)$.

13 For all $\theta$ in $\Theta$ and $g$ in $H^1$, the Fourier coefficients of $\mathcal{E}(g, \theta) - \dot{g}$ in the basis $(\varphi_\ell)_{\ell \ge 1}$ are $e_\ell(g, \theta) =$
14 $\langle \mathcal{E}(g, \theta) - \dot{g}, \varphi_\ell \rangle = \langle \mathcal{E}(g, \theta), \varphi_\ell \rangle + \langle g, \dot{\varphi}_\ell \rangle$, and we introduce the vectors in $\mathbb{R}^L$ $e_L(g, \theta) = (e_\ell(g, \theta))_{\ell=1..L}$
15 and $e_L^*(\theta) = (e_\ell(\phi^*, \theta))_{\ell=1..L}$. Finally, our estimator is defined by the minimization of the quadratic form
16 $Q_{n,L}(\theta) = |e_L(\hat{\phi}, \theta)|^2$:

$$
\hat{\theta}_{n,L} = \arg \min_{\theta \in \Theta} Q_{n,L}(\theta). \tag{2.9}
$$

17 $\hat{\theta}_{n,L}$ is the parameter that “almost” vanishes the first $L$ Fourier coefficients in the orthogonal decompo-
18 sition of $H_0^1 = \mathcal{F} \oplus \mathcal{F}^\perp$:

$$
\mathcal{E}(g, \theta) - \dot{g} = \mathbf{E}_L(g, \theta) + \mathbf{R}_L(g, \theta) + \mathbf{E}_{\mathcal{F}}^{\perp}(g, \theta)
$$

19 with $\mathbf{E}_L(g, \theta) = \sum_{\ell=1}^{L} e_\ell(g, \theta) \varphi_\ell$, $\mathbf{R}_L(g, \theta) = \sum_{\ell>L} e_\ell(g, \theta) \varphi_\ell$ and $\mathbf{E}_{\mathcal{F}}^{\perp}(g, \theta) \in \mathcal{F}^{\perp}$.

20 The function $\mathbf{E}_{\mathcal{F}}^{\perp}(\phi^*, \theta)$ represents the behavior of $\mathcal{E}(g, \theta) - \dot{g}$ at the boundaries of the interval
21 $[0, 1]$. As $\hat{\phi}$ approaches $\phi^*$ asymptotically in supremum norm, the objective function $Q_{n,L}(\theta)$ is close
22 to $Q_{L}^{*}(\theta) = \|\mathbf{E}_{L}(\phi^{*}, \theta)\|_{L^{2}}^{2}$. The discriminative power of $Q_{L}^{*}(\theta)$ can be analyzed locally around $\hat{\theta}$, as
<table><caption>Table 3: MSE, ARE & Tr(V(θ̂)) for Parameter estimation for Ricatti Equation with known T_r</caption><thead><tr><th>×10<sup>-2</sup></th><th colspan="4">MSE</th><th colspan="4">ARE</th><th colspan="3">Tr(V(θ̂))</th></tr><tr><th>(n, σ)</th><th>TS</th><th>OC</th><th>OCopt</th><th>NLS</th><th>TS</th><th>OC</th><th>OCopt</th><th>NLS</th><th>OC</th><th>OCopt</th><th>NLS</th></tr></thead><tbody><tr><td>(400, 0.2)</td><td>0.18</td><td>0.27</td><td>0.51</td><td>0.58</td><td>21.75</td><td>16.47</td><td>26.10</td><td>18.54</td><td>1.76</td><td>1.44</td><td>0.10</td></tr><tr><td>(400, 0.4)</td><td>0.78</td><td>1.21</td><td>1.32</td><td>0.94</td><td>45.03</td><td>31.98</td><td>42.82</td><td>36.25</td><td>2.56</td><td>1.80</td><td>0.38</td></tr><tr><td>(200, 0.2)</td><td>0.33</td><td>0.87</td><td>1.04</td><td>0.57</td><td>27.09</td><td>26.95</td><td>39.63</td><td>25.81</td><td>2.85</td><td>2.07</td><td>0.25</td></tr><tr><td>(200, 0.4)</td><td>1.12</td><td>2.69</td><td>2.62</td><td>1.12</td><td>53.99</td><td>43.17</td><td>54.04</td><td>45.10</td><td>5.64</td><td>2.98</td><td>0.98</td></tr><tr><td>(50, 0.2)</td><td>1.03</td><td>1.30</td><td>2.08</td><td>1.54</td><td>45.53</td><td>41.50</td><td>49.98</td><td>49.12</td><td>4.70</td><td>5.39</td><td>1.00</td></tr><tr><td>(50, 0.4)</td><td>3.80</td><td>4.43</td><td>5.06</td><td>3.94</td><td>77.08</td><td>73.69</td><td>73.30</td><td>70.96</td><td>8.89</td><td>8.76</td><td>4.08</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Figure 5: Coverage Probabilities for the 95% Confidence Ellipse for Ricatti model with known T_r</figcaption></figure>
and the application of proposition 4.3 gives directly the root-n consistency and normality of $\hat{\theta}_{n,L}$ and the asymptotic variance $\lim_n \sqrt{nV}(\hat{\theta}_{n,L})$ is directly derived from the asymptotic variance of the orthogonal conditions i.e it equals $M_L^*V_{e,L}(\theta_L^*)M_L^{\top}$. The latter can be estimated by the plugin estimates $\tilde{V}_{e,L}(\hat{\theta}_{n,L})$ and $M_L^*$ is estimated by using the observed Jacobian matrix $J_{\theta,L}(\hat{\phi}, \hat{\theta}_{n,L})$ (following the estimator of proposition 4.2).

Finally, we can include the bias $r_L = \theta_L^* - \theta^*$, given in theorem 1, in the asymptotics by assuming that we can chose $L = L(n)$ such that $r_n = o(n^{-1/2})$. From lemma 5 (and theorem 1), we can have an estimate of the bias due to the truncation of the Fourier expansion and $r_L = o(1/L)$ i.e $r_n = o(n^{-1/2})$. Then, the influence of the bias is negligible as

$$ \begin{aligned} \hat{\theta}_{n,L(n)} - \theta^* &= \hat{\theta}_{n,L(n)} - \theta_L^*(n) + \theta_L^*(n) - \theta^* \\ &= O_P(n^{-1/2}) + o(n^{-1/2}) \\ &= O_P(n^{-1/2}) \end{aligned} $$

and asymptotic behavior (13) can hold.

□

# 5 From optimal Weighting to a practical algorithm

Theorem 4.1 is of practical interest as it provides a closed-form expression for the asymptotic variance of $\hat{\theta}_{n,L}$ for general ODE. We can make a parallel between the orthogonal condition estimator and the weighted nonlinear least-squares:

$$ \hat{\theta}_c^{WLS} = \arg \min \sum_{i=1}^{n} w(t_i) |y_i - \phi(t_i, \theta_c)|^2 $$

where $w(\cdot)$ is a (positive) weight function and $\theta_c = (\phi_0, \theta)$. $\hat{\theta}_c^{WLS}$ is consistent and $\sqrt{n}(\hat{\theta}_c^{WLS} - \theta_c^*) \rightsquigarrow N(0, V^{AWLS})$ under classical regularity assumptions on $f$, see [5]. The asymptotic variance $V^{AWLS}$ is directly computed from the sensitivity equations, and the optimal weight function $w(\cdot)$ is proportional to the variance function $\sigma^2(\cdot)$, meaning that the unweighted least-squares estimator is optimal in the homoscedastic case. It is clear that $V_L^* \neq V^{AWLS}$, which means that we are sure that the orthogonal condition estimator is not an efficient estimator (at least in the Gaussian case, because NLS and MLE are the same in that case). A striking difference between $V_L^*$ and $V^{AWLS}$ is that the least squares involves the Jacobian of the solution w.r.t. the initial values and parameters, whereas the orthogonal conditions involve the Jacobian of the vector field. It is then hard to compare these two matrices in generality and to evaluate the loss of efficiency. Nevertheless, we can compare the influence of a weight matrix for the minimization of the orthogonal conditions. Indeed, if we introduce a positive definite matrix $W$ in $\mathbb{R}^{L \times L}$, we can define the weighted criterion

$$ Q_{n,L}^{W}(\theta) = e_{L}(\hat{\phi}, \theta)^{\top} W e_{L}(\hat{\phi}, \theta) $$

and the corresponding estimator

$$ \hat{\theta}_{n,L}^{W} = \arg \min_{\theta \in \Theta} Q_{n,L}^{W}(\theta). $$

The results of theorems 3 and 4 are then still true under straightforward adaptations. In particular, $\hat{\theta}_{n,L}^W$ is consistent and asymptotically normal with asymptotic variance $M_{W,L}^*V_{e,L}^*M_{W,L}^{\top}$ where $M_{W,L}^* = [J_L^{*\top}WJ_L^{*\top}]^{-1}J_L^{*\top}W$. Eventually, we can then ask for the best weighting matrix $W$ giving the smallest asymptotic variance. This is a classical result for Generalized Moment Estimators that we recall in the following proposition (see for instance section 3.6 in [2]):

**Proposition 5.1. Optimal weighting matrix**

The minimal asymptotic variance for $\hat{\theta}_{n,L}^W$ is obtained with $W^{opt} = V_{e,L}^{*-1}$ and

$$ V_{L}^{opt}(\theta^{*}) = (J_{\theta,L}^{*\top} V_{e,L}^{*-1} J_{\theta,L}^{*})^{-1} \quad (14) $$

Even if we have an homoscedastic model, we have an interest in using a weighted estimator. Interestingly, the problem of choosing the best weighting matrix is directly related to the choice of the best set of test functions $\phi_1, \dots, \phi_L$. Indeed, the diagonalization of $V_{e,L}^* = U\Lambda U^\top$ permits the introduction of the eigenvalues
<figure><img src="image_3.png" /><figcaption>Figure 1.2: Blowfly Data, collected by Nicholson</figcaption></figure>

1 of the data, see figure 1.2. Nicholson’s dataset is now a classical benchmark for evaluating time series
2 algorithms due its intrinsic complexity. Nevertheless, the following DDE is commonly acknowledged as
3 a correct model [11, 23]:

$$
\dot{N} = PN(t - \tau) \exp(-N(t - \tau)/N_0) - \delta N(t) \quad (1.3)
$$

4 whose parameter fitting (of P, N₀, δ) remains delicate. In particular classical NLS are difficult to use
5 in this setting as the initial condition, which is a function defined on [-τ, 0], is unknown. Alternative
6 solutions, such as Gradient Matching or Bayesian Methods (based on ABC, [41]) give reliable estimates
7 that reproduce the observed dynamics without estimation of the initial condition. These aforementioned
8 methods use particular statistics or functions of the model that provides high-level information on the
9 parameters. The Orthogonal Conditions estimator has a similar approach for dealing with the estimation
10 of Differential Equations.
11

11
1 it behaves approximately as the quadratic form $Q_L^*(\theta) \approx (\theta - \theta^*)^\top \mathbf{J}_{\theta,L}^{*\top} \mathbf{J}_{\theta,L}^* (\theta - \theta^*)$ where $\mathbf{J}_{\theta,L}^*$ is the
2 matrix in $\mathbb{R}^{L \times p}$ with entries $\int_0^1 f_{\theta_j}(t, \phi^*(t), \theta^*) \varphi_\ell(t) dt$, for $j = 1, \dots, p$, $\ell = 1, \dots, L$.

3 Consistency of the Orthogonal Conditions estimator

4 In order to obtain precise results with closed-form expression for the bias and variance estimators,
5 we consider series estimators, i.e. estimators expressed as $\hat{\phi}_j = \sum_{k=1}^K \hat{c}_{k,j} p_{kK} = \hat{c}_j p^K$, where $p^K =
6 (p_{1K}, \dots, p_{kK})$ is a vector of approximating functions and the coefficients $\hat{c}_j = (\hat{c}_{k,j})_{k=1..K}$ are computed
7 by least squares. For notational simplicity, we use the same functions (and the same number $K$) for
8 estimating $\phi_1^*$ and $\phi_2^*$. We denote $P^K = (p_{kK}(t_i))_{1\le i,k \le n,K}$ the design matrix and $\mathbf{Y}_j = (y_{i,j})_{i=1..n}$ the
9 vectors of observations. Hence, the estimated coefficients $\hat{c}_j = (P^{K\top} P^K)^\dagger P^{K\top} \mathbf{Y}_j$ (where $\dagger$ denotes a
10 generalized inverse) gives rise to the so-called hat matrix $\mathbf{H} = P^K (P^{K\top} P^K)^\dagger P^{K\top}$ and the vector of
11 smoothed observations is $\hat{\phi}_j = \mathbf{H} \mathbf{Y}_j$, $j=1,2$. One can typically think of regression splines, [22]. We
12 introduce now the conditions required for the definition and consistency of our estimator.

13 **Condition C1 (a)** $\Theta$ is a compact set of $\mathbb{R}^p$ and $\theta^*$ is an interior point of $\Theta$, $\mathcal{X}$ is an open subset of $\mathbb{R}^2$; (b) $(t, x) \mapsto f(t, x, \theta^*)$ is $L^2$-Lipschitz and $L^2$-Caratheodory.

15 **Condition C2 (a)** $(Y_i, t_i)$ are i.i.d. with variance $V(Y|T=t) = \Sigma_\epsilon = \sigma^2 I_2$; (b) For every $K$, there
16 is a nonsingular constant matrix $B$ such that for $P^K = B_p^K(t)$; (i) the smallest eigenvalue of
17 $E[P^K(T)P^K(T)^\top]$ is bounded away from zero uniformly in $K$ and (ii) there is a sequence of
18 constants $\zeta_0(K)$ satisfying $\sup_t |P^K(t)| \le \zeta_0(K)$ and $K = K(n)$ such that $\zeta_0(K)^2 K/n \to 0$ as
19 $n \to \infty$; (c) There are $\alpha, \mathbf{c}_{1,K}, \mathbf{c}_{2,K}$ such that $\|\phi_j^* - p^K \mathbf{c}_{j,K}\|_\infty = \sup_{[0,1]} |\phi_j^*(t) - p^K(t)^\top \mathbf{c}_{j,K}| = O(K^{-\alpha})$.

21 **Condition C3** There exists $D > 0$, such that the $D$-neighborhood of the solution range $\mathcal{D} = \{x \in \mathbb{R}^2 | \exists t \in [0, 1], |x - \phi^*(t)| < D\}$ is included in $\mathcal{X}$ and $f$ is $C^2$ in $(x, \theta)$ on $\mathcal{D} \times \Theta$ for $t$ in $[0, 1]$ a.e. Moreover, the derivatives of $f$ w.r.t $x$ and $\theta$ (with obvious notations) $f_x, f_\theta, f_{xx}, f_{x\theta}$ and $f_{\theta\theta}$ are $L^2$ uniformly bounded on $\mathcal{D} \times \Theta$ by $L^2$ functions $\bar{h}_x, \bar{h}_\theta, \bar{h}_{x\theta}, \bar{h}_{xx}$ and $\bar{h}_{\theta\theta}$ (respectively).

25 **Condition C4** Let $(\varphi_\ell)_{\ell \ge 1}$ be an orthonormal sequence of functions in $H_0^1$.
2.1.2 Unknown change-point time $T_r$

<table><caption>Table 4: MSE ARE & Sum Empirical Variance for Parameter estimation for Ricatti with unknown $T_r$</caption><thead><tr><th>×10<sup>-2</sup></th><th colspan="2">MSE(α̂)</th><th colspan="2">MSE(ĉ)</th><th colspan="2">MSE(d̂)</th><th colspan="2">MSE(T̂<sub>r</sub>)</th></tr><tr><th>(n, σ)</th><th>OC</th><th>OCopt</th><th>OC</th><th>OCopt</th><th>OC</th><th>OCopt</th><th>OC</th><th>OCopt</th></tr></thead><tbody><tr><td>(400, 0.2)</td><td>0.09</td><td>0.07</td><td>0.00</td><td>0.00</td><td>2.54</td><td>2.35</td><td>1.39</td><td>1.15</td></tr><tr><td>(400, 0.4)</td><td>0.29</td><td>0.26</td><td>0.01</td><td>0.01</td><td>4.27</td><td>3.63</td><td>3.54</td><td>2.95</td></tr><tr><td>(200, 0.2)</td><td>0.21</td><td>0.21</td><td>0.00</td><td>0.00</td><td>4.08</td><td>4.44</td><td>3.18</td><td>3.21</td></tr><tr><td>(200, 0.4)</td><td>0.61</td><td>0.55</td><td>0.01</td><td>0.01</td><td>11.96</td><td>12.80</td><td>6.93</td><td>7.88</td></tr><tr><td>(50, 0.4)</td><td>0.64</td><td>0.48</td><td>0.02</td><td>0.02</td><td>11.20</td><td>14.06</td><td>14.25</td><td>15.92</td></tr><tr><td>(50, 0.4)</td><td>0.77</td><td>0.89</td><td>0.01</td><td>0.02</td><td>17.18</td><td>18.63</td><td>19.40</td><td>18.43</td></tr></tbody></table><table><thead><tr><th>×10<sup>-2</sup></th><th colspan="2">MSE</th><th colspan="2">ARE</th><th colspan="2">Tr(V(θ̂))</th></tr><tr><th>(n, σ)</th><th>OC</th><th>OCopt</th><th>OC</th><th>OCopt</th><th>OC</th><th>OCopt</th></tr></thead><tbody><tr><td>(400, 0.2)</td><td>4.01</td><td>3.57</td><td>30.90</td><td>27.67</td><td>3.97</td><td>3.53</td></tr><tr><td>(400, 0.4)</td><td>8.11</td><td>6.85</td><td>56.12</td><td>54.42</td><td>8.02</td><td>6.62</td></tr><tr><td>(200, 0.2)</td><td>7.47</td><td>7.87</td><td>44.56</td><td>45.11</td><td>7.35</td><td>7.77</td></tr><tr><td>(200, 0.4)</td><td>19.51</td><td>21.23</td><td>78.91</td><td>73.47</td><td>18.94</td><td>20.82</td></tr><tr><td>(50, 0.2)</td><td>26.10</td><td>30.48</td><td>84.33</td><td>81.99</td><td>5.14</td><td>8.61</td></tr><tr><td>(50, 0.4)</td><td>37.36</td><td>37.96</td><td>98.67</td><td>106.39</td><td>9.49</td><td>12.90</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 6: Coverage Probabilities for the 95% Confidence Ellipse for Ricatti model with unknown $T_r$</figcaption></figure>
$\Lambda = \text{diag}(\lambda_1, \dots, \lambda_L)$ (with $\lambda_1 \ge \dots \ge \lambda_L \ge 0$) and $U = (u_1|\dots|u_L)$ is an orthogonal matrix. If $\Lambda$ is nonsingular, $Q_{n,L}^{\text{opt}}(\theta) \doteq Q_{n,L}^{W^{\text{opt}}}(\theta) = \tilde{e}_L(\hat{\phi}, \theta)^{\top}\Lambda^{-1}\tilde{e}_L(\hat{\phi}, \theta)$, with $\tilde{e}_L(\hat{\phi}, \theta) = U^{\top}e_L(\hat{\phi}, \theta)$. By linearity, the new orthogonal conditions $\tilde{e}_L(\hat{\phi}, \theta)$ can be written as orthogonal conditions with the test functions $\psi_1, \dots, \psi_L$ derived from the eigenvectors of $V_{e,L}^*$ as $\psi_\ell = \sum_{k=1}^L u_{k\ell}\phi_k$. The use of the optimal weighting matrix is then equivalent to choose the best test functions in $\mathcal{F}_L$ and the diagonalization show that some care must be taken (typically $L$ have to be reasonably small to avoid $\lambda_L \approx 0$).

The optimal weighting matrix $W^{\text{opt}}$ given in proposition 4.2 is dependent on the true unknown parameter: obviously, this makes this result hard to use in practice. We suggest to approximate $V_{e,L}^{*-1}$ by computing a sequence of weighted estimates, as described in Algorithm 1.

<caption>**Algorithm 1** Iteratively Reweighted Orthogonal Conditions - IRWOC</caption>

**Require:** $\hat{\phi}, \theta_{\text{init}} \in \Theta, \epsilon > 0$, Compute the unweighted estimator

$$ \hat{\theta}_{n,L}^{(0)} = \arg \min_{\theta} Q_{n,L}(\theta). $$

Compute the asymptotic covariance of $e_L(\hat{\phi}, \hat{\theta}_{n,L}^{(0)})$ with equation in Proposition 4.3.

**while** $|\hat{\theta}_{n,L}^{(k)} - \hat{\theta}_{n,L}^{(k-1)}| > \epsilon$ **do**

$$ \left\{ \begin{array}{l} W^{(k)} = \left( V_{e,L}(\hat{\theta}_{n,L}^{(k)}) \right)^{-1} \\ \hat{\theta}_{n,L}^{(k+1)} = \arg \min_{\theta} Q_{n,L}^{W^{(k)}}(\theta) \end{array} \right. . \qquad (15) $$

**end while**

This is known as iterated Generalized Method of Moments (GMM), and consists in successive re-weighting, by using the consistent estimator of the variance $V_{e,L}(\hat{\theta}_{n,L})$, [2]. There is no theoretical guarantee for the (numerical) convergence of this algorithm but if the model is correct, this sequence of iterations tends to ameliorate the quality of the GMM estimator. Possibly, during this stage, it could be necessary to select appropriately $L$, as $V_{e,L}(\hat{\theta}_{n,L}^{(k)})$ can be close to be nonsingular at an iteration $k$.

## References

[1] R.A. DeVore and G.G. Lorentz. *Constructive Approximation*, volume 303 of *Grundlehren der mathematischen Wissenschaften*. Springer-Verlag, 1993.

[2] A.R. Hall. *Generalized Method of Moments*. Oxford University Press, 2005.

[3] W. K. Newey. Convergence rates and asymptotic normality for series estimators. *Journal of Econometrics*, 79:147-168, 1997.

[4] D. O'Regan. *Existence theory for nonlinear ordinary differential equations*. Mathematics and its applications. Kluwer, 1997.

[5] L. Pronzato. Optimal experimental design and some related control problems. *Automatica*, 44:303-325, 2008.

[6] A.W. van der Vaart. *Asymptotic Statistics*. Cambridge Series in Statistical and Probabilities Mathematics. Cambridge University Press, 1998.
# 2 Differential Equation Model and Gradient Matching

## 2.1 ODE models and Gradient Matching

For ease of readability, we focus on a two-dimensional system of ODEs. In our case, as there is no computational and theoretical differences between the situation $d = 2$ and $d > 2$, there is no lack of generality by this assumption. We consider noisy observations $Y_1, \dots, Y_n \in \mathbb{R}^2$ of the function $\phi^*$ measured at random times $t_1 < \dots < t_n \in [0, 1]$:

$$Y_i = \phi^*(t_i) + \epsilon_i \qquad (2.1)$$

where $\epsilon_1, \dots, \epsilon_n$ are i.i.d with $E(\epsilon_i) = 0$ and $V(\epsilon_i) = \sigma^2 I_2$. We suppose that the regression function $\phi^*$ belongs to the Sobolev space $H^1 = \{u \in L^2([0, 1]) \mid \dot{u} \in L^2([0, 1])\}$, and $\phi^*$ is a solution to the parametrized Ordinary Differential Equation (1.1), i.e. there exists a true parameter $\theta^* \in \Theta \subset \mathbb{R}^p$ such that for $t \in [0, 1]$ almost everywhere (a.e.)

$$\dot{\phi}^*(t) = \mathbf{f}(t, \phi^*(t), \theta^*) \qquad (2.2)$$

where $\mathbf{f} = (f_1, f_2)$ is a vector field from $[0, 1] \times \mathcal{X} \times \Theta$ to $\mathbb{R}^2$, where $\mathcal{X} \subset \mathbb{R}^2$.

The statistical problem can be seen as a noisy version of a parametrized Multipoint Boundary-Value Problem (MBVP, [4]). MBVP deals with the existence, uniqueness and computation of a solution $\phi^*$ to equation (1.1), with general boundary conditions $\phi^*(t_1) = y_1, \dots, \phi^*(t_n) = y_n, n \ge 2$. Obviously, MBVP is a much more difficult problem than the classical Initial Value Problem although some theoretical results do exist in some restricted cases ([3, 27] and references therein). On the computational side, numerous algorithms such as collocation, multiple shooting,... have been proposed to solve general Boundary Value Problems, [2]. Among them, the 2 points Boundary Value Problem (BVP) where $G(\phi^*(0), \phi^*(1)) = 0$ with $G$ a given function, is one of the most common and important one, as it arises in numerous applications (physics, control theory, ...). We emphasize that a convenient way to deal theoretically and computationally with BVP, in particular linear second order differential ODEs, is not based on an adaptation of the IVP theory, but it rather involves elaborated concepts from functional
1 **Condition C5** $\theta^*$ is the unique global minimizer of $Q_{\mathcal{F}}^*$ and $\inf_{|\theta-\theta^*|>\epsilon} Q_{\mathcal{F}}^*(\theta) > 0$.

2 **Condition C6** There exists $L_0$ such that for $L \ge L_0$, $J_{\theta,L}(g, \theta)$ is full rank in a neighborhood of $(\phi^*, \theta^*)$.

3 Condition C1 gives the existence and uniqueness of a solution $\phi^*$ in $H^1$ to the IVP for $\theta = \theta^*$ and $x(0) = \phi^*(0)$. If $f$ is continuous in $t$ and $x$, then the derivative $\dot{\phi}^*(t) = f(t, \phi^*(t), \theta^*)$ can be defined on $]0, 1[$ and is also continuous.

6 Under condition C2 (satisfied among others by regression splines with $\zeta_0(K) = \sqrt{K}$), it is known that the series estimator $\hat{\phi}_j$ are consistent estimators of $\phi_j^*$ for usual norms, in particular $\|\hat{\phi}_j - \phi_j^*\|_{\infty} = O_P(\zeta_0(K)(\sqrt{K/n} + K^{-\alpha}))$ (theorem 1, [16]). If $\phi^*$ is $C^s$ and we use splines then $\alpha = s$ and $\|\hat{\phi} - \phi^*\|_{\infty} = O_P(K/\sqrt{n} + K^{1/2-s})$.

11 Condition C3 is here to control the continuity and regularity of the function $\mathcal{E}$ involved in the inverse problem. Moreover, it provides uniform control needed for stochastic convergence.

14 Condition C4 is a sufficient condition for deriving independent conditions $C_\ell(\theta)$, and normalization is useful only to avoid giving implicitly more weight to a condition w.r.t. the other conditions. A principled choice of weights and orthogonal functions will be discussed in section 4.4.

19 Condition C5 is standard in M-estimation [24], but can be hard to check in practice. Indeed, the parametric identifiability of ODE models can be hard to show, even for small systems [14]. Moreover, the natural criterion for estimating $\theta$ and for identifiability analysis is

$$Q^*(\theta) = \|\mathcal{E}(\phi^*, \theta) - \mathcal{E}(\phi^*, \theta^*)\|_{L^2}^2$$

22 but $\|\mathbf{E}_{\mathcal{F}}^{\perp}(\phi^*, \theta)\|_{L^2}^2$ is withdrawn and we use the quadratic form $Q_{\mathcal{F}}^*(\theta)$ in order to avoid boundary effects.

23 This is needed in order to get a parametric rate of convergence, as in the original two-step criterion (2.5).

24 As a consequence, we lose a piece of information brought by the trajectory $t \mapsto \phi^*(t)$ and we have to be sure that the parameter $\theta$ has a low influence on $\|\mathbf{E}_{\mathcal{F}}^{\perp}(\phi^*, \theta)\|_{L^2}^2$. A favorable case is that it is almost
## 2.2 FitzHugh-Nagumo equation

The FitzHugh-Nagumo is a nonlinear two-dimensional ODE introduced for modeling neurons. For well-chosen sets of parameters and initial conditions, it exhibits a periodic behavior, with typical oscillations corresponding to a limit cycle.

$$ \begin{cases} \dot{V} = c \left(V - \frac{V^3}{3} + R\right) \\ \dot{R} = -\frac{1}{c} (V - a + bR) \end{cases} \qquad (1) $$

The true parameters are $a^* = b^* = 0.2$ and $c^* = 3$ and $x_0 = (V_0, R_0) = (-1, 1)$, and are taken from Ramsay et al (2007) where it was introduced as a benchmark for parameter estimation in ODEs. Due to the periodicity of the FitzHugh-Nagumo solution for this parameter set we use the sine basis for the test function and we choose the best number of orthogonal condition between 5 or 6.

<table><thead><tr><th rowspan="2">×10<sup>-2</sup><br>(n, σ)</th><th colspan="4">MSE(̂a)</th><th colspan="4">MSE(̂b)</th><th colspan="4">MSE(̂c)</th></tr><tr><th>TS</th><th>OC</th><th>OCopt</th><th>NLS</th><th>TS</th><th>OC</th><th>OCopt</th><th>NLS</th><th>TS</th><th>OC</th><th>OCopt</th><th>NLS</th></tr></thead><tbody><tr><td>(400, 0.15)</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.00</td><td>0.09</td><td>0.08</td><td>0.07</td><td>0.05</td><td>3.54</td><td>0.94</td><td>0.12</td><td>0.01</td></tr><tr><td>(400, 0.3)</td><td>0.04</td><td>0.04</td><td>0.03</td><td>0.01</td><td>0.36</td><td>0.31</td><td>0.26</td><td>0.17</td><td>23.29</td><td>4.70</td><td>0.49</td><td>0.04</td></tr><tr><td>(200, 0.15)</td><td>0.02</td><td>0.02</td><td>0.01</td><td>0.01</td><td>0.20</td><td>0.18</td><td>0.14</td><td>0.08</td><td>9.80</td><td>2.01</td><td>0.25</td><td>0.02</td></tr><tr><td>(200, 0.3)</td><td>0.06</td><td>0.08</td><td>0.05</td><td>0.01</td><td>0.72</td><td>0.68</td><td>0.62</td><td>0.33</td><td>49.58</td><td>9.97</td><td>0.97</td><td>0.07</td></tr><tr><td>(50, 0.15)</td><td>0.06</td><td>0.07</td><td>0.11</td><td>0.02</td><td>0.77</td><td>0.73</td><td>1.03</td><td>0.29</td><td>168.01</td><td>36.54</td><td>3.07</td><td>0.06</td></tr><tr><td>(50, 0.3)</td><td>0.25</td><td>0.28</td><td>0.32</td><td>0.06</td><td>1.96</td><td>2.14</td><td>3.13</td><td>1.10</td><td>173.38</td><td>58.76</td><td>4.47</td><td>0.25</td></tr></tbody></table>

<table><caption>Table 5: MSE, ARE & Sum Asymptotic Variance for Parameter estimation for FitzHugh-Nagumo model</caption><thead><tr><th rowspan="2">×10<sup>-2</sup><br>(n, σ)</th><th colspan="4">MSE</th><th colspan="4">ARE</th><th colspan="3">Tr(V(θ))</th></tr><tr><th>TS</th><th>OC</th><th>OCopt</th><th>NLS</th><th>TS</th><th>OC</th><th>OCopt</th><th>NLS</th><th>OC</th><th>OCopt</th><th>NLS</th></tr></thead><tbody><tr><td>(400, 0.15)</td><td>3.64</td><td>1.04</td><td>0.20</td><td>0.07</td><td>21.60</td><td>18.12</td><td>15.02</td><td>11.20</td><td>1.09</td><td>0.17</td><td>0.02</td></tr><tr><td>(400, 0.3)</td><td>23.69</td><td>5.05</td><td>0.77</td><td>0.22</td><td>46.77</td><td>36.06</td><td>28.69</td><td>20.30</td><td>3.99</td><td>0.68</td><td>0.09</td></tr><tr><td>(200, 0.15)</td><td>10.02</td><td>2.21</td><td>0.40</td><td>0.11</td><td>32.73</td><td>26.12</td><td>20.98</td><td>14.64</td><td>1.95</td><td>0.3</td><td>0.05</td></tr><tr><td>(200, 0.3)</td><td>50.36</td><td>10.73</td><td>1.64</td><td>0.41</td><td>66.98</td><td>52.63</td><td>42.82</td><td>28.76</td><td>6.62</td><td>1.29</td><td>0.16</td></tr><tr><td>(50, 0.15)</td><td>168.84</td><td>37.35</td><td>4.20</td><td>0.37</td><td>80.47</td><td>58.13</td><td>56.39</td><td>26.75</td><td>4.16</td><td>0.77</td><td>0.18</td></tr><tr><td>(50, 0.3)</td><td>175.59</td><td>61.17</td><td>7.91</td><td>1.41</td><td>122.27</td><td>107.11</td><td>103.35</td><td>53.40</td><td>9.82</td><td>3.29</td><td>0.70</td></tr></tbody></table>
analysis such as weak derivative, variational formulation and Sobolev spaces [10]. If we denote the inner product of $L^2$ as $\forall \varphi, \psi \in L^2([0, 1])$, $\langle \varphi, \psi \rangle = \int_0^1 \varphi(t)\psi(t)dt$, the weak derivative of the function $g$ in $H^1$ is not defined point-wise but as the function $\dot{g} \in L^2$ satisfying $\langle \dot{g}, \varphi \rangle = -\langle g, \dot{\varphi} \rangle$, for all function $\varphi$ in $C^1$ with support included in $]0, 1[$ (denoted $C_C^1([0, 1])$). Of course, if $t \mapsto \phi(t, x_0, \theta)$ is a $C^1$ function on $]0, 1[$, the classical derivative $\dot{\phi}$ is also the weak derivative. We introduce then the (weak) variational formulation of the ODE (1.1): a weak solution $g$ to (1.1) is a function in $H^1$ such that $\forall \varphi \in C_C^1([0, 1])$

$$ \int_0^1 \mathbf{f}(t, g(t), \theta)\varphi(t)dt + \int_0^1 g(t)\dot{\varphi}(t)dt = 0 \qquad (2.3) $$

This variational formulation is the key of the Finite Elements Method which is the reference approach for solving Boundary Value Problems and Partial Differential Equations, [6]. In the case of ODEs, this formulation is not well used for computing solutions, because the geometry of the (1-D) interval $]0, 1[$ is simple, and it is easy to build a spline approximation by collocation that solves approximately the ODE. Nevertheless, the characterization (2.3) is useful for the statistical inference task, as it enables to give necessary conditions for a good estimate. In particular, we emphasize that we do not solve the ODE, but we want to identify a parameter $\theta$ indexing the vector field $\mathbf{f}$. Hence, we develop a new algorithmic approach, different from the one used for solving the direct problem.

## 2.2 Definition

We define a new gradient matching estimator based on (2.3): starting from a nonparametric estimator $\hat{\phi}$, computed from the observations $(t_i, y_i)$, $i = 1, \dots, n$, we want to find the parameter $\theta$ that minimizes the discrepancy between the parametric derivative $t \mapsto \mathbf{f}(t, \hat{\phi}(t), \theta)$ and a nonparametric estimate of the derivative, e.g. $\hat{\phi}$. A classical discrepancy measure is the $L^2$ distance, that gives rise to the two-step estimator $\hat{\theta}^{TS}$ defined as $\hat{\theta}^{TS} = \arg \min_{\theta \in \Theta} R_{n,w}(\theta)$ where

$$ R_{n,w}(\theta) = \int_0^1 |\dot{\hat{\phi}}(t) - \mathbf{f}(t, \hat{\phi}(t), \theta)|^2 w(t)dt. \qquad (2.4) $$

This estimator is consistent for several usual nonparametric estimators [8, 19, 14], but the use of a positive weight function $w$ vanishing at the boundaries ($w(0) = w(1) = 0$) is needed to get the classical
constant on $\Theta$, so that $Q^*$ and $Q_{\mathcal{F}}^*$ are essentially the same functions, with the same global minimum and the same discriminating power.

Finally, Condition C6 is about the influence of the truncation. We use only the first $L$ Fourier coefficients of $\mathcal{E}(g, \theta) - \dot{g}$ to identify the parameter $\theta$, but this might not be sufficient to discriminate between two parameters $\theta$ and $\theta'$. In a way, we perform dimension reduction but we need to be sure that we have an exact recovery when $L$ goes to infinity. More precisely, we expect that the minimum of $|e_L^*(\theta)|^2$ found under condition C5 is also a minima of $Q_{\mathcal{F}}^*(\theta) = \|\mathcal{E}_{\mathcal{F}}(\phi^*, \theta)\|_{L^2}^2$. We have to be sure that $\theta^*$ is a global and isolated minima of $Q_{\mathcal{F}}^*(\theta)$. More generally, we can introduce the Jacobian matrices $\mathbf{J}_{\theta,L}(g, \theta)$ in $\mathbb{R}^{L \times p}$ with entries $\int_0^1 f_{\theta_j}(t, g(t), \theta)\varphi_\ell(t)dt$ and $\mathbf{J}_{x,L}(g, \theta)$ in $\mathbb{R}^{L \times d}$ with entries $\int_0^1 f_{x_i}(t, g(t), \theta)\varphi_\ell(t)dt$. For this reason, we suppose that $\mathbf{J}_{\theta,L}^*$ is full rank, so that $Q_L^*(\theta)$ is locally strictly convex, with a unique local minimum $\theta^*$.

**Theorem 1.** If conditions C1 to C6 are satisfied, then

$$ \hat{\theta}_{n,L} - (\theta^* + r_L) = O_P(1) $$

and the bias $r_L$ tends to zero as $L \to \infty$. If we use the sine basis and if $\mathcal{E}(\phi^*, \theta)$ is in $H^1$ for all $\theta$, then $r_L = o(1/L)$.

## 4 Asymptotics

We give a precise description of the asymptotics of $\hat{\theta}_{n,L}$ (rate, variance and normality), by exploiting the well-known properties of series estimators. We consider the linear case, then we extend the obtained results to general nonlinear ODEs. We show in a preliminary step that the asymptotics of $\hat{\theta}_{n,L} - \theta_L^*$ are directly related to the behavior of $e_L(\hat{\phi}, \theta^*)$, which is a classical feature of Moment Estimators.
<figure><img src="image_7.png" /><figcaption>Figure 7: Ellipse set for parameter estimation for FitzHugh-Nagumo model</figcaption></figure>

Results are presented in table 5 and Figure 7. We can see that NLS estimators is the best among all. Among the non-parametric estimators, the optimal one is OC opt, both in term of MSE and ARE. A parameter by parameter study shows us IRWOC algorithms dramatically improve the accuracy of the estimation of *c*.
parametric root-n rate. The importance of the weight function w for the asymptotics of $\hat{\theta}^{TS}$ is assessed by theorem 3.1 in [8]. Indeed, if w does not vanish at the boundaries, then $\hat{\theta}^{TS}$ does not have a root-n rate, because the asymptotics is then dominated by the nonparametric estimates $\hat{\phi}(0)$ and $\hat{\phi}(1)$. The usefulness of such weighting function is well acknowledged in nonparametric or semiparametric estimation. For instance, the so-called weighted average derivative is based on a similar weight function in order to get estimators with parametric rate in partial index models [25].

The use of a nonparametric proxy (instead of a solution to be computed) gives the opportunity to consider parameter estimation in $f_1$ and in $f_2$ separately. For this reason and ease of readability, we consider only the estimation of the parameter $\theta_1$ when $f$ can be written $f(t, x, \theta) = (f_1(t, x, \theta_1), f_2(t, x, \theta_2))^{\top}$ and $\theta = (\theta_1, \theta_2)^{\top}$ ($\theta_i \in \mathbb{R}^{p_i}$ and $p_1 + p_2 = p$). The joint estimation of $\theta = (\theta_1, \theta_2)^{\top}$ can be done by stacking the observations into a single column: there is no consequence on the asymptotics, but the estimator covariance matrix has to be slightly modified in order to take into account the correlations between the two equations $f_1$ and $f_2$. Having said that, we write simply $f = f_1$ and $\theta = \theta_1$ and we consider only one equation $\dot{x}_1 = f(t, x, \theta)$. We use a nonparametric estimator $\hat{\phi} = (\hat{\phi}_1, \hat{\phi}_2)$ of $\phi^* : [0, 1] \to \mathbb{R}^2$.

Starting from (2.3), a reasonable estimator $\hat{\theta}$ should satisfy the weak formulation

$$ \forall \varphi \in C_C^1([0, 1]), \quad \int_0^1 f(t, \hat{\phi}(t), \hat{\theta}) \varphi(t) dt + \int_0^1 \hat{\phi}_1(t) \dot{\varphi}(t) dt = 0. \qquad (2.5) $$

The vector space $C_C^1([0, 1])$ is not tractable for variational formulation, and one prefers Hilbert space with a structure related to $L^2$. In our case, we use $H_0^1 = \{h \in H^1|h(0) = h(1) = 0\}$ which has a simple description within $L^2$: an orthonormal basis is given by the sine functions $t \mapsto \sqrt{2}\sin(\ell\pi t)$, $\ell \ge 1$ and we have

$$ H_0^1 = \left\{ \sum_{\ell=1}^{\infty} a_{\ell} \sqrt{2} \sin(\ell \pi t) \left| \sum_{\ell=1}^{\infty} \ell^2 a_{\ell}^2 < \infty \right. \right\} \qquad (2.6) $$

Hence, it suffices to consider a countable number of orthogonal conditions (2.5) defined, for instance, with the test functions $\varphi_{\ell} = \sqrt{2}\sin(\ell\pi t)$, $\forall \ell \ge 1$:

$$ C_{\ell}(\theta) : \int_{0}^{1} f(t, \hat{\phi}(t), \hat{\theta}) \varphi_{\ell}(t) dt + \int_{0}^{1} \hat{\phi}(t) \dot{\varphi}_{\ell}(t) dt = 0. \qquad (2.7) $$

More generally, we consider a family of orthonormal functions $\varphi_{\ell} \in H_0^1$, with $\ell \ge 1$, and we introduce
## 4.1 Asymptotic representation for $\hat{\theta}_n - \theta$

2. From the definition (2.9) of $\hat{\theta}_{n,L}$ and differentiability of $f$, the first order optimality condition is

$$ \mathbf{J}_{\theta,L} (\hat{\phi}, \hat{\theta}_{n,L})^{\top} \mathbf{e}_L (\hat{\phi}, \hat{\theta}_{n,L}) = 0 \quad (4.1) $$

3. from which we derive an asymptotic representation for $\hat{\theta}_{n,L}$, by linearizing $\mathbf{e}_L(\hat{\phi}, \hat{\theta}_{n,L})$ around $\theta_L^*$. We

4. need to introduce the matrix-valued function defined on $\mathcal{D} \times \theta$ such that $\mathbf{M}_L(g, \theta) = [\mathbf{J}_{\theta,L}(g, \theta)^{\top} \mathbf{J}_{\theta,L}(g, \theta)]^{-1}$

5. $\mathbf{J}_{\theta,L}(g, \theta)^{\top}$, and proposition 2 shows that $\mathbf{M}_L(\hat{\phi}, \hat{\theta}_{n,L})$ is also a consistent estimator of $\mathbf{M}_L^*$.

6. **Proposition 2.** If conditions C1-C6 are satisfied, then

$$ \left[ \mathbf{J}_{\theta,L} (\hat{\phi}, \hat{\theta}_{n,L})^{\top} \tilde{\mathbf{J}}_L \right]^{-1} \mathbf{J}_{\theta,L} (\hat{\phi}, \hat{\theta}_{n,L})^{\top} \xrightarrow{P} \mathbf{M}_L^* = [\mathbf{J}_{\theta,L}^{*\top} \mathbf{J}_{\theta,L}^{*}]^{-1} \mathbf{J}_{\theta,L}^{*\top} \quad (4.2) $$

7. where The matrix $\tilde{\mathbf{J}}_L$ is the Jacobian $\mathbf{J}_{\theta,L}$ evaluated at a point $\tilde{\theta}$ between $\theta^*$ and $\hat{\theta}_{n,L}$. Moreover, we have

$$ \hat{\theta}_{n,L} - \theta_L^* = -\mathbf{M}_L^* \mathbf{e}_L(\hat{\phi}, \theta^*) + o_P(1). \quad (4.3) $$

## 4.2 Linear differential equations

9. We consider the parametrized linear ODE defined as

$$ \left\{ \begin{aligned} \dot{x}_1 &= a(t, \theta_1)x_1 + b(t, \theta_1)x_2 \\ \dot{x}_2 &= c(t, \theta_2)x_1 + d(t, \theta_2)x_2 \end{aligned} \right. \quad (4.4) $$

10. Since the ODE is linear, conditions IVP(a) and IVP(b) are satisfied as soon as the functions $a(\cdot, \theta^*)$, $b(\cdot, \theta^*)$, $c(\cdot, \theta^*)$, $d(\cdot, \theta^*)$ are in $L^2$, regardless of the initial conditions, [17]. We focus only on the estimation of the parameter $\theta = \theta_1$ involved in the first equation $\dot{x}_1 = a(t, \theta)x_1 + b(t, \theta)x_2$ and we suppose that we have two series estimators $\hat{\phi}_1 = \mathbf{p}_K^{\top}\hat{\mathbf{c}}_1$ and $\hat{\phi}_2 = \mathbf{p}_K^{\top}\hat{\mathbf{c}}_2$ satisfying condition C2. The orthogonal conditions are simple linear functionals of the estimators $e_\ell(\hat{\phi}, \theta) = \langle \hat{\phi}_1, \dot{\varphi}_\ell + a(\cdot, \theta)\varphi_\ell \rangle + \langle \hat{\phi}_2, b(\cdot, \theta)\varphi_\ell \rangle$. Hence the asymptotic behavior of the empirical orthogonal conditions relies on the plug-in properties of $\hat{\phi}_1$ and $\hat{\phi}_2$ into the linear forms $T_\rho: x \mapsto \int_0^1 \rho(t)x(t)dt$ where $\rho$ is a smooth function. Moreover, the
# 3 Real data analysis

## 3.1 Influenza virus growth and migration model

<table><thead><tr><th></th><th>θ̂<sup>OC</sup><sub>3</sub></th><th>θ̂<sup>OC</sup><sub>4</sub></th><th>θ̂<sup>NLS</sup></th><th>θ̂<sup>ref</sup></th></tr></thead><tbody><tr><th>ρ<sub>m</sub></th><td>2.9e-5</td><td>2.7e-5</td><td>1.5e-5</td><td>1.6e-5</td></tr><tr><th>ρ<sub>s</sub></th><td>4.1e-5</td><td>4.7e-5</td><td>4.1e-5</td><td>4.5e-5</td></tr><tr><th>δ<sub>l</sub></th><td>2.0</td><td>3.4</td><td>3.7</td><td>3.96</td></tr><tr><th>γ<sub>ms</sub></th><td>0.39</td><td>0.35</td><td>0.15</td><td>0.157</td></tr><tr><th>γ<sub>sl</sub></th><td>0.72</td><td>0.81</td><td>0.47</td><td>0.49</td></tr><tr><th>RMSE</th><td>13.5</td><td>13.9</td><td>9.0</td><td>9.5</td></tr></tbody></table>

<table><caption>Table 6: Estimates, RMSE and the 95% confidence intervals for different L and estimators.</caption><thead><tr><th rowspan="2"></th><th colspan="2">θ̂<sup>OC</sup><sub>3</sub></th><th colspan="2">θ̂<sup>OC</sup><sub>4</sub></th><th colspan="2">θ̂<sup>NLS</sup></th></tr><tr><th>Low. Bound</th><th>Up. Bound</th><th>Low. Bound</th><th>Up. Bound</th><th>Low. Bound</th><th>Up. Bound</th></tr></thead><tbody><tr><th>ρ<sub>m</sub></th><td>2.1e-5</td><td>3.7e-5</td><td>1.9e-5</td><td>3.4e-5</td><td>0.7e-0.5</td><td>2.4e-0.5</td></tr><tr><th>ρ<sub>s</sub></th><td>0.7e-5</td><td>7.4e-5</td><td>0.9e-5</td><td>8.4e-5</td><td>3.4e-0.5</td><td>4.8e-0.5</td></tr><tr><th>δ<sub>l</sub></th><td>-1.11</td><td>5.21</td><td>-0.28</td><td>7.21</td><td>2.59</td><td>4.93</td></tr><tr><th>γ<sub>ms</sub></th><td>0.27</td><td>0.50</td><td>0.24</td><td>0.46</td><td>0.03</td><td>0.26</td></tr><tr><th>γ<sub>sl</sub></th><td>-0.10</td><td>1.55</td><td>-0.14</td><td>1.76</td><td>0.39</td><td>0.55</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 8: Influenza: Estimated curves for $X_1$ (red), $X_2$ (green), $X_3$ (blue); $\times$: observations, $\square$: solution for $\hat{\theta}_1^{OC}$, $\circ$: solution for $\hat{\theta}_2^{OC}$, solid line: solution with $\hat{\theta}^{NLS}$.</figcaption></figure>
1. the vector space $\mathcal{F} = \overline{\operatorname{span}\{\varphi_\ell, \ell \ge 1\}}$. The vector space $\mathcal{F}$ may not be necessarily dense in $H_0^1$, as
2. the functions $\varphi_\ell$ could be chosen for computational tractability or because of a natural interpretation
3. (for instance B-splines, polynomials, wavelets, ad-hoc functions, ...). For this reason, we introduce the
4. orthogonal decomposition of $H_0^1 = \mathcal{F} \oplus \mathcal{F}^\perp$, where $\mathcal{F}^\perp = \{g \in H_0^1 | \langle g, \varphi \rangle = 0, \varphi \in \mathcal{F}\}$, and we can have
5. $\mathcal{F} \neq H_0^1$. In general, an estimator $\hat{\theta}$ satisfying $C_\ell(\hat{\theta})$ for $\ell \ge 1$ also approximately satisfies (2.5). However
6. in practice, we will use a finite set of orthogonal constraints defined by $L$ test functions ($L > p$).

7. In order to discuss the influence of the choice of $\mathcal{F}$ and of finite dimensional subspace spanned by
8. $\varphi_1, \dots, \varphi_L$ we introduce the nonlinear operator $\mathcal{E}: (g, \theta) \mapsto \mathcal{E}(g, \theta)$, such that $t \mapsto \mathcal{E}(g, \theta)(t) =$
9. $f(t, g(t), \theta)$.

10. For all $\theta$ in $\Theta$ and $g$ in $H^1$, the Fourier coefficients of $\mathcal{E}(g, \theta) - \dot{g}$ in the basis $(\varphi_\ell)_{\ell \ge 1}$ are $e_\ell(g, \theta) =$
11. $\langle \mathcal{E}(g, \theta) - \dot{g}, \varphi_\ell \rangle = \langle \mathcal{E}(g, \theta), \varphi_\ell \rangle + \langle g, \dot{\varphi}_\ell \rangle$, and we introduce the vectors in $\mathbb{R}^L$ $e_L(g, \theta) = (e_\ell(g, \theta))_{\ell=1..L}$ and
12. $e_L^*(\theta) = (e_\ell(\phi^*, \theta))_{\ell=1..L}$. Finally, our estimator is defined by minimizing the quadratic form $Q_{n,L}(\theta) =$
13. $|e_L(\hat{\phi}, \theta)|^2$:

$$ \hat{\theta}_{n,L} = \arg \min_{\theta \in \Theta} Q_{n,L}(\theta). \qquad (2.8) $$

14. $\hat{\theta}_{n,L}$ is the parameter that “almost” vanishes the first $L$ Fourier coefficients in the orthogonal decomposition of $H_0^1 = \mathcal{F} \oplus \mathcal{F}^\perp$:

$$ \mathcal{E}(g, \theta) - \dot{g} = \mathbf{E}_L(g, \theta) + \mathbf{R}_L(g, \theta) + \mathbf{E}_{\mathcal{F}}^{\perp}(g, \theta) $$

16. with $\mathbf{E}_L(g, \theta) = \sum_{\ell=1}^L e_\ell(g, \theta) \varphi_\ell$, $\mathbf{R}_L(g, \theta) = \sum_{\ell>L} e_\ell(g, \theta) \varphi_\ell$ and $\mathbf{E}_{\mathcal{F}}^{\perp}(g, \theta) \in \mathcal{F}^{\perp}$.

17. The function $\mathbf{E}_{\mathcal{F}}^{\perp}(\phi^*, \theta)$ represents the behavior of $\mathcal{E}(g, \theta) - \dot{g}$ at the boundaries of the interval $[0, 1]$. As $\hat{\phi}$ approaches $\phi^*$ asymptotically in supremum norm, the objective function $Q_{n,L}(\theta)$ is close to $Q_L^*(\theta) = \|\mathbf{E}_L(\phi^*, \theta)\|_{L^2}^2$. The discriminative power of $Q_L^*(\theta)$ can be analyzed locally around its global minimum $\theta_L^*$, as it behaves approximately as the quadratic form $Q_L^*(\theta) \approx (\theta - \theta_L^*)^\top \mathbf{J}_{\theta,L}^* \mathbf{J}_{\theta,L} (\theta - \theta_L^*)$ where $\mathbf{J}_{\theta,L}^*$ is the matrix in $\mathbb{R}^{L \times p}$ with entries $\int_0^1 f_{\theta_j}(t, \phi^*(t), \theta_L^*) \varphi_\ell(t) dt$, for $j = 1, \dots, p$, $\ell = 1, \dots, L$.
1 linearity of series estimator makes the orthogonal conditions $\mathbf{e}_L(\hat{\phi}, \theta)$ easy to compute as

$$ \mathbf{e}_L(\hat{\phi}, \theta) = \mathbf{A}(\theta)\hat{\mathbf{c}}_1 + \mathbf{B}(\theta)\hat{\mathbf{c}}_2 \quad (4.5) $$

2 where $\mathbf{A}(\theta)$ and $\mathbf{B}(\theta)$ are matrices in $\mathbb{R}^{L \times K}$ with entries $A_{\ell,k}(\theta) = \int_0^1 (a(t, \theta)\varphi_{\ell}(t) + \dot{\varphi}_{\ell}(t)) p_{kK}(t)dt$

3 and $B_{\ell,k}(\theta) = \int_0^1 (b(t, \theta)\varphi_{\ell}(t)) p_{kK}(t)dt$. The gradient of $\mathbf{e}_L(\hat{\phi}, \theta)$ is $\mathbf{J}_{\theta,L}(\hat{\phi}, \theta) = \partial_{\theta}\mathbf{A}(\theta)\hat{\mathbf{c}}_1 + \partial_{\theta}\mathbf{B}(\theta)\hat{\mathbf{c}}_2$

4 where $\partial_{\theta}\mathbf{A}(\theta)$ and $\partial_{\theta}\mathbf{B}(\theta)$ are straightforwardly computed by permuting differentiation and integration.

5 Although $\mathbf{e}_L(\hat{\phi}, \theta)$ depends linearly on the observations, we have to take care of the asymptotics as we

6 are in a nonparametric framework and $K$ grows with $n$. The behavior of linear functionals $T_{\rho}(\hat{\phi})$ for

7 several nonparametric estimators (kernel regression, series estimators, orthogonal series) is well known

8 [1, 2, 7, 16], and in generality it can be shown that such linear forms can be estimated with the classical

9 root-$n$ rate and that they are asymptotically normal under quite general conditions. In the particular

10 case of series estimators, we rely on theorem 3 of [16] that ensures the root-$n$ consistency and the

11 asymptotic normality of the plugged-in estimators $T_{\rho}(\hat{\phi}_j)$, $j = 1, 2$ under almost minimal conditions.

12 We will give in the next section the precise assumptions required for root-$n$ consistency of linear and

13 nonlinear functional of the series estimator. Moreover, the variance of $\hat{\theta}_{n,L}$ has a remarkable expression

$$ V_{e,L}(\theta) = V(\mathbf{e}_L(\hat{\phi}, \theta)) = \mathbf{A}(\theta)V(\hat{\mathbf{c}}_1)\mathbf{A}(\theta)^{\top} + \mathbf{B}(\theta)V(\hat{\mathbf{c}}_2)\mathbf{B}(\theta)^{\top}. \quad (4.6) $$

14 We remark that there is no covariance term between $\hat{\mathbf{c}}_1$ and $\hat{\mathbf{c}}_2$ since we assume that $V(Y|T=t)$ is

15 diagonal (assumption C2), but in all generality, we should add $2\mathbf{A}(\theta)\mathrm{cov}(\hat{\mathbf{c}}_1, \hat{\mathbf{c}}_2)\mathbf{B}(\theta)^{\top}$. We can use the

16 classical estimates of the variance of $\hat{\mathbf{c}}_1$ and $\hat{\mathbf{c}}_2$ to compute an estimate of $V_{e,L}(\theta)$

$$ \widehat{V_{e,L}(\theta)} = \mathbf{A}(\theta)\widehat{V}(\hat{\mathbf{c}}_1)\mathbf{A}(\theta)^{\top} + \mathbf{B}(\theta)\widehat{V}(\hat{\mathbf{c}}_2)\mathbf{B}(\theta)^{\top} \quad (4.7) $$

17 Thanks to proposition 2, we can estimate the asymptotic variance of the estimator $\hat{\theta}_{n,L}$ with the consistent estimator $\hat{\mathbf{M}}_L = \mathbf{M}_L(\hat{\phi}, \hat{\theta}_{n,L})$ and we estimate $V(\hat{\theta}_{n,L})$ by $V(\widehat{\hat{\theta}}_{n,L}) = \hat{\mathbf{M}}_L V(\widehat{\mathbf{e}}_L(\hat{\phi}, \hat{\theta}_{n,L})) \hat{\mathbf{M}}_L^{\top}$.

18 From the asymptotic normality of the plug-in estimate, we can derive confidence balls or confidence

19
<figure><img src="image_9.png" /><figcaption>Figure 9: Influenza: Estimated curves for $X_1$ (red), $X_2$ (green), $X_3$ (blue); □ solution obtained with OC+NLS, ◦ solution obtained with $\tilde{\theta}^{ref}$.</figcaption></figure>
## 2.3 Boundary Conditions and Construction of Orthogonal Conditions

The construction of the orthogonal conditions $e_\ell(\theta)$ exposed in the previous section is generic and can be proposed for numerous types of Differential Equations, in particular for Ordinary and Delay Differential Equations. Moreover, similar orthogonal conditions could be also derived for solutions of PDEs with a relevant set of test functions $\varphi$, but this extension is beyond the scope of the present paper. A process for deriving "regular" orthogonal conditions, (i.e that gives rise to root-$n$ consistent estimator, as it is shown in section 4) is to use conditions $C_\ell(\theta)$ with an integral expression $\int_0^1 h_\ell(t, \hat{\phi}(t), \theta) dt$. The function $h_\ell: (t, x, \theta) \to \mathbb{R}$ must be smooth and must satisfy the remarkable identity $\int_0^1 h_\ell(t, \phi^*(t), \theta^*) = 0$. The variational formulation generates functions $h_\ell(t, x, \theta) = (f(t, x, \theta) \varphi_\ell(t) - \dot{\varphi}_\ell(t)x)$ whereas the classical Gradient Matching considers a single function $h(t, x, y, \theta) = \|f(t, x, \theta) - y\|^2 \varphi(t)$, and the variable $y$ is evaluated along the derivative $\dot{\phi}(t)$. The asymptotic analysis shows that the dependency in $y$ can be removed and that $h'$ behaves in fact as a function $h(t, x, \theta)$.

The OC framework then generalizes the classical TS estimator and gives ways to ameliorate it. Among other, the use of the boundary vanishing function $\varphi$ implies an information loss close to the boundaries. This loss can be sensible in terms of estimation quality, and should be avoided when the boundary values are known. For instance, for an IVP with known initial condition $\phi(0) = \phi_0$, we can derive an orthogonal condition that takes into account the knowledge of $\phi_0$. By direct computation, we have

$$ \int_0^1 h(t, \phi(t), \theta)dt = \int_0^1 f(t, \phi(t), \theta)\varphi(t)dt - [\phi(1)\varphi(1) - \phi(0)\varphi(0)] \\ + \int_0^1 \phi(t)\dot{\varphi}(t)dt. $$

If $\phi(1)$ is unknown, but $\phi(0)$ is known, it suffices to take $\varphi$ such that $\varphi(1) = 0$ and $\varphi(0) \neq 0$. The orthogonal condition still have the same expression $h(t, x, \theta)$. The same adaptation can be done when boundary values of the derivative are known (called Neumann's condition), for instance $\dot{\phi}(1) = \phi'_1$ is known. Indeed, the ODE gives a relationship between the second order derivative $\ddot{\phi}$ and the state $\phi$, as $\ddot{\phi}(t) = \partial_x f(t, \phi, \theta)f(t, \phi, \theta)$. By choosing $\varphi$ such that $\varphi(0) = 0$ and by Integration By Part, the following
1 with level $1-\alpha$. For instance, for each parameter $\theta_i$, $i = 1, \dots, p$:

$$IC(\theta_i; 1 - \alpha) = \left[ \left( \hat{\theta}_{n,L} \right)_i \pm q_{1-\frac{\alpha}{2}} \hat{V} \left( \hat{\theta}_{n,L} \right)_{ii}^{1/2} \right]$$

2 where $q_{1-\alpha/2}$ is the quantile of order $1 - \frac{\alpha}{2}$ of a standard Gaussian distribution. Nevertheless, we recall

3 that these confidence intervals might be affected by the bias of $\hat{\theta}_{n,L}$ depending on $L$.

4

5 **4.3 Nonlinear differential equations**

6 We give here general results for the asymptotics of $e_\ell(\hat{\phi}, \theta)$ when the functional is linear or not in $\hat{\phi}$.

7 In [16], the root-n consistency and asymptotic normality is obtained if the functional $g \mapsto e_\ell(g, \theta)$ has

8 a continuous Fréchet derivative $De_\ell(g, \theta)$ with respect to the norm $\|\cdot\|_\infty$. If $x \mapsto f(t, x, \theta)$ is twice

9 continuously differentiable for $t \in [0, 1]$ a.e. in $x$ and $\theta$ in $\Theta$, then we can compute easily its Fréchet derivative for $g \in H^1$ in the uniform ball $\|g - \phi^*\|_\infty \le D$. For all $h \in H^1$ such $\|g + h - \phi^*\|_\infty \le D$, we

10 have

$$e_{\ell}(g+h, \theta) - e_{\ell}(g, \theta) = \langle f_x(\cdot, g, \theta)h, \varphi_{\ell} \rangle + \langle h, \dot{\varphi} \rangle + \langle h^{\top} f_{xx}(\cdot, \tilde{g}, \theta)h, \varphi_{\ell} \rangle$$

12 by a Taylor expansion around $g$. As in the linear case, we introduce the tangent linear operator

13 $A_g(\theta) : u \mapsto \dot{u} - a_g(t, \theta)u$ with $a_g(t, \theta) = f_{x_1}(t, g(t), \theta)$ and the function $b_g(t, \theta) = f_{x_2}(t, g(t), \theta)$.

14 For all $\theta$, the Fréchet derivative of $e_\ell(g, \theta)$ (w.r.t to the uniform norm) is the linear operator $h = (h_1, h_2) \mapsto De_\ell(g, \theta).h = \langle h_1, \dot{\varphi}_\ell + a_g(t, \theta)\varphi_\ell \rangle + \langle h_2, b_g(\cdot, \theta)\varphi_\ell \rangle$ and satisfies for all $\theta \in \Theta$

$$|e_{\ell}(g+h, \theta) - e_{\ell}(g, \theta) - De_{\ell}(g, \theta).h| \le C \|h\|_{\infty}^{2}$$

16 because $f_{xx}$ is uniformly dominated on $D \times \Theta$. Moreover, for all $\epsilon$ (with $0 < \epsilon < D$), for all $g, g'$ such

17 that $\|g - \phi^*\|_\infty, \|g' - \phi^*\|_\infty \le \epsilon$, we have

$$\begin{aligned} |De_{\ell}(g, \theta).h - De_{\ell}(g', \theta).h| & \le \int_{0}^{1} h(t)^{\top} f_{xx}(t, \tilde{g}(t), \theta) (g(t) - g'(t)) \varphi_{\ell}(t) dt \\ & \le C \|h\|_{\infty} \|g - g'\|_{\infty} \end{aligned}$$
## 3.2 Blowfly models

<table><thead><tr><th></th><th>L = 11</th><th>L = 9</th><th>L = 12</th></tr></thead><tbody><tr><th>P</th><td>7.81</td><td>7.52</td><td>7.91</td></tr><tr><th>N<sub>0</sub></th><td>381.8</td><td>385.9</td><td>377.7</td></tr><tr><th>δ</th><td>0.154</td><td>0.153</td><td>0.154</td></tr><tr><th>RSSE</th><td>1.7136e+03</td><td>1.7557e+03</td><td>1.7990e+03</td></tr></tbody></table>

<table><caption>Table 7: Estimates, RSSE and 95% confidence intervals for different L</caption><thead><tr><th></th><th colspan="2">L = 11</th><th colspan="2">L = 9</th><th colspan="2">L = 12</th></tr><tr><th>O.C</th><th>Low. Bound</th><th>Up. Bound</th><th>Low. Bound</th><th>Up. Bound</th><th>Low. Bound</th><th>Up. Bound</th></tr></thead><tbody><tr><th>P</th><td>5.80</td><td>9.81</td><td>5.64</td><td>9.40</td><td>5.0416</td><td>10.77</td></tr><tr><th>N<sub>0</sub></th><td>303.62</td><td>459.94</td><td>306.59</td><td>465.38</td><td>289.36</td><td>465.98</td></tr><tr><th>δ</th><td>0.10</td><td>0.20</td><td>0.11</td><td>0.19</td><td>0.10</td><td>0.20</td></tr></tbody></table>
1 identity

$$ \varphi(1)\phi'_{1} = \int_{0}^{1} \partial_{x}f(t, \phi, \theta)f(t, \phi, \theta)\varphi(t)dt + \int_{0}^{1} f(t, \phi, \theta)\dot{\varphi}(t)dt $$

2 gives a new condition that exploits the behavior of the solution at the boundary. Obviously, these 3 conditions can be successfully used if the nonparametric proxy satisfies the boundary conditions of 4 interest. At the contrary, it seems rather difficult to integrate such information about the boundary 5 within the criterion $R_{n,w}(\theta)$. The orthogonal conditions introduced in the previous section are a direct 6 exploitation of the ODE model, and the introduction of the space $\mathcal{F}$ is a way to deal with the problem of 7 the choice of the number of conditions and their type. Nevertheless, it would be useful to introduce model 8 specific conditions $h(t, \phi(t), \theta)$ which are known to have a vanishing integral for $\theta = \theta^*$. Our estimator 9 can be thought as a Generalized Method of Moments estimator, but where Moments do characterize 10 curves and not probability distributions. A similar idea has been developed recently in the context of 11 functional data analysis [18].

## 3 Consistency of the Orthogonal Conditions estimator

12 In order to obtain precise results with closed-form expression for the bias and variance estimators, 13 we consider series estimators, i.e. estimators expressed as $\hat{\phi}_j = \sum_{k=1}^{K} \hat{c}_{k,j} p_{kK} = \hat{c}_j p^K$, where $p^K = 15 (p_{1K}, \dots, p_{kK})$ is a vector of approximating functions and the coefficients $\hat{c}_j = (\hat{c}_{k,j})_{k=1..K}$ are computed 16 by least squares. For notational simplicity, we use the same functions (and the same number K) for 17 estimating $\phi_1^*$ and $\phi_2^*$. We denote $P^K = (p_{kK}(t_i))_{1\le i,k \le n,K}$ the design matrix and $\mathbf{Y}_j = (y_{i,j})_{i=1..n}$ the 18 vectors of observations. Hence, the estimated coefficients $\hat{c}_j = (P^{K\top} P^K)^\dagger P^{K\top} \mathbf{Y}_j$ (where $\dagger$ denotes a 19 generalized inverse) gives rise to the so-called hat matrix $H = P^K (P^{K\top} P^K)^\dagger P^{K\top}$ and the vector of 20 smoothed observations is $\hat{\phi}_j = H \mathbf{Y}_j$, $j = 1, 2$. One can typically think of regression splines, [32]. We 21 introduce now the conditions required for the definition and consistency of our estimator.

22 **Condition C1:** (a) $\Theta$ is a compact set of $\mathbb{R}^p$ and $\theta^*$ is an interior point of $\Theta$, $\mathcal{X}$ is an open subset of $\mathbb{R}^2$; (b) $(t, x) \mapsto f(t, x, \theta^*)$ is $L^2$-Lipschitz and $L^2$-Caratheodory (see Supplementary Material I, section 1).
1 with $C$, a constant independent of $\theta$, $\epsilon$ and $g, g'$ (because $f_{xx}$ is uniformly dominated).

2 As in the linear case, we need to evaluate $De_{\ell}(g, \theta)$ on the basis $\mathbf{p}^K$. We denote $\mathbf{A}(g, \theta)$ and $\mathbf{B}(g, \theta)$

3 the matrices in $\mathbb{R}^{L \times K}$ with entries $\int_0^1 a_g(t, \theta) \varphi_{\ell}(t) p_{kK} dt$ and $\int_0^1 b_g(t, \theta) \varphi_{\ell}(t) p_{kK} dt$ (respectively) and we

4 have the approximation

$$ \mathbf{e}_L(\hat{\phi}, \theta) = \mathbf{e}_L(\phi^*, \theta) + \mathbf{A}(\phi^*, \theta)\hat{\mathbf{c}}_1 + \mathbf{B}(\phi^*, \theta)\hat{\mathbf{c}}_2 + O(\|h\|_{\infty}^2). \quad (4.8) $$

5 We can derive the asymptotic variance of $\mathbf{e}_L(\hat{\phi}, \theta)$ from (4.8)

$$ V_{e,L}(\theta) = \mathbf{A}(\phi^*, \theta)V(\hat{\mathbf{c}}_1)\mathbf{A}(\phi^*, \theta)^{\top} + \mathbf{B}(\phi^*, \theta)V(\hat{\mathbf{c}}_2)\mathbf{B}(\phi^*, \theta)^{\top} \quad (4.9) $$

6 and we can get an estimate $\hat{V}_{e,L}(\theta)$ from the data as in the linear case. We need the two additional

7 conditions for deriving the root-$n$ rate, that comes from [16] (Condition C8 is essential for root-$n$ rate):

8 **Condition C7 (a)** The times $T_1, \dots, T_n$ have a density $\pi$ w.r.t. Lebesgue measure such $0 < c < \pi <$

9 $C < \infty$; (b) $E[\epsilon^4] < \infty$.

10 **Condition C8** There exists $\tilde{\beta}_K$ in $\mathbb{R}^K$ with $E_T\left[\left|f_x(T, \phi^*(T), \theta) \varphi_{\ell}(T) - \tilde{\beta}_K^\top p^K(T)\right|^2\right] \to 0$.

11 **Theorem 3.** Under conditions C1 to C8 and if $f$ is a linear ODE or, $f$ is a nonlinear ODE and

12 $\frac{\zeta_0(K)^4 K^2}{n} \to 0$ then for all $\theta \in \Theta$, $\sqrt{n}(\mathbf{e}_L(\hat{\phi}, \theta) - \mathbf{e}_L(\phi^*, \theta)) \rightsquigarrow N(0, V_{e,L}(\theta))$ with $V_{e,L}(\theta)$ given by

13 (4.9), and we denote $\mathbf{V}_{e,L}^* = V_{e,L}(\theta^*)$. Moreover, $\hat{\theta}_{n,L}$ is such that

$$ \sqrt{n}(\hat{\theta}_{n,L} - \theta_L^*) \rightsquigarrow N(0, \mathbf{V}_L^*) \quad (4.10) $$

14 with

$$ \mathbf{V}_L^* = \mathbf{M}_L^* \mathbf{V}_{e,L}^* \mathbf{M}_L^{*\top}. \quad (4.11) $$

15 The asymptotic variance can be estimated as $\hat{\mathbf{M}}_L V_{e,L}(\hat{\theta}_{n,L}) \hat{\mathbf{M}}_L^\top \xrightarrow{P} \mathbf{V}_L^*$. In particular, if we use re-

16 gression splines and $t \mapsto f(t, \phi^*(t), \theta)$ is $C^s$ on $[0, 1]$ with $s \ge 3$, then (4.10) holds with $K$ such that

17 $\sqrt{n}K^{-s} \to 0$ and $n^{-1}K^4 \to 0$.
1 **Condition C2:** (a) $(Y_i, t_i)$ are i.i.d. with variance $V(Y|T=t) = \Sigma_\epsilon = \sigma^2 I_2$; (b) For every $K$, there
2 is a nonsingular constant matrix $B$ such that for $P^K = B_p^K(t)$; (i) the smallest eigenvalue of
3 $E[P^K(T)P^K(T)^\top]$ is bounded away from zero uniformly in $K$ and (ii) there is a sequence of
4 constants $\zeta_0(K)$ satisfying $\sup_t |P^K(t)| \le \zeta_0(K)$ and $K = K(n)$ such that $\zeta_0(K)^2 K/n \to 0$ as
5 $n \to \infty$; (c) There are $\alpha, c_{1,K}, c_{2,K}$ such that $\|\phi_j^* - p^K c_{j,K}\|_\infty = \sup_{[0,1]} |\phi_j^*(t) - p^K(t)^\top c_{j,K}| = O(K^{-\alpha})$.

7 **Condition C3:** There exists $D > 0$, such that the $D$-neighborhood of the solution range $\mathcal{D} = \{x \in \mathbb{R}^2 | \exists t \in [0, 1], |x - \phi^*(t)| < D\}$ is included in $\mathcal{X}$ and $f$ is $C^2$ in $(x, \theta)$ on $\mathcal{D} \times \Theta$ for $t$ in $[0, 1]$ a.e.
9 Moreover, the derivatives of $f$ w.r.t $x$ and $\theta$ (with obvious notations) $f_x, f_\theta, f_{xx}, f_{x\theta}$ and $f_{\theta\theta}$ are
10 $L^2$ uniformly bounded on $\mathcal{D} \times \Theta$ by $L^2$ functions $\bar{h}_x, \bar{h}_\theta, \bar{h}_{x\theta}, \bar{h}_{xx}$ and $\bar{h}_{\theta\theta}$ (respectively).

11 **Condition C4:** Let $(\varphi_\ell)_{\ell \ge 1}$ be an orthonormal sequence of $C^1$ functions in $H_0^1$.

12 **Condition C5:** $\theta^*$ is the unique global minimizer of $Q_F^*$ and $\inf_{|\theta-\theta^*|>\epsilon} Q_F^*(\theta) > 0$.

13 **Condition C6:** There exists $L_0$ such that for $L \ge L_0$, $J_{\theta,L}(g, \theta)$ is full rank in a neighborhood of $(\phi^*, \theta^*)$.
14

15 Condition **C1** gives the existence and uniqueness of a solution $\phi^*$ in $H^1$ to the IVP for $\theta = \theta^*$ and
16 $x(0) = \phi^*(0)$. If $f$ is continuous in $t$ and $x$, then the derivative $\dot{\phi}^*(t) = f(t, \phi^*(t), \theta^*)$ can be defined on
17 $]0, 1[$ and is also continuous, see appendix A. More generally, **C1** does apply when there is a discontin-
18uous input variable, such as in the Ricatti example described in section 1.2.1.

19 Under condition **C2** (satisfied among others by regression splines with $\zeta_0(K) = \sqrt{K}$), it is known
20 that the series estimator $\hat{\phi}_j$ are consistent estimators of $\phi_j^*$ for usual norms, in particular $\|\hat{\phi}_j - \phi_j^*\|_\infty =$
21 $O_P\left(\zeta_0(K)\left(\sqrt{K/n} + K^{-\alpha}\right)\right)$ (theorem 1, [24]). If $\phi^*$ is $C^s$ and we use splines then $\alpha = s$ and $\|\hat{\phi} - \phi^*\|_\infty =$
22 $O_P\left(K/\sqrt{n} + K^{1/2-s}\right)$.

24 Condition **C3** is here to control the continuity and regularity of the function $\mathcal{E}$ involved in the inverse problem. Moreover, it provides uniform control needed for stochastic convergence.

27
1. Theorem 3 is a direct application of theorem 3 in [16] that claims the root-n consistency and asymptotic normality of general (nonlinear) plug-in estimators.

3.4.4 From optimal weighting to a practical estimation algorithm

4. Theorem (3) is of practical interest as it provides a closed-form expression for the asymptotic variance of $\hat{\theta}_{n,L}$ for general ODE. We can make a parallel between the orthogonal condition estimator and the weighted nonlinear least-squares:

7. where $w(\cdot)$ is a (positive) weight function and $\theta_c = (x_0, \theta)$. $\hat{\theta}_c^{WLS}$ is consistent and $\sqrt{n}(\hat{\theta}_c^{WLS} - \theta_c^*) \rightsquigarrow N(0, V^{AWLS})$ under classical regularity assumptions on $f$, see [19]. The asymptotic variance $V^{AWLS}$ is directly computed from the sensitivity equations, and the optimal weight function $w(\cdot)$ is proportional to the variance function $\sigma^2(\cdot)$, meaning that the unweighted least-squares estimator is optimal in the homoscedastic case. It is clear that $V_L^* \neq V^{AWLS}$, which means that we are sure that the orthogonal condition estimator is not an efficient estimator (at least in the Gaussian case, because NLS and MLE are the same in that case). A striking difference between $V_L^*$ and $V^{AWLS}$ is that the least squares involves the Jacobian of the solution w.r.t. the initial values and parameters, whereas the orthogonal conditions involve the Jacobian of the vector field. It is then hard to compare these two matrices in generality and to evaluate the loss of efficiency. Nevertheless, we can compare the influence of a weight matrix for the minimization of the orthogonal conditions. Indeed, if we introduce a positive definite matrix $W$ in $\mathbb{R}^{L \times L}$, we can define the weighted criterion

19. and the corresponding estimator

$$ \hat{\theta}_{n,L}^W = \arg \min_{\theta \in \Theta} Q_{n,L}^W(\theta). $$
Condition C4 is a sufficient condition for deriving independent conditions $C_\ell(\theta)$, and normalization is useful only to avoid giving implicitly more weight to a condition w.r.t. the other conditions.

Condition C5 means that $\theta^*$ is a global and isolated minima of $Q_{\mathcal{F}}^*(\theta)$, which is standard in M-estimation [37], but can be hard to check in practice. Indeed, the parametric identifiability of ODE models can be hard to show, even for small systems. No general and practical results do exist for assessing the identifiability of an ODE model [21]: it is useful to discriminate between ODE identifiability, statistical identifiability and practical identifiability. The latter being the most useful but almost impossible to check a priori. The essential meaning of condition C5 is that the addition of more and more orthogonal conditions should lead to a perfect and univocal estimation of the true parameter. From our experience and by numerical computations, we can check that $Q_L^*(\theta)$ has a unique minima in $\theta^*$ in a region of interest, for $L$ big enough (usually $L \simeq 2 \times p$). The natural criterion for estimating $\theta$ and for identifiability analysis is

$$Q^*(\theta) = \|\mathcal{E}(\phi^*, \theta) - \mathcal{E}(\phi^*, \theta^*)\|_{L^2}^2$$

but $\|\mathcal{E}_{\mathcal{F}}^{\perp}(\phi^*, \theta)\|_{L^2}^2$ is withdrawn and we use the quadratic form $Q_{\mathcal{F}}^*(\theta)$ in order to avoid boundary effects. This is needed in order to get a parametric rate of convergence, as in the original two-step criterion (2.4). As a consequence, we lose a piece of information brought by the trajectory $t \mapsto \phi^*(t)$ and we have to be sure that the parameter $\theta$ has a low influence on $\|\mathcal{E}_{\mathcal{F}}^{\perp}(\phi^*, \theta)\|_{L^2}^2$. A favorable case is that it is almost constant on $\Theta$, so that $Q^*$ and $Q_{\mathcal{F}}^*$ are essentially the same functions, with the same global minimum and the same discriminating power. In practice, we can check that C5 is approximately satisfied by computing numerically the criterion $\mathcal{E}_{L'}(\phi(\cdot, \hat{\theta}_{n,L}), \theta)$, in a neighborhood of $\hat{\theta}_{n,L}$, for $L' \ge L$.

Finally, Condition C6 is about the influence of the number of test functions used. We use only the first $L$ Fourier coefficients of $\mathcal{E}(g, \theta) - \dot{g}$ to identify the parameter $\theta$, but this might not be sufficient to discriminate between two parameters $\theta$ and $\theta'$. In a way, we perform dimension reduction but we need to be sure that we have an exact recovery when $L$ goes to infinity: we expect that the global minimum $\theta_L^*$ of $|e_L^*(\theta)|^2$ is close to the global minimum $\theta^*$ of $Q_{\mathcal{F}}^*(\theta) = \|\mathcal{E}_{\mathcal{F}}(\phi^*, \theta)\|_{L^2}^2$ (found under condition C5). We introduce the Jacobian matrices $\mathbf{J}_{\theta,L}(g, \theta)$ in $\mathbb{R}^{L \times p}$ with entries $\int_0^1 f_{\theta_j}(t, g(t), \theta)\varphi_\ell(t)dt$ and $\mathbf{J}_{x,L}(g, \theta)$ in $\mathbb{R}^{L \times d}$ with entries $\int_0^1 f_{x_i}(t, g(t), \theta)\varphi_\ell(t)dt$. For this reason, we suppose that $\mathbf{J}_{\theta,L}^*$ is full rank, so that $Q_L^*(\theta)$ is locally strictly convex, with a unique local minimum $\theta_L^*$.
The results of theorems 1 and 3 are then still true under straightforward adaptations. In particular, $\hat{\theta}_{n,L}^W$ is consistent and asymptotically normal with asymptotic variance $M_{W,L}^*V_{e,L}^*M_{W,L}^{*\top}$ where $M_{W,L}^* = [J_L^{*\top}WJ_L^*]^{-1}J_L^{*\top}W$. Eventually, we can then ask for the best weighting matrix $W$ giving the smallest asymptotic variance. This is a classical result for Generalized Moment Estimators that we recall in the following proposition (see for instance section 3.6 in [9]):

### Proposition 4. Optimal weighting matrix

The minimal asymptotic variance for $\hat{\theta}_{n,L}^W$ is obtained with $W^{opt} = V_{e,L}^{*-1}$ and

$$V_{L}^{opt}(\theta^{*}) = (J_{\theta,L}^{*\top}V_{e,L}^{*-1}J_{\theta,L}^{*})^{-1} \quad (4.12)$$

Even if we have an homoscedastic model, we have an interest in using a weighted estimator. Interestingly, the problem of choosing the best weighting matrix is directly related to the choice of the best set of test functions $\varphi_1, \dots, \varphi_L$. Indeed, the diagonalization of $V_{e,L}^* = U\Lambda U^\top$ permits the introduction of the eigenvalues $\Lambda = \text{diag}(\lambda_1, \dots, \lambda_L)$ (with $\lambda_1 \ge \dots \ge \lambda_L \ge 0$) and $U = (u_1|\dots|u_L)$ is an orthogonal matrix. If $\Lambda$ is nonsingular, $Q_{n,L}^{opt}(\theta) \doteq Q_{n,L}^{W^{opt}}(\theta) = \tilde{e}_L(\hat{\phi}, \theta)^\top \Lambda^{-1} \tilde{e}_L(\hat{\phi}, \theta)$, with $\tilde{e}_L(\hat{\phi}, \theta) = U^\top e_L(\hat{\phi}, \theta)$. By linearity, the new orthogonal conditions $\tilde{e}_L(\hat{\phi}, \theta)$ can be written as orthogonal conditions with the test functions $\psi_1, \dots, \psi_L$ derived from the eigenvectors of $V_{e,L}^*$ as $\psi_\ell = \sum_{k=1}^L u_{k\ell} \varphi_k$. The use of the optimal weighting matrix is then equivalent to choose the best test functions in $\mathcal{F}_L$ and the diagonalization shows that some care must be taken (typically $L$ have to be reasonably small to avoid $\lambda_L \approx 0$).

Finally, as the optimal weight $W^{opt}$ depends on the true unknown parameter, the practical use of the previous analysis is hard to apply. We suggest to approximate $V_{e,L}^{*-1}$ by computing a sequence of weighted estimates, as described in Algorithm 1.

This is known as iterated Generalized Method of Moments (GMM), and consists in successive reweighting, by using the consistent estimator of the variance $\tilde{V}_{e,L}(\hat{\theta}_{n,L})$, [9]. There is no theoretical guarantee for the (numerical) convergence of this algorithm but if the model is correct, this sequence of iterations tends to ameliorate the quality of the GMM estimator. Possibly, during this stage, it could be necessary to select appropriately $L$, as $\tilde{V}_{e,L}(\hat{\theta}_{n,L}^{(k)})$ can be close to be nonsingular at an iteration $k$.
1 The Jacobian matrix introduced in condition **C6** is classical in sensitivity analysis (in ODE models). Usu-
2 ally, the sensitivity matrix used is the Jacobian of the least squares criterion (similar to $\mathbf{J}_{\theta,L}(\phi(\cdot, \hat{\theta}), \theta)$);
3 it enables to check a posteriori the identifiability of the parameter $\theta$. Conversely, local non-identifiable
4 parameter (sloppy parameters, [15]) can be detected in that case.

5 Theorem 3.1. If conditions C1 to C6 are satisfied, then

$$
\hat{\theta}_{n,L} - \theta_{L}^{*} = O_{P}(1)
$$

6 and the bias $\mathbf{B}_L = \theta_L^* - \theta^*$ tends to zero as $L \to \infty$.

7 In particular, if we use the sine basis and if $\mathcal{E}(\phi^*, \theta)$ is in $H^1$ for all $\theta$, then $\mathbf{B}_L = o(\frac{1}{L})$.

8 **Remark 3.1.** The convergence rate of the bias $\mathbf{B}_L$ can be refined according to the test functions $\phi_l$: if
9 we use B-splines, the bias is controlled by the meshsize $\Delta = \max_{j>1}(\tau_j - \tau_{j-1})$ of the sequence of knots
10 $\tau_j$, $j = 1, \dots, L$ defining the spline spaces, see section 6 in [34].

11 **Remark 3.2.** In practice, we have $\mathbf{B}_L = 0$ for medium-size $L$, around $2 \times d \times p$.

12 **4 Asymptotics**

13 We give a precise description of the asymptotics of $\hat{\theta}_{n,L}$ (rate, variance and normality), by exploiting the
14 well-known properties of series estimators. We consider the linear case, then we extend the obtained
15 results to general nonlinear ODEs. We show in a preliminary step that the asymptotics of $\hat{\theta}_{n,L} - \theta_L^*$ are
16 directly related to the behavior of $\mathbf{e}_L(\hat{\phi}, \theta^*)$, which is a classical feature of Moment Estimators.

17 **4.1 Asymptotic representation for** $\hat{\theta}_n - \theta_L^*$

18 From the definition (2.8) of $\hat{\theta}_{n,L}$ and differentiability of $f$, the first order optimality condition is

$$
\mathbf{J}_{\theta,L} (\hat{\phi}, \hat{\theta}_{n,L})^{\top} \mathbf{e}_L (\hat{\phi}, \hat{\theta}_{n,L}) = 0 \quad (4.1)
$$

19 from which we derive an asymptotic representation for $\hat{\theta}_{n,L}$, by linearizing $\mathbf{e}_L(\hat{\phi}, \hat{\theta}_{n,L})$ around $\theta_L^*$. We
20 need to introduce the matrix-valued function defined on $\mathcal{D} \times \theta$ such that $\mathbf{M}_L(g, \theta) = [\mathbf{J}_{\theta,L}(g, \theta)^{\top} \mathbf{J}_{\theta,L}(g, \theta)]^{-1}$
<caption>**Algorithm 1** Iteratively Reweighted Orthogonal Conditions</caption>

**Require:** $\hat{\phi}, \theta_{init} \in \Theta, \epsilon > 0$, 
Compute the unweighted estimator

$$ \hat{\theta}_{n,L}^{(0)} = \arg \min_{\theta} Q_{n,L}(\theta). $$

Compute the asymptotic covariance of $e_L(\hat{\phi}, \hat{\theta}_{n,L}^{(0)})$ with equation (4.7).

**while** $|\hat{\theta}_{n,L}^{(k)} - \hat{\theta}_{n,L}^{(k-1)}| > \epsilon$ **do**

$$ \left\{ \begin{aligned} W^{(k)} &= \left( V_{e,L}(\hat{\theta}_{n,L}^{(k)}) \right)^{-1} \\ \hat{\theta}_{n,L}^{(k+1)} &= \arg \min_{\theta} Q_{n,L}^{W(k)}(\theta) \end{aligned} \right. . \qquad (4.13) $$

**end while**

## 1 5 Experiments

We apply the Orthogonal Condition (OC) estimators to first order ODEs linear and nonlinear in state.

We consider situations where the Orthogonal Conditions estimator has a closed-form, or when the

classical NLS estimator or Two-Step estimator may collapse. We compare systematically the NLS

estimator $\hat{\theta}^{NLS}$ and Two-Step Estimator (TS) $\hat{\theta}^{TS}$ to the OC estimator $\hat{\theta}_{n,L}$ and the weighted OC

estimator $\hat{\theta}_{n,L}^W$ for varying sample sizes ($n = 400, 200, 50$) and varying noise level (high and small). This

gives a reasonable picture of the robustness, sensitivity and efficiency of the different estimators. For

the NLS estimator, we compute a one-step ameliorated estimator by using $\hat{\theta}_{n,L}^W$ as a starting point. The

estimator obtained should be efficient and should give an indication about the best attainable variance,

and the efficiency loss of the Gradient Matching methods.

Indeed, we consider a Gaussian and homoscedastic noise, such that the NLS should be nearly effi-

cient. This optimality is asymptotic and the case $n = 200$ or $n = 400$ should give indications on the

efficiency loss, whereas the small size case ($n = 50$) gives some relevant information on the complex

small sample case, where the asymptotic approximation cannot be assessed. In particular, the NLS does

not necessarily provide the best estimator.

The nonparametric estimate for the gradient matching estimators (TS and OC) is the same regression

spline, decomposed on B-splines with a uniform knot sequence $\xi_k, k = 1, \dots, K$. For each dataset (and

each dimension), the number of knots is selected by minimizing the GCV criterion, [22]. For the plain
1. $\mathbf{J}_{\theta,L} (g, \theta)^{\top}$, and proposition 4.1 shows that $\mathbf{M}_L(\hat{\phi}, \hat{\theta}_{n,L})$ is also a consistent estimator of $\mathbf{M}_L^*$.

2. **Proposition 4.1.** If conditions C1-C6 are satisfied, then

$$ \left[ \mathbf{J}_{\theta,L} (\hat{\phi}, \hat{\theta}_{n,L})^{\top} \tilde{\mathbf{J}}_L \right]^{-1} \mathbf{J}_{\theta,L} (\hat{\phi}, \hat{\theta}_{n,L})^{\top} \xrightarrow{P} \mathbf{M}_L^* = \left[ \mathbf{J}_{\theta,L}^{*\top} \mathbf{J}_{\theta,L}^{*} \right]^{-1} \mathbf{J}_{\theta,L}^{*\top} \quad (4.2) $$

3. where the matrix $\tilde{\mathbf{J}}_L$ is the Jacobian $\mathbf{J}_{\theta,L}$ evaluated at a point $\tilde{\theta}$ between $\theta^*$ and $\hat{\theta}_{n,L}$. Moreover, we have

$$ \hat{\theta}_{n,L} - \theta_L^* = -\mathbf{M}_L^* \mathbf{e}_L(\hat{\phi}, \theta_L^*) + o_P(1). \quad (4.3) $$

4. **4.2 Linear differential equations**

5. We consider the parametrized linear ODE defined as

$$ \begin{cases} \dot{x}_1 = a(t, \theta_1)x_1 + b(t, \theta_1)x_2 \\ \dot{x}_2 = c(t, \theta_2)x_1 + d(t, \theta_2)x_2 \end{cases} \quad (4.4) $$

6. where $a(\cdot, \theta)$, $b(\cdot, \theta)$, $c(\cdot, \theta)$, $d(\cdot, \theta)$ are in $L^2$. We focus only on the estimation of the parameter $\theta = \theta_1$

7. involved in the first equation $\dot{x}_1 = a(t, \theta)x_1 + b(t, \theta)x_2$ and we suppose that we have two series estimators

8. $\hat{\phi}_1 = \mathbf{p}_K^{\top} \hat{\mathbf{c}}_1$ and $\hat{\phi}_2 = \mathbf{p}_K^{\top} \hat{\mathbf{c}}_2$ satisfying condition C2. The orthogonal conditions are simple linear

9. functionals of the estimators $e_\ell(\hat{\phi}, \theta) = \langle \hat{\phi}_1, \dot{\varphi}_\ell + a(\cdot, \theta)\varphi_\ell \rangle + \langle \hat{\phi}_2, b(\cdot, \theta)\varphi_\ell \rangle$. Hence the asymptotic

10. behavior of the empirical orthogonal conditions relies on the plug-in properties of $\hat{\phi}_1$ and $\hat{\phi}_2$ into the

11. linear forms $T_\rho : x \mapsto \int_0^1 \rho(t)x(t)dt$ where $\rho$ is a smooth function. Moreover, the linearity of series

12. estimator makes the orthogonal conditions $\mathbf{e}_L(\hat{\phi}, \theta)$ easy to compute as

$$ \mathbf{e}_L(\hat{\phi}, \theta) = \mathbf{A}(\theta)\hat{\mathbf{c}}_1 + \mathbf{B}(\theta)\hat{\mathbf{c}}_2 \quad (4.5) $$

13. where $\mathbf{A}(\theta)$ and $\mathbf{B}(\theta)$ are matrices in $\mathbb{R}^{L \times K}$ with entries $A_{\ell,k}(\theta) = \int_0^1 (a(t, \theta)\varphi_\ell(t) + \dot{\varphi}_\ell(t)) p_{kK}(t)dt$

14. and $B_{\ell,k}(\theta) = \int_0^1 (b(t, \theta)\varphi_\ell(t)) p_{kK}(t)dt$. The gradient of $\mathbf{e}_L(\hat{\phi}, \theta)$ is $\mathbf{J}_{\theta,L}(\hat{\phi}, \theta) = \partial_\theta \mathbf{A}(\theta)\hat{\mathbf{c}}_1 + \partial_\theta \mathbf{B}(\theta)\hat{\mathbf{c}}_2$

15. where $\partial_\theta \mathbf{A}(\theta)$ and $\partial_\theta \mathbf{B}(\theta)$ are straightforwardly computed by permuting differentiation and integration.

16. Although $\mathbf{e}_L(\hat{\phi}, \theta)$ depends linearly on the observations, we have to take care of the asymptotics as we
TS estimator, we use the same weight function as in [4], and we consider two variations for the OC estimator depending on the weight matrix: a uniform weight matrix $W = I_L$, and the optimal weighting matrix $W^{opt}$ computed as described in Algorithm 1. The Orthogonal Conditions are defined with the sine basis and the default number of conditions is $L = 2 \times p \times d$.

Finally, we compute the bias and variance of the different estimators, the confidence intervals and their coverage probabilities by Monte Carlo simulations (based on $N_{MC} = 100$ independent draws).

## 5.1 Linear ODE

The classical linear constant ODE can be written as $\dot{x} = Ax + h(t)$ with constant coefficients $A$ with $A^\top = (A_1|\dots|A_d)$ and $h = (h_1(t), \dots, h_d(t))$. For each $i$, $\dot{x}_i = \sum_{k=1}^d a_{ik}x_k + h_i = A_i^\top x + h_i$. The weak form gives rise to a linear matrix equality for $i = 1, \dots, d$ as $Y_i^\varphi = X^\varphi A_i + H_i^\varphi$: the matrix involved is $X^\varphi$ is $d \times L$ with entries $\langle x_k, \varphi_\ell \rangle$, $Y_i^\varphi$ and the vectors $H_i^\varphi$ are in $\mathbb{R}^L$ with respectively entries equal to $-\langle x_i, \dot{\varphi}_\ell \rangle$ and $\langle h_i, \varphi_\ell \rangle$. For illustration, we consider the $\alpha$-pinene (biochemical) network used in [15] for the comparison of several global optimization algorithms:

$$ \left\{ \begin{aligned} \dot{x}_1 &= -(\theta_1 + \theta_2)x_1 \\ \dot{x}_2 &= \theta_1 x_1 \\ \dot{x}_3 &= \theta_2 x_1 - (\theta_3 + \theta_4)x_3 + \theta_6 x_5 \\ \dot{x}_4 &= \theta_3 x_3 \\ \dot{x}_5 &= \theta_4 x_3 + \theta_5 x_5 \end{aligned} \right. \qquad (5.1) $$

The initial condition is $x_0 = (100, 0, 0, 0, 0)$ and is considered unknown. In that case, the system converges to a stationary point. The global quality can be analyzed from the Mean Square Errors in table (5.1): with large sample case and small noise the estimators are equivalent (with error around 0.1), but in the small sample case ($n = 50$), the best estimators are TS and OCopt, and the worst is NLS. NLS and OC are roughly equivalent, and from this example, it is clear that there is a benefit in using relevant weight $W$ in the OC conditions. Despite the truncation, the optimal weighting enables to recover some information with respect to TS, as it identifies the most relevant directions (conditions) to use for measuring the model discrepancy. We do not detail the mean and variance for each parameter,
1 are in a nonparametric framework and K grows with n. The behavior of linear functionals $T_\rho(\hat{\phi})$ for
2 several nonparametric estimators (kernel regression, series estimators, orthogonal series) is well known
3 $[1, 5, 13, 24]$, and in generality it can be shown that such linear forms can be estimated with the classical
4 root-n rate and that they are asymptotically normal under quite general conditions. In the particular
5 case of series estimators, we rely on theorem 3 of [24] that ensures the root-n consistency and the
6 asymptotic normality of the plugged-in estimators $T_\rho(\hat{\phi}_j)$, $j = 1, 2$ under almost minimal conditions.
7 We will give in the next section the precise assumptions required for root-n consistency of linear and
8 nonlinear functional of the series estimator. Moreover, the variance of $\hat{\theta}_{n,L}$ has a remarkable expression

$$V_{e,L}(\theta) = V(\mathbf{e}_L(\hat{\phi}, \theta)) = \mathbf{A}(\theta)V(\hat{\mathbf{c}}_1)\mathbf{A}(\theta)^{\top} + \mathbf{B}(\theta)V(\hat{\mathbf{c}}_2)\mathbf{B}(\theta)^{\top}. \quad (4.6)$$

9 We remark that there is no covariance term between $\hat{\mathbf{c}}_1$ and $\hat{\mathbf{c}}_2$ since we assume that $V(Y|T=t)$ is
10 diagonal (assumption C2), but in all generality, we should add $2\mathbf{A}(\theta)\text{cov}(\hat{\mathbf{c}}_1, \hat{\mathbf{c}}_2)\mathbf{B}(\theta)^{\top}$. We can use the
11 classical estimates of the variance of $\hat{\mathbf{c}}_1$ and $\hat{\mathbf{c}}_2$ to compute an estimate of $V_{e,L}(\theta)$

$$\widehat{V_{e,L}(\theta)} = \mathbf{A}(\theta)\widehat{V(\hat{\mathbf{c}}_1)}\mathbf{A}(\theta)^{\top} + \mathbf{B}(\theta)\widehat{V(\hat{\mathbf{c}}_2)}\mathbf{B}(\theta)^{\top} \quad (4.7)$$

12 Thanks to proposition 4.1, we can estimate the asymptotic variance of the estimator $\hat{\theta}_{n,L}$ with the con-
13 sistent estimator $\hat{\mathbf{M}}_L = \mathbf{M}_L(\hat{\phi}, \hat{\theta}_{n,L})$ and we estimate $V(\hat{\theta}_{n,L})$ by $V(\widehat{\hat{\theta}_{n,L}}) = \hat{\mathbf{M}}_L V(\widehat{\mathbf{e}}_L(\hat{\phi}, \hat{\theta}_{n,L})) \hat{\mathbf{M}}_L^{\top}$.
14 From the asymptotic normality of the plug-in estimate, we can derive confidence balls with level $1 - \alpha$.
15 For instance, for each parameter $\theta_i$, $i = 1, \dots, p$:

$$IC(\theta_i; 1 - \alpha) = \left[ \left( \hat{\theta}_{n,L} \right)_i \pm q_{1-\frac{\alpha}{2}} \widehat{V(\hat{\theta}_{n,L})}_{ii}^{1/2} \right]$$

16 where $q_{1-\alpha/2}$ is the quantile of order $1 - \frac{\alpha}{2}$ of a standard Gaussian distribution. Nevertheless, we recall
17 that these confidence intervals might be affected by the bias of $\hat{\theta}_{n,L}$ depending on L.

18
<table><caption>Table 5.1: Mean Square Error for $\theta$ for Alpha-Pinene ODE (5.1)</caption><thead><tr><th>(n; σ)</th><th>θ̂<sub>TS</sub></th><th>θ̂<sub>OC</sub></th><th>θ̂<sub>OCopt</sub></th><th>θ̂<sub>NLS</sub></th></tr></thead><tbody><tr><td>(400; 50)</td><td>0.1055</td><td>0.0851</td><td>0.0907</td><td>0.1171</td></tr><tr><td>(400; 100)</td><td>0.1325</td><td>0.2170</td><td>0.1372</td><td>0.2443</td></tr><tr><td>(200; 50)</td><td>0.1139</td><td>0.2021</td><td>0.1183</td><td>0.1964</td></tr><tr><td>(200; 100)</td><td>0.1655</td><td>0.3553</td><td>0.2615</td><td>0.3515</td></tr><tr><td>(50; 50)</td><td>0.2034</td><td>0.5684</td><td>0.3698</td><td>0.4494</td></tr><tr><td>(50; 100)</td><td>0.5289</td><td>0.7130</td><td>0.5473</td><td>0.8042</td></tr></tbody></table>

<table><caption>Table 5.2: Coverage Probability of 95% Asymptotic Confidence Interval for Alpha-Pinene ODE (5.1)</caption><thead><tr><th rowspan="2">(n; σ)</th><th colspan="3">θ₁</th><th colspan="3">θ₂</th><th colspan="3">θ₃</th></tr><tr><th>OC</th><th>OCopt</th><th>NLS</th><th>OC</th><th>OCopt</th><th>NLS</th><th>OC</th><th>OCopt</th><th>NLS</th></tr></thead><tbody><tr><td>(400; 50)</td><td>0.9600</td><td>0.9700</td><td>0.9588</td><td>0.9900</td><td>1.0000</td><td>0.9278</td><td>0.9600</td><td>0.9600</td><td>0.9691</td></tr><tr><td>(400; 100)</td><td>0.9700</td><td>0.9796</td><td>0.9072</td><td>0.9800</td><td>0.9898</td><td>0.8866</td><td>0.9700</td><td>0.9694</td><td>0.9381</td></tr><tr><td>(200; 50)</td><td>0.9800</td><td>0.9800</td><td>0.8866</td><td>1.0000</td><td>1.0000</td><td>0.9175</td><td>0.9600</td><td>0.9600</td><td>0.9278</td></tr><tr><td>(200; 100)</td><td>0.9381</td><td>0.9579</td><td>0.9333</td><td>0.9897</td><td>0.9789</td><td>0.8889</td><td>0.9691</td><td>0.9684</td><td>0.9444</td></tr><tr><td>(50; 50)</td><td>0.9474</td><td>0.9579</td><td>0.8795</td><td>0.9684</td><td>0.9789</td><td>0.8434</td><td>0.9579</td><td>0.9579</td><td>0.9398</td></tr><tr><td>(50; 100)</td><td>0.8681</td><td>0.9176</td><td>0.7571</td><td>0.9341</td><td>0.9647</td><td>0.8286</td><td>0.8462</td><td>0.8824</td><td>0.8143</td></tr></tbody></table>

but we consider only the coverage probabilities of the confidence intervals given in tables (5.2, 5.3).

We see that the asymptotic approximation obtained for the parameters of is really good for $\theta_1, \theta_2, \theta_3$,

and are degraded for $\theta_4, \theta_5, \theta_6$. Finally, we see that the behavior of OCopt is similar to NLS, and better

than simple OC.

5

6

<table><caption>Table 5.3: Coverage Probability of 95% Asymptotic Confidence Interval for Alpha-Pinene ODE (5.1)</caption><thead><tr><th rowspan="2">(n; σ)</th><th colspan="3">θ₄</th><th colspan="3">θ₅</th><th colspan="3">θ₆</th></tr><tr><th>OC</th><th>OCopt</th><th>NLS</th><th>OC</th><th>OCopt</th><th>NLS</th><th>OC</th><th>OCopt</th><th>NLS</th></tr></thead><tbody><tr><td>(400; 50)</td><td>0.7800</td><td>0.8400</td><td>0.9278</td><td>0.7700</td><td>0.8400</td><td>0.9072</td><td>0.8200</td><td>0.8500</td><td>0.9175</td></tr><tr><td>(400; 100)</td><td>0.6300</td><td>0.6939</td><td>0.8763</td><td>0.6400</td><td>0.7245</td><td>0.8866</td><td>0.6600</td><td>0.7449</td><td>0.8763</td></tr><tr><td>(200; 50)</td><td>0.6800</td><td>0.7500</td><td>0.9278</td><td>0.6700</td><td>0.7500</td><td>0.9278</td><td>0.7000</td><td>0.7800</td><td>0.9278</td></tr><tr><td>(200; 100)</td><td>0.5464</td><td>0.6526</td><td>0.8889</td><td>0.5773</td><td>0.6842</td><td>0.8889</td><td>0.7113</td><td>0.8000</td><td>0.8778</td></tr><tr><td>(50; 50)</td><td>0.5474</td><td>0.6737</td><td>0.8795</td><td>0.7263</td><td>0.8632</td><td>0.8795</td><td>0.7158</td><td>0.7895</td><td>0.9036</td></tr><tr><td>(50; 100)</td><td>0.6154</td><td>0.7294</td><td>0.7714</td><td>0.8352</td><td>0.9765</td><td>0.8143</td><td>0.6813</td><td>0.7412</td><td>0.7857</td></tr></tbody></table>
## 4.3 Nonlinear differential equations

We give here general results for the asymptotics of $e_\ell(\hat{\phi}, \theta)$ when the functional is linear or not in $\hat{\phi}$.

In [24], the root-$n$ consistency and asymptotic normality is obtained if the functional $g \mapsto e_\ell(g, \theta)$ has

a continuous Fréchet derivative $De_\ell(g, \theta)$ with respect to the norm $\|\cdot\|_\infty$. If $x \mapsto f(t, x, \theta)$ is twice

continuously differentiable for $t \in [0, 1]$ a.e. in $x$ and $\theta$ in $\Theta$, then we can compute easily its Fréchet

derivative for $g \in H^1$ in the uniform ball $\|g - \phi^*\|_\infty \le D$. For all $h \in H^1$ such $\|g + h - \phi^*\|_\infty \le D$, we

have

$$e_\ell(g + h, \theta) - e_\ell(g, \theta) = \langle f_x(\cdot, g, \theta) h, \varphi_\ell \rangle + \langle h, \dot{\varphi} \rangle + \langle h^\top f_{xx}(\cdot, \tilde{g}, \theta) h, \varphi_\ell \rangle \quad (4.8)$$

by a Taylor expansion around $g$. As in the linear case, we introduce the tangent linear operator

$A_g(\theta) : u \mapsto \dot{u} - a_g(t, \theta)u$ with $a_g(t, \theta) = f_{x_1}(t, g(t), \theta)$ and the function $b_g(t, \theta) = f_{x_2}(t, g(t), \theta)$.

For all $\theta$, the Fréchet derivative of $e_\ell(g, \theta)$ (w.r.t to the uniform norm) is the linear operator $h =$

$(h_1, h_2) \mapsto De_\ell(g, \theta).h = \langle h_1, \dot{\varphi}_\ell + a_g(t, \theta)\varphi_\ell \rangle + \langle h_2, b_g(\cdot, \theta)\varphi_\ell \rangle$ and satisfies for all $\theta \in \Theta$

$$|e_\ell(g + h, \theta) - e_\ell(g, \theta) - De_\ell(g, \theta).h| \le C \|h\|_\infty^2$$

because $f_{xx}$ is uniformly dominated on $D \times \Theta$. Moreover, for all $\epsilon$ (with $0 < \epsilon < D$), for all $g, g'$ such

that $\|g - \phi^*\|_\infty, \|g' - \phi^*\|_\infty \le \epsilon$, we have

$$\begin{aligned} |De_\ell(g, \theta).h - De_\ell(g', \theta).h| &\le \int_0^1 h(t)^\top f_{xx}(t, \tilde{g}(t), \theta) (g(t) - g'(t)) \varphi_\ell(t) dt \\ &\le C \|h\|_\infty \|g - g'\|_\infty \end{aligned}$$

with $C$, a constant independent of $\theta, \epsilon$ and $g, g'$ (because $f_{xx}$ is uniformly dominated).

As in the linear case, we need to evaluate $De_\ell(g, \theta)$ on the basis $p^K$. We denote $A(g, \theta)$ and $B(g, \theta)$

the matrices in $\mathbb{R}^{L \times K}$ with entries $\int_0^1 a_g(t, \theta)\varphi_\ell(t)p_{kK}dt$ and $\int_0^1 b_g(t, \theta)\varphi_\ell(t)p_{kK}dt$ (respectively) and we

have the approximation

$$\mathbf{e}_L(\hat{\phi}, \theta) = \mathbf{e}_L(\phi^*, \theta) + \mathbf{A}(\phi^*, \theta)\hat{\mathbf{c}}_1 + \mathbf{B}(\phi^*, \theta)\hat{\mathbf{c}}_2 + O(\|h\|_\infty^2). \quad (4.9)$$
## 1 5.2 Nonlinear ODEs

### 2 5.2.1 Ricatti equation

3 Ricatti ODE is a family of quadratic ODE that arises commonly in Control Theory. In this example, 4 we consider the simple constant coefficient Ricatti equation

$$ \dot{x} = ax^2 + b\sqrt{t} \qquad (5.2) $$

5 with $a = 0.12$, $b = 0.2$ and $x_0 = -1$, for $t \in [0, 5]$. Because of the squared term $x^2$, the solution can 6 explode on $[0, 5]$, typically if $a$ is too big. This is a problem for NLS, because the parameter optimization 7 can give rise to irrelevant parameter candidate with diverging solution before $T = 5$. This problem does 8 not appear for OC estimator, because we use only the weak form: $\forall \varphi \in C_C^1([0, 1])$,

$$ \langle x, \dot{\varphi} \rangle + a \langle x^2, \varphi \rangle + b \langle \sqrt{t}, \varphi \rangle = 0. \qquad (5.3) $$

9 As the Ricatti equation is linear in the parameter, the OC estimator is computed by solving a similar 10 linear system as in the previous linear example. The essential difference is that the matrix $X^\varphi$ involves 11 nonlinear functionals of $\hat{\phi}$. The Mean Square Errors show that NLS and TS provide better estimation 12 than OC and OC opt. For all the estimators, the bias is similar (and $b$ is notably harder to estimate 13 than $a$), but the estimated asymptotic variance is smaller for NLS, and is smaller than for OC and OC 14 opt. As a consequence, the coverage probability are better for OC and OCopt than NLS (in particular 15 for $b$, see tables 5.6, 5.7). This is classical short-come of NLS to give overoptimistic variance estimators.

<table><caption>Table 5.4: Mean Square Error estimated by Monte Carlo for $a^* = 0.12$ for Ricatti ODE (5.2)</caption><thead><tr><th>(n; σ)</th><th>â<sub>TS</sub></th><th>â<sub>OC</sub></th><th>â<sub>OCopt</sub></th><th>â<sub>NLS</sub></th></tr></thead><tbody><tr><td>(400; 0.2)</td><td>0.0007</td><td>0.0061</td><td>0.0060</td><td>0.0007</td></tr><tr><td>(400; 0.4)</td><td>0.0028</td><td>0.0119</td><td>0.0115</td><td>0.0029</td></tr><tr><td>(200; 0.2)</td><td>0.0015</td><td>0.0105</td><td>0.0106</td><td>0.0016</td></tr><tr><td>(200; 0.4)</td><td>0.0047</td><td>0.1215</td><td>0.0646</td><td>0.0055</td></tr><tr><td>(50; 0.2)</td><td>0.0062</td><td>0.0350</td><td>0.0349</td><td>0.0077</td></tr><tr><td>(50; 0.4)</td><td>0.0213</td><td>0.1755</td><td>0.0831</td><td>0.0310</td></tr></tbody></table>
1 We can derive the asymptotic variance of $\mathbf{e}_L(\hat{\phi}, \theta)$ from (4.9)

$$V_{e,L}(\theta) = \mathbf{A}(\phi^*, \theta)V(\hat{\mathbf{c}}_1)\mathbf{A}(\phi^*, \theta)^{\top} + \mathbf{B}(\phi^*, \theta)V(\hat{\mathbf{c}}_2)\mathbf{B}(\phi^*, \theta)^{\top} \quad (4.10)$$

2 and we can get an estimate $\widehat{V}_{e,L}(\theta)$ from the data as in the linear case.

3

4 In order to assess the previous discussion and for deriving the root-$n$ rate of our estimator, we introduce the following two conditions:

6 **Condition C7: (a)** The times $T_1, \dots, T_n$ have a density $\pi$ w.r.t. Lebesgue measure such $0 < c < \pi < C < \infty$; (b) $E[\epsilon^4] < \infty$.

7

8 **Condition C8:** For $\ell = 1, \dots, L$, $\theta \in \Theta$, there exists $\tilde{\beta}_{K_\ell}$ in $\mathbb{R}^{K_\ell}$ with $\| \frac{f_x(\cdot, \phi^*, \theta)\varphi_\ell + \dot{\varphi}_\ell}{\pi} - \tilde{\beta}_{K_\ell}^\top p^{K_\ell} \|_{L^2} \to 0$.

9 Conditions C7 and C8 are similar to the assumptions given in [24]. Condition C8 is here to ensure that the Fréchet derivative $De_\ell(\phi^*, \theta)$ that drives the asymptotic rate of $e_\ell(g, \theta)$ (see equation 4.8) can be well approximated in the basis $\mathbf{p}^K$ as the nonparametric proxy. Then the linearized nonlinear functional of the nonparametric estimator is well approximated by a linear combination of the regression coefficients. When we use B-splines with uniform knot sequence, condition C8 can be replaced by the simpler condition C9:

10

11

12

13

14

15 **Condition C9: (a)** The series estimator is a regression spline with a uniform knot sequence $(\tau_{1,K}, \dots, \tau_{N_K,K})$ defining the spline basis $\mathbf{p}^K$ satisfies $\max_i |\tau_{i+1,K} - \tau_{i,K}| \to 0$ as $K \to \infty$; (b) For all $\theta \in \Theta$, for $\ell = 1 \dots L$, $v_\ell: t \mapsto \frac{f_x(t, \phi^*(t), \theta)\varphi_\ell(t) + \dot{\varphi}_\ell(t)}{\pi(t)}$ is $C^1$.

16

17

18 **Theorem 4.1.** If either the following conditions are satisfied:

19 $\mathbf{p}^K$ is a general series estimators Under conditions C1-C8 and if $f$ is a linear vector field or, $f$ is a nonlinear vector field and $K$ is chosen such that $\frac{\zeta_0(K)^4 K^2}{n} \to 0$

20

21 $\mathbf{p}^K$ is a uniform knot splines Under conditions C1-C2(a), C3-C7, C9 and if $f$ is a linear vector field and $\frac{K^2}{n} \to 0$, or $f$ is a nonlinear vector field and $\frac{K^4}{n} \to 0$

22
<table><caption>Table 5.5: Mean Square Error estimated by Monte Carlo for $b^* = 0.2$ for Ricatti ODE (5.2)</caption><thead><tr><th>(n; σ)</th><th>b̂<sub>TS</sub></th><th>b̂<sub>OC</sub></th><th>b̂<sub>OCopt</sub></th><th>b̂<sub>NLS</sub></th></tr></thead><tbody><tr><td>(400; 0.2)</td><td>0.0000</td><td>0.0001</td><td>0.0001</td><td>0.0000</td></tr><tr><td>(400; 0.4)</td><td>0.0000</td><td>0.0001</td><td>0.0001</td><td>0.0000</td></tr><tr><td>(200; 0.2)</td><td>0.0000</td><td>0.0001</td><td>0.0001</td><td>0.0000</td></tr><tr><td>(200; 0.4)</td><td>0.0001</td><td>0.0011</td><td>0.0007</td><td>0.0001</td></tr><tr><td>(50; 0.2)</td><td>0.0001</td><td>0.0003</td><td>0.0003</td><td>0.0001</td></tr><tr><td>(50; 0.4)</td><td>0.0003</td><td>0.0025</td><td>0.0011</td><td>0.0004</td></tr></tbody></table>

<table><caption>Table 5.6: Estimated Variance and Coverage Probabilities for $a$, for Ricatti ODE (5.2) with $a^* = 0.12$.</caption><thead><tr><th rowspan="2">(n; σ)</th><th colspan="3">Estimated Asymptotic Variance</th><th colspan="3">Coverage Probability (95%)</th></tr><tr><th>V(̂a<sub>OC</sub>)</th><th>V(̂a<sub>OCopt</sub>)</th><th>V(̂a<sub>NLS</sub>)</th><th>OC</th><th>OC opt</th><th>NLS</th></tr></thead><tbody><tr><td>(400; 0.2)</td><td>0.0047</td><td>0.0047</td><td>0.0007</td><td>0.95</td><td>0.93</td><td>0.89</td></tr><tr><td>(400; 0.4)</td><td>0.0147</td><td>0.0145</td><td>0.0028</td><td>0.95</td><td>0.97</td><td>0.90</td></tr><tr><td>(200; 0.2)</td><td>0.0120</td><td>0.0120</td><td>0.0014</td><td>0.94</td><td>0.93</td><td>0.80</td></tr><tr><td>(200; 0.4)</td><td>0.6987</td><td>0.0797</td><td>0.0054</td><td>0.93</td><td>0.93</td><td>0.90</td></tr><tr><td>(50; 0.2)</td><td>0.0658</td><td>0.2945</td><td>0.0150</td><td>0.89</td><td>0.87</td><td>0.80</td></tr><tr><td>(50; 0.4)</td><td>2.9971</td><td>0.1584</td><td>0.0231</td><td>0.89</td><td>0.84</td><td>0.85</td></tr></tbody></table>

<table><caption>Table 5.7: Estimated Variance and Coverage Probabilities for $a$, for ODE (5.2) with $b^* = 0.2$.</caption><thead><tr><th rowspan="2">(n; σ)</th><th colspan="3">Estimated Asymptotic Variance</th><th colspan="3">Coverage Probability (95%)</th></tr><tr><th>V̂(b̂<sub>OC</sub>)</th><th>V̂(b̂<sub>OCopt</sub>)</th><th>V̂(b̂<sub>NLS</sub>)</th><th>OC</th><th>OC opt</th><th>NLS</th></tr></thead><tbody><tr><td>(400; 0.2)</td><td>0.0001</td><td>0.0001</td><td>0.0000</td><td>0.96</td><td>0.97</td><td>0.82</td></tr><tr><td>(400; 0.4)</td><td>0.0002</td><td>0.0002</td><td>0.0000</td><td>0.97</td><td>0.98</td><td>0.89</td></tr><tr><td>(200; 0.2)</td><td>0.0001</td><td>0.0001</td><td>0.0000</td><td>0.97</td><td>0.97</td><td>0.84</td></tr><tr><td>(200; 0.4)</td><td>0.0058</td><td>0.0009</td><td>0.0001</td><td>0.97</td><td>0.97</td><td>0.89</td></tr><tr><td>(50; 0.2)</td><td>0.0006</td><td>0.0032</td><td>0.0106</td><td>0.96</td><td>0.96</td><td>0.82</td></tr><tr><td>(50; 0.4)</td><td>0.0363</td><td>0.0020</td><td>0.0004</td><td>0.95</td><td>0.95</td><td>0.88</td></tr></tbody></table>

### 5.2.2 FitzHugh-Nagumo

The FitzHugh-Nagumo is a nonlinear two-dimensional ODE introduced for modeling neurons. For well-chosen sets of parameters and initial conditions, it exhibits a periodic behavior, with typical oscillations
1 Then $\hat{\theta}_{n,L}$ is such that

$$ \sqrt{n} (\hat{\theta}_{n,L} - \theta_{L}^{*}) \rightsquigarrow N(0, \mathbf{V}_{L}^{*}) \qquad (4.11) $$

2 with

$$ \mathbf{V}_{L}^{*} = \mathbf{M}_{L}^{*} \mathbf{V}_{e,L}^{*} \mathbf{M}_{L}^{*\top}. \qquad (4.12) $$

3 where $\mathbf{V}_{e,L}^{*} = V_{e,L}(\theta_{L}^{*})$. The asymptotic variance can be estimated by $\hat{\mathbf{M}}_{L} \mathbf{V}_{e,L}(\hat{\theta}_{n,L}) \hat{\mathbf{M}}_{L}^{\top} \xrightarrow{P} \mathbf{V}_{L}^{*}$. In 4 particular, if we use regression splines and $t \mapsto f(t, \phi^{*}(t), \theta)$ is $C^s$ on $[0, 1]$ with $s \ge 3$, then (4.11) holds 5 with $K$ such that $\sqrt{n}K^{-s} \to 0$ and $n^{-1}K^4 \to 0$.

6 Moreover, if $L = L(n) \to \infty, n \to \infty$ is chosen such that the bias $\mathbf{B}_{L(n)} = O(n^{-1/2})$, then we have

$$ \hat{\theta}_{n,L(n)} - \theta^{*} = O_{P}(n^{-1/2}). \qquad (4.13) $$

7 In particular, this is the case when the test functions $\phi_{\ell}$ are the sine basis, and $L(n) = O(n^{\alpha})$ with $\alpha > 1/2$.

8 This theorem is a direct application of theorem 3 in [24] that claims the root-$n$ consistency and 9 asymptotic normality of general (nonlinear) plug-in estimators. The main steps of the proof are given 10 in Supplementary Material I.

## 5 Experiments

### 5.1 Description of the setting

We compare the NLS estimator $\hat{\theta}^{NLS}$, the Two-Step Estimator (TS) $\hat{\theta}^{TS}$ and the OC estimator $\hat{\theta}^{OC}$ for 14 varying sample sizes ($n = 400, 200, 50$) and varying noise levels (high and small). We consider 3 different 15 ODEs with different mathematical structure: the $\alpha$-pinene ODE (linear in state and in parameter), the 16 Ricatti ODE (nonlinear in state, linear in parameter) and the FitzHugh-Nagumo ODE (nonlinear in 17 state and in parameter). Experiments on these three different models provide a good idea of the behavior 18 of the different estimators in terms of the robustness, consistency and efficiency. It helps also in assessing 19 the quality of the linear approximation for the asymptotics (in particular for the computation of the 20
1 corresponding to a limit cycle

$$ egin{cases} \dot{V} = c \left(V - \frac{V^3}{3} + R\right) \\ \dot{R} = -\frac{1}{c} (V - a + bR) \end{cases} \qquad (5.4) $$

2 The true parameters are $a^* = b^* = 0.2$ and $c^* = 3$ and $x_0 = (V_0, R_0) = (-1, 1)$, and are taken from [21] where it was introduced as a benchmark for parameter estimation in ODEs. The noise levels we consider are $\sigma = 0.15, 0.3$. In all the settings we consider (large or small sample size), the NLS provides the smallest MSE, and it is clear that OC improves on TS, and that the optimal weight provides an even better estimator OCopt, see table 5.8.

7 The detailed analysis of the bias show that the NLS have the smallest bias, and that the parameter $c$ is the hardest to estimate. The most essential feature is that in most of the case, the asymptotic variance computed by the NLS is smaller than for OC, which gives significant differences between the coverage probabilities for $a, b, c$ in particular in the small sample case. The confidence sets obtained by OC are rather credible, with coverage probability around 90% in almost any case (except in table 5.11). As the choice of the optimal weighted estimator tends to reduce the variance, we get narrower confidence sets with lower coverage probability.

<table><caption>Table 5.8: Mean Square Error computed by Monte Carlo for FitzHugh-Nagumo ODE (5.4).</caption><thead><tr><th>×10<sup>-3</sup></th><th colspan="4">a</th><th colspan="4">b</th><th colspan="4">c</th></tr><tr><th>(n; σ)</th><th>TS</th><th>OC</th><th>OCopt</th><th>NLS</th><th>TS</th><th>OC</th><th>OCopt</th><th>NLS</th><th>TS</th><th>OC</th><th>OCopt</th><th>NLS</th></tr></thead><tbody><tr><td>(400; 0.15)</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.0</td><td>0.8</td><td>3.2</td><td>1.3</td><td>0.5</td><td>32.4</td><td>13.2</td><td>4.7</td><td>0.2</td></tr><tr><td>(400; 0.3)</td><td>0.3</td><td>0.6</td><td>0.4</td><td>0.1</td><td>4.1</td><td>23.5</td><td>14.2</td><td>1.6</td><td>236.1</td><td>41.6</td><td>18.2</td><td>0.3</td></tr><tr><td>(200; 0.15)</td><td>0.2</td><td>0.4</td><td>0.2</td><td>0.0</td><td>2.0</td><td>23.7</td><td>2.4</td><td>0.9</td><td>110.4</td><td>29.2</td><td>9.1</td><td>0.5</td></tr><tr><td>(200; 0.3)</td><td>0.6</td><td>1.0</td><td>0.6</td><td>0.2</td><td>5.8</td><td>26.4</td><td>9.1</td><td>3.7</td><td>519.1</td><td>82.0</td><td>34.0</td><td>1.8</td></tr><tr><td>(50; 0.15)</td><td>0.3</td><td>0.6</td><td>0.3</td><td>0.1</td><td>3.1</td><td>36.7</td><td>9.0</td><td>1.6</td><td>228.4</td><td>97.1</td><td>35.1</td><td>1.3</td></tr><tr><td>(50; 0.3)</td><td>1.3</td><td>1.6</td><td>1.0</td><td>0.3</td><td>12.4</td><td>51.1</td><td>44.1</td><td>7.4</td><td>893.0</td><td>796.9</td><td>404.8</td><td>3.8</td></tr></tbody></table>

14

15

16

17
covariance matrix).

In the simulations, the noise is homoscedastic and Gaussian, so that the NLS are asymptotically efficient. Hence, the settings $n = 200$ or $n = 400$ indicates the efficiency loss of the Gradient Matching estimators whereas the small size setting ($n = 50$) gives some information on the small sample case, where the asymptotic approximations cannot be assessed.

As the standard reference method, the Sum of Squared Errors (SSE) is minimized by a Levenberg-Marquardt algorithm using 20 starting points centered around the true parameter value $\theta^*$, and we retain the best minimum. The solution of the ODE is computed by a Runge-Kutta algorithm of order 4, implemented in the Matlab function ode45. Hence, we expect that we obtain the true NLS estimator, and that the estimated variance is the true best one.

The Gradient Matching estimators (TS and OC) use the same regression spline, decomposed on a B-spline basis with a uniform knots sequence $\xi_k, k = 1, \dots, K$. For each dataset (and each dimension), the number of knots is selected by minimizing the GCV criterion, [32]. For the plain TS estimator, we use a piecewise affine weight function with $w(0) = w(1) = 0$, as in [8].

The Orthogonal Conditions are defined with the sine basis or B-Splines basis. We have to face with the practical problem of finding the best number of conditions $L$, that depends on the model and on $\hat{\phi}$. In each setting, we have fixed a minimum and a maximum number of conditions $L_{min}$ and $L_{max} \le 2 \times d \times p$ and we select the OC estimator $\hat{\theta}_{n,L}$ that gives the smallest prediction error (i.e that minimizes the SSE):

$$ \hat{\theta}^{OC} = \arg \min_{L_{min} \le L \le L_{max}} \sum_{i=1}^{n} \| y_i - \phi(t_i, \hat{\phi}_0, \hat{\theta}_{n,L}) \|^2 $$

where $\hat{\phi}_0 = \hat{\phi}(0)$ is the nonparametric estimate of the initial condition.

We use Monte Carlo simulations, based on $N_{MC} = 500$ independent draws for comparing the estimators. We compute their Mean Squared Errors $\|\hat{\theta} - \theta^*\|^2$. The accuracy of the estimator is roughly estimated by the trace of the covariance matrices of the estimators, denoted $\mathrm{Tr}(V(\hat{\theta}))$. Moreover, the reliability of the estimates (and asymptotic approximation) is evaluated with the coverage probabilities of the 95% confidence ellipse (except in the case of TS because there is no closed-form for asymptotic variance). For the NLS, the asymptotic variance is computed via the Matlab function nlinfit. A more detailed analysis of the experiments (including coverage probabilities of confidence sets) are given in
<table><caption>Table 5.9: Estimated Variance and Coverage Probabilities for a, for FitzHugh-Nagumo ODE (5.4)</caption><thead><tr><th rowspan="2">(n; σ)</th><th colspan="3">Estimated Asymptotic Variance</th><th colspan="3">Coverage Probability</th></tr><tr><th>V̂ (â<sub>OC</sub>)</th><th>V̂ (â<sub>OCopt</sub>)</th><th>V̂ (â<sub>NLS</sub>)</th><th>OC</th><th>OC opt</th><th>NLS</th></tr></thead><tbody><tr><td>(400; 0.15)</td><td>0.0001</td><td>0.0001</td><td>0.0000</td><td>0.9600</td><td>0.9200</td><td>0.7600</td></tr><tr><td>(400; 0.3)</td><td>0.0011</td><td>0.0002</td><td>0.0000</td><td>0.9394</td><td>0.8750</td><td>0.7292</td></tr><tr><td>(200; 0.15)</td><td>0.0004</td><td>0.0001</td><td>0.0000</td><td>0.8800</td><td>0.8163</td><td>0.7551</td></tr><tr><td>(200; 0.3)</td><td>0.0063</td><td>0.0004</td><td>0.0001</td><td>0.9375</td><td>0.8780</td><td>0.6098</td></tr><tr><td>(50; 0.15)</td><td>0.0577</td><td>0.0002</td><td>0.0000</td><td>0.9400</td><td>0.8876</td><td>0.6180</td></tr><tr><td>(50; 0.3)</td><td>0.7871</td><td>0.0008</td><td>0.0002</td><td>0.9043</td><td>0.9143</td><td>0.6377</td></tr></tbody></table>

<table><caption>Table 5.10: Estimated Variance and Coverage Probabilities for b, for FitzHugh-Nagumo ODE (5.4)</caption><thead><tr><th rowspan="2">(n; σ)</th><th colspan="3">Estimated Asymptotic Variance</th><th colspan="3">Coverage Probability</th></tr><tr><th>V̂ (b̂<sub>OC</sub>)</th><th>V̂ (b̂<sub>OCopt</sub>)</th><th>V̂ (b̂<sub>NLS</sub>)</th><th>OC</th><th>OC opt</th><th>NLS</th></tr></thead><tbody><tr><td>(400; 0.15)</td><td>0.0039</td><td>0.0007</td><td>0.0002</td><td>0.9500</td><td>0.9000</td><td>0.6400</td></tr><tr><td>(400; 0.3)</td><td>0.1516</td><td>0.0058</td><td>0.0008</td><td>0.8889</td><td>0.8333</td><td>0.5521</td></tr><tr><td>(200; 0.15)</td><td>0.0568</td><td>0.0018</td><td>0.0003</td><td>0.9400</td><td>0.8980</td><td>0.5408</td></tr><tr><td>(200; 0.3)</td><td>2.6567</td><td>0.0088</td><td>0.0012</td><td>0.8958</td><td>0.9146</td><td>0.5000</td></tr><tr><td>(50; 0.15)</td><td>39.2617</td><td>0.0048</td><td>0.0006</td><td>0.9000</td><td>0.8315</td><td>0.6180</td></tr><tr><td>(50; 0.3)</td><td>116.2820</td><td>0.0299</td><td>0.0024</td><td>0.8936</td><td>0.8857</td><td>0.4928</td></tr></tbody></table>

# 1 6 Discussion

The experiments done in the previous section show that the Orthogonal Conditions estimator is a reliable estimator that can compete with nonlinear least squares. Indeed, the OC estimators are used as starting point for NLS, so that the NLS improved estimator must be better than OC. When we consider the MSE, it is often the case (except for alpha-pinene), which justifies the common practice of using Gradient Matching for determining interesting initial values for NLS. Nevertheless, one can note that this improvement is rather limited with respect to OC (or OCopt). In particular, this shows that the projection and the use of limited number L of test functions $\phi_{\ell}$ is sufficient for estimating parameters, and that the bias $r_L$ introduced in theorem 1 is in fact quite small. When performing simulations, it is not easy to identify the influence of the quality of the nonparametric proxy in the global quality of a Gradient Matching Estimator (typically OC), in particular when the sample size is small. Hence, some additional work should be done for the construction of good (or the best) nonparametric proxy to use in a Gradient Matching approach in order to control the influence of this step in estimation.

Moreover, the asymptotic variance and the Gaussian approximation have proved to be effective
## 2 5.2 α-pinene

A linear ODE with constant coefficients is written $\dot{x} = Ax$, where $A^T = (A_1|\dots|A_d)$. For $i = 1, \dots, d$, the weak formulation gives the identity $Y_i^\varphi = X^\varphi A_i$ to be satisfied, where $X^\varphi$ is a $d \times L$ matrix with entries $\langle x_k, \varphi_\ell \rangle$ and $Y_i^\varphi$ is a vector in $\mathbb{R}^L$ with entries equal to $-\langle x_i, \dot{\varphi}_\ell \rangle$. For illustration, we consider the α-pinene ODE used in [22] for the comparison of several global optimization algorithms:

$$ \left\{ \begin{aligned} \dot{x}_1 &= -(\theta_1 + \theta_2)x_1 \\ \dot{x}_2 &= \theta_1 x_1 \\ \dot{x}_3 &= \theta_2 x_1 - (\theta_3 + \theta_4)x_3 + \theta_5 x_5 \\ \dot{x}_4 &= \theta_3 x_3 \\ \dot{x}_5 &= \theta_4 x_3 - \theta_5 x_5 \end{aligned} \right. \qquad (5.1) $$

The true parameter to be estimated from a completely observed trajectory on $[0, 100]$ is $\theta^* = (\theta_1, \theta_2, \theta_3, \theta_4, \theta_5)^T$.

As this ODE is linear and time-invariant, we have a closed-form for the solution $\phi^*(t, \theta, \phi_0) = e^{tA}\phi_0$ that can be directly used for the computation of the NLS estimator.

The test functions used for the OC estimators are B-Splines (with uniform knots sequence) $\varphi_\ell$, $\ell = 1, \dots, L$ with compact support included in $]0, 20[$. We consider a varying number of conditions $L$, i.e $2 \le L \le 15$. Finally, we have two settings for the estimation of $\theta$: when the initial condition $\phi_0$ is known (and equal to $(100, 0, 0, 0, 0)^T$ as in [31]), and when $\phi_0$ is unknown and needs to be estimated (for NLS).

### 5.2.1 Known initial condition

For the OC and TS estimator, we constrain the spline estimator $\hat{\phi}$ to satisfy the condition $\hat{\phi}(0) = \phi_0$ (by adding a linear constraint to the classical least-squares minimization). Moreover, following section 2.3, we integrate the knowledge of the initial condition by adding a test function $\varphi_0$ which is a B-spline with $\varphi_0(0) \ne 0$. Hence, we define 2 different OC estimators, respectively, $\hat{\theta}^{OC,0}$ and $\hat{\theta}^{OC}$ that uses or not (resp.) the knowledge of the initial condition.
<table><caption>Table 5.11: Estimated Variance and Coverage Probabilities for c, for FitzHugh-Nagumo (5.4)</caption><thead><tr><th rowspan="2">(n; σ)</th><th colspan="3">Estimated Asymptotic Variance</th><th colspan="3">Coverage Probability</th></tr><tr><th>V̂ (ĉ<sub>OC</sub>)</th><th>V̂ (ĉ<sub>OCopt</sub>)</th><th>V̂ (ĉ<sub>NLS</sub>)</th><th>OC</th><th>OC opt</th><th>NLS</th></tr></thead><tbody><tr><td>(400; 0.15)</td><td>0.0122</td><td>0.0012</td><td>0.0000</td><td>0.9300</td><td>0.7300</td><td>0.6100</td></tr><tr><td>(400; 0.3)</td><td>0.0827</td><td>0.0058</td><td>0.0002</td><td>0.9495</td><td>0.7500</td><td>0.7083</td></tr><tr><td>(200; 0.15)</td><td>0.0401</td><td>0.0024</td><td>0.0001</td><td>0.9500</td><td>0.6531</td><td>0.6735</td></tr><tr><td>(200; 0.3)</td><td>0.5445</td><td>0.0122</td><td>0.0003</td><td>0.7812</td><td>0.7195</td><td>0.5488</td></tr><tr><td>(50; 0.15)</td><td>15.8363</td><td>0.0060</td><td>0.0002</td><td>0.8300</td><td>0.5730</td><td>0.6067</td></tr><tr><td>(50; 0.3)</td><td>3.1898</td><td>0.1295</td><td>0.0008</td><td>0.6809</td><td>0.8000</td><td>0.6667</td></tr></tbody></table>

1 for quantifying the uncertainty about the OC parameters, through classical confidence which is often
delicate for Gradient Matching estimators. From the simulation, it is clear that these estimation are over-
optimistic, which gives too narrow confidence sets: this is a classical problem for nonlinear regression,
where NLS tends also to give unreliable small confidence sets. As the variance $V_L^*$ is a bit bigger,
we can have wider and possibly more relevant confidence sets. Conversely, the OC estimator does not
use errors-in-variables technics as in [12], that can provide huge variances and confidence sets. Hence,
the question of the “best” set estimation in ODEs is still hard to close, as the classical linearization
approaches used for asymptotic approximation can be limited in Differential Equations.

9 It should be noticed that the use of an optimal weight in OC can significantly improve the estimator
10 (by reducing the bias and the variance); but some care must be taken in the iterations as the estimated
11 variance $V(e_L)$ can be singular. A simple device to avoid numerical divergence is to reduce the number
12 of orthogonal conditions.

13 We have shown that there is a real interest in defining more elaborated criterion in Gradient Matching
14 estimators for statistical tractability. In particular, variational characterization is quite new in statis-
tical estimation, and gives computationally and statistically relevant estimators with a broad field of
15 applicability that still needs to be explored.
16

17 Acknowledgments

18 This work was funded by two ANR projects GD2S (ANR-05-MMSA-0013-01) and ODESSA (ANR-
19 09-SYSC-009-01) and received a partial support from the Analysis of Object Oriented Data program
20 (2010-2011) in Statistical and Applied Mathematical Sciences Institute (SAMSI), USA.
<table><caption>Table 5.1: MSE, Asymptotic Variance for α-pinene model with known Initial Condition</caption><thead><tr><th>×10<sup>-2</sup></th><th colspan="4">MSE</th><th colspan="3">Tr (V(θ̂))</th></tr><tr><th>(n, σ)</th><th>TS</th><th>OC</th><th>OC,0</th><th>NLS</th><th>OC</th><th>OC,0</th><th>NLS</th></tr></thead><tbody><tr><td>(400, 3)</td><td>0.72</td><td>0.05</td><td>0.04</td><td>0.02</td><td>0.04</td><td>0.04</td><td>0.02</td></tr><tr><td>(400, 8)</td><td>2.28</td><td>0.22</td><td>0.25</td><td>0.10</td><td>0.95</td><td>1.20</td><td>0.12</td></tr><tr><td>(200, 3)</td><td>1.19</td><td>0.27</td><td>0.30</td><td>0.03</td><td>0.09</td><td>0.13</td><td>0.03</td></tr><tr><td>(200, 8)</td><td>2.95</td><td>0.44</td><td>0.37</td><td>0.18</td><td>2.66</td><td>2.68</td><td>0.27</td></tr><tr><td>(50, 3)</td><td>2.39</td><td>0.27</td><td>0.26</td><td>0.16</td><td>1.37</td><td>1.58</td><td>0.16</td></tr><tr><td>(50, 8)</td><td>4.54</td><td>1.03</td><td>0.93</td><td>0.68</td><td>7.96</td><td>7.27</td><td>1.68</td></tr></tbody></table>

### 1.5.2.2 Unknown initial condition

In this case, the NLS needs to estimate the initial condition as well, whereas it is not needed for Gradient

Matching estimators and we have the same estimates (for $\hat{\theta}^{TS}$ and $\hat{\theta}^{OC}$) as in the previous section. In

this setting, we consider another OC estimator that uses information about the other boundary $T = 100$.

Indeed, we know that the α-pinene network converges to a stationary point, that is almost reached at

time $T = 100$. Hence the boundary condition $\dot{\phi}^{*}(100) = 0$ can be used for estimation (section 2.3): if $\phi_{1}$

is a test function with $\phi_{1}(100) \neq 0$, we have $A^{2} < \phi^{*}, \phi_{1} > +A < \phi^{*}, \dot{\phi}_{1} > = 0$. This gives an additional

condition to be satisfied for the OC estimator, which is denoted as $\hat{\theta}^{OC,1}$, see section 2.3.

<table><caption>Table 5.2: MSE, Asymptotic Variance for α-pinene model with unknown initial conditions</caption><thead><tr><th>×10<sup>-2</sup></th><th colspan="4">MSE</th><th colspan="3">Tr (V(θ̂))</th></tr><tr><th>(n, σ)</th><th>TS</th><th>OC</th><th>OC,1</th><th>NLS</th><th>OC</th><th>OC,1</th><th>NLS</th></tr></thead><tbody><tr><td>(400, 3)</td><td>0.25</td><td>0.11</td><td>0.11</td><td>0.07</td><td>0.10</td><td>0.10</td><td>0.06</td></tr><tr><td>(400, 8)</td><td>1.07</td><td>0.85</td><td>0.56</td><td>0.50</td><td>1.06</td><td>0.82</td><td>0.61</td></tr><tr><td>(200, 3)</td><td>0.6</td><td>0.37</td><td>0.23</td><td>0.14</td><td>0.25</td><td>0.20</td><td>0.14</td></tr><tr><td>(200, 8)</td><td>1.64</td><td>1.42</td><td>0.83</td><td>1.34</td><td>2.36</td><td>1.64</td><td>1.54</td></tr><tr><td>(50, 3)</td><td>1.33</td><td>1.31</td><td>0.80</td><td>0.69</td><td>1.63</td><td>1.02</td><td>0.76</td></tr><tr><td>(50, 8)</td><td>3.64</td><td>2.11</td><td>1.79</td><td>1.96</td><td>5.34</td><td>2.20</td><td>4.38</td></tr></tbody></table>

9

## 10.5.3 Ricatti Equation

The true ODE is $\dot{\phi} = a\phi^2 + c\sqrt{t} - d'1_{[T_r;14]}$, with $a^* = 0.11$, $c^* = 0.09$, $d^* = 2$ and $\phi_0 = -1$, for

$t \in [0, 14]$. For all $\phi$ in $C^1$ with $\phi(0) = \phi(14) = 0$, we have $\langle \phi, \dot{\phi} \rangle + a \langle \phi^2, \phi \rangle + c \langle \sqrt{t}, \phi \rangle > 0$
1 A Proofs

2 Theorem 1

3 The classical proof for the consistency of an M-estimator $\hat{\theta}_{n,L} = \arg \min_{\theta \in \Theta} Q_{n,L}(\theta)$ such as the orthogonal conditions estimator relies basically on two stages: the uniform convergence of $Q_{n,L}(\theta)$ towards $Q_{L}^{*}(\theta) = \sum_{\ell=1}^{L} |e_{\ell}(\phi^{*}, \theta)|^2$ for $\theta$ in $\Theta$, and the fact that the true parameter $\hat{\theta}_{n,L}$ is a unique isolated global maximum (by theorem 5.7 in [24]). In a first step, we assume that $Q_{L}^{*}(\theta)$ has a unique global minimum denoted $\theta_{L}^{*}$, and we will show that this is indeed the case, and that $\theta_{L}^{*}$ is not very far from the true parameter $\theta^{*}$.

9 We have to show that $\sup_{\theta} |Q_{n,L}(\theta) - Q_{L}^{*}(\theta)| \to 0$ as $n \to \infty$. From the simple additive expression of $Q_{n,L}(\theta)$, we see that it suffices to show the uniform convergence of $|e_{\ell}(\hat{\phi}, \theta)|^2$ to $|e_{\ell}(\phi^{*}, \theta)|^2$. From the inequality $|a^2 - b^2| \le |a - b| (|a| + |b|)$, we have

$$ \left| |e_{\ell}(\hat{\phi}, \theta)|^2 - |e_{\ell}(\phi^{*}, \theta)|^2 \right| \le |e_{\ell}(\hat{\phi}, \theta) - e_{\ell}(\phi^{*}, \theta)| \left( |e_{\ell}(\hat{\phi}, \theta)| + |e_{\ell}(\phi^{*}, \theta)| \right). $$

12 As $\hat{\phi}$ is such that for $j = 1,2$ we have $\|\hat{\phi}_j - \phi_j^*\|_{\infty} = O_P(\zeta_0(K) (\sqrt{K/n} + K^{-\alpha}))$ this means that $\hat{\phi}(t) \in \mathcal{D} = \{x \in \mathbb{R}^2 | \exists t \in [0,1], |x - \phi^*(t)| \le D\}$ with a probability tending to 1. This shows that with a probability tending to 1, we have $|f(t, \hat{\phi}(t), \theta)| \le h(t, \theta) \le H(t)$ (because $f$ is uniformly $L^2$-Caratheodory), moreover $\hat{\phi}$ is also bounded by $M+D > 0$ (because $\phi^*$ is bounded as a continuous function on $[0,1]$) and

$$ \begin{aligned} |e_{\ell}(\hat{\phi}, \theta)| + |e_{\ell}(\phi^{*}, \theta)| &\le |\langle \mathcal{E}(\hat{\phi}, \theta), \varphi_{\ell} \rangle| + |\langle \hat{\phi}, \dot{\varphi}_{\ell} \rangle| + |\langle \mathcal{E}(\phi^{*}, \theta), \varphi_{\ell} \rangle| + |\langle \phi^{*}, \dot{\varphi}_{\ell} \rangle| \\ &\le 2 \|H\|_{L^2} + \|\dot{\varphi}_{\ell}\|_{L^2} (2M + D) \end{aligned} $$

17 Hence $|e_{\ell}(\hat{\phi}, \theta)| + |e_{\ell}(\phi^{*}, \theta)|$ is uniformly bounded in probability, so that the uniform convergence of the criterion boils down to the uniform convergence of $|e_{\ell}(\hat{\phi}, \theta) - e_{\ell}(\phi^{*}, \theta)|$.

19 We can re-write $e_{\ell}(\hat{\phi}, \theta) - e_{\ell}(\phi^{*}, \theta) = \langle \mathcal{E}(\hat{\phi}, \theta) - \mathcal{E}(\phi^{*}, \theta), \varphi_{\ell} \rangle + \langle \hat{\phi} - \phi^{*}, \dot{\varphi}_{\ell} \rangle$ and it is clear that the second right-hand side term converges uniformly to zero in probability for all $\ell \ge 1$. Consequently, we just have to check that $\mathcal{E}(\hat{\phi}, \theta) - \mathcal{E}(\phi^{*}, \theta)$ converges uniformly in $\theta$ to 0 in probability.
1. $-d'(\tilde{\varphi}(14) - \tilde{\varphi}(T_r)) = 0$ where $\tilde{\varphi}$ is the antiderivative of $\varphi$.

2. When $T_r$ is known, we use a cubic B-splines basis with 3 knots at $T_r$, meaning that $\hat{\phi}$ can have a discontinuous derivative at time $T_r$ (hence the curve estimation from noisy data is pretty correct at $T_r$).

3. The curve is mainly flat for $t \in [0, T_r]$ and after $T_r$, one can observe a linear behavior: 3 knots are used to estimate the curve, and their positions are selected manually.

4. When $T_r$ is unknown, it is required to estimate $\theta = (a, c, d', T_r)$. The OC is no more linear in parameters, but $\hat{\theta}^{OC}$ can be computed by solving the general nonlinear program. The Two-Step estimator fails to estimate $T_r$ because the derivative of the solution is badly estimated when $T_r$ is unknown. OC estimators still give reliable estimates as it uses only $\hat{\phi}$ in the criterion. Some care has to be taken for the knots selection because of unknown $T_r$: when $n = 200, 400$ we use a uniform grid of 15 knots on $[0, 14]$. For $n = 50$, we have used 8 knots uniformly located on $[0, 14]$. Nevertheless, the nonparametric estimates are too rough to obtaining any correct estimate $\hat{\theta}^{TS}$.

5. Concerning NLS, we were not able to solve the optimization problem and we cannot give Monte Carlo statistics for the evaluation of NLS. NLS collapses in practice because the optimization problem is hard (severely ill-posed problem). Indeed, the Levenberg-Marquardt algorithm becomes very sensitive to initial conditions and gives different solutions for very close starting values, even in the neighborhood of the true value $\theta^*$. Moreover, we have to face with the problem of explosion of the solutions during the optimization process. In particular, this problem is very delicate because we have to chose $(a, c)$ so that the (potential) explosion of the solution can be balanced by a proper choice of $d'$ and $T_r$. Probably, NLS would benefit from a specific optimization algorithm that could exploit the particular properties of the ODE, but this is out of the scope of the paper.

22

23
1 First of all, we remark that $g \mapsto \langle \mathcal{E}(g, \theta), \varphi_\ell \rangle$ is a continuous function from $B_\infty(\phi^*, D) = \{g \in C([0, 1]) \mid \|g - \phi^*\|_\infty \le D\}$ to $\mathbb{R}$ (w.r.t the sup-norm), because

$$ \begin{aligned} |\langle \mathcal{E}(g, \theta), \varphi_\ell \rangle - \langle \mathcal{E}(g', \theta), \varphi_\ell \rangle| & \le \langle a(\cdot, \theta) |g - g'|, \varphi_\ell \rangle \\ & \le \langle a(\cdot, \theta), \varphi_\ell \rangle \|g - g'\| \end{aligned} $$

3 with $a(\cdot, \theta) \in L^2$ for all $\theta$, as the vector field is $L^2$-Lipschitz. By the continuous mapping theorem, 4 we get the point-wise convergence of $\langle \mathcal{E}(\hat{\phi}, \theta) - \mathcal{E}(\phi^*, \theta), \varphi_\ell \rangle$, and the hard part consist in the uniform 5 convergence. This can be proven by controlling the oscillations in $\theta$ of the process $\mathcal{E} = (\mathcal{E}(\theta))_{\theta \in \Theta} = \langle \langle \mathcal{E}(\hat{\phi}, \theta) - \mathcal{E}(\phi^*, \theta), \varphi_\ell \rangle \rangle_{\theta \in \Theta}$ and by using theorem 18.11 in [24], as convergence in probability towards 6 a constant is equivalent to weak convergence. In order to show that the process $\mathcal{E}$ converges weakly to 7 0 in the space $C(\Theta)$ of continuous functions on $\Theta$ equipped with the supremum norm, we have to check 8 that for all $k$ and all $\theta_1, \dots, \theta_k$ in $\Theta$ 9

$$ (\mathcal{E}(\theta_1), \dots, \mathcal{E}(\theta_k)) \rightsquigarrow (0, \dots, 0) \qquad (A.1) $$

10 and that for all $\epsilon, \alpha > 0$ there exists a partition of $\Theta_1, \dots, \Theta_K$ of $\Theta$ such that

$$ \limsup_{n \to \infty} P \left( \sup_k \sup_{\theta, \theta' \in \Theta_k} |\mathcal{E}(\theta) - \mathcal{E}(\theta')| \ge \alpha \right) \le \epsilon. \qquad (A.2) $$

11 The first condition (A.1) is a direct consequence of the point-wise convergence in probability of $\mathcal{E}(\theta)$ in $\Theta$. Concerning the second condition, from we have for $t$ in $[0, 1]$ a.e. 12

$$ \begin{aligned} \left| f(t, \hat{\phi}(t), \theta) - f(t, \hat{\phi}(t), \theta') \right| & \le \sum_{i=1}^{p} \left| f_{\theta_i}(t, \hat{\phi}(t), \tilde{\theta}_n(t)) \right| \left| \theta_i - \theta_i' \right| \\ & \le \bar{a}'(t) \left| \theta - \theta' \right| \end{aligned} $$

13 with a probability tending to 1 ($\tilde{\theta}_n(t)$ being a parameter between $\theta$ and $\theta'$). As a consequence the
<table><caption>Table 5.3: MSE, Tr(V(θ̂)) for Parameter estimation for Ricatti Equation with known T_r</caption><thead><tr><th>×10<sup>-2</sup></th><th colspan="3">MSE</th><th colspan="2">Tr(V(θ̂))</th></tr><tr><th>(n, σ)</th><th>TS</th><th>OC</th><th>NLS</th><th>OC</th><th>NLS</th></tr></thead><tbody><tr><td>(400, 0.2)</td><td>0.18</td><td>0.27</td><td>0.58</td><td>1.76</td><td>0.10</td></tr><tr><td>(400, 0.4)</td><td>0.78</td><td>1.21</td><td>0.94</td><td>2.56</td><td>0.38</td></tr><tr><td>(200, 0.2)</td><td>0.33</td><td>0.87</td><td>0.57</td><td>2.85</td><td>0.25</td></tr><tr><td>(200, 0.4)</td><td>1.12</td><td>2.69</td><td>1.12</td><td>5.64</td><td>0.98</td></tr><tr><td>(50, 0.2)</td><td>1.03</td><td>1.30</td><td>1.54</td><td>4.70</td><td>1.00</td></tr><tr><td>(50, 0.4)</td><td>3.80</td><td>4.43</td><td>3.94</td><td>8.89</td><td>4.08</td></tr></tbody></table>

<table><caption>Table 5.4: MSE, Sum Empirical Variance for Parameter estimation for Ricatti with unknown T_r</caption><thead><tr><th>×10<sup>-2</sup></th><th>MSE(â)</th><th>MSE(ĉ)</th><th>MSE(d̂)</th><th>MSE(T̂<sub>r</sub>)</th></tr><tr><th>(n, σ)</th><th>OC</th><th>OC</th><th>OC</th><th>OC</th></tr></thead><tbody><tr><td>(400, 0.2)</td><td>0.09</td><td>0.00</td><td>2.54</td><td>1.39</td></tr><tr><td>(400, 0.4)</td><td>0.29</td><td>0.01</td><td>4.27</td><td>3.54</td></tr><tr><td>(200, 0.2)</td><td>0.21</td><td>0.00</td><td>4.08</td><td>3.18</td></tr><tr><td>(200, 0.4)</td><td>0.61</td><td>0.01</td><td>11.96</td><td>6.93</td></tr><tr><td>(50, 0.4)</td><td>0.64</td><td>0.02</td><td>11.20</td><td>14.25</td></tr><tr><td>(50, 0.4)</td><td>0.77</td><td>0.01</td><td>17.18</td><td>19.40</td></tr></tbody></table><table><thead><tr><th>×10<sup>-2</sup></th><th>MSE</th><th>Tr(V(θ̂))</th></tr><tr><th>(n, σ)</th><th>OC</th><th>OC</th></tr></thead><tbody><tr><td>(400, 0.2)</td><td>4.01</td><td>3.97</td></tr><tr><td>(400, 0.4)</td><td>8.11</td><td>8.02</td></tr><tr><td>(200, 0.2)</td><td>7.47</td><td>7.35</td></tr><tr><td>(200, 0.4)</td><td>19.51</td><td>18.94</td></tr><tr><td>(50, 0.2)</td><td>26.10</td><td>5.14</td></tr><tr><td>(50, 0.4)</td><td>37.36</td><td>9.49</td></tr></tbody></table>

# 1 6 Real data analysis

## 2 6.1 Influenza virus growth and migration model

3 We consider the ODE model introduced in Wu et. al [42] for the growth and migration of influenza

4 virus-specific effector CD8+ T cells, among lymph node ($T_E^m$), spleen ($T_E^s$), and lung ($T_E^l$) of mice. After

5 a model selection process, it turns out that the following model