The first two moments of $\int_0^1 B_s^2 \, ds$ I was trying to solve the following problem from Continuous Martingales and Brownian Motion by Daniel Revuz and Marc Yor, but got my solution back as the answer for variance was wrong. I have already recomputed it several times and spelled out everything explicitly without any shortcuts, but still do not see what I am missing. I would appreciate any help.
Problem (Chapter I, Exercise 1.15)
Let $B_t$ be a standard linear Brownian motion. Compute the first two moments of
\begin{equation*}
X = \int_0^1 B_s^2 \, ds.
\end{equation*}
Solution
The expected value is
\begin{equation*}
  E(X) = \int_0^1 E \left( B_s^2 \right) \,ds = \frac{1}{2},
\end{equation*}
and the variance is
\begin{align*}
  & \text{Var}(X) = E\left(\left(\int_0^1 B_s^2 ds\right)^2\right) - \frac{1}{4} = E\left( \int_0^1 \int_0^1 (B_s B_t)^2 \,ds\, dt \right) - \frac{1}{4} \\
  & = E\left( \int_0^1 \int_0^t (B_s B_t)^2 \,ds \,dt \right) + E\left( \int_0^1 \int_t^1 (B_s B_t)^2 \,ds\, dt \right) - \frac{1}{4} \\
  & = E\left( \int_0^1 \int_0^t (B_s ((B_t - B_s) + B_s))^2 \,ds\, dt \right) \\
  & \qquad {} + E\left( \int_0^1 \int_t^1 (B_t ((B_s - B_t) + B_t))^2 \,ds\, dt \right) - \frac{1}{4} \\
  & = E\left( \int_0^1 \int_0^t \left[ B^2_s(B_t - B_s)^2 + B_s^4 + 2 B_s^3(B_t - B_s) \right] \,ds\, dt \right) \\
  & \qquad {} + E\left( \int_0^1 \int_t^1 \left[ B^2_t(B_s - B_t)^2 + B_t^4 + 2 B_t^3(B_s - B_t) \right] \,ds \,dt \right) - \frac{1}{4} \\
  & = \int_0^1 \int_0^t \left(s(t - s) + 3 s^2\right) \,ds\, dt + \int_0^1 \int_t^1 \left(t(s - t) + 3 t^2\right) \,ds\, dt - \frac{1}{4} \\
  & = \int_0^1 \int_0^t \left(st + 2 s^2\right) \,ds\, dt + \int_0^1 \int_t^1 \left(ts + 2 t^2\right) \,ds\, dt - \frac{1}{4} \\
  & = \int_0^1 t \int_0^t s \,ds\, dt + 2 \int_0^1 \int_0^t s^2 \,ds\, dt + \int_0^1 t \int_t^1 s \,ds\, dt + 2 \int_0^1 t^2 \int_t^1 \,ds\, dt - \frac{1}{4} \\
  & = \frac{1}{2} \int_0^1 t^3 dt + \frac{2}{3} \int_0^1 t^3 dt + \frac{1}{2} \int_0^1 t (1 - t^2) dt + 2 \int_0^1 t^2 (1 - t) dt - \frac{1}{4} \\
  & = \frac{1}{8} + \frac{1}{6} + \frac{1}{2} \left( \frac{1}{2} - \frac{1}{4} \right) + 2 \left( \frac{1}{3} - \frac{1}{4} \right) - \frac{1}{4} \\
  & = \frac{7}{24} + \frac{7}{24} - \frac{6}{24} = \frac{1}{3}.
\end{align*}
Thank you!
Regards,
Ivan

It turned out that my solution was actually correct.
In addition, here is a little empirical proof in MATLAB:
m = 10000;  % number of paths
n = 10000;  % number of integration segments
dt = 1 / n; % length of one segment

B = cumsum(normrnd(0, sqrt(dt), m, n), 2); % draw m paths of B
X = dt * sum(B.^2, 2); % integrate

E = mean(X)          % expectation
V = mean((X - E).^2) % variance

Running the code, we get
E =
    0.4981
V =
    0.3248

