```

References

Load packages

library(assertthat)
library(tidyverse)

Simulation

## Simulation function
Simulate <- function(init, P, n_steps) {

    assert_that(length(init) == nrow(P))
    assert_that(length(init) == ncol(P))
    assert_that(n_steps >= 1)

    out <- matrix(as.numeric(NA),
                  nrow = n_steps,
                  ncol = length(init))

    ## Initial state
    out[1,] <- init

    for (i in seq_len(n_steps)[-1]) {
        out[i,] <- out[i-1,] %*% P
    }

    out <- as.data.frame(out)
    colnames(out) <- paste0("state", seq_len(ncol(out)))
    out$i <- seq_len(n_steps)
    out
}

## Some transition matrix
P <- matrix(c(0.1, 0.7, 0.1,0.1,
              0.3, 0.5, 0.2,0,
              0.4, 0.3, 0.3,0,
              0.1, 0.1, 0.2, 0.6), byrow = TRUE, nrow = 4)
P
##      [,1] [,2] [,3] [,4]
## [1,]  0.1  0.7  0.1  0.1
## [2,]  0.3  0.5  0.2  0.0
## [3,]  0.4  0.3  0.3  0.0
## [4,]  0.1  0.1  0.2  0.6
## Rows must sum to 1
rowSums(P)
## [1] 1 1 1 1
## Generate data from 3 initial states
df1 <- Simulate(init = c(1,0,0,0),
                P = P,
                n_steps = 10)
df1$init <- "c(1,0,0,0)"
df2 <- Simulate(init = c(0,1,0,0),
                P = P,
                n_steps = 10)
df2$init <- "c(0,1,0,0)"
df3 <- Simulate(init = c(0,0,1,0),
                P = P,
                n_steps = 10)
df3$init <- "c(0,0,1,0)"
df4 <- Simulate(init = c(0,0,0,1),
                P = P,
                n_steps = 10)
df4$init <- "c(0,0,0,1)"

## Clean data for plotting
df_long <- gather(data = rbind(df1, df2, df3, df4),
                  key = state, value = prob,
                  -init, -i)

## Plot
ggplot(data = df_long,
       mapping = aes(x = i, y = prob, color = state, group = state)) +
    geom_point() +
    geom_line() +
    facet_grid(init ~ .) +
    theme_bw() + theme(legend.key = element_blank())

Regardless of the initial state, they all converge to the same equilibrium.

## Check equilibrium
df1
##       state1    state2    state3     state4  i       init
## 1  1.0000000 0.0000000 0.0000000 0.00000000  1 c(1,0,0,0)
## 2  0.1000000 0.7000000 0.1000000 0.10000000  2 c(1,0,0,0)
## 3  0.2700000 0.4600000 0.2000000 0.07000000  3 c(1,0,0,0)
## 4  0.2520000 0.4860000 0.1930000 0.06900000  4 c(1,0,0,0)
## 5  0.2551000 0.4842000 0.1941000 0.06660000  5 c(1,0,0,0)
## 6  0.2550700 0.4855600 0.1939000 0.06547000  6 c(1,0,0,0)
## 7  0.2552820 0.4860460 0.1938830 0.06478900  7 c(1,0,0,0)
## 8  0.2553741 0.4863642 0.1938601 0.06440160  8 c(1,0,0,0)
## 9  0.2554309 0.4865422 0.1938486 0.06417837  9 c(1,0,0,0)
## 10 0.2554630 0.4866451 0.1938418 0.06405011 10 c(1,0,0,0)
## Examine eigen vectors
(res_eigen <- eigen(t(P)))
## $values
## [1]  1.00000000  0.57487218 -0.12886655  0.05399437
## 
## $vectors
##           [,1]        [,2]       [,3]        [,4]
## [1,] 0.4357025 -0.19223810  0.5873163  0.26080813
## [2,] 0.8300883 -0.61328480 -0.7633564 -0.77994640
## [3,] 0.3305329  0.04048207  0.2566195  0.56690484
## [4,] 0.1089256  0.76504083 -0.0805794 -0.04776656
## Eigenvector corresponding to eigenvalue = 1
(eigenvector <- res_eigen$vectors[,1])
## [1] 0.4357025 0.8300883 0.3305329 0.1089256
## Scale so that it sums to 1
(eigenvector <- eigenvector / sum(eigenvector))
## [1] 0.25550661 0.48678414 0.19383260 0.06387665

The stationary distribution \(\Pi\) (row vector) meets the following equality involving the transition matrix \(\mathbf{P}\).

\(\Pi ~ \mathbf{P} = \Pi\)

That is,

\(\mathbf{P}^T ~ \Pi^T = (1)\Pi^T\)

This means the stationary vector is an eigenvector of matrix \(\mathbf{P}^T\), scaled so that sum of elements is 1.

Simulation (pathological example with block diagonal P)

## Some transition matrix
P <- matrix(c(0.3, 0.7, 0.0, 0.0,
              0.1, 0.9, 0.0, 0.0,
              0.0, 0.0, 0.3, 0.7,
              0.0, 0.0, 0.4, 0.6), byrow = TRUE, nrow = 4)
P
##      [,1] [,2] [,3] [,4]
## [1,]  0.3  0.7  0.0  0.0
## [2,]  0.1  0.9  0.0  0.0
## [3,]  0.0  0.0  0.3  0.7
## [4,]  0.0  0.0  0.4  0.6
## Rows must sum to 1
rowSums(P)
## [1] 1 1 1 1
## Generate data from 3 initial states
df1 <- Simulate(init = c(1,0,0,0),
                P = P,
                n_steps = 10)
df1$init <- "c(1,0,0,0)"
df2 <- Simulate(init = c(0,1,0,0),
                P = P,
                n_steps = 10)
df2$init <- "c(0,1,0,0)"
df3 <- Simulate(init = c(0,0,1,0),
                P = P,
                n_steps = 10)
df3$init <- "c(0,0,1,0)"
df4 <- Simulate(init = c(0,0,0,1),
                P = P,
                n_steps = 10)
df4$init <- "c(0,0,0,1)"

## Clean data for plotting
df_long <- gather(data = rbind(df1, df2, df3, df4),
                  key = state, value = prob,
                  -init, -i)

## Plot
ggplot(data = df_long,
       mapping = aes(x = i, y = prob, color = state, group = state)) +
    geom_point() +
    geom_line() +
    facet_grid(init ~ .) +
    theme_bw() + theme(legend.key = element_blank())

Probability cannot travel across the blocks. Thus, not all initial states result in the same equilibrium.

eigen(t(P))
## $values
## [1]  1.0  1.0  0.2 -0.1
## 
## $vectors
##            [,1]       [,2]       [,3]       [,4]
## [1,] -0.1414214  0.0000000 -0.7071068  0.0000000
## [2,] -0.9899495  0.0000000  0.7071068  0.0000000
## [3,]  0.0000000 -0.4961389  0.0000000 -0.7071068
## [4,]  0.0000000 -0.8682431  0.0000000  0.7071068

In this case, the \(\mathbf{P}^T\) has multiple eigenvalues that are 1, resulting in multiple eigenvectors corresponding to them.