A single-qubit state can be parametrized with real $\theta$ and $\phi$ as follows: $$|\psi(\theta, \phi)\rangle = \cos\frac{\theta}{2}|0\rangle + e^{i \phi} \sin\frac{\theta}{2}|1\rangle.$$
I would like to know how to parameterize an arbitrary two-qubit state in a similar way.
TL;DR: The logic behind the equation $|\psi\rangle=\cos\frac{\theta}{2}|0\rangle+e^{i\phi}\sin\frac{\theta}{2}|1\rangle$ can be iterated to obtain a real parameterization for states of arbitrary finite dimension.
We parameterize an expansion of $|\psi\rangle$ in basis $|0\rangle,\dots,|N\rangle$ in $N$ iterations. The $k$th iteration for $k=1,\dots,N$ consists of two steps. First, we introduce $\theta_k\in[0,\pi]$ to apportion the unit norm between the magnitudes of the amplitudes of $|k-1\rangle$ and the projection $|\psi_{k\dots N}\rangle$ of $|\psi\rangle$ onto the subspace spanned by $|k\rangle,\dots,|N\rangle$. Next, we stick a relative phase $\alpha_k\in[0,2\pi)$ on $|k\rangle$. Thus, the $k$th iteration rewrites $|\psi_{k-1\dots N}\rangle$ as $$ |\psi_{k-1\dots N}\rangle=\cos\frac{\theta_k}{2}|k-1\rangle + e^{i\alpha_k}\sin\frac{\theta_k}{2}|\psi_{k\dots N}\rangle.\tag1 $$ The final result can be simplified by setting $\phi_k:=\alpha_1+\dots+\alpha_k\mod 2\pi$.
Each step introduces a single real parameter, for a total of $2N$ real parameters, as expected for an $N+1$ dimensional system.
For a two-qubit state the above procedure yields $$ \begin{align} |\psi\rangle&=\cos\frac{\theta_1}{2}|00\rangle\\ &+ e^{i\phi_1}\sin\frac{\theta_1}{2}\cos\frac{\theta_2}{2}|01\rangle\\ &+ e^{i\phi_2}\sin\frac{\theta_1}{2}\sin\frac{\theta_2}{2}\cos\frac{\theta_3}{2}|10\rangle\\ &+ e^{i\phi_3}\sin\frac{\theta_1}{2}\sin\frac{\theta_2}{2}\sin\frac{\theta_3}{2}|11\rangle \end{align}\tag2 $$ where $\theta_1,\theta_2,\theta_3\in[0,\pi]$ and $\phi_1,\phi_2,\phi_3\in[0,2\pi)$.
This isn't the only way to parameterize a two-qubit state. An alternative approach uses Schmidt decomposition and single-qubit parameterization.