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.