Do the moments of the reciprocal normal distribution exist?

Question

The following question is based on this question: Reciprocal of a normal variable with non-zero mean and small variance

To summarize the main information from that question:

$X$ is a normally distributed random variable:

$$X \sim \mathcal{N}(\mu,\sigma^2)$$

Then $Y = 1/X$ has the following probability density function (see wiki):

$$f(y) = \frac{1}{y^2\sqrt{2\sigma^2\pi}}\,\exp\left(-\frac{(\frac{1}{y} - \mu)^2}{2 \sigma^2}\right)$$

This distribution of $Y$ does not have moments since (stackExchange):

$$\int_{-\infty}^{+\infty}|x|f(x)\,dx = \infty$$

An intuitive explanation of this is that the distributions tails are too heavy and consequently the law of large number fails. The more samples that are drawn and averaged the less stable this average is. The non-zero probability density that $X = 0$ means that $Y$ will not have finite moments since there is a non-zero probability that $Y = \infty$.

However, I am interested in the case where $X$ is normally distributed like above, but in addition $a<X<b$ with $a>0$ and $a, b$ both finite numbers (so a truncated normal distribution).

Would the moments of the distribution of $Y$ exist for that case? Intuitively, I would expect they do because the probability that $X=0$ is zero in this case, as apposed to the original question. However, I do not have the necessary math skills to rigourously work out the integral.

Can anybody tell me whether (a) the (first and second) moments exist and (b) if so, what they are?

Answer 1

We have $E|\frac 1 {X^{n}}| \leq \frac 1 {a^{n}} <\infty$, so all moments exist. I think explicit computation of the moments is not possible.

Answer 2

The moments exist because you have a bounded distribution on the positive real interval $\left[\frac1b,\frac1a\right]$: you are going to find $\frac1b < E[Y] < \frac1a$ and $\frac1{b^2} < E[Y^2] < \frac1{a^2}$.

It is easy enough to use numerical methods to find the moments to arbitrary accuracy. For example with R:

momentsreciptruncnormal <- function(a, b, mu=0, sigmasq=1, cases=10^4){
  values <- 1 / qnorm(pnorm(a, mu, sqrt(sigmasq)) + ppoints(cases) *
            (pnorm(b, mu, sqrt(sigmasq)) - pnorm(a, mu, sqrt(sigmasq))), 
             mu, sqrt(sigmasq))
  return(c(mean(values), mean(values^2)) 
  }

and as an example of the first and second moments with arbitrary $a,b,\mu,\sigma^2$:

momentsreciptruncnormal(a=1, b=2, mu=3, sigmasq=4)
# 0.6783127 0.4787007