What exactly do delta method estimates of moments for $1/\bar X_n$, $\bar X_n\sim\mathcal N(\mu,\sigma^2/n)$ approximate? (not as simple as you think)

Question

Let me start with the excerpt out of Casella & Berger's Statistical Inference (2nd edition, pg. 470) that inspired this question.

---

Definition 10.1.7 For an estimator $T_n$, if $\lim_{n\to\infty}k_n\mathrm{Var}T_n=\tau^2<\infty$, where $\{k_n\}$ is a sequence of constants, then $\tau^2$ is called the limiting variance or limit of the variances.

Example 10.1.8 (Limiting variances) For the mean $\bar X_n$ of $n$ iid normal observations with $\mathrm EX=\mu$ and $\mathrm{Var}\,X=\sigma^2$, if we take $T_n=\bar X_n$, then $\lim n\mathrm{Var}\bar X_n=\sigma^2$ is the limiting variance of $T_n$.

But a troubling thing happens if, for example, we are instead interested in estimating $1/\mu$ using $1/\bar X_n$. If we now take $T_n=1/\bar X_n$, we find that the variance is $\mathrm{Var}\,T_n=\infty$, so the limit of the variances is infinity. But recall Example 5.5.23, where we said that the "approximate" mean and variance of $1/\bar X_n$ are $$ \mathrm E\left(\frac{1}{\bar X_n}\right)\approx\frac{1}{\mu}, $$ $$ \mathrm{Var}\left(\frac{1}{\bar X_n}\right)\approx\left(\frac{1}{\mu}\right)^4\mathrm{Var}\bar X_n, $$ and thus by this second calculation the variance is $\mathrm{Var}\,T_n\approx\frac{\sigma^2}{n\mu^4}<\infty$.

This example points out the problems of using the limit of the variances as a large sample measure. Of course the exact finite sample variance of $1/\bar X$ is $\infty$. However, if $\mu\neq 0$, the region where $1/\bar X$ gets very large has probability going to $0$. So the second approximation in Example 10.1.8 is more realistic (as well as being much more useful). It is this second approach to calculating large sample variances that we adopt.

---

Now, I disagree that $\mathrm{Var}\left(1/\bar X_n\right)=\infty$. While I do agree that $\mathrm E\left(1/\bar X_n^2\right)=\infty$, the first negative moment $\mathrm E\left(1/\bar X_n\right)$ is clearly undefined and thus so is the variance. Given that both $\mathrm E\left(1/\bar X_n\right)$ and $\mathrm{Var}\left(1/\bar X_n\right)$ are strictly undefined, I have the following question.

Question: What exactly do the delta method approximations for $\mathrm E\left(1/\bar X_n\right)$ and $\mathrm{Var}\left(1/\bar X_n\right)$ approximate if the moments of $1/\bar X_n$ are undefined?

My thoughts:

Let $\bar X_n\sim\mathcal N(\mu,\sigma^2/n)$ and define for $t\in\Bbb R$ $$ \mathcal H[f_{\bar X_n}](t)=\lim_{\epsilon\to 0^+}\int_{\Bbb R\setminus(t-\epsilon,t+\epsilon)}\frac{f_{\bar X_n}(x)}{x-t}\,\mathrm dx, $$ which is the Hilbert transform of the density function for $\bar X_n$. Since the Hilbert transform commutes with derivatives, i.e. $$ \mathcal H[\partial_t^k u]=\partial_t^k\mathcal H[u], $$ the expression $\mathcal H[f_{\bar X_n}](t)$ represents a sort of generating function for the negative moments of $\bar X_n$, which do not exist in the traditional sense. We define $$ \mathrm E\bar X_n^{-k}:=\frac{1}{(k-1)!}\partial_t^{k-1}\mathcal H[f_{\bar X_n}](t)\Big|_{t=0}. $$ For our particular example $$ \mathcal H[f_{\bar X_n}](t)=\frac{\sqrt 2}{\sigma/\sqrt n}\mathcal D\left(\frac{\mu-t}{\sqrt 2\,\sigma/\sqrt n}\right), $$ with $\mathcal{D}(z)=e^{-z^{2}}\int_{0}^{z}e^{t^{2}}\,\mathrm{d}t$ being Dawson's integral; thus $$ \mathrm E\bar X_n^{-1}=\frac{\sqrt 2}{\sigma/\sqrt n}\mathcal D\left(\frac{\mu}{\sqrt 2\,\sigma/\sqrt n}\right) $$ and $$ \mathrm{Var}\bar X_n^{-1}:=\mathrm E\bar X_n^{-2}-(\mathrm E\bar X_n^{-1})^2 $$ with $$ \mathrm E\bar X_n^{-2}=\frac{\sqrt 2\,\mu}{(\sigma/\sqrt n)^3}\mathcal D\left(\frac{\mu}{\sqrt 2\,\sigma/\sqrt n}\right)-\frac{1}{(\sigma/\sqrt n)^2}. $$ Now, the Dawson integral admits the following asymptotic expansion for $x\to\infty$ $$ \mathcal D(x/\sqrt 2)\sim\frac{1}{\sqrt 2\,x}\sum_{k=0}^\infty (2k-1)!!\frac{1}{x^{2k}}, $$ so letting $n\to\infty$ we have for the first moment $$ \begin{align} \mathrm E\bar X_n^{-1} &\sim\frac{1}{\mu}\sum_{k=0}^\infty (2k-1)!!\left(\frac{\sigma/\sqrt n}{\mu}\right)^{2 k}\\ &=\frac{1}{\mu}+\mathcal O\left(\frac{1}{n}\right), \end{align} $$ which is the first-order delta method approximation for the mean of $1/\bar X_n$. Likewise, for the second moment we find $$ \mathrm E\bar X_n^{-2} \sim\frac{1}{(\sigma/\sqrt n)^2}\sum_{k=1}^\infty (2k-1)!!\left(\frac{\sigma/\sqrt n}{\mu}\right)^{2 k}, $$ which upon combining with the expansion for $\mathrm E\bar X_n^{-1}$ gives the asymptotic approximation for $n\to\infty$ $$ \begin{align} \mathrm{Var}\bar X_n^{-1} &=\frac{1}{\mu^2}+3\frac{(\sigma/\sqrt n)^2}{\mu^4}+\mathcal O\left(\frac{1}{n^2}\right)-\left(\frac{1}{\mu^2}+2\frac{(\sigma/\sqrt n)^2}{\mu^4}+\mathcal O\left(\frac{1}{n^2}\right)\right)\\ &=\frac{(\sigma/\sqrt n)^2}{\mu^4}+\mathcal O\left(\frac{1}{n^2}\right), \end{align} $$ which is the first-order delta method approximation for the variance of $1/\bar X_n$.

So in this context, the delta method moment "approximations" are equal to the first term in the asymptotic expansions for our exact generalized moments in the case $n\to\infty$.

To ask my question from a different perspective: If we are to assign any value to the moments of $1/\bar X_n$, should we instead assign these generalized values for $E\bar X_n^{-k}$ with the delta method moment estimates simply being an approximation of these "exact" moments?

Indeed, the generalized moments seem to provide a more accurate picture of the moments for $1/\bar X_n$ whenever $|\mu/(\sigma/\sqrt n)|\to\infty$ ($n$ does not necessarily have to be large), which implies $f_{\bar X_n}(0)$ is vanishingly small.

Consider the case $\mu=6$, $\sigma=1$, and $n=1$ so that $\bar X_n=X\sim\mathcal N(\mu,\sigma^2)$. I simulated $X$ in MATLAB with the following code:

mu = 6;
sigma = 1;
Y = 1./normrnd(mu,sigma,1e8,1);

EY = sqrt(2)/sigma*dawson(mu/(sqrt(2)*sigma));
[EY 1/mu mean(Y)]

EY2 = sqrt(2)*mu/sigma^3*dawson(mu/(sqrt(2)*sigma))-1/sigma^2;
[EY2 1/mu^2 mean(Y.^2)]

EY3 = (sqrt(2)*(mu^2-sigma^2)*dawson(mu/(sqrt(2)*sigma))-mu*sigma)/(2*sigma^5);
[EY3 1/mu^3 mean(Y.^3)]

VarY = EY2-(EY)^2;
[VarY sigma^2/mu^4 var(Y)]

fY = @(y) normpdf(1./y,mu,sigma)./y.^2;
ax = linspace(0.05,0.35,512);

figure
hold on;
histogram(Y,linspace(0.05,0.35,40),'Normalization','pdf')
plot(ax,fY(ax),'Color',[0 0 0],'LineWidth',2)

Here are the results showing that the generalized moments more accurately reflect the sample statistics.

$$ \begin{array}{cccc} &\text{generalized} &\text{delta} &\text{sample}\\ \mathrm EX^{-1} &0.1718 &0.1667 &0.1718 \\ \mathrm EX^{-2} &0.0305 &0.0278 &0.0305 \\ \mathrm EX^{-3} &0.0056 &0.0046 &0.0056 \\ \mathrm{Var}X^{-1} &0.0010 &0.0008 &0.0010 \end{array} $$

Answer 1

Too long for comment.

The asymptotic variances are more useful than the variances because some estimators do not have variance for any finite sample and yet the do have asymptotic variance. Consider $Y_1$ and $Y_2$ two independent binomial distributions. The odds ratio is $\theta = \dfrac{\frac{P(Y_2=1)}{P(Y_2=0)}}{\frac{P(Y_1=1)}{P(Y_1=0)}}.$ Write now a (2,2) table of observations as $(n_{ij})$ with $i=1,2,$ denoting $Y_1$ or $Y_2,$ and $j = 1$ denotes failure, while $j = 2$ denotes success. The most natural estimate is $\hat\theta = \dfrac{n_{11}n_{22}}{n_{12}n_{21}}$ which takes the values $\infty$ if $n_{12}n_{21} = 0,$ so its variance does not exists. However, for very large samples, $n_{12}n_{21} = 0$ has an exponentially small probability and this event will not occur in practice and so for all intents and purposes we should be able to calculate the variance. In this sense, the sample variance (not as defined by Casella and Berger but by the "delta method") is the one that we use in practice. The odds ratio can also be written as $\theta = \dfrac{P(Y_2 = 1)}{P(Y_1 = 1)} \dfrac{P(Y_1 = 0)}{P(Y_2 = 0)}.$ When $Y_1$ and $Y_2$ model something rare (so that "success" is a rare event, as in effectiveness of new treatment and we are comparing two treatments, or number of people who use drugs and we are comparing two social factors or number of car accidents than result in death comparing seat-bealt use vs not) an estimator of interest is the relative risk $r = \dfrac{P(Y_2 = 1)}{P(Y_1 = 1)}$ which cannot be estimated usually due to ethical or cost contraints (it can be estimated in clinical trials but not in the social-factor/drug-use scenario or the fatal-accident/seat-belt use), since the probabilty of "failure" is rare for both $Y_1$ and $Y_2$ the odds ratio is very close to the relative risk since $\dfrac{P(Y_1 = 0)}{P(Y_2 = 0)} \approx 1$ and in practice it is not rare that this fraction is $\in [0.97, 1.3].$ (By the way, it can be shown that the odds ratio is an estimator that is invariant under changes of scales, swaping of columns and or of rows, so it is really invariant under retrospective or prospective analyses.)