modelSim.R

#!/usr/bin/env Rscript
 
myargs <- commandArgs(TRUE )
 
START.I <- myargs[[1]]
END.I <- myargs[[2]]
 
DEBUG <- FALSE
 
# Stephen J. Barr
#
 
cat(paste("Running model from: ", START.I, "to:", END.I, "\n"))
MODLIST <- as.list(START.I:END.I)
 
 
 
########################################################3
if(DEBUG) {
  for(i in 1:length(MODLIST)) {
    cat(paste("\tMODEL", MODLIST[[i]], "\n"))
  }
}
########################################################3
 
modelSim <- function(seed, MODNUM) {
  MODFILE <- paste("model.", MODNUM, ".Robj", sep="")
  load(MODFILE) # this brings in the resSave object
 
  policyMat <- resSave[[2]]$POLICYFN
 
  # PARAMETERS
  P <- matrix(c(.6,.4,.4,.6), nrow=2, byrow=TRUE)
  T <- 50
  initialProbVec <- c(.5,1) # this is cumulative
 
  # GENERATE RANDOM u for SIM
  stateSize <- nrow(P);
  u <- runif(T);
 
  # get the probability matrix in the form of cumulative sums by row
  Pcum <- t(apply(P,1,cumsum))
 
  Z <- rep(0,T);
  KLD <- rep(0,T);
 
  # initital state probability
  ltVec <- which(u[1] < initialProbVec)
  Z[1] <- ltVec[1]
  KLD[1] <- ceiling(runif(1)*ncol(policyMat))
 
    # do the rest of the simulation
  for(t in 2:T) {
    prevState <- Z[t-1];
 
    # draw current Z
    ltVec <- which(u[t] < Pcum[prevState,])
    Z[t] <- ltVec[1]
 
    # using policy function and current shock, get new capital stock
    KLD[t] <- policyMat[Z[t],KLD[t-1]]
  } # end for t in 2:T
 
  return(list(KLDVAR=KLD, Z=Z))
} # end model sim
 
 
outerRun <- function(MODNUM, seedlist) {
  cat(paste("outerRun called on MN:", unlist(MODNUM), "seeds:", unlist(seedlist), "\n"))
  cat(length(seedlist))
  outerres <- lapply(seedlist, modelSim, MODNUM)
  return(outerres)
}
 
 
require(multicore)
 
 
#########################################################3
if(DEBUG) {
  seedlist <- c(5678910, 11121314)
  mlshort <- MODLIST[1:2]
  DBGRESULTS <- mclapply(mlshort, outerRun, seedlist)
 
} else {
#########################################################3
  seedlistfull <- as.list(30000:39999)
  RESULTS <- mclapply(MODLIST, outerRun, seedlistfull)
  rname <- paste("RESULTS-", START.I, "-", END.I, ".Robj", sep="")
  save(RESULTS, file=rname)
}