Given a dataset with two variables $X$ and $Y$, with each observation independent of the others, test the null hypothesis $$\mu_X = \mu_Y\\ \sigma^2_X = \sigma^2_Y$$ using a Wald test.
This question comes from a textbook, but my answer isn't going to be handed in. Hence, I'm posting this question here because I'd like my work to be checked. Our class was given a dataset to use for this question, and I use Stata. I'm not sure of the best way to solve this question—the main struggle I had was with figuring out how to find the variance-covariance matrix for the Wald statistic, as I didn't think it was possible to create a MLE random variable for $\widehat{\mu_X} - \widehat{\mu_Y}$ and $\widehat{\sigma_X^2} - \widehat{\sigma^2_Y}$, in order to find their variances and covariances.
To implement the Wald test, I construct the vector of restrictions such that it is predicted to equal $\boldsymbol 0$. In this case, that is $$\left[\begin{matrix} {{\mu_X - \mu_Y}} \\ \sigma^2_X - \sigma^2_Y \end{matrix}\right] = \left[\begin{matrix} 0 \\ 0 \end{matrix}\right] $$ For notational convenience, label this $\boldsymbol v$. Then $\boldsymbol{{\hat v'}}\left(\Sigma_{\boldsymbol{{\hat v}}}\right)^{-1}\boldsymbol{{\hat v}}$ is the Wald Statistic. The covariance matrix is \begin{align*} \Sigma_{\boldsymbol{{\hat v}}} &= (\boldsymbol{{\hat v}} - \mathrm E[\boldsymbol{{\hat v}}])(\boldsymbol{{\hat v}} - \mathrm E[\boldsymbol{{\hat v}}])' \\ &=\left[\begin{matrix} \mathrm{Var}(\widehat{\mu_X} - \widehat{\mu_Y}) & \mathrm{Cov}\left(\widehat{\mu_X} - \widehat{\mu_Y}, \widehat{\sigma^2_X} - \widehat{\sigma^2_Y}\right) \\ \mathrm{Cov}\left(\widehat{\mu_X} - \widehat{\mu_Y}, \widehat{\sigma^2_X} - \widehat{\sigma^2_Y}\right) & \mathrm{Var}\left(\widehat{\sigma^2_X} - \widehat{\sigma^2_Y}\right) \end{matrix}\right] \end{align*}
Let's look at each element of the matrix in turn: \begin{align*} &&\mathrm{Var}\left(\widehat{\mu_X} - \widehat{\mu_Y}\right) &= \mathrm{Var}\left(\frac 1n\sum_{i=1}^n x_i - \frac 1n\sum_{i=1}^n y_i\right) \\ && &= \mathrm{Var}\left(\frac 1n\left(\sum_{i=1}^n(x_i - y_i)\right)\right) \\ && &=\frac1{n^2}\left(\sum_{i=1}^n \mathrm{Var}(x_i) + \sum_{i=1}^n\mathrm{Var}(y_i) - \sum_{i=1}^n 2\mathrm{Cov}(x_i, y_i)\right) \\ && &= \frac1n\left(\sigma^2_X + \sigma^2_Y - 2\mathrm{Cov}(X, Y)\right) \end{align*} where the observations (individuals) are independent of each other, so $\mathrm{Cov}(x_i, y_j) = 0 $ for $i\neq j$. This scalar will be computed in Stata. Easy enough, but trying to use this same method of deriving formulae for the other entries in terms of $\mu_X, \mu_Y, \sigma_X^2$ and $\sigma_Y^2$ seems to be difficult. (In order to do that, I would use the expected value formula for a variance $\mathrm{Var}(W) = \mathrm E[W^2] - \mathrm E[W]^2$, and substitute $x_i = \mu_X + \varepsilon_{X,i}$ (also with $y_i$) to use that fact that $\mathrm E[\varepsilon_{X,i}^2] = \sigma_X^2$ (also with $\varepsilon_{Y, i}^2$).)
I chose a simpler method of generating variables for $(X-\mu_X)^2$ and $(Y-\mu_Y)^2$ in Stata, in order to compute the covariances according to the formulae written below.
The second and third entries are $$% \mathrm{Cov}\left(\widehat{\mu_X} - \widehat{\mu_Y}, \widehat{\sigma^2_X} - \widehat{\sigma^2_Y}\right) = \mathrm{Cov}\left(\widehat{\mu_X}, \widehat{\sigma^2_X}\right) - \mathrm{Cov}\left(\widehat{\mu_Y}, \widehat{\sigma^2_X}\right) -\mathrm{Cov}\left(\widehat{\mu_X}, \widehat{\sigma^2_Y}\right) + \mathrm{Cov}\left(\widehat{\mu_Y}, \widehat{\sigma^2_Y}\right)% $$ Where \begin{align*} \mathrm{Cov}\left(\widehat{\mu_X}, \widehat{\sigma^2_X}\right) &= \mathrm{Cov}\left(\frac 1n\sum_{i=1}^n x_i, \frac 1n\sum_{i=1}^n (x_i - \mu)^2\right) \\ &= \frac 1{n^2}\sum_{i=1}^n\sum_{j=1}^n\mathrm{Cov}\left(x_i, (x_j - \mu_X)^2\right)\\ &\text{where $x_i$ and $(x_j - \mu)^2$ are not correlated when $i\neq j$} \\ &= \frac 1n\mathrm{Cov}\left(X, (X - \mu_X)^2\right) \end{align*} Similarly, \begin{multline*} \mathrm{Cov}\left(\widehat{\mu_X} - \widehat{\mu_Y}, \widehat{\sigma^2_X} - \widehat{\sigma^2_Y}\right) \\ =\frac 1n\bigg(\mathrm{Cov}\left(X, (X - \mu_X)^2\right) - \mathrm{Cov}\left(Y, (X-\mu_X)^2\right) - \mathrm{Cov}\left(X, (Y-\mu_Y)^2\right) + \mathrm{Cov}\left(Y, (Y - \mu_Y)^2\right)\bigg) \end{multline*}
Lastly, $\mathrm{Var}\left(\widehat{\sigma^2_X} - \widehat{\sigma^2_Y}\right)$ \begin{align*} &= \mathrm{Var}\left(\frac 1n\left(\sum_{i=1}^n(x_i - \mu_X)^2 - \sum_{i=1}^n(y_i - \mu_Y)^2\right)\right) \\ &= \frac 1{n^2}\left(\sum_{i=1}^n\mathrm{Var}\left((x_i - \mu_X)^2\right) + \sum_{i=1}^n\mathrm{Var}\left((y_i - \mu_Y)^2\right) - \sum_{i=1}^n\mathrm{Cov}\Big((x_i - \mu_X)^2, (y_i - \mu_Y)^2\Big)\right) \\ &= \frac 1n\bigg(\mathrm{Var}\left((X - \mu_X)^2\right) + \mathrm{Var}\left((Y - \mu_Y)^2\right) - \mathrm{Cov}\Big((X - \mu_X)^2, (Y - \mu_Y)^2\Big)\bigg) \end{align*}
For each of these formulae for the entries of the matrix, substituting $\mu_X = \widehat{\mu_X},\ \mu_Y = \widehat{\mu_Y},\ \sigma_X^2 = \widehat{\sigma^2_X},$ and $\sigma_Y^2 = \widehat{\sigma^2_Y}$ is asymptotically valid. Hence, my Stata code is this.
I thought of using the Fisher Information matrix $I(\boldsymbol{{\hat v}}) \approx \Sigma_{\boldsymbol{{\hat v}}}^{-1}$ according to the Cramer-Rao Lower Bound, but I don't know how to differentiate the bivariate loglikelihood w.r.t. an expression (i.e. $\mu_X - \mu_Y$). I have found the bivariate loglikelihood function under the assumption of normality.
