Properties of a continuous-time semi-Markov process as t -> \infty

Question

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.

Answer 1

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

Answer 2

If you condition on the (random) number of steps $n$ made by the time $t$, you'll see that the distribution is just a mixture of biased random walks with different number of steps. Each particular random walk is close to normal distribution and you'll get the whole thing asymptotically normal if you manage to show that the distribution of the number of steps is close to normal with variance of order $t$ (then the convolution of two normals is normal again). This is (almost) necessary and sufficient.

Now, let's see when it happens. We have just one clock now and we need to know the distribution of the number $n$ of ticks by time $t$. What we expect is that it is normal with the mean $t/ET$ where $T$ is the random time between ticks. Now, to say that $n>N$ is the same as to say that in $N$ ticks we are below $t$. The distribution of the time of the $N$'th tick is approximateny normal with variance $NVT$ where $VT=ET^2-(ET)^2$ is the variance of $T$ if and only if $ET^2<+\infty$. The rest is a straightforward computation. It shows that, in the case of finite second moment of waiting time $T$, the final distribution is approximately the sum of two independent normals $A+B$ where $A$ has mean $\frac{2p-1}{ET}t$ and variance $\frac{4p(1-p)}{ET}t$ and $B$ has mean $0$ and variance $\frac{V(T)(2p-1)^2}{(ET)^3}t$, i.e., a normal with mean $\frac{2p-1}{ET}t$ and variance $\frac{4p(1-p)(ET)^2+V(T)(2p-1)^2}{(ET)^3}t$ (provided I haven't made any stupid errors in arithmetic, of course).