Advanced Search

Posted By

stevejb on 03/24/11

Tagged

Versions (?)

Last Edited at 03/24/11 09:43am

Statistics

Viewed 159 times
Favorited by 0 user(s)

simulateMoments.R

/ Published in: R
Save to your folder(s)

Expand | Embed | Plain Text

Copy this code and paste it in your HTML 
  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
  62.     } # end for j
  64.   } # end for N
  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)
  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)})
  78.   momentVec <- c(simRes, coefficients(partEreg)[2], sdvec, sdc, thirdmoment, thirdmoment.d)
  80.   
  81.   return(momentVec)
  82. }
  92. ###########################################################3
  93. #
  94. # MAIN:
  95. #
  96. #
  99. # I have to do bit of regeneration here, but should be alright
  100. #
  101. #
  103. # the results are in NEILDIAMOND
  104. # the model list is in MODELLIST.Robj
  105. #
  107. # all I need to do is the following
  108. load("../data/MODELLIST.Robj")
  109. load("../data/RESULTLIST.Robj")
  111. # there should be 834 elements in each of these lists
  112. #
  113. #
  118. require(doMC)
  119. registerDoMC()
  121. momentmat <- matrix(nrow=15, ncol=length(modellist))
  124. mm <- foreach(d = 1:length(modellist), .combine=cbind) %dopar% {
  126.   cat(paste("Starting model:", d, "\n"))
  127.   MODEL.TMP <- modellist[[d]]
  128.   RESULTS.TMP <- NEILDIAMOND[[d]]
  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)
  135.   mymoments <- getMomentsFromSim(RESULTS.TMP, moment.par)
  136.   momentmat[,d] <- momentVec
  137.   cat(paste("FINISHED model:", d, "!!! I love Jo !!\n"))
  138. }
  140. filename <- paste(format(Sys.time(), "%s"), "-momentsmat.Robj", sep="")
  141. save(momentmat, file=filename)

Report this snippet


Comments

RSS Icon Subscribe to comments

You need to login to post a comment.