I hesitated a bit whether to use another answer window or to edit the old one but finally decided in favor of a new window. If moderators think it is a bad idea, they are welcome to merge.

Also, it is 6:30AM and it promises to be quite a busy day, so I'll tell you what and how to count but will not attempt to do the entire computation myself.

What we are dealing with here is a Markov chain controlled random walk. The general idea is the following. You have a nice (mixing, etc.) finite state discrete Markov chain with states $X_i$ and transition probabilities $p_{ij}$ and the unique equilibrium state $\pi_j$. At each tick of its internal clock one of the transitions occurs and, if we go from $x_i$ to $x_j$, then we add a new independent copy of $\xi_{ij}$ where $\xi_{ij}$ is some fixed set of random variables with values in $\mathbb R^2=(n,t)$ (I'll assume them all nice too). What we want to find is the distribution of the crossing place $n$ of the level $t$ for large $t$.

Now, the markov chain will quickly go to the equilibrium distribution, so the probability of the transition $x_i\to x_j$ is just $P_{ij}=\pi_i p_{ij}$. Look at the moment $N$ on the internal clock. The random walk needs just to know how many steps of each kind have been made to determine the position. Of course, the answer for the mean is $\bar N_{ij}=P_{ij}N$. Unfortunately, we also need to know the deviations. Each $\xi_{ij}=E\xi_{ij}+\eta_{ij}$ where $\eta$'s have $0$ mean. Thus, we know the average position: $(X,T)=N\sum_{i,j}P_{ij}E\xi_{ij}=N\bar\xi=N(\bar n,\bar t)$.

The variance comes from $\eta$'s and from the deviations of $N_{ij}$. The "cross" effect is small, so we can think that the result is the sum of two independent variables $\Xi_1=\sum N_{ij}(E\xi_{ij}-\bar\xi)=\sum N_{ij}\bar\mu_{ij}$ and $\Xi_2=\sum_{i,j}\bar N_{ij}\eta_{ij}$. The second one is a piece of cake: you just add up the variance matrices $V\eta_{ij}$ with given coefficients. In the first case, since the steps are very weakly dependent, we still know that, after dividing by $\sqrt N$, we get close to *some* normal law, but we also need to find the covariance matrix of that law. Honestly writing the scalar products between steps at distance $m$, we get $N$ times
$$
\sum_{i,j}\pi_i p_{ij}\mu_{ij}^T\mu_{ij}+\sum_{m\ge 0}\sum_{i,j,k,l}\pi_ip_{ij}p^{(m)}_{jk}p_{kl}\mu_{ij}\mu_{kl}^T
$$

$$
=\sum_{i,j}\pi_i p_{ij}\mu_{ij}\mu_{ij}^T+
\sum_{i,j,k,l}\pi_ip_{ij}Q_{jk}p_{kl}[\mu_{ij}^T\mu_{kl}+\mu_{kl}^T\mu_{ij}]
$$
where $Q$ is the inverse matrix to $I-P^*$ and $P^*$ is the transition matrix from which the projection to the equilibrium state is removed.

Now, once we have the expectation and the full covariance matrix after $N$ steps, to find the crossing is easy. Choose $N$ so that the mean has the right value of $t$. Then $n$ will be at the mean. You still have the "deviation" gaussian with mean $0$ and known covariance matrix to add. Project it to the line $t=0$ along the general direction of themovement given by the expectations (a fixed linear operator) and compute the variance of this projection.

All quantities in question can be expressed in terms of your functions but I have no time to do it right now. You have just 2 states ($+$ and $-$), so the matrices are easy.

Hope, it helps for now. I'll continue later if you do not finish it today yourself :).

Continuation

What's below should be more than enough for your mathematician friend to check my argument and for you to write a short code that computes the answer from your functions.

First, let's make sure that we speak the same language. I assume that $p_{++}(t)$ means the probability that the process left the state + by the time $t$ with the first jump in the positive direction and similarly for 3 other functions. I also assume that all my vectors are rows (I changed the formulae above to accomodate for that and corrected one idiotic typo). So, we multiply rows by matrices in the eigenvalue problems, etc.

Now, the underlying Markov control has the transition matrix
$$
\mathcal P=\begin{pmatrix}p_{++}(\infty)&p_{+-}(\infty)\\
pp_{-+}(\infty)&p_{--}(\infty)\end{pmatrix}=\begin{pmatrix}a & b\\
c & d\end{pmatrix}
$$
The second eigenvalue is $a+d-1=1-b-c$, the equilibrium state is
$$
\pi=(\pi_+,\pi_-)=(\frac c{b+c},\frac b{b+c})
$$
The transition matrix with the equilibrium part removed is
$$
\mathcal P^*=(1-b-c)\begin{pmatrix}\frac b{b+c} & -\frac b{b+c}\\
-\frac c{b+c} & \frac c{b+c}\end{pmatrix}
$$
so
$$
I-\mathcal P^*=\begin{pmatrix}b+\frac c{b+c} & -b+\frac b{b+c}\\
-c+\frac c{b+c} & c+\frac b{b+c}\end{pmatrix}
$$

Thus $\operatorname{det}(I-\mathcal P^*)=b+c$ and

$$
Q=(I-\mathcal P^*)^{-1}=\frac 1{b+c}\begin{pmatrix}c+\frac b{b+c} & b-\frac b{b+c}\\
c-\frac c{b+c} & b+\frac c{b+c}\end{pmatrix}
$$

Now, I'll identify $+$ and $+1$ and similarly with $-$. Also, the matrix indices are $+$ and $-$ in this order (so $Q_{++}$ is the left top corner).

If we know that we went from state $\varepsilon$ to state $\delta$, then we added the random "space-time" shift $\xi_{\varepsilon\delta}=(\delta,T_{\varepsilon\delta})$ where $P(T_{\varepsilon\delta}<t)=\frac{p_{\varepsilon\delta}(t)}{p_{\varepsilon\delta}(\infty)}$, (if I understood you right, we always end in the $+$ state when jumping to the right and in the $-$ state when jumping to the left) so
$$
ET_{\varepsilon\delta}=\int_0^\infty(1-\frac{p_{\varepsilon\delta}(t)}{p_{\varepsilon\delta}(\infty)})dt
\qquad
ET_{\varepsilon\delta}^2=\int_0^\infty 2t(1-\frac{p_{\varepsilon\delta}(t)}{p_{\varepsilon\delta}(\infty)})dt
$$
Now, we can find
$$
\bar\xi=\sum_{\varepsilon\delta}\pi_\varepsilon p_{\varepsilon\delta}
\xi_{\varepsilon\delta}=(\bar n,\bar T)
$$
and

The mean $\frac{\bar n}{\bar T}t$.

For the variance, we need
$$
\bar\xi_{\varepsilon\delta}=(\delta,ET_{\varepsilon\delta})
$$
and the "deiations"
$$
\eta_{\varepsilon\delta}=(0,T_{\varepsilon\delta}-ET_{\varepsilon\delta})
$$

$$
\mu_{\varepsilon\delta}=\bar\xi_{\varepsilon\delta}-\bar\xi
$$

We have
$$
V\eta_{\varepsilon\delta}=\begin{pmatrix}0 & 0 \\ 0 & VT_{\varepsilon\delta} \end{pmatrix}
$$
where $VT_{\varepsilon\delta}=ET_{\varepsilon\delta}^2-(ET_{\varepsilon\delta})^2$.
This gives the first component of the "unit step variance"
$$
V_1=\sum_{\varepsilon\delta}\pi_\varepsilon p_{\varepsilon\delta} \begin{pmatrix}0&0\\ 0&VT_{\varepsilon\delta} \end{pmatrix}
$$

The second component is gvien by
$$
V_2=\sum_{\varepsilon\delta} \pi_\varepsilon p_{\varepsilon\delta}
\mu_{\varepsilon\delta}^T\mu_{\varepsilon\delta}
$$

$$
+\sum_{\varepsilon,\delta,\rho,\sigma\in\{\pm\}}
\pi_\varepsilon p_{\varepsilon\delta} Q_{\delta\rho} p_{\rho\sigma}
[\mu_{\varepsilon\delta}^T\mu_{\rho\sigma}+
\mu_{\rho\sigma}^T\mu_{\varepsilon\delta}]
$$

The total "unit step variance" matrix is $V=V_1+V_2$. The projection operator is just $(n,t)\mapsto n-\frac{\bar n}{\bar T}t$, so the

final answer for the variance is $(1,-\frac{\bar n}{\bar T})V(1,-\frac{\bar n}{\bar T})^T\frac t{\bar T}$

I hope I haven't screwed anything up but it would be nice if you check this against simulations before believing it 100%. To justify convergence, it is enough to have moments of power $2+\varepsilon$ (I'm not sure if the second moments are enough because the CLT for slightly dependent random variables is a bit harder than for the independent ones and some estimates in the proof I know require extra leeway, but you have exponential tails anyway, so this shouldn't concern you).

If you have any questions (from seeing a disagreement with simulations to just "what is that symbol?"), feel free to ask. If you find that everything works, as said, let me know and I'll put it out of my head :)

`$\mathbb{Z} \times \{+,-\}$`

. $\endgroup$