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.

  • $\begingroup$ Does the state space consist of all integers between and including -n and n? $\endgroup$ Commented Dec 18, 2010 at 9:03
  • $\begingroup$ It looks like the state space is $\mathbb{Z}\times\lbrace +,-\rbrace$. That is, it knows the level it is at and whether the previous transition was up or down. As a Markov process, it has space $\mathbb{Z}\times\lbrace +,-\rbrace\times\mathbb{R}^+$ as it also remembers how long it has been at its current level for. I expect that whether it really does converge to a normal depends on the form for $p_\pm(t)$. $\endgroup$ Commented Dec 18, 2010 at 14:20
  • $\begingroup$ Yes, the state space is $\mathbb{Z} \times \{+,-\}$. $\endgroup$
    – shoyer
    Commented Dec 18, 2010 at 22:14

2 Answers 2


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


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

  • $\begingroup$ This is wonderful -- It would be great if you could fill out a few more of the details, but I think this is enough that one of my mathematician friends could help me figure it out. $\endgroup$
    – shoyer
    Commented Dec 20, 2010 at 2:12
  • $\begingroup$ I did. I apologize for the delay but to drive 500 miles and to type an answer are "two things too hard to do together" :) $\endgroup$
    – fedja
    Commented Dec 20, 2010 at 20:50
  • $\begingroup$ Okay, I'm in the process of comparing to my Monte-Carlo results :) $\endgroup$
    – shoyer
    Commented Dec 21, 2010 at 7:12
  • $\begingroup$ The predicted mean and variance from this answer match the results of my Monte-Carlo simulations! I checked with several different processes of the same form, and in each case I found values deviated by less than 2% from my simulations using 10 minutes worth of CPU cycles. I will add a comment with the arXiv link once my article is up (hopefully within 1-2 months). Fedja, please contact me via my website if you would like to be recognized in the acknowledgments of my manuscript by your full name. $\endgroup$
    – shoyer
    Commented Dec 21, 2010 at 10:33
  • $\begingroup$ 2% deviation is to be expected: a few steps were done neglecting terms of order $\sqrt N$ in the asymptotics linear in $N$, and the relative deviation of the estimator of the variance is of order $M^{-1/2}$ itself where $M$ is the number of trials, so, for finite times, I wouldn't be surprised with the relative error of order $N^{-1/2}+M^{-1/2}$ with some coefficient between 1 and 10 and I doubt you were able to run much more than 100000 simulations with 100000 steps in each. No "recognition" is necessary beyond mentioning MO (it was just a computation), but I'm glad it worked :). $\endgroup$
    – fedja
    Commented Dec 21, 2010 at 17:04

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

  • $\begingroup$ Thanks fedja -- What does $p$ mean in your answer? The transition probabilities $p_{++}(t)$, $p_{+-}(t)$, $p_{-+}(t)$, $p_{--}(t)$ are each different functions of $t$. To clarify, they cannot be written in the form $p_{++}(t) = p f(t)$ and $p_{+-}(t) = (1-p) f(t)$. Also, the functions $p_{-+}(t)$ and $p_{--}(t)$ are different from $p_{++}(t)$ and $p_{+-}(t)$, which I realize was not clear in my question! Oops! $\endgroup$
    – shoyer
    Commented Dec 18, 2010 at 22:25
  • $\begingroup$ Ah, OK. I thought you have just normal integers and go right with probability $p$ and left with $(1-p)$ but the waiting time is random. So, I was quite confused, as it turns out. I'll think of it a bit more now, when the question is clear. I'll try to adjust everything to how you want it, but it'll take some time :). $\endgroup$
    – fedja
    Commented Dec 19, 2010 at 1:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.