# 
#
# 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 following
load("../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)

simulateMoments.R

R