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)
Thanks for your help!
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.)
But for a time-inhomogeneous CTMC, the holding times are not exactly exponentially distributed. So you get stuck doing something kind of like what you did here where you discretize time and work in subintervals. In some cases you may want to resolve the distribution to a very fine scale, so fine that the rates don't change much on the scale that you're trying to resolve. In this case you can try to discretize time until you find an interval in which an event occurs and then neglect the temporal variation in the rates on that subinterval. In this case the jump will occur at a time that is uniformly distributed within the subinterval. This amounts to assuming that your discretization of time is so fine that your process can be approximated by a Poisson process within each subinterval, which is not always a great assumption but it is certainly better than assuming the jump occurs at an endpoint.