I will focus on the reason of the factor $1/2$ and leave aside the estimation things.
The exact understanding should be : if a scalar Gaussian random variable (rv) is circular symmetric, its real and imaginary parts must be uncorrelated (this is equivalent to independence if they are assumed jointly Gaussian) and identically distributed with zero mean. Thus, your Matlab code is correct for a rv $\sim \mathcal{CN}(0,1)$.
The story behind is that a complex random variable (rv) is simply a vector of two real random variables. A vector of $n$ complex rv is indeed a vector of $2n$ real rv.
You are talking about the case $n=1$, 1 complex Gaussian rv $Z = Z_r + jZ_i$ or a vector of 2 real Gaussian rv $[Z_r,Z_i]^T$. As with real Gaussian rv which is described by its variance, a real vector $[Z_r,Z_i]^T$ must be described by its covariance matrix
$$\mathbb{E}\left\lbrace [Z_r,Z_i]^T \times [Z_r,Z_i] \right\rbrace = \mathbb{E}\left\lbrace\begin{bmatrix}
Z_r^2 & Z_rZ_i \\
Z_iZ_r & Z_i^2 \\
\end{bmatrix}\right\rbrace$$
Take a look at the variance $\mathbb{E}\left\lbrace ZZ^H \right\rbrace = \mathbb{E}\left\lbrace Z_r^2+Z_i^2 \right\rbrace$ and pseudo-variance $\mathbb{E}\left\lbrace ZZ^T \right\rbrace = \mathbb{E}\left\lbrace Z_r^2-Z_i^2+j2Z_rZ_i \right\rbrace$, the covariance matrix of real vector (or complex scalar) is fully described by the variance and pseudo-variance. Threfore, you need both variance and pseudo-variance to characterize a complex Gaussian rv (with the prior condition that its real and imaginary parts are jointly Gaussian).
Now we use the circular symmetry property : $e^{j\phi}Z$ has the same probability distribution as $Z$ for all real $\phi$. This leads to $\mathbb{E}\left\lbrace e^{j\phi}Z(e^{j\phi}Z)^T \right\rbrace = \mathbb{E}\left\lbrace e^{j2\phi}ZZ^T \right\rbrace = \mathbb{E}\left\lbrace ZZ^T \right\rbrace$ for all $\phi$ thus $E[ZZ^T] = 0$ and variance $\mathbb{E}\left\lbrace ZZ^H \right\rbrace$ is sufficient statistic for $Z$. Note that $\mathbb{E}\left\lbrace ZZ^T \right\rbrace = \mathbb{E}\left\lbrace Z_r^2-Z_i^2+j2Z_rZ_i \right\rbrace \implies \mathbb{E}\left\lbrace Z_r^2 \right\rbrace = \mathbb{E}\left\lbrace Z_i^2 \right\rbrace$ and $\mathbb{E}\left\lbrace Z_rZ_i \right\rbrace=0$, the real and imaginary parts are uncorrelated then independent (because they are Gaussian), with the same variance, this is the reason of the factor $1/2$.
To sum up, your code is correct because you are estimating circular symmetric complex Gaussian rv. The jointly Gaussian assumption between real and imaginary parts must be used. If this is not about circular symmetric rv (or real random vector with two elements), you must calculate the pseudo-variance.
For more details and to understand the formula of the wikipedia article, you can read this article Circular Symmetric Gaussian R.Gallager.
