Variance of a function of the sample variance I'm looking for the sampling standard deviation of $\hat\sigma^\gamma$, where $\hat\sigma$ is a sample standard deviation.  For simplicity, lets do the sample variance of the sample variance and take the roots later.  
Before doing any analytical math, I did the following simulation:
N = 10
mu=0
sig=1
gam = 2.5
p = rnorm(1000,mu,sig)
SDgam = c()
sdsd = c()
for (i in 1:4000){
    s = sample(p,N)
    SDgam[i] = sd(s)^gam
    par(mfrow=c(1,2))
    hist(SDgam)
    abline(v = mean(SDgam),col='red')
    sdsd[i] = sd(SDgam)
    if (i>1){
        plot(1:i,sdsd,cex=0)
        lines(1:i,sdsd,cex=0)
    }
}


The distribution is quite skewed when $\gamma=2.5$, which seems natural enough, and it looks like convergence is a bit slow, which also seems natural enough.
I start with the definition of the variance
$$
var(\hat\sigma^{2\gamma}) = E(\hat\sigma^{4\gamma}) - E(\hat\sigma^{2\gamma})^2
$$
then expand it out
$$
var(\hat\sigma^{2\gamma}) = \int\hat\sigma^{4\gamma}f(\hat\sigma^2)d\hat\sigma^2 - \left(\int\hat\sigma^{2\gamma}f(\hat\sigma^2)d\hat\sigma^2\right)^2
$$
I can use that fact that $\frac{(n-1)\hat\sigma^2}{\sigma^2}\sim\chi^2_{n-1}$ to rewrite the $f(\hat\sigma^2)$ from above as
$$
f(\hat\sigma^2) = \frac{\hat\sigma^{2(n/2-1)}e^{-\hat\sigma^2/2}\sigma^2}{2^{n/2}\Gamma(n/2)(n-1)}
$$
But it isn't clear to me what that gets me or whether it is helpful:
$$
var(\hat\sigma^{2\gamma}) = \int\hat\sigma^{4\gamma}\frac{\hat\sigma^{2(n/2-1)}e^{-\hat\sigma^2/2}\sigma^2}{2^{n/2}\Gamma(n/2)(n-1)}d\hat\sigma^2 - \left(\int\hat\sigma^{2\gamma}\frac{\hat\sigma^{2(n/2-1)}e^{-\hat\sigma^2/2}\sigma^2}{2^{n/2}\Gamma(n/2)(n-1)}d\hat\sigma^2\right)^2
$$
How would I go about getting a closed form here?  I want the variance of the sample variance as a function of the sample variance, $N$, and $\gamma$.  Apologies if I'm missing something obvious -- I haven't had to do analytical math in a while.
 A: When faced with a non-linear function of a random variable, we usually apply  the Delta Method (sometimes referred to as "error propagation" when the purpose is to obtain the variance) to obtain an approximate estimation.  
If we have a random variable $X$ and a non-linear function of it $h(X)$, then we apply a first-order Taylor expansion around a fixed point $x_0$:
$$h(X) \approx  h(x_0) + h'(x_0)\cdot (X-x_0)$$
Keeping in mind that $x_0$ is treated as a fixed number, consider
$${\rm Var}[h(X)] \approx {\rm Var}\big[h(x_0) + h'(x_0)\cdot (X-x_0)\big]$$
$$={\rm Var}\big[h'(x_0)\cdot (X-x_0)\big] = \left[h'(x_0)\right]^2{\rm Var}\big(X\big)$$
In your case, $h(X) = X^{\gamma} \implies [h'(x_0)]^2 = [\gamma x_0^{\gamma-1}]^2$, and your random variable is the sample standard deviation. Use the obtained estimate of it as your $x_0$ and also the estimated variance of the sample standard deviation for ${\rm Var}\big(X\big)$ (except if you know a priori the true value of this variance of course).
A: The Question


*

*Let $p$ denote the sample variance. Find $Var(p^\gamma)$  


Exact symbolic solutions can be attained for positive integer-valued $\gamma$. The answer will, of course, depend on what formula you are using for sample variance $p$.  Such problems are known as moments of moments problems and can be solved more easily by working with power sum  $s_r$ notation.
Solution
I will assume that by sample variance $p$, you are referring to the unbiased estimator of population variance $\mu_2$, i.e. the $2^{\text{nd}}$  h-statistic, namely:
$$ p   \quad=\quad \frac{1}{n-1}\sum _{i=1}^n \left(X_i-\bar{X}\right){}^2 \quad =  \quad \frac{n s_2-s_1^2}{n (n-1)} \quad \text{where} \quad s_r=\sum _{i=1}^n X_i^r$$
You seek $Var(p^\gamma)$. Note that the variance operator is the $2^{nd}$ Central Moment of $p^\gamma$. 


*

*In the case of $\gamma = 2$, the solution $\text{Var}(p^2)$ is:



where I am using the CentralMomentToCentral function from the mathStatica package for Mathematica to perform the calculation. The solution is expressed in terms of the central moments $\mu_i$ of the population of $X$.


*

*In the case of $\gamma = 3$, the solution $\text{Var}(p^3)$ is the $2^{\text{nd}}$ central moment of $p^3$:





The Normal case
The above results are completely general and hold for any distributions whose moments exist. In the special case of a Normal parent,  i.e. $X \sim N(\mu, \sigma^2)$, the first 12 central moments $\mu_i$ of the population are:
$\text{mulis}=\left\{\mu _2\to \sigma ^2,\mu _3\to 0,\mu _4\to 3 \sigma ^4, \mu _5\to 0, \mu _6\to 15 \sigma ^6, \mu _7\to 0, \mu _8\to 105 \sigma ^8, \mu _9\to 0,\mu _{10}\to 945 \sigma ^{10}, \mu _{11}\to 0, \mu _{12}\to 10395 \sigma ^{12}\right\}$
and the above solutions simplify to:

and

If your interest is in non-integer $\gamma$ (e.g.$\gamma = 2.5$) and need recourse to approximate solutions, the above exact solutions will provide a helpful benchmark to check how well any approximate method is working.
