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!
 A: 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. 
