  Stochastic volatility models based on Gaussian processes, like fractional
Brownian motion, are able to reproduce important stylized facts of financial
markets such as rich autocorrelation structures, persistence and roughness of
sample paths. This is made possible by virtue of the flexibility introduced in
the choice of the covariance function of the Gaussian process. The price to pay
is that, in general, such models are no longer Markovian nor semimartingales,
which limits their practical use. We derive, in two different ways, an explicit
analytic expression for the joint characteristic function of the log-price and
its integrated variance in general Gaussian stochastic volatility models. Such
analytic expression can be approximated by closed form matrix expressions. This
opens the door to fast approximation of the joint density and pricing of
derivatives on both the stock and its realized variance using Fourier inversion
techniques. In the context of rough volatility modeling, our results apply to
the (rough) fractional Stein--Stein model and provide the first analytic
formulae for option pricing known to date, generalizing that of Stein--Stein,
Sch{\"o}bel-Zhu and a special case of Heston.
