I've made this community wiki so as not to take any credit for @Procrastinator or anyone else's solution, but this question does have an answer.
The first element to consider is that if the true distribution's $\alpha$ is approaching 2 from above, the true variance approaches infinity. Therefore, even if the variance is finite, as $\alpha$ approaches 2 from above, the "true" variance grows similar to $\lim_{x\rightarrow0}\frac{1}{x}$ (hyperbolically?) and the number of simulations needed to capture that extraordinary growth may well be intractable within any feasible time scale (like a human life). The Central Limit Theorem only holds in cases of finite variance, and as the variance gets out of control, its power weakens. For all we know, the third moment is already infinite and the Berry-Esseen theorem no longer holds, and even though the variance is finite, the rate of convergence is completely unknown and almost certainly glacial. The 9,999 observations in the question simply weren't enough to capture the tail.
A second important point is that
$$
\sigma^2=\frac{\beta^2}{(\alpha−1)^2(\alpha−2)}
$$
is the true distribution given know $\alpha$ and $\beta$. When using method of moments, we are making the inference in the opposite direction. We have known empirical $\bar{x}$ and $\bar{s}^2$ from which we want to estimate an $\hat{\alpha}$ and $\hat{\beta}$. Running through the algebra, if I have not made an error, what is actually being solved is:
$$
\begin{aligned}
\hat{\alpha} &= \frac{\bar{x}^2}{\bar{s}^2}+2\\
&=\frac{1}{CV^2} + 2\\
\hat{\beta} &= \left(\frac{\bar{x}^2}{\bar{s}^2}+1\right)\bar{x}\\
&=\left(\frac{1}{CV^2} + 1\right)\bar{x}
\end{aligned}
$$
While the Law of Large Numbers says these estimates will asymptotically approach their true values, once again, the rate of convergence may be on the order of the heat-death of the universe, and this is more likely as the parameters approach their hard limits. For an inverse gamma with a true $\alpha \leq 2$ there is no finite number of simulations whose variance would be captured by MoM. Practically, as the function has a $+2$ in it, but more importantly, as the true variance is infinite, the sum of squared errors does not converge. Similarly, For an inverse gamma with a true $\alpha \leq 1$ there is no finite number of simulations whose mean would be captured by MoM. Practically, note the $+1$. More importantly, in this case, as more observations are generated, larger and larger values are seen and the mean diverges.
And why are inverse gamma distributions used so often as the prior for variances?
Because they're conjugate - at least for $\sigma^2$ in normal models.
How should it work?
You seem to have some confusion between Jeffreys' prior and (more-or-less) "uninformative" conjugate priors more generally. If you are taking a Jeffreys' prior, there's no choice involved: $p(\theta) \propto \sqrt{I(\theta)}\,$
where $I(.)$ is the Fisher information.
http://en.wikipedia.org/wiki/Jeffreys_prior
Best Answer
When estimating the variance of normal distribution with known mean $\mu$ using conjugate inverse gamma prior, the model is
$$ X \sim \mathrm{Normal}(\mu, \sigma^2) \\ \sigma^2 \sim \mathrm{IG}(\alpha, \beta) $$
So we assume an inverse gamma prior parametrized by $\alpha$ and $\beta$ for $\sigma^2$. When we observe some data $x_1,\dots,x_n$ we can update our prior, to get posterior distribution of $\sigma^2$ and since inverse gamma is a conjugate prior, the update has a simple closed-form solution. The updated parameters become
$$ \alpha' = \alpha + \frac{n}{2}, \qquad \beta' = \frac{\sum_{i=1}^n (x_i - \mu)^2}{2} $$
and the posterior distribution of $\sigma^2$ is inverse gamma parametrized by $\alpha'$ and $\beta'$. So by updating the $\alpha$ and $\beta$ parameters you already estimated the distribution of $\sigma^2$. If you want to obtain a point estimate of $\sigma^2$, then you can calculate mean $\tfrac{\beta'}{\alpha'-1}$ or mode $\tfrac{\beta'}{\alpha'+1}$ of it's posterior distribution. What you were trying to do, is you were trying to calculate the variance of $\sigma^2$, that is certainly not it's point estimate.
Check also Bayesian updating with conjugate priors using the closed form expressions and Bayesian updating with new data for similar questions.