# 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?

• In a Kalman filter, you typically just assume that they are. Tests are available for multivariate normality, though I don't know of one that takes as given that each marginal is Gaussian but tests for the joint. – Dougal Jul 15 '15 at 23:58
• One can't prove that things are bivariate Gaussian. Sometimes we can see variables aren't consistent with it, but looking consistent with it doesn't mean it that it is the case; ultimately it's at best a plausible assumption to which (if we're lucky) our inferential procedures may not be especially sensitive. Not that failure to reject a hypothesis of bivariate Gaussian doesn't imply that it is actually the case, only that our sample size wasn't large enough to detect a difference from it with the test we used. – Glen_b -Reinstate Monica Jul 16 '15 at 2:34

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.