R/bootstrap.R

In jpburgard/mixtureVAR: mixtureVAR

# #' bootstrap mixtureVAR model
# #'
# #' @param mod the model to be bootstrapped
# #' @param nboot number of bootstrap iterations
# #' @param eps precision default to 1e-8
# #'
# #' @return return
# #' @export
# Boot_FMVAR = function(mod,nboot=499,eps=1e-8)
# {
#   Y.orig <- mod$Y.orig
#   maxLag <- max(mod$P)
#   maxt   <- length(mod$Tau[[1]])
#   PDiff  <- lapply(mod$P, function(pp) maxLag - pp)
#
#   boot_errors_l      <- mapply(function(tau,e)t(tau * t(e)),mod$Tau,mod$e,SIMPLIFY = FALSE)
#   boot_errors        <- Reduce('+', boot_errors_l)
#   boot_errors_list   <- lapply(1:nboot,function(leer) t(apply(boot_errors,1,sample,size=maxt,replace=TRUE)))
#
#   Y_synth    <- mapply(ysynthfun,mod$X,mod$Theta,mod$Tau,PDiff,SIMPLIFY = FALSE)
#   Y_synth    <- Reduce('+',Y_synth)
#   Y_boot     <- mapply(function(err) {
#     Y_synth+err
#   } , boot_errors_list, SIMPLIFY = FALSE)
#
#   FMVAR_boot <- lapply(Y_boot,function(yboot) {
#     Yboot <- cbind(Y.orig[,1:maxLag],yboot)
#     boot  <- try(mixtureVAR::FMVAR(Yboot, mod$P, mod$EXO.orig, nreps=1, Theta_start=mod$Theta,
#                        multinommod = mod$multinommod,
#                        maxiter = mod$maxiter,
#                        eps=eps),silent = TRUE)
#     cat(".")
#     boot
#   })
#   FMVAR_boot <- FMVAR_boot[!sapply(FMVAR_boot,inherits,"try-error")]
#   cat("\n")
#   return(FMVAR_boot)
# }

#' bootstrap mixtureVAR model
#'
#' Allows to perform a bootstrap on an estimated mixtureVAR model. This bootstrap is necessary for computing the impulse response functions.
#'
#' @param mod the model to be bootstrapped
#' @param nboot number of bootstrap iterations
#' @param eps precision default to 1e-8
#'
#' @return return
#'
#' @export
Boot_mixtureVAR <- function(mod,nboot=499,eps=1e-8)
{
  Y.orig <- mod$Y.orig
  maxLag <- max(mod$P)
  maxt   <- length(mod$Tau[[1]])
  PDiff  <- lapply(mod$P, function(pp) maxLag - pp)

  boot_errors_l      <- mapply(function(tau,e)t(tau * t(e)),mod$Tau,mod$e,SIMPLIFY = FALSE)
  boot_errors        <- Reduce('+', boot_errors_l)
  boot_errors_list   <- lapply(1:nboot,function(leer) t(apply(boot_errors,1,sample,size=maxt,replace=TRUE)))

  Y_synth    <- mapply(ysynthfun,mod$X,mod$Theta,mod$Tau,PDiff,SIMPLIFY = FALSE)
  Y_synth    <- Reduce('+',Y_synth)
  Y_boot     <- mapply(function(err) {
    Y_synth+err
  } , boot_errors_list, SIMPLIFY = FALSE)

  mixtureVAR_boot <- lapply(Y_boot,function(yboot) {
    Yboot <- lapply(Y.orig,function(y)cbind(y[,1:maxLag],yboot))
    boot  <- try(mixtureVAR::mixtureVAR(YFormula = Yboot, LatentFormula = mod$LatentFormula,
                        LatentLag = mod$LatentLag,
                        EXOFormula =  mod$EXO.orig, P = mod$P, data=mod$data,
                        nreps=1, Theta_start=mod$Theta,maxiter = mod$maxiter,
                       eps=eps),silent = TRUE)
    cat(".")
    boot
  })
  mixtureVAR_boot <- mixtureVAR_boot[!sapply(mixtureVAR_boot,inherits,"try-error")]
  cat("\n")
  return(mixtureVAR_boot)
}





#' plot bootstrapped irf's of mixtureVAR
#'
#' @param BOOTAUD bootstrap generated by Boot_FMVAR
#' @param mod the model to be bootstrapped
#' @param impulse which variable should be the impulse variable
#' @param n.ahead timepoint ahead for the irf
#' @param bands lower and upper quantile to be displayed
#' @param type type of plot used
#' @param cumulative TRUE/FALSE Should cumulative irf be produced
#'
#' @export
plotBootIrf <- function(BOOTAUD,mod,impulse,response="all",n.ahead=20,
                        bands=c(0.05,0.95),type=c("mean","median","both","p-value"),
                        relation="same", grey=FALSE,modnames=NULL, cumulative=FALSE){
  if(length(bands)>2) stop("bands must have length 2: lower and upper quantile for band")
  if(!impulse%in%rownames(mod$Theta[[1]])) stop("impulse must be a variable of the model")
  mod1irf  <- irf(mod,n.ahead=n.ahead,cumulative=cumulative)
  mod1irfmat<-abind(do.call(mapply,
                            args = list(FUN=abind,SIMPLIFY=FALSE,along=3,mod1irf[[1]])),along=4)
  if(!is.null(modnames) & length(modnames)==dim(mod1irfmat)[4])
    model = modnames
  else
    model = paste0("Model",1:dim(mod1irfmat)[4])
  DIF1 <- mod1irfmat[,,,1]-mod1irfmat[,,,2]
  mod1irfmat <- abind::abind(mod1irfmat,DIF1,along=4)
  dimnames(mod1irfmat) <- list(time=paste0("t:",(1:dim(mod1irfmat)[1]-1)),
                               variable = dimnames(mod1irfmat)[[2]],
                               svariable = dimnames(mod1irfmat)[[3]],
                               model = c(modnames,"difference"))

  BOOTIRF      <- lapply(BOOTAUD,irf,n.ahead=n.ahead,cumulative=cumulative)
  BOOTIRF2     <- t(lapply(BOOTIRF,"[[",1))
  BOOTIRFfinal <- abind(lapply(1:length(BOOTIRF2[[1]]),function(i)
    abind(do.call(mapply,
                  args = list(FUN=abind,SIMPLIFY=FALSE,along=3,lapply(BOOTIRF2,"[[",i))),along=4)
  ),along=5)
  DIF <- BOOTIRFfinal[,,,,1]-BOOTIRFfinal[,,,,2]
  BOOTIRFfinal <- abind(BOOTIRFfinal,DIF,along=5)
  dimnames(BOOTIRFfinal) <- list(time=paste0("t:",(1:dim(BOOTIRFfinal)[1]-1)),
                                 variable = dimnames(BOOTIRFfinal)[[2]],
                                 svariable = dimnames(BOOTIRFfinal)[[3]],
                                 boot=paste0("B:",1:dim(BOOTIRFfinal)[4]),
                                 model=c(modnames,"difference"))
  # bands <- sort(unique(c(bands,.5))) # change for article
  bands                  <- c(0.025,.16,.5,0.84,.975)
  BOOTIRFfinalQ          <- apply(BOOTIRFfinal,c(1:3,5),quantile,probs=bands)
  LATQ                   <- expand.grid(dimnames(BOOTIRFfinalQ)[-1])
  LATQ$min               <- as.vector(apply(BOOTIRFfinal,c(1:3,5),min))
  LATQ$max               <- as.vector(apply(BOOTIRFfinal,c(1:3,5),max))
  LATQ$lower1            <- as.vector(BOOTIRFfinalQ[1,,,,])
  LATQ$lower2            <- as.vector(BOOTIRFfinalQ[2,,,,])
  LATQ$median            <- as.vector(BOOTIRFfinalQ[3,,,,])
  LATQ$upper2            <- as.vector(BOOTIRFfinalQ[4,,,,])
  LATQ$upper1            <- as.vector(BOOTIRFfinalQ[5,,,,])
#  LATQ$origirf           <- as.vector(mod1irfmat)
  LATQ$mean              <- as.vector(apply(BOOTIRFfinal,c(1:3,5),mean))
  LATQ$pvalue            <- as.vector(apply(BOOTIRFfinal,c(1:3,5),function(a){
    meana <- mean(a<0)
    ifelse(mean(a)>=0,meana,1-meana)
    })) # only relevant for model="difference"
  #LATQ[LATQ$time=="t:0" & LATQ$variable == LATQ$svariable,"mean"]
  SummaryOutput <<- LATQ
  cat("Variable >> SummaryOutput << created\n")

  if(response[1]=="all") response <- unique(LATQ$variable)
  if(grey){
    mmcolors <- c("black","darkgrey")
    mmlty    <- c(1,1)
    mmpch    <- c(1,4)
  } else {
    mmcolors <- c("blue","black")
    mmlty    <- c(1,1)
    mmpch    <- c(1,4)
  }
  if(type=="both"){
    PLOT <- xyplot( mean ~ time| model + variable , data=LATQ, subset = svariable==impulse & variable%in%response & model!="difference",
                    upper = LATQ$upper, lower = LATQ$lower, median = LATQ$median, prepanel = prepanel.default.xyplot2,
                    panel = function(x, y, upper, lower, subscripts, median, ...){
                      #panel.superpose(x, y, panel.groups = my.panel.bands, type='l', col='gray',...)
                      upper  <- upper[subscripts]
                      lower  <- lower[subscripts]
                      median <- median[subscripts]
                      panel.polygon(c(x, rev(x)), c(upper, rev(lower)), fill=.5, col="lightgrey", border=FALSE,...)
                      panel.abline(h=0,col="red")
                      panel.xyplot(x, median, type='b', cex=0.4, lty=1, col=mmcolors[1], pch = mmpch[2],...)
                      panel.xyplot(x, y, type='b', cex=0.4, lty=1, col=mmcolors[1], pch = mmpch[1],...)
                    },ylab="",
                    #main=paste("The impulse variable is",impulse,"with bands [",paste(range(bands),collapse = ";"),"]"),
                    key=list(space="top", columns = 2,lines=list(col=mmcolors, type="b", pch=mmpch, lty=mmlty, lwd=2),
                             text=list(c("bootirf mean","bootirf median"))
                    ), strip=strip.custom(  ),
                    strip.left=TRUE)
  } else if(type=="mean"){
    PLOT <- xyplot( mean ~ time| model + variable , data=LATQ, subset = svariable==impulse& variable%in%response & model!="difference",
                    upper1 = LATQ$upper1, lower1 = LATQ$lower1,upper2 = LATQ$upper2, lower2 = LATQ$lower2, prepanel = prepanel.default.xyplot2,
                    panel = function(x, y, upper1, lower1, upper2, lower2, subscripts, ...){
                      upper1 <- upper1[subscripts]
                      lower1 <- lower1[subscripts]
                      upper2 <- upper2[subscripts]
                      lower2 <- lower2[subscripts]
                      panel.polygon(c(x, rev(x)), c(upper1, rev(lower1)), fill=.5, col="lightgrey", border=FALSE, ... )
                      panel.polygon(c(x, rev(x)), c(upper2, rev(lower2)), fill=.5, col="grey", border=FALSE, ... )
                      panel.abline(h=0, col="red")
                      panel.xyplot(x, y, type='b', cex=0.2, lty=1, col=mmcolors[1], pch = mmpch[1], ... )
                    },ylab="",xlab="",
                    #main=paste("The impulse variable is",impulse,"with bands [",paste(range(bands),collapse = ";"),"]"),
                    #key=list(space="top", columns = 1,lines=list(col=mmcolors[1], type="b",  pch=mmpch[1], lty=mmlty[1], lwd=2),text=list(c("bootirf mean"))),
                    strip=strip.custom(  ),
                    strip.left=TRUE)
  } else if(type=="median"){
    PLOT <- xyplot( median ~ time| model + variable , data=LATQ, subset = svariable==impulse& variable%in%response & model!="difference",
                    upper = LATQ$upper, lower = LATQ$lower, prepanel = prepanel.default.xyplot2,
                    panel = function(x, y, upper, lower, subscripts, ...){
                      upper <- upper[subscripts]
                      lower <- lower[subscripts]
                      panel.polygon(c(x, rev(x)), c(upper, rev(lower)), fill=.5, col="lightgrey", border=FALSE,...)
                      panel.abline(h=0,col="red")
                      panel.xyplot(x, y, type='b', cex=0.4, lty=1, col=mmcolors[1], pch = mmpch[1],...)
                    },ylab="",xlab="",
                    #main=paste("The impulse variable is",impulse,"with bands [",paste(range(bands),collapse = ";"),"]"),
                    #key=list(space="top", columns = 1,lines=list(col=mmcolors[1], type="b", pch=mmpch[1], lty=mmlty[1], lwd=2),text=list(c("bootirf median"))), strip=strip.custom(  ),
                    strip.left=TRUE)
  } else if(type=="p-value"){
    PLOT <- xyplot( pvalue ~ time| variable , data=LATQ, subset = svariable==impulse& variable%in%response & model=="difference",
                    #prepanel = prepanel.default.xyplot2,
                    panel = function(x, y, ...){
                      panel.abline(h=c(0.05,0.1),col="lightgrey",lwd=0.9)
                      panel.xyplot(x, y, type='b', cex=0.2, lty=1, col=mmcolors[1], pch = mmpch[1], ... )
                    },ylab="",xlab="",ylim = c(0,1), strip=strip.custom(factor.levels=c("Industr. Prod.","Inflation")))
  } else stop("plot type not defined: Please use either 'mean','median' or 'both' as argument for type.")
  #library(latticeExtra)
  at.seq <- seq(0,n.ahead,by = 3)+1
  at.lab <- at.seq-1
  at.lab[!((1:length(at.lab))%%2)] <- ""
  at.lab[1] <- ""
  if(type!="p-value"){
  PLOT <- update(PLOT, scales = list(y=list(relation=relation,alternating=3),
                                     x=list(alternating=1,at=at.seq,labels=at.lab, rot = 0, tck=c(1,0))),
                 par.settings=list(strip.background = list( col="transparent" ),
                                   strip.border     = list( col = NA )))
  if(length(response)==1 & response!="all") PLOT <- useOuterStrips(PLOT,strip.left=FALSE) else PLOT <- useOuterStrips(PLOT)
  } else {
    PLOT <- update(PLOT, scales = list(y=list(relation=relation,alternating=3),
                                       x=list(alternating=1,at=at.seq,labels=at.lab, rot = 0, tck=c(1,0))),
                   par.settings=list(strip.background = list( col="transparent" ),
                                     strip.border     = list( col = NA ),
                                     strip.names=c(F,F),
                                     strip.levels=c(T,F)))
    #PLOT <- useOuterStrips(PLOT,strip.left=FALSE)
    PLOT <- update(PLOT, layout=c(2,1))
  }
PLOT
}

prepanel.default.xyplot2 <-
  function (x, y, type, subscripts, groups = NULL, upper1, lower1, ...)
  {
    if (any(!is.na(x)) && any(!is.na(y))) {
      ord <- order(as.numeric(x))
      if (!is.null(groups)) {
        gg <- groups[subscripts]
        dx <- unlist(lapply(split(as.numeric(x)[ord], gg[ord]),
                            diff))
        dy <- unlist(lapply(split(as.numeric(y)[ord], gg[ord]),
                            diff))
      }
      else {
        dx <- diff(as.numeric(x[ord]))
        dy <- diff(as.numeric(y[ord]))
      }
      list(xlim = lattice:::scale.limits(x), ylim = lattice:::scale.limits(c(upper1[subscripts],y[subscripts],lower1[subscripts])),
           dx = dx, dy = dy, xat = if (is.factor(x)) sort(unique(as.numeric(x))) else NULL,
           yat = if (is.factor(y)) sort(unique(as.numeric(y))) else NULL)
    }
    else prepanel.null()
  }