Here is another description, using only a single q-bit, that shows how the imaginary part of $|\psi\rangle$ plays a central role in capturing the state and its dynamics.
A spin-$\frac{1}{2}$ system is described by a two-state Hilbert space, so a general state takes the form
$$|\psi\rangle=\left[\begin{array}{c}
\alpha \\
\beta
\end{array}
\right],$$
with $|\alpha|^{2}+|\beta|^{2}=1$, so that the state is normalized. With two complex numbers $\alpha$ and $\beta$, subject to one equation $|\alpha|^{2}+|\beta|^{2}=1$, there are three (real) parameters in this $|\psi\rangle$. However, one of those parameters is superfluous: the overall phase. If we define a new spinor $|\psi'\rangle$ equal to $e^{i\chi}|\psi\rangle$, then $|\psi\rangle$ and $|\psi'\rangle$ describe the same physical state. This means, for instance, that any expectation value $\langle\psi'|{\cal O}|\psi'\rangle$ remains equal to $\langle\psi|{\cal O}|\psi\rangle$, since
$$\langle\psi'|{\cal O}|\psi'\rangle=\langle\psi|e^{-i\chi}{\cal O}e^{i\chi}|\psi'\rangle=\langle\psi|{\cal O}|\psi\rangle,$$
since $e^{i\chi}$ is just a number and so commutes with the operator ${\cal O}$.
So the overall phase of this (or any other wave function) does not contain any physical information. That leaves two real degrees of freedom that parameterize the physical state $|\psi\rangle$; these may be encapsulated in the single complex number $\alpha/\beta$ which can take any complex value (or be infinite, if $\beta=0$). The quantum-mechanical spin-$\frac{1}{2}$ can describe the same dynamics as a classical vector $\vec{S}$ of fixed length, with the two parameters specifying the state mappable onto the spherical angles $\theta$ and $\phi$ describing the direction of $\vec{S}$.
The connection is made through taking the expectation values of the three Cartesian components of the spin operator,
$$\left\langle\vec{S}\right\rangle=\frac{\hbar}{2}\left[\langle\sigma_{1}\rangle\hat{\imath}+
\langle\sigma_{2}\rangle\hat{\jmath}+\langle\sigma_{3}\rangle\hat{k}\right].$$
The expectation values of the Pauli matrices are
$$\begin{array}{rcl}
\langle\sigma_{1}\rangle=\langle\psi|\sigma_{1}|\psi\rangle & = &
[\alpha^{*}\,\beta^{*}]
\left[\begin{array}{cc}
0 & 1 \\
1 & 0
\end{array}
\right]
\left[\begin{array}{c}
\alpha \\
\beta
\end{array}
\right]=\beta^{*}\alpha+\alpha^{*}\beta=2\,\Re\{\alpha^{*}\beta\}\\
\langle\sigma_{2}\rangle=\langle\psi|\sigma_{2}|\psi\rangle & = &
[\alpha^{*}\,\beta^{*}]
\left[\begin{array}{cc}
0 & -i \\
i & 0
\end{array}
\right]
\left[\begin{array}{c}
\alpha \\
\beta
\end{array}
\right]=i\beta^{*}\alpha-i\alpha^{*}\beta=2\,\Im\{\alpha^{*}\beta\}\\
\langle\sigma_{3}\rangle=\langle\psi|\sigma_{3}|\psi\rangle & = &
[\alpha^{*}\,\beta^{*}]
\left[\begin{array}{cc}
1 & 0 \\
0 & -1
\end{array}
\right]
\left[\begin{array}{c}
\alpha \\
\beta
\end{array}
\right]=
|\alpha|^{2}-|\beta|^{2}
\end{array}.$$
Note that the square of the length of the expectation value, $\langle\sigma_{1}\rangle^{2}+\langle\sigma_{2}\rangle^{2}+\langle\sigma_{3}\rangle^{2}$ (which is not the same as the expectation value of the square, $\langle\sigma_{1}^{2}\rangle+\langle\sigma_{2}^{2}\rangle+\langle\sigma_{3}^{2}
\rangle$) is normalized so that
$$\sqrt{\left\langle\vec{S}\right\rangle^{2}}=\frac{\hbar}{2}.$$
However, to have a fully accurate representation of the full space of possible directions evidently requires the use of complex numbers; if $\alpha$ and $\beta$ are both strictly real, then we will inevitably have $\left\langle S_{y}\right\rangle=0$.
It is straightforward to verify that if the direction of the vector $\left\langle\vec{S}\right\rangle$ is described by the polar angle $\theta$ and azimuthal angle $\phi$, meaning
$$\left\langle\vec{S}\right\rangle=\hat{n}=\frac{\hbar}{2}\left(\sin\theta\cos\phi\,\hat{\imath}+\sin\theta\sin\phi\,\hat{\jmath}+
\cos\theta\,\hat{k}\right)=\frac{\hbar}{2}\hat{n}$$
(in terms of a unit vector $\hat{n}$), then $|\psi\rangle$ takes the form
$$|\psi\rangle=\left[\begin{array}{c}
\cos(\theta/2) \\
e^{i\phi}\sin(\theta/2)
\end{array}
\right],$$
(up to a total phase, as always). Moreover, this spinor is also an eigenstate of the projection of the total spin $\vec{S}$ along the $\hat{n}$-direction:
$$\left(\vec{S}\cdot\hat{n}\right)|\psi\rangle=+\frac{\hbar}{2}|\psi\rangle.$$
Finally, when we introduce a time dependence—with, say, the Hamiltonian $H=-\mu B_{0}\,\sigma_{3}$ for a spin carrying magnetic moment $\mu$ in the presence of a magnetic field $\vec{B}=B_{0}\hat{k}$ along the $z$-direction—even an initially real spinor wave function will develop an imaginary part as it evolves in time. Under the action of
$$H=\mu B_{0}\left[\begin{array}{cc}
1 & 0 \\
0 & -1
\end{array}
\right],$$
an initial state
$$|\psi(0)\rangle=\left[\begin{array}{c}
\alpha \\
\beta
\end{array}
\right],$$
will evolve according to
$$|\psi(t)\rangle=\left[\begin{array}{c}
e^{-i\omega t/2}\alpha \\
e^{i\omega t/2}\beta
\end{array}\right].$$
The frequency here is $\omega=2\mu B_{0}/\hbar$. The presence of the imaginary phases is what causes the state to evolve in time. If the initial spinor is strictly real, then the expectation value of the spin will evolve in time according to
$$\left\langle\vec{S}\right\rangle(t)=\frac{\hbar}{2}\left(\sin\theta\cos\omega t\,
\hat{\imath}+\sin\theta\sin\omega t\,\hat{\jmath}+\cos\theta\,{\hat k}\right),$$
which is just the same as the precession of a classical magnetic moment $\vec{\mu}$ in an external magnetic field.
