I want to calculate expected hitting times (or first passage times) in a continuous-time birth-death process.
Context Let $Q$ be the instant transition rate matrix and $S = \{0, 1, \ldots\}$ the state space. Fix $k \in S$ and let $t_i = \mathbb E[T_k \mid X_0 = i]$ be the hitting time of the state $k$ if the process starts from the state $i \in S$. It can be shown (e.g. here [PDF], page 55) that the vector $\mathbf{t} = {(t_0, t_1, \ldots)}^{\mathsf{T}}$ solves the linear system $$Q_{(-k)} \mathbf{t} = {(-1, -1, \ldots)}^{\mathsf{T}}$$ where $Q_{(-k)}$ is the matrix $Q$ without the $k$-th row or column (so it's still square).
Question This is a very useful result when the state space is finite, as the linear system is trivial to solve. In my case, though, $Q$ is infinite but has a convenient tridiagonal form: $$Q = \begin{pmatrix} -\lambda & \lambda\\ \mu & -(\mu + \lambda) & \lambda\\ & 2\mu & (-2\mu + \lambda) & \lambda\\ && \ddots & \ddots & \ddots \end{pmatrix}$$
Alternatively, we may specify $Q = (q_{ij})$, $i, j \ge 0$ where $$q_{ij} = \begin{cases} -(i\mu + \lambda) & j = i\\ \lambda & j = i + 1\\ i\mu & j = i - 1\\ 0 & \text{otherwise} \end{cases}$$
I was wondering if there is a way to obtain closed forms for the elements $t_i$ in this infinite case. I suppose the fact that the system matrix $Q_{(-k)}$ is still tridiagonal could make this problem tractable.
If that's not possible, is there a way to find a closed form for a particular element? In fact, I'm only interested in $t_0$, but I couldn't come up with any other way.
EDIT: I just realized that in order to find $t_0$ it's not necessary to solve an infinite system. This is due to the fact that the time to reach state $k$ in the original (infinite) chain, is the same as in a chain with $k + 1$ states, since it's not possible to hit a state larger than $k$ without hitting $k$ first. This allows us to reduce the computation to a finite system.
I managed to obtain $t_0$ in closed form.
The solution $\mathbf{t}$ is given by $$\mathbf{t} = Q_{(-k)}^{-1} {(-1, -1, \ldots)}^{\mathsf{T}}$$ Since we want to find $t_0$, it's sufficient to calculate the first row of $Q_{(-k)}^{-1}$. Maybe there's a way to do this analytically (by using the fact that $Q_{(-k)}$ is tridiagonal perhaps?), but I couldn't do it. Therefore I calculated the values of $Q_{(-k)}^{-1}$ with SymPy and reverse-engineered a formula. For example, the first values of $t_0$ for $k \in \{1, \ldots, 4\}$ are $$\begin{array}{ll} \hline k & t_0\\ \hline 1 & \dfrac{1}{\lambda}\\[3ex] 2 & \dfrac{2\lambda + \mu}{\lambda^2}\\[3ex] 3 & \dfrac{3\lambda^2 + 3\lambda\mu + 2\mu^2}{\lambda^3}\\[3ex] 4 & \dfrac{4\lambda^3 + 6\lambda^2\mu + 8\lambda\mu^2 + 6\mu^3}{\lambda^4}\\[1ex] \hline \end{array}$$
I verified up to $k = 17$ that the value of $t_0$ is given by $$t_0 = \frac{1}{\lambda^k} \sum_{i = 0}^{k - 1} \frac{k!}{(i + 1)(k - i - 1)!} \lambda^{k - i - 1} \mu^i.$$
This value is consistent with numerical simulations.
I would be very grateful if someone could chime in to remedy my lack of mathematical skill and offer some suggestions regarding the computation of the first row of $Q_{(-k)}^{-1}$. The way I followed, while yielding the (correct) result, is certainly not very mathematical.