# Stationary distribution for time-inhomogeneous Markov process

I have a two state, discrete time, time-inhomogeneous Markov process with transition matrix defined by $$T_i=\begin{pmatrix} 1-p_i\alpha & p_i\alpha \\ p_i\beta& 1-p_i\beta \end{pmatrix}$$ where $p_i$ is a random variable and $\alpha,\beta$ are both defined on $[0,1]$. In the application I have in mind, it would be reasonable to assume that the $p_i$ are iid in a probability distribution on $[0,1)$ or (probably) a small subset.

I have done some numerical experiments in SAGE and I believe that this process has a steady state/equilibrium/stationary distribution (whatever you want to call it), but I am having trouble finding literature on the existence of steady states for time inhomogeneous Markov processes. If it helps, I am reasonably sure that one could also apply the adjectives irreducible, aperiodic, and possibly ergodic to this Markov process.

I have found a theorem that says that a finite-state, irreducible, aperiodic Markov process has a unique stationary distribution (which is equal to its limiting distribution). What is not clear (to me) is whether this theorem is still true in a time-inhomogeneous setting.

Does anyone know of a theorem I could reference that gives the existence of this solution? I am out of my area of expertise, and in this case I would just need the reference. Alternatively, if you know of a more hack-and-slash way to determine if an equilibrium exists in this case (eigenvalues? limits of products of random matrices?), I'm willing to do some computations!

• This is secretly a time-homogeneous Markov process if the p_i are iid. – guest Jun 22 '15 at 20:23
• Consider a time-inhomogeneous Markov process $X_t$ on a finite state space. Let $Q(t)$ denote the generator, and let $P(s,t)$ denote the corresponding transition kernel, i.e. $P(s,t) = U^{-1}(s)U(t)$, where the Markov propagator is $$U(t) := \mathcal{TO}^* \exp \int_0^t Q(s) \ ds,$$ and $\mathcal{TO}^*$ indicates the formal adjoint or reverse time-ordering operator. Now an initial distribution $\pi(0)$ is propagated as $\pi(t) = \pi(0)U(t)$. (See, e.g., Kleinrock, L. Queueing Systems, vol. 1. Wiley (1975).) – Steve Huntsman Jun 22 '15 at 21:18
• You might also find the multiplicative ergodic theorem helpful. – Steve Huntsman Jun 22 '15 at 21:28
• If (in addition to $0 < \alpha, \beta < 1$) $p_i$ are bounded away from zero and one, then a subsequence of the products of matrices (starting at any $i$) will converge to a rank one matrix, hence the system will have unique harmonic function---which is what you want (I think). In fact, this can still happen if $p_i$ isn't likely to get too close to zero. A stationary distribution is then plausible–transpositions don't occur, so telescoping convergence can probably be replaced by convergence. Two extreme harm fns (so no stationary distns) occur if $p_i$ is summable (not iid). – David Handelman Jun 22 '15 at 22:51
• @c.williams This isn't an answer or reference or proof and it doesn't directly speak to time homogeneity, but one thing to think about is what happens if you draw a probability according to a known distribution (with a known expectation -- hint hint) and you subsequently flip one coin controlled by that probability. Can this procedure be simplified? Another question is what stationary distribution are you seeing empirically and why is it proportional to (beta, alpha) regardless of the iid distribution on p (possibly subject to some fine print about bounding)? – guest Jun 24 '15 at 14:03

This answer is just spelling out what guest already said: If $p_k$ is iid over time $k \in \{0, 1, 2, \ldots\}$ then your system is equivalent to a discrete time homogeneous Markov chain with a fixed transition probability matrix: $$P =\begin{pmatrix} 1-E[p]\alpha & E[p]\alpha \\ E[p]\beta& 1-E[p]\beta \end{pmatrix}$$ where $E[p]$ is the expectation of $p_1$. For a proof, you can compute the upper-left component of the matrix as follows (assuming $p$ has a density $f(p)$ for simplicity): \begin{align} Pr[X_{k+1}=1|X_k=1] &= \int_{p=0}^1 Pr[X_{k+1}=1|X_k=1, p_k=p]f(p)dp \\ &=\int_{p=0}^1 (1-p\alpha)f(p)dp \\ &= 1-E[p]\alpha \end{align} and so this answer is the same for all time slots $k$. You can compute similar values for the other three entries of the matrix.
Assuming that $\alpha E[p] \in (0,1)$ and $\beta E[p] \in (0,1)$, a steady state distribution exists: $$\lim_{k\rightarrow\infty} Pr[X_k=1] = \frac{\beta}{\alpha+\beta} \: \: , \: \: \lim_{k\rightarrow\infty} Pr[X_k=2] = \frac{\alpha}{\alpha+\beta}$$ In particular, the specific value of $E[p]$ (and the distribution of $p$) have no influence on the steady state distribution. However, the value of $E[p]$ affects the convergence time to this distribution: It scales the time we spend in each state.