Simulating death processes

Let $\{N(t)\}$ be a death process with the death rate $\mu$ and initial population size $N(0)=n_0$. It is characterized by the following differential-difference equations for the probability mass function $p_n(t)=Pr(N(t)=n)$:

$$\begin{array}{}\text{(Eq:DeathN)}& \frac{d{p}_n(t)}{dt}& =& \mu (n+1)p_{n+1}(t)-\mu n{p}_n(t),(0\le n\le n_0-1)\text{(Eq:DeathN0)}& \frac{d{p}_{n_0}(t)}{dt}& =& -\mu n_0{p}_{n_0}(t)\end{array}$$

with the initial condition

$$p_n(0)={\delta }_{n,n_0}.$$

See also: Death process

Here, we want to simulate them numerically to obtain some concrete sample paths like this one:

We can do so just like we did for the birth process.

That is, we first discretize the time variable so that we consider time steps with some small interval $\delta t$: $N(0),N(\delta t),N(2\delta t),N(3\delta t),\cdots ,N(k\delta t),\cdots$. 

Next, we exploit the Markov property of the birth process: Every death is the first death since the last one. That is, given $N(t)=n_t$ at time $t$, we may regard this as the "initial" condition for the next step: $N(t+\delta t)=n_t$ or $n_t-1$. We can solve (Eq:Death0) (with $n_0$ replaced with $n_t$) to find $p_{n_t}(t+\delta t)$ with the "initial" condition $p_n(t)={\delta }_{n,n_t}$. A bit of exercise gives

$$p_{n_t}(t+\delta t)=e^{-\mu n_t\delta t},$$

which is the probability that $N(t+\delta t)=n_t$ (i.e., no death during $\delta t$). Accordingly, we have the probability that $N(t+\delta t)=n_t+1$ (i.e., one death during $\delta t$) as

$$p_{n_t-1}(t+\delta t)=1-e^{-\mu n_t\delta t}.$$

Thus, the algorithm for simulating the birth process is the following:

1. Set $N(0)=n_0$; set k = 0.
2. Generate a uniformly distributed random number $r\in [0,1)$.
3. If $r<1-e^{-\mu \delta tN(k\delta t)}$, then set $N((k+1)\delta t)=N(k\delta t)-1$; Otherwise, set $N((k+1)\delta t)=N(k\delta t)$.
4. Update $k:=k+1$, go to Step 2, and repeat.