We consider stochastic differential equations of the form dY(t) = V(Y-t)dX(t) + V-0(Y-t)dt driven by a multi-dimensional Gaussian process. Under the assumption that the vector fields V-0 and V = (V-1,...,V-d) satisfy Hormander's bracket condition, we demonstrate that Y-t admits a smooth density for any t is an element of (0, T], provided the driving noise satisfies certain nondegeneracy assumptions. Our analysis relies on relies on an interplay of rough path theory, Malliavin calculus and the theory of Gaussian processes. Our result applies to a broad range of examples including fractional Brownian motion with Hurst parameter H > 1/4, the Ornstein-Uhlenbeck process and the Brownian bridge returning after time T.