R/bootstrap.R

In Sondre91/STGARCH: STGARCH

#' Running boostrap experiments
#'
#' This function is used to run boostrap experiments.
#' @param true_param Numeric, Parameter vector, consiting of \eqn{\alpha_0} and \eqn{\alpha_1}. Default
#'     value draws random parameter vector well within the stationary region.
#' @param N Integer, number of Monte Carlo repetitions of bootstrap experiment.
#' @param M Integer, number of Boostrap replicas made per Monte Carlo repetition.
#' @param burnin Integer, size of temporal burnin
#' @param cutoff Integer, number of spatial points to be "cut off" in each spatial direction in order
#' to make the simulated process non-circular. If \code{cutoff=0}, the process is circular.
#' @param dim Numeric vector of length 3, indicating dimension size on form spatial1 x spatial2 x temporal
#' @param ncores Integer, number of parallel processors to be used.
#' @param n.iterations Integer, maximum number of iterations in the maximazation routine.
#' @param print Logical, indicating if reports should be printed.
#' @return The function returns a list consisting of summary statistics and the raw data estimates
#' @import doParallel
#' @import foreach
#' @export


running.bootstrap<-function(true_param=c(p1<-sample(1:50,1)/100,
                               sample(1:floor((.9-p1)*100/9),1)/100),
                  N=100, M=100, burnin=200,
                  cutoff=100, dim=c(30,30,200),
                  ncores=doParallel::detectCores()-1, n.iterations=30, print=TRUE){
  result<-matrix(NA_real_, ncol=4, nrow=N)
  W <- neighbourmatrix(dim[1], prod(dim[1:2]))
  tot.time<-system.time({
  for(i in 1:N){
    time<-system.time({
    result[i,]<-bootstrap.est(true_param=true_param, M=M, burnin=burnin,
                              cutoff=cutoff, dim=dim,
                              ncores=ncores, W = W, n.iterations=n.iterations)
    })
    if(print)
      cat("Finished with ",i, "iterations. It took ",time[3]," seconds.\n")
  }
  #Computing statistics
  rep.t.par<-rep(true_param,2)
  means<-apply(result[,1:4],2,mean)
  MSE<-rowMeans(apply(result[,1:4],1,function(x) return((x-rep.t.par)^2)))
  bias<-rowMeans(apply(result[,1:4],1,function(x) return((x-rep.t.par))))
  variance<-apply(result[,1:4],2,var)

  # Structuring output
  result.df<-data.frame(cbind(rep.t.par,means, bias, MSE, variance))
  rownames(result.df)<-c("est.0","est.1","bias.corr.0","bias.corr.1")
  names(result.df)<-c("true.param","mean.estimates", "mean.bias",
                      "MSE", "variance")
  retur_list<-list()
  retur_list$result.df<-result.df
  retur_list$raw<-result
  })
  if(print)
    cat("###################################################\n",
        "Finished with all in ", tot.time[3]/3600, " hours\n")
  return(retur_list)
}

#' Boostrap estimation function
#'
#' This function simulates a dataset considered to be the "truly observed" data, which are then estimated.
#' We imaging we do not know the true parameter, but want to estimate the bias and correct our estimate.
#' We therefore simulate \code{M} new datasets using the first estimate as the true parameters and estimate
#' the parameters for these as well. Here we know the true parameter value and can therefore estimate the bias.
#' We then use this bias estimate to correct the original estimate.
#'
#' @param true_param Numeric, Parameter vector, consiting of \eqn{\alpha_0} and \eqn{\alpha_1}.
#' @param M Integer, number of bootstrap replicas to be made
#' @param burnin Integer, size of temporal burnin
#' @param cutoff Integer, number of spatial points to be "cut off" in each spatial direction in order
#' to make the simulated process non-circular. If \code{cutoff=0}, the process is circular.
#' @param dim Numeric vector of length 3, indicating dimension size on form spatial1 x spatial2 x temporal
#' @param ncores Integer, number of parallel processors to be used.
#' @param W Neighbourhood matrix of size \code{prod(dim[1:2])x prod(dim[1:2])}
#' @param n.iterations Integer, maximum number of iterations in the maximazation routine.
#' @return Vector of the original estimates and the parametric bootstrap bias corrected estimates.
#' @import doParallel
#' @import foreach
#' @export
bootstrap.est<-function(true_param, M=100, burnin=200,
                        cutoff=100,dim=c(30,30,200),
                        ncores=32, W = neighbourmatrix(dim[1],prod(dim[1:2])), n.iterations=30){
    #Simulating "original" dataset
    data<- simulation_2d_arch(param=true_param, dim=dim, burnin=burnin, cutoff = cutoff)
    param.est<-est(X = data, n.iterations = n.iterations, W = W)
    doParallel::registerDoParallel(cores=ncores)
    value<-foreach::foreach(j = 1:M) %dopar%{
    # Simulating data
      data.boot<- simulation_2d_arch(param=param.est, dim=dim, burnin=burnin, cutoff = cutoff)
    # Estimation by maximum likelihood
      boot.est<-est(X = data.boot,n.iterations = n.iterations, W = W)
      names(boot.est)<-c("param1", "param2")
      boot.est
    }
    theta<-t(sapply(value, function(x) return(x)))
    theta_mean<-apply(theta,2,mean) # calculating mean of bootstrap estimates
    bias_corrected<-2*param.est-theta_mean # bias corrected estimate
  result<-c(param.est, bias_corrected) # vector to be returned
  names(result)<-c("est.0", "est.1","bias.corr.0","bias.corr.1")
  return(result)
}