The random variable $1+z_1+z_1z_2+\dots$ appears in many contexts and was shown by Kesten to exhibit a heavy tail distribution. We consider natural extensions of this variable and its associated recursion to $N \times N$ matrices either real symmetric $β=1$ or complex Hermitian $β=2$. In the continuum limit of this recursion, we show that the matrix distribution converges to the inverse-Wishart ensemble of random matrices. The full dynamics is solved using a mapping to $N$ fermions in a Morse potential, which are non-interacting for $β=2$. At finite $N$ the distribution of eigenvalues exhibits heavy tails, generalizing Kesten's results in the scalar case. The density of fermions in this potential is studied for large $N$, and the power-law tail of the eigenvalue distribution is related to the properties of the so-called determinantal Bessel process which describes the hard edge universality of random matrices. For the discrete matrix recursion, using free probability in the large $N$ limit, we obtain a self-consistent equation for the stationary distribution. The relation of our results to recent works of Rider and Valkó, Grabsch and Texier, as well as Ossipov, is discussed.