Variance of a particular stochastic process Given the following stochastic process:
$$X_t = \epsilon_t\epsilon_{t-1}~~~~,~\epsilon_t \sim N(0, \sigma^2)$$
Which is the variance of the process $X_t^2$?
Firstly I defined the process $X_t$ as:
$$W_t = (\epsilon_t\epsilon_{t-1})^2= \epsilon_t^2\epsilon_{t-1}^2$$ 
Then its expected value as:
$$E[W_t]= E[\epsilon_t^2\epsilon_{t-1}^2]=E[\epsilon_t^2]E[\epsilon_{t-1}^2]= \sigma^4$$
And the variance as:
$$\text{Var}[W_t]=E[(\epsilon_t^2\epsilon_{t-1}^2-E[X_t])^2]$$
$$=E[\epsilon_t^4\epsilon_{t-1}^4]-2E[\epsilon_t^2\epsilon_{t-1}^2]\sigma^4+\sigma^8 = \sigma^8 - 2\sigma^8+\sigma^8=0$$
I am not able to understand where is my error.
 A: \begin{eqnarray*}
Var(X_t^2) &:=& Var(W_t) \\
&=& E[W_t^2] - E[W_t]^2 \\
&=& E[X_t^4] - E[X_t^2]^2 \\
&=& E[(\varepsilon_t\varepsilon_{t-1})^4] - E[(\varepsilon_t\varepsilon_{t-1})^2]^2
\end{eqnarray*}
Recall that the product of two standard normal random variables (that is, $N(0,1)$) is distributed according to a $\chi^2_1$ distribution. Then,
\begin{eqnarray*}
\varepsilon_t\varepsilon_t &=& \sigma^2\frac{\varepsilon_t}{\sigma}\frac{\varepsilon_t}{\sigma} \\
&:=& \sigma^2Y \\
\Rightarrow Y &\sim& \chi^2_1 \\
\Rightarrow \sigma^2Y &\sim& \Gamma\left(\frac{1}{2},2\sigma^2\right) \\
\Rightarrow \varepsilon_t\varepsilon_t &\sim& \Gamma\left(\frac{1}{2},2\sigma^2\right)
\end{eqnarray*}
The second implication in the above work is taken from this link. Note that the exact same holds for $\varepsilon_{t-1}\varepsilon_{t-1}$ because $\varepsilon_t$ and $\varepsilon_{t-1}$ follow the same distribution. Thus we can write $X^4$ as $(\varepsilon_t\varepsilon_t)^2(\varepsilon_{t-1}\varepsilon_{t-1})^2$, where $X$ follows a $\Gamma\left(\frac{1}{2},2\sigma^2\right)$ distribution.
For any random variable $X$, $E[X^n] = M^{(n)}_X(0)$, where $M^{(n)}_X(0)$ is the $n$th derivative of the moment generating function (mgf) of $X$, evaluated at 0. (This is Theorem 2.3.7 in Casella and Berger's Statistical Inference.) The moment generating function for a $\Gamma(\alpha,\beta)$ variable is $M_X(t)=\left(\frac{1}{1-\beta t}\right)^\alpha$. Then...
\begin{eqnarray*}
E[X_t^2] &=& M^{(2)}_X(0) \\
&=& \frac{d^2}{dt^2}\left(\frac{1}{1-\beta t}\right)^\alpha\biggr\rvert_{t=0} \\
&=& \frac{d^2}{dt^2}(1-\beta t)^{-\alpha}\biggr\rvert_{t=0} \\
&=& \frac{d}{dt}(-\alpha\cdot(1-\beta t)^{-\alpha-1}\cdot-\beta)\biggr\rvert_{t=0} \\
&=& \frac{d}{dt}(\alpha\beta(1-\beta t)^{-\alpha-1})\biggr\rvert_{t=0} \\
&=& \left((-\alpha-1)\cdot\alpha\beta(1-\beta t)^{-\alpha-2}\cdot-\beta\right)\biggr\rvert_{t=0} \\
&=& \left(\beta(\alpha+1)\alpha\beta(1-\beta t)^{-\alpha-2}\right)\biggr\rvert_{t=0} \\
&=& \left(\frac{\beta(\alpha+1)\alpha\beta}{(1-\beta t)^{\alpha+2}}\right)\biggr\rvert_{t=0} \\
&=& \left(\beta^2\alpha(\alpha+1)\right)
\end{eqnarray*}
You can find $E[X^4]$ in the same manner or by using formulas available online. Be careful, however, as some Web sites will have different parameterizations of the Gamma distribution.
\begin{eqnarray*}
E[X^4] &=& M^{(4)}_X(0) \\
&=& \beta^4\alpha(\alpha+1)(\alpha+2)(\alpha+3)
\end{eqnarray*}
From the top, we had: 
\begin{eqnarray*}
Var(X_t^2) &=& E[X_t^4] - E[X_t^2]^2 \\
&=& \beta^4\alpha(\alpha+1)(\alpha+2)(\alpha+3) - \left(\beta^2\alpha(\alpha+1)\right)^2 \\
&=& \beta^4\alpha(\alpha+1)(\alpha+2)(\alpha+3) - \beta^4\alpha^2(\alpha+1)^2 \\
Var(X_t^2) &=& \beta^4\alpha(\alpha+1)\left[(\alpha+2)(\alpha+3)-\alpha(\alpha+1)\right]
\end{eqnarray*}
A: The mistake of the OP lies in the last line of the post: in reality
$$E[\epsilon_t^4\epsilon_{t-1}^4] \neq \sigma^8$$
because, under assumed independence,
$$E[\epsilon_t^4\epsilon_{t-1}^4] = E[\epsilon_t^4]\cdot E[\epsilon_{t-1}^4]$$
and
$$E[\epsilon_t^4] = E[(\epsilon_t^2)^2] > [E(\epsilon_t^2)]^2 = (\sigma^2)^2 = \sigma^4$$
due to Jensen's Inequality.
