Alexandra Chronopoulou

Samy Tindel

ON INFERENCE FOR FRACTIONAL DIFFERENTIAL EQUATIONS

Keywords: 2010 Mathematics Subject Classification. Primary 60H35; Secondary 60H07, 60H10, 65C30, 62M09 Fractional Brownian motion, Stochastic differential equations, Malliavin calculus, Inference for stochastic processes

Based on Malliavin calculus tools and approximation results, we show how to compute a maximum likelihood type estimator for a rather general differential equation driven by a fractional Brownian motion with Hurst parameter H > 1/2. Rates of convergence for the approximation task are provided, and numerical experiments show that our procedure leads to good results in terms of estimation.

Introduction

In this introduction, we first try to motivate our problem and outline our results. We also argue that only a part of the question can be dealt with in a single paper. We briefly sketch a possible program for the remaining tasks in a second part of the introduction.

1.1. Motivations and outline of the results. The inference problem for diffusion processes is now a fairly well understood problem. In particular, during the last two decades, several advances have allowed to tackle the problem of inference based on discretely observed diffusions [START_REF] Durham | Numerical techniques for maximum likelihood estimation of continuous-time diffusion processes[END_REF][START_REF] Pedersen | Consistency and asymptotic normality of an approximate maximum likelihood estimator for discretely observed diffusion processes[END_REF][START_REF] Sorensen | Parametric inference for discretely sampled stochastic differential equations[END_REF], which is of special practical interest.

More specifically, consider a family of stochastic differential equations of the form

Y t = a + t 0 µ(Y s ; θ) ds + d l=1 t 0 σ l (Y s ; θ) dB l s , t ∈ [0, T ], (1)
where a ∈ R m , µ(•; θ) : R m → R m and σ(•; θ) : R m → R m,d are smooth enough functions, B is a d-dimensional Brownian motion and θ is a parameter varying in a subset Θ ⊂ R q . If one wishes to identify θ from a set of discrete observations of Y , most of the methods which can be found in the literature are based on (or are closely linked to) the maximum likelihood principle. Indeed, if B is a Brownian motion and Y is observed at some equally distant instants t i = iτ for i = 0, . . . , n, then the log-likelihood of a sample (Y t 1 , . . . , Y tn ) can be expressed as

ℓ n (θ) = n i=1 ln p τ, Y t i-1 , Y t i ; θ , (2)
where p stands for the transition semi-group of the diffusion Y . If Y enjoys some ergodic properties, with invariant measure ν θ 0 under P θ 0 , then we get a.s.-lim

n→∞ 1 n ℓ n (θ) = E θ 0 [p (τ, Z 1 , Z 2 ; θ)] J θ 0 (θ), (3)
where Z 1 ∼ ν θ 0 and L(Z 2 | Z 1 ) = p(τ, Z 1 , • ; θ). Furthermore, it can be shown in a general context that θ → J θ 0 (θ) admits a maximum at θ = θ 0 . This opens the way to a MLE analysis which is similar to the one performed in the case of i.i.d observations, at least theoretically. However, in many interesting cases, the transition semi-group p is not amenable to explicit computations, and thus expression [START_REF] Baudoin | Gradient Bounds for Solutions of Stochastic Differential Equations Driven by Fractional Brownian Motions[END_REF] has to be approximated in some sense. The most common approach, advocated for instance in [START_REF] Pedersen | Consistency and asymptotic normality of an approximate maximum likelihood estimator for discretely observed diffusion processes[END_REF], is based on a linearization of each p(τ, Y t i-1 , Y t i ; θ), which transforms it into a Gaussian density

N Y t i-1 + µ(Y t i-1 ; θ) τ, σσ * (Y t i-1 ; θ) τ .
This linearization procedure is equivalent to the approximation of equation ( 1) by an Euler (first order) numerical scheme. Refinements of this procedure, based on Milstein type discretizations, are proposed in [START_REF] Durham | Numerical techniques for maximum likelihood estimation of continuous-time diffusion processes[END_REF].

Some special situations can be treated differently (and often more efficiently): for instance, in case of a constant diffusion coefficient, the continuous time likelihood can be computed explicitly by means of Girsanov's theorem. When the dimension of the driving Brownian motion B is d = 1, one can also apply Itô's formula in order to be back to an equation with constant diffusion coefficient, or use Doss-Sousman representation of solutions to [START_REF] Baudoin | Operators associated with a stochastic differential equation driven by fractional Brownian motions[END_REF]. Let us also mention that statistical inference for SDEs driven by Lévy processes is currently intensively investigated, with financial motivations in mind.

The current article is concerned with the estimation problem for equations of the form [START_REF] Baudoin | Operators associated with a stochastic differential equation driven by fractional Brownian motions[END_REF], when the driving process B is a fractional Brownian motion. Let us recall that a fractional Brownian motion B with Hurst parameter H ∈ (0, 1), defined on a complete probability space (Ω, F , P), is a d-dimensional centered Gaussian process. Its law is thus characterized by its covariance function, which is given by

E B i t B j s = 1 2 t 2H + s 2H -|t -s| 2H 1 (i=j) , s, t ∈ R + .
The variance of the increments of B is then given by

E B i t -B i s 2 = |t -s| 2H , s, t ∈ R + , i = 1, . . . , d,
and this implies that almost surely the fBm paths are γ-Hölder continuous for any γ < H. Furthermore, for H = 1/2, fBm coincides with the usual Brownian motion, converting the family {B H ; H ∈ (0, 1)} into the most natural generalization of this classical process.

In the last decade, some important advances have allowed to solve [START_REF] Nualart | Differential equations driven by fractional Brownian motion[END_REF][START_REF] Zähle | Integration with respect to fractal functions and stochastic calculus I[END_REF] and understand [START_REF] Hu | Differential equations driven by Hölder continuous functions of order greater than 1/2[END_REF][START_REF] Nualart | Malliavin calculus for stochastic differential equations driven by a fractional Brownian motion[END_REF] differential systems driven by fBm for H ∈ (1/2, 1). The rough paths machinery also allows to handle fBm with H ∈ (1/4, 1/2), as nicely explained in [START_REF] Friz | Multidimensional stochastic processes as rough paths[END_REF][START_REF] Gubinelli | Controlling rough paths[END_REF][START_REF] Lejay | An Introduction to Rough Paths[END_REF][START_REF] Lyons | System control and rough paths[END_REF]. However, the irregular situation H ∈ (1/4, 1/2) is not amenable to useful moments estimates for the solution Y to (1) together with its Jacobian (that is the derivative with respect to the initial condition). This is why we concentrate, in the sequel, on the simpler case H > 1/2 for our estimation problem. In any case, many real world noisy systems are currently modeled by equations like [START_REF] Baudoin | Operators associated with a stochastic differential equation driven by fractional Brownian motions[END_REF] driven by fBm, and this is particularly present in the Biophysics literature, as assessed by [START_REF] Kou | Generalized Langevin equation with fractional Gaussian noise: subdiffusion within a single protein molecule[END_REF][START_REF] Odde | Stochastic dynamics of the nerve growth cone and its microtubules during neurite outgrowth[END_REF], or for Finance oriented applications as in [START_REF] Cheridito | Arbitrage in fractional Brownian motion models[END_REF][START_REF] Guasoni | No arbitrage under transaction costs, with fractional Brownian motion and beyond[END_REF][START_REF] Hu | Fractional white noise calculus and applications to finance, Infnite dimentional analysis[END_REF][START_REF] Hu | Optimal sonsumption and portfolio in a Black-Scholes market driven by fractional Brownian motion, Infnite dimentional analysis[END_REF][START_REF] Rogers | Arbitrage with fractional Brownian motion[END_REF][START_REF] Willinger | Stock market prices and long-range dependence[END_REF]. This leads to a demand for rigorous estimation procedures for SDEs driven by fractional Brownian motion, which is the object of our paper.

Concerns about the inference problem for fractional diffusion processes started a decade ago with the analysis of fractional Ornstein-Uhlenbeck processes [START_REF] Kleptsyna | Statistical analysis of the fractional Ornstein-Uhlenbeck type process[END_REF]. Then a more recent representative set of references on the topic includes [START_REF] Papavasiliou | Parameter estimation for rough differential equations[END_REF][START_REF] Tudor | Statistical aspects of the fractional stochastic calculus[END_REF]. More specifically, [START_REF] Tudor | Statistical aspects of the fractional stochastic calculus[END_REF] handles the case of a one-dimensional equation of the form

Y t = a + θ t 0 µ(Y s ) ds + B t , t ∈ [0, T ], ( 4
)
where µ is regular enough, and where B is a fBm with H ∈ (0, 1). The simple dependence on the parameter θ and the fact that an additive noise is considered enables the use of Girsanov's transform in order to get an exact expression for the MLE. Convergence of the estimator is then obtained through an extensive use of Malliavin calculus.

As far as [START_REF] Papavasiliou | Parameter estimation for rough differential equations[END_REF] is concerned, it is focused on the case of a polynomial equation, for which the exact moments of the solution can be computed. The estimator relies then on a generalization of the moment method, which tries to fit empirical moments of the solution with their theoretical value. The range of application of this method is however confined to specific situations, for the following reasons:

• It assumes that N independent runs of equation ( 1) can be obtained, which is usually not the case. • It hinges on multiple integrals computations, which are time consuming and are avoided in most numerical schemes.

As can be seen from this brief review, parameter estimation for rough equations is still in its infancy. We shall also argue that it is a hard problem.

Indeed, if one wishes to transpose the MLE methods used for diffusion processes to the fBm context, an equivalent of the log-likelihood functions (2) should first be produced. But the covariance structure of B is quite complex and the attempts to put the law of Y defined by (1) into a semigroup setting are cumbersome, as illustrated by [START_REF] Baudoin | Operators associated with a stochastic differential equation driven by fractional Brownian motions[END_REF][START_REF] Hairer | Ergodicity theory of SDEs with extrinsic memory[END_REF][START_REF] Neuenkirch | Trees and asymptotic developments for fractional diffusion processes[END_REF]. We have thus decided to consider a highly simplified version of the log-likelihood. Namely, still assuming that Y is observed at a discrete set of instants 0

< t 1 < • • • < t n < ∞, set ℓ n (θ) = n i=1 ln (f (t i , Y t i ; θ)) , (5)
where we suppose that under P θ the random variable Y t i admits a density z → f (t i , z; θ). Notice that in case of an elliptic diffusion coefficient σ the density f (t i , •; θ) is strictly positive, and thus expression (5) makes sense by a straightforward application of [START_REF] Friz | Multidimensional stochastic processes as rough paths[END_REF]Proposition 19.6]. However, the successful replication of the strategy implemented for Brownian diffusions (that we have tried to summarize above) relies on some highly non trivial questions: existence of an invariant measure for equation [START_REF] Baudoin | Operators associated with a stochastic differential equation driven by fractional Brownian motions[END_REF], rate of convergence to this invariant measure, convergence of expressions like (5), characterization of the limit in terms of θ as in [START_REF] Baudoin | Gaussian bounds for the density of solutions of stochastic differential equations driven by fractional Brownian motions[END_REF], to mention just a few. We shall come back to these considerations in the next section, but let us insist at this point on the fact that all those questions would fit into a research program over several years. Our aim in this paper is in a sense simpler: we assume that quantities like (5) are meaningful for estimation purposes. Then we shall implement a method which enables to compute ℓ n (θ) and optimize it in θ, producing thus a pseudo MLE estimator. We focused first on this specific aspect of the problem for the following reasons:

(1) From a statistical point of view, it is obviously essential to obtain a computationally efficient estimation procedure. This will allow us for instance to evaluate numerically the accuracy of our method. (2) The procedure itself is nontrivial, and requires the use of advanced stochastic analysis tools: probabilistic representation of the density, Malliavin type integration by parts, Stratonovich-Skorohod correction terms, discretization of systems of pathwise stochastic differential equations for instance. We have thus decided to tackle the implementation issues first. If it turns out to be satisfying, we shall then try to proceed to a full justification of our method.

Let us also mention that it might not be clear to the reader that ℓ n (θ) can be meaningful in terms of statistical estimation, since it only involves evaluations at single points Y t i . However our numerical experiments indicate that this quantity behaves nicely for our purposes. Moreover, it will become clear from the forthcoming computations that our methodology can be extended to handle quantities like

ln (θ) := n i=1 ln f (t i , t i+1 , Y t i , Y t i+1 ; θ) ,
where f (s, t, x, z; θ) stands for the density of the couple (Y s , Y t ). This kind of pseudo log-likelihood is obviously closer in spirit to the diffusion case. Densities of tuples could also be considered at the price of technical complications.

Let us now try to give a flavor of the kind of result we shall obtain in this article, in a very loose form: Theorem 1.1. Consider Equation (1) driven by a d-dimensional fractional Brownian motion B with Hurst parameter H > 1/2. Assume µ and σ are smooth enough coefficients, and that σσ * is strictly elliptic. For a sequence of times t 0 < • • • < t n < ∞, let y t i , i = 1, . . . , n be the corresponding observations. Then: (i) The gradient of the log-likelihood function admits the following probabilistic representation:

∇ l ℓ n (θ) = n i=1 V i (θ) W i (θ) , with W i (θ) = E 1 (Yt i (θ)>yt i ) H (1,...,m) Y t i (θ) (6)
where H (j 1 ,...,jn) (Y t i (θ)) is an expression involving Malliavin derivatives an Skorohod integrals of Y (θ). A similar expression is also available for V i (θ).

(ii) A computational procedure is constructed in order to obtain H (1,...,m) (Y t i (θ)) in a suitable way.

(iii) When Y t is replaced by its Euler scheme approximation with step T /M and expected values in [START_REF] Dalang | Potential theory for hyperbolic SPDEs[END_REF] are approximated thanks to N Monte Carlo steps, we show that

• N can be chosen in function of M in an optimal way (see Proposition 4.7).

• The corresponding approximation of ∇ l ℓ n (θ) converges to the real one with rate n -(2γ-1) for any 1/2 < γ < H.

All those results are stated in a more rigorous way in the remainder of the article.

Here is how our article is structured: we give some preliminaries and notations on Young and Malliavin calculus for fractional Brownian motion at Section 2. The probabilistic representation for the log-likelihood is given at Section 3. Discretization procedures are designed at Section 4, and finally numerical examples are given at Section 5.

1.2.

Remaining open problems. We emphasized above the fact that only a part of the problem at stake was going to be solved in the current article. We now briefly sketch the remaining tasks to be treated.

The most important obstacle in order to fully justify our methodology is to get a suitable convergence theorem for ℓ n (θ)/n, where ℓ n (θ) is defined by [START_REF] Cheridito | Arbitrage in fractional Brownian motion models[END_REF]. In a natural way, this should be based on some strong ergodicity properties for Y t . After a glance at the literature on ergodicity for fractional systems, one can distinguish two cases: [START_REF] Hairer | Ergodicity of stochastic differential equations driven by fractional Brownian motion[END_REF], with a (presumably non optimal) rate of convergence t -1/8 . (ii) For a general smooth and elliptic coefficient σ, only the uniqueness of the invariant measure is shown in [START_REF] Hairer | Ergodicity theory of SDEs with extrinsic memory[END_REF], with an interesting extension to the hypoelliptic case in [START_REF] Hairer | Ergodicity of hypoelliptic SDEs driven by fractional Brownian motion[END_REF]. Nothing is known about the convergence of L(Y t ), not to mention rates. This brief review already indicates that the convergence to invariant measures is still quite mysterious for fractional differential equations, at least for a non constant coefficient σ. Moreover, recall that if ν(θ) stands for the invariant measure corresponding to the system with coefficients µ(•; θ), σ(•; θ), we also wish to retrieve some information on the dependence θ → ν(θ) (See [START_REF] Hairer | A simple framework to justify linear response theory[END_REF] for some partial results in this direction).

(i) When σ(•; θ) is constant, the convergence of L(Y t ) as t → ∞ is established in
Let us mention another concrete problem: even in the case of a constant σ, the convergence of L(Y t ) to an invariant measure ν(θ) is proven in [START_REF] Hairer | Ergodicity of stochastic differential equations driven by fractional Brownian motion[END_REF] in the total variation sense. In terms of the density p(t, x; θ) of Y t , it means that p(t, •; θ) converges to the density of ν in L 1 topology. However, in order to get a limit for ℓ n (θ)/n, one expects to use at least a convergence in some Sobolev space W α,p for α, p large enough.

One possibility in order to get this sharper convergence is to bound first the density p(t, •; θ) in another Sobolev space W α ′ ,p ′ and then to use interpolation theory. It seems thus sufficient to obtain Gaussian bounds on p(t, •; θ), uniformly in t. In case of Brownian diffusions, these Gaussian bounds are obtained by analytic tools, thanks to the Markov property. This method being obviously not available for systems driven by fBm, a possible inspiration is contained in the upper Gaussian bounds for the stochastic wave equation which can be found in [START_REF] Dalang | Potential theory for hyperbolic SPDEs[END_REF]. The latter technical results stem from an intensive use of Malliavin calculus, which should also be invoked in our case, and notice the recent efforts [START_REF] Baudoin | Gradient Bounds for Solutions of Stochastic Differential Equations Driven by Fractional Brownian Motions[END_REF][START_REF] Baudoin | Gaussian bounds for the density of solutions of stochastic differential equations driven by fractional Brownian motions[END_REF] in this direction.

Finally, let us mention that it seems possible to produce some reasonable convergent parametric estimators for equations driven by fBm in a rather general context. Among the methods which can be adapted from the diffusion case with the current stochastic analysis techniques, let us mention the least square estimator of [START_REF] Kasonga | Maximum likelihood theory for large interacting systems[END_REF], as well as the local asymptotic normality property shown in [START_REF] Gobet | Local asymptotic mixed normality property for elliptic diffusion: a Malliavin calculus approach[END_REF]. However, it seems obvious that the road to a complete picture of parameter estimation for stochastic equations driven by fBm is still hard and long. We hope to complete it in some subsequent communications.

Preliminaries and notations

As mentioned in the introduction, we are concerned with equations driven by a ddimensional fractional Brownian motion B. We recall here some basic facts about the way to solve those equations, and some Malliavin calculus tools which will be needed later on. Let us introduce first some general notation for Hölder type spaces: Notation 2.1. We will denote by C α (V ) the set of V -valued α-Hölder functions for any α ∈ (0, 1), and by C n b (U; V ) the set of n times differentiable functions, bounded together with all their derivatives, from U to V . In the previous notation, U and V stand for two finite dimensional vector spaces. The state space V can be omitted for notational sake when its value is non ambiguous. When we want to stress the fact that we are working on a finite interval [0, T ], we write C α T (V ) for the space of α-Hölder functions f from [0, T ] to V . The corresponding Hölder norms shall be denoted by f α,T .

2.1. Differential equations driven by fBm. Recall that the equation we are interested in is of the form [START_REF] Baudoin | Operators associated with a stochastic differential equation driven by fractional Brownian motions[END_REF]. Before stating the assumptions on our coefficients we need an additional notation:

Notation 2.2. For n, p ≥ 1, a function f ∈ C p (R n ; R) and any tuple (i 1 , . . . i p ) ∈ {1, . . . , d} p , we set ∂ i 1 ...ip f for ∂ p f ∂x i 1 ...∂x ip .
Similarly, consider a function g θ ∈ C p (Θ; R), for n, p ≥ 1 and a vector of parameters θ ∈ Θ ⊂ R q . For any tuple (i 1 , . . . i p ) ∈ {1, . . . , q} p , we set ∇ i 1 ...ip g i θ for

∂ p g i θ ∂θ i 1 ...∂θ ip , where i = 1, . . . , n.

Using this notation, we work under the following set of assumptions:

∇ l i 1 •••i l µ(•; θ) ∞ + ∇ l i 1 •••i l σ(•; θ) ∞ < ∞.
When equation ( 1) is driven by a fBm with Hurst parameter H > 1/2 it can be solved, thanks to a fixed point argument, with the stochastic integral interpreted in the (pathwise) Young sense (see e.g. [START_REF] Gubinelli | Controlling rough paths[END_REF]). Let us recall that Young's integral can be defined in the following way:

Proposition 2.4. Let f ∈ C γ , g ∈ C κ with γ + κ > 1, and 0 ≤ s ≤ t ≤ 1.
Then the integral t s g ξ df ξ is well-defined as limit of Riemann sums along partitions of [s, t]. Moreover, the following estimation is fulfilled:

t s g ξ df ξ ≤ C f γ g κ |t -s| γ , (7)
where the constant C only depends on γ and κ. A sharper estimate is also available:

t s g ξ df ξ ≤ |g s | f γ |t -s| γ + c γ,κ f γ g κ |t -s| γ+κ . (8)
With this definition in mind and under assumptions 2.3, we can solve our differential system of interest, and the following moments bounds are proven in [START_REF] Friz | Multidimensional stochastic processes as rough paths[END_REF][START_REF] Hu | Differential equations driven by Hölder continuous functions of order greater than 1/2[END_REF]: Proposition 2.5. Consider a fBm B with Hurst parameter H > 1/2. Then:

(1) Under Hypothesis 2.3, equation [START_REF] Baudoin | Operators associated with a stochastic differential equation driven by fractional Brownian motions[END_REF] driven by B admits a unique β-Hölder continuous solution Y , for any β < H.

(2) Furthermore,

Y T,β ≤ |a| + c f,T B 1/β β,T .
(3) If we denote by Y a the solution to (1) with initial condition a, then

Y b -Y a T,β ≤ |b -a| exp c f,T B 1/β β,T .
(4) If we only assume that f has linear growth, with ∇f, ∇ 2 f bounded, the following estimate holds true:

sup t∈[0,T ] |Y t | ≤ (1 + |a|) exp c f,T B 1/β β,T .
Remark 2.6. The framework of fractional integrals is used in [START_REF] Hu | Differential equations driven by Hölder continuous functions of order greater than 1/2[END_REF] in order to define integrals with respect to B. It is however easily seen to be equivalent to the Young setting we have chosen to work with.

Some differential calculus rules for processes controlled by fBm will also be useful in the sequel: andc, ĉ ∈ C T (R) (all these assumptions are understood in the almost sure sense). Define two processes z, ẑ on [0, T ] by

Proposition 2.7. Let B be a d-dimensional fBm with Hurst parameter H > 1/2. Con- sider a, â ∈ R, b, b ∈ C α T (R d ) with α + H > 1,
z t = a + d j=1 t 0 b j u dB j u + t 0 c u du, and ẑt = â + d j=1 t 0 bj u dB j u + t 0 ĉu du.
Then for t ∈ [0, T ], one can decompose the product z t ẑt into

z t ẑt = a â + n j=1 t 0 ẑu b j u + z u bj u dB j u + t 0 [z u ĉu + ẑu c u ] du,
where all the integrals with respect to B are understood in the Young sense.

The proof of this elementary and classical result is omitted here. See [28, Proposition 2.8] for the proof of a similar rule.

Malliavin calculus techniques.

Our representation of the density for the solution to (1) obviously relies on Malliavin calculus tools that we proceed now to recall. As already mentioned in the introduction, on a finite interval [0, T ] and for some fixed H ∈ (1/2, 1), we consider (Ω, F , P ) the canonical probability space associated with a fractional Brownian motion with Hurst parameter H. That is, Ω = C 0 ([0, T ]; R d ) is the Banach space of continuous functions vanishing at 0 equipped with the supremum norm, F is the Borel sigma-algebra and P is the unique probability measure on Ω such that the canonical process B = {B t , t ∈ [0, T ]} is a d-dimensional fractional Brownian motion with Hurst parameter H. Remind that this means that B has d independent coordinates, each one being a centered Gaussian process with covariance R H (t, s) = 1 2 (s 2H + t 2H -|t -s| 2H ). 2.2.1. Functional spaces. Let E be the space of d-dimensional elementary functions on [0, T ]:

E = f = (f 1 , . . . , f d ); f j = n j -1 i=0 a j i 1 [t j i ,t j i+1 ) , 0 = t 0 < t j 1 < • • • < t j n j -1 < t j n j = T, for j = 1, . . . , d . (9)
We call H the completion of E with respect to the semi-inner product

f, g H = d i=1 f i , g i H 0 , where 1 [0,t] , 1 [0,s] H 0 := R(s, t), s, t ∈ [0, T ].
Then, one constructs an isometry

K * H : H → L 2 ([0, 1]; R d ) such that K * H 1 [0,t 1 ] , . . . , 1 [0,t d ] = 1 [0,t 1 ] K H (t 1 , •), . . . , 1 [0,t d ] K H (t d , •)
, where the kernel K H is given by

K H (t, s) = c H s 1 2 -H t s (u -s) H-3 2 u H-1 2 du
and verifies that R H (t, s) = s∧t 0 K H (t, r)K H (s, r) dr, for some constant c H . Moreover, let us observe that K * H can be represented in the following form:

for ϕ = (ϕ 1 , . . . , ϕ d ) ∈ H, we have K * H ϕ K * H ϕ = K * H ϕ 1 , . . . , K * H ϕ d , where [K * H ϕ i ] t = 1 t ϕ i r ∂ r K H (r, t) dr.

Malliavin derivatives.

Let us start by defining the Wiener integral with respect to B: for any element f in E whose expression is given as in ( 9), we define the Wiener integral of f with respect to B as

B(f ) := d j=1 n j -1 i=0 a j i (B j t j i+1 -B j t j i ) .
We also denote this integral as T 0 f t dB t , since it coincides with a pathwise integral with respect to B.

For θ : R → R, and j ∈ {1, . . . , d}, denote by θ [j] the function with values in R d having all the coordinates equal to zero except the j-th coordinate that equals to θ. It is readily seen that E B 1

[j] [0,s) B 1 [k] [0,t) = δ j,k R s,t .
This definition can be extended by linearity and closure to elements of H, and we obtain the relation E [B(f ) B(g)] = f, g H , valid for any couple of elements f, g ∈ H. In particular, B(•) defines an isometric map from H into a subspace of L 2 (Ω).

We can now proceed to the definition of Malliavin derivatives. With this notation 2.2 in hand, let us consider S be the family of smooth functionals F of the form

F = f (B(h 1 ), . . . , B(h n )), (10)
where h 1 , . . . , h n ∈ H, n ≥ 1, and f is a smooth function with polynomial growth, together with all its derivatives. Then, the Malliavin derivative of such a functional F is the H-valued random variable defined by

DF = n i=1 ∂ i f (B(h 1 ), . . . , B(h n )) h i .
For all p > 1, it is known that the operator D is closable from L p (Ω) into L p (Ω; H) (see e.g. [32, Section 1]). We will still denote by D the closure of this operator, whose domain is usually denoted by D 1,p and is defined as the completion of S with respect to the norm

F 1,p := (E(|F | p ) + E( DF p H )) 1
p . It should also be noticed that partial Malliavin derivatives with respect to each component B j of B will be invoked: they are defined, for a functional F of the form (10) and j = 1, . . . , d, as

D j F = n i=1 ∂ i f (B(h 1 ), . . . , B(h n ))h [j] i ,
and then extended by closure arguments again. We refer to [32, Section 1] for the definition of higher derivatives and Sobolev spaces D k,p for k > 1. Another essential object related to those derivatives is the so-called Malliavin matrix of a R m -valued random variable F ∈ D 1,2 , defined by

γ F = DF i , DF j 1≤i,j≤m . (11)
2.2.3. Skorohod integrals. We will denote by δ the adjoint of the operator D (also referred to as the divergence operator). This operator is closed and its domain, denoted by Dom(δ), is the set of H-valued square integrable random variables u ∈ L 2 (Ω; H) such that

|E [ DF, u H ] | ≤ C F 2 , for all F ∈ D 1,2
, where C is some constant depending on u. Moreover, for u ∈ Dom(δ), δ(u) is the element of L 2 (Ω) characterized by the duality relationship:

E [F δ(u)] = E [ DF, u H ] , for any F ∈ D 1,2 . (12)
The quantity δ(u) is usually called Skorohod integral of the process u. Skorohod integrals are obviously analytic objects, not suitable for easy numerical implementations. However, they can be related to the Young type integrals introduced at Proposition 2.4. For this, we need to define another functional space as follows: Notation 2.8. We call |H| the space of measurable functions ϕ : [0, T ] → R d such that

ϕ 2 |H| := c H 1 0 1 0 |ϕ r ||ϕ u ||r -u| 2H-2 drdu < +∞,
where c H = H(2H -1), and we denote by •, • |H| the associated inner product. We also write D k,p (|H|) for the space of D k,p functionals with values in |H|.

The following proposition is then a slight extension of [32, Proposition 5.2.3]:

Proposition 2.9. Let {u ij t , t ∈ [0, 1]}, for i = 1, . . . , m and j = 1, . . . , d, be a stochastic process in D 1,2 (|H|) such that d j=1 1 0 1 0 |D j s u ij t | |t -s| 2H-2 dsdt < +∞ a.s. ( 13
)
We also assume that almost surely, u has β-Hölder paths with β +H > 1. Then the Young integral d j=1 T 0 u ij t dB j t exists and for all i = 1, . . . , m can be written as

d j=1 T 0 u ij t dB j t = δ(u i ) + d j=1 T 0 T 0 D j s u ij t |t -s| 2H-2 dsdt,
where δ(u) stands for the Skorohod integral of u.

Probabilistic expression for the log-likelihood

Recall that we are focusing on equation (1) driven by a d-dimensional fBm B, and that we have chosen to use expression ( 5) as a substitute to the log-likelihood function. We have thus reduced the initial maximization problem to the solution of ∇ l ℓ n (θ) = 0. This will be performed numerically by means of a root approximation algorithm.

Observe first that in order to define [START_REF] Cheridito | Arbitrage in fractional Brownian motion models[END_REF], the density of Y t (θ) must exist for any t > 0. Let us thus recall the classical setting (given in [START_REF] Hu | Differential equations driven by Hölder continuous functions of order greater than 1/2[END_REF]) under which Y t admits a smooth density: Hypothesis 3.1. Let µ and σ be coefficients satisfying Hypothesis 2.3. For ξ ∈ R m and θ ∈ Θ, set α(ξ) = σ(ξ, θ)σ * (ξ, θ). Then we assume that (i) For any k ≥ 0 and j 1 , . . . , j k ∈ {1, . . . , m} we have

sup θ∈Θ 2 l=0 1≤p 1 ,...,p l ≤q ∇ l p 1 •••p l ∂ k j 1 ,...,j k µ(•; θ) ∞ + ∇ l p 1 •••p l ∂ k j 1 ,...,j k σ(•; θ) ∞ ≤ c k ,
for a strictly positive constant c k .

(ii) There exists a strictly positive constant

ε such that α(ξ; θ)η, η R m ≥ ε|η| 2 R m for any couple of vectors η, ξ ∈ R m , uniformly in θ ∈ Θ.
Then the density result for Y t can be read as follows: Theorem 3.2. Consider the stochastic differential equation ( 1) with initial condition a ∈ R m . Assume Hypothesis 3.1 is satisfied. Then, for any t > 0 and θ ∈ Θ, the law of Y t (θ) admits a C ∞ density, denoted by f (t, •; θ), with respect to Lebesgue's measure.

In the sequel, we shall suppose that the density f (t, •; θ) exists without further mention, the aim of this section being to produce a probabilistic representation of f (t, •; θ) for computational purposes. To this aim, we shall first give the equations governing the Malliavin derivatives of the processes Y (θ) and ∇Y (θ), and then use a stochastic analysis formula in order to represent our log-likelihood. We separate these tasks in two different subsections.

3.1. Some Malliavin derivatives. This section is devoted to a series of preliminary lemmas which will enable to formulate our probabilistic representation of f (t, •; θ). Let us first introduce a notation which will prevail until the end of the paper: Notation 3.3. For a set of indices or coordinates (k 1 , . . . , k r ) of length r ≥ 1 and 1 ≤ j ≤ r, we denote by (k 1 , . . . , ǩj , . . . , k r ) the set of indices or coordinates of length r -1 where k j has been omitted.

We now give a general expression for the higher order derivatives of Y t , borrowed from [START_REF] Nualart | Malliavin calculus for stochastic differential equations driven by a fractional Brownian motion[END_REF].

Lemma 3.4. Assume Hypothesis 2.3 and 3.1 hold true. For n ≥ 1 and

(i 1 , . . . , i n ) ∈ {1, . . . , d} n , denote by D i 1 ,...,in Y i t (θ) the n th Malliavin derivative of Y i t (θ) with respect to the coordinates B i 1 , . . . , B in of B. Then D i 1 ,...,in Y i t (θ)
, considered as an element of H ⊗n , satisfies the following linear equation:

for t ≥ r 1 ∨ • • • ∨ r n , D i 1 ,...,in r 1 ,...,rn Y i t (θ) = n p=1 α i ip,i 1 .
..,ǐp,...,in (r p ; r 1 , . . . , řp , . . . , r n ; θ)

+ t r 1 ∨•••∨rn β i i 1 ,...,in (s; r 1 , . . . , r n ; θ) ds + d l=1 t r 1 ∨•••∨rn α i l,i 1 ,...,in (s; r 1 , . . . , r n ; θ) dB l s , (14)
where

α i j,i 1 ,...,in (s; r 1 , . . . , r n ; θ) = m k 1 ,...,kν =1 ∂ ν k 1 ...kν σ ij (Y s (θ); θ) D i(I 1 ) r(I 1 ) Y k 1 s (θ) . . . D i(Iν ) r(Iν ) Y kν s (θ)
β i i 1 ,...,in (s; r 1 , . . . , r n ; θ) = m k 1 ,...,kν =1 ∂ ν k 1 ...kν µ i (Y s (θ); θ) D i(I 1 ) r(I 1 ) Y k 1 s (θ) . . . D i(Iν ) r(Iν ) Y kν s (θ).
In the expressions above, the first sums are extended to the set of all partitions I 1 , . . . , I ν of {1, . . . , n} and for any subset K = {i 1 , . . . , i η } of {1, . . . , n} we set D i(K) r(K) for the derivative operator D i 1 ,...,iη r 1 ,...,rη . Notice that D i 1 ,...,in r 1 ,...,rn Y i t (θ) = 0 whenever t < r 1 ∨ • • • ∨ r n . The formulas above might seem intricate. The following example illustrate their use in a simple enough situation:

Example 3.5. The second order derivative D 1,3 r 1 ,r 2 Y 2 t (θ) can be computed as:

D 1,3 r 1 ,r 2 Y 2 t (θ) = α 2 1,3 (r 1 , r 2 ; θ) + α 2 3,1 (r 2 , r 1 ; θ) + t r 1 ∨r 2 β 2 1,3 (s, r 1 , r 2 ; θ) ds + d l=1 t r 1 ∨r 2 α 2 l,1,3 (s, r 1 , r 2 ; θ)dB l s ,
where

α 2 1,3 (r 1 , r 2 ; θ) = m k=1 ∂ k σ 21 (Y r 2 (θ); θ) D 3 r 2 Y k r 1 (θ), α 2 3,1 (r 2 , r 1 ; θ) = m k=1 ∂ k σ 23 (Y r 1 (θ); θ) D 1 r 1 Y k r 2 (θ)
and

β 2 1,3 (s, r 1 , r 2 ; θ) = ∂ 2 kk µ 2 (Y s (θ); θ)D 1,3 r 1 ,r 2 Y k s (θ) + ∂ 2 k 1 k 2 µ 2 (Y s (θ); θ)D 1 r 1 Y k 1 s (θ)D 3 r 2 Y k 2 s (θ), α 2 l,1,3 (s, r 1 , r 2 ; θ) = ∂ 2 kk σ 2l (Y s (θ); θ)D 1,3 r 1 ,r 2 Y k s (θ) + ∂ 2 k 1 k 2 σ 2l (Y s (θ); θ)D 1 r 1 Y k 1 s (θ)D 3 r 2 Y k 2 s (θ),
where we have used the convention of summation over repeated indices.

Our formula for the log-likelihood will also involve some derivatives of the process Y (θ) with respect to the parameter θ. The existence of this derivative is assessed below: Proposition 3.6. Under the same hypothesis as for Lemma 3.4, the random variable Y i t (θ) is a smooth function of θ for any t ≥ 0. We denote by ∇ l Y i t (θ) the derivative of Y i t (θ) with respect to the l th element of the vector of parameters θ. This process satisfies the following SDE:

∇ l Y i t (θ) = t 0 [∂ i µ i (Y u (θ); θ)∇ l Y i u (θ) + ∇ l µ i (Y u (θ); θ)]du + d j=1 t 0 [∂σ ij (Y u (θ); θ)∇ l Y i u (θ) + ∇ l σ ij (Y u (θ); θ)]dB j u .
Proof. The proof goes exactly along the same lines as for [34, Proposition 4], and the details are left to the reader.

We shall also need some equations governing the Malliavin derivatives of ∇ l Y (θ). This is the aim of the following lemma: Lemma 3.7. For any l ∈ {1, . . . , q} and n ≥ 1, the process ∇ l D i 1 ,...,in Y (θ) is n-times differentiable in the Malliavin calculus sense. Moreover, taking up the notations of Lemma 3.4, the process ∇ l D i 1 ,...,in Y i t (θ) satisfies the following linear equation:

for t ≥ r 1 ∨ • • • ∨ r n , ∇ l D i 1 ,...,in r 1 ,...,rn Y i t (θ) = n p=1
αi,l ip,i 1 ...,ǐp,...,n (r ip , r 1 , . . . , řp , . . . , r n ; θ)

+ t r 1 ∨•••∨rn βi,l i 1 ,...,in (s; r 1 , . . . , r n ; θ) ds + d l=1 t r 1 ∨•••∨rn αi,l l,i 1 ,...,in (s; r 1 , . . . , r n ; θ) dB l s ,
where αi,l j,i 1 ,...,in = ∇ l α i j,i 1 ,...,in and βi,l j,i 1 ,...,in = ∇ l β i i 1 ,...,in . More specifically, βi,p j,i 1 ,...,in is defined recursively by βi,p i 1 ,...,in (s; r 1 , . . . , r n ; θ)

= I 1 ∪...∪Iν m k 1 ,...,kν =1 ∇ p [∂ ν k 1 ...kν µ i (Y s (θ); θ)] D i(I 1 ) r(I 1 ) Y k 1 s (θ) • • • D i(Iν ) r(Iν ) Y kν s (θ) + ∂ ν k 1 ...kν µ i (Y s (θ); θ) ν p=1 ∇ p D i(Ip) r(Ip) Y kp s (θ) D i(I 1 ) r(I 1 ) Y k 1 s (θ) • • • D ǐ(Ip) ř(Ip) Y ǩp s (θ) • • • D i(Iν ) r(Iν ) Y kν s (θ) ,
where we have set

∇ p [∂ ν k 1 ...kν µ i (Y s (θ); θ)] = ∇ p ∂ ν k 1 ...kν µ i (Y s (θ); θ) + ∂∂ ν k 1 ...kν µ i (Y s (θ); θ)∇ p Y s (θ).
Notice that the same kind of equation (skipped here for sake of conciseness) holds true for the coefficients αi,l j,i 1 ,...,in . The next object we need for our calculations is the inverse of the Malliavin matrix γ Yt(θ) of Y t (θ). Recall that according to [START_REF] Friz | Multidimensional stochastic processes as rough paths[END_REF], the Malliavin matrix of Y t (θ) is defined by

γ t (θ) := γ Yt(θ) = D • Y i t (θ), D • Y j t (θ) 1≤i,j≤m , (15)
where we have set γ t (θ) := γ Yt(θ) for notational sake in the computations below. We shall now compute γ -1 t (θ) as the solution to a SDE: Proposition 3.8. The matrix valued process γ -1 t (θ) is the unique solution to the following linear equation in η:

η t (θ) = α-1 0 (Y t (θ); θ) - d l=1 t 0 [η u (θ) αl (Y u (θ); θ) + αT l (Y u (θ); θ)η u ]dB l u - t 0 [η u (θ) β(Y u (θ); θ) + βT (Y u (θ); θ)η u (θ)]du, (16)
with

α0 (Y t (θ); θ) = m j=1 t 0 t 0 σ ij (Y r (θ); θ)σ i ′ j (Y r ′ (θ); θ) |r -r ′ | 2H-2 dr dr ′ , i, i ′ = 1, . . . , m
and where the other coefficients α and β are defined by

αl (Y u (θ); θ) = ∂ k σ i ′ l (Y u (θ); θ) 1≤i ′ ,k≤m and β(Y u (θ); θ) = ∂ k µ i ′ (Y u (θ); θ) 1≤i ′ ,k≤m . (17)
Proof. The proof of this fact is an adaptation of [START_REF] Hu | Differential equations driven by Hölder continuous functions of order greater than 1/2[END_REF]Theorem 7] to the case of a SDE with drift. We include it here for sake of completeness, and we drop the dependence of Y on θ for notational sake in the computations below.

Let us start by invoking Proposition 2.7 and equation ( 14) in order to compute the product of two first-order Malliavin derivatives:

D j r Y i t D j r ′ Y i ′ t = σ ij (Y r )σ i ′ j (Y r ′ ) + (18)
+ m k=1 t 0 d l=1 ∂ k σ il (Y u ) D j r ′ Y i ′ u D j r Y k u + ∂ k σ i ′ l (Y u ) D j r Y i u D j r ′ Y k u dB l u + t 0 ∂ k µ i (Y u ) D j r ′ Y i ′ u D j r Y k u + ∂ k µ i ′ (Y u ) D j r Y i u D j r ′ Y k u du .
Moreover, recall that γ t is defined by [START_REF] Hairer | Ergodicity of stochastic differential equations driven by fractional Brownian motion[END_REF]. Thus, the covariance matrix becomes

γ ii ′ t = d j=1 D j Y i t , D j Y i ′ t H = c H d j=1 t 0 t 0 D j r Y i t (θ) D j r ′ Y i ′ t (θ) |r -r ′ | 2H-2 dr dr ′ .
Plugging [START_REF] Hairer | Ergodicity of hypoelliptic SDEs driven by fractional Brownian motion[END_REF] into this relation, we end up with the following equation for γ ii ′ :

γ ii ′ t = αii ′ 0 + d l=1 t 0 m k=1 ∂ k σ i l (Y u ) γ i ′ k u + ∂ k σ i ′ l (Y u ) γ ik u dB l u + t 0 m k=1 ∂ k µ i (Y u ) γ i ′ k u + ∂ k µ i ′ (Y u ) γ i k u du.
Using our notation [START_REF] Hairer | Ergodicity theory of SDEs with extrinsic memory[END_REF] and matrix product rules, we obtain that γ t is solution to:

γ t = d l=1 t 0 ( αl (Y u )γ u + γ u αT l (Y u ))dB l u + t 0 ( β(Y u )γ u + γ u βT (Y u ))du.
Consider now η solution to [START_REF] Hairer | A simple framework to justify linear response theory[END_REF]. Applying again Proposition 2.7, it is readily checked that γ t η t = Id for any t ∈ [0, T ], which ends the proof.

Remark 3.9. Gathering equation ( 16) and Proposition 2.5, it is easily seen that for any t > 0 and θ ∈ Θ, Y t (θ) is a non degenerate random variable in the sense given at [32, Definition 2.1.2]: we have det(γ -1 t ) ∈ L p (Ω) for any p > 1.

Now that we have derived an equation for η = γ -1 , an equation for the Malliavin derivative of η is also available: Proposition 3.10. For any l ∈ {1, . . . , q} and n ≥ 1, the process η t = γ -1 t is n-time differentiable in the Malliavin calculus sense. Moreover, the process D i 1 ,...,in η t satisfies the following equation: for

t ≥ r 1 ∨ • • • ∨ r n , D i 1 ,...,in r 1 ,...,rn η ij t (θ) = - n k 1 =1 k 1 k 2 =1 (D i 1 ,...,i k 2 r 1 ,...,r k 2 α-1 0 D i 1 ,i 2 ,...,i k 1 -k 2 r 1 ,r 2 ,...,r k 1 -k 2 α0 D i 1 ,...,i n-k 1 r 1 ,...,r n-k 1 α-1 0 ) ij - d ℓ=1 t r 1 ∨...∨rn C ij ℓ,i 1 ,...,in (s; r 1 , . . . , r n ; θ)dB ℓ s - t r 1 ∨...∨rn
A ij i 1 ,...,in (s; r 1 , . . . , r n ; θ)ds,

where

A ij i 1 ,...,in (s; r 1 , . . . , r n ; θ) = m k 1 ,...,kν m k=1 [∂ ν k 1 ,...,kν ( β(Y u (θ); θ)) kj D i(I 1 ) r(I 1 ) η ik s (θ) . . . D i(Iν ) r(Iν ) η ik s (θ) D i(I 1 ) r(I 1 ) Y i s (θ) . . . D i(Iν ) r(Iν ) Y i s (θ)] +[∂ ν k 1 ,...,kν ( β(Y u (θ); θ)) ik D i(I 1 ) r(I 1 ) η kj s (θ) . . . D i(Iν ) r(Iν ) η kj s (θ) D i(I 1 ) r(I 1 ) Y j s (θ) . . . D i(Iν ) r(Iν ) Y j s (θ)]
and the same kind of equation holds for C ij l,i 1 ,...,in (s; r 1 , . . . , r n ; θ), with the coefficients β replaced by α l .

Proof. The proof of this proposition is based on Lemma 3.4 and the fact that

dA -1 λ dλ = -A -1 λ dA λ dλ A -1 λ .
Finally, one can also differentiate η with respect to our standing parameter θ, which yields: Lemma 3.11. The derivative of the inverse of the Malliavin matrix η t with respect to θ satisfies the following SDE

∇ l η t (θ) = ∇ l α-1 0 - d ℓ=1 t 0 {∇ l η u (θ) αℓ (Y u (θ); θ) + η u (θ)∇ l [ αℓ (Y u (θ); θ)] +∇ l [ αT ℓ (Y u (θ); θ)]η u (θ) + αT ℓ (Y u (θ); θ)∇ l η u (θ)}dB ℓ u - t 0 {∇ l η u (θ) β(Y u (θ); θ) +η u (θ)∇ l [ β(Y u (θ); θ)] + ∇ l [ βT (Y u (θ); θ)]η u (θ) + βT (Y u (θ); θ)∇ l η u (θ)}du, where ∇ l [ βℓ (Y u )] = ∂ βℓ (Y u )∇ l Y u + ∇ l βℓ (Y u ) and ∇ l [ αℓ (Y u )] = ∂ αℓ (Y u )∇ l Y u + ∇ l αℓ (Y u ).
3.2. Probabilistic representation of the likelihood. We have chosen to represent the log-likelihood of our sample thanks to the following formula borrowed from the stochastic analysis literature:

Proposition 3.12. Let F be a R m -valued non degenerate random variable (see Remark 3.9 for references on this concept), and let f be the density of F . For n ≥ 1 and (j 1 , . . . , j n ) ∈ {1, . . . , m} n , let H (j 1 ,...,jn) (F ) be defined recursively by H (j 1 ) (F ) = m j=1 δ((γ -1 F ) j 1 j DF j ) and

H (j 1 ,...,jn) (F ) = m j=1 δ γ -1 F jnj DF j H (j 1 ,...,j n-1 ) (F ) , (19)
where the Skorohod operator δ is defined at Section 2.2.3. Then one can write

f (x) = E 1 (F >x) H (1,...,m) (F ) = E (F -x) + H (1,...,m,1,...,m) (F ) , (20)
where 1 (F >x) := m i=1 1 (F i >x i ) and (Fx) + := m i=1 (F ix i ) + . Proof. The first formula is a direct application of [START_REF] Nualart | Malliavin Calculus and Related Topics[END_REF]Proposition 2.1.5]. The second one is obtained along the same lines, integrating by parts m additional times with respect to the first one.

The formula above can obviously be applied to Y t (θ) for any strictly positive t, since we have noticed at Remark 3.9 that Y t (θ) is a non-degenerate random variable. However, the expression of H (j 1 ,...,jn) (Y t (θ)) given by ( 19) is written in terms of Skorohod integrals, which are not amenable to numerical computations. We will thus recast this expression in terms of Young integrals plus some correction terms: Proposition 3.13. Under Hypothesis 2.3 and 3.1, let us define Q pji st := (γ -1 s ) pj D i s Y j t (θ) for 0 ≤ s < t ≤ T , p, j ∈ {1, . . . , m} and i ∈ {1, . . . , d}. Consider p ∈ {1, . . . , m} and a real valued random variable G which is smooth in the Malliavin calculus sense. Set

U p (G) = m i=1 d j=1 G t 0 Q pji st dB i s -c H m i=1 d j=1 t 0 t 0 D i s GQ pji rt |r -s| 2H-2 drds, (21)
where the integral with respect to B is understood in the Young sense. Then the quantities H (j 1 ,...,jn) (Y t (θ)) defined at Proposition 3.12 can be expressed as

H (j 1 ,...,jn) (Y t (θ)) = m j=1 U jn • • • • • U j 1 Y j t (θ) . (22)
Proof. It is an immediate consequence of Proposition 2.9, since we have noticed in our Remark 3.9 that Y t (θ) is a non-degenerate random variable.

The previous proposition is still not sufficient to warranty an effective computation of the log-likelihood. Indeed, the right hand side of (21) contains terms of the form D s [GQ pji rt ], which should be given in a more explicit form. This is the content of our next proposition.

Proposition 3.14. Set H (j 1 ,...,jn) (Y t (θ)) := K j 1 ...jn . Then the term D s [K j 1 ...jn Q pji rt ] in ( 21) can be computed inductively as follows:

(i) We have D s [K j 1 ...jn Q pji rt ] = D s K j 1 ...jn Q pji rt + K j 1 .
..jn D s Q pji rt , and D s Q pji rt is computed by invoking Proposition 3.10 for the derivative of γ -1 t and Lemma 3.4 for the derivative of Y t (θ). We are thus left with the computation of D s K j 1 ...jn .

(ii) Assume now that we can compute nr Malliavin derivatives of K j 1 ...jr . Notice that this condition is met for r = 0, since Y t (θ) itself can be differentiated n times in an explicit way according to Lemma 3.4 again. Then for any j 1 , . . . , j r+1 and k ≤ nr -1, the quantity K j 1 ...j r+1 can be differentiated k times, with a Malliavin derivative given by

D i 1 ,...,i k ρ 1 ...ρ k K j 1 ...j r+1 = k ℓ=1 D i 1 ,..., ǐℓ ...,i k ρ 1 ...ρ ℓ ...ρ k (K j 1 ...jr Q pji ρ ℓ t ) + d j=1 t 0 D i 1 ,...,i k ρ 1 ...ρ k (K j 1 ...jr Q pji st )dB j s -c H t 0 t 0 D k+1 r 1 ρ 1 ...ρ k (K j 1 ...jr Q pji r 2 t )|r 1 -r 2 | 2H-2 dr 1 dr 2 . ( 23
)
Proof. We focus on the induction step (ii), the other one being straightforward: for a smooth random variable W , one easily gets by induction that

D i 1 ,...,ip r 1 ...rp δ(W ) = p ℓ=1 D i 1 ,..., ǐℓ ...,ip r 1 ...ř ℓ ...rp W r ℓ + δ(D i 1 ,...,ip r 1 ...rp W ). ( 24
)
Suppose we know the nr Malliavin derivatives for

U jr • • • • • U j 1 (F ) := K j 1 ...jr . Recall moreover that K j 1 ...j r+1 = U j r+1 (K jr...j 1 ) = δ(K jr...j 1 Q •t )
Applying directly relation [START_REF] Kohatsu-Higa | Lower bounds for densities of uniformly elliptic random variables on Wiener space[END_REF] we thus get, for k ≤ m -1:

D i 1 ,...,i k ρ 1 ...ρ k δ(K jr...j 1 Q •t ) = k ℓ=1 D i 1 ,..., ǐℓ ,...,i k-1 ρ 1 ...ρ ℓ ...ρ k (K jr...j 1 Q ρ ℓ t ) + δ(D i 1 ,...,i k ρ 1 ...ρ k (K jr...j 1 Q •t )).
Our formula ( 23) is now obtained by applying Proposition 2.9 to the Skorohod integral δ(D k ρ 1 ...ρ k (K jr...j 1 Q •t )) above.

Example 3.15. As an illustration of the proposition above, we compute

U 2 • U 1 (F ) for F = Y i t , i ∈ {1, . . . , m} and our d-dimensional fBm B. Write first U 1 (Y i t ) = δ(Y i t (γ -1 ) 1j 1 D j 1 Y i t )
, and since this quantity has to be expressed in a suitable way for numerical approximations, we have

U 1 (Y i t ) = d j 1 =1 Y i t t 0 Q 1ij 1 ut dB j 1 u -c H d j 1 =1 t 0 t 0 D j 1 u 1 [Y i t Q 1ij 1 u 2 t ]|u 1 -u 2 | 2H-2 du 1 du 2 ,
where Q is defined at Proposition 3.13 and where the first integral in the right hand side is understood in the Young sense. In order to compute the second one, we have to compute Malliavin derivatives. This is done through Lemma 3.4 for Y and Proposition 3.10 for Q.

We now have to differentiate U 1 (Y i t ): the derivation rules for Skorohod integrals immediately yield

D j 2 u 2 [U 1 (Y i t )] = d j 2 =1 Y i t Q 2ij 2 u 2 t + d j 2 =1 δ(D j 2 u 2 Y i t Q 2ij 2 •t ).
Once again, the Skorohod integral above is not suitable for numerical approximations. Write thus

D j 2 u 2 [U 1 (Y i t )] = d j 2 =1 Y i t Q 2ij 2 u 2 t + d j 2 =1 t 0 D j 2 u 2 [Y i t Q 2ij 2 rt ]dB j 2 r -c H d j 2 =1 t 0 t 0 D j 2 u 2 D j 1 u 1 [Y i t Q 2ij 2 u 2 t ]|u 2 -u 1 | 2H-2 du 1 du 2 ,
and compute the Malliavin derivatives of the products Y Q thanks to Lemma 3.4 for Y and Proposition 3.10 for Q. Once this is done, just write

U 2 (U 1 (Y i t )) = δ(U 1 (Y i t )Q i •t ) = d j 2 =1 U 1 (Y i t ) t 0 Q 2ij 2 ut dB j 2 u -c H d j 2 =1 t 0 t 0 D j 2 u 2 [U 1 (Y i t )Q 2ij 2 u 1 t ]|u 2 -u 1 | 2H-2 du 1 du 2 .
In order to give our formula for the derivative of the log-likelihood, we still need to compute the derivative with respect to θ of H (j 1 ,...,jn) (Y t (θ)). For this we state the following lemma Lemma 3.16. The derivative with respect to θ of U p (Y t (θ)) can be written as

∇ l U p (Y i t (θ)) = d j=1 [∇ l Y i t (θ) t 0 Q pij st (θ) dB j s + Y i t (θ) t 0 ∇ l [Q pij st (θ)] dB j s ] -c H d j=1 t 0 t 0 ∇ l [D j s Y i t (θ)Q pij rt (θ)]|r -s| 2H-2 drds,
where ∇ l Y i t (θ) is computed according to Proposition 3.6 and ∇ l [D j s Y i t ] is given by Lemma 3.7. As far as ∇ l [Q pj st (θ)] is concerned, it is obtained through the following equation:

∇ l [Q pj st (θ)] = ∇ l η pj s (θ) D s Y j t (θ) + η pj s (θ)∇ l [D s Y j t (θ)],
where the expression for ∇ l η pj s (θ) is a consequence of Lemma 3.11.

We are now ready to state our probabilistic expression for the log-likelihood function [START_REF] Cheridito | Arbitrage in fractional Brownian motion models[END_REF].

Theorem 3.17. Assume Hypothesis 2.3 and 3.1 hold true. Let y t i , i = 1, . . . , n be the observation arriving at time t i . Let also Y t i be the solution to the SDE (1) at time t i . Then, the gradient of the log-likelihood function admits the following probabilistic representation:

∇ l ℓ n (θ) = n i=1 V i (θ) W i (θ) , with W i (θ) = E 1 (Yt i (θ)>yt i ) H (1,...,m) Y t i (θ) ( 25
)
and

V i (θ) = E ∇ l Y t i (θ) 1 (Yt i (θ)>yt i ) H (1,...,m,1,...,m) Y t i (θ) + Y t i (θ) -y t i + ∇ l H (1,...,m,1,...,m) Y t i (θ) , ( 26
)
where (i) H (j 1 ,...,jn) (Y t i (θ)) is given recursively by ( 22) and computed at Proposition 3.14

(ii) ∇ l Y t i (θ)
(θ) = n i=1 ln(f (t i , y t i ; θ)). Thus ∇ l ℓ n (θ) = n i=1 ∇ l f (t i , y t i ; θ) f (t i , y t i ; θ) := n i=1 V i (θ) W i (θ) .
Now W i (θ) can be expressed like ( 25) by a direct application of ( 20), first relation. As far as V i (θ) is concerned, write

f (t i , y t i ; θ) = E (Y t i (θ) -y t i ) + H (1,...,m,1,...,m) (Y t i (θ)) ,
according to the second relation in [START_REF] Hu | Fractional white noise calculus and applications to finance, Infnite dimentional analysis[END_REF]. By using standard arguments, one is allowed to differentiate this expression within the expectation, which directly yields [START_REF] Kushner | Stochastic approximation algorithms and applications[END_REF].

Discretization of the log-likelihood

The expression of the log-likelihood that we derived in Proposition 3.17 is a fraction of two expectations that do not have explicit formulas even in the one-dimensional case. In addition, our goal is to find the root of this non-explicit expression, the ML estimator, which is an even harder task. To solve this problem in practice we first use a stochastic approximation algorithm in order to find the root of ∇ l ℓ n (θ). In each iteration of the algorithm we compute the value of the expression using Monte-Carlo (MC) simulations. For each Monte-Carlo simulation, since we do not have available an exact way of simulating the kernels of the expectation, we use an Euler approximation scheme. More specifically, we simulate using Euler approximation terms such as Y t , DY t , which are solutions to fractional stochastic differential equations.

Therefore, in our approach we have three types of error in the computation of the MLE: the error of the stochastic approximation algorithm, the Monte-Carlo error and the discretization bias introduced by the Euler approximation for the stochastic differential equations. Our aim here is to combine the Monte Carlo and Euler approximations in an optimal way in order to get a global error bound for the computation of ∇ l ℓ n (θ).

4.1.

Pathwise convergence of the Euler scheme. The Euler scheme is the main source of error in our computations. There is always a trade-off between the number of Euler steps and the number of simulations, but what is usually computationally costly is the number of Euler steps. This is even worse when we deal with fractional SDEs, since the rate of convergence depends on H and the closer the value of H to 1/2, the more steps are required for the simulation.

In this section, we compute the magnitude of the discretization error we introduce. We measure the bias of the Euler scheme via the root mean square error. That is, we want to estimate the quantity sup

τ ∈[0,T ] (E|Y τ (θ) -Ȳ M τ (θ)| 2 ) 1/2
, where Y t (θ) is the solution to the SDE (1) and Ȳ M τ (θ) is the Euler approximation of Y τ (θ) given on the grid {τ k ; k ≤ M} by

Ȳ M τ k+1 (θ) = Ȳ M τ k (θ) + µ( Ȳ M τ k (θ); θ)(τ k+1 -τ k ) + d j=1 σ j ( Ȳ M τ k (θ); θ)δB M,j τ k τ k+1 , (27)
in which we denote δB M,j τ k τ k+1 = B M,j τ k+1 -B M,j τ k and τ k = kT M for k = 0, . . . , M -1. Notice that those estimates can be found in [START_REF] Deya | A Milstein-type scheme without Lévy area terms for SDEs driven by fractional Brownian motion[END_REF][START_REF] Friz | Multidimensional stochastic processes as rough paths[END_REF][START_REF] Mishura | The rate of convergence for Euler approximations of solutions of stochastic differential equations driven by fractional Brownian motion[END_REF]. We include their proof here because it is simple enough, and also because they can be easily generalized to the case of a linear equation. This latter case is of special interest for us, since it corresponds to Malliavin derivatives, and is not included in the aforementioned references. Notation 4.1. For simplicity, in this section we write Y := Y (θ). Proposition 4.2. Let T > 0 and recall that Ȳ M is defined by equation [START_REF] Lejay | An Introduction to Rough Paths[END_REF]. Then, there exists a random variable C with finite L p moments such that for all γ < H and H > 1/2 we have

Y t -Ȳ γ,T ≤ C T M 1-2γ (28)
Consequently, we obtain that the MSE is of order O(M 1-2γ ).

Proof. In order to prove [START_REF] León | Malliavin calculus for fractional delay equations[END_REF] we apply techniques of the classical numerical analysis for the flow of an ordinary differential equation driven by a smooth path. Namely, the exact flow of (1) is given by Φ(y; s, t) := Y t , where Y t is the unique solution of (1) when t ∈ [s, T ] and the initial condition is Y s = y. Introduce also the numerical flow

Ψ(y; τ k , τ k+1 ) := y + µ(y)(τ k+1 -τ k ) + d j=1 σ j (y)δB M,j τ k τ k+1 , (29)
where τ k = kT M , k = 0, . . . , M -1. Thus, we can write that

Ȳ M τ k+1 = Ψ Ȳ M τ k ; τ k , τ k+1 , k = 0, . . . , M -1 Y M 0 = α.
For q > k we also have that

Ψ(y; τ k , τ q ) := Ψ(•; τ q-1 , τ q ) • Ψ(•; τ q-2 , τ q-1 ) • . . . • Ψ(y; τ k , τ k+1 ).
The one-step error computes as

r k = Φ(y; τ k , τ k+1 ) -Ψ(y; τ k , τ k+1 ) = τ k+1 τ k µ(Y s ) -µ(y) ds + τ k+1 τ k σ(Y s ) -σ(y) dB s (30)
Furthermore, since Y ∈ C γ and B ∈ C γ for γ > 1/2, using [START_REF] Dasgupta | Arbitrage opportunities for a class of Gladyshev processes[END_REF] we have

τ k+1 τ k σ(Y s ) -σ(y) dB s ≤ c γ ∂σ ∞ Y γ B γ T M 2γ ≤ c γ,σ ∂σ ∞ B 1/γ γ B γ T M 2γ ,
where we used the fact that Y γ ≤ c σ B

1/γ γ (see Proposition 2.5). Similarly, for the drift part we have

τ k+1 τ k µ(Y s ) -µ(y) ds ≤ c γ ∂µ ∞ Y γ T M γ+1 ≤ c γ,µ ∂µ ∞ B 1/γ γ T M γ+1 .
Therefore, the one-step error [START_REF] Mishura | The rate of convergence for Euler approximations of solutions of stochastic differential equations driven by fractional Brownian motion[END_REF] satisfies

|r k | ≤ c µ,σ B 1+1/γ γ T M 2γ . (31)
Now, we can write the classical decomposition of the error in terms of the exact and numerical flow. Since Ȳ

M τ k = Φ( Ȳ M τ k ; τ k , τ k ) and Y τ k = Φ( Ȳ M τ 0 ; τ 0 , τ k ) we have Ȳ M τq -Y τq = Φ( Ȳτ 0 ; τ 0 , τ k ) -Φ( Ȳτq ; τ q , τ q ) = q-1 k=0 Φ( Ȳ M τ k ; τ k , τ q ) -Φ( Ȳτ k+1 ; τ k+1 , τ q ) . ( 32
) Since Φ Ȳ M τ k ; τ k , τ q = Φ Φ( Ȳ M τ k ; τ k , τ k+1 ); τ k+1 , τ q we obtain Φ( Ȳ M τ k ; τ k , τ q ) -Φ( Ȳ M τ k+1 ; τ k+1 , τ q ) = Φ Φ( Ȳ M τ k ; τ k , τ q ); τ k+1 , τ q -Φ( Ȳ M τ k+1 ; τ k+1 , τ q ) ≤ C T ( B γ ) |Φ( Ȳ M τ k ; τ k , τ k+1 ) -Ȳ M τ k+1 |,
where we have used the fact that

|Φ(α; t, s) -Φ(β; t, s) ≤ C T ( B γ )|α -β|,
where C T is a subexponential function (see Proposition 2.5 again). Moreover, owing to relation [START_REF] Neuenkirch | Trees and asymptotic developments for fractional diffusion processes[END_REF],

|Φ( Ȳ M τ k ; τ k , τ q ) -Ȳ M τ k+1 | = |r k | ≤ c µ,σ B 1+1/γ γ T M 2γ . (33)
Therefore, replacing [START_REF] Nualart | Differential equations driven by fractional Brownian motion[END_REF] in [START_REF] Nualart | Malliavin Calculus and Related Topics[END_REF] for any q ≤ n we obtain

| Ȳ M τq -Y τq | ≤ c µ,σ B 1+1/γ γ q-1 k=0 T M 2γ
Let us push forward this analysis to Hölder type norms on the grid 0 ≤ τ 1 < . . . < τ n = T . We have for

q ≥ p δ Y -Ȳ M τpτq = Φ(Y τp ; τ p , τ q ) -Y τp -Ψ( Ȳ M τp ; τ p , τ q ) -Ȳ n τp = Φ(Y τp ; τ p , τ q ) -Y τp -Φ( Ȳ M τp ; τ p , τ q ) -Ȳ M τp -Ψ( Ȳ M τp ; τ p , τ q ) -Φ( Ȳ M τp ; τ p , τ q ) = Φ(Y τp ; τ p , τ q ) -Φ( Ȳ M τp ; τ p , τ q ) -Y τp -Ȳ M τp -Ψ( Ȳ M τp ; τ p , τ q ) -Φ( Ȳ M τp ; τ p , τ q ) .
Similar to the calculations leading to [START_REF] Nualart | Differential equations driven by fractional Brownian motion[END_REF] we obtain

Ψ( Ȳ M τp ; τ p , τ q ) -Φ( Ȳ M τp ; τ p , τ q ) ≤ c µ,σ B 1+1/γ γ q-1 k=p T M 2γ .
Moreover, owing to Proposition 2.5 part (2), observe that

Φ(Y τp ; τ p , τ q ) -Φ( Ȳ M τp ; τ p , τ q ) -Y τp -Ȳ M τp |τ q -τ p | γ ≤ c( B γ ) |Y τp -Ȳ M τp |.
Consequently, we have that for 0

≤ p < q ≤ M δ Y -Ȳ M τpτq ≤ c ′ ( B 1+1/γ γ ) q-1 k=p T M 2γ + |τ q -τ p | γ q k=0 T M 2γ which easily yields that sup p,q=0,1,...,M -1,p =q δ Y -Ȳ M τpτq |τ p -τ q | γ ≤ c( B γ ) M 1-2 γ .
By "lifting" this error estimate to [0, T ] and since |t -

s| ≤ T /M, Y t -Ȳ γ,∞,T ≤ C M 1-2 γ , (34)
which concludes the first part of the proof.

Regarding the order of the Mean Square Error, it suffices to note that the constant C has finite L p moments.

As mentioned before, an elaboration of Proposition 4.2 is needed in the sequel. Indeed, in the expression of the log-likelihood in Proposition 3.17 we need to discretize more complicated quantities of the underlying process, such as [START_REF] Gubinelli | Controlling rough paths[END_REF] or [START_REF] Hairer | A simple framework to justify linear response theory[END_REF]. To this aim, let us notice first that all those equations can be written under the following generic form:

Z t = α + t 0 ξ 2 u Z u du + t 0 ξ 1,j u Z t dB j u , (35)
where ξ 1 , ξ 2 are stochastic processes with bounded moments of any order. The corresponding Euler discretization is

ZM τ k = ZM τ k + ξ 2 τ k ZM τ k (τ k+1 -τ k ) + d j=1 ξ 1,j τ k Zτ k δB j,M τ k τ k+1 , (36)
and we give first an approximation result in this general context: Proposition 4.3. Let T > 0, and consider the R q -valued solution Z to equation [START_REF] Odde | Stochastic dynamics of the nerve growth cone and its microtubules during neurite outgrowth[END_REF], where α ∈ R q , ξ 2 , ξ 1,j ∈ R q,q and we suppose that ξ 2 γ and ξ 1,j γ belong to L p (Ω) for any value of p ≥ 1. Let ZM be defined by equation [START_REF] Pedersen | Consistency and asymptotic normality of an approximate maximum likelihood estimator for discretely observed diffusion processes[END_REF]. Then, there exists a random variable C

′ with L p finite moments, such that for all γ < H and H > 1/2 we have

Z -Z γ,T ≤ C ′ T M 1-2γ (37)
Consequently, we obtain that the Mean Square Error is of order O(M 1-2γ ).

Proof. We follow a similar approach as in the previous proposition. Thus, the exact flow is equal to Φ(ζ; s, t) := Z t , where Z t is the unique solution of equation ( 35) when t ∈ [s, T ] and the initial condition is Z s = ζ. Consider also the numerical flow

Ψ(ζ; τ k , τ k+1 ) := ζ + ξ 2 u ζ(τ k+1 -τ k ) + d j=1 ξ 1,j u ζδB j,M τ k τ k+1 ,
where τ k = kT /M, n = 0, . . . , M -1. Thus, we have

ZM τ k+1 = Ψ( ZM τ k ; τ k+1 , τ k ), k = 0, . . . , M -1 ZM 0 = α.
In this case, the one-step error can be written as

r k = Φ(ζ; τ k , τ k+1 ) -Ψ(ζ; τ k , τ k+1 ) = τ k+1 τ k ξ 2 u (Z s -ζ)du + τ k+1 τ k ξ 1 u (Z s -ζ)dB u
We now treat each term separately. Therefore, using the fact that Z γ ≤ exp(c B 1/γ γ ), which is recalled at Proposition 2.5 point (4) in a slightly different context, we have that

τ k+1 τ k ξ 1 s (Z s -ζ)dB s ≤ c γ Zξ 1 γ B γ T M 2γ ≤ c γ exp( B 1/γ γ ) B γ T M 2γ .
Similarly, we also have

τ k+1 τ k ξ 2 s (Z s -ζ)ds ≤ c γ Zξ 2 γ B γ T M 2γ ≤ c γ exp( B 1/γ γ ) B γ T M 2γ .
Therefore, the one-step error satisfies the following inequality

|r k | ≤ c γ exp( B 1/γ γ ) B γ T M 2γ .
Along the same lines as for Proposition 4.2, the decomposition of the error in terms of the exact and numerical flow becomes ZM τq -Z τq = Φ( ZM τq ; τ q , τ q ) -Φ( ZM

τ 0 ; τ 0 , τ k ) = q-1 k=0 Φ( ZM τ k+1 ; τ k+1 , τ q ) -Φ( ZM τ k ; τ k , τ q ) ,
and the same inequalities allowing to go from [START_REF] Nualart | Malliavin Calculus and Related Topics[END_REF] to [START_REF] Nualart | Differential equations driven by fractional Brownian motion[END_REF] yield

| Zτq -Z τq | ≤ c γ T M 2γ .
The claim of the proposition follows now as in Proposition 4.2.

We now use the previous proposition in order to approximate the kernels of the expectations in ∇ l ℓ n (θ). Let us first introduce the following notation: Notation 4.4. Let W i (θ), V i (θ) as in ( 25) and ( 26) respectively and define w i (θ) and v i (θ) as

w i (θ) = 1 (Yt i (θ)>yt i ) H (1,...,m) Y t i (θ) (38)
v i (θ) = ∇ l Y t i (θ) 1 (Yt i (θ)>yt i ) H (1,.,m,1,.,m) + Y t i (θ) -y t i + ∇ l H (1,.,m,1,.,m) . (39)
Let also wM i and vM i to be the Euler discretized versions of ( 38) and (39) using 4.3, and set

W M i (θ) = E[ wM i ] and V M i (θ) = E[v M i ].
Our convergence result for ∇ l ℓ n (θ) can be read as follows:

Theorem 4.5. Recall from Theorem 3.17 that ∇ l ℓ n (θ) can be decomposed as (θ) . Then the following approximation result holds true:

∇ l ℓ n (θ) = n i=1 V i (θ) W i
V i (θ) -V M i (θ) + W i (θ) -W M i (θ) ≤ c M 2γ-1 ,
for a strictly positive constant c.

Proof. We focus on the bound for |V i (θ) -V M i (θ)|, the other one being very similar. Now, applying Proposition 4.3 to the particular case of the equations governing Malliavin derivatives, we easily get

v t -v γ,T ≤ C 2 M 1-2γ ,
for an integrable random variable C 2 . The proof is now easily finished by invoking the inequality

V i (θ) -V M i (θ) ≤ E [ v t -v γ,T ] .
Remark 4.6. We have given two separate approximations for V i (θ) and W i (θ). In order to fully estimate (V i (θ)/W i (θ)) -( V M i (θ)/ W M i (θ)), one should also prove that W i (θ) is bounded away from 0. This requires a lower bound for densities of differential equations driven by fractional Broawnian motion, which are out of the scope of the current article. 4.2. Efficiency of the Monte Carlo simulation. In this section we aim to study the computational tradeoff between the length of a time period in the Euler discretization (i.e. 1/M) and the number of Monte Carlo simulations of the sample path (i.e. N). In order to do so we consider wM i and vM i as above. Recall that, given t units of computer time, the Monte-Carlo estimators for W i (θ) and V i (θ) can be written as

1 c 1 (t, 1 M ) c 1 (t, 1 M ) k=1 w M i,k , 1 c 2 (t, 1 M ) c 2 (t, 1 M ) k=1 v M i,k
where {w M i,ℓ ; ℓ ≥ 1} (resp. {v M i,ℓ ; ℓ ≥ 1}) is a sequence of i.i.d. copies of w M i (resp. of v M i ), and c 1 (t, 1 M ), c 2 (t, 1 M ) are the maximal number of runs one is allowed to consider with t units of computer time. Using the result by [START_REF] Durham | Numerical techniques for maximum likelihood estimation of continuous-time diffusion processes[END_REF] we can state the following proposition: Proposition 4.7. Let N be the number of Monte Carlo simulations and M the number of steps of the Euler scheme, then the tradeoff between N and M for computing W i (θ) (and similarly

V i (θ)) is N ≍ M γ 2γ-1 -3
, for all 1/2 < γ < H and γ = T m(d + 1), where T is the time horizon, m the dimension of the observed process and d the dimension of the noise process.

Proof. We discuss the proof only for W i , by following exactly the same steps we can obtain the same result for V i .

We only need to check that our process w satisfies the conditions of Theorem 1 in [START_REF] Durham | Numerical techniques for maximum likelihood estimation of continuous-time diffusion processes[END_REF].

(i) We can easily see that the discretized wM t i converges uniformly to w t i . (ii) In addition, we have bounded moments of w t i , thus E[ W 2

t i ] → E[w 2 t i ]. (iii) From Proposition 4.5 we have that the rate of convergence of the Euler scheme of wM t i is M 1-2γ , for 1/2 < γ < H. (iv) The computer time required to generate wM t i is given by τ (1/M), which satisfies:

τ (1/M) = T m(d + 1)M = γM
where T is the length of the time period, m is the dimension of the SDE, d is the dimension of the fBm and M is the number of Euler steps. By applying Theorem 1 (by [START_REF] Durham | Numerical techniques for maximum likelihood estimation of continuous-time diffusion processes[END_REF]) the optimal rule for choosing the number of Monte-Carlo simulations and the number of Euler steps is chosen such that the asymptotic error is minimized. Therefore, for t the total budget of computer time, as t increases, then the Euler step should converge to zero with order 1-2γ γ+2-4γ or equivalently:

1 M ≍ t 1-2γ γ+2-4γ thus t ≍ M -γ+2-4γ 1-2γ .
But the number of operations needed for an arbitrary Monte Carlo simulation t 0 is equal to γMN. Thus, we finally obtain that N ≍ M -γ+2-4γ 1-2γ -1 .

∇l ℓ n (θ) = Vi Ŵi := 1 N N k=1 vM i,k 1 N N k=1 wM i,k , (40)
where wM 1,k , wM 2,k , . . . and vM 1,k , vM 2,k , . . . are iid copies of wM i and vM i respectively. Our aim in this section is to give a global bound for the mean square error obtained by approximating ∇ l ℓ n (θ) by ∇l ℓ n (θ). Proof. We discuss the idea of the proof for the W i term first:

E Ŵi -W i 2 = E 1 N N k=1 wM i,k -E[w i (θ)] 2 = E 1 N N k=1 wM i,k - 1 N N k=1 w i,k + 1 N N k=1 w i,k -E[w i (θ)] 2 .
Thanks now to the independence property between Monte Carlo runs, we get

E Ŵi -W i 2 ≤ 2 N N k=1 E( wM i,k -w i,k ) 2 + 2 E 1 N N k=1 w i,k -E[w i (θ)] 2 = 1 N N k=1 (Euler MSE) 2 + (Monte Carlo MSE) 2 ≍ (M 1-2γ ) 2 + 1 N ,
and thus

MSE Ŵi -W i ≍ (M 1-2γ ) 2 + 1 N .
Now, if we use Proposition 4.7, i.e. N ≍ M -γ+2-4γ 1-2γ -1 , for all 1/2 < γ < H, and γ = T m(d + 1), where T is the time horizon, m the dimension of the observed process and d the dimension of the noise process, we have

MSE Ŵi -W i ≍ M 2-4γ + M γ 1-2γ +3 ≍ M 1-2γ ,
since the first is the dominant term above. Following the same procedure, we can show that MSE( Vi -V i ) ≍ M 1-2γ and thus the claim of the proposition follows easily. Remark 4.9. In Proposition 4.8 the rate of convergence is independent of the dimension of the problem, i.e. it is independent of the parameter γ = T m(d + 1).

Numerical Examples

In this section our aim is to investigate the performance of the suggested maximum likelihood method in practice. We study the one-dimensional fractional Ornstein-Uhlenbeck process, a linear two-dimensional system of fractional SDEs and then some real data given by a financial time series. Before presenting our results, we first discuss some technical issues raised by the algorithmic implementation of our method.

The goal is to find the root of the quantity ∇ l ℓ n (θ) with respect to θ. We can divide this procedure in two parts. The first part consists in computing the root of the log-likelihood using a stochastic approximation algorithm. This is a stochastic optimization technique firstly introduced by Robbins and Monro (1951) that is used when only noisy observations of the function are available. In our case it is appropriate, since we want to solve

∇ l ℓ n (θ) = 0,
where ∇ l ℓ n (θ) is given by Theorem 3.17 and has to be approximated by ∇l ℓ n (θ). Thus, the recursive procedure is of the following form θk+1 = θka k ∇l ℓ n ( θk ). [START_REF] Tudor | Statistical aspects of the fractional stochastic calculus[END_REF] where ∇l ℓ n is the estimate of ∇ l ℓ n at the k-th iteration based on the observations and a k is a sequence of real numbers such that ∞ k=1 a k = ∞ and ∞ k=1 a 2 k < ∞. Under appropriate conditions (see for example [START_REF] Blum | Multidimensional stochastic approximation methods[END_REF]), the iteration in (41) converges to θ almost surely. The step sizes satisfy a k > 0 and the way that we choose them can be found in [START_REF] Kushner | Stochastic approximation algorithms and applications[END_REF].

The second part consists of the computation of ∇l ℓ n ( θk ) at each step of the stochastic approximation algorithm. Thus, for a given value of θ k (the one computed at the k-th iteration) we want to compute ∇l ℓ n (θ k ) when we are given n discrete observations of the process: y t i , i = 1, . . . , n. Here, we describe the main idea of the algorithm we use for only one step. Thus, assume that we are at [t i-1 , t i ], and at time t i we obtain the i-th observation. We want to compute W i (θ) and V i (θ) according to expressions [START_REF] Kou | Generalized Langevin equation with fractional Gaussian noise: subdiffusion within a single protein molecule[END_REF] and ( 26) respectively. To compute the expectations we use simple Monte-Carlo simulations.Therefore, we discretize the time interval into N steps

t i-1 = s 0 < s 1 < • • • < s N = t i .
From each simulated path (apart from that of fBm) we only need to keep the terminal value which is the value of the process at time t i . The algorithm is the following (1) Simulate N values of fBm in the interval [t i-1 , t i ] using for example the circulant matrix method (any exact -preferably-simulation technique can be used). (2) Using the simulated values from step 1 and an Euler scheme for the SDE (1), simulate the value of the process at time t i . For example, for k = 0, . . . , N 8) Using the previous steps, numerical integration for the double integral and numerical integration for the stochastic integral we compute U p (Y t i (θ)) as defined in Proposition 3.12. (9) Recursively compute H (1,...,m) (Y t i (θ)) as given in [START_REF] Hu | Differential equations driven by Hölder continuous functions of order greater than 1/2[END_REF]. [START_REF] Durham | Numerical techniques for maximum likelihood estimation of continuous-time diffusion processes[END_REF] Combine steps 3 and 9 to obtain the kernel W i (θ). [START_REF] Friz | Multidimensional stochastic processes as rough paths[END_REF] We repeat steps 1 through 10 N times and we average these values to obtain estimate for the expectation W i (θ). Using a similar procedure we can obtain an estimate for the expectation V i (θ). Finally, for each i = 1, . . . , n we compute V i (θ)/W i (θ) and sum over i to obtain the desired value of the log-likelihood at θ k .

Ȳ M s k = Ȳ M s k-1 + µ( Ȳ M s k )(s k -s k-1 ) + d j=1 σ (j) ( Ȳ M s k-1 )(B (j) s k -B (j) s k-1 ). ( 3
We have completed the study of our numerical approximation of the log-likelihood, and are now ready for the analysis of some numerical examples. 5.1. Fractional Ornstein-Uhlenbeck process. Though our method can be applied to highly nonlinear contexts, we focus here on some linear situations, which allow easier comparisons with existing methods or exact computations. Let us first study the onedimensional fractional Ornstein-Uhlenbeck process, i.e.

dY t = -λY t dt + dB t , (42)
where the solution is given Y t (λ) = t 0 e -λ(t-s) dB s (notice the existence of an explicit solution here). In this case our methodology is quite simplified. The log-likelihood can be written as follows:

∂ λ ℓ(λ; y) = n i=1 E ∂ λ Y t (λ) 1 (Yt(λ)>y) H (1,1) (λ) + Y t (λ) -y + ∂ λ H (1,1) (λ) E 1 (Yt(λ)>y) H (1) Y t (λ), 1 . The Malliavin derivative of Y t (λ) satisfies the following ODE D s Y t (λ) = 1 -λ t s D s Y u (λ)du, with solution D s Y t (λ) = e -λ t 1 {s≤t} . The corresponding norm is D • Y t (λ) 2 = c H t s t s e -λ(u+v) |u -v| 2H-2 dudv.
The higher order derivatives of Y t (λ) are equal to zero. Therefore, Remark 5.1. The value of H used for the simulation of the process is 0.6. The number of observations is n = 50, the number of Euler steps is M = 500, the number of stochastic approximation steps is K = 50 and the number of MC simulations N = 500.

H (1) Y t (λ) = 1 D • Y t (λ) 2
∂ λ H (1,1) (λ) = 1 D • Y t (λ) 8 D • Y t (λ)
5.2. Two-dimensional fractional SDE. In this section we study the following system of fractional OU processes:

dY (1) t = -αY (2)
t dt + βdB

(1) t dY

(2) t

= -βY

(1)

t dt + βdB (2)
t .

In this case, the computations are more involved even though the SDEs are linear functions of Y . Furthermore, the parameter we want to estimate is two-dimensional as well (θ = (α, β) T ), which complicated the optimization procedure. Therefore, instead of computing only one derivative, we need to compute both derivatives with respect to α and β and then compute the solution of the system of two equations The covariance matrix γ t is given by ( D • Y i t , D • Y j t ) 1≤i,j≤2 . The inverse of the covariance matrix satisfies the following SDE

γ -1 t = - t 0 [γ -1 u M + M T γ -1 u ]du,
where M = 0 α β 0 Now, it remains to compute the quantities H (1,2) and H (1,2,1,2) . This can be done using the recursive formulas in Proposition 3.12, but we need to keep in mind that higher order derivatives of Y are equal to zero, thus they will be simplified. Indeed, Of course, recall that Q pj st = (γ -1 s ) pj D s Y j t . In practice, these quantities are computed recursively. The last step is to compute the derivative of H (1,2,1,2) (Y t ) with respect to α and β, which in this case is not as complicated and compute the MLEs using the algorithm discussed in the previous section. The table below summarizes our results, and we have plotted the corresponding histograms in Figure 1 Remark 5.2. The value of H used for the simulation of the process is 0.6. The number of observations is n = 50, the number of Euler steps is N = 500, the number of stochastic approximation steps is K = 50 and the number of MC simulations M = 500.

H (1) (Y t ) =
5.3. Application to financial data. One of the most popular applications of fractional SDEs is in finance. Hu and Oksendal, [START_REF] Hu | Fractional white noise calculus and applications to finance, Infnite dimentional analysis[END_REF], introduced the fractional Black-Scholes model in order to account for inconsistencies of the existing models in practice. More specifically, the stock price is described therein by a fractional geometric Brownian motion with Hurst parameter 1/2 < H < 1. The choice of this model is based on empirical studies that displayed the presence of long-range dependence on stock prices, for example in Willinger, Taqqu and Teverovsky, [START_REF] Willinger | Stock market prices and long-range dependence[END_REF]. However, the presence of fractional Brownian motion in the model allows for arbitrage in the general setting. It has been shown that arbitrage opportunities can be avoided in a number of ways, for example the reader can refer to Rogers [START_REF] Rogers | Arbitrage with fractional Brownian motion[END_REF], Dasgupta and Kallianpur [START_REF] Dasgupta | Arbitrage opportunities for a class of Gladyshev processes[END_REF] and Cheridito [START_REF] Cheridito | Arbitrage in fractional Brownian motion models[END_REF]. We choose to model the stock price as as follows:

Histogram of alpha

dS t = µS t dt + σdB t , ( 44
)
where B is a fractional Brownian motion with Hurst index 1/2 < H < 1. For this SDE (as well as for a more general class of fractional SDEs) Guasoni, [START_REF] Guasoni | No arbitrage under transaction costs, with fractional Brownian motion and beyond[END_REF], proved that there is no arbitrage when transaction costs are present. Our goal is to estimate the unknown parameters µ and σ based on daily observations of the S&P 500 index (data from June 2010 until December 2010). Since the Hurst parameter is piece-wise constant, we devide the data in three groups (of 50 daily observations each) and we compute for each one the Hurst index using the Rescaled-Range (R/S) statistic. We obtain that for the first group of data Ĥ1 = 0.59, for the second Ĥ2 = 0.63 and for the third one Ĥ3 = 0.61. For the different groups, we apply our maximum likelihood approach in order to estimate µ and σ. The estimates are summarized in the following table :  Estim. Parameters Group 1: Ĥ1 = 0.59 Group 2: Ĥ2 = 0.63 Group 3: Ĥ3 = 0.61 μ 0.015 (0.0123) 0.019 (0.0144) 0.011 (0.0214) σ 0.352 (0.058) 0.339 (0.046) 0.341 (0.024)

Remark 5.3. The volatility is computed in years. In addition, during this period of time the historical volatility is around 0.38, which is coherent with our own estimation.

Hypothesis 2 . 3 . 2 l=0

 232 For any θ ∈ Θ, we assume that µ(•; θ) : R m → R m and σ(•; θ) : R m → R m,d are C 2 b coefficients. Furthermore, we have sup θ∈Θ 1≤i 1 ,...,i l ≤q

4. 3 .

 3 Discretization of the score function. Consider the following discretized version of the score function, i.e. ∇ l ℓ n (θ):

Proposition 4 . 8 .

 48 The discretized score function ∇l ℓ n (θ) converges to the continuous score function ∇ l ℓ n (θ) with rate of convergence of order M -(2γ-1) , where 1/2 < γ < H and M is the number of Euler steps used in the discretization.

e

  -λ(u+v) dB u dB vc H D • Y t (λ) -2 .The derivative with respect to the unknown parameter λ satisfies∂ λ Y t (λ) = -t 0 Y s (λ) -λ∂ λ Y s (λ)ds with solution ∂ λ Y t (λ) = t 0 (ts)e -λ(t-s) dB s .The last term we need to compute is:

∇[E[ 1 ( 1 ×

 11 α ℓ(α, β; y) = 0, ∇ β ℓ(α, β; y) = 0, where∇ l ℓ(α, β; y) = n i=1 Yt(α,β)>y) H (1,2) (Y t (α, β))] -{E[∇ l Y t (α, β) 1 (Yt(α,β)>y) H (1,2,1,2) (α, β) + (Y t (α, β)y) + ∇ l H (1,2,1,2) (α, β)]}and l = α or β. The Malliavin derivative of Y t computes as follows:

Figure 1 .

 1 Figure 1. Empirical Distribution of the estimators for α and β.

  is given by Proposition 3.6 (iii) ∇ l H (1,...,m,1,...,m) is obtained by applying Lemma 3.16.Proof. Recall that under Hypothesis 2.3 and 3.1, Y t (θ) admits a C ∞ density f (t, •; θ) for any t > 0 and θ ∈ Θ. Moreover, we have defined ℓ n (θ) as ℓ n

  ) Using step 2 and the observation at time t i , compute the indicator function 1 (Yt i (θ)>yt i ) . (4) Using step 1 and an Euler scheme simulate D t i Y i τ , as given in Lemma 3.3 for n = 1 -first Malliavin derivative-. (5) Using step 1 and an Euler scheme simulate η kj t i , k, j = 1, . . . , m, as given in Proposition 3.7. (6) Steps 4 and 5 are used to compute Q pj st i , p ∈ {1, . . . , m}, j ∈ {1, . . . , d} as defined in Propositions 3.12 and 3.14. (7) Simulate the Malliavin derivative of the product D s [Y t Q pj rt ]. (

  + v)e -λ(u+v) dB u dB v -2c H D • Y t (λ) 2 + v)e -λ(u+v) |u -v| 2H-2 dudv + v)e -λ(u+v) |u -v| 2H-2 dudv D • Y t (λ) 4 .Now, we compute the MLE following the algorithm we described above. The results we obtained are summarized in the following table:

			t	t
			4 -(u t s r t
	s -(u -r c 2 H t s t r -(u True λ MLE λ Standard Error
	0.5	0.497	0.00369
	4	3.861	0.00127

  .

	Parameter True Value MLE Standard Error
	α	2	2.003	0.0518
	β	4	3.987	0.0157

S. Tindel is partially supported by the ANR grant ECRU. Both authors are part of the BIGS (Biology, Genetics and Statistics) team at INRIA Nancy.