Determining if two correlated Gaussian RVs are jointly Gaussian One common way we can find the posterior distribution of Gaussian parameters for a Kalman filter on Gaussian observations is by first computing the covariance between the parameters and forecast given data until the previous time point and then conditioning on the observed value. In other words, for the model
$$\begin{align*}
y_t &= \mathbf{F}_t'\boldsymbol{\theta}_t + \nu_t, \quad \nu_t \sim \mathcal{N}(0,v_t)\\
\boldsymbol{\theta}_t &= \boldsymbol{\theta}_{t-1} + \boldsymbol{\omega}_t, \quad \boldsymbol{\omega}_t \sim \mathcal{N}(0, \mathbf{W}_t)
\end{align*}$$
we can compute the covariance between $\boldsymbol{\theta}_t$ and $y_t$ given data up to time $t-1$ and then condition on an observation $y_t^*$ to update the distribution $\boldsymbol{\theta}_t$ given data up to time $t$.
As shown here, in general, two correlated variables that are marginally Gaussian may not be jointly Gaussian. My question is, how do we know that the joint distribution of $\boldsymbol{\theta}_t$ and $y_t$ is Gaussian? Unless we establish this, we can't simply condition on $y_t^*$ and use the conditional Gaussian equations to update the distribution of $\boldsymbol{\theta}_t$.
Although it's possible to come up with counterexamples to show that two marginally Gaussian variables that are correlated may not be jointly Gaussian, in general, how can you know when they are?
 A: For the standard Kalman filter we must assume $\boldsymbol{\theta}_t$ and $y_t$ are jointly Gaussian.  So in short, we know they are joint Gaussian because we assumed it in the first place.  This assumption is made for convenience in estimation and I do not see any obvious reason why you would want to relax this assumption given a linear, Gaussian, state space model.
The more general question you pose, given two Gaussian distributed variables, how do I test that they are distributed joint Gaussian?, is a more interesting question.  In this setting, it is generally easier to specify a particular class of alternative distributions and compare model fit rather than conduct formal hypothesis testing.  For example, let $X_1,..X_n \stackrel{iid}{\sim}N(\mu_x,\sigma_x)$ and $Y_1,..,Y_n\stackrel{iid}{\sim}N(\mu_y,\sigma_y);$ it is then much easier to compare the fit of a bivariate normal model to the fit of, say  a bivariate Clayton copula then to test the hypothesis
$H_0:$ $(\mathbf{X,Y})$ is bivariate normal vs $H_a:$ $(\mathbf{X,Y})$ is not bivariate normal
However, you may still be interested in such a hypothesis.  In this case, residual diagnostic testing can be used to asses the validity of the joint normality assumption given a set of parameters, which comes very close to the above hypothesis test.  
To do this we begin with a likelihood.  If we assume $(\mathbf{X,Y})$ are bivariate normal than we can estimate the univariate parameters $\mu_x,\mu_y,\sigma_x,\sigma_y$ and the correlation coefficient, $\rho$, via maximum likelihood.
$$
L(\mu_x,\mu_y,\sigma_x,\sigma_y,\rho|\mathbf{X,Y})=\prod_{i=1}^n\boldsymbol{\phi}(X_i,Y_i|\mu_x,\mu_y,\sigma_x,\sigma_y,\rho)
$$ 
where $\boldsymbol{\phi}$ is the bivariate normal pdf.  
Given the maximum likelihood estimates $\hat \mu_x,\hat\mu_y,\hat\sigma_x,\hat\sigma_y,\hat\rho$, we can generate a generalized residual for each $i$ by first calculating $u_i=\boldsymbol{\Phi}(X_i,Y_i|\hat \mu_x,\hat\mu_y,\hat\sigma_x,\hat\sigma_y,\hat\rho)$ where $\boldsymbol{\Phi}$ is the bivariate normal cdf.  By the probability integral transform, it should be the case and that $u_i \sim uniform(0,1)$ and in turn, the generalized residual $r_i=\Phi^{-1}(u_i)$, where $\Phi$ is the univariate standard normal cdf, should be distributed standard normal.
You can then test $r_1,..r_n$ for standard normality using the KS-test, Jarque-Bera test or whatever else you may prefer.  The idea is that if  $(\mathbf{X,Y})$ are truly bivariate normal with correctly estimated parameters, then $r_1,..r_n$ should be distributed standard normal.  If we reject normality in $r_1,..r_n$, than we in-effect reject that $(X_i,Y_i) \sim \mathcal{N}(\hat \mu_x,\hat\mu_y,\hat\sigma_x,\hat\sigma_y,\hat\rho)$. 
One nuance is that the hypothesis test is conditional on the estimated parameters so we would not be testing the "unconditional" hypothesis
$H_0:$ $(\mathbf{X,Y})$ is bivariate normal vs $H_a:$ $(\mathbf{X,Y})$ is not bivariate normal  
exactly.  In a Bayesian setting one could probably find a way to integrate over parameter uncertainty, but given a sufficient amount of data the difference between the conditional and unconditional test becomes negligible.
