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
stemming from Wishart distributions. 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