Every transition matrix on a nonempty finite state space has a stationary distribution

Fix the state $x_1$ and let $\mu$ be its point mass. For every $N\ge 1$ define the row vector ${\nu }_N:=\frac{1}{N}{\sum }_{k=0}^{N-1}\mu p^{(k)},$ where $p^{(k)}$ is the $k$-th matrix power. Its entries are finite nonnegative sums of products of the entries of $p$, hence ${\nu }_N(y)\in [0,1]$ for every $y\in E$; and using [F1] twice, ${\sum }_{y\in E}{\nu }_N(y)=\frac{1}{N}{\sum }_{k=0}^{N-1}{\sum }_{y\in E}\mu p^{(k)}(y)=\frac{1}{N}{\sum }_{k=0}^{N-1}1=1$, since a point mass is a probability vector and rows of every $p^{(k)}$ sum to one. So every ${\nu }_N$ lies in $S:=\{\rho :E\to [0,1]:{\sum }_y\rho (y)=1\}$.

For every $N\ge 1$ the exact telescoping identity ${\nu }_{N}p-{\nu }_N=\frac{1}{N}{\sum }_{k=0}^{N-1}(\mu p^{(k+1)}-\mu p^{(k)})=\frac{1}{N}(\mu p^{(N)}-\mu )$ holds coordinatewise as an identity of finite real sums, its entries being differences of numbers in $[0,1]$, so no infinite sum is rearranged.

There is a strictly increasing sequence $N_j$ and a vector $\pi :E\to [0,1]$ such that ${\nu }_{N_j}(y)\to \pi (y)$ for every $y\in E$. Enumerate $E=\{x_1,\dots ,x_n\}$ and argue by induction on the number $i$ of coordinates already handled: the $i$-th coordinate sequence along the subsequence produced so far is bounded in $[0,1]$, so [F3] supplies a further strictly increasing subsequence on which it converges; after finitely many successive subsequence choices, every coordinate converges. Finite induction on these existential choices requires no choice axiom.

Consequently ${\nu }_{N}p-{\nu }_N\to 0$ coordinatewise as $N\to \mathrm{\infty }$: each entry of $\frac{1}{N}(\mu p^{(N)}-\mu )$ has absolute value at most $\frac{2}{N}$, since every entry of $\mu p^{(N)}$ and of $\mu$ lies in $[0,1]$.

The limit $\pi$ is a probability vector: $\pi (y)\ge 0$ for every $y$ because a limit of nonnegative numbers is nonnegative; and ${\sum }_{y\in E}\pi (y)=1$ because a finite sum of convergent sequences converges to the sum of the limits, applied to the constant sums $1$ from step 2.1 and step 1.1.

Passing to the subsequence of step 2.1, ${\nu }_{N_j}\to \pi$ coordinatewise and hence ${\nu }_{N_j}p\to \pi p$ coordinatewise, because each entry of the finite matrix product is a finite sum ${\sum }_{y\in E}{\nu }_{N_j}(y)p(y,z)$ of finitely many convergent sequences. By step 2.2 the left side of ${\nu }_{N}p-{\nu }_N$ tends to $0$ along $N_j$, so $\pi p-\pi =0$, that is, $\pi p=\pi$.

By step 3.1, $\pi$ is a probability vector and by step 3.2 it satisfies $\pi p=\pi$; [F2] then identifies it as an invariant probability distribution for $p$. The state space was required nonempty so that $x_1$ exists; the empty matrix has no probability vector, so $E=\mathrm{\varnothing }$ is excluded by the hypothesis. If $n=1$, then $p(x_1,x_1)=1$ and $\pi ={\delta }_{x_1}$ is invariant, consistent with the construction. The argument uses only finite enumerations, finite sums and the subsequence theorem; no irreducibility, aperiodicity, recurrence, product-space path law or choice principle is used, and the conclusion is a one-way existence assertion rather than an equivalence.