I'm trying to find the mean and correlation functions of $ X_t:=W_t^2 $, where $W_t$ is a standard Wiener process ($\sigma=1,E\{W_t^2\}=t$ for $t\geq 0$).
Unless I'm missing something large here, the mean is very easy to compute.
$m_x(t) = E[X_t] = E[W_t^2] = \sigma^2 t = t$, for $t\geq 0$
Where I run into trouble is the correlation function. I've omitted some intermediate steps for time's sake that I could include if asked.
$$R_x(t_1,t_2)=E[X_{t_1}X_{t_2}]=E[W_{t_1}^2 W_{t_2}^2]=E[W_{t_1}W_{t_1}W_{t_2}W_{t_2}] = E[W_{t_1}^4] + (\sigma^2t_1)(\sigma^2(t_2-t_1))=E[W_{t_1}^4]+t_1(t_2-t_1) $$
I run into trouble calculating $E[W_{t_1}^4]$. I saw a solution somewhere online that led me to believe it should be $3\sigma^4t_1^2$, but I have no idea how to arrive at this answer.
I've been trying to use the method for solving $E[V_t^2]$, where $V(t)$ is integrated white noise, $V(t)=\int_0^tX_\tau d\tau$ to give an idea...
$$E[V_t^2] = E[(\int_0^tX_\tau d\tau)(\int_0^tX_\theta d\theta)]=\int_{0}^{t} \int_{0}^{t} R_{X}(\tau-\theta) d \tau d \theta=\int_{0}^{t} \int_{0}^{t} \sigma^{2} \delta(\tau-\theta) d \tau d \theta=\sigma^{2} t$$
... but I seriously doubt it's as easy as squaring the result.
A nudge in the right direction for solving $E[W_{t_1}^4]$ would be much appreciated (assuming it's correct in the first place)
Ito's Lemma: If $dX_t=\mu(t,X_t)dt+\sigma(t,X_t)dW_t$ and $f(t,x)$ is twice continuously differentiable then $$ df(t,X_t)=\left[f_t(t,X_t)+\mu(t,X_t)f_x(t,X_t)+\frac{1}{2}\sigma^2(t,X_t)f_{xx}(t,X_t)\right]dt+\Big[\sigma(t,X_t) f_{x}(t,X_t)\Big]dW_t $$ So, in your case take $dX_t=0dt+1dW_t=dW_t$ and $f(t,x)=x^4$. Then you have $\mu(t,x)=0$, $\sigma(t,x)=1$, $f_t=0$, $f_x=4x^3$ and $f_{xx}=12x^2$. Hence
\begin{eqnarray*} dW_t^4&=&df(t,W_t)\\ &=&\left[0+0f_x(t,W_t)+\frac{1}{2} 1^2f_{xx}(t,W_t)\right]dt+\Big[1f_{x}(t,W_t)\Big]dW_t\\ &=& \Big[6W_t^2\Big]dt+\Big[4W_t^3\Big]dW_t \end{eqnarray*}
So, we have $$ W_t^4 = 6\int_0^tW_s^2ds + 4\int_0^tW_s^3dW_s $$ which means that \begin{eqnarray*} E[W_t^4] &=& 6E\left[\int_0^tW_s^2ds\right] + 4E\left[\int_0^tW_s^3dW_s\right]\\ &=& 6\int_0^tE\left[W_s^2\right]ds + 0\\ &=& 6\int_0^t s ds\\ &=& 3t^2 \end{eqnarray*}
Hence \begin{eqnarray*} E\Big[W^2_{t+h}W^2_t\Big] &=& E\Big[E\Big[W^2_{t+h}W^2_t\Big|{\cal F}_t\Big]\Big]\\ &=& E\Big[E\Big[W^2_{t+h}\Big|{\cal F}_t\Big]W^2_t\Big]\\ &=& E\Big[E\Big[(W_{t+h}-W_t+W_t)^2\Big|{\cal F}_t\Big]W^2_t\Big]\\ &=& E\Big[E\Big[(W_{t+h}-W_t)^2+2(W_{t+h}-W_t)W_t+W_t^2\Big|{\cal F}_t\Big]W_t^2\Big]\\ &=& E\Big[\Big(E\Big[(W_{t+h}-W_t)^2\Big|{\cal F}_t\Big]+2E\Big[(W_{t+h}-W_t)W_t\Big|{\cal F}_t\Big]+E\Big[W_t^2\Big|{\cal F}_t\Big]\Big)W_t^2\Big]\\ &=& E\Big[\Big(h+2E\Big[W_{t+h}-W_t\Big|{\cal F}_t\Big]W_t+W_t^2\Big)W_t^2\Big]\\ &=& E\Big[(h+0+W_t^2)W_t^2\Big]\\ &=& hE\Big[W^2_t\Big]+E\Big[W_t^4\Big]\\ &=& ht+3t^2 \end{eqnarray*} where $\cal{F}$ is the filtration generated by the Brownian motion.