# First entrance time distribution CTMC

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\}$$?

• You stated that this was a continuous-time Markov chain. What's the distribution of the jump times? Probably something like $\operatorname{Exponential}(1)$? Commented Dec 14, 2019 at 16:11

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...

• Many thanks to you! I will post what i did, but asssume that i am wrong. Maybe you can help me finding the mistakes! Commented Dec 20, 2019 at 7:35
• UPDATE: I found here people.bath.ac.uk/maspm/Prob2B-Part2.pdf also a nice explanation for a similar way via laplace transform Commented Dec 22, 2019 at 13:50
• @stievie Yeah I initially tried doing this via MGFs/Laplace transforms/resolvents but using characteristic functions let me readily appeal to residue calculus in the end so I went with that route. I can sketch more details if you are interested (either the computation of the CF or the resulting integral for the PDF) but it will take a little time. Commented Dec 22, 2019 at 22:14
• Oh that would be great! I am really happy and thankful for the result but also your enormous effort! Since I am not an markov chain expert this would help me a lot seeing more details about finding and computing the CF! Commented Dec 23, 2019 at 21:03

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 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

• Yeah sorry it is difficult to say. Seems like the end computation will be pretty unwieldy even if this works out. Also, don’t you need to force $Y$ to take the same paths of $Z$ up until hitting $V$ instead of just taking any CTMC with the reduced $Q$ matrix? Commented Dec 21, 2019 at 21:09
• This is one of my problems. What do you think about those events in $P(T>t\mid N(t)=2n)$; are they really independent? I don't really know since $S_n$ depends on $T_D^i$... Commented Dec 22, 2019 at 8:59