Merge pull request #303 from LAMPSPUC/small_improvements

guilhermebodin · web-flow · commit 6a825f78f7a8 · 2022-01-17T12:58:15.000-03:00
add some improvements
diff --git 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.5.22"
+version = "0.6.0"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
@@ -19,8 +19,8 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
-Distributions = "0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24"
-MatrixEquations = "1"
+Distributions = "0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25"
+MatrixEquations = "1, 2"
 Optim = "0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 1.2"
 OrderedCollections = "1"
 Polynomials = "1, 2"
diff --git a/README.md b/README.md
@@ -26,7 +26,7 @@ model = LocalLevel(y)
 
 fit!(model)
 
-results(model)
+print_results(model)
 
 forecast(model, 10)
 
@@ -105,6 +105,7 @@ Quick examples on automatic forecasting. When performing automatic forecasting
 users should provide the seasonal period if there is one.
 ```julia
 model = auto_ets(log_air_passengers; seasonal = 12)
+model = auto_arima(log_air_passengers; seasonal = 12)
 ```
 
 ## Contributing
diff --git a/docs/src/examples.md b/docs/src/examples.md
@@ -147,7 +147,7 @@ air_passengers = CSV.File(StateSpaceModels.AIR_PASSENGERS) |> DataFrame
 log_air_passengers = log.(air_passengers.passengers)
 model = SARIMA(log_air_passengers; order = (0, 1, 1), seasonal_order = (0, 1, 1, 12))
 fit!(model)
-results(model)
+print_results(model)
 ```
 
 To make a forecast of 24 steps ahead of the model the command is:
diff --git a/docs/src/manual.md b/docs/src/manual.md
@@ -19,7 +19,7 @@ model = LocalLevel(y)
 
 fit!(model)
 
-results(model)
+print_results(model)
 
 forec = forecast(model, 10)
 
diff --git a/src/StateSpaceModels.jl b/src/StateSpaceModels.jl
@@ -119,7 +119,7 @@ export kalman_smoother
 export loglike
 export num_states
 export number_hyperparameters
-export results
+export print_results
 export set_initial_hyperparameters!
 export simulate
 export simulate_scenarios
diff --git a/src/cross_validation.jl b/src/cross_validation.jl
@@ -58,7 +58,7 @@ function cross_validation(model::StateSpaceModel, steps_ahead::Int, start_idx::I
         y_to_fit = model.system.y[1:start_idx - 1 + i]
         y_to_verify = model.system.y[start_idx + i:start_idx - 1 + i + steps_ahead]
         model_to_fit = reinstantiate(model, y_to_fit)
-        fit!(model_to_fit; filter=filter, optimizer=optimizer)
+        fit!(model_to_fit; filter=filter, optimizer=optimizer, save_hyperparameter_distribution=false)
         forec = forecast(model_to_fit, steps_ahead; filter=filter)
         scenarios = simulate_scenarios(model_to_fit, steps_ahead, n_scenarios; filter=filter)
         expected_value_vector = forecast_expected_value(forec)[:]
diff --git a/src/filters/multivariate_kalman_filter.jl b/src/filters/multivariate_kalman_filter.jl
@@ -154,7 +154,7 @@ function update_Ptt!(kalman_state::MultivariateKalmanState{Fl}, Z::Matrix{Fl}) w
 end
 
 function repeat_P_in_Ptt!(kalman_state::MultivariateKalmanState{Fl}) where Fl
-    for i in axes(kalman_state.P, 1), j in axes(kalman_state.P, 2)
+    for j in axes(kalman_state.P, 1), i in axes(kalman_state.P, 2)
         kalman_state.Ptt[i, j] = kalman_state.P[i, j]
     end
     return kalman_state
diff --git a/src/filters/sparse_univariate_kalman_filter.jl b/src/filters/sparse_univariate_kalman_filter.jl
@@ -163,7 +163,7 @@ function update_Ptt!(kalman_state::SparseUnivariateKalmanState{Fl}) where Fl
 end
 
 function repeat_P_in_Ptt!(kalman_state::SparseUnivariateKalmanState{Fl}) where Fl
-    for i in axes(kalman_state.P, 1), j in axes(kalman_state.P, 2)
+    for j in axes(kalman_state.P, 1), i in axes(kalman_state.P, 2)
         kalman_state.Ptt[i, j] = kalman_state.P[i, j]
     end
     return kalman_state
diff --git a/src/filters/univariate_kalman_filter.jl b/src/filters/univariate_kalman_filter.jl
@@ -131,7 +131,7 @@ end
 
 function update_att!(kalman_state::UnivariateKalmanState{Fl}, Z::Vector{Fl}) where Fl
     copyto!(kalman_state.att, kalman_state.a)
-    @inbounds for i in axes(kalman_state.P, 1), j in axes(kalman_state.P, 2)
+    @inbounds for j in axes(kalman_state.P, 1), i in axes(kalman_state.P, 2)
         kalman_state.att[i] +=
             (kalman_state.v / kalman_state.F) * kalman_state.P[i, j] * Z[j]
     end
@@ -162,7 +162,7 @@ function update_Ptt!(kalman_state::UnivariateKalmanState{Fl}) where Fl
 end
 
 function repeat_P_in_Ptt!(kalman_state::UnivariateKalmanState{Fl}) where Fl
-    for i in axes(kalman_state.P, 1), j in axes(kalman_state.P, 2)
+    for j in axes(kalman_state.P, 1), i in axes(kalman_state.P, 2)
         kalman_state.Ptt[i, j] = kalman_state.P[i, j]
     end
     return kalman_state
diff --git a/src/fit.jl b/src/fit.jl
@@ -2,7 +2,8 @@
     fit!(
         model::StateSpaceModel;
         filter::KalmanFilter=default_filter(model),
-        optimizer::Optimizer=Optimizer(Optim.LBFGS())
+        optimizer::Optimizer=Optimizer(Optim.LBFGS()),
+        save_hyperparameter_distribution::Bool=true
     )
 
 Estimate the state-space model parameters via maximum likelihood. The resulting optimal
@@ -29,6 +30,7 @@ function fit!(
     model::StateSpaceModel;
     filter::KalmanFilter=default_filter(model),
     optimizer::Optimizer=default_optimizer(model),
+    save_hyperparameter_distribution::Bool=true
 )
     isfitted(model) && return model
     @assert has_fit_methods(typeof(model))
@@ -45,12 +47,29 @@ function fit!(
     opt_loglikelihood = -opt.minimum * size(model.system.y, 1)
     opt_hyperparameters = opt.minimizer
     update_model_hyperparameters!(model, opt_hyperparameters)
-    # TODO
-    # I leaned that this is not a good way to compute the covariance matrix of paarameters
-    # we should investigate other methods
-    numerical_hessian = Optim.hessian!(func, opt_hyperparameters)
-    std_err = diag(pinv(numerical_hessian))
-    fill_results!(model, opt_loglikelihood, std_err)
+
+    if save_hyperparameter_distribution
+        numerical_hessian = Optim.hessian!(func, opt_hyperparameters)
+        try 
+            std_err = numerical_hessian |> pinv |> diag .|> sqrt
+            fill_results!(model, opt_loglikelihood, std_err)
+        catch
+            @warn(
+                "The optimization process converged but the Hessian matrix is not positive definite. " *
+                "This means that StateSpaceModels.jl cannot estimate the distribution of the hyperparameters " *
+                "If you are interested in estimates of the distribution of ther hyperparameters we advise you to" *
+                "change the optimization algorithm by using the kwarg fit(...; optimizer = "*
+                "Optimizer(StateSpaceModels.Optim.THE_METHOD_OF_YOUR_CHOICE()))" * 
+                "The list of possible algorithms can be found on this link https://julianlsolvers.github.io/Optim.jl/stable/#" * 
+                "Otherwise you can simply skip this proccess by using fit(...; save_hyperparameter_distribution=false) "
+            )
+            std_err = fill(NaN, number_hyperparameters(model))
+            fill_results!(model, opt_loglikelihood, std_err)
+        end
+    else
+        std_err = fill(NaN, number_hyperparameters(model))
+        fill_results!(model, opt_loglikelihood, std_err)
+    end
     return model
 end
 
@@ -121,12 +140,7 @@ function Results{Fl}() where Fl
     return Results{Fl}("", CoefficientTable{Fl}(), Fl(NaN), Fl(NaN), Fl(NaN), Fl(NaN), 0, 0)
 end
 
-"""
-    results(model::StateSpaceModel)
-
-Query the results of the optimization called by `fit!`.
-"""
-results(model::StateSpaceModel) = model.results
+get_results(model::StateSpaceModel) = model.results
 function Base.isempty(results::Results)
     return isempty(results.coef_table) &&
            isnan(results.llk) &&
diff --git a/src/models/sarima.jl b/src/models/sarima.jl
@@ -761,7 +761,7 @@ function fit_candidate_models!(candidate_models::Vector{SARIMA}, show_trace::Boo
     non_converged_models = Int[]
     for (i, model) in enumerate(candidate_models)
         try 
-            fit!(model)
+            fit!(model; save_hyperparameter_distribution=false)
             if isnan(model.results.llk)
                 show_trace && println(model, " - diverged")
                 push!(non_converged_models, i)
diff --git a/src/prints.jl b/src/prints.jl
@@ -1,35 +1,39 @@
-function Base.show(io::IO, results::Results{Fl}) where Fl
+"""
+    print_results(model::StateSpaceModel)
+
+Query the results of the optimization called by `fit!`.
+"""
+function print_results(model::StateSpaceModel)
+    results = get_results(model)
     if !isempty(results)
-        println(io, "                             Results                           ")
-        println(io, "===============================================================" * adjust_charachters("=", results.coef_table))
-        println(io, "Model:                        ", results.model_name)
-        println(io, "Number of observations:       ", results.num_observations)
-        println(io, "Number of unknown parameters: ", results.num_hyperparameters)
-        println(io, "Log-likelihood:               ", pretty_number(results.llk))
-        println(io, "AIC:                          ", pretty_number(results.aic))
-        println(io, "AICc:                         ", pretty_number(results.aicc))
-        println(io, "BIC:                          ", pretty_number(results.bic))
-        print_coef_table(io, results.coef_table)
+        println("                             Results                           ")
+        println("===============================================================" * adjust_charachters("=", results.coef_table))
+        println("Model:                        ", results.model_name)
+        println("Number of observations:       ", results.num_observations)
+        println("Number of unknown parameters: ", results.num_hyperparameters)
+        println("Log-likelihood:               ", pretty_number(results.llk))
+        println("AIC:                          ", pretty_number(results.aic))
+        println("AICc:                         ", pretty_number(results.aicc))
+        println("BIC:                          ", pretty_number(results.bic))
+        print_coef_table(results.coef_table)
     else
-        # Up to discussion of what should be the behaviour
-        error("You must call fit! to get results.")
+        println("Results are empty.")
     end
-    return nothing
 end
 
 function Base.show(io::IO, model::StateSpaceModel)
     print(io, model_name(model))
     return nothing
 end
 
-function print_coef_table(io::IO, coef_table::CoefficientTable{Fl}) where Fl
+function print_coef_table(coef_table::CoefficientTable{Fl}) where Fl
     # The repeats adjust the space for large names
-    println(io, "---------------------------------------------------------------" * adjust_charachters("-", coef_table))
-    println(io, "Parameter" * adjust_charachters(" ", coef_table) * "      Estimate      Std.Error      z stat      p-value")
+    println("---------------------------------------------------------------" * adjust_charachters("-", coef_table))
+    println("Parameter" * adjust_charachters(" ", coef_table) * "      Estimate      Std.Error      z stat      p-value")
     offset = 1
     for i in 1:length(coef_table)
         p_c, p_std, p_z_stat, p_p_val = print_coef(coef_table, offset)
-        println(io, build_print(p_c, p_std, p_z_stat, p_p_val, offset, coef_table))
+        println(build_print(p_c, p_std, p_z_stat, p_p_val, offset, coef_table))
         offset += 1
     end
     return nothing
diff --git a/test/models/linear_regression.jl b/test/models/linear_regression.jl
@@ -6,7 +6,7 @@
         X = rand(100, 5)
         y = rand(100)
         model = LinearRegression(X, y)
-        fit!(model)
+        fit!(model; save_hyperparameter_distribution=false)
         β = X \ y
         for i in 1:size(X, 2)
             @test get_constrained_value(model, "β_$i") .≈ β[i] atol = 1e-5
@@ -27,7 +27,7 @@
         10.0
     ]
     model = LinearRegression(X, y)
-    fit!(model)
+    fit!(model; save_hyperparameter_distribution=false)
     forec = forecast(model, [1 6.0])
     forec.expected_value[1] == [11.0]
 end
diff --git a/test/models/locallevel.jl b/test/models/locallevel.jl
@@ -14,7 +14,6 @@ using CSV, DataFrames
     @test_throws ErrorException get_filtered_state_variance(model)
     @test_throws ErrorException get_predictive_state(model)
     @test_throws ErrorException get_predictive_state_variance(model)
-    @test_throws ErrorException show(stdout, results(model))
 
     fit!(model)
     @test loglike(model) ≈ -632.5376 atol = 1e-5 rtol = 1e-5
@@ -71,7 +70,7 @@ using CSV, DataFrames
 
     filter = kalman_filter(model)
     smoother = kalman_smoother(model)
-    @test filter.Ptt[end] == smoother.V[end] # by construction
+    @test filter.Ptt[end] ≈ smoother.V[end] atol = 1e-7 # by construction
 
     for t in 2:(length(model.system.y) - 1)
         @test filter.Ptt[t][1] > smoother.V[t][1] # by construction
diff --git a/test/prints.jl b/test/prints.jl
@@ -5,7 +5,7 @@
     show(model)
     println("")
     fit!(model)
-    show(results(model))
+    print_results(model)
     println("")
     # Specific show methods
     model = SARIMA(random_series)
diff --git a/test/visualization/forecast.jl b/test/visualization/forecast.jl
@@ -2,7 +2,7 @@
     air_passengers = CSV.File(StateSpaceModels.AIR_PASSENGERS) |> DataFrame
     log_air_passengers = log.(air_passengers.passengers)
     model = BasicStructural(log_air_passengers, 12)
-    fit!(model)
+    fit!(model; save_hyperparameter_distribution=false)
     # forecasting
     forec = forecast(model, 12)
     rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), model, forec)
