Amount of time spent in each state by a Markov chain with matrix of transition times

Question

When we have a simple Markov matrix (rows sum to 1) $P$ where $P_{i,j}$ gives the probability of making a transition from state $i$ to state $j$, We can get the steady state probabilities of being in each of the states ($\pi$) by solving the following system of equations -

$$P \pi = \pi$$

Now, let's say we have a matrix of transition times as well $T$, where $T_{i,j}$ gives the time it takes to make the transition from state $i$ to state $j$. When $T$ has all entries of the same magnitude, it is clear that $\pi$ would not change. If not however, $T$ will also have some impact on $\pi$. Any ideas on how to find the new $\pi$?

Answer 1

I think I figured it out. I don't have a concrete proof, but strong numerical verification (through simulations) for a variety of matrices (which is good enough for me on most days). @Math1000 had the right idea of taking the $\pi$ matrix and modifying it with the transition times matrix.

$\pi$ represents the proportion of times the system enters each state. However, once the system enters state $i$, it spends on average not one unit of time there now. No, it spends $t_i = \Sigma_j P_{i,j}T_{i,j}$. In this way, we can get the $t$ vector. Then, simply element-wise multiply $t$ with $\pi$ and renormalize the result (to make it sum to one). And that is it.

The Python code I used to verify this is below. It returns two arrays, the simulation result and the closed form calculated through the method outlined. The two are always very close regardless of the $P$ and $T$ matrices input to the routine. The code is not very clean but then again, not too much is going on here.

'''
Given the transition probabilities and transition times matrices, outputs the proportion of time spent in each state via simulation and closed form.
The two should always be close.
'''
def SteadyStateMonteCarlo(
  p = np.matrix([
          [0,.2,.4,.4],
          [.2,0,.4,.4],
          [.2,.3,0,.5],
          [.3,.4,.3,0]
          ]),
    t = np.matrix([
          [2,2,2,1],
          [1,1,5,1],
          [1,1,1,1],
          [1,1,1,1]
          ]),
    starting_state = 2
):
  states = np.zeros(p.shape[0])
  states_w_times = np.zeros(p.shape[0])
  curr_state = starting_state
  for i in range(10000):
    next_state = np.random.choice(p.shape[0], p=np.array(p[curr_state,])[0])
    states[curr_state] = states[curr_state] + 1
    states_w_times[curr_state] = states_w_times[curr_state] + t[curr_state, next_state]
    curr_state = next_state
  p_times_t = np.array(np.sum(np.multiply(p,t),axis=1).T)[0]
  # Solve for pi by finding the null space of (P-I)
  pis = np.linalg.svd(p-np.eye(p.shape[0]))[0][:,p.shape[0]-1]
  pis = pis/sum(pis)
  res = [states_w_times/sum(states_w_times), pis, p_times_t]
  props1 = np.multiply(res[1].T,res[2])/sum(np.multiply(res[1].T,res[2]).T)
  props2 = res[0]
  return [props1, props2]

Answer 2

Let $X:=\{X_n:n=0,1,\ldots\}$ be the original Markov chain and $Y:=\{Y_n:n=0,1\ldots\}$ the process obtained by adding the transition times. $Y$ is not a Markov chain, because at time $n$, the $(n+1)^{\mathrm{th}}$ transition time depends on the $(n+1)^{\mathrm{th}}$ state. So we cannot take the approach of finding a stationary distribution for a transition matrix.

Let $E$ be the state space of $X$, then $Y$ takes values in $E\times\mathbb N$. That is, for each $n$ we have $Y_n=(X_n, J_n)$, where $J_n$ is the $n^{th}$ jump time. For example, if $X_0 = i_0$, $X_1 = i_1$, and $T_{i_0,i_1} = 2$, then $Y_0 = (i_0,2)$ and $$ \mathbb P(Y_1 = (i_1, k)) = \sum_{j\in E} kP_{i_1,j}\cdot \mathsf 1_k(T_{i_1,j}). $$ $Y$ is a Markov renewal process, a generalization of Markov chains (and of renewal processes!). For each $n$, let $$\nu_{i,n}=\left(\frac1{1 + \sum_{m=0}^n J_m}\right)\sum_{m=0}^nJ_m\cdot \mathsf 1_i(X_m)$$ be the time spent in state $i$ after $n$ jumps. We want to find $\nu_i :=\lim_{n\to\infty} \nu_{i,n}$.

I believe (though it remains to be proven), that $$ \nu_i = \pi_i \cdot \left(\frac{\sum_{j\in E} T_{i,j}}{\sum_{k\in E}\sum_{l\in E} T_{k,l}} \right). $$ That is, the limiting fraction of time $Y$ is in state $i$ is given by the limiting fraction of time $X$ is in state $i$ multiplied by the ratio of the mean transition time from state $i$ to the mean transition time between any pair of states.