The discounted stock price satisfies
$$
dX(t) = \big(\mu(t) - r(t)\big)X(t) dt + \sigma_S(t) X(t) dW_S^{\mathbb P}(t).
$$
The Girsanov density for $X$ is
$$
Z(T) = \exp\left\{\int_0^T \nu(t)dW^{\mathbb P}_S(t) - \frac12\int_0^T \nu(t)^2 dt \right\}
$$
with $\nu(t) = (r(t)-\mu(t))/\sigma_S(t)$. The problem, as usual, is to establish that $\mathbb{E}[Z(t)]=1$. Since $r$ has Gaussian distribution, it is not possible to proceed with Novikov or Kazamaki conditions, as usual. However, when $\mu$, $\theta$ and $\sigma_S$ are reasonable enough, it is possible to use the less known Beneŝ condition as follows. The solution to the equation for interest rate is
$$
r(t) = r(0)e^{-at} + \int_0^t e^{a(s-t)}\big(\theta(s) ds + \sigma_r dW^{\mathbb{P}}_r(s)\big) \\
= r(0)e^{-at} + \int_0^t e^{a(s-t)}\theta(s) ds + \sigma_r \int_0^t e^{a(s-t)} d\big(\rho W^{\mathbb{P}}_S(s) +\gamma(s)\big) \\
= \Gamma(t) + \Lambda(t),
$$
where $$\Gamma(t) = r(0)e^{-at} +\int_0^t e^{a(s-t)}\theta(s) ds + \int_0^t e^{a(t-s)}d\gamma(s),\\
\Lambda(t) = \sigma_r \rho\left( W^{\mathbb{P}}_S(t) + a \int_0^t e^{a(t-s)}W^{\mathbb{P}}_S(s)ds\right),
$$ and $\gamma$ is, up to a constant, a Wiener process independent of $W^{\mathbb{P}}_S$. Then, assuming $\mu$, $\theta$ and $\sigma_S$ to be nice (say, all bounded, and $\inf|\sigma_S|>0$), we have $|\nu(t)|\le C(\gamma)+ C\max_{s\in [0,t]} |W_S^{\mathbb P}(s)|$, where $C(\gamma)$ is a finite random variable depending on $\gamma$. Then we have, see e.g. here,
$$
\mathbb{E}\left[Z(T) \mid \gamma\right]=1,\tag{1}
$$
whence $\mathbb{E}\left[Z(T)\right]=1$, as required.
As @michael wrote, there is no uniqueness (unless $|\rho| = 1$ or $\sigma_r = 0$). Financially, the explanation is that the interest rate process is non-traded, thus there is no completeness. Mathematically, in view of $(1)$, for any $\gamma$-measurable positive random variable $L(\gamma)$, which integrates to $1$, $Z'(T) = L(\gamma) Z(T)$ defines a density of equivalent martingale measure.
