Let $X_1,\dots,X_n$ be iid $\Gamma(p,1/\lambda)$ with density $g_\theta (x) = \frac{1}{\Gamma(p)} \lambda^p x^{p-1} e^{-\lambda x}$, $x>0$, $\theta = (p,\lambda)$, $p > 0$, $\lambda > 0$. Show that the moment estimates, $\hat{\theta}$, are asymptotically bi-variate normal and give their asymptotic mean and variance-covariance matrix.
First off, I calculated that $\mathbb{E}X = \frac{p}{\lambda}$ and $\mathbb{E}X^2 = \frac{p(p+1)}{\lambda^2}$. So, setting $\bar{X} = \frac{p}{\lambda}$ and $\frac{1}{n}\sum X_i^2 = \frac{p(p+1)}{\lambda^2}$, I got the moment estimates $$\hat{\lambda} = \frac{\bar{X}}{\frac{1}{n} \sum (X_i - \bar{X})^2},$$ and $$\hat{p} = \frac{\bar{X}^2}{\frac{1}{n} \sum (X_i - \bar{X})^2}.$$ Now... I have no idea where to go from here.