Why does quantum mechanics involve complex numbers?

First, I am using this helpful expository paper by Ricardo Karam. He gives four reasons, labeled (A) through (D). I want to discuss his third reason, (C), based on the Stern-Gerlach experiment, and my understanding of it.

First, there is an experiment which does the following; I'll be as vague as possible to avoid getting details wrong. A particle is sent through various boxes which observe its spin in the $x$, $y$ and $z$ directions, with the "equal filtering probability" property, i.e. that if a pure state for one observable is then observed by any other, it has equal probability of $1/2$ of being positive or negative spin.

Let's introduce some notation. First, the eigenvectors for the $z$-observation $S_z$ are $\lvert + \rangle$ and $\lvert - \rangle$, and we take this to be our standard basis. After performing a $x$-observation $S_x$ (or a $y$-observation $S_y$), we obtain the states $\lvert + \rangle_x$ and $\lvert - \rangle_x$ (respectively $\lvert + \rangle_y, \lvert - \rangle_y$). The matrices for the observables in this standard basis are
$A_z = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}, \;\;\; A_x = \begin{pmatrix} \vert & \vert \\ \lvert + \rangle_x & \lvert - \rangle_x \\ \vert & \vert \end{pmatrix}, \;\;\; A_y = \begin{pmatrix} \vert & \vert \\ \lvert + \rangle_y & \lvert - \rangle_y \\ \vert & \vert \end{pmatrix}$.
For convenience, we record that a solution (up to various symmetries) for the problem is: $A_z = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}, \;\;\; A_x = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}, \;\;\; A_y = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ i & -i \end{pmatrix}$. We observe that these matrices are unitary.

Because I'm not too familiar with bra-ket notation, I need to make a digression. Vectors of the form $\langle \bullet \rvert$ are covectors, e.g. ${}_x\langle + \rvert$ is the first column of the conjugate transpose $A_x^*$. What is the meaning of these covectors? By definition, to pair with $\langle + \vert$ and $\langle - \vert$ gives the proportion of particles which have positive or negative $z$-spin. We generalize: ${}_y\langle +\vert - \rangle_x$ is a number whose square magnitude is the probability that a particle which has been observed to have negative $x$-spin will then be observed to have positive $y$-spin. Note that we compute the quantity ${}_y\langle +\vert - \rangle_x$ by taking the corresponding entry of the matrix $A_y^* A_x$.

Why square magnitude and not magnitude? All our matrices will have the property that $A^* A = I$, i.e. they are unitary, because repeating an observation should result in the same thing. The columns of unitary matrices are have norm one, which satisfy the property that the sum of the squares of the magnitudes equals one. Thus the correct quantity to interpret as a probability (i.e. the thing that sums to one) are the squares.

Anyway, the point is, the "equal probability" property means that any pair of eigenvectors (from different observables) should pair to value $1/2$ (and the fact that repeated measurements give the same value means for a single observable the eigenvalues must be orthogonal). If they were all real eigenvectors, this means that their squared dot product is $1/2$, i.e. the angle between the vectors is $45^\circ$. But this is impossible to do with three pairs of eigenvectors in $\mathbb{R}^2$, so the eigenvectors must be complex. Moreover, the matrices have complex entries with respect to the standard basis which we have fixed.

Finally, it's worth constrasting this situation with polarizing filters; I saw an argument on Twitter at some point about whether polarization is a quantum phenomenon or not. I'm a bit confused by this. On the one hand, one can clearly model light polarization in terms of projections, but it appears there is both a quantum and a classical model for polarization of photons. I want to say this is due to the fact that photon polarization only involves two observables, but I haven't really checked this.