stresstest.Rcheck/00_pkg_src/stresstest/R/stress_interface.R

In gamalamboy/stresstest: Implement the ETI Stress Testing Simulation

#' Low-level interface to retrieve and convert .csv files into a zoo object with an arbitrary transformation.
#' 
#' This function takes an input path, reads the file with \code{read.csv},
#' and creates a \code{zoo} object with the specified time sequence.
#' 
#' @details The function only accepts data with either numeric, logical, or
#' integer values. It also adds an extra column "i0" to be used for
#' constructing the model intercepts. If there are more end years than data,
#' the function pads the additional rows with NA's. This interface is called by
#' the forecasting functions when creating the Bayesian samples.
#' 
#' @param inputfile a character string indicating the path of the file to be
#' read.
#' @param start the start date. Should be inserted as in \code{ts}.
#' @param end the end date.
#' @param freq the frequency for the time interval.
#' @param transf an optional transformation to be given to the data before
#' creating the zoo object.
#' 
#' @return a zoo object.
#' 
#' @rdname DataInput
#' 
#' @export
.DataInput <- function(inputfile, start = 1992, end = 2016, freq = 1, 
transf = NULL)
{

    data <- utils::read.csv( inputfile, stringsAsFactors = FALSE,
    na.strings = "")
    
    if( end - start > nrow(data))
    {
        m <- matrix( NA, nrow = end - start - nrow(data) + 1, 
        ncol = ncol(data))
        colnames(m) <- colnames(data)
        data <- rbind( data, m)
    }
    
    data <- cbind( data, i0 = 1)
    
    cl <- vapply( data, class, character(1))
    
    if( any( !cl %in% c("numeric", "logical", "integer")))
    {
        warning("Non-Numeric Data Found!")
    }
    
    index1 <- seq( start, end, freq)
    zut <- zoo::zoo( data[,colnames(data) != "year"], order.by = index1)
    
    if( !is.null( transf))
    {
        zut <- transf(zut)
    }
    
    zut
}

#' Extract formula list from text file
#' 
#' This is an input function that retrieves the model from a text file,
#' and converts it to a list of formula objects.
#' 
#' @details as the function to interpret the model structure
#' (\code{Ztprops}) is slightly different from the base R functionality, a
#' particular class was given to ensure that the intercepts and terms are
#' interpreted consistently in other functions. This is planned to be
#' homogeneized with the R base interface sometime soon, by implementing
#' these functions as methods.
#' 
#' @param forfile a character string, indicating the formula path to be
#' read.
#' 
#' @return a list with additional class \code{ForList}.
#' 
#' @rdname FormulaInput
#' @export
.FormulaInput <- function(forfile = NULL)
{
    if( !file.exists(forfile)) stop( gsub( "\n", "", "Please input a Text 
    File with the appropriate model structure for evaluation."))
    
    dt1 <- gsub( "=", "~", scan( forfile, what = character(), sep = "\n"))
    
    dt1 <- lapply( dt1, function(y) 
    {
        tmp <- stats::as.formula(y)
        if( attr( stats::terms.formula( tmp), "intercept") == 1)
        {
            tmp <- stats::update.formula( tmp, ~ . + i0)
        }
        
        class(tmp) <- c("formula", "ForList")        
        tmp 
    })
    
    names(dt1) <- vapply( dt1, function(y) all.vars(y)[[1L]], "")
    class(dt1) <- c("list", "ForList")
    dt1
}

.centerT <- function(phrase, char = "#", termlen = 60, pad = 5)
{
    
    phrase_wrap <- do.call(c, strsplit(phrase, "\n"))
    
    
    lv <- list()
    
    for(y in 1:length(phrase_wrap))
    {
        ln <- floor( (termlen - nchar(phrase_wrap[[y]]) - 2) / 2)
        
        padv <- paste( rep( char, ln),collapse = "")
        cat( paste(padv, phrase_wrap[[y]], padv, sep = " "), "\n")
    }
    lv
}

#' Interface for specifying the paths for the variables
#' 
#' This function is run at startup, when the environment is being prepared. It
#' basically presents the user with a menu screen to specify each of the
#' required files for the model to work.
#' 
#' @details this function allows the user to select the terminal length, as
#' well as the level of "padding" the centered text has. The idea is that
#' on different monitors, this may need to be changed.
#' 
#' @param ... parameters to be passed to .centerT: \describe{
#' \item{char}{ the character used for developing the menu structure. Could be
#' useful in specialized consoles.}
#' \item{termlen}{the length of the interface in characters. Defaults to 60.}
#' \item{pad}{the level of padding that will be added to larger text lines
#' that are split in two.}
#' }
#' @return the list of file paths, as required by the model functions.
#' @export
SelectPaths <- function(...)
{

    compnames <- c("forfile", "inputfile", "priorbeta", "priorvbeta", 
    "workdir", "fundir", "resultsdir")
    
    components <- rep("", 5)
    
    compselect <- c("the formula file:\n", "the data file:\n",
    "the file containing the parameter prior means (if any):\n",
    "the file containing the parameter prior variances (if any):\n",
    "the work directory:\n")
    
    compselect <- paste( " ", 1:length(compselect), ") ", compselect, 
    sep = "")
    
    MainMenu <- function()
    {
        .centerT( "Path Selection Menu", ...)
        .centerT( "Please Select the Path You Wish to Specify:", char = "#")
        cat( "\n", compselect, "To quit, type \"6\".\n", sep = "")
    }
    
    MainMenu()
    
    selection <- readline()
    seln <- as.numeric(selection)
    tmp <- ""
    
    while( selection != "6" | components[5] == "")
    {
        while( !seln %in% 1:6)
        {
            cat("Selection must be a number between 1 and 6.\n")
            selection <- readline()
            seln <- as.numeric(selection)
        }
        
        if( seln %in% 1:4)
        {
            tmp <- tcltk::tk_choose.files(multi = FALSE)
        } else tmp <- tcltk::tk_choose.dir( caption = 
        "Please Select the Work Directory")
        
        if( length(tmp) == 0 | is.na(tmp)) tmp <- ""
        
        components[seln] <- tmp
        
        if( components[seln] != "")
        {
            compselect[seln] <- gsub(":\\\n$", ": [Selected]\n", compselect[
            seln])
        }
        
        if( selection == "6" & components[5] == "")
        {
            cat("Cannot close until work directory is specified.\n")
        }
        
        system("clear")
        
        MainMenu()
        selection <- readline()
        seln <- as.numeric(selection)
    }
    
    components[[6]] <- paste( components[[5]], "R", sep = "/")
    components[[7]] <- paste( components[[5]], "/Stress Scenarios ",
    format( Sys.Date(), "%Y-%m-%d"), sep = "")
    
    names(components) <- compnames
    
    return(as.list(components))
}

#' Select the prior parameters
#' 
#' This function calls the internal functions \code{.priorBeta},
#' \code{.priorVbeta}, and \code{.priorVar} to create the prior values,
#' whilst assigning them their proper formal arguments under a \code{tryCatch}
#' handler. Errors are stored in the global environment.
#' 
#' @param data the model data, used when constructing the parameter covariance
#' matrix.
#' @param ztprop1 the model structure (a \code{Ztprops} object): required for
#' all three priors.
#' @param ... additional parameters required by the priors.
#' 
#' @return a list with additional class \code{PriorList}, containing:
#' \enumerate{
#' \item the prior parameter means.
#' \item the prior parameter precision matrix.
#' \item the model covariance and precision matrices (\eqn{\Sigma} and
#' \eqn{\Sigma^{-1}})
#' }
#' 
#' @rdname SelectPriors
#' @export
.SelectPriors <- function( data, ztprop1, ...)
{
#     if(!inherits( ztprop1, "Ztprops")) stop( 
#     "formula_list must be of class \"Ztprops\"")
    
    plist <- c( list( data = data, ztprop1 = ztprop1), list(...))
    
    funs <- list( .priorBeta, .priorVbeta, .priorVar)
    names(funs) <- paste("Prior", c("Beta", "Vbeta", "Var"), sep = "")
    
    ilist <- lapply( funs, function(v)
    {
        pv <- stats::na.omit( match( methods::formalArgs(v), names(plist)))
        do.call(v, plist[pv])
        
    })
    
    names(ilist) <- names(funs)
    
    ilist$InvPriorVbeta <- solve( ilist$PriorVbeta)
    ilist$InvVar <- solve( ilist$PriorVar)
    
    ilist <- ilist[ c("PriorBeta", "InvPriorVbeta", "PriorVar", "InvVar")]
    
    class( ilist) <- append( class(ilist), "PriorList")
    ilist
}

#' Select prior parameters from file
#' 
#' This is a convenience function for specifying priors from a file and
#' preparing the \code{ilist} object for the Gibbs sampler.
#' 
#' @param PATHS the PATHS variable; see 'SelectPaths'.
#' @param select a vector with entries \emph{str}, \emph{scenariob}, and
#' \emph{scenariov}, specifying the structure and mean/variance scenarios
#' from which to construct the \code{ilist} object.
#' @param ... additional parameters to be passed to \code{.DataInput}.
#' 
#' @return a list with additional class \code{PriorList}, containing:
#' \enumerate{
#' \item the prior parameter mean.
#' \item the prior parameter precision matrix.
#' \item the model's covariance and precision matrices (\emph{PriorVar} and
#' \emph{InvVar}).
#' }
#' 
#' @export
PriorsFromFile <- function(PATHS, select, ...)
{
    zt1 <- Ztprops( .FormulaInput( PATHS$forfile))
    dats <- .DataInput( inputfile = PATHS$inputfile, ...)
    .pbeta <- utils::read.csv( PATHS$priorbeta, stringsAsFactors = FALSE)
    .pvbeta <- utils::read.csv( PATHS$priorvbeta, stringsAsFactors = FALSE)
    

    base_file <- do.call(rbind, lapply(zt1, function(x) x$props))
    
    pbeta <- .pbeta[ .pbeta$structure == select[["str"]], -1]
    pvbeta <- .pvbeta[ .pvbeta$structure == select[["str"]], -1]
    
    pbeta <- merge( base_file, pbeta, sort = FALSE)
    
    bval <- stats::setNames( pbeta[, select[["scenariob"]]], pbeta$terms)
    bval <- .priorBeta(zt1, bval)
    
    pvb <- as.list(pvbeta[ pvbeta$Parameter != "var", select[["scenariov"]]])
    
    names(pvb) <- pvbeta$Parameter[ pvbeta$Parameter != "var"]
    
    pvb <- .priorVbeta( data = dats, ztprop1 = zt1, params = pvb)
    
    vl <- pvbeta[ pvbeta$Parameter == "var", select[["scenariov"]] ]
    
    ilist <- list( PriorBeta = bval, PriorVbeta = pvb,
    PriorVar = .priorVar(  zt1, value = vl))
    
    ilist$InvPriorVbeta <- solve( ilist$PriorVbeta)
    ilist$InvVar <- solve( ilist$PriorVar)
    
    ilist <- ilist[ c("PriorBeta", "InvPriorVbeta", "PriorVar", "InvVar")]
    
    class( ilist) <- append( class(ilist), "PriorList")
    ilist

    
}

#' High-Level interface to create the necessary objects for the Gibbs sampler
#' 
#' This function calls a series of helper functions to create the model
#' objects to run the Gibbs sampler.
#' 
#' @details This function is the main workhorse to 'set the stage' for the
#' Gibbs sampler. All required objects are placed inside a list and returned
#' as inputs.
#' 
#' @param PATHS a list of file paths. This function is intended to be used
#' alongside a GUI that allows the User to select the file paths for each
#' component. The following paths are required in the list: 
#' \describe{
#' \item{forfile}{the formula file.}
#' \item{fundir}{the directory containing the functions.}
#' \item{workdir}{the directory where the results and intermediate inputs are 
#' to be stored.}
#' \item{inputfile}{path to the basic input file (a \code{.csv}).}
#' \item{resultsdir}{the path of the results directory, where the sampler
#' output is to be stored.}
#' \item{priorbeta}{the path of the prior parameters to be used in the
#' simulation.}
#' \item{vbeta}{the path for the prior parameters to the covariance matrix of
#' the parameters.}
#' }
#' @param priors an optional \code{ilist} object. If NULL, the .SelectPriors 
#' interface will be called to enter the priors; otherwise, it will use the
#' \code{ilist} supplied by the user.
#' 
#' @param ... additional parameters to be passed to \code{.DataInput} and
#' \code{ .SelectPriors}.
#' 
#' @return a list with the following components: \describe{
#' \item{dats}{the data set (a \code{zoo} object).}
#' \item{form}{a \code{ForList} object containing the formulas.}
#' \item{Zprop}{a \code{Ztprops} object containing the model structure.}
#' \item{Z}{the model matrices in \code{Zlist} form}
#' \item{Y}{the model dependent vectors in \code{Ylist} form}
#' \item{ilist}{the model prior list.}}
#' 
#' @export
CreateData <- function(PATHS, priors = NULL, ...)
{
    
    dats <- .DataInput( inputfile = PATHS$inputfile, ...)

    form <- .FormulaInput( PATHS$forfile)
    Zprop <- Ztprops(form)
    
    Zlst <- create.Z( dats, Zprop)
    
    Ylst <- create.Y( dats, Zprop)
    
    if( is.null(priors)) ilist <- .SelectPriors( dats, Zprop, ...) else { 
    ilist <- priors }
    
#     save( ilist, file = paste( PATHS$resultsdir, "/prior_params", length(
#     list.files( PATHS$resultsdir)), ".rds", sep = ""))
    
    list( dats = dats, form = form, Zprop = Zprop, Z = Zlst, Y = Ylst, 
    ilist = ilist)
}

#' Run a Gibbs Simulation
#' 
#' This is the main user function, designed to obtain posterior samples from
#' the Normal-Wishart VAR model.
#' 
#' @details The main underlying functions are \code{CreateData}, which is used
#' to set the environment, and Gibbs, which runs the Gibbs sampler. As a
#' safety measure, one cannot run the Gibbs sampler without pre-specifying the
#' \emph{base_objs} argument; the function will return a pre-made 
#' \emph{base_objs} for revision, and must be run with this environment for
#' it to begin the sampling process.
#' 
#' @param PATHS the PATHS variable to be used.
#' @param iters the number of Gibbs iterations.
#' @param base_objs an optional parameter providing the environment from 
#' \code{CreateData}. If NULL, the \code{dataentry} interface for .SelectPriors
#' is activated.
#' @param burn the number of burn-in samples to discard from the Gibbs Sampler.
#' @param ... additional parameters to be passed to \code{CreateData}.
#' 
#' @return if base_objs is NULL, it will return the base environment for
#' revision. If not, then a \code{Gibbs} object.
#' @export
RunScenario <- function( PATHS, iters = 40000, base_objs = NULL, burn = 10000,
                         ...)
{
    if( is.null(base_objs))
    {
        return( CreateData(PATHS, ...))
    } else 
    {
        res1 <- Gibbs( base_objs$Z, base_objs$Y, base_objs$ilist, iters)
        rownames(res1$Beta) <- rownames(base_objs$ilist$PriorBeta)
        dimnames(res1$Vbeta)[1:2] <- dimnames(base_objs$ilist$InvPriorVbeta)
        dimnames(res1$Var)[1:2] <- dimnames(res1$InvVar) <- dimnames(
        base_objs$ilist$PriorVar)
        
        BurnIn(res1, burn)
    }
}

#' Eliminate burn-in samples from Gibbs
#' 
#' This function returns the Gibbs object minus a randomized
#' burn period.
#' 
#' @param Gibbs an object of class \code{Gibbs}.
#' @param burn the burn-in sample size.
#' 
#' @return a Gibbs object without the burn-in samples.
BurnIn <- function(Gibbs, burn)
{
    GibbsC <- c( list( Gibbs$Beta[,-c(1:burn)]),
    lapply( Gibbs[-1], function(x) x[,,-c(1:burn)]))
    
    names(GibbsC) <- names( Gibbs)
    
    GibbsC
    
}