Asymptotic Distribution of Vector of Method of Moments Estimators? 
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?
 A: 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}$. 
