Deriving the Distribution of Markov Chain Times

Question

I am interested in learning how to derive the probability distributions for the Time to Absorption in Markov Chains (Discrete and Continuous).

In the past, I have usually done one of the following:

• Use the Fundamental Matrix to directly determine the Expected Time to Absorption of a Discrete Time Markov Chain
• Repeatedly simulate time to absorption for a Markov Chain and plot the distribution of these times (this works for both continuous and discrete)

Now, I am trying to learn how to derive the exact probability distribution for these absorption times.

Loosely speaking, I am aware that Absorption Time Distributions are based on Phase Type Distributions (mixture distributions of exponential/geometric distributions, i.e. weighted sums of exponential/geometric distributions) - but I am not sure how to derive the exact formulas for the probability distributions of absorption times.

Doing some research online, I came across the following references:

1) (Distribution for the Absorption Time of a) Discrete Markov Chain: https://math.stackexchange.com/questions/3281928/what-is-the-distribution-of-time-to-absorption-for-an-absorbing-markov-chain :

$$p(k) = \tau T^{k-1} T_0 $$

Where:

• the last row of the $n\times n$ transition matrix corresponds to the absorbing state

• $T$ be the upper left $(n-1)\times(n-1)$ block.

• $\tau$ is the initial distribution

• $T_0$ is the vector formed by the first $n-1$ entries of the last column of the transition matrix

2) (Distribution for the Absorption Time of a) Continuous Markov Chain: https://www.louisaslett.com/PhaseType/:

$$\begin{array}{rcl} F_X(x) &=& 1 - \boldsymbol{\pi}^\mathrm{T} \exp\{x \mathbf{S}\} \mathbf{e}\\ f_X(x) &=& \boldsymbol{\pi}^\mathrm{T} \exp\{x \mathbf{S}\} \mathbf{s} \end{array}$$

Where:

• $\pi$ is the initital state distribution
• $\mathbf{T} = \left( \begin{array}{cc} \mathbf{S} & \mathbf{s} \\ \mathbf{0}^\mathrm{T} & 0 \end{array} \right)$
• $S$ is a $n$ x $n$ matrix of transition rates

Can anyone show me how I can derive these probability distributions from first principles? Once I have these probability distributions, I can then work towards finding the formulas for the variances.

I am in the process of trying to derive them myself and I can share my progress if needed.

Answer 1

I can show you how to derive the discrete-space result, and the extension to the continuous-space result is essentially a limiting version of this.

---

Distribution of absorption time: The formula you give in your question is for the discrete phase-type distribution corresponding to a Markov chain with a finite number of states and an absorbing state. Without further loss of generality, suppose we have $n$ states in the Markov chain with a single absorbing state (so $n-1$ transient states). Varying your notation slightly, we can write the $n \times n$ transition probability matrix for such a chain as:

$$\mathbf{p} = \begin{bmatrix} \mathbf{T} & \mathbf{T}_0 \\ \mathbf{0}^\text{T} & 1 \end{bmatrix} = \begin{bmatrix} \mathbf{T} & (\mathbf{I}-\mathbf{T}) \mathbf{1} \\ \mathbf{0}^\text{T} & 1 \end{bmatrix},$$

where $\mathbf{T}$ is the $(n-1) \times (n-1)$ upper-left block for the transition probabilities between the transient states and $\mathbf{0}$ and $\mathbf{1}$ are $(n-1) \times 1$ column-vectors of zeroes and ones respectively. (The remaining vector $\mathbf{T}_0 \equiv (\mathbf{I}-\mathbf{T}) \mathbf{1}$ is the $(n-1) \times 1$ upper-right block for the transition probabilities from the transient states to the absorbing state.) You can see from this form of transition matrix that the first $n-1$ states are the transient states and the last state is the absorbing state.

Suppose we let $X_t$ denote the state of the Markov chain at time $t=0,1,2,...$. As is standard in Markov chain analysis, we can use the law of total probability to write the relevant multi-step transition probabilities in terms of the powers of the transition probability matrix, which are given by:

$$\begin{align} \mathbf{p}^k &= \begin{bmatrix} \mathbf{T} & (\mathbf{I}-\mathbf{T}) \mathbf{1} \\ \mathbf{0}^\text{T} & 1 \end{bmatrix}^k \\[6pt] &= \begin{bmatrix} \mathbf{T}^k & (\mathbf{I}-\mathbf{T}^k) \mathbf{1} \\ \mathbf{0}^\text{T} & 1 \end{bmatrix}. \\[6pt] \end{align}$$

Suppose we denote the initial probabilities for the transient states of the chain as the row vector $\boldsymbol{\tau} = [\tau_1,...,\tau_{n-1}]$ where $\tau_i \equiv \mathbb{P}(X_0 = i)$, and we then have the full initial probability vector:

$$\boldsymbol{\tau}_* = \begin{bmatrix} \boldsymbol{\tau} & & 1- \boldsymbol{\tau} \mathbf{1} \end{bmatrix}$$

We can use the law of total probability to get the marginal state probabilities:

$$\begin{align} \begin{bmatrix} \mathbb{P}(X_{k} = 1) & \cdots & \mathbb{P}(X_{k} = n) \end{bmatrix} &= \sum_{i=1}^n \mathbb{P}(X_{k} = j| X_0 = i) \cdot \mathbb{P}(X_0 = i) \\[6pt] &= \begin{bmatrix} \boldsymbol{\tau} & & 1- \boldsymbol{\tau} \mathbf{1} \end{bmatrix} \begin{bmatrix} \mathbf{T}^k & (\mathbf{I}-\mathbf{T}^k) \mathbf{1} \\ \mathbf{0}^\text{T} & 1 \end{bmatrix} \\[6pt] &= \begin{bmatrix} \boldsymbol{\tau} \mathbf{T}^k & & \boldsymbol{\tau} (\mathbf{I}-\mathbf{T}^k) \mathbf{1} + (1- \boldsymbol{\tau} \mathbf{1}) \end{bmatrix} \\[12pt] &= \begin{bmatrix} \boldsymbol{\tau} \mathbf{T}^k & & \boldsymbol{\tau} \mathbf{1} - \boldsymbol{\tau} \mathbf{T}^k \mathbf{1} + 1- \boldsymbol{\tau} \mathbf{1} \end{bmatrix} \\[12pt] &= \begin{bmatrix} \boldsymbol{\tau} \mathbf{T}^k & & 1 - \boldsymbol{\tau} \mathbf{T}^k \mathbf{1} \end{bmatrix} \\[6pt] \end{align}$$

In particular, the probability that absorption has occurred is:

$$\mathbb{P}(X_{k} = n) = 1 - \boldsymbol{\tau} \mathbf{T}^k \mathbf{1}.$$

Now, the absorption time for the Markov chain is the hitting time for the absorbing state $n$, which we can denote as $K = \min \{ k =0,1,2,... | X_k = n \}$. We can easily establish the following logical equivalence:

$$X_k = n \quad \quad \iff \quad \quad K \leqslant k.$$

Using the initial probability vector we have $\mathbb{P}(K = 0) = 1- \boldsymbol{\tau} \mathbf{1}$ which is the probability that the chin starts in the absorbing state. The probability of first absorption at any later time $k>0$ is then given by:

$$\begin{align} \mathbb{P}(K = k) &= \mathbb{P}(K \leqslant k) - \mathbb{P}(K \leqslant k-1) \\[6pt] &= \mathbb{P}(X_{k} = n) - \mathbb{P}(X_{k-1} = n) \\[6pt] &= (1 - \boldsymbol{\tau} \mathbf{T}^k \mathbf{1}) - (1 - \boldsymbol{\tau} \mathbf{T}^{k-1} \mathbf{1}) \\[6pt] &= \boldsymbol{\tau} (\mathbf{T}^{k-1} - \mathbf{T}^k) \mathbf{1} \\[6pt] &= \boldsymbol{\tau} \mathbf{T}^{k-1} (\mathbf{I} - \mathbf{T}) \mathbf{1} \\[6pt] &= \boldsymbol{\tau} \mathbf{T}^{k-1} \mathbf{T}_0. \\[6pt] \end{align}$$

This is the result asserted in the formula you gave in your question for the discrete case. The extension to the continuous case requires some further work; it can be established either through taking relevant limits to the discrete case or directly from the rules for continuous-state Markov chains.

---

Moments of absorption time: The moments of the absorption time can be computed using a recursive method. The $r$th raw moment of the absorption time is:

$$\begin{align} \mathbb{E}(K^r) &= \sum_{k=0}^\infty k^r \cdot \mathbb{P}(K = k) \\[6pt] &= (1- \boldsymbol{\tau} \mathbf{1}) \cdot \mathbb{I}(r=0) + \sum_{k=1}^\infty k^r \cdot \boldsymbol{\tau} \mathbf{T}^{k-1} \mathbf{T}_0 \\[6pt] &= (1- \boldsymbol{\tau} \mathbf{1}) \cdot \mathbb{I}(r=0) + \boldsymbol{\tau} \Bigg[ \sum_{k=1}^\infty k^r \cdot \mathbf{T}^{k-1} \Bigg] \mathbf{T}_0 \\[12pt] &= (1- \boldsymbol{\tau} \mathbf{1}) \cdot \mathbb{I}(r=0) + \boldsymbol{\tau} \mathbf{M}_r, \\[6pt] \end{align}$$

where we use the quantities:

$$\mathbf{M}_r \equiv \Bigg[ \sum_{k=1}^\infty k^r \cdot \mathbf{T}^{k-1} \Bigg] \mathbf{T}_0.$$

The values of $\mathbf{M}_r$ can be computed recursively in $r$ using the method shown in Dayar (2005), which involves LU factorisation of the matrix $\mathbf{I}-\mathbf{T}$.

Answer 2

I tried to solve this myself:

1) Discrete Time Markov Chain

A) First Step Analysis:

Let us consider a discrete-time Markov chain with state space $S = \{0, 1, 2, ..., N\}$ and transition probabilities $P_{ij} = P(X_{n+1} = j | X_n = i)$, where $X_n$ is the state at time $n$.

The concept of Absorption Time is closely related to the Hitting Time - i.e. when is the first time when the Markov Chain enters a certain state, conditional on it being in some state at some initial time.

Let's denote $T_i$ as the hitting time of state $N$, given that we started from state $i$, i.e., the first time the process reaches state $N$ when it starts from state $i$. We want to find $E[T_i]$, the expected hitting time. Or more generally, we can calculate this for starting in any state with respect to some initial distribution:

$$\pi = ( \pi_0, ..., \pi_{N-1})$$

It makes sense to write this as weighted sum of expectations [$\pi_i$ is the initial probability of starting in state $i$, i.e. $\pi_i$ is an element in the $\pi = ( \pi_0, ..., \pi_{N-1})$ vector].

The expected hitting time to state $j$, considering the initial distribution, can be calculated as (need to solve this system of equations):

$$E[T_j] = \sum_{i=1}^{n} \pi_i (1 + \sum_{k \neq j} P_{ik} E_k[T_j])$$

where:

• $\pi_i$ is the initial probability of being in state $i$,
• $P_{ik}$ is the transition probability from state $i$ to state $k$, and
• $E_k[T_j]$ is the expected hitting time to state $j$ starting from state $k$.
• note: if you want to do the analysis conditional on starting in a specific state (e.g. state 1), this is equivalent to an initial distribution with probability 1 of being in that state and all other probabilities are 0, i.e. $\pi = (1,0,0,...)$

I think the absorption probabilities can be expressed in the same way. The expected absorption probability to state $j$, considering the initial distribution, can be calculated as (need to solve this system of equations):

$$\phi_j = \sum_{i=1}^{n} \pi_i (P_{ij} + \sum_{k \neq j} P_{ik} \phi_k)$$

where:

• $\pi_i$ is the initial probability of being in state $i$,
• $P_{ij}$ is the transition probability from state $i$ to state $j$, and
• $\phi_k$ is the absorption probability to state $j$ starting from state $k$.

B) Fundamental Matrix Approach:

In a Discrete Time Markov Chain, all transition matrices (with an absorption state) can be re-written in the following form. For example, suppose we have a 5 state Markov Chain where the 5th state is an Absorbing State:

$$ P = \begin{bmatrix} p_{11} & p_{12} & p_{13} & p_{14} & p_{15} \\ p_{21} & p_{22} & p_{23} & p_{24} & p_{25} \\ p_{31} & p_{32} & p_{33} & p_{34} & p_{35} \\ p_{41} & p_{42} & p_{43} & p_{44} & p_{45} \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix} $$

We can re-write this matrix as:

$$ P = \begin{bmatrix} \color{red}{Q} & \color{blue}{R} \\ \color{green}{0} & I \end{bmatrix} = \begin{bmatrix} p_{11} & p_{12} & p_{13} & p_{14} & p_{15} \\ p_{21} & p_{22} & p_{23} & p_{24} & p_{25} \\ p_{31} & p_{32} & p_{33} & p_{34} & p_{35} \\ p_{41} & p_{42} & p_{43} & p_{44} & p_{45} \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} p_{11} & p_{12} & p_{13} & p_{14} & | & p_{15} \\ p_{21} & p_{22} & p_{23} & p_{24} & | & p_{25} \\ p_{31} & p_{32} & p_{33} & p_{34} & | & p_{35} \\ p_{41} & p_{42} & p_{43} & p_{44} & | & p_{45} \\ - & - & - & - & - & - \\ 0 & 0 & 0 & 0 & | & 1 \end{bmatrix} = \begin{bmatrix} \color{red}{p_{11}} & \color{red}{p_{12}} & \color{red}{p_{13}} & \color{red}{p_{14}} & \color{blue}{p_{15}} \\ \color{red}{p_{21}} & \color{red}{p_{22}} & \color{red}{p_{23}} & \color{red}{p_{24}} & \color{blue}{p_{25}} \\ \color{red}{p_{31}} & \color{red}{p_{32}} & \color{red}{p_{33}} & \color{red}{p_{34}} & \color{blue}{p_{35}} \\ \color{red}{p_{41}} & \color{red}{p_{42}} & \color{red}{p_{43}} & \color{red}{p_{44}} & \color{blue}{p_{45}} \\ \color{green}{0} & \color{green}{0} & \color{green}{0} & \color{green}{0} & \color{black}{1} \end{bmatrix} $$

• The red entries represent the $Q$ matrix (transitions between transient states).
• The blue entries represent the $R$ matrix (transitions from transient states to the absorbing state).
• The green entries represent the $0$ matrix (no transitions from the absorbing state to a transient state).
• The black entry represents the $I$ matrix (once the process enters the absorbing state, it stays there with probability

Based on this, we can observe the following pattern:

$$ P^2 = P \cdot P = \begin{bmatrix} Q & R \\ 0 & I \end{bmatrix} \cdot \begin{bmatrix} Q & R \\ 0 & I \end{bmatrix} = \begin{bmatrix} Q^2 + R \cdot 0 & Q \cdot R + R \cdot I \\ 0 \cdot Q + I \cdot 0 & 0 \cdot R + I \cdot I \end{bmatrix} = \begin{bmatrix} Q^2 & Q \cdot R + R \\ 0 & I \end{bmatrix} $$

$$ P^3 = P^2 \cdot P = \begin{bmatrix} Q^2 & Q \cdot R + R \\ 0 & I \end{bmatrix} \cdot \begin{bmatrix} Q & R \\ 0 & I \end{bmatrix} = \begin{bmatrix} (Q^2 \cdot Q) + ((Q \cdot R + R) \cdot 0) & (Q^2 \cdot R) + ((Q \cdot R + R) \cdot I) \\ (0 \cdot Q) + (I \cdot 0) & (0 \cdot R) + (I \cdot I) \end{bmatrix} = \begin{bmatrix} Q^3 & Q^2 \cdot R + Q \cdot R + R \\ 0 & I \end{bmatrix} $$

$$ P^4 = P^3 \cdot P = \begin{bmatrix} Q^3 & Q^2 \cdot R + Q \cdot R + R \\ 0 & I \end{bmatrix} \cdot \begin{bmatrix} Q & R \\ 0 & I \end{bmatrix} = \begin{bmatrix} (Q^3 \cdot Q) + ((Q^2 \cdot R + Q \cdot R + R) \cdot 0) & (Q^3 \cdot R) + ((Q^2 \cdot R + Q \cdot R + R) \cdot I) \\ (0 \cdot Q) + (I \cdot 0) & (0 \cdot R) + (I \cdot I) \end{bmatrix} = \begin{bmatrix} Q^4 & Q^3 \cdot R + Q^2 \cdot R + Q \cdot R + R \\ 0 & I \end{bmatrix} $$

$$ P^k = \begin{bmatrix} Q^k & (I + Q + Q^2 + .... + Q^{k-1}) R \\ 0 & I \end{bmatrix} $$

In the $P^k$ case, we can see that :

• $\lim_{{k \to \infty}} Q^k = 0$ . This is because $Q$ is a matrix of probabilities (i.e. numbers between 0 and 1), thus repeatedly raising these numbers to increasing powers will bring them closer to 0.

• $I + Q + Q^2 + .... + Q^{k-1}$ is a Geometric Series:

• $\text{If } 0 \leq x < 1, \text{ then } \sum_{n=0}^{\infty} x^n = 1 + x + x^2 + x^3 + ... = \frac{1}{1-x}$

• Thus, $I + Q + Q^2 + .... + Q^{k-1}$ can be written as : $(I-Q)^{-1} = \frac{I}{I-Q}$

• Let's call define $N = (I-Q)^{-1}$

We can define the "Fundamental Matrix" as : $$N = (I - Q)^{-1}$$ .

This matrix tells us the number of steps that the Markov Chain spends in transient states, i.e. each element in the Funfamental Matrix $N$ can be thought of as $N_{i,j} = E(T_{i,j})$

Now, we can again write (these formulas make sense practically, but I not am completely sure why they are true):

• Expected Time until Absorption with Initial Distribution:

$$t = \sum_{i=1}^{n} \pi_i t_i = \sum_{i=1}^{n} \pi_i \sum_{j=1}^{n} N_{ij}$$

where $t_i$ is the expected time until absorption from state $i$ and $N_{ij}$ is the $(i, j)$-th entry of the fundamental matrix $N$.

• Absorption Probabilities with Initial Distribution

$$\phi_{j} = \sum_{i=1}^{n} \pi_i \phi_{ij} = \sum_{i=1}^{n} \pi_i \sum_{k=1}^{n} N_{ik} R_{kj}$$

where $\phi_{ij}$ is the probability that the chain will be absorbed in state $j$ given that it started in state $i$, $N_{ik}$ is the $(i, k)$-th entry of $N$, and $R_{kj}$ is the $(k, j)$-th entry of $R$.

Although I have not verified this myself, I think that expected times and probabilities obtained from the Fundamental Matrix approach will be equivalent to the times and probabilities obtained from the First Step Analysis (?). The main advantage being that the Fundamental matrix approach is easier to calculate compared to First Step Analysis (i.e. computing $ (I-Q)^{-1}$ is easier than solving the system of linear equations from First Step Analysis )

C) Distribution of Absorption Times

This is a very "hand-wavy" way I thought of to explain why the Distribution of Absorption Times: Absorption depends on:

• the initial conditions of the Markov Chain
• the probabilities of going from the initial conditions to the non-absorbing part of the Markov Chain (and remaining in the non-absorbing part of the Markov Chain)
• the probabilities of moving from the non-absorbing part of the Markov Chain to the absorbing part of the Markov Chain

Thus, if we want to find out the probability of being absorbed in $k$ steps - this is equivalent to spending $k-1$ steps in the non-absorbing part of the Markov Chain ($Q$) and $1$ step in the absorbing part of the Markov Chain ($R$). Using the notation from the Fundamental Matrix, we can write the probability distribution of being absorbed in $k$ steps (i.e. $k = k -1 + 1$):

$$\sum_{k=1}^{n} P(K=k) = \pi_0 \cdot \sum_{k=1}^{n} \cdot Q^{k-1} \cdot R^1$$

And we can see that the above formula is in the same form as the formula provided here https://math.stackexchange.com/questions/3281928/what-is-the-distribution-of-time-to-absorption-for-an-absorbing-markov-chain (i.e. Discrete Phase Type Distribution). I am not sure how to calculate the variance of this distribution yet.

2) Continuous Time Markov Chains

I will try to keep this short. Define the basics equations of a Continuous Time Markov Chain (Here $Q$ is not referring to the same thing as in the Fundamental Matrix used in the Discrete Markov Chain above).

• Rate Matrix (note: for all absorbing states $i$, the associated rate $q_{i,}$ will be 0 by definition - this is because there is no emission from an absorbing state):

$$ Q_{ij} = \begin{bmatrix} q_{11} & q_{12} & \cdots & q_{1m} \\ q_{21} & q_{22} & \cdots & q_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ q_{m1} & q_{m2} & \cdots & q_{mm} \end{bmatrix} $$ $$q_{ii} = -\sum_{\substack{j=1 \\ j \neq i}}^{m} q_{ij}$$

• Kolmogorov Equation $$ \lim_{{t \to 0}} \frac{{P_{ij}(t)}}{t} = \begin{cases} q_i & \text{if } i = j \\ q_{ij} & \text{if } i \neq j \end{cases} $$

$$\frac{dP(t)}{dt} = P(t)Q$$ $$P(t) = P(t=0) \cdot e^{Qt}$$ $$\exp(Qt) = \sum_{n=0}^{\infty} \frac{(Qt)^n}{n!} = I + Qt + \frac{1}{2!} Q^2t^2 + \frac{1}{3!} Q^3t^3 + \ldots$$

• Embedded Discrete Jump Process:

Holding Rates: $\lambda_i = \lambda_1, \lambda_2, ... \lambda_m$, such that $\lambda_i = \sum_{j} q_{ij}$

(Exponential) Distribution of Holding Times (i.e. Memoryless Property): $f_T(t) = \lambda_i e^{-\lambda_i t}$

Jump Process Probabilities: $P_{ij} = \frac{q_{ij}}{\sum_{j} q_{ij}}$

As an example, for a 4 state Continuous Time Markov Chain with one absorbing state, we can write:

$$ Q = \begin{bmatrix} -q_{11} & q_{12} & q_{13} & q_{14} \\ q_{21} & -q_{22} & q_{23} & q_{24} \\ q_{31} & q_{32} & -q_{33} & q_{34} \\ 0 & 0 & 0 & 0 \end{bmatrix} = \left[\begin{array}{ccc|c} -q_{11} & q_{12} & q_{13} & q_{14} \\ q_{21} & -q_{22} & q_{23} & q_{24} \\ q_{31} & q_{32} & -q_{33} & q_{34} \\ \hline 0 & 0 & 0 & 0 \end{array}\right]$$

$$T = \begin{bmatrix} -q_{11} & q_{12} & q_{13} \\ q_{21} & -q_{22} & q_{23} \\ q_{31} & q_{32} & -q_{33} \\ \end{bmatrix}$$

$$r = -T \cdot \mathbf{1} = \begin{bmatrix} q_{14} \\ q_{24} \\ q_{34} \\ \end{bmatrix}$$

$$ Q = \begin{bmatrix} T_{n \times n} & r_{n \times 1} \\ 0 & 0 \end{bmatrix} $$

As is typically done in Continuous Time Markov Chains, we want to use the Rate Matrix $Q$ to determine the time dependent probabilities of transitioning between states (via the Kolmogorov Equation).

Using the expansion $\exp(Qt) = \sum_{n=0}^{\infty} \frac{(Qt)^n}{n!} = I + Qt + \frac{1}{2!} Q^2t^2 + \frac{1}{3!} Q^3t^3 + \ldots$, we should be able to compactly write (I haven't verified this, its from the youtube video in the references) the Kolmogorov Equation as ($\pi_0$ is the initial distribution):

$$ \pi(t) = \pi_0 \left[ I + \begin{bmatrix} Tt & rt \\ 0 & 0 \end{bmatrix} + \frac{1}{2!} \begin{bmatrix} (Tt)^2 + Ttrt & 0 \\ 0 & 0 \end{bmatrix} + \frac{1}{3!} \begin{bmatrix} (Tt)^3 & (Tt)^2 rt \\ 0 & 0 \end{bmatrix} + \ldots \right] $$

$$ \pi_(t) = \pi_0 \left\{ I_{n \times n} + \begin{bmatrix} \left(Tt + \frac{(Tt)^2}{2!} + \frac{(Tt)^3}{3!} + \ldots\right) & \left(rt + Tt \cdot rt + (Tt)^2 \cdot rt + \ldots\right) \\ 0 & 0 \end{bmatrix} \right\} $$

$$\pi(t) = \pi_0 \left\{ \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix}_{n \times n} + \begin{bmatrix} \left(Tt + \frac{(Tt)^2}{2!} + \frac{(Tt)^3}{3!} + \ldots\right) & \left(rt + Tt \cdot rt + (Tt)^2 \cdot rt + \ldots\right) \\ 0 & 0 \end{bmatrix} \right\} $$

$$ \pi_(t) = \pi_0 \left\{ \begin{bmatrix} \left(1+ Tt + \frac{(Tt)^2}{2!} + \frac{(Tt)^3}{3!} + \ldots\right) & \left(rt + Tt \cdot rt + (Tt)^2 \cdot rt + \ldots\right) \\ 0 & 0 \end{bmatrix} \right\} $$

$$ \pi_(t) = \pi_0 e^{Qt} = \pi_0 \cdot \begin{bmatrix} \sum_{k=0}^{\infty} \frac{(Tt)^k}{k!} & \sum_{k=1}^{\infty} \frac{T^{(k-1)} (t^k) r}{k!} \\ 0 & 1 \end{bmatrix} $$

$$e^{Tt} = \sum_{k=0}^{\infty} \frac{(Tt)^k}{k!}$$

Again, the top left of this matrix refers to transition rates within the transient states, the top right of this matrix refers to transition rates between transient to absorbing states, and the bottom left/right refer to the activity of the Markov Chain once within the absorbing states (i.e. 0 emission rates)

Finally, let's define the initial distribution $\pi_0 = [\alpha, 0]$, where $\alpha = (\alpha_1, \alpha_2, ...., \alpha_m)$ itself is a sub-vector for the probabilities of being in the set of all $m$ number of non-absorbing state at time=0. After some more manipulation (have not fully confirmed this yet), we can define:

The probability of being in any state at time=t (note: this vector should sum to 1, since it contains probabilities): $$\pi(t) = \begin{bmatrix} \alpha e^{Tt} \quad , \quad 1 - \alpha e^{Tt} \end{bmatrix}$$

Thus, the probability distribution of being absorbed at time = t $$\Pi(t) = 1 - \alpha e^{Tt}$$

However, the variance still remains unknown to me.

Conclusion: Thus after much struggle and heartache, I have tried to understand parts of the analysis needed to derive the probability distributions of the absorption times in Discrete Time Markov Chains and Continuous Time Markov Chains. I am open to any feedback/criticism and would love to hear all thoughts and ideas.

References:

• https://math.stackexchange.com/questions/3780106/expected-time-until-absorption-and-variance-of-time-until-absorption-for-absorbi
• https://escholarship.org/content/qt52h7q78h/qt52h7q78h_noSplash_a425672630bd475d97bfa84204576f55.pdf?t=pdajxt
• Spread of number of steps to reach absorbing state in markov chain
• https://www.youtube.com/watch?v=oUwrCxqOQVU

Answer 3

Would love to also understand how closed-form solutions for the variances of these distributions can also be calculated (i.e. non simulations) if possible.

Below is an example of Markov chain with 5 states plus a 6th absorbing state.

The transition matrix for such terminating Markov chain can be described as consisting of several parts

$$T = \begin{bmatrix} {\color{gray} S} & {\color{red} {s}} \\ 0 & {\color{blue} 1} \end{bmatrix}$$

• The transitions between the different non-terminating states. (In gray)
• The transitions towards the absorbing/terminating state. (In red)
• The transitions from the absorbing state to itselve. (In blue)

Let's consider the state-vector $X(k)$ at step $k$, which describes the fractions of the different states that are being occupied. Note that we have $X(k) = X(k-1) \cdot T$, ie. you can compute the $k$-th state vector from the $k-1$-th state vector by multiplying with the transition matrix. By induction we can express it in terms of the begin state $\pi = X(0)$ multiplied by a power of the transition matrix (repeated multiplications)

$$X(k) = \pi T^k$$

The upper left block of $T^k$ will be of the form $S^k$ and is the part of the transitions between the non-terminating states. If we consider a reduced states-vector with only the non-terminating states, $X'(k)$, which is the states-vector $X(k)$ without the terminating state, then these are

$$X'(k) = \pi' {\color{gray}S}^k$$

And from these states $X'(k)$ a fraction ${\color{red}s}$ will terminate so the probability to end in the $k$-th step is

$$P(K_{terminate} = k) = X'(k-1) \cdot {\color{red}s} = \pi' {\color{gray}S}^{k-1} {\color{red}s}$$

---

To compute the raw moments you can use

$$E[K^n] = \sum_{k=1}^\infty k^n P(K = k)$$

E.g.

$$E[K] = \sum_{k=1}^\infty k P(K = k) = \pi' \left[\sum_{k=1}^\infty k S^{k-1}\right] s$$

$$E[K^2] = \sum_{k=1}^\infty k^2 P(K = k) = \pi' \left[\sum_{k=1}^\infty k^2 S^{k-1}\right] s$$

and substitute $$\sum_{k=1}^\infty k S^{k-1} = (S-I)^{-1}(S-I)^{-1}$$

$$\sum_{k=1}^\infty k^2 S^{k-1} = -(S+I)(S-I)^{-1}(S-I)^{-1}(S-I)^{-1}$$

Alternatively, the expressions are easier when we use

$$\begin{array}{rcl} E[K] &=& \sum_{k=0}^\infty P(K > k)\\ E[K^2] &=& \sum_{k=0}^\infty 2k P(K > k)\\ \end{array}$$

and $P(K > k) = \pi' S^{k} \cdot \mathbf{1}$ where $\mathbf{1}$ is a column vector with ones, which effectively sums up the fractions in all the non-terminateds states and that is equal to the probability of not being yet in the terminating state.

Here is an R computer code that demonstrates that these computations give the same result as a simulation (with some small difference because the simulations are only approximate)

### just some funny matrix for transitions
S = matrix( c(0.1,0.3,0.2,0.1,0.2,
              0.3,0.2,0.1,0.1,0.1,
              0.1,0.1,0.2,0.2,0.1,
              0.2,0.1,0.2,0.1,0.3,
              0.1,0.3,0.1,0.1,0.2), 5, byrow = TRUE)

T = cbind(S,1-rowSums(S))
T = rbind(T,c(0,0,0,0,0,1))

T

start = c(0.2,0.2,0.2,0.2,0.2,0)

simulate = function() {
    t = 0
    state = sample(1:6,1, p = start)

    ### simulate transitions untill termination
    while (state < 6) {
       state = sample(1:6,1, p = T[state,])
       t = t+1
    }
    return(t)
}


#### simulations

set.seed(1)
t = replicate(10^5,simulate())

           ### results
mean(t)    # 5.44687
mean(t^2)  # 53.97489

### computations

SI1 = solve(S-diag(c(1,1,1,1,1)))
s = 1-rowSums(S)
EX = t(start[-6]) %*% SI1 %*% SI1 %*% s
EX2 = -t(start[-6]) %*% (S+diag(c(1,1,1,1,1))) %*% SI1 %*% SI1 %*% SI1 %*% s

           ### results
EX         # 5.450171
EX2        # 53.81013