Subsampling

1 Subsampling Suppose Xi , i = 1, . . . , n is an i.i.d. sequence of random variables with distribution P . Let θ(P ) be some real-valued parameter of interest, and let θ̂n = θ̂n (X1 , . . . , Xn ) be some estimate of θ(P ). It is natural to assume in the i.i.d. setting that θ̂n is symmetric in its arguments, so we will make this assumption below, but it is not necessary. Consider the root Rn = Rn (X1 , . . . , Xn , θ(P )) = τn (θ̂n − θ(P )) , where τn is some normalizing sequence. Often, τn = √ n, but we do not wish to assume this. Let Jn (x, P ) denote the distribution of Rn ; that is, Jn (x, P ) = Pr{ τn (θ̂n − θ(P )) ≤ x} . We wish to estimate Jn (x, P ) so we can make inferences about θ(P ). For example, we would like to estimate quantiles of Jn (x, P ), so we can construct confidence sets for θ(P ). Unfortunately, we do not know P , and, as a result, we do not know Jn (x, P ). The bootstrap solved this problem simply by replacing the unknown P with an estimate P̂n . In the case of i.i.d. data, a typical choice of P̂n is the empirical distribution of the Xi , i = 1, . . . , n. For this approach to work, we essentially required that Jn (x, P ) when viewed as a function of P was continuous in a certain neighborhood of P . We will now explore an alternative to the bootstrap known as subsampling that does not impose this requirement. Subsampling is originally due to Politis and Romano (1994). In order to motivate the idea behind subsampling, consider the following thought experiment. Suppose for the time being that θ(P ) is known. Suppose that, instead of n i.i.d. observations from P , we had a very, very large number of i.i.d. observations from P . For concreteness, suppose Xi , i = 1, . . . , m is an i.i.d. sequence of random variables with distribution P with m = nk for some very big k. We could then estimate Jn (x, P ) by looking at the empirical distribution of τn (θ̂n (Xn(j−1)+1 , . . . , Xnj ) − θ(P )), j = 1, . . . , k . 1 This is an i.i.d. sequence of random variables with distribution Jn (x, P ). Therefore, by the Glivenko-Cantelli theorem, we know that this empirical distribution is a good estimate of Jn (x, P ), at least for large k. In fact, with a simple trick, we could show that it is even possible to improve upon this estimate by using all possible sets of data of size n from the m observations, not just those that are disjoint; that is, estimate Jn (x, P ) with the empirical distribution of the τn (θ̂n,j   m − θ(P )), j = 1, . . . , , n where θ̂n,j is the estimate of θ(P ) computed using the jth set of data of size n from the original m observations. In practice m = n, so, even if we knew θ(P ), this idea won’t work. The key idea behind subsampling is the following simple observation: replace n with some smaller number b that is much smaller than n. We would then expect τb (θ̂b,j   n − θ(P )), j = 1, . . . , , b where θ̂b,j is the estimate of θ(P ) computed using the jth set of data of size b from the original n observations, to be a good estimate of Jb (x, P ), at
least if nb is large. Of course, we are interested in Jn (x, P ), not Jb (x, P ). We therefore need some way to force Jn (x, P ) and Jb (x, P ) to be close to one another. To ensure this, it suffices to assume that Jn (x, P ) → J(x, P ). Therefore, Jb (x, P ) and Jn (x, P ) are both close to J(x, P ), and thus close to one another as well, at least for large b and n. In order to ensure that
both b and nb are large, at least asymptotically, it suffices to assume that b → ∞, but b/n → 0. This procedure is still not feasible because in practice we typically do not know θ(P ). But we can replace θ(P ) with θ̂n . This would cause no problems if τb (θn − θ(P )) = τb τn (θn − θ(P )) τn were small. Since τn (θn − θ(P )) = OP (1), it is enough to assume that 2 τb /τn → 0. Typically, this assumption will be implied by the assumption that b/n → 0, but it may not be, so we need to assume it separately. The next theorem formalizes the above discussion. Theorem 1.1 Let Xi , i = 1, . . . , n be an i.i.d. sequence of random variables with distribution P . Let θ(P ) be a real-valued parameter of interest, and let θ̂n = θ̂n (X1 , . . . , Xn ) be some estimate of θ(P ). Assume that θ̂n is symmetric in its arguments. Let Jn (x, P ) be the distribution of the root τn (θ̂n − θ(P )). Suppose Jn (x, P ) converges in distribution to J(x, P ). Let b = bn > 0 be a sequence of integers such that b → ∞, b/n → 0, and τb /τn → 0. Index by
i = 1, . . . , Nn = nb the different subsets of data of size b, and let θ̂b,i be θ̂b evaluated on the ith subset of data. Define 1 X I{τb (θ̂b,i − θ̂n ) ≤ x} . Ln (x) = Nn 1≤i≤Nn P (i) Ln (x) → J(x, P ) for all continuity points of J(x, P ). (ii) If J(x, P ) is continuous at c(1 − α) = J −1 (1 − α, P ), then ĉn (1 − α) = P L−1 n (1 − α) → c(1 − α, P ). Proof: (i) Let Un (x) = 1 Nn X I{τb (θ̂b,i − θ(P )) ≤ x} . 1≤i≤Nn Notice Un (x) only differs from Ln (x) in that θ̂n is replaced with θ(P ). Intuitively, we expect Un (x) and Ln (x) to be close under our assumptions, so P we will first show that Un (x) → J(x, P ) for all continuity points of J(x, P ). Note that Un (x) is an average of the Nn terms I{τb (θ̂b,i − θ(P )) ≤ x} and E[I{τb (θ̂b,i − θ(P )) ≤ x}] = Jb (x, P ) . P Therefore, we might expect Un (x) − Jb (x, P ) → 0. If this were true, then P Un (x) → J(x, P ) for all continuity points of J(x, P ), since Jb (x, P ) → J(x, P ) for all continuity points of J(x, P ) (this follows from the assumptions that b → ∞ and Jn (x, P ) converges in distribution to J(x, P )). 3 The object Un (x) is an example of a U -statistic, and the fact that Un (x)− P Jb (x, P ) → 0 actually follows from the theory of U -statistics, but we will try to give a more direct proof. To this end, first consider the behavior of 1 X Ūn (x) = I{τb (θ̂b (Xb(i−1)+1 , . . . , Xbk ) − θ(P )) ≤ x} , k 1≤i≤k where k = kn = bn/bc. Since this is an average of i.i.d. random variables, it P is easy to see that Ūn (x) − Jb (x, P ) → 0. Concretely, let  > 0 be given and use Chebychev’s inequality to write Pr{|Ūn (x) − Jb (x, P )| > } ≤ Var(Ūn (x)) . 2 But, Var(Ūn (x)) = Var(I{τb (θ̂b (Xb(i−1)+1 , . . . , Xbk ) − θ(P )) ≤ x}) →0. k P Thus, Ūn (x) − Jb (x, P ) → 0. By the same argument, we see that in order to prove that Un (x) − P Jb (x, P ) → 0, it is enough to prove that Var(Un (x)) → 0. Intuitively, we might expect the Var(Un (x)) ≤ Var(Ūn (x)). In order to formalize this idea, note that we may write Un (x) = E[Ūn (x)|X(1) , . . . , X(n) ] . To see this, note that X X k Ūn (x, Xπ(1) , . . . , Xπ(n) ) = kb!(n − b)! I{τn (θ̂b,i − θ(P )) ≤ x} 1≤i≤Nn π∈Sn   n = kb!(n − b)! Un (x) = kn!Un (x) , b where Sn is the set of all permutations π of 1, . . . , n. Therefore, 1 X Un (x) = Ūn (x, Xπ(1) , . . . , Xπ(n) ) = E[Un (x)|X(1) , . . . , X(n) ] . n! π∈Sn The last equality follows from the fact that since the data is i.i.d., conditional on X(1) , . . . , X(n) , each of the n! orderings are equally likely. 4 Now the desired result follows because Var(Un (x)) = E[(Un (x) − Jb (x, P ))2 ] = E[E[Ūn (x) − Jb (x, P )|X(1) , . . . , X(n) ]2 ] ≤ E[E[(Ūn (x) − Jb (x, P ))2 |X(1) , . . . , X(n) ]] = Var(Ūn (x)) . P We have therefore shown that Un (x) → J(x, P ) for all continuity points of J(x, P ). We now use this fact to establish the desired result, that is, P Ln (x) → J(x, P ) for all continuity points of J(x, P ) To this end, let x be a continuity point of J(x, P ) and note that Ln (x) = 1 Nn X I{τb (θ̂b,i − θ(P )) − τb (θ̂n − θ(P )) ≤ x} . 1≤i≤Nn For  > 0, let En = {|τb (θ̂n − θ(P ))| < } . Since I{τb (θ̂b,i − θ(P )) − τb (θ̂n − θ(P )) ≤ x} = I{τb (θ̂b,i − θ(P )) ≤ x + τb (θ̂n − θ(P ))} , it follows that when En is true, I{τb (θ̂b,i − θ(P )) − τb (θ̂n − θ(P )) ≤ x} ≤ I{τb (θ̂b,i − θ(P )) ≤ x + } I{τb (θ̂b,i − θ(P )) − τb (θ̂n − θ(P )) ≤ x} ≥ I{τb (θ̂b,i − θ(P )) ≤ x − } . Since τb /τn → 0 and Jn (x, P ) converges in distribution to J(x, P )En has probability tending to 1. Therefore, with probability tending to 1, Un (x − ) ≤ Ln (x) ≤ Un (x + ) . If  > 0 is such that x −  and x +  are continuity points of J(x, P ), then P Un (x − ) → J(x − , P ) P Un (x + ) → J(x + , P ) 5 Therefore, with probability tending to 1, for any δ > 0, we have that δ 2 δ Un (x + ) < J(x + , P ) + . 2 Un (x − ) > J(x − , P ) − Choose  > 0 so that x −  and x +  are continuity points of J(x, P ) and δ 2 δ J(x + , P ) + 2 J(x − , P ) − > J(x, P ) − δ < J(x, P ) + δ . This is possible because J(x, P ) is continuous at x. Putting this all together, we have that with probability tending to 1, J(x, P ) − δ < Ln (x) < J(x, P ) + δ . Since the choice of δ > 0 was arbitrary, the desired result follows. (ii) This follows from an argument almost identical to the proof of Lemma 2.1 in the lecture notes on the bootstrap. In practice, Nn is too large to actually compute Ln (x), so what one would do is randomly sample B of the Nn possible data sets of size b and just use B in place of Nn when computing Ln (x). Provided B = Bn → ∞, all the conclusions of the theorem remain valid. This approximation step is similar in spirit to approximating the bootstrap distribution Jn (x, P̂n ) using simulations from P̂n . It is in fact possible to show that P sup |Un (x) − Jb (x, P )| → 0 . x∈R Remarkably, this convergence is in fact even uniform in P over any set of distributions for the observed data! This fact can be used to analyze when the coverage probability of confidence sets constructed using subsampling converges to the nominal level not just for a fixed distribution of the observed data, but uniformly over a set of distributions for the observed data. 6 See Romano and Shaikh (2006) and Andrews and Guggenberger (2007) for details. Essentially, all we required was that Jn (x, P ) converged in distribution to a limit distribution J(x, P ), whereas for the bootstrap we required this and additionally that Jn (x, P ) was continuous in a certain sense. Showing continuity of Jn (x, P ) was very problem specific. We also saw an example where Jn (x, P ) → J(x, P ), but this continuity failed (e.g., the extreme order statistic). Subsampling would have no problems handling the extreme order statistic. Typically, when both the bootstrap and subsampling are valid, the bootstrap works better in the sense of higher-order asymptotics (see the lecture notes on the bootstrap), but subsampling is more generally valid. There is a variant of the bootstrap known as the m-out-of-n bootstrap. Instead of using Jn (x, P̂n ) to approximate Jn (x, P ), one uses Jm (x, P̂n ) where m is much smaller than n. If one assumes that m2 /n → 0, then all the conclusions of the theorem remain valid with Jm (x, P̂n ) in place of Ln (x). This follows because if m2 /n → 0, then (i) m/n → 0 and (ii) with probability tending to 1, the approximation to Jm (x, P̂n ) is the same as the approximation to Ln (x) because the probability of drawing all distinct observations tends to 1. To see this, note that this probability is simply equal to Since 1 − n(n − 1)(n − 2) · · · (n − b + 1) = nb i n Y (1 − 1≤i≤b−1 i ). n ≥ 1 − nb , we have that 2 Y 1≤i≤b−1 b i b b (1 − ) ≥ (1 − ) = (1 − n )b . n n b If b2 /n → 0, then for every  > 0 we have that b2 /n <  for all n sufficiently large. Therefore, (1 − b2 n b  ) > (1 − )b → exp(−) . b b 7 By choosing  > 0 sufficiently small, we see that the desired probability converges to 1. 8