# Calculating confidence intervals via bootstrap on dependent observations

The bootstrap, in its standard form, can be used to calculate confidence intervals of estimated statistics provided that observations are iid. I. Visser et al. in "Confidence Intervals for Hidden Markov Model Parameters," used a parametric bootstrap to calculate CIs for HMM parameters. However, when we fit an HMM on an observation sequence, we have already assumed the observations are dependent (in contrast to mixture models).

I have two questions:

1. What does the iid assumption do with the bootstrap?
2. Can we ignore the iid requirement in a parametric bootstrap?

Visser et al. method is briefly as follows:

1. Assume we have an observation sequence $Y=o_1,o_2,...,o_n$ resulted from sampling an HMM with real but unknown set of parameters $\theta=\theta_1,\theta_2,...,\theta_l$.
2. The parameters can be estimated using the EM algorithm: $\hat{\theta}=\hat{\theta}_1,\hat{\theta}_2,...,\hat{\theta}_l$
3. Use the estimated HMM to generate a bootstrap sample of size $n$: $Y^*=o^*_1,o^*_2,...,o^*_n$
4. Estimate parameters of the HMM according to the bootstrap sample: $\hat{\theta}^*=\hat{\theta}^*_1,\hat{\theta}^*_2,...,\hat{\theta}^*_l$
5. Repeat the steps 3 and 4 for $B$ times (e.g. $B$=1000) resulting in $B$ bootstrap estimations: $\hat{\theta}^*(1),\hat{\theta}^*(2),...,\hat{\theta}^*(B)$
6. Calculate the CI of each estimated parameter $\hat{\theta}_i$ using the distribution of $\hat{\theta}^*_i$ in bootstrap estimations.

Notes (my findings):

1. The percentiles method should be used to calculate CIs in order to have correct coverage (normality is a bad assumption).
2. Bias of the bootstrap distribution should be corrected. Meaning that the distribution mean of $\hat{\theta}^*_i$ should be shifted to $\hat{\theta}_i$
-
First question in other words: What is the effect of iid assumption on the bootstrap? Is it a simplifying assumption that can be removed by following a more complex algorithm, or formula? –  Sadeghd Aug 13 '11 at 17:37

Short answers: 1. It simplifies it. (Frankly, I did not get the question). 2. No, you can never ignore it, as lack of i.i.d. has immediate consequences on the variances of whatever you are estimating.

Medium answer: Pretty much the central issue with the bootstrap is, 'Does the proposed procedure reproduce the features of the data?'. Violation of the i.i.d. assumption is a big deal: your data are dependent, you (most likely) have less information in your data than you would have in an i.i.d. sample of the same size, and if you run a naive bootstrap (resample the individual observations), the standard errors you get from it will be too small. The proposed procedure circumvents the problem of lack of independence by capturing (or at least attempting to capture) the dependence in the model structure and parameters. If successful, each bootstrap sample would reproduce the features of the data, as needed.

Long answer: There are multiple layers of assumptions concerning the bootstrap, and even in the simplest possible case (i.i.d. data, estimation of the mean), you have to make at least three assumptions: (1) the statistic of interest is a smooth function of the data (true in the case of the mean, not so true even in the case of percentiles, totally off with say nearest neighbor matching estimators); (2) the distribution from which you bootstrap is "close" to the population distribution (works OK in the case of i.i.d. data; may not work OK in the case of dependent data, where you essentially have only one trajectory = one observation in the case of time series, and you have to invoke additional assumptions like stationarity and mixing to strecth this single observation into a quasi-population); (3) your Monte Carlo bootstrap sampling is a good enough approximation to the complete bootstrap with all possible subsamples (the inaccuracy from using Monte Carlo vs. the complete bootstrap is much less than the uncertainty you are trying to capture). In the case of the parametric bootstrap, you also make an assumption that (4) your model perfectly explains all the features of the data.

As a warning of what could go wrong with (4), think about regression with heteroskedastic errors: $y=x\beta + \epsilon$, Var$[\epsilon] = \exp[ x\gamma]$, say. If you fit an OLS model and resample the residuals as if they were i.i.d., you will get a wrong answer (some sort of $\bar\sigma^2 (X'X)^{-1}$ where $\bar\sigma^2$ is the average $1/n \sum_i \exp[x_i \gamma]$, instead of the appropriate $(X'X)^{-1} \sum \exp[x_i \gamma] x_i x_i' (X'X)^{-1}$). So if you wanted to have a fully parametric bootstrap solution, you would've have to fit the model for heteroskedasticity along with the model for the mean. And if you suspect serial or other sort of correlation, you would've have to fit the model for that, too. (See, the non-parametric distribution-free flavor of the bootstrap is pretty much gone for now, as you have replaced the voice of the data with the synthesized voice of your model.)

The method you described works around the i.i.d. assumption by creating a whole new sample. The greatest problem with the dependent data bootstrap is to create the sample that would have the dependence patterns that would be sufficiently close to those in the original data. With time series, you could use block bootstraps; with clustered data, you bootstrap the whole clusters; with heteroskedastic regression, you have to with wild bootstraps (which is a better idea than the bootstrap of residuals, even if you have fitted a heteroskedasticty model to it). In the block bootstrap, you have to make an educated guess (or, in other words, have good reasons to believe) that distant parts of time series are approximately independent, so that all of the correlation structure is captured by the adjacent 5 or 10 observtations that form the block. So instead of resampling observations one by one, which totally ignores the correlation structure of the time-series, you resampling them in blocks, hoping that this would respect the correlation structure. The parametric bootstrap you referred to says: "Rather than fiddling with the data and assembling the new dolls from the pieces of the old ones, why don't I just stamp the whole molded Barbie for you instead? I've figured out what kind of Barbies you like, and I promise I will make you one you'd like, too."

In case of the parametric bootstrap you described, you have to be pretty damn sure that your HMM model fit is pretty much perfect, otherwise your parametric bootstrap may lead to incorrect results (Barbies that cannot move their arms). Think about the above heteroskedastic regession example; or think about fitting an AR(1) model to AR(5) data: whatever you do with the parametrically simulated data, they won't have the structure the original data used to have.

Edit: as Sadeghd clarified his question, I can respond to that, as well. There is a humongous variety of the bootstrap procedures, each addressing the particular quirk in either the statistic, the sample size, the dependence, or whatever an issue with the bootstrap could be. There is no single way to address dependence, for instance. (I've worked with survey bootstraps, there are about 8 different procedures, although some are mostly of methodological rather than practical interest; and some are clearly inferior in that they are only applicable in special, not easily generalizable, cases.) For a general discussion of issues you could face with the bootstrap, see Canty, Davison, Hinkley and Ventura (2006). Bootstrap diagnostics and remedies. The Canadian Journal of Statistics, 34 (1), 5-27.

-
Just to add a bit to your statement about having less information when you have dependent clusters of data (in the Medium section), I believe this is true where there is positive intraclass correlation within a cluster, but the opposite is true when there is negative intraclass correlation. Of course, it seems that in most real data applications intraclass correlations are positive. –  Macro Aug 13 '11 at 17:03
@Macro: you are certainly right on both counts (that this is technically possible, and that it is practically irrelevant). The same will be true if you estimate the mean level of an AR(1) process with a negative correlation, but again I am at a loss thinking of real processes that might have this feature. Unlike the positive autocorrelation that is self-reproducible at different time scales, the negative correlation has to disappear if you double the length of your reference period. (The business cycles data, like US GDP, have negative correlations at the lag length of about three years.) –  StasK Aug 13 '11 at 17:25
Thanks for your detailed answer. I concluded that the parametric re-sampling may diminish the effect of dependence. However, the parametric distribution must be, to a good extent, representative of the true population, and the dependency patterns be regenerated in re-sampling. –  Sadeghd Aug 13 '11 at 18:08