Include explanatory tests init

src/components_dynamics.jl

@@ -51,7 +51,7 @@ Checks if the provided autoregressive (AR) dictionary indicates the presence of
 ## Returns
 - `true` if there is at least one component with autoregressive structure, `false` otherwise.
 """
-has_AR(ar::Union{Dict{Int64, Int64}, Dict{Int64, Vector{Int64}}, Dict{Int64, Bool}, Dict{Int64, Any}}) = !all(isequal.(typeof.(values(ar)), Bool))
+has_AR(ar::Union{Dict{Int64, Int64}, Dict{Int64, Vector{Int64}}, Dict{Int64, Bool}, Dict{Int64, Any}, Dict{Int64, Integer}}) = !all(isequal.(typeof.(values(ar)), Bool))
 
 
 """
@@ -111,7 +111,7 @@ Checks if the specified parameter has an autoregressive (AR) component in the gi
 ## Returns
 - `true` if the specified parameter has an autoregressive (AR) component, `false` otherwise.
 """
-has_AR(ar::Union{Dict{Int64, Int64}, Dict{Int64, Vector{Int64}}, Dict{Int64, Bool}, Dict{Int64, Any}}, param::Int64) = typeof(ar[param]) == Bool || ar[param] == 0 ? false : true
+has_AR(ar::Union{Dict{Int64, Int64}, Dict{Int64, Vector{Int64}}, Dict{Int64, Bool}, Dict{Int64, Any}, Dict{Int64, Integer}}, param::Int64) = typeof(ar[param]) == Bool || ar[param] == 0 ? false : true
 
 """
 # get_AR_order(ar::Union{Dict{Int64, Int64}, Dict{Int64, Vector{Int64}}, Dict{Int64, Bool}, Dict{Int64, Any}})
@@ -129,23 +129,25 @@ Extracts the autoregressive (AR) orders for each parameter from the given dictio
 """
 
 # To do: Troquei o retorno da função caso não tenha AR para nothing. Lembrar de trocar isso nos if's de outras funções.
-function get_AR_order(ar::Union{Dict{Int64, Int64}, Dict{Int64, Vector{Int64}}, Dict{Int64, Bool}, Dict{Int64, Any}})
+function get_AR_order(ar::Union{Dict{Int64, Int64}, Dict{Int64, Vector{Int64}}, Dict{Int64, Bool}, Dict{Int64, Any}, Dict{Int64, Integer}})
 
     num_params = length(ar)
     order      = Union{Vector{Int64}, Vector{Nothing}}[]
 
     for i in 1:num_params
-        if typeof(ar[i]) == Int64
+
+        if ar[i] == 0 || isnothing(ar[i])
+            push!(order, [nothing])
+        elseif typeof(ar[i]) == Int64
             push!(order, Int64.(collect(1:ar[i])))
         elseif typeof(ar[i]) == Vector{Int64}
             push!(order, Int64.(ar[i]))
-        else
-            push!(order, [nothing])
         end
     end
     return order
 end
 
+
 """
 # add_AR!(model::Ml, s::Vector{Fl}, T::Int64, ar::Union{Dict{Int64, Int64}, Dict{Int64, Vector{Int64}}, Dict{Int64, Bool}, Dict{Int64, Any}}) where {Ml, Fl}
 
@@ -160,12 +162,12 @@ Incorporate the autoregressive (AR) component into the dynamics of the specified
 ## Returns
 - Modifies the input model by adding autoregressive (AR) dynamics.
 """
-function add_AR!(model::Ml, s::Vector{Fl}, T::Int64, ar::Union{Dict{Int64, Int64}, Dict{Int64, Vector{Int64}}, Dict{Int64, Bool}, Dict{Int64, Any}}) where {Ml, Fl}
+function add_AR!(model::Ml, s::Vector{Fl}, T::Int64, ar::Union{Dict{Int64, Int64}, Dict{Int64, Vector{Int64}}, Dict{Int64, Bool}, Dict{Int64, Any}, Dict{Int64, Integer}}) where {Ml, Fl}
 
     idx_params = findall(i -> i != 0, ar) # Time-varying parameters with autoregressive dynamic
     order      = get_AR_order(ar)
 
-    max_order     = maximum(vcat(order...)) # Maximum lag in the model
+    max_order     = maximum(filter(x -> !isnothing(x), vcat(order...))) # Maximum lag in the model
     unique_orders = filter(x -> !isnothing(x), unique(vcat(order...)))#[findall(i -> i != 0.0, vcat(order...))]
 
     @variable(model, AR[1:T, idx_params])

src/forecast.jl

@@ -61,7 +61,7 @@ function get_dict_hyperparams_and_fitted_components_with_forecast(gas_model::GAS
     end
 
     dict_hyperparams_and_fitted_components["ar"]["value"]  = zeros(num_params, T_fitted + steps_ahead, num_scenarios)
-    dict_hyperparams_and_fitted_components["ar"]["ϕ"]      = zeros(maximum(vcat(order...)), num_params)
+    dict_hyperparams_and_fitted_components["ar"]["ϕ"]      = zeros(maximum(filter(x -> !isnothing(x), vcat(order...))), num_params)
     dict_hyperparams_and_fitted_components["ar"]["κ"]      = zeros(num_params)
 
     for i in 1:num_params

src/initialization.jl

@@ -1,53 +1,142 @@
-function get_initial_values(y::Vector{Float64}, X::Union{Matrix{Float64}, Missing}, has_level::Bool, has_slope::Bool, has_seasonality::Bool, seasonal_period::Union{Missing, Int64}, stochastic::Bool, order::Union{Vector{Int64}, Vector{Nothing}}, max_order::Int64)
+# Fit and return the predictive state of a StateSpaceModel
+function fit_and_get_preditive_state(model::M) where M
+    StateSpaceModels.fit!(model)
+    return StateSpaceModels.get_predictive_state(model)
+end
 
-    #T = length(y)
-    has_explanatories = !ismissing(X) ? true : false
+# Case without explanatory
+function define_state_space_model(y::Vector{Float64}, has_level::Bool, has_slope::Bool, 
+                                has_seasonality::Bool, seasonal_period::Union{Missing, Int64}, stochastic::Bool)
 
-    if has_level || has_slope || has_seasonality
+    T                   = length(y)
+    initial_level       = zeros(T)
+    initial_slope       = zeros(T)
+    initial_seasonality = zeros(T)
+    initial_γ           = zeros(1)
+    initial_γ_star      = zeros(1)
+    if has_seasonality
+        if !has_level && !has_slope
+            ss_model = SeasonalNaive(y, seasonal_period)
+            StateSpaceModels.fit!(ss_model)
+            initial_seasonality = ss_model.y .- vcat(zeros(seasonal_period), ss_model.residuals)
+            res = ss_model.residuals
+            initial_γ, initial_γ_star = fit_harmonics(initial_seasonality, seasonal_period, stochastic)
+        else
+            #Basic structural
+            ss_model                  = BasicStructural(y, seasonal_period)
+            pred_state                = fit_and_get_preditive_state(ss_model)
+            for t in 1:T
+                initial_seasonality[t] = -sum(pred_state[t+1, end-(seasonal_period-2):end])
+            end
+            initial_γ, initial_γ_star = fit_harmonics(initial_seasonality, seasonal_period, stochastic)
 
-        if has_level && !has_slope
-            trend = "local level"
-        elseif has_level && has_slope
-            trend = "local linear trend"
+            if has_level && has_slope
+                initial_level = pred_state[2:end,1]
+                initial_slope = pred_state[2:end,2]
+            elseif has_level && !has_slope
+                initial_level = pred_state[2:end,1] + pred_state[2:end,2]
+            end
+            res = StateSpaceModels.get_innovations(ss_model)[:, 1]
+        end
+    else
+        if has_level && has_slope
+            #Local Linear Trend
+            ss_model   = LocalLinearTrend(y)
+            pred_state = fit_and_get_preditive_state(ss_model)
+            initial_level = pred_state[2:end,1]
+            initial_slope = pred_state[2:end,2]
+        elseif has_level && !has_slope
+            # Local Level
+            ss_model   = LocalLevel(y)
+            pred_state = fit_and_get_preditive_state(ss_model)
+            initial_level = pred_state[2:end,1]
         end
+        res = StateSpaceModels.get_innovations(ss_model)[:, 1]
+    end
+
+    return Dict("level" => initial_level,"slope" => initial_slope,"seasonality" => initial_seasonality,
+            "γ" => initial_γ,"γ_star" => initial_γ_star,"explanatory" => missing,"res" => res)
+end
 
-        if has_seasonality && stochastic
-            seasonal = "stochastic "*string(seasonal_period)
-        elseif has_seasonality && !stochastic
-            seasonal = "deterministic "*string(seasonal_period)
+# Case with explanatory
+function define_state_space_model(y::Vector{Float64}, X::Union{Matrix{Float64}, Missing}, has_level::Bool, has_slope::Bool, 
+                                has_seasonality::Bool, seasonal_period::Union{Missing, Int64}, stochastic::Bool)
+
+    T                   = length(y)
+    N                   = size(X, 2)
+    initial_level       = zeros(T)
+    initial_slope       = zeros(T)
+    initial_seasonality = zeros(T)
+    initial_γ           = zeros(1)
+    initial_γ_star      = zeros(1)
+    explanatory_coefs   = zeros(N)
+
+    # Be aware that if we select just seasonal as the mean component with explanatories, 
+    ## there wont be a initial value for the explanatories coefs (it will be zero)
+    if has_seasonality
+        if !has_level && !has_slope
+            ss_model = SeasonalNaive(y, seasonal_period)
+            StateSpaceModels.fit!(ss_model)
+            initial_seasonality = ss_model.y .- vcat(zeros(seasonal_period), ss_model.residuals)
+            explanatory_coefs   = zeros(N)
+            res = ss_model.residuals
+            initial_γ, initial_γ_star = fit_harmonics(initial_seasonality, seasonal_period, stochastic)
+        else
+         #Basic structural
+            ss_model   = BasicStructuralExplanatory(y, seasonal_period, X)
+            pred_state = fit_and_get_preditive_state(ss_model)
+            for t in 1:T
+                initial_seasonality[t] = -sum(pred_state[t+1, end-(seasonal_period+N-2):end-N])
+            end
+            initial_γ, initial_γ_star = fit_harmonics(initial_seasonality, seasonal_period, stochastic)
+
+            if has_level && has_slope
+                initial_level = pred_state[2:end,1]
+                initial_slope = pred_state[2:end,2]
+            elseif has_level && !has_slope
+                initial_level = pred_state[2:end,1] + pred_state[2:end,2]
+            end
+            res = StateSpaceModels.get_innovations(ss_model)[:, 1]
         end
-
-        state_space_model = StateSpaceModels.UnobservedComponents(Float64.(y) ; trend = trend ,seasonal = seasonal) 
-        # if has_explanatories
-        StateSpaceModels.fit!(state_space_model)
+    else
+        if has_level && has_slope 
+            # Since there is no LocalLinearTrendExplanatory ...
+            ss_model   = BasicStructuralExplanatory(y, 12, X)
+            pred_state = fit_and_get_preditive_state(ss_model)
+            initial_level = pred_state[2:end,1]
+            initial_slope = pred_state[2:end,2]
+        elseif has_level && !has_slope
+            # Local Level
+            ss_model   = LocalLevelExplanatory(y, X)
+            pred_state = fit_and_get_preditive_state(ss_model)
+            initial_level = pred_state[2:end,1]
+        end
+        res = StateSpaceModels.get_innovations(ss_model)[:, 1]
+    end
+
+    explanatory_coefs = ss_model.hyperparameters.constrained_values[end-N+1:end]
 
-        #pred_state = StateSpaceModels.get_predictive_state(state_space_model)
-
-        # inov = StateSpaceModels.get_innovations(state_space_model)
-        # fit_1 = y .- vec(inov)
-        # fit_2 = sum(a[:, i] for i in 1:13)[2:end] .+ inov[:, 1]
+    return Dict("level" => initial_level,"slope" => initial_slope,"seasonality" => initial_seasonality,
+                "γ" => initial_γ,"γ_star" => initial_γ_star,"explanatory" => explanatory_coefs,"res" => res)
+end
 
-        # plot(fit_1)
-        # plot!(fit_2)
 
+function get_initial_values(y::Vector{Float64}, X::Union{Matrix{Float64}, Missing}, has_level::Bool, has_slope::Bool, has_seasonality::Bool, seasonal_period::Union{Missing, Int64}, stochastic::Bool, order::Union{Vector{Int64}, Vector{Nothing}}, max_order::Int64)
 
-        # kf = kalman_filter(state_space_model)
+    #T = length(y)
+    has_explanatories = !ismissing(X) ? true : false
 
-        # p1 = plot(y[15:end], label = "serie") #serie
-        # p2 = plot(a[15:end, 1], label="level") #level
-        # p3 = plot(a[15:end, 2], label = "slope") #slope
-        # p4 = plot(s[15:end], label = "seasonal")
-        # plot(p1,p2,p3,p4, layout=(4,1))
+    if has_level || has_slope || has_seasonality
 
-        # plot(state_space_model, kf)
-        # s = zeros(length(y))
-        # for t in  1:144
-        #     s[t] = -sum(a[t, 3:end])
-        # end
+        if has_explanatories
+            ss_components = define_state_space_model(y, X, has_level, has_slope, has_seasonality, seasonal_period, stochastic)
+        else
+            ss_components = define_state_space_model(y, has_level, has_slope, has_seasonality, seasonal_period, stochastic)
+        end
 
+        println("------",order)
         if !isnothing(order[1])
-            #res = y .- output.fit
-            res = StateSpaceModels.get_innovations(state_space_model)[:, 1]
+            res = ss_components["res"]
             fit_ar_model, ar_coefs, ar_intercept = fit_AR_model(res, order)
 
             initial_ar = fit_ar_model
@@ -73,14 +162,12 @@
         initial_intercept = ar_intercept
     end
 
-    pred_state = StateSpaceModels.get_predictive_state(state_space_model)
-
     if has_slope && has_level
-        initial_rws   = pred_state[2:end,1]#output.components["level"]["values"]
-        initial_slope = pred_state[2:end,2]#output.components["slope"]["values"]
+        initial_rws   = ss_components["level"]
+        initial_slope = ss_components["slope"]
         initial_rw    = zeros(length(y))
     elseif !has_slope && has_level
-        initial_rw = pred_state[2:end,1]#output.components["level"]["values"]
+        initial_rw = ss_components["level"]
         initial_rws  = zeros(length(y))
         initial_slope = zeros(length(y))
     elseif !has_slope && !has_level
@@ -90,47 +177,15 @@
     end
 
     if has_seasonality
-        initial_seasonality = zeros(length(y))
-        for t in  1:length(y)
-            initial_seasonality[t] = -sum(pred_state[t+1, end-(seasonal_period-2):end])
-        end
-        # initial_seasonality = output.components["seasonality"]["values"]
-        initial_γ, initial_γ_star = fit_harmonics(initial_seasonality, seasonal_period, stochastic)
+        initial_seasonality = ss_components["seasonality"]
+        initial_γ           = ss_components["γ"]
+        initial_γ_star      = ss_components["γ_star"]
     else
         initial_seasonality = zeros(length(y))
-        initial_γ = zeros(1)
-        initial_γ_star = zeros(1)
+        initial_γ           = zeros(1)
+        initial_γ_star      = zeros(1)
     end
 
-    # initial_values = Dict{String}{Any}()
-    # # initial_values["param"] = y
-
-    # if has_level && has_slope
-    #     initial_values["rws"] = Dict{String}{Any}()
-    #     initial_values["rws"]["values"] = output.components["level"]["values"]
-    #     #initial_values["level"]["κ"] =  haskey(output.variances,"ξ") ? output.variances["ξ"] > 0.0 ? sqrt(output.variances["ξ"]) : 0.02 : 0.02
-    #     initial_values["rws"]["κ"] = 0.02
-
-    #     initial_values["slope"] = Dict{String}{Any}()
-    #     initial_values["slope"]["values"] = has_slope ? output.components["slope"]["values"] : zeros(T)
-    #     #initial_values["slope"]["κ"] = haskey(output.variances, "ζ") ? output.variances["ζ"] > 0.0 ? sqrt(output.variances["ζ"]) : 0.02 : 0.02
-    #     initial_values["slope"]["κ"] = 0.02
-    # end
-
-    # if has_level && !has_slope
-    #     initial_values["rw"] = Dict{String}{Any}()
-    #     initial_values["rw"]["values"] = has_slope ? output.components["level"]["values"] : zeros(T)
-    #     #initial_values["slope"]["κ"] = haskey(output.variances, "ζ") ? output.variances["ζ"] > 0.0 ? sqrt(output.variances["ζ"]) : 0.02 : 0.02
-    #     initial_values["rw"]["κ"] = 0.02
-    # end
-
-    # if has_seasonality
-    #     initial_values["seasonality"] = Dict{String}{Any}()
-    #     initial_values["seasonality"]["values"] = has_seasonality ? output.components["seasonality"]["values"] : zeros(T)
-    #     #initial_values["seasonality"]["κ"] =  haskey(output.variances, "ω") ? output.variances["ω"] > 0.0 ? sqrt(output.variances["ω"]) : 0.02 : 0.02
-    #     initial_values["seasonality"]["κ"] = 0.02
-    # end
-
     selected_explanatories = missing
 
     initial_values                          = Dict()
@@ -156,18 +211,7 @@
     initial_values["ar"]["κ"]               = 0.02
 
     if has_explanatories
-
-        explanatory_idx = collect(output.components["explanatory"]["idx"])
-        # selected_explanatories = []
-        # for i in eachindex(explanatory_idx)
-        #     if explanatory_idx[i] in output.selected_variables
-        #         push!(selected_explanatories, i)
-        #     end
-        # end
-
-        initial_values["explanatories"] = output.coefs[explanatory_idx]#[selected_explanatories]
-
-        #X = X[:, selected_explanatories]
+        initial_values["explanatories"] = ss_components["explanatory"]
     end
 
     return initial_values#, X, selected_explanatories
@@ -183,7 +227,7 @@
     idx_fixed_params        = setdiff(1:num_params, idx_time_varying_params)
     T                       = length(y)
     order                   = get_AR_order(ar)
-    max_order               = has_AR(ar) ? maximum(vcat(order...)) : 0
+    max_order               = has_AR(ar) ? max_order = maximum(filter(x -> !isnothing(x), vcat(order...))) : 0
 
     initial_params = get_initial_params(y, time_varying_params, dist, seasonality)
 
@@ -323,7 +367,7 @@
     T             = length(fitted_params["param_1"])
 
     order                   = get_AR_order(ar)
-    max_order               = has_AR(ar) ? maximum(vcat(order...)) : 0#maximum(vcat(order...))
+    max_order               = has_AR(ar) ? maximum(filter(x -> !isnothing(x), vcat(order...))) : 0
 
     output_initial_values       = Dict()
     initial_time_varying_params = zeros(length(fitted_params["param_1"]), num_params) 

src/structures.jl

@@ -30,7 +30,7 @@ mutable struct GASModel
     d::Union{Float64, Missing}
     random_walk::Dict{Int64, Bool} 
     random_walk_slope::Dict{Int64, Bool} 
-    ar::Union{Dict{Int64, Int64}, Dict{Int64, Vector{Int64}}, Dict{Int64, Bool}, Dict{Int64, Any}}
+    ar::Union{Dict{Int64, Int64}, Dict{Int64, Vector{Int64}}, Dict{Int64, Bool}, Dict{Int64, Any}, Dict{Int64, Integer}}
     seasonality::Union{Dict{Int64, Int64}, Dict{Int64, Bool}}
     robust::Bool
     stochastic::Bool
@@ -40,7 +40,7 @@ mutable struct GASModel
                         d::Union{Float64, Missing},
                         random_walk::Dict{Int64, Bool}, 
                         random_walk_slope::Dict{Int64, Bool}, 
-                        ar::Union{Dict{Int64, Int64}, Dict{Int64, Vector{Int64}}, Dict{Int64, Bool}, Dict{Int64, Any}},
+                        ar::Union{Dict{Int64, Int64}, Dict{Int64, Vector{Int64}}, Dict{Int64, Bool}, Dict{Int64, Any}, Dict{Int64, Integer}},
                         seasonality::Union{Dict{Int64, Int64}, Dict{Int64, Bool}},
                         robust::Bool, stochastic::Bool)
 
@@ -77,6 +77,8 @@ mutable struct GASModel
                 end
             end
         end
+        println(ar)
+        println(seasonality)
 
         return new(dist, time_varying_params, d, random_walk, random_walk_slope, ar, seasonality, robust, stochastic)   
     end

src/utils.jl

@@ -406,8 +406,8 @@ function fit_AR_model(y::Vector{Fl}, order::Union{Vector{Int64}, Vector{Nothing}
     @variable(model, ϕ[order])
     @variable(model, y_hat[1:T])
 
-    println(order)
-    println(typeof(order))
+    # println(order)
+    # println(typeof(order))
     @constraint(model, [t = max_order+1:T], y_hat[t] == c + sum(ϕ[i]*y[t - i] for i in order))
 
     @objective(model, Min, sum((y .- y_hat).^2))