A Functional Central Limit Theorem for a Class of Interacting Markov Chain Monte Carlo Methods

We present a functional central limit theorem for a new class of interacting Markov chain Monte Carlo algorithms. These stochastic algorithms have been recently introduced to solve non-linear measure-valued equations. We provide an original theoretical analysis based on semigroup tech-niques on distribution spaces and ﬂuctuation theorems for self-interacting random ﬁelds. Addi-tionally we also present a series of sharp mean error bounds in terms of the semigroup associated with the ﬁrst order expansion of the limiting measure-valued process. We illustrate our results in the context of Feynman-Kac semigroups.


Nonlinear Measure-Valued Equations and Self Interacting Markov chains
Let (S (l) ) l≥0 be a sequence of measurable spaces endowed with some σ-fields ( (l) ) l≥0 , and for every l ≥ 0, denote by (S (l) ) the set of all probability measures on S (l) . Let also (π (l) ) l≥0 be a sequence of probability measures on S (l) satisfying a nonlinear equation of the following form Φ l (π (l−1) ) = π (l) (1.1) for some mappings Φ l : (S (l−1) ) → (S (l) ).
Being able to solve these equations numerically has numerous applications in nonlinear filtering, global optimization, Bayesian statistics and physics as it allow us to approximate any sequence of fixed "target" probability distributions (π (l) ) l≥0 . For example, in a nonlinear filtering framework π (l) corresponds to the posterior distribution of the state of an unobserved dynamic model at time l given the observations collected from time 0 to time l. In an optimization framework, π (l) could correspond to a sequence of annealed versions of a distribution π that we are interested in maximizing. When Φ l is a Feynman-Kac transformation [2], there has been much work on mean field interacting particle interpretations of measure-valued equations of the form (1.1); see for example [2; 6]. An alternative iterative approach to solve these equations named Interacting Markov chain Monte Carlo methods (i-MCMC) has recently been proposed in [1; 3].
Let us give a rough description of i-MCMC methods. At level l = 0, we run a standard MCMC chain say X (0) n with a prescribed simple target distribution π (0) . Then, we use the occupation measure of the chain X (0) to design a second MCMC chain say X (1) n with a more complex target limiting measure π (1) . These two mechanisms are combined online so that the pair process interacts with the occupation measure of the system at each time. More formally, the elementary transition of the chain X (1) at a current time n depends on the occupation measure of the chain X (0) p from the origin p = 0, up to the current time p = n. This strategy can be extended to design a series of MCMC samplers (X (l) ) l≥0 targeting a prescribed sequence of distributions (π (l) ) l≥0 of increasing complexity. These i-MCMC samplers are self interacting Markov chains as the Markov chain X [m] n = (X (l) n ) 0≤l≤m associated with a fixed series of m levels, evolves with elementary transitions which depend on the occupation measure of the whole system from the origin up to the current time.
The theoretical analysis of these algorithms is more complex than the one of mean field interacting particle algorithms developed in [2]. In a pair of recent articles [1; 3], we initiated the theoretical analysis of i-MCMC algorithms and we provided a variety of convergence results including exponential estimates and an uniform convergence theorem with respect to the index l. However, we did not present any central limit theorem. The main purpose of this paper is to provide the fluctuation analysis of the occupation measures of a specific class of i-MCMC algorithms around its limiting values π (l) , as the time parameter n tends to infinity. We also make no restriction on the state spaces and we illustrate these models in the context of Feynman-Kac semigroups.
The rest of the paper is organized as follows: the i-MCMC methodology is described in section 1.2. For the convenience of the reader, we have collected in the short preliminary section 1.3 some of the main notation and conventions used in the paper. We also recall some regularity properties of integral operators used in this article. Section 2 is devoted to the two main results of this work. We provide a multivariate functional central limit theorem together with some sharp non asymptotic mean error bounds in terms of the semigroup associated with a first order expansion of the mappings Φ l . To motivate our approach, we illustrate in section 3 these results in the context of Feynman-Kac semigroups. The last four sections are essentially concerned with the proof of the two theorems presented in section 2. Our analysis is based on two main ingredients; namely, a multivariate functional central limit theorem for local interaction fields and a multilevel expansion of the occupation measures of the i-MCMC model at a given level l in terms of the occupation measures at the levels with lower indexes. These two results are presented respectively in section 4 and section 5. The final section 6 and section 7 are devoted respectively to the proofs of the functional central limit theorem and the sharp non asymptotic mean error bounds theorem. In appendix, we have collected the proofs of a series of technical lemmas.

Description of the i-MCMC methodology
To introduce the i-MCMC methodology, we first recall how the traditional Metropolis-Hastings algorithm proceeds. This algorithm samples a Markov chain with the following transition kernel where P(x, d y) is a Markov transition on some measurable space S, and G : S 2 → + . Let π be a probability measure on S, and let (π × P) 1 and (π × P) 2 be the pair of measures on (S × S) defined by (π × P) 1 (d(x, y)) = π(d x) P(x, d y) = (π × P) 2 (d( y, x)).
We further assume that (π × P) 2 ≪ (π × P) 1 and the corresponding Radon-Nykodim ratio is given by In this case, π is an invariant measure of the time homogeneous Markov chain with transition M . One well-known difficulty with this algorithm is the choice of P.
Next, assume π = Φ(η) is related to an auxiliary probability measure η on some possibly different measurable space S ′ with some mapping Φ : We also suppose that there exists an MCMC algorithm on S ′ whose occupation measure, say η n , approximates the measure η as the time parameter n increases. At each current time, say n, we would like to choose a Metropolis-Hastings transition M n associated to a pair (G n , P n ) with invariant measure Φ(η n ). The rationale behind this is that the resulting chain will behave asymptotically as a Markov chain with invariant distribution Φ(η) = π, as soon as η n is a good approximation of η.
The choice of (G n , P n ) is far from being unique and we present a few solutions in [3]. A natural and simple one consists of selecting P n (x, d y) = Φ(η n )(d y) which ensures that G n is equal to one. This is the class of i-MCMC algorithms studied in this paper. Iterating this strategy, we design a sequence of i-MCMC chains X (k) = (X (k) n ) n≥0 on every level set S (k) whose occupation measures η (k) n approximate the solution π (k) of the system (1.1) for every k ≥ 0.
To describe more precisely the resulting algorithm, we need a few notation. We introduce a sequence of initial probability measures ν k on S (k) , with k ≥ 0, and we denote by the sequence of occupation measures of the chain X (k) at level k from the origin up to time n, where δ x stands for the Dirac measure at a given point x ∈ S (k) . The i-MCMC algorithm proposed here is defined inductively on the level parameter k. First, we shall suppose that ν 0 = π (0) and we let X (0) = (X (0) n ) n≥0 be a collection of independent random variables with common distribution ν 0 = π (0) . For every k ≥ 1, and given a realization of the chain X (k) , the (k + 1)-th level chain X (k+1) = (X (k+1) n ) n≥0 is a collection of independent random variables with distributions given, for any n ≥ 0, by the following formula where we use the convention that we have Moreover, d(x 0 , . . . , x n ) = d x 0 × . . . × d x n stands for an infinitesimal neighborhood of a generic path sequence (x 0 , . . . , x n ) ∈ (S (k+1) ) (n+1) .

Notation and conventions
We denote respectively by (E), 0 (E), (E), and (E), the set of all finite signed measures on some measurable space (E, ), the convex subset of measures with null mass, the set of all probability measures, and the Banach space of all bounded measurable functions f on E. We equip (E) with the uniform norm f = sup x∈E | f (x)|. We also denote by 1 (E) ⊂ (E) the unit ball of functions f ∈ (E) with f ≤ 1, and by Osc 1 (E), the convex set of -measurable functions f with oscillations less than one, which means that We denote by the Lebesgue integral of a function f ∈ (E) with respect to µ ∈ (E). We slightly abuse the notation, and we denote by µ(A) = µ(1 A ) the measure of a measurable subset A ∈ .
We recall that a bounded integral operator M from a measurable space (E, ) into an auxiliary measurable space (F, ), is an operator f → M ( f ) from (F ) into (E) such that the functions are -measurable and bounded for any f ∈ (F ). A bounded integral operator M from a measurable space (E, ) into an auxiliary measurable space (F, ) also generates a dual operator µ → µM from (E) into (F ) defined by (µM )( f ) = µ(M ( f )). We denote by the norm of the operator f → M ( f ) and we equip the Banach space (E) with the corresponding total variation norm µ = sup We also denote by β(M ) the Dobrushin coefficient of a bounded integral operator M defined as When M has a constant mass M (1)(x) = M (1)( y), for any (x, y) ∈ E, the operator µ → µM maps 0 (E) into 0 (F ), and β(M ) coincides with the norm of this operator. We equip the sets of sequence of distributions (E) with the uniform total variation distance defined for all η = We extend a given bounded integral operator µ ∈ (E) → µM ∈ (F ) into a mapping Sometimes, we slightly abuse the notation and we denote by ν instead of (ν) n≥0 the constant distribution flow equal to a given measure ν ∈ (E). Finally we denote by c(k) with k ∈ , a generic constant whose values may change from line to line, but they only depend on the parameter k. For any pair of integers 0 ≤ m ≤ n, we denote by (n) m := n!/(n − m)! the number of one to one mappings from {1, . . . , m} into {1, . . . , n}. Finally, we shall use the usual conventions = 0 and = 1.

Fluctuations theorems
This section is mainly concerned with the statement of the two main theorems presented in this article. These results are based on a first order weak regularity condition on the mappings (Φ l ) appearing in (1.1). We assume that the mappings Φ l : (S (l−1) ) → (S (l) ) satisfy the following first order decomposition for any l ≥ 1 In the formula above, l,η is a collection of bounded integral operators from S (l−1) into S (l) , indexed by the set of probability measures η ∈ (S (l−1) ) and Ξ l (µ, η) is a collection of remainder signed measures on S (l) indexed by the set of probability measures µ, η ∈ (S (l−1) ). We also require that for some integral operator Ξ l from (S (l) ) into the set 2 (S (l−1) ) of all tensor product functions This weak regularity condition is satisfied in a variety of models including Feynman-Kac semigroups discussed in the next section. We also mention that, under weaker assumptions on the mappings Φ l , we already proved in [1, Theorem 1] and [3, Theorem 1.1.] that for every l ≥ 0 and any function f ∈ (S (l) ), we have the almost sure convergence result The article [3, Theorem 1.1.] also provides exponential inequalities and uniform estimates with respect to the index l.
In order to describe precisely the fluctuations of the empirical measures η (l) n around their limiting values π (l) , we need to introduce some notation. We denote by D l the first integral operator l,π (l−1) associated with the target measure π (l−1) , and we set D k,l with 0 ≤ k ≤ l for the corresponding semigroup. More formally, we have D l = l,π (l−1) and for all 1 ≤ k ≤ l, For k > l, we use the convention D k,l = I d , the identity operator. The reader may find in section 3 an explicit functional representation of these semigroups in the context of Feynman-Kac models.
We now present the functional central limit theorem describing the fluctuations of the i-MCMC process around the solution of equation n which are i.i.d. samples from π (0) when k = 0 and generated according to (1.4) otherwise.
Theorem 2.1. For every k ≥ 0, the sequence of random fields (U (k) n ) n≥0 on (S (k) ) defined by U (k) n := n η (k) n − π (k) converges in law, as n tends to infinity and in the sense of finite dimensional distributions, to a sequence of Gaussian random fields U (k) on (S (k) ) given by the following formula where V (l) is a collection of independent and centered Gaussian fields with a covariance function given for any ( f , g) ∈ (S (l) ) 2 by We now recall that if Z is a Gaussian (0, 1) random variable, then for any m ≥ 1, where Γ stands for the standard Gamma function. Consequently, in view of Theorem 2.1, the reader should be convinced that the estimates presented in the next theorem are sharp with respect to the parameters m and k.

Theorem 2.2.
For any k, n ≥ 0, any f ∈ Osc 1 (S (k) ), and any integer m ≥ 1, we have the non asymptotic mean error bounds with the collection of constants a(m) given by a(2m) 2m = (2m) m and a(2m and for some constants b(m) (respectively c(k)) whose values only depend on the parameter m (respectively k).

Feynman-Kac semigroups
In this section, we shall illustrate the fluctuation results presented in this paper with the Feynman-Kac mappings (Φ l ) given for all l ≥ 0 and all (µ, f ) ∈ ( (S (l) ) × (S (l+1) )) by where G l is a positive potential function on S (l) and L l+1 is a Markov transition from S (l) to S (l+1) . In this scenario, the solution of the measure-valued equation (1.1) is given by where (Y l ) l≥0 stands for a Markov chain taking values in the state spaces (S (l) ) l≥0 , with initial distribution π (0) and Markov transitions (L l ) l≥1 .
These probabilistic models arise in a variety of applications including nonlinear filtering and rare event analysis as well as spectral analysis of Schrödinger type operators and directed polymer analysis. These Feynman-Kac distributions are quite complex mathematical objects. For instance, the reference Markov chain may represent the paths from the origin up to the current time of an auxiliary sequence of random variables Y ′ l taking values in some state spaces E ′ l .
To be more precise, we have To get an intuitive understanding of i-MCMC algorithms in this context, we note that for the Feynman-Kac mappings (3.1) we have for all f ∈ (S (k) ) It follows that each random state X (k) p is sampled according to two separate genetic type mechanisms. First, we select randomly one state X (k−1) q at level (k − 1), with a probability proportional to . Second, we evolve randomly from this state according to the exploration transition L k . This i-MCMC algorithm model can be interpreted as a spatial branching and interacting process. In this interpretation, the k-th chain duplicates samples with large potential values, at the expense of samples with small potential values. The selected offspring evolve randomly from the state space S (k−1) to the state space S (k) at the next level. The same description for path space models (3.2) coincides with the evolution of genealogical tree based i-MCMC algorithms.
For this class of Feynman-Kac models, we observe that the decomposition (2.1) is satisfied with the first order integral operator D l,η defined for all f ∈ (S (l) ) by and the remainder measures Ξ l (µ, η) given by We can observe that the regularity condition Indeed, it is easy to check that for any function f ∈ 1 (S (l) ), Finally, we mention that the semigroup D k,l introduced above can be explicitly described in terms of the semigroup The Dobrushin coefficient β(D k,l ) can also be computed in terms of these semigroups. We have with the Markov integral operator P k,l given, for all 1 ≤ k ≤ l, by Therefore, it follows that

Introduction
As for mean field interacting particle models, see section 9.2 in [2], the local errors induced by the i-MCMC algorithm can be interpreted as small disturbances of the measure-valued equation (1.1).
To be more precise, we need a few additional definitions. (S (l) ) defined for all n ≥ 0 and all l ≥ 0 by For n = 0, we use the convention Φ l η In our context, the i-MCMC process at level l consists in a series of conditionally independent random variables X (l) n with different distributions Φ l (η By definition of the i-MCMC algorithm, we have (∆ (l) n ( f )) = 0 for any f ∈ (S (l) ). In addition, we have Recalling (2.5), as η (l) n ( f ) converges almost surely to π (l) ( f ), as n → ∞, for every l ≥ 0 and any function f ∈ (S (l) ), we deduce that Φ l (η (l−1) n )( f ) converges almost surely to π (l) ( f ), as n → ∞, which implies, via Cesaro mean convergence that, as n → ∞, This shows that the random disturbances ∆ (l) n ( f ) are unbiased with finite variance. It is possible to obtain a much stronger result by applying the traditional central limit theorem for triangular arrays; see for instance Theorem 4, page 543 of [11]. More precisely, we find that ∆ (l) n ( f ) converges in law as n → ∞ to a centered Gaussian random variable with variance given by (4.2). If we rewrite ∆ (l) n in a different way, we have the following decomposition Equation (4.3) shows that the evolution equation of the sequence of occupation measures η (l) n can be interpreted as a stochastic perturbation of the following measure-valued equation initialized using conditions µ (0) n = η (0) n . Note that the constant distribution flow µ (l) n = π (l) associated with the solution of (1.1) also satisfies (4.4).
Although the above discussion gives some insight on the asymptotic normal behavior of the local errors accumulated by the i-MCMC algorithm, it does not give directly the precise fluctuations of the sequence of measures (η (l) n ) n≥0 around their limiting values π (l) . The analysis of these fluctuations is based on the way the semigroup associated with the evolution (4.4) propagates these local perturbation random fields. We can observe that the one step transformation of the equation (4.4) is an averaging of a nonlinear transformation Φ l of the complete sequence of measures µ (l−1) p from the origin p = 0 up to the current time n. Iterating these transformations, we end up with a complex semigroup on the sequence of measures between successive levels. We refer the interested reader to [3] for an "explicit" calculation of these semigroups in the context of Feynman-Kac models. In the further developments of section 5, we shall derive a first order expansion of these semigroups. These multilevel decompositions express the propagations of the local perturbation errors in terms of weighted random fields. The next section is devoted to the fluctuations analysis of a general class of weighted random fields associated with the i-MCMC algorithm. We present an admissible class of array type weight functions for which the corresponding weighted perturbation random fields behave as a collection of independent and centered Gaussian fields.

A martingale convergence theorem
This section is concerned with the fluctuations analysis of weighted random fields associated with the local perturbation random fields (4.1) discussed in section 4.1. Our theorem is basically stated in terms of the convergence of a sequence of martingales. To describe these processes, we need the following definitions. We observe that the standard sequence w n (p) = 1/ n belongs to , with the identity function ̟(ǫ) = ǫ.
Denote by f = ( f l ) l≥0 ∈ l≥0 (S (l) ) d a collection of d-valued functions, and for every 1 ≤ i ≤ d, let f i be the i-th coordinate collection of functions f i = ( f i l ) l≥0 ∈ l≥0 (S (l) ). We consider the dvalued and (n) -martingale (n) ( f ) = ( (n) ( f i )) 1≤i≤d defined for any l ≥ 0 and any 1 ≤ i ≤ d, We can extend this discrete generation process to the half line + = [0, ∞) by setting, for all t ∈ + , We can observe that ( (n) t ( f )) is a càdlàg martingale with respect to the filtration associated with the σ-fields (n) t generated by the random variables (X (k) p ) 0≤p≤n , with k ≤ ⌊t⌋, and X ⌊t⌋+1 p with 0 ≤ p < ⌊{t}(n + 1)⌋.
Before establishing the proof of this convergence result, we illustrate some direct consequences of this theorem.
On the one hand, it follows from this convergence that the discrete generation martingales ( (n) l ( f )), where l ∈ , converge in law as n → ∞ to an d -valued and discrete Gaussian martingale l ( f ) = ( l ( f i )) 1≤i≤d with bracket given, for any l ≥ 0 and any 1 ≤ i, j ≤ d, by On the other hand, we can fix a parameter l and introduce a collection of functions g In other words, all the functions g j k are null except on the diagonal k = j. By construction, for every 0 ≤ k ≤ l, we have We can deduce from the above martingale convergence theorem that the random fields (V ( j) n ) with 0 ≤ j ≤ l, converge in law, as n tends to infinity and in the sense of finite dimensional distributions, to a sequence of independent and centered Gaussian fields V ( j 0≤ j≤l such that, for any 1 ≤ i, i ′ ≤ d j , We can achieve this discussion by the following corollary of Theorem 4.4.

Corollary 4.5. For every collection of weight functions (w (l) ) ∈
, the sequence of random fields (V (l) n ) l≥0 converges in law, as n tends to infinity and in the sense of finite dimensional distributions, to a sequence of independent and centered Gaussian fields V (l) l≥0 with covariance function given for any l ≥ 0 and any ( f , g) ∈ (S (l) ) 2 by V (l) ( f )V (l) (g) = π (l) ( f − π (l) ( f )) (g − π (l) (g)) .
We now come back to the proof of Theorem 4.4.

Proof of Theorem 4.4:
For every 1 ≤ i ≤ d and any 0 ≤ k ≤ l such that f i k ∈ (S (k) ), the martingale increments are given by Our goal is to make use of the central limit theorem for triangular arrays of d -valued martingales given in [8] page 437. We first rewrite the martingale where for every 0 ≤ i ≤ l(n + 1) + n, with i = k(n + 1) + p for some 0 ≤ k ≤ l, and 0 ≤ p ≤ n We further denote by (n) i the σ-field generated by the random variables X (k) p for any pair of parameters (k, p) such that k(n + 1) + p ≤ i. By construction, for any sequence of functions f = ( f l ) l≥0 and g = (g l ) l≥0 ∈ l≥0 (S (l) ) and for every 0 ≤ i ≤ l(n + 1) + n, with i = k(n + 1) + p for some 0 ≤ k ≤ l, and 0 ≤ p ≤ n, we find that with the local covariance function This implies that Recalling (2.5), we have that η (k) n (h) converges almost surely to π (k) (h), as n → ∞, for every k ≥ 0 and any function h ∈ (S (k) ). Under our regularity conditions on the sequence of mappings (Φ k ), we find that Φ k (η (k−1) n )(h) converges almost surely to π (k) (h), as n → ∞, for any function h ∈ (S (k) ). It immediately implies the almost sure convergence .
Under our assumptions on the weight array functions (w (k) n (p)), we can deduce that almost surely Consequently, it leads to the almost sure convergence In order to go one step further, we introduce the sequence of d -valued and continuous time càdlàg random processes given, for any f = ( f l ) l≥0 ∈ l≥0 (S (l) ) d and t ∈ + , by It is not hard to see that for any l ∈ , and any 0 ≤ k ≤ n, and ⌊t(n + 1)⌋ − ⌊t⌋(n + 1) = ⌊{t}(n + 1)⌋ = k where {t} = t − ⌊t⌋ stands for the fractional part of t ∈ + . Using the fact that ⌊t(n + 1)⌋ + n = (⌊t⌋(n + 1) + n) + ⌊{t}(n + 1)⌋, we obtain the decomposition This readily implies that Furthermore, using the same calculations as above, for any pair of integers 1 ≤ j, j ′ ≤ d, we find that Under our assumptions on the weight array functions (w (k) n (p)), we find that Hereafter, we have for every 0 ≤ i ≤ l(n + 1) + n In addition, we also have lim n→∞ sup 0≤k≤l sup 0≤p≤n w (k) n (p) = 0.
Consequently, we conclude that the conditional Lindeberg condition is satisfied. Hence, the dvalued martingale ( (n) t ( f )) converge in law, as n → ∞, to a continuous martingale t ( f ) with predictable bracket given, for any air of indexes 1 ≤ j, j ′ ≤ d, by which completes the proof of Theorem 4.4.

A multilevel expansion formula
We present in this section a multilevel decomposition of the sequence of occupation measures (η (k) n ) around their limiting value π (k) . In section 4.1, we have shown that the sequence (η (k) n ) satisfies a stochastic perturbation equation (4.3) of the measure-valued equation (4.4). In this interpretation, the forthcoming developments provide a first order expansion of the semigroup associated to equation (4.3). This result will be used in various places in the rest of the paper. In section 6, we combine these first order developments with the local fluctuation analysis presented in section 4 to derive a "two line" proof of the functional central limit theorem 2.1. In section 7, we use these decompositions to prove the non asymptotic mean error bounds presented in Theorem 2.2. The forthcoming multilevel expansion is expressed in terms of the following time averaging operators. Definition 5.1. For any l ≥ 0, denote by the mapping : η ∈ (S (l) ) → (η) = ( n (η)) n≥0 ∈ (S (l) ) defined for any sequence of measures η = (η n ) ∈ (S (l) ) , and any n ≥ 0, by We also denote by k = • k−1 the k-th iterate of the mapping .
We are now in position to state the following pivotal multilevel expansion.

Proposition 5.2.
For every k ≥ 0, we have the following multilevel expansion with a sequence of signed random measures Ξ (k) = (Ξ (k) n ) such that, for any m ≥ 1, for some constants b(m) (respectively c(k)) whose values only depend on the parameter m (respectively k).

Proof of Proposition 5.2:
Recall that we use the convention D 1,0 = I d . We prove the mean error bounds by induction on the integer parameter k. For k = 0, we readily find that Hence, (5.2) is satisfied with the sequence of null measures Ξ (0) n = 0. We further assume that the expansion is satisfied at rank k. We shall make use of the following decomposition It follows from our assumptions on the mappings (Φ k ) that for all µ, η ∈ (S (k) ), Using the fact that On the other hand, using the first order expansion of the mapping Φ k+1 , we have the decomposition ) be a tensor product function associated with a subset I ⊂ , a pair of functions (h 1 i , h 2 i ) i∈I ∈ ( (S (k) ) 2 ) I and a sequence of numbers (a i ) i∈I ∈ I . Using the generalized Minkowski integral inequality, we find that . Moreover, via the mean error bounds established in [3, Theorem 1.1.], it follows that for any i ∈ I and j = 1, 2 for some finite constant a(m) whose values only depend on the parameter m. These estimates readily implies that for some finite constant b(m) whose values only depend on the parameter m. This implies that for any f ∈ 1 (S (k) ) To take the final step, we use the induction hypothesis to check that This yields the decomposition n (δ (k+1) ) + In the last assertion, we have used the convention that D k+2,k+1 = I d stands for the identity operator. We also note that The end of the proof of the result at rank (k + 1) is now a consequence of decomposition (5.4) in conjunction with the estimates (5.5) and (5.6). This completes the proof of Proposition 5.2.

Proof of the central limit theorem
This section is mainly concerned with the proof of the functional central limit theorem presented in section 2.1. We first need to express the time averaging semigroup k introduced in definition 5.1 in terms of the following weighted summations k n (η) = 1 n + 1 0≤p≤n s (k) n (p) η p with the weight array functions s (k) n = (s (k) n (p)) 0≤p≤n defined by the initial condition s (1) n (p) = 1 and, for any k ≥ 1 and, for all n ≥ 0 and 0 ≤ p ≤ n, by the recursive equation The proof of Theorem 2.1 is a direct consequence of the following proposition whose proof is postponed to the end of the section.
We shall now proceed to the "two line" proof of Theorem 2.1.
Proof of Theorem 2.1: Fix a parameter k ≥ 0 and for any 0 ≤ l ≤ k, let W (l) be the distribution flow mappings associated with the weight functions (w ((k−l)+1) ) defined in Proposition 6.1 and given by By construction, we have for any 0 ≤ l ≤ k Using the multilevel expansion presented in Proposition 5.2, we obtain that where we recall that, for any k ≥ 0, V (k) n = W (k) n (δ (k) ). Finally, Theorem 2.1 is a direct consequence of Corollary 4.5 together with the mean error bound (5.3) with m = 1 and (6.1).
The proof of Proposition 6.1 relies on two pivotal lemmas. Their proofs follow elementary techniques but require tedious calculations which are given in appendix.

Lemma 6.2.
For any k ≥ 0, and any 0 ≤ p ≤ n, we have the formula where the remainder sequence (r (k) n (p)) is such that r (0) n (p) = r (1) n (p) = 0 and, for all k ≥ 2 and 0 ≤ p ≤ n, we have with |ρ (k) (n)| ≤ c(k) log (n + 1) k . In addition, for any increasing sequence of integers (κ n ) such that κ n ≤ n and κ n /n → κ > 0, as n → ∞, we have We are now in position to prove Proposition 6.1.
We achieve the proof of (6.3) by a direct application of Lemma 6.3. From the above discussion, we also find that 0≤p≤κ n n | ≤ |r (k) n |. We deduce from Lemma 6.3 that for any increasing sequence of integers (κ n ) such that κ n ≤ n and κ n /n → κ > 0, lim n→∞ 1 n 0≤p≤κ n s (k) n (p) 2 = (2(k − 1))!

Proof of the sharp m -estimates
This short section is concerned with the proof of non asymptotic estimates presented in Theorem 2.2. We start with a pivotal technical lemma.
In summary, we have proved that with some remainder sequence satisfying r (k) n (p) ≤ c(k) (log (n + 1)) k−2 p≤q≤n 1 (q + 1) 2 which completes the proof of Lemma 6.2.
We observe that the first inequality in the left hand side also holds true for p = 0. Moreover, using the fact that for any p ≥ 2, we have Finally, the second assertion of Lemma 6.3 follows from (7.2) and (7.3).