Suppose I have an estimator $B\in\mathbb{R}^m$ converging to $\beta$, such that $$ \sqrt{n}(B-\beta)\rightarrow\mathcal{N}(0,\Sigma). $$
I am interested in a quantity $\mathbf{h}(B):\ \mathbb{R}^m\rightarrow\mathbb{R}^p$, and would like to use delta method to approximate the asymptotic distribution.
Let $\nabla\mathbf{h}(\beta)$ be a matrix $\left[\frac{\partial h_k(\beta)}{\partial \beta_j}\right]_{jk}$. When $p=1$, multivariate Delta method would provide $$ \sqrt{n}(\mathbf{h}(B)-\mathbf{h}(\beta))\rightarrow\mathcal{N}(0,\nabla\mathbf{h}(\beta)^T\cdot\Sigma\cdot\nabla\mathbf{h}(\beta)). $$ Can we say that it is the same for the case when $p>1$? When I actually tried this numerically, I find that the covariance $\nabla\mathbf{h}(\beta)^T\cdot\Sigma\cdot\nabla\mathbf{h}(\beta)$ is not invertible.
