A stochastic matrix has a stationary distribution

The simplex is invariant. Let $\Delta=\{\pi\in\mathbb{R}^n:\pi_i\ge0,\sum\pi_i=1\}$, a compact convex set. If $\pi\in\Delta$ then $\pi P\ge0$ and $\sum_j(\pi P)_j=\sum_i\pi_i\sum_jP_{ij}=1$, so $\pi P\in\Delta$.

Cesàro averages. Let $\pi_0=(\frac1n,\dots,\frac1n)$ and
$$m_N=\frac1N\sum_{k=0}^{N-1}\pi_0P^k\in\Delta$$
(an average of points of the convex set $\Delta$). Then
$$m_NP-m_N=\frac1N\Big(\sum_{k=1}^{N}\pi_0P^k-\sum_{k=0}^{N-1}\pi_0P^k\Big)=\frac{\pi_0P^N-\pi_0}N.$$
Since $\pi_0P^N,\pi_0\in\Delta$ have $1$-norm $1$, $\|m_NP-m_N\|_1\le\frac2N\to0$.

Compactness. $\Delta$ is compact, so a subsequence $m_{N_j}$ converges to some $\pi\in\Delta$. By continuity of $x\mapsto xP-x$,
$$\pi P-\pi=\lim_j(m_{N_j}P-m_{N_j})=0.$$

Example. $P=\begin{pmatrix}0&1\\1&0\end{pmatrix}$: $P^k$ does not converge, but the averages $m_N\to(\frac12,\frac12)$, which is stationary.

Remark. This is the Markov-Kakutani / Brouwer-type argument. The vector $\pi$ is generally not unique (e.g. $P=I$) and $\pi_0P^k$ need not converge, which is why Cesàro averaging is needed; under the positivity assumption of the next problem, uniqueness and convergence hold.