add error for stochastic singularity

DynareJulia · Jan 6, 2024 · 2a6aa16 · 2a6aa16
1 parent 741d8c9
commit 2a6aa16
Show file tree

Hide file tree

Showing 2 changed files with 11 additions and 2 deletions.
diff --git a/src/kalman_filter.jl b/src/kalman_filter.jl
@@ -147,6 +147,7 @@ function kalman_filter!(Y::AbstractArray{X},
     vP = view(P, :, :, 1)
     copy!(ws.oldP, vP)
     cholHset = false
+    Falwayssingular = true
     while t <= last
         pattern = data_pattern[t]
         ndata = length(pattern)
@@ -194,13 +195,17 @@ function kalman_filter!(Y::AbstractArray{X},
         info = get_cholF!(vcholF, vF)
         if rcond(vcholF) < ws.kalman_tol || info != 0
             # F is near singular
+            #=
             if !cholHset
                 get_cholF!(vcholH, vvH)
                 cholHset = true
             end
             ws.lik[t] = ndata*l2pi + extended_univariate_step!(vatt, va1, vPtt, vP1, Y, t, c, ws.Zsmall, vvH, d, T, ws.QQ, va, vP, ws.kalman_tol, ws, pattern)
             t += 1
             continue
+            =#
+        else
+            Falwayssingular = false
         end 
         # iFv = inv(F)*v
         get_iFv!(viFv, vcholF, vv)
@@ -214,6 +219,9 @@ function kalman_filter!(Y::AbstractArray{X},
         end
         t += 1
     end
+    if Falwayssingular
+        error("The covariance matrix of the one-step ahead forecast error for the observed variables (F) is singular. This a case of stochastic singularity.")
+    end
     vlik = view(ws.lik, start + presample:last)
     return -0.5*sum(vlik)
 end

diff --git a/src/kalman_smoother.jl b/src/kalman_smoother.jl
@@ -143,6 +143,7 @@ function kalman_smoother!(Y::AbstractArray{V},
         #        ws.csmall .= view(vc, pattern)
         vZsmall = get_vZsmall(ws.Zsmall, ws.iZsmall, Z, pattern, ndata, ny, t)
         vH = changeH ? view(H, pattern, pattern, t) : view(H, pattern, pattern)
+        hasmeasurementerror = !all(vH .== 0)
         vT = changeT ? view(T, :, :, t) : view(T, :, :)
         va = changeA ? view(a, :, t) : view(a, :)
         vatt = changeAtt ? view(att, :, t) : view(att, :)
@@ -163,13 +164,13 @@ function kalman_smoother!(Y::AbstractArray{V},
         # r_{t-1} = Z_t'*iF_t*v_t + L_t'r_t (DK 4.44)
         update_r!(ws.r, vZsmall, viFv, ws.L, ws.r_1)
         if (length(alphah) > 0 ||
-            length(epsilonh) > 0 ||
+            hasmeasurementerror ||
             length(etah) > 0)
             # N_{t-1} = Z_t'iF_t*Z_t + L_t'N_t*L_t (DK 4.44)
             get_iFZ!(viFZ, vcholF, vZsmall)
             update_N!(ws.N, vZsmall, viFZ, ws.L, ws.N_1, ws.PTmp)
         end
-        if length(epsilonh) > 0
+        if hasmeasurementerror
             vepsilonh = view(epsilonh, pattern, t)
             # epsilon_t = H*(iF_t*v_t - KDK_t'*r_t) (DK 4.69)
             vtmp1 = view(ws.tmp_ny, 1:ndata)