I'm trying to calculate the variance of the inverse gamma distribution using the method of movements. According to wikipedia the variance should be:
$$\sigma^2 =\frac{\beta^2}{(\alpha-1)^2(\alpha-2)}$$
Where $\alpha$ is the shape and $\beta$ is the scale of the inverse gamma distribution. When trying this out in R it works reasonably well when $\alpha$ is not to close to 2. For example:
library(MCMCpack) # for the rinvgamma function
a <- 10
b 100
# The variance according to the method of movements
b^2/((a-1)^2*(a-2))
## 15.4321
# The variance by generating inverse gamma distributed random numbers and
# calculating the sample variance
var(rinvgamma(n=9999, shape=a, scale=b)) #
## 15.84388
But when $\alpha$ gets close to 2 the method of movements doesn't seem to work anymore. In the following example the sample variance is much smaller than the method of movements variance:
a <- 2.2
b <- 100
# The variance according to the method of movements
b^2/((a-1)^2*(a-2))
## 34722.22
# The variance by generating inverse gamma distributed random numbers and
# calculating the sample variance
var(rinvgamma(n=9999, shape=a, scale=b)) #
##14479.56
Why doesn't the method of movements work? Am I doing something wrong that can be fixed or is there some other way that I can calculate the variance of an inverse gamma distribution?
Best Answer
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.