Let $f:\Omega\rightarrow\mathbb{R}$ be a Borel-measurable function, $X$ a random variable with values in $\Omega$ and $X_i\in\mathbb{R},i\in\mathbb{N}$ realizations of $X$.
Literature
In their book on Monte Carlo Methods (Simulation and the Monte Carlo method, Second Edition) Rubinstein and Kroese give in Section 4.2.1 an approximate confidence interval for the mean $E[f(X)]$: $$\left(\hat{\mu}\pm z_{1-\frac{\alpha}{2}}\frac{\hat{\sigma_f}}{\sqrt{N}}\right),$$ where we write $N$ for the number of Monte Carlo samples, $z_{1-\frac{\alpha}{2}}$ for the $1-\frac{\alpha}{2}$ quantile of the standard normal distribution, $\hat{\mu}= \frac{1}{N}\sum_{i=1}^Nf(X_i)$ for the unbiased estimator of the mean and $\hat{\sigma}_f=\sqrt{\frac{1}{N-1}\sum_{i=1}^{N}(f(X_i)-\hat{\mu})^2}$ for the empirical standard deviation derived from the unbiased estimator of the variance.
Extension to variance
In my opinion we can derive in a similar way a confidence interval for the variance $\operatorname{var}[f(X)]$ if we define $g(X):=(f(X)-\hat{\mu})^2$ and repeat the construction. This leads to the interval $$\left(\hat{\sigma}_f^2\pm z_{1-\frac{\alpha}{2}}\frac{\hat{\sigma}_g}{\sqrt{N}}\right),$$ where $\hat{\sigma}_g=\frac{1}{N-1}\sum_{i=1}^N\left(g(X)-\frac{1}{N}\sum_{i=1}^Ng(X_i)\right)^2$.
Question
How can I transfer this to give a confidence interval for the standard deviation? Surely I can't just take the square root of the boundaries. Any literature, help and thoughts are welcome. Thanks.
The approach is different and yields a confidence interval of the form $$\left(\sqrt{\frac{(N-1)\hat{\sigma}^2_f}{\chi^2_{\alpha/2}}},\sqrt{\frac{(N-1)\hat{\sigma}^2_f}{\chi^2_{1-\alpha/2}}}\right),$$ where $\chi^2_\alpha$ denotes the $\alpha$-quantile of the $\chi^2$ distribution with N-1 degrees of freedom.