0

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
}
Matheus Lacerda
  • 5,983
  • 11
  • 29
  • 45
Christian
  • 3
  • 4

1 Answers1

0

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)
FAMG
  • 395
  • 1
  • 9
  • OK, the key step appears to be saving the simulation before plotting. Didn't realize I had to do this. Thank you for your help. – Christian Oct 28 '17 at 21:39
  • You are welcome. Would be nice if you could mark my answer as solving your problem. Thx! – FAMG Oct 29 '17 at 08:50