stresstest.Rcheck/00_pkg_src/stresstest/blah/stress_exec_2017.R

In gamalamboy/stresstest: Implement the ETI Stress Testing Simulation

STR <- "str1"
SCE <- 1

SCES <- expand.grid( c("Base", "Optimist", "Pessimist", "Risk"), 
c("Hi_Var", "Low_Var"), stringsAsFactors = FALSE, KEEP.OUT.ATTRS = FALSE)
SCES <- as.matrix(SCES)
colnames(SCES) <- c("scenariob", "scenariov")

PATHS <- "./blah"
PATHS <- list(
forfile = paste0( PATHS, "/eqs_orig/eqs_econ_str1"),
fundir = PATHS,
workdir = PATHS,
inputfile = paste(PATHS, "stress.csv", sep = "/"),
resultsdir = paste( PATHS, "/Stress Scenarios ", format( Sys.Date(), 
"%Y-%m-%d"), sep = ""),
priorbeta = paste( PATHS, "/ilist_beta20170707.csv", sep = ""),
priorvbeta = paste( PATHS, "/ilist_vbeta.csv", sep = ""))

# source(paste( PATHS$workdir, "/stress_interface.R", sep = ""))
# execF <- interfaceF$EnvPrep(PATHS)

set.seed( 1000)
library(devtools)
load_all()

# Rcpp::sourceCpp("~/Documents/ETI/stress/src/draft_bayes.cpp", 
# rebuild = TRUE, cleanupCacheDir = TRUE)

# library(stress, lib.loc="~/.R")

# execF <- EnvPrep()

pt <- diff

# for( i in 1:nrow(SCES))
# {
#     ilist <- PriorsFromFile( select = c(str = STR, SCES[1,]))
if( !exists("boj0"))
{
    boj0 <- RunScenario( PATHS, base_objs = NULL)
    boj <- RunScenario( PATHS, base_objs = boj0, transf = pt)
}

# source(paste0(PATHS$workdir, "/fixed.R"))
# 
zp <- Ztprops( boj0$form)
# pv4 <- .create.zmt( boj0$dats, pv3)
Z <- create.Z( boj0$dats, zp)
Y <- create.Y( boj0$dats, zp)

Dv <- .DataInput( PATHS$inputfile, end = 2020, transf = pt)

# pv <- with(boj0, Ofore(Z[,,1], ilist$PriorBeta, ilist$PriorVar))
pv2 <- ForecastGibbs( Dv, 2017:2020, boj0$Zprop, boj$Beta, 
                      boj$Var)

# 
#     pv <- St( boj0$ilist$PriorBeta, boj0$ilist$PriorVar, Y, Z)
# pv2 <- St( boj0$ilist$PriorBeta, boj0$ilist$PriorVar, Y, Z)
# pv2c <- RWishart( 35, pv2)
# 
# pv1 <- Vbeta( boj0$ilist$InvPriorVbeta, pv2c, Z)
# 
# pv0 <- Beta( boj0$ilist$PriorBeta, boj0$ilist$InvPriorVbeta, 
# pv1, pv2c, Y, Z)
# 
# pvg <- Gibbs( Z, Y, boj0$ilist, 20)

# UpdateZ <- function( data, ztprop1, ofore)
# {
#     yr <- max(index(data)) + 1
#     
#     exovars <- setdiff( ztprop1$vars_indep, ztprop1$vars_dep)
#     
#     return(exovars)
#     
#     base1 <- matrix(NA, 1, ncol(data), dimnames = list(yr, colnames(data)))
#     base1[,colnames(ofore), drop = FALSE] <- coredata(ofore)
# }

# pvu <- UpdateZ( boj0$dats, pv3, Z[,,1] %*% pvg$Beta[,1])

#     
#     boj <- RunScenario( PATHS, base_objs = boj0, transf = pt, iters = 1)
#     boj0 <- Gibbs( boj0$Z, boj0$Y, boj0$ilist, 1)
#     boj <- NoCorSamp( boj, burn = 15000)
#     save( boj, file = paste0(PATHS$resultsdir, "/samp_", STR, "_", i, ".rda"))
# }
# # q(save = "no")

# if( !exists("lfore"))
# {
#     lfore <- lapply( c("opt", "base", "pes", "risk"), function(v)
#     {
#        lfore <- ForecastGibbs(tv, boj0$Zprop, boj, 2017:2020, 
#        pt)
#     })
#     
#     lfore <- array( unlist(lfore), dim = dim( lfore[[1]]) * c(1,1,4), 
#     dimnames = c( dimnames( lfore[[1]])[1:2], list(NULL)))
# }




# 
# ppb <- paste("priors/priorlist1_", STR, ".rda", sep = "")
# if( !dir.exists( paste( PATHS$resultsdir, "priors", sep = "/"))) dir.create(
# paste( PATHS$resultsdir, "priors", sep = "/"))
# 
# if( ppb %in% list.files( PATHS$resultsdir))
# {
#     load(paste( PATHS$resultsdir,ppb, sep = "/"))
# } else
# {
#     .pbeta <- read.csv( PATHS$priorbeta, stringsAsFactors = FALSE)
#     .pvbeta <- read.csv( PATHS$priorvbeta, stringsAsFactors = FALSE)
#     priorlist <- expand.grid( c("Base", "Optimist", "Pessimist", "Risk"), 
#     c("Hi_Var", "Low_Var"))
# 
#     colnames(priorlist) <- c("scenariob", "scenariov")
# 
#     priorlist1 <- apply( priorlist, 1, function(x) as.list( c( str = STR, x)))
# 
#     priorlist1 <- lapply( priorlist1, function(x) interfaceF$PriorsFromFile( 
#     PATHS, select = x, transf = pt))
#     
#     names( priorlist1) <- apply( priorlist, 1, function(x) paste( x, STR,
#     collapse = "_", sep = "_"))
#     
#     save( priorlist1, file = paste( PATHS$resultsdir, ppb, sep = "/"))
# }
# 
# if( !exists( "kb"))
# {
#     boj0 <- RunScenario( PATHS, base_objs = NULL, priors = priorlist1[[SCE]])
#     boj <- RunScenario( PATHS, base_objs = boj0, transf = pt, iters = 40000)
#     boj <- interfaceF$NoCorSamp( boj, burn = 15000)
# 
#     kb <- ForecastGibbs( inputfile = PATHS$inputfile, ztprop1 = boj0$Zprop,
#     Gibbs = boj, period = 2017:2020, transf = pt, mean = FALSE)
#     
# }
# 
# save( kb, file = paste( PATHS$resultsdir, "/vars_filtered_", STR, "_", SCE, 
# ".rda", sep = ""))
# 
# print( dim(kb))
# quit( save = "no")



# kb <- apply(kb, 2, cumsum)

# write.csv( kb, "tmp_results.csv", row.names = FALSE)


# ppb <- merge( ppv, ppb[ ppb$Structure == STR,], by = c("vars_dep", 
# "vars_indep", "terms"),sort = FALSE)




# if(!exists( "sce1")) 
# {
#     sce1 <- RunScenario(PATHS, iters = 100, transf = function(x) diff(
#     log(x + 10)))
# #     sce1 <- interfaceF$NoCorSamp( sce1)
# }



# Convergence Diagnostics

# interfaceF$BetaConvPlots( sce1[,"Beta"], base_objs$Zprop)

# if( !exists("irfm"))
# {
#     ca <- with( base_objs, MakeComps( Zprop, sce1[,"Beta"], FALSE))
# #     cb <- with( base_objs, MakeComps( Zprop, sce1[,"Beta"], TRUE))
#     irfm <- IRFMats( 10, base_objs$Zprop, ca)
# # }
# 
# par1 <- list( s = 10, q = "mean", ztprop1 = base_objs$Zprop, Gibbs = sce1)
# # if(!exists("if1")) 
# if1 <- do.call( IRF, c( par1, exog = FALSE))




# hh <- lapply( rownames(irfm)[[1]], function(x) IRF( 0.2, 0.5, "vspend_tourists", x, 
# irfm, sce1, FALSE))
# conv1 <- BetaConvPlots( 
# bm <- execF$Stability( zp1, rowMeans( do.call( cbind, sce1$Beta)))
# ba <- execF$MakeA( zp1, sce1$Beta[[100]])

# if(!exists( "debl"))
# {
#     debl <- list()
#     for( i in 1:3)
#     {
#         PATHS$forfile <- paste("eqs_orig/eqs_econ_str", i, sep = "")
#         debl[[i]] <- CreateData( PATHS, priors = TRUE)
#     }
# }
# 
# if( !exists( "GibbsL"))
# {
#     GibbsL <- lapply( 1:3, function(x) with(debl[[x]], Gibbs(Z,Y, ilist, 
#     10000)))
# }