You must put:
$$\sigma^2=\frac{1}{N-1}\sum_{i=1}^N(y_i-\hat{y}_i)\overline{(y_i-\hat{y}_i)}$$
that is, replacing your code by:
for i = 1: N
dy(i,1) = (y(i)-yhat(i))*conj(y(i)-yhat(i));
end
sigma2 = 1/(N-1)*sum(dy);
The variance is the second moment of a probability distribution, which can be defined for any metric space in the sense of the definition over the random variable (integral of a density), or in any vector space, in the sense of the estimator of the random variable (sum of n samples).
For $X$ a random variable $\in \mathbb{R}^n$, with probability density function $p(x):I_X \subset \mathbb{R}^n \rightarrow \mathbb{R}$, the second moment $\in \mathbb{R}$ is:
$$\sigma^2=\mathbf{E}\{(X-\mu)^2\}=\int_{I_X} (x-\mu)^2 p(x)dx$$
which is the variance for $n=1$ and the covariance matrix in general.
For other fields $\mathbb{F}$, complex and even real numbers, the metric $d(x,y): \mathbb{F}x\mathbb{F}\rightarrow \mathbb{R}$ can redefine the moment definition
$$\sigma^2_{\mathbb{F}}=\mathbf{E}_{\mathbb{F}}\{(X-\mu)^2\}=\int_{I_X} d(x,\mu)^2 p(x)dx$$
For example, for the complex field, the metric function can be the standard distance from the hermitian inner product, and thus we have the variance definition:
$$d(x,y)=\sqrt{(x-y)\overline{(x-y)}}$$
$$\sigma^2_{\mathbb{C}}=\mathbf{E}_{\mathbb{C}}\{(X-\mu)^2\}=\int_{I_X} (x-y)\overline{(x-y)} p(x)dx$$
But the standard product, without the conjugated term, though can be a valuable figure (named as the "complementary variance") is not a second moment:
$${\sigma^2}_{\mathbb{C}comp}=\int_{I_X} (x-y)^2 p(x)dx$$
Of course, you always could be more "reasonable" and prefer to drop the complex analysis and move onto a $\mathbb{C}^2$ cartesian plane, and start taking covariances instead of a single number which dont dive you the full distribution information of your stochastics.