[R] GLS models - bootstrapping

[Christian Kamenik]

Dear Lillian,

I tried to estimate parameters for time series regression using time 
series bootstrapping as described on page 434 in Davison & Hinkley 
(1997) - bootstrap methods and their application. This approach is based 
on an AR process (ARIMA model) with a regression term (compare also with 
page 414 in Venable & Ripley (2002) - modern applied statistics with S) 
I rewrote the code for R (this comes without any warranty):

fit <- function( data )
{ X <- cbind(rep(1,100),data$activ)
   para <- list( X=X,data=data)
   assign("para",para)
   d <- arima(x=para$data$temp,order=c(1,0,0),xreg=para$X)
   res <- d$residuals
   res <- res[!is.na(res)]
   list(paras=c(d$model$ar,d$reg.coef,sqrt(d$sigma2)),
        res=res-mean(res),fit=X %*% d$reg.coef)
}
beaver.args <- fit( beaver )
white.noise <- function( n.sim, ts) sample(ts,size=n.sim,replace=T)
beaver.gen <- function( ts, n.sim, ran.args )
{ tsb <- ran.args$res
   fit <- ran.args$fit
   coeff <- ran.args$paras
   ts$temp <- fit + coeff[4]*arima.sim( model=list(ar=coeff[1]),
                             n=n.sim,rand.gen=white.noise,ts=tsb )
   ts }
new.beaver <- beaver.gen( beaver, 100, beaver.args )

beaver.fun <- function(ts) fit(ts)$paras
beaver.boot <- tsboot( beaver, beaver.fun, R=99,sim="model",
                n.sim=100,ran.gen=beaver.gen,ran.args=beaver.args)
names(beaver.boot)
beaver.boot$t0
beaver.boot$t[1:10,]

Maybe there is a more elegant way for doing this. Anyway, boot.ci should 
give you confidence intervals.

Let me know how you are doing.

Best, Christian




> From: Lillian Sandeman <l.sandeman>
> Date: Mon, 2 Oct 2006 13:59:09 +0100 (BST)
> 
> Hello,
> 
> I am have fitted GLS models to time series data. Now I wish to bootstrap
> this data to produce confidence intervals for the model.
> 
> However, because this is time series data, normal bootstrapping is not
> applicable. Secondly, 'tsboot' appears to only be useful for ar models -
> and does not seem to be applicable to GLS models.
> 
> I have written code in R to randomly sample blocks of the data (as in
> Davison & Hinkley's book - bootstrap methods and their application) and
> use this resampling to re-run the model, but this does not seem to be the
> correct approach since Confidence Intervals produced do not show the
> underlying pattern (cycles) in the data [even when block length is
> increased, it only picks up a little of this variation].
> 
> Any help as to how to proceed with this would be greatly appreciated, as I
> cannot find anything applicable on the R pages. Alternatively, if there
> is another method to proceed with this (other than bootstrapping), I would
> also be happy to try it.
> 
> Thankyou,
> 
> Lillian.