Steady States in Markov Chains

Below, we analyze the long term behavior of system under three policies using a stochastic model to make better decisions on what policy to employ.

First let us begin by loading relevant libraries:

library(expm)

Warning messages:
1: In formals(fun) : argument is not a function
2: In formals(fun) : argument is not a function
3: In formals(fun) : argument is not a function

library(reshape2)

In this example, we consider a situation where we have a viral infection spreading through a community. We have three possible policies:

• Policy 0: Baseline case, we do nothing
• Policy I: Education + early screening
• Policy II: Late stage treatment focus

We have three different transition matrices corresponding to the probability of transitioning between states in the three different situations.

Warning message:
In formals(fun) : argument is not a function

P.control <- matrix(c(0.85, 0.1,    0.05,   0,  0,
0.2,    0.5,    0.2,    0.1,    0,
0,  0.1,    0.6,    0.2,    0.1,
0,  0.05,   0.15,   0.7,    0.1,
0,  0,  0.1,    0.3,    0.6), nrow=5, byrow=TRUE)


P.policy.i <- matrix(c(0.9, 0.08,   0.02,   0,  0,
0.15,   0.6,    0.15,   0.1,    0,
0,  0.05,   0.65,   0.25,   0.05,
0,  0.02,   0.1,    0.8,    0.08,
0,  0,  0.05,   0.2,    0.75), nrow=5, byrow=TRUE)

P.policy.ii <- matrix(c(0.85,   0.1,    0.05,   0,  0,
0.2,    0.5,    0.2,    0.1,    0,
0,  0.05,   0.6,    0.25,   0.1,
0,  0.02,   0.1,    0.75,   0.13,
0,  0,  0.1,    0.25,   0.65), nrow=5, byrow=TRUE)

We can track the evolution of the system under each policy by multiplying the transition matrix raised to the power of the iteration by the initial distribution (\(\pi_0\)):

\(\pi_i = \pi_0 P^i\)

The steady state will be reached at \(\pi_\infty\) when \(\pi_\infty P = \pi_\infty\). In this example, we simply calculate the distribution of the system after sufficiently enough time has passed.

To find the steady state probability we have to find:

\(\pi_\infty(P-I)=0\)

Below we define the initial population:

Initial.pop <- c(300000, 100000, 50000, 30000, 20000)

Define parameters for solution:

num.itr <- 100

Baseline

control <- matrix(nrow=num.itr, ncol=5)

for(i in 1:num.itr){
  control[i,] <- Initial.pop %*% (P.control %^% i)
}


control_df <- as.data.frame(control)
colnames(control_df) <- c('Healthy', 'At-risk', 'Infected', 'Recovered', 'Chronic Condition')
control_df$Iteration <- 1:num.itr


control_long <- melt(control_df, id.vars = "Iteration", variable.name = "State", value.name = "Population")

ggplot(control_long, aes(x = Iteration, y = Population, color = State, group = State)) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  theme_minimal() +
  labs(title = "Distribution of population - baseline",
       x = "Iteration",
       y = "Population",
       color = "State") +
  theme(text = element_text(size=14))+
  scale_color_manual(values = c('Healthy'='#27445D', 'At-risk'='#780C28', 'Infected'='#71BBB2', 'Recovered'='#3F4F44', 'Chronic Condition'='#ADB2D4'))

NA
NA

n <- nrow(P.control)
A <- t(P.control - diag(n)) 
A <- rbind(A, rep(1, n))  
b <- c(rep(0, n), 1)  


steady_state <- ginv(A) %*% b 
steady_state <- steady_state / sum(steady_state)  


print(steady_state*500000)

          [,1]
[1,]  75520.83
[2,]  56640.62
[3,] 121093.75
[4,] 173177.08
[5,]  73567.71

Great, now let us plot the average cost for the population over time.

ggplot()+geom_line(aes(x = 1:length(control%*%c(0, 200, 8000, 500, 25000)), y = control%*%c(0, 200, 8000, 500, 25000)), color = '#ADB2D4', size=1.5) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  theme_minimal() +
  labs(title = "Cost over time",
       x = "Iteration",
       y = "Cost") +
  theme(text = element_text(size=14))

Okay, let’s see how our system behaves under the two other policies.

Policy II

policy.i <- matrix(nrow=num.itr, ncol=5)

for(i in 1:num.itr){
  policy.i[i,] <- Initial.pop %*% (P.policy.i %^% i)
}


policy.i_df <- as.data.frame(policy.i)
colnames(policy.i_df) <- c('Healthy', 'At-risk', 'Infected', 'Recovered', 'Chronic Condition')
policy.i_df$Iteration <- 1:num.itr


policy.i_long <- melt(policy.i_df, id.vars = "Iteration", variable.name = "State", value.name = "Population")

ggplot(policy.i_long, aes(x = Iteration, y = Population, color = State, group = State)) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  theme_minimal() +
  labs(title = "Distribution of population - policy i",
       x = "Iteration",
       y = "Population",
       color = "State") +
  theme(text = element_text(size=14))+
  scale_color_manual(values = c('Healthy'='#27445D', 'At-risk'='#780C28', 'Infected'='#71BBB2', 'Recovered'='#3F4F44', 'Chronic Condition'='#ADB2D4'))

NA
NA

Let’s evaluate steady states

n <- nrow(P.policy.i)
A <- t(P.policy.i - diag(n)) 
A <- rbind(A, rep(1, n))  
b <- c(rep(0, n), 1)  


steady_state <- ginv(A) %*% b 
steady_state <- steady_state / sum(steady_state)  


print(steady_state*500000)

          [,1]
[1,]  50111.36
[2,]  33407.57
[3,]  95662.32
[4,] 228550.22
[5,]  92268.53

Average cost over time for policy i:

ggplot()+geom_line(aes(x = 1:length(policy.i%*%c(0, 200, 8000, 500, 25000)), y = policy.i%*%c(0, 200, 8000, 500, 25000)), color = '#ADB2D4', size=1.5) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  theme_minimal() +
  labs(title = "Cost over time",
       x = "Iteration",
       y = "Cost") +
  theme(text = element_text(size=14))

Policy ii

policy.ii <- matrix(nrow=num.itr, ncol=5)

for(i in 1:num.itr){
  policy.ii[i,] <- Initial.pop %*% (P.policy.ii %^% i)
}


policy.ii_df <- as.data.frame(policy.ii)
colnames(policy.ii_df) <- c('Healthy', 'At-risk', 'Infected', 'Recovered', 'Chronic Condition')
policy.ii_df$Iteration <- 1:num.itr


policy.ii_long <- melt(policy.ii_df, id.vars = "Iteration", variable.name = "State", value.name = "Population")

ggplot(policy.ii_long, aes(x = Iteration, y = Population, color = State, group = State)) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  theme_minimal() +
  labs(title = "Distribution of population - policy ii",
       x = "Iteration",
       y = "Population",
       color = "State") +
  theme(text = element_text(size=14))+
  scale_color_manual(values = c('Healthy'='#27445D', 'At-risk'='#780C28', 'Infected'='#71BBB2', 'Recovered'='#3F4F44', 'Chronic Condition'='#ADB2D4'))

NA
NA

Steady states:

n <- nrow(P.policy.ii)
A <- t(P.policy.ii - diag(n)) 
A <- rbind(A, rep(1, n))  
b <- c(rep(0, n), 1)  


steady_state <- ginv(A) %*% b 
steady_state <- steady_state / sum(steady_state)  


print(steady_state*500000)

          [,1]
[1,]  34843.21
[2,]  26132.40
[3,] 101742.16
[4,] 224738.68
[5,] 112543.55

Average cost over time for policy ii:

ggplot()+geom_line(aes(x = 1:length(policy.i%*%c(0, 200, 8000, 500, 25000)), y = policy.i%*%c(0, 200, 8000, 500, 25000)), color = '#ADB2D4', size=1.5) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  theme_minimal() +
  labs(title = "Cost over time",
       x = "Iteration",
       y = "Cost") +
  theme(text = element_text(size=14))

As demonstrated all three Markov chains are ergodic and reach a steady state distribution after some time.

Comparison of costs under different policies

df <- data.frame(
  Iteration = 1:length(control %*% c(0, 200, 8000, 500, 25000)),
  Control = as.vector(control %*% c(0, 200, 8000, 500, 25000)),
  Policy_I = as.vector(policy.i %*% c(0, 200, 8000, 500, 25000)),
  Policy_II = as.vector(policy.ii %*% c(0, 200, 8000, 500, 25000))
)

# Convert data to long format for ggplot
library(reshape2)
df_long <- melt(df, id.vars = "Iteration", variable.name = "Policy", value.name = "Cost")

# Plot with legend
ggplot(df_long, aes(x = Iteration, y = Cost, color = Policy, group = Policy)) +
  geom_line(size=1.5) +
  geom_point(size = 2) +
  theme_minimal() +
  labs(title = "Cost over Time",
       x = "Iteration",
       y = "Cost",
       color = "Policy") +  # Legend label
  scale_color_manual(values = c("Control" = "#ADB2D4", "Policy_I" = "darkred", "Policy_II" = "#27445D")) +
  theme(text = element_text(size=14))

It seems a little bit counter intuitive and strange to see that long term it appears that the control case with no intervention is less expensive compared to the other two policies long-term. However, upon further inspection we realize that having patients with chronic condition is the most costly. Interestingly, under the two policies the probability of staying in the chronic condition state starting from the chronic condition state is higher than the control case. The chronic condition cases appear to be much lower under the control condition compared to the other two policies. However, the number of infected individuals remain much higher in the control case compared to the two policies. There is a trade-off to this and that is the fact that the number of at-risk people is significantly higher in the control case compared to the other conditions involving intervention.

Mean First Passage Time (MFPT)

find.MFPT <- function(P) {
  n <- nrow(P)
  MFPT <- matrix(0, n, n)  
  
  for (i in 1:n) {
    for (j in 1:n) {
      if (i != j) {
        
        A <- diag(n) - P
        A[j, ] <- 0  
        A[j, j] <- 1  
        b <- rep(1, n)
        b[j] <- 0 
        
    
        m <- solve(A, b)
        MFPT[i, j] <- m[i]  
      }
    }
  }
  
  return(MFPT)
}

Let’s find the overall accessibility of each state using the MFPT matrix and compare the different conditions:

print(paste0("Control: ", colSums(find.MFPT(P.control))))

[1] "Control: 184.137931034483" "Control: 69.6551724137931" "Control: 36.1290322580645"
[4] "Control: 37.2180451127819" "Control: 76.1946902654867"

print(paste0("Policy I: ", colSums(find.MFPT(P.policy.i))))

[1] "Policy I: 464.444444444447" "Policy I: 158.666666666667"
[3] "Policy I: 58.4811529933481" "Policy I: 43.7470997679815"
[5] "Policy I: 99.6206896551725"

print(paste0("Policy II: ", colSums(find.MFPT(P.policy.ii))))

[1] "Policy II: 439.000000000002" "Policy II: 164.666666666667"
[3] "Policy II: 41.5582191780822" "Policy II: 33.6821705426357"
[5] "Policy II: 58.0185758513932"

It seems like based on this analysis, the healthy state is much easier to reach in the control case than the two other conditions.