I'm working on the following birth-death process $$X\rightleftharpoons A$$ in which the birth and death (transition) rates are $t^+=k_2a$ and $t^-=k_1n$ respectively. $n$ is the concentration of $X$ and $a$ is a constant concentration of $A$.
I wrote down a continuous time Monte Carlo code in MATLAB that shows the number of particles as a function of time, $n(t)$. What I'm trying to do now is find the first passage time distribution to get to $n=0$ starting from some initial number of particles $n(0)=N$ and I'm not sure on how to do this.
What I tried doing is the following: I ran the program for many different realizations and in each realization I found the time at which $n=0$ for the first time and saved this realization and time. I was hoping I could extract some information from this on the distribution but I don't know how to continue on from here. I don't know if it's even the right approach to this question...
Help would be very much appreciated.
Let $X_t$ denote a continuous time Markov chain on state space $\mathbb{Z}_+:=\{0,1,2,\dotsc\}$ and transitions $n\to n+1$ at rate $\lambda$ and $n\to n-1$ at rate $\mu n$. That is $q_{n n+1}=\lambda$, $q_{nn-1}=n\mu$, where $q_{ij}$ are the entries of the transition rate matrix, $Q$ and assume $X_0=N$ a.s. Note, the average holding time in state $i\in\mathbb{Z}_+$ is $1/q_i$ where $q_i:=-q_{ii}=-(\lambda+i\mu)$. We also assume $0$ is an absorption state, so that once we enter that state we cannot leave (no births occur at zero population). Please indulge me in the change of notation.
Now set $T^j=\inf\{t>0 : X_t=j\}$, the hitting or passage time of state $j$. Then the time of extinction is just $T^0$ (here subscripts are not powers, of course). A first step to extract some information about the distribution is to compute the mean extinction time first. As the standard theory goes, we can compute $\mathbb{E}(T^0)$ by first computing $k_j:=\mathbb{E}_j(T^0):=\mathbb{E}(T^0 | X_0=j)$ for every positive integer $j$. In principle we can always compute $k_j$ by solving the corresponding difference equation: $$k_j=\frac{1}{j\mu+\lambda}+\frac{\lambda}{j\mu+\lambda} k_{j+1}+\frac{j\mu}{j\mu+\lambda} k_{j-1},$$ with "boundary" condition $k_0=0$. All of this and more is covered in J.R. Norris' Markov Chains. Of course, as I said this is just for computing the mean hitting time of a given state. The distribution is another beast but is obtainable sometimes. It is actually not too difficult to instead compute $\mathbb{P}(X_t=0)$ in some cases. For example when both rates are proportional to the population level and the initial population is $1$, i.e. $q_{n n+1}=n\lambda$, $q_{nn-1}=n\mu$, $X_0=1$, then one finds $$\mathbb{P}(X_t=0)=\frac{\mu(e^{\mu t}-e^{\lambda t})}{\mu e^{\mu t}-\lambda e^{\lambda t}},$$ for any given time $t\geq 0$.