Solving a linear system represented by an infinite matrix (tridiagonal) 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.
 A: 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.
