I am interested in calculating properties of a continuous-time random walk problem which I believe is a type of semi-Markov process.


I have states of the form $n\_\pm \in \mathbb{Z} \times \{ +, -\}$. For a state $n\_\pm$, I have time-dependent transition probabilities $p\_{\pm+}(t)$ and $p\_{\pm-}(t)$ for jumping to the states $(n+1)\_-$ and $(n-1)\_+$ respectively, where $t$ is the time than the walker has remained in the state. To be clear, each of the functions $p\_{++}(t)$, $p\_{+-}(t)$, $p\_{-+}(t)$ and $p\_{--}(t)$ is different and they do not share the same time-dependence.


The time between jumps is almost exponentially distributed as in a Markov process, but not quite (there are deviations at small times). I am interested in the long-time distribution of this process as $t \to \infty$. 


My Monte-Carlo simulations suggest that the probability distribution converges to a normal distribution over $n$ with parameters of the form $\mu = v t$ and $\sigma^2 = 2 D t$. I would like to be able to calculate $v$ and $D$ directly from my (numerical) transition probabilities without extracting them from such a simulation.


I apologize in advance for any informal language -- I am a physicist.



Properties of a continuous-time semi-Markov process as t -> \infty
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 :)


