Marie Bernhart
email: marie-externe.bernhart@edf.fr

Huyên Pham
email: pham@math.jussieu.fr

Peter Tankov
email: peter.tankov@polytechnique.org

Xavier Warin
email: xavier.warin@edf.fr

Swing Options Valuation: a BSDE with Constrained Jumps Approach

Keywords: Backward stochastic differential equations with constrained jumps, Impulse control problems, Swing options, Monte Carlo methods

We introduce a new probabilistic method for solving a class of impulse control problems based on their representations as Backward Stochastic Differential Equations (BSDEs for short) with constrained jumps. As an example, our method is used for pricing Swing options. We deal with the jump constraint by a penalization procedure and apply a discrete-time backward scheme to the resulting penalized BSDE with jumps. We study the convergence of this numerical method, with respect to the main approximation parameters: the jump intensity λ, the penalization parameter p > 0 and the time step. In particular, we obtain a convergence rate of the error due to penalization of order (λp) α-1 2 , ∀α ∈ 0, 1 2 . Combining this approach with Monte Carlo techniques, we then work out the valuation problem of (normalized) Swing options in the Black and Scholes framework. We present numerical tests and compare our results with a classical iteration method.

Introduction

In this report, we introduce a new probabilistic method for solving impulse control problems based on their representations as Backward Stochastic Differential Equations (BSDEs for short) with constrained jumps. As an example, our method is used for pricing Swing options in the Black and Scholes framework.

BSDEs provide alternative characterizations of the solution to multiple-obstacle, optimal switching (see among others [START_REF] Hamadène | On the Starting and Stopping Problem: Application in reversible investments[END_REF][START_REF] Carmona | Pricing Asset Scheduling Flexibility Using Optimal Switching[END_REF][START_REF] Hu | Multi-dimensional BSDE with oblique reflection and optimal switching[END_REF][START_REF] Porchet | Problems of Valuation and Organization in Energy Markets[END_REF][START_REF] Djehiche | Stochastic Impulse Control of Non-Markovian Process[END_REF]) and more generally impulse control problems: Kharroubi et al. [START_REF] Kharroubi | Backward SDEs with constrained jumps and Quasi-Variational Inequalities[END_REF] recently introduced a family of BSDEs with constrained jumps providing a representation of the solution to such problems. A challenging question is that of the numerical approximation of this kind of BSDEs with constrained jumps.

A discrete-time backward scheme for solving BSDEs with jumps (without constraint) has been introduced by Bouchard and Elie [START_REF] Bouchard | Discrete-time approximation of decoupled forwardbackward SDE with jumps[END_REF]. In our case, the main difficulty comes from the constraint, which concerns the jump component of the solution. These BSDEs do not a priori involve any Skorohod type minimality condition. In consequence, classical approaches by projected schemes (discretely reflected backward schemes) used for example by [START_REF] Bouchard | Discrete-time approximation for continuously and discretely reflected BSDEs[END_REF] and [START_REF] Chassagneux | Discrete-time Approximation of Multidimensional BSDEs with oblique reflections[END_REF] cannot be used.

We consider a penalization procedure to deal with the constraint on jumps and provide a convergence rate of the penalized solution to the exact solution. This allows us to establish a convergence rate of the error between the solution of the considered impulse control problem and the numerical approximation given by the discrete-time solution to the penalized BSDE with jumps, as the penalization coefficient and the number of time steps go to infinity.

The rest of the report is structured as follows: in Section 2, we set the considered impulse control problem in the mathematical framework of BSDEs with constrained jumps. We present in Section 3 our penalization approach and provide a global convergence rate of our approximation. In Section 4, our method is used for pricing multi-exercise options, so-called (normalized) Swing options. This multiple optimal stopping time problem leads to a particularly degenerate three-dimensional impulse control problem. We combine our BSDE-based approach with Monte Carlo techniques and deal with Swing options with a small maximal number of exercises rights, due to large computational times. We compare our pricing results with those obtained by a classical iteration-based approach proposed for example by [START_REF] Carmona | Optimal multiple stopping and valuation of Swing options[END_REF].

BSDE Representation for Impulse Control Problems

Let T be a given time horizon. We work in a complete probability space (Ω, F , P), on which is defined a d-dimensional Brownian motion W and a Poisson process N with intensity λ > 0. We denote by F = (F t ) t≥0 , the augmentation of the natural filtration generated by W and N , by F W = (F W t ) t≥0 the one generated by W , and by P, the σ-algebra of predictable sub-sets of Ω × [0, T ].

Notation

Throughout this report, the euclidean norm defined on R d or on R will be indiscriminately denoted by | • |. The matrix transposition is denoted by ⊥. In addition, unless specified otherwise, C will denote a strictly positive constant depending only on Lipschitz constants of the coefficients of the problem, see assumptions Besides, we shall use the standard notations:

• S 2 , the set of real-valued càdlàg adapted processes Y = (Y t ) 0≤t≤T such that

Y S 2 := E sup 0≤t≤T |Y t | 2 1 2
< ∞ .

• A 2 , the sub-set of S 2 such that

A 2 := K ∈ S 2 : (K t ) 0≤t≤T nondecreasing , K 0 = 0 . • L 2 F ([0, T ]
), the set of real-valued adapted processes (φ t ) 0≤t≤T such that

E T 0 |φ t | 2 dt < ∞ .
• L 2 (W ), the set of real-valued P-measurable processes Z = (Z t ) s≤t≤r such that

Z L 2 (W ) := E T 0 |Z t | 2 dt 1 2
< ∞ .

• L 2 (N ), the set of real-valued P-measurable processes V = (V t ) s≤t≤r such that

V L 2 (N ) := E T 0 |V t | 2 λdt 1 2
< ∞ .

• V denotes the set of P-measurable essentially bounded processes, valued in (0, ∞) and V p = {ν p ∈ V : ν p t ≤ p a.s.}.

A Class of Impulse Control Problems

We consider the class of impulse control problems whose value function is defined by:

v(t, x) = sup u=(τ k ) k≥1 ∈U (t,T ] E   g(X t,x,u T ) + T t f (X t,x,u s )ds + k≥1 t<τ k ≤T κ(X t,x,u τ - k )    . ( 1
)
An impulse strategy u = (τ k ) k≥1 is said to be admissible for problem (1) (and belongs to U (t,T ] ) if it is a non-decreasing sequence of F W -stopping times valued in (t, T ] (we set by convention τ 0 = t) such that, if

n u (t,T ] := ♯ {k ≥ 1 : t < τ k ≤ T }
denotes the (random) number of interventions of the strategy u before time T , then

E n u (t,T ] 2 < C, (2)
for some universal constant C > 0. The controlled state variable X t,x,u is a càdlàg process with dynamics

X t,x,u s = x + s t b(X t,x,u r )dr + s t σ(X t,x,u r )dW r + t<τ k ≤s γ(X t,x,u τ - k ), ∀s ≥ t. (3)
Between two successive intervention times τ k and τ k+1 , the state variable evolves as a diffusion process and the controller makes an integral profit f . At each decided intervention time τ k , he gives an impulse to the system: the state process jumps with a size

X u τ k -X u τ - k = γ(X u τ - k
) and he obtains the intervention gain κ.

We consider standard assumptions on the coefficients of the problem: A straightforward computation using (H), ( 2) and Gronwall's lemma shows that

∀(t, x) ∈ [0, T ] × R d , sup u∈U (t,T ] E sup t≤s≤T X t,x,u s 2 < ∞. (4)
Finally, we will assume the existence of an optimal strategy u * = (τ * k ) k≥1 ∈ U (t,T ] to problem (1). We refer for example to [START_REF] Bensoussan | Impulse control and quasi-variational inequalities[END_REF] and [START_REF] Øksendal | Applied stochastic control of jump diffusions[END_REF] in the infinite horizon case, for specific conditions on the coefficients of the problem which ensures such an existence.

Link to BSDEs with Constrained Jumps

Let us consider the BSDE with constrained jumps

Y t = g(X T ) + T t f (X s )ds - T t Z s dW s - T t V s dN s + T t dK s , ∀0 ≤ t ≤ T V t + κ(X t -) ≤ 0, ∀0 ≤ t ≤ T (5)
where X is the (uncontrolled) jump diffusion process with dynamics

dX t = b(X t )dt + σ(X t )dW t + γ(X t -)dN t . (6)
Under (H), this SDE admits an unique solution in S 2 and it is shown in [START_REF] Kharroubi | Backward SDEs with constrained jumps and Quasi-Variational Inequalities[END_REF] that under the additional assumption (H 1 ) given below, (5) admits a unique minimal

solution (Y, Z, V, K) ∈ S 2 × L 2 (W ) × L 2 (N ) × A 2 with K predictable.
The solution (Y, Z, V, K) is said to be minimal if and only if it has the smallest component Y in the (infinite) class of solutions to [START_REF] Bouchard | Discrete-Time Approximation and Monte-Carlo Simulation of Backward Stochastic Differential Equations[END_REF]. (Y t ) t≥0 is called the value process and jumps with a size

V t = Y t -Y t -. (H 1 ) There exists a solution ( Ỹ , Z, K) ∈ S 2 × L 2 (W ) × A 2 to Y t = g(X T ) + T t f (X s )ds - T t Z s dW s + T t κ(X s -)dN s + T t dK s .
(H ′ 1 ) (H 1 ) holds and Ỹt = ṽ(t, X t ), ∀0 ≤ t ≤ T for some ṽ with linear growth.

(H 2 ) There exists a non negative function ϕ ∈ C 2 (R d ) and a constant ρ > 0 s.t.

Lϕ + f ≤ ρϕ, ϕ -Hϕ > 0, ϕ ≥ g, lim |x|→∞ ϕ(x) 1+|x| = ∞,
in which L is the local component of the generator of the process X and H, the intervention operator:

Lv(t, x) = b(x) • D x v(t, x) + 1 2 Tr σσ ⊥ (x)D 2 x v(t, x) , Hv(t, x) = v(t, x + γ(x)) + κ(x) .
Let (Y t,x s , Z t,x s , V t,x s , K t,x s ) t≤s≤T be the solution to (5) when X ≡ (X t,x s ) t≤s≤T is the solution starting at x in t to SDE [START_REF] Bouchard | Monte-Carlo valorisation of American options: facts and new algorithms to improve existing methods[END_REF]. Under assumptions (H), (H ′ 1 ) and (H 2 ), [START_REF] Kharroubi | Backward SDEs with constrained jumps and Quasi-Variational Inequalities[END_REF] show that the solution to impulse control problem (1) coincides with initial value of component Y t,x :

Y t,x t = v(t, x) (7)
and is equal to the (unique) solution with linear growth to quasi-variationnal inequality

min -∂v ∂t (t, x) -Lv(t, x) -f (t, x) ; v(t, x) -Hv(t, x)} = 0, ∀(t, x) ∈ [0, T ) × R d , min {v(T -, x) -g(x); v(T -, x) -Hv(T -, x)} = 0, ∀x ∈ R d .
(8) Let us mention that in general, the terminal condition v(T -, •) = g is irrelevant, because of the possible discontinuity of Y in T -due to constraints: the relaxed terminal condition in [START_REF] Carmona | Pricing Asset Scheduling Flexibility Using Optimal Switching[END_REF] expresses the possibility of a jump at time T -. Remark 1. For a better intuition, the following interpretation to solution (Y, Z, V, K) holds when assuming v ∈ C 1,2 ([0, T ], R d ):

∀0 ≤ t ≤ T, Y t = v(t, X t ) Z t = σ(t, X t -)D x v(t, X t -) V t = v(t, X t -+ γ(X t -)) -v(t, X t -) = Hv(t, X t -) -v(t, X t -) -κ(X t -) K t = t 0 -∂v ∂t -Lv -f (s, X s )ds.
The constraint in [START_REF] Bouchard | Discrete-Time Approximation and Monte-Carlo Simulation of Backward Stochastic Differential Equations[END_REF] means thus that the obstacle condition is satisfied, namely v(t, X t -) -Hv(t, X t -) ≥ 0.

Convergence of the Numerical Approximation by Penalization

It does not seem possible to use the minimality condition of the solution to BSDE with constrained jumps (5) directly in a numerical scheme. We thus propose an approach by penalization of the jump constraint. The penalized constraint is introduced in the BSDE driver: when the constraint is fulfilled, this penalization term disappears, and otherwise penalizes the driver with an exploding factor p.

In Theorem 1, we provide an explicit rate of convergence of our approximation with respect to the parameters introduced: namely, the jump intensity λ, the penalization coefficient p and the time step. Such an error estimate is essential for numerical purposes (understanding of the numerical impact of those parameters) and allows to adjust in practice the fineness of the time grid in relation to (λ, p).

Approximation by Penalization

Given a parameter value p > 0, the penalized BSDE is:

Y p t = g(X T ) + T t f (X s ) + p (V p s + κ(X s -)) + λ ds (9)
- T t Z p s dW s - T t V p s dN s , ∀0 ≤ t ≤ T which admits an unique solution (Y p , Z p , V p ) ∈ S 2 × L 2 (W ) × L 2 (N )
from the classical theory of BSDEs with jumps. In addition, the sequence of penalized solutions 5) as p goes to infinity, see [START_REF] Kharroubi | Backward SDEs with constrained jumps and Quasi-Variational Inequalities[END_REF]. Besides, the convergence of (Y p ) p to Y is monotone and increasing.

(Y p , Z p , V p ) p tends in L 2 F ([0, T ]) × L 2 (W ) × L 2 (N ) to the minimal solution (Y, Z, V ) to (
Let (Y p,t,x , Z p,t,x , V p,t,x ) be the solution to (9) when X ≡ (X t,x s ) t≤s≤T . We consider the following error introduced by this penalization procedure:

E p := sup 0≤t≤T v(t, x) -Y p,t,x t . ( 10
)
For any t < η ≤ T , let us introduce

v η T (t, x) = sup u=(τ k ) k≥1 ∈U (t,T -η] E   g(X t,x,u T ) + T t f (X t,x,u s )ds + k≥1 t<τ k ≤T κ(X t,x,u τ - k )    (11)
which corresponds to initial problem (1) restricted to the sub-set of strategies taking values in (t, Tη]. We shall denote by

u η * = (τ η * k ) k≥1 an η 1 2
-optimal strategy to problem [START_REF] Cont | Financial modelling with jump processes[END_REF] (the existence of an optimal strategy is not ensured) and by n η * be the number of impulses in strategy u η * that is:

n η * := ♯ k ≥ 1 : t < τ η * k ≤ T -η .
We will use the following additional assumptions:

(H n ) There exists some n ∈ N * such that ∀j ≥ n, P (n η * ≥ j) ≤ l(j)
for some map l such that l(j) ≤ e -Cj for some some constant C > 0.

(H * ) There exists a map h such that h

(ε) = O ε→0 (ε 1
2 ) and

∀ε > 0, P min k≥1 τ η * k+1 -τ η * k ≤ ε ≤ h(ε).
Remark 2 (Assumptions (H n ) and (H * )). Both assumptions (H n ) and (H * ) are directly satisfied for the problem of Swing options valuation since the number of exercises right is almost surely bounded by some n max and there is some fixed time delay δ > 0 between two consecutive interventions. More generally, (H n ) is intuitively satisfied as soon as the controlled state variable is constrained almost surely and admits jumps of constant sign, see Example 1. (H * ) is verified for sufficiently smooth problems, see for example the case of optimal forest management studied in [START_REF] Bernhart | Modelization and valuation methods of gas contracts[END_REF].

Example 1. Let us assume that the state variable defined in ( 3) is such that • b is uniformly bounded and σ > 0 constant,

• for some constant c > 0, sup

x∈R d γ(x) ≤ -c,
and that the optimal strategy u * implies X u * T ≥ 0 a.s. Then a straightforward computation shows that the (random) number n * (0,T ] of optimal impulses before time T satisfies, for any a > 0,

P n * (0,T ] > n = O e -an
as n → +∞.

Proposition 1. Let assumptions (H), (H n ), (H * ), (H ′ 1 ) and (H 2 ) be satisfied. Then the penalization error in [START_REF] Chassagneux | Discrete-time Approximation of Multidimensional BSDEs with oblique reflections[END_REF] admits the following bound as p goes to infinity:

E p ≤ C n C n (λp) 1 2 -α , ∀α ∈ 0, 1 2 .
for some constants C > 0 and C > 1, which do not depend either on λ, p, n or α.

Proof. We provide the main arguments of the proof and refer the reader to [START_REF] Bernhart | Modelization and valuation methods of gas contracts[END_REF] for more details. The main idea comes from the following explicit functional representation available for Y p,t,x , see [START_REF] Kharroubi | Backward SDEs with constrained jumps and Quasi-Variational Inequalities[END_REF]:

Y p,t,x t = ess sup ν p ∈V p E ν p g(X t,x T ) + T t f (X t,x s )ds + T t κ(X t,x s -)dN s , (12)
where E ν p denotes the expectation under the probability measure P ν p equivalent to P on (Ω, F T ) with Radon-Nikodym density

dP ν p dP FT = e -T 0 (ν p s -1)λds e T 0 ln(ν p s )dNs .
The specificity of such a change of measure is that it impacts only the jump parts of the processes: under P ν p , the Brownian motion W remains unchanged whereas N has a (stochastic) intensity (λν p s ) s≥0 . We shall denote by N p , the doubly stochastic Poisson process (Cox process) with intensity (λν p s ) s≥0 under P, by (τ p k ) k≥1 the sequence of its jump dates and by X p the solution to (6) driven by N p .

In view of ( 11) and ( 12), we introduce a convenient measure change, which intuitively forces the penalized solution to jump as soon as possible after that an optimal impulse occurs:

∀s ≥ 0, ν p s = p if k≥1 1 {τ η * k <s} = N p s , 0 else. ( 13
)
Notice that ν p is a P-measurable process bounded by p a.s. By definition of the counting process N p ,

∀s ≥ 0, P   τ p k -τ η * k > s | j≥1 1 {τ p j ≤τ η * k } = k -1   = e -λps .
In other words, the increment τ p kτ η * k has an exponential distribution with parameter (λp), conditionally to the fact that X p has jumped one time less than X u η * . This allows us to compute an estimate of the distance

sup 0≤t≤T v η T (t, x) -Y p,t,x t ,
which holds for E p by sending η to 0 together with a continuity argument of the value function in its maturity variable.

Convergence Rate of the Numerical Scheme

Given a regular time grid π = {t 0 = 0, t 1 , . . . , t N = T }, we assume that the solution X to (6) can be simulated on π either perfectly or by using an Euler scheme and denote its discrete-time version by X π . Along the lines of [START_REF] Bouchard | Discrete-time approximation of decoupled forwardbackward SDE with jumps[END_REF], we consider the following backward discrete-time scheme for numerically solving the penalized BSDE with jumps (9),

                       Y p,π tN = g(X π tN ) ∀t n ∈ π, t n < T : V p,π tn = 1 λ∆tn+1 E Y p,π tn+1 ∆ Ñtn+1 |F tn Z p,π tn = 1 ∆tn+1 E Y p,π tn+1 ∆W tn+1 |F tn Y p,π tn = E Y p,π tn+1 |F tn + f (X π tn ) + p V p,π tn + κ(X π tn ) + -V p,π tn λ ∆t n+1 (14)
where ∆t n+1 = t n+1t n , ∆W tn+1 is the Brownian increment on [t n , t n+1 ] and ∆ Ñtn+1 the compensated version of the Poisson increment ∆N tn+1 on [t n , t n+1 ).

We consider the classical discretization error between the continuous-time solution (Y p , Z p , V p ) in [START_REF] Carmona | Optimal multiple stopping and valuation of Swing options[END_REF] and its discrete-time approximation (Y p,π , Z p,π , V p,π ) in [START_REF] Djehiche | Popier: The Finite Horizon Optimal Multiple Switching Problem[END_REF], that is

E π (Y p ) := max n<N -1 sup tn≤t≤tn+1 E Y p t -Y p,π tn 2 1 2 E π (Z p ) := N -1 n=0 tn+1 tn E Z p t -Z p,π tn 2 dt 1 2 E π (V p ) := N -1 n=0 tn+1 tn E V p t -V p,π tn 2 λdt 1 2
.

Because of the lack of first order regularity of the driver of the penalized BSDE

f p (x, v) := f (x) + p(v + κ(x)) + -v λ, ∀(x, v) ∈ R d × R, (15)
classical regularization arguments for the FBSDE coefficients and Malliavin differentiation representations allow us to provide an explicit convergence rate of order |π| 1 2 for errors E π (Y p ) and E π (V p ), but only of order |π| 1 4 for error E π (Z p ), see Proposition 2.

The impact of the penalization coefficient p on the convergence of backward discrete-time schemes is well known in practice, even if there does not exist, to our best knowledge, any explicit computation in the literature (see for example the numerical experiments of [START_REF] Lemor | Approximation par projections et simulations de Monte-Carlo des équations différentielles stochastiques rétrogrades[END_REF] for the resolution by penalization of a BSDE with one reflecting barrier). Basically, as p increases at fixed discrete-time step, the quantity f p (•)∆t n+1 explodes, leading to a numerical explosion of the approximate values Y p,π , see [START_REF] Djehiche | Popier: The Finite Horizon Optimal Multiple Switching Problem[END_REF].

We show rigorously in Proposition 2 that the discretization error grows exponentially with (λp 2 ). This is due to the linear dependence in λ and p of the BSDE driver f p , see [START_REF] Elie | Contrôle stochastique et méthodes numériques en finance mathématique[END_REF], and estimate computations based on the use of Gronwall's lemma.

Proposition 2. Assume (H). Then, as soon as

|π| = O 1 λp 2 , ( 16
)
we get the following bounds as p goes to infinity:

E π (Y p ) ≤ C (1 + λ) 2 λp Cλp 2 |π| 1 2 E π (V p ) ≤ C (1 + λ) 2 λ 3 2 p 2 Cλp 2 |π| 1 2
.

Under (H ′ ), there exists a version of Z p such that,

E π (Z p ) ≤ C (1 + λ) 2 λ 5 2 p 3 Cλp 2 |π| 1 4
, for some constants C > 0 and C > 1, which do not depend either on λ, p or |π|.

Proof. This follows from the same arguments as [START_REF] Elie | Contrôle stochastique et méthodes numériques en finance mathématique[END_REF]: computations using Itô and Gronwall's lemma and regularization and Malliavin differentiation arguments applied to the penalized BSDE with jumps [START_REF] Carmona | Optimal multiple stopping and valuation of Swing options[END_REF]. We refer the reader to [START_REF] Bernhart | Modelization and valuation methods of gas contracts[END_REF] for a detailed proof.

Propositions 1 and 2 enable us to establish a global convergence rate of the error introduced by our approximation by penalization.

Theorem 1. Let the assumptions of Proposition 1 be satisfied. Then

E p + E π (Y p ) ≤ C 1 (λp) 1 2 -α + (1 + λ) 2 λp Cλp 2 |π| 1 2 , ∀α ∈ 0, 1 2 (17)
for some constants C > 0 and C > 1, which do not depend either on λ, p, |π|, n or α. Thus, for a sufficiently small time step |π| with respect to λ and p, the global error is such that

[E p + E π (Y p )] * = O 1 (λp) 1 2 -α , ∀α ∈ 0, 1 2 .
Remark 3 (Global convergence rate). At fixed time step |π|, the convergence rate strongly deteriorates as λ or p increases, see [START_REF] Gobet | A regression-based Monte Carlo method to solve backward stochastic differential equations[END_REF]. The numerical method is more sensitive to p than to λ according to ( 16) and ( 17). ( 16) constitutes a necessary condition for the convergence of the backward discrete-time scheme. In practice, the penalization parameter will need to be chosen relatively small and the time step |π| very small to avoid multiple jump times on each time step (otherwise, this introduces a bias).

Application to Swing Options Valuation

In this section, our method is applied for the valuation of Swing options in the Black and Scholes framework. This constitutes a multiple optimal stopping time problem which can be reformulated as a particularly degenerate three-dimensional impulse control problem. We have been able to achieve convergence for a small maximal number of exercises rights (n max ≤ 2) due to the slow computational speed of our method.

Swing Options Valuation as an Impulse Control Problem

We consider a (normalized) Swing option: the holder of the option is given a maximal number of exercise rights, say n max ≥ 1, and has the opportunity to sell whenever he wants over a time period [0, T ] an underlying asset against a fixed strike price. For some fixed strike price K, we shall denote by φ(s) = (Ks) + the reward function corresponding to the profit made at each exercise date and by S the underlying asset spot price. We concentrate here on the risk-neutral Black and Scholes framework in which r > 0 is a constant interest rate and the spot price process is defined by

S t = s exp (r - 1 2 σ 2 )t + σW t , ∀t ≥ 0, (18)
where σ > 0 denotes the volatility coefficient and s an initial price value.

A delay δ > 0 between two consecutive exercise dates is introduced. Indeed, without any delay, the optimal strategy would consist in n max simultaneous exercise at an unique optimal date (so that this option is equivalent to n max identical American options). The value of such an option at time 0 can be written as the solution to the following multiple optimal stopping time problem

v(nmax) (s) = sup u=(τ k ) k≥1 ∈U δ,nmax (0,T ] E   k≥1 e -rτ k φ(S τ k )   (19)
in which a strategy u = (τ k ) k≥1 is said to be admissible and belongs to U δ (0,T ] if and only if it is a vector of F W -stopping times valued in (0, T ] with maximal lenght n max which satisfied the constraint on delay, that is

∀k ≥ 1, τ k+1 -τ k ≥ δ.
As a multiple optimal stopping time problem, this problem can be formulated as a particular (and strongly degenerate) impulse control problem. An impulse control corresponds to a sequence of exercise dates and the intervention gain is written as the payoff function φ multiplied by an indicator function which allows to satisfy both the constraint on the number of exercise rights and the constraint on delay between exercise dates. Namely:

v(s) = sup u=(τ k ) k≥1 E    k≥1 τ k ≤T e -rτ k φ(S τ k )1 Θ k-1 τ - k ≥δ ∩ Q u τ - k <nmax    , (20)
where u = (τ k ) k≥1 is the sequence of F W -stopping times valued in (0, T ] and two additional state variables which both are controlled and discontinuous (càdlàg) processes are introduced:

• Q u counts the number of exercise rights used before considered time

Q u 0 = 0, Q u t = ♯ {k ≥ 1, τ k ≤ t} , ∀t ≥ 0.
• Θ u t := Θ k t = inf {tτ k , τ k ≤ t} corresponds to the delay between t and last exercise date

Θ k t = t -τ k , ∀τ k ≤ t < τ k+1 , Θ k τ k+1 = 0, ∀k ≥ 0,
where by convention Θ u 0 = Θ 0 0 = 0. Obviously, the problem ( 19) is equivalent to the impulse control problem [START_REF] Hu | Multi-dimensional BSDE with oblique reflection and optimal switching[END_REF]:

∀s ∈ R, v(nmax) (s) = v(s).

Penalized BSDE Associated to Swing Option Pricing Problem

Let us introduce the uncontrolled variable (Q, Θ) defined by

Q t = N t , ∀t ≥ 0, Θ t = t - t 0 Θ s -dN s , (21)
where we recall that N is a Poisson process with intensity λ > 0. For a penalization coefficient p > 0, the penalized BSDE with jumps associated to problem [START_REF] Hu | Multi-dimensional BSDE with oblique reflection and optimal switching[END_REF] is

Y p t = κ(S T , Q T -, Θ T -) - T t rY p u du + p T t (V p u + κ(S u , Q u -, Θ u -)) + λdu - T t Z p u dW u - T t V p u dN u , ∀0 ≤ t ≤ T. ( 22
)
in which the intervention gain is defined by

κ(s, q, θ) := φ(s)1 {(θ≥δ)∩(q≤nmax-1)} , ∀(s, q, θ) ∈ R × N × R + . ( 23
)

Numerical Valuation Algorithm

The three-dimensional process (S, Q, Θ) can be exactly computed on the time grid π. We shall denote by (S π , Q π , Θ π ) its version on π. In particular, the pure jump processes Q and Θ can be computed without approximation error, recall [START_REF] Kharroubi | Backward SDEs with constrained jumps and Quasi-Variational Inequalities[END_REF], from a simulated trajectory of a Poisson process with intensity λ, see for example [START_REF] Cont | Financial modelling with jump processes[END_REF].

As the driver of the penalized BSDE with jumps [START_REF] Vath | A model of portfolio selection under liquidity risk and price impact[END_REF] does not depend on Z p , it is sufficient to compute (Y p,π , V p,π ) on π, which can be made by the backward recursive scheme:

                 Y p,π tN = κ(S π tN , Q π tN , Θ π tN ) ∀t n ∈ π, t n < T : V p,π tn = 1 λ∆tn+1 E Y p,π tn+1 ∆ Ñtn+1 |F tn Y p,π tn = 1 1+r∆tn+1 E Y p,π tn+1 |F tn + p V p,π tn + κ(S π tn , Q π tn , Θ π tn ) + -V p,π tn λ∆t n+1 (24)
and in our setting:

E [•|F tn ] = E •|(S π tn , Q π tn , Θ π tn ) . ( 25
)

Monte Carlo-based Resolution

We compute estimators of the conditional expectations ( 25) by a classical least squares Monte Carlo technique. We use least squares regressions on adaptative local basis functions, see [START_REF] Bouchard | Monte-Carlo valorisation of American options: facts and new algorithms to improve existing methods[END_REF]. Since such an approach is only relevant for real-valued variables (the regression basis functions have compact support), it cannot handle the integer-valued variable Q π .

We thus deal with this variable explicitly. Namely, we simulate

M ≥ 1 i.i.d. paths of (S π , Q π , Θ π ) (S π,(m) , Q π,(m) , Θ π,(m) ), ∀m ≤ M,
and at each time step t n < T of the backward recursion, the set of M Monte Carlo samples is separated in n max + 1 sub-sets corresponding to the samples on which Q π tn = 0, 1, . . . , n max -1 and Q π tn ≥ n max . Then, we just need to estimate n max conditional expectation operators, namely

E •|(S π tn , Q π tn = q, Θ π tn ) , ∀q ∈ {0, 1, . . . , n max -1}
by using the corresponding Monte Carlo samples, since (Y p,π tn , V p,π tn ) = (0, 0) when

Q π tn ≥ n max , see Remark 4. Remark 4. Let us show that Q π tn ≥ n max implies Y p,π tn = V p,π tn = 0. We have E Y p,π tn+1 ∆ Ñtn+1 |(S π tn , Q π tn ≥ n max , Θ π tn ) = E E Y p,π tn+1 |Q π tn ≥ n max = 0 ∆ Ñtn+1 |(S π tn , Q π tn ≥ n max , Θ π tn ) = 0,
by definition of κ in ( 23) and [START_REF] Ludkovski | Optimal Switching with Applications to Energy Tolling Agreements[END_REF]. In the same way

E Y p,π tn+1 |(S π tn , Q π tn ≥ n max , Θ π tn ) = 0
which allows to conlude.

The Monte Carlo-based algorithm is then the following:

I. Initialization: Y p,π,(m) tN = φ(S π,(m) tN )1 Θ π,(m) t N ≥δ ∩ Q π,(m) t N ≤nmax-1
, ∀m ≤ M.

II. Computation backward in time of (V p,π,(m) , Y p,π,(m) ) on each sample m ≤ M .

For n = N -1, . . . , 0, set

M q tn := m = 1, . . . , M : Q π,(m) tn = q , ∀q ≤ n max -1, M nmax tn := m = 1, . . . , M : Q π, (m) tn
≥ n max .

Then:

1. For any m ∈ M nmax tn , V p,π,(m) tn = Y p,π, (m) tn
= 0.

2. Set q := n max -1.

3. If q ≥ 0, for any m ∈ M q tn , the conditional expectations estimators

ε V,q,(m) tn ≈ E Y p,π tn+1 ∆ Ñtn+1 | S π, (m) tn , Θ π,(m) tn ε
Y,q,(m) tn ≈ E Y p,π tn+1 | S π, (m) tn , Θ π,(m)
}. Then,            V p,π,(m) tn = 1 λ∆tn+1 ε V,q,(m) tn Y p,π,(m) tn = 1 1+r∆tn+1 ε Y,q,(m) tn + p V p,π,(m) tn + φ(S π,(m) tn )1 Θ π,(m) tn ≥δ + -V p,π,(m) tn λ∆t n+1 .
4. q := q -1 and go to 3. III. At time t 0 (M 0 t0 = {1, . . . , M } and M q t0 = ∅, ∀q ≥ 1) the Swing option price estimator is given by Y p,π t0 such that

   V p,π t0 = 1 λ∆t1 1 M M m=1 Y p,π,(m) t1 ∆ Ñ (m) t1 Y p,π t0 = 1 1+r∆t1 1 M M m=1 Y p,π,(m) t1 + p(V p,π t0 ) + -V p,π t0 λ∆t 1 .
Let us highlight some features of the above-presented Monte Carlo procedure. At each backward induction date t n < T , we have to estimate in worst cases 2 × n max conditional expectations, performed on each subset M q tn , q ≤ n max -1. When n max increases, much more Monte Carlo samples are needed as each least squares regression requires a sufficient number of samples.

In addition, the number of local basis functions b q tn has to be adapted to the number of Monte Carlo samples used for the least squares regression, namely card(M q tn ). We thus introduce a dynamic choice for b q tn : it is fixed proportionally to card(M q tn ) for any q ≤ n max -1 and t n < T .

Remark 5 (Statistical error of our method). A control of the statistical error introduced by the least squares Monte Carlo approach is provided in Gobet et al. [START_REF] Gobet | Rate of convergence of empirical regression method for solving generalized BSDE[END_REF] (see [START_REF] Lemor | Approximation par projections et simulations de Monte-Carlo des équations différentielles stochastiques rétrogrades[END_REF] for further details). By extension, this applies to BSDE with jumps (see [START_REF] Elie | Contrôle stochastique et méthodes numériques en finance mathématique[END_REF]) to control the error on the jump component V p,π and thus ensures that the least squares Monte Carlo error tends to 0 as the number of samples M and the number of basis functions b tends to +∞.

A Benchmark Method Based on Iteration

The classical method to value such a Swing option, recall formulation [START_REF] Hamadène | On the Starting and Stopping Problem: Application in reversible investments[END_REF], is based on an iteration over the number of exercise rights, see for example [START_REF] Carmona | Optimal multiple stopping and valuation of Swing options[END_REF]. The dynamic programming principle provides a direct link between the solution v (j) to the same problem as [START_REF] Hamadène | On the Starting and Stopping Problem: Application in reversible investments[END_REF] but with at most j ≤ n max exercise rights and the solution v (j-1) with at most (j -1) exercise rights.

The value of the Swing option with 0 exercise right is obviously zero v (0) = 0. Then, we compute the sequence of values of Swing options with j exercise rights v (j) , j = 1, . . . , n max according the backward recursion scheme:

                     v (j) (t N , s) = φ(s) ∀t n ∈ π, T -δ < t n < T : v (j) (t n , s) = max φ(s) ; e -r∆tn+1 E (tn,s) v (j) (t n+1 , S π tn+1 ) ∀t n ∈ π, t n ≤ T -δ : v (j) (t n , s) = max φ(s) + e -rδ E (tn,s) v (j-

Pricing Results

We consider put options with maturity T = 1 year and a strike price K = 100. The Black and Scholes parameters, see [START_REF] Gobet | Rate of convergence of empirical regression method for solving generalized BSDE[END_REF], are r = 0.05, σ = 0.3 and s = 100.

Special Case of American Options: n max = 1

In the single-exercise case, the additional variable Θ disappears (there is no delay constraint). This helps simplify the Monte Carlo procedure described in Paragraph 4.2. In particular, the algorithm implies only two sequences of samples subsets, whether one jump of N occurs before T on the considered path or not: that is ((M 1 ) tn ) n≤N -1 and ((M 0 ) tn ) n≤N -1 .

In our numerical experiments, we find out that increasing too much λ makes the variance of the Monte Carlo procedure explode. It would be necessary to increase the number of Monte Carlo samples, which leads to prohibitive computational times (each pricing result presented below was obtained after a computation between 6 and 8 hours). For the same reason (exploding behavior of the penalized BSDE driver), we restrict our numerical experiments to penalization parameters ≤ 5.

The benchmark price for the American put option is 9.88 (by a binomial approach or classical Monte Carlo). We report in Table 1 the price given by our method when varying λ and the number of time steps N for a penalization parameter equal to 5. We used 20 million of Monte Carlo paths.

In all the experiments that we performed in this simple case, we numerically observed that the limiting prices of our method (with respect to N ) are below the benchmark value: this is due to penalization. For each considered value of λ, we retrieve a convergence in the number of time steps N of our method. As p = 5, approximate prices are converged from N = 160, so that we restrict ourselves to N ≤ 320 time steps as p = 10. The limiting values are still below the benchmark but accurate option prices (relative error less than 1%) are obtained with a penalization coefficient p equal to 10 and N = 160. See also Table 2 in which the (signed) relative error to the benchmark is given in brackets. Besides, we observe a monotone convergence in λ of our approximate method.

We should point out that fine-tuning the parameters of the algorithm is difficult. As already mentioned, since the number of Monte Carlo paths is different in each set of sample paths M q tn , q = 0, 1, 2, the number of basis functions used for the least squares regressions has to be dynamically adapted. And when increasing much more the jump intensity λ, more Monte Carlo samples would be necessary. For such a Swing option, the running time is much longer because the conditional expectations are computed by regression with respect to the bidimensional state variable (S π , Θ π ). The computation of one option price takes at least 15 hours in above cases (when N ≥ 80). In comparison, the benchmark method takes less than 5 minutes. Besides, the complexity of our method increases with n max , leading to untractable computational times for bigger values of n max , see Remark 6.

On this particular case of Swing options valuation, it seems that our method is less competitive than the classical approach. This is without any doubt due to the strong degeneracy of such a problem in our impulse control context: the valuation problem is 3-dimensional and involves an additional integer-valued state variable Q representing the number of exercise rights used at any considered time.

However, our method works and the numerical results that we obtain are consistent with the theoretical convergence rate given in Theorem 1. One can expect that our method would work better on less degenerate problems.

Remark 6 (Dealing with more exercise rights). The computational time of our method intuitively increases linearly with the number of exercise rights n max . Indeed, at each time step of the backward induction procedure, the number of conditional expectation estimations is proportional to n max . Besides, when multiplying by 2 the number of exercise rights, it would require, at least, a double number of Monte Carlo samples for a same accuracy of the computation of conditional expectation estimators.

Let us mention that the computational time of the benchmark method using iteration increases linearly as function of the the maximal number of exercise rights as well.

  (H) and (H ′ ) below, and constants T , |b(0)|, |σ(0)|, |γ(0)|, |f (0)|, |κ(0)| and |g(0)|.

  (H) b : R d → R d , σ : R d → R d×d and γ : R d → R d are Lipschitz continuous and γ is uniformly bounded. f : R d → R, κ : R d → R and g : R d → R are Lipschitz continous. (H ′ ) The maps b, σ, γ, and g belong to C 1 b (R d ) and have Lipschitz continuous derivatives.

tn are approximated by least squares regression of Y

  with b := b q tn basis functions {ψ 1 , . . . , ψ b

4. 3 . 2 , 2 10 , 3 10 .

 321010 Figures 1, 2 and 3 the corresponding approximate prices when varying λ and N for a penalization parameter equal to 5 and 10 (we used 40 million of Monte Carlo paths and N = 20, 40, 80, 160, 320, 640).

Figure 1 :

 1 Figure 1: Approximate prices of a Swing option with 2 exercise rights and δ = 1 10 , with p = 5 (left) and p = 10 (right)

Figure 2 :

 2 Figure 2: Approximate prices of a Swing option with 2 exercise rights and δ = 2 10 , with p = 5 (left) and p = 10 (right)

17

 17

Figure 3 :

 3 Figure 3: Approximate prices of a Swing option with 2 exercise rights and δ = 3 10 , with p = 5 (left) and p = 10 (right)

Table 1 :

 1 Approximate prices of an American option with p = 5 where E (tn,s) [•] := E •|S π tn = s . We use the same least squares Monte Carlo regression-based method for approximating the conditional expectations operators as above.

	1) (t n + δ, S π tn+δ ) ;
	e -r∆tn+1 E (tn,s) v (j) (t n+1 , S π tn+1 )

Table 2 :

 2 Prices of a Swing option with 2 exercise rights (limiting values with N = 160)

		H H H H H H p λ 3	4	5
	δ = 1 10	5 10	18.80 (-2.44%) 18.95 (-1.66%) 19.00 (-1.40%) 19.16 (-0.57%) 19.24 (-0.16%) 19.27 (0.00%)
	δ = 2 10	5 10	18.23 (-2.88%) 18.37 (-2.13%) 18.43 (-1.81%) 18.59 (-0.96%) 18.65 (-0.64%) 18.69 (-0.43%)
	δ = 3 10	5 10	17.62 (-2.97%) 17.76 (-2.20%) 17.81 (-1.93%) 17.99 (-0.94%) 18.04 (-0.66%) 18.05 (-0.61%)


) Maier
email: ouanezar@tsi.enst.fr

Performance analysis of an off-line BMI for the control of a robotic digit (1,3)

Keywords: BMI, neurobotics, motor cortex, EMG, neural networks, computational neuroscience

. Here we report the performance analysis based on the same data and robot, the statistical results including aspects of learning/training, prediction, as well as control.

Introduction

Recent developments in BMI technologies have shown how cortical activities can be used to control the position and/or the state of different types of external devices in real time. These efforts have focused mainly on controlling a computer cursor [START_REF] Wu | Neural decoding of cursor motion using a Kalman filter[END_REF], [START_REF] Mijail | Brain-machine interface: Instant neural control of a movement signal[END_REF], the state of on/off switches or the decoding and control of reach trajectories of a robotic arm, sometimes equipped with a simple 1 DoF gripper [START_REF] Ouanezar | A BMI application for the control of a robotic digit[END_REF][START_REF] Velliste | Cortical control of a prosthetic arm for self-feeding[END_REF][START_REF] Carmena | Learning to control a brain-machine interface for reaching and grasping by primates[END_REF][START_REF] Taylor | Direct cortical control of 3d neuroprosthetic devices[END_REF]. However, studies dedicated to the decoding of cortical signals for the eventual control of a dexterous hand with multiple degrees of freedom are rare. The goal would be to control individuated and/or combined finger movements. First attempts of decoding finger movements have only been undertaken very recently [START_REF] Georgopoulos | Neural coding of finger and wrist movements[END_REF][START_REF] Aggarwal | Asynchronous decoding of dexterous finger movements using M1 neurons[END_REF][START_REF] Acharya | Decoding individuated finger movements using volumeconstrained neuronal ensembles in the M1 hand area[END_REF][START_REF] Aggarwal | Cortical decoding of individual finger and wrist kinematics for an upper-limb neuroprosthesis[END_REF][START_REF] Shin | Neural decoding of finger movements using Skellam-based maximum-likelihood decoding[END_REF]. None of these previous results have so far been applied to robot artifacts. Many studies on the cortical activity of the motor cortex in behaving non-human primates (e.g. [START_REF] Maier | Contribution of the monkey corticomotoneuronal system to the control of force in precision grip[END_REF]) have shown that neurons in the primary motor cortex (M1) code for some aspects of finger movements, such as force. In the aim of providing a proof-of-concept that a robot finger of multiple DoF can be controlled by intra-cortical signals, we first set up an experimental site (robotic digit) and implemented a kinematic control of its movement (i.e. positional control). To mimic the monkey finger, we opted for an anthropomorphic robot finger of 4 DoF (designed and manufactured by Shadow Robot Company, UK) actuated by McKibben-like pneumatic artificial muscles in an antagonist actuation scheme. This robot digit is equipped with hall effect position sensors and with fingertip quantum tunnel tactile sensors. The robot digit was then mounted on an experimental set-up, which contained additional mechanical parts and sensors to imitate the monkey task (precision grip). Finally, a hybrid motor control scheme composed of an artificial neural network and PID modules has been implemented [START_REF] Ouanezar | A BMI application for the control of a robotic digit[END_REF]. We then performed a statistical analysis on the decoding of previously recorded cortical data for controlling finger position. This is reported in the following.

Method

The sensory-motor control of an anthropomorphic robot finger using intra-cortical signals is composed of three main steps:

• Recording the data.

• Decoding the neural data to predict grip motion (BMI part). • Controlling the robot.

Data recording

The biological signal consists of the activity of 33 corticomotoneuronal (CM) cells recorded in M1 of a macaque monkey performing a precision grip task. Recordings were obtained in 7 separate behavioral sessions. In addition, EMG signals were recorded from up to eight intrinsic hand and forearm muscles. These recordings were obtained from the Institute of Neurology (UCL, London, UK, courtesy of RN Lemon).

Decoding spiking data

We attempted to decode two different signals from the CM cell spike trains: (i) the time-varying endpoint position of the index finger during precision grip. The position was recorded via a spring-loaded lever in the experimental setup. (ii) the EMG of extrinsic or intrinsic hand muscles recorded during this same task.

We used a neural network approach, a Time Delayed Multi-Layer Perceptron (TDMLP) to decode and thus to predict either the finger position or the EMG from the CM cell spike trains. The input vector consisted of a binary sequence corresponding to the presence/absence of spikes of the CM cell (down-sampled to 250 Hz) over a given period, i.e. the duration of the sliding input window, which varied between 25 ms and 400 ms.

The TDMLP, a three-layer feed-forward neural network was implemented with the Matlab neural network tool box. The input layer consists of D × NC units (D is the number of bins in the sliding window, and NC the number of CM cells used) which represent the spiking response recorded between time t-D and t. Sigmoid activation functions (range 0-1) were used in the hidden layer, which consisted of 10 units. The one output unit, also using a sigmoid activation function, provides a positive scalar that corresponds to the predicted position (or EMG) at time t. This is illustrated in Fig. 3.

Hidden layer units and the output unit also received bias inputs. Bias and connection weights are adjusted during the learning process, based on the gradient descent error to minimize the sum of the squared errors between the desired output and the network output. The number of units in the hidden layer and the particular learning algorithm have been chosen by the technique of crossvalidation.

Learning of the input-output transformation occurred through the use of multiple trials (min: 15, max: 92 trials) depending on the recording session. The number of learning epochs was fixed to 1000, each epoch containing all trials. The Nguyen-Widrow algorithm was used for weight initialization.

Controlling the robot digit

The third and final part of the BMI based motor control chain consists in the reproduction of estimated index finger movement pushing on one of the manipulandum's lever, shown in Fig. 2, using the 4 DoF robot finger. On the finger there are four joints and three human-sized phalanges. Both middle and proximal (MCP) joints are actuated in the flexion/extension plane by a pair of antagonist artificial muscles. The MCP joint, also responsible for the abduction/adduction motion of the proximal phalange, is in addition driven by one muscle with an opposing return spring. Similar to the human finger the distal joint is coupled to the middle joint such that the angle of the middle joint is always greater than or equal to the angle of the distal joint. The joint angles relative to the estimated fingertip trajectory were calculated by an inverse kinematics model of the robot finger taking into account the coupled motion of the distal phalange. Finally, we used a classical PID controller, the last link of the control chain, to ensure the tracking of the desired trajectory as shown in Figure 4.

Results

The performance of the TDMLP-based prediction was systematically analysed by varying the size D of the sliding window and the number of cortical cells NC used as inputs. This was done with a TDMLP with 10 units in the hidden layer, a number of hidden units found to be optimal across different sizes of input windows, irrespective of the type of the output signal, i.e. position or EMG. Performance (Per) was calculated as follows: 4 shows the performance of the TDMLP in estimating the recorded position for different sizes of sliding input windows, as well as for two modes of training: (i) a prediction based on the spike train of a single CM cell, and then averaged across cells in a given recording session. This was called "cell-by-cell" prediction. (ii) a prediction based on spike trains of multiple, simultaneously recorded CM cells. This was called prediction of "combined cells". The number of combined cells varied between 2 and 6 depending on the session. Fig. 4 shows the grand average across three recording sessions. The network behaved similarly for predicting position or EMG. Three results were found: (i) In both cases the performance increased continuously as a function of increasing input window size (from 25 to 400 ms). Saturation was observed at a window size of about 400 ms. (ii) Performance was significantly better when using multiple simultaneously recorded spike trains compared to prediction based on single cells. (iii) Independent of input size and number of cells, the prediction of position was better than that of EMG.

∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ = = = = - -- - = = = =                                 = = = =
How the observed performance translates to trial-bytrial prediction is shown in Fig. 6. Recorded (blue) and estimated (red) position signals for three representative trials, based on a TDMLP with a 400 ms input window of 4 concurrent CM cells. Also shown is the prediction of two corresponding EMG signals for Flexor Digitorum Superficialis (FDS) and Extensor Digitorum Communis (EDC).

Recorded and estimated EMG of EDC time[s]

Norm. Ampl. Since each CM cells facilitates a given set of target muscles, it was of interest to investigate whether prediction of EMG varied as a function of the muscle field of the input cells. We hypothesised that the prediction of a target EMG would be better when based on CM cells that facilitated that particular muscle. Presence of facilitation was estimated by post-spikefacilitation (PSF) of spike-triggered averaging of the EMG [START_REF] Maier | Contribution of the monkey corticomotoneuronal system to the control of force in precision grip[END_REF]. Fig. 7. Prediction of FDS EMG with CM cells having FDS as target muscle compared to a prediction based on CM cells without facilitation of FDS (same 3 sessions as in Fig. 5).

Fig. 7 shows the performance of predicting (over 3 sessions) the FDS EMG as a function of the window size. Prediction based on CM cells that facilitated the FDS EMG (green triangles) tends to be better than those based on CM cells that facilitated other muscles than the FDS (yellow triangles).

Finally, we used the predicted position in order to control the robot digit. A comparison between the performance achieved by the TDMLP and that realized by the robot digit shows a systematic lower performance for the robot. The decrease in performance, measured over 25 trails, is in the order of 10% to 20% (Fig. 8). However, this still allows for a functional reproduction of the monkey's digit movement.

Discussion

With few exceptions [START_REF] Georgopoulos | Neural coding of finger and wrist movements[END_REF][START_REF] Aggarwal | Asynchronous decoding of dexterous finger movements using M1 neurons[END_REF][START_REF] Acharya | Decoding individuated finger movements using volumeconstrained neuronal ensembles in the M1 hand area[END_REF][START_REF] Aggarwal | Cortical decoding of individual finger and wrist kinematics for an upper-limb neuroprosthesis[END_REF][START_REF] Shin | Neural decoding of finger movements using Skellam-based maximum-likelihood decoding[END_REF], most studies on the control of the upper limb through data from a cortical invasive BMI have focused on reaching [e.g. 1,2,5,6], not on hand or finger movements. Our results on the control of a single finger demonstrate that neural activity from cells recorded in the M1 hand area contain sufficient information to enable decoding of an asynchronous finger movement.

Several differences distinguish our approach from that of Georgopoulos et al. [START_REF] Georgopoulos | Neural coding of finger and wrist movements[END_REF] or from those of Schieber and Thakor [START_REF] Aggarwal | Asynchronous decoding of dexterous finger movements using M1 neurons[END_REF][START_REF] Acharya | Decoding individuated finger movements using volumeconstrained neuronal ensembles in the M1 hand area[END_REF][START_REF] Aggarwal | Cortical decoding of individual finger and wrist kinematics for an upper-limb neuroprosthesis[END_REF][START_REF] Shin | Neural decoding of finger movements using Skellam-based maximum-likelihood decoding[END_REF]. Both of these groups also worked on decoding of digit movements.

First, we used a neural network approach to obtain a transfer function between spike trains and finger trajectory, whereas they used a coding based on the population vector [START_REF] Georgopoulos | Neural coding of finger and wrist movements[END_REF], on non-linear filters [START_REF] Aggarwal | Asynchronous decoding of dexterous finger movements using M1 neurons[END_REF], or on a maximum-likelihood scheme [START_REF] Shin | Neural decoding of finger movements using Skellam-based maximum-likelihood decoding[END_REF].

Second, they decoded a binary movement direction (flexion or extension) together with a decoding of one of the five digits to be moved [START_REF] Georgopoulos | Neural coding of finger and wrist movements[END_REF]. Essentially, their decoder worked as a movement classifier [START_REF] Acharya | Decoding individuated finger movements using volumeconstrained neuronal ensembles in the M1 hand area[END_REF][START_REF] Shin | Neural decoding of finger movements using Skellam-based maximum-likelihood decoding[END_REF]. In our case, decoding concerned a single finger whose endpoint trajectory was a time-varying scalar. We showed that the accuracy of the prediction increased as a function of the size (duration) of the input window. This was shown previously for hand trajectories [START_REF] Paninski | Superlinear population encoding of dynamic hand trajectory in primary motor cortex[END_REF].

Third, they used the activity of up to 50 simultaneously recorded but unidentified neurons in the M1 hand area as inputs, whereas we used few (up to 5) but identified CM cells. We showed that the accuracy of the prediction increased as a function of the number of simultaneously recorded CM cells acting as inputs. Such an increase in performance has been shown previously for unidentified cells in the case of classifying digit movements [START_REF] Aggarwal | Asynchronous decoding of dexterous finger movements using M1 neurons[END_REF][START_REF] Shin | Neural decoding of finger movements using Skellam-based maximum-likelihood decoding[END_REF], as well as in decoding of the arm trajectory [START_REF] Carmena | Learning to control a brain-machine interface for reaching and grasping by primates[END_REF][START_REF] Sanchez | Interpreting neural activity through linear and nonlinear models for brain machine interfaces[END_REF].

We also predicted EMG activity. Although EMG is an inherently more noisy signal than position, prediction accuracy reached up to 60%, for windows of ≥ 200 ms and multiple CM cells. Slightly higher levels of accuracy (70-75%), but based on many more and unidentified cells, have been obtained previously for EMGs during arm movements [START_REF] Georgopoulos | Neural coding of finger and wrist movements[END_REF]. Furthermore, we show a trend for improved accuracy if the prediction is based on CM cells that facilitate the target muscle EMG.

Finally, the TDMLP-based prediction of finger position was used to replicate the recorded movement with a robotic digit. Our robot differs from those of other BMI applications. We used a robot finger actuated by artificial muscles in an antagonist actuation scheme. Three properties make this robot biologically more plausible than conventional robots: i) one can establish a correspondence between a particular anatomical muscle and its artificial counterpart in the robot. ii) The dynamic control -future work in this project -needs a model of the actuators. It has been demonstrated that actuators of this type are biologically highly plausible [START_REF] Tondu | The McKibben muscle and its use in actuating robot-arms showing similarities with human arm behavior[END_REF]. iii) Finally, an analogy between the EMG signal and the pressurized air used to contract the artificial muscle looks promising. The observed performance difference between estimated and realized robot movement are most likely due to non-resolved compatibility issues in the inverse kinematics model as well as in the PID controller.

In conclusion, decoding of spike trains from simultaneously recorded CM cells in M1 allows for an off-line reproduction of a digit movement during precision grip.

Fig. 1 .

 1 Fig. 1. Left: Schematic diagram of the manipulandum used for the precision grip task. Right: example of the recorded signal during the precision grip task. From top to bottom: force of the index finger and

  Figure taken from[START_REF] Maier | Contribution of the monkey corticomotoneuronal system to the control of force in precision grip[END_REF].

Fig. 2 .

 2 Fig. 2. The control of a robotic digit of 4 DoF using activity of multiple CM cells.

Fig. 3 .

 3 Fig. 3. Schematic of the TDMLP and its input-output transformation.

  Fig.4shows the performance of the TDMLP in estimating the recorded position for different sizes of sliding input windows, as well as for two modes of training: (i) a prediction based on the spike train of a single CM cell, and then averaged across cells in a given recording session. This was called "cell-by-cell" prediction. (ii) a prediction based on spike trains of multiple, simultaneously recorded CM cells. This was called prediction of "combined cells". The number of combined cells varied between 2 and 6 depending on the session. Fig.4shows the grand average across three recording sessions.

Fig. 4 .

 4 Fig. 4. Performance of prediction for fingertip position as a function of window size and number of cortical cells. Open triangles: "cell-bycell" estimation. Filled triangles: "combined cell" estimation. Grand average across 3 sessions. Similar results were found for the remaining 4 sessions.

Fig. 5

 5 Fig. 5 shows the performance of a similar network for predicting the EMG signal of Abductor Pollicis Longus (AbPL) muscle.

Fig. 5 .

 5 Fig. 5. Performance of prediction for EMG activity of AbPL muscle as function of window size and number of cortical cells. Same conventions as in Fig. 4.

Fig. 6 .

 6 Fig. 6. The recorded (blue) and estimated (red) finger trajectories (top) and EMG activity of FDS and EDC for three successive trials.

Fig. 8 .

 8 Fig. 8. The performance comparison between the TDMLP-estimated and reproduced fingertip endpoint positions over 25 trials recorded during one session.

Acknowledgements

This study was in part funded through a grant by the 'Fondation pour la Recherche Médicale' (FRM), France. DBC20080713368.


Bruno Cessac

Juan-Carlos Vasquez

Hassan Nasser

Horacio Rostro-González

Thierry Viéville

Adrian Palacios

PARAMETRIC ESTIMATION OF SPIKE TRAIN STATISTICS BY GIBBS DISTRIBUTIONS : AN APPLICATION TO BIO-INSPIRED AND EXPERIMETAL DATA

Keywords: Spike trains statistics, Gibbs measures, Numerical implementations, MEA recordings

, in order to characterizing the statistics of multi-unit spike trains. We present this here with the aim to analyze and modeling synthetic data, especially bio-inspired simulated data e.g. from Virtual Retina [2], but also experimental data Multi-Electrode-Array(MEA) recordings from retina obtained by Adrian Palacios. We remark that Gibbs distribution allow us to estimate the spike statistics, given a design choice, but also to compare different models, thus answering comparative questions about the neural code.

Context.

Processing and encoding of information in neural dynamics is a very active research field [START_REF] Rieke | Spikes, Exploring the Neural Code[END_REF], although still much of the role of neural assemblies and their internal interactions remains unknown [START_REF] Pouget | Paying attention to correlated neural activity[END_REF]. The simultaneously recording of the activity of groups of neurons (up to several hundreds) over a dense configuration, supply today a critical database to unravel the role of specific neural assemblies.

A recent popular approach has been applied by [START_REF] Schneidman | Weak pairwise correlations imply string correlated network states in a neural population[END_REF]. They have shown that a model taking into account pairwise synchronizations between neurons in a small assembly (10-40 retinal ganglion cells) describes most (90%) of the correlation structure and of the mutual information of the block activity, and performs much better than a non homogeneous Poissonian model. Analogous results were presented the same year by [START_REF] Shlens | The structure of multi-neuron firing patterns in primate retina[END_REF]. The model used by both teams is based on a probability distribution known as the Gibbs distribution of the Ising model and comes from statistical physics. We develop here a numerical method allowing to parametrically estimate and compare, models for the statistics of simulated multi-cell-spike trains. Our method is not limited to firing rates models, pairwise synchronizations as [START_REF] Schneidman | Weak pairwise correlations imply string correlated network states in a neural population[END_REF][START_REF] Shlens | The structure of multi-neuron firing patterns in primate retina[END_REF][START_REF] Shlens | The structure of large-scale synchronized firing in primate retina[END_REF] or 1-step time pairwise correla-INRIA Cortex , UHP et Nancy 2 et CREA Ecole Polytechnique tions models as [START_REF] Marre | Prediction of spatiotemporal patterns of neural activity from pairwise correlations[END_REF], but deals with general form of Gibbs distributions, with parametric potentials corresponding to a spike n-uplets expansion with multi-units and multi-times terms. The method is exact (in the sense that is does not involve heuristic minimization techniques) since it is based on spectral properties of a transition matrix associated with a Markov chain.

Theory

Spike trains and Raster Plots

Lets denote by i = 1 . . . N the neuron index, and assume we have discretized the time in bins of size ∆t so t denotes the bin number, then to each membrane potential value corresponding to each bin V i (t) we associate a binary variable ω i (t) = Z(V i (t)), where Z(x) = 1 whenever x ≥ θ and Z(x) = 0 otherwise and θ is the firing threshold. The "spiking pattern" of the neural network at time t is the vector ω(t) = (ω i (t)) N i=1 : it tells us which neurons are firing at time t, (ω i (t) = 1) and which neurons are not firing at time t (ω i (t) = 0). We denote by ω t s the sequence or spike block ω(s) . . . ω(t). A bi-infinite sequence ω = {ω(t)} +∞ t=-∞ of spiking patterns is called a "raster plot". In practice, raster plots are obviously finite sequences of spiking patterns In the next sections we give our main results whose details can be found in [START_REF] Vasquez | Entropy-based parametric estimation of spike train statistics[END_REF].

Empirical averaging

As a starting point, one computes experimental averages allowing us to estimate the average value of some prescribed observable φ(ω), i.e. a function which associates to a raster plot a real number. Typical examples of observables are φ(ω) = ω i (0) in which case the average of φ is the firing rate of neuron i; φ(ω) = ω i (0)ω j (0) then the average of φ measures the probability of spike coincidence for neuron j and i; φ(ω) = ω i (τ )ω j (0) then the average of φ measures the probability of the event "neuron j fires and neuron i fires τ time step later" (or sooner according to sign(τ )).

The average of φ is computed, from the experimental raster plots. Let us assume that it is sufficient to focus on time average for a single raster plot, via the time-empirical average:

π (T ) ω (φ) = 1 T T t=1 φ(σ t ω). (1)
Note that in [START_REF] Vasquez | Entropy-based parametric estimation of spike train statistics[END_REF] we have used the shift σ t for the time evolution of the raster plot.

Inferring statistics with fixed observables

We are seeking a probability distribution ν which is compatible with these empirical average i.e. for all l = 1 . . . L,

ν(φ l ) = π (T ) ω (φ l ) = C l . (2)
This is a minimal, but insufficient requirement, since one can construct infinitely many probability distributions satisfying these constraints. One has therefore to add more structure. We propose two additional requirements: Stationarity and using Gibbs measures as parametric statistical models.

Stationarity The first one corresponds to assuming that data are extracted from a stationary dynamics i.e. that statistics are invariant under time translation, during the observation window.On mathematical grounds, this corresponds to requiring that ν is invariant. Indeed, by definition, an invariant measure ν is ergodic if for any typical raster plot (i.e. ν-almost every raster plot) ω, and for any (measurable) observable φ, the empirical average π

(T ) ω (φ) converges, as T → +∞ to the theoretical expectation of φ with respect to ν. Namely, lim T →+∞ π (T ) ω (φ) = ν(φ),
for ν-almost every ω. This means that the empirical average better and better approximates the expected value as T → +∞ Gibbs distributions The second requirement states that one only wants to take into account the constraints (2) without adding additional assumptions. Statistical physics naturally proposes a canonical way to construct such a statistical model: "Maximizing the statistical entropy under the constraints (2)" [START_REF] Jaynes | Information theory and statistical mechanics[END_REF]. In the present context, this amounts to solving a variational principle:

P [ψ] = sup ν∈m (inv) h [ν] + ν [ψ] , (3)
where m (inv) is the set of invariant measures and h the Kolmogorov-Sinai entropy or entropy rate. The term ψ is given by:

ψ = L l=1 λ l φ l , (4)
where the λ l 's are adjustable Lagrange multipliers. It is called a "potential. Typically, in the case considered here, the observables φ l and therefore the potential ψ depends on the raster plot over a finite time horizon R, i.e. ψ(ω) = ψ(ω(0), . . . , ω(R -1)). One speaks of "range-R potentials". They constitute Markovian approximations of more general (infinite range) potentials occurring e.g. in IF models with noise, where the degree of accuracy of the approximation can be mathematically controlled : the accuracy increases exponentially fast as R growths (see [START_REF] Cessac | A discrete time neural network model with spiking neurons ii. dynamics with noise[END_REF] for details).

A probability ν ψ which realizes the supremum (3), i.e.

P [ψ] = h [ν ψ ] + ν ψ [ψ] , (5)
is called, in this context, an "equilibrium state". Typically, equilibrium states are also Gibbs states or Gibbs distributions. The Gibbs distribution property means that one can find some constants P (ψ), c 1 , c 2 with 0 < c 1 ≤ 1 ≤ c 2 such that for all n ≥ 0 and for all ω ∈ X:

c 1 ≤ µ ψ ω 0 -n exp -(n + 1)P (ψ) + 0 k=-n ψ(T k ω) ≤ c 2 .
where T is the right shift over the set of infinite sequences X i.e. (T ω)(t) = ω(t -1), t ≤ 0. This means that the probability of (large) sequences behaves like

µ ψ ω 0 -n ∼ P 0 k=-n ψ(T k ω)
exp (n+1)P (ψ) . The term P (ψ), called the topological pressure in this context is the formal analog of a thermodynamic potential (free energy). It is a generating function for the cumulants of ψ.

Numerical Estimation of Gibbs Distributions

P (ψ) and µ ψ can be found by a spectral approach, being respectively the (unique) largest eigenvalue and related left eigenfunction of the Ruelle-Perron-Frobenius operator

L ψ f (ω) = ω ′ :σω ′ =ω ψ(ω ′ )f (ω ′ ), acting on C(X, Ê),
the set of continuous real functions on X.

Ruelle-Perron-Frobenius operator

The Ruelle-Perron-Frobenius (RPF) operator was introduced by Ruelle in [START_REF] Ruelle | Statistical mechanics of a one-dimensional lattice gas[END_REF]. In the present case of range-R potentials it acts on vectors v in H as:

(L ψ v) w = w ′ ∈2 N R e ψ w ′ v w ′ G w ′ ,w , (6)
where w = σw ′ , σ denoting again the time evolution. Thus, here L ψ is a positive 2 N R × 2 N R matrix (while it acts on functional spaces in the infinite range case). We note L ψ (w ′ , w) the entries of L ψ . The matrix G is called the grammar. It encodes the essential fact that the underlying dynamics is not able to produce all possible raster plots in its definition:

G w ′ ,w = 1, if the transition w ′ → w is admissible; 0, otherwise. (7
) Note therefore that L ψ is sparse.

In the present paper we make the assumption that the underlying (and hidden) dynamics is such that the L ψ matrix is primitive, i.e. ∃n > 0, s.t. ∀w, w ′ L n ψ (w ′ , w) > 0. This assumption holds for Integrate and Fire models with noise and is likely to hold for more general neural networks models where noise renders dynamics ergodic and mixing [START_REF] Cessac | A discrete time neural network model with spiking neurons ii. dynamics with noise[END_REF]. where . is the scalar product in H. The remaining part of the spectrum is located in a disk in the complex plane, of radius strictly lower than s ψ . Moreover, for all v in H,

1 s n ψ L n ψ v → b ψ b ψ . v, (8)
The Gibbs distribution is

ν ψ = b ψ b ψ , (9)
i.e. the probability of a spin block w of length R is

ν ψ (w) = b ψ ,w b ψ,w ,
where b ψ ,w is the w-th component of b ψ .

Computing averages

Since ν ψ [φ l ] = w ν ψ [w]φ l (w) one obtains using (9):

ν ψ [φ l ] = w∈H b ψ ,w φ l (w)b ψ,w . (10)
This provides a fast way to compute ν ψ [φ l ] and will be used to tune the λ l 's so that the theoretical averages ν ψ [φ l ] matches the empirical average C l according to eq. ( 16).

Entropy

An exact expression for the Kolmogorov-Sinai entropy can be readily obtained from eq. ( 5) giving:

h [ν ψ ] = log(s ψ ) - l λ * l ν ψ [φ l ] . (11)

Comparing several Gibbs statistical models

The choice of a potential (4), i.e. the choice of a set of observables, fixes a statistical model for the statistics of spike trains. Clearly, there are many choices of potentials and one needs to propose a criterion to compare them. The Kullback-Leibler divergence,

d(µ, ν) = lim sup n→∞ 1 n [ω] n-1 0 µ [ω] n-1 0 log   µ [ω] n-1 0 ν [ω] n-1 0   , (12)
where µ and ν are two invariant probability measures, provides some notion of asymmetric "distance" between µ and ν. Minimizing this divergence, i.e. the conditional entropy h[µ|ν], corresponds to minimizing "what is not explained in the empirical measure µ by the theoretical measure ν".

The computation of d(µ, ν) is numerically delicate but, in the present context, the following holds. For µ an invariant measure and ν ψ a Gibbs measure with a potential ψ, both defined on the same set of sequences Σ, one has [START_REF] Keller | Equilibrium States in Ergodic Theory[END_REF][START_REF] Chazottes | Pressure and Equilibrium States in Ergodic Theory, chapter Ergodic Theory[END_REF]:

d (µ, ν ψ ) = P [ψ] -µ(ψ) -h(µ).
This suggests us 2 to use this relation to compare different statistical models (corresponding to different potentials) by choosing the one which minimizes the quantity:

d(π (T ) ω , ν ψ ) = P [ψ] -π (T ) ω (ψ) -h(π (T ) ω ). ( 13
)
The advantage is that this quantity can be numerically estimated, since the topological pressure is known from the Ruelle-Perron-Frobenius theorem, while π

(T ) ω (ψ) is di- rectly computable. Since π (T )
ω is fixed by the experimental raster plot, h(π (T ) ω ) is fixed independently of the statistical model. Thus, comparing two statistical models corresponding to 2 potentials ψ 1 , ψ 2 amounts to comparing

P [ψ 1 ] -π (T ) ω (ψ 1 ) and P [ψ 2 ] -π (T ) ω (ψ 2 ). Introducing h [ψ] = P [ψ] -π (T ) ω (ψ) = P [ψ] - l λ l π (T ) ω (φ l ), (14)
the comparison of two statistical models, i.e. knowing if ψ 2 is significantly "better" that ψ 1 , reduces to verify:

h [ψ 2 ] ≪ h [ψ 1 ] , (15)
easily computable at the implementation level, as developed below. Note that h has the dimension of entropy. Since we compare entropies, which units are bits of information, defined in base 2, the previous comparison units is well-defined.

Estimating the potential parameters

The final step of the estimation procedure is to find the parameters λ such that the Gibbs measure fits at best with the empirical measure. We have discussed why minimizing ( 13) seems the best choice in this context. Since h(π

(T )
ω ) is a constant with respect to λ, it is equivalent to minimize h [ψ λ ] eq. ( 14) while

∂ h ∂λ l = ∂P [ψ λ ] ∂λ l -π (T ) ω (φ l ) = ν ψ (φ l ) -π (T ) ω (φ l ), (16)
where ν ψ (φ l ) is given by [START_REF] Cessac | A discrete time neural network model with spiking neurons ii. dynamics with noise[END_REF]. Equivalently, we are looking for a Gibbs distribution ν ψ such that

∂P [ψ λ ] ∂λ l = π (T ) ω (φ l ) which expresses that π (T ) ω
is tangent to P at ψ λ [START_REF] Keller | Equilibrium States in Ergodic Theory[END_REF].

Applications

This analysis opens up the possibility of developing efficients algorithms to estimate at best the statistic of spike 2 Although π (T ) ω is not invariant. This approximation becomes exact as T → +∞. trains from experimental data, using several guess potential and selecting the one which minimizes the KL divergence [START_REF] Vasquez | Entropy-based parametric estimation of spike train statistics[END_REF]. The idea is to start from a parametric form of potential (Ψ), of range R, and to compute the empirical average of all monomials φ l from data (say, an experimental raster ω (exp) ). Then, one adjust the parameters λ l by minimizing the KL divergence between the Gibbs measure µ ψ (R) λ and the empirical measure attached to the experimental raster ω (exp) .

The complete computational strategy and the description of the implemented numerical algorithms is described in [START_REF] Vasquez | Entropy-based parametric estimation of spike train statistics[END_REF] and is freely available as a C++ source at http://enas.gforge.inria.fr.

Preliminary data analysis

As a first step, we observe the outlook of the correlation in these different data types. To do so, we estimate the correlation between a couple of neurons for several various raster length. The correlation between two neuron is then:

C i,j = |µ(ω i (t)ω j (t)) -µ(ω i (t))µ(ω j (t))| (17)
Thus, C i,j is equal to how many time the two neurons i and j fired together minus the product of their firing rate.

The computed correlation between non-correlated data e.g between two neurons whose raster plots follow a Bernoulli distribution decreases to zero exponentially with a decay constant ∼ 1/ √ T as described by large deviations theory. Figure 1 using data obtained from virtual retina software ( [START_REF] Wohrer | Virtual retina: a biological retina model and simulator, with contrast gain control[END_REF]) shows the same behavior as in randomly generated data (Bernoulli distribution, figure not shown) as independently of distance between cells. We perform the same analysis for real data acquired with and MEA chip (60 channels) on retinal ganglion cell with no stimuli (spontaneous regime). Figure 2 shows results for correlation between several couples of neurons, at different distances proving the presence of effective (i.e. non-zero asymptotically) correlations

Detailed Results for data

We expect to provide detailed results in time for the conference. In particular, MEA recordings obtained by Adrien Palacios are currently analyzed within this approach.

Discussion

In this paper we have addressed the question of characterizing the spike train statistics of a group of neurons in the stationary case, thanks to the framework of thermodynamic formalism in ergodic theory (see [START_REF] Cessac | How gibbs distribution may naturally arise from synaptic adaptation mechanisms: a model based argumentation[END_REF] for the related theory). We have shown that the Jaynes method, based on an a priori choice of a "guess" potential, with finite range, amounts to approximate the exact probability distribution by the Gibbs distribution of a Markov chain. The degree of approximation can be controlled by the Kullback-Leibler divergence. We will show detailed results on different type of data, in time for the conference With respect to the state of the art, this method allows us to consider either more neurons with non-trivial statistics (e.g. beyond rate models and even models with correlation), thus targeting models with complex spike pattern. This method is in a sense the next step after Ising models, known as being able to represent a large but limited part of the encoded information (e.g. [START_REF] Schneidman | Weak pairwise correlations imply string correlated network states in a neural population[END_REF][START_REF] Mezard | Constraint satisfaction problems and neural networks: A statistical physics perspective[END_REF]). Another very important difference with respect to other current methods is that we perform the explicit variational optimization of a well defined quantity, i.e., the KL-divergence between the observed and estimated distributions. The method proposed here does not rely on Monte Carlo Markov Chain methods but on a spectral computation based on the RPF operator, providing exact formula, while the spectral characteristics are easily obtained from standard numerical methods.

A step further, the non-trivial but very precious virtue of the method is that it allows us to efficiently compare models. We thus not only estimate the optimal parameters of a model, but can also determine among a set of models which model is the most pertinent. This means for instance that we can determine if either only rates, or rates and correlations matters, for a given piece of data. Another example is to detect if a given spike pattern is significant, with respect to a model not taking this pattern into account. These elements push the state of the art regarding statistical analysis of spike train a step further.

One weakness of the present work is that it only considers stationary dynamics, where e.g. the external current I i is independent of time. However, real neural systems are submitted to non static stimuli, and transients play a crucial role, but considering non stationary requires to handle time dependent Gibbs measures. In the realm of ergodic theory applied to non equilibrium statistical physics, Ruelle has introduced the notion of time-dependent SRB measure [START_REF] Ruelle | Smooth dynamics and new theoretical ideas in nonequilibrium statistical mechanics[END_REF]. A similar approach could be used here, at least formally.

  Then L ψ obeys the Perron-Frobenius theorem 1 : Theorem 1 L ψ has a unique maximal eigenvalue s ψ = e P [ψ] associated with a right eigenvector b ψ and a left eigenvector b ψ such that L ψ b ψ = s ψ b ψ , and b ψ L ψ = s ψ b ψ . Those vectors can be chosen such that b ψ .b ψ = 1

Figure 1 .

 1 Figure 1. Correlation of spiking activity vs the sample size, computed from data obtained from Virtual Retina ([2]) for several couples of neurons situated at different distances. The distance is measured in pixel and the virtual retina is configured such with 1 neuron/pixel.

Figure 2 .

 2 Figure 2. Correlation between ganglion cells in real MEA measurements. The figure shows the correlation for three couples of neurons, at several distances. The time scale corresponds to the first 30 sec. of the measured neurons .

This theorem has been generalized by Ruelle to infinite range potentials under some regularity conditions[START_REF] Ruelle | Statistical Mechanics: Rigorous results[END_REF][START_REF] Ruelle | Thermodynamic formalism[END_REF].