Solution (source code)

= Solution
{c}

<Independent increments> give $B_s\sim N(0,s)$ independent of $B_t-B_s\sim N(0,t-s)$. Their <joint probability density function> at $(x,v)$ is the product
$$
\boxed{\frac1{2\pi\sqrt{s(t-s)}}\exp\left(-\frac{x^2}{2s}-\frac{v^2}{2(t-s)}\right).}
$$
The <change of variables> $(x,v)\mapsto(x,y)=(x,x+v)$ has <Jacobian determinant> one. Therefore the <joint probability density function> of $(B_s,B_t)$ is
$$
\boxed{\frac1{2\pi\sqrt{s(t-s)}}\exp\left(-\frac{x^2}{2s}-\frac{(y-x)^2}{2(t-s)}\right).}
$$
In particular these coordinates are correlated, with <covariance> $s$, although the first pair has independent coordinates.