Salima El

Florian Pelgrin

PARAMETRIC INFERENCE OF HIDDEN DISCRETE-TIME DIFFUSION PROCESSES BY DECONVOLUTION

We study a parametric approach for hidden discrete-time diffusion models based on contrast minimization and deconvolution. This approach leads to estimate a large class of stochastic models with nonlinear drift and nonlinear diffusion. It can be applied, for example, for ecological and financial state space models. After proving consistency and asymptotic normality of the estimator, leading to asymptotic confidence intervals, we provide a thorough numerical study, which compares many classical methods used in practice (Monte Carlo Expectation Maximization Likelihood estimator and Bayesian estimators) to estimate stochastic volatility models. We prove that our estimator clearly outperforms the Maximum Likelihood Estimator in term of computing time, but also most of the other methods. We also show that this contrast method is the most stable and also does not need any tuning parameter.

Introduction

This paper is motivated by the parametric estimation of hidden stochastic models of the form:

Y i = X i + ε i X i+1 = b θ 0 (X i ) + σ θ 0 (X i )η i+1 , (1)
where one observes Y 1 ,• • • ,Y n , and where the random variables ε i , η i and X i are unobserved. Notably (X i ) i≥0 is a strictly stationary, ergodic process that depends on two measurable functions b θ 0 and σ θ 0 and its stationary density is f θ 0 , where θ 0 belongs to Θ ⊂ R p . The functions b θ 0 , σ θ 0 and f θ 0 are known up to a finite dimensional parameter, θ 0 , and the dependence with respect to θ 0 is not required to be the same in b θ 0 and σ θ 0 . Finally, the innovations (η i ) i≥0 and the errors (ε i ) i≥0 are independent and identically distributed (i.i.d.) random variables, the distribution of the innovations being known for identifiability of the model.

In this work, we propose to estimate the parameters of the two functions b θ 0 and σ θ 0 driving the dynamics of the hidden variables (X i ) i≥0 . The principle of the estimation method goes as follows.

Taking that the stationary density, f θ 0 , is known up to the finite dimensional parameter θ 0 , our Mestimator consists in optimizing a contrast function that exploits a Fourier deconvolution strategy in a parametric framework. In so doing, we exploit a "Nadaraya-Watson approach" in the sense that we estimate b θ 0 (respectively, b 2 θ 0 + σ 2 θ 0 ) as ratio of an estimator of l θ 0 = b θ 0 f θ 0 (respectively, l θ 0 = (b 2 θ 0 + σ 2 θ 0 )f θ 0 ) and an estimator of f θ 0 . Notably we provide an analytical expression of the contrast function for a well-known example and characterize their main properties. Moreover we show that this deconvolution-based estimator is consistent and asymptotically normally distributed for αmixing processes which leads to obtain confidence intervals in practice for many processes. Finally, our Monte Carlo simulations show that our approach gives good results and is fast computing. All the results are illustrated on the famous stochastic volatility model with discrete time version of CIR 1 (Cox Ingersoll Ross, see [START_REF] Cox | A theory of the term structure of interest rates[END_REF]) process for the volatility and are compared with many others methods used in the literature to estimate this model (Monte Carlo Expectation Maximisation, Sequential Monte Carlo).

Our approach extends the previous work [El Kolei, 2013] where parametric estimation of models of type of (1) is handled for constant volatility function (∀x, σ θ 0 (x) = σ θ 0 ) and where the estimator proposed by the author is not adapted for stochastic process with nonlinear diffusion. As in the previous work [El Kolei, 2013], this approach is the extension to the parametric framework of the work [START_REF] Comte | Adaptive estimation of the dynamics of a discrete time stochastic volatility model[END_REF] where the authors propose a non-parametric estimation procedure in the case of discrete time stochastic models of the form of (1). 1 Our aim consists in showing how their procedure can be extended to a parametric framework and further by obtaining confidence intervals which are useful in practice.

Applications. This class of parametric models includes, among others, the autoregressive model with measurement errors, the autoregressive stochastic volatility model ( [Taylor, 2005]), the discrete time versions of well-known diffusion processes in finance ( [START_REF] Hull | Pricing interest-rate-derivative securities[END_REF]], [Heston, 1993]) and some families of stochastic processes: Vasicek, CIR, modified CIR and hyperbolic processes (see [START_REF] Cox | A theory of the term structure of interest rates[END_REF] and [START_REF] Genon-Catalot | Parameter estimation for discretely observed stochastic volatility models[END_REF]).

Here, we focus on a stochastic volatility model of the form:

R i+1 = exp X i+1 2 ξ i+1 , X i+1 = X i + κ(µ -X i )∆ + σ √ ∆X i η i+1 (2)
where ξ i and η i are centered gaussian random variables and ∆ the sampling interval. Hence, the unobserved variance process X i is driven by a mean reverting stochastic process which was introduced in [START_REF] Cox | A theory of the term structure of interest rates[END_REF] to model the short term interest rates.

By applying a log-transformation

Y i+1 = log(R 2 i+1 )-E[log(ξ 2 i+1 )] and ε i+1 = log(ξ 2 i+1 )-E[log(ξ 2 i+1 )
], the SV model is a particular version of (1) since it can be written as

Y i = X i + ε i X i+1 = X i + κ(µ -X i )∆ + σ √ ∆X i η i+1 (3)
where ε i follows a log chi-squared distribution.

From a practical point of view, the observed component (Y i ) i≥1 stands for the log-return of an asset price while the unobserved component stands for the volatility of this asset. The parameter κ is the positive mean reverting parameter, µ is the positive long run parameter and σ the positive volatility of the stochastic volatility process (X i ) i≥1 .

Organization of this paper. The paper is organized as follows. Section 2 presents the notations and the model assumptions. Section 3 defines the deconvolution-based M-estimator and states all of the theoretical properties. Some Monte Carlo simulations are discussed in Section 4 and some concluding remarks are provided in the last section. All the proofs can be found in Appendix 5.

1 See also Comte et Taupin in [START_REF] Comte | Adaptive estimation in a nonparametric regression model with errors-in-variables[END_REF].

General setting and assumptions

In this section, we introduce some preliminary main notations and provide the assumptions of the model (1).

Notations

Subsequently, for any function v : R → R, we denote by v * the Fourier transform of the function v: v * (t) = e itx v(x)dx, by ||v|| 2 its L 2 (R)-norm, ||v|| ∞ its supremum norm, •, • stands for the scalar product in L 2 (R) and " " for the usual convolution product. Moreover, for any integrable and squareintegrable functions u, u 1 , and u 2 : we have (u

* ) * (x) = 2πu(-x)and u 1 , u 2 = 1 2π u * 1 , u * 2 . Finally, A denotes the Euclidean norm of a matrix A, Y i = (Y i , Y i+1
) and y i = (y i , y i+1 ), P n (respectively, P) the empirical (respectively, theoretical) expectation, that is, for any stochastic variable:

P n (X) = 1 n n i=1 X i (respectively, P(X) = E[X]
). Regarding the partial derivatives, for any function h θ , ∇ θ h θ is the vector of the partial derivatives of h θ with respect to (w.r.t) θ and ∇ 2 θ h θ is the Hessian matrix of h θ w.r.t θ.

Assumptions

We consider the hidden discrete-time diffusion model (1). The assumptions are the following. A0 θ 0 belongs to the interior Θ 0 of a compact set Θ, θ 0 ∈ Θ ⊂ R p . A1 The errors (ε i ) i≥0 are independent and identically distributed centered random variables with finit variance, E ε 2 1 = s 2 . The density of ε 1 , f ε , belongs to L 2 (R), and for all x ∈ R, f * ε (x) = 0. A2 The innovations (η i ) i≥0 are independent and identically distributed centered random variables with unit variance E η 2 1 = 1 and E η 3 1 = 0. A3 The X i 's are strictly stationary and ergodic with invariant density f θ 0 . A4 The sequences (X i ) i≥0 and (ε i ) i≥0 are independent. The sequence (ε i ) i≥0 and (η i ) i≥0 are independent. A5 On Θ 0 , the functions θ → b θ and θ → σ θ admit continuous derivatives with respect to θ up to order 2. A6 The function to estimate ), is twice continuously differentiable w.r.t θ ∈ Θ for any x and measurable w.r.t x for all θ in Θ. Each element of ∇ θ l θ and ∇ 2 θ l θ belongs to L 1 (R) ∩ L 2 (R). A7 The application θ → Pm θ admits a unique minimum and its Hessian matrix, denoted by V θ , is non-singular in θ 0 .

l θ := b 2 θ + σ 2 θ f θ belongs to L 1 (R) ∩ L 2 (R
The compactness assumption A0 might be relaxed by assuming that θ 0 is an element of the interior of a convex parameter space Θ ∈ R p . In this case, the statistical properties of the M-estimator can be proved in the light of convex optimization arguments. Assumptions A1-A3 are quite standard when considering estimation for convolution models. On the other hand, Assumption A3 implies that if (X i ) i≥0 is an ergodic process then (Y i ) i≥0 is stationary and ergodic since it is the sum of an ergodic process and an i.i.d. noise process [START_REF] Dedecker | Weak dependence: with examples and applications[END_REF]). Consequently Y i = (Y i , Y i+1 ) inherits the ergodicity property. According to Assumption A4 the unknown density g θ 0 of the Y i 's is defined to be

f θ 0 f ε . It turns out that g * θ 0 = f * θ 0 f * ε and thus f * θ 0 = g * θ 0 /f * ε .
Assumption A5 ensures some smoothness for the drift and diffusion functions. Assumption A6 is also quite usual in the literature and serves for the construction and for asymptotic properties of our estimator.

3 Parametric deconvolution estimator

The contrast function

Definition 1 (Theoretical and empirical contrast functions). For any square integrable real function v, we set

u v (x) = 1 2π v * (-x) f * ε (x)
,

where we recall that f * ε is the Fourier transform of the density of the observation noise. Let ϕ : R → R be given by

(1) ϕ : x → x, if σ θ 0 is a constant function of the hidden variable (2) ϕ : x → x 2 -s 2 , if σ θ 0 is not a constant function of the hidden variable where s 2 = E[ε 2
1 ], and let us define the mapping m as m

• (•) : (θ, y i ) ∈ (Θ × R 2 ) → m θ (y i ) = ||l θ || 2 -2ϕ (y i+1 ) u * l θ (y i ).
Then, under Assumptions A1 up to A7, the contrast function is defined by:

E [m θ (Y 1 )] := l θ 2 -2E ϕ(Y 2 )u * l θ (Y 1 ) , (4)
and its empirical couterparts is given by

P n m θ = 1 n n i=1 m θ (Y i ). ( 5
)
Remark 1. As said in the introduction, the case where the diffusion function σ θ 0 is a constant function of the hidden variable has already been studied in [El Kolei, 2013]. Therefore, from now on, we focus on the case 2 in Definition 1 and we refer to the aforementioned paper for the case 1 in Definition 1.

Definition 2 (Minimal contrast estimator). Suppose that Assumptions A0-A7 hold true then, the minimum-contrast estimator θ n is defined as any solution of

θ n = arg min θ∈Θ P n m θ . (6)
The existence of our estimator can be deduce from regularity properties of the function l θ and compactness argument of the parameter space. See Appendix 5.2.1.

Remark 2. In this paper we consider the situation in which the observation noise variance is known. This assumption, which is often not satisfied in practice, is necessary for the identifiability of the model and so is a standard assumption for state-space models given in (1). There is some restrictions on the distribution of the innovations in the Nadaraya-Watson approach. It is known that the rate of convergence for estimating the function l θ is related to the rate of decreasing of f * ε . In particular, the smoother f ε is, the slower the rate of convergence for estimating is. This rate of convergence can be improved by assuming some additional regularity conditions on l θ (see [START_REF] Comte | Adaptive estimation of the dynamics of a discrete time stochastic volatility model[END_REF] and [START_REF] Comte | Penalized contrast estimator for adaptive density deconvolution[END_REF]).

Remark 3. Let us explain the choice of the contrast function and how the strategy of deconvolution works. The convergence of P n m θ to Pm θ = E [m θ (Y 1 )] as n tends to the infinity is provided by the Ergodicity Theorem. Moreover, the limit E [m θ (Y 1 )] of the contrast function can be explicitly computed. Using (1) and Assumptions A1-A3, standard computations (see Appendix 5.1) lead to

E [m θ (Y 1 )] = l θ 2 -2 l θ , l θ 0 = l θ -l θ 0 2 -l θ 0 2 ,
which is, obviously, minimal at point θ = θ 0 .

Asymptotic properties

In this section we first show that our estimator is weakly consistent and asymptotically normally distributed for mixing processes. To this aim, we further assume that for ϕ defined in (2) of Definition 1 the following assumptions hold true:

A8 (Local dominance): E sup θ∈Θ ϕ(Y i+1 )u * l θ (Y i ) < ∞. A9 (Moment condition): For some δ > 0, E ϕ(Y i+1 )u * ∇ θ l θ (Y i ) 2+δ < ∞.
A10 (Hessian Local dominance): For some neighbourhood U of θ 0 :

E sup θ∈U ϕ(Y i+1 )u * ∇ 2 θ l θ (Y i ) < ∞ .

Asymptotic properties of the estimator: consistency and normality

The first result regards the (weak) consistency of our estimator.

Theorem 1. Consider the model (1) under the assumptions A0-A8, the estimator θn defined by ( 6) is weakly consistent:

θ n -→ θ 0 as n → ∞ in P θ 0 -probability.
Sketch of proof. The main idea for proving the consistency of a M-estimator comes from the following observation: if P n m θ converges to Pm θ in probability, and if the true parameter solves the limit minimization problem, then, the limit of the argminimum θ n is θ 0 . By using an argument of uniform convergence in probability and by compactness of the parameter space, we show that the argminimum of the limit is the limit of the argminimum. A standard method to prove the uniform convergence is to use the Uniform Law of Large Numbers (see Lemma 1 in Appendix 5.2.2). Combining these arguments with the dominance argument (A8) give the consistency of our estimator, and then, the Theorem 1. For further details see Appendix 5.2.2.

The second result states our estimator is √ n-consistent and asymptotically normally distributed. Besides, we give in Corollary 1 the different terms of the asymptotic variance-covariance matrix. For the CLT, we need some mixing properties (we refer the reader to [START_REF] Dedecker | Weak dependence: with examples and applications[END_REF] for a complete review of mixing processes). Hence, in the following, we further assume that:

A11 The stochastic process X i is α-mixing.
Theorem 2. Consider the model (1) under the assumptions A0-A7, and suppose that the conditions A8-A11 hold true. Then θn defined by ( 6) is a √ n-consistent estimator of θ 0 which satisfies:

√ n( θ n -θ 0 ) L → N (0, Σ(θ 0 )) .
Sketch of proof. The asymptotic normality follows essentially from Central Limit Theorem for mixing processes (see [Jones, 2004]). Thanks to the consistency, the proof is based on a moment condition of the Jacobian vector of the function m θ (y) and on a local dominance condition of its Hessian matrix.

For further details, see Appendix 5.2.3.

The following corollary gives an expression of the variance-covariance matrix Σ(θ 0 ) of Theorem 2 for the practical implementation:

Corollary 1. Under our assumptions, the variance-covariance matrix Σ(θ 0 ) is given by:

Σ(θ 0 ) = V -1 θ 0 Ω(θ 0 ) = V -1 θ 0   Ω 0 (θ 0 ) + 2 +∞ j=2 Ω j-1 (θ 0 )   V -1 θ 0 , with Ω 0 (θ 0 ) = 4 E ϕ(Y 2 )u * ∇ θ l θ (Y 1 ) ϕ(Y 2 )u * ∇ θ l θ (Y 1 ) -E b 2 θ 0 (X 1 ) + σ 2 θ 0 (X 1 ) ∇ θ l θ (X 1 ) E b 2 θ 0 (X 1 ) + σ 2 θ 0 (X 1 ) ∇ θ l θ (X 1 ) Ω j-1 (θ 0 ) = 4 E b 2 θ 0 (X 1 ) + σ 2 θ 0 (X 1 ) ∇ θ l θ (X 1 ) b 2 θ 0 (X j ) + σ 2 θ 0 (X j ) ∇ θ l θ (X j ) -E b 2 θ 0 (X 1 ) + σ 2 θ 0 (X 1 ) ∇ θ l θ (X 1 ) E b 2 θ 0 (X 1 ) + σ 2 θ 0 (X 1 ) ∇ θ l θ (X 1 ) .
and the gradient ∇ θ l θ is taken at point θ = θ 0 ; furthermore, the Hessian matrix V θ 0 is given by:

[V θ 0 ] j,k 1≤j,k≤r = 2 ∂l θ ∂θ k , ∂l θ ∂θ j j,k at point θ = θ 0 .
Proof. See Appendix 5.2.4 for further details.

4

Application: the CIR process

We consider the following stochastic volatility model

Y i = X i + ε i X i+1 = X i + κ(µ -X i )∆ + σ √ ∆X i η i+1 (7)
where ε i follows a log chi-squared distribution and η i a gaussian distribution. This can be also seen as a discrete version of the so called Heston model with independent noise between the log-returns Y i and its volatility X i . Here, we further assume that the variance process X i is greater than zero. To ensure this condition, we make the following assumption:

F a := 2κθ σ 2 ≥ 1 and c := 2κ σ 2 > 0, which is known as the Feller condition (see [START_REF] Cox | A theory of the term structure of interest rates[END_REF]) and implies that the variance process X i is ergodic and ρ-mixing. Furthermore, the stationary distribution f θ for this process is the gamma distribution γ(a, c) (see [START_REF] Genon-Catalot | Parameter estimation for discretely observed stochastic volatility models[END_REF]).

Minimum contrast estimator

In this case, the functions b • , σ • and l • are given by:

b θ (x) = (1 -κ)x + κθ, σ θ (x) = σ √ x and l θ = b 2 θ (x) + σ 2 θ (x) γ(a, c) ∀x ∈ R * + with θ = (κ, µ, σ).
Using the Fourier transform of the Gamma and the log chi-squared density, we have

f * θ (x) = 1 - ix c -a and f * ε (x) = 1 √ π 2 ix Γ( 1 2 + ix) exp(-iCx),
with C the expectation of the logarithm of a chi-squared random variable, i.e. C = -1.27 (see [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]] and Appendix 5.3 for the expression of the Fourier transform). Next, the Fourier transform of l θ is given by

l * θ (x) = -α 1 -a c 2 (a + 1) 1 - ix c -a-2 + iα 2 a c 1 - ix c -a-1 + α 3 1 - ix c -a
.

with

α 1 = (1 -κ) 2 , α 2 = 2 (1 -κ) κθ + σ 2 , α 3 = (κθ) 2 .
Finally, the L 2 -norm of l θ is given by:

l θ 2 2 = α 2 1 2 -(2a+3) c -3 Γ(2a + 3) Γ 2 (a) + 2α 1 α 2 2 -(2a+2) c -2 Γ(2a + 2) Γ 2 (a) + 2α 1 α 3 + α 2 2 2 -(2a+1) c -1 Γ(2a + 1) Γ 2 (a) + α 2 α 3 2 -(2a) Γ(2a) Γ 2 (a) +α 2 3 2 -2a+1 c Γ(2a -1) Γ 2 (a) .
where Γ corresponds to the Gamma function given by Γ(z) = R + t z-1 exp(-t)dt.

Proof. See Appendix 5.3.

Hence, the M-estimator solves:

θ n = arg min θ∈Θ l θ 2 2 - 2 n n i=1 Y i+1 u * l θ (Y i ) (8)
where:

u l θ (y) = √ π     α 1 a c 2 (a + 1) 1 -iy c -a-2 + iα 2 a c 1 -iy c -a-1 + α 3 1 -iy c -a 2 iβy Γ( 1 2 + iβy) exp(-i Cy)    

Others methods

Particle filters : EKF, APFS, APF and KSAPF. For the comparison with our contrast estimator given in (8), we use the following methods: the Extended Kalman Filter (EKF), the Auxiliary Particle filter (APF), the Auxiliary Particle filter with static parameter (APFS) and the Kernel Smoothing Auxiliary Particle filter (KSAPF). We refer the reader to [START_REF] Reif | Stochastic stability of the discrete-time extended Kalman filter[END_REF], [START_REF] Pitt | Filtering via simulation: auxiliary particle filters[END_REF], [START_REF] Doucet | [END_REF] and [START_REF] Liu | Combined parameter and state estimation in simulation-based filtering[END_REF]] for a complete revue of these methods.

In order to estimate the parameters with these methods, for the EKF, APF and KSAPF estimators we use the Kitagawa and al.'s approach (see [START_REF] Doucet | [END_REF] chapter 10 p.189) in which the parameters are supposed time-varying: θ i+1 = θ i + G i+1 where G i+1 is a centered Gaussian random with a variance matrix Q supposed to be known. Hence, we consider the augmented state vector Xi+1 = (X i+1 , θ i+1 ) where X i+1 is the hidden state variable and θ i+1 the unknown vector of parameters. Furthermore, for the numerical part we call APFS an Auxiliary Particle filter without time-varying parameters.

The MCEM. From a theoretical point of view, the MLE is asymptotically efficient. However, in practice since the states (X 1 • • • , X n ) are unobservable and since the model ( 7) is non Gaussian, the likelihood is untractable. We have to use numerical methods to approximate it. In this section, we illustrate the MCEM estimator which consists in approximating the likelihood and applying the Expectation-Maximisation algorithm introduced by Dempster [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] to find the parameter θ.

Numerical Results

In this section we present some Monte Carlo simulations using the model ( 7). For the analysis we consider the following "true parameter" θ 0 = (κ 0 , µ 0 , σ 2 0 ) = (4, 0.03, 0.4) which is consistent with empirical applications of daily data (see [Do, 2005]). Thus, we have sampled the trajectory of the X i , and conditionally to the trajectory, we have sampled the variables Y i with a variance noise s 2 ε = 0.1 2 The numerical illustration goes as follows: we work with a number of observations n equal to 1000, we first compare all the methods proposed in term of computing time, then we run N = 100 estimates for each method and we compare the performance of our estimator with others methods by computing the Mean Square Error (MSE) defined as:

2 For the simulation (see Appendix 5.3) we take f

* ε (x) = 1 √ π 2 iβx Γ( 1 2 + iβx) exp(-i Cx) with β = 2s 2 ε /π 2 and C = βC where C = -1.27. MSE = 1 N p i=1 N j=1 ( θi j -θ i 0 ) 2 , (9)
where p corresponds to the dimension of the vector of parameters.

Then, we illustrate the statistical properties of our contrast estimator, by computing the confidence intervals for different number of observations. Finally, we study the influence of the signal noise ratio since the performance of our estimator depends on the regularity of l * θ and f * ε (see Remark 2).

For particles methods, we take a number of particles M equal to 5000. Note that for the Bayesian procedure (APF, APFS and KSAPF) we need a prior on θ, and this only at the first step. The prior for θ is taken to be the Uniform law and conditionally to θ the distribution of X 1 is its stationary law:

p(θ) = U(3, 5) × U(0.02, 0.04) × U(0.3, 0.5) f θ (X 1 ) = γ (θ)
For the KSAPF, we take a bandwidth h = 0.1 and for the APF we take a matrix Q = [10 -3 , 0.1.10 -4 , 10 -4 ]I 3 with I 3 the identity matrix in R 3 (see Section 4.3.1 for the definition of the matrix Q). This comparison illustrates the numerical complexity of the MCEM. Therefore, in the following, we only compare our contrast estimator with the particles estimators.

Computing time

Parameter estimates

We illustrate by boxplots the different estimates (see Figures [1] up to [3]). In Table [2], we compute the MSE defined in (9) for each method and with a number of MC equal to N = 100 and the CPU for a number of observations n = 1000 (see Table [1]).

We note that for all parameters, the EKF estimator is very bad since the stochastic volatility model ( 7) is strongly nonlinear, and its corresponding boxplots have the largest dispersion meaning that this filter is not stable and not appropriated to estimate this model. Among particle filters, the KSAPF and the APF are the best estimators although the dispersion is huge for the mean reversion parameter κ and the volatility parameter σ. Besides, the APFS is less efficient than the others particle filters. Our estimator is stable and performs the others when one compare the MSE. From a computational point of view, all particles filters have an equivalent CPU. Our contrast estimator is fast and its implementation is straightforward and the MSE is the smallest (see Table [1]).

Confidence Interval of the minimum contrast estimator

To illustrate the statistical properties of our contrast estimator, we compute the confidence intervals with the confidence level 1 -α equal to 0.95 for N = 100 estimators. The coverage corresponds to the number of times for which the true parameter θ i 0 , i = 1, • • • , p belongs to the confidence interval. The results are illustrated in Figure [4]. We note that the coverage converges to 95% for a small number of observations and as expected, the confidence interval decreases with the number of observations. Note that of course a MLE confidence interval would be smaller since the MLE is efficient but the corresponding computing time would be huge (see Table [1]).

Ratio signal-noise for the contrast estimator

We denote by r = [s 2 ε /σ 2 ] the ratio signal-noise and in Table (3) we compare the MSE for different r and different number of observations n for the contrast estimator. We note that the MSE decreases with the number of observations and is smaller for small ratio-signal-noise. As explained in Section 3.1 and see [START_REF] Comte | Adaptive estimation of the dynamics of a discrete time stochastic volatility model[END_REF] for more details, the rate of convergence of our approach depends on the regularity of the noise density f ε . And, in particular, the smoother the noises are, the slower the rate of convergence is. For the CIR model, the density of the noises and the function l θ are ordinary smooth, so we are in a favourable case.

Summary and Conclusions

In this paper we have proposed a new method to estimate hidden nonlinear diffusion process. This method is based on a deconvolution strategy and leads to consistent and asymptotically normal estimator. We have numerically studied the performance of our estimator for the CIR process widely used in many domains and we were able to construct confidence interval (see Figure [4]). As the boxplots [1] up to [3] show, only Contrast, APF, and KSAPF estimators are comparable. Indeed EKF and APFS estimators are biased and their MSE are bad, especially for the EKF method since the CIR process is nonlinear. Furthermore, if one compares the MSE of the particle filters, the KSAPF estimator is the best method. Among particles filters, it is clearly known that the APFS is less efficient than the APF filter since the parameters are not time-varying and so the only randomness is made at the first step by the prior law and not in each propagation step. Then, the Contrast, APF, and KSAPF methods lead to unbiased and not so much varying estimator. We emphasize that our estimator performs the others in a MSE aspect (see Table [1]). Most importantly, our estimator can be constructed without any arbitrary parameters choice, is straightforward to implement, fast and allows to construct confidence interval.

Appendix

The contrast function

The procedure: The limit E [m θ (Y 1 )] of the contrast function can be explicitly computed. Using (1) and Assumptions A1-A3, we obtain:

E (Y 2 2 -s 2 ε )u * l θ (Y 1 ) = E X 2 2 + 2X 2 ε 2 + ε 2 2 -s 2 ε u * l θ (Y 1 ) = E X 2 2 u * l θ (Y 1 ) by assumption A1 = E b 2 θ 0 (X 1 ) + σ 2 θ 0 (X 1 )η 2 2 + 2b θ 0 (X 1 )σ θ 0 (X 1 )η 2 u * l θ (Y 1 ) = E b 2 θ 0 (X 1 ) + σ 2 θ 0 (X 1 ) u * l θ (Y 1
) by assumption A2, Using Fubini's Theorem and (1), it follows that:

E b 2 θ 0 (X 1 ) + σ 2 θ 0 (X 1 ) u * l θ (Y 1 ) = E b 2 θ 0 (X 1 ) + σ 2 θ 0 (X 1 ) e iY 1 y u l θ (z)dz = E b 2 θ 0 (X 1 ) + σ 2 θ 0 (X 1 ) 1 2π 1 f * ε (z) e iY 1 z (l θ (-z)) * dy = 1 2π E b 2 θ 0 (X 1 ) + σ 2 θ 0 (X 1 ) e i(X 1 +ε 1 )z 1 f * ε (z) (l θ (-z)) * dz = 1 2π E e iε 1 z f * ε (z) E b 2 θ 0 (X 1 ) + σ 2 θ 0 (X 1 ) e iX 1 z (l θ (-z)) * dy = 1 2π E b 2 θ 0 (X 1 ) + σ 2 θ 0 (X 1 ) e iX 1 z (l θ (-z)) * dz = 1 2π E b 2 θ 0 (X 1 ) + σ 2 θ 0 (X 1 ) ((l θ (-X 1 )) * ) * = E b 2 θ 0 (X 1 ) + σ 2 θ 0 (X 1 ) l θ (X 1 ) . = b 2 θ 0 (x) + σ 2 θ 0 (x) f θ 0 (x) b 2 θ (x) + σ 2 θ (x) f θ (x)dx = l θ , l θ 0 .
Then,

E [m θ (Y 1 )] = l θ 2 -2 l θ , l θ 0 = l θ -l θ 0 2 -l θ 0 2 .

Proofs

For the reader convenience we split the proof of Theorems 1 and 2 into three parts: in Subsection 5.2.1, we give the proof of the existence of our contrast estimator defined in (6). In Subsection 5.2.2, we prove the consistency, that is, the Theorem 1. Then, we prove the asymptotic normality of our estimator in Subsection 5.2.3, that is, the Theorem 2. The Subsection 5.2.4 is devoted to Corollary 1.

Recall from Remark 1 that we only made the proof for the function ϕ defined by (2) in Definition 1 and we refer to [El Kolei, 2013] for the proof in the case (1) of Definition 1.

Proof of the existence and measurability of the M-Estimator

By assumption, the function θ → l θ 2 2 is continuous. Moreover, l * θ and then u *

l θ (x) = 1 2π e ixy l * θ (-y) f * ε (y) dy are continuous w.r.t θ. In particular, the function m θ (y i ) = l θ 2 2 -2ϕ(y i+1 )u * l θ (y i ) is continuous w.r.t θ, for ϕ : x → x 2 -s 2 ε . Hence, the function P n m θ = 1 n n i=1 m θ (Y i ) is continuous w.r.t.
θ belonging to the compact subset Θ. So, there exists θ that belongs to Θ such that:

inf θ∈Θ P n m θ = P n m θ.

Proof of the Consistency

For the consistency of our estimator, we need to use the uniform convergence given in the following Lemma. Let us consider the following quantities:

P n h θ = 1 n n i=1 h θ (Y i ); P n S θ = 1 n n i=1 ∇ θ h θ (Y i ) and P n H θ = 1 n n i=1 ∇ 2 θ h θ (Y i )
where h θ is real function from Θ × Y with value in R.

Lemma 1. Uniform Law of Large Numbers (ULLN)(see [START_REF] Newey | Large sample estimation and hypothesis testing[END_REF] for the proof.)

Let (Y i ) be an ergodic stationary process and suppose that:

(1) For all y ∈ R, θ → h θ (y) is continuous and for all θ ∈ Θ, y → h θ (y) is measurable.

(2) There exists a function (called the dominating function) s : R → R such that for all θ ∈ Θ |h θ (y)| ≤ s(y) and E[s(Y 1 )] < ∞. Then:

sup θ∈Θ |P n h θ -Ph θ | → 0 in probability as n → ∞.
Moreover, Ph θ is a continuous function of θ.

By assumption for all x, θ → l θ (x) is continuous and for all θ, x → l θ (x) is measurable which implies the continuity and the measurability of the function P n m θ on the compact subset Θ. Furthermore, the local dominance assumption (A8)

implies that E [sup θ∈Θ |m θ (Y i )|] is finite. Indeed, |m θ (y i )| = l θ 2 2 -2ϕ(y i+1 )u * l θ (y i ) ≤ l θ 2 2 + 2 ϕ(y i+1 )u * l θ (y i ) .
with ϕ the function defined in (2) in Definition 1.

As l θ 2 2 is continuous on the compact subset Θ, sup θ∈Θ l θ 2 2 is finite. Therefore, E [sup θ∈Θ |m θ (Y i )|] is finite if E sup θ∈Θ ϕ(Y i+1 )u * l θ (Y i ) is finite.
Lemma 1 gives us the uniform convergence in probability of the contrast function: for any ε > 0,

lim n→+∞ P sup θ∈Θ |P n m θ -Pm θ | ≤ ε = 1.
Combining the uniform convergence with Theorem 2.1 p. 2121 chapter 36 in [START_REF] Hansen | robert f. engle and daniel l. mcfadden[END_REF] yields the weak (convergence in probability) consistency of the estimator.

Proof of the asymptotic normality

Consider the model (1) under the assumptions A0-A7. The proof of the asymptotic normality results from assumptions A8-A11 is a straighforward application of [Hayashi, 2000] (Propostion 7.8. p. 472 and [Jones, 2004]). Furthermore, for the CLT, we need to recall some mixing properties (we refer the reader to [START_REF] Dedecker | Weak dependence: with examples and applications[END_REF] for a complete revue of mixing processes).

Let K θ (x, dy) be a Markov transition kernel on a general space (X , B(X ))) and let K n θ (x, dy) denotes the n step Markov transition corresponding to K θ . Then, for k ∈ N, x ∈ X and a measurable set A:

K n θ (x, dy) = P θ (X n+k ∈ A|X k = x)
Let M : R → R be a nonnegative function and γ : N → Z + be a nonnegative decreasing function such that:

(C) ||K n θ (x, .) -f θ (.)|| V T ≤ M (x)γ(n) where ||.|| V T denotes the total variation norm. Remark 4. X is geometrically ergodic if (C) holds with γ(n) = t n for some t < 1. X is uniform ergodic if (C) holds with M bounded and γ(n) = t n for some t < 1. X is polynomial ergodic of order m where m ≥ 0 if (C) holds with γ(n) = n -m .
The proof of the asymptotic normality is based on the following Lemma: Lemma 2. [see [Hayashi, 2000] and [Jones, 2004] for the proof.] Suppose that the conditions of the consistency hold. Suppose further that:

(1) Y i is α-mixing.

(2) (Moment condition): for some δ > 0 and for each j ∈ {1, • • • , p} :

E ∂m θ (Y 1 ) ∂θ j 2+δ < ∞ . (3) Assumption (C) holds such that E[M (X 1 )] < ∞ and γ(n) satisfies n γ(n) δ 2+δ < ∞ with δ defined in (2).
(4) (Hessian Local condition): for some neighbourhood U of θ 0 and for j, k ∈ {1, • • • , p}

E sup θ∈U ∂ 2 m θ (Y 1 ) ∂θ j ∂θ k < ∞.
Then, θ n defined in Eq.( 6) is asymptotically normal with asymptotic covariance matrix given by:

Σ(θ 0 ) = V -1 θ 0 Ω(θ 0 )V -1 θ 0
where V θ 0 is the Hessian of the application Pm θ given in Eq.(4).

Proof. It just remains to check that the conditions (2) and (4) of Lemma 2 hold under our assumptions.

Moment condition: as the function l θ is twice continuously differentiable w.r.t θ, for all y i ∈ R 2 , the application m θ (y i ) :

θ ∈ Θ → m θ (y i ) = ||l θ || 2 2 -2ϕ(y i+1 )u * l θ (y i
) is twice continuously differentiable for all θ ∈ Θ and its first derivatives are given by:

∇ θ m θ (y i ) = ∇ θ ||l θ || 2 2 -2ϕ(y i+1 )∇ θ u * l θ (y i ).
By assumption, for each j ∈ {1, • • • , p}, ∂l θ ∂θ j ∈ L 1 (R), therefore one can apply the Lebesgue differentiation Theorem and Fubini's Theorem to obtain :

∇ θ m θ (y i ) = ∇ θ ||l θ || 2 2 -2ϕ(y i+1 )u * ∇ θ l θ (y i ) . (10)
Then, for some δ > 0:

|∇ θ m θ (y i )| 2+δ = ∇ θ ||l θ || 2 2 -2ϕ(y i+1 )u * ∇ θ l θ (y i ) 2+δ ≤ C 1 ∇ θ ||l θ || 2 2 2+δ + C 2 ϕ(y i+1 )u * ∇ θ l θ (y i ) 2+δ , (11)
where C 1 and C 2 are two positive constants. By assumption, the function

||l θ || 2 2 is twice continuously differentiable w.r.t θ. Hence, θ → ∇ θ ||l θ || 2
2 is continuous on the compact subset Θ and the first term of equation ( 11) is finite. The second term is finite by the moment assumption (A9).

Hessian Local dominance: for j, k ∈ {1, • • • , p}, ∂ 2 l θ ∂θ j ∂θ k ∈ L 1 (R)
, the Lebesgue differentiation Theorem gives:

∇ 2 θ m θ (y i ) = ∇ 2 θ ||l θ || 2 2 -2ϕ(y i+1 )u * ∇ 2
θ l θ (y i ), and, for some neighbourhood U of θ 0 :

E sup θ∈U ∇ 2 θ m θ (Y i ) ≤ sup θ∈U ∇ 2 θ ||l θ || 2 2 + 2E sup θ∈U ϕ(Y i+1 )u * ∇ 2 θ l θ (Y i ) .
The first term of the above equation is finite by continuity and compactness argument and the second term is finite by the Hessian local dominance assumption (A10).

Proof of Corollary 1

By replacing ∇ θ m θ (Y 1 ) by its expression (10), we have:

Ω 0 (θ) = Var ∇ θ ||l θ || 2 2 -2ϕ(Y 2 )u * ∇ θ l θ (Y 1 ) = 4Var ϕ(Y 2 )u * ∇ θ l θ (Y 1 ) = 4 E ϕ(Y 2 ) 2 u * ∇ θ l θ (Y 1 ) u * ∇ θ l θ (Y 1 ) -E ϕ(Y 2 )u * ∇ θ l θ (Y 1 ) E ϕ(Y 2 )u * ∇ θ l θ (Y 1
) . Furthermore, by Eq.( 1) and by independence of the centered noise (ε 2 ) and (η 2 ), we have:

E ϕ(Y 2 )u * ∇ θ l θ (Y 1 ) = E (b 2 θ 0 + σ 2 θ 0 )(X 1 )u * ∇ θ l θ (Y 1 ) .
Using Fubini's Theorem and Eq.( 1) we obtain:

E (b 2 θ 0 + σ 2 θ 0 )(X 1 )u * ∇ θ l θ (Y 1 ) = E (b 2 θ 0 + σ 2 θ 0 )(X 1 ) e iY 1 z u ∇ θ l θ (z)dz = E (b 2 θ 0 + σ 2 θ 0 )(X 1 ) 1 2π 1 f * ε (z) e iY 1 z (∇ θ l θ ) * (-z)dz = 1 2π E (b 2 θ 0 + σ 2 θ 0 )(X 1 )e i(X 1 +ε 1 )z 1 f * ε (z) (∇ θ l θ ) * (-z)dz = 1 2π E e iε 1 z f * ε (z) E (b 2 θ 0 + σ 2 θ 0 )(X 1 )e iX 1 z (∇ θ l θ ) * (-z)dz, so that E (b 2 θ 0 + σ 2 θ 0 )(X 1 )u * ∇ θ l θ (Y 1 ) = 1 2π E (b 2 θ 0 + σ 2 θ 0 )(X 1 )e i(X 1 +ε 1 )z 1 f * ε (z) (∇ θ l θ ) * (-z)dz = 1 2π E e iε 1 z f * ε (z) E (b 2 θ 0 + σ 2 θ 0 )(X 1 )e iX 1 z (∇ θ l θ ) * (-z)dz = 1 2π E (b 2 θ 0 + σ 2 θ 0 )(X 1 ) e iX 1 z (∇ θ l θ ) * (-z)dz = 1 2π E (b 2 θ 0 + σ 2 θ 0 )(X 1 ) ((∇ θ l θ ) * (-X 1 )) * = E (b 2 θ 0 + σ 2 θ 0 )(X 1 )∇ θ l θ (X 1 ) . (12)
Hence, Ω 0 (θ) = 4 (P 2 -P 1 ) ,

where

P 1 = E (b 2 θ 0 + σ 2 θ 0 )(X 1 )∇ θ l θ (X 1 ) E (b 2 θ 0 + σ 2 θ 0 )(X 1 )∇ θ l θ (X 1 ) , P 2 = E ϕ(Y 2 ) 2 u * ∇ θ l θ (Y 1 ) u * ∇ θ l θ (Y 1 ) .
Calculus of the covariance matrix of Corollary (1): By replacing (∇ θ m θ (Y 1 )) by its expression (10) we have:

Ω j-1 (θ) = Cov ∇ θ ||l θ || 2 2 -2ϕ(Y 2 )u * ∇ θ l θ (Y 1 ), ∇ θ ||l θ || 2 2 -2ϕ(Y j+1 )u * ∇ θ l θ (Y j ) , = 4Cov ϕ(Y 2 )u * ∇ θ l θ (Y 1 ), ϕ(Y j+1 )u * ∇ θ l θ (Y j ) , = 4 E ϕ(Y 2 )u * ∇ θ l θ (Y 1 )ϕ(Y j+1 )u * ∇ θ l θ (Y j ) -E ϕ(Y 2 )u * ∇ θ l θ (Y 1 ) E ϕ(Y j+1 )u * ∇ θ l θ (Y j ) .
By using Eq.( 12) and the stationary property of the Y i , one can replace the second term of the above equation by: Finally, we obtain: Ω

(θ) = Ω 0 (θ)+2 ∞ j>1 Ω j-1 (θ) with Ω 0 (θ) = 4 (P 2 -P 1 ) and Ω j-1 (θ) = 4 Cj-1 -P 1 .
Expression of the Hessian matrix V θ : We have:

Pm θ = ||l θ || 2 2 -2 l θ , l θ 0 . (17)
For all θ in Θ, the application θ → Pm θ is twice differentiable w.r.t θ on the compact subset Θ. And for j ∈ {1, • • • , p}:

∂Pm ∂θ j (θ) = 2 ∂l θ ∂θ j , l θ -2 ∂l θ ∂θ j , l θ 0 = 2 ∂l θ ∂θ j , l θ -l θ 0 ,
= 0 at the point θ 0 , and for j, k ∈ {1, • • • , p}:

∂ 2 Pm ∂θ j ∂θ k (θ) = 2 ∂ 2 l θ ∂θ j θ k , l θ -l θ 0 + ∂l θ ∂θ k , ∂l θ ∂θ j j,k = 2 ∂l θ ∂θ k , ∂l θ ∂θ j j,k
at the point θ 0 .

5.3 M-estimator using the example in Section 4

Expression of f * ε : consider the random variable ε = ε-C √ V with ε = log(X 2 ) where X is standard Gaussian random variable, C = E[log(X 2 )] and V = V[log(X 2 )].
The Fourier transform of ε is given by:

E [exp (iεy)] = exp - iC √ V y E [exp (iεy)] = exp - iC √ V y E X 2iy √ V = exp - iC √ V y 1 √ 2π +∞ -∞ x 2iy √ V exp - x 2 2 dx
Using a change of variable z = x 2 2 , we get:

E [exp (iεy)] = exp - iC √ V y 2 iy √ V √ π +∞ 0 z iy √ V -1 2 e -z dz ≡ exp - iC √ V y 2 iy √ V √ π Γ 1 2 + iy √ V . Then f * ε = exp - iC √ V y 2 iy √ V √ π +∞ 0 z iy √ V -1 2 e -z dz ≡ exp - iC √ V y 2 iy √ V √ π Γ 1 2 + iy √ V .
The CIR process: taking that η i+1 ∼ N (0, 1) and ε i has a (log-) Chi-squared probability density function, if the Feller's condition holds true (a = 2κθ σ 2 ≥ 1) and c = 2κ σ 2 > 0, then the volatility process X i is stationary ergodic and ρ -mixing. The stationary distribution f θ is the gamma distribution γ(a, c) (see [START_REF] Genon-Catalot | Parameter estimation for discretely observed stochastic volatility models[END_REF]). On the other hand, the functions b θ , σ θ and l θ are given by: b θ (x) = (1 -κ)x + κθ and σ θ (x) = σ √ x l θ (x) = b 2 θ (x) + σ 2 θ (x) γ(a, c) =

(1 -κ) 2 x 2 + 2x (1 -κ) κθ + σ 2 + (κθ) 2 f θ (x), = α 1 x 2 + α 2 x + α 3 f θ (x).

where θ = (κ, µ, σ) and α 1 = (1 -κ) 2 , α 2 = 2 (1 -κ) κθ + σ 2 , α 3 = (κθ) 2 .

Therefore

l * θ (t) = E[e itX l θ (X)] = α 1 E[X 2 e itX ] + α 2 E[Xe itX ] + α 3 E[e itX ] with X ∼ Γ(a, c) = -α 1 ∂ 2 f * θ ∂t 2 (t) -iα 2 ∂f * θ ∂t (t) + α 3 f * θ (t).
After replacing f * θ (t) by 1 -it c -a , we obtain:

l * θ (x) = -α 1 -a c 2 (a + 1) 1 - ix c -a-2 + iα 2 a c 1 - ix c -a-1 + α 3 1 - ix c -a
.

It follows that for all θ, the squared norme of l θ is given by:

l θ 2 = b 2 θ + σ 2 θ (x) 2 Γ 2 (a, c)dx = β 1 x 4 + β 2 x 3 + β 3 x 2 + β 4 x + β 5 Γ 2 (a, c)dx,
where β 1 = α 2 1 , β 2 = 2α 1 α 2 , β 3 = 2α 1 α 3 + α 2 2 , β 4 = 2α 2 α 3 , β 3 = α 2 3 . Finally, using the non-centered moments of a Gamma-distributed random variable, E[X r ] = Γ(a+r) Γ(a)c r , we get: and the expression of the contrast function ( 8) is obtained. It is worth noting that the function u * l θ (y) defined in Definition 1 must be approximated numerically by using standard quadrature methods or Fast Fourier Transform.

β 1 x 4 Γ 2 (a, c)dx = β 1 x 4 c 2a Γ 2 (

Proof of Theorems 1 and 2 for the CIR process

Mixing property: under the Feller's condition, the volatility process X t is ρ-mixing and so αmixing by using the strong Markov property.

Regularity conditions: for the CIR process, the function l θ is given by the following polynomial function x → (α 1 x 2 + α 2 x + α 3 )f θ (x) with α 1 = (1 -κ) 2 , α 2 = 2 (1 -κ) κθ + σ 2 , α 3 = (κθ) 2 and for all x this function is smooth w.r.t. θ ∈ Θ. Hence, it remains to prove the moment condition and the local dominance to apply Theorem 2.

Since the function l θ is polynomial w.r.t θ belonging to the compact subset Θ, all the derivatives exist and in particular sup θ∈Θ l θ and sup θ∈Θ ∇ 2 θ l θ are finite. Furthermore, by combining the compactness argument and as the Fourier transform f * satisfies (see [START_REF] Fan | Nonparametric function estimation involving errors-invariables[END_REF]):

|f * ε (x)| = √ 2 exp - π 2 |x| 1 + O 1 |x| , |x| → ∞,
which means that f ε is ordinary-smooth in its terminology, we obtain:

           E sup θ∈Θ ϕ(Y 2 )u * l θ (Y 1 ) < ∞ E ϕ(Y 2 )u * ∇ θ l θ (Y 1 ) 2+δ < ∞ for some δ > 0, E sup θ∈U ϕ(Y 2 )u * ∇ 2 θ l θ (Y 1 ) < ∞
for some neighbourhood U of θ 0 .

'Figure 1 .

 1 Figure 1. Boxplot of the parameter κ. True value equal to 4.

Figure 2 .

 2 Figure 2. Boxplot of the parameter µ. True value equal to 0.03.

Figure 3 .

 3 Figure 3. Boxplot of the parameter σ 2 . True value equal to 0.4.

Figure 4 .

 4 Figure 4. Coverage with respect to the number of observations n = 500 up to 3000 for N = 100 estimators.

  = β 1 2 -(2a+3) c -3 Γ(2a + 3) Γ 2 (a)andβ 2 x 3 Γ 2 (a, c)dx = β 2 2 -(2a+2) c -2 Γ(2a + 2) Γ 2 (a) Γ(2a + 2, 2c)dx = β 2 2 -(2a+2) c -2 Γ(2a + 2) Γ 2 (a) , β 3 x 2 Γ 2 (a, c)dx = β 3 2 -(2a+1) c -1 Γ(2a + 1) Γ 2 (a) Γ(2a + 1, 2c)dx = β 3 2 -(2a+1) c -1 Γ(2a + 1) Γ 2 (a) , β 4 xΓ 2 (a, c)dx = β 4 2 -(2a) Γ(2a) Γ 2 (a) , β 5 Γ 2 (a, c)dx = β 5 2 -2a+1 c Γ(2a -1) Γ 2 (a).

Table 1 .

 1 Comparison of the CPU for Particle filters estimators, MCEM estimator and Contrast estimator.

	APF	KSAPF	APFS	EKF Contrast MCEM
	CPU (sec) 105.1695 93.8846 192.2166 0.2	20.4074 217430

Table 3 .

 3 Ratio signal-noise for the estimation of the CIR model

		Mean(μ n ) Mean(κ n ) Mean(σ 2 n ) MSE
	n = 500 and r = 0.1	0.0315	3.88	0.401	0.14
	n = 500 and r = 1	0.0303	3.89	0.405	0.16
	n = 1000 and r = 0.1	0.0312	3.76	0.401	0.11
	n = 1000 and r = 1	0.0308	3.83	0.41	0.18


Alain Gellé

Andreas M Läuchli

Brijesh Kumar

Frédéric Mila

A 2D quantum antiferromagnet with a four-fold degenerate valence-bond-solid ground state

Keywords: numbers: 75.10.Jm, 75.30.Kz, 75.50.Ee, 75.40.Mg

We study the competition between antiferromagnetic order and valence bond solid formation in a two-dimensional frustrated spin-1/2 model. The J1-J2 model on the square lattice is further frustrated by introducing products of three-spin projectors which stabilize four dimer-product states as degenerate ground states. These four states are reminiscent of the dimerized singlet ground state of the Shastry-Sutherland model. Using exact diagonalizations, we study the evolution of the ground state by varying the ratio of interactions. For a large range of parameters (J2 0.25J1), our model shows a direct transition between the valence-bond-solid phase and the collinear antiferromagnetic phase. For small values of J2, several intermediate phases appear which are also analyzed.

I. INTRODUCTION

In the last several years, a large amount of work has been devoted to the study of quantum systems with frustrated magnetic interactions due to their propensity to present spin liquid states. These phases indeed present interesting low energy properties, starting with the absence of magnetic order even at zero temperature. Recently, a lot of attention has been focused particularly on the Valence Bond Solids (VBS), and more precisely on the transition between these states and a magnetically ordered phase. In this context, it has been suggested that this transition could be a second order transition that does not belong to the Landau-Ginsburg paradigm 1,2 .

In a VBS state, spins are coupled in pairs forming singlet states, evocating valence bonds. These pairwise singlets are themselves arranged in a periodic pattern. Such states break the translation symmetry, and it is possible to define an order parameter quantifying the singlet long range order. The question on the nature of quantum phase transition that separates this state from a longrange ordered magnetic phase is very interesting as well as current. Supposing a second order transition, according to Landau-Wilson paradigm [START_REF] Landau | Statistical Physics[END_REF] , the order parameter of both phases should vanish precisely at the transition. It seems more likely that the two parameters will not vanish exactly at the same point, leading either to a first order transition, or to two second order transitions separated by an intermediate phase. Considering spin-1/2 on square lattice, Senthil and coworkers have recently suggested that there could exist second order transitions that are not described by Landau-Ginsburg theory 1,2 . The transition could instead be described by means of fractional degrees of freedom, namely spinons. These spinons become deconfined at the transition point, called the deconfined quantum critical point.

A priori, the spin-1/2 J 1 -J 2 antiferromagnetic model on the square lattice seems the simplest choice to inves-tigate the relevance of this theory. In the two limiting cases where either J 1 or J 2 is very large compared to the other, the system presents antiferromagnetic order. In the former case it corresponds to the usual Néel order, and in the latter case the ground state has a collinear antiferromagnetic order, corresponding to the Néel order on two sublattices (that are obtained by connecting second neighbor sites). It is generally accepted that this model presents an intermediate spin disordered phase in the range 0.4 < J 2 /J 1 < 0.6, which might break the translational symmetry [START_REF] Schulz | [END_REF]5,6,7,8,9,10 . It is still unclear if the transition between the Néel state and the spin disordered phase is a deconfined quantum critical point or a simple first order transition, although recent works are in favor of a weak first order transition 11,12 .

The difficulty in understanding the J 1 -J 2 model comes from the intermediate phase whose nature is subject to discussions 13 . Some studies show a four-spin plaquette order 9,10 , while others are in favor of a columnar dimer order 6,7,8 . In the latter case, not only the translational symmetry, but also the rotational symmetry is broken. However, if this phase develops dimer-dimer correlations, it is far away from being the simple direct product of dimer singlet wave functions. Indeed, if that was the case, the ground state would present a strong signal of long range order which is not observed for the intermediate phase of the J 1 -J 2 model. Exact diagonalization studies show that the dimer-dimer correlation is rapidly decreasing with distance 6,14 , and that, if the dimer-dimer order exists, it should be rather small 14 . The analysis of the properties of the phase transition is therefore not so easy in the J 1 -J 2 model. Recently, a new model including ring exchange has been proposed 15 to explore the possibility of non Landau-Ginsburg phase transitions. It has however been shown that this model presents a first order transition 15 .

In the present paper, we propose a new quantum spin-1/2 model with frustrated antiferromagnetic interactions. Interestingly, for a simple choice of the interac-tion parameters, this model has an exactly solvable singlet ground state. The exact ground state is a pure direct product of dimers, arranged in the same pattern as that in the ground state of the Shastry-Sutherland model 16 . However, unlike in the Shastry-Sutherland model, the Valence Bond Solid (VBS) ground state of our model presents a case of spontaneous symmetry breaking, and is four-fold degenerate 17,18 . Using exact diagonalization, we investigate the evolution of the ground state as the interaction parameters are varied away from the exactly solvable VBS case, particularly towards magnetically ordered phases. In the following section, we describe the model and discuss its main properties and its relation to the J 1 -J 2 model. In section III-A, we study the competition between the VBS phase and the antiferromagnetic collinear phase. In section III-B, we analyze the competition with the usual Néel antiferromagnetic phase. Finally, the last section is devoted to conclusions and perspectives. Consider a system of spin-1/2 on the square lattice. The model we consider contains only two and four spins interactions which can be conveniently described using a six-site rectangular plaquette (see Fig. 1). It is also convenient to distinguish between the two sublattices A and B of the square lattice. The six sites of a plaquette thus contain 3 sites belonging to the A sublattice and other 3 sites belonging to the B sublattice. For each set of three spins, we consider the spin projector on the quartet (spin 3/2) state:

P A i,j,k = | S i + S j + S k | 2 - 3 4 (1)
P B l,m,n = | S l + S m + S n | 2 - 3 4 (2)
where A and B refer to the two sublattices, and i, j, k, l, m, n are the sites of the plaquettes as depicted in Fig 1 . The interaction we consider is obtained by tak-ing the product of the two projectors of a plaquette:

H 0 = [i,j,k,l,m,n] 1 4 P A i,j,k P B l,m,n (3)
where the sum runs over all horizontal and vertical plaquettes, and [i, j, k, l, m, n] are the sites of one plaquette. One should notice that the value of a projector is always positive except when the three spins are in a doublet state (that is when it becomes zero). The latter condition is fulfilled, in particular, if any two of the three spins form a singlet. These two spins could either be two first neighbor or two second neighbor spins of a sub-lattice, and correspond respectively to two second or two third neighbor spins of the original square lattice. Since the interaction in the model is the product of two projectors, it is only necessary to cancel one of these projectors to minimize the corresponding term. For this purpose, it is thus sufficient to have one dimer on the plaquette, either between second or third neighbor spins. Therefore, it is obvious that if a valence bond configuration has one such dimer on every plaquette, it will form a zero energy ground state of H 0 . The idea of stabilizing valence bond states by spin projectors is also contained in the Majumdar-Ghosh model on the linear chain 19 . Recently, Batista and Trugman studied a model with such projectors on four-site square plaquettes 15 , but their Hamiltonian has a highly degenerate ground state since a great number of dimer patterns minimize all projectors. In the present model, the ground state is obtained when dimers are arranged in the pattern given in Fig. 2 and is much less degenerate (see also the discussion in the Appendix). This arrangement evokes the Shastry-Sutherland model (and we will further refer to this phase as the SS-VBS phase), but contrary to the SS model, our Hamiltonian does not break the translational symmetry of the square lattice. The ground state presents a spontaneous symmetry breaking, and is fourfold degenerate. The other three of SS-VBS states can be obtained by translation of the dimer pattern by lattice vectors. Interestingly, the four SS-VBS states are also exact zero energy eigenstates (although not the ground state) of the nearest-neighbor Heisenberg model on the square lattice, which motivated us to seek for a model on the square lattice with four SS-VBS states as the exact ground state 20 .

The spin projectors are two-spin operators which can be rewritten in terms of the exchange couplings between the first and second neighbor spins of a sublattice:

P A i,j,k = 2( S i • S j + S i • S k + S j • S k ) + 3 2 (4)
Therefore, the product of two projectors generates twoand four-spin operators in H 0 . As the projector is invariant under rotation and involves spins on a given sublattice, it will commute with the total spin of this sublattice, and obviously, with the total spin of the other sublattice. H 0 will thus conserve the total spin on each sublattice. This extra symmetry may introduce a further factor two in the degeneracy of excited states.

It is interesting to notice that this Hamiltonian is somehow related to the J 1 -J 2 model. Let us assume that a state can be expressed as a direct product of the wavefunction of spins on A sublattice |Ψ A and of spins on B sublattice |Ψ B :

|Ψ = |Ψ A ⊗ |Ψ B (5)
In that case, it is possible to define a Hamiltonian for one sublattice. Let us, for instance, choose the A sublattice:

H A = [i,j,k,l,m,n] 1 4 P A i,j,k Ψ B |P B l,m,n |Ψ B (6)
This Hamiltonian corresponds, for the A sublattice, to a model with first and second neighbor interactions whose amplitudes are modulated by the local spin state of the B sub-lattice. It is easy to show that, in the case where the mean value of the projector on B sublattice is homogeneous, the Hamiltonian can be expressed as:

H A = P B <i,j>A S i • S j + 1 2 P B ≪i,j≫A S i • S j (7)
where < i, j > A are the first neighbor couples of spins of A sublattice, and ≪ i, j ≫ A the second neighbor ones, P B is the uniform value of Ψ B |P B l,m,n |Ψ B . This Hamiltonian is precisely the J 1 -J 2 model on A sublattice at the J 2 /J 1 = 0.5 point, which is located in the intermediate spin liquid phase of the model. Clearly, if we consider the ground state of our model, we do not expect the mean value of the P B projector to be homogeneous. Nevertheless, this comparison introduces a simple picture of the model, in which the local spin correlations self-consistently modulate the J 1 -J 2 model on each sublattice. This process allows the system to stabilize a pure direct product of dimers as ground state. This valence bond ground state is different in nature from the one expected in the intermediate phase of J 1 -J 2 model, since the dimer order on each sub-lattice is not columnar but staggered 21 (see Fig. 2). The existence of this exact ground state makes this model an interesting candidate to explore the transition between the VBS phase and the antiferromagnetically ordered state. In order to drive the system into such phases, we consider the following Hamiltonian:

H = J(1 -γ)(1 -δ) <i,j> S i • S j +J γ(1 -δ) ≪i,j≫ S i • S j + J δ H 0 (8)
where < i, j > and ≪ i, j ≫ are respectively first and second neighbor spins of the original (square) lattice, γ and δ are dimensionless parameters, and J is the energy scale (assumed to be positive). When δ is equal to 1, one just retrieves the H 0 model, while when δ = 0 is zero, the Hamiltonian is simply the J 1 -J 2 model, with γ being equal to J 2 /(J 1 + J 2 ). A sketch of the phase diagram of the Hamiltonian ( 8) is shown in Fig. 3. It contains three phases where the order has been unambiguously identified denoted by N (antiferromagnetic (π, π) Néel order), col. AF (collinear (π, 0)/(0, π) order with Néel order on the two sublattices) and SS-VBS (Valence-Bond order with Shastry-Sutherland arrangement). Between these phases, there is a region where correlations change very significantly, defining possible phase transitions between three phases denoted by A, B and C. These phases are discussed in the last part of the next section.

III. NUMERICAL RESULTS

We study the evolution of eigenstates of this model by an exact diagonalization technique based on the Lanczos algorithm. Finite clusters of N = 16, 20 and 32 sites were used. Larger clusters were not accessible because of the complexity of the interaction. The first noticeable result concerns the δ = 1 case, when H = JH 0 . We checked numerically that the four SS ground states presented in previous section are the only ground states of the model. This is always true except for the 16 site cluster, for which two more ground states are also present, as shown in Fig. 4. These two states are ground states only for the 16 site cluster, due to the very short loops which wrap around the boundaries of the cluster (see also the Appendix). Indeed, in this case a third neighbor dimer between the sites (x, y) and (x, y+2) also represents a dimer between the sites (x, y + 2) and (x, y + 4) because of the periodic boundary conditions. This is no longer true for larger clusters and by extension in the thermodynamic limit.

In the following two subsections we map out the phase diagram of the model, with emphasis on two lines: In subsection III A we study first the case γ = 1, which connects to the case of J 1 = 0, J 2 = J once δ = 0. In subsection III B the case γ = 0 is presented, which connects the unfrustrated square lattice Heisenberg antiferromagnet to the fourfold dimerized ground state.

A. SS-VBS versus Collinear Order (γ = 1)

In the present section, we consider the case where γ = 1. According to the previous discussion on the relation between our model and the J 1 -J 2 model, J 2 corresponds to a first neighbor interaction on each sublattice which, for large enough values, will lead to a Néel state on each sublattice.

For δ = 1, the SS-VBS states are exact ground states even for finite clusters. Note that these four states are nonorthogonal on a finite cluster: they have a finite overlap that decreases exponentially with the cluster size 22 . However, since they are linearly independent, the ground state is indeed four-fold degenerate on a finite cluster for δ = 1. This degeneracy is lifted however by the J 2 interactions that appear for δ < 1. Fig. 5 represents the energy differences between the ground state and the lowest state of some symmetry sectors, obtained for the 32 site cluster. The energy differences between the four lowest singlet states stay relatively small in the range 0.8 < δ ≤ 1, while the spin gap progressively decreases. Below 0.8 the energy difference between singlet states rapidly increases, while the spin gap remains approximatively constant suggesting that the system is in the antiferromagnetic phase that will present a zero spin gap in the thermodynamic limit. In order to further investigate the transition, we calculate the Q dependent magnetic susceptibility for each cluster 25 :

M 2 N ( Q) = 1 N (N + 2) i,j S i S j e i Q( rj -ri) (9)
where r i denotes the position of i th spin, <> the mean value in the ground state and N the number of sites in the cluster. The evolution of the Q 1 = (π, 0) magnetic susceptibility relevant for antiferromagnetic collinear order is presented in Fig. 6. The extrapolation to the thermodynamic limit of the corresponding sublattice magnetization 4 has been performed using the finite size scaling predicted by non-linear sigma model studies 23,24 :

M 2 N ( Q 1 ) = 1 8 m 2 0 ( Q 1 ) + const. √ N (10)
The fitted value of m 0 ( Q 1 ) is also shown in Fig. 6. The extrapolated magnetization stays large up to δ ≃ 0.7, which confirms that the collinear phase is stable in this range of parameter. It then rapidly drops, and vanishes around δ ≈ 0.75. The VBS phase, which is expected for larger values of δ, is characterized by long range dimer-dimer correlation. This long range correlation can be considered as the order parameter of this symmetry-breaking phase. In order to determine the stability of the phase, we compute the following dimer-dimer correlations in the ground state: (11) where 0 stands for the origin, and r 1 and r 2 can be either equal to (1, 1), (1, -1), (-1, 1) or to (-1, -1). As expected, for δ ∼ 1, values of D( r) that corresponds to the dimers of the Shastry-Sutherland pattern are quite large. For a given cluster, and close to δ = 1, fluctuations of these values are very small, of the order of a few percent of the average value of these correlations. Since we are interested in the value of D( r) for r going to infinity, we only considered among these correlations the one obtained for the largest r value (D N ( r m )). These correlations, shown in Fig. 7, are quite small in the antiferromagnetic phase, and rapidly increase at around δ ∼ 0.8. Interestingly enough, the curves cross at δ = 0.76, very close to the point were the antiferromagnetic order vanishes. This behaviour is consistent with a first order transition, with an order parameter scaling down to zero with the cluster size below a critical value, and scaling up to a finite value above. However, with only two sizes available (the 16 site cluster turns out to be rather pathological with essentially δ independent correlations), this information should be taken with care, and a definitive identification of the nature of the phase transition requires further investigation.

D N ( r) = ( S 0 • S r1 )( S r • S r+ r2 ) -S 0 S r1 S 0 S r2

B. Competition with Néel order (γ = 0)

We now consider the case of γ = 0, which corresponds to the competition between the first neighbor coupling and the six sites plaquette interaction. The interest in the model obtained for γ = 0 comes from the fact that the SS-VBS states remain eigenstates for all values of δ with an energy equal to zero. It follows that the transition at which the valence bond state vanishes is necessarily a level crossing and is therefore first order.

It should be emphasized here that, in the previous section, both SS-VBS and the collinear order phases presented antiferromagnetic correlations between second neighbor spins. Therefore, a direct transition between them seems to intuitively exist as if the spin-liquid phase (SS-VBS) knows, through the spin-spin correlations present in it, to which ordered phase it must go. This is no more the case for the Néel order which is stabilized by the first neighbor coupling, and presents antiferromagnetic correlations between first neighbor spins, and ferromagnetic correlations between second neighbor spins. According to this simple consideration, a direct transition between the SS-VBS phase and the Néel phase seems therefore unlikely. We start the discussion by presenting in Fig. 8 the evolution of the energies of some of the lowest eigenstates obtained for the 32-site cluster, taking the lowest fully symmetric k = (0, 0) level as the energy reference. Based on this figure and Fig. 10, we are lead to identify three different phases, tentatively labeled "A", "B" and "C" in addition to the well characterized (π, π) Néel phase at δ = 0 and the fourfold degenerate Shastry-Sutherland states for δ 0.4. We discuss each of these three phases in the following. Phase "A" -Starting at small values of δ, we note that the plaquette interactions favor antiferromagnetic correlations between second neighbor spins. Therefore one could expect the evolution of the system to be similar to the J 1 -J 2 model for small values of J 2 . The symmetry of the four lowest states near δ = 0.15 with momenta (0, 0) (two states), (π, 0) and (0, π), are indeed compatible with the hypothesis of a four-fold degenerate ground state with a translational symmetry breaking. Fig. 9 shows the correlations defined in Eqn. (11), this time between dimers made of first neighbor spins obtained for δ = 0.15. These correlations are relatively large, and show a columnar ordering of the dimers. Nevertheless, one should note that, as for the J 1 -J 2 model, it is difficult to determine whether this phase presents columnar dimer order or plaquette order in the thermodynamic limit. So we believe the phase "A" is formed of some sort of valence bond solid with dimers on nearest-neighbor sites.

Phase "C" -At the other end of the δ axis, starting from large values of δ, one can see that the ground state stays exactly fourfold degenerate down to δ ≃ 0.4, below which a level crossing occurs. Near δ = 0.35 many levels are very close in energy. Some of these states are even found to be lower in energy than the fully symmetric k = (0, 0) state -which is again the ground state for δ ≤ 0.3, but this may well be a finite size effect on this particular sample.

In the case γ = 1 studied in section III A, we en-countered a direct transition between a collinear (π, 0) Néel ordered phase and the Shastry-Sutherland type VBS state. In order to test this scenario here, we determined the static spin structure factors for different momenta in Fig. 10. Indeed the (π, 0) components are strongest around δ ∼ 0.35. In order to shed further light on the presence of magnetic long range order we study the evolution of the collinear magnetic order as a function of γ for a fixed value of δ = 0.35. Using samples of 20 and 32 sites we obtain in Fig. 11 a finite size scaling which shows that the magnetic order is lost at a finite value of γ ∼ 0.1, i.e. the point δ = 0.35, γ = 0 does not sustain magnetic long range order. Looking at real-space spin correlations, the second neighbor correlation is much larger than all other spin correlations, so that the possibility of a dimer VBS state needs to be checked. We therefore also computed the real-space dimer-dimer correlations which are presented in Fig. 12-(a). They present an ordering reminiscent of the SS-VBS phase, although correlations are smaller than in the pure phase. For the purpose of comparison the same correlations in the pure SS-VBS phase are displayed in Fig. 12-(b). We recall at this stage that in this range of parameters, several states are very close in energy (see Fig. 8), all singlets. We also observe that they all present the same kind of dimer-dimer correlations. From our results it is difficult to determine if these dimer-dimer correlations are short ranged or long ranged. However it seems as though this phase could be best visualized by a condensation of singlet excitations above the SS ground states, and not by a simple level crossing into a new ground state of completely different character, as it seems to happen in the original Shastry-Sutherland model 26 . Phase "B" -Upon close inspection of Fig. 10 it is rather natural to suppose that there exists a third phase around 0.2 δ 0.3 sandwiched between the phases "A" and "C" which displays enhanced magnetic correlations for Q = (π, π/2) and Q = (π/2, π/2). Unfortunately the presence of correlations at these wavevectors renders the study of the system more difficult, since these momenta are not present on the 20 site cluster. It is therefore not possible to perform a finite size scaling of the Q dependent magnetic susceptibility. Nevertheless, one should note that there is also an enhancement of these components on the 16 site cluster. At that stage, it is difficult to characterize this phase "B" and to know if these large components correspond to a long range spin order, or to some more exotic phase. N ( Q1) square magnetic susceptibility for 20 and 32 sites clusters and its extrapolation to thermodynamic limit M∞( Q1) 2 = m 2 0 ( Q1)/8, calculated for the lowest state with momenta k = (0, 0) and highest symmetry.

IV. CONCLUSION

We introduced a new model with frustrated interactions which provides an interesting case of competition between antiferromagnetic orders and a valence bond solid order. We have shown that, for some values of the interaction parameters, the four fold degenerate VBS ground state is exactly a direct product of the dimer singlet wave functions. This model is thus an interesting candidate for investigating the possibility of a newly proposed scenario of quantum phase transition. It indeed presents, for a large range of parameters (namely γ > 0.2), a direct transition between the SS-VBS phase and a collinear antiferromagnetic phase. Further investigations are needed to determine if the transition is first order, of if it could correspond to the non-Landau-Ginsburg transition proposed by Senthil and coworkers 1,2 .

For smaller values of the second neighbor interaction (i.e. small values of γ), the Néel and SS-VBS phases are separated by an intermediate region where different types of correlations dominate depending on the value of δ. It is likely that these correlations are the trace of exotic intermediate phases, but further work is clearly needed to fully characterize these phases and the nature of the transitions between them. Let us now consider one diagonal dimer. The remaining sites of the square plaquette on which this dimer sits have to be part of a dimer. But since two dimers cannot be on the same rectangular plaquette, the only possibility is that these dimers are perpendicular to the first one. This is precisely the prescription to construct a Shastry-Sutherland state. The freedom to chose the position and orientation of the first dimer leads to four different states.

The exact diagonalization calculations on the 20 site and 32 site clusters presented in this paper support this assertion.

Note that with periodic boundary conditions of length 4, a dimer constructed from third neighbours satisfies four rectangles, which leads to two additional ground states on the 16-site cluster.

FIG. 1 :

 1 FIG. 1: Representation of a 3 by 2 horizontal rectangular plaquette. A vertical 2 by 3 plaquette can be obtained by a 90 • rotation of the plaquette. Gray plain and dashed lines represent the nearest-neighor bonds of the A and B sub-lattices.

FIG. 2 :

 2 FIG. 2: Representation of one of the four ground states of model H0. Black lines indicate the pairs of spins that form singlet states. The three other ground states can be obtained by translation of the dimer pattern by lattice vectors.

FIG. 3 :

 3 FIG.3:(color online) Sketch of the phase diagram of the model(8). Phase transitions obtained by exact diagonalizations are indicated by full dots. The lines are guide to the eyes. The nature of the different phases is described in the following sections.

FIG. 4 :

 4 FIG.4: Representation of one of the two additional ground states for the 16 site cluster. Full lines connect the 16 sites of the cluster, dashed lines symbolize the boundary conditions. The other additional ground state can be obtained by translation of the dimer pattern by lattice vectors.

FIG. 5 :

 5 FIG.5:(color online) Energy differences between the ground state and some of the lower states of the 32 site cluster. The symmetry of the different states is indicated on the figure. States with (π, 0) -A1 and (0, π) -A1 symmetry are triplet states, others are singlet states.

FIG. 6 :

 6 FIG.6:(color online) MN ( Q1) 2 square magnetic susceptibility for 20 and 32 sites clusters and its extrapolation to thermodynamic limit M∞( Q1) 2 = m 2 0 ( Q1)/8.

FIG. 7 :

 7 FIG.7:(color online) Dimer-dimer correlation Dm obtained for the 20 and 32 sites cluster. This value correspond for one cluster to the correlation between the two dimers of the Shastry-Sutherland pattern separated by the largest distance.

FIG. 8 :

 8 FIG.8:(color online) Energy differences between the lower fully symmetric state with k = (0, 0) momenta and some of the lowest states of the 32 site cluster. The symmetry of the different states is indicated on the figure.

FIG. 9 :

 9 FIG.9:(color online) Correlations between first neighbor spin dimers on the 32 sites cluster for γ = 0 and δ = 0.15. Positive values are represented by (blue) plain lines, and negative values by (red) dashed lines. The thickness is proportional to the amplitude of the correlation as depicted on the figure. The (black) dimer in the lowest corner of the cluster is the reference dimer.

FIG. 10 :

 10 FIG.10:(color online) Q dependent magnetic susceptibility obtained for the 32 sites cluster.

FIG. 11 :

 11 FIG.11:(color online) M 2 N ( Q1) square magnetic susceptibility for 20 and 32 sites clusters and its extrapolation to thermodynamic limit M∞( Q1) 2 = m 2 0 ( Q1)/8, calculated for the lowest state with momenta k = (0, 0) and highest symmetry.

FIG. 12 :

 12 FIG.12:(color online) Representation of correlations between second neighbor spin dimer on the 32 sites cluster calculated for the lowest state with momenta k = (0, 0) and highest symmetry, (a) for γ = 0 and δ = 0.375, (b) for γ = 0 and δ > 0.4. Conventions are the same than on Fig.9

Acknowledgments

We acknowledge the support of the Swiss National Fund and of MaNEP. The exact diagonalization computations have been enabled by allocation of resources on the IBM Regatta machines of CSCS Manno (Switzerland).

APPENDIX

The Hamiltonian H 0 (Eq. 3 of section II) has four zero energy SS-VBS singlet configurations forming an exact ground state. Proving (in a mathematically rigorous way) that these four SS-VBS states are the only states in the ground state and there exists no fifth state, is a non-trivial and hard task. We are not going to attempt it here. Historically, for the Majumdar-Ghosh model, which is considerably simpler as compared to H 0 , it was already very hard to prove the exact two-fold degeneracy of the ground state (which was eventually shown by AKLT). However, it was much easier to show that there are two dimerized singlet configurations which form the exact ground state of the Majumdar-Ghosh model, and to argue that other dimer-singlet configurations, generated by the allowed variations, will not be the eigenstates. We will do a similar exercise for H 0 , showing that the four SS-VBS state are the only allowed dimer configurations in the ground state.

The block Hamiltonian, h 6 = P A P B , of a six-site plaquette is the basic building block of H 0 . Since P A and P B are the spin projectors, the lowest energy of h 6 is zero. This corresponds to either P A or P B or both becoming zero in a given spin configuration of the block. This happens when the three A sublattice spins in a sixsite plaquette form a total spin=1/2 state, or the same thing happens for B sublattice spins or for both. One way, in which this can be achieved, is by forming exactly one singlet bond out of three A or B sublattice spins of a plaquette. Thus a simple rule emerges for constructing the dimerized ground state of H 0 . If a dimer configuration on the full square lattice is such that, on every six-site plaquette, there exists at least one singlet bond (dimer) between only A or only B sublattice spins, then all the plaquette Hamiltonians can be simultaneously satisfied (that is, every h 6 is in its ground state), and such a configuration will be an exact zero energy ground state of H 0 . Now, the number of rectangular plaquettes is equal to two times the number of sites, hence four times the number of dimers of any dimer covering of the lattice. Since a dimer belongs at most to four rectangles, to satisfy all rectangles simultaneously, each dimer should belong exactly to four rectangles, and each rectangle should contain a single dimer. Since a dimer constructed from third neighbours belongs to only two rectangles, such dimers have to be rejected, and one should only use diagonal dimers.


Marcelo Mandelli
email: mandelli@ieee.org

Luciano Ost
email: luciano.ost@leicester.ac.uk

Gilles Sassatelli
email: sassatelli@lirmm.fr

Fernando Gehm Moraes

Fernando Moraes
email: moraes@pucrs.br

Trading-off System Load and Communication in Mapping Heuristics for Improving NoC-Based MPSoCs Reliability

Keywords: MPSoC, NoC, dynamic mapping, lifetime, reliability

come   L'archive ouverte pluridisciplinaire

Introduction and Related Work

Semiconductor technology evolution has enabled to reduce transistors' size and, consequently, to integrate multiple cores in the same chip. MPSoCs are expected to integrate thousands of cores by the end of the decade [START_REF]International Technology Roadmap for Semiconductors[END_REF], providing a high computational performance. However, such systems became more susceptible to failures due to transistor shrinking, which can result in high power density and temperature [START_REF] Meyer | Cost-effective and yield optimization for NoC-based MPSoCs[END_REF][3] [START_REF] Wang | System-level reliability exploration framework for heterogeneous MPSoC[END_REF]. In NoC-based MPSoCs, it is also important to consider the NoC links reliability. A higher data volume transferred through the NoC may induce links to fail, which may make processors unusable. Unusable processors induce mapping of applications onto other system processors, increasing their workload and, consequently, reducing their lifetime. Furthermore, links failures compromise system performance by increasing NoC congestion, latency and routing algorithm computation.

Such systems are likely to execute several applications in parallel. The workload imposed by such applications and its distribution impact considerably on system reliability [START_REF] Chandra | Quantifying workload dependent reliability in embedded processors[END_REF] [START_REF] Huang | Lifetime reliability-aware task allocation and scheduling for MPSoC platforms[END_REF]. Workload imbalance can generate hotspots zones (i.e. peaks of power dissipation) and consequent thermal implications, leading to unreliable system operation. Thus, task mapping has to consider the wear state of processors to better balance applications workload over the system [START_REF] Huang | Lifetime reliability-aware task allocation and scheduling for MPSoC platforms[END_REF].

The literature presents different task mapping heuristics to improve system lifetime reliability. Huang et al. [START_REF] Huang | Lifetime reliability-aware task allocation and scheduling for MPSoC platforms[END_REF] propose a task mapping allied to scheduling strategies based on simulated annealing, which also considers applications deadlines.

Hartman et al. [START_REF] Hartman | Lifetime improvement through runtime wear-based task mapping[END_REF] present a run-time task mapping that uses sensors to capture wear patterns in the system, aiming balancing the processors' loads. Das et al. [START_REF] Das | Reliability-driven task mapping for lifetime extension of networks-on-chip based multiprocessor systems[END_REF] propose a wearbased task mapping technique that generates mapping solutions at design-time, aiming to satisfy application deadlines and to maximize lifetime of cores and NoC links. This solution is improved in [START_REF] Das | Temperature aware energy-reliability trade-offs for mapping of throughput-constrained applications on multimedia MPSoCs[END_REF] by including an offline DVFS technique defining voltage and frequency levels of all cores. Bolchini et al. [START_REF] Bolchini | Run-time mapping for reliable manycores based on energy/performance trade-offs[END_REF] present a run-time task mapping technique aiming to improve system reliability under energy and performance constraints. For this purpose, this work combines the technique of Das et al. [START_REF] Das | Reliability-driven task mapping for lifetime extension of networks-on-chip based multiprocessor systems[END_REF] with a remap strategy. Chantem et al. [START_REF] Chantem | Enhancing multicore reliability through wear compensation in online assignment and scheduling[END_REF] propose both dynamic task mapping and scheduling. Authors adopt the LTF-based algorithm [START_REF] Chen | Approximation Algorithms for Multiprocessor Energy-Efficient Scheduling of Periodic Real-Time Tasks with Uncertain Task Execution Time[END_REF] for dynamic task mapping. A scheduling technique is executed periodically, depending on system wear state conditions, trying to compensate uneven core wear states.

Most reviewed works use static mapping techniques [START_REF] Huang | Lifetime reliability-aware task allocation and scheduling for MPSoC platforms[END_REF][8][9] [START_REF] Bolchini | Run-time mapping for reliable manycores based on energy/performance trade-offs[END_REF], which are not suitable to deal with the dynamic aspects of the system. Among fully dynamic mapping heuristics [3][7], only [START_REF] Hartman | Lifetime improvement through runtime wear-based task mapping[END_REF] considers NoC links in the mapping decision. Embedded silicon sensors are used in [START_REF] Hartman | Lifetime improvement through runtime wear-based task mapping[END_REF], while our approach is software-based, which makes our system more flexible while consuming less energy and area when compared to Hartman's approach.

The main contribution of this work is a fully dynamic mapping heuristic that makes a trade-off between processors' workload and NoC communication volume, improving systems reliability.

MPSoC Architecture

This work adopts a homogeneous MPSoC architecture divided into clusters [START_REF] Mandelli | Enhancing Performance of MPSoCs through Distributed Resource Management[END_REF] as presented in Figure 1. The processing elements (PEs) are interconnected through a NoC. A PE contains a processor, a local memory (scratchpad memory, without caches), a DMA module and a router. The NoC assumes a 2D-mesh topology, credit-based flow control, round-robin arbitration, and XY routing algorithm. An external memory, called task repository, stores all applications, which are loaded into the system at runtime.

PEs may act as: Global Manager Processor (GMP), Local Manager Processor (LMP) and Slave Processor (SP). SPs are responsible for executing user's applications, running a simple operating system (OS) that allows task communication and multitask execution. An SP local memory is divided into k pages, also called resources. A task executing in the system is mapped in one resource. Therefore, an SP may execute k simultaneous tasks. If a resource does not have a task mapped on it, it is considered free or available. LMPs are responsible for cluster control, such as executing task mapping and reclustering. GMP is a single processor in the system that accumulates the functions of an LMP and the overall system management functions, including application-to-cluster mapping, external devices access (e.g. task repository). All clusters have a fixed size, defined at design-time. At runtime, a reclustering process may modify the cluster size, according to the application mapping. When an application is mapped in a given cluster, its LMP checks the availability of SPs, and if necessary, SPs are borrowed from neighbor clusters. When the application finishes its execution, the borrowed SPs are released to the original cluster.

Mapping Heuristics

This work models an application as a directed graph GApp = (T, E), where each vertex t i ∈ T represents an application task and each directed weighted edge e ij ∈ E represents a communication dependence between tasks t i and t j . The weight of an edge e ij is denoted by comm ij , representing the total data communication volume transferred between application tasks t i and t j . The set T is divided in two subsets iT and niT, where (iT U niT) = T. The subset iT contains the application initial tasks, which are those without dependences to other tasks; and the subset niT contains non-initial tasks. A task t i ∈ T contains: • a set called communicating task list, defined as

C i = {(t j , comm ij ); (t k , comm ik ); … (t n , comm in )}
, where each element is a tuple containing a task t j that communicates with t i and the volume comm ij transferred between t i and this task. • a load value Ld i , which represents the energy consumption [START_REF] Rosa | Fast Energy Evaluation of Embedded Applications for Manycore Systems[END_REF] related to the execution of this task on an SP. Each application is stored in the task repository with a header called application description. The application description defines the application size, the application initial tasks and information about application tasks, such as tasks' communicating task list and load.

This work adopts a mapping process divided into three steps, which may vary according to the employed mapping heuristic:

• cluster selection: whenever a new application is required to be mapped, the system alerts the GMP, which verifies if the system has available resources to map the entire application. In not, the application is scheduled to be mapped later. Otherwise, the GMP selects a cluster to map the required application, according to the employed heuristic. Then, the required application description is sent to the LMP of the selected cluster. • initial task mapping: the LMP of the selected cluster receives the application description. Then, the LMP maps the application initial tasks inside the cluster, according to the used heuristic. The mapping of initial tasks will start the application execution. • non-initial task mapping: the mapping of non-initial tasks occurs whenever a given task t i needs to communicate with a non-mapped task t j . In this situation, the mapping of t j is requested by t i to the LMP of the cluster. LMP maps the task t j on the SP selected by the adopted heuristic. Three mapping heuristics are evaluated in this work. Section 3.1 presents the reference mapping heuristic proposed in [START_REF] Mandelli | Multi-Task Dynamic Mapping onto NoC-based MPSoCs[END_REF], which main cost function is to reduce communication volume transferred through the NoC. Section 3.2 proposes a heuristic aiming to balance processors' load in the system, and 3.3 proposes a heuristic that aims to make a trade-off between communication volume reduction and processor load balance. All heuristics are detailed according to the three-step mapping process previously discussed.

Premap-DN

The Premap-DN heuristic reduces the communication energy through the NoC by approximating communicating tasks that exchange a high communication volume. This heuristic acts as follows in the three mapping steps:

• cluster selection: selects the cluster selected_cluster with the large number of available resources to map an application.

• initial task mapping: this heuristic evaluates the SPs inside the selected_cluster, selecting the SP with the largest region_free. The function region_free(sp i , n_hops) returns the total number of available resources of the set containing sp i and all SPs up to n_hops hops from sp i , as detailed in Figure 2. • non-initial task mapping: when a task t i is required to be mapped, the heuristic creates a list with all communicating tasks with t i already mapped onto the SPs within the cluster. Next, a bounding box rectangle is defined, containing all mapped communicating tasks. When there are more than one mapped communicating tasks, the bounding box is increased by one hop offering a larger search space to map t i . Figure 3 illustrates the mapping search space when one (a) or more communicating tasks (b) are mapped in the cluster. The adoption of a bounding box aims to reduce the distance among communicating tasks. In the case of just one mapped communicating task t j , this heuristic will map t i on a SP as close as possible to the one where t j is mapped. Otherwise, the heuristic will map t j onto the SP with the lowest communication energy cost according to the volume based energy model proposed by Hu et al [START_REF] Hu | An Efficient Power-Aware Optimization for Task Scheduling on NoC-based Many-core System[END_REF]. In both cases, if there are no available SPs inside the bounding box, it is increased by one hop. This heuristic also uses a premap method with the goal to group tasks that exchanges a high volume of communication in the same SP. The premap defines tasks' addresses before the request to map them. The premap is invoked when the SP is available, and a task t i is assigned to it. In such a case, the nonmapped tasks communicating with t i are assigned to the same processor, according to the communication volume.

2 0 1 1 1 2 2 3 3 2 3 2 3

Proposed Load (L) heuristic

This heuristic uses the task load Ld to distribute tasks and balance processor load of the system. For this purpose, this heuristic assumes that each sp i ∈ SP has a total computational load (TL i ), which corresponds to the sum of loads (Ld) of all executed tasks and the tasks that are currently being executed on this processor. Whenever a task is mapped onto sp i , the TL i is updated. The three steps of the mapping process for this heuristic are:

• cluster selection: this heuristic first computes the total computational load load(c i ) of each cluster c i by summing the TL of each SP inside the cluster. Then, the cluster with the smallest load(c i ) is selected. • initial task mapping: this heuristic evaluates all SPs inside the cluster and maps an initial task on the SP with lowest load Ld. • Non-initial task mapping: in this step, this heuristic acts as in the initial task mapping, choosing a SP inside the cluster with the lowest Ld.

Proposed Load-Communication (LC) heuristic

This heuristic makes a trade-off between processor load and reduction of communication volume. Thus, this heuristic also assumes each sp i ∈ SP has a total computational load (TL i ).

• cluster selection: this heuristic selects the cluster with the smallest load(c i ), as the previous heuristic. • initial task mapping: this heuristic divides the initial task process into two phases. The first phase selects the SP selected_sp with the smallest region_load to map an initial task. The region_load(sp i , n_hops) returns the average TL from the set containing an sp i and all SPs up to n_hops hops from sp i , as detailed in Figure 4. If there is more than one initial task, it is created a set with all SPs up to n_hops from the selected_sp, selecting the SP of this set with the smallest load. • Non-initial task mapping: this heuristic uses a similar bounding box search method used in the Premap-DN. The only difference is that the bounding box is increased by one hop in both cases: when there is one and when there are more than one communicating tasks mapped in the cluster. Figure 5 illustrates the mapping search space in the cluster.

In both cases, this heuristic selects the SP inside the bounding box with the lowest TL.

Results

This Section employs an OVP [START_REF] Ovp | [END_REF] NoC-based MPSoC platform model to compare the mapping heuristics in terms of processor load distribution, communication volume, and total execution time. For this, purpose, five applications are used as benchmarks:

• DTW -Digital Time Warping (DTW), with 10 tasks;

• MPEG decoder, with 5 tasks;

• DJK -Dijkstra, with 6 tasks; • SYN1, synthetic application, with 12 tasks, which emulates the communication behavior of an MPEG4 full decoder; • SYN2, synthetic application, with 12 tasks, that emulates the communication behavior of VOP (Video Object Plane) decoder application. Six different scenarios are evaluated, varying the number of applications and using an 8x8 MPSoC instance with 4x4 cluster size, as shown in Table 1.

Load Distribution

Figure 6 presents the number of instructions executed by each SP, for the three heuristics. Note Figure 6(a) illustrates the presence of hotspots (yellow squares), i.e., some PEs execute a larger number of instructions, while others execute a small number of instructions. Figure 6(b) and Figure 6(c) present a uniform load distribution for the LC and L heuristics. Such different behavior highlights the benefits of using heuristics that consider the load in the mapping process to avoid hotspots, therefore increasing the system lifetime.

Table 2 presents the average number of instructions executed per SP and its standard deviation for all evaluated scenarios. Results show that the average number of executed instructions per SP, for the same scenario, has small variations for the different heuristics. This is an expected result, since the workload is the same.

Most relevant result in Table 2 is the standard deviation value. A small standard deviation value represents a better load distribution among the SPs. The L and LC heuristics reduces the standard deviation in average 68% and 72% respectively when compared to the PREMAP-DN heuristic. The reason explaining this result is the fact that the PREMAP-DN heuristic does not take into account the tasks' load in the mapping decision. In addition, the PREMAP-DN heuristic tends to group tasks together in the same SP, eliminating communication through the NoC. Furthermore, the L heuristic has an average reduction of 13% compared to LC, since it has a larger search space searching inside the whole cluster.

Communication Volume

Table 3 shows the total communication volume in flits transferred through the NoC for each scenario. The PREMAP-DN heuristic presents the large reduction in the communication volume, being used as the reference mapping. The LC and L heuristics increase the communication volume in 81 and 154% (average values), respectively, when compared to the PREMAP heuristic. This happens since these heuristics have as main function to distribute tasks along the SPs, inducing the use of the NoC. Otherwise, as explained before, the PREMAP-DN heuristic tends to map tasks together in the same SP. In this case, there is an intra-SP communication between tasks. Moreover, the LC heuristic reduces the total communication volume in average 29% compared to the L heuristic. The LC heuristic considers the distance of the communicating tasks in the mapping decision, explaining the obtained results. Overloaded processors

Execution Time

Table 4 compares the execution time of each scenario. Results show the LC heuristic reduces the total execution time by 13% and 8% (average values) when compared to the PREMAP-DN and L heuristics. This result comes from the load distribution and communication energy trade-off used by the LC heuristic.

The PREMAP-DN increases the execution time because it increases the processor sharing (communicating tasks are mapped in the same SP). Distributing tasks and using as much as possible SPs of the system at the same time reduces the execution time, minimizing the processor sharing. The L heuristic increases the total execution time for two reasons. First, due to the larger search space, the heuristic takes longer to execute. Second, by not considering the NoC traffic in the mapping cost function, a large communication volume transferred through the NoC induces congestion and increases latency.

Reliability

The reliability of a given system is proportional to its use [START_REF] Chandra | Quantifying workload dependent reliability in embedded processors[END_REF], i.e. processors with the highest usage have more probability to suffer failures. The same happens with NoC links: an increased use of a given link will reduce its lifetime. The PREMAP-DN heuristic compromises processors reliability due to its behavior that tries to group tasks at the same SP to reduce communication energy. At the same way, the L heuristic compromises links reliability by increasing communication volume in 154% and 29% compared to the PREMAP-DN and LC heuristics, respectively.

Conclusion and Future Work

This work proposed a task mapping heuristic that makes a trade-off between processors load and communication volume -LC heuristic. This heuristic is compared with two heuristics: one reduces the communication volume and energy, and another evenly distributes the system workload. Results show the proposed heuristic reduces communication volume without compromising processors workload balance. Further, the proposed heuristic reduces the total execution time in average 10% compared to the other heuristics.

Future works aim to incorporate a reliability model in the MPSoC architecture to estimate the system mean-time-tofailure (MTTF).

Figure 1 :

 1 Figure 1: 6x6 MPSoC instance, with four 3x3 clusters [12].

Figure 2 :

 2 Figure 2: Hypothetical example of region_free(sp 2,2 , 2) in a 5x5 cluster, being sp i the central SP. The numbers inside at each SP corresponds to the number of available resources. The region load(sp 2,2 , 2) is then 25.

Figure 3 :

 3 Figure 3: (a) search space when one communicating tasks is already mapped (t 1 ); (b) search space when more than one communicating task is already mapped (t 1 and t 2 ). Thicker lines correspond to the original bounding box, dashed lines to the bounding box increased by one hop.

Figure 4 :

 4 Figure 4: Hypothetical example of region load in a 7x7 cluster, being sp i the central SP, and n_hops=3 (the number of SPs inside the region is 25): region load(sp 3,3 , 3) = 4100/25=164. The numbers inside at each SP corresponds to the processor load.

Figure 5 :

 5 Figure5: (a) search space when one communicating tasks is already mapped (t 1 ); (b) search space when more than one communicating task is already mapped (t 1 and t 2 ). Thicker lines correspond to the original bounding box, dashed lines to the bounding box increased by one hop.

Figure 6 :

 6 Figure 6: Thousands of instructions executed by each PE, for the 3 heuristics, scenario 6. Black square are manager PEs, not considered in the evaluation.

Table 1 :

 1 Evaluated scenarios.

	Scenario	Applications	Number of Applications	Number of tasks
	1	8 x DTW, 8 x MPEG, 8 x DJK, 8 x SYN1, 8 x SYN2	40	360
	2	8 x DJK, 11 x SYN1, 11 x SYN2	30	312
	3	21 x MPEG, 12 x DJK, 5 x SYN1, 12 x SYN2	50	381
	4	4 x DTW, 6 x MPEG, 14 x DJK, 7 x SYN1, 13 x SYN2	44	394
	5	8 x DTW, 11 x MPEG, 10 x DJK, 11 x SYN1, 9 x SYN2	49	435
	6	12 x DTW, 5 x MPEG, 8 x DJK, 13 x SYN1, 12 x SYN2	50	493

Table 2 :

 2 Average number of instructions per SP and standard deviation on the number of executed instructions per SP (thousands of instructions).

		PREMAP-DN	LC		L
	Scenario Avg.1 2,812	3,165	2,790 1,153	2,783	873
	2	2,402	3,402	2,413	777	2,422	805
	3	2,326	2,662	2,348	923	2,363	868
	4	3,055	3,497	3,026 1,266	3,064 1,115
	5	3,392	3,888	3,412 1,260	3,433 1,081
	6	3,914	4,189	3,814 1,174	3,829	836

/SP Std. Dev. Avg./SP Std. Dev. Avg./SP Std. Dev.

Table 3 :

 3 Total communication volume (thousands of flits) and the difference compared to the reference mapping.

	Scenario PREMAP-DN	LC	L
		Reference	Diff. %	Diff. %
	1	1,771	3,048 +72.15 4,128 +133.12
	2	1,417	3,070 +73.34 4,134 +133.42
	3	1,982	2,515 +42.05 3,444 +94.49
	4	1,956	2,761 +55.91 4,043 +128.32
	5	2,143	3,803 +114.78 5,637 +218.34
	6	2,528	4,057 +129.08 5,688 +221.22

Table 4 :

 4 Execution time (thousands of clock cycles).

	Scenario PREMAP-DN	LC	L
	1	18,520	16,514	18,854
	2	18,939	16,893	18,230
	3	17,083	14,866	17,409
	4	23,923	18,933	18,977
	5	21,139	18,927	20,322
	6	25,998	22,382	24,008

Acknowledgements

The Authors would like to thank Imperas Software Ltd. (www.imperas.com) and Open Virtual Platforms (www.OVPworld.org) for their support and access to their models and simulator. Fernando Moraes is supported by CNPq projects 472126/2013-0 and 302625/2012-7, and CAPES project CAPES-COFECUB 708/11.