I am trying to understand stochastic calculus and got stuck calculating the following.
I need the distribution of a zero bond under the black model, so I am deriving the variance using the second moment of $\int_{t}^{t+h} r_s\,ds.$
Given the differential equation for $r_t$ (with $W$ standard Brownian motion): $$dr_t \;=\; \mu dt + \sigma dWt^Q,$$
I want to calculate
$$\mathbb{E}_t^Q\left[\,\left(\int_{t}^{t+h} r_s\,ds\right)^2\,\right].$$
I have gotten this far:
$$=\;\mathbb{E}_t^Q\left(\int_{t}^{t+h} r_s\,ds\,\int_{t}^{t+h}r_u\,du\right)$$
$$=\;\int_{t}^{t+h}\int_{t}^{t+h} \big(\mathbb{E}_t^Q r_s r_u\big)\,ds\,du.$$
Now I know that $r_s$ and $r_u$ are normally distributed, but where to go from here?
Could you maybe use $Var[\int_t^{t+h}r_sds]=E[(\int_t^{t+h}r_sds)^2]-(E[\int_t^{t+h}r_sds])^2 $. Then if $r_s$ is Normal, then the integral itself follows Normal distribution with $\mu=\int_t^{t+h}E[r(s)]ds$ and $\sigma^2=\int_t^{t+h} \int_t^{t+h}E[r_s r_t]dtds$
Then you yo must solve the SDE and plug everything in