Merge pull request #20 from LAMPSPUC/add_gamma

Add gamma
LAMPSPUC · Sep 3, 2024 · 693c2d3 · 693c2d3
2 parents 9301d9d + 8ca507e
commit 693c2d3
Show file tree

Hide file tree

Showing 11 changed files with 507 additions and 98 deletions.
diff --git a/.gitignore b/.gitignore
@@ -1,10 +1,16 @@
 test_gas_alves.jl
-Manifest.toml
+
 forecast_example.jl
 residuals_diagnostic_example.jl
 test_numeric_integration.jl
 test_m3.jl
 test_dynamic_comb.jl
 inovations_ssmodels.jl
 test_component_expression.jl
-results_variables_expressions.csv
+test_gas.jl
+
+# Manifest file
+Manifest.toml
+
+# Files with csv extension
+*.csv
diff --git a/src/UnobservedComponentsGAS.jl b/src/UnobservedComponentsGAS.jl
@@ -23,6 +23,7 @@ module UnobservedComponentsGAS
     include("distributions/normal.jl")
     include("distributions/t_location_scale.jl")
     include("distributions/log_normal.jl")
+    include("distributions/gamma.jl")
     include("initialization.jl")
     include("fit.jl")
     include("utils.jl")
@@ -32,17 +33,22 @@ module UnobservedComponentsGAS
     include("residuals_diagnostics.jl")
 
     const DICT_CODE = Dict(1 => "Normal",
-                           2 => "tLocationScale" )
+                           2 => "tLocationScale",
+                           3 => "Gamma")
 
     const DICT_SCORE = Dict("Normal"         => score_normal,
-                            "tLocationScale" => score_tlocationscale)
+                            "tLocationScale" => score_tlocationscale,
+                            "Gamma"          => score_gamma)
 
     const DICT_FISHER_INFORMATION = Dict("Normal"         => fisher_information_normal,
-                                         "tLocationScale" => fisher_information_tlocationscale)
+                                         "tLocationScale" => fisher_information_tlocationscale,
+                                         "Gamma"          => fisher_information_gamma)
 
     const DICT_LOGPDF = Dict("Normal"         => logpdf_normal,
-                             "tLocationScale" => logpdf_tlocationscale)
+                             "tLocationScale" => logpdf_tlocationscale,
+                             "Gamma"          => logpdf_gamma)
 
     const DICT_CDF = Dict("Normal"         => cdf_normal,
-                          "tLocationScale" => cdf_tlocationscale)
+                          "tLocationScale" => cdf_tlocationscale,
+                          "Gamma"          => cdf_gamma)
 end 
diff --git a/src/components_dynamics.jl b/src/components_dynamics.jl
@@ -221,7 +221,8 @@ Adds AutoRegressive (AR) components to the optimization model.
 ### Returns
 - Modifies the provided optimization model by adding AR components.
 """
-function add_AR!(model::Ml, s::Vector{Fl}, T::Int64, ar::Union{Int64, Vector{Int64}, Vector{Missing}, Vector{Union{Int64, Missing}}, Missing, Vector{Vector{Int64}}, Vector{Union{Vector{Int64}, Missing}}}) where {Ml, Fl}
+function add_AR!(model::Ml, s::Vector{Fl}, T::Int64, ar::Union{Int64, Vector{Int64}, Vector{Missing}, Vector{Union{Int64, Missing}}, Missing, Vector{Vector{Int64}}, Vector{Union{Vector{Int64}, Missing}}};
+                κ_min::Int64 = 0, κ_max::Int64 = 2) where {Ml, Fl}
 
     idx_params = findall(i -> !ismissing(i), ar) # Time-varying parameters with autoregressive dynamic
     order      = get_AR_order(ar)
@@ -233,7 +234,7 @@ function add_AR!(model::Ml, s::Vector{Fl}, T::Int64, ar::Union{Int64, Vector{Int
     @variable(model, ϕ[1:max_order, idx_params])
     @variable(model, κ_AR[idx_params])
 
-    @constraint(model, [i in idx_params], -2 ≤ κ_AR[i] ≤ 2)
+    @constraint(model, [i in idx_params], κ_min ≤ κ_AR[i] ≤ κ_max)
 
     # Revisar essas restrições com o Alves !!
     for i in unique_orders
@@ -262,7 +263,8 @@ Incorporate the random walk with slope component into the dynamics of the specif
 ## Returns
 - Modifies the input model by adding random walk with slope components.
 """
-function add_random_walk_slope!(model::Ml, s::Vector{Fl}, T::Int64, random_walk_slope::Dict{Int64, Bool}) where {Ml, Fl}
+function add_random_walk_slope!(model::Ml, s::Vector{Fl}, T::Int64, random_walk_slope::Dict{Int64, Bool};
+                                κ_min::Int64 = 0, κ_max::Int64 = 2) where {Ml, Fl}
 
     idx_params = findall(i -> i == true, random_walk_slope) #Time-varying parameters with the random walk with slope dynamic
 
@@ -273,8 +275,8 @@ function add_random_walk_slope!(model::Ml, s::Vector{Fl}, T::Int64, random_walk_
 
     @constraint(model, [t = 2:T, j in idx_params], b[t, j]   == b[t - 1, j] + κ_b[j] * s[j][t])
     @constraint(model, [t = 2:T, j in idx_params], RWS[t, j] == RWS[t - 1, j] + b[t - 1, j] + κ_RWS[j] * s[j][t])
-    @constraint(model, [j in idx_params], -2 ≤ κ_RWS[j] ≤ 2)
-    @constraint(model, [j in idx_params], -2 ≤ κ_b[j] ≤ 2)
+    @constraint(model, [j in idx_params], κ_min ≤ κ_RWS[j] ≤ κ_max)
+    @constraint(model, [j in idx_params], κ_min ≤ κ_b[j] ≤ κ_max)
 end
 
 """
@@ -291,15 +293,16 @@ Incorporate the random walk component into the dynamics of the specified paramet
 ## Returns
 - Modifies the input model by adding random walk components.
 """
-function add_random_walk!(model::Ml, s::Vector{Fl}, T::Int64, random_walk::Dict{Int64, Bool}) where {Ml, Fl}
+function add_random_walk!(model::Ml, s::Vector{Fl}, T::Int64, random_walk::Dict{Int64, Bool};
+                        κ_min::Int64 = 0, κ_max::Int64 = 2) where {Ml, Fl}
 
     idx_params = findall(i -> i == true, random_walk) # Time-varying parameters with the random walk dynamic
 
     @variable(model, RW[1:T, idx_params])
     @variable(model, κ_RW[idx_params])
 
     @constraint(model, [t = 2:T, j in idx_params], RW[t, j] == RW[t-1, j] + κ_RW[j] * s[j][t])
-    @constraint(model, [j in idx_params], -2 ≤ κ_RW[j] ≤ 2)
+    @constraint(model, [j in idx_params], κ_min ≤ κ_RW[j] ≤ κ_max)
 end
 
 """
@@ -316,7 +319,8 @@ Adds AutoRegressive (AR(1)) components to the optimization model.
 ### Returns
 - Modifies the provided optimization model by adding AR(1) components.
 """
-function add_ar1!(model::Ml, s::Vector{Fl}, T::Int64, ar1::Dict{Int64, Bool})  where {Fl, Ml}
+function add_ar1!(model::Ml, s::Vector{Fl}, T::Int64, ar1::Dict{Int64, Bool};
+                κ_min::Int64 = 0, κ_max::Int64 = 2)  where {Fl, Ml}
 
     idx_params = findall(i -> i == true, ar1)
 
@@ -325,7 +329,7 @@ function add_ar1!(model::Ml, s::Vector{Fl}, T::Int64, ar1::Dict{Int64, Bool})  w
     @variable(model, κ_AR1_LEVEL[idx_params])
     @variable(model, ω_AR1_LEVEL[idx_params])
 
-    @constraint(model, [i in idx_params], -2 ≤ κ_AR1_LEVEL[i] ≤ 2)
+    @constraint(model, [i in idx_params], κ_min ≤ κ_AR1_LEVEL[i] ≤ κ_max)
     @constraint(model, [i in idx_params], -0.95 <= ϕ_AR1_LEVEL[i] <= 0.95)
 
     @constraint(model, [t = 2:T, j in idx_params], AR1_LEVEL[t, j] == ω_AR1_LEVEL[j] + ϕ_AR1_LEVEL[j] * AR1_LEVEL[t-1, j] + κ_AR1_LEVEL[j] * s[j][t])
@@ -346,29 +350,31 @@ Adds level components to the optimization model based on the provided level info
 ### Returns
 - Modifies the provided optimization model by adding level components.
 """
-function add_level!(model::Ml, s::Vector{Fl}, T::Int64, level::Vector{String}) where {Fl, Ml}
+function add_level!(model::Ml, s::Vector{Fl}, T::Int64, level::Vector{String};
+                    κ_min::Int64 = 0, κ_max::Int64 = 2) where {Fl, Ml}
+
     if "random walk" ∈ level 
         random_walk = Dict{Int64, Bool}()
         for i in 1:length(level)
             level[i] == "random walk" ? random_walk[i] = true : random_walk[i] = false
         end
-        add_random_walk!(model, s, T, random_walk)
+        add_random_walk!(model, s, T, random_walk; κ_min = κ_min, κ_max = κ_max)
     end
 
     if "random walk slope" ∈ level 
         random_walk_slope = Dict{Int64, Bool}()
         for i in 1:length(level)
             level[i] == "random walk slope" ? random_walk_slope[i] = true : random_walk_slope[i] = false
         end
-        add_random_walk_slope!(model, s, T, random_walk_slope)
+        add_random_walk_slope!(model, s, T, random_walk_slope; κ_min = κ_min, κ_max = κ_max)
     end
 
     if "ar(1)" ∈ level 
         ar1 = Dict{Int64, Bool}()
         for i in 1:length(level)
             level[i] == "ar(1)" ? ar1[i] = true : ar1[i] = false
         end
-        add_ar1!(model, s, T, ar1)
+        add_ar1!(model, s, T, ar1; κ_min = κ_min, κ_max = κ_max)
     end
 end
 
@@ -418,7 +424,8 @@ Adds trigonometric seasonality components to the optimization model based on the
 ### Returns
 - Modifies the provided optimization model by adding trigonometric seasonality components.
 """
-function add_trigonometric_seasonality!(model::Ml, s::Vector{Fl}, T::Int64, seasonality::Vector{String}) where {Ml, Fl}
+function add_trigonometric_seasonality!(model::Ml, s::Vector{Fl}, T::Int64, seasonality::Vector{String};
+                                        κ_min::Int64 = 0, κ_max::Int64 = 2) where {Ml, Fl}
 
     seasonality_dict, stochastic, stochastic_params = get_seasonality_dict_and_stochastic(seasonality)
 
@@ -435,7 +442,7 @@ function add_trigonometric_seasonality!(model::Ml, s::Vector{Fl}, T::Int64, seas
 
 
         @variable(model, κ_S[idx_params_stochastic])
-        @constraint(model, [i in idx_params_stochastic], -2 ≤ κ_S[i] ≤ 2)    
+        @constraint(model, [i in idx_params_stochastic], κ_min  ≤ κ_S[i] ≤ κ_max)    
         #JuMP.fix.(model[:κ_S][idx_params_deterministic], 1e-4)
 
         @variable(model, γ_sto[1:unique_num_harmonic, 1:T, idx_params_stochastic])
@@ -471,36 +478,6 @@ function add_trigonometric_seasonality!(model::Ml, s::Vector{Fl}, T::Int64, seas
 
     @expression(model, S[t=1:T, j in idx_params], S_aux[t, j])
 
-    # for i in idx_params
-    #     if idx_params[i] ∈ stochastic_params
-
-
-    #     elseif idx_params[i] ∈ idx_params_deterministic
-    #     end
-    # end
-
-    # if stochastic
-    #     @variable(model, κ_S[idx_params])
-    #     @constraint(model, [i in idx_params], 1e-4 ≤ κ_S[i])    
-    #     JuMP.fix.(model[:κ_S][idx_params_deterministic], 1e-4)
-
-    #     @variable(model, γ[1:unique_num_harmonic, 1:T, idx_params])
-    #     @variable(model, γ_star[1:unique_num_harmonic, 1:T, idx_params])
-
-    #     @constraint(model, [i = 1:unique_num_harmonic, t = 2:T, j in idx_params], γ[i, t, j] == γ[i, t-1, j] * cos(2*π*i / seasonal_period[j]) + 
-    #                                                                                 γ_star[i,t-1, j]*sin(2*π*i / seasonal_period[j]) + κ_S[j] * s[j][t])
-    #     @constraint(model, [i = 1:unique_num_harmonic, t = 2:T, j in idx_params], γ_star[i, t, j] == -γ[i, t-1, j] * sin(2*π*i / seasonal_period[j]) + 
-    #                                                                                 γ_star[i,t-1, j]*cos(2*π*i / seasonal_period[j]) + κ_S[j] * s[j][t])
-
-    #     @expression(model, S[t = 1:T, j in idx_params], sum(γ[i, t, j]  for i in 1:unique_num_harmonic))
-    # else
-    #     @variable(model, γ[1:unique_num_harmonic, idx_params])
-    #     @variable(model, γ_star[1:unique_num_harmonic, idx_params])
-
-    #     @expression(model, S[t = 1:T, j in idx_params], sum(γ[i, j]*cos(2 * π * i * t/seasonal_period[j]) + 
-    #                                         γ_star[i, j] * sin(2 * π * i* t/seasonal_period[j]) for i in 1:unique_num_harmonic))
-    # end
-
 end
 
 """
@@ -517,20 +494,21 @@ Incorporates various components into the specified model based on the configurat
 ## Returns
 - Modifies the input model by adding components such as random walk, random walk slope, autoregressive (AR), and trigonometric seasonality based on the configurations provided in the `GASModel`.
 """
-function include_components!(model::Ml, s::Vector{Fl}, gas_model::GASModel, T::Int64) where {Ml, Fl}
+function include_components!(model::Ml, s::Vector{Fl}, gas_model::GASModel, T::Int64;
+                            κ_min::Int64 = 0, κ_max::Int64 = 2) where {Ml, Fl}
 
     @unpack dist, time_varying_params, d, level, seasonality, ar = gas_model
 
     if has_level(level)
-        add_level!(model, s, T, level)
+        add_level!(model, s, T, level; κ_min = κ_min, κ_max = κ_max)
     end
 
     if has_AR(ar)
-        add_AR!(model, s, T, ar)
+        add_AR!(model, s, T, ar; κ_min = κ_min, κ_max = κ_max)
     end
 
     if has_seasonality(seasonality)
-        add_trigonometric_seasonality!(model, s, T, seasonality)
+        add_trigonometric_seasonality!(model, s, T, seasonality; κ_min = κ_min, κ_max = κ_max)
     end
 end