It is well known that the integral $$\int_0^t W_s ds,$$ where $(W_s)_s$ is a Brownian motion, can be derived using Ito's Lemma. More precisely, Ito's lemma on $d(tW_t)$ implies that $$d(tW_t) = tdW_t + W_t dt.$$ Therefore, $$\int_0^t W_s ds = tW_t - \int_0^t sdW_s.$$ Its mean and variance can be obtained from this expression. This leads to my question below.
Question: Given a positive integer $n,$ what is the mean and variance $$\int_0^t (W_s)^n ds?$$
Calculation above is for $n=1.$
For the mean you can use Fubini's Theorem to change the order of integration $$ E\int_0^t (W_s)^n ds = \int_0^t E(W_s)^n ds$$ Then we can use the fact that $W_s \sim N(0,\sqrt{s})$ to obtain \begin{align*} E (W_s)^{2k+1} &= 0, k=0,1,2,... \tag*{(odd n)} \\ E (W_s)^{2k} &= (2k-1)!!s^k, k=1,2,... \tag*{(even n)} \end{align*} Therefore, \begin{align*} E\int_0^t (W_s)^{2k+1} ds &= 0, k=0,1,2,... \tag*{(odd n)} \\ E\int_0^t (W_s)^{2k}ds &= (2k-1)!! \int_0^t s^k ds = (2k-1)!! \frac{t^{k+1}}{k+1}, k=1,2,... \tag*{(even n)} \end{align*}
