How to discretely stochastically simulate a continuous-time Markov chain? A continuous-time markov chain describes a continuously varying process, such that future state only depends on the current state.
A sampling of a continuous markov chain can be described in terms of how long it stays in each state. For example, with the three states $A$, $B$, and $C$, we can imagine a sampling to be
Stay in $A$ for $1$ second, stay in $B$ for $.7$ seconds, stay in $C$ for $3$ seconds, stay in $B$ for $2$ seconds, stay in $D$ forever.
My question is, how does one stochastically simulate this? Namely, what stochastic algorithm can produce a series of events of the form $(\text{Chain State}, \text{time})$, such that the probability of a certain sequence of events matches the probability of that behavior in the continuous markov chain?
 A: Say your generator matrix is:
$$ G =
        \begin{bmatrix}
        -\lambda_0 & \lambda_{01} & \lambda_{02}  & \cdots \\
        \lambda_{10} & -\lambda_{1} & \lambda_{12}  & \cdots \\
        \lambda_{20} & \lambda_{21} & -\lambda_{2} & \cdots \\
        \cdots \\
        \end{bmatrix}
$$
where, as usual, for each $i,\;$ $\lambda_i = \sum\limits_{k\neq i}\lambda_{ik}$.
Each pair generated $(S_n,T_n) = (j,t)$, for the $n^{th}$ transition depends on the previous state $S_{i-1}=i$ as follows:
$T_n\sim Exp(\lambda_i)$ so draw $t$ randomly from that distribution.
Independently to that, $P(S_n=k)=\lambda_k/\lambda_i$ for $k\neq i$ so draw $u$, say, randomly from a Uniform$(0,1)$ distribution. Choose the $j$ value for which
$$\dfrac{\sum\limits_{k\neq i}^{j-1} \lambda_k}{\lambda_i} \lt u \leq \dfrac{\sum\limits_{k\neq i}^{j} \lambda_k}{\lambda_i}.$$
Essentially, for each row/state $i$, you have a partition of the $(0,1)$ interval with sub-interval sizes proportional to the values $\lambda_{ij},\; j\neq i$.
