Convergencia de una cadena de Markov finita ergódica

P : matriz de transición de un paso

suppressPackageStartupMessages(library(expm))

El paquete expm permite elevar una matriz a una potencia dada con la función %^%

P = matrix(c(0.4, 0.3, 0.1, 0.1, 0.3, 0.2, 0.4, 0.4, 0.2, 0.4, 0.2, 0.3, 0.1, 
    0.1, 0.3, 0.2), nrow = 4)

Estimando P^∞

P

##      [,1] [,2] [,3] [,4]
## [1,]  0.4  0.3  0.2  0.1
## [2,]  0.3  0.2  0.4  0.1
## [3,]  0.1  0.4  0.2  0.3
## [4,]  0.1  0.4  0.3  0.2

P %^% 10

##        [,1]  [,2]   [,3]   [,4]
## [1,] 0.2326 0.314 0.2801 0.1734
## [2,] 0.2326 0.314 0.2801 0.1734
## [3,] 0.2326 0.314 0.2801 0.1734
## [4,] 0.2326 0.314 0.2801 0.1734

Simulación de los estados del proceso.

nsim = 1e+05
x = numeric(nsim)

x[1] = 4  # inicio estado

for (i in 2:nsim) {
    if (x[i - 1] == 1) 
        x[i] = sample(1:4, 1, prob = P[1, ])
    if (x[i - 1] == 2) 
        x[i] = sample(1:4, 1, prob = P[2, ])
    if (x[i - 1] == 3) 
        x[i] = sample(1:4, 1, prob = P[3, ])
    if (x[i - 1] == 4) 
        x[i] = sample(1:4, 1, prob = P[4, ])
}

head(x, 20)

##  [1] 4 4 1 1 3 2 2 3 3 3 4 2 3 4 2 1 3 4 3 3

table(x)

## x
##     1     2     3     4 
## 23093 31213 28200 17494

table(x)/nsim

## x
##      1      2      3      4 
## 0.2309 0.3121 0.2820 0.1749