Computer Intensive Methods: Appendix

R.C.H. Cheng School of Mathematics University of Southampton Southampton SO17 1BJ UK

June 2005

1 The Bootstrap

1.1 The Bootstrap Concept

To start with we shall just consider samples Y = (Y_{1}, Y_{2}, ..., Y_{n}) and statistics s(Y) calculated from them. The Y_{i} may themselves be vectors, and they are not necessarily independent, though often they will be. We shall assume initially however that the Y_{i} are mutually independent so that Y is a random sample. (The case of several samples is considered in Section (5) and the case when the Y_{i} are dependent will be discussed in Section (7).)

Bootstrapping is simply a numerical way of answering the question: ‘What is G(.), the probability distribution of a statistic s(Y), calculated from a sample Y = (Y_{1}, Y_{2}, ..., Y_{n})?’

If we knew F (.), the distribution of Y, the above question, ‘What is G(.)?’, is easily answered numerically. We simply generate a number, B, of independent samples of Y: y_{1}, y_{2}, ... y_{B}, and calculate the statistic s_{j} = s(y_{j}) from each sample. Then the empirical distribution function (EDF) of the sample s = (s_{1}, s_{2}, ..., s_{B}),

G_{n}(s | s) = \frac{# of s_{j} ≤ s}{B},

will converge pointwise with probability one to G(s). [Here and in what follows we shall use the notation G_{n}(.) to denote the EDF of a sample of size n drawn from the distribution G(.). Where the precise sample s used to construct G_{n} needs to be made explicit then we use the notation G_{n}(.| s) or G_{n}(s| s)] The basic sampling process is thus:

The Basic Sampling Process

For j = 1 to B

For i = 1 to n

Draw y_{ij} from F (.)
Next i
Calculate s_{j} = s(y_{j})

Next j

Form G_{n}(.| s) The problem is that we do not know F (.) usually. However we do have the EDF, F_{n}(.| y), formed from the sample y. The bootstrap idea is to use F_{n}(.| y) instead of F (.) in the above basic process. To obtain a sample we draw values, not from F (.), but from F_{n}(.| y). This is equivalent to drawing a sample of the same size n, with replacement, from the original set of y’s. We call such a sample a bootstrap sample to distinguish it from the original sample, and write it as y^{∗} = (y_{1}^{∗}, y_{2}^{∗}, ..., y_{n}^{∗}). As in the basic process, B such bootstrap samples are drawn, and the statistic is calculated from these each of these samples: s_{j}^{∗} = s(y_{j}^{∗}). The EDF, G_{n}(. | s^{∗}), of these bootstrap statistics s_{j}^{∗}, is our estimate of G(.). The bootstrap process is as follows

The Bootstrap Sampling Process

Given a random sample y = (y_{1}, y_{2}, ...y_{n}) from F (.)

Form the EDF F_{n}(.|y)

For j = 1 to B

For i = 1 to n

Draw y_{ij}^{∗} from F_{n}(.|y)
Next i
Calculate s_{j}^{∗} = s(y_{j}^{∗})

Next j

Form G_{n}(.|s^{∗}) The underlying hope and expectation is that the bootstrap process will under many conditions reproduce the behaviour of the original process.

1.2 Basic Method

We shall assume that the objective is to use the random sample y = {y_{1}, y_{2}, ..., y_{n}} to estimate a parameter η = η(F (.)) that describes some quantity of interest of the distribution of Y. Typical examples are the mean

η = ∫_{−∞}^{∞} ydF (y)

or some other moment of Y, or a quantile η = q(p), defined by

∫_{−∞}^{q(p)} dF (y) = p.

The statistic s(Y) can then be regarded as any appropriate quantity that we may care to use for estimating η. We shall emphasise this view by using the alternative but entirely equivalent notation η̂(Y) ≡ s(Y), with the ‘hat’ indicating an extimated quantity. Thus, when η is the mean, an obvious statistic is the sample mean

s(Y) = Ȳ .

It is important to realise that we do not have to use the sample version as the estimator, assuming we can calculate it at all. In the present example we might use a trimmed sample mean, or even the median rather than the sample mean. What we shall do is regard the sample version of the parameter as being the quantity of interest in the bootstrap process, i.e. the parameter of interest in the bootstrap process is

η^{∗} = η(F_{n}(.))

Thus we think of the bootstrap process as one where we are generating bootstrap samples y_{j}^{∗} from which bootstrap estimates s_{j}^{∗} = s(y_{j}^{∗}) are obtained that estimate the bootstrap parameter η^{∗}.

With this viewpoint, the bootstrap process is a prime example of the so-called plug-in approach, being a precise analogue of the original process with the only difference being that the known F_{n}(.) is plugged-in for, i.e. replaces, the unknown F (.). The bootstrap principle is that the bootstrap process reproduces the properties of the original basic process.

Useful quantities that are obtained from the bootstrap process are the sample mean and variance of the bootstrap estimates:

s̄^{∗} = \frac{1}{B} ∑_{j=1}^{B} s_{j}^{∗},

σ^{∗2} = \frac{1}{B − 1} ∑_{j=1}^{B}(s_{j}^{∗} − s̄^{∗})^{2}.

We also have an estimate of the bias of the mean of the bootstrap estimates: b^{∗} = s̄^{∗} − η^{∗}. A bootstrap bias adjusted estimate for η is thus

η̆ = η̂ − b^{∗}.

If the bootstrap expectation satisfies E^{∗}(s_{j}^{∗}) = η̂ (here E^{∗} denotes the conditional expectation E[s_{j}^{∗}| F_{n}(.|y)], where F_{n}(.) is treated as being the parent distribution), then s̄^{∗} → η̂ in probability as B → ∞, and we have

η̆ = η̂ − b^{∗} = η̂ − (s̄^{∗} − η^{∗}) → η̂ − (η̂ − η^{∗}) = η^{∗}.

Thus, in this case, adjusting the original estimator for bias using the bootstrap bias is equivalent to using the bootstrap parameter η^{∗} as the initial estimate.

1.3 The Double Bootstrap and Bias Correction

The bootstrap bias adjustment might not remove all the bias. At the expense of significantly increased computation, a second application of the bootstrap can be made. This is known as the double bootstrap, and consists of an outer bootstrap that is exactly the structure of the initial bootstrap, but where each step within this outer bootstrap is itself a (inner) bootstrap. The precise calculation is as follows.

The Double Bootstrap (for Bias Correction)

Given a random sample y = (y_{1}, y_{2}, ...y_{n}) from F (.)

Outer Bootstrap:

For j = 1 to B

For i = 1 to n

Draw y_{ij}^{∗} from F_{n}(.|y)

Next i

Let η_{j}^{∗∗} = η(F_{n}(.|y_{j}^{∗})) Inner Bootstrap: For k = 1 to B

For i = 1 to n
Draw y_{ijk}^{∗∗} from F_{n}(.|y_{j}^{∗})
Next i
Calculate s_{jk}^{∗∗} = s(y_{jk}^{∗∗})

Next k

Let s̄_{j·}^{∗∗} = \frac{1}{B} ∑_{k=1}^{B} s_{jk}^{∗∗} Let b_{j}^{∗∗} = s̄_{j·}^{∗∗} − η_{j}^{∗∗} Let η̆_{j}^{∗} = η̂_{j}^{∗} − b_{j}^{∗∗} End of Inner Bootstrap

Next j

Let ¯ η̆_{·}^{∗} = \frac{1}{B} ∑_{j=1}^{B} η̆_{j}^{∗}

Let b^{∗} = ¯ η̆_{·}^{∗} − η^{∗}

Let η̆ = η̂ − b^{∗}

End of Outer Bootstrap Whilst this can lead to a reduction in the bias, the effect on the variance of the estimator is not so clear. The main obvious penalty is the increase in computing cost. If a full numerical method is used to carry out the double bootstrap then the computation increases quadratically with B.

The full double bootstrap is not always necessary. It may be possible to avoid the inner bootstrap by using an estimator of η that is already bias adjusted. For example if it is thought that the original estimator might seriously biased then one might replace it immediately with η^{∗} and then use bootstrap bias adjustment on η^{∗} using the single bootstrap.

The single bootstrap yields an estimate of the variance of η̂^{∗}. When the bias correction is small this same variance will of course estimate the variance of η̆ as well. However when the bias correction is large then the double bootstrap has to be used

1.4 Parametric Bootstrap

We now consider parametric bootstrapping. Again we focus on a random sample y = (y_{1}, y_{2}, ..., y_{n}) drawn from a distribution F (y, θ), with the difference that the form of F is known but which depends on a vector θ of parameters that is unknown. We denote the unknown true value of θ by θ_{0}. It is then natural to carry out bootstrap sampling, not from the EDF F_{n}(y | y), but from the distribution F (y, θ̂), where θ̂ is an estimate of θ_{0}.

An obvious choice for θ̂ is the maximum likelihood (ML) estimator because it has an asymptotic distribution that is easy to calculate. We define the log-likelihood by

L(θ, y) = ∑_{i=1}^{n} log f(y_{i}; θ).  (1)

where y = (y_{1}, y_{2}, ..., y_{n}). The method of ML is to estimate θ_{0} by maximizing L, treated as a function of θ, subject to θ ∈ Θ, where Θ is some allowed region of the parameter space. for the remainder of the subsection θ̂ will represent the ML estimator.

Regularity conditions which ensure that standard asymptotic theory applies are given in Wilks (1962) or Cox and Hinkley (1974) for example. When they apply θ̂ can found by solving \frac{∂L(θ)}{∂θ} |θ̂ = 0. Moreover

−n^{−1} \frac{∂^{2}L(θ)}{∂θ^{2}} |_{θ_{0}} → I(θ_{0}),  (2)

a positive definite matrix, and the asymptotic distribution of θ̂, as n → ∞, is

n^{1/2}(θ̂ − θ_{0}) is N(0, V (θ_{0}))  (3)

where V (θ_{0}) = [I(θ_{0})]^{−1}. From (2) a reasonable approximation for I(θ_{0}) is

I(θ̂) ≃ −n^{−1} \frac{∂^{2}L(θ)}{∂θ^{2}} |θ̂ .  (4)

More generally suppose that η is a smooth function of θ, ie η = η(θ). The ML estimate of η is then simply

η̂ = η(θ̂),

where θ̂ is the ML estimate of θ. The asymptotic normality of η̂ follows from that of θ̂, assuming that η(θ) can be approximated by a linear Taylor series in θ. This method of obtaining the asymptotic behaviour of η̂ is called the delta method. We have

n^{1/2}(η̂ − η_{0}) ∼ N(0, σ^{2}(θ_{0}))  (5)

with

σ^{2}(θ_{0}) = g^{T} (θ_{0})V (θ_{0})g(θ_{0})  (6)

where

g(θ_{0}) = \frac{∂η(θ)}{∂θ}|_{θ_{0}} .  (7)

An approximation for the distribution of η̂ is given by (5) with θ̂ used in place of the unknown θ_{0}.

The above asymptotic theory is firmly established and much used in practice. For this reason parametric bootstrapping is not especially needed when parametric models are fitted to large samples. However, for small samples Monte-Carlo simulation is a well known method for examining the behaviour of estimators and fitted distributions in parametric problems . Such methods are actually parametric bootstrapping under a different name and as such have been known and used for a long time. In this sense parametric bootstrapping significantly predates the direct non-parametric bootstrap. Parametric bootstrapping lacks the computational novelty of direct sampling in that resampling is not involved. Perhaps it is because of this that the significance of parametric bootstrapping was not fully appreciated, and its theoretical importance not recognized until the landmark paper of Efron (1979).

If we view bootstrapping as a numerical sampling approach involving models that may include parameters, then its potential is significantly extended.

2 Percentiles and Confidence Intervals

2.1 Percentiles

We shall denote the estimate of a percentile q_{p}, with percentile probability p, by q̂_{p}(y).

The obvious choice for estimating a percentile non-parametrically is to use the corresponding percentile of the EDF. There is a certain ambiguity in doing this as the EDF is a step function. Thus all percentiles in the interval (\frac{k−1}{n}, \frac{k}{n}) are estimated by the observed order statistic y_{(k)}. Conversely the percentile \frac{k}{n} is estimated by any point y for which F_{n}(y) = k/n. Such ambiguities are removed if the EDF is smoothed. One way is to use a kernel estimator of the density but in most simulations this is perhaps unnecessarily elaborate.

The simplest smoothing procedure is to note that if Y_{(k)} is the kth order statistic (as opposed to the observed value y_{(k)}) then E[F (Y_{(k)})] = k/(n+1) and to use this value to estimate the value of F at the observed order statistic y_{(k)}. We can now simply interpolate between these estimated points (y_{(k)}, k/(n + 1)) for k = 1, 2, ..., n. The range can be extended to (0, 0) (using the line segment joining (0, 0) and (y_{(1)}, 1/(n + 1)) if Y is known to be a positive random variable. If we denote this smoothed version of F_{n}(.) by F̃_{n}(.) then an obvious estimator for q̂_{p}(y) is

q̂_{p}(y) = F̃_{n}^{−1}(p)

Estimating F (.) in the range (y_{(n)}, ∞), or equivalently q_{p} for p > n/(n+1) is not advisable unless the tail behaviour of F is known.

The bootstrap analogue of q̂_{p}(y) is obtained in the usual way. We draw y^{∗}, a bootstrap sample from y, and construct F̃_{n}(.| y^{∗}), the smoothed version of the bootstrap EDF. [Where there is no ambiguity we shall write F_{n}^{∗}(.) for F_{n}(.| y^{∗}) and F̃_{n}^{∗}(.) for F̃_{n}(.| y^{∗}).] We can now calculate q̂_{p}^{∗} as q̂_{p}^{∗} = F̃_{n}^{∗−1}(p). The bootstrap parameter here is η^{∗} = q̂_{p}(y).

In the parametric case, things are much easier. Quantiles are simply estimated from the fitted distribution F (., θ̂).

2.2 Confidence Intervals by Direct Bootstrapping

We now consider the construction of confidence intervals for a parameter of interest, η, and consider to what extent the corresponding bootstrap process. can be used to supply a confidence interval for this original parameter η. The key requirement is that the distribution of the bootstrap difference η̂^{∗} − η^{∗} should be close to that of η̂ − η. Roughly speaking we need

P^{∗}(\sqrt{n}(η̂^{∗} − η^{∗}) ≤ y) − P (\sqrt{n}(η̂ − η) ≤ y) → 0  (8)

for any y as n → ∞. (As before with E^{∗}, P^{∗} denotes the conditional probability evaluation treating F_{n}(.| y) as the parent distribution.) More precise statements of this convergence (8) and other results are given in Section (3).

We can proceed as follows. We generate B bootstrap samples y_{j}^{∗}, j = 1, 2, ..., B and corresponding bootstrap estimates η̂_{j}^{∗}, j = 1, 2, ..., B. Then we select an appropriate confidence level (1 − 2α) (we use 2α rather than α so that α corresponds to each of the two tail probabilities in two-sided intervals) and find values a^{∗} and b^{∗} for which

P^{∗}(a^{∗} ≤ η̂^{∗} − η^{∗} ≤ b^{∗}) ≃ 1 − 2α.

We then appeal to (8) to replace η̂^{∗} − η^{∗} by η̂ − η and obtain

P (a^{∗} ≤ η̂ − η ≤ b^{∗}) ≃ 1 − 2α

which on inversion gives the approximate (1 − 2α) confidence interval

η̂ − b^{∗} ≤ η ≤ η̂ − a^{∗}  (9)

for η.

One way of obtaining a^{∗} and b^{∗} is as follows. Let the bootstrap estimates of the α and (1 − α) quantiles obtained from the smoothed EDF of the η̂_{j}^{∗}’s be q̂_{α}(η̂^{∗}) and q̂_{1−α}(η̂^{∗}). These are estimates of the quantiles of the distribution of η̂^{∗}; so the corresponding quantiles of η̂^{∗} − η^{∗} are

a^{∗} = q̂_{α}(η̂^{∗}) − η^{∗} and b^{∗} = q̂_{1−α}(η̂^{∗}) − η^{∗}.  (10)

Substituting these into (9) and noting that η^{∗} = η̂ gives what is normally called the basic bootstrap confidence interval.

2η̂ − q̂_{1−α}(η̂^{∗}) ≤ η ≤ 2η̂ − q̂_{α}(η̂^{∗}).  (11a)

This confidence interval is also known as the hybrid bootstrap confidence interval as we have essentially replaced percentage points of the unknown original EDF by those of the smoothed bootstrap EDF. There is substantial evidence that this basic bootstrap can be significantly improved (in terms of giving coverage that is closer to the nominal level) if we use a pivotal quantity, like a studentized mean, rather than focusing just on the difference η̂ − η. The reason is that a pivotal quantity is less dependent (ideally not dependent at all) ) on the form of the original distribution. This would mean that the confidence interval is less influenced by any difference there may be in the distributions of η̂^{∗} − η^{∗} and η̂ − η. We consider this next.

2.3 Studentization

Suppose that σ̂^{2}(y), an estimate of the variance of η̂, can be calculated from the sample y. This will be the case for η̂ = ȳ say, when we have σ̂^{2}(y) = s^{2}, the sample variance.

As before suppose we can find a^{∗} and b^{∗} such that

P^{∗}(a^{∗} ≤ (η̂^{∗} − η^{∗})/σ̂^{∗} ≤ b^{∗}) ≃ 1 − 2α.

If we then appeal to a similar result to (8) and replace (η̂^{∗} − η^{∗})/σ̂^{∗} by (η̂ − η)/σ̂, we obtain

P (σ̂a^{∗} ≤ η̂ − η ≤ σ̂b^{∗}) ≃ 1 − 2α

which on inversion gives the approximate (1 − 2α) confidence interval

η̂ − σ̂b^{∗} ≤ η ≤ η̂ − σ̂a^{∗}.

We can now calculate a^{∗} and b^{∗} by drawing B bootstrap versions of z = (η̂ − η)/σ̂ : z_{j}^{∗} = (η̂_{j}^{∗} − η^{∗})/σ̂_{j}^{∗}, j = 1, 2, ..., B. Let q̂_{α}(z^{∗}) and q̂_{1−α}(z^{∗}) be the quantiles obtained from the EDF of these z_{j}^{∗}. The confidence interval for η is now

(η̂_{α}^{Stud}, η̂_{1−α}^{Stud}) = (η̂ − σ̂q̂_{1−α}(z^{∗}), η̂ − σ̂q̂_{α}(z^{∗})).  (12)

This is usually called the studentized bootstrap confidence interval.

Studentized bootstrap intervals can be readily used with the parametric model F (y, θ) of Section (1.4). Suppose that y is a random sample of size n drawn from the distribution F (y, θ_{0}), and that η_{0} = η(θ_{0}) is the quantity of interest. We can estimate η_{0} using η̂ = η(θ̂) where θ̂ is the MLE of θ_{0}.When n is large we can use the asymptotic approximation

n^{1/2}(η̂ − η_{0})/σ(θ̂) ∼ N(0, 1).

When n is not large it is worth employing bootstrapping. We draw B bootstrap samples y_{j}^{∗} j = 1, 2, ..., B from the fitted distribution F (y, θ̂). From each sample we obtain the bootstrap ML estimator θ̂_{j}^{∗} and bootstrap studentized quantity z_{j}^{∗} = n^{1/2}(η̂_{j}^{∗} − η^{∗})/σ(θ̂_{j}^{∗}). The quantiles q̂_{α}(z^{∗}) and q̂_{1−α}(z^{∗}) can be obtained from the EDF of these z_{j}^{∗} and the studentized interval (12) constructed.

In the non-parametric case an estimate of the variance of η̂ may not be immediately available. One possibility is to use the double bootstrap as given previously (where s = η̂) with the inner loop taking the form Inner Bootstrap: For k = 1 to B

For i = 1 to n

Draw y_{ijk}^{∗∗} from F_{n}(.|y_{j}^{∗})

Next i
Let s_{jk}^{∗∗} = s(y_{jk}^{∗∗})

Next k

Let η̂_{j}^{∗} = s(y_{j}^{∗})

Let s̄_{j·}^{∗∗} = \frac{1}{B} ∑_{k=1}^{B} s_{jk}^{∗∗}

Let v_{j}^{∗∗} = \frac{1}{B−1} ∑_{k=1}^{B}(s_{jk}^{∗∗} − s̄_{j·}^{∗∗})^{2} Let z_{j}^{∗} = n^{1/2}(η̂_{j}^{∗} − η^{∗})/(v_{j}^{∗∗})^{1/2} End of Inner Loop

If the double bootstrap is considered computationally too expensive, then an alternative using influence functions can be used provided our statistic is expressible as functional of the EDF, i.e. s(y) = η[F_{n}(.|y)]. Such statistics are termed statistical functions, and were introduced by von Mises (1947). We then assume that the relevant functional η = η[F (.)] has a linearised form (akin to a first order linear generalisation of a Taylor series): η(G(.)) ≃ η(F (.)) + ∫ L_{η}(y|F (.))dG(y). Here

L_{η}(y|F (.)) = \frac{∂η\{(1 − ε)F (.) + εH_{y}(.)\}}{∂ε} = lim_{ε→0} ε^{−1}[η\{(1 − ε)F (.) + εH_{y}(.)\} − η(F (.))], (13)

the derivative of η at F, is called the influence function of η, and H_{y}(x) is the Heaviside unit step function with jump from 0 to 1 at x = y. The sample approximation

l(y|y) = L_{η}(y|F_{n}(.|y))  (14)

is called the empirical influence function. An analogous argument to that used in the delta method yields the non-parametric analogue of (6) as

Var[η{F (.)}] = n^{−1} ∫ L_{η}^{2}{y|F (.)}dF (y),

with sample version

Var[η{F_{n}(.|y)}] = n^{−2} ∑_{j=1}^{n} l^{2}(y_{j}|y),

where the l(y_{j}|y) are the empirical function values evaluated at the observed values y_{j}. In practice these values have to be evaluated numerically using (13) with, typically, ε = (100n)^{−1}.

The ease with which the empirical influence function can be evaluated is rather problem dependent, and the potential complexity of this calculation detracts from this approach. One simple approximation is the jacknife approximation

l_{jack}(y_{j}|y) = (n − 1)(η̂(y) − η̂(y_{∖j}))

where η̂(y_{∖j}) is the estimate calculated from the sample y but with the observation y_{j} omitted. (See Davison and Hinkley 1997, problem 2.18.)

2.4 Percentile Methods

Suppose now that the distribution of η̂ is symmetric, or more generally, there is some transformation w = h(η̂) for which the distribution is symmetric. Examination of the form of the bootstrap percentiles (10), (see for example Hjorth, 1994, Section 6.6.2) shows that they do not depend on the explicit form of h(.) at all. Instead we find that they can, with a change of sign, be swapped, ie a^{∗} is replaced by −b^{∗} and b^{∗} by −a^{∗}. Then (9) becomes

(η̂_{α}^{P er}, η̂_{1−α}^{P er}) = (q̂_{α}(η̂^{∗}), q̂_{1−α}(η̂^{∗})).  (15)

This (1 − 2α) confidence interval is known as the bootstrap percentile interval. This confidence interval is easy to implement but can be improved. We give two variants called the bias corrected (BC) and accelerated bias corrected (BCa) percentile methods.

The basic idea in the BCa method is that there is a studentized transformation g for which

\frac{g(η̂) − g(η)}{1 + ag(η)} + b ∼ N(0, 1).  (16)

The BC method is simply the special case where a = 0.

If we calculate a one sided confidence interval for η under the assumption (16) we find that if we set

β_{1} = Φ (b + \frac{b + z_{α}}{1 − a(b + z_{α})})

where Φ(.) is the standard normal distribution function, then with confidence (1 − α)

q̂_{β_{1}}(η̂^{∗}) < η

where q̂_{β}(η̂^{∗}) is the bootstrap percentile of the η̂^{∗} with probability β. In a similar way a two sided 1 − 2α confidence interval is

q̂_{β_{1}}(η̂^{∗}) < η < q̂_{β_{2}}(η̂^{∗})

where

β_{2} = Φ (b + \frac{b + z_{1−α}}{1 − a(b + z_{1−α})}) .

The bias parameter b is obtained as

b = Φ^{−1}{F̆_{n}(η̂ | η̂^{∗})}

where F̃_{n}(η̂ | η̂^{∗}) is the smoothed EDF of the bootstrap sample η̂^{∗} evaluated at the original estimated value η̂. If we now set a = 0 this gives the BC method. We denote the resulting two-sided version of the confidence interval by (η̂_{α}^{BC}, η̂_{1−α}^{BC}).

The acceleration parameter a is arguably less easily calculated. Efron (1987) suggests the approximation (which is one-sixth the estimated standard skewness of the linear approximation to η) :

a = \frac{1}{6} \frac{∑_{i=1}^{n} l^{3}(y_{i}|y)}{[∑_{i=1}^{n} l^{2}(y_{i}|y)]^{3/2}}

where l(y_{i}|y) are the values of the empirical influence function (14) evaluated at the observations y_{i}. We denote the two-sided version of the confidence interval given by this BCa method by (η̂_{α}^{BCa}, η̂_{1−α}^{BCa}).

3 Theory

Like many statistical methods, understanding of the practical usefulness of the bootstrap as well as its limitations has been built up gradually with experience through applications. This practical experience has been underpinned by a growing asymptotic theory which provides a basis for understanding when bootstrapping will work and when it will not. A truly general theory rapidly becomes very technical and is still incomplete. A detailed treatment is not attempted here, but some of the more accessible and useful results are summarized.

Bootstrapping relies on sampling from the EDF F_{n}(.|y) reproducing the behaviour of sampling from the original distribution F (.). We therefore need convergence of F_{n}(.|y) to F (.). This is confirmed by the Glivenko-Cantelli theorem which guarantees strong convergence, i.e.

sup_{y} |F_{n}(y|y) − F (y)| → 0 with probability 1

as n → ∞. Though reassuring this does not throw any direct light on the bootstrap process.

To investigate the bootstrap process we look instead at the equally well-known result that the process

Z_{n}(y) = \sqrt{n}(F_{n}(y) − F (y))

converges in probability to the Brownian bridge B(F (y)) as n → ∞. [A Brownian bridge B(t), 0 < t < 1, is a Gaussian process for which B(0) = B(1) = 0, E[B(u)] = 0 and E[B(u)B(v)] = u(1 − v) for 0 < u < v < 1.] Bickel and Freedman (1981) give the following bootstrap version of this. Theorem Let Y_{1}, Y_{2}, ... be a sequence of independent observations from F (.), with Y_{n} comprising the first n values. Then for almost all such sequences (i.e. with probability 1) the bootstrap process

Z_{n}^{∗}(y) = \sqrt{n}(F_{n}^{∗}(y|Y_{n}) − F_{n}(y))

converges in probability to the Brownian bridge B(F (y)) as n → ∞.

An interesting immediate application of this result to the calculation of confidence bands for F (y) is given by Bickel and Freedman. We set a confidence level 1 − α and, from the bootstrap process, we then select c_{n}(α) so that, as n → ∞,

P^{∗}{sup_{y} |Z_{n}^{∗}(y)| ≤ c_{n}(α)} → 1 − α.

Then, as Z_{n}^{∗}(y) and Z_{n}(y) converge to the same process B(F (y)), we have also that

P^{∗}{sup_{y} |Z_{n}(y)| ≤ c_{n}(α)} → 1 − α.

Thus an asymptotically correct (1 − α) confidence band for F (y) is

F_{n}(y) ± \frac{c_{n}(α)}{\sqrt{n}}.

We now consider the distribution of an estimator η̂ of a parameter η. We need the distributions of the bootstrap estimator η̂^{∗} and of the difference η̂^{∗} − η^{∗} to reproduce respectively the distributions of η̂ (≡ η^{∗}) and of η̂ − η. The most general results obtained to date assume that the statistic can be expressed as s = η(F_{n}) where both the EDF F_{n} and the underlying distribution F belong to a space F of distributions. Moreover such results require s to be differentiable in such a way that a linear approximation can be formed with, as first order terms, the previously described influence functions, i.e.

η(F_{n}) = η(F) + n^{−1} ∑_{j=1}^{n} l(y_{j}|y) + o_{p}(n^{−1/2}).

Then it is often possible to establish the analogous result for the bootstrap version

η(F_{n}^{∗}) = η(F) + n^{−1} ∑_{j=1}^{n} l(y_{j}^{∗}|y^{∗}) + o_{p}(n^{−1/2}),

and, provided Var[η(F)] is finite, then, in probability, the sequence y_{1}, y_{2}, ... is such that

sup_{s} |P^{∗}(η(F_{n}^{∗}) ≤ s) − P (η(F_{n}) ≤ s)| → 0

as n → ∞. Here we have written P^{∗}(η(F_{n}^{∗}) ≤ s) for P (η(F_{n}^{∗}) ≤ s | F_{n}). If, in addition, s is continuously differentiable, then we can usually strengthen the result to almost sure convergence.

A detailed discussion is given by Shao and Tu (1995) who give examples of η̂(F) that satisfy such conditions. Convergence has been studied in detail for the particular case of means and percentiles. In this case a more accessible approach using the Berry-Esséen inequality is possible, and this we discuss in the next section.

3.1 Convergence Rates

Efron (1979) considered the finite case where Y takes just a finite set of values which we can take to be {1, 2, ..., m} with probabilities p =(p_{1}, p_{2}, ..., p_{m})^{T} . We shall not discuss this particular case, a good summary for which is provided by Hjorth (1994).

We shall instead focus on the continuous case where η = μ ≡ E(Y), with η̂ = Ȳ_{n} = ∑_{i=1}^{n} Y_{i}/n. We also write V ar(Y) = σ^{2}, ρ = E(|Y |^{3}). The well-known Berry-Esséen theorem (see Serfling, 1980, for example) then provides a powerful result for investigating bootstrap convergence. To avoid long formulas we shall write:

H_{n}(y) = P (\sqrt{n}(Ȳ_{n} − μ) ≤ y)

for the (unknown) CDF of interest and

H_{n}^{∗}(y) = P^{∗}(\sqrt{n}(Ȳ_{n}^{∗} − Ȳ_{n}) ≤ y)

for the bootstrap version of this. We write

H̃_{n}(y) = P (\sqrt{n}\frac{(Ȳ_{n} − μ)}{σ} ≤ y)

for the (unknown) standardized CDF of interest,

H̆_{n}(y) = P (\sqrt{n}\frac{(Ȳ_{n} − μ)}{σ̂_{n}} ≤ y)

for the studentized CDF, and

H̃_{n}^{∗}(y) = P^{∗}(\sqrt{n}\frac{(Ȳ_{n}^{∗} − Ȳ_{n})}{σ̂_{n}} ≤ y)

for its bootstrap studentized version. We write also Φ(y) for the CDF of the standard normal distribution.

The Berry-Esséen theorem states that if Y_{1}, Y_{2}, ... are independent identically distributed random variables and ρ = E[|Y_{i}^{3}|] < ∞, then, for all n

sup_{y} |H_{n}(y) − Φ(\frac{y}{σ})| < K \frac{ρ}{σ^{3}\sqrt{n}}.

The value of the constant given by Berry and Esséen, K = \frac{33}{4}, has been improved and reduced to .7975 (van Beeck,1972). Thus the theorem gives a bound on the effectiveness of the normal distribution in approximating H_{n}(.).

A typical application of the theorem is in proving the following. Theorem If ρ = E[|Y_{i}^{3}|] < ∞ then, as n → ∞,

sup_{y} |H_{n}^{∗}(y) − H_{n}(y)| → 0  (17)

for almost all sequences y_{1}, y_{2}, ...of independent observations drawn from F (.). Proof Outline. By the Berry-Esséen theorem we have

|P (H_{n}(y) − Φ(\frac{y}{σ})| < K \frac{ρ}{σ^{3}\sqrt{n}}

and

|P^{∗}(H_{n}^{∗}(y) − Φ(\frac{y}{σ̂_{n}})| < K \begin{matrix}ρ & ∗ n \\ σ̂ & 3 n & \sqrt{n} \end{matrix}.

By the strong law of large numbers μ^{∗} ≡ ȳ_{n} → μ, σ̂_{n} → σ, ρ_{n}^{∗} → ρ all with probability 1. This shows that the two probabilities in (17) converge to the same normal distribution from which the result follows easily. □

In fact the above result does not depend on ρ < ∞; the condition σ^{2} < ∞ is actually enough (see Bickel and Freedman, 1981 or Singh 1981.) However the condition ρ < ∞ and explicit use of the Berry-Esséen theorem is needed to establish the rate of convergence more precisely. In the remainder of this section we summarize results obtained by Singh (1981) and Bickel and Freedman (1981).

We need one further definition. A random variable has distribution F (.) that is lattice if there is zero probability of it taking values outside the discrete points y_{j} = c + jd j = 0, 1, 2, ... where c and d are constants. Theorem For almost all sequences Y_{1}, Y_{2}, ...of independent observations drawn from F (.):

(i) If E(Y^{4}) < ∞, then

lim sup_{n→∞}[\begin{matrix}√ & \overline{n} \\ √ & \overline{log log n} \end{matrix} sup_{y} |H_{n}^{∗}(y) − H_{n}(y)|] = \frac{\sqrt{V ar[(Y − μ)^{2}]}}{2σ^{2}\sqrt{πe}}

(ii) If E[|Y |^{3}] < ∞ and F (.) is lattice then

lim sup_{n→∞}[\sqrt{n} sup_{y} |H̃_{n}^{∗}(y) − H̃_{n}(y)|] = \frac{d}{\sqrt{2π}σ}

(iii) If E[|Y |^{3}] < ∞ and F (.) is non-lattice then

\sqrt{n} sup_{y} |H̃_{n}^{∗}(y) − H̃_{n}(y)| → 0

The result (i) shows that when E(Y^{4}) < ∞ then the convergence rate of H_{n}^{∗} is O(\sqrt{log log n/n}). This is actually the same as that for approximating \sqrt{n}(Ȳ_{n} − μ)/σ̂_{n}, the studentized original Ȳ by Z, a standard normal, i.e. approximating H_{n}(y) by Φ(y/σ̂_{n}). So in this case we do no better using the bootstrap than standard normal theory.

Results (ii) and (iii) show that when E[|Y |^{3}] < ∞, then the convergence rate of H̃_{n}^{∗}(.) to H̃_{n}(.) is at least O(n^{−1/2}) and is o(n^{−1/2}) when F is non-lattice. Thus bootstrapping does better than the usual normal approximation in this case. This is the key theoretical result that underpins much of bootstrapping.

We consider now quantile estimation. Perhaps the most informative general result is that due to Singh (1981). Let η = q(p) = F^{−1}(p); where the conditions of the theorem ensure that the inverse is uniquely defined. We can estimate η from the EDF using η̂ = F_{n}^{−1}(p) = sup{y : F_{n}(y) ≤ p\} and define η̂^{∗} = F_{n}^{∗−1}(p) from y^{∗} as usual. We take H_{n} = P (\sqrt{n}(η̂ − η) ≤ y) and H_{n}^{∗} = P^{∗}(\sqrt{n}(η̂^{∗} − η̂) ≤ y) Theorem If F has a bounded second derivative in a neighborhood of η and f(η) = dF/dy|_{η} > 0, then

lim sup_{n→∞} \frac{n^{1/4} sup_{y} |H_{n}^{∗}(y) − H_{n}(y)|}{\sqrt{log log n}} = c_{p,F} with probability 1

This result shows that H_{n}^{∗}(y) converges at rate O(n^{−1/4} log log n) to H_{n}(y). Now for a percentile estimator we have from the Berry-Esséen theorem

sup_{y} |H_{n}(y) − Φ(\frac{y}{τ})| = O(n^{−1/2})

where τ = \sqrt{p(1 − p)}/f(η). Thus if we were to use a normal approximation for H_{n}(y) we would use Φ(y/τ̂) where τ̂ is an estimate of τ. Whether this is better than H_{n}^{∗}(y) thus depends on the converence rate of τ̂, and the matter is not a clear one.

3.2 Asymptotic Accuracy of EDF’s

Edgeworth expansions are asymptotic series aimed at improving the normal approximation by introducing additional terms that try to correct for effects of skewness, kurtosis and higher moment which slow the rate of convergence to normality.

For fixed n, asymptotic series usually diverge as more and more terms are included. However for a fixed number of terms, k say, the series converges as n → ∞. The rate of convergence is usually of smaller order of than the last included term. We shall only consider the special case η̂ = Ȳ_{n}. Here, the general Edgeworth expansion is a power series in n^{−1/2} and has the form:

H̃_{n}(y) = Φ(y) + n^{−1/2}p^{(1)}(y)ϕ(y) + ... + n^{−k/2}p^{(k)}(y)ϕ(y) + o(n^{−k/2})  (18)

where ϕ(y) = (2π)^{−1/2} exp(−y^{2}/2) is the standard normal density, and p_{j} is a polynomial of degree 3j − 1. We have explicitly

p^{(1)}(y) = −\frac{1}{6}κ_{3}(y^{2} − 1)

and

p^{(2)}(y) = −\{\frac{1}{24}κ_{4}(y^{2} − 3) + \frac{1}{72}κ_{3}^{2}(y^{4} − 10y^{2} + 15)\}

where κ_{3} and κ_{4} are the skewness and kurtosis of F (.). The term involving p^{(1)} corrects for the main effect of skewness, the term involving p^{(2)} corrects for the main effect of kurtosis and for the secondary effect of skewness.

Often the remainder term o(n^{−k/2}) can be replaced by O(n^{−(k+1)/2}) when the Edgeworth expansion is said to be (k + 1)th order accurate. Usually inclusion of more than one or two correction terms becomes counter-productive as the coefficients associated with the powers of n^{−1/2} rapidly become large with k.

When E(|Y |^{3}) < ∞ and Y is non-lattice then one-term Edgeworth expansions for both H_{n}(y) and H_{n}^{∗}(y) exist and a comparison (see Shao and Tu, 1995, Section 3.3.3) shows that H_{n}^{∗}(y) has smaller asymptotic mean square error than Φ(y/σ̂_{n}) unless the skewness is zero. In fact comparison of H_{n}^{∗}(y) and the one-term Edgeworth expansion estimator

H_{n}^{EDG}(y) = Φ(z) + n^{−1/2}p_{n}^{(1)}(z | y_{n})ϕ(z)

where z = σ̂_{n}^{−1}y, and p_{n}^{(1)}(z | y_{n}) is the polynomial p^{(1)}(z) with the moments of F (.) replaced by the sample moments calculated from y_{n}, shows that both have the same asymptotic mean square error.

If we considered studentized versions the bootstrap does even better. Under appropriate moment conditions both H̆_{n}(y) and H̃_{n}^{∗}(y) have two-term Edgeworth expansions and a comparison of these shows that H̃_{n}^{∗}(y) has a smaller asymptotic mean square error than the corresponding one-term Edgeworth estimator, though they have the same convergence rate.

The normal, bootstrap and one-term Edgeworth expansion estimators of the standardized distribution have been compared by Hall (1988) using the criterion of asymptotic relative error.

When E(|Y |^{3}) < ∞ then the bootstrap does better than the normal approximation, however it may or may not do better than the one term Edgeworth expansion (see Shao and Tu, 1995).

When E(|Y |^{3}) = ∞, the situation is more complicated and depends on the tail behaviour of F (.). When the tail is thin the bootstrap can be worse than the normal approximation.

In estimating tail behaviour the bootstrap is comparable to the one term Edgeworth expansion except in the extreme of the tail.

3.3 Asymptotic Accuracy of Confidence Intervals

The above analysis focuses on distribution functions, and does not give the whole picture. It is helpful to consider also the coverage accuracy of confidence intervals. We shall write the basic confidence limit that we seek as η̂_{α} defined by

Pr(η ≤ η̂_{α}) = α

and the normal, basic, percentile, studentized bootstrap, BC and BCa approximations as η̂_{α}^{Norm}, η̂_{α}^{Boot}, η̂_{α,}^{P er} η̂_{α}^{Stu}, η̂_{α}^{BC} and η̂_{α}^{BCa} respectively. We summarise the results given in Shao and Tu (1995). These apply when the parameter of interest is a smooth function of the mean, η = η(ȳ_{n}). Then an analysis analogous to that used for EDF’s can be carried out, but now relying on an expansion of the quantile, that is the inverse of Edgeworth series, called the Cornish-Fisher expansion.

Let Φ(z_{α}) = α. Then, under appropriate moment conditions η̂_{α} has an expansion of the form

q(α) = z_{α} + n^{−1/2}q^{(1)}(z_{α}) + n^{−1}q^{(2)}(z_{α}) + o(n^{−1})

where comparison with the Edgeworth exapnasion shows that

q^{(1)}(y) = −p^{(1)}(y)

and

q^{(2)}(y) = p^{(1)}(y)p^{(1)′}(y) − \frac{1}{2}p^{(1)}(y)^{2} − p^{(2)}(y).

Under appropriate conditions we find that analogous expansions exist for the quantile approximations listed above. We find in particular that

η̂_{α}^{Boot} − η̂_{α} = O_{p}(n^{−1})

and in general

Pr(η̂_{α}^{Boot} ≤ η) = 1 − α + O(n^{−1/2}).

However the two tailed version is second-order accurate:

Pr(η̂_{α}^{Boot} ≤ η ≤ η̂_{1−α}^{Boot}) = 1 − 2α + O(n^{−1}).

For symmetric distributions, these results apply to the percentile approximation as well, for example:

η̂_{α}^{P er} − η̂_{α} = O_{p}(n^{−1}).

The bootstrap BC method turns out to perform no better than the percentile limit, in terms of convergence rate, but the constant factor is smaller so it is marginally to be preferred.

Studentization is definitely better with

η̂_{α}^{Stu} − η̂_{α} = O_{p}(n^{−3/2}) and η̂_{α}^{BCa} − η̂_{α} = O_{p}(n^{−3/2})

and

Pr(η̂_{α}^{Stu} ≤ η) = 1 − α + O(n^{−1}) and Pr(η̂_{α}^{BCa} ≤ η) = 1 − α + O(n^{−1}).

It follows that the two-sided intervals for both limits are also both second-order accurate.

Finally we note that

η̂_{α}^{Boot} − η̂_{α}^{Norm} = O_{p}(n^{−1})

so that the usual normal approximation and the basic bootstrap behave similarly.

3.4 Failure of Bootstrapping

It should be clear from the previous three subsections that, for bootstrapping to work well, regularity conditions are required on both the distribution F and also on the statistic of interest.

More explicitly, bootstrapping is sensitive to the tail behaviour of F; convergence of H_{n}^{∗} usually requires moment conditions on F that are more stringent than those needed for convergence of H_{n}.

Also the statistic s(y) has to be suitably smooth in an appropriate sense.

Finally it is possible for convergence of the bootstrap to be sensitive to the method used in carrying out the bootstrapping.

An example of the first situation is inconsistency of s(y), when it is simply the variance, even when the asymptotic variance is finite (see Ghosh et al., 1984). This can occur if F (.) is fat-tailed and does not have appropriate moments. The problem then arises because the bootstrap s(y^{∗}) can take exceptionally large values.

An example of the second situation is where y is a random sample, from U(0, b) say, and we wish to consider y_{(n)}, the largest order statistic. Then a natural statistic to consider is s(y) = n(b − y_{(n)})/b, which has a limiting standard exponential distribution as n → ∞. The bootstrap version is then s(y^{∗} | y) = n(y_{(n)} − y_{(n)}^{∗})/y_{(n)}. But

P^{∗}(s(y^{∗} | y) = 0) = P^{∗}(y_{(n)}^{∗} = y_{(n)}) = 1 − (1 − n^{−1})^{n} → 1 − e^{−1}.

Thus H_{n}^{∗} does not tend to H_{n} as n → ∞. This result applies to any given order statistic y_{(n−k)} where k is fixed. However the problem does not arise for a given quantile y_{(pn)} where p is fixed with 0 < p < 1.

We shall see an example of the last situation, where convergence depends not on the problem but on the bootstraaping method employed, when we consider model selection in Section (4.6).

4 Monte-Carlo/Simulation Models

4.1 Direct Models

Monte-Carlo or simulation models are direct models in the sense that they attempt to mimic the behaviour of an actual system by simulating the different objects of a system and the (typically) dynamic relationships between them. (The two names are interchangeable, but Monte-Carlo is the name more favoured by statisticians, whereas simulation is more used by the operational researchers.) A simple example is a single server queue we try to capture the actual arrival patterns of customers and how they are then processed by the server. The quantities that we try to analyse are summary characteristics such as average queue length and waiting times. We can therefore think of a set of n simulation runs as yielding observations

y_{j} = y(u_{j}, v_{j}, θ̂(w), x_{j}) j = 1, 2, ..., n.  (19)

where the y_{j} depend on a number of quantities that we now explain.

Firstly u_{j} denotes the stream of uniformly distributed U(0, 1) random numbers used in the jth run. Typically the uniforms are not used directly, but are transformed into random variables drawn from distributions other than the uniform.

Next comes v_{j}. This represents a sequence of inputs that are random, but that has been generated independently of the simulation model. A typical instance is where v_{j} comprises sampled real observations taken from some real system, separate from the system being modelled, but on which y_{j} depends. Such a sampled real process is sometimes called a trace. We shall simply think of v_{j} as being just a sample of observations.

In addition there may be further quantities which may affect y. These are denoted by x and θ. There is no essential difference between x and θ in the way that they influence y. They are simply variables on which y depends. However we make a distinction in supposing that x are decision variables, i.e. quantities that can be selected by the simulator. However we also include in x parameters whose values are known, and that might affect y, but which are not selectable by the simulator. The reason for adopting this convention is that it then allows us to assume θ to be those parameter values, whose values are not known, and so have to be estimated. We will therefore assume that in addition to v, there exists input data w that is used exclusively to estimate θ and it is the estimated values, θ̂(w), that are used in the simulation runs.

A simple example is a simulation model of a multiserver queue, where y is the observed average queue length over a given period, v might be a previously observed set of interarrival times, x might be the (scalar) number of servers and θ the service rates of servers. Here we treat x as being selectable by the simulator so that it is a design variable, θ may not be known and have to be estimated from available sampled real service times, w.

Resampling might already feature in carrying out the runs yielding the y_{j} of (19). For example variability in the v will contribute to the variability of y. We can therefore automatically build this effect into the simulation runs simply by using resampling from an initial given sequence v to construct resampled sequences v_{j}^{∗} for use in the actual runs.

Similarly in the case (19), because θ̂(w) is estimated, we can allow for its randomness affecting the distribution of y by using independent values of θ_{j} in each simulation run. When θ̂(w) is the ML estimate we can draw sample values of θ :

θ_{j}^{∗} ∼ N(θ̂, I^{−1}(θ̂))  (20)

from the asymptotic normal distribution (3) with I(θ̂) calculated from (4). An alternative is to use produce bootstrap samples w_{j}^{∗} from w, and to use the bootstrap estimate

θ_{j}^{∗} = θ̂(w_{j}^{∗})  (21)

for θ in the jth run.

4.2 Metamodels

We now consider a very different approach where the behaviour of a system is represented by a metamodel. This is simply where we regard the behaviour of y as representable by a statistical regression model. The runs are then not represented as (19) but instead as

y_{j} = η(x_{j}; β) + ε_{j} j = 1, 2, ..., n.

where η(x_{j}; β) is called the regression function. Its form is usually known, though in the problem of model selection, one has to choose from a number of possible competing models. The errors ε_{j} are assumed to be independent with E(ε_{i}) = 0. We shall consider mainly the case where all the error variances are equal V ar(ε_{i}) = σ^{2}, but will consider the heteroscedastic case too. The β are unknown coefficients that have to be estimated. We consider linear metamodels first.

4.3 Linear Metamodels

In the linear metamodel, observations are assumed to take the form

Y_{j} = x_{j}^{T} β + ε_{j} j = 1, 2, ..., n  (22)

where the errors are assumed to be independent and identically distributed with E(ε_{i}) = 0 and V ar(ε_{i}) = σ^{2}, and the β are unknown and have to be estimated. Let

X = (x_{1}, ..., x_{n})^{T} , Y = (Y_{1}, ..., Y_{n})^{T}

and

β̂ = (X^{T} X)^{−1}X^{T} Y

be the usual least squares estimate (equivalent to maximum likelihood when the ε_{i} are normally distributed). The fitted values are

Ŷ = HY  (23)

where

H = X(X^{T} X)^{−1}X^{T}  (24)

is the well-known known ’hat’ matrix.

The raw residuals are r = (I − H)Y. One way of resampling is to sample from these residuals as if they were the unknown errors. This gives the bootstrap sample

Y_{j}^{∗} = x_{j}^{T} β̂ + r_{j}^{∗} j = 1, 2, ..., n.  (25)

We shall call this residual sampling. However we have that E(r) = 0 and V ar(r) = σ^{2}(I − H). Thus a better bootstrapping method is to sample from the adjusted residuals

e_{j} = r_{j}/(1 − h_{jj})^{1/2}  (26)

where h_{jj} is the jth diagonal entry of H (h_{jj} commonly known as the leverage of the jth point). This form of model-based resampling yields bootstrap samples of form:

Y_{j}^{∗} = x_{j}^{T} β̂ + e_{j}^{∗} j = 1, 2, ..., n.  (27a)

We shall call this adjusted residual sampling.

An alternative is to use case sampling and to sample from the (Y_{j}, x_{j}) pairs directly. The main problem with this method is that it introduces extra variation into the data. In the conventional situation X is regarded as fixed, but using case sampling replaces this by a bootstrap X^{∗} which will be variable. This is not so much a problem when the number of parameters, p say, is small compared with n.

A partial correction can often be made. For example suppose we are interested in the distribution of c^{T} β̂ where c is a constant vector. Then V ar[c^{T} (β̂ − β)] = σ^{2}c^{T} (X^{T} X)^{−1}c. We can then study the behaviour of

c^{T} (β̂ − β)/(c^{T} (X^{T} X)^{−1}c).

indirectly by considering the bootstrap version

c^{T} (β̂^{∗} − β^{∗})/(c^{T} (X^{∗T} X^{∗})^{−1}c)

(where β^{∗} = β̂).

A common variation is when the errors are heteroscedastic so that their variances are unequal. In this case resampling can be carried out as

Y_{j}^{∗} = x_{j}^{T} β̂ + e_{j}t_{j}^{∗} j = 1, 2, ..., n

with e_{i} as defined in (26) and with t_{j}^{∗} independent and random with E(t_{j}^{∗}) = 0, V ar(t_{j}^{∗}) = 1. Wu (1986) suggests sampling the t_{j}^{∗} with replacement from the set

(r_{j} − r̄)/[n^{−1} ∑_{i=1}^{n}(r_{i} − r̄)^{2}]^{1/2} j = 1, 2, ..., n.

4.4 NonLinear Metamodels

The previous ideas apply with little change to nonlinear models. For additive errors we have

Y_{j} = η(x_{j}; β) + ε_{j} j = 1, 2, ..., n.  (28)

where η(x_{j}; β) is not necessarily linear in β. Again β̂ can be obtained by ML say and residuals formed: r_{j} = Y_{j} − η(x_{j}; β̂) j = 1, 2, ..., n. This time there is no obvious equivalent of standardization so the bootstrap sample is simply

Y_{j}^{∗} = η(x_{j}; β̂) + r_{j}^{∗} j = 1, 2, ..., n.

However the method extends easily to more complex nonlinear situations. For example it may be that we have observations (y_{j}, x_{j}) where

y_{j} ∼ P oisson(μ_{j}), μ_{j} = exp(x_{j}^{T} β).

Again β̂ can be obtained by ML and a bootstrap sample can then be formed by drawing a sample as

y_{j}^{∗} ∼ P oisson(μ_{j}^{∗}), μ_{j}^{∗} = exp(x_{j}^{T} β̂).

A clear introduction to resampling using generalized linear models is given by Hjorth (1994).

4.5 Uses of Metamodels

We end by discussing some examples of why we might wish to carry out bootstrap resampling of a regression metamodel. We give three examples.

The first is when the regression function η(x; β) represents a performance index and we might wish to construct a confidence interval for it at some given value of x, x_{0} say. Suppose that we draw B sets of bootstrap observations of the form (27a)

{Y_{j}^{(i)∗} = η(x_{j}; β̂) + e_{j}^{(i)∗}, j = 1, 2, ..., n}, i = 1, 2, ..., B,

and from each set we calculate a bootstrap estimator β̂^{(i)∗} of β̂ using the ML method say. Then the corresponding set of bootstrap regression function values η(x_{0}, β̂^{(i)∗}) can be used to calculate a bootstrap confidence interval. For example we could use (15) with η̂^{∗} in that formula having

components η̂_{i}^{∗} = η(x_{0}, β̂^{(i)∗}).

The second example is when we wish to find an optimal setting for x. We consider the simplest case where

Y = η(x; β) + ε,

with x scalar and where

η(x; β) = β_{0} + β_{1}x + β_{2}x^{2}.

Then the optimal setting is given by dη/dx = 0, i.e. at x_{opt} = −β_{1}/(2β_{2}). Using ML to estimate this yields

x̂_{opt} = −β̂_{1}/(2β̂_{2}).

Again a simple confidence interval for the unknown x_{opt} is given by (15), only this time, with η̂^{∗}

in that formula having components η̂_{i}^{∗} = −β̂_{1}^{(i)∗}/(2β̂_{2}^{(i)∗}). This is an example where bootstrapping furnishes a relatively easy answer to what in classical statistics is a problem that is not straightforward.

The final example concerns the identification of important factors. Suppose we are interested in identifying those coefficients in the regression function η(x; β) for which |β_{i}| > β_{0i} where β_{0i} > 0 are given constants. Let S = {β_{i} : |β_{i}| > β_{0i}} denote the important set of coefficients. The obvious estimate is to select those coefficients β_{i} for which |β̂_{i}| > β_{0i},i.e. Ŝ = {β_{i} : |β̂_{i}| > β_{0i}}. Bootstrapping is a simple way of assessing the stability of this choice. We generate B bootstrap samples (y_{j}^{(k)∗}, x_{j}^{(k)∗}) j = 1, 2, ..., n, k = 1, 2, ..., B, using either residual sampling or case sampling say and fit the regression model to each bootstrap sample to give bootstrap estimates β̂^{(k)∗}, k = 1, 2, ..., B. We then calculate Ŝ^{(k)∗} = {β_{i} : |β̂_{i}^{(k)∗}| > β_{0i}}, k = 1, 2, ..., K. Now assume that B is sufficiently large so that each distinct bootstrap important set that has been obtained occurs a reasonable number of times. Then a (1−α) confidence region for the unknown true important set can be constructed by selecting bootstrap important sets in decreasing order of their observed frequency of occurrence until (1 − α)100% of the Ŝ^{(k)∗} have been chosen.

Related to identifying important factors in metamodels is the rather more difficult problem of metamodel selection. We consider this in the next section.

4.6 Metamodel Comparison and Selection

Metamodel comparison and selection is a difficult subject. The main trouble is that models that we wish to compare may be of different functional complexity. The statistical properties associated with quantities derived from different models may therefore be hard to put on a common scale on which they can be compared.

A reasonably satisfactory approach is based on the use of cross-validation technique. We suppose that the initial data has the form

S = {(y_{j}, x_{j}), j = 1, 2, ..., n}

with

y_{j} = η(x_{j}, β) + ε_{j}, j = 1, 2, ..., n.

Suppose that our fitted regression is

ŷ(x) = η(x, β̂).

In its basic form cross-validation is used to assess how effective the fitted regression is for predicting the response at some new design point x_{new}. Rather than explicitly choosing such a new point, a simple idea is the leave-one-out method where we drop an observed point (y_{j}, x_{j}) from the set of all observed points, fit the metamodel to the remaining points S_{−j} = S ∖(y_{j}, x_{j}) giving the fitted regression as

ŷ_{−j}(x) = η(x, β̂_{−j}),

and then look at the squared difference between the omitted y_{j} and the value of y at x_{j}, as predicted by the fitted model; i.e.

(y_{j} − ŷ_{−j}(x_{j}))^{2}.

If we do this even handedly by leaving out each observation in turn we have as an estimate of cross-validation prediction error:

L̂_{CV} = n^{−1} ∑_{j=1}^{n}(y_{j} − ŷ_{−j}(x_{j}))^{2}.  (29)

It turns out that, for the linear regression model, this formula simplifies elegantly to one where we only need comparisons with the one model fitted to all the original observations:

L̂_{CV} = n^{−1} ∑_{j=1}^{n} \frac{(y_{j} − ŷ(x_{j}))^{2}}{1 − h_{jj}}

where

ŷ(x_{j}) = x_{j}^{T} β̂,

and h_{jj} is the jth main diagonal entry in the hat-matrix (24), as before.

The distributional property of L̂_{CV} can be obtained in the usual way by bootstrapping to get B bootstrap samples. If we use case resampling say then

S^{(i)∗} = {(y_{j}^{(i)∗}, x_{j}^{(i)∗}), j = 1, 2, ..., n}

where each (y_{j}^{(i)∗}, x_{j}^{(i)∗}) is an observation drawn at random from the set S.

We turn now to model selection. For simplicity we shall only consider the linear case where y = x^{T} β + ε, though with an obvious adjustment, the discussion does generalize. Let β be of dimension p. If we are uncertain as to which design variables are really needed in the model, then we would have in principle up to 2^{p} models to choose from with \binom{p}{q} (sub)models where there were precisely q covariates selected. We denote a typical submodel by M. The estimate of prediction error using (29) is

L̂_{CV} (M) = n^{−1} ∑_{j=1}^{n} \frac{(y_{j} − ŷ_{M}(x_{j}))^{2}}{1 − h_{Mjj}} .

It turns out that this measure does not work as well as it might even when n is large. (This is an example previously referred to in subsection (3.4) where inconsistency arises not in the problem itself but in the bootstrap sampling method.) A much better variant is not to leave out just one observation at a time but instead to split the observations into two: a training set, and an assessment set with respectively n_{t} = n − m and n_{a} = m observations in each. We shall moreover select not one but K such pairs and denote the sets as S_{t,k} and S_{a,k}, k = 1, 2, ..., K.

We shall write β̂_{Mk} for the coefficients of model M fitted to the kth training set, and write ŷ_{Mjk} = x_{Mj}^{T}β̂_{Mk} for the value of x_{j}^{T} β predicted by this model M at x_{j}. Then L̂_{CV} (M) becomes

L̂_{CV} (M) = K^{−1} ∑_{k=1}^{K} m^{−1} ∑_{j∈S_{a,k}}(y_{j} − ŷ_{Mjk})^{2}.  (30)

We use the same set of K pairs for all models being compared. Provided n − m → ∞ and m/n → 1 as n → ∞ then selecting M to minimize (30) will yield a consistent estimator for the correct model. When n is small it may not be possible to select m large enough in which case Davison and Hinkley (1997) suggest taking m/n ≃ 2/3.

The variability of L̂_{CV} (M) can be estimated by bootstrapping in the usual way.

5 Bootstrap Comparisons

In this Section we consider problems where there are a number of different samples to compare. We suppose that we have m data sets

y_{i} = (y_{i1}, y_{i2}, ..., y_{i,n_{i}}), i = 1, 2, ..., m  (31)

the ith being of size n_{i} runs. This situation covers two different scenarios in particular.

The first is where we are considering the question: Does the output of the simulation model accurately represent the output of the system being modelled? Here we have just two data sets, so that m = 2, and one data set comprises observations of a real system whilst the other comprises the observed output from a simulation model. Thus this is a question of validation.

The second is when we have a number of different simulation models to assess. Here validation against real data is not the issue. The different models simply represent different systems, and their validity is not in doubt. Instead we are just interested in assessing how differently the systems behave, as represented by the output of their simulation models. Typically we might be interested in which system produces the greatest average output, or least average delay. Alternatively we might wish to compare output variability.

A similar methodology can be applied to either problem. We discuss this in the remainder of this section.

5.1 Goodness-of-Fit and Validation

A discussion of validation involving trace-driven simulation is described in Kleijnen et al. ( 2000, 2001). The basic idea used there can be generalized. We treat the problem as essentially one of goodness-of-fit.

A simple starting point is where we have n simulation outputs y_{j}, j = 1, 2, ..., n, to compare with m observations y_{j}^{Real} , j = 1, 2, ..., m of a real system and we wish to know if the two samples are identically distributed. A good goodness-of-fit statistic is the two sample Cramer - von Mises statistic (see Anderson, 1962):

W^{2} = ∫ [F_{n}(y) − F_{m}^{Real} (y)]^{2}dH_{n+m}(y)  (32)

where F_{n}, F_{m}^{Real} and H_{n+m} are respectively the EDF’s of the simulated, real and combined samples. The asymptotic distribution of W^{2} is known when the samples are drawn from continuous distributions, but it is easy to set up a bootstrap version of the test that can handle discrete distributions just as well.

We can use a bootstrapping method as follows. We obtain B pairs of bootstrap samples (y^{(i)∗}, y^{Real (i)∗}), i = 1, 2, ..., B, each sample in the pair being obtained by resampling with replacement from y and y^{Real} combined but with y^{(i)∗} the same sample size as y, and y^{Real (i)∗} the same sample size as y^{Real} . Denote the EDFs of the samples of each pair by F_{n}^{(i)∗} and F_{n}^{Real (i)∗}, i = 1, 2, ..., B. We can now calculate B bootstrap two sample Cramer - von Mises statistics from each pair:

W^{(i)∗2} = ∫ [F_{n}^{(i)∗}(y) − F_{m}^{Real (i)∗}(y)]^{2}dH_{n+m}^{(i)∗}(y).  (33)

Under the null assumption that y and y^{Real} are drawn from the same distribution it follows that each W^{(i)∗2} is just a bootstrap version of W^{2}. We can thus calculate a critical p-value from the EDF of the W^{(i)∗2}, i = 1, 2, ..., B. The null hypothesis that y and y^{Real} are drawn from the same distribution can be tested simply by checking if the original W^{2} exceeds this p-value or not.

The validation method just described is easily extended to allow the representational accuracy of a number of different competing simulation models to be compared against (the same) real data. We order the Cramer von-Mises goodness of fit statistics (32) obtained from comparing the real data with the simulated output of each of the models. We can assess the reliability of the comparison by bootstrapping as in (33) to obtain bootstrap distributions of each such Cramer von-Mises statistic and looking at the degree of overlap of amongst these distributions.

When the W^{2} statistic shows the data samples y_{j} to be significantly different, there is an interesting decomposition of W^{2} that allows the nature of the differences to be more clearly identified . Durbin and Knott (1972) show that in the one sample case where F_{n}(.) is being compared with a uniform null distribution, so that W^{2} = n ∫ [F_{n}(y) − y]^{2}dy, then W^{2} has the orthogonal representation

W^{2} = ∑_{j=1}^{∞}(jπ)^{−2}z_{nj}^{2}

where the z_{nj} are scaled versions of the coefficients in the Fourier decomposition of the process\sqrt{n} (F_{n}(y) − y). The z_{nj} are stochastic of course, depending on F_{n}(.), but they are independent and, under the null, have identical distributions. It can be shown that z_{n1}, z_{n2} and z_{n3} are dependent essentially exclusively on deviations respectively of the mean, variance and skewness of F () from that of the null. In other words these components z_{n1}, z_{n2} and z_{n3} can be used to provide convenient statistical tests for these differences. Cheng and Jones (2000, 2004) describe a generalisation of the results of Durbin and Knott (1972), which they term EDFIT statistics, suitable for comparison of data of the form (31) arising in the simulation context. Their formulation uses ranks rather than the original. Though this leads to some loss of power the method does then have the advantage of enabling critical values of tests to be easily obtained by bootstrapping, and moreover in an exact way.

5.2 Comparison of Different Systems

We now consider the situation of (31) where the data sets are outputs from different simulation models only, and no comparison is made with real data. Here we are interested simply in making comparisons between different models. The discussion of the previous subsection goes through with no essential change. The only difference is that all samples have the same status; there is no sample being singled out as special by being drawn from a real system. The EDFIT approach can therefore be used directly, with bootstrapping providing critical null values.

A more direct approach, not using goodness of fit ideas, is possible. We focus on examining differences between the means, ȳ_{i}, i = 1, 2, .., m, of the m samples of (31). Comparison of other sample statistics be handled in the same way. If we suppose that largest is best we can rank the means as ȳ_{(1)} > ȳ_{(2)} > ... > ȳ_{(m)}. The question then is how stable is this order? Bootstrapping provides an easy answer. We generate B bootstrap sets of observations

y_{i}^{(k)∗} = (y_{i1}^{(k)∗}, y_{i2}^{(k)∗}, ..., y_{i,n_{i}}^{(k)∗}), i = 1, 2, ..., m, k = 1, 2, ..., B  (34)

where each sample y_{i}^{(k)∗}, k = 1, 2, ..., B, is obtained by sampling with replacement from y_{i}. We then order the means in each bootstrapped set: ȳ_{(1)}^{(k)∗} > ȳ_{(2)}^{(k)∗} > ... > ȳ_{(m)}^{(k)∗}, k = 1, 2, ..., B. The frequency count of how many times the mean of a given sample comes out on top is a measure its relative merit. A more comprehensive picture is provided by setting up a two-way table showing the number of times ȳ_{i}^{(k)∗} > ȳ_{j}^{(k)∗} out of the B bootstrapped sets, for each possible pair 1 ≤ i, j ≤ m.

The parametric form of the bootstrapping procedure is equally easy to implement, though it is not clear there is much advantage to be had in doing so. Suppose for example that the ith sample y_{i} is drawn from the distribution F^{(i)}(., θ_{i}). Usually the F^{(i)}(., θ_{i}) will have the same functional form, for example all normal distributions, but the procedure works equally well when the F^{(i)} are functionally different. Let the mean of the F^{(i)} distribution be μ^{(i)}(θ_{i}). Then we can estimate, by ML say, the parameters of each distribution from their corresponding sample. Let these estimates be θ̂_{i}. We can then carry out parametric bootstrapping by forming (34) but this time with each sample y_{i}^{(k)∗}, k = 1, 2, ..., B, obtained by sampling from the fitted distribution F^{(i)}(., θ̂_{i}). The analysis then proceeds as before.

6 Bayesian Models

Much of the current popularity of Bayesian methods comes about because increased computing power makes possible Bayesian analysis by numerical procedures, most notably Markov Chain Monte Carlo (MCMC). See Gilks et al. (1996).This allows numerical sampling to be carried out in much the same way as is done in bootstrap and other resampling procedures. We indicate the commonality by considering just one or two situations where Bayesian methods and resampling overlap.

We again focus on a random sample w = (w_{1}, w_{2}, ..., w_{n}) drawn from a distribution F (w, θ). As in the situation where we considered ML estimation, the form of F is known but it depends on a vector θ of parameters. In the ML estimation situation θ_{0},the true value of θ, was assumed unknown. In the Bayesian case a prior distribution is assumed known. We consider just the continuous case and denote the pdf of the prior by π(θ). The main step in Bayesian analysis is to construct the posterior distribution π(θ | w) which shows how the sample w, which depends θ, has modified the prior. The Bayesian formula is

π(θ | w) = \frac{p(w | θ)π(θ)}{∫ p(w | θ)π(θ)dθ} .  (35a)

The difficult part is the evaluation of the normalizing integral ∫ p(w | θ)π(θ)dθ. MCMC has proved remarkably successful in providing a powerful numerical means of doing this.

In the context of simulation a Bayesian approach is useful in assessing uncertainty concerning input parameter values used in a simulation model. We can adopt the viewpoint given in (19) and regard the runs as taking the form:

y_{j} = y(u_{j}, v_{j}, θ, x_{j}) j = 1, 2, ..., n.  (36)

only now θ is assumed to have a Bayesian prior distribution π(θ). Note that in this situation we do not have the equivalent of data w from which to calculate a posterior distribution for θ. What we can do is to treat the prior π(θ) as inducing a prior distribution on y = y(u, v, θ, x). In this sense we can think of y itself as having a prior which can then be estimated by sampling θ from its prior π(θ), yielding values θ_{j} j = 1, 2, ..., n, and then running the model to produce a set of observations

y_{j} = y(u_{j}, v_{j}, θ_{j}, x_{j}) j = 1, 2, ..., n.  (37)

The EDF of these y_{j} then estimates this prior distribution of y. This process clearly is the analogue to the classical situation where θ is estimated from data w, as is supposed in (19), and we then allowed for this variability in the y_{j} by replacing θ̂(w) by θ_{j}^{∗} calculated from (20) or (21).

More interestingly we can take this a step further if there are real observations available of the y process itself. We do not attempt a fully general formulation but give an example to indicate the possibilities. Suppose that

y_{j}^{Real} = y^{Real} (x_{j}), j = 1, 2, ..., n

comprise n observations on the real system being modelled by the simulation. The precise relationship between the simulation output y = y(u, v, θ, x) and y^{Real} (x) is not known. However we might assume that

y(u, v, θ, x) ∼ N(y^{Real} (x), σ^{2})

where σ^{2} is an additional parameter that will also be treated as Bayesian in the sense of having a prior. We can express our great prior uncertainty about σ by assuming a reference prior distribution, ρ(σ), for it. Then, by (35a), the posterior distribution of (θ, σ) is proportional to

σ^{−n}π(θ)ρ(σ) exp\{−\frac{1}{2σ^{2}} ∑_{j=1}^{n}[y(u, v, θ, x_{j}) − y^{Real} (x_{j})]^{2}\}.

The posterior distribution can thus be obtained by MCMC methods for example. An interesting application of this problem occurs in epidemiological modelling. Suppose that y is a measure of the progress of an epidemic that is dependent on factors θ for which there are previous measurements or for which there is expert information. A Bayesian approach is then a natural way of incorporating this prior information. We need also to have a good epidemic model for producing a simulated y. Thus reverse use of simulation as indicated above has allowed this prior information to be updated.

7 Time Series Output

Bootstrapping of time series is a well studied problem. In simulation the most likely use of such procedures is to generate correlated input for a model. As usual the parametric form is relatively easy to explain and implement and we discuss this first.

7.1 Residual Sampling

We consider the case where the time-series is an autoregressive model. Here the residual sampling method used to construct a bootstrap metamodel applies with little change. Suppose we have the autorgressive model

Y_{t} = ∑_{j=1}^{p} a_{j}Y_{t−j} + ε_{t}

where the ε_{t} are independently and identically distributed, commonly called the innovations. Suppose that y_{1}, y_{2}, ..., y_{n} are a series drawn from this model. Then we can estimate the a_{j} by least squares say, yielding the estimates â_{j}, j = 1, 2, ..., p, and form the residuals

r_{t} = y_{t} − ∑_{j=1}^{p} â_{j}y_{t−j}, t = p + 1, p + 2, ..., n.

We can then form the bootstrap time-series as

y_{t}^{∗} = ∑_{j=1}^{p} â_{j}y_{t−j}^{∗} + r_{t}^{∗}, t = 1, 2, ..., n

by sampling the r_{t}^{∗} from the EDF of the residuals {r_{t}}. We need y_{1}^{∗}, y_{2}^{∗}, ..., y_{p}^{∗} to initiate the process, but if we assume that the observed series is stationary, it is probably easiest to simply initiate the series with some arbitrary starting values, possibly the original y_{1}, y_{2}, ..., y_{p}, then run the bootstrapping sufficiently long until initial effect are negligible and collect the actual y_{t}^{∗} from that point on. Freedman (1984) gives conditions for the asymptotic validity of this procedure, Basawa et al. (1989) extending these results to the case of nonlinear time-series.

7.2 Block Sampling

For time-series the analogue of case sampling is block sampling. We cannot sample individual observations y_{t} because this loses the correlation between observations. If the series is long enough then we can take n = bl and think of the series as comprising b blocks each of length l. We write the ith block as y_{i} = (y_{l(i−1)+1}, y_{l(i−1)+2}, ..., y_{(li)}) i = 1, 2, ..., b. Bootstrapping is done by sampling blocks with replacement from this set of b blocks, retaining the order of the observations in each block when writing down the individual observations of the bootstrapped series. A balance needs to be struck between having block lengths long enough to retain the correlation properties between neighboring observations, and having enough blocks to measure the inherent variation of the series. A typical compromise is to use say b = l = n^{1/2} so that both quantities tend to infinity as n → ∞.

A major weakness of block sampling is the loss of correlation incurred by the random sampling of blocks. This loss of correlation is called whitening. It is especially serious when the statistic of interest involves correlations of high lag. The crude block sampling just described may be quite ineffective if the size of blocks is not large enough, because calculation involves quantities which straddle blocks and which are then not sufficiently correlated because of the whitening. There are many variants of block sampling aimed at reducing the effect of whitening in specific situations. A good example is estimation of the lag m covariance

c_{m} = \frac{1}{n − m} ∑_{t=1}^{n−m}(y_{t} − ȳ)(y_{t+m} − ȳ).

Here we can define a two-dimensional process

z_{t} = ( \begin{matrix}z_{1t} \\ z_{2t} \end{matrix} ) = ( \begin{matrix}y_{t} \\ y_{t+m} \end{matrix} ) , t = 1, 2, .., n − m

with z̄_{i} = (n − m)^{−1} ∑_{t=1}^{n−m} z_{it} and think of c_{1} as a statistic of this process

c_{1} = \frac{1}{n − 1} ∑_{t=1}^{n−1}(z_{1t} − z̄_{1})(z_{2t} − z̄_{2}).

We can then obtain bootstrap z_{t}^{∗} by sampling with replacement from the set {z_{t}}.The bootstrap lag m covariance then clearly substantially retains the covariance of the original series as we have, in effect, bootstrap sampled the terms (y_{t} − ȳ)(y_{t+m} − ȳ) appearing in the formula giving c_{m}. Generalizations of this technique are known as block of blocks sampling.

7.3 Spectral resampling

Residual and block sampling are time domain techniques. An alternative approach is to sample in the frequency domain. A big advantage is that spectral increments are uncorrelated, and for Gaussian processes this strengthens to independent increments.

Suppose we have n = 2m + 1 observations

y_{t}, t = −m, −m + 1, ..., m − 1, m

for which there is a continuous spectral density S(ω) and define the frequencies ω_{k} = 2πk/n, −m ≤ k ≤ m. Then the observations have the spectral representation

y_{t} = ∑_{k=−m}^{m} a_{k}e^{iω_{k}t}

where

a_{k} = n^{−1} ∑_{t=−m}^{m} y_{t}e^{−iω_{k}t}.

In this section i = \sqrt{−1}. For a Gaussian process the real and purely imaginary components of the a_{k} are independent, and the a_{k} are asymptotically so.

Norgaard (1992) gives two possible ways of obtaining bootstrap samples of the a_{k.} The simplest version is to draw a_{k}^{∗} at random from one of the twelve elements

(±a_{k−1}, ±a_{k}, ±a_{k+1}, ±ia_{k−1}, ±ia_{k}, ±ia_{k+1})  (38)

when 0 < k < n and to draw a_{n}^{∗} at random from one of the twelve elements

(±a_{n−1}, ±a_{n}, ±a_{n}, ±ia_{n−1}, ±ia_{n}, ±ia_{n}).  (39)

(note the repetition of some elements in this latter case). In both cases we set

a_{−k}^{∗} = \overline{a_{k}^{∗}}

i.e. a_{−k}^{∗} is the complex conjugate of a_{k}^{∗}. The value of a_{0}^{∗} needs to be real. The simplest case is when E(Y_{t}) = 0 when we can select a_{0}^{∗} at random as one of the four elements (±a_{0}, ±|a_{1}|).

A variant which seems to be significantly superior is called the extended circle (EC) method by Norgaard (1992). Here we select |a|_{k}^{∗} as a random element of the set (|a_{k−1}|, |a_{k}|, |a_{k+1}|) and then give this a random spin in both the real and imaginary directions:

a_{k}^{∗} = |a|_{k}^{∗}(cos ϕ_{1k} + i sin ϕ_{2k}), k = 1, 2, ..., n

where ϕ_{1k} and ϕ_{2k} are independent U(−π, π) angles. We have a_{−k}^{∗} = \overline{a_{k}^{∗}} as before.

8 Final Comment

There are many topics not covered in this appendix that could well have been included. We have for instance, not touched on the closely related technique of jackknifing. Other topics we have made only cursory mention, such as cross validation and when bootstrapping fails.

References

[1] Anderson, T. W. 1962. On the distribution of the two-sample Cramer-von Mises criterion. Ann. of Math. Statist., 33, 1148-1159.

[2] Bickel, P. J. and D. A. Freedman. (1981). Some asymptotic theory for the bootstrap. Ann. Statist., 9, 1196-1217.

[3] Cheng, R. C. H. and O. D. Jones. 2000. Analysis of simulation factorial experiments by EDF resample statistics. Proceedings of the 2000 Winter Simulation Conference, eds J. A. Joines, R. Barton, P. Fishwick, K. Kang. Piscataway, New Jersey: IEEE, 697-703.

[4] Cheng, R. C. H. and Jones, O.D. 2004. Analysis of distributions in factorial experiments. Statistica Sinica. To appear.

[5] Cheng, R. C. H. 1995. Bootstrap methods in computer simulation experiments. In Proceedings of 1995 Winter Simulation Conference, eds W. R. Lilegdon, D. Goldsman, C. Alexopoulos and K. Kang. Piscataway, New Jersey: IEEE, 171-177.

[6] Chernick, M. R. 1999. Bootstrap Methods. New York: John Wiley.

[7] Cox, D. R. and D. V. Hinkley. 1974. Theoretical Statistics. London: Chapman & Hall.

[8] Davison, A. C., D. V. Hinkley and E. Schechtman. 1986. Efficient bootstrap simulations. Biometrika, 73, 555-566.

[9] Davison, A. C. and D. V. Hinkley. 1997. Bootstrap Methods and Their Application. Cambridge, England: Cambridge University Press.

[10] Durbin, J. and M. Knott. 1972. Components of Cramer-von Mises statistics I. J. Royal Statist. Soc., 34, 290-307.

[11] Efron, B. 1979. Bootstrap methods: another look at the jackknife. Ann. Statist., 7, 1-26.

[12] Efron, B. 1987. Better bootstrap confidence intervals. J. Amer. Statist Assoc. 82, 171-185.

[13] Efron, B. and R. J. Tibshirani. 1993. An Introduction to the Bootstrap. New York and London: Chapman and Hall.

[14] Gilks, W. R., S. Richardson and D. J. Spiegelhalter. 1996. Markov Chain Monte Carlo in Practice. London: Chapman & Hall.

[15] Hall, P. 1988. Rate of convergence in bootstrap approximations. Ann. Prob., 16, 1165-1685.

[16] Hall, P.and Pittelkow, Y.E. 1990, "Simultaneous bootstrap confidence bands in regression." J. Statist. Comput. Simul. 37, 99-113.

[17] Hall, P. 1992. The Bootstrap and Edgeworth Expansion. New York: Springer-Verlag.

[18] Hjorth, J. S. U. 1994. Computer Intensive Statistical Methods. London: Chapman & Hall.

[19] Johns, M. V. Jr. 1988. Importance sampling for bootstrap confidence itervals. J. Amer. Statist. Assoc., 83, 709-714.

[20] Kleijnen, J.P.C., R. C. H. Cheng and B. Bettonvil. 2000. Validation of trace-driven simulation models: more on bootstrap tests. In Proceedings of the 2000 Winter Simulation Conference, eds J. A. Joines, R. R. Barton, K. Kang and P. A. Fishwick, IEEE Piscataway, New Jersey, 882- 892

[21] Kleijnen, J. P. C., R. C. H. Cheng and B. Bettonvil. 2001 Validation of Trace Driven Simulation Models: Bootstrap Tests, Management Science, 47, 1533-1538.

[22] Norgaard, A. 1992. Resampling stochastic processes using a bootstrap approach. In Bootstrapping and Related Topics. Lecture Notes in Economics and Mathematical Systems, Springer, 376, 181-185.

[23] Serfling, R. J. 1980. Approximation Theorems of Mathematical Statistics. New York: Wiley.

[24] Shao, J. and D. Tu. 1995. The Jackknife and the Bootstrap. New York: Springer.

[25] Singh, K. 1981. On the asymptotic accuracy of Efron’s bootstrap. Ann. Statist., 9, 1187-1195.

[26] van Beeck, P. 1972. An application of Fourier methods to the problem of sharpening the Berry-Esséen inequality. Z. Wahrscheinlichkeitstheorie und. Verw. Gebiete, 15, 279-290.

[27] von Mises, R. 1947. On the asymptotic distribution of differentiable statistical functions. Ann. Math. Statist., 18, 309-348.

[28] Wilks, S. S. 1962. Mathematical Statistics. New York: Wiley.

[29] Wu, C. F. J. 1986. Jackknife, bootstrap and other resampling methods in regression analysis. Ann. Statist., 14, 1261-1295.