First entrance time distribution CTMC

Question

Let $Z$ be a continuous time markov chain on $S=\{A,D,V\}$ with $Q$-matrix \begin{pmatrix}-(1+d)& d& 1\\a & -a & 0\\0 & 0 & 0 \end{pmatrix} and $Z(0)=A$. How to find a closed expression for the distribution of the time $$T=\inf\{t\geq 0: Z(t)=V\}$$?

Answer 1

We first note that, at least, the expected hitting time $\mathbb{E}(T_V)=(a+d)/a$ is easily computed by standard Markov chain theory. Now to compute the density, assuming it exists, we shall try to find the CF and then invert it.

If we let $J_1$ be the first jump time, and condition on $Z_0=A$, by studying $$\phi(u):=\mathbb{E}_A(e^{iuT_V})=\frac{d}{d+1} \mathbb{E}_A(e^{iu T_V}\mid Z_{J_1}=D) + \frac {1}{d+1} \mathbb{E}_A(e^{iu T_V}\mid Z_{J_1}=V),$$

we can find the characteristic function of $T_V$ is given by, $$\phi(u)=\frac{a-iu}{a-u^2-iu(1+a+d)}.$$ This is done by exploiting properties of the specific $Q$-matrix given (and the associated jump-chain matrix) and using the strong Markov property on the second jump time but I am omitting details for brevity and it would provide a good exercise for readers (as well as to check this result!).

Thus, by taking the (probability theorist's) inverse Fourier transform, the density of $T_V$, denoted by $f(t)$, is given by, $$f(t)=\frac{1}{2\pi} \int_{\mathbb{R}} \frac{e^{-itu}(a-iu)}{a-u^2-iu(1+a+d)} \, du.$$ For $a=d$, I've empirically deduced that $$f(t)=\frac 1{2c_0} e^{0.5(\sqrt{c_0}-c_1)t} \left[(c_0+\sqrt{c_0}) e^{-\sqrt{c_0}t}+c_0-\sqrt{c_0} \right] \cdot \mathbb{1}_{t>0},$$ where $c_0=(a+d)^2+1=4a^2+1$ and $c_1=1+a+d=1+2a$, but I do not believe this holds when $a\neq d$.

Here is a plot of empirical density function, an exponential PDF with rate $a/(a+d)$, and the above density function for when $a=d=\pi$, with the empirical density function estimated by R's density() function from $N=10000$ simulations of the Markov Chain for time length $t=100$, recording how long it takes to reach $V$ (it is always reached in each simulation due the length of $[0,t]$):

For a sanity check, we can compute sample moments from the simulation, quadrature approximations using the density function, and exact expressions from expanding $\phi(u)$ and/or taking derivatives and scaling by appropriate powers of $i$:

$$\begin{array}{|c|c|c|c|} \hline k & sample & quadrature & CF \\ \hline 1 & 1.961350 & 1.998869 & 2 \\ \hline 2 & 8.318155 & 8.608140 & 8.63662 \\ \hline 3 & 52.761806 & 55.567803 & 56.24736 \\ \hline 4 & 441.649876 & 472.203482 & 488.60583 \\ \hline \end{array}$$

where values in the third column are computed from the exact expressions at $a=d=\pi$, $$\mathbb{E}(T_V)=\frac{a+d}a$$ $$\mathbb{E}(T_V^2)=2!\frac{(a^2+2ad +d^2+d)}{a^2}$$ $$\mathbb{E}(T_V^3)=3!\frac{(a^3+3a^2 d+ad(3d+2)+d(d+1)^2)}{a^3}$$ $$\mathbb{E}(T_V^4)=4!\frac{(a^4 + 4a^3d+ 3a^2d(2d+1)+2ad(2d^2+3d+1)+d(d+1)^3)}{a^4}$$ etc, and these formulas are for any $a,d>0$. All mistakes are my own, and I can provide more details (e.g. on the details of computing the CF) later with some time. I would be interested in the computation of the integral if anyone else has it, since those details eluded me.

---

EDIT: Update 12/19/2019

Using an adopted result from Ahlfors' text to compute integrals $\int_{-\infty}^\infty R(u)e^{-itu}du$ for rational functions $R$ and $t>0$ via residue calculus. This gives (after a bit of work) that the density of the hitting time of $V$ is $$f(t)=\frac{1}{2 c_0} \left[A e^{-0.5t(c-\sqrt{c_0})} + B e^{-0.5t(c+\sqrt{c_0})} \right]\mathbb{1}_{t>0},$$ where $A=c_0-(c-2a)\sqrt{c_0},\, B=c_0+(c-2a)\sqrt{c_0},\, c_0 = c^2-4a,\, c=1+a+d$ and $a,d>0$, for all $t\in \mathbb{R}$. One can check that $f$ is a proper PDF: that it is non-negative and $\int_0^\infty f(t)=1$. The details of this computation and the justification of the modified result from Ahlfors' book probably merit another question, so I am omitting them...

Answer 2

EDIT 20.12.2019 (12/20/19)

I tried the following. For simplicity, let's assume $a=d=1$.

Just consider the markov chain $Y$ on $\{A,D\}$ with Generator $$Q=\begin{pmatrix} -d & d\\ a & -a\end{pmatrix}$$

We consider the number of jumps $N(t)$, where $N(0)=0$. We know then that $N(t)$ is poisson distributed with rate $1$. Distinguish two cases:

$$P(T>t)=\sum_{n=0}^{\infty}\left(P(T>t\mid N(t)=2n)P(N(t)=2n)+P(T>t\mid N(t)=2n+1)P(N(t)=2n+1)\right)$$

so we compute $P(T>t\mid N(t)=2n), P(T>t\mid N(t)=2n+1)$.

If $N(t)=2n+1$, then $Y(t)=D$. By the construction of the orgininal markov chain $Z$, it is not possible to jump into $V$. Hence $P(T>t\mid N(t)=2n+1)=1$.

If $N(t)=2n$, then $Y(t)=A$. Then $T$ has not occured in the past if in after every time point when $Y$ just has jumped in $A$ there was an Exp($1$)-time $T_V^i$ (for jumping in $V$) larger than an Exp($1$)-time $T_D^i$ for jumping in $D$. Since $N(t)=2n$, this has happend $n$ times independently.

Further, before the last jump into $A$ there occured $n$ phases in $A$ and also $n$ phases in $D$ which durations are $2n$ independent Exp($1$) random variables. This sum is now denoted by $S_n$ and is Erlang(2n)-distributed. Now, since the last jumping point into $A$ there is again an Exp($1$)-time $T_V$ for the possibility in jumping into $V$. So $T$ has not occured if $S_n+T_V\leq T$.

Then we get $$P(T>t\mid N(t)=2n)=P(t<S_n+T_V,\bigcap_{i=1}^n\{T_D^i\leq T_V^i\})$$ Now $T_D^i\leq T_V^i=1/2$ by property of Exp-Distribution and $S_n$ is independent from $T_V$ so $S_n+T_V$ is Erlang(2n+1)-distributed, denoted by the random variable $E$. By independency,

$$P(T>t\mid N(t)=2n)=P(t<E,\bigcap_{i=1}^n\{T_D^i\leq T_V^i\})=P(t<E)\frac{1}{2}^n$$

We can finally proceed by plugging in the poisson probabilties. What do you think about this?