Randomized Dimension Reduction for Monte Carlo Simulations
Nabil Kahale





























dcu







exampleExample[section]
theoremTheorem[section]

conjectureConjecture[section]
lemmaLemma[section]
propositionProposition[section]
remarkRemark[section]
corollaryCorollary[section]
definitionDefinition[section]
equationsection

Randomized Dimension Reduction for Monte Carlo Simulations
    Nabil Kahalé
ESCP Europe, Labex ReFi, Big data research center,  75011 Paris,
France; e-mail: nkahale@escpeurope.eu.

    December 30, 2023
=========================================================================================================================






 


















 [, ](),a,
 
 
 24pt
  
 and













 
                 




Kahalé
Randomized Dimension Reduction
Randomized Dimension Reduction for Monte Carlo Simulations

Nabil Kahalé

ESCP Europe, Labex Réfi, Big data research center, 75011 Paris, France, nkahale@escpeurope.eu 
 

dimension reduction,   variance reduction,  effective dimension, Markov chains, Monte Carlo methods

informs2014 



We present  a new unbiased algorithm that  estimates the expected value of f(U) via Monte
Carlo  simulation, where U is a vector of d independent random variables, and f is
a function of d variables. We assume that f does not depend equally on all its arguments.  Under certain conditions we prove that, for the same computational cost, the variance of  our estimator is lower than the variance of the standard Monte Carlo estimator by a factor of order  d.  Our method can be used to obtain  a low-variance unbiased estimator  for the expectation  of a function of the state of  a Markov chain at a given time-step. We study applications to  volatility forecasting and time-varying queues. Numerical experiments show that our algorithm dramatically improves upon the standard Monte Carlo method for large values of d, and is highly resilient to discontinuities.
     

Keywords: 
dimension reduction;   variance reduction;  effective dimension; Markov chains; Monte Carlo methods
Randomized Dimension Reduction for Monte Carlo Simulations
    Nabil Kahalé
ESCP Europe, Labex ReFi, Big data research center,  75011 Paris,
France; e-mail: nkahale@escpeurope.eu.

    December 30, 2023
=========================================================================================================================




§ INTRODUCTION

Markov chains arise in a variety of fields such as finance, queuing  theory,  and social networks.   While much research has been devoted to the study of steady-states of Markov chains, several practical applications rely on the transient behavior of Markov chains. For example, the volatility of an index can be modelled as a Markov chain using the GARCH model <cit.>. Financial institutions conducting  stress tests may need to estimate the probability that the volatility  exceeds a given level in a few years from now. Also, due to  the nature of human activity, queuing systems  in  areas such as health-care, manufacturing, telecommunication and transportation networks, have often time-varying features and do not have a steady-state. For instance, empirical data show significant daily variation in  traffic in wide-area networks <cit.> and vehicular flow on roads <cit.>. Estimating the expected delay of packets in    a wide-area network at a specific time of the day (12pm, say) could be used to dimension such  networks. Similarly, estimating the velocity of cars in a region at 6pm could be used to design  transportation       
networks. In the same vein, consider the problem of estimating the queue-length at the end of a business day
in a call center that operates with fixed hours. In such call centers, knowing how many
calls would still need to be answered at 5pm could be an important metric that would
be needed in estimating their staffing requirements.    Methods to determine appropriate
staffing levels in call centers and other
many-server queueing systems with time-varying
arrival rates have been designed in <cit.>. Also, approximation tools have been developed to study  time-varying queues (see <cit.> and references therein).  However, in many situations, there are no  analytical tools, except Monte Carlo simulation,  to study accurately systems modeled by a  Markov chain.   A drawback of Monte Carlo simulation is its high computation cost. This motivates the need to design efficient simulation tools to study the transient behavior of  Markov chains, with or without time-varying features.

This paper gives a new  unbiased algorithm to estimate
E(f(U)),   where      U=(U_1,…,U_d) is a vector of d independent  random variables
    U_1,…,U_d
  taking values in a measurable space F, and  f  is a real-valued Borel-measurable function on F^d such that 
f(U) is square-integrable. For instance,  F can be equal to ℝ or to any vector space over ℝ. Under certain conditions, we show that  our  algorithm yields substantial lower variance than the standard Monte Carlo method for the same computational effort.  Our techniques can be used to efficiently estimate the expected value of a function of the state of a  Markov chain at a given time-step d, for a class of  Markov chains driven by independent random variables. An alternative  algorithm for Markov chains estimation, based on Quasi-Monte Carlo  sequences, that substantially improves upon standard Monte Carlo in certain numerical  examples, is given in <cit.>, with 
bounds on  the variance proven for special situations where the state space of the chain is a subset of the real
numbers. 

In a standard Monte Carlo scheme, E(f(U))    is estimated by simulating  n  independent vectors in F^d having the same distribution as U,   and taking the average of f over the n vectors.  In the related  Quasi-Monte Carlo method (see <cit.>),  f is evaluated at a predetermined deterministic sequence of points.
  In several  applications, the efficiency of   Quasi-Monte Carlo algorithms  can be improved by reordering the U_i's and/or making a change of variables, so that the value of  f(U) depends  mainly on the first few  U_i's.
For instance, the Brownian bridge construction and principal components analysis have been
used <cit.>  to reduce the error in the valuation of  financial derivatives via Quasi-Monte Carlo methods (see <cit.>
for related results). The relative importance of the first variables can  formally be measured by calculating the effective dimension in the truncation sense, a concept defined in <cit.>: when the first variables are important, the effective dimension in the truncation sense  is low in comparison to the nominal dimension.  It is proven in <cit.>  that Quasi-Monte Carlo methods are  effective for a class of functions where the importance of U_i decreases with i. <cit.>  apply Quasi-Monte Carlo methods to   queueing simulation    and option pricing, and examine their connection
to the effective dimension.  The truncation dimension and a related notion, the effective dimension in the superposition sense, are studied in <cit.>. It is shown in <cit.>  that the Brownian bridge and/or principal components analysis algorithms substantially reduce the truncation dimension of certain financial instruments. Alternative linear transformations have been
proposed in <cit.>  to reduce the effective dimension of financial derivatives and improve the performance of Quasi-Monte Carlo methods. 

Other previously known variance reduction techniques have exploited the  importance of certain variables  or states. 
For instance,     
stratified sampling along    important directions
 is used in pricing  path-dependent options <cit.>.  Importance sampling methods aim to increase the number of  samples that hit an important set via a change of measure technique <cit.>.   
When d=2 and  f(U_1,U_2) is more influenced by U_1 than by U_2,   and the expected time to generate U_1 is much lower than the expected time to generate U_2, the splitting   technique <cit.> simulates several independent copies of U_1 for each copy of U_2.  <cit.> give
the variance of the splitting estimator and the optimal  number of copies of U_1,
and show that the splitting technique is related to the conditional  Monte Carlo method. Multilevel splitting techniques are often  used for variance reduction in  the estimation of rare event probabilities <cit.>. The idea is to  split each path that reaches an important region into a number of subpaths in order to produce more paths that hit the rare event set. The rare event probability is then evaluated via a telescoping product. <cit.> analyse  multilevel splitting techniques that estimate  functionals of Markov chains with a discrete state space and   of  ergodic Markov chains in their steady state.   <cit.> analyse the performance of multilevel splitting techniques for rare event estimation and give, under certain conditions,  the optimal degree of splitting as the probability of the event goes to 0. Multilevel splitting methods have had many applications, such as the estimation of network reliability <cit.> and of  rare events in Jackson networks <cit.>. Multilevel splitting techniques for rare event simulation with finite time constraints are analysed in <cit.>. A comprehensive survey  on multilevel splitting techniques with applications to rare event simulations, sampling from complicated distributions, Monte Carlo counting, and randomized optimization, can be found in <cit.>. 

     Another technique, the multilevel Monte Carlo 
(MLMC) method introduced in <cit.>,  which relies on  low dimensional approximations of the function to be estimated,  dramatically reduces the computational complexity of estimating an expected value arising
from a stochastic differential equation.  Related randomized multilevel methods that produce  unbiased estimators for equilibrium expectations of functionals defined on homogeneous Markov chains have been provided in <cit.>. These methods apply to    the class of positive Harris recurrent Markov chains, and to chains that are contracting on average.   It is shown in <cit.> that similar randomized multilevel methods can be used to efficiently compute  unbiased estimators for expectations of functionals of solutions to stochastic differential equations. The MLMC method has had numerous other applications (e.g., <cit.>). 
  
 




 

  The basic idea behind our algorithm is that, if  f  does not depend  equally on all its arguments, the standard Monte Carlo method can be inefficient because it simulates all d arguments of f at each iteration. Assuming that  the expected time needed to simulate f(U) is  of order  d and that the variance of  f(U) is upper and lower-bounded by constants, the expected time needed to achieve  variance  ϵ^2 by standard Monte Carlo simulation is  Θ(dϵ^-2). In contrast, our algorithm simulates  at each iteration  a random subset  of  arguments of f, and reuses the remaining arguments from the previous iteration.
Under certain conditions,
we show that  our algorithm estimates E(f(U)) with variance  ϵ^2 in O(d+ϵ^-2) expected time. We also establish central limit theorems
on the statistical error of our algorithm.
When d=2, our method is very similar to the splitting technique in <cit.>. Our approach can thus be viewed as a   multidimensional version of this technique. However, in contrast with
existing multilevel splitting algorithms
where splitting decisions   typically depend on the   current state, in our method, the  arguments of f to be redrawn  are independent of  previously generated copies of  U.  
   
 

In order to  optimize  the performance of our estimator, we minimize the asymptotic product of the variance and expected running time, in the same spirit as stratified sampling <cit.>, the splitting technique <cit.>),   MLMC <cit.>, and related methods <cit.>.
   This minimization is performed via a new geometric algorithm that solves in O(d) time a d-dimensional optimisation problem. Our geometric algorithm is of independent interest and can be used to solve an optimization problem of the same type that was  solved in <cit.> in   O(d^3) time. We are not aware of other previous algorithms that solve this problem.
 This work extends the research in <cit.>, where the variance properties of the randomized estimator presented in this paper were announced without  proof.

Our method has the following features:



  * Under certain conditions,
it estimates E(f(U)) with variance  ϵ^2 in O(d+ϵ^-2) expected time. We are not aware of any previous method that achieves, under the same conditions, such a tradeoff between the  expected running time and accuracy. In contrast with stratified sampling, which can be performed in practice  only along a small number of dimensions <cit.>, our method 
is targeted at high-dimensional problems. The efficiency of our method typically increases with d, even though it can be used in principle for any  d≥2. 
  * It is easy to implement, does not make any continuity assumptions on f, nor does it    require  a detailed knowledge of the structure of f or U.
In contrast,  importance sampling,  multilevel    splitting  and multilevel Monte Carlo methods can achieve substantial variance reduction by exploiting the  structure of the simulated model. The standard Quasi-Monte Carlo method does not    necessitate  a detailed knowledge of the model structure, but makes regularity assumptions on the function to be integrated, and its efficiency does not  increase with d. 

        
  
The
rest of the paper is organized as follows.   
<ref> presents our generic randomized dimension reduction  algorithm and analyses its performance. <ref>  describes the aforementioned geometric algorithm and gives a numerical  implementation of the randomized dimension reduction algorithm.    <ref>  provides applications to Markov chains. <ref> presents and analyses a deterministic version of our algorithm. <ref>   compares our algorithm to a class of MLMC algorithms. 
 <ref> gives numerical simulations. <ref> contains concluding remarks. Most proofs are contained in the appendix or the on-line supplement.
The connection between our approach and the ANOVA decomposition and truncation dimension is studied in the appendix.  The appendixon-line supplement explores further the relation between our method,   the splitting technique, and the conditional Monte Carlo method, and contains more numerical simulations. 
   

    

§ THE  GENERIC RANDOMIZED DIMENSION REDUCTION ALGORITHM
 


 §.§ The algorithm description

We assume that all random variables in this paper are defined on the same probability space (Ω,ℱ,ℙ). Our algorithm
   estimates E(f(U)) by performing n iterations, where n is
   an arbitrary positive integer. The algorithm samples more often the first arguments of f than the last ones. It implicitly assumes that, roughly speaking, the importance  of the i-th argument of  f  decreases with i.   In many Markov chain examples, the last random variables are more important than the first ones, but our algorithm  can still be used efficiently after re-ordering the random variables, as described in detail in <ref>. A general  approach to rank input variables according to their importance is described in <cit.>, but we will not   use such an approach in our examples.  
 

Let 
A={(q_0,…,q_d-1)∈ℝ^d:1=q_0≥q_1≥⋯≥q_d-1>0}.
Throughout the paper,      q=(q_0,…,q_d-1) denotes
an element of A.  Our   generic algorithm takes such a vector      q as parameter. Let (N_k), 
k≥ 1,  be    a sequence of  independent random integers in [1,d] such that ℙ(N_k>i)=q_i for 0≤ i≤ d-1 and
k≥ 1.  The algorithm simulates n copies V^(1),…,V^(n) of U and consists of
the following steps:


  * First iteration.
Simulate a vector V^(1) that has the same distribution as U and calculate f(V^(1)).


  * Loop. In iteration k+1, where 1≤ k≤ n-1, let V^(k+1) be the vector obtained from  V^(k) by redrawing the first N_k components of V^(k), and keeping the remaining components unchanged.
Calculate f(V^(k+1)).

  * Output the average  of f(V^(1)),…,f(V^(n)). 


More formally, 
 consider a sequence (U^(k)), 
k≥ 1,  of independent  copies of U such that the two sequences  (N_k), 
k≥ 1, and  (U^(k)), 
k≥ 1, are independent.   Define the sequence (V^(k)),  k≥1,  in F^d    as follows: V^(1)=U^(1) and, for k≥1, the first N_k components of V^(k+1) are the same as the corresponding components of U^(k+1), and the remaining components of V^(k+1)  are the same as the corresponding components of V^(k).   The algorithm then outputs 
f_n≜f(V^(1))+…+f(V^(n))/n.

Note that f_n is an unbiased estimator of E(f(U)) since V^(k)  d  = U for   1≤ k≤ n. 


 §.§ Performance analysis

For ease of presentation, we  ignore the time needed to generate N_k and the running time of the third step of the algorithm. For 1≤ i≤ d, let t_i  be  the expected time needed to generate V^(k+1) and calculate f(V^(k+1)) when N_k=i.  Equivalently,  t_i is  the expected time needed to perform Step 2 of the algorithm when N_k=i. Thus, t_i  is  the expected time needed to re-draw the first i components of U and  recalculate  f(U), and  t_d is  the expected
time needed to simulate f(U).
By convention,  t_0=0. We  assume for simplicity that t_i is a strictly increasing function of i. In many examples (see  <ref> and <ref>), it can be shown that t_i=O(i).         As ℙ(N_k=i)=q_i-1-q_i for 1≤ i≤ d and k≥1, where q_d=0, the expected running time of a single iteration of our algorithm, excluding the first one, is equal to T, where

    T≜∑^d_i=1(q_i-1-q_i)t_i=∑^d-1_i=0q_i(t_i+1-t_i).
  
For 0≤ i≤ d, define  
    C(i)≜(E(f(U)|U_i+1,…,U_d)).
Thus, C(0)= (f(U)), while C(d)=0, and we can interpret C(i) as the variance captured by the last d-i components
of U. Note that if f depends only on its first i arguments, then f(U) is independent of (U_i+1,…,U_d), and
so 
E(f(U)|U_i+1,…,U_d)=E(f(U)),
 which implies that C(i)=0. More generally, if the last  d-i arguments of f are not important, the conditional expectation E(f(U)|U_i+1,…,U_d) is “almost” constant, and its  variance    C(i)  should be  small. Thus C(i)/C(0) can be used to measure the importance of
the last d-i components
of U. As shown in the appendix, when  U is uniformly distributed on the d-dimensional unit cube [0,1]^d, this ratio coincides with a global sensitivity index for the subset {i+1,…,d},  defined in <cit.>
 in terms  of the ANOVA decomposition of f. Proposition <ref> below shows that (C(i)), 0≤ i≤ d,  is always a decreasing sequence,  and gives an alternative   expression for C(i),  which can be viewed as a  variant of Theorem 2 of <cit.>.   

The sequence (C(i)), 0≤ i≤ d, is decreasing. If U'_1,… ,U'_i are  random variables such that  U'_j  d  = U_j  for 1≤ j≤ i, and   U'_1,… ,U'_i,  U are independent, then   

    C(i)=(f(U),f(U'_1,… ,U'_i,U_i+1,…,U_d)).


Theorem <ref> below establishes a formal relationship between the variance of f_n and the C(i)'s. Let  ν^* be the vector of ℝ^d+1 with  ν^*_0=C(0) and ν^*_i=2C(i) for 1≤ i≤ d. 
For n≥1,

    n(f_n)≤∑^d-1_i=0ν^*_i-ν^*_i+1/q_i.
Furthermore, the LHS of (<ref>) converges to its RHS as n goes to infinity. 

As, for ν=(ν_0,…,ν_d)∈ℝ^d×{0},   
    ∑^d-1_i=0ν _i-ν _i+1/q_i=ν _0+∑^d-1_i=1ν _i(1/q_i-1/q_i-1),
the RHS of (<ref>) is a weighted combination of the C(i)'s, with positive weights. Thus, the smaller  the C(i)'s, the smaller  the RHS of (<ref>). Furthermore, as C(i) is the variance of   the conditional expectation E(f(U)|U_i+1,…,U_d), which can be considered as  a smoothed version of f(U), we expect our algorithm to be resilient to discontinuities of f. 

Denote by ℝ_+  the set of nonnegative real numbers.
For q∈ A, and ϑ=(ϑ_0,…,ϑ_d)∈{0}×ℝ_+^d, and ν=(ν_0,…,ν_d)∈ℝ_+^d×{0},     set  
   
    R(q;ϑ,ν)=(∑^d-1_i=0q_i(ϑ_i+1-ϑ_i))(∑^d-1_i=0ν _i-ν _i+1/q_i).

The expected time needed to perform n iterations of the algorithm, including  the first one, is  T_n=(n-1)T+t_d. Theorem <ref> and (<ref>) imply that    T_n(f_n) converges to R(q;t,ν^*) as n goes to infinity, where t=(t_0,…,t_d).  By (<ref>),
 R(q;ϑ,ν) is an increasing function with respect to 
 ν, i.e. 
 R(q;ϑ,ν)≤ R(q;ϑ ,ν') for ν≤ν', where the symbol ≤ between vectors represents componentwise inequality. Similarly,  it is   easy to see that  R(q;ϑ ,ν) is increasing with respect to ϑ. Let   T^tot(q,ϵ) be   the total expected time  it takes for our algorithm to   guarantee that Std(f_n)≤ϵ.     Corollary <ref> below gives an upper bound on   T^tot(q,ϵ)  in terms of R(q;t,ν^*). It also  implies that,  if R(q;t,ν^*) is upper bounded by a constant independent of d,  and    (f(U))=Θ(1),   then our algorithm outperforms the standard Monte Carlo algorithm by a factor of order t_d. More precisely, running our algorithm for     n=⌈ t_dT^-1⌉ iterations has the same expected cost, up to a constant, as a single iteration of the standard Monte Carlo method, but produces an unbiased estimator of  E(f(U)) with O(1/t_d) variance.           
 

For ϵ>0,

    T^tot(q,ϵ)≤ t_d+ R(q;t,ν^*)ϵ^-2.
 Furthermore, if   n=⌈ t_dT^-1⌉,  the expected running time of n iterations of the algorithm is at most 2t_d, and
    (f_n)≤R(q;t,ν^*)/t_d.
   

Proof. Theorem <ref> and (<ref>) imply that n(f_n)T≤ R(q;t,ν^*). Thus,  Std(f_n)≤ϵ for n=⌈ R(q;t,ν^*)T^-1ϵ^-2⌉. The expected time needed to calculate f_n is  T_n, which is upper-bounded by t_d+ R(q;t,ν^*)ϵ^-2 since n-1≤  R(q;t,ν^*)T^-1ϵ^-2. Hence (<ref>).
On the other hand, if    n=⌈ t_dT^-1⌉,
then T_n≤2t_d since  (n-1)T≤ t_d, and (<ref>) holds since nT≥ t_d. 

Theorem <ref> below establishes a central limit theorem on f_n. It also establishes a central limit theorem on the estimate of  E(f(U)) that can be obtained with a  computational budget c, using the framework described by <cit.>. Denote by Ñ (c)  the number of iterations generated  by our algorithm incunits
of computation time. In other words,   Ñ (c) is the maximum integer n such that f_n is calculated   within c  time (with f_0≜0). As the time
to calculate f_n is random,   Ñ (c)   is a random integer. Let ⇒  denote weak convergence (see <cit.>).

As   n→∞, 

    √(n)(f_ n-E(f(U))) ⇒ N(0,σ^2 ),

whereσ^2=∑^d-1_i=0ν^*_i-ν^*_i+1/q_i.

Furthermore, as c→∞,

    √(c)(f_Ñ(c)-E(f(U))) ⇒ N(0,R(q;t, ν^*)).
 In light of above, we will use R(q;t,ν^*)  to  measure the performance of our algorithm. The smaller the C(i)'s and t_i's, the smaller  R(q;t,ν^*), and the better the performance
of our algorithm.  Proposition <ref>  below  shows that  C(i) is small if f is well-approximated  by a function of its first i arguments.      

For  1≤ i≤ d, if   f_i is a  measurable function from F^i  to ℝ such that f_i(U_1,…,U_i) is square-integrable, then 

    C(i)≤(f(U)-f_i(U_1,…,U_i)).


 §.§ Explicit and semi-explicit distributions
     
An optimal choice for q is a one that minimizes R(q;t,ν^*). A numerical algorithm that performs such minimization is presented in <ref>.  
This subsection gives explicit or semi-explicit choices for q,  with corresponding
upper-bounds on  R(q;t,ν^*).
  
Proposition <ref> below gives upper bounds on R(q;t,ν^*)    if    t_i=O(i)  and (C(i))  decreases at a sufficiently high rate. It implies in particular that, if    t_i=O(i)   and C(i)=O((i+1)^γ) with  γ<-1, then   T^tot(q,ϵ)=O(d+ ϵ^-2). 

Assume that d≥2 and there are constants c and c' and γ<0 independent of d such that t_i≤ ci and C(i)≤ c' (i+1)^γ for 0≤ i≤ d. Then, for q_i=(i+1)^(γ-1)/2,    0≤ i≤ d-1,  there     is a  constant c_1 independent of d such that 
    R(q;t,ν^*) ≤ c_1,    γ<-1, 
    
     c_1ln^2(d),    γ=-1, 
    
     c_1  d^γ+1,      -1<γ<0.
 and 

    T^tot(q,ϵ) ≤ c_2(d+ϵ^-2),    γ<-1, 
    
     c_2(d+ln^2(d)ϵ^-2),    γ=-1, 
    
     c_2(d+ d^γ+1ϵ^-2),     -1<γ<0.

Below is a simple example where C(0)=1 and the C(i)'s do not meet the conditions of
Proposition <ref>.  
Suppose  that F=ℝ and that  U_1,…,U_d are square-integrable real-valued  random variables with unit variance. Assume  that  f(x_1,…,x_d)=d^-1/2(∑^d_j=1x_j) for (x_1,…,x_d)∈ℝ^d.  As  f depends equally on its arguments,  our algorithm does not improve upon the standard Monte Carlo method.    Since  
E(f(U)|U_i+1,…,U_d)=d^-1/2(E(U_1+⋯+U_i)+U_i+1+⋯+U_d),
C(i)=(d-i)/d. The conditional variance C(i) decreases very slowly i since C(d/2) has the same order of magnitude as C(0). Thus the C(i)'s do not meet the conditions of
Proposition <ref>.  


When  upper-bounds on the C(i)'s and t_i's satisfying a convexity condition are known, Proposition <ref> below gives an explicit vector q together
with an upper bound on R(q;t,ν^*). 

Assume that t_i≤ϑ_i for 0≤ i≤ d, where ϑ_0,…,ϑ_d is a strictly increasing sequence with ϑ_0=0. Assume further that     ν_0,…,ν_d-1 are positive real numbers such that  C(i)≤ν_i for 0≤ i≤ d-1, and  that the sequence 

    θ_i=ν_i+1-ν_i/ϑ_i+1-ϑ_i,
0≤ i≤ d-1, is  increasing (by convention, ν_d=0). Then, for  q_i=√(θ_i/θ_0),   0≤ i≤ d-1,

    R(q;t,ν^*)≤ 2(∑_i=0^d-1√((ν_i-ν_i+1)(ϑ_i+1-ϑ_i)))^2.
Proof. 

 We first observe that θ_i≤θ_d-1<0 for 0≤ i≤ d-1. Thus q  is well-defined and belongs to A.  Let ϑ=(ϑ_0,…,ϑ_d). As t≤ϑ and  ν^*≤2ν, and since  R(q;.,.) is increasing with respect to its second and third  arguments, we have  R(q;t,ν^*)≤ R(q;ϑ,2ν).    This concludes the proof. 

 Proposition <ref>  below  yields an upper bound on  √(R(q;t,ν^*)) in terms of a weighted sum of the square roots of  the C(i)'s, for a semi-explicit vector q.      
      Assume  that C(d-1)>0. If, for 0≤ i≤ d-1, 
    q_i=√(t_1C(i)/t_i+1C(0)),
then 

    R(q;t,ν^*)
    ≤8(∑^d-1_i=0(√(t_i+1)-√(t_i))√(C(i)))^2.

Proposition <ref> below  gives an explicit  distribution which is   optimal up to a logarithmic factor, without  requiring any prior knowledge on the C(i)'s.       

For  any  q∈ A, 

    R(q;t,ν^*)≥∑^d-1_i=0C(i)(t_i+1-t_i).
  Furthermore, if q_i=t_1/t_i+1 for 0≤ i≤ d-1,  then 
    R(q;t,ν^*)≤ 2(1+ln (t_d/t_1))∑^d-1_i=0C(i)(t_i+1-t_i).


 §.§ A Lipschitz function example

Assume  that F=ℝ and that  U_1,…,U_d are square-integrable real-valued  random variables, with σ_1≥⋯≥σ_d>0,
where σ_i is the standard deviation of U_i. Assume also that  f(x_1,…,x_d)=g(∑^d_j=1x_j) for (x_1,…,x_d)∈ℝ^d, where g is a real-valued 1-Lipschitz function on ℝ that can be calculated in constant time. For instance,  f(x_1,…,x_d)=max(∑^d_j=1x_j-K,0), where K is a constant, satisfies this condition. Assume further that  each U_i   can  be simulated in constant time.  For 1≤ k≤ n, let S_k be the sum of all components of V^(k). Thus S_k+1 can be calculated recursively in O(N_k) time by adding to S_k  the first  N_k components of V^(k+1) and subtracting  the first  N_k components of V^(k). Hence  t_i≤ ci, for some constant c. 

In order to bound the C(i)'s, we show that f(U) can be approximated by f_i(U_1,…,U_i), where  f_i(x_1,…,x_i)=g(∑^i_j=1x_j+∑^d_j=i+1E(U_j))  for (x_1,…,x_i)∈ℝ^i. Let ||Z||=√(E(Z^2)) for a real-valued random variable Z. By Proposition <ref>, 
    C(i)   ≤   ||f(U)-f_i(U_1,…,U_i)||^2
       ≤    ||∑^d_j=i+1(U_j-E(U_j))||^2
       =   (∑^d_j=i+1U_j)
       =   ∑^d_j=i+1σ_j^2.
 
The second equation follows from the assumption that g is  1-Lipschitz. By applying Proposition <ref>, with ϑ_i=ci and ν_i=∑^d_j=i+1σ_j^2, and setting      q_i=σ_i+1/σ_1,  0≤ i≤ d-1, we infer that

    R(q;t,ν^*)≤2c(∑^d_i=1σ_i)^2.
  Thus, if  σ_i=O(i^γ), with γ<-1, then R(q;t,ν^*)=O(1) and T^tot(q,ϵ)=O(d+ ϵ^-2).
Also, by (<ref>) and a standard calculation, C(i)=O((i+1)^2γ+1)  for 0≤ i≤ d. As 2γ+1<-1, Proposition <ref> is also applicable in this case. 

  §.§.§ A polynomial counter-example
Assume that 
f(x_1,…,x_d)=∑^d_j=1x_jx_1^j.
 In general, in a standard implementation,  t_1=Ω(d),   since calculating f(U) after redrawing U_1  typically takes  Θ(d) steps. Thus, T=Ω(d), for any q∈ A. On the other hand, by (<ref>), the RHS of (<ref>) is lower-bounded by (f(U)). Thus, R(q;t,ν^*)=Ω(d(f(U)), which corresponds to to the time-variance product of the standard Monte Carlo algorithm. Therefore, a  standard implementation of the Randomized Dimension Reduction algorithm does not outperform the standard Monte Carlo algorithm in this example.  

§ THE OPTIMAL DISTRIBUTION

We now seek to calculate a vector q that minimizes R(q;t,ν^*). Given a vector  ν   in ℝ^d×{0} whose first d components are positive,  Theorem <ref> below gives a geometric algorithm that finds in O(d) time a vector q^* that minimizes R(q;t,ν) under the constraint that q∈ A. Note that the vector q^* depends on ν.
   In <cit.>, a dynamic programming algorithm that  calculates such a vector q^* in  O(d^3) time   has been described. 

Let ν'=(ν'_0,…,ν'_d)∈ℝ^d+1 be such that the set  {(t_i,ν'_i):0≤ i≤ d}
forms the lower hull of the set  {(t_i,ν_i):0≤ i≤ d}. In other words,  ν' is the supremum  of all sequences in ℝ^d+1 such that ν'≤ν and  the sequence   (θ_i) is  increasing, where  

    θ_i=ν'_i+1-ν'_i/t_i+1-t_i,
0≤ i≤ d-1.  For instance, if d=6, with t_i=i and ν=(20,21,13,8,7,2,0), then  ν'=(20,16,12,8,5,2,0), as illustrated in Fig. <ref>. <ref> shows how to calculate ν' in O(d) time. 


Let   ν    be a vector in ℝ^d×{0} whose first d components are positive. For 0≤ i≤ d-1,  set  q^*_i=√(θ_i/θ_0), where θ_i is given by (<ref>), and let q^*=(q^*_0,…,q^*_d-1). Then q^*=min_q∈ AR(q;t,ν), and 

    R(q^*;t,ν)=(∑_i=0^d-1√((ν'_i-ν'_i+1)(t_i+1-t_i))
    )^2.


 §.§ Lower hull calculation

   Given ν,   the following algorithm, due to <cit.>,  first generates recursively a subset B(j) of  {1,…,d}, 2≤ j≤ d, then calculates ν' via B(d). The algorithm runs in O(d) time.
  *  
Set B(2)={1,2}.

  *  For j=3 to d, denote by i_1<⋯<i_m  the elements of B(j-1).
Let k be the largest element of   {2,…,m} such that (t_i_k,ν_i_k) lies below the segment
[(t_i_k-1,ν_i_k-1),(t_j,ν_j)], if such k exists, otherwise let k=1. Set  B(j)={i_1,…,i_k,j}. 

  * 
For i=1 to d, let i' and i” be two elements of B(d) with i'≤ i≤ i”. Set ν'_i so that (t_i,ν'_i) lies on the segment [(t_i',ν_i'),(t_i”,ν_i”)]. 

 §.§ Estimating the C(i)'s

The calculation of a vector q^*∈ A  that minimizes  R(q;t,ν^*) requires the knowledge of the C(i)'s. Proposition <ref> below can be used
to estimate  C(i)    via Monte Carlo simulation. Assuming that f(U) can be approximated by a function of its first i arguments, we expect that both components of the product in the RHS of (<ref>)  to be small, on average. Thus, (<ref>) can be considered as a “control variate” version of (<ref>), and should yield a more accurate estimate of C(i) via Monte Carlo simulation for large values of i.

Assume that U'_1,… ,U'_d, and  U”_i+1,…,U”_d, are  random variables such that U'_j  d  = U_j  for 1≤ j≤ d, and  U^”_j  d  = U_j  for i+1≤ j≤ d, and   U'_1,… ,U'_d,  U,  U^”_i+1,…,U^”_d are independent. Then   

    C(i)=E((f(U)-f(U_1,… ,U_i,U'_i+1,…,U'_d)) 
    
     (f(U'_1,… ,U'_i,U_i+1,…,U_d)-f(U'_1,… ,U'_i,U^”_i+1,…,U^”_d)).


 §.§ Numerical algorithm
Building upon the previously discussed elements, the algorithm that we have used for our numerical experiments is as follows. It  constructs a vector (ν_0,…,ν_d) and uses it as a proxy for ν^*.   

  *  
For i=0 to d-1, if  i+1 is a power of 2, estimate  C(i)  by Monte Carlo simulation with 1000 samples via Proposition <ref>. 
  *  Set   ν_0=C_(0), and    ν_d=0.  For 1≤ i≤ d-1, let ν_i=2C(j),  where j is the largest index in [0,i] such that j+1 is a power of 2. 

  *  For i=d-1 down to 1, set ν_i←max(ν_i,ν_i+1). Set  ν_0←max(ν_0,ν_1/2).

  *  Let (ν'_0,…,ν'_d)∈ℝ^d+1 be such that the set  {(t_i,ν'_i):0≤ i≤ d}
forms the lower hull of the set  {(t_i,ν_i):0≤ i≤ d}. For 0≤ i≤ d-1,  set  q_i=√(θ_i/θ_0), where θ_i is given by (<ref>). 
  *   Calculate T via (<ref>).
For 0≤ i≤ d-1, set 
q_i←min(1,max(q_i,T/t_i+1ln(t_d/t_1))).

  *  Run Steps 1 through 3 of the generic randomized dimension reduction algorithm of <ref> using q.

The purpose of Steps 3 and 5 is to reduce the impact on  q of statistical errors that arise in Step 1.
Because of statistical errors,  Step 1 may underestimate or overestimate the C(i)'s. Step 3 guarantees that the ν_i's are non-negative, so that the q_i's can be calculated in Step 4. Step 5 yields a cap on the RHS of (<ref>) by ensuring that the q_i's are not too small.   Using  (<ref>) and the proof of Proposition <ref>, and assuming that t_d≥2t_1, it can be shown that Step 5  increases T  by at most a constant multiplicative factor. 
An alternative way to implement our algorithm
is to skip Steps 1 through 5 and  run the generic algorithm with  q_i=t_1/t_i+1 for 0≤ i ≤ d-1.   By Proposition <ref>, the resulting
vector q is optimal up to a logarithmic factor. 


§ APPLICATIONS TO MARKOV CHAINS

In queueing systems, the performance metrics at a specific time instant
are heavily dependent on the last busy cycle, i.e., the events that occurred after the
queue was empty for the last time. Thus, the performance metrics depend
a lot more on the last random variables driving the system than on the initial ones.
Nevertheless, we can apply our algorithm to queueing systems by  using a time-reversal transformation inspired from <cit.>. More generally, using such a time-reversal transformation, this section shows that our algorithm can efficiently estimate the expected value of a function of the  state of a Markov chain at time-step d, for a class of Markov chains driven by independent random variables.

 Let (X_m), 0≤ m≤ d, be a Markov chain with state-space
F' and deterministic initial value X_0.
Assume  that  there are independent  random variables Y_i,  0≤ i≤ d-1, that take values in  F, and measurable
functions g_i from F'× F to F' such
that X_i+1=g_i(X_i,Y_i) for 0≤ i≤ d-1.     We want to estimate  E(g(X_d))
for a given positive integer d, where g is a deterministic real-valued
measurable function on F' such that g(X_d) is square-integrable.
  For 1≤ i≤ d, set U_i=Y_d-i. It can be shown by  induction that  X_d=G_i(U_1,…,U_i,X_d-i), where G_i, 0≤ i≤ d, is a measurable function from F^i×
F' to F', and so there is a real-valued measurable function f on F^d
with g(X_d)=f(U_1,…,U_d). We can thus  use our  randomized dimension reduction algorithm to estimate E(g(X_d)). Recall that, in iteration k+1 in Step 2 of the generic algorithm  of  <ref>, conditioning on N_k=i, the first i arguments of f are re-drawn, and the remaining arguments are unchanged. This is equivalent to re-drawing the last i random variables driving the Markov chain, and keeping the first d-i variables unchanged. In light of above, the generic randomized dimension reduction  algorithm for Markov chains estimation takes as parameter a vector q∈ A  and consists of
the following steps:

  *  
First iteration.
Generate recursively  X_0,…,X_d.
Calculate g(X_d).

  *  Loop. In iteration k+1, where 1≤ k≤ n-1, keep  X_0,…,X_d-N_k
unchanged, and calculate recursively X_d-N_k+1,…,X_d by re-drawing  Y_d-N_k,…,Y_d-1, where N_k is a random integer in [1,d] such that ℙ(N_k>i)=q_i. Calculate g(X_d).

  *  Output the average  of g over the n copies of X_d generated in the first two steps.



We assume that g and the g_i's
can be calculated in constant time, and that the expected time needed to simulate each Y_i is upper-bounded by a constant independent of d. Thus, given N_k,   the expected time needed to perform iteration k+1   is O(N_k). Hence t_i≤ ci, for some constant c independent of d. Proposition <ref> below shows that, roughly speaking, C(i) is small if  X_d-i and  X_d are “almost” independent. 

For 0≤ i≤ d, we have C(i)= (E(g(X_d)|X_d-i)).
By Proposition <ref>, if there are constants c'>0 and γ<-1 independent of d such that C(i)≤ c'(i+1)^γ for 0≤ i≤ d-1, then  R(q;t,ν^*) is upper-bounded by a constant independent of d, where q_i=(i+1)^(γ-1)/2 for 0≤ i≤ d-1.
The analysis in <cit.>, combined with Proposition <ref>, suggests that   C(i)   decreases exponentially with i  for a variety of Markov chains.  

For x∈𝔽', and 0≤ i≤ d, let 

    X_i,x=G_i(U_1,…,U_i,x).
In other words,  X_i,x is the state of the chain at time-step d  if  the chain is at state x at time-step d-i.  Intuitively, we  expect X_i,x to be close to X_d for large i if X_d depends mainly on the last Y_j's. By Proposition <ref>,  if  g(X_i,x) is square-integrable,
 
    C(i)≤||g(X_d)-g(X_i,x)||^2.

In the following examples,  we prove   that under certain conditions,   R(q;t,ν^*) is upper-bounded by a constant independent of d  for an explicit vector q∈ A, and so T^tot(q,ϵ)= O(d+ ϵ^-2). 


 §.§ GARCH volatility model

In the  GARCH(1,1) volatility model (see <cit.>), the  variance  X_i  of an index return between day i and day i+1, as estimated at the end of day  i, satisfies the following recursion:

    X_i+1= w +α X_iY_i^2+β X_i,
i≥0, where w, α and β are positive constants with α+β<1, and Y_i,  i≥0, are independent standard Gaussian random variables. The variable Y_i is known at the end of day i+1. At the end of day 0, given  X_0≥0, a positive integer d and a real number z, we want to estimate ℙ(X_d> z). In this example, F=F'=ℝ, and g_i(x,y)= w +α xy^2+β x, with   g(u)=1{u> z} for u∈ℝ. Proposition <ref> below shows that C(i) decreases exponentially with i. 

There is a constant κ independent of d such that
C(i)≤κ(α+β)^i/2 for 0≤ i≤ d-1.

By applying Proposition <ref> with  ϑ_i=ci and  ν_i=κ(α+β)^i/2,  and setting q_i=(α+β)^i/4 for 0≤ i≤ d-1, we infer that R(q;t,ν^*) is  upper-bounded by  a constant independent of d.


 §.§ G_t/D/1 queue

Consider a queue  where   customers arrive at time-step
i, 1≤ i ≤ d, and are served by a single server in order of arrival. 
Service times are all equal to 1.  Assume the system starts empty at time-step
0, and that   A_i    customers arrive
at time-step i, 0≤ i≤ d, where  A_0=0 and the A_i's are independent square-integrable random variables. Let X_i be the number of customers waiting in the queue  at time-step i.  Then  X_0=0  and (X_i) satisfies the Lindley equation 
    X_i+1=(X_i+Y_i)^+,
for 0≤ i≤ d-1, with Y_i=A_i+1-1.   We want to estimate E(X_d). In this example, g is the identity function, F=F'=ℝ, and g_i(x,y)=(x+y)^+. Proposition <ref> below shows that C(i) decreases exponentially with i under certain conditions on the  service times. 
If  there are constants γ>0 and κ<1 independent of d such that

    E(e^γ Y_i)≤κ
  for 0≤ i≤ d-1, then  C(i)≤γ'κ^i for 0≤ i≤ d-1, where γ' is a constant independent of d. 
  By applying Proposition <ref> with ϑ_i=ci and ν_i=γ'κ^i,   and setting  q_i=κ^i/2 for 0≤ i≤ d-1, we conclude that, under the assumption of Proposition <ref>, R(q;t,ν^*) is  upper-bounded by  a constant independent of d. The assumption in Proposition <ref> can be justified as follows. Given i∈[0,d-1], if   E(A_i)<1 and the function h(γ)= E(e^γ Y_i) is bounded on a neighborhood of 0, then h'(0)=E(Y_i)<0. As   h(0)=1, there is  γ>0   such that  h(γ)<1, and (<ref>)  holds for  κ=h(γ).    The assumption in Proposition <ref> says that  γ and κ can be chosen independently of i and of d. 


 §.§ M_t/GI/1 queue 

Consider  a  M_t/GI/1 queue
 where   customers are served by a single server in order of arrival. 
We assume that   customers arrive
according to a Poisson process with positive and continuous time-varying rate λ_s≤λ^*, where λ^* is a fixed positive real number.  The service times are assumed to be  i.i.d. and  independent of the arrival times. Assume that the system starts empty at time 0. For simplicity, we assume that the number of customers that arrive in any bounded time interval is finite (rather than finite with probability 1). Consider a  customer present in the system at a given time s. If the customer   has been served for a period of length τ,  its remaining service time is equal to its service time minus τ, and if the customer   is in the queue,   its remaining service time is equal to its service time. The residual work W_s at time s is defined as the sum of remaining service times of customers present in the system at   s.   We want to estimate the expectation of W_θ, where  θ is a fixed  time.   Let   d=⌈λ^^*θ⌉,  and assume that d≥2. For 0≤ i≤ d, let X_i=W_iθ/d be the  residual work at time iθ/d. For 0≤ i≤ d-1, let Y_i be the vector that consists of  arrival and  service times of customers  that arrive during the interval (iθ/d,(i+1)θ/d].   
In this example, g is the identity function, F is  equal to the set of real-valued sequences with finite support, and   F'=ℝ.   Let 0≤ s<s'. If no costumers arrive in (s,s']  then W_s'=(W_s-s'+s)^+. On the other hand, if no costumers arrive in (s,s')  and a customer with service time S arrives at s', then W_s'=S+(W_s-s'+s)^+.  Thus, given the set of arrival and service times of customers that arrive in (s,s'], we can calculate iteratively W_s'   from W_s. This implies that  X_i+1 is a deterministic measurable function of X_i and Y_i, for 0≤ i≤ d-1. Proposition <ref> below shows that C(i) decreases exponentially with i under certain conditions on the arrival and service times.

For 0≤ s≤θ, let  Z_θ(s) be the cumulative service time of costumers that arrive in [s,θ]. Assume there are constants γ>0 and κ<1 independent of d such that, for 0≤ s≤ s'≤θ and s'-s≤1/λ^*,

    E(e^γ (Z_θ(s')-Z_θ(s)-1/λ^*))≤κ.
 Then C(i)≤γ'κ^i/2 for 0≤ i≤ d-1, where γ' is a constant independent of d.     By applying Proposition <ref> with ϑ_i=ci and ν_i=γ'κ^i/2, and setting q_i=κ^i/4 for 0≤ i≤ d-1, we conclude that, under the assumption of Proposition <ref>, R(q;t,ν^*) is  upper-bounded by  a constant independent of d. 
The assumption in Proposition <ref> can be justified as follows. For 0≤ s≤ s'≤θ and s'-s≤1/λ^*,
  the cumulative service times of customers that arrive in [s,s') is  Z_θ(s')-Z_θ(s). If E(Z_θ(s')-Z_θ(s))<s'-s   and   h(γ)= E(e^γ (Z_θ(s')-Z_θ(s)-1/λ^*)) is bounded on a neighborhood of 0, then h'(0)<0. Thus  h(γ)<1 for some γ>0 and (<ref>)  holds for  κ=h(γ).    The assumption in Proposition <ref> says that  γ and κ can be chosen independently of d, s and s'. 

 

 §.§ Time-dependent contracting Markov chains

In <cit.>, an algorithm that gives unbiased estimators  for equilibrium expectations associated with real-valued functionals defined on time-homogeneous contracting Markov Chains is described. This subsection defines the related notion of time-dependent contracting Markov Chain and shows that our algorithm ...
 

§ DETERMINISTIC DIMENSION REDUCTION

This section studies a deterministic dimension reduction algorithm that performs the same steps as the generic randomized dimension reduction algorithm of <ref>, but uses a deterministic integral sequence (N_k),   k≥  1,   taking values in [1,d], to estimate E(f(U)). As for the randomized algorithm, denote by   f_n the output of the  deterministic dimension reduction  algorithm, and by T_n its expected running time, where  n is the number of iterations, including the first one. The sequence  (N_k),   k≥  1, may depend on n.  As the algorithm generates n copies of U, 
the random variable  f_n is an unbiased estimator of E(f(U)). 

 Assume that C(d-1)>0 and let q̂=min_q∈ AR(q;t,C), where C denotes the vector (C(0),…,C(d)). The existence of q̂ follows from Theorem <ref>. Define the integers μ_0,…,μ_d-1 recursively as follows. Let μ_0=1 and, for 1≤ i≤ d-1,
let  μ_i be the largest multiple of μ_i-1 in the interval [0,1/q̂_i], i.e.μ_i=μ_i-1⌊1/μ_i-1q̂_i⌋.
 It can be shown by induction that μ_i is well-defined and positive. Let q̅   be the vector in A defined by q̅_i=1/μ_i,  for 0≤ i≤ d-1. As ⌊ x⌋≤ x<2⌊ x⌋ for x≥1, we have μ_i≤1/q̂_i<2μ_i. It follows that, for 0≤ i≤ d-1,  
    q̂_i≤q̅_i< 2q̂_i.
 Define the sequence  (N̅_k),   k≥ 1,    as follows: 
    N̅_k=max{i∈[1,d]:k is a multiple of  μ_i-1}.
  As μ_0=1, such  i always exists. For 0≤ i≤ d-1 and  k≥1, if N̅_k=j with j>i,  then  k is a multiple of μ_j-1, and so    k is a multiple of μ_i, since  μ_j-1/μ_i is an integer. Conversely, if     k is a multiple of μ_i then, by construction, N̅_k>i. Hence

    N̅_k>i ⇔ k≡ 0μ_i.
   Given i∈[0,d-1], the inequality N̅_k>i occurs once as k ranges in a set of μ_i   consecutive positive integers. The sequence   (N̅_k) can thus be considered as a deterministic counterpart to the random sequence   (N_k) generated by the randomized dimension reduction algorithm when q=q̅. 

Theorem <ref> 
below gives a lower bound on the performance of the deterministic dimension reduction algorithm for any sequence (N_k), k≥ 1, and analyses the algorithm when N_k=N̅_k for k≥ 1.   

 For n≥1 and any deterministic sequence (N_k), k≥ 1, 
    T_n(f_n)≥ R(q̂;t,C).
If    N_k=N̅_k for k≥ 1 then, for  n≥1, 
    T_n=t_d+∑^d-1_i=0⌊(n-1)q̅_i⌋(t_i+1-t_i),

and 

    n(f_n)≤∑^d-1_i=0C(i)-C(i+1)/q̅_i.

Furthermore, the LHS of (<ref>) converges to its RHS as n goes to infinity.  
Using again the framework of <cit.>, we measure the performance of an estimator via the work-normalized
variance, i.e. the product of the variance
and expected running time. If N_k=N̅_k for k≥1 then, by Theorem <ref>, 
    T_n(f_n)→ R(q̅;t,C)

as n goes to infinity. Furthermore, it follows from (<ref>) and (<ref>) that R(q̅;t,C)≤2R(q̂;t,C).
Thus, up to a factor of 2, the sequence (N̅_k) asymptotically  minimizes the work-normalized variance of the deterministic dimension reduction algorithm. Moreover, by definition of q̂, and since C≤ν^*, 

    R(q̂;t,C)
       ≤    R(q^*;t,C)
       ≤    R(q^*;t,ν^*).

    R(q̂;t,C)
    ≤ R(q^*;t,C)
    ≤ R(q^*;t,ν^*),
 
where  q^*=min_q∈ AR(q;t,ν^*). Hence R(q̅;t,C)≤2R(q^*;t,ν^*).
Similarly, as ν^*≤2C,
R(q^*;t,ν^*)≤ R(q̅;t,ν^*)≤ 2R(q̅;t,C).
 Thus, the asymptotic work-normalized variances of the randomized dimension reduction algorithm, with q=q^*,  and of the deterministic dimension reduction algorithm, with N_k=N̅_k for  k≥ 1,   are within a factor of 2 from each other. 

Proposition <ref> below shows that if, after generating the first copy of U, we generate  the next n(q_0-q_1) samples by only changing the first component of the U in the previous
iteration, and the next  n(q_1-q_2) samples by only changing the first two components of the U in the previous
iteration, and so on, the resulting algorithm is asymptotically less efficient than standard Monte Carlo.    
 Assume that C(d-1)>0. Let q∈ A with q_d-1<1. If  N_k=i for   1≤ i≤ d and integer k∈(n(1-q_i-1),n(1-q_i)], then n(f_n)→∞ as n goes to infinity.


§ COMPARISON WITH A CLASS OF  MULTILEVEL ALGORITHMS

We  compare our method to a class of MLMC algorithms, adapted from <cit.>,  that   efficiently  estimate E(f(U)) under the assumption that f is approximated,  in the L^2 sense, by functions of its first arguments. Under conditions described in <ref>, we prove that, up to a constant, the  randomized dimension reduction algorithm is at least as  efficient as this class of   MLMC  algorithms. It should be stressed, however, that there may exist other MLMC algorithms that estimate E(f(U)) more efficiently than  the class of MLMC algorithms described below. 


 §.§ The  MLMC algorithms description and analysis


 Let L be a positive integer and let (m_l), 0≤ l≤ L, be a strictly  increasing integral sequence, with m_0=0 and m_L=d. For 1≤ l≤ L, let  ϕ_l be a square-integrable random variable equal to a deterministic measurable function of U_1,…,U_m_l, with  ϕ_L=f(U). The  ϕ _l's are chosen so that,  as l increases, ϕ_l gets closer to f(U), in the L^2 sense.    For instance, L  can be proportional to ln(d),  the m_l's can increase exponentially with l, and     ϕ_l could equal  f(U_1,…,U_m_l,x,…,x_d-m_l), for some  x∈ F. For   1≤ l≤ L,   let  ϕ̂_l be the average of n_l independent copies of ϕ_l-ϕ_l-1 (with ϕ_0≜0), where n_l is a positive integer to be specified later. Assume that  the estimators   ϕ̂_1,…,ϕ̂_L are independent. As 
    E(f(U))=∑^L_l=1E(ϕ_l-ϕ_l-1),
ϕ̂=∑^L_l=1ϕ̂_l is an unbiased estimator of E(f(U)). Following the analysis in <cit.>, 

    (ϕ̂)=∑^L_l=1V_l/n_l,
where V_l≜(ϕ_l-ϕ_l-1) for 1≤ l≤ L.   The expected  time needed to simulate ϕ̂ is  T_ML≜∑^L_l=1n_lt̂_l, where    t̂_l  is the expected time needed to simulate ϕ_l-ϕ_l-1.   
As the variance of the average of n i.i.d. square-integrable random variables is proportional to 1/n,   for ϵ>0, we need ⌈(ϕ̂)ϵ^-2⌉  independent samples   of  ϕ̂   to achieve an estimator variance at most  ϵ^2. Thus the total expected time T^MLMC(ϵ) needed for the MLMC algorithm to estimate E(f(U)) with variance at most ϵ^2       satisfies the relation
   
    T^MLMC(ϵ)=Θ(T_ML+T_ML(ϕ̂)ϵ^-2).
  
The first term in the RHS of (<ref>) accounts for the fact that   ϕ̂  is simulated at least once. As shown in <cit.>, the work-normalized variance  T_ML(ϕ̂) is minimized when the n_l's are proportional to √(V_l/t̂_l) (ignoring the integrality constraints on
the n_l's), in which case 

    T_ML(ϕ̂)=(∑^L_l=1√(V_lt̂_l))^2.
 In line with <cit.>, if t̂_l=O(2^l) and ||ϕ_l-ϕ_L||^2=O(2^β l), with β<-1, where the constants behind the O-notation do not depend on d, then T_ML(ϕ̂) is upper-bounded by a constant independent of d.
This can be shown by observing that 
V_l≤||ϕ_l-ϕ_l-1||^2≤(||ϕ_l-ϕ_L||+||ϕ_l-1-ϕ_L||)^2.
 
Theorem <ref> below  shows that, under certain conditions,  the randomized dimension reduction method is, up to a  multiplicative constant,    at least as  efficient as the class of MLMC methods described above.  Indeed, under the assumptions of Theorem <ref>, by (<ref>),   T^tot(q,ϵ)=O(d+T_ML(ϕ̂)ϵ^-2). On the other hand,   T_ML≥t̂_L≥ĉd since m_L=d. Thus (<ref>) implies that T^MLMC(ϵ)≥ c'(d+T_ML(ϕ̂)ϵ^-2), for some constant c'.    

Assume that there are constants c and ĉ independent of d such that  t_i≤ ci for 1≤ i≤ d, and   t̂_l≥ĉm_l for 1≤ l≤ L, and that C(d-1)>0. Then  R(q;t,ν^*)≤ (32c/ĉ)T_ML(ϕ̂)  if, for 0≤ i≤ d-1,  

    q_i=√(C(i)/(i+1)C(0)).


§ NUMERICAL EXPERIMENTS

Our simulation experiments, using the examples in <ref>, 
were implemented  in the C++ programming language. The randomized dimension
reduction algorithm   (RDR)  was
implemented as described in <ref>. The deterministic dimension
reduction algorithm  (DDR) was implemented similarly with  N_k=N̅_k.
In both algorithms,     n  was chosen so that the expected total number of  simulations of the U_i's  in iterations 2 through n  is approximately   10d. The actual total number of simulations of the U_i's, denoted by “Cost” in   our computer experiments, is about 11d because it includes the d simulations of the first iteration.    

We have  implemented the multilevel algorithm (MLMC) described in <ref>, with L=⌊log_2(d)⌋+1,  and m_l=⌊2^l-Ld⌋ for 1≤ l≤ L, and   ϕ_l=f(U_1,…,U_m_l,X_0,…,X_0_d-m_l). The V_l's were estimated by Monte Carlo simulation with 1000 samples, and the  n_l's   were scaled up so that  the actual total number of simulations   of the U_i's   is about 11d. We have also implemented
a  randomized quasi-Monte
Carlo method (QMC) with a random shift <cit.>.   Our implementation  uses  the C++ program available  at 
http://web.maths.unsw.edu.au/˜fkuo/sobol<http://web.maths.unsw.edu.au/ fkuo/sobol>
 to generate d-dimensional Sobol sequences of length n=4096. For
practical reasons linked to computing time and storage cost, the QMC algorithm
was tested for d up to 10^4.



In Tables <ref> through <ref> and in the appendixon-line supplement, the variable Std  refers to the standard deviation of f_n for the RDR and DDR algorithms,  to the standard deviation of ϕ̂ for the MLMC algorithm, and to the
standard deviation of the Quasi-Monte Carlo estimator for the QMC algorithm.
 The variable Std and a 90% confidence interval for E(f(U)) were estimated using 1000 independent runs of these two algorithms. For the RDR algorithm, a 90% confidence interval for 
the variable Cost was reported as well. The  variance reduction factor VRF is defined asVRF=d(f(U))/Cost×Std^2.

We estimated (f(U)) by  using 10000 independent samples of U.   


 §.§ GARCH volatility model

Table <ref> shows results of our simulations of the GARCH volatility  model for estimating 
ℙ(X_d>  z),
with z=4.4×10^-5, X_0=10^-4,  α=0.06, β=0.9, and w=1.76×10^-6. As expected, the variable Cost is about 11d for the  RDR, DDR, and MLMC  algorithms. For these algorithms, the variable Cost × Std^2 is roughly independent of d, and  the variance reduction factors are roughly proportional to d.   In contrast, for the QMC algorithm, the variable Cost × Std^2 is roughly proportional to d, and the variance reduction factors are roughly constant. The  90% confidence interval  of
the RDR algorithm running time has a  negligible length in comparison
to the running time.  
The  RDR algorithm outperforms the MLMC algorithm by about a factor of 10, and the QMC
algorithm by a factor ranging from 5 to 20.
In all our numerical experiments, the DDR algorithm outperforms the RDR algorithm by a factor between
1 and 2.

 §.§ G_t/D/1 queue

Assume that  A_i has a Poisson distribution with time-varying rate λ_i=0.75+0.5cos(π i/50),  for 1≤ i≤ d, (recall that A_0=0).
  These parameters are taken from <cit.>. Table <ref> estimates E(X_d), and   Table <ref> gives VRFs in the estimation of  ℙ(X_d> z), for selected values of z. Once again, for  the  RDR, DDR, and MLMC algorithms, the variable Cost × Std^2 is roughly independent of d, and the variance reduction factors are roughly proportional to d. The VRFs of the  RDR  and DDR  algorithms in Table <ref>  are greater than or equal to the corresponding VRFs in Table <ref>, which confirms the resiliency of these algorithms to discontinuities of g. In contrast, the VRFs of the MLMC algorithm in Table <ref> are lower than the corresponding VRFs  in Table <ref>. The RDR algorithm outperforms the MLMC
algorithm by  a factor  ranging from 1 to 2  in Table <ref>,
and a factor ranging from 2 to 17  in Table <ref>. Table <ref> estimates E(X_d) for shifted values of d. The results in
Table <ref> are similar to those of Table <ref>,
but the values of E(X_d)   in Table <ref>
are significantly smaller  than those in Table <ref>. This
can be explained by observing that λ_d is maximized (resp. minimized)
at the values of d listed in Table <ref> (resp. Table <ref>).
   
  
 
 


 §.§ M_t/GI/1 queue

Assume that  λ_s=0.75+0.5cos(π s/50) for s≥0. These parameters are taken from <cit.>. Assume further that, for j≥1, the service time S_j for the j-th customer has a Pareto distribution with  ℙ(S_j≥ z)=(1+z/α)^-3   for z≥0, for some constant α>0.
  A simple calculation shows that E(S_j)=α/2.  In our simulations, we have set d=⌈θ⌉. Table <ref> gives  our simulation results   for estimating
ℙ(W_θ> 1) when α=2, and Table <ref> lists VRFs for estimating
ℙ(W_θ> 1) for selected values of α.
Here again, for  the  RDR, DDR, and MLMC algorithms, the variable Cost × Std^2 is roughly independent of d, and the VRFs are roughly proportional to d. The RDR algorithm outperforms the MLMC algorithm by a factor ranging from      1 to  10,
 depending on the value of α.
In Table <ref>,   the  RDR, DDR, and MLMC algorithms become less efficient as α increases. This can be explained by noting that, as  α increases, the length of the last busy cycle 
increases as well, which renders W_θ more dependent on the first Y_i's. 


§ CONCLUSION

We have described  a randomized dimension reduction algorithm that  estimates E(f(U)) via Monte
Carlo  simulation, assuming that f does not depend equally on all its arguments.  We formally prove
that under some conditions, in order to achieve an estimator variance  ϵ^2, our algorithm requires
O(d + ϵ^-2) computations as opposed to O(d ϵ^-2) under the standard Monte Carlo method. Our algorithm can be used to efficiently estimate the expected value of a function of the state of a  Markov chain at time-step d, for a class of Markov chains driven by random variables.  The
numerical implementation of our algorithm uses a new geometric procedure  of independent interest that solves  in O(d) time a d-dimensional optimisation problem that was  previously solved in O(d^3) time.    
We have argued intuitively that our method is resilient to discontinuities of f, and have described and analysed a deterministic version of our algorithm. Our numerical experiments confirm  that our approach highly outperforms the  standard Monte Carlo  method for large values of d, and show its high resilience to discontinuities. Whether our approach can be combined with the Quasi-Monte Carlo method to produce a provably efficient estimator is left for future work.  


§ ACKNOWLEDGMENTS
 This research has been presented at the  9th NIPS Workshop on Optimization for Machine Learning, Barcelona, December 2016, the Stochastic  Methods in Finance seminar at École des Ponts Paris-Tech, March 2018, and the 23rd International Symposium on Mathematical Programming,  Bordeaux, July 2018.    The author thanks  Bernard Lapeyre, seminar and conference participants for helpful conversations. He is grateful to three anonymous referees,  an anonymous associate editor, and Baris Ata (department editor), for insightful comments and suggestions. This work was achieved through the Laboratory of Excellence on Financial Regulation (Labex ReFi) under the reference ANR-10-LABX-0095. It benefitted from a French government support managed by the National Research Agency (ANR).
 


   





This can  explained by the fact that the performance of our algorithm (resp. multilevel algorithms) actually relies on the weak dependence of (resp. f(U)) on U_i+1,…,U_d, as i gets large, as shown formally in    Sections <ref> and <ref>. Furthermore, while smoothing techniques are often needed to efficiently implement multilevel <cit.> and QMC algorithms (see <cit.> and references therein) in the presence of discontinuities, no continuity assumptions are made in the analysis of our algorithm, and we did not  need to smooth f in our numerical experiments. This can be explained by the fact that  the performance of our algorithm relies on  properties of conditional expectations of f, which can be considered as  smoothed versions of f. Note that it is not always possible to smooth f explicitly, especially if F  is multi-dimensional.














§ PROOF OF PROPOSITION <REF>

For 0≤ i≤ d, let 
    f^(i)=E(f(U)|U_i+1,…,U_d).
 Thus, C(i)=(f^(i)). By the tower law, for 0≤ i≤ d-1, 
E(f^(i)|U_i+2,…,U_d)=f^(i+1),
 and so 
    C(i+1)=(E(f^(i)|U_i+2,…,U_d)).

As the   variance decreases by taking the conditional expectation,   it follows that C(i+1)≤ C(i), as desired. 

We now prove (<ref>).  Let W=(U'_1,… ,U'_i,U_i+1,…,U_d).
 Since U and W are conditionally independent given U_i+1,…,U_d, and E(f(W)|U_i+1,…,U_d)=f^(i),

E(f(U)f(W)|U_i+1,…,U_d))=(f^(i))^2.

Hence, by the tower law,
E(f(U)f(W))=E((f^(i))^2).

On the other hand, using the tower law once again, 
E(f(U))=E(f(W))=E(f^(i)),

and so the RHS of (<ref>) is equal to (f^(i)),
as required.


§ PROOF OF THEOREM <REF>

We first prove Lemma <ref> below, which follows by classical calculations (see e.g. <cit.>).
Let (Z_k), k≥1, be an homogeneous stationary Markov chain in ℝ^d, and let g be a real-valued  Borel-measurable function on ℝ^d such that g(Z_1) is square-integrable, and a_j=(g(Z_1),g(Z_1+j)) is non-negative for j≥0. Assume that  ∑^∞_j=1a_j is finite. Then 

    n^-1(∑^n_m=1g(Z_m))≤ a_0+2∑^∞_j=1a_j.
Furthermore, the LHS of (<ref>) converges to its RHS as n goes to infinity. 
Proof.
Since (Z_k), k≥1, is homogeneous and stationary, (g(Z_m)g(Z_m+j))=a_j for m≥1 and j≥0. Thus,  
    (∑^n_m=1g(Z_m))
       =   ∑^n_m=1(g(Z_m))+2∑_1≤ m<m+j≤ n(g(Z_m)g(Z_m+j))
       =   na_0+2∑^n_j=1(n-j) a_j.
Hence  

    n^-1(∑^n_m=1g(Z_m))= a_0+2∑^n_j=1n-j/na_j,
which implies (<ref>).  Using Lebesgue's dominated convergence theorem concludes the proof.

We now prove Theorem <ref>.
Since the C(i)'s and the variance of f_n remain unchanged if we
add a constant to f, we can assume without loss of generality that E(f(U))=0.   For 0≤ i≤ d-1, set p_i=1-q_i. Let  m<k be two integers in [1,n]. For 0≤ i≤ d-1, ℙ(max_m≤ j<kN_j≤ i)= p_i^k-m
and so, for 1≤ i≤ d-1, 
    ℙ(max_m≤ j<kN_j= i)= p_i^k-m-p_i-1^k-m.

 For i∈ [1,d], conditional on the event  max_m≤ j<kN_j= i, the first i  components of V^(m) and V^(k) are independent, and  the  last d-i  components of V^(m) and V^(k) are  the same. This is because, if N_j=i, with m≤ j<k,  the vector V^(m)  and the first  i  components of V^(j+1) are independent.   Thus,  by Proposition <ref>,
    E(f(V^(m))f(V^(k))|max_m≤ j<kN_j=i)=C(i).
 As   C(d)=0,
it follows from Bayes' formula that 

    E(f(V^(m))f(V^(k)))   =   ∑_i=1^d-1ℙ(max_m≤ j<kN_j= i)C(i)
       =   ∑^d-1_i=1(p_i^k-m-p_i-1^k-m)C(i).

Let a_j=(f(V^(1)),f(V^(1+j))). Thus, for j>0,  
    a_j=∑^d-1_i=1(p_i^j-p_i-1^j)C(i)
is non-negative, and 
    ∑_j=1^∞ a_j   =   ∑^d-1_i=1(p_i/1-p_i-p_i-1/1-p_i-1)C(i)
       =   ∑^d-1_i=1(1/1-p_i-1/1-p_i-1)C(i)
       =   -C(1)+∑^d-1_i=1C(i)-C(i+1)/q_i
 is finite. Since a_0=C(0), it follows that

    a_0+2∑^∞_j=1a_j=C(0)-2C(1)+2∑^d-1_i=1C(i)-C(i+1)/q_i.

We conclude the proof using Lemma <ref>.


§ PROOF OF PROPOSITION <REF>

 Let U'_1,… ,U'_i be random variables  satisfying the conditions of  Proposition <ref>, and W=(U'_1,…,U'_i,U_i+1,…,U_d). Since (U_1,…,U_i) and W are independent, (f_i(U_1,…,U_i),f(W))=0.
 Similarly, (f_i(U_1,…,U_i),f_i(U'_1,…,U'_i))=(f(U),f_i(U'_1,…,U'_i))=0.
 Thus, by Proposition <ref> and  bilinearity of the covariance,

    C(i)   =   (f(U)-f_i(U_1,…,U_i),f(W)-f_i(U'_1,…,U'_i))
       ≤   Std(f(U)-f_i(U_1,…,U_i)) Std(f(W)-f_i(U'_1,…,U'_i))
       =   (f(U)-f_i(U_1,…,U_i)).
 
The last equation follows by observing that f(W)-f_i(U'_1,…,U'_i) d=f(U)-f_i(U_1,…,U_i).


§ PROOFS OF PROPOSITIONS <REF> AND <REF>
 
We first prove the following.

Let ν=(ν_0,…,ν_d) be an element of ℝ^d×{0}. Assume that ν_0,…,ν_d-1 are positive, and that the sequence (ν_i/t_i+1),  0≤ i≤ d-1, is decreasing. For 0≤ i≤ d-1, set q_i=√((t_1ν_i)/(ν_0 t_i+1)).
Then 
R(q;t,ν)≤4(∑^d-1_i=0√(ν_i)(√(t_i+1)-√(t_i)))^2.
Proof.
Applying the inequality x-y≤2√(x)(√(x)-√(y)), which holds for   x≥0 and y≥0, to x=ν_i and y=ν_i+1 yields

    ∑^d-1_i=0ν_i-ν_i+1/q_i   ≤    2√(ν_0/t_1)∑_i=0^d-1(√(ν_i)-√(ν_i+1))√(t_i+1)
       =   2√(ν_0/t_1)∑_i=0^d-1√(ν_i)(√(t_i+1)-√(t_i)).
 
Similarly, since t_i+1-t_i≤2√(t_i+1)(√(t_i+1)-√(t_i)), ∑^d-1_i=0q_i(t_i+1-t_i)≤2√(t_1/ν_0)∑^d-1_i=0√(ν_i)(√(t_i+1)-√(t_i)).
Taking the product completes the proof.
We now prove Proposition <ref>. Set  t'=(t'_0,…,t'_d),
with t'_i=ci, 0≤ i≤ d, and let ν=(ν_0,…,ν_d-1,0), with   ν_i=c' (i+1)^γ, 0≤ i ≤ d-1.  Thus   R(q;t,ν^*)≤ R(q;t',2ν) since  t≤ t' and ν^*≤ 2ν.  
By Proposition <ref>, 
R(q;t',2ν)≤8cc'(∑^d-1_i=0(i+1)^γ/2(√(i+1)-√(i)))^2.
 As √(i+1)-√(i)≤(i+1)^-1/2 for i≥0,  it follows
that   

    R(q;t,ν^*)≤8 cc' (∑_i=1^di^(γ-1)/2)^2.
    
The inequality ∑_i=1^di^(γ-1)/2≤1+∫^d_1x^(γ-1)/2 dx
implies that 
    ∑^d-1_i=0i^(γ-1)/2≤ 1+2/(γ+1),    γ<-1, 
    
     1+ln(d),    γ=-1, 
    
     2  d^(γ+1)/2/(γ+1),     -1<γ<0.

We conclude that there is a constant c_1 such that (<ref>) holds. This concludes the proof of Proposition <ref>.   
Corollary <ref> then   implies the existence of a constant c_2 such that (<ref>) holds. 

We now prove Proposition <ref>.
 By applying Proposition <ref> to  ν=(2C(0),…,2C(d-1),0), it follows that  R(q;t,ν) is upper-bounded by the RHS of (<ref>). Since R(q;t,ν^*)≤ R(q;t,ν), this implies (<ref>).


§ PROOF OF PROPOSITION <REF>

Since R(q;t,ν) is increasing with respect to ν, 
    R(q;t,ν^*)   ≥    R(q;t,C(0),…,C(d))
       =   (∑^d-1_i=0C(i)-C(i+1)/q_i)(∑^d-1_i=0q_i(t_i+1-t_i)) .

However, as  ∑^d-1_j=0q_j
(t_j+1-t_j)≥ q_it_i+1 for 0≤ i≤ d-1, and since (C(i)) is a decreasing sequence,

    (∑^d-1_i=0C(i)-C(i+1)/q_i)(∑^d-1_j=0q_j(t_j+1-t_j))
    ≥∑^d-1_i=0(C(i)-C(i+1))t_i+1.

This implies (<ref>) since 

    ∑^d-1_i=0(C(i)-C(i+1))t_i+1
    =∑^d-1_i=0C(i)(t_i+1-t_i).



Assume now that q_i=t_1/t_i+1 for 0≤ i≤ d-1. Since R(q;t,ν) is increasing with respect to ν, 

    R(q;t,ν^*)   ≤    R(q;t,2C(0),…,2C(d))
       =   2(∑^d-1_i=0C(i)-C(i+1)/q_i)(∑^d-1_i=0q_i(t_i+1-t_i))
       =   2(∑^d-1_i=0(C(i)-C(i+1))t_i+1)(∑^d-1_i=0t_i+1-t_i/t_i+1) .

Furthermore, as (t_i+1-t_i)/t_i+1≤ln (t_i+1/t_i) for 1≤ i≤ d-1,∑^d-1_i=0t_i+1-t_i/t_i+1≤1+ln (t_d/t_1).
  
Using (<ref>) once again yields (<ref>). 

§ RELATION WITH THE ANOVA DECOMPOSITION AND THE TRUNCATION DIMENSION
This section assumes that f is a square-integrable function on [0,1]^d, and that each   U_i is  uniformly distributed on [0,1], with (f(U))>0. Consider a decomposition of   f in the following form:
    f=∑_Y⊆{1,…,d}f_Y,
 where f_Y is a measurable function on [0,1]^d and  f_Y(u) depends on u only through  (u_j)_j∈ Y,  for u=(u_1,…,u_d)∈[0,1]^d. For instance, f_∅ is a constant, f_{j}(u) is a function of u_j, and f_{j,k}(u) is a function of (u_j,u_k). The relation (<ref>) is called ANOVA representation of f if, for Y⊆{1,…,d},  any vector u∈[0,1]^d, and any  j∈ Y,

    ∫^1_0 f_Y(u_1,…,u_j-1,x,u_j+1,…,u_d) dx=0.

It can be shown <cit.> that there is a unique ANOVA representation of f, that  the f_Y's are square-integrable, and that   (f(U))=∑_Y⊆{1,…,d}σ_Y^2,

where σ_Y is the standard deviation of f_Y(U). Furthermore,  for 0≤ i≤ d-1,    

    E(f(U)|U_i+1,…,U_d)=∑_Y⊆{i+1,…,d}f_Y(U),
 and the covariance between f_Y(U) and f_Y'(U) is null if Y≠ Y'. Hence,   for 0≤ i≤ d-1, 

    C(i)=∑_Y⊆{i+1,…,d}σ_Y^2.
It follows that C(i) is equal to the variance corresponding to the subset {i+1,…,d}, as defined in <cit.>.
 Thus, C(i)/C(0) is equal to the global sensitivity index S_{i+1,…,d} for the subset {i+1,…,d}  (see <cit.>).
  
Proposition <ref> below  relates the performance of our algorithm to the truncation dimension d_t of f,
  defined in <cit.> as
    d_t≜∑_Y⊆{1,…,d},Y≠∅max(Y)σ_Y^2/(f(U)).

Under the conditions in Proposition <ref>, if  t_d=Θ(d) and d_t is upper bounded by a constant independent of d, the asymptotic (as n goes to infinity) work-normalized variance  of our algorithm
is O(ln(d)(f(U))), whereas the  work-normalized variance  of the standard
Monte Carlo algorithm is  Θ(d(f(U))).    

Assume that there is a real number c such that
t_i≤ ci for 1≤ i≤ d, and that q_i=1/(i+1) for 0≤ i≤ d-1. Then  
    R(q;t,ν^*)≤2c(1+ln(d))d_t(f(U)).
Proof.
For  0≤ i≤ d-1, we can rewrite (<ref>) as 

    C(i)=∑_Y⊆{1,…,d},Y≠∅1{i<min(Y)}σ_Y^2.
Hence, 
    ∑ ^d-1_i=0C(i)
       =   ∑_Y⊆{1,…,d},Y≠∅∑^d-1_i=01{i<min(Y)}σ_Y^2
       =   ∑_Y⊆{1,…,d},Y≠∅min(Y)σ_Y^2
       ≤    d_t(f(U)).

Let t'_i=ci for 0≤ i≤ d, and t'=(t'_0,…,t'_d). The proof of (<ref>) shows that this relation still holds if
t is replaced by any increasing sequence of length d+1 starting at
0. Replacing t by t' implies that  
    R(q;t',ν^*)
       ≤   2c(1+ln(d))∑ ^d-1_i=0C(i)
       ≤    2c(1+ln(d))d_t(f(U)).

Moreover, R(q;t,ν^*)≤ R(q;t',ν^*)  since t≤ t'. This completes the proof.


§ PROOF OF THEOREM <REF>

We use the  following proposition, whose proof follows immediately from (<ref>).If   ν∈ℝ^d×{0} and  ν'∈ℝ^d×{0}
are such that ν'≤ν,  and  q∈ A, then R(q;t,ν')≤ R(q;t,ν),
 with equality if  ν_0=ν'_0 and, for 1≤ i≤ d-1,  (ν_i-ν'_i)(q_i-1- q_i)=0. 
 By definition of
the lower hull,  (θ_i), 0≤ i≤ d-1, is an increasing sequence. Furthermore, θ_d-1<0 since it is equal to the slope of a segment joining (t_i,ν_i) to  (t_d,0), for some i∈[0,d-1]. Hence θ_i<0 for   0≤ i≤ d-1, and so q^* is well defined and belongs to A. Furthermore,   (ν'_i), 0≤ i≤ d, is a decreasing sequence, and ν'_d=ν_d=0. On the other hand, by (<ref>),

    R(q^*;t,ν')=(∑_i=0^d-1√((ν'_i-ν'_i+1)(t_i+1-t_i))
    )^2.
 Since, by the Cauchy-Schwartz inequality, for all non-negative sequences (x_i) and (y_i),
(∑^d-1_i=0√(x_iy_i))^2≤(∑^d-1_i=0x_i)(∑^d-1_i=0y_i),
it follows that R(q^*;t,ν')≤ R(q;t,ν') for q∈ A. Furthermore, by Proposition <ref>,   R(q;t,ν')≤ R(q;t,ν), and so  R(q^*;t,ν')≤ R(q;t,ν). On the other hand,  (ν_i-ν'_i)(q^*_i-1- q^*_i)=0 for 1≤ i≤ d-1. This is because, if  ν_i≠ν'_i, then the point (t_i,ν_i) does not belong to the lower hull of the set  {(t_j,ν_j):0≤ i≤ d}. Hence  (t_i,ν'_i) belongs to the segment ( (t_i-1,ν'_i-1),(t_i+1,ν'_i+1), which implies that θ_i-1= θ_i and q^*_i-1= q^*_i.   Thus, as   ν_0=ν'_0,   Proposition <ref> shows that    R(q^*;t,ν)= R(q^*;t,ν').  This implies (<ref>)
and that  R(q^*;t,ν)≤ R(q;t,ν) for q∈ A, as desired.


§ PROOF OF PROPOSITION <REF>

The random vectors W^(1)=(U_1,… ,U_i,U'_i+1,…,U'_d), W^(2)=(U'_1,… ,U'_i,U_i+1,…,U_d), and W^(3)=(U'_1,… ,U'_i,U^”_i+1,…,U^”_d) have the same distribution as U.
Hence 
E((f(U)-f(W^(1)))(f(W^(2))-f(W^(3)))=(f(U)-f(W^(1)),f(W^(2))-f(W^(3))).
As W^(1) and W^(2) are independent,   (f(W^(1)),f(W^(2))=0.
 Similarly,    (f(W^(1)),f(W^(3))=(f(U),f(W^(3))=0.
 By bilinearity of the covariance, it follows that 
E((f(U)-f(W^(1)))(f(W^(2))-f(W^(3)))=(f(U),f(W^(2)).

We conclude the proof using Proposition <ref>.


§ PROOF OF PROPOSITION <REF>

We first prove the following Markov property.

Let   i∈[0,d]. If H is a bounded random variable which is measurable with respect to the σ-algebra generated by X_i,Y_i,…,Y_d-1,  then  

    E(H|Y_0,…,Y_i-1)=E(H|X_i).
Proof.
Let ℋ be the vector space of bounded real-valued random variables H satisfying (<ref>). Clearly, the constant random variables belong to  ℋ. Let (H_m), m≥0, be an increasing sequence of  positive elements of    ℋ such that H=sup_m≥0H_m is bounded. For m≥0, 
    E(H_m|Y_0,…,Y_i-1)=E(H_m|X_i).
 By the conditional  Lebesgue dominated convergence theorem <cit.>, the LHS (resp. RHS) of (<ref>)   converges to E(H|Y_0,…,Y_i-1)  (resp. E(H|X_i)) as m goes to infinity, and so H∈ℋ.

 Let 𝒢  (resp.  𝒢')  be the set of  bounded real-valued  random variables which are measurable with respect to the σ-algebra generated by X_i (resp. (Y_i,…,Y_d-1)),  and let 𝒞 be the set of random variables of the form GG', with     G∈𝒢 and G'∈𝒢'.
For      G∈𝒢 and G'∈𝒢', 
    E(GG'|Y_0,…,Y_i-1)
       =   GE(G'|Y_0,…,Y_i-1)
    
       =   GE(G').
 
The first equation holds since X_i is a measurable function of Y_0,…,Y_i-1, which implies that G  is measurable with respect to the σ-algebra generated by Y_0,…,Y_i-1. The second equation follows from the independence of  (Y_i,…,Y_d-1) and (Y_0,…,Y_i-1). Similarly, since  (Y_i,…,Y_d-1) and  X_i are independent,

    E(GG'|X_i)=GE(G'|X_i)=GE(G').

Thus, GG'∈ℋ,  and so 𝒞⊆ℋ. As 𝒞 is closed under pointwise multiplication, by the monotone class theorem <cit.>, ℋ contains all bounded random variables which are measurable with respect to the σ-algebra generated by the elements of 𝒞. Since Y_i,…,Y_d-1 and X_i belong to 𝒞, it follows that  ℋ contains all bounded random variables which are measurable with respect to the σ-algebra generated by X_i,Y_i,…,Y_d-1. This completes the proof.

  
  As g(X_d) is a measurable function of (X_i,Y_i,…,Y_d-1), for any integer m, the random variable min(m,g^+(X_d)) belongs to ℋ. By the conditional  Lebesgue dominated convergence theorem, taking the limit as m goes to infinity implies that

    E(g^+(X_d)|Y_0,…,Y_i-1)=E(g^+(X_d)|X_i).
 A similarly equation holds for g^-(X_d).  Thus, 
    E(g(X_d)|Y_0,…,Y_i-1)=E(g(X_d)|X_i).

Replacing i with d-i implies that  

    E(g(X_d)|U_i+1,…,U_d)=E(g(X_d)|X_d-i).
Taking the variance of both sides concludes the proof.


§ PROOF OF PROPOSITION <REF>

For 0≤ i≤ d-1, set Z_i=α Y_i^2+β, so that X_i+1= w+Z_iX_i.  It can be shown by induction on i that  X_d= w+ w Z_d-1Z'+Z” for 2≤ i≤ d,  where 
 
    Z'=∑^i-1_j=1∏^j_k=2Z_d-k

  and

    Z”=Z_d-i⋯ Z_d-1X_d-i.

   By convention,  the product over an empty set is equal to 1. Since X_i,0 is equal to the state of the Markov chain at step d if X_d-i=0, it follows that X_i,0= w+ w Z_d-1Z'. Note that E(Z”)=(α+β)^iE(X_d-i) as E(Y_i^2)=1. By <cit.>, for 0≤ m≤ d, 
    E(X_m)≤max(X_0, w /1-α-β),
 and so 
E(Z”)≤max(X_0, w /1-α-β)(α+β)^i.
 Since the density of Y_d-1 is upper-bounded by 1/2,
for γ≤γ',

    ℙ(γ≤ Z_d-1≤γ')
       =   2 ℙ(√((γ-β)^+/α)≤ Y_d-1≤√((γ'-β)^+/α))
       ≤   √((γ'-β)^+/α) - √((γ-β)^+/α)
       ≤   √(γ'-γ/α).
 
The last equation follows from the inequality √(y')-√(y)≤√(y'-y), which holds for 0≤ y≤ y'. On the other hand, as X_i,0≤ X_d, 
g(X_d)-g(X_i,0) =  1     X_i,0≤ z<X_d, 
 
0    Thus, by (<ref>),
    C(i)   ≤   ||g(X_d)-g(X_i,0)||^2
       =   ℙ(X_i,0≤ z< X_d).
But 

    ℙ(X_i,0≤ z< X_d)
       =   ℙ (z-Z”- w/ w Z'<Z_d-1≤z- w/ w Z')
       =    E(ℙ (z-Z”- w/ w Z'<Z_d-1≤z- w/ w Z'|Z',Z”))
       ≤   E(√(Z”/α w Z'))
       ≤   κ(α+β)^i/2 ,
 
 where κ is a constant  that depends only on w, X_0, α and β. The third equation follows from (<ref>) and the independence of Z_d-1 and (Z',Z”). The last equation follows from Jensen's inequality and the inequality Z'≥1.
Hence C(i)≤κ(α+β)^i/2 for 2≤ i≤ d. This inequality also holds for 0≤ i≤1 by replacing κ with max(κ,1/(α+β)).

If Y and Y' are independent and X is a measurable function of  Y, and X' is a measurable function of (X,Y'), then E(X'|Y)=E(X'|X).


§ PROOF OF PROPOSITION <REF>

By classical calculations (see, e.g., <cit.>), it can be shown by induction on d that X_d=max _0≤ j≤ dS_j, where  S_j=∑^d-1_k=d-jY_k for 1≤  j≤ d, with S_0=0. For  0≤ i≤ d-1, since X_i,0  is the number of customers in the queue at time-step d if there are no costumers in the queue at time-step d-i, it can be shown by induction on d that X_i,0=max _0≤ j≤ iS_j.  Hence 
X_d=max(X_i,0,max _i+1≤ j≤ dS_j).
 Thus 

    X_d-X_i,0≤max _i+1≤ j≤ dS_j^+,
and so 
    C(i)   ≤   ||X_d-X_i,0||^2
       ≤   ∑ _j=i+1^d||S_j^+||^2.

On the other hand, since the Y_i's are independent,   E(e^γ S_j)≤κ^j for 0≤ j≤ d. Furthermore, as (x^+)^2≤2e^x for x∈ℝ, γ^2(S_j^+)^2≤ 2e^γ S_j. Taking expectations implies that γ^2||S_j^+||^2≤ 2κ^j, and so C(i)≤γ'κ^i, where γ'=2γ^-2/(1-κ). 

§ PROOF OF PROPOSITION <REF>

It can be shown by induction on the number of arrivals in (θ',θ]  (see also <cit.>)  that, if the system is empty at time θ'∈[0,θ], 

    W_θ=sup_s∈[θ',θ]{Z_θ(s)-(θ -s)}.

 Hence 

    X_d=sup_s∈[0,θ]{Z_θ(s)-(θ -s)}.
 Also, for 0≤ i≤ d, since X_i,0 is the residual work at time  θ if the residual work at time (d-i)θ/d is 0, 
    X_i,0=sup_s∈[(d-i)θ/d,θ]{Z_θ(s)-(θ -s)}.
   Thus,
by a calculation similar to the proof of Proposition <ref>,
    X_d-X_i,0≤sup_0≤ s≤(d-i) θ/d{Z_θ(s)-(θ -s)}^+.
For non-negative integer j, set 
S_j=Z_θ((θ-j/λ ^*)^+)-j-1/λ ^*.
Fix  s∈[0,(d-i) θ/d], and let j=⌈(θ-s)λ^*⌉. As θ-j/λ ^*≤ s and Z_θ is a decreasing function, 

Z_θ(s)≤ Z_θ((θ-j/λ ^*)^+).
 Since θ-s≥(j-1)/λ ^*, this implies that
Z_θ(s)-(θ -s)≤ S_j. But  2j≥ i as  j≥ iλ^*θ/d and λ^*θ≥ d/2. Hence 
    X_d-X_i,0≤sup_j≥ i/2S_j^+.

By calculations similar to the proof of Proposition <ref>, it  follows that
C(i)≤∑ _j≥ i/2||S_j^+||^2.
On the other hand, (<ref>) implies that   E(e^γ (S_j+1-S_j))≤κ for j≥0, and so  E(e^γ S_j)≤ e^γ /λ^*κ^j. We conclude the proof in a way similar to the proof of Proposition <ref>.

Supplementary Material

§ PROOF OF THEOREM <REF>

 If 𝒞 is a collection of random variables, denote by σ[𝒞] the  σ-algebra generated by the elements of 𝒞. Let (ξ_k), k≥1, be a stationary sequence of real-valued random variables. For k≥1, define the σ-algebras   ℱ_k=σ[ξ_n:1≤ n≤ k] and   ℱ'_k=σ[ξ_n: n≥ k].  Following <cit.>,
for n≥1, let φ_n=sup_k≥1sup{ℙ (B'|B)-ℙ (B'):B∈ℱ_k, ℙ(B)>0, B'∈ℱ'_k+n}.

The sequence (φ_n) is used to study the mixing proprieties of the sequence  (ξ_k). Theorem <ref> below establishes a functional central limit theorem on the partial sums of   (ξ_k)
under certain conditions on (φ_n).
Theorem <ref>  follows immediately from    Theorem 19.2 and the discussion on p. 203 of <cit.>. Let D_∞ denote the set of real-valued functions on the interval [0,∞)  that are right-continuous and have left-hand limits, endowed with the Skorokhod topology (see <cit.>).   
 Assume that  (ξ_k), k≥1, is stationary, that    ξ_1 is square-integrable, and that 

    ∑_n=0^∞√(φ_n)<∞.
 
Then the series 
    σ̅^2=(ξ_1)+2∑^∞_k=2(ξ_1,ξ_k)

converges absolutely. If σ̅>0  and S_n=∑^n_k=1(ξ_k-E(ξ_k)) then, as  n→∞, 
    S_⌊ ns⌋/√(n)⇒σ̅B_s
in the sense of D_∞, where B_s is a standard Brownian motion. 
We show that the conditions of Theorem <ref> hold when ξ_k=f(V^(k)) for k≥1.
As (V^(k)), k≥1, is a stationary Markov chain, the sequence (f(V^(k))), k≥1, is stationary. We first prove the following lemma, which follows intuitively from the fact that if all components of  the copy of U  are redrawn at a certain iteration, future copies of U are independent of past ones.
 For positive integers n and k, if B∈ℱ_k and B'∈ℱ'_k+n, then B and B'∩ I are independent, where 
I=I_k,n={ω∈Ω:∃ j∈[k,k+n-1],N_j(ω)=d}.
Proof.
   Let 𝒢_k=σ[N_j,U^(j+1):j≥ k].  
  Fix   l≥ k+n. By construction, if ω∈ I, there is j∈[k,k+n-1] is such that  V^(j+1)(ω)=U^(j+1)(ω).  Hence,  if   ω∈ I, the vector  V^(l) does not depend V^(k). More precisely, it can be shown by induction that there is a measurable function H:F^d(l-k)×ℝ^(l-k)→𝔽^d such that, for  ω∈ I,  

    V^(l)(ω)=H(U^(k+1)(ω),…,U^(l)(ω),N_k(ω),…,N_l-1(ω)).

For ω∈Ω, let
 
    G(ω)=f(H(U^(k+1)(ω),…,U^(l)(ω),N_k(ω),…,N_l-1(ω))).
   Thus,     f(V^(l)(ω))=G(ω) for  ω∈ I. As  the maps f and H are measurable, and since the random variables U^(k+1),…,U^(l) and N_k,…,N_l-1 are measurable with respect  to 𝒢_k,  the random variable G is measurable with respect  to 𝒢_k as well. 
  For  z∈ℝ,  let B_z={ω∈Ω:f(V^(l)(ω))<z}. Since  I∈𝒢_k and
B_z∩ I={ω∈Ω:G(ω)<z}∩ I,
  we have B_z∩ I∈𝒢_k. Thus B_z∈𝒢'_k for any real number z, where𝒢'_k={B∈ℱ:B∩ I∈𝒢_k}.
It is easy to check that  𝒢'_k is a σ-algebra.  We conclude that  f(V^(l)) is measurable with respect to 𝒢'_k, for any   l≥ k+n, and so   ℱ'_k+n⊆𝒢'_k.
Consequently, if B∈ℱ_k and B'∈ℱ'_k+n, then  B'∈𝒢'_k, and so B'∩ I∈𝒢_k. Since  ℱ_k and 𝒢_k are independent, so are   B and B'∩ I. 

We now prove the theorem. For positive integers n and k, let B∈ℱ_k with ℙ(B)>0 and B'∈ℱ'_k+n, and define I as in Lemma <ref>. Thus,ℙ(B'∩ I|B)=ℙ(B'∩ I).

 Since   ℙ( B')≤ℙ(B'∩ I)+ℙ( I^c),
 
it follows that
ℙ(B')≤ℙ( B'|B)+ℙ( I^c).
 
Replacing B' with its complement and noting that   ℙ( I^c)=(1-q_d-1)^n, we conclude that

|ℙ( B'|B)-ℙ(B')|≤ (1-q_d-1)^n.
  Hence  φ_n≤  (1-q_d-1)^n, and so (<ref>) holds. Thus the conclusions of Theorem <ref> hold for the sequence (f(V^(k))), k≥1. It follows from (<ref>) and (<ref>) that σ̅=σ.  
By (<ref>),  as  n→∞,

 
    ⌊ ns⌋(f_⌊ ns⌋-E(f(U)))/√(n)⇒σ̅B_s
 
 in D_∞. Setting s=1, which can be justified by applying the  continuous mapping theorem <cit.>  with the projection map <cit.>,   implies (<ref>). 
Let τ_k denote the
random amount time to generate V^(k) and calculate f(V^(k)). Thus E(τ_k)=T. By the strong law of large numbers, with probability 1, 
    1/n∑^n_k=1τ_k→ T
 as n goes to infinity.  By (<ref>), (<ref>) and <cit.>, it follows that 
    √(c)(f_Ñ(c)-E(f(U))) ⇒ N(0,Tσ^2 ),
 as c→∞. Since Tσ^2 = R(q;t, ν^*), this implies (<ref>).



§ PROOF OF THEOREM <REF> AND OF PROPOSITION <REF>
For convenience, set N_0=N̅_0=d. For 0≤ i≤ d-1, let
    𝒮_i=𝒮_i(n)={k∈[0,n-1]:N_k>i}.
 Let  u_i(0)=0,u_i(1),…,u_i(|𝒮_i|-1) be the elements of 𝒮_i  sorted in increasing order.  Set  u_i(|𝒮_i|)=n and let Q(i)=Q(i,n)=∑^|𝒮_i|_l=1(u_i(l)-u_i(l-1))^2. Note that 𝒮_0={0,…,n-1} and Q(0)=n. The following lemma gives an alternative characterization of Q(i).

For 0≤ i≤ d-1,
∑_1≤ m<k≤ n1{max_m≤ j<kN_j≤ i}=1/2(Q(i)-n).
Proof.
Given integers m and k with 1≤ m<k≤ n, the condition max_m≤ j<kN_j≤ i holds if and only if [m,k)∩𝒮_i=∅.
It is thus equivalent to the existence of an integer l∈[1,|𝒮_i|] such that u_i(l-1)<m<k≤ u_i(l).
Hence ∑_1≤ m<k≤ n1{max_m≤ j<kN_j≤ i}=1/2∑^|𝒮_i|_l=1(u_i(l)-u_i(l-1))(u_i(l)-u_i(l-1)-1).
 
We conclude the proof by noting that ∑^|𝒮_i|_l=1(u_i(l)-u_i(l-1))=n.
      
Lemma <ref> below relates the variance of f_n to the Q(i)'s and C(i)'s.

For n≥1,

    (f_n)=n^-2∑_i=0^d-1Q(i)(C(i)-C(i+1)).
Proof.
As in the  proof of Theorem <ref>, we  assume without loss of generality that E(f(U))=0. Let  m<k be two integers in [1,n]. By arguments similar to those that lead to (<ref>), 
E(f(V^(m))f(V^(k)))=C(max_m≤ j<kN_j).

Since, for 0≤ l≤ d, 
    C(l)   =   ∑_i=l^d-1(C(i)-C(i+1))
       =   ∑_i=0^d-11{l≤ i}(C(i)-C(i+1)),
it follows that 
    E(f(V^(m))f(V^(k)))=∑_i=0^d-11{max_m≤ j<kN_j≤ i}(C(i)-C(i+1)).

As
 
    (f_n)=n^-2(∑_m=1^nE((f(V^(m)))^2)+2∑_1≤ m<k≤ nE(f(V^(m))f(V^(k)))),

and E((f(V^(m)))^2)=C(0) for 1≤ m≤ n, we conclude that 
    (f_n)=n^-2(nC(0)+2∑_i=0^d-1∑_1≤ m<k≤ n1{max_m≤ j<kN_j≤ i}(C(i)-C(i+1))).

Thus, by Lemma <ref>,



    (f_n)=n^-2(nC(0)+∑_i=0^d-1(Q(i)-n)(C(i)-C(i+1))),
which, after some simplifications, completes the proof.
     
We now prove Theorem <ref>. By the Cauchy-Schwartz inequality, for 0≤ i≤ d-1 and  any sequence (x_l), 1≤ l≤𝒮_i, 
(∑^|𝒮_i|_l=1x_l)^2≤|𝒮_i|(∑^|𝒮_i|_l=1x_l^2).

Replacing x_l with u_i(l)-u_i(l-1) shows that n^2≤|𝒮_i|Q(i). By Lemma <ref>, it follows that 
    (f_n)≥∑_i=0^d-1C(i)-C(i+1)/|𝒮_i|.

 For 0≤ k≤ n-1, the expected running time of iteration k+1 is t_i if N_k=i. Furthermore, for i∈[1,d], there are |𝒮_i-1|-|𝒮_i|  integers k in [0,n-1] with N_k=i, where 𝒮_d≜∅, and so

    T_n=∑^d_i=1(|𝒮_i-1|-|𝒮_i|)t_i=∑^d-1_i=0|𝒮_i|(t_i+1-t_i).
Thus, 
    T_n(f_n)≥ R(q;t,C),
where q_i=|𝒮_i|/n for 0≤ i≤ d-1. This implies (<ref>). 

Assume now that N_k=N̅_k for k≥1.
Then, by (<ref>), for 0≤ i≤ d-1,
    𝒮_i={k∈[0,n-1]:  k  is a multiple of μ_i}.
  Thus  |𝒮_i|=1+⌊(n-1)q̅_i⌋ which, by (<ref>), implies (<ref>).   Since  u_i(l)=lμ_i for 0≤ l≤ |𝒮_i|-1,  

    Q(i)
       =   μ_i∑^|𝒮_i|-1_l=1(u_i(l)-u_i(l-1))
    +(n-max(𝒮_i))^2
       =   μ_imax𝒮_i+(n-max(𝒮_i))^2.
 
As 0≤ n-max(𝒮_i)≤μ_i, it follows that μ_i( n-μ_i)≤ Q(i)≤μ_in.
By Lemma <ref>, this implies (<ref>) and shows that, as n goes to infinity, the LHS of (<ref>) converges to its RHS. 


We now prove Proposition <ref>. Assume that (N_k) satisfies the conditions in the proposition, and that nq_d-1>1.  Then 𝒮_d-1 consists of the  integers belonging to {0}∪(n-nq_d-1,n), and so Q(d-1)≥ n^2(1-q_d-1)^2. By Lemma <ref>, it follows that    (f_n)≥ (1-q_d-1)^2C(d-1), which completes the proof.


§ PROOF OF THEOREM <REF>

We first prove the following lemma.   
 
 Let (ν_i), 0≤ i≤ d, be a decreasing sequence such that ν_m_l≤(ϕ_L-ϕ_l) for 0≤ l≤ L, with ν_d=0. Then ∑_i=0^d-1√(ν_i/i+1)≤2∑^L_l=1√(m_lV_l).
Proof.
Since the sequence (ν_i) is decreasing and (i+1)^-1/2≤2(√(i+1)-√(i)) for   i≥0,  ∑_i=j^k-1√(ν_i/i+1)≤2(√(k)-√(j))√(ν_j),
 for 0≤ j≤ k≤ d. Hence,

    ∑_i=0^d-1√(ν_i/i+1)   =   ∑^L-1_l=0∑_i=m_l^m_l+1-1√(ν_i/i+1)
       ≤   2∑^L-1_l=0(√(m_l+1)-√(m_l))√(ν_m_l)
       ≤   2∑^L-1_l=0(√(m_l+1)-√(m_l))Std(ϕ_L-ϕ_l)
       =   2∑^L_l=1√(m_l)(Std(ϕ_L-ϕ_l-1)-Std(ϕ_L-ϕ_l))
    
       ≤   2∑^L_l=1√(m_lV_l),
 where the last equation follows by sub-linearity of the standard deviation. 
We now prove Theorem <ref>.
Since ϕ_l is square-integrable and is a measurable function of  U_1,…,U_m_l,
by  Proposition <ref>, C(m_l)≤(ϕ_L-ϕ_l) for 0≤ l≤ L. By Proposition <ref>, the sequence (C(i)), 0≤ i≤ d, is decreasing, and so it satisfies the conditions of Lemma <ref>. Thus,  ∑_i=0^d-1√(C(i)/i+1)≤2∑^L_l=1√(m_lV_l).
Furthermore, by  Proposition <ref>, 

    R(q;t,ν^*)
    ≤8c(∑^d-1_i=0(√(i+1)-√(i))√(C(i)))^2.
 Since √(i+1)-√(i)≤(i+1)^-1/2, it follows that  

    R(q;t,ν^*)
    ≤32 c(∑^L_l=1√(m_lV_l))^2.
 Using (<ref>) concludes
the proof.


§ RELATION WITH  SPLITTING AND CONDITIONAL MONTE CARLO 

Like the splitting algorithm described in <cit.>
when d=2,
our method samples more often important random variables.    This section explores further the relation between our method and the splitting and conditional Monte Carlo methods.  Fix n≥1 and assume for simplicity that the sequence (N_k) is deterministic. We first analyse the relation between our method and the conditional Monte Carlo method in the general case, then show that
 the generic dimension reduction algorithm for Markov chains estimation can  be efficiently cast as a  splitting algorithm.


 §.§ Relation with conditional Monte Carlo

Define  𝒮_i  via (<ref>) for 0≤ i≤ d-1, and set  𝒮_d={0}.  For 0≤ i≤ d and m∈𝒮_i,   let 
    𝒮_i,m={k∈[m,n-1]:(m,k]∩𝒮_i=∅},
and
    f_i,m=1/|𝒮_i,m|∑^_k∈𝒮_i,mf(V^(k+1)).

Thus, if (m,m') are consecutive elements of 𝒮_i∪{n}, then   𝒮_i,m={m,…, m'-1}.  In particular, for 0≤ m≤ n-1, we have 𝒮_0,m={m},   and so  f_m,0=f(V^(m+1)). Similarly, 𝒮_d,0={0,…,n-1} and f_d,0=f_n. For 0≤ i≤ d,    m∈𝒮_i,   and k ∈𝒮_i,m, we have  N_l≤ i for m< l≤ k, and so  the  last d-i  components of V^(k+1) and V^(m+1) are  the same. In other words, 
    V^(k+1)_j=V_j^(m+1),    for  0≤ i<j≤  d  and  k ∈𝒮_i,m.
 
 We can thus view f_i,m  as a discrete analog to   E(f(V^(m+1))|V_i+1^(m+1),…,V_d^(m+1)).
The following proposition shows that the random variables f_i,m, for  0≤ i≤ d and m∈𝒮_i,  can be calculated inductively via (<ref>). This is reminiscent of the conditional Monte Carlo method, where an expectation  is estimated via an average of conditional expectations.  

For 1≤ i≤ d and m∈𝒮_i,  

    f_i,m=∑_k∈𝒮_i-1∩𝒮_i,m|𝒮_i-1,k|/|𝒮_i,m|f_i-1,k,
 and

    |𝒮_i,m|=∑_k∈𝒮_i-1∩𝒮_i,m|𝒮_i-1,k|.
Proof.
We  first show that 
    𝒮_i,m=⋃_k∈𝒮_i-1∩𝒮_i,m𝒮_i-1,k.
 If  k∈𝒮_i-1∩𝒮_i,m and k'∈𝒮_i-1,k,
we have (m,k]∩𝒮_i=∅ and  (k,k']∩𝒮_i-1=∅.
As 𝒮_i⊆𝒮_i-1 and m≤ k≤ k', this implies that  (m,k']∩𝒮_i=∅, and so   k'∈𝒮_i,m.    Thus  
    ⋃_k∈𝒮_i-1∩𝒮_i,m𝒮_i-1,k⊆𝒮_i,m.
    Conversely, given k'∈𝒮_i,m, let
    k=max([0,k']∩𝒮_i-1).
 Since m∈[0,k']∩𝒮_i-1, the integer k is well-defined and m≤ k.     Since k≤ k' and (m,k']∩𝒮_i=∅, we have (m,k]∩𝒮_i=∅. Hence k∈𝒮_i,m, and so  k∈𝒮_i-1∩𝒮_i,m.  Furthermore,      (k,k']∩𝒮_i-1=∅ by (<ref>).  Thus  k'∈𝒮_i-1,k, and so  
   
    𝒮_i,m⊆⋃_k∈𝒮_i-1∩𝒮_i,m𝒮_i-1,k.
This implies (<ref>).  Moreover,  if k∈𝒮_i-1 and j∈𝒮_i-1,k, then (k,j]∩𝒮_i-1=∅, and so k=max([0,j]∩𝒮_i-1). Thus,  if k and k' are distinct elements of 𝒮_i-1,
the sets 𝒮_i-1,k and 𝒮_i-1,k' are disjoint.
Together with (<ref>), this immediately implies (<ref>) and (<ref>).            


 §.§ The Markov chains case

Using the same notation and assumptions as in <ref>, we  show
how to cast the generic dimension reduction algorithm for Markov chains estimation as a  splitting algorithm. For each integer k in  [1,n], define the Markov chain  (X_i^(k)), 0≤ i≤ d, by induction on i as follows: X_0^(k)=X_0 and X^(k)_i+1=g_i(X^(k)_i,V^(k)_d-i) for 0≤ i≤ d-1. Then it can be shown by induction that X_d^(k)=G_d(V^(k),X_0), and that g(X_d^(k))=f(V^(k)). The generic dimension reduction algorithm for Markov chains estimation
described in <ref>  thus outputs the average of g(X_d^(1)),…,g(X_d^(n)).
Furthermore, it can be shown by induction that, for 0≤ i ≤ d and 1≤ k≤ n, the random variable X_i^(k) is a deterministic function of (V^(k)_j), d-i<j≤ d. On the other hand, by (<ref>), for 0≤ i≤ d and  0≤ k≤ n-1,   if m=max([0,k]∩𝒮_d-i) then k∈𝒮_d-i,m. The integer m is well-defined since 0∈𝒮_d-i. By (<ref>), the last i components of   V^(k+1) and V^(m+1) are  the same, and so X^(k+1)_i=X_i^(m+1). However, for 0≤ i≤ d-1 and k∈𝒮_d-i-1, we have V^(k+1)_d-i=U^(k+1)_d-i since N_k≥ d-i, and so 
    X^(k+1)_i+1=g_i(X^(m+1)_i,U^(k+1)_d-i),  with m=max([0,k]∩𝒮_d-i).
We can  calculate  X_i^(k+1) by induction on i via (<ref>)  for  0≤ i≤ d and  k∈𝒮_d-i. Recalling that 𝒮_0={0,…,n-1}, this allows us to simulate X_d^(1),…,X_d^(n). The generic dimension reduction algorithm for Markov chains estimation described in <ref> can thus be rewritten as follows: 
  *  Generate   N_k for 1≤ k≤ n-1 and calculate
the sets 𝒮_i, 0≤ i≤ d. 

  *  Set X_0^(1)=X_0.

  *  For i=0,…,d-1 and k∈𝒮_d-i-1, sample a copy U^(k+1)_d-i of Y_i and calculate  X_i+1^(k+1) via (<ref>). 
  *  Output the average  of g(X_d^(1)),…,g(X_d^(n)).

Step 3 generates one or several copies of X_i+1 for each copy of X_i, and so the above algorithm may be viewed as a splitting algorithm. Assume now that  N_k=N̅_k for k≥1.
 Step 1 can then be implemented via (<ref>).   In (<ref>),  m=0 if i=0 and, by (<ref>),  m=⌊ k/μ_d-i⌋μ_d-i if 1≤ i≤ d-1.  The expected running times of this algorithm and of the algorithm in <ref>  are  within a constant from each other as they are  both 
 proportional to the number of sampled copies of the U_i's. 

§ FURTHER NUMERICAL EXPERIMENTS


 §.§ Comparison with Quasi-Monte Carlo

Consider the G_t/D/1 queue with the same parameters as in <ref>. Table <ref> estimates E(X_d), and   Table <ref> gives VRFs in the estimation of  ℙ(X_d> z), for selected values of z. Our numerical results for the RDR and DDR algorithms are similar to those of <ref>.
Once again, for  the RDR and DDR algorithms, the variable Cost × Std^2 is roughly independent of d, and the variance reduction factors are roughly proportional to d. In contrast, for the QMC algorithm, the variable Cost × Std^2 is roughly proportional to d, and the variance reduction factors are roughly constant. The VRFs of the RDR and DDR algorithms in Table <ref> are, in general, greater than or equal to the corresponding VRFs in Table <ref>, which confirms the resiliency of these algorithms to discontinuities of g. In contrast, the VRFs of the QMC algorithm  in Table <ref> are lower than the corresponding VRFs  in Table <ref>. The RDR and DDR algorithms  outperform the QMC
algorithm. The VRFs of the QMC algorithm in Tables <ref> and <ref> are of the same order of magnitude as those obtained
by <cit.>, who
have reduced the variance in the simulation of a M/M/1 queue by a  factor
ranging between 5 and 10 via lattice rules. 
   
 

 §.§ G_t/D/1 queue with time-varying
amplitude

Consider a G_t/D/1 queue where  A_i has a Poisson distribution with time-varying rate

    λ_i=(1-1/ln (i+2))(0.75+0.5cos(π i/50)),
   for 1≤ i≤ d (recall that A_0=0).
Thus, up to a time-varying multiplicative factor,   λ_i  has
the same expression as in <ref>. Table <ref> estimates E(X_d), and   Table <ref> gives VRFs in the estimation of  ℙ(X_d> z), for selected values of z. Our numerical results are similar to those of <ref>, except
that there is a big difference in the means at different times.
 Once again, for  the RDR, DDR, and MLMC algorithms, the variable Cost × Std^2 is roughly independent of d, and the variance reduction factors are roughly proportional to d. The VRFs of the RDR and DDR algorithms in Table <ref>  are, in most cases, greater than or equal to the corresponding VRFs in Table <ref>, which confirms the resiliency of these algorithms to discontinuities of g. In contrast, the VRFs of the MLMC algorithm in Table <ref> are lower than the corresponding VRFs  in Table <ref>. The RDR algorithm outperforms the MLMC
algorithm by  a factor  ranging from 1 to 2  in Table <ref>,
and a factor ranging from 2 to 14  in Table <ref>.



 §.§ A multi-frequency G_t/D/1 queue

Consider a G_t/D/1 queue where  A_i has a Poisson distribution with time-varying rate

    λ_i=0.75+0.2cos(π i/50)+0.1cos(π i/5000)+0.05cos(π i/500000),
   for 1≤ i≤ d, with A_0=0.
Table <ref> estimates E(X_d). 
 Once again, for  the RDR, DDR, and MLMC algorithms, the variable Cost × Std^2 is roughly independent of d, and the variance reduction factors are roughly proportional to d.    The RDR algorithm outperforms the MLMC
algorithm by  about a factor  of 1.5.



 §.§ Comparison with a long-run  average estimator

In several examples, when d is large, the Markov chain (X_m), 0≤
m≤ d, has some notion
of stationarity, and so E(g(X_d))  can be estimated  via a suitable
long-run average. A drawback of such an estimator is that it is   biased.
We compare below a long-run average  estimator with the RDR and DDR algorithms. In Tables <ref> through <ref>, the RDR and DDR algorithms  use 10^9/d  samples, while the long-run algorithm uses 11×10^9/d  samples. Thus, for
each    of the three algorithms and each d,  the  total number of simulations of the U_is throughout the independent samples is roughly 11×10^9. The bias of the long-run average  estimator is calculated by taking the  difference
with the DDR estimator.
 

For the GARCH volatility example of <ref>, a natural long-run average  estimator for (X_d > z) is 

    1/d∑^d_i=11{X_i>z}.

Table <ref>  compares this estimator with the RDR and DDR estimators.
 The work-normalized variance of each of the three algorithms is roughly independent
of d. The work-normalized variance of the long-run average  estimator  is smaller than those of the RDR and DDR algorithms by about
a factor of 3 and 2, respectively.  The bias of the long-run average  estimator decreases
as  d increases, but   is
much larger than the standard deviations of the long-run average  and DDR estimators.


For the G_t/D/1 queue example of <ref> (resp. <ref>), the arrival rate is periodic (resp. almost periodic) with  period  100, and so a long-run average   estimator for E(X_d)
is 

    1/⌈ d/100⌉∑^⌈ d/100⌉-1_i=0X_d-i*100.

Tables <ref> and <ref> compares this estimator with the RDR and DDR estimators. The work-normalized variance of the long-run average  estimator is larger than those of the RDR and DDR estimators by 
a factor ranging from 1.5 and 2.2. The bias of the long-run algorithm
is not reported in Table <ref> because it is not statistically significant. In Table <ref>, however, the bias of the long-run average  estimator is
much larger than the standard deviations of the long-run average  and DDR estimators.





Consider now the  M_t/GI/1 queue example of <ref>, and assume that θ is an integer. In our experiments, we have set d=θ and X_i=W_i for 0≤ i≤ d,  and so a long-run average  estimator for 
 ℙ(W_θ> 1) is 

    1/⌈ d/100⌉∑^⌈ d/100⌉-1_i=01{X_d-i*100>1}.

Table <ref> compares this estimator with the RDR and DDR estimators. The work-normalized variance of the long-run average  estimator is larger than those of the RDR and DDR estimators by about
a factor of 1.3 and 2, respectively. Here again, the bias of the long-run algorithm
is not reported because it is not statistically significant.
 

In summary, for large d and Markov chains with periodic features, E(g(X_d)) can be  estimated via a suitable biased long-run average  estimator. The order of magnitude of the bias  depends on
the application and on the value of  d. The bias of the long-run average  estimator is  difficult to evaluate without using an alternative estimator, though. In the examples in this subsection, the work-normalized variances of the long-run average, RDR and DDR estimators have the same order of magnitude. In the  G_t/D/1 queue example of <ref>, however, the arrival rate is periodic with  period  10^6. This example does not seem to admit a suitable long-run  average estimator for the values of d listed in Table <ref>.