Consider at dynamic linear model where $$ \theta_{1} \sim N(\mu_{1}, W_{1}), $$ $$ \theta_{i}=G\theta_{i-1} + w_{i}, w_{i}\sim N(0,W), $$ $$ Y_{i} = F\theta_{i} + v_{i}, v_{i}\sim N(0,V) $$ and $ \theta_{1}, w_{i}, v_{i} $ all independent random vectors. Let $ \theta_{0:t} : = (\theta_{t}, \theta_{t-1},\ldots, \theta_{0}) $ and $ Y_{1:t}:= (Y_{t},Y_{t-1},\ldots, Y_{1})$.
A generel result from multivariate normal distribution is as follows: Suppose $ Y|X=x \sim N(Ax,V)$ and $ X \sim N(\mu, \Sigma) $. Then $ Y\sim(A\mu, A\Sigma A^{T} + V) $.
Question 1: How do I see, that the joint density of $ (\theta_{0:t}, Y_{1:t}) $ is Gaussian? From the definitions above, we have $Y_{t}|\theta_{t}$ is normal and $\theta_{t}|\theta_{t-1}$ is normal. From the general result the joint densities of $(Y_{t},\theta_{t})$ and $(\theta_t, \theta_{t-1})$ is also normal. In fact, since a dynamic linear model is also a state space model, it can be shown (and I have) that $$p(\theta_{0:t}, y_{1:t}) = p(\theta_0)\prod_{j=1}^{t}p(y_j | \theta_j)p(\theta_j | \theta_{j-1}).$$ But I can not complete the puzzle; i.e why is $ (\theta_{0:t}, Y_{1:t}) $ is Gaussian?
Question 2: Assume we know $$ \theta_{n-1}\mid Y_{1:n-1} = y_{1:n-1} \sim N(\hat \theta_{n-1}, \Sigma_{n-1}).$$ Now, how do I apply the generel result stated above, to see that $$ \theta_{n}\mid Y_{1:n-1} = y_{1:n-1} \sim N(G\hat \theta_{n-1}, G\Sigma_{n-1}G^{T} + W).$$ Of course we have $E[\theta_n | Y_{1:n-1}] = G\hat \theta_{n-1}$ and $Var[\theta_n | Y_{1:n-1}] = G\Sigma_{n-1}G^{T} + W$. But I am not sure why $\theta_{n}\mid Y_{1:n-1}$ is normally distributed. My guess is, that I can apply the result about linear combinations of Gaussians is again Gaussian.
I have come up with an answer myself.
We know that $ y_{t}\mid \theta_{t} $ is Gaussian, and since $ \theta_{t} $ is Gaussian the joint density of $ (y_{t},\theta_{t}) $ is also Gaussian (By general result). By conditional independence $ \theta_{t}\mid y_{t-1}, \theta_{t-1} $ is Gaussian (it has same distribution as $\theta_{t}\mid \theta_{t-1}$), and therefore the joint density of $ (\theta_{t}, y_{t-1}, \theta_{t-1}) $ is Gaussian. Proceeding this way, the joint density of $ (\theta_{0:t}, y_{1:t}) $ is seen to be Gaussian. This implies that every sub vector and conditional distributions are Gaussian.