3
$\begingroup$

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)
  # Calculate updates
  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. enter image description here

The mathematical calculations are (sorry for a changed notation): enter image description here enter image description here

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: enter image description here The diagonal of the second calculated H_t I get: enter image description here

$\endgroup$
  • $\begingroup$ 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)? $\endgroup$ – Mark L. Stone May 26 '18 at 1:42
  • $\begingroup$ What is t(Gamma)? Where is the closing parenthesis for the first character in: (sigma_updates[, , t] %*% t(Gamma) %*% $\endgroup$ – Mark L. Stone May 26 '18 at 1:53
  • $\begingroup$ 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. $\endgroup$ – Dorthe May 26 '18 at 2:06
  • $\begingroup$ 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? $\endgroup$ – Mark L. Stone May 26 '18 at 2:09
  • $\begingroup$ Hope it makes sense. It is already in the first iteration, that negative variances occur. $\endgroup$ – Dorthe May 26 '18 at 2:53
1
$\begingroup$

You appear to be getting roundoff error level negative variances. Subtraction is a dangerous thing. There is a perhaps a better "root cause" fixing up which can be done, by modifying the updating equation.

But in the absence of that, one quick "fix", which I'm not saying is a good thing to do here, is to adjust the covariance matrix on each iteration to be symmetric positive semidefinite, with some minimum eigenvalue value, eigmin (a very small positive number). Do this by first symmetrizing, then eigenvalue eigenvector decomposiition, then replacing eigenvalues below eigmin with eigmin, then reconstituting the matrix. Symmetrize if needed at the end.

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.

$\endgroup$
  • $\begingroup$ @Dorthe Does my answer help you or is it off the mark? $\endgroup$ – Mark L. Stone May 30 '18 at 0:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.