# esempio catena con tre stati (Andrieu at al.) 
T <- rbind(c(0,1,0), c(0., 0.1, 0.9), c(0.6, 0.4, 0))
mu <- c(0.5, 0.2, 0.3)

mu %*% T
mu %*% T %*% T


matrix.pow <- function(a, n) {
  if (n == 1) return(a)
  for (i in 2:n) a <- a %*% a
  return(a)
}

matrix.pow(T,1)
matrix.pow(T,10)
c(1, 0, 0) %*% matrix.pow(T,10)
c(0, 0, 1) %*% matrix.pow(T,10)
c(0.33, 0.33, 0.34) %*% matrix.pow(T,10)

# attenzione: ci sono problemi di instabilita' per n > 50

# autovettori
library(panel)
# per istallare il pacchetto eseguire il comando
# install.packages("panel")
ris <- eddcmp(t(T)) 
p <- t(ris$evectors[,3]/sum(ris$evectors[,3]))
#check 
p %*% T