Predrag Pilipovic
email: predrag@math.ku.dk

Adeline Samson
email: adeline.leclercq-samson@univ-grenoble-alpes.fr

Susanne Ditlevsen
email: susanne@math.ku.dk

PARAMETER ESTIMATION IN NONLINEAR MULTIVARIATE STOCHASTIC DIFFERENTIAL EQUATIONS BASED ON SPLITTING SCHEMES A PREPRINT

Keywords: Asymptotic normality, Consistency, L p convergence, Splitting schemes, Stochastic differential equations, Stochastic Lorenz system

Surprisingly, general estimators for nonlinear continuous time models based on stochastic differential equations are yet lacking. Most applications still use the Euler-Maruyama discretization, despite many proofs of its bias. More sophisticated methods, such as Kessler's Gaussian approximation, Ozak's Local Linearization, Aït-Sahalia's Hermite expansions, or MCMC methods, lack a straightforward implementation, do not scale well with increasing model dimension or can be numerically unstable. We propose two efficient and easy-to-implement likelihood-based estimators based on the Lie-Trotter (LT) and the Strang (S) splitting schemes. We prove that S has L p convergence rate of order 1, a property already known for LT. We show that the estimators are consistent and asymptotically efficient under the less restrictive one-sided Lipschitz assumption. A numerical study on the 3-dimensional stochastic Lorenz system complements our theoretical findings. The simulation shows that the S estimator performs the best when measured on precision and computational speed compared to the state-of-the-art.

Introduction

Stochastic differential equations (SDEs) are popular models for physical, biological, and socio-economic processes. Some recent applications include tipping points in the climate (Ditlevsen and Ditlevsen, 2023), the spread of COVID-19 [START_REF] Arnst | A hybrid stochastic model and its Bayesian identification for infectious disease screening in a university campus with application to massive COVID-19 screening at the University of Liège[END_REF][START_REF] Kareem | A Stochastic Differential Equations Model for Internal COVID-19 Dynamics[END_REF], animal movements [START_REF] Michelot | The Langevin diffusion as a continuous-time model of animal movement and habitat selection[END_REF][START_REF] Michelot | Varying-Coefficient Stochastic Differential Equations with Applications in Ecology[END_REF] and cryptocurrency rates [START_REF] Dipple | Using correlated stochastic differential equations to forecast cryptocurrency rates and social media activities[END_REF]. The advantage of SDEs is their ability to capture and quantify the randomness of the underlying dynamics. They are especially applicable when the dynamics are not entirely understood, and the unknown parts act as random. The following parametric form is common for an SDE model with additive noise: dX t = F (X t ; β) dt + Σ dW t , X 0 = x 0 .

(1)

We want to estimate the underlying drift parameter β and diffusion parameter Σ based on discrete observations of X t . The transition density is necessary for likelihood-based estimators and, thus, a closed-form solution to (1). However, the transition density is only available for a few SDEs, including the Ornstein-Uhlenbeck (OU) process, which has a linear drift function F. Extensive literature exists on MCMC methods for the nonlinear case [START_REF] Fuchs | Inference for Diffusion Processes with Applications in Life Sciences[END_REF][START_REF] Chopin | An Introduction to Sequential Monte Carlo[END_REF] however, these are often computationally intensive and do not always converge to the correct values for complex models. Thus, we need a valid approximation of the transition density to perform likelihood-based statistical inference.

The most straightforward discretization scheme is the Euler-Maruyama (EM) [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF]. Its main advantage is the easy-to-implement and intuitive Gaussian transition density. Both frequentist and Bayesian approaches extensively employ EM across theoretical and applied studies. However, the EM-based estimator has many disadvantages. First, it exhibits pronounced bias as the discretization step increases (see [START_REF] Florens-Zmirou | Approximate discrete-time schemes for statistics of diffusion processes[END_REF] for a theoretical study, or [START_REF] Gloaguen | Stochastic differential equation based on a multimodal potential to model movement data in ecology[END_REF], [START_REF] Gu | Parameter Estimation for Multivariate Nonlinear Stochastic Differential Equation Models: A Comparison Study[END_REF] for applied studies). Second, [START_REF] Hutzenthaler | Strong and weak divergence in finite time of Euler's method for stochastic differential equations with non-globally Lipschitz continuous coefficients[END_REF] showed that it is not mean-square convergent when the drift function F of (1) grows super-linearly. Consequently, we should avoid EM for models with polynomial drift. Third, it often fails to preserve important structural properties, such as hypoellipticity, geometric ergodicity, and amplitudes, frequencies, and phases of oscillatory processes [START_REF] Buckwar | A splitting method for SDEs with locally Lipschitz drift: Illustration on the FitzHugh-Nagumo model[END_REF].

Some pioneering papers on likelihood-based SDE estimators are [START_REF] Dacunha-Castelle | Estimation of the coefficients of a diffusion from discrete observations[END_REF]; [START_REF]On Estimating the Diffusion Coefficient[END_REF]; [START_REF] Florens-Zmirou | Approximate discrete-time schemes for statistics of diffusion processes[END_REF]; [START_REF] Genon-Catalot | On the estimation of the diffusion coefficient for multi-dimensional diffusion processes[END_REF]; [START_REF] Kessler | Estimation of an Ergodic Diffusion from Discrete Observations[END_REF]. The first two only estimate the diffusion parameter. [START_REF] Florens-Zmirou | Approximate discrete-time schemes for statistics of diffusion processes[END_REF] used EM to estimate both parameters and derived asymptotic properties. Genon-Catalot and Jacod (1993) generalized to higher dimensions, non-equidistant discretization step, and a generic form of the objective function, however only estimating the diffusion parameter. [START_REF] Kessler | Estimation of an Ergodic Diffusion from Discrete Observations[END_REF] proposed an estimator (denoted K) approximating the unknown transition density with a Gaussian density using the true conditional mean and covariance, or approximations thereof using the infinitesimal generator. He proved consistency and asymptotic normality under the commonly used, but too restrictive, global Lipschitz assumption on the drift function F.

A competitive likelihood-based approach relies on local linearization (LL), initially proposed by [START_REF] Ozaki | Statistical Identification of Storage Models with Application to Stochastic Hydrology[END_REF] and later extended by [START_REF] Ozaki | A Bridge between Nonlinear Time Series Models and Nonlinear Stochastic Dynamical Systems: A Local Linearization Approach[END_REF]; [START_REF] Shoji | Estimation for nonlinear stochastic differential equations by a local linearization method[END_REF]. They approximated the drift between two consecutive observations by a linear function. In the case of additive noise, this corresponds to an OU process with a known Gaussian transition density. Thus, the likelihood approximation is a product of Gaussian densities. [START_REF] Shoji | Approximation of Continuous Time Stochastic Processes by a Local Linearization Method[END_REF] proved that LL discretization is one-step consistent and L p convergent with order 1.5. [START_REF] Shoji | A note on convergence rate of a linearization method for the discretization of stochastic differential equations[END_REF], [START_REF] Jimenez | A weak Local Linearization scheme for stochastic differential equations with multiplicative noise[END_REF] extended the theory of LL for SDEs with multiplicative noise. Simulation studies show the superiority of the LL estimator compared to other estimators [START_REF] Shoji | Estimation for nonlinear stochastic differential equations by a local linearization method[END_REF][START_REF] Hurn | Seeing the Wood for the Trees: A Critical Evaluation of Methods to Estimate the Parameters of Stochastic Differential Equations[END_REF][START_REF] Gloaguen | Stochastic differential equation based on a multimodal potential to model movement data in ecology[END_REF][START_REF] Gu | Parameter Estimation for Multivariate Nonlinear Stochastic Differential Equation Models: A Comparison Study[END_REF].

Until recently, the implementation of the LL estimator was numerically ill-conditioned due to the possible singularity of the Jacobian matrix of the drift function F. However, [START_REF] Gu | Parameter Estimation for Multivariate Nonlinear Stochastic Differential Equation Models: A Comparison Study[END_REF] proposed an efficient implementation that overcomes this. The main disadvantage of the LL method is its slow computational speed. [START_REF] Aït-Sahalia | Maximum Likelihood Estimation of Discretely Sampled Diffusions: A Closed-form Approximation Approach[END_REF] proposed Hermite expansions (HE) to approximate the transition density, focusing on univariate time-homogeneous diffusions. This method, widely utilized in finance, was later extended to both reducible and irreducible multivariate diffusions [START_REF] Aït-Sahalia | Closed-form likelihood expansions for multivariate diffusions[END_REF]. [START_REF] Chang | On the approximate maximum likelihood estimation for diffusion processes[END_REF] found conditions under which the HE estimator has the same asymptotic distribution as the exact maximum likelihood estimator (MLE). [START_REF] Choi | Closed-form likelihood expansions for multivariate time-inhomogeneous diffusions[END_REF][START_REF] Choi | Explicit form of approximate transition probability density functions of diffusion processes[END_REF] further broadened the technique to time-inhomogeneous settings. [START_REF] Picchini | Practical estimation of high dimensional stochastic differential mixedeffects models[END_REF] used the method for multidimensional diffusions with random effects. When an SDE is irreducible, Aït-Sahalia (2008) applied Kolmogorov's backward and forward equations to develop a small-time expansion of the diffusion probability densities. [START_REF] Yang | A new delta expansion for multivariate diffusions via the Itô-Taylor expansion[END_REF] introduced a delta expansion method, using Itô-Taylor expansions to derive analytical approximations of the transition densities of multivariate diffusions inspired by [START_REF] Aït-Sahalia | Maximum Likelihood Estimation of Discretely Sampled Diffusions: A Closed-form Approximation Approach[END_REF]. While Aït-Sahalia's approach allows for a broad class of drift and diffusion functions, the implementation can be complex. To our knowledge, there have not been any applications to models with more than four dimensions. Furthermore, computing coefficients even up to order two can be challenging, while higher-order approximations are often necessary for non-linear models. [START_REF] Hurn | Seeing the Wood for the Trees: A Critical Evaluation of Methods to Estimate the Parameters of Stochastic Differential Equations[END_REF] implemented HE up to third order in univariate cases, emphasizing the importance of symbolic computation tools like Mathematica or Maple. Their survey concluded that while LL is the best among discrete maximum likelihood estimators, HE is the preferred overall choice. They highlighted that the HE proposed by [START_REF] Aït-Sahalia | Maximum Likelihood Estimation of Discretely Sampled Diffusions: A Closed-form Approximation Approach[END_REF] has the best trade-off between speed and accuracy, proving more feasible than LL in most financial applications. This finding aligns with the newer review study from [START_REF] López-Pérez | Parametric Estimation of Diffusion Processes: A Review and Comparative Study[END_REF]. However, LL's broad applicability contrasts with the limitations of Hermite expansions, particularly for high-dimensional multivariate models exceeding three dimensions.

Apart from the above-mentioned general methods, there are some specific setups. [START_REF] Sørensen | Small-diffusion asymptotics for discretely sampled stochastic differential equations[END_REF] investigated a small-diffusion estimator, [START_REF] Ditlevsen | Inference for observations of integrated diffusion processes[END_REF]; [START_REF] Gloter | Parameter Estimation for a Discretely Observed Integrated Diffusion Process[END_REF] worked with integrated diffusion, and [START_REF] Uchida | Adaptive estimation of an ergodic diffusion process based on sampled data[END_REF] used adaptive maximum likelihood estimation. [START_REF] Bibby | Martingale estimation functions for discretely observed diffusion processes[END_REF] and [START_REF] Forman | The Pearson diffusions: A class of statistically tractable diffusion processes[END_REF] explored martingale estimation functions (EF) in one-dimensional diffusions, but they are difficult to extend to multidimensional SDEs. [START_REF] Ditlevsen | Hypoelliptic diffusions: filtering and inference from complete and partial observations[END_REF] used the 1.5 scheme to solve the problem of hypoellipticity when the diffusion matrix is not of full rank.

More recently, contributions from [START_REF] Gloter | Parameter Estimation for a Discretely Observed Integrated Diffusion Process[END_REF]Yoshida (2020, 2021) have extended the research of [START_REF] Uchida | Adaptive estimation of an ergodic diffusion process based on sampled data[END_REF]. [START_REF] Gloter | Adaptive and non-adaptive estimation for degenerate diffusion processes[END_REF] introduced a non-adaptive approach and offered similar analytic asymptotic results as [START_REF] Ditlevsen | Hypoelliptic diffusions: filtering and inference from complete and partial observations[END_REF] without imposing strict limitations on the model class. [START_REF] Iguchi | Parameter Estimation with Increased Precision for Elliptic and Hypo-elliptic Diffusions[END_REF] proposed sampling schemes for elliptic and hypoelliptic models that often result in conditionally non-Gaussian integrals, distinguishing their approach from prior works. As the transition density of their new scheme is typically complex, [START_REF] Iguchi | Parameter Estimation with Increased Precision for Elliptic and Hypo-elliptic Diffusions[END_REF] created a closed-form density expansion using Malliavin calculus. They recommended a transition density scheme that retained second-order precision through prudent truncation of the expansion. This closed-form expansion aligns with the works of [START_REF] Aït-Sahalia | Maximum Likelihood Estimation of Discretely Sampled Diffusions: A Closed-form Approximation Approach[END_REF][START_REF] Aït-Sahalia | Closed-form likelihood expansions for multivariate diffusions[END_REF] and [START_REF] Li | Maximum-likelihood estimation for diffusion processes via closed-form density expansions[END_REF] on elliptic SDEs, although with a different approach. [START_REF] Iguchi | Parameter Estimation with Increased Precision for Elliptic and Hypo-elliptic Diffusions[END_REF] deliver asymptotic results with analytically available rates, beneficial for both elliptic and hypoelliptic models. An essential aspect of any estimator is the practical execution in real-world applications. Although the previously mentioned research contributes significantly to the theoretical development and broadens our understanding of inference for SDEs, its practical implementations tend not to be user-friendly. Except for precomputed models, applications by non-specialists can be challenging. Our main contribution is proposing estimators that are intuitive, easy to implement, computationally efficient, and scalable with increasing dimensions. These characteristics make the estimators accessible to researchers in various applied sciences while maintaining desirable statistical properties. Moreover, these estimators remain competitive with the best state-of-the-art methods, particularly concerning estimation bias and variance.

We propose to use the LT or the S splitting schemes for statistical inference. These numerical approximations were first suggested for ordinary differential equations (ODEs) (see for example, [START_REF] Mclachlan | Splitting methods[END_REF]; [START_REF] Blanes | Splitting and composition methods in the numerical integration of differential equations[END_REF]), but their extension to SDEs is straightforward. A few studies have investigated numerical properties [START_REF] Bensoussan | Approximation of Some Stochastic Differential Equations by the Splitting Up Method[END_REF][START_REF] Ableidinger | A Stochastic Version of the Jansen and Rit Neural Mass Model: Analysis and Numerics[END_REF][START_REF] Ableidinger | Splitting Integrators for the Stochastic Landau-Lifshitz Equation[END_REF][START_REF] Buckwar | A splitting method for SDEs with locally Lipschitz drift: Illustration on the FitzHugh-Nagumo model[END_REF]. [START_REF] Barbu | A Product Formula Approach to Nonlinear Optimal Control Problems[END_REF] applied LT splitting on nonlinear optimal control problems, while [START_REF] Hopkins | Lie-Trotter Product Formulas for Nonlinear Filtering[END_REF] used it for nonlinear filtering. Bou-Rabee and Owhadi (2010); [START_REF] Abdulle | Long Time Accuracy of Lie-Trotter Splitting Methods for Langevin Dynamics[END_REF] used LT splitting to investigate conditions for preserving the measure of the ergodic nonlinear Langevin equations. Recently, [START_REF] Bréhier | Analysis of a Positivity-preserving Splitting Scheme for Some Nonlinear Stochastic Heat Equations[END_REF] showed that LT splitting successfully preserved positivity for a class of nonlinear stochastic heat equations with multiplicative space-time white noise. Additional studies on the application of splitting schemes to SDEs include those by [START_REF] Misawa | A lie algebraic approach to numerical integration of stochastic differential equations[END_REF]; [START_REF] Milstein | Quasi-symplectic methods for Langevin-type equations[END_REF]; [START_REF] Leimkuhler | Molecular dynamics[END_REF]; [START_REF] Alamo | A Technique for Studying Strong and Weak Local Errors of Splitting Stochastic Integrators[END_REF]; [START_REF] Bréhier | Analysis of some splitting schemes for the stochastic Allen-Cahn equation[END_REF].

Regarding statistical applications, to the best of our knowledge, only [START_REF] Buckwar | Spectral density-based and measure-preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs[END_REF]; Ditlevsen et al. (2023) used splitting schemes for parametric inference in combination with Approximate Bayesian Computation, and Ditlevsen and Ditlevsen (2023) used it for prediction of a forthcoming collapse in the climate.

This paper presents five main contributions:

1. We introduce two new efficient, easy-to-implement, and computationally fast estimators for multidimensional nonlinear SDEs. 2. We establish L p convergence of the S splitting scheme. 3. We prove consistency and asymptotic normality of the new estimators under the less restrictive assumption of one-sided Lipschitz. This proof requires innovative approaches. 4. We demonstrate the estimators' performance in a stochastic version of the chaotic Lorenz system, in contrast to prior studies that primarily addressed the deterministic Lorenz system. 5. We compare the new estimators to three discrete maximum likelihood estimators from the literature in a simulation study, comparing the accuracy and computational speed.

The rest of this paper is structured as follows. In Section 2 we introduce the SDE model class and define the splitting schemes and the estimators. In Section 3, we show that the S splitting has better one-step predictions than the LT, and we prove that the S splitting is L p consistent with order 1.5 and L p convergent with order 1.

K up to order J General J fixed: h → 0, N h → ∞, N h p → 0,
for any p ∈ N a one-sided global Lipschitz assumption. We illustrate in Section 6 the theoretical results in a simulation study on a model that is not globally Lipschitz, the 3-dimensional stochastic Lorenz systems. Since the objective functions based on pseudo-likelihoods are multivariate in both data and parameters, we use automatic differentiation (AD) to get faster and more reliable estimators. We compare the precision and speed of the EM, K, LL, LT, and S estimators. We show that the EM and LT estimators become biased before the others with increasing discretization step h and that the LL and S perform the best. However, S is much faster than LL because LL calculates a new covariance matrix for each combination of data points and parameter values.

Notation. We use capital bold letters for random vectors, vector-valued functions, and matrices, while lowercase bold letters denote deterministic vectors. ∥ • ∥ denotes both the L 2 vector norm in R d and the matrix norm induced by the L 2 norm, defined as the square root of the largest eigenvalue. Superscript (i) on a vector denotes the i-th component, while on a matrix it denotes the i-th row. Double subscript ij on a matrix denotes the component in the i-th row and j-th column. If a matrix is a product of more matrices, square brackets with subscripts denote a component inside the matrix. The transpose is denoted by ⊤. Operator Tr(•) returns the trace of a matrix and det(•) the determinant. Sometimes, we denote by [a i ] d i=1 a vector with coordinates a i , and by [b ij ] d i,j=1 a matrix with coordinates b ij , for i, j = 1, . . . , d. We denote with ∂ i g(x) the partial derivative of a generic function g : R d → R with respect to x (i) and ∂ 2 ij g(x) the second partial derivative. The nabla operator ∇ denotes the gradient vector of a function g, ∇g(

x) = [∂ i g(x)] d i=1 . The differential operator D denotes the Jacobian matrix DF(x) = [∂ i F (j) (x)] d i,j=1
, for a vector-valued function

F : R d → R d . H denotes the Hessian matrix of a real-valued function g, H g (x) = [∂ ij g(x)] d i,j=1
. Let R represent a vector (or a matrix) valued function defined on (0, 1) × R d , such that, for some constant C, ∥R(a, x)∥ < aC(1 + ∥x∥) C for all a, x. When denoted R, it is a scalar. ----→ for convergence in probability, distribution, and almost surely, respectively. I d denotes the d-dimensional identity matrix, while 0 d×d is a d-dimensional zero square matrix. For an event E ∈ F, we denote by ⊮ E the indicator function.

The

Problem setup

Let X in (1) be defined on a complete probability space (Ω, F, P θ ) with a complete right-continuous filtration (F t ) t≥0 , and let the d-dimensional Wiener process W = (W t ) t≥0 be adapted to F t . The probability measure P θ is parameterized by the parameter θ = (β, Σ). Rewrite equation (1) as follows: The matrix ΣΣ ⊤ is assumed to be positive definite and determines the variance of the process. Since any square root of ΣΣ ⊤ induces the same distribution, Σ is only identifiable up to equivalence classes. Thus, instead of estimating Σ, we estimate ΣΣ ⊤ . The drift function F in (1) is split up into a linear part given by matrix A and vector b and a nonlinear part given by N. This decomposition is essential for defining the splitting schemes and the objective functions used for estimating θ.

dX t = A(β)(X t -b(β)) dt + N (X t ; β) dt + Σ dW t , X 0 = x 0 , (2)
such that F(x; β) = A(β)(x -b(β)) + N (x; β). Let Θ = Θ β × Θ Σ be
We denote the true parameter value by θ 0 = (β 0 , Σ 0 ) and assume that θ 0 ∈ Θ. Sometimes we write A 0 , b 0 , N 0 (x) and ΣΣ ⊤ 0 instead of A(β 0 ), b(β 0 ), N(x; β 0 ) and Σ 0 Σ ⊤ 0 , when referring to the true parameters. We write A, b, N(x) and ΣΣ ⊤ for any parameter θ. Sometimes, we suppress the parameter to simplify notation. For example, E implicitly refers to E θ .

Remark The drift function F(x) can always be rewritten as A(xb) + N(x) for any A, b by setting N(x) = F(x) -A(xb), including choosing A and b to be zero. In this case, the splitting proposed below will result in a Brownian motion (3) and a nonlinear ODE (4).

Remark We assume additive noise, meaning that the diffusion matrix does not depend on the current state. While this assumption is natural in some applications, it can be restrictive in others. The proposed methodology could potentially be extended to reducible diffusions by applying the Lamperti transform to obtain a unit diffusion coefficient, as demonstrated by [START_REF] Aït-Sahalia | Closed-form likelihood expansions for multivariate diffusions[END_REF]. However, if the transform depends on the parameter, estimation is not straightforward. In this paper, we only consider additive noise.

Assumptions

The main assumption is that (2) has a unique strong solution X= (X t ) t∈[0,T ] , adapted to (F t ) t∈[0,T ] , which follows from the following first two assumptions (Theorem 2 in [START_REF] Alyushina | Euler Polygonal Lines for Itô Equations with Monotone Coefficients[END_REF], Theorem 1 in [START_REF] Krylov | A Simple Proof of the Existence of a Solution of Itô's Equation with Monotone Coefficients[END_REF], Theorem 3.5 in [START_REF] Mao | Stochastic differential equations and applications[END_REF]). We need the last three assumptions to prove the properties of the estimators.

(A1) Function N is twice continuously differentiable with respect to x and θ, i.e., N ∈ C 2 . Additionally, it is one-sided globally Lipschitz continuous with respect to x on R d × Θ β , i.e., there exists a constant C > 0 such that:

(x -y) ⊤ (N(x; β) -N(y; β)) ≤ C∥x -y∥ 2 , ∀x, y ∈ R d .
(A2) Function N grows at most polynomially in x, uniformly in θ, i.e., there exist constants C > 0 and χ ≥ 1 such that:

∥N (x; β) -N (y; β) ∥ 2 ≤ C 1 + ∥x∥ 2χ-2 + ∥y∥ 2χ-2 ∥x -y∥ 2 , ∀x, y ∈ R d .
Additionally, its derivatives are of polynomial growth in x, uniformly in θ.

(A3) The solution X of SDE (1) has invariant probability ν 0 (dx).

(A4) ΣΣ ⊤ is invertible on Θ Σ . (A5) Function F is identifiable in β, i.e., if F(x, β 1 ) = F(x, β 2 ) for all x ∈ R d , then β 1 = β 2 .
Assumption (A3) is required for the ergodic theorem to ensure convergence in distribution. Assumption (A4) implies the ellipticity of model (1), which is not needed for the S estimator. On the contrary, the EM estimator breaks down in hypoelliptic models. We will treat the hypoelliptic case in a separate paper where the proofs are more involved. Assumption (A5) ensures the identifiability of the parameter.

Assume a sample (X t k ) N k=0 ≡ X 0:t N from (2) at time steps 0 = t 0 < t 1 < • • • < t N = T , which we, for notational simplicity, assume equidistant with step size h = t k -t k-1 .

Moments

Assumption (A1) ensures finiteness of the moments of the solution X [START_REF] Tretyakov | A Fundamental Mean-Square Convergence Theorem for SDEs with Locally Lipschitz Coefficients and Its Applications[END_REF], i.e.,

E[ sup t∈[0,T ] ∥X t ∥ 2p ] < C(1 + ∥x 0 ∥ 2p ), ∀ p ≥ 1.
Furthermore, we need the infinitesimal generator L of (1) defined on sufficiently smooth functions g : R d × Θ → R given by:

L θ0 g (x; θ) = F (x; β 0 ) ⊤ ∇g (x; θ) + 1 2 Tr(ΣΣ ⊤ 0 H g (x; θ)).
The moments of SDE (1) are expanded using the following lemma (Lemma 1.10 in [START_REF] Sørensen | Estimating functions for diffusion-type processes, chapter 1[END_REF]).

Lemma 2.1 Let Assumptions (A1)-(A2) hold. Let X be a solution of (1). Let g ∈ C (2l+2) be of polynomial growth and p ≥ 2. Then,

E θ0 [g(X t k ; θ) | F t k-1 ] = l j=0 h j j! L j θ0 g(X t k-1 ; θ) + R(h l+1 , X t k-1 ).
We need terms up to order R(h 3 , X t k-1 ). After applying the generator L θ on g(x) = x (i) , the previous Lemma yields:

E[X (i) t k | X t k-1 = x] = x (i) + hF (i) (x) + h 2 2 (F(x) ⊤ ∇F (i) (x) + 1 2 Tr(ΣΣ ⊤ H F (i) (x))) + R(h 3 , x).

Splitting Schemes

Consider the following splitting of (2):

dX [1] t = A(X [1] t -b) dt + Σ dW t , X [1] 0 = x 0 , (3)
dX [2] t = N(X [2] t ) dt, X [2] 0 = x 0 . (4)
The solution of equation ( 3) is an OU process given by the following h-flow:

X [1] t k = Φ [1] h (X [1] t k-1 ) = e Ah X [1] t k-1 +(I -e Ah )b + ξ h,k , (5)
where ξ h,k i.i.d

∼ N d (0, Ω h ) for k = 1, . . . , N [START_REF] Vatiwutipong | Alternative way to derive the distribution of the multivariate Ornstein-Uhlenbeck process[END_REF]. The covariance matrix Ω h and the conditional mean of the OU process (5) are provided by:

Ω h = h 0 e A(h-u) ΣΣ ⊤ e A ⊤ (h-u) du = hΣΣ ⊤ + h 2 2 (AΣΣ ⊤ + ΣΣ ⊤ A ⊤ ) + R(h, x 0 ), (6)
µ h (x; β) := e A(β)h x + (I -e A(β)h )b(β).
(7) Assumptions (A1) and (A2) ensure the existence and uniqueness of the solution of (4) (Theorem 1.2.17 in [START_REF] Humphries | Deterministic and random dynamical systems: theory and numerics[END_REF]). Thus, there exists a unique function f h : R d × Θ β → R d , for h ≥ 0, such that:

X [2] t k = Φ [2] h (X [2] t k-1 ) = f h (X [2] t k-1 ; β).
(8) For all β ∈ Θ β , the time flow f h fulfills the following semi-group properties:

f 0 (x; β) = x, f t+s (x; β) = f t (f s (x; β); β), t, s ≥ 0. ( 9
)
Remark Since only one-sided Lipschitz continuity is assumed, the solution to (4) might not exist for all h < 0 and all x 0 ∈ R d , implying that the inverse f -1 h might not exist. If it exists, then f -1 h = f -h . For the S estimator, we need a well-defined inverse. This is not an issue when N is globally Lipschitz .

We, therefore, introduce the following and last assumption.

(A6) Function f -1 h (x; β) is defined asymptotically, for all x ∈ R d , β ∈ Θ β , when h → 0. Before defining the splitting schemes, we present a useful proposition for expending the nonlinear solution f h (Section 1.8 in [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF]).

Proposition 2.2 Let Assumptions (A1)-(A2) hold. When h → 0, the h-flow of (4) is f h (x) = x + hN(x) + h 2 2 (DN(x)) N(x) + R(h 3 , x).
Now, we introduce the two most common splitting approximations, which serve as the main building blocks for the proposed estimators.

Definition 2.3 Let Assumptions (A1) and (A2) hold. The Lie-Trotter and Strang splitting approximations of the solution of (2) are given by:

X [LT] t k := Φ [LT] h (X [LT] t k-1 ) = (Φ [1] h • Φ [2] h )(X [LT] t k-1 ) = µ h (f h (X [LT] t k-1 )) + ξ h,k , (10)

X

[S]

t k := Φ [S] h (X [S] t k-1 ) = (Φ [2] h/2 • Φ [1] h • Φ [2] h/2 )(X [S] t k-1 ) = f h/2 (µ h (f h/2 (X [S] t k-1 )) + ξ h,k ). (11)
Remark The order of composition in the splitting schemes is not unique. Changing the order in the S splitting leads to a sum of 2 independent random variables, one Gaussian and one non-Gaussian, whose likelihood is not trivial. Thus, we only use the splitting (11). The reversed order in the LT splitting can be treated the same way as the S splitting.

Remark Splitting the drift F(x) into a linear and a nonlinear part is not unique. However, all theorems and properties, particularly consistency and asymptotic normality of the estimators, hold for any splitting choice. Yet, for fixed step size h and sample size N , certain splittings perform better than others. In this paper, we present two general and intuitive strategies. The first applies when the system has a fixed point; here, the linear part of the splitting is the linearization around the fixed point. The linear OU performs accurately near the fixed point, with the nonlinear part correcting for nonlinear deviations. Simulations consistently show this approach to perform best. Another strategy is to linearize around the measured average value for each coordinate. An in-depth analysis of the splitting strategies for a specific example is provided in Section 2.5.

Remark Overall trajectories of the S and LT splittings coincide up to the first h/2 and the last h/2 move of the flow Φ [2] h/2 . Indeed, when applied k times, the S splitting can be written as:

(Φ [S] h ) k (x 0 ) = (Φ [2] h/2 • (Φ [LT] h ) k • Φ [2] -h/2 )(x 0
). Thus, it is natural that LT and S have the same order of L p convergence. We prove this in Section 3. However, the LT and S trajectories differ in their output points ( 10) and ( 11). The S splitting outputs the middle points of the smooth steps of the deterministic flow (8), while the LT splitting outputs the stochastic increments in the rough steps. We conjecture that this is one of the reasons why the S splitting has superior statistical properties.

Estimators

In this section, we first introduce two new estimators, LT and S, given a sample X 0:t N . Subsequently, we provide a brief overview of the EM, K, and LL estimators, which will be compared in the simulation study.

Splitting estimators

The LT scheme (10) follows a Gaussian distribution. Consequently, the objective function corresponds to (twice) the negative pseudo-log-likelihood:

L [LT] (X 0:t N ; θ) θ = N log(det Ω h (θ)) + N k=1 (X t k -µ h (f h (X t k-1 ; β); β)) ⊤ Ω h (θ) -1 (X t k -µ h (f h (X t k-1 ; β); β)). ( 12
)
The S splitting ( 11) is a nonlinear transformation of the Gaussian random variable

µ h (f h/2 (X t k-1 ; β); β) + ξ h,k .
We first define:

Z t k (β) := f -1 h/2 (X t k ; β) -µ h (f h/2 (X t k ; β); β). (13)
Afterwards, we apply a change of variables to derive the following objective function:

L [S] (X 0:t N ; θ) θ = N log(det Ω h (θ)) + N k=1 Z t k (β) ⊤ Ω h (θ) -1 Z t k (β) -2 N k=1 log | det Df -1 h/2 (X t k ; β)|. ( 14
)
The last term is due to the nonlinear transformation and is an extra term that does not appear in commonly used pseudo-likelihoods.

The inverse function f -1 h may not exist for all parameters in the search domain of the optimization algorithm. However, this problem it can often be solved numerically. When f -1 h is well defined, we use the identity 14) to increase the speed and numerical stability. Finally, we define the estimators as:

-log | det Df -1 h (x; β) | = log | det Df h (x; β) | in (
θ [k] N := arg min θ L [k] (X 0:t N ; θ) , k ∈ {LT, S}. (15)

Euler-Maruyama

The EM method uses first-order Taylor expansion of (1):

X [EM] t k := X [EM] t k-1 + hF(X [EM] t k-1 ; β) + ξ [EM] h,k , (16)
where ξ [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF]. The transition density p

[EM] h,k i.i.d. ∼ N d (0, hΣΣ ⊤ ) for k = 1, . . . , N (
[EM] (X t k | X t k-1 ; θ)
is Gaussian, so the pseudo-likelihood follows trivially.

Kessler

The K estimator uses Gaussian transition densities p [K] (X t k | X t k-1 ; θ) with the true mean and covariance of the solution X [START_REF] Kessler | Estimation of an Ergodic Diffusion from Discrete Observations[END_REF]. When the moments are unknown, they are approximated using the infinitesimal generator (Lemma 2.1). We implement the estimator K based on the 2nd-order approximation:

X [K] t k := X [K] t k-1 + hF(X [K] t k-1 ; β) + ξ [K] h,k (X [K] t k-1 ) + h 2 2 DF(X [K] t k-1 ; β)F(X [K] t k-1 ; β) + 1 2 [Tr(ΣΣ ⊤ H F (i) (X [K] t k-1 ; β))] d i=1 , (17)
where ξ

[K] h,k (X [K] t k-1 ) ∼ N d (0, Ω [K] h,k (θ)),
and

Ω [K] h,k (θ) = hΣΣ ⊤ + h 2 2 (DF(X [K] t k-1 ; β)ΣΣ ⊤ + ΣΣ ⊤ D ⊤ F(X [K] t k-1 ; β)
). The covariance matrix is not constant which makes the algorithm slower for a larger sample size.

Ozaki's local linearization

Ozaki's LL method approximates the drift of (1) between consecutive observations by a linear function [START_REF] Jimenez | Simulation of Stochastic Differential Equations Through the Local Linearization Method. A Comparative Study[END_REF]. The LL method consists of the following steps:

(1) Perform LL of the drift F in each time interval [t, t + h) by the Itô-Taylor series;

(2) Compute the analytic solution of the resulting linear SDE.

The approximation becomes:

X [LL] t k := X [LL] t k-1 + Φ [LL] h (X [LL] t k-1 ; θ) + ξ [LL] h,k (X [LL] t k-1 ), (18)
where

Ω [LL] h,k (θ) = h 0 e DF(X [LL] t k-1 ;β)(h-u) ΣΣ ⊤ e DF(X [LL] t k-1 ;β) ⊤ (h-u) du and ξ [LL] h,k (X [LL] t k-1 ) ∼ N d (0, Ω [LL] h,k (θ)). More- over, Φ [LL] h (x; θ) = R h,0 (DF(x; β))F(x; β) + (hR h,0 (DF(x; β)) -R h,1 (DF(x; β)))M(x; θ), R h,i (DF(x; β)) = h 0 exp(DF(x; β)u)u i du, i = 0, 1, M(x; θ) = 1 2 (Tr H 1 (x; θ), Tr H 2 (x; θ), ..., Tr H d (x; θ)) ⊤ , H k (x; θ) = [ΣΣ ⊤ ] ij ∂ 2 F (k) ∂x (i) ∂x (j) (x) d i,j=1
.

Building on the approach by [START_REF] Gu | Parameter Estimation for Multivariate Nonlinear Stochastic Differential Equation Models: A Comparison Study[END_REF], we can efficiently compute R h,i and Ω

[LL] h,k (θ) using the following procedure. To begin, let us define three block matrices:

P 1 (x) = 0 d×d I d 0 d×d DF(x; β) , P 2 (x) = -DF(x; β) I d 0 d×d 0 d×d 0 d×d I d 0 d×d 0 d×d 0 d×d , P 3 (x) = DF(x; β) ΣΣ ⊤ 0 d×d -DF(x; β) ⊤ . (19)
Then, we compute the matrix exponential of matrices hP 1 (x) and hP 2 (x):

exp(hP 1 (x)) = ⋆ R h,0 (DF(x; β)) 0 d×d ⋆ , exp(hP 2 (x)) = ⋆ ⋆ B R h,1 (DF(x; β)) 0 d×d ⋆ ⋆ 0 d×d 0 d×d ⋆ .
Starting with the first matrix, we derive R h,0 (DF(x; β)). Then, we compute R h,1 (DF(x; β)) using the formula

R h,1 (DF(x; β)) = exp(hDF(x; β))B R h,1 (DF(x; β)).
The terms marked with ⋆ symbols can be disregarded.

Finally, we obtain Ω [LL]

h,k (θ) from the matrix exponential:

exp(hP 3 (x)) = B Ω h,k (DF(x; β); θ) C Ω h,k (DF(x; β); θ) 0 d×d ⋆ , Ω [LL] h,k (θ) = C Ω h,k (DF(x; β); θ)B Ω h,k (DF(x; β); θ) ⊤ .
Thus, we have a Gaussian density p [LL] (X t k | X t k-1 ; θ) and standard likelihood inference. Like in the case of K, the covariance matrix Ω

[LL]

h,k (θ) depends on the previous state X

[LL] t k-1 , which is a major downside since it is harder to implement and slower to run due to the computation of N -1 covariance matrices. Unlike K, LL does not use Taylor expansions of the approximated drift and covariance matrix, so the influence of the sample size N on computational times is much stronger. For details on the derivations of the previous formulas, see [START_REF] Gu | Parameter Estimation for Multivariate Nonlinear Stochastic Differential Equation Models: A Comparison Study[END_REF].

An example: the stochastic Lorenz system

The Lorenz system is a 3D system introduced by [START_REF] Lorenz | Deterministic Nonperiodic Flow[END_REF] to model atmospheric convection. The model is originally deterministic exhibiting deterministic chaos. It means that tiny differences in initial conditions lead to unpredictable and widely diverging trajectories. The Lorenz system evolves around two strange attractors. It means that the trajectories remain within some bounded region, while points that start in close proximity may eventually separate by arbitrary distances as time progresses. [START_REF] Hilborn | Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers[END_REF]. We add noise to include unmodelled forces and randomness in the Lorenz system . The stochastic Lorenz system is given by:

dX t = p(Y t -X t ) dt + σ 1 dW (1) t , dY t = (rX t -Y t -X t Z t ) dt + σ 2 dW (2) t , dZ t = (X t Y t -cZ t ) dt + σ 3 dW (3) t .
(20)

The variables X t , Y t , and Z t represent convective intensity, and horizontal and vertical temperature differences, respectively. Parameters p, r, and c denote the Prandtl number, the Rayleigh number, and a geometric factor, respectively [START_REF] Tabor | Chaos and Integrability in Nonlinear Dynamics: An Introduction[END_REF]. [START_REF] Lorenz | Deterministic Nonperiodic Flow[END_REF] used the values p = 10, r = 28 and c = 8/3, yielding chaotic behavior.

The system does not fulfill the global or the one-sided Lipschitz condition because it is a second-order polynomial [START_REF] Humphries | Runge-Kutta Methods for Dissipative and Gradient Dynamical Systems[END_REF]. However, it has a unique global solution and an invariant probability [START_REF] Keller | Attractors and Bifurcations of the Stochastic Lorenz System[END_REF]. Thus, all assumptions (A2)-(A5), except (A1) hold. Even so, we show in Section 6 that the estimators work.

Different approaches for estimating parameters in the Lorenz system have been proposed, mostly in the deterministic case. [START_REF] Zhuang | Parameter Estimation of Lorenz Chaotic System Based on a Hybrid Jaya-Powell Algorithm[END_REF] and [START_REF] Lazzús | Parameter estimation of Lorenz chaotic system using a hybrid swarm intelligence algorithm[END_REF] used sophisticated optimization algorithms to achieve better precision. Dubois et al. (2020) and [START_REF] Ann | Parameter Estimation of Lorenz Attractor: A Combined Deep Neural Network and K-Means Clustering Approach[END_REF] used deep neural networks in combination with other machine learning algorithms. [START_REF] Ozaki | The Role of the Likelihood Function in the Estimation of Chaos Models[END_REF] used Kalman filtering based on LL on the stochastic Lorenz system. Figure 1 shows an example trajectory of the stochastic Lorenz system. The trajectory was generated by subsampling from an EM simulation, such that N = 10000 and h = 0.05, with parameter values p = 10, r = 28, c = 8/3, σ 2 1 = 1, σ 2 2 = 2 and σ 2 3 = 1.5. Even if the trajectory had not been stochastic, the unpredictable jumps in the first row of Figure 1 would still have been there due to the chaotic behavior .

We suggest to split SDE (20) by choosing the OU part (3) as the linearization around one of the two fixed points (x ⋆ , y ⋆ , z ⋆ ) = (± c(r -1), ± c(r -1), r -1). For simplicity, we exclude the fixed point (0, 0, 0) since X and Y spend little time around this point, see Figure 1. Specifically, we apply a mixture of two splittings, linearizing around ( c(r -1), c(r -1), r -1) when X > 0 and around (-c(r -1), -c(r -1), r -1) when X < 0. We denote these estimators by LT mix and S mix . The splitting is given by:

A mix = -p p 0 1 -1 -x ⋆ y ⋆ x ⋆ -c , b mix = x ⋆ y ⋆ z ⋆ , N mix (x, y, z) = 0 -(x -x ⋆ )(z -z ⋆ ) (x -x ⋆ )(y -y ⋆ )
.

The OU process is mean-reverting towards b mix = (x ⋆ , y ⋆ , z ⋆ ). The nonlinear solution is

f mix,h (x, y, z) = x (y -y ⋆ ) cos(h(x -x ⋆ )) -(z -z ⋆ ) sin(h(x -x ⋆ )) + y ⋆ (y -y ⋆ ) sin(h(x -x ⋆ )) + (z -z ⋆ ) cos(h(x -x ⋆ )) + z ⋆ .
The solution is a composition of a 3D rotation and translation of (y, z) around the fixed point. The inverse always exists, and thus, Assumption (A6) holds. Moreover, det Df -1 mix,h (x, y, z) = 1. The mixing strategy does not increase the complexity of the implementation significantly, and it is straightforward to incorporate into the existing framework. Thus, this splitting strategy is convenient when the model has several fixed points.

An alternative splitting linearizes around the average of the observations. Let (µ x , µ x , µ z ) be the average of the data, where we put µ x = µ y since the difference of their averages is small, around 10 -3 . We denote these estimators by LT avg and S avg .The splitting is given by:

A avg = -p p 0 r -µ z -1 -µ x µ x µ x -c , b avg = µ x µ x µ z , N avg (x, y, z) =   0 -(x -µ x )(z -µ z ) + (r -1 -µ z )µ x (x -µ x )(y -µ x ) + µ 2 x -cµ z   .
The nonlinear solution is:

f avg,h (x, y, z) =    µ x µ x + cµz-µ 2 x x-µx µ z + µx(r-1-µz) x-µx    +    x -µ x (y -µ x - cµz-µ 2 x x-µx ) cos(h(x -µ x )) -(z -µ z -µx(r-1-µz) x-µx ) sin(h(x -µ x )) (y -µ x - cµz-µ 2 x x-µx ) sin(h(x -µ x )) + (z -µ z -µx(r-1-µz) x-µx ) cos(h(x -µ x ))    , where we define f avg,h (µ x , y, z) = (µ x , y + hµ x (r -1 -µ z ), z + hµ 2 x -cµ z ) ⊤ . Again, det Df -1 avg,h (x, y, z) = 1.
3 Order of one-step predictions and L p convergence

In this Section, we investigate L p convergence of the splitting schemes and the order of the one-step predictions. Theorem 2.1 in Tretyakov and Zhang (2013) extends Milstein's fundamental theorem on L p convergence for global Lipschitz coefficients [START_REF] Milstein | A Theorem on the Order of Convergence of Mean-Square Approximations of Solutions of Systems of Stochastic Differential Equations[END_REF] to Assumptions (A1) and (A2). This theorem provides the theoretical underpinning for our approach, drawing on the key concepts of L p consistency and boundedness of moments.

Definition 3.1 (L p consistency of a numerical scheme) The one-step approximation Φ h of the solution X is L p consistent, p ≥ 1, of order q 2 -1/2≥ 0, if for k = 1, . . . , N , and some q 1 ≥ q 2 + 1/2:

∥E[X t k -Φ h (X t k-1 ) | X t k-1 = x]∥ = R(h q1 , x), (E[∥X t k -Φ h (X t k-1 )∥ 2p | X t k-1 = x]) 1 2p = R(h q2 , x),
Definition 3.2 (Bounded moments of a numerical scheme) A numerical approximation X of the solution X has bounded moments, if for all p ≥ 1,there exists constant C > 0, such that, for k = 1, . . . , N :

E[∥ X t k ∥ 2p ] ≤ C(1 + ∥x 0 ∥ 2p ).
The following theorem (Theorem 2.1 in Tretyakov and Zhang ( 2013)) gives sufficient conditions for L p convergence of a numerical scheme in a one-sided Lipschitz framework.

Theorem 3.3 (L p convergence of a numerical scheme) Let Assumptions (A1) and (A2) hold, and let X t k be a numerical approximation of the solution X t k of (1) at time t k . If

(1) The one-step approximation

X t k = Φ h ( X t k-1 ) is L p consistent of order q 2 -1/2; and
(2) X has bounded moments, then, the numerical method X is L p convergent, p ≥ 1, of order q 2 -1/2, i.e., for k = 1, . . . , N , it holds:

(E[∥X t k -X t k ∥ 2p ]) 1 2p = R(h q2-1/2 , x 0 ).

Lie-Trotter splitting

We first show that the one-step LT approximation is of order R(h 2 , x 0 ) in mean. The following proposition is proved in the Supplementary Material [START_REF] Pilipovic | Supplement to "Parameter Estimation in Nonlinear Multivariate Stochastic Differential Equations Based on Splitting Schemes[END_REF] for scheme (10), as well as for the reversed order of composition. We demonstrate that the order of one-step prediction can not be improved unless the drift F is linear.

Proposition 3.4 (One-step prediction of LT splitting) Let Assumptions (A1) and (A2) hold, let X be the solution to SDE (1) and let Φ

[LT] h be the LT approximation (10). Then, for k = 1, . . . , N , it holds:

∥E[X t k -Φ [LT] h (X t k-1 ) | X t k-1 = x]∥ = R(h 2 , X t k-1 ).
L p convergence of the LT splitting scheme is established in Theorem 2 in [START_REF] Buckwar | A splitting method for SDEs with locally Lipschitz drift: Illustration on the FitzHugh-Nagumo model[END_REF], which we repeat here for convenience.

Theorem 3.5 (L p convergence of the LT splitting) Let Assumptions (A1) and (A2) hold, let X [LT] be the LT approximation defined in (10), and let X be the solution of (1). Then, there exists C ≥ 1 such that for all p ≥ 2, and k = 1, . . . , N , it holds:

(E[∥X t k -X [LT] t k ∥ p ]) 1 p = R(h, x 0 ).
Now, we investigate the same properties for the S splitting.

Strang splitting

The following proposition states that the S splitting ( 11) has higher order one-step predictions than the LT splitting (10).

The proof can be found in Supplementary Material [START_REF] Pilipovic | Supplement to "Parameter Estimation in Nonlinear Multivariate Stochastic Differential Equations Based on Splitting Schemes[END_REF].

Proposition 3.6 Let Assumptions (A1) and (A2) hold, let X be the solution to (1), and let

Φ [S]
h be the S splitting approximation (11). Then, for k = 1, . . . , N , it holds:

∥E[X t k -Φ [S] h (X t k-1 ) | X t k-1 = x]∥ = R(h 3 , X t k-1 ). (21)
Remark Even though LT and S have the same order of L p convergence, the crucial difference is in the one-step prediction. The approximated transition density between two consecutive data points depends on the one-step approximation. Thus, the objective function based on pseudo-likelihood from the S splitting is more precise than the one from the LT.

To prove L p convergence of the S splitting scheme for (1) with one-sided Lipschitz drift, we follow the same procedure as in [START_REF] Buckwar | A splitting method for SDEs with locally Lipschitz drift: Illustration on the FitzHugh-Nagumo model[END_REF]. The proof of the following theorem is in Section 7.1.

Theorem 3.7 (L p convergence of S splitting) Let Assumptions (A1), (A2) and (A6) hold, let X [S] be the S splitting defined in (11), and let X be the solution of (1). Then, there exists C ≥ 1 such that for all p ≥ 2 and k = 1, . . . , N , it holds:

(E[∥X t k -X [S] t k ∥ p ]) 1 p = R(h, x 0 ).
Before we move to parameter estimation, we prove a useful corollary.

Corollary 3.8 Let all assumptions from Theorem 3.7 hold. Then,

(E[∥Z t k -ξ h,k ∥ p ]) 1/p = R(h, x 0 ).
Proof From the definition of Z t k in (13), it is enough to prove that:

(E[∥f -1 h/2 (X t k ) -µ h (f h/2 (X t k-1 )) -ξ h,k ∥ p ]) 1/p = R(h, x 0 ). From (11) we have that ξ h,k = f -1 h/2 (X [S] t k ) -µ h (f h/2 (X [S] t k-1 )). Then, E[∥f -1 h/2 (X t k ) -µ h (f h/2 (X t k-1 )) -ξ h,k ∥ p ] 1/p ≤ C(E[∥f -1 h/2 (X t k ) -f -1 h/2 (X [S] t k )∥ p ] + E[∥f h/2 (X t k-1 ) -f h/2 (X [S] t k-1 )∥ p ]) 1/p ≤ C(E[∥X t k -X [S] t k ∥ p ] + E[∥X t k-1 -X [S] t k-1 ∥ p ]) 1/p + R(h, x 0
). We used Proposition 2.2, together with the fact that X and X [S] have finite moments and f h/2 and f -1 h/2 grow polynomially. The result follows from L p convergence of the S splitting scheme, Theorem 3.7.

Auxiliary properties

This paper centers around proving the properties of the S estimator. There are two reasons for this. First, most numerical properties in the literature are proved only for LT splitting because proofs for S splitting are more involved. Here, we establish both the numerical properties of the S splitting as well as the properties of the estimator. Second, the S splitting introduces a new pseudo-likelihood that differs from the standard Gaussian pseudo-likelihoods. Consequently, standard tools, like those proposed by [START_REF] Kessler | Estimation of an Ergodic Diffusion from Discrete Observations[END_REF], do not directly apply.

The asymptotic properties of the LT estimator are the same as for the S estimator. However, the following auxiliary properties will be stated and proved only for the S estimator. They can be reformulated for the LT estimator following the same logic.

Before presenting the central results for the estimator, we establish the groundwork with two essential lemmas that rely on the model assumptions. Lemma 4.1 (Lemma 6 in [START_REF] Kessler | Estimation of an Ergodic Diffusion from Discrete Observations[END_REF]) deals with the p-th moments of the solution of the SDE increments and also provides a moment bound of a polynomial map of the solution. The proof of this lemma, presented in Section 7.2, differs from that in [START_REF] Kessler | Estimation of an Ergodic Diffusion from Discrete Observations[END_REF] due to our relaxation of the global Lipschitz assumption of the drift F. Instead, we use a one-sided Lipschitz condition for the drift function F in conjunction with the generalized Grönwall's inequality (Lemma 2.3 in [START_REF] Tian | Nonlinear integral inequality with power and its application in delay integro-differential equations[END_REF] stated in Supplementary Material [START_REF] Pilipovic | Supplement to "Parameter Estimation in Nonlinear Multivariate Stochastic Differential Equations Based on Splitting Schemes[END_REF]) to establish the result. Lemma 4.2 (Lemma 8 in [START_REF] Kessler | Estimation of an Ergodic Diffusion from Discrete Observations[END_REF], Lemma 2 in Sørensen and Uchida ( 2003)) constitutes a central ergodic property that is essential for establishing the asymptotic behavior of the estimator. The proof when the drift F is one-sided Lipschitz is identical to the one presented in [START_REF] Kessler | Estimation of an Ergodic Diffusion from Discrete Observations[END_REF], particularly when combined with Lemma 4.1.

Lemma 4.1 Let Assumptions (A1) and (A2) hold. Let X be the solution of (1). For t k ≥ t ≥ t k-1 , where h = t k -t k-1 < 1, the following two statements hold.

(1) For p ≥ 1, there exists C p > 0 that depends on p, such that:

E[∥X t -X t k-1 ∥ p | F t k-1 ] ≤ C p (t -t k-1 ) p/2 (1 + ∥X t k-1 ∥) Cp .
(2) If g : R d × Θ → R is of polynomial growth in x uniformly in θ, then there exist constants C and C t-t k-1 that depends on t -t t k-1 , such that:

E[|g(X t ; θ)| | F t k-1 ] ≤ C t-t k-1 (1 + ∥X t k-1 ∥) C .
Lemma 4.2 Let Assumptions (A1), (A2) and (A3) hold, and let X be the solution to (1). Let g : R d × Θ → R be a differentiable function with respect to x and θ with derivative of polynomial growth in x, uniformly in θ.

If h → 0 and N h → ∞, then, 1 N N k=1 g (X t k , θ) P θ 0 -----→ N h→∞ h→0 g (x, θ) dν 0 (x),
uniformly in θ.

Lastly, we state the moment bounds needed for the estimator asymptotics. The proof is in Supplementary Material [START_REF] Pilipovic | Supplement to "Parameter Estimation in Nonlinear Multivariate Stochastic Differential Equations Based on Splitting Schemes[END_REF].

Proposition 4.3 (Moment Bounds) Let Assumptions (A1), (A2) and (A6) hold. Let X be the solution of (1), and Z t k as defined in (13). Let g(x; β) be a generic function with derivatives of polynomial growth, and β ∈ Θ β . Then, for k = 1, . . . , N , the following moment bounds hold:

(i) E θ0 [Z t k (β 0 ) | X t k-1 = x] = R(h 3 , X t k-1 ) (ii) E θ0 [Z t k (β 0 )g(X t k ; β) ⊤ | X t k = x] = h 2 (ΣΣ ⊤ 0 D ⊤ g(x; β) + Dg(x; β)ΣΣ ⊤ 0 ) + R(h 2 , X t k-1 ); (iii) E θ0 [Z t k (β 0 )Z t k (β 0 ) ⊤ | X t k-1 = x] = hΣΣ ⊤ 0 + R(h 2 , X t k-1 ).

Asymptotics

The estimators θN are defined in (15). However, we do not need the full objective functions ( 12) and ( 14) to prove consistency and asymptotic normality. It is enough to approximate the covariance matrix Ω h up to the second order by hΣΣ ⊤ + h 2 2 (AΣΣ ⊤ + ΣΣ ⊤ A ⊤ ) (see equation ( 6)). Indeed, after applying Taylor series on the inverse of Ω h , we get:

Ω h (θ) -1 = 1 h (ΣΣ ⊤ ) -1 (I + h 2 (A(β) + ΣΣ ⊤ A(β) ⊤ (ΣΣ ⊤ ) -1 ) -1 ) + R(h, x 0 ) = 1 h (ΣΣ ⊤ ) -1 (I - h 2 (A(β) + ΣΣ ⊤ A(β) ⊤ (ΣΣ ⊤ ) -1 ) + R(h, x 0 ) = 1 h (ΣΣ ⊤ ) -1 - 1 2 ((ΣΣ ⊤ ) -1 A(β) + A(β) ⊤ (ΣΣ ⊤ ) -1 ) + R(h, x 0 ).
Similarly, we approximate the log-determinant as:

log det Ω h (θ)= log det(hΣΣ ⊤ + h 2 2 (A(β)ΣΣ ⊤ + ΣΣ ⊤ A(β) ⊤ )) + R(h 2 , x 0 ) θ = log det ΣΣ ⊤ + log det(I + h 2 (A(β) + ΣΣ ⊤ A(β) ⊤ (ΣΣ ⊤ ) -1 )) + R(h 2 , x 0 ) = log det ΣΣ ⊤ + h 2 Tr(A(β) + ΣΣ ⊤ A(β) ⊤ (ΣΣ ⊤ ) -1 ) + R(h 2 , x 0 ) = log det ΣΣ ⊤ + h Tr A(β) + R(h 2 , x 0 ).
Using the same approximation we obtain:

2 log | det Df h/2 (x; β) |= 2 log | det(I + h 2 DN(x; β))| = 2 log |1 + h 2 Tr DN(x; β)| + R(h, x) = h Tr DN(x; β) + R(h 2 , x 0 )
Subsequently, retaining terms up to order R(N h 2 , x 0 ) from objective functions ( 12) and ( 14), we establish the approximate objective functions:

L [LT] N (θ) := N log det ΣΣ ⊤ +N h Tr A(β) + 1 h N k=1 (X t k -µ h (f h (X t k-1 ; β); β)) ⊤ (ΣΣ ⊤ ) -1 (X t k -µ h (f h (X t k-1 ; β); β)) (22) - N k=1 (X t k -µ h (f h (X t k-1 ; β); β)) ⊤ (ΣΣ ⊤ ) -1 A(β)(X t k -µ h (f h (X t k-1 ; β); β)) L [S] N (θ) := N log det ΣΣ ⊤ +N h Tr A(β) + 1 h N k=1 Z t k (β) ⊤ (ΣΣ ⊤ ) -1 Z t k (β) - N k=1 Z t k (β) ⊤ (ΣΣ ⊤ ) -1 A(β)Z t k (β) + h N k=1 Tr DN(X t k ; β). ( 23
)
Unlike other likelihood-based methods, such as [START_REF] Kessler | Estimation of an Ergodic Diffusion from Discrete Observations[END_REF], [START_REF] Aït-Sahalia | Maximum Likelihood Estimation of Discretely Sampled Diffusions: A Closed-form Approximation Approach[END_REF][START_REF] Aït-Sahalia | Closed-form likelihood expansions for multivariate diffusions[END_REF], [START_REF] Choi | Closed-form likelihood expansions for multivariate time-inhomogeneous diffusions[END_REF][START_REF] Choi | Explicit form of approximate transition probability density functions of diffusion processes[END_REF], [START_REF] Yang | A new delta expansion for multivariate diffusions via the Itô-Taylor expansion[END_REF], our estimators do not involve expansions. The objective functions are formulated in simple terms without hyperparameters, such as the order of the expansions. Hence, our approach is robust and user-friendly, as we directly employ ( 12) and ( 14) without requiring approximations. However, we leverage the approximations ( 22) and ( 23) for the mathematical analysis and the proofs.

Consistency

Now, we state the consistency of βN and ΣΣ

P θ 0 --→ β 0 , ΣΣ ⊤ N P θ 0 --→ ΣΣ ⊤ 0 .

Asymptotic normality

In this Section, we state the asymptotic normality of the estimator. First, we need some preliminaries. Let ρ > 0 and B ρ (θ 0 ) = {θ ∈ Θ | ∥θ -θ 0 ∥ ≤ ρ} be a ball around θ 0 . Since θ 0 ∈ Θ, for sufficiently small ρ > 0, B ρ (θ 0 ) ∈ Θ. Let L N be one of the two objective functions ( 22) or ( 23). For θN ∈ B ρ (θ 0 ), the mean value theorem yields:

1 0 H L N (θ 0 + t( θN -θ 0 )) dt ( θN -θ 0 ) = -∇L N (θ 0 ) . ( 24
) With ς := vech(ΣΣ ⊤ ) = ([ΣΣ ⊤ ] 11 , [ΣΣ ⊤ ] 12 , [ΣΣ ⊤ ] 22 , ..., [ΣΣ ⊤ ] 1d , ..., [ΣΣ ⊤ ] dd )
, we half-vectorize ΣΣ ⊤ to avoid working with tensors when computing derivatives with respect to ΣΣ ⊤ . Since ΣΣ ⊤ is a symmetric d × d matrix, ς is of dimension s = d(d + 1)/2. For a diagonal matrix, instead of a half-vectorization, we use ς := diag(ΣΣ ⊤ ). Define:

C N (θ) := 1 N h ∂ ββ L N (θ) 1 N √ h ∂ βς L N (θ) 1 N √ h ∂ βς L N (θ) 1 N ∂ ςς L N (θ) , ( 25
)
s N := √ N h( βN -β 0 ) √ N (ς N -ς 0 ) , λ N :=    - 1 √ N h ∂ β L N (θ 0 ) - 1 √ N ∂ ς L N (θ 0 )    , (26)
and

D N := 1 0 C N (θ 0 + t( θN -θ 0 )) dt. Then, (24) is equivalent to D N s N = λ N . Let: C(θ 0 ) := C β (θ 0 ) 0 r×s 0 s×r C ς (θ 0 ) , (27)
where:

[C β (θ 0 )] i1,i2 := (∂ βi 1 F 0 (x)) ⊤ (ΣΣ ⊤ 0 ) -1 (∂ βi 2 F 0 (x)) dν 0 (x), 1 ≤ i 1 , i 2 ≤ r, [C ς (θ 0 )] j1,j2 := 1 2 Tr((∂ς j1 ΣΣ ⊤ 0 )(ΣΣ ⊤ 0 ) -1 (∂ς j2 ΣΣ ⊤ 0 )(ΣΣ ⊤ 0 ) -1 ), 1 ≤ j 1 , j 2 ≤ s.
Now, we state the theorem for asymptotic normality, whose proof is in Section 7.4.

Theorem 5.2 Let Assumptions (A1)-(A6) hold, X be the solution of (1), and θ N = ( β N , ς N ) be the estimator that minimizes one of the objective functions (22) or (23

). If θ 0 ∈ Θ, C(θ 0 ) is positive definite, h → 0, N h → ∞, and N h 2 → 0, then, √ N h( βN -β 0 ) √ N (ς N -ς 0 ) d -→ N (0, C -1 (θ 0 )), (28)
under P θ0 .

The estimator of the diffusion parameter converges faster than the estimator of the drift parameter. [START_REF] Gobet | LAN property for ergodic diffusions with discrete observations[END_REF] showed that for a discretely sampled SDE model, the optimal convergence rates for the drift and diffusion parameters are 1/ √ N h and 1/ √ N , respectively. Thus, our estimators reach optimal rates. Moreover, the estimators are asymptotically efficient since C is the Fisher information matrix for the corresponding continuous-time diffusion (see [START_REF] Kessler | Estimation of an Ergodic Diffusion from Discrete Observations[END_REF], [START_REF] Gobet | LAN property for ergodic diffusions with discrete observations[END_REF]). Finally, since the asymptotic correlation is zero between the drift and diffusion estimators , they are asymptotically independent.

6 Simulation study

Estimators used in the study

The EM transition distribution (16) for the Lorenz system (20) is:

X t k Y t k Z t k | X t k-1 Y t k-1 Z t k-1 = x y z ∼ N   x + hp(y -x) y + h(rx -y -xz) z + h(xy -cz) ,   hσ 2 1 0 0 0 hσ 2 2 0 0 0 hσ 2 3     .
We do not write the closed-form distributions for the K (17) and LL (18) estimators, but we use the corresponding formulas to implement the likelihoods. We implement the two splitting strategies proposed in Section 2.5, leading to four estimators: LT mix , LT avg , S mix , and S avg . To further speed up computation time, we use the same trick for calculating Ω h in (6) as the one suggested by [START_REF] Gu | Parameter Estimation for Multivariate Nonlinear Stochastic Differential Equation Models: A Comparison Study[END_REF] for calculating Ω [LL] h . Namely, for the splitting schemes, we adapt P 3 from (19) accordingly:

P 3 = A ΣΣ ⊤ 0 d×d -A ⊤ .
Therefore, computing exp(hP 3 ) circumvents the need for evaluating the integral in Ω h (6), following the approach described in Section 2.4.4.

Trajectory simulation

To simulate sample paths, we use the EM discretization with a step size of h sim = 0.0001, which is small enough for the EM discretization to perform well. Then, we sub-sample the trajectory to get a larger time step h, decreasing discretization errors. We perform M = 1000 Monte Carlo repetitions.

Optimization in R

To optimize the objective functions we use the R package torch [START_REF] Falbel | torch: Tensors and Neural Networks with 'GPU' Acceleration[END_REF], which uses AD instead of the traditional finite differentiation used in optim. The two main advantages of AD are precision and speed. Finite differentiation is subject to floating point precision errors and is slow in high dimensions [START_REF] Baydin | Automatic Differentiation in Machine Learning: A Survey[END_REF]. Conversely, AD is exact and fast and thus used in numerous applications, such as MLE or training neural networks.

We tried all available optimizers in the torch package and chose the resilient backpropagation algorithm optim_rprop based on [START_REF] Riedmiller | RPROP -A Fast Adaptive Learning Algorithm[END_REF]. It performed faster than the rest and was more precise in finding the global minimum. We used the default hyperparameters and set the optimization iterations to 200. We chose the precision of 10 -5 between the updated and the parameters from the previous iteration as the convergence criteria. For starting values, we used (0.1, 0.1, 0.1, 0.1, 0.1, 0.1). All estimators converged after approximately 80 iterations.

Comparing criteria

We compare seven estimators based on their precision and speed. For the precision, we compute the absolute relative error (ARE) for each component θ(i) N of the estimator θN :

ARE( θ(i) N ) = 1 M M r=1 | θ(i) N,r -θ (i) 0,r | θ (i) 0,r .
For S and LL, we compare the distributions of θNθ 0 to investigate the precision more closely.

The running times are calculated using the tictoc package in R, measured from the start of the optimization step until the convergence criterion is met. To avoid the influence of running time outliers, we compute the median over M repetitions.

Results

In Figure 2, AREs are shown as a function of the discretization step h. For a clearer comparison, we use the log scale on the y axis. While most estimators work well for a step size no greater than 0.01, only LL, S mix , and S avg perform well for h = 0.05. The LT avg is not competitive even for h = 0.005. The performance of LT mix varies, sometimes approaching the performance of K, while other times performing similarly to EM. Thus, LT mix is not a good choice for this specific model. The bias of EM starts to show for h = 0.01 escalating for h = 0.05. The largest bias appears in the diffusion parameters, which is due to the poor approximation of Ω EM h . K is less biased than EM except for p and r when h = 0.05. Note how some parameters are estimated better for larger h when N is fixed. This is due to a longer observation interval T = N h, reflecting the √ N h rate of convergence.

Since S mix , S avg , and LL perform the best with large time steps, we zoom in on their distributions in Figure 3. To make the figure clearer, we removed some outliers for σ 2 1 and σ 2 2 . This did not change the shape of the distributions, it only truncated the tails. The three estimators perform similarly, especially for small h. For h = 0.05, the drift parameters are underestimated by approximately 5 -10%, while the diffusion parameters are overestimated by up to 20%. Both S estimators exhibited superior performance compared to LL, except for the parameters p and σ 2 1 . While the LL and S estimators perform similarly in terms of precision, Figure 4 shows the superiority of the S estimators over LL in computational costs. The LL becomes increasingly computationally expensive for increasing N because it calculates N covariance matrices for each parameter value. The second slowest estimator is S mix , followed by LT mix , S avg , K, LT avg , and, finally, EM is the fastest. The speed of EM is almost constant in N. Additionally, it seems that the running times do not depend on h. Thus, we recommend using the S estimators, especially for large N .

Figures 5 and6 show that the theoretical results hold for the S mix and LT mix estimators. We compare how the distributions of θNθ 0 change with sample size N and step size h. With increasing N , the variance decreases, whereas the mean does not change. For that, we need smaller h. To obtain negligible bias for LT mix , we need a step size smaller than h = 0.005 . However, S mix is practically unbiased up to h = 0.01. This shows that LT estimators might not be a good choice in practice, while S estimators are.

The solid black lines in Figures 5 and6 represent the theoretical asymptotic distributions for each parameter computed from (28). For the Lorenz system (20), the precision matrix ( 27) is given by:

C(θ 0 ) = diag (y -x) 2 σ 2 1,0 dν 0 (x), x 2 σ 2 2,0 dν 0 (x), z 2 σ 2 3,0 dν 0 (x), 1 2σ 4 1,0 , 1 2σ 4 2,0 , 1 2σ 4 3,0 .
The integrals are approximated by taking the mean over all data points and all Monte Carlo repetitions.

Some outliers of σ2

2 are removed from Figures 5 and 6 by truncating the tails. The proof of Proposition 3.6 is in Supplementary Material [START_REF] Pilipovic | Supplement to "Parameter Estimation in Nonlinear Multivariate Stochastic Differential Equations Based on Splitting Schemes[END_REF]. Here, we present the proof of L p convergence stated in Theorem 3.7. 2). To prove condition (1), we need to prove the following property:

(E[∥X t k -Φ [S] h (X t k-1 )∥ p | X t k-1 = x]) 1 p = R(h q2 , x),
where q 2 = 3/2. We start with

∥X t k -Φ [S] h (X t k-1 )∥ p = ∥X t k -X t k-1 -hF(X t k-1 ) -ξ h,k + R(h 3/2 , X t k-1 )∥ p .

For more details on the expansion of Φ [S]

h , see Supplementary Material [START_REF] Pilipovic | Supplement to "Parameter Estimation in Nonlinear Multivariate Stochastic Differential Equations Based on Splitting Schemes[END_REF]. We approximate

ξ h,k = t k t k-1 e A(t k -s) Σ dW s by: ξ h,k = t k t k-1 (I + (t k -s)A)Σ dW s + R(h 2 , X t k-1 ) = Σ(W t k -W t k-1 ) + AΣ t k t k-1 (t k -s) dW s + R(h 2 , X t k-1 ).
Using the fact that

t k t k-1 (t k -s) dW s ∼ N (0, h 3 3 I), we deduce that ξ h,k = Σ(W t k -W t k-1 ) + R(h 3/2 , X t k-1
). Then, Hölder's inequality yields:

∥X t k -X t k-1 -hF(X t k-1 ) -Σ(W t k -W t k-1 )∥ p ≤ h p-1 t k t k-1 ∥(F(X s ) -F(X t k-1 ))∥ p ds.
Assumption (A2), the integral norm inequality, Cauchy-Schwartz, and Hölder's inequalities, together with the mean value theorem yield:

E[∥X t k -Φ [S] h (X t k-1 )∥ p | X t k-1 = x] ≤ C(E[h p-1 t k t k-1 ∥F(X s ) -F(X t k-1 )∥ p ds | X t k-1 = x]) ≤C(h p-1 t k t k-1 E[∥X s -X t k-1 ∥ p ∥ 1 0 D x F(X s -u(X s -X t k-1 )) du∥ p | X t k-1 = x] ds) ≤ C h p-1 t k t k-1 (E[∥X s -X t k-1 ∥ 2p | X t k-1 = x]) 1 2 (E[∥ 1 0 D x F(X s -u(X s -X t k-1 )) du∥ 2p | X t k-1 = x]) 1 2 ds ≤ C(h p-1 t k t k-1 h p 2 ds) = R(h 3p/2 , x).
In the last line, we used Lemma 4.1. This proves condition (1) of Theorem 3.3. Now, we prove condition (2). We use ( 5) and ( 11) to write X

[S]

t k = f h/2 (e Ah (f h/2 (X [S] t k-1 ) -X [1] t k-1 ) + X [1] t k ). Define R t k := e Ah (f h/2 (X [S] t k ) -X [1]
t k ), and use the associativity (9) to get

R t k = e Ah (f h (R t k-1 + X [1] t k ) -X [1] t k ).
The proof of the boundness of the moments of R t k is the same as in Lemma 2 in [START_REF] Buckwar | A splitting method for SDEs with locally Lipschitz drift: Illustration on the FitzHugh-Nagumo model[END_REF]. Finally, we have X

[S]

t k = f -1 h/2 (e -Ah R t k + X [1]
t k ). Since f -1 h/2 grows polynomially and X

[1]

t k has finite moments, X

[S]

t k must have finite moments too. This concludes the proof.

Proof of Lemma 4.1

Proof of Lemma 4.1 We first prove (1). In the following, C 1 and C 2 denote constants. We use the triangular inequality and Hölder's inequality to obtain:

∥X t -X t k-1 ∥ p ≤ 2 p-1 (∥ t t k-1 F(X s ; θ) ds∥ p + ∥Σ(W t -W t k-1 )∥ p ) ≤ 2 p-1 (( t t k-1 C 1 (1 + ∥X s ∥) C1 ds) p + ∥Σ(W t -W t k-1 )∥ p ) ≤ 2 p-1 C p 1 ( t t k-1 (1 + ∥X s -X t k-1 ∥ + ∥X t k-1 ∥) C1 ds) p + 2 p-1 ∥Σ(W t -W t k-1 )∥ p ≤ 2 C1+2p-3 C p 1 (t -t k-1 ) p-1 ( t t k-1 ∥X s -X t k-1 ∥ pC1 ds + (t -t k-1 ) p (1 + ∥X t k-1 ∥) pC1 ) + 2 p-1 ∥Σ(W t -W t k-1 )∥ p .
In the second inequality, we used the polynomial growth (A2) of F. Furthermore, for some constant C 2 that depends on p, we have

E[∥Σ(W t -W t k-1 )∥ p | F t k-1 ] = (t -t t k-1 ) p/2 C 2 (p).
Then, for h < 1, there exists a constant C p that depends on p, such that:

C p (t -t k-1 ) 2p-1 (1 + ∥X t k-1 ∥) Cp + C p (t -t t k-1 ) p/2 ≤ C p (t -t k-1 ) p/2 (1 + ∥X t k-1 ∥) Cp .
The last inequality holds because the term of order p/2 is dominating when

t -t k-1 < 1. Denote m(t) = E[∥X t - X t k-1 ∥ p | F t k-1 ].
Then, we have:

m(t) ≤ C p (t -t k-1 ) p/2 (1 + ∥X t k-1 ∥) Cp + C p t t k-1 m C1 (s) ds. ( 29
)
Now, we apply the generalized Grönwall's inequality (Lemma 2.3 in [START_REF] Tian | Nonlinear integral inequality with power and its application in delay integro-differential equations[END_REF], stated in Supplementary Material [START_REF] Pilipovic | Supplement to "Parameter Estimation in Nonlinear Multivariate Stochastic Differential Equations Based on Splitting Schemes[END_REF]) on ( 29). Since we consider a super-linear growth, we can assume that there exist C 1 > 1 and C p > 0, such that:

m(t) ≤ C p (t -t k-1 ) p/2 (1 + ∥X t k-1 ∥) Cp + (κ 1-C1 (t) -(C 1 -1)2 C1-1 C p (t -t k-1 )) 1 1-C 1 ≤ C p (t -t k-1 ) p/2 (1 + ∥X t k-1 ∥) Cp + Cκ(t), (30)
where

κ(t) = C p (t -t k-1 ) C1p/2+1 (1 + ∥X t k-1 ∥) Cp .
The bound C in inequality (30) makes sense, because the term:

(1 -(C 1 -1)2 C1-1 C p (t -t k-1 )κ 1 1-C 1 (t)) 1 1-C 1
is positive by Lemma 2.3 from [START_REF] Tian | Nonlinear integral inequality with power and its application in delay integro-differential equations[END_REF]. Additionally, the same term reaches its maximum value of 1, for t = t k-1 . The constant C in (30) includes some terms that depend on t -t k-1 . However, these terms will not change the dominating term of κ(t) since h < 1. Finally, the terms in κ(t) are of order p/2, thus for large enough constant

C p , it holds m(t) ≤ C p (t -t k-1 ) p/2 (1 + ∥X t k-1 ∥) Cp .
To prove (2), we use that g is of polynomial growth:

E[|g(X t ; θ)| | F t k-1 ] ≤ C 1 E[(1 + ∥X t k-1 ∥ + ∥X t -X t k-1 ∥) C1 | F t k-1 ] ≤ C 2 (1 + ∥X t k-1 ∥ C1 + E[∥X t -X t k-1 ∥ C1 | F t k-1 ]).
Now, we apply the first part of the lemma, to get:

E[|g(X t ; θ)| | F t k-1 ] ≤ C 2 (1 + ∥X t k-1 ∥ C1 + C ′ t-t k-1 (1 + ∥X t k-1 ∥) C3 ) ≤ C t-t k-1 (1 + ∥X t k-1 ∥) C . That concludes the proof.

Proof of consistency of the estimator

The proof of consistency consists in studying the convergence of the objective function that defines the estimators. The objective function L N (β, ς) (23) can be decomposed into sums of martingale triangular arrays. We thus first state a lemma that proves the convergence of each triangular array involved in the objective function. Then, we will focus on the proof of consistency. The proof of the Lemma is in Supplementary Material [START_REF] Pilipovic | Supplement to "Parameter Estimation in Nonlinear Multivariate Stochastic Differential Equations Based on Splitting Schemes[END_REF].

Lemma 7.1 Let Assumptions (A1)-(A6) hold, and X be the solution of (1). Let g, g 1 , g 2 : R d × Θ × Θ → R be differentiable functions with respect to x and θ, with derivatives of polynomial growth in x, uniformly in θ. If h → 0 and N h → ∞, then:

1. 1 N h N k=1 Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 Z t k (β 0 ) P θ 0 -----→ N h→∞ h→0 Tr((ΣΣ ⊤ ) -1 ΣΣ ⊤ 0 ); 2. h N N k=1 g(X t k-1 ; β 0 , β) ⊤ (ΣΣ ⊤ ) -1 g(X t k-1 ; β 0 , β) P θ 0 -----→ N h→∞ h→0 0; 3. 1 N N k=1 Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 g(X t k-1 ; β 0 , β) P θ 0 -----→ N h→∞ h→0 0; 4. 1 N h N k=1 Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 g(X t k-1 ; β 0 , β) P θ 0 -----→ N h→∞ h→0 0; 5. 1 N N k=1 Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 g(X t k ; β 0 , β) P θ 0 -----→ N h→∞ h→0 0; 6. 1 N h N k=1 Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 g(X t k ; β 0 , β) P θ 0 -----→ N h→∞ h→0
Tr(Dg(x; β 0 , β)ΣΣ ⊤ 0 (ΣΣ ⊤ ) -1 ) dν 0 (x);

7. h N N k=1 g 1 (X t k-1 ; β 0 , β) ⊤ (ΣΣ ⊤ ) -1 g 2 (X t k ; β 0 , β) P θ 0 -----→ N h→∞ h→0 0,
uniformly in θ.

Proof of Theorem 5.1 To establish consistency, we follow the proof of Theorem 1 in [START_REF] Kessler | Estimation of an Ergodic Diffusion from Discrete Observations[END_REF] and study the limit of L

[S]

N (β, ς) from ( 23), rescaled by the correct rate of convergence. More precisely, the consistency of the diffusion parameter is proved by studying the limit of 1 N L

[S]

N (β, ς), while the consistency of the drift parameter is proved by studying the limit of 1 N h (L

[S] N (β, ς) -L [S]
N (β 0 , ς)). We start with the consistency of the diffusion parameter ς. We need to prove that:

1 N L [S] N (β, ς) → log(det(ΣΣ ⊤ )) + Tr((ΣΣ ⊤ ) -1 ΣΣ ⊤ 0 ) =: G 1 (ς, ς 0 ), (31)
in P θ0 , for N h → ∞, h → 0, uniformly in θ. To study the limit, we first decompose 1 N L

[S]

N (β, ς) as follows:

1 N L [S] N (β, ς) = log det ΣΣ ⊤ + T 1 + T 2 + T 3 + 2(T 4 + T 5 + T 6 ) + R(h, x 0 ). ( 32
)
The terms T 1 , . . . , T 6 are derived from the quadratic form in ( 23) by adding and subtracting the corresponding terms with β 0 , followed by rearrangements, resulting in the following expressions:

T 1 := 1 N h N k=1 Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 Z t k (β 0 ), T 2 := 1 N h N k=1 (f -1 h/2,k (β) -f -1 h/2,k (β 0 )) ⊤ (ΣΣ ⊤ ) -1 (f -1 h/2,k (β) -f -1 h/2,k (β 0 )), T 3 := 1 N h N k=1 (µ h,k-1 (β 0 ) -µ h,k-1 (β)) ⊤ (ΣΣ ⊤ ) -1 (µ h,k-1 (β 0 ) -µ h,k-1 (β)), T 4 := 1 N h N k=1 Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 (µ h,k-1 (β 0 ) -µ h,k-1 (β)), T 5 := 1 N h N k=1 (f -1 h/2,k (β) -f -1 h/2,k (β 0 )) ⊤ (ΣΣ ⊤ ) -1 (µ h,k-1 (β 0 ) -µ h,k-1 (β)), T 6 := 1 N h N k=1 (f -1 h/2,k (β) -f -1 h/2,k (β 0 )) ⊤ (ΣΣ ⊤ ) -1 Z t k (β 0 ). Previously, we defined f -1 h/2,k (β) := f -1 h/2 (X t k ; β) and µ h,k-1 (β) := µ h (f h/2 (X t k-1 ; β); β)
. These terms will also play a significant role in proving the asymptotic normality.

The first term of ( 32) is a constant. Properties 1, 2, 3, 5, and 7 from Lemma 7.1 give the following limits T 1 → Tr((ΣΣ ⊤ ) -1 ΣΣ ⊤ 0 ) and for l = 2, 3, ..., 6, T l → 0, uniformly in θ. The convergence in probability is equivalent to the existence of a subsequence converging almost surely. Thus, the convergence in ( 31) is almost sure for a subsequence ( β N l , ς N l ). This implies:

ς N l P θ 0 -a.s. ------→ N h→∞ h→0 ς ∞ . The compactness of Θ implies that ( β N l , ς N l ) converges to a limit (β ∞ , ς ∞ ) almost surely. By continuity of the mapping ς → G 1 (ς, ς 0 ) we have 1 N l L [S] N l ( βN l , ς N l ) → G 1 (ς ⊤ ∞ , ς 0 ), in P θ0 , for N h → ∞, h → 0, uniformly in θ. By the definition of the estimator, G 1 (ς ∞ , ς 0 ) ≤ G 1 (ς 0 , ς 0 ). We also have: G 1 (ς ∞ , ς 0 ) ≥ G 1 (ς 0 , ς 0 ) ⇔ log(det(ΣΣ ⊤ ∞ )) + Tr((ΣΣ ⊤ ∞ ) -1 ΣΣ ⊤ 0 ) ≥ log(det(ΣΣ ⊤ 0 )) + Tr(I d ) ⇔ Tr((ΣΣ ⊤ ∞ ) -1 ΣΣ ⊤ 0 ) -log(det((ΣΣ ⊤ ∞ ) -1 ΣΣ ⊤ 0 )) ≥ d ⇔ d i=1 λ i -log d i=1 λ i ≥ d i=1 1 ⇔ d i=1 (λ i -1 -log λ i ) ≥ 0,
where λ i are the eigenvalues of (ΣΣ ⊤ ∞ ) -1 ΣΣ ⊤ 0 , which is a positive definite matrix. The last inequality follows since for any positive x, log x ≤ x -1. Thus, G 1 (ς ∞ , ς 0 ) = G 1 (ς 0 .ς 0 ). Then, all the eigenvalues λ i must be equal to 1, hence, ΣΣ ⊤ ∞ = ΣΣ ⊤ 0 . We proved that a convergent subsequence of ς N tends to ς 0 almost surely. From there, the consistency of the estimator of the diffusion coefficient follows.

We now focus on the consistency of the drift parameter. The objective is to prove that the following limit in P θ0 , for N h → ∞, h → 0, uniformly with respect to θ:

1 N h (L [S] N (β, ς) -L [S] N (β 0 , ς)) → G 2 (β 0 , ς 0 , β, ς), (33)
where:

G 2 (β 0 , ς 0 , β, ς) := (F 0 (x) -F(x)) ⊤ (ΣΣ ⊤ ) -1 (F 0 (x) -F(x)) dν 0 (x) + Tr(D(N 0 (x) -N(x))(ΣΣ ⊤ 0 (ΣΣ ⊤ ) -1 -I)) dν 0 (x).
To prove it, we decompose 1 N h (L

[S] N (β, ς) -L [S]
N (β 0 , ς)) as follows:

1 N h (L [S] N (β, ς) -L [S] N (β 0 , ς)) = Tr(A(β) -A(β 0 )) + 1 h (T 2 + T 3 + 2(T 4 + T 5 + T 6 )) + 1 N h N k=1 (Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 A(β 0 )Z t k (β 0 ) -Z t k (β) ⊤ (ΣΣ ⊤ ) -1 A(β)Z t k (β)) (34)
+ 1 N N k=1 Tr D(N(X t k ; β) -N(X t k ; β 0 )) + R(h, x 0 ). The term 1 N h N k=1 (Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 A(β 0 )Z t k (β 0 ) -Z t k (β) ⊤ (ΣΣ ⊤ ) -1 A(β)Z t k (β)) converges to Tr(A(β 0 ) -A(β))
, which thus cancels out with the first term in (34). Lemma 4.2 provides the uniform convergence of 1 h T 2 with respect to θ:

1 h T 2 = 1 4N N k=1 (N 0 (X t k ) -N(X t k )) ⊤ (ΣΣ ⊤ ) -1 (N 0 (X t k ) -N(X t k )) + R(h, x 0 ) → 1 4 (N 0 (x) -N(x)) ⊤ (ΣΣ ⊤ ) -1 (N 0 (x) -N(x)) dν 0 (x).
The limit of 1 h T 3 computes analogously. To prove 1 h T 4 → 0, we use Lemma 9 in Genon-Catalot and Jacod (1993) and Property 4 from Lemma 7.1. Lemma 4.2 yields:

1 h T 5 P θ 0 -----→ N h→∞ h→0 1 4 (N 0 (x) -N(x)) ⊤ (ΣΣ ⊤ ) -1 (N 0 (x) -N(x)) dν 0 (x) + 1 2 (A 0 (x-b 0 ) -A(x-b)) ⊤ (ΣΣ ⊤ ) -1 (N 0 (x) -N(x)) dν 0 (x).
Finally,

1 h T 6 → 1 2 Tr(D(N 0 (x) -N(x)) ⊤ ΣΣ ⊤ 0 (ΣΣ ⊤ ) -1
) dν 0 (x) uniformly in θ, by Property 6 of Lemma 7.1. Lemma 4.2 gives:

1 N N k=1 Tr D(N(X t k ) -N 0 (X t k )) P θ 0 -----→ N h→∞ h→0 Tr D(N(x) -N 0 (x)) dν 0 (x),
uniformly in θ. This proves (33). Then, there exists a subsequence N l such that ( β N l , ς N l ) converges to a limit (β ∞ , ς ∞ ), almost surely. By continuity of the mapping (β, ς) → G 2 (β 0 , ς 0 , β, ς), for N l h → ∞, h → 0, we have the following convergence in P θ0 :

1

N l h (L [S] N l ( β N l , ς N l ) -L [S] N l (β 0 , ς N l )) → G 2 (β 0 , ς 0 , β ∞ , ς ∞ ). Then, G 2 (β 0 , ς 0 , β ∞ , ς ∞ ) ≥ 0 since ΣΣ ⊤ ∞ = ΣΣ ⊤ 0 .
On the other hand, by the definition of the estimator

L [S] N l ( β N l , ς N l ) -L [S]
N l (β 0 , ς N l ) ≤ 0. Thus, the identifiability assumption (A5) concludes the proof for the S estimator.

To prove the same statement for the LT estimator, the representation of the objective function (32) has to be adapted. In the LT case, this representation is straightforward. There is no extra logarithmic term and only three instead of six auxiliary T terms are used. This is due to the Gaussian transition density in the LT approximation.

Proof of asymptotic normality of the estimator

Proof of Theorem 5.2 According to Theorem 1 in [START_REF] Kessler | Estimation of an Ergodic Diffusion from Discrete Observations[END_REF] or Theorem 1 in [START_REF] Sørensen | Small-diffusion asymptotics for discretely sampled stochastic differential equations[END_REF], Lemmas 7.2 and 7.3 below are enough for establishing the asymptotic normality of θN . Here, we only present the outline of the proof. For more details, see proof of Theorem 1 in [START_REF] Sørensen | Small-diffusion asymptotics for discretely sampled stochastic differential equations[END_REF].

Lemma 7.2 Let C N (θ 0 ) and C(θ 0 ) be as defined in (25) and (27), respectively. If h → 0, N h → ∞, and ρ N → 0, then:

C N (θ 0 ) P θ 0 --→ 2C(θ 0 ), sup ∥θ∥≤ρ N ∥C N (θ 0 + θ) -C N (θ 0 )∥ P θ 0 --→ 0. Lemma 7.3 Let λ N be as defined (26). If h → 0, N h → ∞ and N h 2 → 0, then: λ N d -→ N (0, 4C(θ 0 )), under P θ0 .
Lemma 7.2 states that C N (θ 0 ) approaches 2C(θ 0 ) as h → 0 and N h → ∞. Moreover, the difference between C N (θ 0 + θ) and C N (θ 0 ) approaches zero when θ approaches θ 0 , within a distance specified by balls B ρ N (θ 0 ), where ρ N → 0. To ensure the asymptotic normality of θN , Lemma 7.2 is employed to restrict the term ∥D N -C N (θ 0 )∥ when θN ∈ Θ ∩ B ρ N (θ 0 ) as follows:

∥D N -C N (θ 0 )∥⊮ { θN ∈Θ∩Bρ N (θ0)} ⩽ sup θ∈Bρ N (θ0) ∥C N (θ) -C N (θ 0 )∥ P θ 0 -----→ N h→∞ h→0 0
Applying again Lemma 7.2 on the previous line, we get D N → 2C(θ 0 ) in P θ0 , as h → 0 and N h → ∞.

Lemma 7.3 establishes the convergence in distribution of λ N to N (0, 4C(θ 0 )), under P θ0 , as h → 0 and N h → ∞. This result provides the groundwork for the asymptotic normality of θN . Indeed, consider the set D N composed of instances where D N is invertible. The probability, under θ 0 , of D N occurring approaches 1, as h → 0 and N h → ∞. This implies that D N is almost surely invertible in this limit. Furthermore, we define E N as the intersection of { θN ∈ Θ} and D N . Then, it can be shown that

⊮ E N → 1 in P θ0 when h → 0 and N h → ∞. For E N := D N on E N , we have E N → 2C(θ 0 ) in P θ0 as h → 0 and N h → ∞. Given that s N ⊮ E N = E -1 N D N s N ⊮ E N = E -1 N λ N ⊮ E N and according to Lemma 7.3, s N ⊮ E N → N (0, C(θ 0 ) -1 ) in distribution as h → 0, N h → ∞ and N h 2 → 0.
In conclusion, under P θ0 , as h → 0, N h → ∞ and N h 2 → 0, s N ⊮ E N is shown to converge in distribution to N (0, C(θ 0 ) -1 ). The asymptotic normality for θN is, thus, confirmed due to the convergence of ⊮ E N → 1.

Proof of Lemma 7.2 To prove the first part of the lemma, we aim to represent C N (θ 0 ) from the objective function ( 14). In doing so, we again employ the approximation (23), focusing solely on the terms that do not converge to zero as N h → ∞ and h → 0. We start as in the approximation (34) and compute the corresponding derivatives to obtain the first block matrix of C N (25). We begin with

∂ βi 1 βi 2 L [S] N (β, ς): 1 N h ∂ βi 1 βi 2 L [S] N (β, ς)= ∂ βi 1 βi 2 Tr A(β) + 1 N N k=1 ∂ βi 1 βi 2 Tr DN(X t k ; β) + ∂ βi 1 βi 2 1 h T 2 (β 0 , β, ς) + T 3 (β 0 , β, ς) + 2(T 4 (β 0 , β, ς) + T 5 (β 0 , β, ς) + T 6 (β 0 , β, ς)) - 1 N h N k=1 ∂ βi 1 βi 2 (Z t k (β) ⊤ (ΣΣ ⊤ ) -1 A(β)Z t k (β))) + R(h, x 0 ).
To determine the convergence of each of the previous terms, we use the definitions of the sums T i s and approximate each T i using Proposition 2.2 and the Taylor expansion of the function µ h . As we apply the derivatives ∂β i1 β i2 , the order of h in each sum increases since terms of order R(1, x 0 ) are constant with respect to β. Finally, when evaluating

1 N h ∂ βi 1 βi 2 L [S]
N (β, ς) at θ = θ 0 , numerous terms will cancel out due to differences of the type g(β 0 ; X t k , X t k-1 ) -g(β; X t k , X t k-1 ). Using the results from Lemma 7.1 and the proof of Theorem 5.1, we get the following limits:

∂ βi 1 βi 2 1 h T 2 (β 0 , β, ς 0 ) β=β0 P θ 0 --→ 1 2 (∂ βi 1 N 0 (x)) ⊤ (ΣΣ ⊤ 0 ) -1 ∂ βi 2 N 0 (x) dν 0 (x), ∂ βi 1 βi 2 1 h T 3 (β 0 , β, ς 0 ) β=β0 P θ 0 --→ 1 2 (∂ βi 1 N 0 (x) + 2∂ βi 1 A 0 (x-b 0 )) ⊤ (ΣΣ ⊤ 0 ) -1 (∂ βi 2 N 0 (x) + 2∂ βi 2 A 0 (x-b 0 )) dν 0 (x), ∂ βi 1 βi 2 1 h T 5 (β 0 , β, ς 0 ) β=β0 P θ 0 --→ 1 2 (∂ βi 1 F 0 (x)) ⊤ (ΣΣ ⊤ 0 ) -1 ∂ βi 2 N 0 (x) dν 0 (x) + 1 2 (∂ βi 2 A 0 (x-b 0 )) ⊤ (ΣΣ ⊤ 0 ) -1 ∂ βi 1 N 0 (x) dν 0 (x), ∂ βi 1 βi 2 1 h T 6 (β 0 , β, ς 0 ) β=β0 P θ 0 --→ - 1 2 Tr(D∂ βi 1 βi 2 N 0 (x)) dν 0 (x),
for N h → ∞, h → 0. Since 1 h T 4 → 0, the partial derivatives go to zero too. From Lemma 4.2, for N h → ∞, h → 0, we have:

1 N N k=1 ∂ βi 1 βi 2 Tr DN(X t k ; β) P θ 0 --→ Tr(D∂ βi 1 βi 2 N 0 (x)) dν 0 (x). Term 1 N h N k=1 ∂ βi 1 βi 2 (Z t k (β) ⊤ (ΣΣ ⊤ ) -1 A(β)Z t k (β))
, evaluated in θ = θ 0 , has only one term of order h:

1 N h N k=1 Z t k (β 0 ) ⊤ (ΣΣ ⊤ 0 ) -1 ∂ βi 1 βi 2 A(β 0 )Z t k (β 0 ), which converges to ∂ βi 1 βi 2 Tr A(β 0 ) (Property 1 Lemma 7.1). Thus, 1 N h ∂ βi 1 βi 2 L [S] N (β, ς 0 )| β=β0 → 2 (∂ βi 2 F 0 (x)) ⊤ (ΣΣ ⊤ 0 ) -1 ∂ βi 2 F 0 (x) dν 0 (x), in P θ0 for N h → ∞, h → 0. Now, we prove 1 N √ h ∂ βς L [S] N (β, ς)| β=β0,ς=ς0 → 0, in P θ0 for N h → ∞, h → 0.
For a constant C h , depending on h, l = 2, 3, ..., 6, and generic functions g, g 1 , the following term is at most of order R(h, x 0 ):

∂ βi T l (β, ς) = C h N k=1 (g(β 0 ; X t k , X t k-1 ) -g(β; X t k , X t k-1 )) ⊤ (ΣΣ ⊤ ) -1 g 1 (β; X t k , X t k-1 ), Then, term ∂ βς L [S]
N (β, ς) still contains g(β 0 ; X t k , X t k-1 ) -g(β; X t k , X t k-1 ) which is 0 for β = β 0 . Moreover, the term 1

N N k=1 ∂ βς (Z t k (β) ⊤ (ΣΣ ⊤ ) -1 A(β)Z t k (β)) is at most of order R(h, x 0 ).
Thus,

1 N √ h ∂ βς L [S] N (β, ς)| β=β0,ς=ς0 = 0. Finally, we compute 1 N ∂ ςj 1 ςj 2 L [S]
N (β, ς). As before, it holds 1 N ∂ ςj 1 ςj 2 T l (β, ς)| β=β0,ς=ς0 → 0, for l = 2, 3, ..., 6. Similarly, we see that

1 N N k=1 Z t k (β 0 ) ⊤ ∂ ςj 1 ςj 2 (ΣΣ ⊤ ) -1 A(β 0 )Z t k (β 0
) is at most of order R(h, x 0 ). So, we need to compute the following second derivatives ∂ ςj 1 ςj 2 log(det ΣΣ ⊤ ) and

∂ ςj 1 ςj 2 1 N h N k=1 Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 Z t k (β 0 ).
The first one yields:

∂ ςj 1 ςj 2 log(det ΣΣ ⊤ ) = Tr((ΣΣ ⊤ ) -1 ∂ ςj 1 ςj 2 ΣΣ ⊤ ) -Tr((ΣΣ ⊤ ) -1 (∂ ςj 1 ΣΣ ⊤ )(ΣΣ ⊤ ) -1 ∂ ςj 2 ΣΣ ⊤ ).
On the other hand, we have:

∂ ςj 1 ςj 2 1 N h N k=1 Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 Z t k (β 0 ) = - 1 N h N k=1 Tr(Z t k (β 0 )Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 (∂ ςj 1 ςj 2 ΣΣ ⊤ )(ΣΣ ⊤ ) -1 ) + 1 N h N k=1 Tr(Z t k (β 0 )Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 (∂ ςj 1 ΣΣ ⊤ )(ΣΣ ⊤ ) -1 (∂ ςj 2 ΣΣ ⊤ )(ΣΣ ⊤ ) -1 ) + 1 N h N k=1 Tr(Z t k (β 0 )Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 (∂ ςj 2 ΣΣ ⊤ )(ΣΣ ⊤ ) -1 (∂ ςj 1 ΣΣ ⊤ )(ΣΣ ⊤ ) -1 ).
Then, from Property 1 of Lemma 7.1, we get:

∂ ςj 1 ςj 2 1 N h N k=1 Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 Z t k (β 0 ) ς=ς0 P θ 0 -----→ N h→∞ h→0 2 Tr((ΣΣ ⊤ 0 ) -1 (∂ ςj 1 ΣΣ ⊤ 0 )(ΣΣ ⊤ 0 ) -1 ∂ ςj 2 ΣΣ ⊤ 0 ) -Tr((ΣΣ ⊤ 0 ) -1 ∂ ςj 1 ςj 2 ΣΣ ⊤ 0 ). Thus, 1 N ∂ ςj 1 ςj 2 L [S] N (β, ς)| β=β0,ς=ς0 → Tr((ΣΣ ⊤ 0 ) -1 (∂ ςj 1 ΣΣ ⊤ 0 )(ΣΣ ⊤ 0 ) -1 ∂ ςj 2 ΣΣ ⊤ 0 )
. Since all the limits used in this proof are uniform in θ, the first part of the lemma is proved. The second part is trivial, because all limits are continuous in θ.

Proof of Lemma 7.3 First, we compute the first derivatives. We start with:

∂ βi L [S] N (β, ς) = -2 N k=1 Tr(Df h/2,k (β)D x ∂ βi f -1 h/2,k (β)) + 2 h N k=1 (f -1 h/2,k (β) -µ h,k-1 (β)) ⊤ (ΣΣ ⊤ ) -1 (∂ βi f -1 h/2,k (β) -∂ βi µ h,k-1 (β)).
The first derivative with respect to ς is:

∂ ςj L [S] N (β, ς) = N ∂ ςj log det(ΣΣ ⊤ ) + 1 h ∂ ςj N k=1 (f -1 h/2,k (β) -µ h,k-1 (β)) ⊤ (ΣΣ ⊤ ) -1 (f -1 h/2,k (β) -µ h,k-1 (β)) = - 1 h N k=1 Tr (f -1 h/2,k (β) -µ h,k-1 (β))(f -1 h/2,k (β) -µ h,k-1 (β)) ⊤ (ΣΣ ⊤ ) -1 (∂ ςj ΣΣ ⊤ )(ΣΣ ⊤ ) -1 + Tr((ΣΣ ⊤ ) -1 ∂ ςj ΣΣ ⊤ ) .
Define:

η (i) N,k (θ) := 2 √ N h Tr(Df h/2,k (β)D x ∂ βi f -1 h/2,k (β)) (35) - 2 √ N hh Z t k (β) ⊤ (ΣΣ ⊤ ) -1 ∂ βi (f -1 h/2,k (β) -µ h,k-1 (β)) ζ (j) N,k (θ) := 1 √ N h Tr(Z t k (β)Z t k (β) ⊤ (ΣΣ ⊤ ) -1 (∂ς j ΣΣ ⊤ )(ΣΣ ⊤ ) -1 ) (36) - 1 √ N Tr((ΣΣ ⊤ ) -1 ∂ς j ΣΣ ⊤ ),
and rewrite

λ N as λ N = N k=1 [η (1) N,k (θ 0 ), . . . , η (r) N,k (θ 0 ), ζ (1)
N,k (θ 0 ), . . . , ζ (s)
N,k (θ 0 )] ⊤ . Now, by Proposition 3.1 from [START_REF] Crimaldi | Convergence results for multivariate martingales[END_REF], it is sufficient to prove Lemma 7.4. For more details, see Supplementary Material [START_REF] Pilipovic | Supplement to "Parameter Estimation in Nonlinear Multivariate Stochastic Differential Equations Based on Splitting Schemes[END_REF].

Lemma 7.4 Let η (i) N,k (θ) and ζ (j)
N,k (θ) be defined as in (35) and (36), respectively. If h → 0, N h → ∞, and N h 2 → 0, then for and all i, i 1 , i 2 = 1, 2, ..., r, and j, j 1 , j 2 = 1, 2, ..., s, it holds:

[(i)] 1. E θ0 [sup 1≤k≤N |η (i) N,k (θ 0 )|] -→ 0, and E θ0 [sup 1≤k≤N |ζ (j) N,k (θ 0 )|] -→ 0; 2. N k=1 E θ0 [η (i) N,k (θ 0 ) | X t k-1 ] P θ 0 --→ 0, and N k=1 E θ0 [ζ (j) N,k (θ 0 ) | X t k-1 ] P θ 0 --→ 0; 3. N k=1 E θ0 [η (i1) N,k (θ 0 ) | X t k-1 ]E θ0 [η (i2) N,k (θ 0 ) | X t k-1 ] P θ 0 --→ 0; 4. N k=1 E θ0 [ζ (j1) N,k (θ 0 ) | X t k-1 ]E θ0 [ζ (j2) N,k (θ 0 ) | X t k-1 ] P θ 0 --→ 0; 5. N k=1 E θ0 [η (i) N,k (θ 0 ) | X t k-1 ]E θ0 [ζ (j) N,k (θ 0 ) | X t k-1 ] P θ 0 --→ 0; 6. N k=1 E θ0 [η (i1) N,k (θ 0 )η (i2) N,k (θ 0 ) | X t k-1 ] P θ 0 --→ 4[C β (θ 0 )] i1i2 ; 7. N k=1 E θ0 [ζ (j1) N,k (θ 0 )ζ (j2) N,k (θ 0 ) | X t k-1 ] P θ 0 --→ 4[C ς (θ 0 )] j1j2 ; 8. N k=1 E θ0 [η (i) N,k (θ 0 )ζ (j) N,k (θ 0 ) | X t k-1 ] P θ 0 --→ 0; 9. N k=1 E θ0 [(η (i1) N,k (θ 0 )η (i2) N,k (θ 0 )) 2 | X t k-1 ] P θ 0 --→ 0; 10. N k=1 E θ0 [(ζ (j1) N,k (θ 0 )ζ (j2) N,k (θ 0 )) 2 | X t k-1 ] P θ 0 --→ 0; 11. N k=1 E θ0 [(η (i) N,k (θ 0 )ζ (j) N,k (θ 0 ) 2 ) | X t k-1 ] P θ 0 --→ 0.
The proof of the previous Lemma is technical and is shown in Supplementary Material [START_REF] Pilipovic | Supplement to "Parameter Estimation in Nonlinear Multivariate Stochastic Differential Equations Based on Splitting Schemes[END_REF].

Conclusion

We proposed two new estimators for nonlinear multivariate SDEs. These estimators are based on splitting schemes, a numerical approximation that preserves all important properties of the model. It was known that the LT splitting scheme has L p convergence rate of order 1. We proved that the same holds for the S splitting. This result was expected because the overall trajectories of the S and LT splittings coincide up to the first h/2 and the last h/2 move of the flow Φ [2] h/2 . Nonetheless, S splitting is more precise in one-step predictions, which is crucial for the estimators because the objective function consists of densities between consecutive data points. Therefore, the obtained S estimator is less biased than the LT.

We proved that both estimators have optimal convergence rates for discretized observations of the SDEs. These rates are

√

N for the diffusion parameter and √ N h for the drift parameter. We also showed that the asymptotic variance of the estimators is the inverse of the Fisher information for the continuous time model. Thus, the estimators are efficient.

In the simulation study conducted with the stochastic Lorenz system, we the superior performance of the S estimators. We compared seven estimators based on different discretization schemes. Estimators based on Ozaki's LL and the S splitting schemes demonstrated the highest precision. However, the running time of LL is notably influenced by the sample size N , unlike the S estimator, which experiences a more gradual increase in runtime with larger N . This property makes the S estimator a more appropriate choice for large sample sizes. The LT, EM, and K estimators perform well for small h, but for larger h the bias increases, especially for the diffusion parameters in the EM case.

While the proposed estimators are versatile, they come with certain limitations. These include assumptions like the presence of additive noise and equidistant observations. However, under specific conditions, the Lamperti transformation can relax the constraint of additive noise. Equidistant observations can easily be relaxed due to the continuous-time formulation. Furthermore, we assumed that the diffusion parameter ΣΣ ⊤ is invertible. However, there are applications where models with degenerate noise naturally arise, like second-order differential equations. We will thoroughly investigate these cases in another paper with more involved proofs.

S1 Supplementary Material

This section provides proofs for all propositions, lemmas, and theorems. References to equations and sections that do not begin with "S" refer to the main paper. If not stated, we assume the parameters are the true ones and the expectations are taken under the probability measure. Occasionally, we omit explicit parameter notation to enhance clarity. For instance, E implicitly denotes E θ .

In Section S1.1, we provide the proof for the Lie-Trotter splitting, while Section S1.2 contains the proofs for the Strang splitting. Additionally, the proofs of moment bounds are detailed in Section S1.3. These three sections rely on manipulating the Itô-Taylor expansion. The foundational properties necessary for subsequent proofs are outlined in Section S1.4. These properties encompass Grönwall's and Rosenthal's inequalities, as well as Central Limit Theorems for a sum of triangular arrays. The proof for Lemma 7.1 can be found in Section S1.5, while the section addressing asymptotic normality is presented in Section S1.6.

S1.1 Proof for the Lie-Trotter splitting

Proof of Proposition 3.4 To establish the proposition, we compare the actual first moment of the solution to SDE (1), as obtained from Lemma 2.1, with the moment derived through Taylor expansion of the LT approximation. First, we prove the proposition for Lie-Trotter splitting as defined in the paper. By performing the Taylor expansion of

E[Φ [LT] h (x)] = µ h (f h/2 ((x)) =e Ah f h (x)+(I -e Ah )
b around h = 0, using Proposition 2.2, we arrive at:

µ h (f h/2 ((x)) = x + h(A(x-b) + N(x)) + h 2 2 (A 2 (x-b) + 2AN(x) + (DN(x))N(x)) + R(h 3 , x). (S1)
The coefficient of h in (S1) is F(x), which aligns with the coefficient of h in the theoretical moment of the solution to (1) as provided in Lemma 2.1. However, in Lemma 2.1, Σ appears in the coefficient of h 2 , while it does not appear in (S1). Consequently, to achieve the order of convergence R(h 3 , x), we need to make the following unrealistic assumption.

(SA)

d i=1 d j=1 [ΣΣ ⊤ ] ij ∂ 2 ij F (i) (x) = 0, for all k = 1, . . . , d.
Upon comparing expression (S1) with the true moments of the SDE solution under Assumption (SA), we arrive at (DF(x))N(x) = (DN(x))F(x) to ensure equality of the coefficient at order h 2 . However, the last equation holds true for all x ∈ R d only when N is linear. Therefore, achieving the order R(h 3 , x) one-step convergence is feasible only if SDE (1) is linear.

We now aim to show that changing the composition order within the LT splitting does not affect the one-step convergence order. To demonstrate this, we define the reversed Lie-Trotter splitting:

X [LT]⋆ t k := Φ [LT]⋆ h (X [LT]⋆ t k-1 ) = (Φ [2] h • Φ [1] h )(X [LT]⋆ t k-1 ) = f h (µ h (X [LT]⋆ t k-1 ) + ξ h,k ). We compute E[f h (µ h (X t k-1 ) + ξ h,k ) | X t k-1 = x], which is equivalent to calculating E[f h (X [1] t k ) | X [1] t k-1 = x] = E[f h (µ h (X [1] t k-1 ) + ξ h,k ) | X [1] t k-1 = x]. The infinitesimal generator L [1] for SDE (3) is defined on the class of sufficiently smooth functions g : R d → R by L [1] g(x) = (A(x-b)) ⊤ ∂g(x) ∂x + 1 2 Tr(ΣΣ ⊤ H g (x)
). This yields:

E[g(X [1] t k ) | X [1] t k-1 = x] = g(x) + hL [1] g(x) + h 2 2 L 2 [1] g(x) + R(h 3 , x). (S2) We apply (S2) to g(x) = f (i) h (x). For calculating L [1] f (i) h (x) and L 2 [1] f (i) h (x)
, we use the Taylor expansion of f h (x) around h = 0, as provided in Proposition 2.2. The partial derivatives are ∂ j f

(i) h (x) = δ i j + h∂ j N (i) (x) + R(h 2 , x) and ∂ 2 jk f (i) h (x) = h∂ 2 jk N (i) (x) + R(h 2 , x). Since L [1] f (i)
h (x) is multiplied by h in (S2), we only need to calculate it up to order R(h, x). We have L

[1] f (i) h (x) = (A(x-b)) (i) + h(A(x-b)) ⊤ ∇N (i) (x) + h 2 Tr(ΣΣ ⊤ H N (i) (x)) + R(h 2 , x). Similarly, we have L 2 [1] f (i) h (x) = (A(x-b)) ⊤ ∇(A(x-b)) (i) + R(h, x) = (A(x-b)) ⊤ A (i) + R(h, x). Thus, E[f (i) h (X [1] t k-1 ) | X [1] t k-1 = x] = x (i) + hN (i) (x) + h 2 2 (N(x)) ⊤ ∇N (i) (x) (S3) + h(A(x-b)) (i) + h 2 (A(x-b)) ⊤ ∇N (i) (x) + h 2 2 Tr(ΣΣ ⊤ H N (i) (x)) + h 2 2 (A(x-b)) ⊤ A (i) + R(h 3 , x) = x (i) + hF (i) (x) + h 2 2 ((F(x)) ⊤ (∇N (i) (x)) + (A(x-b)) ⊤ ∇F (i) (x) + Tr(ΣΣ ⊤ H N (i) (x))) + R(h 3 , x).
Lastly, we compute:

E[ξ ⊤ h H N (i) (X)ξ h | X = x] = E[tr(ξ ⊤ h H N (i) (X)ξ h ) | X = x] = tr(H N (i) (X)E[ξ h ξ ⊤ h ]) = d j,k=1 ∂ 2 jk N (i) (x)[var(ξ h )] jk = d j,k=1 ∂ 2 jk F (i) (x)[Ω h ] jk .
We use the approximation of the variance of the random vector

ξ h to get E[N(X + Q h (X) + ξ h ) | X = x] = N(x) + (DN(x))Q h (x) + h 2 [ d j,k=1 [ΣΣ ⊤ ] jk ∂ 2 jk F (i) (x)] d i=1 + R(h 2 , x).
Taking the expectation of (S4) and incorporating the previous equation completes the proof.

S1.3 Proofs of the Moment Bounds

Before proving the moment bounds, we first demonstrate in Lemma S1.2 how the infinitesimal generator L operates on a product of two functions. Lemma S1.2 Let L be the infinitesimal generator defined in the main text of SDE (1). For sufficiently smooth functions α, β : R d → R, it holds:

L(α(x)β(x)) = α(x)Lβ(x) + β(x)Lα(x) + 1 2 Tr(ΣΣ ⊤ (∇α(x)∇ ⊤ β(x) + ∇β(x)∇ ⊤ α(x))).
Proof We use the generator L and the product rule to get:

L(α(x)β(x)) = F(x) ⊤ α(x)∇β(x) + F(x) ⊤ β(x)∇α(x) + 1 2 Tr(ΣΣ ⊤ (α(x)H β (x) + β(x)H α (x))) + 1 2 Tr(ΣΣ ⊤ (∇α(x)∇ ⊤ β(x) + ∇β(x)∇ ⊤ α(x))) = α(x)Lβ(x) + β(x)Lα(x) + 1 2 Tr(ΣΣ ⊤ (∇α(x)∇ ⊤ β(x) + ∇β(x)∇ ⊤ α(x))).
This concludes the proof.

Proof of Proposition 4.3 Proof of (i). Lemma S1.1 yields:

E[f -1 h/2 (X t k ) -µ h (f h/2 (X t k-1 )) | X t k-1 = x] = E[f -1 h/2 (X t k ) | X t k-1 = x] -µ h (f h/2 (x)) = E[f -1 h/2 (X t k ) | X t k-1 = x] -f -1 h/2 (x) -hF(x) - h 2 2 AF(x) + R(h 3 , x).
Now, we use the infinitesimal generator L to evaluate the expectation in the last line where the generator L is applied to a vector-valued function. We have:

E[f -1 h/2 (X t k ) | X t k-1 = x] = f -1 h/2 (x) + hLf -1 h/2 (x) + h 2 2 L 2 f -1 h/2 (x) + R(h 3 , x).
We use f -1 h/2 (x) = f -h/2 (x) and Proposition 2.2 to get:

Lf -1 h/2 (x) = Lx - h 2 LN(x) + R(h 2 , x) = F(x) - h 2 LN(x) + R(h 2 , x), L 2 f -1 h/2 (x) = LA(x-b) + LN(x) + R(h, x) = AF(x) + LN(x) + R(h, x). It follows that E[f -1 h/2 (X t k ) -µ h (f h/2 (X t k-1 )) | X t k-1 = x] = R(h 3 , x).
Proof of (ii). In this proof, we distinguish the true parameters θ 0 from a generic parameter θ. We start with the expansions of f -1 h and µ h :

E θ0 [(f -1 h/2 (X t k ; β 0 ) -µ h (f h/2 (X t k-1 ; β 0 ); β 0 ))g(X t k ; β) ⊤ | X t k-1 = x] = E θ0 [X t k g(X t k ; β) ⊤ | X t k-1 = x] - h 2 E θ0 [N(X t k ; β 0 )g(X t k ; β) ⊤ | X t k-1 = x] -xE θ0 [g(X t k ; β) ⊤ | X t k-1 = x] - h 2 (2A 0 (x-b 0 ) + N 0 (x))E θ0 [g(X t k ; β) ⊤ | X t k-1 = x] + R(h 2 , x) = xg(x; β) ⊤ + hL θ0 (xg(x; β) ⊤ ) - h 2 N 0 (x)g(x; β) ⊤ -xg(x; β) ⊤ -hxL θ0 g(x; β) ⊤ -hA 0 (x-b 0 )g(x; β) ⊤ - h 2 N 0 (x)g(x; β) ⊤ + R(h 2 , x) = hL θ0 (xg(x; β) ⊤ ) -hxL θ0 g(x; β) ⊤ -hF 0 (x)g(x; β) ⊤ + R(h 2 , x).
Lastly, Lemma S1.2 and the definition of L θ0 yield:

L θ0 (xg(x; β) ⊤ ) = xL θ0 g(x; β) ⊤ + (L θ0 x)g(x; β) ⊤ + 1 2 (ΣΣ ⊤ 0 D ⊤ g(x; β) + Dg(x; β)ΣΣ ⊤ 0 ) = xL θ0 g(x; β) ⊤ + F(x; β 0 )g(x; β) ⊤ + 1 2 (ΣΣ ⊤ 0 D ⊤ g(x; β) + Dg(x; β)ΣΣ ⊤ 0 ).
Proof of (iii). We introduce g(X t k ; β 0 ) = f -1 h/2 (X t k ; β 0 ) and use (ii) to show:

E θ0 [(f -1 h/2 (X t k ; β 0 ) -µ h (f h/2 (X t k-1 ; β 0 ); β 0 ))(f -1 h/2 (X t k ; β 0 ) -µ h (f h/2 (X t k-1 ; β 0 ); β 0 )) ⊤ | X t k-1 = x] = h 2 (ΣΣ ⊤ 0 D ⊤ g(x; β 0 ) + Dg(x; β 0 )ΣΣ ⊤ 0 ) -E θ0 [f -1 h/2 (X t k ; β 0 ) -µ h (f h/2 (X t k-1 ; β 0 ); β 0 ) | X t k-1 = x]µ h (f h/2 (x; β 0 ); β 0 ) ⊤ + R(h 2 , x
). The result follows from property (i) and Dg(x; β 0 ) = I + R(h, x).

S1.4 Auxiliary properties

In this section, we revisit crucial properties essential for establishing the consistency and asymptotic normality of the proposed estimators. To begin, we invoke Lemma 2.3 from Tian and Fan (2020) as Lemma S1.3, which was used in proving Lemma 4.1. This lemma offers a generalization of the Grönwall's inequality. Furthermore, Lemma 9 in Genon-Catalot and Jacod (1993) provides conditions for the convergence of a sum of a triangular array and is recalled as Lemma S1.4. Lemmas S1.5 and S1.6 give sufficient conditions for uniform convergence. The former is sourced from Proposition A1 in [START_REF] Gloter | Parameter Estimation for a Discretely Observed Integrated Diffusion Process[END_REF], while the latter comes from Lemma 3.1 from [START_REF] Yoshida | Asymptotic behavior of M-estimator and related random field for diffusion process[END_REF]. On occasions, Lemma S1.5 might not suffice, warranting the use of Lemma S1.6. Theorem S1.7 is a helpful tool for assessing the conditions of these two lemmas is the Rosenthal's inequality for martingales (Theorem 2.12 in [START_REF] Hall | Martingale Limit Theory and Its Application[END_REF]).

Lastly, Theorem S1.8 presents a special case of the central limit theorem for multivariate martingale triangular arrays (Proposition 3.1 from [START_REF] Crimaldi | Convergence results for multivariate martingales[END_REF]). This theorem is pivotal for proving the asymptotic normality of the proposed estimators. (S7)

Proof of Lemma 7.1 Proof of 1. As previously discussed, we introduce a martingale array that corresponds to the limit outlined in point 1. We then utilize Lemma S1.4 to facilitate our analysis. We denote Y N k (β 0 , ς) := 1 N h Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 Z t k (β 0 ). We have:

N k=1 E θ0 [Y N k (β 0 , ς) | X t k-1 ] = 1 N h N k=1 E θ0 [Tr(Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 Z t k (β 0 )) | X t k-1 ] = 1 N h N k=1 Tr((ΣΣ ⊤ ) -1 E θ0 [Z t k (β 0 )Z t k (β 0 ) ⊤ | X t k-1 ]) = 1 N h N k=1
Tr((ΣΣ ⊤ ) -1 hΣΣ ⊤ 0 + R(h 2 , X t k-1 ))

P θ 0 -----→ N h→∞ h→0
Tr((ΣΣ ⊤ ) -1 ΣΣ ⊤ 0 ).

To use the result of Lemma S1.4, we need to prove that covariance of Y N k (β 0 , ς) goes to zero. To achieve this, we leverage Corollary 3.8 and recall that if ρ is a Gaussian random vector ρ ∼ N (0, Π), then E[(ρ T Mρ) 2 ] = 2 Tr((MΠ) 2 ) + (Tr(MΠ)) 2 . This leads to:

N k=1 E θ0 [Y N k (β 0 , ς) 2 | X t k-1 ] = 1 N 2 h 2 N k=1 (E θ0 [(ξ ⊤ h,k (ΣΣ ⊤ ) -1 ξ h,k ) 2 | X t k-1 ] + R(h 3/2 , X t k-1 )) = 1 N h 1 N N k=1
(2 Tr((ΣΣ ⊤ ) -1 Σ 0 Σ ⊤ 0 ) 2 + (Tr((ΣΣ ⊤ ) -1 Σ 0 Σ ⊤ 0 )) 2 + R(h 1/2 , X t k-1 ))

P θ 0 --→ 0,
for N h → ∞, h → 0. Then, by Lemma S1.4 1

N h N k=1 Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 Z t k (β 0 ) P θ 0
--→ Tr((ΣΣ ⊤ ) -1 ΣΣ ⊤ 0 ), for N h → ∞, h → 0. To establish the uniformity of the limits with respect to ς, we turn to Lemma S1.5 and introduce sets Θ ςj such that ς = (ς 1 , ς 2 , . . . , ς s ) ∈ Θ ς1 × Θ ς2 × • • • × Θ ςs = Θ ς . Then it is enough to show that for all j = 1, . . . , s, it holds:

sup N ∈N E θ0 [ sup ςj ∈Θς j |∂ ςj 1 N h N k=1 Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 Z t k (β 0 )|] < ∞. (S9)
We use the well-known rule of matrix differentiation ∂ X (a ⊤ X -1 a) = -X -1 aa ⊤ X -1 , where a is a vector and X is a symmetric matrix, to get:

∂ x (i)
Tr(a ⊤ C -1 (x)a) = -Tr(C -1 (x)aa ⊤ C -1 (x)∂ x (i) C(x)) = -Tr(aa ⊤ C -1 (x)(∂ x (i) C(x))C -1 (x)).

We omit writing β 0 for ease of notation. Then, by using the trace bound, the norm inequality, and Assumption (A4), we can deduce that:

sup N ∈N E θ0 [ sup ςj ∈Θς j |∂ ςj 1 N h N k=1 Z ⊤ t k (ΣΣ ⊤ ) -1 Z t k |] ≤ sup N ∈N E θ0 [ 1 N h N k=1 sup ςj ∈Θς j |∂ ςj Tr(Z ⊤ t k (ΣΣ ⊤ ) -1 Z t k )|] ≤ sup N ∈N E θ0 [ 1 N h N k=1 Tr(Z t k Z ⊤ t k ) sup ςj ∈Θς j ∥(ΣΣ ⊤ ) -1 (∂ ςj ΣΣ ⊤ )(ΣΣ ⊤ ) -1 ∥] ≤ sup N ∈N E θ0 [ 1 N h N k=1 Tr(Z t k Z ⊤ t k ) sup ςj ∈Θς j ∥(ΣΣ ⊤ ) -1 ∥ 2 ∥∂ ςj ΣΣ ⊤ ∥] ≤ C sup N ∈N E θ0 [ 1 N h N k=1 Tr(Z t k Z ⊤ t k )] = C sup N ∈N E θ0 [E θ0 [ 1 N h N k=1 Tr(Z t k Z ⊤ t k ) | X t k-1 ]] = C sup N ∈N E θ0 [ 1 N h N k=1
Tr(hΣΣ ⊤ 0 + R(h 2 , X t k-1 ))] < ∞.

Proof of 2. We use Lemma 4.2 to deduce:

1 N N k=1
g(X t k-1 ; β 0 , β) ⊤ (ΣΣ ⊤ ) -1 g(X t k-1 ; β 0 , β)

P θ 0
--→ g(x; β 0 , β) ⊤ (ΣΣ ⊤ ) -1 g(x; β 0 , β) dν 0 (x), uniformly in θ, for N h → ∞, h → 0. Then we use the bound of g to conclude the proof of 2.

Proof of 3. For Y N k (β 0 , θ) := 1 N Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 g(X t k-1 ; β 0 , β), the limit of

N k=1 E θ0 [Y N k (β 0 , θ) | X t k-1 ] rewrites as: N k=1 E θ0 [Y N k (β 0 , θ) | X t k-1 ] = 1 N N k=1 E θ0 [Tr(Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 g(X t k-1 ; β 0 , β)) | X t k-1 ] = 1 N N k=1
Tr((ΣΣ ⊤ ) -1 g(X t k-1 ;

β 0 , β)E θ0 [Z t k (β 0 ) ⊤ | X t k-1 ]) = 1 N N k=1
R(h 3 , X t k-1 )

P θ 0 --→ 0,
for N h → ∞, h → 0 . Then, we study the limit of

N k=1 E θ0 [Y N k (β 0 , θ) 2 | X t k-1 ]: N k=1 E θ0 [Y N k (β 0 , θ) 2 | X t k-1 ] = 1 N 2 N k=1 E θ0 [g(X t k-1 ; β 0 , β) ⊤ (ΣΣ ⊤ ) -1 Z t k (β 0 )Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 g(X t k-1 ; β 0 , β) | X t k-1 ] = 1 N 2 N k=1 g(X t k-1 ; β 0 , β) ⊤ (ΣΣ ⊤ ) -1 E θ0 [Z t k (β 0 )Z t k (β 0 ) ⊤ | X t k-1 ](ΣΣ ⊤ ) -1 g(X t k-1 ; β 0 , β) = 1 N N k=1 R( h N , X t k-1 ) P θ 0 --→ 0,
for N h → ∞, h → 0. Lemma S1.4 yields that 1 N N k=1 Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 g(X t k-1 ; β 0 , β)

P θ 0
--→ 0, for N h → ∞, h → 0. To show the uniformity of the limits with respect to θ, we leverage Lemma S1.6. It is sufficient to demonstrate the existence of constants p ≥ l > r + s and C > 0 such that for all θ, θ 1 and θ 2 it holds:

E θ0 [| N k=1 Y N k (β 0 , θ)| p ] ≤ C, (S10) E θ0 [| N k=1 (Y N k (β 0 , θ 1 ) -Y N k (β 0 , θ 2 ))| p ] ≤ C∥θ 1 -θ 2 ∥ l . (S11)
We begin by considering equation (S10). Based on the definition of Z t k (β 0 ) and the assumptions made about N, as well as the fact that h < 1, there exist constants C 1 and C 2 such that:

∥Z t k (β 0 )∥ p ≤ ∥X t k -X t k-1 ∥ p + C 1 h p (1 + ∥X t k ∥) C1 + C 2 h p (1 + ∥X t k-1 ∥) C2 , (S12)
Then, Lemma 4.1 yields:

E θ0 [∥Z t k (β 0 )∥ p | X t k-1 ] ≤ Ch p/2 (1 + ∥X t k-1 ∥) C . (S13)
Subsequently, we use the norm inequality, (S13) and both statements of Lemma 4.1 to get:

E θ0 [| N k=1 Y N k (β 0 , θ)| p ] ≤ N p-1 N k=1 E θ0 [|Y N k (β 0 , θ)| p ] = 1 N N k=1 E θ0 [E θ0 [|Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 g(X t k-1 ; β 0 , β)| p | X t k-1 ]] ≤ 1 N N k=1 E θ0 [E θ0 [∥Z t k (β 0 )∥ p | X t k-1 ]∥(ΣΣ ⊤ ) -1 ∥ p ∥g(X t k-1 ; β 0 , β)∥ p ] ≤ 1 N • N • C. (S14)
This completes the proof of (S10). Now, we focus on (S11). We use the triangular inequality and the Hölder's inequality to derive:

E θ0 [| N k=1 (Y N k (β 0 , θ 1 ) -Y N k (β 0 , θ 2 ))| p ] ≤ 2 p-1 N N k=1 E θ0 [|Z t k (β 0 ) ⊤ (Σ 1 Σ ⊤ 1 ) -1 (g(X t k-1 ; β 1 , β 0 ) -g(X t k-1 ; β 2 , β 0 ))| p ] (S15) + 2 p-1 N N k=1 E θ0 [|Z t k (β 0 ) ⊤ ((Σ 1 Σ ⊤ 1 ) -1 -(Σ 2 Σ ⊤ 2 ) -1 )g(X t k-1 ; β 2 , β 0 )| p ]. (S16)
First, we study sum (S15). We use the mean value theorem and the triangular inequalities to get:

1 N N k=1 E θ0 [|Z t k (β 0 ) ⊤ (Σ 1 Σ ⊤ 1 ) -1 (g(X t k-1 ; β 1 , β 0 ) -g(X t k-1 ; β 2 , β 0 ))| p ] ≤ 1 N N k=1 E θ0 [E θ0 [∥Z t k (β 0 )∥ p | X t k-1 ]∥(Σ 1 Σ ⊤ 1 ) -1 ∥ p ∥g(X t k-1 ; β 1 , β 0 ) -g(X t k-1 ; β 2 , β 0 )∥ p ] ≤ 1 N N k=1 E θ0 [C p (1 + ∥X t k-1 ∥) Cp ∥β 1 -β 2 ∥ p ∥ 1 0 D β g(X t k-1 ; β 2 + t(β 1 -β 2 ), β 0 ) dt∥ p ] ≤ C∥β 1 -β 2 ∥ p . (S17)
To bound sum (S16), we introduce the following multivariate matrix-valued function G(ς) := (ΣΣ ⊤ ) -1 . Then, we use the inequality between the operator 2-norm and Frobenius norm, and the definition of the Frobenius norm to get:

∥G(ς 1 ) -G(ς 2 )∥ ≤ ( d i,j=1 ∥G ij (ς 1 ) -G ij (ς 2 )∥ 2 ) 1 2 .
Now, apply the mean value theorem on each G ij and Assumption (A4) to get:

∥G(ς 1 ) -G(ς 2 )∥ ≤ ( d i,j=1 ∥ς 1 -ς 2 ∥ 2 ∥ t 0 ∇ ς G ij (ς 2 + t(ς 1 -σ 2 )) dt∥ 2 ) 1 2 ≤ C∥ς 1 -ς 2 ∥.
Finally, combining the previous results, we conclude that:

E θ0 [| N k=1 (Y N k (β 0 , θ 1 ) -Y N k (β 0 , θ 2 ))| p ] ≤ C(∥β 1 -β 2 ∥ p + ∥ς 1 -ς 2 ∥ p ) ≤ C(∥β 1 -β 2 ∥ 2 + ∥ς 1 -ς 2 ∥ 2 ) p/2 = C∥θ 1 -θ 2 ∥ p ,
for p ≥ 2. This concludes the proof of 3.

Proof of 4. For Y N k (β 0 , θ) := 1 N h Z t k (β 0 ) ⊤ (ΣΣ ⊤ ) -1 g(X t k-1 ; β 0 , β), we repeat the same derivations as in the proof of 3. to show that the limit of

N k=1 E θ0 [Y N k (β 0 , θ) | X t k-1 ] satisfies: N k=1 E θ0 [Y N k (β 0 , θ) | X t k-1 ] = 1 N h N k=1
Tr((ΣΣ ⊤ ) -1 g(X t k-1 ;

β 0 , β)E θ0 [Z t k (β 0 ) ⊤ | X t k-1 ]) = 1 N N k=1
R(h 2 , X t k-1 )

P θ 0 --→ 0,

  Kronecker delta function is denoted by δ j i . For an open set A, the bar A indicates closure. We use θ = to indicate equality up to an additive constant that does not depend on θ. We write P

  the parameter space with Θ β and Θ Σ being two open convex bounded subsets of R r and R d×d , respectively. Functions F, N : R d × Θ β → R d are locally Lipschitz, and A, b are defined on Θ β and take values in R d×d and R d , respectively. Parameter matrix Σ takes values in R d×d .

Figure 1 :

 1 Figure 1: An example trajectory of the stochastic Lorenz system (20) starting at (0, 1, 0) for N = 10000 and h = 0.005. The first row shows the evolution of the individual components X, Y , and Z. The second row shows the evolution of component pairs: (Y, Z), (X, Z) and (X, Y ). Parameters are p = 10, r = 28, c = 8/3, σ 2 1 = 1, σ 2 2 = 2 and σ 2 3 = 1.5.

⊤N.

  The proof of Theorem 5.1 is in Section 7.3. Theorem 5.1 Let Assumptions (A1)-(A6) hold, X be the solution of (1), and θ N = ( β N , ΣΣ ⊤ N ) be the estimator that minimizes one of objective functions (22) or (23). If h → 0 and N h → ∞, then, βN

Figure 2 :

 2 Figure 2: Comparing the absolute relative error (ARE) as a function of increasing h for seven different estimators in the stochastic Lorenz system. The estimators are obtained for sample sizes of N = 10000. The y-axis is on a log scale.

Figure 3 :

 3 Figure 3: Comparison of the normalized distributions of ( θNθ 0 ) ⊘ θ 0 (where ⊘ is the element-wise division) in the Lorenz system for the S mix , S avg , and LL estimators for N = 10000. Each column represents one parameter, and each row represents one value of the discretization step h. A black dot with a vertical bar in each violin plot represents the mean and the standard deviation.

Figure 4 :

 4 Figure 4: Running times as a function of N for different estimators of the Lorenz system. Each column shows one value of h. On the x-axis is the sample size N , and on the y-axis is the running time in seconds.

Figure 5 :

 5 Figure 5: Comparing distributions of θNθ 0 for the S mix estimator with theoretical asymptotic distributions (28) for each parameter (columns), for h = 0.01 and N ∈ {1000, 5000, 10000} (colors). The black lines correspond to the theoretical asymptotic distributions computed from data and true parameters for N = 10000 and h = 0.01.

Figure 6 :

 6 Figure6: Comparing distributions of θNθ 0 for the LT mix estimator with theoretical asymptotic distributions (28) for each parameter (columns), for h ∈ {0.005, 0.01} (rows) and N ∈ {1000, 5000, 10000} (colors). The black lines correspond to the theoretical asymptotic distributions computed from data and true parameters for N = 10000 and corresponding h.

  Lemma S1.3 (Generalized Grönwall's inequality, Lemma 2.3 in Tian and Fan (2020)) Let p > 1 and b > 0 be constants, and let a : (0, +∞) → (0, +∞) be a continuous function. Ifu(t) ≤ a(t) + b t 0 u p (s) ds, then u(t) ≤ a(t) + (κ 1-p (t) -(p -1)2 p-1 bt)1 1-p and κ 1-p (t) > (p -1)2 p-1 bt, where κ(t) := 2 p-1 b t 0 a p (s) ds.

Table 1

 1 provides a comprehensive overview of estimator properties, finite sample performance, and required model assumptions for the most prominent state-of-the-art methods. While asymptotic properties might be similar in most cases, the finite sample properties are often different. The table also includes the Lie-Trotter (LT) and the Strang (S) splitting estimators, which we propose in this paper. The comparison encompasses four key characteristics: (1) Diffusion coefficient allowed in the model class, distinguishing between additive and general noise; (2) Asymptotic

regime, the conditions needed to prove the asymptotic properties; (3) Implementation, assessing the complexity of implementation, dependence on model dimension and parameter optimization time; and (4) Finite sample properties, evaluating performance for fixed sample size N and discretization step size h.

This Section presents the simulation study of the Lorenz system, illustrating the theory and comparing the proposed estimators with other likelihood-based estimators from the literature. We briefly recall the estimators, describe the simulation process and the optimization in the programming language R (R Core Team, 2022), and present and analyze the results.

Acknowledgement

This work has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 956107, "Economic Policy in Complex Environments (EPOC)"; Novo Nordisk Foundation NNF20OC0062958; and Independent Research Fund Denmark | Natural Sciences 9040-00215B.

Using that F (i) (x) = (A(x-b)) (i) + N (i) (x), ∂F (i) (x) ∂x = (A (i) ) ⊤ + ∇N (i) (x) and H F (i) (x) = H N (i) (x), the expectation of the true process rewrites as:

The final equation coincides with equation (S3) only up to order R(h, x). Despite the reversed LT splitting having the term with ΣΣ ⊤ at the order h 2 , the coefficients do not match. Thus, to obtain order R(h 2 , x), the condition (N(x)) ⊤ ∇F (i) (x) -1 2 Tr(ΣΣ ⊤ H N (i) (x)) = (F(x)) ⊤ ∇N (i) (x), must hold for all i = 1, . . . , d. Given Assumption (SA), the condition for achieving a higher one-step convergence order remains equivalent to the case of the original LT splitting.

S1.2 Proof for the Strang Splitting

We continue employing the Taylor expansion to establish the numerical properties of the Strang splitting approximation. To begin, we introduce a helpful Lemma S1.1 regarding the approximation of the composition of the mean function µ h and the nonlinear solution f h/2 . Lemma S1.1 expends µ h (f h/2 (x)) around h = 0, in various ways, each retaining the crucial terms necessary for the subsequent proofs.

Lemma S1.1 For the mean function µ h and the nonlinear solution f h/2 the following three identities hold:

Proof We prove only the first two identities, as the last one follows the same reasoning. Utilizing the definition of µ h , its Taylor expansion, and the expansion of f h/2 , we obtain:

, which concludes the first part. For the second part, Proposition 2.2 gives x). This leads to:

Proof of Proposition 3.6 We begin by introducing a new function of x, arising from the third property of Lemma S1.1:

Then, for a generic random vector X we use Proposition 2.2 and Lemma S1.1 to write:

Consequently, we expand:

The term

, while the terms with only one ξ h have zero means. Thus,

Lemma S1.4 (Lemma 9 in Genon-Catalot and Jacod (1993)) Let (X N k ) N ∈N,1≤k≤N be a triangular array with each row N adapted to a filtration (G N k ) 1≤k≤N , and let U be a random variable. If

Lemma S1.5 (Proposition A1 in [START_REF] Gloter | Parameter Estimation for a Discretely Observed Integrated Diffusion Process[END_REF]) Let S N (ω, θ) be a sequence of measurable real-valued functions defined on Ω × Θ, where (Ω, F, P) is a probability space, and Θ is product of compact intervals of R. We assume that S N (•, θ) converges to a constant C in probability for all θ ∈ Θ; and that there exists an open neighbourhood of Θ on which S N (ω, •) is continuously differentiable for all ω ∈ Ω. Furthermore, we suppose that:

Then, S N (θ)

C uniformly in θ.

Lemma S1.6 (Lemma 3.1 in [START_REF] Yoshida | Asymptotic behavior of M-estimator and related random field for diffusion process[END_REF]) Let F ⊂ R d be a convex compact set, and let {ξ N (θ); θ ∈ F }, be a family of real-valued random processes for N ∈ N. If there exist constants p ≥ l > d and C > 0 such that for all θ, θ 1 and θ 2 , it holds:

(1)

Theorem S1.7 (Rosenthal's inequality, Theorem 2.12 in [START_REF] Hall | Martingale Limit Theory and Its Application[END_REF]) Let (X N k ) N ∈N,1≤k≤N be a triangular array with each row N adapted to a filtration (G N k ) 1≤k≤N and let:

be a martingale array. Then, for all p ∈ [2, ∞) there exist constants C 1 , C 2 such that:

Theorem S1.8 (Proposition 3.1. in [START_REF] Crimaldi | Convergence results for multivariate martingales[END_REF]) Let (X N,k ) N ∈N,1≤k≤N be a triangular array of d-dimensional random vectors, such that, for each N , the finite sequence (X N,k ) 1≤k≤N is a martingale difference array with respect to a given filtration (G N k ) 1≤k≤N such that:

(2)

U, for some non-random positive semi-definite matrix U,

Remark Instead of using the second condition of Theorem S1.8, Lemma S1.6 yields that it is sufficient to prove that, for all i, j = 1, ..., d, it holds:

Remark For a martingale difference array the conditional expectations need to be zero almost surely, i.e:

does not fulfil the previous condition. Hence, similar to the approach in Corollary 2.6 of [START_REF] Mcleish | Dependent Central Limit Theorems and Invariance Principles[END_REF], we need the following two additional conditions on (X N,k ) N ∈N,1≤k≤N :

] satisfies conditions of the previous theorem. To prove that the first condition is satisfied, we write:

We used the Doob's inequality for the last submartingale. To demonstrate the second condition we fix i, j to get:

The first term goes to U ij , and the last term goes to zero. To prove that middle terms also vanish, we use the following inequalities:

S1.5 Proof of Lemma 7.1

Lemma 7.1 plays a central role in demonstrating the consistency and asymptotic normality of the proposed estimators.

The lemma deals with the uniform convergence of multiple triangular arrays, and proving various aspects of it involves a range of technical tools and methods. Different parts of Lemma 7.1 require distinct strategies to establish appropriate bounds, which can be intricate. Once these bounds are established, we leverage the properties discussed in the preceding section.

For instance, when establishing point-wise convergence, we primarily rely on Lemma S1.4. On the other hand, for proving uniform convergence, we utilize both Lemma S1.5 and Lemma S1.6. Throughout the proof of Lemma 7.1, a recurring theme is to interpret quadratic forms as traces and exploit the cyclic property inherent to them. Additionally, we employ fundamental mathematical tools like the mean value theorem, the Cauchy-Schwartz inequality, and Hölder's inequality in various instances.

Furthermore, there are occasions where we require inequality for norms, particularly the Frobenius norm. To address this, we introduce the Frobenius inner product of matrices M 1 and M 2 in R n×m as ⟨M 1 , M 2 ⟩ F := Tr(M ⊤ 1 M 2 ). Leveraging Hölder's inequality on Frobenius norm provides us with the following bound for the trace of a matrix product:

for h → 0. Similarly, we deduce that:

for N h → ∞. To prove uniform convergence, we use Lemma S1.6 along with Rosenthal's inequality from Theorem S1.7, resulting in:

The first term is bounded because of (S18). To bound the second term on the right-hand side, we use (S14). Then, for N h → ∞ and h → 0 and p > 2 it holds:

To conclude the proof of uniform convergence, we once again apply Rosenthal's inequality to get:

To bound the first term in (S19), we follow the reasoning from (S17) and start with:

Then, the rest is the same. Similarly, to bound the second term in (S19), we repeat derivations from (S17) to get:

Finally, (S18) and conclusions after (S17) complete the proof of 4.

Proof of 5. We introduce Y

Hölder's inequality, together with Cauchy-Schwartz inequality, Lemma 4.1 and (S13), yield:

for N h → ∞, h → 0. To prove the uniform convergence, we use Lemma S1.6. Again, it is enough to prove (S10) and (S11). Repeating the same steps as in the proof of (S14) leads to (S10). Similarly, to prove (S11) we repeat the same steps as in (S17) using Hölder's inequality, Cauchy-Schwartz inequality, and Lemma 4.1 with (S13).

Proof of 6. We introduce

for N h → ∞, h → 0. On the other hand,

which follows from derivations in (S20). To prove uniform convergence, we repeat the same approach as in the previous two proofs. Proof of 7. First, we use the fact that E[g(X

Again, the proofs of (S10) and (S11) are the same as in property 3, with a distinction of rewriting:

S1.6 Proof of the Asymptotic Normality

In this section, we distinguish between the true parameter θ 0 and a generic parameter θ. To complete the proof of asymptotic normality, we need to prove Lemma 7.4. The proof of this lemma is technical and involves bounding the sums of triangular arrays in such a way that the bound converges to zero in probability P θ0 as h → 0, N h → ∞, and N h 2 → 0. Unlike in the previous proof, this time we do not require uniform convergence.

Proof of Lemma 7.4 We begin by expanding η

k to differentiate between terms that vanish and those that do not in the limits:

Proof of (i). Let us begin by examining the limit of the expectation of sup 1≤k≤N |η

In equation (S21), all the involved functions are bounded, and the term with the largest order is R(

). The remaining terms converge to zero. Moreover, terms with coefficients 1 √ N h take the form Z t k (β 0 ) ⊤ (ΣΣ ⊤ 0 ) -1 g, where g is a vector-valued function of either X t k-1 or X t k . Their expected values are bounded by R(h, X t k-1 ) at most. Thus, the dominant order becomes R( h/N , X t k-1 ), which indeed converges to zero. We proceed to analyze the limit of the expectation of sup 1≤k≤N |ζ

N,k (θ 0 ), as defined in the paper, has an order R(1/ √ N h 2 , X t k-1 ). Upon calculating its expected value, we obtain an order of R(h, X t k-1 ). This concludes the proof of (i).

To establish limits (ii)-(v), we need to calculate the expectations of η

Similarly, from:

). Then, combining the previous, we get:

Now, we prove limit (vi). Here, we focus on the terms of order 1/

N,k which are the only ones that will not converge to zero when multiplying η

In the previous calculations, we introduced a new notation ψ

N,k (θ 0 ) and again focus only on the terms with coefficient 1/N h:

In the previous equation, we must show that the sum of expectations of all the terms except the first converges to zero. We only prove this for the second row; the rest follows analogously. Due to the definition of ψ i , it is clear that

), for all p ≥ 1. Then, we use property (S13) to obtain:

Finally, we use Lemma 4.2 to get:

(Tr(∂ βi 2 F(X t k-1 ; β 0 ) ⊤ (ΣΣ ⊤ 0 ) -1 ∂ βi 1 F(X t k-1 ; β 0 )) + R( √ h, X t k-1 ))

To prove (vii) we use Corollary 3.8:

Now, we use the expectation of a product of two quadratic forms of normally distributed random vectors (see for example Section 2 in [START_REF] Kumar | Expectation of Product of Quadratic Forms[END_REF]) to get:

This proves (vii). We omit the proofs of (viii)-(xi) since they follow the same pattern. Namely, we find the leading term and ensure it goes to zero. For the expectations of squares, we can apply the same approach with a product of two quadratic forms.