If we denote:
$a_k=E_{\theta}(X^k)$, $a_{nk}=\frac{1}{n}\sum_{i=1}^n X_{i}^k$. Further assume $\frac{1}{n}\sum_{i=1}^n X_{i}^{2k}<\infty$.
$X_1,X_2,\cdots,X_n\;(iid)\sim f(x;\theta),\; \theta \in \Theta$.
$g(\theta)$ is a real valued function defined on $\Theta$ and can be rewritten as $g(\theta)=G(a_1,a_2,\cdots,a_k)$ where $G$ is a function that ensures $\frac{\partial G}{\partial a_i}$ exist and continuous.
$\hat{g}(\vec{X})=\hat{g}(X_1,X_2,\cdots,X_n)=G(a_{n1},a_{n2},\cdots,a_{nk})$ is the moment estimator of $\hat{g}(\theta)$
Then
$\hat{g}(\vec{X})$ is a consistent asymptotically normal (CAN) estimator for $g(\theta)$, i.e.
$$ \begin{cases} \hat{g}(\vec{X}) \text{ is a weak consistent estimator for } g(\theta) \qquad (1)\\ \sqrt{n}\cdot\left( \hat{g}(\vec{X})-G(a_1,a_2,\cdots,a_k) \right) \rightarrow^{d} N(0,b^2)\quad (2) \end{cases} $$
where in (2) $b^2=\mathbf{d}^{T}\mathbf{B}\mathbf{d}$, here $\mathbf{d}$ is a vector where $d_i=\frac{\partial G(a_1,a_2,\cdots, a_k)}{a_i}$, $\mathbf{B}$ is a matrix where $B_{ij}=a_{i+j}-a_{i}a_{j}$. And $\rightarrow^{d}$ means convergence in distribution.
Right now I know how to prove (1) since moment estimators can easily be shown as strong consistent estimators. However, the proof of (2) troubles me and I don't have idea for it (especially those matrix forms).
Should this part be proven with Delta Method? How to?
Really appreciate it if you can stop and take a glance of my question.