# How can I efficiently approximate the stationary distribution of an infinite CTMC with a sparse rate matrix?

I am looking for methods to approximate the stationary distribution of an infinite CTMC with a sparse rate matrix. Each row and column of the rate matrix has a finite number of non-zero elements. There is this paper Bounding the equilibrium distribution of Markov population models by Dayar et al. that identifies a region of the state space where most of the probability of the stationary distribution is concentrated and then finds upper and lower bounds on the stationary probabilities for the states in that region.

I wondered whether there was something similar for other classes of CTMC or something a bit more general.

A lot of what I've found has been applicable to various classes of queue systems but my case doesn't fall within that.

• you mention that "my case doesn't fall within that", but you do not really tell us what your case is; I guess you'll get a better response if you can be more specific. May 31, 2016 at 14:05

Since each row $r$ has a finite number of non-zero elements, say $N_r \subset \{1, 2, \dots\}$, you could use the Gillespie algorithm to simulate the continuous time Markov chain. A stationary distribution can then be approximated by the normalized occupancy times. If the matrix is irreducible, this will be the stationary distribution.

Here is an outline of how this would work:

Gillespie algorithm

1. Given current state $X_t = r$, simulate an exponential random variable of rate $\lambda_r = \sum_{c \in N_r, c\neq r} q_{r,c}$. This corresponds to the time $t' = t + \lambda_r$ where the current state jumps to a different one.
2. Pick a new distinct state $X'_{t'}$ from a categorical random variable with probability $P(X_{t'} = c|X_{t'-} = r) = q_{r,c}/\lambda_r$ for $c \neq r$.
3. Repeat until a trajectory of fixed length $T$ has been produced.

Occupancy time

The approximation of a stationary distribution is then given by

$$\pi(x) \approx \frac{\text{length}\{t \in [0, T) : X_t = x\}}{T}.$$

• This Monte-Carlo technique is also known as the Doob-Gillespie or Stochastic Simulation Algorithm. If the process $\{ X_t \}$ has a stationary distribution and satisfies a strong LLN then: $\pi(x) = \lim_{t \to \infty} \frac{1}{t} \int_0^t \mathbf{1}_{\{X_s = x\}} ds$. This quantity basically gives the fraction of time spent by the (infinitely long) trajectory $\{ X_t \}$ in state $x$. Oct 2, 2016 at 13:51