Search code examples
rbayesianmarkov-chainseconomics

Histogram of MC simulation


I am trying to compare histograms of a Markov chain (MC) simulation and actual data. I have tried to run the simulation using the code below, which I don't fully understand. R seems to have accepted the code, but I don't know how to run the histograms... For background, the data are Expansions and Contractions of the US economy (found here: http://www.nber.org/cycles.html). I have set up the transition matrix between these two states as P where columns sum to 1, and the transitions between states are calculated as "transitions / months in each state". I think n corresponds to transitions here, but I could be wrong...

P <- matrix(c(0.74961, 0.57291, 0.25039, 0.42709),2,2)
P <- t(P)
colSums(P)
n <- 33

MC.sim <- function(n,P) {
    sim<-c()
    m <- ncol(P)    
    sim[1] <- sample(1:m,1) 
    for(i in 2:n){
        newstate <- sample(1:m,1,prob=P[,sim[i-1]])
        sim[i] <- newstate
    }
    sim
}

Solution

  • As stated in the comments you could use hist():

    P <- matrix(c(0.74961, 0.57291, 0.25039, 0.42709),2,2)
    P <- t(P)
    colSums(P)
    n <- 33
    
    MC.sim <- function(n,P) {
      sim<-c()
      m <- ncol(P)    
      sim[1] <- sample(1:m,1) 
      for(i in 2:n){
        newstate <- sample(1:m,1,prob=P[,sim[i-1]])
        sim[i] <- newstate
      }
      return(sim)
    }
    
    # Save results of simulation
    temp <- MC.sim(n,P)
    
    #plot basic histogram
    hist(temp)
    

    I did not check your simulation. But it only generates values of 1 or 2.

    You could also use ggplot():

    # Save results of simulation
    temp <- MC.sim(n,P)
    
    # Make data frame for ggplot
    t <- as.data.frame(temp)
    names(t) <- "Sim"
    
    p <- ggplot(t, aes(x = Sim))
    p + geom_histogram(binwidth = 0.5)