Covariance of the Kalman Filter innovation I am trying to fully understand the derivation of the covariance of the innovation vector, however I am stucked conceptually at a point. I will show you my reasoning and where I am stuck (if someone find something wrong I kindly ask for corrections):
Let's write the true measurement model and the measurement model after the state time prediction from epoch $k-1$ to $k$:
$z_k = H_k x_k + \epsilon_k~~~~~~~~~~\epsilon_k\sim N(0,R_k)$  
and
$z_k = H_k \hat{x}^-_k + \upsilon_k$  where
$\upsilon $ accounts for the $measurement$ and $prediction$ noise such that $\sim N(0,S_k)$
Therefore the KF innovation vector (assumed column) is the following:
$\gamma = z_k-H_k \hat{x}^-_k~~~~~~~~~~\gamma_k\sim N(0,S_k)$
In this sense, the covariance of $\gamma$ is expressed as:
$S_k= E[(\gamma - \underbrace{E[\gamma]}_0)(\gamma - \underbrace{E[\gamma]}_0)^T ]$
$S_k = E[\gamma\gamma^T] = E[(z_k-H_k \hat{x}^-_k)(z_k-H_k \hat{x}^-_k)^T]$ 
$S_k = E[z_kz_k^T-\color{blue}{z_k(H_k \hat{x}^-_k)^T}-\color{blue}{H_k \hat{x}^-_kz_k^T} +H_k {\hat{x}^-_k}{\hat{x}^-_k}^TH^{T}_k]$
Now here I don't really know how to get rid of the blue terms (if is possible). Because if I do, then I will have:
$S_k = E[z_kz_k^T] + H_kE[{\hat{x}^-_k}{\hat{x}^-_k}^T]H^{T}_k$
$S_k = H_kP^{-}_kH^{T}_k + R_k$ which is the final expression (rearranged)
Any opinion will be highly appreciated.
 A: Sorry to reply to this question so late.
The issue isn't actually about the cancellation of the terms, but rather the conceptual set up of the problem.  To re-frame the problem, let the state vector at time $k$ be given by $x_k$, and the observation model be given as
$y_k = H_k x_k + v_k$, 
where $v_k \sim N(0,R_k)$.  Remember, we are modelling the distribution of the random variable $x_k$, but are only given the measurements $y_k$ which we assume to be given by the above model.  Our mean state is one of the parameters which defines the distribution of $x_k$, let the mean state be defined to be $\widehat{x}_k$.  Now, the innovation at time $k$ is defined as
$\gamma_k = y_k - H_k \widehat{x}_k$;
note that our measurement of the mean is not artificially perturbed as in the definition of $z_k$ above.  In fact, in our estimation of the mean state, and its measurement, we always use the deterministic components of the model.  The covariance Riccati equation incorporates the model noise and the observational noise covariances, but this is only encoded in the mean state via the Kalman gain.
Now consider the covariance of $\gamma_k$ as defined above.  In particular, let us define the forecast error to be
$\epsilon_k = x_k - \widehat{x}_k$,
with covariance defined 
$P_k = \mathbb{E}[\epsilon_k \epsilon_k^{\text T}]$.
We may re-write the innovation $\gamma_k$ as
\begin{align}
\gamma_k &= H_k ( x_k - \widehat{x}_k) + v_k \\
& = H_k \epsilon_k + v_k
\end{align}.
Take the covariance of the above, using the independence and zero mean of each of the terms to recover the desired equation.
