Using $\delta$-method to "estimate" undefined moments of a random variable? I posted this over on MSE without much luck. Not sure if posting here is considered cross-posting but I can remove it if it is.
Let $X\sim\mathcal N(\sqrt 2,1/x^2)$. The expected value $\mathsf EX^{-1}$ is undefined; however, we can assign it a value via the Cauchy principal value interpretation
$$
\mathsf{E}X^{-1}\overset{\text{p.v.}}{=}\lim_{\epsilon\nearrow 0}\left(\int_{-\infty}^{-\epsilon}+\int_\epsilon^\infty\right)\frac{1}{t}f_X(t)\,\mathrm dt=\sqrt 2 x\mathcal D(x),
$$
where $\mathcal D(x):=e^{-x^2}\int_0^xe^{t^2}\,\mathrm dt$ is Dawson's integral. Notice that as $x\to\infty$, $\mathsf{Var}X\to 0$ and so I wondered what would happen if I were to apply the $\delta$-method to estimating $\mathsf EX^{-1}$. Let $g(X)=1/X$. Expanding $g$ in a Taylor series about $\mathsf EX=\sqrt 2$ gives
$$
g(X)=\frac{1}{\sqrt 2}\sum_{k=0}^\infty\left(\frac{\sqrt 2-X}{\sqrt 2}\right)^k
$$
and so evaluating the expected value simply requires knowing central moments of the normal distribution. We find
$$
\mathsf Eg(X)
%=\frac{1}{\sqrt 2}\sum_{k=0}^\infty\mathsf E\left(\frac{\sqrt 2-X}{\sqrt 2}\right)^k
%=\frac{1}{\sqrt 2}\sum_{k=0}^\infty\mathsf E\left(\frac{X-\sqrt 2}{\sqrt 2}\right)^{2k}
=\frac{1}{\sqrt 2}\sum_{k=0}^\infty\frac{(2k-1)!!}{(2x^2)^k}.
$$
Looking back at our expression for $\mathsf EX^{-1}$ it stands to reason that if we then divide $\mathsf Eg(X)$ by $\sqrt 2 x$ that we obtain an estimate for $\mathcal D(x)$ at large $x$, namely,
$$
\mathcal D(x)\sim\frac{1}{2x}\sum_{k=0}^\infty\frac{(2k-1)!!}{(2x^2)^k}=\frac{1}{2x}+\frac{1}{4x^3}+\mathcal O(x^{-5}).
$$
Comparing these first two terms with eqn. (9) here seems to indicate the above expression is indeed an asymptotic expansion of $\mathcal D(x)$ for large $x$ and thus $\mathsf E g(X)$ provides us with an asymptotic expansion for the Cauchy principal value of  $\mathsf EX^{-1}$ as $x\to\infty$.
There is no reason to stop here as we could go further and apply the same approach to "estimate" $\mathsf{Var}X^{-1}$ yielding
$$
\mathsf{Var}X^{-1}\sim
%\mathsf Eg^2(X)-(\mathsf Eg(X))^2
\frac{1}{2}\sum_{k=0}^\infty\frac{(2k+1)(2k-1)!!}{(2x^2)^k}-\frac{1}{2}\left(\sum_{k=0}^\infty\frac{(2k-1)!!}{(2x^2)^k}\right)^2
=\frac{1}{4x^2}+\frac{1}{x^4}+\mathcal O(x^{-6}).
$$

My question fundamentally has to do with what exactly the $\delta$-method is estimating when we apply it to estimating moments of random variables for which the moments to not exist? In this specific case, we could claim (at least I think) that the $\delta$-method gave us an asymptotic expansion for the Cauchy principal value of $\mathsf EX^{-1}$. However, for the second moment, the Cauchy principal value interpretation of the integral would give us $\mathsf EX^{-2}\overset{\text{p.v.}}{=}\infty$ whereas the $\delta$-method gave us a finite expression. So what the heck did the $\delta$-method estimate?

 A: $\newcommand\vp\varepsilon$

*

*Of course, your divergent series for $EX^{-1}$ and $EX^{-2}$ should be understood as asymptotic expansions. The delta method practically never involves series; it involves asymptotic expansions instead, with an appropriately controlled remainder.


*Take any natural $k$. With the modifications discussed above, your delta method estimates a lot of things. In particular, for $X\sim N(a,1/x^2)$ with a real $a>0$ and $x\to\infty$ and for any $\vp\in(0,a/2)$, it estimates
$$EX^{-k}1(|X|>\vp)=I_{k,a}+J_{k,a}(\vp),$$
where
$$I_{k,a}:=EX^{-k}1(|X-a|<a/2),$$
$$J_{k,a}(\vp):=EX^{-k}1(|X|>\vp,|X-a|>a/2).$$
We have
$$|J_{k,a}(\vp)|\le\vp^{-k}P(|X-a|>a/2)=2\vp^{-k}P(Z>ax/2)=o(x^{-m})\tag{1}$$
for any natural $m$, where $Z\sim N(0,1)$.
Using the series
$$x^{-k}=\sum_{n=0}^\infty \binom{-k}{n}a^{-k-n}(x-a)^n$$
for $x$ with $|x-a|<a/2$, we get the asymptotic expansion
$$
\begin{aligned}I_{k,a}&\sim\sum_{n=0}^\infty \binom{-k}{n}a^{-k-n} E(X-a)^n1(|X-a|<a/2) \\ 
&\sim\sum_{n=0}^\infty \binom{-k}{n}a^{-k-n} E(X-a)^n+o(x^{-m})
\end{aligned}$$
for any natural $m$ (cf. (1)). So,
$$EX^{-k}1(|X|>\vp)\sim\sum_{n=0}^\infty \binom{-k}{n}a^{-k-n} E(X-a)^n.$$
In particular,
$$EX^{-k}1(|X|>\vp)=\frac1{a^k}\Big(1+\frac{(k+1) k}{2 a^2 x^2}+\frac{(k+3) (k+2) (k+1) k}{8 a^4 x^4}\Big)+O(x^{-6}).$$
These results will hold if, instead of fixing $\vp$, we allow $\vp$ to go to $0$, but not overly too fast: in particular,
the following very mild requirement is enough: $\vp>e^{-ba^2x^2/k}$ for some real $b\in(0,1/2)$ and all large enough $x>0$.

(For $k=1$, $EX^{-k}$ exists in $\mathbb R$ in the principal value sense -- that is, as $\lim_{\vp\downarrow0}EX^{-1}1(|X|>\vp)$, but this does not hold for any other natural $k$. This distinction has little, if anything, to do with the delta method.)
A: Too long for comment. As pointed out by @Iosif, $\mathsf EX^{-k}$ does not exist in the Cauchy principal value sense for $k=2,3,\dots$; however, we may regularize the high order negative moments as
$$
\mathsf E_\mathcal PX^{-n}:=\frac{1}{(n-1)!}\partial_t^{n-1}\mathsf E_\mathcal P(X-t)^{-1}\Big|_{t=0},
$$
with $\mathsf E_\mathcal P(X-t)^{-1}$ being defined by the Cauchy principal value.
Likewise, we may define the higher order negative moments in the sense of the "delta method" by making use of the series expansion
$$
\mathsf E_\delta X^{-n}:=\sum_{k=0}^\infty \binom{-n}{k}\mu^{-n-k}\mathsf E(X-\mu)^k.
$$
Using @Iosif's analysis one can show that the delta method moments $\mathsf E_\delta$ yield asymptotic expansions for the regularized moments $\mathsf E_\mathcal P$ when $|\mu/\sigma|\to\infty$.
For example, consider $X\sim\mathcal N(\mu,\sigma^2)$ so that
$$
\mathsf E_\mathcal PX^{-2}=\frac{1}{\sigma^2}\sqrt 2\frac{\mu}{\sigma}\mathcal D\left(\frac{\mu/\sigma}{\sqrt 2}\right)-\frac{1}{\sigma^2}
$$
and as $|\mu/\sigma|\to\infty$ we may utilize the asymptotic expansion of the Dawson integral to write
$$
\mathsf E_\mathcal PX^{-2}\sim\frac{1}{\sigma^2}\sum_{k=1}^\infty(2k-1)!!\left(\frac{\sigma}{\mu}\right)^{2k}.
$$
Now using our delta method we find
$$
\begin{align}
\mathsf E_\delta X^{-2}
&=\sum_{k=0}^\infty (-1)^k(k+1)\mu^{-2-k}\mathsf E(X-\mu)^k\\
&=\sum_{k=0}^\infty (2k+1)\mu^{-2-2k}\sigma^{2k}(2k-1)!!\\
&=\frac{1}{\sigma^2}\sum_{k=1}^\infty(2k-1)!!\left(\frac{\sigma}{\mu}\right)^{2k},
\end{align}
$$
which is precisely our asymptotic expansion for $\mathsf E_\mathcal PX^{-2}$ when $|\mu/\sigma|\to\infty$.
