Negative variances in Kalman smoother (FFBS)

I have implemented the forward-filtering-backwards-sampling (ffbs) algorithm. It consists of kalman filtering forward in time (to obtain mean and sigma). Then it uses these values and the Kalman smoother backwards in time to get mean and Sigma for all time points.

My problem is that the Kalman smoother seems to get negative variances. Can this be the case - and if so, how does it affect my results?

The code is made in R (and Sigma_predictions and Sigma_upates are obtained from the Kalman filter)

# ---- KALMAN / FORWARD FILTERING ----

# Kalman filter iterations

# Kalman loop starts
for(i in 1:big_t){ # Loop through each time point

# In the first iteration the inittial value is set
thisIsTheFirstValue <- (i == 1)
if(thisIsTheFirstValue){
phi_upd <- mean_0
sigma_upd <- sigma_0
}

# Prediction step
phi_pre <- (Gamma %*% phi_upd) + Theta
sigma_pre <- (Gamma %*% sigma_upd %*% t(Gamma)) + Psi_cov

# Update step
v <- y[, i] - (alpha + B %*% phi_pre) # y_t - f_t
big_F <- (B %*% sigma_pre %*% t(B)) + epsilon_cov # epsilon_cov = sigma * Idnetity matrix
inv_F <- solve(big_F, tol = 1e-40)
phi_upd <- phi_pre + (sigma_pre %*% t(B) %*% inv_F %*% v)
sigma_upd <- sigma_pre - (sigma_pre %*% t(B) %*% inv_F %*%  B %*% sigma_pre)

# save values in the arrays (used later in backward sampling step)
phis[, i] <- phi_upd
phi_predictions[ , , i] <- phi_pre
sigma_predictions[ , , i] <- sigma_pre
phi_updates[ , , i] <- phi_upd
sigma_updates[ , , i] <- sigma_upd

} # end for loop

# ---- BACKWARDS SAMPLING ----

# Find h_t and H_t by backwards sampling
loop_seq <- seq(from = (big_t - 1), to = 1, by = -1)
for(t in loop_seq){

# In the first iteration the inittial value is set
thisIsTheSecondLastTimePoint <- (t == (big_t-1))
if(thisIsTheSecondLastTimePoint){
phi_next <- matrix(rmvnorm(1, mean = phi_updates[ , , big_t], sigma = (sigma_updates[ , ,  big_t])), nrow = (p+1), ncol = 1)
# Save phi drawing
phi_array[ , , big_t] <- phi_next
}

# Calculate mean value dep. on this value
inv_sig_pred <- solve(sigma_predictions[ , , t+1])
mean_t <- phi_updates[, , t] + (sigma_updates[, , t] %*% t(Gamma) %*% inv_sig_pred %*% (phi_array[ , , (t+1)] - phi_predictions[ , , t+1]))
sigma_t <- sigma_updates[ , , t] - (sigma_updates[, , t] %*% t(Gamma) %*% inv_sig_pred %*% Gamma %*% sigma_updates[, , t])
# Make the drawing
phi_next <- t(rmvnorm(1, mean = mean_t, sigma = sigma_t))

# Save phi drawings
phi_array[, , t] <- phi_next
}


It is a state-space system, I am considering. lambda is set to 0.5. All is set due to an article, I am following.

The mathematical calculations are (sorry for a changed notation):

The problems, I am facing are problems with diagonal elements being negative in the H_t-matrix. The diagonal of the first calculated H_t I get: The diagonal of the second calculated H_t I get:

• Is t(Gamma) > 1? I don't know what your algorithm is supposed to do, so I don't know whether your problem is due to algorithm deficiency or implementation error. Is sigma_t actually a covariance matrix (not just 1 by 1 variance)? May 26, 2018 at 1:42
• What is t(Gamma)? Where is the closing parenthesis for the first character in: (sigma_updates[, , t] %*% t(Gamma) %*% May 26, 2018 at 1:53
• I think it is there (or am I wrong?). The line was just broken into two lines. It should be seen clearer now, as I set them to one. May 26, 2018 at 2:06
• Show us mathematically what the covariance update formulas are, and specifically what covariance is coming out with negative (variances?, meaning diagonal elements? or eigenvalues?). Are there some number of updates before the negative variances occur? May 26, 2018 at 2:09
• Hope it makes sense. It is already in the first iteration, that negative variances occur. May 26, 2018 at 2:53

Symmetrization: $H = 0.5*(H + H^T)$
Let V = eigenvector matrix of H, D = diagonal; matrix of eigenvalues of H. Therefore, $H = VDV^{-1}$. Replace diagonal elements of D below eigmin by eigmin, resulting in $D_{adjusrted}$. Now form $VD_{adjusted}V^{-1}$, and symmetrize if necessary.