This question has been migrated from the MSE.
Background/Motivation:
We have $Z=X/Y$ where $X$ and $Y$ are independent and $X\sim\mathcal N(\mu,\sigma^2)$. The density of $Y$ is not important here. We can write the distribution and density functions of $Z$ in terms of expected values w.r.t. $Y$ as $$ F_Z(z)=\mathsf E\Phi\left(\frac{z|Y|-\operatorname{sign}(Y)\mu}{\sigma}\right) $$ and $$ f_Z(z)=\mathsf E\left(\frac{|Y|}{\sigma}\phi\left(\frac{zY-\mu}{\sigma}\right)\right), $$ where $\Phi(\cdot)$ and $\phi(\cdot)$ represent the standard normal cdf and pdf, respectively. This leads to unbiased Monte Carlo estimators of the distribution and density functions. For example, given a sample $Y_1,\dots,Y_n$ we can estimate the distribution function $F_Z$ at the point $z$ via $$ \hat F_Z(z)=\frac{1}{n}\sum_{k=1}^n\Phi\left(\frac{z|Y_k|-\operatorname{sign}(Y_k)\mu}{\sigma}\right) $$ I am interested in evaluating the variance of these estimators as a function of $z$, i.e. $\mathsf{Var}(\hat F_Z)(z)$ and $\mathsf{Var}(\hat f_Z)(z)$.
Approach:
It turns out in my application $\sigma\ll\mathsf{Var}Y$ so much so that $X$ looks nearly constant in comparison to $Y$. As such, taking limit $\sigma\to 0$ in the above expressions still gives good approximations to the cdf/pdf of $Z$. For example, taking the limit $\sigma\to0$ in the expression for the cdf we make use of the fact that the normal cdf tends to a step function giving the approximation $$ F_Z(z)\approx\mathsf E(\mathbf 1_{z|Y|-\operatorname{sign}(Y)\mu>0}), $$ and so we have the corresponding MC estimator $$ \hat F_Z(z)\approx\frac{1}{n}\sum_{k=1}^n\mathbf 1_{z|Y_k|-\operatorname{sign}(Y_k)\mu>0}. $$ This approximation is very convenient because $\mathbf 1_{z|Y|-\operatorname{sign}(Y)\mu>0}$ is Bernoulli distributed with success probability $p=\mathsf E(\mathbf 1_{z|Y|-\operatorname{sign}(Y)\mu>0})\approx F_Z(z)$, that is we have the distributional approximation $\mathbf 1_{z|Y|-\operatorname{sign}(Y)\mu>0}\sim\operatorname{Binomial}(1,F_Z(z))$. As such we obtain the approximation $$ \mathsf{Var}(\hat F_Z)(z)\approx\frac{F_Z(z)(1-F_Z(z))}{n}. $$ I performed simulations to estimate $\mathsf{Var}(\hat F_Z)(z)$ and compared the estimates to this approximation which showed excellent agreement. However, I am unable to see how to extend this idea to estimate $\mathsf{Var}(\hat f_Z)(z)$.
Given $Y$ we note that $$ |Y|\frac{1}{\sigma}\phi\left(\frac{zY-\mu}{\sigma}\right) =\frac{1}{\sqrt{2\pi}\sigma/|Y|}\exp\left(-\frac{(zY-\mu)^2}{2\sigma^2}\right) =\frac{1}{\sqrt{2\pi}\sigma/|Y|}\exp\left(-\frac{(z-\mu/Y)^2}{2(\sigma/|Y|)^2}\right), $$ which is a normal density with mean $\mu/Y$ and variance $(\sigma/|Y|)^2$. So taking the limit $\sigma\to 0$ in the expression for $f_Z$ gives $$ f_Z(z)\approx \mathsf E(\delta(z-\mu/Y)) $$ and the corresponding "estimator" $$ \hat f_Z(z)\approx\frac{1}{n}\sum_{k=1}^n\delta(z-\mu/Y_k). $$ But here is where I run into trouble. In calculating the variance $\mathsf{Var}(\hat f_Z)(z)$ using this approximation we would have to evaluate $\mathsf E\delta^2(z-\mu/Y)$, which I do not know what to do with. How do I proceed? Why did this approach work for approximating $\mathsf{Var}(\hat F_Z)(z)$ but runs into problems in estimating $\mathsf{Var}(\hat f_Z)(z)$?