Why do sequences $(x_n)$ of real numbers defined by $x_{n+1}=1+\frac 1{x_n}$ nearly always converge to the golden ratio?

The map $x\mapsto 1+{1\over x}$ defined on ${\mathbb R}\setminus\{0\}$ can be extended to a Moebius transformation of the Riemann sphere $\bar{\mathbb C}:={\mathbb C}\cup\{\infty\}$: $$T:\quad \bar{\mathbb C}\to \bar{\mathbb C},\qquad z\mapsto {z+1\over z},\quad T(0)=\infty,\quad T(\infty)=1\ .\tag{1}$$ Its fixed points are $\alpha:={1\over2}(1+\sqrt{5})$ and $\beta={1\over2}(1-\sqrt{5})$, obtained by solving the equation $z^2-z-1=0$.

We now introduce a new complex coordinate $w$ on $\bar {\mathbb C}$, related to $z$ via $$w=\phi(z):={z-\alpha\over z-\beta},\qquad{\rm resp.}\qquad z=\phi^{-1}(w):={\alpha-\beta w\over 1-w}\ .$$ The fixed points now are $w=0$ and $w=\infty$. In fact, in terms of the new coordiate $w$ the transformation $T$ appears as $\hat T=\phi\circ T\circ\phi^{-1}$, and computes to $$\hat T:\quad \bar{\mathbb C}\to \bar{\mathbb C},\qquad w\mapsto{\beta\over\alpha}w,\quad \hat T(0)=0,\quad \hat T(\infty)=\infty\ .$$ Since $${\beta\over\alpha}=-{3-\sqrt{5}\over2}\doteq-0.382$$ we can infer that the fixed point $0$ is attracting with basin of attraction all of ${\mathbb C}$, while $\infty$ is repelling. This allows to conclude that in the original setting all initial points $x_0\ne \beta$ lead to $\lim_{n\to\infty} x_n=\alpha$ (assuming the "exception handling" described in $(1)$).

You can proceed as follows. It will suffice to show that there is an index $r$ such that $x_r> 0$.

If $x_0\in A_0=(0,\infty)$, then $r=0$ and we are done.

If $x_0\in B_0=(-\infty, -1)$, then $r=1$ and we are done.

So we may assume without loss that $x_0\in (-1,0)$.

If $x_0\in A_1=\left(-\frac{1}{2},0\right)$, then $x_1\in B_0$, $r=2$ and we are done. So we may assume without loss that $x_0\in\left(-1,-\frac{1}{2}\right)$.

If $x_0\in B_1=\left(-1,-\frac{2}{3}\right)$, then $x_1 \in A_1$, $r=3$ and we are done. So we may assume without loss that $x_0\in\left(-\frac{2}{3},-\frac{1}{2}\right)$.

Continuing this way, we obtain (for $n\geq 1$) the two families $A_n=\left(-u_{n+1},-u_{n}\right), B_n=\left(\frac{-1}{u_{n}+1},\frac{-1}{u_{n+1}+1}\right)$ where $(u_n)$ is defined by $u_1=0$ and $u_{n+1}=\frac{u_n+1}{u_n+2}$. It is easy to check that $u_n$ stays in $(0,1)$, is increasing and converges to $\frac{-1+\sqrt{5}}{2}$.

Equally easily, we have $f(B_n)\subseteq A_n$ and $f(A_n)\subseteq B_{n-1}$, whence $r=2n$ whenever $x_0\in A_n$ and $r=2n+1$ whenever $x_0\in B_n$.

So the only case that's left is when $x_0$ is not in any of the $A_n$ or $B_n$. This means that $x_0$ is at one of the endpoints of $A_n,B_n (n\geq 1)$, i.e. $x_0$ is either one of your $a_k$'s or is $\frac 1 2 - \frac {\sqrt 5}2$.

Edit by the OP (ahorn):

We have seen from Ewan's answer that we can consider $x_0>0$ without loss of generality. What follows is an attempt to solidify the claim that $(x_n)$ converges in this case. Taking advice from Ewan, let $g=f\circ f$ where $f(x):=1+\frac 1 x$. That is, $g(x)=2-\frac 1{x+1}$ which is an increasing function.

Let $x_0\in(0,\phi)$, where $\phi=\frac 1 2 +\frac{\sqrt 5} 2$.

Since $x<g(x)<\phi=\sup\{g(x)\ |\ x\in(0,\phi) \}$ when $x\in(0,\phi)$, $$ x_{2n}<g(x_{2n})=x_{2n+2}<\phi=\sup\{x_{2k}\} $$ where $n, k\in \Bbb N$. So $(x_{2n})$ is an increasing sequence that converges to $\phi$.

$x_0\in(0,\phi)\implies x_{1}\in(\phi,\infty)$.

Since $x>g(x)>\phi=\inf\{g(x)\ |\ x\in(\phi,\infty) \}$ when $x\in(\phi, \infty)$, $$ x_{2n+1}>g(x_{2n+1})=x_{2n+3}>\phi=\inf\{x_{2k+1}\} $$ so $(x_{2n+1})$ is a decreasing sequence that converges to $\phi$.

Now, $x_n\in(0,\phi)\implies x_{n+1}\in(\phi,\infty)$ and $x_n\in(\phi, \infty)\implies x_{n+1}\in(0, \phi)$, so $x_0\in(0,\phi)$ was chosen without loss of generality.

Suppose that for any $\epsilon>0$, $2k\geq N_1\implies|x_{2k}-\phi|<\epsilon$ and $2k+1\geq N_2\implies|x_{2k+1}-\phi|<\epsilon$. Let $N=\max\{N_1, N_2\}$ so that $n\geq N \implies |x_n-\phi|<\epsilon.$