A method I've seen suggested (e.g. p 446 of this text) for estimating the Fisher information matrix (FIM) is by computing the sampling covariance matrix of the scores. That is, $$ \hat{\mathcal{I}}_n = \frac{1}{n} \sum_{i=1}^n (y_i - \bar{y}) (y_i - \bar{y})^T, $$ where $$ y_i = \nabla \log f(x_i; \hat{\theta}_n), \qquad i=1,\ldots,n $$ is the score function evaluated at the MLE estimates $\hat{\theta}_n$.
For normally-distributed data $X \sim \mathcal{N}(\mu, \sigma^2)$, the analytical FIM for the parameter vector $\theta := (\mu, \sigma)$ is given by $$ \mathcal{I}(\mu,\sigma) = \frac{1}{\sigma^2} \begin{pmatrix} 1 & 0 \\ 0 & 2 \end{pmatrix}. $$ However, when I estimate the FIM using the sample covariance technique for normally-distributed pseudorandom numbers, the error is quite large. E.g., for $n = 100$ and $X_i \sim \mathcal{N}(0,1)$, for one particular simulation I got $$ \hat{\mathcal{I}}_n = \begin{pmatrix} 0.9536 & 0.0332 \\ 0.0332 & 2.5354 \end{pmatrix}, \qquad \mathcal{I} = \begin{pmatrix} 1 & 0 \\ 0 & 2 \end{pmatrix}. $$ For $n = 1000$, for one particular simulation I got $$ \hat{\mathcal{I}}_n = \begin{pmatrix} 1.0089 & -0.1460 \\ -0.1460 & 2.1624 \end{pmatrix}, $$ where you can see the error in the covariances actually got much worse (this was just for one particular simulation, of course, so other times it was better).
My question is, is there a rule of thumb for using this technique, regarding the sample size $n$? I'm use the "n large enough" meaning $n \geq 30$ for CLT results, but it seems this doesn't apply here.
