V Konakov
email: vkonakov@hse.ru

A Kozhina
email: aakozhina@hse.ru

S Menozzi

Stability of Densities for Perturbed Diffusions and Markov Chains

Keywords: 1991 Mathematics Subject Classification. Primary 60H10; Secondary 65C30. November 25, 2016 Diffusion Processes, Markov Chains, Parametrix, Hölder Coefficients, bounded drifts

  L'archive ouverte pluridisciplinaire

Introduction

Setting.

For a fixed given deterministic final horizon T > 0, let us consider the following multidimensional SDE:

dX t = b(t, X t )dt + σ(t, X t )dW t , t ∈ [0, T ],
(1.1)

where b : [0, T ] × R d → R d , σ : [0, T ] × R d → R d ⊗ R d
are bounded coefficients that are measurable in time and Hölder continuous in space (this last condition will be possibly relaxed for the drift term b) and W is a Brownian motion on some filtered probability space (Ω, F , (F t ) t≥0 , P). Also, a(t, x) := σσ * (t, x) is assumed to be uniformly elliptic. In particular those assumptions guarantee that (1.1) admits a unique weak solution, see e.g. Bass and Perkins [START_REF] Bass | A new technique for proving uniqueness for martingale problems[END_REF], [START_REF] Menozzi | Parametrix techniques and martingale problems for some degenerate Kolmogorov equations[END_REF] from which the uniqueness to the martingale problem for the associated generator can be derived under the current assumptions.

We now introduce, for a given parameter ε > 0, a perturbed version of (1.1) with dynamics:

dX (ε) t = b ε (t, X (ε)
t )dt + σ ε (t, X (ε)
t )dW t , t ∈ [0, T ],

(1.2)

where b ε : [0, T ] × R d → R d , σ ε : [0, T ] × R d → R d ⊗ R d
satisfy at least the same assumptions as b, σ and are in some sense meant to be close to b, σ when ε is small. It is known that, under the previous assumptions, the density of the processes (X t ) t≥0 , (X

t ) t≥0 exists in positive time and satisfies some Gaussian bounds, see e.g Aronson [START_REF] Aronson | The fundamental solution of a linear parabolic equation containing a small parameter[END_REF], Sheu [START_REF] Sheu | Some estimates of the transition density of a nondegenerate diffusion Markov process[END_REF] or [START_REF] Delarue | Density estimates for a random noise propagating through a chain of differential equations[END_REF] for extensions to some degenerate cases.

The goal of this work is to investigate how the closeness of (b ε , σ ε ) and (b, σ) is reflected on the respective densities of the associated processes. Important applications can for instance be found in mathematical finance. If the dynamics of (1.1) models the evolution of the (log-)price of a financial asset, it is often very useful to know how a perturbation of the volatility σ impacts the density, and therefore the associated option prices (see also Corielli et al. [CFP10] or Benhamou et al. [START_REF] Benhamou | Expansion formulas for European options in a local volatility model[END_REF] for related problems).

In the framework of parameter estimation it can be useful, having at hand estimators (b ε , σ ε ) of the true parameters (b, σ) and some controls for the differences |bb ε |, |σσ ε | in a suitable sense, to quantify the difference p εp of the densities corresponding respectively to the dynamics with the estimated parameters and the one of the model.

Another important application includes the case of mollification by spatial convolution. This specific kind of perturbation is useful to investigate the error between the densities of a non-degenerate diffusion of type (1.1) with Hölder coefficients (or with piecewise smooth bounded drift) and its Euler scheme. In this framework, some explicit convergence results can be found in [START_REF] Konakov | Weak error for the Euler scheme of a diffusion with non-regular coefficients[END_REF].

More generally, this situation can appear in every applicative field for which the diffusion coefficient might be misspecified.

The previously mentioned Gaussian bounds on the density are derived through the so-called parametrix expansion which will be the crux of our approach. Roughly speaking, it consists in approximating the process by a proxy which has a known density, here a Gaussian one, and then in investigating the difference through the Kolmogorov equations. Various approaches to the parametrix expansion exist, see e.g. Il'in et al. [START_REF] Il'in | Second-order linear equations of parabolic type[END_REF], Friedman [START_REF] Friedman | Partial Differential Equations of Parabolic Type[END_REF] and McKean and Singer [START_REF] Mckean | Curvature and the eigenvalues of the Laplacian[END_REF]. The latter approach will be the one used in this work since it appears to be the most adapted to handle coefficients with no a priori smoothness in time and can also be directly extended to the discrete case for Markov chain approximations of equations (1.1) and (1.2). Let us mention in this setting the works of Konakov et al., see [START_REF] Konakov | Local limit theorems for transition densities of Markov chains converging to diffusions[END_REF], [START_REF] Konakov | Edgeworth type expansions for Euler schemes for stochastic differential equations[END_REF].

Our stability results will also apply to two Markov chains with respective dynamics:

Y t k+1 = Y t k + b(t k , Y t k )h + σ(t k , Y t k ) √ hξ k+1 , Y 0 = x, Y (ε) t k+1 = Y (ε) t k + b ε (t k , Y (ε)
t k )h + σ ε (t k , Y (ε)
t k ) √ hξ k+1 , Y (ε) 0 = x, (1.3)
where h > 0 is a given time step, for which we denote for all k ≥ 0, t k := kh and the (ξ k ) k≥1 are centered i.i.d. random variables satisfying some integrability conditions. Again, the key tool will be the parametrix representation for the densities of the chains and the Gaussian local limit theorem.

Assumptions and Main Results.

Let us introduce the following assumptions. Below, the parameter ε > 0 is fixed and the constants appearing in the assumptions do not depend on ε.

(A1) (Boundedness of the coefficients). The components of the vector-valued functions b(t, x), b ε (t, x) and the matrix-functions σ(t, x), σ ε (t, x) are bounded measurable. Specifically, there exist constants K 1 , K 2 > 0 s. (A2) (Uniform Ellipticity). The matrices a := σσ * , a ε := σ ε σ * ε are uniformly elliptic, i.e. there exists

Λ ≥ 1, ∀(t, x, ξ) ∈ [0, T ] × (R d ) 2 , Λ -1 |ξ| 2 ≤ a(t, x)ξ, ξ ≤ Λ|ξ| 2 , Λ -1 |ξ| 2 ≤ a ε (t, x)ξ, ξ ≤ Λ|ξ| 2 .
(A3) (Hölder continuity in space). For some γ ∈ (0, 1] , κ < ∞, for all t ∈ [0, T ],

|σ(t, x) -σ(t, y)| + |σ ε (t, x) -σ ε (t, y)| ≤ κ |x -y| γ .
Observe that the last condition also readily gives, thanks to the boundedness of σ, σ ε that a, a ε are also uniformly γ-Hölder continuous.

For a given ε > 0, we say that assumption (A) holds when conditions (A1)-(A3) are in force. Let us now introduce, under (A), the quantities that will bound the difference of the densities in our main results below. Set for ε > 0:

∆ ε,b,∞ := sup (t,x)∈[0,T ]×R d |b(t, x) -b ε (t, x)| , ∀q ∈ (1, +∞), ∆ ε,b,q := sup t∈[0,T ] b(t, .) -b ε (t, .) L q (R d ) .
Since σ, σ ε are both γ-Hölder continuous, see (A3) we also define

∆ ε,σ,γ := sup u∈[0,T ] |σ(u, .) -σ ε (u, .)| γ , (1.4)
where for γ ∈ (0, 1], |.| γ stands for the usual Hölder norm in space on C γ b (R d , R d ⊗ R d ) (space of Hölder continuous bounded functions, see e.g. Krylov [START_REF] Krylov | Lectures on elliptic and parabolic equations in Hölder spaces[END_REF]) i.e. :

|f | γ := sup x∈R d |f (x)| + [f ] γ , [f ] γ := sup x =y,(x,y)∈(R d ) 2 |f (x) -f (y)| |x -y| γ .
The previous control in particular implies for all (u, x, y)

∈ [0, T ] × (R d ) 2 : |a(u, x) -a(u, y) -a ε (u, x) + a ε (u, y)| ≤ 2(K 2 + κ)∆ ε,σ,γ |x -y| γ .
We eventually set for q ∈ (1, +∞], ∆ ε,γ,q := ∆ ε,σ,γ + ∆ ε,b,q , (1.5) which will be the key quantity governing the error in our results. We will denote, from now on, by C a constant depending on the parameters appearing in (A) and T . We reserve the notation c for constants that only depend on (A) but not on T . The values of C, c may change from line to line and do not depend on the considered parameter ε. Also, for given integers i, j ∈ N s.t. i < j, we will denote by [[i, j]] the set {i, i + 1, • • • , j}.

We are now in position to state our main results.

Theorem 1 (Stability Control for diffusions). Fix ε > 0 and a final deterministic time horizon T > 0. Under (A) and for q > d, there exist constants C := C(q) ≥ 1 and c ∈ (0, 1] s.t. for all 0 ≤ s < t ≤ T, (x, y) ∈ (R d ) 2 :

p c (t -s, y -x) -1 |(p -p ε )(s, t, x, y)| ≤ C∆ ε,γ,q , (1.6)
where p(s, t, x, .), p ε (s, t, x, .) respectively stand for the transition densities at time t of equations (1.1), (1.2) starting from x at time s. Also, we denote for a given c > 0 and for all u > 0,

z ∈ R d , p c (u, z) := c d/2 (2πu) d/2 exp(-c |z| 2 2u ).
Remark 1 (About the constants). We mention that the constant C := C(q) in (1.6) explodes when q ↓ d and is decreasing in q (see Lemmas 2 and 3 below for details). In particular, it can be chosen uniformly as soon as q ≥ q 0 > d.

Before stating our results for Markov Chains we introduce two kinds of innovations in (1.3). Namely: (I G ) The i.i.d. random variables (ξ k ) k≥1 are Gaussian, with law N (0, I d ). In that case the dynamics in (1.3) corresponds to the Euler discretization of equations (1.1) and (1.2). (I P,M ) For a given integer M > 2d + 5 + γ, the innovations (ξ k ) k≥1 are centered and have C 5 density f ξ which has, together with its derivatives up to order 5, at least polynomial decay of order M . Namely, for all z ∈ R d and multi-index ν, |ν| ≤ 5:

|D ν f ξ (z)| ≤ CQ M (z), (1.7)
where we denote for all

r > d, z ∈ R d , Q r (z) := c r 1 (1+|z|) r , R d dzQ r (z) = 1.
Theorem 2 (Stability Control for Markov Chains). Fix ε > 0 and a final deterministic time horizon T > 0.

For h = T /N, N ∈ N * := N \ {0}, we set for i ∈ N, t i := ih. Under (A), assuming that either (I G ) or (I P,M ) holds, and for q > d there exist

C := C(q) ≥ 1, c ∈ (0, 1] s.t. for all 0 ≤ t i < t j ≤ T, (x, y) ∈ (R d ) 2 : χ c (t j -t i , y -x) -1 |(p h -p h ε )(t i , t j , x, y)| ≤ C∆ ε,γ,q , (1.8)
where p h (t i , t j , x, .), p h ε (t i , t j , x, .) respectively stand for the transition densities at time t j of the Markov Chains Y and Y (ε) in (1.3) starting from x at time t i . Also:

-If (I G ) holds: χ c (t j -t i , y -x) := p c (t j -t i , y -x),
with p c as in Theorem 1.

-If (I P,M ) holds:

χ c (t j -t i , y -x) := c d (t j -t i ) d/2 Q M-(d+5+γ) |y -x| (t j -t i ) 1/2 /c .
The behavior of C := C(q) is similar to what is described in Remark 1.

1.3. On Some Related Applications.

1.3.1. Model Sensitivity for Option Prices.

Assume for instance that the (log)-price of a financial asset is given by the dynamics in (1.1). Under suitable assumptions the price of an option on that asset writes at time t and when X t = x as E[f (exp(X t,x T ))] up to an additional discounting factor. In the previous expression f is the pay-off function. For a rather large class of pay-offs, say measurable functions with polynomial growth, including irregular ones, Theorem 1 allows to specifically quantify how a perturbation of the coefficients impacts the option prices. Precisely for a given ε > 0, under (A):

|E ε (t, T, x, f )| := |E[f (exp(X t,x T ))] -E[f (exp(X (ε),t,x T ))]| ≤ C∆ ε,γ,q R d f (exp(y))p c (T -t, x, y)dy.
This previous control can be as well exploited to investigate perturbations of a model which provides some closed formulas, e.g. a perturbation of the Black and Scholes model that would include a stochastic volatility taking for instance σ ε (x) = σ + εψ(x) for some bounded γ-Hölder continuous function ψ and ε small enough.

In that case, assuming that the drift is known and unperturbed, we have ∆ ε,γ,∞ = |σ ε -σ| γ = ε|ψ| γ .

In connection with this application, we can quote the work of Corielli et al. [START_REF] Corielli | Parametrix approximation of diffusion transition densities[END_REF] who give estimates on option prices through parametrix expansions truncating the series. Some of their results, see e.g. their Theorem 3.1, can be related to a perturbation analysis since they obtain an approximation of an option price for a local volatility model in terms of the Black-Scholes price and a correction term corresponding to the first order term in the parametrix series. A more probabilistic approach to similar problems can be found in Benhamou et al. [START_REF] Benhamou | Expansion formulas for European options in a local volatility model[END_REF]. However, none of the indicated works indeed deals with the global perturbation analysis we perform here.

Weak Error Analysis

It is well known that if the coefficients b, σ in (1.1) are smooth and a satisfies the non-degeneracy condition (A2), then the weak error on the densities for the approximation by the Euler scheme is well controlled. Precisely, for a given time step h > 0, let us set for i ∈ N, t i := ih. Introduce now, for a given starting time t i , the Euler scheme

X h ti = x, ∀j > i, X h tj+1 = X h tj + b(t j , X h tj )h + σ(t j , X h tj )(W tj+1 -W tj
) and denote by p h (t i , t j , x, .) its density at time t j . The dynamics of the Euler scheme clearly enters the scheme (1.3). It can be derived from Konakov and Mammen [START_REF] Konakov | Edgeworth type expansions for Euler schemes for stochastic differential equations[END_REF] (see also Bally and Talay [START_REF] Bally | The law of the Euler scheme for stochastic differential equations, II. Convergence rate of the density[END_REF] for an extension to the hypoelliptic setting) that:

|pp h |(t i , t j , x, y) ≤ Chp c (t jt i , yx). If the coefficients in (1.1) are not smooth, it is then possible to use a mollification procedure, taking for

x ∈ R d , b ε (t, x) := b(t, .) ⋆ ρ ε (x), σ ε (t, x) := σ(t, .) ⋆ ρ ε (x) with ρ ε (x) := ε -d ρ(x/ε) and ρ ∈ C ∞ (R d , R + ), R d ρ(x)dx = 1, supp(ρ) ⊂ K for some compact set K of R d .
For the mollifying kernel ρ ε , one then easily checks that for γ-Hölder continuous in space coefficients b, σ there exists C s.t.

sup t∈[0,T ] |b(t, .) -b ε (t, .)| ∞ ≤ Cε γ , sup t∈[0,T ] |σ(t, .) -σ ε (t, .)| η ≤ Cε γ-η , η ∈ (0, γ).
(1.9)

The important aspect is that we lose a bit with respect to the sup norm when investigating the Hölder norm.

We then have by Theorems 1 and 2 and their proof, that, for γ-Hölder continuous in space coefficients b, σ and taking q = ∞, there exist c, C s.t. for all 0 ≤ s < t ≤ T, 0 ≤ t i < t j ≤ T, (x, y) ∈ (R d ) 2 :

|(p -p ε )(s, t, x, y)| ≤ CC η ε γ-η p c (t -s, y -x), |(p h -p h ε )(t i , t j , x, y)| ≤ CC η ε γ-η p c (t j -t i , y -x),
where the constant C η explodes when η tends to 0.

To investigate the global weak error (pp h )(t i , t j , x, y) = {(pp ε ) + (p εp h ε ) + (p h εp h )}(t i , t j , x, y), it therefore remains to analyze the contribution (p εp h ε )(t i , t j , x, y). The results of [START_REF] Konakov | Edgeworth type expansions for Euler schemes for stochastic differential equations[END_REF] indeed apply but yield

|(p ε -p h ε )(t i , t j , x, y)| ≤ C ε hp c (t j -t i , y -x)
where C ε is explosive when ε goes to zero. The global error thus writes:

|(p -p h )(t i , t j , x, y)| ≤ C{C η ε γ-η + C ε h}p c (t j -t i , y -x)
, and a balance is needed to derive a global error bound. This is precisely the analysis we perform in [START_REF] Konakov | Weak error for the Euler scheme of a diffusion with non-regular coefficients[END_REF]. In this work, we extend to densities (up to a slowly growing factor) the results previously obtained by Mikulevičius and Platen [START_REF] Mikulevičius | Rate of convergence of the Euler approximation for diffusion processes[END_REF] on the weak error, i.e. they showed

|E[f (X T ) -f (X h T )]| ≤ Ch γ/2 provided f ∈ C 2+γ b (R d , R
). Precisely, we obtain through a suitable analysis of the constants C η , C ε , which respectively depend on behavior of the parametrix series and of the derivatives of the heat kernel with mollified coefficients, that |p-p h |(t i , t j , x, y) ≤ Ch γ/2-ψ(h) p c (t jt i , yx) for a function ψ(h) going to 0 as h → 0 (which is induced by the previous loss of η in (1.9)). In the quoted work, we also obtain some error bounds for piecewise smooth drifts having a countable set of discontinuities. This part explicitly requires the stability result of Theorems 1, Theorems 2 for q < +∞. The idea being that the difference between the piece-wise smooth drift and its smooth approximation (actually the mollification procedure is only required around the points of discontinuity), is well controlled in L q norm, q < +∞.

Extension to some Kinetic Models

The results of Theorems 1 and 2 should extend without additional difficulties to the case of degenerate diffusions of the form:

dX 1 t = b(t, X t )dt + σ(t, X t )dW t , dX 2 t = X 1 t dt,
(1.10) denoting X t = (X 1 t , X 2 t ), under the same previous assumptions on b, σ when we consider perturbations of the non-degenerate components, i.e. for a given ε > 0, X

(ε) t = (X 1,(ε) t , X 2,(ε) t
) where:

dX 1,(ε) t = b ε (t, X (ε) t )dt + σ ε (t, X (ε) t )dW t , dX 2,(ε) t = X 1,(ε) t dt.
(1.11) Indeed, under (A), the required parametrix expansions of the densities associated with the solutions of equation (1.10), (1.11) have been established in [START_REF] Konakov | Explicit parametrix and local limit theorems for some degenerate diffusion processes[END_REF].

A posteriori Controls in Parameter Estimation.

Let us consider to illustrate this application a parametrized family of diffusions of the form:

dX t = b(t, X t )dt + σ(ϑ, t, X t )dW t , (1.12)
where ϑ ∈ Θ ⊂ R d , the coefficients b, σ are smooth, bounded and the non-degeneracy condition (A2) holds.

A natural practical problem consists in estimating the true parameter ϑ from an observed discrete sample (X

t n i ) i∈[[0,n]] where the {(t i ) i∈[[0,n]] } form a partition of the observation interval, i.e. if T = 1, 0 = t n 0 < t n 1 < • < t n n = 1.
Introducing the contrast:

U n (ϑ) := 1 n n i=1 log det(a(ϑ, t n i-1 , X t n i-1 )) + a -1 (ϑ, t n i-1 , X t n i-1 )X n i , X n i , ∀i ∈ [[1, n]], X n i := X t n i -X t n i-1 t n i -t n i-1
, and denoting by θn the corresponding minimizer, it was shown by Genon-Catalot and Jacod [GCJ93] that under P ϑ , √ n( θnϑ) converges in law towards a mixed normal variable S which is, conditionally to F 1 := σ{(X s ) s∈[0,1] }, centered and Gaussian. For a precise expression of the covariance which depends on the whole path of (X t ) t∈[0,1] we refer to Theorem 3 and its proof in [START_REF] Genon-Catalot | On the estimation of the diffusion coefficient for multi-dimensional diffusion processes[END_REF].

This means that, when n is large, conditionally to F 1 , we have on a subset Ω ⊂ Ω which has high probability,

that | θn -ϑ| ≤ C √ n for a certain threshold C. Setting ε n = n -1/2 , σ εn (t, x) := σ( θn , t, x
) and with a slight abuse of notation σ(t, x) := σ(ϑ, t, x), one gets that, on Ω:

|σ(t, x) -σ εn (t, x) -(σ(t, y) -σ εn (t, y))| ≤ |x -y| ∧ Cn -1/2 ⇒ |σ -σ ε | η ≤ (Cn -1/2 ) 1-η , η ∈ (0, 1].
We can then invoke our Theorem 1 to compare the densities of the diffusions with the estimated parameter and the exact one in (1.12).

The paper is organized as follows. We recall in Section 2 some basic facts about parametrix expansions for the densities of diffusions and Markov Chains. We then detail in Section 3 how to perform a stability analysis of the parametrix expansions in order to derive the results of Theorems 1 and 2.

Derivation of formal series expansion for densities

Parametrix Representation of the Density for Diffusions

In the following, for given (s, x) ∈ R + × R d , we use the standard Markov notation (X s,x t ) t≥s to denote the solution of (1.1) starting from x at time s.

Assume that (X s,x t ) t≥s has for all t > s a smooth density p(s, t, x, .) (which is the case if additionally to (A) the coefficients are smooth see e.g. Friedman [START_REF] Friedman | Partial Differential Equations of Parabolic Type[END_REF]). We would like to estimate this density at a given point y ∈ R d . To this end, we introduce the following Gaussian inhomogeneous process with spatial variable frozen at y. For all (s, x) ∈ [0, T ] × R d , t ≥ s we set:

Xy t = x + t s σ(u, y)dW u .
Its density py readily satisfies the Kolmogorov backward equation:

   ∂ u py (u, t, z, y) + Ly u py (u, t, z, y) = 0, s ≤ u < t, z ∈ R d , py (u, t, ., y) → u↑t δ y (.), (2.1)
where for all ϕ ∈ C 2 0 (R d , R) (twice continuously differentiable functions with compact support) and z ∈ R d :

Ly u ϕ(z) = 1 2 Tr σσ * (u, y)D 2 z ϕ(z) ,
stands for the generator of Xy at time u.

On the other hand, since we have assumed the density of X to be smooth, it must satisfy the Kolmogorov forward equation (see e.g. Dynkin [START_REF] Dynkin | Markov Processes[END_REF]). For a given starting point

x ∈ R d at time s,    ∂ u p(s, u, x, z) -L * u p(s, u, x, z) = 0, s < u ≤ t, z ∈ R d , p(s, u, x, .) → u↓s δ x (.), (2.2)
where L * u stands for the formal adjoint (which is again well defined if the coefficients in (1.1) are smooth) of the generator of (1.1) which for all ϕ ∈ C 2 0 (R d , R), z ∈ R d writes:

L u ϕ(z) = 1 2 Tr σσ * (u, z)D 2 z ϕ(z) + b(u, z), D z ϕ(z) .
Equations (2.1), (2.2) yield the formal expansion below which is initially due to McKean and Singer [START_REF] Mckean | Curvature and the eigenvalues of the Laplacian[END_REF].

(p -py )(s, t, x, y) = t s du∂ u R d dzp(s, u, x, z)p y (u, t, z, y) = t s du R d dz (∂ u p(s, u, x, z)p y (u, t, z, y) + p(s, u, x, z)∂ u py (u, t, z, y)) = t s du R d dz L * u p(s, u, x, z)p y (u, t, z, y) -p(s, u, x, z) Ly u py (u, t, z, y) = t s du R d dzp(s, u, x, z)(L u -Ly u )p y (u, t, z, y), (2.3)
using the Dirac convergence for the first equality, equations (2.2) and (2.1) for the third one. We eventually take the adjoint for the last equality. Note carefully that the differentiation under the integral is also here formal since we would need to justify that it can actually be performed using some growth properties of the density and its derivatives which we a priori do not know.

Let us now introduce the notation

f ⊗ g(s, t, x, y) = t s du R d dzf (s, u, x, z)g(u, t, z, y)
for the time-space convolution and let us define p(s, t, x, y) := py (s, t, x, y), that is in p(s, t, x, y) we consider the density of the frozen process at the final point and observe it at that specific point. We now introduce the parametrix kernel: H(s, t, x, y) := (L s -Ls )p(s, t, x, y) := (L s -Ly s )p y (s, t, x, y).

(2.4)

With those notations equation (2.3) rewrites:

(p -p)(s, t, x, y) = p ⊗ H(s, t, x, y).
From this expression, the idea then consists in iterating this procedure for p(s, u, x, z) in (2.3) introducing the density of a process with frozen characteristics in z which is here the integration variable. This yields to iterated convolutions of the kernel and leads to the formal expansion:

p(s, t, x, y) = ∞ r=0 p ⊗ H (r) (s, t, x, y), (2.5)
where p ⊗ H (0) = p, H (r) = H ⊗ H (r-1) , r ≥ 1. Obtaining estimates on p from the formal expression (2.5) requires to have good controls on the right-hand side. The remarkable property of this formal expansion is now that the right-hand-side of (2.5) only involves controls on Gaussian densities which in particular will provide, associated with our assumption (A) a smoothing in time property for the kernel H.

Proposition 1. Under the sole assumption (A), for t > s, the density of X x,s t solving (1.1) exists and can be written as in (2.5).

Proof. The proof can already be derived from a sensitivity argument. We first introduce two parametrix series of the form (2.5). Namely, p(s, t, x, y) := p(s, t, x, y)

+ ∞ r=1 p ⊗ H (r) (s, t, x, y) (2.6)
and

p ε (s, t, x, y) := p ε (s, t, x, y) + ∞ r=1 pε ⊗ H (r) ε (s, t, x, y). (2.7)
Let us point out that, at this stage, p and p ε are defined as sums of series. The purpose is then to identify those sums with the densities of the processes X s,x t , X (ε),s,x t at point y. The convergence of the series (2.6) and (2.7) is in some sense standard (see e.g. [START_REF] Menozzi | Parametrix techniques and martingale problems for some degenerate Kolmogorov equations[END_REF] or Friedman [START_REF] Friedman | Partial Differential Equations of Parabolic Type[END_REF]) under (A). We recall for the sake of completeness the key steps for (2.6).

From direct computations, there exist c 1 ≥ 1, c ∈ (0, 1] s.t. for all T > 0 and all multi-index α, |α| ≤ 8,

∀0 ≤ u < t ≤ T, (z, y) ∈ (R d ) 2 , |D α z p(u, t, z, y)| ≤ c 1 (t -u) |α|/2 p c (t -u, y -z), (2.8)
where

p c (t -u, y -z) = c d/2 (2π(t -u)) d/2 exp - c 2 |y -z| 2 t -u ,
stands for the usual Gaussian density in R d with 0 mean and covariance (tu)c -1 I d . From (2.8), the boundedness of the drift and the Hölder continuity in space of the diffusion matrix we readily get that there exists

c 1 ≥ 1, c ∈ (0, 1], |H(u, t, z, y)| ≤ c 1 (1 ∨ T (1-γ)/2 ) (t -u) 1-γ/2 p c (t -u, z -y).
(2.9)

Now the key point is that the control (2.9) yields an integrable singularity giving a smoothing effect in time once integrated in space in the time-space convolutions appearing in (2.6) and (2.7). It follows by induction that:

| p ⊗ H (r) (s, t, x, y)| ≤ ((1 ∨ T (1-γ)/2 )c 1 ) r+1 r i=1 B( γ 2 , 1 + (i -1) γ 2 )p c (t -s, y -x)(t -s) rγ 2 = ((1 ∨ T (1-γ)/2 )c 1 ) r+1 Γ( γ 2 ) r Γ(1 + r γ 2 ) p c (t -s, y -x)(t -s) rγ 2 , (2.10) where for a, b > 0, B(a, b) = 1 0 t -1+a (1 -t) -1+b
dt stands for the β-function, and using as well the identity B(a, b) = Γ(a)Γ(b) Γ(a+b) for the last equality. These bounds readily yield the convergence of the series as well as a Gaussian upper-bound. Namely

p(s, t, x, y) ≤ c 1 exp((1 ∨ T (1-γ)/2 )c 1 [(t -s) γ/2 ])p c (t -s, y -x).
(2.11)

An important application of the stability of the perturbation consists in considering coefficients b

ε := b ⋆ ζ ε , σ ε := σ ⋆ ζ ε in (2.7)
, where ζ ε is a mollifier in time and space. For mollified coefficients which satisfy the non-degeneracy assumption (A2), the existence and smoothness of the density of X (ε) is a standard property (see e.g. Friedman [START_REF] Friedman | Partial Differential Equations of Parabolic Type[END_REF], [START_REF] Friedman | Stochastic differential equations[END_REF]). It is also well known, see [START_REF] Konakov | Local limit theorems for transition densities of Markov chains converging to diffusions[END_REF], that this density p ε admits the series representation given in (2.7). Observe now carefully that the previous Gaussian bounds also hold for p ε uniformly in ε and independently of the mollifying procedure. This therefore gives that p ε (s, t, x, y) -→ ε→0 p(s, t, x, y),

(2.12) boundedly and uniformly. Thus, for every continuous bounded function f we derive from the bounded convergence theorem and (2.11) that for all 0 ≤ s < t ≤ T, x ∈ R d :

E[f (X (ε),s,x t )] = R d f (y)p ε (s, t, x, y)dy -→ ε→0 R d f (y)p(s, t, x, y)dy.
(2.13)

Taking f = 1 gives that R d p(s, t, x, y)dy = 1 and the uniform convergence in (2.12) gives that p(s, t, x, .) is non negative. We therefore derive that p(s, t, x, •) is a probability density on R d .

On the other hand, under (A), we can derive from Theorem 11.3.

4 of [SV79] that (X ε s ) s∈[0,T ] (law) =⇒ ε→0 (X s ) s∈[0,T ] .
This gives that for any bounded continuous function f :

E[f (X (ε),s,x t )] -→ ε→0 E[f (X s,x t )].
This convergence and (2.13) then yield that the random variable X s,x t admits p(s, t, x, •) as density. We can thus now conclude that the processes X, X (ε) in (1.1), (1.2) have transition densities given by the sum of the series (2.6), (2.7).

Parametrix for Markov Chains

One of the main advantages of the formal expansion in (2.5) is that it has a direct discrete counterpart in the Markov chain setting. Indeed, denote by (Y ti,x tj ) j≥i the Markov chain with dynamics (1.3) starting from x at time t i . Observe first that if the innovations (ξ k ) k≥1 have a density then so does the chain at time t k .

Let us now introduce its generator at time t i , i.e. for all ϕ ∈ C 2 0 (R d , R), x ∈ R d :

L h ti ϕ(x) := h -1 E[ϕ(Y ti,x ti+1 ) -ϕ(x)].
In order to give a representation of the density of p h (t i , t j , x, y) of Y ti,x tj at point y for j > i, we introduce similarly to the continuous case, the Markov chain (or inhomogeneous random walk) with coefficients frozen in space at y. For given (t i , x) ∈ [0, T ] × R d , t j ≥ t i we set:

Ỹ ti,x,y tj := x + h 1/2 j-1 k=i σ(t k , y)ξ k+1 ,
and denote its density by ph,y (t i , t j , x, .). Its generator at time

t i writes for all ϕ ∈ C 2 0 (R d , R), x ∈ R d : Lh,y ti ϕ(x) = h -1 E[ϕ( Ỹ ti,x,y ti+1 ) -ϕ(x)].
Using the notation ph (t i , t j , x, y) := ph,y (t i , t j , x, y), we introduce now for 0 ≤ i < j ≤ N the parametrix kernel: y). Analogously to Lemma 3.6 in [START_REF] Konakov | Local limit theorems for transition densities of Markov chains converging to diffusions[END_REF], which follows from a direct algebraic manipulation, we derive the following representation for the density which can be viewed as the Markov chain analogue of Proposition 1.

H h (t i , t j , x, y) := (L h ti -Lh,y ti )p h (t i + h, t j , x,
Proposition 2 (Parametrix Expansion for the Markov Chain). Assume (A) is in force. Then, for 0 ≤ t i < t j ≤ T ,

p h (t i , t j , x, y) = j-i r=0 ph ⊗ h H h,(r) (t i , t j , x, y),
where the discrete time convolution type operator ⊗ h is defined by

f ⊗ h g(t i , t j , x, y) = j-i-1 k=0 h R d f (t i , t i+k , x, z)g(t i+k , t j , z, y)dz.
Also g ⊗ h H h,(0) = g and for all r ≥ 1, H h,(r) := H h ⊗ h H h,(r-1) denotes the r-fold discrete convolution of the kernel H h .

Stability of Parametrix Series.

We will now investigate more specifically the sensitivity of the density w.r.t. the coefficients through the difference of the series. For a given fixed parameter ε, under (A) the densities p(s, t, x, .), p ε (s, t, x, •) at time t of the processes in (1.1), (1.2) starting from x at time s both admit a parametrix expansion of the previous type.

Stability for Diffusions: Proof of Theorem 1

Let us consider the difference between two parametrix expansions:

|p(s, t, x, y) -p ε (s, t, x, y)| = | ∞ r=0 p ⊗ H (r) (s, t, x, y) - ∞ r=0 pε ⊗ H (r) ε (s, t, x, y)| ≤ |(p -pε )(s, t, x, y)| + | ∞ r=1 p ⊗ H (r) (s, t, x, y) - ∞ r=1 pε ⊗ H (r) ε (s, t, x, y)|. (3.1)
The strategy to study the above difference, using some well known properties of the Gaussian kernels and their derivatives recalled in (2.8), consists in first studying the difference of the main terms.

We have the following Lemma.

Lemma 1 (Difference of the first terms and their derivatives). Under (A), there exist c 1 ≥ 1, c ∈ (0, 1] s.t. for all 0 ≤ s < t, (x, y) ∈ (R d ) 2 and all multi-index α, |α| ≤ 4,

|D α x p(s, t, x, y) -D α x pε (s, t, x, y)| ≤ c 1 (t -s) |α|/2 ∆ ε,σ,γ p c (t -s, y -x).
Proof. Let us first consider |α| = 0 and introduce some notations. Set:

Σ(s, t, y) := t s a(u, y)du, Σ ε (s, t, y) := t s a ε (u, y)du. (3.2)
Let us now identify the columns of the matrices Σ(s, t, y), Σ ε (s, t, y) with d-dimensional column vectors, i.e. for Σ(s, t, y):

Σ(s, t, y) = Σ 1 Σ 2 • • • Σ d (s, t, y).
We now rewrite p(s, t, x, y) and pε (s, t, x, y) in terms of vectors Θ(s, t, y),

Θ ε (s, t, y) ∈ R d 2 : p(s, t, x, y) = f x,y (Θ(s, t, y)), Θ(s, t, y) = ((Σ 1 ) * , • • • , (Σ d ) * ) * (s, t, y), pε (s, t, x, y) = f x,y (Θ ε (s, t, y)), Θ ε (s, t, y) = ((Σ 1 ε ) * , • • • , (Σ d ε ) * ) * (s, t, y),
where (•) * stands for the transpose and

f x,y : R d 2 → R Γ → f x,y (Γ) = 1 (2π) d/2 det(Γ 1:d ) 1/2 exp - 1 2 (Γ 1:d ) -1 (y -x), y -x , (3.3)
where Γ :=

     Γ 1 Γ 2 . . . Γ d      and each (Γ i ) i∈[[1,d]]
belongs to R d . Also, we have denoted:

Γ 1:d := Γ 1 Γ 2 • • • Γ d , the d × d matrix formed with the entries (Γ i ) i∈[[1,d]]
, each entry being viewed as a column.

The multidimensional Taylor expansion now gives:

|(p -pε )(s, t, x, y)| = |f x,y (Θ(s, t, y)) -f x,y (Θ ε (s, t, y))| = |ν|=1 D ν f x,y (Θ(s, t, y)){(Θ ε -Θ)(s, t, y)} ν +2 |ν|=2 {(Θ ε -Θ)(s, t, y)} ν ν! 1 0 (1 -δ)D ν f x,y ([Θ + δ(Θ ε -Θ)](s, t, y))dδ , (3.4)
where for a multi-index

ν := (ν 1 , • • • , ν d 2 ) ∈ N d 2 , we denote by |ν| := d 2 i=1 ν i the length of the multi-index, ν! = d 2 i=1 ν i ! and for h ∈ R d 2 , h ν := d 2
i=1 h νi i (with the convention that 0 0 = 1). With these notations, from (3.2), (3.3), (3.4) we get:

|f x,y (Θ(s, t, y)) -f x,y (Θ ε (s, t, y))| ≤ c |ν|=1 |D ν f x,y (Θ(s, t, y))|∆ ε,σ,γ (t -s) + ∆ 2 ε,σ,γ (t -s) 2 max δ∈[0,1] |ν|=2 |D ν f x,y ([Θ + δ(Θ ε -Θ)](s, t, y))| , (3.5)
recalling that ∆ ε,σ,γ has previously been defined in (1.4) 1 . Since f x,y in (3.3) is a Gaussian density in the parameters x, y, we recall from Cramer and Leadbetter [CL04] (see eq. (2.10.3) therein), that for all Γ ∈ R d 2 and any multi index ν, |ν| ≤ 2:

D ν f x,y (Γ) = 1 2 |ν|   d 2 i=1 ( ∂ 2 ∂ x ⌊ i-1 d ⌋+1 ∂ x i-⌊ i-1 d ⌋d ) νi f x,y (Γ)   ,
where ⌊•⌋ stands for the integer part. Hence, taking from (3.5), for all δ ∈ [0, 1], Γ ε,δ (s, t, y) := [Θ + δ(Θ ε -Θ)](s, t, y) yields, thanks to the non-degeneracy conditions (see equation (2.8)):

|D ν f x,y (Γ ε,δ (s, t, y))| ≤ c1 (t -s) |ν| f x,y cΓ ε,δ (s, t, y) ≤ c1 (t -s) |ν| p c(t -s, y -x), (3.6)
for some c1 ≥ 1, c ∈ (0, 1]. Thus, from (3.3), (3.4), equations (3.5) and (3.6) give:

|p(s, t, x, y) -pε (s, t, x, y)| ≤ c1 ∆ ε,σ,γ p c(t -s, y -x).
Up to a modification of c1 , c or c 1 , c in (2.8) we can assume that the statement of the lemma and (2.8) hold with the same constants c 1 , c. The bounds for the derivatives are established similarly using the controls of (2.8). This concludes the proof.

Remark 2. Observe from equation (3.4) that the previous Lemma still holds with ∆ ε,σ,γ replaced by ∆ ε,σ,∞ := sup t∈[0,T ] |σ(t, .)σ ε (t, .)| ∞ . The Hölder norm is required to control the differences of the parametrix kernels.

The previous lemma quantifies how close are the main parts of the expansions. To proceed we need to consider the difference between the one-step convolutions. Combining the estimates of Lemmas 1 and 2 below will yield by induction the result stated in Theorem 1.

1 Actually only the sup norm term in ∆ε,σ,γ could be considered for this part of the anlalysis.

Lemma 2 (Control of the one-step convolution). For all 0 ≤ s < t ≤ T, (x, y) ∈ (R d ) 2 and for q ∈ (d, +∞]:

|p ⊗ H (1) (s, t, x, y) -pε ⊗ H (1) ε (s, t, x, y)| ≤ c 2 1 p c (t -s, y -x) 2(1 ∨ T (1-γ)/2 ) 2 [∆ ε,σ,γ + I q=+∞ ∆ ε,b,+∞ ]B(1, γ 2 )(t -s) γ 2 +I q∈(d,+∞) ∆ ε,b,q B( 1 2 + α(q), α(q))(t -s) α(q) , (3.7)
where c 1 , c are as in Lemma 1 and for q ∈ (d, +∞] we set α(q) = 1 2 (1 -d q ). The above control then yields for a fixed q ∈ (d, +∞]:

|p ⊗ H (1) (s, t, x, y) -pε ⊗ H (1) ε (s, t, x, y)| ≤ 2 C2 ∆ ε,γ,q p c (t -s, y -x)(t -s) γ 2 ∧α(q) (B(1, γ 2 ) ∨ B( 1 2 + α(q), α(q))), C = c 1 (1 ∨ T (1-γ)/2 ), (3.8)
which will be useful for the iteration (see Lemma 3).

Proof. Let us write:

|p ⊗ H (1) (s, t, x, y) -pε ⊗ H (1) ε (s, t, x, y)| ≤ |p -pε | ⊗ |H|(s, t, x, y) + pε ⊗ |H -H ε |(s, t, x, y) =: I + II. (3.9)
From Lemma 1 and (2.9) we readily get for all q ∈ (d, +∞]:

|p -pε | ⊗ |H|(s, t, x, y)| ≤ ((1 ∨ T (1-γ)/2 )c 1 ) 2 ∆ ε,γ,q p c (t -s, y -x)B(1, γ 2 )(t -s) γ 2 .
(3.10) Now we will establish that for all 0 ≤ u < t ≤ T, (z, y) ∈ (R d ) 2 and q = +∞:

|(H -H ε )(u, t, z, y)| ≤ ∆ ε,γ,∞ (1 ∨ T (1-γ)/2 )c 1 (t -u) 1-γ 2 p c (t -u, y -z).
(3.11) Equations (3.11) and (2.8) give that II can be handled as I which yields the result for q = +∞. It therefore remains to prove (3.11). Let us write with the notations of (3.3):

(H -H ε )(u, t, z, y) := 1 2 Tr (a(u, z) -a(u, y))D 2 z f z,y Θ(u, t, y) + b(u, z), D z f z,y Θ(u, t, y) - 1 2 Tr (a ε (u, z) -a ε (u, y))D 2 z f z,y Θ ε (u, t, y) + b ε (u, z), D z f z,y Θ ε (u, t, y) .
Thus,

(H -H ε )(u, t, z, y) = 1 2 Tr (a(u, z) -a(u, y)){D 2 z f z,y Θ(u, t, y) -D 2 z f z,y Θ ε (u, t, y) } -Tr [(a ε (u, z) -a ε (u, y) -(a(u, z) -a(u, y))]D 2 z f z,y Θ ε (u, t, y) + b(u, z), {D z f z,y Θ(u, t, y) -D z f z,y Θ ε (u, t, y) } -(b ε (u, z) -b(u, z)), D z f z,y Θ ε (u, t, y) . (3.12)
Observe now that, similarly to (3.6) one has for all i ∈ {1, 2}:

|D i z f z,y Θ(u, t, y) | + |D i z f z,y Θ ε (u, t, y) | ≤ c1 (t -u) i/2 p c(t -u, y -z), |D i z f z,y Θ(u, t, y) -D i z f z,y Θ ε (u, t, y) | ≤ c1 ∆ ε,σ,γ (t -u) i/2 p c(t -u, y -z).
Also,

|(a ε (u, z) -a ε (u, y) -(a(u, z) -a(u, y))| ≤ c∆ ε,σ,γ |z -y| γ , |b ε (u, z) -b(u, z)| ≤ c∆ ε,b,∞ .
Thus, provided that c 1 , c have been chosen large and small enough respectively in Lemma 1, the definition in (1.5) gives:

|(H -H ε )(u, t, z, y)| ≤ (1 ∨ T (1-γ)/2 )c 1 ∆ ε,γ,∞ (t -u) 1-γ/2 p c (t -u, y -z).
This establishes (3.11) for q = +∞. For q ∈ (d, +∞) we have to use Hölder's inequality in the time-space convolution involving the difference of the drifts (last term in (3.12)). Set:

D(s, t, x, y) := t s du R d pε (s, u, x, z) (b ε (u, z) -b(u, z)), D z f z,y Θ ε (u, t, y) dz.
Denoting by q the conjugate of q, i.e. q, q > 1, q -1 + q-1 = 1, we get from (2.8) and for q > d that:

|D(s, t, x, y)| ≤ c 2 1 t s du (t -u) 1/2 b(u, .) -b ε (u, .) L q (R d ) R d [p c (u -s, z -x)p c (t -u, y -z)] q dz 1/q ≤ c 2 1 ∆ ε,b,q t s c d (2π) d(1-1 q ) (cq) d/q R d p cq (u -s, z -x)p cq (t -u, y -z)dz 1/q du (u -s) d 2 (1-1 q ) (t -u) 1 2 + d 2 (1-1 q ) ≤ c 2 1 c(t -s) 2π d 2 (1-1 q ) q-d 2 q ∆ ε,b,q p c (t -s, y -x) t s du (u -s) d 2 (1-1 q ) (t -u) 1 2 + d 2 (1-1 q )
. Now, the constraint d < q < +∞ precisely gives that 1 < q < d/(d -1) ⇒ 1 2 + d 2 (1 -1 q ) < 1 so that the last integral is well defined. We therefore derive:

|D(s, t, x, y)| ≤ c 2 1 (t -s) 1 2 -d 2 (1-1 q ) ∆ ε,b,q p c (t -s, y -x)B(1 - d 2 (1 - 1 q ), 1 2 - d 2 (1 - 1 q )).
In the case d < q < +∞, recalling that α(q) = 1 2 (1 -d q ), we eventually get :

pε ⊗ |H -H ε |(s, t, x, y) ≤ c 2 1 p c (t -s, y -x){∆ ε,b,q (t -s) α(q) B( 1 2 + α(q), α(q)) +2∆ ε,σ,γ (1 ∨ T (1-γ)/2 )(t -s) γ/2 B(1, γ/2)}.
(3.13)

The statement now follows in whole generality from (3.9), (3.10), equations (3.11), (2.8) for q = ∞ and (3.13) for d < q < +∞.

The following Lemma associated with Lemmas 1 and 2 allows to complete the proof of Theorem 1.

Lemma 3 (Difference of the iterated kernels). For all 0 ≤ s < t ≤ T, (x, y) ∈ (R d ) 2 and for all q ∈ (d, +∞], r ∈ N:

|(p ⊗ H (r) -pε ⊗ H (r) ε )(s, t, x, y)| ≤ (r + 1)∆ ε,γ,q Cr+1 Γ( γ 2 ∧ α(q)) r Γ(1 + r( γ 2 ∧ α(q))) p c (t -s, y -x)(t -s) r( γ 2 ∧α(q)) . (3.14)
where c is as in Lemma 1 and C as in equation (3.8) of Lemma 2.

Proof. Observe that Lemmas 1 and 2 respectively give (3.14) for r = 0 and r = 1. Let us assume that it holds for a given r ∈ N * and let us prove it for r + 1.

Let us denote for all r ≥ 1, η r (s, t, x, y)

:= |(p ⊗ H (r) -pε ⊗ H (r) ε )(s, t, x, y)|. Write η r+1 (s, t, x, y) ≤ |[p ⊗ H (r) -pε ⊗ H (r) ε ] ⊗ H(s, t, x, y)| + |p ε ⊗ H (r) ε ⊗ (H -H ε )(s, t, x, y)| ≤ η r ⊗ |H|(s, t, x, y) + |p ε ⊗ H (r) ε | ⊗ |(H -H ε )|(s, t, x, y).
Recall now that under (A), the terms |H|(s, t, x, y) and

|p ε ⊗ H (r)
ε | satisfy respectively and uniformly in ε the controls of equations (2.9), (2.10). The result then follows from the proof of Lemma 2 (see equation (3.11) for q = ∞ and (3.13) for q ∈ (d, +∞)) and the induction hypothesis.

Theorem 1 now simply follows from the controls of Lemma 3, the parametrix expansions (2.6) and (2.7) of the densities p, p ε and the asymptotics of the gamma function.

Stability for Markov Chains.

In this Section we prove Theorem 2. The strategy is rather similar to the one of Section 3.1 thanks to the series representation of the densities of the chains given in Proposition 2.

Recall first from Section 2.2 that we have the following representations for the density p h , p h ε of the Markov chains Y, Y (ε) in (1.3). For all 0 ≤ t i < t j ≤ T, (x, y) ∈ (R d ) 2 :

p h (t i , t j , x, y) = j-i r=0 ph ⊗ h H h,(r) (t i , t j , x, y), p h ε (t i , t j , x, y) = j-i r=0 ph ε ⊗ h H h,(r) ε (t i , t j , x, y).

Comparison of the frozen densities

The first key point for the analysis with Markov chains is the following Lemma.

Lemma 4 (Controls and Comparison of the densities and their derivatives).

There exist c, c 1 s.t. for all 0 ≤ t i < t j ≤ T, (x, y) ∈ (R d ) 2 and for all multi-index α, |α| ≤ 4:

|D α x ph (t i , t j , x, y)| + |D α x ph ε (t i , t j , x, y)| ≤ 1 (t j -t i ) |α|/2 ψ c,c1 (t j -t i , y -x), |D α x ph (t i , t j , x, y) -D α x ph ε (t i , t j , x, y)| ≤ ∆ ε,σ,γ (t j -t i ) |α|/2 ψ c,c1 (t j -t i , y -x), where -Under (I G ): ψ c,c1 (t j -t i , y -x) := c 1 p c (t j -t i , y -x),
-Under (I P,M ):

ψ c,c1 (t j -t i , y -x) := c 1 c d (t j -t i ) d/2 Q M-d-5 |y -x| (t j -t i ) 1/2 /c .
Proof. Note first that under (I G ) the statement has already been proved in Lemma 1. We thus assume that (I P,M ) holds. Introduce first the random vectors with zero mean:

Zy k,j := 1 (t j -t k ) 1/2 j-1 l=k σ(t l , y) √ hξ l+1 , Zy,(ε) k,j := 1 (t j -t k ) 1/2 j-1 l=k σ ε (t l , y) √ hξ l+1 .
Denoting by q j-k , q j-k,ε their respective densities, one has:

D α x ph (t k , t j , x, y) = 1 (t j -t k ) (d+|α|)/2 (-1) |α| D α z q j-k (z)| z= y-x (t j -t k ) 1/2 , D α x ph ε (t k , t j , x, y) = 1 (t j -t k ) (d+|α|)/2 (-1) |α| D α z q j-k,ε (z)| z= y-x (t j -t k ) 1/2 . (3.15)
From the Edgeworth expansion of Theorem 19.3 in Bhattacharya and Rao [START_REF] Bhattacharya | Normal approximations and asymptotic expansions[END_REF], for q j-k , q j-k,ε , one readily derives under (A), for |α| = 0 that there exists c 1 s.t. for all 0 ≤ t k < t j ≤ T, (x, y)

∈ (R d ) 2 , ph (t k , t j , x, y) + ph ε (t k , t j , x, y) ≤ c 1 (t j -t k ) d/2 1 1 + |x-y| (tj -t k ) 1/2 m , (3.16)
for all integer m < Md, where we recall that M stands for the initial decay of the density f ξ of the innovations bounded by Q M (see equation (1.7)). We can as well derive similarly to the proof of Theorem 19.3 in [START_REF] Bhattacharya | Normal approximations and asymptotic expansions[END_REF], see also Lemma 3.7 in [START_REF] Konakov | Local limit theorems for transition densities of Markov chains converging to diffusions[END_REF], that for all α, |α| ≤ 4:

|D α x ph (t k , t j , x, y)| + |D α x ph ε (t k , t j , x, y)| ≤ c 1 (t j -t k ) (d+|α|)/2 1 1 + |x-y| (tj -t k ) 1/2 m ,
(3.17) for all m < Md -4. Note indeed that differentiating in D α x the density and the terms of the Edgeworth expansion corresponds to a multiplication of the Fourier transforms involved by ζ α , ζ standing for the Fourier variable. Hence, from our smoothness assumptions in (I P,M ), after obvious modifications, the estimates of Theorem 9.11 and Lemma 14.3 from [START_REF] Bhattacharya | Normal approximations and asymptotic expansions[END_REF] apply for these derivatives. With these bounds, one then simply has to copy the proof of Theorem 19.3. Roughly speaking, taking derivatives deteriorates the concentration of the initial control in (3.16) up to the derivation order. On the other hand, the bound in (3.16) is itself deteriorated w.r.t. the initial concentration condition in (1.7). The key point is that the techniques of Theorem 19.3 in [START_REF] Bhattacharya | Normal approximations and asymptotic expansions[END_REF] actually provide concentration bounds for inhomogeneous sums of random variables with concentration as in (1.7) in terms of the moments of the innovations. To explain the bound in (3.16) let us observe that the m th moment of ξ is finite for m < Md.

Equations (3.16) and (3.17) give the first part of the lemma. Still from the proof of Theorem 19.3 in [START_REF] Bhattacharya | Normal approximations and asymptotic expansions[END_REF], one gets, under (A), that there exists C > 0 s.t. for all multi-indexes ᾱ, |ᾱ| ≤ 4, β, | β| ≤ m ≤ Md -5 for all j > k:

R d |ζ ᾱ| |D β ζ qj-k (ζ)| + |D β ζ qj-k,ε (ζ)| dζ ≤ C, (3.18)
where qj-k (ζ), qj-k,ε (ζ) stand for the respective characteristic functions of the random variables Zy k,j , Zy,(ε)

k,j at point ζ. To investigate the quantity |D α x ph (t k , t j , x, y) -D α x ph ε (t k , t j , x, y)| thanks to (3.15) define now for all α, |α| ≤ 4, β, |β| ≤ m ≤ M -d -5: ∀z ∈ R d , Θ j-k,ε (z) := z β D α z (q j-k (z) -q j-k,ε (z)) , ∀ζ ∈ R d , Θ j-k,ε (ζ) := (-i) |α|+|β| D β ζ (ζ α {q j-k (ζ) -qj-k,ε (ζ)}) . (3.19)
Let us now estimate the difference between the characteristic functions. From the Leibniz formula, we are led to investigate for all multi-indexes β, ᾱ, | β| ≤ |β|, |ᾱ| ≤ |α| quantities of the form:

(i β ) -1 ζ ᾱ(D β ζ qj-k (ζ) -D β ζ qj-k,ε (ζ)) = ζ ᾱE ( Zy k,j ) β exp[iζ • Zy k,j ] -( Zy,(ε) k,j ) β exp[iζ • Zy,(ε) k,j ] .
Assume first that j > k + 1. In that case, set now Zy k,j,1 := Zy k,⌈(j+k)/2⌉ , Zy k,j,2 := Zy k,j -Zy k,j,1 . Denoting similarly Zy,(ε) k,j,1 := Zy,(ε) k,⌈(j+k)/2⌉ , Zy,(ε) k,j,2 := Zy,(ε) k,j -Zy,(ε) k,j,1 for the perturbed process, we get:

(i β ) -1 ζ ᾱ(D β ζ qj-k (ζ) -D β ζ qj-k,ε (ζ)) = ζ ᾱ E ( Zy k,j,1 + Zy k,j,2 ) β exp[iζ • Zy k,j,1 ] exp[iζ • Zy k,j,2 ] -E ( Zy,(ε) k,j,1 + Zy,(ε) k,j,2 ) β exp[iζ • Zy,(ε) k,j,1 ] exp[iζ • Zy,(ε) k,j,2 ] = ζ ᾱ l,|l|≤| β| C l β E ( Zy k,j,1 ) l exp[iζ • Zy k,j,1 ] E ( Zy k,j,2 ) β-l exp[iζ • Zy k,j,2 ] - l,|l|≤| β| C l β E ( Zy,(ε) k,j,1 ) l exp[iζ • Zy,(ε) k,j,1 ] E ( Zy,(ε) k,j,2 ) β-l exp[iζ • Zy,(ε) k,j,2 ] = ζ ᾱ l,|l|≤| β| C l β E ( Zy k,j,1 ) l exp[iζ • Zy k,j,1 ] -E ( Zy,(ε) k,j,1 ) l exp[iζ • Zy,(ε) k,j,1 ] E ( Zy k,j,2 ) β-l exp[iζ • Zy k,j,2 ] + E ( Zy,(ε) k,j,1 ) l exp[iζ • Zy,(ε) k,j,1 ] E ( Zy k,j,2 ) β-l exp[iζ • Zy k,j,2 ] -E ( Zy,(ε) k,j,2 ) β-l exp[iζ • Zy,(ε) k,j,2
, where in the above expression we considered the binomial expansion for multi-indexes denoting by

C l β := β! ( β-l)!l!
with the corresponding definitions for factorials (see the proof of Lemma 1). Introduce now, for a multi-index l, |l| ∈ [[0, | β|]], the functions:

Ψ ᾱ, β-l 1 (ζ) := ζ ᾱE ( Zy k,j,2 ) β-l exp[iζ • Zy k,j,2 ] , Ψ ᾱ,l 2 (ζ) := ζ ᾱE ( Zy,(ε) k,j,1 ) l exp[iζ • Zy,(ε) k,j,1 ] ,
and

E 1,l (ζ) := E ( Zy k,j,1 ) l exp[iζ • Zy k,j,1 ] -E ( Zy,(ε) k,j,1 ) l exp[iζ • Zy,(ε) k,j,1 ] , E 2, β-l (ζ) := E ( Zy k,j,2 ) β-l exp[iζ • Zy k,j,2 ] -E ( Zy,(ε) k,j,2 ) β-l exp[iζ • Zy,(ε) k,j,2 ] .
Thus, we can rewrite from the previous computations:

(i β ) -1 ζ ᾱ(D β ζ qj-k (ζ) -D β ζ qj-k,ε (ζ)) = l,|l|≤| β| C l β (E 1,l Ψ ᾱ, β-l 1 )(ζ) + (E 2, β-l Ψ ᾱ,l 2 )(ζ) . (3.20)
Recall from (3.18) that we already have integrability for the contributions Ψ ᾱ, β-l 1 (ζ) and Ψ ᾱ,l 2 (ζ). Let us thus proceed with the control of E 1,l (ζ), E 2, β-l (ζ). We only give details for E 1,l (ζ), the contribution E 2, β-l can be handled similarly. We also consider |l| ≥ 2, since the cases |l| < 2 can be handled more directly. Write:

|E 1,l (ζ)| ≤ E[|( Zy k,j,1 ) l -( Zy,(ε) k,j,1 ) l |] + E[|( Zy,(ε) k,j,1 ) l || exp(iζ • Zy k,j,1 ) -exp(iζ • Zy,(ε) k,j,1 )|] ≤ C E[| Zy k,j,1 - Zy,(ε) k,j,1 |(| Zy k,j,1 | |l|-1 + | Zy,(ε) k,j,1 | |l|-1 )] + E[| Zy,(ε) k,j,1 | |l| |ζ|| Zy k,j,1 - Zy,(ε) k,j,1 |] .
Apply now Hölder's inequality with p 1 = |l|, q 1 = |l|/(|l|-1) for the first term and p 2 = (|l|+1)/|l|, q 2 = |l|+1 for the second one so that all the contribution appear with the same power (in order to equilibrate the constraints concerning the intregrability conditions). One gets:

|E 1,l (ζ)| ≤ C E[| Zy k,j,1 - Zy,(ε) k,j,1 | |l| ] 1/|l| {E[| Zy k,j,1 | |l| ] (|l|-1)/|l| + E[| Zy,(ε) k,j,1 | |l| ] (|l|-1)/|l| } + |ζ|E[| Zy,(ε) k,j,1 | |l|+1 ] |l|/(|l|+1) E[| Zy k,j,1 - Zy,(ε) k,j,1 | |l|+1 ] 1/(|l|+1) . (3.21)
The point is now to prove, since we have assumed

m ≤ M -d -5 ⇐⇒ m + 1 ≤ M -d -4, that there exists c s.t. for all r ≤ m + 1, E[| Zy k,j,1 - Zy,(ε) k,j,1 | r ] 1/r ≤ c∆ ε,σ,γ , E[| Zy k,j,1 | r ] 1/r + E[| Zy,(ε) k,j,1 | r ] 1/r ≤ c. (3.22)
Let us establish the point for the difference, the other bounds can be derived similarly. Define for all i ∈

[[k, j]], Mi := √ h i-1 r=k (σ -σ ε )(t r , y)ξ r+1 . The process ( Mi ) i∈[[k,j]] is a square integrable martingale (in discrete time, w.r.t. F i := Σ(ξ r , r ≤ i), Σ-

field generated by the innovation up to the current time). Its quadratic variation writes [ M ]

i = h i-1 r=k |(σ -σ ε )(t r , y)| 2 |ξ r+1 | 2
and the Burkholder-Davies-Gundy inequalities, see e.g. Shiryaev [START_REF] Shiryaev | Probability, Second Edition[END_REF], give for all r ≤ Md -4:

E[ sup i∈[[k,j]] | Mi | r ] ≤ c r E[[ M ] r/2 j ] = c r h r/2 E[( j-1 i=k |(σ -σ ε )(t i , y)| 2 |ξ i+1 | 2 ) r/2 ].
(3.23)

If r = 2 one readily gets:

E[| Zy k,j,1 - Zy,(ε) k,j,1 | 2 ] ≤ c 2 (t j -t k ) E[ sup i∈[[k,j]] | Mi | 2 ] ≤ c 2 h (t j -t k ) ∆ 2 ε,σ,γ j-1 i=k E[|ξ i+1 | 2 ] ≤ c2 ∆ 2 ε,σ,γ .
Let us thus assume r > 2 and derive from (3.23)

E[| Zy k,j,1 - Zy,(ε) k,j,1 | r ] ≤ c r (t j -t k ) r/2 E[ sup i∈[[k,j]] | Mi | r ] ≤ c r h r/2 (t j -t k ) r/2 E[( j-1 i=k |(σ -σ ε )(t i , y)| r |ξ i+1 | r )( j-1 i=k 1) r/2(1-2/r) ],
applying Hölder's inequality for the counting measure with p = r/2, q = r/(r -2) for the last inequality. This finally gives:

E[| Zy k,j,1 - Zy,(ε) k,j,1 | r ] ≤ c r h r/2 (t j -t k ) r/2 (j -k) r/2-1 ∆ r ε,σ,γ j-1 i=k E[|ξ i+1 | r ] ≤ cr ∆ r ε,σ,γ .
Since we have assumed r ≤ m + 1 ≤ Md -4, this gives the first control in (3.22). The other one readily follows replacing σσ ε by σ or σ ε .

From equations (3.21), (3.22) and similar controls for E 2, β-l (ζ) we finally derive:

|E 1,l (ζ)| + |E 2, β-l (ζ)| ≤ C 1 ∆ ε,σ,γ (1 + |ζ|).
As a result we have from (3.19) and (3.20):

|D β ζ (ζ α (q j-k (ζ) -D β ζ qj-k,ε (ζ)))| ≤ C∆ ε,σ,γ β, | β| ≤ |β| ᾱ = α -(β -β). l,|l|≤| β| (|Ψ ᾱ, β-l 1 (ζ)| + |Ψ ᾱ,l 2 (ζ)|)(1 + |ζ|) .
We finally derive from (3.19) and (3.18) (which thanks to the smoothness assumption on Q M in (I P,M ) holds as well for a multi-index ᾱ, |ᾱ| = 5):

|Θ j-k,ε (z)| ≤ 1 (2π) d R d | Θj-k,ε (ζ)|dζ ≤ c∆ ε,σ,γ .
(3.24)

From (3.15) this concludes the proof for j > k + 1. If j = k + 1 the previous arguments can be simplified and lead to the same results.

Comparison of the parametrix kernels

This step is crucial and actually the key to the result for the Markov chains. We focus for simplicity on the case q = +∞, for which pointwise controls for the differences between the drift coefficients are available, and which already emphasizes all the difficulties. The case q ∈ (d, +∞) for the drifts could be handled as in Lemma 2, using similar Hölder inequalities.

We actually have the following Lemma.

Lemma 5 (Control of the One-Step Convolution for the Chain.). There exists c 1 , c s.t. for all 0 ≤ t k < t j ≤ T, (z, y) ∈ (R d ) 2 :

|(H h -H h ε )(t k , t j , z, y)| ≤ ∆ ε,γ,∞ (t j -t k ) 1-γ/2 Φ c,c1 (t j -t k , z -y), with -Φ c,c1 (t j -t k , z -y) = ψ c,c1 (t j -t k , z -y) under (I G ). -Φ c,c1 (t j -t k , z -y) = ψ c,c1 (t j -t k , z -y) 1 + |z-y| (tj -t k ) 1/2 γ , under (I P,M ),
where ψ c,c1 is defined according to the assumptions on the innovations in Lemma 4.

Proof. The case k = j + 1 involves directly differences of densities and could be treated more directly than the case k > j + 1. We thus focus on the latter. Introduce for k ∈ [[0, N ]], (x, w) ∈ (R d ) 2 the one step transitions:

T h (t k , x, w) := b(t k , x)h + h 1/2 σ(t k , x)w, T h ε (t k , x, w) := b ε (t k , x)h + h 1/2 σ ε (t k , x)w, T h 0 (t k , x, w) := h 1/2 σ(t k , x)w, T h ε,0 (t k , x, w) := h 1/2 σ ε (t k , x)w.
(3.25)

From the definition of H h , H h ε , recalling that f ξ stands for the density of the innovation, the difference of the kernels writes:

(H h -H h ε )(t k , t j , z, y) = h -1 R d
dwf ξ (w) ph (t k+1 , t j , z + T h (t k , z, w), y)ph (t k+1 , t j , z + T h 0 (t k , y, w), y)

ph ε (t k+1 , t j , z + T h ε (t k , z, w), y)ph ε (t k+1 , t j , z + T h 0,ε (t k , y, w), y) .

(3.26)

Let us now perform a Taylor expansion at order 2 with integral rest. To this end, let us first introduce for λ ∈ [0, 1] the mappings:

ϕ h λ : R d × R d -→ R (T 1 , T 2 ) -→ Tr D 2 z ph (t k+1 , t j , z + λT 1 , y)[T 2 T * 2 ] , ϕ h λ,ε : R d × R d -→ R (T 1 , T 2 ) -→ Tr D 2 z ph ε (t k+1 , t j , z + λT 1 , y)[T 2 T * 2 ] , (3.27)
where T 2 is viewed as a column vector and T * 2 denotes its transpose. Recalling as well that ξ is centered we get:

∆H h,ε (t k , t j , z, y) := (H h -H h ε )(t k , t j , z, y) = D z ph (t k+1 , t j , z, y), b(t k , z) -D z ph ε (t k+1 , t j , z, y), b ε (t k , z) +h -1 R d dwf ξ (w) 1 0 dλ(1 -λ) × ϕ h λ (T h (t k , z, w), T h (t k , z, w)) -ϕ h λ (T h 0 (t k , y, w), T h 0 (t k , y, w)) -ϕ h λ,ε (T h ε (t k , z, w), T h ε (t k , z, w)) -ϕ h λ,ε (T h 0,ε (t k , y, w), T h 0,ε (t k , y, w)) =: (∆ 1 H h,ε + ∆ 2 H h,ε )(t k , t j , z, y), (3.28)
where for i ∈ {1, 2}, ∆ i H h,ε is associated with the terms of order i. The idea is now to make ∆ ε,γ,∞ appear explicitly. The term ∆ 1 H h,ε is the easiest to handle. We can indeed readily write:

∆ 1 H h,ε (t k , t j , z, y) = D z ph (t k+1 , t j , z, y), [b(t k , z) -b ε (t k , z)] -(D z ph ε -D z ph )(t k+1 , t j , z, y), b ε (t k , z) .
From Assumption (A3), equation (1.5) and Lemma 4 we derive for q = +∞:

|∆ 1 H h,ε (t k , t j , z, y)| ≤ C∆ ε,γ,∞ (t j -t k ) 1/2 ψ c,c1 (t j -t k , y -z).
(3.29)

The term ∆ 2 H h,ε is trickier to handle. Define to this end:

∆ϕ h,ε λ (t k , z, y, w) := ϕ h λ (T h (t k , z, w), T h (t k , z, w)) -ϕ h λ (T h 0 (t k , y, w), T h 0 (t k , y, w)) -ϕ h λ,ε (T h ε (t k , z, w), T h ε (t k , z, w)) -ϕ h λ,ε (T h 0,ε (t k , y, w), T h 0,ε (t k , y, w)) .
Let us then decompose:

∆ϕ h,ε λ (t k , z, y, w) := ϕ h λ (T h (t k , z, w), T h (t k , z, w)) -ϕ h λ (T h (t k , z, w), T h 0 (t k , y, w)) -ϕ h λ,ε (T h ε (t k , z, w), T h ε (t k , z, w)) -ϕ h λ,ε (T h ε (t k , z, w), T h 0,ε (t k , y, w)) + ϕ h λ (T h (t k , z, w), T h 0 (t k , y, w)) -ϕ h λ (T h 0 (t k , y, w), T h 0 (t k , y, w)) -ϕ h λ,ε (T h ε (t k , y, w), T h 0,ε (t k , y, w)) -ϕ h λ,ε (T h 0,ε (t k , z, w), T h 0,ε (t k , y, w)) =: (∆ 1 ϕ h,ε λ + ∆ 2 ϕ h,ε λ )(t k , z, y, w), (3.30)
and write from (3.28):

∆ 2 H h,ε (t k , t j , z, y) = h -1 R d dwf ξ (w) 1 0 dλ(1 -λ)(∆ 1 ϕ h,ε λ + ∆ 2 ϕ h,ε λ )(t k , z, y, w) =: (∆ 21 H h,ε + ∆ 22 H h,ε )(t k , t j , z, y), (3.31)
for the associated contributions in ∆ 2 H h,ε . Again, we have to consider these two terms separately.

Term ∆ 21 H h,ε . We first write from (3.30):

∆ 1 ϕ h,ε λ (t k , z, y, w) = ϕ h λ (T h (t k , z, w), T h (t k , z, w)) -ϕ h λ (T h (t k , z, w), T h 0 (t k , y, w)) - ϕ h λ (T h (t k , z, w), T h ε (t k , z, w)) -ϕ h λ (T h (t k , z, w), T h 0,ε (t k , y, w)) + ϕ h λ (T h (t k , z, w), T h ε (t k , z, w)) -ϕ h λ (T h (t k , z, w), T h 0,ε (t k , y, w)) -ϕ h λ (T h ε (t k , z, w), T h ε (t k , z, w)) -ϕ h λ (T h ε (t k , z, w), T h 0,ε (t k , y, w)) -ϕ h λ,ε (T h ε (t k , z, w), T h ε (t k , z, w)) -ϕ h λ,ε (T h ε (t k , z, w), T h 0,ε (t k , y, w)) -ϕ h λ (T h ε (t k , z, w), T h ε (t k , z, w)) -ϕ h λ (T h ε (t k , z, w), T h 0,ε (t k , y, w)) =: 3 i=1 ∆ 1i ϕ h,ε λ (t k , z, y, w). (3.32)
We now state some useful controls for the analysis. Namely, setting:

D(t k , z, y, w) := T h (t k , z, w)T h (t k , z, w) * -T h 0 (t k , y, w)T h 0 (t k , y, w) * , D ε (t k , z, y, w) := T h ε (t k , z, w)T h ε (t k , z, w) * -T h 0,ε (t k , y, w)T h 0,ε (t k , y, w) * ,
we have from (A3) and equation (1.5) for q = +∞ :

(|D| + |D ε |)(t k , z, y, w) ≤ c(h 2 + h 3/2 |w| + h(1 ∧ |z -y|) γ |w| 2 ), |D -D ε |(t k , z, y, w) ≤ c∆ ε,γ,∞ (h 2 + h 3/2 |w| + h(1 ∧ |z -y|) γ |w| 2 ). (3.33)
From the definition of ϕ h λ in (3.27), equation (3.32), the control (3.33) and Lemma 4, we get:

|∆ 11 ϕ h,ε λ |(t k , z, y, w) ≤ c∆ ε,γ,∞ ψ c,c1 (t j -t k , y -(z + λT h (t k , z, w))) (t j -t k ) (h 2 + h 3/2 |w| + h(1 ∧ |z -y|) γ |w| 2 ). (3.34)
We would similarly get from Lemma 4 and (3.33):

|∆ 13 ϕ h,ε λ |(t k , z, y, w) ≤ c∆ ε,γ,∞ ψ c,c1 (t j -t k , y -(z + λT h ε (t k , z, w))) (t j -t k ) (h 2 + h 3/2 |w| + h(1 ∧ |z -y|) γ |w| 2 ), |∆ 12 ϕ h,ε λ |(t k , z, y, w) ≤ ψ c,c1 (t j -t k , y -(z + θλT h (t k , z, w) + (1 -θ)λT h ε (t k , z, w))) (t j -t k ) 3/2 ×|(T h -T h ε )(t k , z, w)||D ε |(t k , z, y, w) ≤ c∆ ε,γ,∞ ψ c,c1 (t j -t k , y -(z + θλT h (t k , z, w) + (1 -θ)λT h ε (t k , z, w))) (t j -t k ) 3/2 ×(h 2 + h 3/2 |w| + h(1 ∧ |z -y|) γ |w| 2 )(h + h 1/2 |w|), (3.35)
for some θ ∈ (0, 1), using as well (3.25) and (1.5) for the last inequality. The point is now to get rid of the transitions appearing in the function ψ c,c1 . We separate here the two assumptions at hand.

-Under (I G ), it suffices to remark that by the convexity inequality |zy -

Θ| 2 ≥ 1 2 |z -y| 2 -|Θ| 2 , for all Θ ∈ R d : ψ c,c1 (t j -t k , y -z -Θ) ≤ c 1 c d/2 (2π(t j -t k )) d/2 exp - c 4 |z -y| 2 t j -t k exp c 2 |Θ| 2 t j -t k .
Now, if Θ is one of the above transitions or linear combination of transitions, we get from (3.25):

ψ c,c1 (t j -t k , y -z -Θ) ≤ c 1 (c/2) d/2 (2π(t j -t k )) d/2 exp - c 4 |z -y| 2 t j -t k exp( c 2 K 2 2 |w| 2 ), (3.36)
up to a modification of c 1 observing that h/(t jt k ) ≤ 1 and with K 2 as in (A 1 ). Since c can be chosen small enough in the previous controls, up to deteriorating the concentration properties in Lemma 4, the last term can be integrated by the standard Gaussian density f ξ appearing in (3.31). We thus derive, from (3.36), (3.34), (3.35) and the definition in (3.32), up to modifications of c, c 1 :

|∆ 1 ϕ h,ε λ |(t k , z, y, w) ≤ ∆ ε,γ,∞ hcψ c,c1 (t j -t k , z -y) exp(c|w| 2 ) 1 + |w| (t j -t k ) 1/2 + |z -y| γ |w| 2 t j -t k ,
which plugged into (3.31) yields up to modifications of c, c, c 1 :

|∆ 21 H h,ε (t k , t j , z, y)| ≤ c ∆ ε,γ,∞ (1 ∨ T (1-γ)/2 )ψ c,c1 (t j -t k , z -y) (t j -t k ) 1-γ/2 .
(3.37)

-Under (I P,M ), we only detail the computations for the off diagonal regime |z -y| ≥ c(t jt k ) 1/2 which is the most delicate to handle. In this case, we have to discuss according to the position of w w.r.

| + |∆ 13 ϕ h,ε λ |)(t k , z, y, w) ≤ c∆ ε,,γ,∞ ψ c,c1 (t j -t k , y -z) (t j -t k ) (h 2 + h 3/2 |w| + h(1 ∧ |z -y|) γ |w| 2 ), |∆ 12 ϕ h,ε λ |(t k , z, y, w) ≤ c∆ ε,γ,∞ ψ c,c1 (t j -t k , y -z) (t j -t k ) 3/2 (h 2 + h 3/2 |w| + h(1 ∧ |z -y|) γ |w| 2 )(h + h 1/2 |w|).
On the other hand, when w ∈ D we use f ξ to make the off-diagonal bound of ψ c,c1 (t jt k , yz) appear. Namely, we can write: We emphasize that in the case of innovations with polynomial decays, the control on the difference of the kernels again induces a loss of concentration of order γ in order to equilibrate the time singularity.

f ξ (w) ≤ c 1 (1 + |w|) M ≤ c 1 (1 + |z-y| h 1/2 ) M-(d+4) 1 (1 + |w|) d+4 ≤ c 1 (1 + |z-y| (tj -t k ) 1/2 ) M-(d+4)
Term ∆ 22 H h,ε . This term can be handled with the same arguments as ∆ 21 H h,ε . For the sake of completeness we anyhow specify how the different contributions appear. Namely, with the notations of (3.30) and (3.31): ∆ 2 ϕ h,ε λ (t k , z, y, w) = 1 0 dµ D T1 ϕ h λ (T h 0 (t k , y, w) + µ(T h (t k , z, w) -T h 0 (t k , y, w)), T h 0 (t k , y, w)), T h (t k , z, w) -T h 0 (t k , y, w)

-D T1 ϕ h λ,ε (T h 0,ε (t k , y, w) + µ(T h ε (t k , z, w) -T h 0,ε (t k , y, w)), T h 0,ε (t k , y, w)), T h ε (t k , z, w) -T h 0,ε (t k , y, w) = 1 0 dµ D T1 ϕ h λ (T h 0 (t k , y, w) + µ(T h (t k , z, w) -T h 0 (t k , y, w)), T h 0 (t k , y, w)), (T h (t k , z, w) -T h 0 (t k , y, w)) -(T h ε (t k , z, w) -T h 0,ε (t k , y, w)) -1 0 dµ D T1 ϕ h λ,ε (T h 0,ε (t k , y, w) + µ(T h ε (t k , z, w) -T h 0,ε (t k , y, w)), T h 0,ε (t k , y, w))

-D T1 ϕ h λ (T h 0 (t k , y, w) + µ(T h (t k , z, w) -T h 0 (t k , y, w)), T h 0 (t k , y, w)) , T h ε (t k , z, w) -T h 0,ε (t k , y, w) =: (∆ 21 ϕ h,ε λ + ∆ 22 ϕ h,ε λ )(t k , z, y, w).

In ∆ 21 ϕ h,ε λ we have sensitivities of order 3 for the density, giving time singularities in (t jt k ) -3/2 , which are again equilibrated by the the multiplicative factor: |T h 0 (t k , y, w)[T h 0 (t k , y, w)] * | ×|(T h (t k , z, w) -T h 0 (t k , y, w)) -(T h ε (t k , z, w) -T h 0,ε (t k , y, w))| ≤ c(h 2 + h 3/2 |w| + h|w| 2 )∆ ε,γ,∞ (h + h 1/2 (1 ∧ |z -y|) γ |w|), where the last inequality is obtained similarly to (3.33) using as well (1.5). The same kind of controls can be established for ∆ 22 ϕ h,ε λ . Anyhow, the analysis of this term leads to investigate the difference of third order derivatives, which finally yields contributions involving derivatives of order four. This is what induces the final concentration loss under (I P,M ), i.e. we need to integrate a term in |w| 4 (see also equation (3.38) in which we performed the splitting of f ξ on the off-diagonal region to integrate a contribution in |w| 3 ).

We can thus claim that |∆H h,ε 22 (t k , t j , z, y)| ≤ ∆ ε,γ,∞ Φ c,c1 (t jt k , zy) (t jt k ) 1-γ/2 .

Plugging the above control and (3.39) (or (3.37) under (I G )) into (3.31) we derive:

|∆H h,ε 2 (t k , t j , z, y)| ≤ ∆ ε,γ,∞ Φ c,c1 (t j -t k , z -y) (t j -t k ) 1-γ/2
, which together with (3.29) and the decomposition (3.28) completes the proof.

From Lemmas 4 and 5 the proof of Theorem 2 is achieved, under (I G ), following the steps of Lemmas 2 and 3, using the Hölder inequalities for the differences of the drift terms for q ∈ (d, +∞).

The point is that we want to justify the following inequality under (I P,M ) and q = +∞:

|(p h ⊗ h H h,(r)ph ε ⊗ h H h,(r) ε )(t i , t j , x, y)| (3.40)

≤ (r + 1)∆ ε,γ,∞ {(1 ∨ T (1-γ)/2 )c 1 } r+1 Γ( γ 2 ) r Γ(1 + r γ 2 ) c d (t j -t i ) d/2 Q M-(d+5+γ)
yx (t jt i ) 1/2 /c (t jt i ) rγ 2 .

The only delicate point, w.r.t. the analysis performed for diffusions, consists in controlling the convolutions of the densities with polynomial decay. To this end, we can adapt a technique used by Kolokoltsov [Kol00] to investigate convolutions of "stable like" densities. Set m := M -(d + 5 + γ) and denote for all 0 ≤ i < j ≤ N, x ∈ R d by q m (t jt i , x) :=

c d (tj -ti) d/2 Q M-(d+5+γ)
x (tj -ti) 1/2 /c the density with polynomial decay appearing in Lemmas 4 and 5. Let us consider for fixed i < k < j, (x, y) ∈ (R d ) 2 the convolution:

I 1 t k (t i , t j , x, y) := R d
dzq m (t kt i , zx)q m (t jt k , yz).

(3.41)

-If |x -y| ≤ c(t jt i ) 1/2 (diagonal regime for the parabolic scaling), it is easily seen that one of the two densities in the integral (3.41) is homogeneous to q m (t jt i , yx). Namely, if (t kt i ) ≥ (t jt i )/2, q m (t kt i , zx) ≤ c d cm (t k -ti) d/2 ≤ 2 d/2 c d cm (tj -ti) d/2 ≤ cq m (t jt i , yx). Thus, I 1 t k (t i , t j , x, y) ≤ cq m (t jt i , yx)

R d
dzq m (t jt k , yz) = cq m (t jt i , yx).

If (t kt i ) < (t jt i )/2, the same operation can be performed taking q m (t jt k , yz) out of the integral, observing again that in that case q m (t jt k , yz) ≤ cq m (t jt i , yx).

  t. yz. With the notations of (A2), introduce D := { w ∈ R d : {Λh} 1/2 | w| ≤ |z -y|/2}. If w ∈ D, then, still from (3.34), (3.35),

	(|∆ 11 ϕ h,ε λ

  1 (1 + |w|) d+4 , (3.38)where the last splitting is performed in order to integrate the contribution in |w| 3 coming from the upper bound for |∆ 12 ϕ h,ε λ | in (3.35). Plugging the above controls in (3.31) yields:|∆H h,ε 21 (t k , t j , z, y)| ≤ ∆ ε,γ,∞ Φ c,c1 (t jt k , zy) (t jt k ) 1-γ/2 .(3.39)

Acknowledgements

The authors would like to thank the two anonymous referees for their careful reading and suggestions that truly improved the previous version of this work.

* The article was prepared within the framework of a subsidy granted to the HSE by the Government of the Russian Federation for the implementation of the Global Competitiveness Program. * * Support by Deutsche Forschungsgemeinschaft through the Research Training Group RTG 1953 is gratefully acknowledged by the second author.

-If |x -y| > c(t jt i ) 1/2 (off-diagonal regime), we introduce

recalling that, under (I P,M ), m > d for the last but one inequality. The same operation could be performed on A 1 .

We have thus established that, there exist c > 1 s.t. for all 0 ≤ i < k < j, (x, y) ∈ (R d ) 2 :

From the controls of Lemma 5 and following the strategy of Lemma 3, we will be led to consider convolutions of the previous type involving Γ functions. The above strategy thus yields (3.40) by induction.


Deqing Zou
email: deqingzou@hust.edu.cn

Hao Qin

Hai Jin

Weizhong Qiang

Zongfen Han

Xueguang Chen

Improving Log-based Fault Diagnosis by Log Classification

In modern computer systems, system event logs have always been the primary source for checking the system status. As computer systems become more complex, such as cloud computing systems, the interaction among software and hardware is increasingly frequently. These components will generate enormous log information, including running reports and fault information. The massive data is a great challenge for analysis with manual method. In this paper, we implement a log management and analysis system, which can assist system administrators to understand the real-time status of the entire system, classify logs into different fault types, and determine the root cause of the faults. In addition, we improve the existing fault correlation analysis method based on the results of system log classification. We apply the log management and analysis system to cloud computing environment for evaluation. The results show that our system can classify fault logs effectively and automatically. By using the proposed system, administrators can easily detect the root cause of faults.

Introduction

With the widespread usage of cloud computing, computer systems are becoming increasingly co mplex and the components within the entire system also become diverse. Once some key parts failed, the whole system would be seriously imp licated due to the frequent interactions and high coupling. Therefo re, an effective fault detection and analysis method can help system administrators to locate the fault and identify the cause, which plays an essential role in large systems management. System and software logs are important sources for diagnosing the system and software faults. Ho wever, for large systems , various components will generate amounts of log informat ion in real t ime. If a fault occurs, it is difficult to extract useful informat ion fro m the system efficiently and locate the fault accurately. Typically, we have to manually extract the useful info rmation fro m vast amounts of data, which would seriously delay the response time of fault recovery. Therefore, a n unified management system for fault log analysis is required, which can automatically identify the fault type and analyze the cause of the fau lts. It will provide a great help for system management.

A lot of studies have been proposed on log-based fault analysis, main ly falling in the following directions . The first is log collection and analysis, which investigate how to effectively and efficiently gather log info rmation [START_REF] Zawoad | SecLaaS: secure logging-as-a-service for cloud forensics[END_REF]. This information will be used to get profile in formation fo r analy zing the system situation [START_REF] Rao | Identifying faults in large-scale distributed systems by filtering noisy error logs[END_REF] and e xtract the feature [START_REF] Yuan | SherLog: error diagnosis by connecting clues from run-time logs[END_REF]. The second part is fault location, which aims at determining the control flo w of software through logs [START_REF] Fu | Execution anomaly detection in distributed systems through unstructured log analysis[END_REF] and uses the source code [START_REF] Xu | Detecting large-scale system problems by mining console logs[END_REF] to locate the position of the faults occurred. The researches mentioned above use the log fault analysis in different scenarios. But few of them focus on the system ad ministrator's perspective to design an integrated fault analysis system. The third part is fault correlation analysis. In the mult i-node environ ment, the fault propagation is a critical problem for fault diagnosis. The area focuses on using log informat ion to determine the connection between different fau lts. Researches use time and spatial correlation to find some connection, few of them consider the meaning of lo gs.

An integrated fault log analysis system will assist admin istrators to perform fault analysis, improving the ad ministrator's ability to respond to the system fault and reduce the time consuming of fault processing. We implement a new fault classification method to assist manager to understand computer system and use fault correlation analysis to locate the root-cause of fault. In this paper, we propose an integrated fault log analysis platform (UiLog) to collect and manage various components logs, storing, filtering, and analy zing logs for administrators to quickly locate fault and analyze the cause of faults.

This platform consists of three components: 1) Fault log collection module is mainly used to collect log data fro m various components. 2) Fault log analysis module classifies log into the identifies fault type in real time. 3) Fault log correlation analysis module collects fault log caused by same root -fault as a tuple and tries to find the root cause of such faults. UiLog have deployed in a practical cloud environment, help ing admin istrator to troubleshoot and find the root cause of fault. Our main contributions are:

We propose a novel classification method for fault logs, using fault keyword matrix to imp rove the accuracy. It reduces the time of d etermin ing the fault type and the workload of manual processing. Moreover, this method can be more convenient to add a new fault type without recalculated.

We improve the existing log correlat ion analysis. It co mbines the result of fault classification and time windows correlat ion analysis. Our method uses the fault type of logs as one factor in determining the size of the time window. It imp roves the accuracy of the log correlation and the location of the fault's root cause.

We illustrate a comprehensive log management system, which can help administrators to quickly grasp the operation status of the system and save troubleshooting time.

The remain ing of this paper is organized as follo ws. Section 2 discusses background and related work. Sect ion 3 outlines the structure of UiLog system and describes the implementation of the system. Sect ion 4 describes the evaluation of our system whereas section 5 presents conclusions and future work.

Related Work

The aim of log information is to extract useful information fro m fault logs. The important techniques of previous work are mainly relied on regular expressions [START_REF] Xu | Detecting large-scale system problems by mining console logs[END_REF]. However, the ru les of regular exp ressions require different knowledge fro m laborious and expert [START_REF] James | Listening to your cluster with LoGS[END_REF]. In addition, the deploy ment of new application and upgrading of system will change these templates of log frequency. It is difficult for designing these regular expressions [START_REF] Jain | Extracting the textual and temporal structure of supercomputing logs[END_REF].

Various studies are looking for how to understand the mean of logs, but it may be useful for detecting faults in cloud computing system [START_REF] Stearley | Bad words: Finding faults in Spirit's syslogs[END_REF]. Many interesting features are displayed by these studies. For instance, Stearley [START_REF] Sandia | Towards informatic analysis of syslogs[END_REF] has a new discovery that only through words cannot detect the fault type fro m logs. The position of each word is a powerful indicator to distinguish different messages.

Researchers have also looked at other way except system logs to diagnose system, such as application console logs [START_REF] Xu | M ining Console Logs for Large-Scale System Problem Detection[END_REF]. However, this technique is limited to application specific anomalies and requires source code [START_REF] Salfner | Error Log Processing for Accurate Failure Prediction[END_REF].

In addition, several techniques and algorithms for auto matic log classific ation have been developed. [START_REF] Park | Proactive self-healing system based on multi-agent technologies[END_REF] attempts to classify different raw logs into a set of categories. Moreover, in [START_REF] Li | An integrated framework on mining logs files for computing system management[END_REF], the authors try to use the modified naive Bayesian model and Hidden Markov Models (HMM) to classify event logs based on the IBM CBE (Co mmon Base Event) format. On the other hand, SLCT [START_REF] Vaarandi | A data clustering algorithm for mining patterns from event logs[END_REF] and Loghound [START_REF] Vaarandi | A breadth-first algorithm for mining frequent patterns from event logs[END_REF] are designed specifically to discover the format of logs and classify ro w log automatically. They use two similar algorith ms, which are useful to extract the template from logs.

Similar to correlation analysis, time and spatial correlation techniques have been applied to a variety of large scale co mputing systems [START_REF] Oliner | What supercomputers say: A study of five system logs[END_REF]. The current trend of this study is to use tuple with a fixed value for time window, such as 5 minutes [START_REF] Pecchia | Improving log-based field failure data analysis of multi-node computing systems[END_REF].

Content-based correlation is also a hot topic. For examp le, [START_REF] Pecchia | Improving log-based field failure data analysis of multi-node computing systems[END_REF] applies the lift data by min ing operator to find frequent event patterns starting fro m log contents and try to isolate accidental patterns.

3

UiLog System

Overview of UiLog

The Unify Log Analysis System (Ui Log) is a fau lt analysis and diagnosis system, which collects the system log informat ion of each component and track logs for statistics.

Through the fault classification, UiLog learns the classification rules fro m training set of art ificial classification. After that, the system determines the fault type in real time according the rule lib rary. In addition, the fault correlation analysis can be deployed when system ad ministrators need to diagnosis fault. Considering the propagation of the fault, UiLog can mining the association between faults generated by the same root-cause and collect these fault logs into the same cluster.

To implement the process, we apply three modules to represent the UiLog: fault log collection, fault log analysis, and fault log correlation analysis. As shown in Fig. 1, the fault log co llect ion module co llects logs fro m the entire target node. It collects software and system logs in all co mponents. The analysis node is responsible for fault analysis and diagnose. In the analysis node, the fault log co llect ion module stores all logs in log information database for analysing.

The process of fault log analysis module is shown in Fig. 1. The fault log collection module will be deployed respectively into Target Node for log collection and Analysis Node for storage. This module will gather software log and system log fro m Target Node and store log informat ion into Log Information database. In the Analysis Node, the Fault log analysis module will ext ract log structure from fau lt log collection module. At first, these structures will be classified into different fault type by administrator. The Fault log analysis module will learn classification rules fro m administrator and store it into fault template database. After that, the fault log analysis module will automatic classify log for log diagnose.

The correlation analysis module will use the result of fault classification and expert knowledge to provide correlation analysis report from artificial analysis.

Fault Log Analysis

The fault log analysis module aims at classifying log int o different log types. We adopt an examp le to illustrate the process: when an ad min istrator analyzes the logs (Table 1) the first step is to determine whether the logs are fault record or not. The most commonly used method is finding keyword. The second st ep is the log analysis.

The system ad ministrator will weed out the details of log informat ion to determine the cause of the fault. The first row in Table 1 shows that "Read socket failed", which indicates a socket problem when reading. The third step is th e log classification. It is easily to determine that the first log belongs to the network fault and the subcategory is the remote network connection fault.

Through the above analysis, we can conclude that if we can deal with the semantic analysis of keywords and determine the fault type in advance, a large load of work can be automatically processed by the log analysis and classification.

As shown in Fig. 2, UiLog log classification method is divided into five steps. 1) Log Pretreat ment. 2) Extracting Invariants. 3) Filtering Template Information. 4) Obtaining Fau lt Keyword Matrix. 5) Classifying Log Information. These steps will gradually ext ract log information fo r fault classification and remove irrelevant content. The input of the whole algorith m is the log data flow obtained by fault log collection module. The output has two parts. One is the log classification ru les (output at Step 4) based on the training set. The rules in our algorithm appear as the fault keyword matrix. The other is fau lt category of every log (output at Step 5). The following sub-sections describe each step of the algorithm in more detail.

1) Log Preprocessing

Log preprocessing contains two parts. The first part is to filter repeated logs generated by the system. Many software faults are insufficient to cause the collapse in co mputer system and the co mponent will persistent send fault message. We can choose an appropriate threshold to filter the duplicate logs.

The second part is the log filter of mean ingless words. We use English Stop Word Table to filter meaningless words. In accessory, cons idering the special of logic, we also filter out the most of adjectives and adverbs.

2) Extracting Invariants

After preprocessing, the next step is to extract the template information. The template of log is used for classification. This step can reduce the s olution space and compress the size of fault keyword mat rix. In addition, the te mplate will be used to match logs for fast classification.

3) Filtering Template Information

Before classification, we use automatic classification approach DBSCAN (Density-Based Spatial Clustering of Applications with Noise) [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] to classify the log for further reducing the manually determining space.

For a classificat ion algorith m, based on Levenshten Distance (LD), we define our Log Levenshten Distance (LLD) to measure the distance between a log and a clu ster. The distance between log A and log B is described in [START_REF] Zawoad | SecLaaS: secure logging-as-a-service for cloud forensics[END_REF].

        2, , LD A B LLD A B length A length B    (1)
wherein LD(A, B) is the original Levenshten distance, length() indicates the length of the log.

4) Obtaining Fault Keyword Matrix

When the system is in the learning stage, this step will learn the result of artificial classification. While in the running stage, this step will modify the classification's rules. Through the step of filter template information , the remain ing nu mbers of raw logs have less than original.

For learning stage, the first requires administrator to classify logs manually. We pre-define a fault catalogue in the cloud environment. Through the last step, the classification has already marked a label to each fault cluster. Now the ad min istrator will first modify the auto matically label, then adjust the result of the class ification. As shown in Table 2, after indicated the category, the result will be stored in the way of "Content: Mark". Label is a number that represents the fault type. Next, Ui Log needs to learn the results of the artificial classification. Here we propose the Fa ult Keyword Matrix (FK-Matrix) for saving the learning result of artificial classification.

The FK-Matrix (mat rix A) is a two-dimensional matrix constructed by the probability of each word appeared in each fault type in the template. It is an m×n matrix. M represents the different number of words in all the sample log, while n represents the number o f fault types, a i,j denotes the probability of the i-th word belong to the j-th catalogue. (Note: a i,j is only a relat ive probability factor, not the true probability.) The follo wing describes how to calculate a i,w . The value of a i,w is used to determine whether a word belongs to a certain type. As the value indicates the frequency of each word in a particu lar type of fault, we consider the probability of the word i in the fault type w as the ratio between the number of word i in type w and the total number o f words in type w. The basic formu la is described in [START_REF] Yuan | SherLog: error diagnosis by connecting clues from run-time logs[END_REF]. P(i, w) represents the probability that the i-th word appears in the fault type w, count(i, w) represents the times that the i-th word appears in the fault type w.

1, 1 1, 2 1, 2, 1 2, 2 2, , 1 , 2 ,
      1 , , , m j count i w P i w count j w    (3)
However, this formu la only considers the distribution of different words in the same fault type, but it ignores the same word between different fault types in the log template. Fo r examp le, if the word i only appears in the type w, then we believe that a log is very likely to belong to the fault type w as long as it contains the word i, even if how many time the word i appears in the type w. Therefore, we can amend the formula by adding a scale factor shown as [START_REF] Fu | Execution anomaly detection in distributed systems through unstructured log analysis[END_REF].

        , , log 1
sum i count i w K i w sum i         (4)
The sum(i) represents the number of times that the word i appears in all the fault type of the template library, namely in [START_REF] Xu | Detecting large-scale system problems by mining console logs[END_REF].

    1 , n i sum i count i t    (5)
K(i, w) indicates the importance of the wo rd i in the type w and it is in inverse proportion to the frequency of word i occurs in other types, i.e. if the occurrence that word i appears in the type w is more than in the other fault types, the word i is mo re important to determine whether the log is belong to type w.

Thereby, we can conclude that the probability coefficient a i,w is calculated as the product between the frequency of words in the fault type and the importance in the entire template, namely in [START_REF] James | Listening to your cluster with LoGS[END_REF].

    , ,, iw a P i w K i w  ( 6
)
Importing the equations ( 3) and ( 4) can obtain equation [START_REF] Jain | Extracting the textual and temporal structure of supercomputing logs[END_REF].

1 ( , ) ( , ) , log 1 1 ( ) ( , ) m j count i w count i w i w sum i count j w a                     (7)
The following describes the learning process by using the FK-matrix. After obtaining the results of manual classification, accord ing to the formula, we can calculate the matrix A by column. Next we will describe how to solve the w-th column as an example.

The program will count the total nu mber of words in type w fro m the entire library template. The amount is named as T(w), which is shown in [START_REF] Stearley | Bad words: Finding faults in Spirit's syslogs[END_REF].

  1 ( , ) m j T w count j w    (8)
To facilitate rev ised and updated the FK-matrix, a new ro w will be added in the FK-matrix to store T(w) for reducing the number of double counting.

For each word i in the fault type w, the program will calculate sum(i). Similar to the T(w), in order to reduce the number of calculations, the additional space will be used to store sum(i). Thus, there is an expanded FK-matrix adding the additional ro w and column for storing statistical information. The final matrix A is shown in [START_REF] Sandia | Towards informatic analysis of syslogs[END_REF].

            1, 1 1, 2 1, 2, 1 2, 2 2, , 1 , 2 , 1 2 12 n
n m m m n a a a sum a a a sum A a a a sum m T T T n          (9)
According to equation [START_REF] Jain | Extracting the textual and temporal structure of supercomputing logs[END_REF], we should calculate count(i, w) before calculating a i,w . Considering the value of count(i, w) will be changed with updating of FK-matrix, count(i, w) will be saved in the FX-matrix instead of a i,w in pract ice. UiLog will calculate a i,w until the occurrence probability of word is needed. Due to related variables has restored in the FK-matrix, it does not add any additional overhead.

After obtaining the FK-matrix through learning period, UiLog can classify log through FK-matrix. In the running period, ad min istrators might modify or add new categories to different temp late for coping changes of system or software environment. At that moment, UiLog will modify FK-matrix.

5) Classifying Log Information

Through the fault keyword matrix, we can easily classify the fau lt log on the system. In the running period of the system, the fault log collection module will send fault logs to the fault log analysis module fro m various components. In this step, UiLog will scan every log message to compute the probability of d ifferent fault types of log, according the Fault Keyword matrix. Suppose the fault log L need to be classified, UiLog will co mpare each word in L with the FK-matrix. It will use an array (array s) to store the probability that every word in L belongs to different fault types. The Algorithm 1 gives the pseudo-code for calculating the probability.

Algorithm 1.

Log classification 1 function Classification()

2
for w in counts of fault types 3 for i ← every word in w 4

if ( i L ) then s[w] ← s[w] + A[w][i]

Fault Log Correlation Analysis

Fault log correlation analysis module diagnoses faults generated by different components of system and software through logs. It will use the result of log analysis to find connection between different fault logs.

In the UiLog, we use the results of our previous fault log classification to improve the trad itional fau lt correlat ion method based on time. We note that different fault type has different time range to affect system. For examp le, a hardware fau lt has a relatively small range of t ime to affect the system, but it has a huge impact on the system within a short time. On the other hand, the time range of the influence of network fau lt is relat ively wide. The associated component will produce a fault report after a long time. Thus, we can use different size of time windows to diagnose different fault types. It can improve the accuracy of judgment homologous fault.

As a practical application of the fault log classification, we collect log informat ion generated by all hosts and virtual machines in StrongCloud within 10 months and analyze the fault happened in this time. Here we first use traditional fault correlation method based on time to diagnosis the fault. It uses the un iform window size. Then we manually analyze each log in the tuple to determine th e different window size for every fault type.

In the specific process, the ad min istrator will point out which log they want to diagnose. UiLog will query the fault keyword matrix for finding the appropriate time window after confirming the fault type of log. Then UiLog will use this time window to diagnose the fault.

Performance Evaluations

In this section, we test the main functions of UiLog system and evaluate the effect including fault log analysis . We construct a fault-tolerance testbed, StrongCloud [19], which is made up of five Inspur's NF5240M 3 servers, each with 24 Intel Xeon E5 -2420 1.90 GHz processors, 32 GB of RAM and four Gigabit Ethernet ports. Each host has a Domain0 with CentOS release 6.3 of kernel 3.7.1-xen.x86_64, and the hypervisor is Xen 4.2.1. UiLog is a sub-system of StrongCloud. The evaluations show the importance of log management.

Log Analysis

For log analysis module, the most important is to test the efficiency of log classification. In this test, UiLog classifies all the fault logs collected by log collection module. The sample data is derived from StrongCloud within one year. The details are shown in Table 3. We first use the data of January and February to train UiLog to obtain the basic Fault Keyword matrix. The test is launched on the data from March 2013.

Fig. 3. Types of logs

To classify logs, we use trigger mechanism and User Defined Function (UDF) in database for analysis log timely. Whenever there is a log from the log collection mo dule, the trigger is activated and UDF function will call log analysis module. The log analysis module will classify the fault log according to the described steps of running period by using the fault keyword matrix. If the type of log cannot be determined, this log will be saved in unprocessed database and UiLog will inform the admin istrator to manually classify this log. After dealing with this unsorted log artificially, UiLog will automat ically learn fro m these new results and modify the fault keyword matrix to improve classification accuracy.

The bar chart (Fig. 3) shows the classification results for the one year fault log of StrongCloud. We divide fault log into 12 different types. Fig. 3 gives the ratio of each type. From the figure we can conclude that the main function of StrongCloud platform is network-related. The net error and web error are sum up to 43% of all faults. In addition, we can find that the platform has been subject to external attack, network authentication faults occupy 15%. The port SSL services or disable services are not commonly used.

Fro m Fig. 3, we can demonstrate that through the fault log classification analysis, UiLog can help administrator to manage computer systems, finding pro blems and bottlenecks in the system to compensate the deficiencies in the system.

Learning Efficiency

The learning efficiency of artificial classification is a vary import indicator for evaluation the new fault classification method of UiLog. We use two months fault logs as a sample for learn ing period. When the system is running, we will classify those non-classification logs manually at the end of every month. UiLog will re-learn these new classification results to improve the Fault Key matrix for better automatic classification result. Fig. 4 illustrates the number of classified fault logs and unclassified fault logs fro m March to October. It shows that there are only 568 fault log types cannot be judged in March after learning compared 2350 defined types. Then UiLog studies the results of new type rules by the end of March. The number of unclassified fault logs had dropped significantly, reaching 208 types on April. After that, the apparent decreasing tendency can be seen during March to June, reaching the lowest point (48 types) in June. That is because of re-learning by every end of month. Through the above experiments and analysis, our fault classification method is effective. It can be sufficiently carried out using the results of the automatic learn ing to automatically determine the type of fault log.

Conclusions and Future Work

With the development of cloud computing, the architecture of system and software are more co mplicated than before. This paper presents an integrated fault log analysis platform UiLog system, help ing administrators to manage the log generated by the all the components of system, mon itoring the running of the system and diagnosing the fault.

To effectively analy ze the logs, we propose a new method to classify log into different catalogs according to the different fau lt types. We use Fault Keyword matrix to accelerate the speed of classificat ion. In addition, we imp rove the fault correlation analysis. We use the result of fault classification to fix the time correlation window to reduce the truncation error and collision error.

Fig. 2 .

 2 Fig. 2. Number of key words appeared in the fault logs

Fig. 4 .

 4 Fig. 4. Number of clustersHowever, Fig.4also presents that the percentage of unclassified fault logs has a growth in July. According the previous analysis, it is due to the new applic ation or software developed in July. This change generates a lot of new fault logs and leads to 324 new fault types. After finishing the re-learning step by the end of July, the percentage of unclassified logs has dropped dramatically and the system become s gradually stabilizing. There are only 3 undefined log typ es.Through the above experiments and analysis, our fault classification method is effective. It can be sufficiently carried out using the results of the automatic learn ing to automatically determine the type of fault log.

Table 1 .

 1 Example of log messages

					Administrator
					Semantics, grammar
					System knowledge
		Raw Logs		Complex logic
					Fault Alarm
	Information extraction	Scan	Collect Filter		Human process
	and classification	Monitor system running			Experience
			Structuralization	Fault template
		System running report	Log classification
	Log				Artificial
	correlation		Correlation analysis	Knowledge library	analysis Manually
					Validation	analysis
					tools
			Analysis report	Analysis report

analysis Fig. 1. Workflow for log analysis and fault diagnosis

Table 2 .

 2 Fault types

	Content of Log			Label	Meaning
	INFO: task * blocked * more * |NUM| seconds. 14	Disk
	udevd (|NUM|): |DIR| is deprecated, please use	11	File
	|DIR| instead.
	pam_succeed_if(*):	error	retrieving	7	Authenticate
	information about user *
	Kernel reported iSCSI connection |NUM| error	6	Drive
	* state
	* received packet with * address * source	5	Network

Table 3 .

 3 Experimental environment

	System	Beginning	Ending	Days Size	Messages
	CentOS	2012-12	2013-12	386	1.7GB	14450302

Acknowledg ments.

This work is supported by National 973 Fundamental Basic Research Program under grant No. 2014CB340600 and National Science Foundation of China under grant No. 61272072.