Posted By

stevejb on 03/24/11


Tagged

R simulation dynamicprogramming


Versions (?)

modelSim.R


 / Published in: R
 

  1. #!/usr/bin/env Rscript
  2.  
  3. myargs <- commandArgs(TRUE )
  4.  
  5. START.I <- myargs[[1]]
  6. END.I <- myargs[[2]]
  7.  
  8. DEBUG <- FALSE
  9.  
  10. # Stephen J. Barr
  11. #
  12.  
  13. cat(paste("Running model from: ", START.I, "to:", END.I, "\n"))
  14. MODLIST <- as.list(START.I:END.I)
  15.  
  16.  
  17.  
  18. ########################################################3
  19. if(DEBUG) {
  20. for(i in 1:length(MODLIST)) {
  21. cat(paste("\tMODEL", MODLIST[[i]], "\n"))
  22. }
  23. }
  24. ########################################################3
  25.  
  26. modelSim <- function(seed, MODNUM) {
  27. MODFILE <- paste("model.", MODNUM, ".Robj", sep="")
  28. load(MODFILE) # this brings in the resSave object
  29.  
  30. policyMat <- resSave[[2]]$POLICYFN
  31.  
  32. # PARAMETERS
  33. P <- matrix(c(.6,.4,.4,.6), nrow=2, byrow=TRUE)
  34. T <- 50
  35. initialProbVec <- c(.5,1) # this is cumulative
  36.  
  37. # GENERATE RANDOM u for SIM
  38. stateSize <- nrow(P);
  39. u <- runif(T);
  40.  
  41. # get the probability matrix in the form of cumulative sums by row
  42. Pcum <- t(apply(P,1,cumsum))
  43.  
  44. Z <- rep(0,T);
  45. KLD <- rep(0,T);
  46.  
  47. # initital state probability
  48. ltVec <- which(u[1] < initialProbVec)
  49. Z[1] <- ltVec[1]
  50. KLD[1] <- ceiling(runif(1)*ncol(policyMat))
  51.  
  52. # do the rest of the simulation
  53. for(t in 2:T) {
  54. prevState <- Z[t-1];
  55.  
  56. # draw current Z
  57. ltVec <- which(u[t] < Pcum[prevState,])
  58. Z[t] <- ltVec[1]
  59.  
  60. # using policy function and current shock, get new capital stock
  61. KLD[t] <- policyMat[Z[t],KLD[t-1]]
  62. } # end for t in 2:T
  63.  
  64. return(list(KLDVAR=KLD, Z=Z))
  65. } # end model sim
  66.  
  67.  
  68. outerRun <- function(MODNUM, seedlist) {
  69. cat(paste("outerRun called on MN:", unlist(MODNUM), "seeds:", unlist(seedlist), "\n"))
  70. cat(length(seedlist))
  71. outerres <- lapply(seedlist, modelSim, MODNUM)
  72. return(outerres)
  73. }
  74.  
  75.  
  76. require(multicore)
  77.  
  78.  
  79. #########################################################3
  80. if(DEBUG) {
  81. seedlist <- c(5678910, 11121314)
  82. mlshort <- MODLIST[1:2]
  83. DBGRESULTS <- mclapply(mlshort, outerRun, seedlist)
  84.  
  85. } else {
  86. #########################################################3
  87. seedlistfull <- as.list(30000:39999)
  88. RESULTS <- mclapply(MODLIST, outerRun, seedlistfull)
  89. rname <- paste("RESULTS-", START.I, "-", END.I, ".Robj", sep="")
  90. save(RESULTS, file=rname)
  91. }

Report this snippet  

You need to login to post a comment.