Variance of Integrated Geometric Brownian Motion

Question

I'm just asking for verification that my derivation is correct, as I can't seem to find this result elsewhere. I'd like to calculate $Var(\int_0^T X(t) dt)$ where $X(t) = X_0e^{(\mu - \frac{1}{2}\sigma^2)t + \sigma W(t)}$, $W(t)$ is a standard Brownian motion and $\mu, \sigma$ and $T$ are constants. By definition $$ \operatorname{Var}\left(\int_0^T X(t) dt\right) = \mathbb{E}\left[\int_0^T\int_0^T X(t)X(s) dt ds\right] - \mathbb{E}\left[\int_0^T X(t) dt\right]^2. $$ Assuming $\mu \neq 0$ the second term is easily calculated to be (neglecting the square) $$ \begin{align*} \mathbb{E}\left[\int_0^T X(t) dt\right] & = \int_0^T \mathbb{E}\left[X(t)\right] dt \\ & = \int_0^T X_0e^{\mu t} dt \\ & = \frac{X_0}{\mu}(e^{\mu T} - 1). \end{align*} $$ The first term is a bit nastier. Assume $\mu \neq 0$, $\mu + \sigma^2 \neq 0$, and $2\mu + \sigma^2 \neq 0$ (these terms appear in denominators after the final integration, which is left out). Then using this result with no restriction on the order of $s,t$ to get from the second to the third line I get: $$ \begin{align*} \mathbb{E}\left[\int_0^T\int_0^T X(t)X(s) dt ds\right] & = \int_0^T\int_0^T \mathbb{E}\left[X(t)X(s) \right] dt ds \\ & = X_0^2\int_0^T\int_0^T e^{(\mu - \frac{1}{2}\sigma^2)(t+s)}\mathbb{E}\left[e^{\sigma(W_t + W_s)} \right] dt ds \\ & = X_0^2\int_0^T\int_0^T e^{(\mu - \frac{1}{2}\sigma^2)(t+s)} e^{\frac{1}{2}\sigma^2(3\min\{s,t\} + \max\{s,t\})} dt ds \\ & = X_0^2\int_0^T\int_0^s e^{(\mu - \frac{1}{2}\sigma^2)(t+s)} e^{\frac{1}{2}\sigma^2(3t + s)} dt \\ & \qquad + \int_s^T e^{(\mu - \frac{1}{2}\sigma^2)(t+s)} e^{\frac{1}{2}\sigma^2(3s + t)} dtds \\ & = X_0^2\int_0^T\left(\int_0^s e^{\mu s} e^{(\mu + \sigma^2)t} dt + \int_s^T e^{\mu t} e^{(\mu + \sigma^2)s} dt\right) ds. \end{align*} $$ The rest is just a tedious integration. If anyone knows where I can find another derivation of this to compare with my result it would be appreciated. Thanks!