Consistency of the maximum likelihood estimator of $\rho$

Suppose $(X_i,Y_i)$, $i=1,2,...,n$ are a random sample from a bivariate Normal distribution with parameters $(\mu,-\mu,\sigma^2,\sigma^2,\rho)$. Obtain the Maximum Likelihood Estimator of $(\mu,\sigma^2,\rho)$ and check if the MLE of $\rho$ is consistent i.e. if $r_n$ be the MLE of $\rho$ based on $n$ observations, then check if $r_n\to^P\rho$

The loglikelihood is $$l=const.-n\log(\sigma^2)-\dfrac{n}{2}\log(1-\rho^2)-\dfrac{1}{2\sigma^2}[\sum_{i=1}^n (x_i-\mu)^2+\sum_{i=1}^n (y_i+\mu)^2-2\rho\sum_{i=1}^n (x_i-\mu)(y_i+\mu)]$$

The MLEs of $(\mu,\sigma^2,\rho)$ satisfy:

$\sigma^2=\dfrac{1}{2n}[\sum_{i=1}^n (x_i-\mu)^2+\sum_{i=1}^n (y_i+\mu)^2-2\rho\sum_{i=1}^n (x_i-\mu)(y_i+\mu)]$

$\sum_{i=1}^n -(x_i-\mu)+\sum_{i=1}^n(y_i+\mu)-2\rho \sum_{i=1}^n (x_i-y_i-2\mu)=0$

$\dfrac{n\rho}{1-\rho^2}+\sum_{i=1}^n \dfrac{(x_i-\mu)(y_i+\mu)}{\sigma^2}=0$

Solving these three equations seems nightmarish! And also, how can one check if MLE of $\rho$ is consistent?

Here is a naive attempt: From the third equation, we see that if we had MLE of $\mu$ and $\sigma^2$ consistent for the corresponding parameters then if MLE of $\rho$ was consistent for $\rho$, we would get, as $n\to\infty$, $\dfrac{\rho}{1-\rho^2}+\rho=0$, clearly not true for every possible value of $\rho$.

This is just a guess. As one can see, I am assuming consistency of MLEs of $\mu$ and $\sigma^2$. Can anyone give a proper solution?

You could denote $Z_i=X_i-Y_i\sim\mathcal N(2\mu,2\sigma^2(1-\rho))$, $W_i=X_i+Y_i\sim\mathcal N(0,2\sigma^2(1+\rho))$, and (it's simplifies things a lot) $Z_i$ and $W_i$ are uncorrelated. After that you may express the log-likelihood function via $Z$, $W$: $$\log L=-n\log(2\pi)-n\log2-n\log(\sigma^2)-\frac n2\log(1-\rho^2)-\frac1{4\sigma^2(1-\rho)}\sum_{i=1}^n(Z_i-2\mu)^2-\frac1{4\sigma^2(1+\rho)}\sum_{i=1}^nW_i^2$$ and get the first order conditions: \begin{aligned} \frac{\partial\log L}{\partial\mu}&=\frac1{\sigma^2(1-\rho)}\sum_{i=1}^n(Z_i-2\mu)=0\\ \frac{\partial\log L}{\partial(\sigma^2)}&=-\frac n{\sigma^2}+\frac1{4(\sigma^2)^2(1-\rho)}\sum_{i=1}^n(Z_i-2\mu)^2+\frac1{4(\sigma^2)^2(1+\rho)}\sum_{i=1}^nW_i^2=0\\ \frac{\partial\log L}{\partial\rho}&=n\frac\rho{1-\rho^2}-\frac1{4\sigma^2(1-\rho)^2}\sum_{i=1}^n(Z_i-2\mu)^2+\frac1{4\sigma^2(1+\rho)^2}\sum_{i=1}^nW_i^2=0 \end{aligned} This system simplifies to the following: \begin{aligned} \mu&=\overline Z/2\\ 2\sigma^2(1+\rho)&=s^2_W\\ 2\sigma^2(1-\rho)&=s^2_Z \end{aligned} where $s^2_Z=\frac1n\sum_{i=1}^n(Z_i-\overline Z)^2$, $s^2_W=\frac1n\sum_{i=1}^nW_i^2$. So, ultimately $$\hat\mu=\overline Z/2,\quad\hat\sigma^2=\frac{s^2_W+s^2_Z}4,\quad\hat\rho=\frac{s^2_W-s^2_Z}{s^2_W+s^2_Z}$$ And since the $s^2_W$ is a consistent estimator for $\mathop{\mathrm{Var}}W=2\sigma^2(1+\rho)$, and $s^2_Z$ is a consistent estimator for $\mathop{\mathrm{Var}}Z=2\sigma^2(1-\rho)$ (both are the sample variances, the difference is that we already know the expected value of $W$), then using the Slutsky's theorem, we can conclude that $\hat\rho$ is a consistent estimator of $$\frac{2\sigma^2(1+\rho)-2\sigma^2(1-\rho)}{2\sigma^2(1+\rho)+2\sigma^2(1-\rho)}=\rho.$$