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?
var(rnorm(10000,100,35000))(less severe, of course). $\endgroup$b^2/((a-1)^2*(a-2))is actually the theoretical variance. The method of moments consists of solving the system $\dfrac{\beta}{\alpha-1}=\bar{x}$ and $\dfrac{\beta^2}{(\alpha-1)^2(\alpha-2)}=\dfrac{1}{n}\sum_{j=1}^n(x_j-\bar{x})^2$ in terms of $(\alpha,\beta)$. After this, you could plug the solutions $(\hat{\alpha},\hat{\beta})$ into the expression for the variance $\dfrac{\beta^2}{(\alpha-1)^2(\alpha-2)}$ in order to obtain an estimator of this quantity. $\endgroup$