R/bivar.R

In mboeck11/BTSM: Bayesian Time Series Models

#' @name bivar
#' @title Bayesian Interacted Vector Autoregression
#' @usage bivar(Data, plag=1,draws=5000,burnin=5000,prior="NG",SV=TRUE,h=0,thin=1,
#'              hyperpara=NULL,eigen=FALSE,
#'              Ex=NULL,cons=FALSE,trend=FALSE,applyfun=NULL,cores=NULL,verbose=TRUE)
#' @param Data Data in matrix form.
#' @param Idata Data for interaction term in matrix form.
#' @param plag number of lags
#' @param draws number of saved draws.
#' @param burnin number of burn-ins.
#' @param prior which prior
#' @param SV If set to \code{TRUE} stochastic volatility is enabled.
#' @param h holdout-sample.
#' @param thin thinning factor
#' @param hyperpara hyperparameter set
#' @param eigen should eigenvalues be computed?
#' @param Ex exogenous variables to add to the model
#' @param cons If set to \code{TRUE} a constant is included.
#' @param trend If set to \code{TRUE} a trend is included.
#' @param applyfun parallelization
#' @param cores number of cores
#' @param verbose verbosity option
#' @export
#' @importFrom MASS ginv
#' @importFrom methods is
#' @importFrom GIGrvg rgig
#' @importFrom parallel parLapply mclapply
#' @importFrom Rcpp evalCpp
#' @importFrom stats is.ts median time ts
#' @importFrom xts is.xts
#' @importFrom zoo coredata
bivar<-function(Data,Idata,plag=1,draws=5000,burnin=5000,prior="NG",SV=TRUE,h=0,thin=1,hyperpara=NULL,eigen=TRUE,Ex=NULL,cons=TRUE,trend=FALSE,applyfun=NULL,cores=NULL,verbose=TRUE){
  start.bivar <- Sys.time()
  #--------------------------------- checks  ------------------------------------------------------#
  if(!is.matrix(Data)){
    stop("Please provide the argument 'Data' either as 'matrix' object.")
  }
  if(!is.matrix(Idata)){
    stop("Please provide the argument 'Idata' either as 'matrix' object.")
  }
  if(!is.null(Ex)){
    if(!is.list(Ex) & !is.matrix(Ex)){
      stop("Please provide the argument 'Ex' either as 'list' or as 'matrix' object.")
    }
  }
  if(!is.numeric(plag)){
    stop("Please specify number of lags as numeric.")
  }
  if(any(is.na(plag))){
    stop("Please specify number of lags.")
  }
  if(length(plag)>1 || plag<1){
    stop("Please specify number of lags accordingly. One lag length parameter for the whole model.")
  }
  if(!is.numeric(draws) | !is.numeric(burnin)){
    stop("Please specify number of draws and burnin as numeric.")
  }
  if(length(draws)>1 || draws<0 || length(burnin)>1 || burnin<0){
    stop("Please specify number of draws and burnin accordingly. One draws and burnin parameter for the whole model.")
  }
  if(!prior%in%c("MN","SSVS","NG")){
    stop("Please choose an available prior specification.")
  }
  #-------------------------- construct arglist ----------------------------------------------------#
  args <- .construct.arglist(bvar)
  if(verbose){
    cat("\nStart estimation of Bayesian Interacted Vector Autoregression.\n\n")
    cat(paste("Prior: ",ifelse(prior=="MN","Minnesota prior",ifelse(prior=="SSVS","Stochastic Search Variable Selection prior",ifelse(prior=="NG","Normal-Gamma prior","Natural conjugate prior"))),".\n",sep=""))
    cat(paste("Lag order: ",plag,"\n",sep=""))
    cat(paste("Stochastic volatility: ", ifelse(SV,"enabled","disabled"),".\n",sep=""))
  }
  #------------------------------ user checks  ---------------------------------------------------#
  # check Data
  if(is.matrix(Data)){
    if(any(is.na(Data))){
      stop("The data you have submitted contains NAs. Please check the data.")
    }
    isTS  <- is.ts(Data)
    isXTS <- is.xts(Data)
    Traw  <- nrow(Data)
    if(isTS || isXTS){
      temp       <- as.character(time(Data))
      years      <- unique(regmatches(temp,regexpr("^[0-9]{4}",temp)))
      months     <- temp
      for(kk in 1:length(years)) months <- gsub(paste(years[kk],"(\\.)?",sep=""),"",months)
      freq       <- length(unique(months))
      months     <- strtrim(months,3)
      startmonth <- ifelse(months[1]=="","01",ifelse(months[1]=="083","02",ifelse(months[1]=="166","03",ifelse(months[1]=="25","04",
                                                                                                               ifelse(months[1]=="333","05",ifelse(months[1]=="416","06",ifelse(months[1]=="5","07",ifelse(months[1]=="583","08",
                                                                                                                                                                                                           ifelse(months[1]=="666","09",ifelse(months[1]=="75","10",ifelse(months[1]=="833","11","12")))))))))))
      timeindex  <- seq.Date(from=as.Date(paste(years[1],"-",startmonth,"-01",sep=""), format="%Y-%m-%d"),
                             by=ifelse(freq==12,"months","quarter"), length.out = Traw)
      Data       <- ts(coredata(Data),  start=c(as.numeric(years[1]),as.numeric(startmonth)),frequency=freq)
      Idata      <- ts(coredata(Idata), start=c(as.numeric(years[1]),as.numeric(startmonth)),frequency=freq)
    }else{
      timeindex  <- seq.Date(from=as.Date("1830-08-01", format="%Y-%m-%d"), by="month", length.out = Traw)
      Data       <- ts(Data,  start=c(1830,8), frequency=12)
      Idata      <- ts(Idata, start=c(1830,8), frequency=12)
    }
    args$time <- timeindex
    args$Traw <- length(timeindex)
  }
  # check interactdion data
  if(nrow(Idata)!=nrow(Data)){
    stop("Please provide Data and Idata with same length.")
  }
  # check truly exogenous variables
  if(!is.null(Ex)){
    if(is.matrix(Ex)){
      if(any(is.na(Ex))){
        stop("The data for exogenous variables you have submitted contains NAs. Please check the data.")
      }
      if(nrow(Ex)!=args$Traw){
        stop("Provided data and truly exogenous data not equally long. Please check.")
      }
    }
  }
  # check thinning factor
  if(thin<1){
    thin_mess <- paste("Thinning factor of ",thin," not possible. Adjusted to ",round(1/thin,2),".\n",sep="")
    thin <- round(1/thin,2)
  }
  if(draws%%thin!=0){
    thin_mess <- paste("Thinning factor of ",thin," no divisor of ",draws," (number of draws to save for posterior analysis).\n",sep="")
    div <- .divisors(draws,thin)
    thin <- min(div[which(abs(div-thin)==min(abs(div-thin)))])
    thin_mess <- paste(thin_mess,"New thinning factor: ", thin,". This means every", ifelse(thin==1,"",ifelse(thin==2,paste(" ",thin,"nd ",sep=""), ifelse(thin==3,paste(" ",thin,"rd ",sep=""),paste(" ",thin,"th ",sep="")))), "draw is saved.\n",sep="")
  }else{
    thin_mess <- paste("Thinning factor: ", thin,". This means every ",ifelse(thin==1,"",ifelse(thin==2,paste(thin,"nd ",sep=""),ifelse(thin==3,paste(thin,"rd ",sep=""),paste(thin,"th ",sep="")))),"draw is saved.\n",sep="")
  }
  if(verbose) cat(thin_mess)
  args$thindraws <- draws/thin
  # set default
  if(verbose) cat("Hyperparameter setup: \n")
  default_hyperpara <- list(c=10, Multiplier=10, # hyperparameter setup for natural conjugate case
                            a_1=0.01,b_1=0.01, prmean=1,# Gamma hyperparameter SIGMA (homoskedastic case) and mean
                            Bsigma=1, a0=25, b0=1.5, bmu=0, Bmu=100^2, # SV hyper parameter
                            shrink1=0.1,shrink2=0.2,shrink3=10^2,shrink4=0.1, # MN
                            tau0=0.1,tau1=3,kappa0=0.1,kappa1=3,p_i=0.5,q_ij=0.5,   # SSVS
                            e_lambda=0.01,d_lambda=0.01,a_start=0.7,sample_A=FALSE,a_log=TRUE,
                            use_R=FALSE) # NG
  paras     <- c("c","a_1","b_1","prmean","Bsigma_sv","a0","b0","bmu","Bmu","shrink1","shrink2","shrink3","shrink4",
                 "shrink4","tau0","tau1","kappa0","kappa1","p_i","q_ij","e_lambda","d_lambda","a_log","a_start","sample_A","use_R")
  if(is.null(hyperpara)){
    if(verbose) cat("\t No hyperparameters are chosen, default setting applied.\n")
  }
  if(!is.null(hyperpara)){
    for(para in names(hyperpara)){
      if(!para%in%paras){
        warning(paste0(para," no valid hyperparameter. Please check.\n"))
        next
      }
      default_hyperpara[para] <- hyperpara[para]
      if(para=="a_start") default_hyperpara["a_log"] <- FALSE
    }
    if(verbose) cat("Default values for chosen hyperparamters overwritten.\n")
  }
  #----------------------------------transform to matirx--------------------------------------------------------#
  Yfull <- as.matrix(Data)
  Dfull <- as.matrix(Idata)
  #---------------------------------hold out sample------------------------------------------------------------#
  args$yfull <- Yfull
  Yraw       <- Yfull[1:(nrow(Yfull)-h),,drop=FALSE]
  Draw       <- Dfull[1:(nrow(Dfull)-h),,drop=FALSE]
  args$time  <- args$time[1:(length(args$time)-h)]
  #------------------------------ prepare applyfun --------------------------------------------------------#
  if(is.null(applyfun)) {
    applyfun <- if(is.null(cores)) {
      lapply
    } else {
      if(.Platform$OS.type == "windows") {
        cl_cores <- parallel::makeCluster(cores)
        on.exit(parallel::stopCluster(cl_cores))
        function(X, FUN, ...) parallel::parLapply(cl = cl_cores, X, FUN, ...)
      } else {
        function(X, FUN, ...) parallel::mclapply(X, FUN, ..., mc.cores =
                                                   cores)
      }
    }
  }
  if(is.null(cores)) {cores <- 1}
  #--------------------------- checking eigenvalues ----------------------------------------------------------#
  if(is.logical(eigen)){
    if(eigen){trim<-1.00}else{trim<-NULL}
  }else{
    trim<-eigen;eigen<-TRUE
  }
  #------------------------------ estimate BVAR ---------------------------------------------------------------#
  if(verbose) cat("\nEstimation of model starts...\n")
  globalpost <- .BIVAR_linear_wrapper(Yraw=Yraw,Draw=Draw,prior=prior,plag=plag,draws=draws,burnin=burnin,cons=cons,trend=trend,SV=SV,thin=thin,default_hyperpara=default_hyperpara,Ex=Ex,applyfun=applyfun,cores=cores,eigen=eigen,trim=trim)
  args$Y <- globalpost$Y
  args$X <- globalpost$X
  args$D <- globalpost$D
  args$thindraws <- ifelse(eigen,sum(globalpost$store$A.eigen<trim),args$thindraws)
  if(verbose) cat(paste0("\nModel yields ", args$thindraws," stable draws."))
  #---------------------- return output ---------------------------------------------------------------------------#
  out  <- structure(list("args"=args,
                         "xglobal"=Yfull,
                         "xint"=Dfull,
                         "post"=globalpost$post,
                         "store"=globalpost$store), class = "bivar")
  end.bivar <- Sys.time()
  diff.bivar <- difftime(end.bivar,start.bivar,units="mins")
  mins.bivar <- round(diff.bivar,0); secs.bivar <- round((diff.bivar-floor(diff.bivar))*60,0)
  if(verbose) cat(paste("\nNeeded time for estimation of bvar: ",mins.bivar," ",ifelse(mins.bivar==1,"min","mins")," ",secs.bivar, " ",ifelse(secs.bivar==1,"second.","seconds.\n"),sep=""))
  return(out)
}

#' @method print bivar
#' @export
#' @importFrom utils object.size
print.bivar<-function(x, ...){
  cat("---------------------------------------------------------------------------------------")
  cat("\n")
  cat("Model Info:")
  cat("\n")
  prior <- ifelse(x$args$prior=="MN","Minnesota Prior",ifelse(x$args$prior=="SSVS","Stochastic Search Variable Selection Prior",ifelse(x$args$prior=="NG","Normal Gamma Prior","not defined.")))
  cat(paste("Prior: ",prior,sep=""))
  cat("\n")
  cat(paste("Nr. of lags: ",x$args$plag,sep=""))
  cat("\n")
  cat(paste("Nr. of posterior draws: ",x$args$draws,"/",x$args$thin,"=",floor(x$args$draws/x$args$thin),sep=""))
  cat("\n")
  cat(paste("Size of BVAR object: ",format(object.size(x),units="MB"),sep=""))
  cat("\n")
  if(x$args$eigen){
    cat(paste("Model has ",sum(x$post$A.eigen<1)," stable draws.",sep=""))
    cat("\n")
    cat("---------------------------------------------------------------------------------------")
  }
  cat("\n")
  cat("Model specification:")
  cat("\n")
  cat(colnames(x$args$yfull))
  cat("\n")
  invisible(x)
}

#' @name logLik
#' @title Extract Log-likelihood of Bayesian VAR
#' @description Extract Log-Likelihood for \code{bivar}.
#' @param object an object of class \code{bivar}.
#' @param ... additional arguments.
#' @param quantile reported quantiles. Default is set to median.
#' @return Returns an vector of dimension \code{q} (number of specified quantiles) of log-likelihoods.
#' @importFrom stats logLik
#' @export
logLik.bivar<-function(object, ..., quantile=.50){
  if(length(quantile)!=1){
    stop("Please provide only one quantile.")
  }
  temp <- object$args$logLik
  xglobal <- object$xglobal
  bigT <- nrow(xglobal)
  bigK <- ncol(xglobal)
  plag <- object$args$plag
  if(is.null(temp)){
    cons      <- object$args$cons
    trend     <- object$args$trend
    thindraws <- object$args$thindraws
    X_large   <- .mlag(xglobal,plag)
    if(cons)  X_large <- cbind(X_large,1)
    if(trend) X_large <- cbind(X_large,seq(1:bigT))
    Y_large   <- xglobal[(plag+1):bigT,,drop=FALSE]
    X_large   <- X_large[(plag+1):bigT,,drop=FALSE]
    A_large   <- object$store$A_store
    S_large   <- object$store$Smed_store
    loglik <- try(loglik_C(Y_in=Y_large,X_in=X_large,A_in=A_large,S_in=S_large,thindraws=thindraws)$loglik,silent=TRUE)

    if(is(loglik,"try-error")){
      out <- -Inf
    }else{
      out <- quantile(loglik,quantile,na.rm=TRUE)
    }
    eval.parent(substitute(object$args$logLik<-out))
  }else{
    out <- temp
  }
  attributes(out) <- list(nall=bigT, nobs=bigT, df=bigK)
  class(out) <- "logLik"
  return(out)
}