For a Kalman filter you are not so much interested in the "stability" of $\hat{x}$ (full state estimation), but in the error between $x$ (the actual state) and $\hat{x}$. Because if this error goes to zero, then $\hat{x}$ will become equal to $x$, which is what we want.
The dynamics of $\hat{x}$ can be written as,
$$
\dot{\hat{x}} = A\,\hat{x} + B\,u + K\,(y - C\,\hat{x} - D\,u), \tag{1}
$$
Defining the error as $e = x - \hat{x}$, and thus its derivative as $\dot{e} = \dot{x} - \dot{\hat{x}}$. Substituting in the expressions for $\dot{x}$, $\dot{\hat{x}}$ and $y$ yields,
$$
\dot{e} = \left(A\,x + B\,u\right) - \left(A\,\hat{x} + B\,u + K\,((C\,x + D\,u) - C\,\hat{x} - D\,u)\right) = \left(A - K\,C\right)e. \tag{2}
$$
So if $A - K\,C$ is a Hurwitz matrix, then the magnitude of $e$ will go down in time and $\hat{x}$ will keep becoming a better estimation of $x$.
You are correct that $K=P\,C^T(D\,R\,D^T)^{-1}$. Now if you define the dynamics of some new variable $\zeta$ as,
$$
\dot{\zeta} = \left(A - P\,C^T(D\,R\,D^T)^{-1}C\right)^T \zeta, \tag{3}
$$
then the stability of this dynamics will be the same as for $e$, since transposing does not change the eigenvalues of a matrix. A quadratic Lyapunov for this system can shown to be,
$$
V(\zeta) = \zeta^T P\,\zeta, \tag{4}
$$
since its time derivative equals,
$$
\dot{V}(\zeta) = \zeta^T P\,\dot{\zeta} + \dot{\zeta}^T P\,\zeta = \zeta^T \left(P\,A^T + A\,P - 2\,P\,C^T(D\,R\,D^T)^{-1}C\,P\right) \zeta. \tag{5}
$$
Using that $P\,A^T + A\,P - \,P\,C^T(D\,R\,D^T)^{-1}C\,P + B^T B = 0$, then
$$
P\,A^T + A\,P - 2\,P\,C^T(D\,R\,D^T)^{-1}C\,P = -B^T B - \,P\,C^T(D\,R\,D^T)^{-1}C\,P, \tag{6}
$$
but $B^T B$ and $P\,C^T(D\,R\,D^T)^{-1}C\,P$ will be (semi-)positive definite and a sum of a positive definite matrix and a semi-positive definite matrix will also yield a positive definite matrix. So the left hand side of equation $(6)$ will be a negative definite matrix and therefore $V(\zeta)$ will be quadratic Lyapunov of equation $(3)$, which implies that the corresponding matrix must be Hurwitz.
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.
Best Answer
A short stability proof of the Kalman filter is given in section 2.4.1 of:
I will give a short summary of what is stated in this dissertation. First of all it is assumed that $Q$ and $R$ (these two matrices can also be time varying) are such that $(F_k,Q)$ is uniformly controllable and $(F_k,H_k^\top R^{-1} H_k)$ is uniformly observable, then the error covariance matrix $P_k$ has a finite positive definite upper and lower bound.
The a priori estimation error $\tilde{x}_k = \hat{x}_k - x_k$ can shown to be stable using the following Lyapunov function
$$ V(\tilde{x}_k) = \tilde{x}_k^\top P_k^{-1} \tilde{x}_k, $$
which is guaranteed to be a positive definite function when the controllability and observability conditions are satisfied. The estimation error dynamics is given by
$$ \tilde{x}_{k+1} = F_k\,(I - K_k\,H_k)\,\tilde{x}_{k} + F_k\,K_k\,v_k - w_k. $$
For the bounded input bounded output (BIBO) stability analysis it is only required to consider the homogenous part of the error dynamics. Doing so yields to following increment of the proposed Lyapunov function
\begin{align} \Delta V(\tilde{x}_k) &= \tilde{x}_{k+1}^\top P_{k+1}^{-1} \tilde{x}_{k+1} - \tilde{x}_k^\top P_k^{-1} \tilde{x}_k, \\ &= \tilde{x}_k^\top \left[(I - K_k\,H_k)^\top F_k^\top P_{k+1}^{-1} F_k (I - K_k\,H_k) - P_k^{-1}\right] \tilde{x}_k, \end{align}
which after some more manipulation can shown to be negative definite. Previously it was shown that $P_k$ is upper and lower bounded and therefore $P_k^{-1}$ should be as well, which together with $\Delta V(\tilde{x}_k)\prec 0$ enables one to show that the Kalman filter is BIBO stable.
For showing that the Kalman filter in stable one only needs that $w_k$ and $v_k$ are bounded, so nowhere is it required that $Q$ and and $R$ are such that $w_k \sim \mathcal{N}(0, Q)$ and $v_k \sim \mathcal{N}(0, R)$. So the stochastic terms can have very different covariance matrices associated with them. Since the Kalman filter is dealing with a stochastic system an more in-depth analysis is needed to show that covariance of the error remains bounded when there is a mismatch between noise covariance matrices the actual noise compared to the ones used in the Kalman filter. Further resources regarding this are mentioned earlier mentioned dissertation.
However, having a mismatch in the noise covariance matrices can result in a far from optimal noise suppression, as is discussed in: