Merge pull request #300 from LAMPSPUC/fix_missing_values

Fix missing values

Author: guilhermebodin

Project.toml

@@ -1,7 +1,7 @@
 name = "StateSpaceModels"
 uuid = "99342f36-827c-5390-97c9-d7f9ee765c78"
 authors = ["raphaelsaavedra <[email protected]>, guilhermebodin <[email protected]>, mariohsouto"]
-version = "0.5.20"
+version = "0.5.21"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

src/models/common.jl

@@ -14,3 +14,20 @@ function lagmat(y::Vector{Fl}, k::Int) where Fl
     end
     return X
 end
+
+function variance_of_valid_observations(y::Vector{Fl}) where Fl
+    return var(filter(!isnan, y))
+end
+
+function mean_of_valid_observations(y::Vector{Fl}) where Fl
+    return mean(filter(!isnan, y))
+end
+
+function assert_zero_missing_values(y::Vector)
+    for el in y
+        if isnan(el)
+            return error("This model does not support missing values.")
+        end
+    end
+    return nothing
+end

src/models/dar.jl

@@ -19,6 +19,9 @@ A Dynamic Autorregressive model is defined by:
     The dynamic autorregressive model does not have the [`forecast`](@ref) method implemented yet. 
     If you wish to perform forecasts with this model please open an issue.
 
+!!! warning Missing values
+    The dynamic autorregressive model currently does not support missing values (`NaN` observations.)
+
 # Example
 ```jldoctest
 julia> model = DAR(randn(100), 2)
@@ -34,6 +37,8 @@ mutable struct DAR <: StateSpaceModel
 
     function DAR(y::Vector{Fl}, lags::Int) where Fl
 
+        assert_zero_missing_values(y)
+
         X = lagmat(y, lags)
         num_observations = size(X, 1)
         first_observations = y[1:lags]
@@ -73,8 +78,8 @@ end
 
 function initial_hyperparameters!(model::DAR)
     Fl = typeof_model_elements(model)
-    observed_variance = var(model.system.y[findall(!isnan, model.system.y)])
-    observed_mean = mean(model.system.y[findall(!isnan, model.system.y)])
+    observed_variance = variance_of_valid_observations(model.system.y)
+    observed_mean = mean_of_valid_observations(model.system.y)
     initial_hyperparameters = Dict{String,Fl}(
         "sigma2_ε" => observed_variance, "intercept" => observed_mean
     )

src/models/exponential_smoothing.jl

@@ -135,8 +135,8 @@ end
 
 function initial_hyperparameters!(model::ExponentialSmoothing)
     Fl = typeof_model_elements(model)
-    observations = model.system.y[findall(!isnan, model.system.y)]
-    observed_variance = var(observations)
+    observations = filter(!isnan, model.system.y)
+    observed_variance = variance_of_valid_observations(model.system.y)
     initial_hyperparameters = Dict{String,Fl}(
         "sigma2" => observed_variance,
         "smoothing_level" => Fl(0.1),

src/models/locallevelcycle.jl

@@ -65,7 +65,7 @@ end
 
 function initial_hyperparameters!(model::LocalLevelCycle)
     Fl = typeof_model_elements(model)
-    observed_variance = var(model.system.y[findall(!isnan, model.system.y)])
+    observed_variance = variance_of_valid_observations(model.system.y)
     initial_hyperparameters = Dict{String,Fl}(
         "sigma2_ε" => observed_variance,
         "sigma2_η" => observed_variance,

src/models/locallevelexplanatory.jl

@@ -68,7 +68,7 @@ end
 
 function initial_hyperparameters!(model::LocalLevelExplanatory)
     Fl = typeof_model_elements(model)
-    observed_variance = var(model.system.y[findall(!isnan, model.system.y)])
+    observed_variance = variance_of_valid_observations(model.system.y)
     initial_hyperparameters = Dict{String,Fl}(
         "sigma2_ε" => observed_variance, "sigma2_η" => observed_variance
     )

src/models/locallineartrend.jl

@@ -50,16 +50,17 @@ end
 function default_filter(model::LocalLinearTrend)
     Fl = typeof_model_elements(model)
     steadystate_tol = Fl(1e-5)
-    a1 = zeros(Fl, 2)
-    P1 = Fl(1e6) .* Matrix{Fl}(I, 2, 2)
-    return UnivariateKalmanFilter(a1, P1, 2, steadystate_tol)
+    a1 = zeros(Fl, num_states(model))
+    P1 = Fl(1e6) .* Matrix{Fl}(I, num_states(model), num_states(model))
+    return UnivariateKalmanFilter(a1, P1, num_states(model), steadystate_tol)
 end
 
 function initial_hyperparameters!(model::LocalLinearTrend)
     Fl = typeof_model_elements(model)
+    observed_variance = variance_of_valid_observations(model.system.y)
     initial_hyperparameters = Dict{String,Fl}(
-        "sigma2_ε" => var(model.system.y),
-        "sigma2_ξ" => var(model.system.y),
+        "sigma2_ε" => observed_variance,
+        "sigma2_ξ" => observed_variance,
         "sigma2_ζ" => one(Fl),
     )
     set_initial_hyperparameters!(model, initial_hyperparameters)

src/models/naive_models.jl

@@ -13,7 +13,7 @@ typeof_model_elements(model::NaiveModel) = eltype(model.y)
 function assert_zero_missing_values(model::NaiveModel)
     for i in 1:length(model.y)
         if isnan(model.y[i])
-            return error("model $(typeof(model)) does not support missing values.)")
+            return error("model $(typeof(model)) does not support missing values.")
         end
     end
     return nothing

src/models/unobserved_components.jl

@@ -455,8 +455,7 @@ end
 
 function initial_hyperparameters!(model::UnobservedComponents)
     Fl = typeof_model_elements(model)
-    y = filter(!isnan, model.system.y)
-    observed_variance = var(y)
+    observed_variance = variance_of_valid_observations(model.system.y)
     # TODO add heuristic for initial hyperparameters
     initial_hyperparameters = Dict{String,Fl}(get_names(model) .=> one(Fl))
     if model.has_irregular

test/models/dar.jl

@@ -9,6 +9,8 @@
     # Runned on Python statsmodels
     @test loglike(model) ≈ -1175.9129 atol = 1e-5 rtol = 1e-5
 
+    @test_throws ErrorException DAR(vcat(rand(10), NaN, rand(10)), 2)
+
     # forecasting
     @test_throws ErrorException forecast(model, 10)
 end