From 1dfc0c25cd3fce52657ae3ce785bd8e8f481b69b Mon Sep 17 00:00:00 2001
From: Marina Dietze <87202177+dietzemarina@users.noreply.github.com>
Date: Thu, 4 Aug 2022 10:56:11 -0300
Subject: [PATCH] fixed cholesky error in simulate (#314)

* fixed cholesky error in simulate

* created cholesky_decomposition function

* increase rounding digits

* updated version in Project.toml
---
 Project.toml                              |  2 +-
 src/models/basicstructural_explanatory.jl |  2 +-
 src/systems.jl                            | 21 +++++++++++++++++++--
 3 files changed, 21 insertions(+), 4 deletions(-)

diff --git a/Project.toml b/Project.toml
index b848e2a..cdb52d5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "StateSpaceModels"
 uuid = "99342f36-827c-5390-97c9-d7f9ee765c78"
 authors = ["raphaelsaavedra <raphael.saavedra93@gmail.com>, guilhermebodin <guilherme.b.moraes@gmail.com>, mariohsouto"]
-version = "0.6.4"
+version = "0.6.5"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
diff --git a/src/models/basicstructural_explanatory.jl b/src/models/basicstructural_explanatory.jl
index b0136ea..ca73ca1 100644
--- a/src/models/basicstructural_explanatory.jl
+++ b/src/models/basicstructural_explanatory.jl
@@ -192,7 +192,7 @@ function simulate(
     alpha = Matrix{Fl}(undef, n + 1, m)
     # Sampling errors
     chol_H = sqrt(sys.H[1])
-    chol_Q = cholesky(sys.Q[1])
+    chol_Q = cholesky_decomposition(sys.Q[1])
     standard_ε = randn(n)
     standard_η = randn(n + 1, size(sys.Q[1], 1))
 
diff --git a/src/systems.jl b/src/systems.jl
index 2431c48..cf3307e 100644
--- a/src/systems.jl
+++ b/src/systems.jl
@@ -312,6 +312,23 @@ function to_multivariate_time_variant(system::LinearUnivariateTimeVariant{Fl}) w
 end
 
 to_multivariate_time_variant(system::LinearMultivariateTimeVariant) = system
+ispossemdef(M::Matrix{Fl}) where Fl = all(i -> i >= 0.0, eigvals(M))
+
+function cholesky_decomposition(M::Matrix{Fl}) where Fl
+    if isposdef(M)
+        return cholesky(M)
+    elseif ispossemdef(M)
+        size_matrix = size(M, 1)
+        chol_M = cholesky(M .+ I(size_matrix) .* floatmin(Fl))
+        chol_M.L[:, :] = round.(chol_M.L; digits = 10)
+        chol_M.U[:, :] = round.(chol_M.U; digits = 10)
+        chol_M.UL[:, :] = round.(chol_M.UL; digits = 10)
+
+        return chol_M
+    else
+        @error("Matrix is not positive definite or semidefinite. Cholesky decomposition cannot be performed.")
+    end
+end
 
 # Functions for simulations
 function simulate(
@@ -325,7 +342,7 @@ function simulate(
     alpha = Matrix{Fl}(undef, n + 1, m)
     # Sampling errors
     chol_H = sqrt(sys.H)
-    chol_Q = cholesky(sys.Q)
+    chol_Q = cholesky_decomposition(sys.Q)
     standard_ε = randn(n)
     standard_η = randn(n + 2, size(sys.Q, 1))
 
@@ -357,7 +374,7 @@ function simulate(
     alpha = Matrix{Fl}(undef, n + 1, m)
     # Sampling errors
     chol_H = cholesky(sys.H)
-    chol_Q = cholesky(sys.Q)
+    chol_Q = cholesky_decomposition(sys.Q)
     standard_ε = randn(n, size(sys.H, 1))
     standard_η = randn(n + 2, size(sys.Q, 1))