## Posted By

stevejb on 03/24/11

# simulateMoments.R

/ Published in: R

`# ## Stephen J. Barr## get moments from the simulated data## ## derivative # require(moments) getMomentsFromSim <- function(mySim, moment.par) {   # there are 50 time periods. Drop first 10  T <- 50  Tsub <- T - 10;  N <- moment.par\$N  STATEMAT <- moment.par\$STATEMAT  ValFn <- moment.par\$VALFN  sumStat <- matrix(nrow=(Tsub*N), ncol=2);  KLDmat <- matrix(nrow=(Tsub*N), ncol=3)  colnames(sumStat) <- c("investment", "tobinsq")  colnames(KLDmat) <- c("K", "L", "D")  KM <- 1  LM <- 2  DM <- 3  delta <- .15   for(i in 1:N) {    #  bellmanSubset <- bellmanEQ[(length(S)*(i-1) + 1):(i*length(S)),]    startIdx <- (Tsub*(i-1) + 1)    endIdx <- Tsub * i    POLICY <- mySim[[i]]\$KLDVAR    ZTMP <- mySim[[i]]\$Z     kprimeTMP <- STATEMAT[POLICY[(T-Tsub+2):T], KM]    kTMP <- STATEMAT[POLICY[(T-Tsub+1):(T-1)], KM]    investTMP <- (kprimeTMP - (1-delta)*kTMP)/kTMP;    sumStat[(startIdx+1):endIdx, 1] <- investTMP;     kt1 <- STATEMAT[POLICY[(T-Tsub):(T-1)], KM]    lt1 <- STATEMAT[POLICY[(T-Tsub):(T-1)], LM]    dt1 <- STATEMAT[POLICY[(T-Tsub):(T-1)], DM]     KLDmat[startIdx:endIdx,KM] <- kt1    KLDmat[startIdx:endIdx,LM] <- lt1    KLDmat[startIdx:endIdx,DM] <- dt1     # calculate Tobin's Q    count <- 0    for(j in (T-Tsub+1):T) {      # val function at j      vfj <- ValFn[ZTMP[j], POLICY[j]]      kj <- STATEMAT[POLICY[j], KM]      tqTmp <- vfj/kj      sumStat[(startIdx + count), 2] <- tqTmp;      count <- count+1     } # end for j   } # end for N   data.filtered <- as.data.frame(na.omit(sumStat))  simRes <- colMeans(data.filtered)  ## colMeans(na.omit(sumStat))  ## investment    tobinsq   ##  0.2099869  4.7489438   partEreg <- lm(investment ~ tobinsq, data=data.filtered)   sdvec <- apply(KLDmat, 2, function(x) {moment(x,order=2)})  sdc <-   apply(KLDmat, 2, function(x) {moment(x,order=2, central=TRUE)})  thirdmoment <- apply(KLDmat, 2, function(x) {moment(x,order=3)})  thirdmoment.d <-   apply(KLDmat, 2, function(x) {moment(x,order=3, central=TRUE)})   momentVec <- c(simRes, coefficients(partEreg)[2], sdvec, sdc, thirdmoment, thirdmoment.d)    return(momentVec)}         ###########################################################3## MAIN:##  # I have to do bit of regeneration here, but should be alright## # the results are in NEILDIAMOND# the model list is in MODELLIST.Robj# # all I need to do is the followingload("../data/MODELLIST.Robj")load("../data/RESULTLIST.Robj") # there should be 834 elements in each of these lists##    require(doMC)registerDoMC() momentmat <- matrix(nrow=15, ncol=length(modellist))  mm <- foreach(d = 1:length(modellist), .combine=cbind) %dopar% {   cat(paste("Starting model:", d, "\n"))  MODEL.TMP <- modellist[[d]]  RESULTS.TMP <- NEILDIAMOND[[d]]   # construct the paramters necessary to get the moments  moment.par <- list(STATEMAT = MODEL.TMP[[1]]\$STATEMATRIX,                     VALFN = MODEL.TMP[[2]]\$VALUEFN,                     N=10000)   mymoments <- getMomentsFromSim(RESULTS.TMP, moment.par)  momentmat[,d] <- momentVec  cat(paste("FINISHED model:", d, "!!! I love Jo !!\n"))} filename <- paste(format(Sys.time(), "%s"), "-momentsmat.Robj", sep="")save(momentmat, file=filename)`