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
}
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)