Posted By

stevejb on 03/24/11


Tagged

simulation


Versions (?)

simulateMoments.R


 / Published in: R
 

  1. #
  2. #
  3. # Stephen J. Barr
  4. #
  5. # get moments from the simulated data
  6. #
  7. #
  8.  
  9. #
  10. # derivative
  11. #
  12.  
  13. require(moments)
  14.  
  15. getMomentsFromSim <- function(mySim, moment.par) {
  16.  
  17. # there are 50 time periods. Drop first 10
  18. T <- 50
  19. Tsub <- T - 10;
  20. N <- moment.par$N
  21. STATEMAT <- moment.par$STATEMAT
  22. ValFn <- moment.par$VALFN
  23. sumStat <- matrix(nrow=(Tsub*N), ncol=2);
  24. KLDmat <- matrix(nrow=(Tsub*N), ncol=3)
  25. colnames(sumStat) <- c("investment", "tobinsq")
  26. colnames(KLDmat) <- c("K", "L", "D")
  27. KM <- 1
  28. LM <- 2
  29. DM <- 3
  30. delta <- .15
  31.  
  32. for(i in 1:N) {
  33. # bellmanSubset <- bellmanEQ[(length(S)*(i-1) + 1):(i*length(S)),]
  34. startIdx <- (Tsub*(i-1) + 1)
  35. endIdx <- Tsub * i
  36. POLICY <- mySim[[i]]$KLDVAR
  37. ZTMP <- mySim[[i]]$Z
  38.  
  39. kprimeTMP <- STATEMAT[POLICY[(T-Tsub+2):T], KM]
  40. kTMP <- STATEMAT[POLICY[(T-Tsub+1):(T-1)], KM]
  41. investTMP <- (kprimeTMP - (1-delta)*kTMP)/kTMP;
  42. sumStat[(startIdx+1):endIdx, 1] <- investTMP;
  43.  
  44. kt1 <- STATEMAT[POLICY[(T-Tsub):(T-1)], KM]
  45. lt1 <- STATEMAT[POLICY[(T-Tsub):(T-1)], LM]
  46. dt1 <- STATEMAT[POLICY[(T-Tsub):(T-1)], DM]
  47.  
  48. KLDmat[startIdx:endIdx,KM] <- kt1
  49. KLDmat[startIdx:endIdx,LM] <- lt1
  50. KLDmat[startIdx:endIdx,DM] <- dt1
  51.  
  52. # calculate Tobin's Q
  53. count <- 0
  54. for(j in (T-Tsub+1):T) {
  55. # val function at j
  56. vfj <- ValFn[ZTMP[j], POLICY[j]]
  57. kj <- STATEMAT[POLICY[j], KM]
  58. tqTmp <- vfj/kj
  59. sumStat[(startIdx + count), 2] <- tqTmp;
  60. count <- count+1
  61.  
  62. } # end for j
  63.  
  64. } # end for N
  65.  
  66. data.filtered <- as.data.frame(na.omit(sumStat))
  67. simRes <- colMeans(data.filtered)
  68. ## colMeans(na.omit(sumStat))
  69. ## investment tobinsq
  70. ## 0.2099869 4.7489438
  71. partEreg <- lm(investment ~ tobinsq, data=data.filtered)
  72.  
  73. sdvec <- apply(KLDmat, 2, function(x) {moment(x,order=2)})
  74. sdc <- apply(KLDmat, 2, function(x) {moment(x,order=2, central=TRUE)})
  75. thirdmoment <- apply(KLDmat, 2, function(x) {moment(x,order=3)})
  76. thirdmoment.d <- apply(KLDmat, 2, function(x) {moment(x,order=3, central=TRUE)})
  77.  
  78. momentVec <- c(simRes, coefficients(partEreg)[2], sdvec, sdc, thirdmoment, thirdmoment.d)
  79.  
  80.  
  81. return(momentVec)
  82. }
  83.  
  84.  
  85.  
  86.  
  87.  
  88.  
  89.  
  90.  
  91.  
  92. ###########################################################3
  93. #
  94. # MAIN:
  95. #
  96. #
  97.  
  98.  
  99. # I have to do bit of regeneration here, but should be alright
  100. #
  101. #
  102.  
  103. # the results are in NEILDIAMOND
  104. # the model list is in MODELLIST.Robj
  105. #
  106.  
  107. # all I need to do is the following
  108. load("../data/MODELLIST.Robj")
  109. load("../data/RESULTLIST.Robj")
  110.  
  111. # there should be 834 elements in each of these lists
  112. #
  113. #
  114.  
  115.  
  116.  
  117.  
  118. require(doMC)
  119. registerDoMC()
  120.  
  121. momentmat <- matrix(nrow=15, ncol=length(modellist))
  122.  
  123.  
  124. mm <- foreach(d = 1:length(modellist), .combine=cbind) %dopar% {
  125.  
  126. cat(paste("Starting model:", d, "\n"))
  127. MODEL.TMP <- modellist[[d]]
  128. RESULTS.TMP <- NEILDIAMOND[[d]]
  129.  
  130. # construct the paramters necessary to get the moments
  131. moment.par <- list(STATEMAT = MODEL.TMP[[1]]$STATEMATRIX,
  132. VALFN = MODEL.TMP[[2]]$VALUEFN,
  133. N=10000)
  134.  
  135. mymoments <- getMomentsFromSim(RESULTS.TMP, moment.par)
  136. momentmat[,d] <- momentVec
  137. cat(paste("FINISHED model:", d, "!!! I love Jo !!\n"))
  138. }
  139.  
  140. filename <- paste(format(Sys.time(), "%s"), "-momentsmat.Robj", sep="")
  141. save(momentmat, file=filename)

Report this snippet  

You need to login to post a comment.