Patrice Bertail

Stéphan Clémençon

Approximate regenerative-block bootstrap for Markov chains: some simulation studies

Keywords: Markov chains, Bootstrap, Regeneration, Small set, Nummelin splitting technique, Simulation PACS: 62E20, 62D05, G2N01

, a novel methodology for bootstrapping general Harris Markov chains has been developed, the (approximate) regenerative block-bootstrap. It is built on the renewal properties of the chain (or of a Nummelin extension of the latter) and has theoretical properties that surpass other existing methods within the Markovian framework. This paper is devoted to discuss practical issues related to the implementation of this specific resampling method and to present various simulation studies for investigating its performance and comparing it to other bootstrap resampling schemes, standing as natural candidates in the Markov setting.

Introduction

Adapting the naive bootstrap method introduced by [START_REF] Efron | Bootstrap methods: another look at the jackknife[END_REF] in the i.i.d. setting, so as to produce asymptotically valid procedures for dependent data, and time series in particular, constitutes an intense field of research, still developing. The principle underlying such generalizations generally consists in resampling whole blocks of observations instead of single observed values in order to mimic the data dependence (see [START_REF] Lahiri | Resampling methods for dependent Data[END_REF] and the references therein). The moving-block bootstrap (MBB) illustrates this idea, it lies in resampling (overlapping or disjoint) data blocks of fixed length to capture the dependence structure of the observations (refer to [START_REF] Bühlmann | Bootstrap for time series[END_REF] and [START_REF] Politis | The impact of bootstrap methods on time series analysis[END_REF] for recent surveys). Although this approach may yield consistent procedures in many weakly dependent settings, it has several important drawbacks. First, stationarity of the observations is usually required by the validity framework of the MBB approach. Furthermore, implementing the MBB method calls for a preliminary estimation of the bias and of the asymptotic variance of the statistic of interest. This makes its application difficult in practice (see [START_REF] Götze | Second order correctness of the blockwise bootstrap for stationary observations[END_REF]). And from a theoretical viewpoint, the rate of convergence of the MBB distribution is slower than the one of the i.i.d. bootstrap: at best it is of order O P (n -3/4 ) under restrictive conditions (stipulating that all moments are finite and that strong mixing coefficients decrease exponentially fast). Eventually, the results obtained highly depend on the choice of the block size. Except in very particular situations (see [START_REF] Götze | Second order correctness of the blockwise bootstrap for stationary observations[END_REF]), no general method for determining the adequate block size has yet been developed. In [START_REF] Hall | On blocking rules for the bootstrap with dependent data[END_REF], [START_REF] Bühlmann | Block length selection in the bootstrap for time series[END_REF] and [START_REF] Politis | Automatic Block-Length Selection for the Dependent Bootstrap[END_REF], data-driven methods for selecting the block size are proposed, aiming at approximating the variance, rather than at estimating directly the studentized distribution in a precise fashion (see [START_REF] Lahiri | Resampling methods for dependent Data[END_REF]).

Various approaches for bootstrapping particular types of stationary Markov chains have recently been proposed in the statistical literature. If a parametric Markovian model is a priori specified (an ARMA model for instance), the problem simply reduces then to sampling randomly estimated centered residuals (see [START_REF] Bose | Edgeworth correction by bootstrap in autoregressions[END_REF]). Following this idea, [START_REF] Bühlmann | Sieve Bootstrap for time series[END_REF] introduced a sieve bootstrap method based on approximating the time series by some AR model with large lagorder. This method presents both promising theoretical results and good practical performance at the same time, but is well suited to linear stationary time series rather than to general Markov chains. Another approach to bootstrapping Markov chains follows the proposal of [START_REF] Athreya | Bootstrapping Markov chains: countable case[END_REF] (see also [START_REF] Rajarshi | Bootstrap in Markov-sequences based on estimates of transition density[END_REF]) in the finite case, which uses a nonparametric estimate of the transition probability for reproducing the markovian structure of the data series (see also the local Markov bootstrap in [START_REF] Paparoditis | The local bootstrap for Markov processes[END_REF] based on similar ideas).

Following the approach proposed by [START_REF] Datta | Regeneration-based bootstrap for Markov chains[END_REF], it has been shown in [START_REF] Bertail | Regenerative Block Bootstrap for Markov Chains[END_REF], that a specific resampling technique for bootstrapping some statistics of regenerative Markov chains offers attracting advantages both regarding to asymptotic second order properties and from a practical viewpoint. This method, the Regenerative Block Bootstrap (RBB), consists in resampling (a random number of) data blocks corresponding to 'cycles' of the observed path (i.e. data segment between consecutive regeneration times) until the length of the reconstructed trajectory is larger than the initial one. In the general (non-regenerative) setting, [START_REF] Bertail | Regenerative Block Bootstrap for Markov Chains[END_REF] proved that this principle can still be successfully applied in numerous situations, by adding a preliminary stage consisting in estimating the distribution of a regenerative extension of the chain. Due to the approximation step, we call this method Approximate Regenerative Block Bootstrap (ARBB).

The purpose of the paper is to describe the mechanic of the (A)RBB method from a practical angle and evaluate empirically its performance in comparison with major competitors. It is organized as follows. In section 2, basics concern-ing the regenerative method in connection with the Markov chain theory and the Nummelin splitting technique are briefly recalled. Section 3 is devoted to the description of the (A)RBB algorithm. Practical issues related to its implementation are discussed as well. The performance of the (A)RBB methodology is then investigated through several applications: simulation studies are carried out in section 4, which show on some (regenerative and pseudo-regenerative) examples of Markov chains, arising from operational research or standard time series analysis, performs, when compared to natural bootstrap method competitors. In section 5, some concluding remarks are collected, together with several lines of further research.

Theoretical Background

In what follows X = (X n ) n∈N denotes an aperiodic Harris Markov chain on a countably generated state space (E, E), with transition probability Π, and initial probability distribution ν (see [START_REF] Revuz | Markov chains[END_REF] for an account of the Markov chain theory). We also denote by P ν (respectively by P x for x in E) the probability measure on the underlying space such that X 0 ∼ ν (resp. X 0 = x), by E ν [.] the P ν -expectation (resp. by E x [.] the P x -expectation) and by I{A} the indicator function of the event A.

We now recall key notions, concerning the regenerative method and its application to the analysis of the behavior of general Harris chains via the Nummelin splitting technique (refer to [START_REF] Nummelin | General irreducible Markov chains and non negative operators[END_REF], for further detail).

Regenerative Markov chains Here we assume that the chain X possesses an accessible atom, i.e. a Harris set A such that for all (x, y) ∈ A 2 , Π(x, .) = Π(y, .). Denote by τ A = τ A (1) = inf {n ≥ 1, X n ∈ A} the hitting time on A, by τ A (j) = inf {n > τ A (j -1), X n ∈ A} for j ≥ 2, the successive return times to A, and by E A [.] the expectation conditioned on X 0 ∈ A.

Regeneration blocks

From the strong Markov property it is immediate that, for any initial distribution ν, the sample paths of the chain may be divided into i.i.d. blocks of random length corresponding to consecutive visits to the atom A

B 1 = X τ A (1)+1 , ..., X τ A (2) , ..., B j = X τ A (j)+1 , ..., X τ A (j+1) , ... (1)
taking their values in the torus T = ∪ ∞ n=1 E n . The τ A (j)'s are thus successive random times at which the chain forgets its past, namely regeneration times. In this regenerative setting, the stochastic stability properties of the chain amount to properties concerning the speed of return time to the atom only. For instance, X is positive recurrent iff E A [τ A ] < ∞ (see Kac's theorem in [START_REF] Meyn | Markov chains and stochastic stability[END_REF]). In such a case the unique invariant probability distribution µ is the Pitman's occupation measure given by:

∀B ∈ E, µ(B) = 1 E A [τ A ] E A [ τ A i=1 I{X i ∈ B}]. (2)
In particular, parameters of interest for a positive Harris Markov chain may be expressed in many cases in terms of regeneration cycles only (mainly those related to the long term behaviour of the process, as indicated by ( 2)). Numerous examples are given in [START_REF] Bertail | Regeneration-based Statistics for Markov Chains[END_REF]).

The regenerative method We point out that first order limit results such as the Law of Large Numbers, the Central Limit Theorem or the Law of Iterated Logarithm for an additive functional i f (X i ) of a regenerative Harris positive chain X may be easily derived by applying the corresponding i.i.d. results to functionals of the i.i.d. regeneration blocks (B j ) j 1 (see [START_REF] Meyn | Markov chains and stochastic stability[END_REF] for such illustrations of the regenerative method introduced by [START_REF] Smith | Regenerative stochastic processes[END_REF]). However, when the matter is to establish higher order limit results (see [START_REF] Bölthausen | The Berry-Esseen Theorem for strongly mixing Harris recurrent Markov Chains[END_REF], [START_REF] Malinovskii | Limit theorems for Harris Markov chains I[END_REF] or [START_REF] Bertail | Edgeworth expansions for suitably normalized sample mean statistics of atomic Markov chains[END_REF] for refinements of the CLT in the Markovian setting), the fact that the data blocks

B 0 = (X 1 , ..., X τ A (1) ), B 1 , ..., B ln-1 , B (n)
ln = (X τ A (ln)+1 , ..., X n ) defined by the l n = 1 k n I{X k ∈ A} regeneration times over a trajectory of finite length n are not independent (the sum of the block lengths is of course n) is essential. The randomness of the number of blocks plays a crucial role in the distribution of any statistic based on a finite sample. This observation lies at the heart of the RBB procedure (see section 3, details are in [START_REF] Bertail | Regenerative Block Bootstrap for Markov Chains[END_REF]).

Regeneration-based statistics

In the time series framework, inference is generally based on a single trajectory X 1 , ..., X n of the Harris chain X. Therefore, in a nonstationary setting the distribution of the non-regenerative blocks, B 0 and B (n) ln , cannot be estimated from a single realization of the chain only (due to their dramatical dependence upon the unknown initial distribution). Furthermore, their contribution to the value of a functional T (X 1 , ..., X n ) of interest is generally significant, leading to first order bias terms in particular (see the discussion in [START_REF] Bertail | Edgeworth expansions for suitably normalized sample mean statistics of atomic Markov chains[END_REF]). Hence, statistics involving B 0 and B (n) ln must be avoided in practice, when estimating sampling distributions is essential (for building confidence intervals for instance, see [START_REF] Bertail | Regeneration-based Statistics for Markov Chains[END_REF]).

As an illustration, when X is a positive recurrent chain with limiting probability distribution µ, consider f : E → R a µ-integrable function. In the nonstationary case, when the matter is to recover the asymptotic mean µ(f

) = f dµ from data X 1 , ..., X n (notice that µ(f ) = E A [τ A ] -1 E A [ 1 i τ A f (X i )]
), together with assessing the accuracy of the estimation, rather than the standard sample mean µ n (f ) = n -1 1 k n f (X k ), it is preferable to use the truncated mean computed using data collected between the first and last regeneration times, namely

µ n (f ) = τ A (ln) k=1+τ A f (X k ) τ A (l n ) -τ A = ln-1 j=1 f (B j ) ln-1 j=1 L(B j ) , (3)
where

f (B j ) = τ A (j+1)
k=1+τ A (j) f (X k ) and L(B j ) = τ A (j + 1)τ A (j) for j 1, with the convention that µ n (f ) = 0 when l n 2. Furthermore, under suitable moment conditions (see the paragraph below) [START_REF] Bertail | Edgeworth expansions for suitably normalized sample mean statistics of atomic Markov chains[END_REF] have shown that the estimator µ n (f ) is asymptotically normal with asymptotic mean µ(f ) and variance 1, when standardized by the following sequence, based on regeneration data blocks as well (see [START_REF] Hobert | On the applicability of regenerative simulation in Markov chain Monte Carlo[END_REF])

σ 2 n (f ) = ln-1 j=1 (f (B j ) -µ n (f )L(B j )) 2 ln-1 j=1 L(B j ) . (4)
Precisely, they proved that σ 2 n (f ) is a strongly consistent and asymptotically normal estimator of the limiting variance

σ 2 f = µ(A)E A [ 1 i τ A f (X i )], with a bias of order O(1/n) as n → ∞.
When implementing the MBB, the choice of the standardization, the bias it induces and the definition of its Bootstrap counterpart are key points to obtain the second order validity of the method. As shown in [START_REF] Bertail | Regenerative Block Bootstrap for Markov Chains[END_REF], the standardization (4), which is specifically tailored for the regenerative setting, does not weaken the performance of the RBB (see also section 3 below), while the standardization of the MBB distribution in the strong mixing case is the main barrier to achieve good performance (as shown by [START_REF] Götze | Second order correctness of the blockwise bootstrap for stationary observations[END_REF]). In most practical situations (except for the very special case of m -dependence), positive moving-block based estimates of the asymptotic variance with such good properties are not available.

Technical assumptions

The assumptions required by the validity framework of the RBB method (see [START_REF] Bertail | Regenerative Block Bootstrap for Markov Chains[END_REF] for further details) are of the following type.

• Regularity conditions: there exists κ ≥ 1 such that

E A [τ κ A ] < ∞ and E ν [τ κ A ] < ∞. • Block-moment conditions: there exists κ ≥ 1 such that E A [( τ S i=1 |f (X i )|) κ ] < ∞ and E ν [( τ S i=1 |f (X i )|) κ ] < ∞.
Using well-known results (see Chapt. 11 in [START_REF] Meyn | Markov chains and stochastic stability[END_REF] and the references therein), such assumptions may be checked in practice by establishing drift criteria of Lyapounov's type for the chain. When considering Edgeworth expansion or second order results, the "block-Cramer conditions" of the following type are additionally required.

• Block-Cramer condition:

lim sup t→∞ | E A [exp(it τ A i=1 f (X i ))] |< 1.
Such conditions are checked on various models in section 4.3 of [START_REF] Bertail | Regeneration-based Statistics for Markov Chains[END_REF] (see also [START_REF] Hobert | On the applicability of regenerative simulation in Markov chain Monte Carlo[END_REF]), including the examples considered in section 4.

Regenerative extension

Now we recall the splitting technique introduced in [START_REF] Nummelin | A splitting technique for Harris recurrent chains[END_REF]. This theoretical construction aims at extending in some sense the probabilistic structure of a general Harris chain, so as to artificially build a regeneration set. It is based on the following notion. A set S ∈ E is small for X if there exist m ∈ N * , a probability measure Φ supported by S, and δ > 0 such that

∀x ∈ S, ∀A ∈ E, Π m (x, A) ≥ δΦ(A), (5)
where Π m denotes the m-th iterate of Π. Roughly speaking, the small sets are the ones on which an iterate of the transition probability is uniformly bounded below. When (5) holds, we shall say that X satisfies the minorization condition M(m, S, δ, Φ). Small sets do exist for irreducible chains1 , a fortiori for Harris chains (any accessible set actually contains small sets, see [START_REF] Jain | Contributions to Doeblin's theory of Markov processes[END_REF]). Suppose that X fulfills M = M(m, S, δ, Φ) for some accessible set S. Take m = 1, even if it entails to replace X by the chain ((X nm , ..., X n(m+1)-1 )) n∈N . The regenerative chain onto which the initial chain X is embedded is constructed by expanding the sample space, so as to define a specific sequence (Y n ) n∈N of independent Bernoulli r.v.'s with parameter δ. This joint distribution P ν,M is obtained by randomizing the transition Π each time the chain X hits S (this happens a.s. since X is Harris). If X n ∈ S and

• if Y n = 1 (which happens with probability δ ∈ ]0, 1[), then draw X n+1 according to Φ, • if Y n = 0, (which happens with probability 1-δ), then draw X n+1 according to (1δ) -1 (Π(X n , .) -δΦ(.)).

For obtaining an insight into this construction, observe simply that, if condition [START_REF] Bertail | Note on the regeneration-based bootstrap for atomic Markov chains[END_REF] holds with m = 1, when X n ∈ S, one may write the distribution of X n+1 conditioned on X n as the following mixture

Π(X n , .) = (1 -δ){(1 -δ) -1 (Π(X n , .) -δΦ(.))} + δΦ(.), (6)
which second component is independent from X n . The bivariate Markov chain X M = ((X n , Y n )) n∈N constructed this way is called the split chain. The key point lies in the fact that S × {1} is then an atom for the split chain X M , the latter inheriting all the communication and stochastic stability properties from X. In particular the blocks constructed from the consecutive times when X M visits S × {1} are independent (if X satisfies M = M(m, S, δ, Φ) for m > 1, the resulting blocks are 1-dependent only, a form of dependence that can also be easily handled). Using this construction, one may enlarge the range of applications of the regenerative method so as to extend all of the results established for atomic chains to general Harris chains. We omit the subscript M in what follows and abusively denote by P ν the extensions of the underlying probability we consider.

On approximating the regenerative extension

Here we assume further that the conditional distributions {Π(x, dy)} x∈E and the initial distribution ν are dominated by a σ-finite measure λ of reference, so that ν(dy) = f (y)λ(dy) and Π(x, dy) = p(x, y)λ(dy) for all x ∈ E. For simplicity's sake, we suppose that condition M is fulfilled with m = 1. This entails that Φ is absolutely continuous with respect to λ too, and that p(x, y) ≥ δφ(y), λ(dy)a.s.

for any x ∈ S, with Φ(dy) = φ(y)dy.

If we were able to generate practically binary random variables Y 1 , ..., Y n , so that X M (n) = ((X 1 , Y 1 ), ..., (X n , Y n )) be a realization of the split chain X M described above, then we could divide the sample path X (n) = (X 1 , ..., X n ) into regeneration blocks, as in §2.1. Therefore, knowledge of Π is required to draw Y 1 , ..., Y n this way. As a matter of fact, the distribution

L (n) (p, S, δ, φ, x (n+1) ) of Y (n) = (Y 1 , ..., Y n ) conditioned on X (n+1) = (x 1 , ..., x n+1
) is the tensor product of Bernoulli distributions given by:

∀β (n) = (β 1 , ..., β n ) ∈ {0, 1} n , ∀x (n+1) = (x 1 , ..., x n+1 ) ∈ E n+1 , P ν (Y (n) = β (n) | X (n+1) = x (n+1) ) = n i=1 P ν (Y i = β i | X i = x i , X i+1 = x i+1 )(8)
with for 1 i n:

if x i / ∈ S, P ν (Y i = β i | X i = x i , X i+1 = x i+1 ) = Ber δ (β i ) , (9)
and if

x i ∈ S, P ν (Y i = 1 | X i = x i , X i+1 = x i+1 ) = δφ(x i+1 )/p(x i , x i+1 ), P ν (Y i = 0 | X i = x i , X i+1 = x i+1 ) = 1 -δφ(x i+1 )/p(x i , x i+1 ). (10)
In short, given X (n+1) , the Y i 's are Bernoulli r.v.'s with parameter δ, unless X has hit the small set S at time i: in this case Y i is drawn from the Bernoulli distribution with parameter δφ(X i+1 )/p(X i , X i+1 ). Our proposition for constructing data blocks relies in approximating this construction by computing first an estimate p n (x, y) of the transition density p(x, y) from data X 1 , ..., X n+1 , and then drawing a random vector ( Y 1 , ..., Y n ) from the distribution

L (n) (p n , S, δ, φ, X (n+1)
), obtained by simply plugging p n in ( 10) and ( 10) (the estimate p n (x, y) may be picked such that p n (x, y) ≥ δφ(y), λ(dy) a.s., and

p n (X i , X i+1 ) > 0, 1 i n).
From a practical viewpoint, it actually suffices to draw the Y i 's only at times i when the chain hits the small set S, Y i indicating then whether the trajectory should be divided at time point i or not (see Fig. 1 for instance). This way, on gets the approximate regeneration blocks B 1 , ..., B ln-

1 with l n = 1 k n I{X k ∈ S, Y k = 1}.
Of course, knowledge of parameters (S, δ, φ) of condition ( 7) is required for this construction. In § 3.2, we shall discuss a practical method for selecting those parameters.

The question of accuracy of this approximation has been addressed in [START_REF] Bertail | Regenerative Block Bootstrap for Markov Chains[END_REF]. Precisely, they established a bound for the deviation between the distribution of ((X i , Y i )) 1 i n and the one of the ((X i , Y i )) 1 i n in the sense of the Mallows distance, which essentially depends on the rate of the uniform convergence of p n (x, y) to p(x, y) over S × S.

The (A)RBB methodology

Now that necessary background material has been reviewed, we turn to describe the following block-resampling method and discuss practical issues encountered for implementing the latter.

The (A)RBB algorithm

Suppose that the finite sample path has been divided into true or approximate regeneration blocks B 1 , ..., B ln-1 . The (approximate) regenerative blockbootstrap algorithm for estimating the sample distribution of some statistic T n = T (B 1 , ..., B ln-1 ) estimating some parameter θ with standardization

σ n = σ(B 1 , ..., B ln-1 ), namely H(x) = P(σ -1 n (T n -θ) x), (11)
is performed in three steps as follows.

(1) Draw sequentially bootstrap data blocks B * 1 , ..., B * k independently from the empirical distribution F n = (l n -1) -1 ln-1 j=1 δ B j of the blocks B 1 , ..., B ln-1 , conditioned on X (n) until the length of the bootstrap data series

l * (k) = k j=1 l(B * j ) is larger than n. Let l * n = inf{k 1, l * (k) > n}.
(2) From the bootstrap data blocks generated at step 1, reconstruct a pseudotrajectory by binding the blocks together, getting the reconstructed (A)RBB sample path

X * (n) = (B * 1 , ..., B * l * n -1 ). ( 12
)
Then compute the (A)RBB statistic and the (A)RBB standardization

T * n = T (X * (n) )andσ * (n) n = σ(X * (n) ). ( 13
)
(3) The (A)RBB distribution is then given by

H (A)RBB (x) = P * (σ * -1 n (T * n -T n ) x | X (n+1) ), (14)
denoting by P * (. | X (n+1) ) the conditional probability given X (n+1) .

A Monte-Carlo approximation to H ARBB (x) may be straightforwardly computed by repeating independently N times the procedure above. Based on Edgeworth expansions proved in [START_REF] Bertail | Edgeworth expansions for suitably normalized sample mean statistics of atomic Markov chains[END_REF], one may show that in the regenerative positive recurrent case, the RBB method inherits the accuracy of the standard i.i.d. bootstrap (see [START_REF] Hall | The Bootstrap and Edgeworth Expansion[END_REF]) up to O Pν (n -1 ) for additive functionals of type n -1 1 k n f (X k ) under weak conditions (see Theorem 3.3 in [START_REF] Bertail | Regenerative Block Bootstrap for Markov Chains[END_REF] for further details). In [START_REF] Bertail | Regeneration-based Statistics for Markov Chains[END_REF] asymptotic validity of the RBB has also been established for more general functionals, including U or V statistics based on regeneration blocks. In the general Harris recurrent case, the ARBB method for bootstrapping Markov chains simply relies in applying the RBB procedure to the data ((X 1 , Y 1 ), ..., (X n , Y n )) as if they were exactly drawn from the atomic chain X M . However, as shown in [START_REF] Bertail | Regenerative Block Bootstrap for Markov Chains[END_REF], even if it requires to use a consistent estimate of the "nuisance parameter" p and the corresponding approximate blocks it induces, this bootstrap method still remains asymptotically valid.

In [START_REF] Bertail | Edgeworth expansions for suitably normalized sample mean statistics of atomic Markov chains[END_REF] (see Prop. 3.1) it is shown that in the nonstationary case (i.e. when the initial law ν differs from µ), the first data block B 0 induces a significant bias, of order O(n -1 ), which cannot be estimated from a single realization X (n) of the chain starting from ν. Practitioners are thus recommended not to use estimators based on the whole trajectory. This fact is known as the burnin (time) problem in the bayesian literature on MCMC algorithms, related to the time from which the 1-dimensional marginal of a (simulated) chain is close enough to the limit distribution µ. When the statistic is built using regeneration blocks only, the first (non-regenerative) block has no impact on the second order properties of the ARBB estimate . However, from a practical viewpoint, it may happen that the size of the first block is large compared to the size n of the whole trajectory (for instance in the case where the expected return time to the (pseudo-)atom when starting with ν is large), the effective sample size for constructing the data blocks and the corresponding statistic is then dramatically reduced. In such a case, for mimicking the distribution of the original statistic, it is preferable, heuristically speaking, to draw sequentially the bootstrap blocks B * 1 , ..., B * k independently from the empirical distribution F n , until l(B 0 ) + k j=1 l(B * j ) is larger than n, taking practically into account the size l(B 0 ) this way (although it does not play any role in the asymptotic behavior, since l(B 0 )/n = O Pν (n -1 ) as n → ∞).

Tuning parameters

In the general (non-regenerative) case, the procedure above may be very sensitive to the choice of the minorization condition parameters (S, δ, Φ). It is essential to pick the latter in a data-driven fashion, so that enough blocks may be obtained for computing meaningful statistics, their accuracy increasing as the mean number of pseudo-regenerative blocks, that is

N n (S) = E ν [ n i=1 I{X i ∈ S, Y i = 1} |X (n+1) ], (15)
for a given realization of the trajectory. Therefore, this is somehow determined by the size of the small set chosen. More precisely, it depends on how often the chain visits the latter in a finite length path) and how sharp is the lower bound in the minorization condition. The trade-off is as follows: as the size of the small set S used for the data blocks construction increases, the number of points of the trajectory that are candidates for determining a 'cut' in the trajectory naturally increases, but, since the uniform lower bound for p(x, y) over S 2 then decreases, the probability of drawing Y i = 1 also decreases (see expression [START_REF] Bose | Edgeworth correction by bootstrap in autoregressions[END_REF]). Thus one may heuristically expect better numerical results for the ARBB, when one implements it by choosing S so as to maximize the expected number of data blocks given the trajectory, namely N n (S) -1.

In the case when the chain takes real values and in lack of any prior information about its structure, a possible data-driven method for selecting the tuning parameters could be as follows. Let S be a collection of borelian sets S (typically compact intervals) and let U S (dy) = φ S (y).λ(dy) denote the uniform distribution on S, where φ S (y) = I{y ∈ S}/λ(S) and λ is the Lebesgue measure on R. For any S ∈ S, we clearly have p(x, y) ≥ δ(S)φ S (y) for all x, y in S, with δ(S) = λ(S). inf (x,y)∈S 2 p(x, y). When δ(S) > 0, the theoretical criterion [START_REF] Bühlmann | Block length selection in the bootstrap for time series[END_REF], that one would ideally seek to maximize over S, can be written as follows

N n (S) = inf (x,y)∈S 2 p(x, y) × n i=1 I{(X i , X i+1 ) ∈ S 2 } p(X i , X i+1 ) . ( 16
)
Observing that N n (S)/n converges P ν -a.s. to λ(S)µ(S) as n → ∞, an alternative criterion to maximize, independent from the data and asymptotically equivalent to N n (S), is given by

N n (S) = n inf (x,y)∈S 2 p(x, y) λ(S)µ(S), (17)
One gets an empirical counterpart of these quantities by replacing the unknown transition density p(x, y) by an estimate p n (x, y) in expression ( 16) or [START_REF] Clémençon | Adaptive estimation of the transition density of a regular Markov chain[END_REF], and µ(S) by the empirical estimator µ n (S) = n -1 1 i n I{X i ∈ S}. Actually from a bootstrap viewpoint, the conditional criterion ( 16) is more pertinent because the rate of convergence of the ARBB distribution is directly related to the effective number of observed regeneration times conditionally to the trajectory. Note furthermore that many nonparametric estimators of the transition density of Harris recurrent chains have been proposed in the literature, among which the standard Nadaraya-Watson estimator

p n (x, y) = n i=1 K(h -1 (x -X i ))K(h -1 (y -X i+1 )) n i=1 K(h -1 (x -X i )) , (18)
computed from a Parzen-Rosenblatt kernel K(x) and a bandwidth h > 0.

In the positive recurrent case, their estimation rates have been established under various smoothness assumptions on the density of the joint distribution µ(dx)Π(x, dy) and the one of µ(dx) (see [START_REF] Athreya | Kernel estimation for real-valued Markov chains[END_REF] or [START_REF] Clémençon | Adaptive estimation of the transition density of a regular Markov chain[END_REF] and the references therein for instance).

Once p n (x, y) is computed, calculate its minimum over sets S of the class S and maximize then the practical empirical criterion over S:

S * = arg max S∈S N n (S) (19)
with

N n (S) = inf (x,y)∈S 2 p n (x, y) × n i=1 I{(X i , X i+1 ) ∈ S 2 } p n (X i , X i+1 ) . ( 20
)
On many examples of real valued chains (see section 4 below), it is possible to check at hand that any compact interval V x 0 (ε) = [x 0ε, x 0 + ε] for a suitably chosen x 0 ∈ R and ε > 0 small enough, is small, choosing φ as the density φ Vx 0 (ε) of the uniform distribution on V x 0 (ε). For practical purpose, one may perform the optimization over ε > 0, while x 0 is kept fixed (see [START_REF] Bertail | Regenerative Block Bootstrap for Markov Chains[END_REF], [START_REF] Bertail | Regeneration-based Statistics for Markov Chains[END_REF]). But both x 0 and ε may be considered as tuning parameters: searching for (x 0 , ε) over a pre-selected grid G = {(x 0 (k), ε(l)), 1 k K, 1 l L} such that inf (x,y)∈Vx 0 (ε) 2 p n (x, y) > 0 for any (x 0 , ε) ∈ G could lead to the following numerically feasible selection rule. For all (x 0 , ε) ∈ G, compute the estimated expected number of approximate pseudo-regenerations:

N n (x 0 , ε) = δ n (x 0 , ε) 2ε n i=1 I{(X i , X i+1 ) ∈ V x 0 (ε) 2 } p n (X i , X i+1 ) , (21)
with

δ n (x 0 , ε) = 2ε. inf (x,y)∈Vx 0 (ε) 2 p n (x, y). Then, pick (x * 0 , ε * ) ∈ G maximiz- ing N n (x 0 , ε) over G, corresponding to the set S * = [x * 0 -ε * , x * 0 + ε * ] and the minorization constant δ * n = δ n (x * 0 , ε * ).
It remains next to construct the approximate pseudo-blocks using S * , δ * n and p n as described in § 2.3. We point out that other approaches may be considered for determining practically small sets and establishing accurate minorization conditions, which conditions do not necessarily involve uniform distributions besides. Refer for instance to [START_REF] Roberts | Quantitative bounds for convergence rates of continuous time Markov processes[END_REF] for Markov diffusion processes.

We end this paragraph by making the following remarks about the practical implementation of the ARBB method. We first emphasize that estimation in specific null recurrent cases (including AR(p) models with unit roots for instance) has been dealt with in [START_REF] Karlsen | Nonparametric estimation in null recurrent time series[END_REF], which established in particular consistency results for the Nadaraya-Watson estimator [START_REF] Datta | Regeneration-based bootstrap for Markov chains[END_REF]. And it is noteworthy that the procedures described above, the approximate Nummelin construction and the ARBB algorithm, are actually still asymptotically valid in this framework, when applied to adequate functions f . However, in the null recurrent case, the choice of the standardization may be cumbersome. Investigating the asymptotic properties of the ARBB at the first order, which corresponds to choosing

σ n = l 1/2 n and σ * n = l * 1/2 n
, may be done using the same approach as in [START_REF] Bertail | Regenerative Block Bootstrap for Markov Chains[END_REF] (the study of second order properties is currently in progress). As indicated by the results in [START_REF] Karlsen | Nonparametric estimation in null recurrent time series[END_REF], accurate estimation of the underlying transition density in the null recurrent case is naturally possible only when a very large data sample is at disposal. In β-null recurrent cases (i.e. when the distribution of the return time to the small set has power tail), [START_REF] Chen | How often does a Harris Recurrent Markov Chain recur?[END_REF] also established deterministic approximations of l n (respectively, of l * n ). To give an insight into the problems encountered in this case, we considered the case of an AR(1) model with a unit root among our simulation studies (see § 4.2): the number l n of regenerations over a trajectory of length n for the split chain being of order n 1/2 , only large sample sizes n enable us then to get enough (pseudo-) regeneration cycles for computing significant statistics with our methodology.

Secondly, a natural question arising from the practical considerations discussed above is to determine whether the use of the preliminary estimate p n and of µ n eventually for selecting S and building the pseudo-blocks affect the second order properties of the resulting ARBB distribution. This seems to be a very difficult problem, since by construction the pseudo-regeneration times and the data blocks B j they induce, all depend on the whole trajectory now, owing to the transition probability estimation step. A possible construction to avoid this theoretical problem consists in using a double splitting trick in a semiparametric sense (see [START_REF] Schick | Sample splitting with Markov Chains[END_REF]). This amounts first to construct the transition density estimator using the first m n observations say (with m n → ∞, m n /n → 0 as n → ∞), then to drop the next q n observations (typically q n << m n , q n → ∞ as n → ∞) for allowing the split chain to regenerate with overwhelming probability, and finally to build the pseudo-blocks B j from the n-m n -q n remaining observations. It is easy to understand (but technical to prove) that these blocks are then, asymptotically i.i.d conditionally to the first m n observations. One may then prove the second order validity of the procedure in both the studentized and unstudentized cases. As shown in [START_REF] Bertail | Approximate Regenerative Block-Bootstrap for Markov Chains: second-order properties[END_REF], this splitting trick entails some loss in the rate of the ARBB distribution, but the latter remains anyway faster than the best rate the MBB may achieve. However, one may argue, as in the semiparametric case, that such a modification of the initial procedure is essentially motivated by our limitations in the analysis of asymptotic properties of the estimators. From our own practical experience, this construction generally deteriorates the finite sample performance of the initial algorithm and estimating p(x, y) from the whole trajectory leads to better numerical results.

Simulation studies

We now give two examples, with a view to illustrate the scope of applications of our methodology. The first example presents a regenerative Markov chain described and studied at greater length in [START_REF] Harrison | The stationary distribution and first exit probabilities of a storage process with general release rule[END_REF] (see also [START_REF] Brockwell | Storage processes with general release rules and additive inputs[END_REF] and [START_REF] Browne | Work-modulated queues with applications to storage processes[END_REF]) for modeling storage systems, regenerative chains being widely used in operations research. We point out that our method also applies to the framework of MCMC (Monte-Carlo Markov Chain) in order to control estimates based on regenerative MCMC trajectories (see [START_REF] Mykland | Regeneration in Markov chain samplers[END_REF]). In consideration of the recent emphasis on nonlinear models in the time series literature, our second example shows to what extent the ARBB method may apply to a general nonlinear AR model. Further, we point out that the principles exposed in this paper are by no means restricted to the markovian setting, but may apply to any process for which a regenerative extension can be constructed and simulated from the data available (see Chapt. 10 in [START_REF] Thorisson | Coupling, Stationarity and Regeneration[END_REF]).

Example 1 : content-dependent storage systems

We consider a general model for storage, evolving through a sequence of input times (T n ) n∈N (with T 0 = 0 by convention), at which the storage system is replenished. Let S n be the amount of input into the storage system at the n th input time T n and C t be the amount of contents of the storage system at time t. When possible, there is withdrawal from the storage system between these input times at the constant rate r and the amount of stored contents that drops in a time period [T, T + ∆T ] since the latter input time is equal to C T -C T +∆T = r∆T , and when the amount of contents reaches zero, it continues to take the value zero until it is replenished at the next input time. If X n denotes the amount of contents immediately before the input time T n (i.e. X n = C Tn -S n ), we have for all n ∈ N,

X n+1 = (X n + S n -r∆T n+1 ) + , (22)
with (x) + = sup (x, 0) , X 0 = 0 by convention and ∆T n = T n -T n-1 for all n ≥ 1. Let K(x, ds) be a transition probability on R + . Assume that, conditionally to X 1 , ..., X n , the amounts of input S 1 , ..., S n are independent from each other and independent from the inter-arrival times ∆T 1 , ..., ∆T n and that the distribution of S i is given by K(X i , .), for 0 i n. Under the further assumption that (∆T n ) n 1 is an i.i.d. sequence with common distribution G, independent from X = (X n ) n∈N , the storage process X is a Markov chain with transition probability Π given by Π(

X n , {0}) = Γ(X n , [X n , ∞[), Π(X n , ]x, ∞[) = Γ(X n , ]-∞, X n -x[
) for all x > 0, where the transition probability Γ is given by the convolution product Γ(

x, ]-∞, y[) = ∞ t=0 ∞ z=0 G(dt)K(x, dz)I{rt -z < y}.
One may check that the chain Π is δ 0 -irreducible as soon as K(x, .) has infinite tail for all x 0. In this case, {0} is an accessible atom for X and it can be shown that it is positive recurrent if and only if there exists b > 0 and a test function

V : R + → [0, ∞] such that V (0) < ∞ and for all x 0 : Π(x, dy)V (y) -V (x) -1 + bI{x = 0}. ( 23
)
The times at which the storage process X reaches the value 0 are thus regeneration times, and allow to define regeneration blocks dividing the sample path, as shown in Figure 1. Figure 2 below shows a reconstructed RBB data series, generated by a sequential sampling of the regeneration blocks (as described in § 3.1), on which RBB statistics may be based.

Simulation results

We simulated two trajectories of respective length n = 100 and n = 200 drawn from this Markov chain with r = 1, K(x, dy) = Exp 3 (dy) and G(dy) = Exp 1 (dy), denoting by Exp λ (dy) the exponential distribution with mean 1/λ > 0, which is a standard M/M/1 model (see [START_REF] Asmussen | Applied Probabilities and Queues[END_REF] for instance). In Fig. 3 below,a Monte-Carlo estimate of the true distribution of the sample mean standardized by its estimated standard error (as defined in ( 4)) computed with 10000 simulated trajectories is compared to the RBB distribution (in both cases, Monte-Carlo approximations of RBB estimates are computed from B = 2000 repetitions of the RBB procedure) and to the gaussian approximation. Note also that in the ideal case where one a priori knows the exact form of the markovian data generating process, one may naturally construct a bootstrap distribution in a parametric fashion by estimating first the parameters of the M/M/1 model, and then simulating bootstrap trajectories based on these estimates. Such an ideal procedure naturally performs very well in practice. In our simulation study, the resulting distribution estimate was actually so close to the true distribution that one could not distinguish one from the other in the plot. Of course, in most applications practitioners have generally no knowledge of the exact form of the underlying Markov model, since this is often one of the major goals of statistical inference.

With the aim of constructing accurate confidence intervals, Table 1 compares the quantile of order γ of the true distribution, the one of the gaussian approximation (both estimated with 10000 simulated trajectories) and the mean of the quantile of order γ of the RBB distribution over 100 repetitions of the RBB procedure in the tail regions.

The left tail is clearly very well estimated, whereas the right tail gives a better approximation than the asymptotic distribution. The gain in term of coverage accuracy is quite enormous in comparison to the asymptotic distribution. For instance at the level 95%, for n = 200, the asymptotic distribution yields a bilateral coverage interval of level 71% only, whereas the RBB distribution yields a level of 92% in our simulation.

Example 2 : General autoregressive models

Consider now the general heteroscedastic autoregressive model

X n+1 = m(X n ) + σ(X n )ε n+1 , n ∈ N, (24)
where m : R → R and σ : R → R * + are measurable functions, (ε n ) n∈N is a i.i.d. sequence of r.v.'s drawn from g(x)dx such that, for all n ∈ N, ε n+1 is independent from the X k 's, k n with E[ε n+1 ] = 0 and var(ε n+1 ) = 1. See [START_REF] Franke | Bootstrap of kernel smoothing in nonlinear time series[END_REF] for some proposals for bootstrapping such models. The transition density of the chain is given by p(x, y) = g((ym(x))/σ(x)), (x, y) ∈ R 2 . Assume further that g, m and σ are continuous functions and there exists x 0 ∈ R such that p(x 0 , x 0 ) > 0. Then, the transition density is uniformly bounded from below over some neighborhood

V x 0 (ε) 2 = [x 0 -ε, x 0 + ε] 2 of (x 0 , x 0 ) in R 2 : there exists δ = δ(ε) ∈]0, 1[ such that, inf (x,y)∈V 2 x 0 p(x, y) δ(2ε) -1 . (25)
Any compact interval V x 0 (ε) is thus a small set for the chain X, which satisfies the minorization condition M(1,

V x 0 (ε), δ, U Vx 0 (ε) )
, where U Vx 0 (ε) denotes the uniform distribution on V x 0 (ε). Hence, in the case when one knows x 0 , ε and δ such that ( 5) holds (this simply amounts to know a uniform lower bound estimate for the probability to return to V x 0 (ε) in one step), one may effectively apply the ARBB methodology to X. In the following, we use the practical criterion N n (x 0 , ε) with x 0 = 0. The choice x 0 = 0 is simply motivated by observing that our temporal simulated data fluctuate around 0. Actually, to our own practical experience, optimizing over x 0 does not really improve the performance of the procedure in this case.

In what follows, we shall compare the performance of the ARBB to the one of some reference competitors for bootstrapping time series. In all our simulations the Markov bootstrap (consisting in generating a Markov chains with an estimated transition probability) has performed always worse than all the other methods (due to the difficulty of estimating accurately the transition probability on the whole real line). We do not present the results for this method to alleviate the graphics and tables.

The sieve bootstrap is specifically tailored for linear time series (see [START_REF] Bühlmann | Sieve Bootstrap for time series[END_REF], [START_REF] Bühlmann | Bootstrap for time series[END_REF]).

The main idea consists in fitting an AR(p) model (eventually with p unknown depending on ) first and then applying a residual based resampling method. The fact that it fully exploits the underlying linear structure explains why it performs very well in this framework. When simulating linear time series, we use it as a benchmark for evaluating the pertinence of the ARBB distribution. Recall also that this method requires a preliminary estimation of the order q of the sieve : for this purpose we choose an AIC criterion of the type AIC(q) = nlog( M SE) + 2q in the sequel. In the linear AR(q) model below, this information criterion enables us to pick the right order of the model. And the resulting sieve bootstrap behaves like a parametric bootstrap method in these cases (see [START_REF] Bose | Edgeworth correction by bootstrap in autoregressions[END_REF]), leading to very good numerical results, as soon as the roots of the AR(q) model are far from the unit circle. In contradistinction, we actually experienced problems in our simulations, when dealing with an AR(1) model with a root close to 1: in such cases, it may happen with high probability that one gets an estimate of the root larger than one, yielding to explosive bootstrap trajectories.

We also compared the ARBB method to the usual MBB. The difficulty for applying the latter method essentially relies in the choice of the block size for estimating the variance and in the choice of the block size for the resampling procedure. As there is actually no reason for these two sizes to be equal, they should be picked separately and the estimator of the variance should be correctly unbiased (see [START_REF] Götze | Second order correctness of the blockwise bootstrap for stationary observations[END_REF]). To our knowledge, the problem of simultaneously calibrating these two quantities has not been treated yet and leads to extremely volatile results. For comparing directly the MBB distribution to the true studentized distribution [START_REF] Brockwell | Storage processes with general release rules and additive inputs[END_REF], we have chosen here to standardize all the distributions by the estimator (4), so as to avoid a deteriorating preliminary variance estimation step. The MBB distribution is also correctly centered (at the bootstrap mean). The block size for the MBB is chosen according to the method of [START_REF] Hall | On blocking rules for the bootstrap with dependent data[END_REF]. It consists in estimating first the MSE of the MBB distribution corresponding to blocks of size l with a subsampling technique for various size values l and then picking the size corresponding to a minimum MSE estimate. This unfortunately requires to select a subsampling size and a plausible pilot size, which are in their turn also difficult to calibrate (see the discussion in Section 7.3 of [START_REF] Lahiri | Resampling methods for dependent Data[END_REF]): here we have chosen n 1/4 as pilot size and b n = n 10/21 as subsampling size (which is close to n 1/2 in our simulations and satisfies the conditions needed for the MBB to be asymptotically valid). When standardized this way, the MBB has performed quite well in most simulations, except notably when data exhibit significant nonlinear features and/or nonstationarity. The reason of this misbehavior arises from the fact that, for some drawing of the fixed size blocks, the jumps between the blocks were so important, that the reconstructed series could not be splitted according to our randomized procedure leading to an invalid estimator of the variance. In these case (too few regenerations), we have eliminated the corresponding MBB simulation. Thus the MBB considered here can be considered as a MBB with a Markovian control ensuring that the MBB reconstructed series has some regeneration properties. Such procedure clearly improved the resulting estimated distributions.

Simulation results

Here are empirical evidences for three specific autoregressive models.

The AR(1) model :

X i+1 = αX i + ε i+1 , i ∈ N, (26)
with i.i.d. ǫ i ∼ N (0, 1), α = 0.8, X 0 = 0 and for a trajectory of length n = 200.

The AR(1) model with ARCH(1) residuals called AR-ARCH model :

X i+1 = αX i + (1 + βX 2 i ) 1/2 ε i+1, i ∈ N, (27)
with i.i.d. ε i ∼ N (0, 1), α = 0.6, β = 0.35, X 0 = 0 and for a trajectory of length n = 200.

The so called ExpAR(1) model

X i+1 = (α 1 + α 2 e -|X i | 2 )X i+1 + ε i+1, i ∈ N, (28)
with i.i.d. ε i ∼ N (0, 1), α 1 = 0.6, α 2 = 0.1, X 0 = 0 and for a trajectory of length n = 200. Such a chain is recurrent positive under the sole assumption that |α 1 | < 1, see [START_REF] Tjøstheim | Non Linear Time series[END_REF]. This highly nonlinear model behaves like a threshold model: when the chain takes large values, this is almost an AR(1) model with coefficient α 1 , whereas for small values, it behaves as an AR(1) model with a larger autoregressive coefficient

α 1 + α 2 .
Here the true distribution of the sample mean is estimated with 10000 simulations. And for a given trajectory, the ARBB distribution is approximated with B = 1000 resamplings of the pseudo-blocks. In a previous simulation work, we experienced that the ARBB distribution obtained may strongly fluctuate, depending on the randomization steps (see §2.3). For a given trajectory, this problem may be avoided by repeating the ARBB procedure several times (50 times in our simulations) and averaging the resulting ARBB distribution estimates. According to our experiments, only a small number of repetitions (leading to different ways of dividing the same trajectory) suffices for smoothing the ARBB distribution.

For the ARBB, the sieve and the MBB methods, the whole procedure has been repeated 1000 times. The small set is selected by maximizing over ε > 0 the empirical criterion N n (0, ε) described above. The main steps of the procedure are summarized in the graph panels shown below.

The first figure in Graph panel 1 shows the Nadaraya-Watson (NW) estimator ( 18), the second one represents N n (0, ε) as ε grows (as well as the smoother empirical criterion [START_REF] Clémençon | Adaptive estimation of the transition density of a regular Markov chain[END_REF], see the dotted line). It clearly allows to identify an optimal value for the size of the small set. In the case of the AR model for instance, this selection rule leads to pick in mean ε = 0.83 and δ = 0.123. Our empirical criterion tends to overestimate very slightly the size of the "optimal" small set (a phenomenon that we have noticed on several occasions in our simulations). The level sets of the NW estimator, the data points (X i , X i+1 ) and the estimated small set are represented in the next graphic. This also shows that the small set chosen may be not that "small" if the transition density is flat around (x 0 , x 0 ) = (0, 0) (in some cases it may be thus preferable to choose x 0 = 0 so as to be in this situation). In the second line of the panel, the figure on the left hand side represents a sample path of the chain and indicates the pseudo-regenerative blocks obtained by applying the randomization rule with Ber(1δ(2ε) -1 / p n (X i , X i+1 )) at times i when (X i , X i+1 ) ∈ V 0 (ε) 2 . The next figure shows how binded blocks form a typical ARBB trajectory. It is noteworthy that such a trajectory presents less artificial "jumps" than a trajectory reconstructed from a classical MBB procedure: by construction, blocks are joined end to end at values belonging to the small set. For comparison purpose, the figure on the right hand side displays a typical realization of a MBB trajectory. Finally, on the last line of the panel, the true distribution (green), the ARBB distribution (black), the sieve bootstrap distribution (gray), the MBB distribution (red dotted line) and the asymptotic gaussian distribution (blue dotted line) are compared.

And the last figure shows the QQ-plots α ∈

[0, 1] → G n (H -1 (α))
, where H is the true distribution and G n denotes one of the approximations: this enables us to discriminate between the various approximations in a sharper fashion, especially in the tail regions.

These results clearly indicate that both the sieve and MBB methods perform very well for linear time series. In this case, the ARBB distribution tends to have larger tails. However, when considering nonlinear models, the advantage of the ARBB method over its rivals plainly come into sight: for moderate sample sizes n, the sieve bootstrap tends to choose a too large value q n for the lag order of the approximate sieve AR( q n ). This problem is less serious for larger sample sizes, as shown in Graph panel 4 (with n = 500). In these situations, the MBB may behave very poorly especially when the non-linearity and the non-stationarity is important: we conjecture that it could be possibly improved by investigating further how to tune optimally the block size, especially for standardized distributions.

Pictures in Graph panels 3 and 4 speak volumes: for both nonlinear models, the true distribution is accurately approximated by the ARBB distribution. Note nevertheless the difference in the size of the "optimal small set" and in the number of pseudo-regenerations between these models. We point out that, though remarkable when compared to the gaussian approximation, the gain in accuracy obtained by applying the ARBB methodology to the EXP-AR model is higher than the one obtained for the AR-ARCH type model. As may be confirmed by other simulations, the ARBB method provides less accurate results for a given (moderate) sample size, as one gets closer to a unit root model (i.e. as α tends to 1): one may get an insight into this phenomenon by simply noticing that the rate of the number of regenerations (respectively, of the number of visits to the small set) then drastically decreases.

Concluding remarks

We finally summarize our empirical findings. We first point out that, in the linear case when roots are much less than 1 in amplitude, the sieve bootstrap clearly surpasses its competitors. But it is noteworthy that both the ARBB and the MBB also provides very good numerical results in this case.

Besides, all these methods seem to break down from a practical viewpoint for an AR(1) model with an autoregressive coefficient α tending to 1 and with a fixed (moderate) sample size: in such a case, too few pseudo-regeneration blocks may be constructed for the ARBB methodology to be practically performant (although it is asymptotically valid). In this respect, the graph of the estimated number of pseudo-regenerations (see Graph panels 1-4) provides a crucial help for diagnosing the success or the failure of the ARBB method.

It is also remarkable that the sieve bootstrap can lead to very bad results in this case, due to the fact that the estimated AR model may have a root larger than 1 (generating then explosive sieve bootstrap trajectories). This strongly advocates the use of preliminary tests or constrained estimation procedures (ensuring that the resulting reconstructed series is asymptotically stationary).

And as may be reported from our simulation results, the advantage of the ARBB over the sieve bootstrap, the MBB and the asymptotic distributions, clearly appears when dealing with nonlinear models even if in some case the MBB can still give some good approximation (see the AR-ARCH(1) case, Graph-panel 2). Even if the lag is chosen very large (in mean 85 for the AR-ARCH(1) model and 21 for the EXP-AR model), the linear sieve method is unable to capture the non-linearities and performs very badly for moderate sample sizes. The MBB also performs poorly in some nonlinear setting for moderate sample sizes, whereas the ARBB provides very accurate approximations of the tail distributions in these examples. It should be mentionned that using a moving-block estimator of the variance leads to even worse results.

In any case, it is recommended to use all the available methods and to compare the results. The ARBB being much more robust it can be used to check whether the other methods are trustable for the data at hand.

As pointed out by a referee, we restricted our study to the case of Markov models of order 1 (see our framework in section 2). However, by vectorizing, any 1-dimensional Markov model of order p classically boils down to a p dimensional Markov model of order 1. Our theoretical work thus applies in this context. However, statistical problems related to the curse of dimensionality may appear, when considering models of large orders (arising in the choice of the small set or the transition density estimation step), Hence, numerical experiments should be carried out for such models, in order to determine whether theoretical results are still supported by empirical evidence. This is beyond the scope of the present paper, but will certainly be the subject of further investigation. Besides, another challenging line of research could consist in determining whether the ARBB performs well when applied to long memory Markov chains. As a matter of fact, this problem is of different nature. Recall that a Markov model of order 1 may naturally have a long memory: if X is a positive recurrent chain with limiting distribution µ and possesses an atom A, the long memory property for the sequence {f (X n )} n∈N is equivalent to the condition E A [( i≤τ A {f (X i )µ(f )}) 2 ] = ∞. No theoretical result for the (A)RBB asymptotic validity in this framework has yet been established, even if one may reasonably expect that it is the case for β-null recurrent chains (which exhibit long range dependence when β < 1, see [START_REF] Karlsen | Nonparametric estimation in null recurrent time series[END_REF]).

Akcnowledgements. We thank the referees for their valuable remarks and suggestions.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Dividing the trajectory of the storage process into data blocks corresponding to the regeneration times.

Fig. 3 .

 3 Fig. 3. Distribution estimates.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Graph panel 1: AR(1) model with α = 0.8, n = 200.

Table 1 :

 1 Comparison of the tails of the true distribution (TD), RBB and gaussian distributions.

	n=	100	200	∞	n=	100	200	∞
	γ%	TD	RBB	TD	RBB	ASY	γ%	TD	RBB	TD	RBB	ASY
	1	-7.733	-7.044	-5.492	-5.588	-2.326	90	1.041	1.032	1.029	1.047	1.282
	2	-6.179	-5.734	-4.607	-4.695	-2.054	91	1.078	1.085	1.083	1.095	1.341
	3	-5.302	-5.014	-4.170	-4.165	-1.881	92	1.125	1.145	1.122	1.150	1.405
	4	-4.816	-4.473	-3.708	-3.757	-1.751	93	1.168	1.207	1.177	1.209	1.476
	5	-4.374	-4.134	-3.430	-3.477	-1.645	94	1.220	1.276	1.236	1.277	1.555
	6	-4.086	-3.853	-3.153	-3.243	-1.555	95	1.287	1.360	1.299	1.356	1.645
	7	-3.795	-3.607	-2.966	-3.045	-1.476	96	1.366	1.453	1.380	1.442	1.751
	8	-3.576	-3.374	-2.771	-2.866	-1.405	97	1.433	1.568	1.479	1.549	1.881
	9	-3.370	-3.157	-2.606	-2.709	-1.341	98	1.540	1.722	1.646	1.685	2.054
	10	-3.184	-2.950	-2.472	-2.560	-1.282	99	1.762	1.970	1.839	1.916	2.326

  Table2below gives the median of the quantiles at several orders γ of the bootstrap distributions over the 1000 replications for each of the three AR models, compared to the true and asymptotic corresponding quantiles.

	n=200	AR(1)				AR-ARCH(1)			EXP-AR(1)
	γ% TD ARBB Sieve MBB TD ARBB Sieve MBB TD ARBB Sieve MBB ASY
	1	-3.51	-3.61	-3.41	-3.42	-3.03	-3.23	-5.26	-3.16	-4.48	-5.23	-5.59	-9.61	-2.33
	2..5	-2.84	-2.78	-2.81	-2.72	-2.41	-2.61	-3.52	-2.52	-3.35	-3.87	-4.76	-6.44	-1.96
	5	-2.23	-2.13	-2.11	-2.10	-1.97	-2.14	-2.85	-2.06	-2.58	-2.79	-3.74	-5.00	-1.65
	10	-1.62	-1.57	-1.65	-1.55	-1.52	-1.59	-2.25	-1.53	-1.83	-1.98	-2.93	-3.51	-1.28
	n=200	AR(1)				AR-ARCH(1)			EXP-AR(1)
	γ% TD ARBB Sieve MBB TD ARBB Sieve MBB TD ARBB Sieve MBB ASY
	90	1.62	1.52	1.61	1.61	1.52	1.33	2.26	1.58	1.80	1.89	2.74	2.26	1.28
	95	2.21	2.08	2.19	2.14	2.01	1.74	3.04	2.07	2.58	2.68	3.89	3.07	1.65
	97.5	2,79	2.71	2.73	2.69	2.37	2.03	3.93	2.44	3.24	3.47	4.79	4.02	1.96
	99	3.46	3.73	3.86	3.48	3.11	2.62	5.97	3.22	4.37	5.36	5.92	6.25	2.33

Table 2 :

 2 Comparison of the tails of the true, ARBB and gaussian distributions for the three models

Recall that a Markov chain X with state space (E, E) and transition Π is irreducible if there exists a positive measure that dominates n≥1 Π n (x.) for all x ∈ E.

Fig. 6. Graph panel 3: EXP-AR(1) model with α 1 = 0.8 and α 2 = 0.5, n = 200.


M Ashour Pour

M Zanardini

A Bacchetti

S Zanoni

An Economic Insight into Additive Manufacturing System Implementation

Keywords: Additive Manufacturing, Economic Assessment, Implementation Guideline

With an ever growing diffusion of Additive Manufacturing (AM) system in industrial and commercial level, as well as the direct and indirect dynamics which are being introduced resulting from its inclusion as a possible production technology on companies' portfolio, the need to reconfigure production system and adapt the production strategy becomes even more relevant than before. There are several studies which have emphasized on the importance of a paradigm shift in order to exploit advantages of AM, not only considering changes within design and functionality of the product, but also concerning AM's impact on the entire value chain (re)configuration. Thus, it is of crucial importance to take into consideration that for this shift to be feasible and manageable, there is a need to include both technical and managerial aspects of manufacturing. This work proposes an economic insight in order to provide a guideline for the proper evaluation of AM system implementation.

Introduction

Considering their evolution rate, Additive Manufacturing (AM) technologies are becoming more interesting for a growing spectrum of industries, and the applications will soon impact the production activities of components and products [START_REF] Wohler | Wohlers Report[END_REF] as well. The authors have found that the scarce available literature is not aligned with companies' requirements. They are focusing mostly on technological aspects, and do not provide an holistic view to ensure the evaluation of all the benefits and limitations of AM. In order to highlight the main limitations, three categories have been identified: (i) The data provided are often not up to date (and generally refer to a particular technology or even to a particular machine); (ii) The majority of the works consider a specific technology, and do not provide a priori analysis for identifying which is the most suitable for the context of application; (iii) As a consequence of point ii, the economic analysis generally does not compare different AM technologies to evaluate which is the best one and thus, no investment decision analysis is considered.

In order to overcome the shortcomings described above, this work provides an holistic guideline; deploying an economic analysis to properly investigate the use of AM technologies within an industrial sector. The main objectives can be summarized as the following:

 To define a guideline to support manufacturing companies to understand whether AM techniques are suitable for their context and products; proposing criteria for implementing an a priori analysis, and providing a list of the (eventual) appropriate technologies and related machines,

 To provide a model to include economic implications resulting from AM application, facilitating a comparison between conventional subtractive techniques.

AM system

In order to have a better understanding of AM, it is necessary to look at it from the system point of view. This is due to the wide range of impacts that is accompanied by its implementation: extending from raw material suppliers and procurement, towards production level, distributors and even customers. Furthermore, a systematic analysis is able to characterize AM in a more comprehensive way. AM is a term applied to a technological class which consists of multiple subsets that make up the technological variations. Each of these technologies could be applied in various industries ranging from electronics to medical, automotive, armament and aerospace.

2.1

AM Characterization One of the most remarkable aspects of AM which has enhanced its position among other manufacturing techniques is the flexibility, not only enabling economical low volume production [START_REF] Mellor | Additive manufacturing: A framework for implementation[END_REF] by eliminating the need for tooling, but also providing designers with a degree of freedom that no more limits functionality in favor of feasibility of the process. This feature provides manufacturers with two remarkable opportunities regarding the design: faster time-to-market, and almost real-time design changes that happen as improvements and optimizations are made to the original design. According to [START_REF] Holmström | Rapid manufacturing in the spare parts supply chain: Alternative approaches to capacity deployment[END_REF], special features of AM would result in the following benefits:  Tool independency, since no tooling is required,  Economic production of small batches becomes feasible,  High flexibility for changing the design of the parts/products,  Optimization for product functionality would be achievable,  Customization of products which are based on individual customers' needs,  High possibilities of wastes elimination during production phase,  Possibility of having simpler and shorter supply chains.

Another study [START_REF] Bianchi | Additive manufacturing : towards a new operations management paradigm ?[END_REF] characterizes AM by highlighting its distinguishing features. High automation and part consolidation which provides the possibility to build parts as a single piece and therefore eliminate the assembly would consequently lead to a great reduction of the labor, storage, handling and logistics costs. Economies of scale are one of the most remarkable properties of mass manufacturing. Manufacturing in large volumes allows for reduction of cost per unit as a result of the fixed-cost proration. However, since AM machinery requires no setups, production in small batches becomes economically feasible and this is a direct result of "economies of one".

Economic inefficiency in large volume production, inability of processing large parts due to the chamber size limitations [START_REF] Huang | Additive manufacturing and its societal impact: a literature review[END_REF], process variability [START_REF] Bianchi | Additive manufacturing : towards a new operations management paradigm ?[END_REF] and lack of consistency among produced parts to ensure mechanical properties of the parts [START_REF] Guo | Additive manufacturing: technology, applications and research needs[END_REF], incompetency of the companies struggling with process automation and digitalization, limited range of raw materials and lack of international standardization are amongst the most important barriers towards considering AM as manufacturing method.

AM technology variations

Different approaches have been introduced to classify AM processes. One classification is based on the raw materials feed. While it classifies processes according to four types of input raw materials, another type of classification is based solely on the working principle of process [START_REF] Guo | Additive manufacturing: technology, applications and research needs[END_REF]. Amongst the others, the most comprehensive classification is presented by the AM subset of American Society for Testing and Materials (ASTM F42), categorizing AM variants into seven broad groups [START_REF] Guo | Additive manufacturing: technology, applications and research needs[END_REF].

The notable areas in which AM has been deployed with high rates of success are currently limited to medical devices, consumer goods (e.g. electronics), aerospace, automotive, jewelry, architecture and defense [START_REF] Mellor | Additive manufacturing: A framework for implementation[END_REF]. Although various studies have considered the issue of energy usage in AM machines, a unified and standard procedure to measure energy consumption is still lacking and there needs to be more data for making comparisons among conventional technologies and AM. However, there are multiple studies [START_REF] Mani | Sustainability Characterization for Additive Manufacturing[END_REF] which show when it comes to the environmental aspects and carbon footprint, AM has a positive impact. Needless to say that a majority of these researches would still pinpoint the focus of their investigations into the lack of detailed information regarding wastes, energy consumption and environmental impacts.

Impactful dynamics of AM

As it was mentioned earlier, AM is a system which is attributed by a variety of dynamics. One of these attributes which directly impacts the value chain is the supply chain management. The ability to redesign products with fewer components and the possibility of manufacturing products near customers' physical location are two valuable opportunities offered by AM [START_REF] Huang | Additive manufacturing and its societal impact: a literature review[END_REF]. This would not only reduce the need for warehousing, transportation and inventory, but it would also make the supply chain simpler by reducing time-to-market and lead-time. Design for Additive Manufacturing (DFAM) is a term which is used to emphasize on the flexibility of AM; meaning that since there are no limitations imposed by the design of the product to reduce its functionality, parts can be redesigned into single components and thus, AM's capabilities would be exploited in a more efficient way. By doing so, a reduction of the materials, energy and natural resources would take place which would eventually result in significant sustainable and economic benefits. An exciting area for AM to implement is in the spare parts supply chain. A thorough investigation [START_REF] Walter | Rapid manufacturing and its impact on supply chain management[END_REF] of spare parts supply in aircraft industry shows that rapid manufacturing (a term used for AM of individual parts/small lot sizes) can be used for low volume production of parts in a centralized location and at the place of consumption, if inventory holding and logistics costs are high in comparison with the production costs. This strategy would keep stock level down and AM capacity utilization high. In another study [START_REF] Khajavi | Additive manufacturing in the spare parts supply chain[END_REF] four scenarios were studied in two dimensions of supply chain configuration (centralized and decentralized) and AM machine technology (current and future technology). One significant outcome of the study showed that with the current maturity of AM in which machines are both capital-and labor-intensive, centralized production is more efficient, while with the evolution of technology in the future, characterized by cheaper and more automated machines, distribution of production would be a better choice for the spare parts supply chain.

Although the lack of comprehensive data to assess sustainability aspects of AM is a big impediment, some researchers have tried to identify the key concepts of AM which are relevant to sustainable manufacturing [START_REF] Mani | Sustainability Characterization for Additive Manufacturing[END_REF]. These are the same advantages that distinguish AM from conventional and traditional manufacturing processes. Considering the current legislation and regulatory laws that exist on the environmental aspects of manufacturing processes, and manufacturers' tendency towards moving to cleaner and more sustainable production, the environmental impacts of AM is part and parcel of any analytic assessment. An analytic model on the evaluation of environmental impacts in AM [START_REF] Bourhis | Sustainable manufacturing : evaluation and modeling of environmental impacts in additive manufacturing To cite this version[END_REF] which considered the whole environmental flows, shows that in order to study the global environmental impacts, not only the electricity power consumption, but also the materials, and fluid consumption need to be taken into account.

State of the art of AM adoption frameworks

One of the first contributions in this research stream comes from [START_REF] Hopkinson | Analysis of rapid manufacturing-using layer manufacturing processes for production[END_REF] that provided a model for cost estimation of AM applications. It analyzes the direct cost of production considering the machine, labor, and material costs, omitting the overhead costs, as well as the energy consumption. One of the most relevant outputs of the work is evaluation of the typical 3D printing cost profile, which is independent from quantity of the production. In the subsequent studies, the results are discussed in detail and then confuted, especially for small production batches.

The use of activity-based costing for the economic analysis of an AM alternative is provided by [START_REF] Ruffo | Cost estimation for rapid manufacturinglaser sintering production for low to medium volumes[END_REF], however the proposed model has strong limitations since it considers only one single technique. Nevertheless, [START_REF] Ruffo | Cost estimation for rapid manufacturinglaser sintering production for low to medium volumes[END_REF] confirmed previous assumptions, meaning that the more production chamber is saturated, the more the unit cost production is reduced. The models provided in the following years [START_REF] Elenora Atzeni | Economics of additive manufacturing for end-usable metal parts[END_REF] try to evaluate cost of 3D printing application in an holistic view, considering a life cycle approach. In these works, the authors encompassed also re-designing activities which are required for a full exploitation of 3D printing capabilities, while incorporating full advantages enabled by AM. The approach considered by [START_REF] Kair | Additive Manufacturing and Production of Metallic Parts in Automotive Industry: A Case Study on Technical[END_REF] is one of the most comprehensive ones. First, they identify a list of possible products that may be revisited by AM, for each of which, they then evaluate the most appropriate technology that matches the firm's requirements, and only at the end of this evaluation process, do the authors develop an economic analysis.

Considering a more consultancy-oriented approach, one has to notice Senvol (included in [START_REF] Wohler | Wohlers Report[END_REF]), the company which experts in AM machinery and applications based in the US. In the paragraph titled "Cost-Benefit Analyses for Final Production Parts", the authors explain applications of their cost evaluation model. Contrary to the previous works cited earlier (e.g. [START_REF] Hopkinson | Analysis of rapid manufacturing-using layer manufacturing processes for production[END_REF]& [START_REF] Ruffo | Cost estimation for rapid manufacturinglaser sintering production for low to medium volumes[END_REF]), and due to the inefficiencies caused by print batches, their model does not provide a constant production cost. Thus, until the printing chamber is not completely saturated, the production cost per part provided is not constant. Considering the assumption that the more the machine is saturated, the lower the final production cost per part, previous scientific works that hypothesize to fully load the printer capacity seem more attractive. This assumption is reasonable, taking into account that (due to the absence of setup costs) a given company could saturate the build chamber with other parts/products and hence produce with a fully saturated chamber.

Evaluation guideline and case study

According to [START_REF] Achillas | A methodological framework for the inclusion of modern additive manufacturing into the production portfolio of a focused factory[END_REF], one of the main requirements of companies with respect to AM is the analytic support in order to evaluate whether or not AM could be suitable for their production processes and products. In this era, academics have to propose guidelines that help "senior management to reconsider whether they will continue using current production technologies, or they could benefit by exploiting the benefits of modern AM technologies". In accordance with this statement, the authors have identified a logical path that a company which is approaching AM for the first time, could follow in order to have a comprehensive evaluation of the AM techniques (Error! Reference source not found.).

Fig. 1. Proposed AM evaluation framework

Reminding the main shortcomings described in the introduction section, the guideline aims to provide a complete assessment of this new paradigm; considering not only an economic evaluation or a context/product analysis, but also a comprehensive assessment of the alternative AM technologies that can be employed. The text provides an a priori analysis of the AM applications, considering the main features of the company context and products properties, and then proposes a technical and economic model to provide data for a quantitative assessment.

5.1

Preliminary assessment The first step that a company should take to understand whether these technologies may bring advantage to its business is a preliminary (qualitative) assessment of the context in which it operates. A framework provided by [START_REF] Conner | Making sense of 3-D printing: creating a map of additive manufacturing products and services[END_REF] encompasses three key attributes: production volume, customization, and complexity of the products.

The cited works do not provide quantitative drivers for an easy evaluation of which products are more promising for AM. Indeed, it is argued that already in this preliminary step a company should consider indicators related both to product and supply chain features for identifying which products (or range of products) could be encompassed in the following evaluation steps. For this reason, four different quantitative drivers are introduced that permit a selection of the most promising products for AM:

𝐶𝑜𝑠𝑡 𝑤𝑒𝑖𝑔ℎ𝑡 𝑖𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦 = 𝑃𝑟𝑜𝑑𝑢𝑐𝑡 𝑜𝑣𝑒𝑟𝑎𝑙𝑙 𝑐𝑜𝑠𝑡 [€] 𝑃𝑟𝑜𝑑𝑢𝑐𝑡 𝑤𝑒𝑖𝑔ℎ𝑡 [𝑘𝑔] (1) 𝐵𝑢𝑦 𝑡𝑜 𝑓𝑙𝑦 = 𝑇𝑜𝑡𝑎𝑙 𝑚𝑎𝑡𝑒𝑟𝑖𝑎𝑙 𝑐𝑜𝑛𝑠𝑢𝑚𝑝𝑡𝑖𝑜𝑛 [𝑘𝑔] 𝑃𝑟𝑜𝑑𝑢𝑐𝑡 𝑤𝑒𝑖𝑔ℎ𝑡 [𝑘𝑔] (2) 𝑀𝑜𝑢𝑙𝑑 𝑐𝑜𝑠𝑡 𝑖𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦 = 𝐶𝑜𝑠𝑡 𝑜𝑓 𝑡ℎ𝑒 𝑚𝑜𝑢𝑙𝑑 𝑎𝑙𝑙𝑜𝑐𝑎𝑡𝑒𝑑 𝑡𝑜 𝑡ℎ𝑒 𝑝𝑟𝑜𝑑𝑢𝑐𝑡 [€] 𝑃𝑟𝑜𝑑𝑢𝑐𝑡 𝑜𝑣𝑒𝑟𝑎𝑙𝑙 𝑐𝑜𝑠𝑡 [€] (3) 𝐶𝑁𝐶 𝑡𝑖𝑚𝑒 𝑖𝑛𝑡𝑒𝑛𝑠𝑖𝑡𝑦 = 𝐶𝑁𝐶 𝑡𝑖𝑚𝑒 𝑐𝑜𝑛𝑠𝑢𝑚𝑝𝑡𝑖𝑜𝑛 [ℎ] 𝑃𝑟𝑜𝑑𝑢𝑐𝑡 𝑣𝑜𝑙𝑢𝑚𝑒 [𝑑𝑚 3 ] (4)
In which, the higher the ratios, the more appropriate AM alternative would be.

Technical assessment

The next evaluation phase takes the input from selected products of the last part, aiming at the evaluation of the technological feasibility to manufacture them through AM. A more quantitative analysis is performed in order to map some relevant product features that have to be considered in this technological assessment; including dimensions, materials, physical and mechanical properties etc. The comparison of these parameters with a machines' database ensures to identify the (eventual) technologies and the related machines that are suitable for the company's products and needs. The output of this step is a list of (technology and) machines that fulfill the company's requirements, along with information about the machinery price as well as the retailers that could provide them.

Economic assessment

Having identified the list of suitable machines, it is possible to perform a preliminary evaluation of the costs incurred by the company. The developed model ensures two different types of analysis: one for evaluating whether products or components made by AM are more cost-effective than the same products or components realized through conventional subtractive techniques (e.g. injection molding or CNC machining), and the other one for evaluating which of the AM technologies that fulfill company's needs is more cost-effective, overcoming a general limitation of the literature.

According to [START_REF] Hopkinson | Analysis of rapid manufacturing-using layer manufacturing processes for production[END_REF], the provided model computes the direct costs of the AM application in terms of machine, materials and workforce. Indeed, thanks to the rigorous data collection, the costs related to the maintenance activities that [START_REF] Ruffo | Cost estimation for rapid manufacturinglaser sintering production for low to medium volumes[END_REF] took into account as indirect cost, is considered to be a direct cost.

5.2

Case study The proposed case study considers a company that has exploited AM since 2001, and reached a high level of expertise, especially in SLA [START_REF] Guo | Additive manufacturing: technology, applications and research needs[END_REF]. The company belongs to the automotive sector (specialized in the production of racing components), and operates following an Engineer-To-Order strategy in a one-of-a-kind-production context.

Initially, the SLA technology is compared with SLS [START_REF] Guo | Additive manufacturing: technology, applications and research needs[END_REF], and due to several reasons, SLA has been selected as the most appropriate choice. First of all, the company considered the liquid material (raw material for SLA) easier to manage (than the powder required for SLS). Secondly, the durability of the products (longer for SLS) was not required, since the prototypes are used no longer than one month. These requirements have been changed in the last years, and so the company needs to reevaluate the available technologies. The company's choices are limited to SLA and SLS technologies, since other AM variants produce products whose technical characteristics wouldn't match with those of the company and its customers. Thus, the company is happy with limiting its range of selection to these two AM variants.

According to the guideline described before, the company context is first considered for a preliminary analysis. After the positive qualitative results, the proper tech-nologies for the specific requirements are identified and then an economic analysis of the selected technologies (and printers) is performed.

Having already discussed with R&D department and wind tunnel managers in order to analyze the main product features as well as the company's specialty, the company is understood to be operating in a one of a kind production sector, where each product is realized upon specific customers' requirements and needs. Considering the quantitative drivers introduced above, (at least) two of them are considered relevant for AM evaluation:  For the majority of the products and prototypes realized, the buy-to-fly ratio is very high (exploiting conventional casting and molding technologies) according to the hollow structure required,  And considering the uniqueness of the products, the mold cost intensity is also very high, in accordance with the allocation of the mold cost to only one product manufactured.

Considering these two drivers, and in response to the high levels of product complexity and customization, as well as the low volume productions, the company context immediately appears to be highly suitable for AM application. Once the evaluation of the preliminary feasibility study of AM techniques is performed for the products, the next step i.e. the technical assessment is started, through which it becomes possible to identify which technologies are able to satisfy company's requirements. Taking into account products' dimensions, surface finish, mechanical properties, and the required production volume (data collected through interviews with R&D manager), the technical database is consulted to exhaust all the available options: not surprisingly, the output provides 15 printers using SLA and SLS technologies for polymer materials. These technologies ensure an appropriate level of product porosity, which is one of the main requirements for the company (even though SLA reaches better performance levels related to this parameter). The technical database provides: 6 printers out of 15 exploit SLA technology (differing only in the chamber's dimensions), while the remaining 9 printers make use of SLS technology for manufacturing parts.

According to the proposed guideline and using the collected data, an economic evaluation is then performed to compute the total production cost for all 15 printers. In the remainder of the paper the data for two analyses are discussed: one related to the SLA technology (the iPro800 printer which is the actual one adopted by the company), in order to validate the model, and one related to the SLS technology, in order to evaluate an alternative scenario (considering the most advanced printer coming from previous step, that is the Spro140 HD by 3D Systems).

In the first analysis in which the company's 3D printing cost structure (given by the interviewed personnel) is compared with the value computed by the proposed model, an initial model validation is obtained: the output approximates the actual cost given by the company with high accuracy (>95%). So, in a preliminary way, it is assumed that the model adequately represents the behavior of the real system for the project objectives. Then the model is modified in a way to estimate the SLS printing process cost structure (according to the 3 cost elements described earlier). As a result, an evaluation of the performances in the as-is and to-be scenarios is obtainable. The comparison between SLS and SLA model outputs, highlights lower overall cost for SLS, with a global saving of more than €300,000 per year (about 25% reduction). In Error! Reference source not found. a comparison of the cost structures for as-is and to-be scenarios is illustrated. Fig. 2. Economic analysis of as-is and to-be scenarios  Material cost is lower for SLS, considering that higher waste rate of this technology is balanced out with cost per kg of raw materials, which is roughly half compared with SLA.  Machine cost is higher for SLA, and constituted the majority of the gap between the overall costs of the two technologies. This is due to the stacking parameter: SLA does not permit to stack up different products on different layers, while SLS (that exploits powder material) ensures to fully saturate the printing chamber with different layers of products. For this reason SLA technology requires more printers to fulfill the annual demand of products than SLS does. In particular, SLA requires 6 printers while SLS needs only 2.  Work cost is similar for the two technologies. SLS requires more time for finishing operations, while SLA requires longer time for the setup activities.

Considering that the evolution of SLA and SLS technologies has enhanced precision and resolution of these technologies onto very high (and almost similar) levels, the decision concerning which of them to implement is closely related to the total manufacturing process cost. Therefore, taking into account the economic assessment provided, and the strategic alternative about the technical performance of the two technologies, the company is now evaluating the possibility to (gradually) change its technology and substitute the old SLA machine with a new SLS machine.

Conclusion

There are two points about AM that draws a lot of attention among concerning researchers in the field. First, the rather young age of AM compared with the traditional and conventional technologies and second, the ongoing process towards full adoption of AM as a viable manufacturing system in the industrial world. However, it must be noted that due to the incomplete maturity and ongoing research, many of AM's aspects including, but not limited, to process measurement and standardization, finish surface quality, throughput rate, raw material selection, still lack enough competence to replace conventional technologies and become a widely accepted manufacturing system. In order to define a more holistic approach to AM, implementation of more case studies and accurate tests deem necessary to provide guidelines which help to identify the threshold value for the abovementioned four drivers to immediately discriminate which products should be subjected to a technical and economic evaluation, and which should be excluded from further analysis.