Simulate a continuous time Markov chain with two states (one absorbing)

Let $$X_t$$, $$t\geq 0$$ be a continuous time Markov chain with the two states: $$\{\ast,\dagger\}$$ being $$\ast$$ alive and $$\dagger$$ dead (absorbing). We start with $$X_s=\ast$$, at some point $$s\geq 0$$, then I would like to simulate $$X_t$$, $$t\geq s$$. Eventually, $$X_t$$ will become $$\dagger$$ and stay there forever. I would like to simulate paths of this process or, in other words, random times until $$X_t$$ reaches $$\dagger$$ and stays. That is, life spans.

I am given the transition probabilities: $$p_{\ast\ast}(s,t)=P(X_t=\ast|X_s=\ast)$$, $$s\leq t$$.

How can I do that? I tried diving the interval in steps and generating Bernoulli trials with probabilities $$p_{\ast\ast}(n,n+1)$$, $$n\geq s$$, then stop when I get $$\dagger$$. But then if I do that with say, $$1000$$ trials, I count those which are above say $$L$$, this should coincide with $$1000p_{\ast\ast}(s,L)$$, but I get a slightly significant difference.

My question is: how to simulate life spans and if it is done like I did, if there is a huge discretization error we have to live with? (I tried discretizing more and more but the difference is still there)

• So $p_{**}(s,t)$ is already the probability that the object is alive at time $t$, which means that $1-p_{**}(s,t)$ is the CDF of the lifetime. So you can use the probability integral transformation to sample from this without actually doing anything with CTMCs (since this single number characterizes the entire trajectory).
– Ian
Oct 5, 2020 at 11:48
• That being said, for a time-homogeneous CTMC you can avoid all discretization error (as opposed to floating point error) by splitting the chain into the "jump chain" which says where each jump takes you (obtained by taking ratios of rates) and the "holding times" which are exactly exponentially distributed. For a time-inhomogeneous CTMC this doesn't work, because the holding times aren't exponentially distributed anymore. Thus you get locked into a procedure more or less like what you did here.
– Ian
Oct 5, 2020 at 11:57
• What you did here can require an annoyingly fine discretization in order to resolve the holding times if the rates in the CTMC don't change all that fast. If the rates don't change too fast, then you can approximate the process as being a Poisson process on each subinterval, and then once you have decided that the object dies in a particular subinterval then the moment of death will be uniformly distributed within that subinterval (within the Poisson process approximation).
– Ian
Oct 5, 2020 at 11:57
• Let me know whether that covers your question, if so I will copy it into an answer.
– Ian
Oct 5, 2020 at 11:58
• Oh of course, yeah, that is a quite straight forward way, just simulating the time $T_s$ with distribution function $F_s(x)=1-p_{\ast\ast}(s,x)$. I thought about it actually, but it is not possible to invert $F_s$ when things are not stationary and $mu$ is complicated. Nevertheless, I used $R$ to invert it numerically and it works, so I would accept it as an answer, thanks!! In general, with many states I assume that what I did by discretizing is the only easy way then right? Discretizing and generating multinomial random variables. Thanks for your help! Oct 6, 2020 at 6:04

In your simple two state case where one of the states is absorbing, the entire trajectory is determined by the absorption time, which is just one random variable. Thus you can just sample from the absorption time distribution. Specifically, $$1-p_{**}(s,t)$$ is the CDF of the time of absorption, which you can sample from by using the probability integral transformation. (Specifically it suffices to solve $$p_{**}(s,t)=U$$ for $$t$$ given $$U$$ uniformly distributed on $$(0,1)$$.) This can be done numerically if you can't actually do the inversion analytically.
As for more general CTMCs, in the time-homogeneous case, you can avoid all discretization error by splitting the chain into the "jump chain" (a DTMC that tells you where each jump takes you, which is constructed by taking ratios of rates) and the "holding times" which are exactly exponentially distributed and independent of one another. The jump chain and holding times can then be simulated; the jump chain doesn't need anything else to be simulated, while knowing the distribution of the $$(n+1)$$th holding time requires you to know where the $$n$$th jump took you. (Thus you can simulate them in lockstep, or simulate the jump chain and then the holding times, or perform both simulations together and allow the holding times to lag behind the jump chain. In theory all three of those implementations will behave the same, though particular realizations will behave differently, even if you use a seeded PRNG for your noise source.)