Manually Performing Calculations in a Continuous Time Markov Chain

Question

As a learning exercise, I am trying to learn how to fit and simulate from Continuous Time Markov Chains.

Suppose I have a Stochastic Process that can assume 3 States: S1, S2 or S3. Lets say I have observe this process making 100 state transitions and I record the time that each transition was made at. I simulated this data using the R programming language:

# Define the possible states
states <- c("S1", "S2", "S3")

# Initialize the sequence of state transitions with the first state
transitions <- "S1"

# Loop to generate 99 more state transitions
for (i in 2:100) {
  # Get the current state
  current_state <- transitions[i-1]
 
  # Generate a new state that is different from the current state
  new_state <- current_state
  while (new_state == current_state) {
    new_state <- sample(states, 1)
  }
 
  # Add the new state to the sequence of state transitions
  transitions <- c(transitions, new_state)
}

# Set the rate parameter of the exponential distribution
lambda <- 0.5

# Simulate 100 observations from the exponential distribution
times <- rexp(100, rate = lambda)

# final data: note that special emphasis was placed to ensure that no two consecutive state transitions occur as this can not happen within a Continuous Time Markov Chain
my_data = data.frame(id = 1:100, transitions, times)

As I understand, there are 4 Steps required in fitting a Continuous Time Markov Chain:

• Step 1: First we need to calculate the probabilities within the Embedded Discrete Time Markov Chain. This contains the probabilities of transitioning from "state i to state j" and can be calculated as the number of times the process went from 'state i to state j' divided by the number of total exits out of 'state i', represented as $$p_{ij} = \frac{n_{ij}}{n_i}$$

• Step 2: Next, we need to calculate the "lambda_i" - this is the "rate of leaving state i". The formula for "lambda_i" is based on the Maximum Likelihood Estimate of the "Rate Parameter" in an Exponential Distribution. This is equal to the number of transitions out of "state i" divided by the total time spent in "state i" prior to exiting: $$\lambda_i = \frac{n_i}{t_i}$$

• Step 3: Next, we need to calculate the entries of the "Rate Matrix" (Q_ij). These entries represent the rates of transitioning between two states "i" and "j". This can be considered as the "rate of leaving "state i" multiplied by the probability of transitioning from "state i to state j". We can write this as follows:

$$Q_{ij} = \begin{cases} \lambda_i p_{ij}, & \text{if } i \neq j \ \sum_{k\neq i} Q_{ik}, & \text{if } i = j \end{cases}$$

• Step 4: Finally, we can calculate the "time dependent probabilities of transitioning from state i to state j" :

$$ p_{ij}(t) = e^{Qt} $$

The above expression is said to be approximated using Numerical Methods such as Taylor Expansions or Pade Approximants.

For the above data - I tried to write a computer program to first fit a Continuous Time Markov Chain to this data, and then to simulate trajectories of this Continuous Time Markov Chain at different time periods:

# Step 1

# Define the possible states
states <- c("S1", "S2", "S3")

# Initialize the count matrix
count_matrix <- matrix(0, nrow = length(states), ncol = length(states), dimnames = list(states, states))

# Loop through the state transitions and update the count matrix
for (i in 2:100) {
  # Get the current and previous states
  current_state <- transitions[i]
  prev_state <- transitions[i-1]

  # Increment the count for the transition from prev_state to current_state
  count_matrix[prev_state, current_state] <- count_matrix[prev_state, current_state] + 1
}

# Calculate the total number of exits from each state
exit_counts <- colSums(count_matrix)

# Calculate the transition probabilities
transition_probs <- count_matrix / exit_counts


# STEP 2

# Calculate lambda_1
lambda_1 <- sum(my_data$times[my_data$transitions == "S1"]) / nrow(my_data[my_data$transitions == "S1", ])

# Calculate lambda_2
lambda_2 <- sum(my_data$times[my_data$transitions == "S2"]) / nrow(my_data[my_data$transitions == "S2", ])

# Calculate lambda_3
lambda_3 <- sum(my_data$times[my_data$transitions == "S3"]) / nrow(my_data[my_data$transitions == "S3", ])


# STEP 3

# Define the rate matrix
rate_matrix <- matrix(0, nrow = length(states), ncol = length(states), dimnames = list(states, states))

# Define the lambda vector
lambda <- c(lambda_1, lambda_2, lambda_3)

# Calculate the rate matrix
for (i in 1:length(states)) {
    for (j in 1:length(states)) {
        if (i != j) {
            rate_matrix[i, j] <- lambda[i] * transition_probs[i, j]
        }
    }
    rate_matrix[i, i] <- -sum(rate_matrix[i, ])
}

Now, I would like to try and simulate trajectories from the Continuous Time Markov Chain at different times (e.g. suppose I start at State 1 with probability = 1):

 library(Matrix)
 library(expm)
   
    p0 <- c(1, 0, 0)
    t <- 5.4
    P_t <- expm(rate_matrix * t)
    p_t <- p0 %*% P_t

We can also observe the long term probabilities of being in any state over a period of time:

library(ggplot2)

# define the time range
t_range <- seq(0, 100, length.out=1000)

# compute the state probabilities at each time point
p_t <- t(sapply(t_range, function(t) p0 %*% expm(rate_matrix * t)))

# create a data frame for the plot
df <- data.frame(t = t_range, p1 = p_t[,1], p2 = p_t[,2], p3 = p_t[,3])

# plot the state probabilities over time
ggplot(df, aes(x = t)) +
  geom_line(aes(y = p1, color = "State 1")) +
  geom_line(aes(y = p2, color = "State 2")) +
  geom_line(aes(y = p3, color = "State 3")) +
  labs(x = "Time", y = "State Probability", color = "State")

My Question: Have I understood and done this correctly?

Thanks!

Note: Stationary Distribution Calculations for the Continuous Time Markov Chain

As I understand, the Stationary Distribution of a Markov Process can be considered to be a "magical distribution" such that if this "magical distribution" happens to be the initial distribution of the Markov Process, then the transition dynamics of the Markov Process will reach an equilibrium (i.e. stabilize).

In a Discrete Time Markov Process, the Stationary Distribution can be written as follows: $$\pi P = \pi$$

In a Continuous Time Markov Process, the Stationary Distribution has a similar implications:

   # Calculate the eigenvectors and eigenvalues of the rate matrix
        eig <- eigen(rate_matrix)
       
        # Extract the eigenvector corresponding to the eigenvalue of zero
        pi <- eig$vectors[, which.min(abs(eig$values))]
       
        # Normalize the eigenvector to obtain the stationary distribution
        pi <- pi / sum(pi)

Answer 1

A few comments on the simulation script:

• Currently you don't allow for alternative dynamics (alternative rate matrices Q).
• You shouldn't hardcode variables if you plan on playing around with them later. If you do, you're making it more difficult to change them by hardcoding them in several places.
• Don't grow a vector. Preallocate it. Growing a vector is very slow. You're destroying and re-creating a vector many times. Creating a vector once and accessing its indexes is much faster.
• It is easier to label your states with numeric variables because the labels will also serve as indexes you can subset things by.
• Maybe add some code to verify Q is a proper rate matrix.

Try something like this:

# variables you change a lot go at the top
numJumps <- 100 
numStates <- 3
Q <- matrix(data = c(-2, 1, 1,
                     1,-2, 1,
                     1, 1,-2), 
            byrow = TRUE,
            nrow = numStates, ncol = numStates)



states <- 1:numStates
trajectory <- vector(mode = "numeric", length = numJumps)
jumpDurations <- vector(mode = "numeric", length = numJumps)


trajectory[1] <- 1
jumpDurations[1] <- 0
for (i in 2:numJumps) {
  
  # Get the current state's row
  currentState <- trajectory[i-1]
  positiveColumns <- Q[currentState,] > 0
  
  # simulate waiting times
  expRvs <- rexp(n = numStates - 1, rate = Q[currentState, positiveColumns])
  
  # store the jump time
  jumpDurations[i] <- min(expRvs)
  
  # based on that amount of time, simulate a new state
  transitionMat <- Matrix::expm( jumpDurations[i] * Q )
  trajectory[i] <- sample(x = states, size = 1, prob = transitionMat[currentState,])

}

plot(cumsum(jumpDurations), trajectory, type = "s")

Regarding the estimation, I recommend writing down the log-likelihood function and stating which algorithm you plan on using before you start writing any code. I’m vaguely aware that there are many techniques to perform maximum likelihood estimation, and it will be easier to discuss if you TeX things up first.