Let $X_1,...,X_n$ be IID from a $GAMMA(\alpha, \beta)$ distribution, where $E(X_1)=\alpha\beta$. Derive the asymptotic distribution of the method of moment estimator $\tilde{\theta}$ of $\theta=(\alpha,\beta),$ that is:
$$\sqrt{n}(\tilde{\theta}-\theta)\to^dW$$
and give the expression of $W$.
In the above problem, both $\tilde{\theta}$ and $\theta$ should be bold to represent vectors. I can calculate the methods of moments estimators, easily; they are:
$$\tilde\alpha=\frac{\bar{x}^2}{\bar{x^2}-\bar{x}^2}$$
and
$$\tilde\beta=\frac{\bar{x^2}-\bar{x}^2}{\bar{x}}$$
where $\bar{x^2}$ represents $\sum x_i^2/n$ and $\bar{x}^2$ represents $(\sum x_i/n)^2$. However, I have no idea how to procde from here. How can I possibly find that limiting distribution with this information?
You can use multivariate Delta-method. Look here (paragraph 5.4) for similar calculations.
By multivariate CLT we can derive the asymptotic distribution of vector of sample moments: $$ \sqrt{n}\left[\pmatrix{\bar x\cr \bar{x^2}}-\pmatrix{\mathbb E[x_1]\cr\mathbb E[x_1^2]}\right]=\sqrt{n}\left[\pmatrix{\bar x\cr \bar{x^2}}-\pmatrix{\alpha\beta\cr\beta^2\alpha(\alpha+1)}\right]\xrightarrow{d} \mathcal N\left(\pmatrix{0\cr 0},\Sigma\right) $$ where $\Sigma$ is a covariance matrix of $\pmatrix{x_1\cr x_1^2}$. Please check all the calculations. $$ \Sigma =\pmatrix{\text{Var}(x_1) & \text{Cov}(x_1,x_1^2)\cr \text{Cov}(x_1,x_1^2) & \text{Var}(x_1^2)} =\pmatrix{\alpha\beta^2 & 2\alpha\beta^3(\alpha+1)\cr 2\alpha\beta^3(\alpha+1) & 2\alpha\beta^4(\alpha+1)(2\alpha+3)}. $$ Next, we have pair of estimates $$ \pmatrix{\tilde\alpha\cr\tilde\beta} = g\left(\bar x,\bar{x^2}\right), $$ where the vector-function $g:\mathbb R^2\to\mathbb R^2$ looks as $$ g(x,y)=\pmatrix{g_1(x,y)\cr g_2(x,y)} = \pmatrix{\frac{x^2}{y-x^2} \cr \frac{y-x^2}{x}}. $$ By multivariate Delta-method, $$ \sqrt{n}(\tilde{\theta}-\theta)=\sqrt{n}\left[g\left(\bar x,\bar{x^2}\right) - g\left(\mathbb E[x_1], \mathbb E[x_1^2]\right)\right] \xrightarrow{d} \mathcal N\left(\pmatrix{0\cr 0},\Sigma^*\right), $$ where $$ \Sigma^*=\left(\nabla_g\left(\mathbb E[x_1], \mathbb E[x_1^2]\right)\right)^T \cdot \Sigma \cdot \left(\nabla_g\left(\mathbb E[x_1], \mathbb E[x_1^2]\right)\right). $$ Here $$ \nabla_g\left(x, y\right) = \pmatrix{\frac{\partial g_1(x,y)}{\partial x} & \frac{\partial g_2(x,y)}{\partial x}\cr \frac{\partial g_1(x,y)}{\partial y} & \frac{\partial g_2(x,y)}{\partial y}} = \pmatrix{\frac{2xy}{(y-x^2)^2}& -\frac{y}{x^2}-1\cr -\frac{x^2}{(y-x^2)^2} & \frac1x }. $$
Substitute expectations instead of variables and get $$ \nabla_g\left(\mathbb E[x_1], \mathbb E[x_1^2]\right) = \pmatrix{ \frac{2(1+\alpha)}{\beta} & -\frac{2\alpha+1}{\alpha} \cr -\frac1{\beta^2} & \frac1{\alpha\beta}} $$
And finally get the covariance matrix of limiting two-dimensional normal distribution: $$ \Sigma^* = \pmatrix{2\alpha(\alpha+1) & -2\beta(\alpha+1) \cr -2\beta(\alpha+1) & \frac{\beta^2(2\alpha+3)}{\alpha}}. $$
At the first sight, this is a valid answer: it is a valid covariance matrix since it is symmetric and positively definite, and the variance of $\tilde \alpha$ does not depend on scaling parameter $\beta$. Note that the distribution of $\tilde \alpha$ does not depend on $\beta$ too since $\frac{x_1}{\beta}$ follows $GAMMA(\alpha,1)$ distribution and $\tilde\alpha$ does not change if we replace $x_i$ with $\frac{x_i}{\beta}$.