fix errors in unit tests and add in sample metrics

LAMPSPUC · Jul 22, 2024 · 96e6d5b · 96e6d5b
1 parent 1de5c12
commit 96e6d5b
Show file tree

Hide file tree

Showing 6 changed files with 91 additions and 59 deletions.
diff --git a/src/utils.jl b/src/utils.jl
@@ -515,3 +515,20 @@ function get_seasonality_dict_and_stochastic(seasonality::Vector{String})
 
     return seasonality_dict, any(stochastic), stochastic
 end
+
+function get_in_sample_statistics(output::Output, y::Vector{Fl}) where Fl
+
+    y_pred = output.fit_in_sample
+    aic    = output.information_criteria["aic"]
+    bic    = output.information_criteria["bic"]
+    aicc   = output.information_criteria["aicc"]
+
+    mse  = mean((y .- y_pred).^2)
+    rmse = sqrt(mse)
+    mae  = mean(abs.(y .- y_pred))
+    mape = mean(abs.((y .- y_pred) ./ y)) * 100
+
+    return Dict("MSE" => mse, "RMSE" => rmse, 
+                "MAE" => mae, "MAPE" => mape,
+                "BIC"  => bic, "AIC"  => aic, "AICc" => aicc)
+end
diff --git a/test/test_fit_forecast_lognormal.jl b/test/test_fit_forecast_lognormal.jl
@@ -50,14 +50,22 @@
     gas_model_lognormal_X         = deepcopy(gas_model_lognormal)
     gas_model_lognormal_X_2params = deepcopy(gas_model_lognormal_2params)
 
+    model_lognormal, parameters_lognormal, initial_values_lognormal                         = UnobservedComponentsGAS.create_model(gas_model_lognormal, y, missing)
+    model_lognormal_2params, parameters_lognormal_2params, initial_values_lognormal_2params = UnobservedComponentsGAS.create_model(gas_model_lognormal_2params, y, missing)
+
+    model_lognormal_X, parameters_lognormal_X, initial_values_lognormal_X                         = UnobservedComponentsGAS.create_model(gas_model_lognormal_X, y, X, missing);
+    model_lognormal_X_2params, parameters_lognormal_X_2params, initial_values_lognormal_X_2params = UnobservedComponentsGAS.create_model(gas_model_lognormal_X_2params, y, X, missing);
+
+    fitted_model_lognormal         = UnobservedComponentsGAS.fit(gas_model_lognormal, y; tol = 5e-2)
+    fitted_model_lognormal_2params = UnobservedComponentsGAS.fit(gas_model_lognormal_2params, y; tol = 5e-2)
+    fitted_model_lognormal_X         = UnobservedComponentsGAS.fit(gas_model_lognormal_X, y, X)
+    fitted_model_lognormal_X_2params = UnobservedComponentsGAS.fit(gas_model_lognormal_X_2params, y, X)
+
+    forecast_lognormal         = UnobservedComponentsGAS.predict(gas_model_lognormal, fitted_model_lognormal, y, steps_ahead, num_scenarious)
+    #forecast_lognormal_X       = UnobservedComponentsGAS.predict(gas_model_lognormal_X, fitted_model_lognormal_X, y, X_lognormal_forec, steps_ahead, num_scenarious)
+    forecast_lognormal_2params = UnobservedComponentsGAS.predict(gas_model_lognormal_2params, fitted_model_lognormal_2params, y, steps_ahead, num_scenarious)
+
     @testset " --- Testing create_model functions" begin
-        # Create model with no explanatory series
-        model_lognormal, parameters_lognormal, initial_values_lognormal                         = UnobservedComponentsGAS.create_model(gas_model_lognormal, y, missing)
-        model_lognormal_2params, parameters_lognormal_2params, initial_values_lognormal_2params = UnobservedComponentsGAS.create_model(gas_model_lognormal_2params, y, missing)
-
-        model_lognormal_X, parameters_lognormal_X, initial_values_lognormal_X                         = UnobservedComponentsGAS.create_model(gas_model_lognormal_X, y, X, missing);
-        model_lognormal_X_2params, parameters_lognormal_X_2params, initial_values_lognormal_X_2params = UnobservedComponentsGAS.create_model(gas_model_lognormal_X_2params, y, X, missing);
-
         @test(size(parameters_lognormal)         == (T,2))
         @test(size(parameters_lognormal_2params) == (T,2))
         @test(typeof(model_lognormal)            == JuMP.Model)
@@ -73,12 +81,7 @@
         @test(test_initial_values_components(initial_values_lognormal_X_2params, rw, rws, ar, seasonality))
     end
 
-    @testset " --- Testing fit functions" begin
-        fitted_model_lognormal         = UnobservedComponentsGAS.fit(gas_model_lognormal, y; tol = 5e-2)
-        fitted_model_lognormal_2params = UnobservedComponentsGAS.fit(gas_model_lognormal_2params, y; tol = 5e-2)
-        fitted_model_lognormal_X         = UnobservedComponentsGAS.fit(gas_model_lognormal_X, y, X)
-        fitted_model_lognormal_X_2params = UnobservedComponentsGAS.fit(gas_model_lognormal_X_2params, y, X)
-
+    @testset " --- Testing fit functions" begin        
         # "Test if termination_status is correct"
         possible_status = ["LOCALLY_SOLVED", "TIME_LIMIT"]
         @test(fitted_model_lognormal.model_status in possible_status)
@@ -108,11 +111,6 @@
     end
 
     @testset " --- Test forecast function ---" begin
-        forecast_lognormal         = UnobservedComponentsGAS.predict(gas_model_lognormal, fitted_model_lognormal, y, steps_ahead, num_scenarious)
-        #forecast_lognormal_X       = UnobservedComponentsGAS.predict(gas_model_lognormal_X, fitted_model_lognormal_X, y, X_lognormal_forec, steps_ahead, num_scenarious)
-        forecast_lognormal_2params = UnobservedComponentsGAS.predict(gas_model_lognormal_2params, fitted_model_lognormal_2params, y, steps_ahead, num_scenarious)
-
-
         @test(isapprox(forecast_lognormal["mean"], vec(mean(forecast_lognormal["scenarios"], dims = 2)); rtol = 1e-3)) 
         @test(size(forecast_lognormal["scenarios"]) == (steps_ahead, num_scenarious))
 

diff --git a/test/test_fit_forecast_normal.jl b/test/test_fit_forecast_normal.jl
@@ -49,14 +49,24 @@
 
     gas_model_normal_X            = deepcopy(gas_model_normal)
     gas_model_normal_X_2params    = deepcopy(gas_model_normal_2params)
+
+    model_normal, parameters_normal, initial_values_normal                         = UnobservedComponentsGAS.create_model(gas_model_normal, y, missing)
+    model_normal_2params, parameters_normal_2params, initial_values_normal_2params = UnobservedComponentsGAS.create_model(gas_model_normal_2params, y, missing)
+
+    model_normal_X, parameters_normal_X, initial_values_normal_X                         = UnobservedComponentsGAS.create_model(gas_model_normal_X, y, X, missing);
+    model_normal_X_2params, parameters_normal_X_2params, initial_values_normal_X_2params = UnobservedComponentsGAS.create_model(gas_model_normal_X_2params, y, X, missing);
 
+    fitted_model_normal         = UnobservedComponentsGAS.fit(gas_model_normal, y; tol = 5e-2)
+    fitted_model_normal_2params = UnobservedComponentsGAS.fit(gas_model_normal_2params, y; tol = 5e-2)
+    fitted_model_normal_X         = UnobservedComponentsGAS.fit(gas_model_normal_X, y, X; tol = 5e-2)
+    fitted_model_normal_X_2params = UnobservedComponentsGAS.fit(gas_model_normal_X_2params, y, X; tol = 5e-2)
+
+    forecast_normal         = UnobservedComponentsGAS.predict(gas_model_normal, fitted_model_normal, y, steps_ahead, num_scenarious)
+    forecast_normal_X       = UnobservedComponentsGAS.predict(gas_model_normal_X, fitted_model_normal_X, y, X_normal_forec, steps_ahead, num_scenarious)
+    forecast_normal_2params = UnobservedComponentsGAS.predict(gas_model_normal_2params, fitted_model_normal_2params, y, steps_ahead, num_scenarious)
+
     @testset "create_model" begin
     # Create model with no explanatory series
-        model_normal, parameters_normal, initial_values_normal                         = UnobservedComponentsGAS.create_model(gas_model_normal, y, missing)
-        model_normal_2params, parameters_normal_2params, initial_values_normal_2params = UnobservedComponentsGAS.create_model(gas_model_normal_2params, y, missing)
-
-        model_normal_X, parameters_normal_X, initial_values_normal_X                         = UnobservedComponentsGAS.create_model(gas_model_normal_X, y, X, missing);
-        model_normal_X_2params, parameters_normal_X_2params, initial_values_normal_X_2params = UnobservedComponentsGAS.create_model(gas_model_normal_X_2params, y, X, missing);
 
         @test(size(parameters_normal)         == (T,2))
         @test(size(parameters_normal_2params) == (T,2))
@@ -75,11 +85,6 @@
 
     #@info(" --- Testing fit functions")
     @testset "fit" begin
-        fitted_model_normal         = UnobservedComponentsGAS.fit(gas_model_normal, y; tol = 5e-2)
-        fitted_model_normal_2params = UnobservedComponentsGAS.fit(gas_model_normal_2params, y; tol = 5e-2)
-        fitted_model_normal_X         = UnobservedComponentsGAS.fit(gas_model_normal_X, y, X; tol = 5e-2)
-        fitted_model_normal_X_2params = UnobservedComponentsGAS.fit(gas_model_normal_X_2params, y, X; tol = 5e-2)
-
         # "Test if termination_status is correct"
         possible_status = ["LOCALLY_SOLVED", "TIME_LIMIT"]
         @test(fitted_model_normal.model_status in possible_status)
@@ -110,11 +115,7 @@
 
     #@info(" --- Test forecast function ---")
     @testset "forecast" begin
-        forecast_normal         = UnobservedComponentsGAS.predict(gas_model_normal, fitted_model_normal, y, steps_ahead, num_scenarious)
-        forecast_normal_X       = UnobservedComponentsGAS.predict(gas_model_normal_X, fitted_model_normal_X, y, X_normal_forec, steps_ahead, num_scenarious)
-        forecast_normal_2params = UnobservedComponentsGAS.predict(gas_model_normal_2params, fitted_model_normal_2params, y, steps_ahead, num_scenarious)
-
-
+
         @test(isapprox(forecast_normal["mean"], vec(mean(forecast_normal["scenarios"], dims = 2)); rtol = 1e-3)) 
         @test(size(forecast_normal["scenarios"]) == (steps_ahead, num_scenarious))
 
@@ -156,6 +157,8 @@
 
     #@info(" --- Test quality of fit - Normal with 2 params")
     @testset "quality of fit - Normal with 2 params" begin
+        y         = time_series[1:end-steps_ahead,5]
+        y_test    = time_series[end-steps_ahead+1:end, 5]
         gas_model = UnobservedComponentsGAS.GASModel(UnobservedComponentsGAS.NormalDistribution(), [true, true],
                                                         0.0, ["random walk slope", "random walk"], ["deterministic 12", "deterministic 12"], [missing, missing])
         fitted_model = UnobservedComponentsGAS.fit(gas_model, y)
@@ -167,6 +170,8 @@
 
     # #@info(" --- Test quality of fit - Normal with robust")
     @testset "quality of fit - Normal with robust" begin
+        y         = time_series[1:end-steps_ahead,5]
+        y_test    = time_series[end-steps_ahead+1:end, 5]
         gas_model = UnobservedComponentsGAS.GASModel(UnobservedComponentsGAS.NormalDistribution(), [true, false],
                                                         1.0, "random walk slope", "deterministic 12", 1)
         fitted_model = UnobservedComponentsGAS.fit(gas_model, y; α = 0.0, robust = true)
@@ -177,6 +182,8 @@
     end
 
     @testset "AR(1) level" begin
+        y         = time_series[1:end-steps_ahead,5]
+        y_test    = time_series[end-steps_ahead+1:end, 5]
         gas_model = UnobservedComponentsGAS.GASModel(UnobservedComponentsGAS.NormalDistribution(), [true, false],
                                                         0.5, ["ar(1)", ""], ["deterministic 12", ""], [missing, missing])
         fitted_model = UnobservedComponentsGAS.fit(gas_model, y)

diff --git a/test/test_fit_forecast_t.jl b/test/test_fit_forecast_t.jl
@@ -50,15 +50,23 @@
 
     gas_model_t_X            = deepcopy(gas_model_t)
     gas_model_t_X_2params    = deepcopy(gas_model_t_2params)
+
+    model_t, parameters_t, initial_values_t                         = UnobservedComponentsGAS.create_model(gas_model_t, y, ν)
+    model_t_2params, parameters_t_2params, initial_values_t_2params = UnobservedComponentsGAS.create_model(gas_model_t_2params, y,  ν)
+
+    model_t_X, parameters_t_X, initial_values_t_X                         = UnobservedComponentsGAS.create_model(gas_model_t_X, y, X,  ν);
+    model_t_X_2params, parameters_t_X_2params, initial_values_t_X_2params = UnobservedComponentsGAS.create_model(gas_model_t_X_2params, y, X,  ν);
+
+    fitted_model_t         = UnobservedComponentsGAS.fit(gas_model_t, y; tol = 5e-2)
+    fitted_model_t_2params = UnobservedComponentsGAS.fit(gas_model_t_2params, y; tol = 5e-2)
+    fitted_model_t_X         = UnobservedComponentsGAS.fit(gas_model_t_X, y, X; tol = 5e-2)
+    fitted_model_t_X_2params = UnobservedComponentsGAS.fit(gas_model_t_X_2params, y, X; tol = 5e-2)
+
+    forecast_t         = UnobservedComponentsGAS.predict(gas_model_t, fitted_model_t, y, steps_ahead, num_scenarious)
+    forecast_t_X       = UnobservedComponentsGAS.predict(gas_model_t_X, fitted_model_t_X, y, X_t_forec, steps_ahead, num_scenarious)
+    forecast_t_2params = UnobservedComponentsGAS.predict(gas_model_t_2params, fitted_model_t_2params, y, steps_ahead, num_scenarious)
 
-    @testset " --- Testing create_model functions" begin
-        # Create model with no explanatory series
-        model_t, parameters_t, initial_values_t                         = UnobservedComponentsGAS.create_model(gas_model_t, y, ν)
-        model_t_2params, parameters_t_2params, initial_values_t_2params = UnobservedComponentsGAS.create_model(gas_model_t_2params, y,  ν)
-
-        model_t_X, parameters_t_X, initial_values_t_X                         = UnobservedComponentsGAS.create_model(gas_model_t_X, y, X,  ν);
-        model_t_X_2params, parameters_t_X_2params, initial_values_t_X_2params = UnobservedComponentsGAS.create_model(gas_model_t_X_2params, y, X,  ν);
-
+    @testset " --- Testing create_model functions" begin        
         @test(size(parameters_t)         == (T,3))
         @test(size(parameters_t_2params) == (T,3))
         @test(typeof(model_t)            == JuMP.Model)
@@ -74,12 +82,7 @@
         @test(test_initial_values_components(initial_values_t_X_2params, rw, rws, ar, seasonality))
     end
 
-    @testset " --- Testing fit functions" begin
-        fitted_model_t         = UnobservedComponentsGAS.fit(gas_model_t, y; tol = 5e-2)
-        fitted_model_t_2params = UnobservedComponentsGAS.fit(gas_model_t_2params, y; tol = 5e-2)
-        fitted_model_t_X         = UnobservedComponentsGAS.fit(gas_model_t_X, y, X; tol = 5e-2)
-        fitted_model_t_X_2params = UnobservedComponentsGAS.fit(gas_model_t_X_2params, y, X; tol = 5e-2)
-
+    @testset " --- Testing fit functions" begin    
         # "Test if termination_status is correct"
         possible_status = ["LOCALLY_SOLVED", "TIME_LIMIT"]
         @test(fitted_model_t.model_status in possible_status)
@@ -109,11 +112,6 @@
     end
 
     @testset " --- Test forecast function ---" begin
-        forecast_t         = UnobservedComponentsGAS.predict(gas_model_t, fitted_model_t, y, steps_ahead, num_scenarious)
-        forecast_t_X       = UnobservedComponentsGAS.predict(gas_model_t_X, fitted_model_t_X, y, X_t_forec, steps_ahead, num_scenarious)
-        forecast_t_2params = UnobservedComponentsGAS.predict(gas_model_t_2params, fitted_model_t_2params, y, steps_ahead, num_scenarious)
-
-
         @test(isapprox(forecast_t["mean"], vec(mean(forecast_t["scenarios"], dims = 2)); rtol = 1e-3)) 
         @test(size(forecast_t["scenarios"]) == (steps_ahead, num_scenarious))
 
@@ -154,7 +152,8 @@
     end
 
     @testset " --- Test quality of fit - tDist with 2 params" begin
-
+        y         = time_series[1:end-steps_ahead,1]
+        y_test    = time_series[end-steps_ahead+1:end, 1]
         gas_model = UnobservedComponentsGAS.GASModel(UnobservedComponentsGAS.tLocationScaleDistribution(), [true, true, false],
                                                         0.0, ["random walk slope", "random walk", ""], ["deterministic 12", "deterministic 12", ""], [missing, missing, missing])
         fitted_model = UnobservedComponentsGAS.fit(gas_model, y)
@@ -166,6 +165,8 @@
 
 
     @testset " --- Test quality of fit - tDist with robust" begin
+        y         = time_series[1:end-steps_ahead,1]
+        y_test    = time_series[end-steps_ahead+1:end, 1]
         gas_model = UnobservedComponentsGAS.GASModel(UnobservedComponentsGAS.tLocationScaleDistribution(), [true, false, false],
                                                         1.0, "random walk slope", "deterministic 12", 1)
         fitted_model = UnobservedComponentsGAS.fit(gas_model, y; α = 0.0, robust = true)

diff --git a/test/test_initialization.jl b/test/test_initialization.jl
@@ -8,10 +8,6 @@
     stochastic = false
 
     @testset "structures" begin
-        params      = 
-        level       = 
-        seasonality = 
-        ar          = 
         gas_model = UnobservedComponentsGAS.GASModel(UnobservedComponentsGAS.NormalDistribution(), [true, false],
                                                      0.5, ["random walk", "ar(1)"], 
                                                     ["deterministic 12", "deterministic 12"],
@@ -21,7 +17,7 @@
         @test(ismissing(gas_model.ar[2]))
         @test(!isempty(gas_model.seasonality[1]))
         @test(!isempty(gas_model.level[1]))
-        @test(ismissing(gas_model.ar[1]))
+        @test(!ismissing(gas_model.ar[1]))
     end
 
 
@@ -356,7 +352,6 @@
         @test(all(initial_values["seasonality"]["γ"] .== initial_values_state_space["γ"]))
         @test(all(initial_values["seasonality"]["γ_star"] .== initial_values_state_space["γ_star"]))
         @test(all(initial_values["seasonality"]["values"] .== zeros(T)))
-        @test(all(initial_values["ar1"]["values"] .== zeros(T)))
         @test(all(initial_values["slope"]["values"] .== zeros(T)))
         @test(all(initial_values["rw"]["values"] .==  zeros(T)))
         @test(all(initial_values["explanatories"] .== initial_values_state_space["explanatory"]))

diff --git a/test/test_residuals_diagnostics.jl b/test/test_residuals_diagnostics.jl
@@ -92,4 +92,18 @@
               haskey(diagnosis_pvalues_std, "ARCH")
     end
 
+    @testset "In-Sample Statistics Tests" begin
+        stats = UnobservedComponentsGAS.get_in_sample_statistics(fitted_model, y)
+
+        @test stats["MSE"] > 0
+        @test stats["RMSE"] > 0
+        @test stats["RMSE"] == sqrt(stats["MSE"])
+        @test stats["MAE"] > 0
+        @test stats["MAPE"] > 0
+        @test stats["MAPE"] < 100
+        @test stats["AIC"] > 0
+        @test stats["BIC"] > 0
+        @test stats["AICc"] > 0
+    end
+
 end