Blue format

Author: andre_ramos

paper_tests/simulation_test/evaluate_models.jl

@@ -14,7 +14,7 @@ function get_SSL_results(
     X_train::Matrix{Fl},
     inf_criteria::String,
     true_β::Vector{Fl},
-) where Fl <: AbstractFloat
+) where {Fl<:AbstractFloat}
     series_result = nothing
 
     model = StateSpaceLearning.StructuralModel(
@@ -78,7 +78,7 @@ function get_SS_res_results(
     X_train::Matrix{Fl},
     inf_criteria::String,
     true_β::Vector{Fl},
-) where Fl <: AbstractFloat
+) where {Fl<:AbstractFloat}
     py"""
     import math
     import statsmodels.api as sm
@@ -141,7 +141,9 @@ function get_SS_res_results(
     return series_result, converged
 end
 
-function get_exogenous_ss_inf_criteria(y_train::Vector{Fl}, X_train::Matrix{Fl}) where Fl <: AbstractFloat
+function get_exogenous_ss_inf_criteria(
+    y_train::Vector{Fl}, X_train::Matrix{Fl}
+) where {Fl<:AbstractFloat}
     py"""
     import math
     import statsmodels.api as sm
@@ -166,7 +168,7 @@ function get_forward_ss(
     X_train::Matrix{Fl},
     inf_criteria::String,
     true_β::Vector{Fl},
-) where Fl <: AbstractFloat
+) where {Fl<:AbstractFloat}
     best_inf_crit = Inf
     current_inf_crit = 0
     coefs = nothing

src/estimation_procedure.jl

@@ -79,7 +79,9 @@ function get_path_information_criteria(
     Lasso_y::Vector{Fl},
     information_criteria::String;
     intercept::Bool=true,
-)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}} where {Fl <: AbstractFloat, Tl <: AbstractFloat}
+)::Tuple{
+    Vector{AbstractFloat},Vector{AbstractFloat}
+} where {Fl<:AbstractFloat,Tl<:AbstractFloat}
     path_size = length(model.lambda)
     T = size(Lasso_X, 1)
     K = count(i -> i != 0, model.betas; dims=1)'
@@ -136,7 +138,9 @@ function fit_glmnet(
     information_criteria::String="aic",
     penalty_factor::Vector{Pl}=ones(size(Lasso_X, 2) - 1),
     intercept::Bool=intercept,
-)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}} where {Fl <: AbstractFloat, Tl <: AbstractFloat, Pl <: AbstractFloat}
+)::Tuple{
+    Vector{AbstractFloat},Vector{AbstractFloat}
+} where {Fl<:AbstractFloat,Tl<:AbstractFloat,Pl<:AbstractFloat}
     model = glmnet(
         Lasso_X,
         Lasso_y;
@@ -188,7 +192,9 @@ function fit_lasso(
     components_indexes::Dict{String,Vector{Int}},
     penalty_factor::Vector{Pl};
     rm_average::Bool=false,
-)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}} where {Fl <: AbstractFloat, Tl <: AbstractFloat, Pl <: AbstractFloat}
+)::Tuple{
+    Vector{AbstractFloat},Vector{AbstractFloat}
+} where {Fl<:AbstractFloat,Tl<:AbstractFloat,Pl<:AbstractFloat}
     outlier_duplicate_columns = get_outlier_duplicate_columns(
         Estimation_X, components_indexes
     )
@@ -278,7 +284,9 @@ function estimation_procedure(
     ϵ::AbstractFloat,
     penalize_exogenous::Bool,
     penalize_initial_states::Bool,
-)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}} where {Fl <: AbstractFloat, Tl <: AbstractFloat}
+)::Tuple{
+    Vector{AbstractFloat},Vector{AbstractFloat}
+} where {Fl<:AbstractFloat,Tl<:AbstractFloat}
     @assert 0 <= α <= 1 "α must be in [0, 1]"
     @assert ϵ > 0 "ϵ must be positive"
 
@@ -398,10 +406,12 @@ function estimation_procedure(
     ϵ::AbstractFloat,
     penalize_exogenous::Bool,
     penalize_initial_states::Bool,
-)::Tuple{Vector{Vector{AbstractFloat}},Vector{Vector{AbstractFloat}}} where {Fl <: AbstractFloat, Tl <: AbstractFloat}
+)::Tuple{
+    Vector{Vector{AbstractFloat}},Vector{Vector{AbstractFloat}}
+} where {Fl<:AbstractFloat,Tl<:AbstractFloat}
     coefs_vec = Vector{AbstractFloat}[]
-    ε_vec     = Vector{AbstractFloat}[]
-    N_series  = size(estimation_y, 2)
+    ε_vec = Vector{AbstractFloat}[]
+    N_series = size(estimation_y, 2)
 
     for i in 1:N_series
         coef_i, ε_i = estimation_procedure(
@@ -418,4 +428,4 @@ function estimation_procedure(
         push!(ε_vec, ε_i)
     end
     return coefs_vec, ε_vec
-end
+end

src/fit_forecast.jl

@@ -50,18 +50,26 @@ function fit!(
 
     residuals_variances = get_variances(model, estimation_ε, coefs, components_indexes)
 
-    T = typeof(model.y) <:Vector ? length(model.y) : size(model.y, 1)
+    T = typeof(model.y) <: Vector ? length(model.y) : size(model.y, 1)
 
-    ε, fitted = get_fit_and_residuals(
-        estimation_ε, coefs, model.X, valid_indexes, T
-    )
+    ε, fitted = get_fit_and_residuals(estimation_ε, coefs, model.X, valid_indexes, T)
 
-    if typeof(model.y) <:Vector
+    if typeof(model.y) <: Vector
         output = Output(coefs, ε, fitted, residuals_variances, valid_indexes, components)
     else
         output = Output[]
         for i in eachindex(coefs)
-            push!(output, Output(coefs[i], ε[i], fitted[i], residuals_variances[i], valid_indexes, components[i]))
+            push!(
+                output,
+                Output(
+                    coefs[i],
+                    ε[i],
+                    fitted[i],
+                    residuals_variances[i],
+                    valid_indexes,
+                    components[i],
+                ),
+            )
         end
     end
     return model.output = output
@@ -85,11 +93,13 @@ function forecast(
     model::StateSpaceLearningModel,
     steps_ahead::Int;
     Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, 0),
-)::Union{Matrix{<:AbstractFloat}, Vector{<:AbstractFloat}} where Fl <: AbstractFloat
-
-    exog_idx = typeof(model.output) == Output ? model.output.components["Exogenous_X"]["Indexes"] : model.output[1].components["Exogenous_X"]["Indexes"]
-    @assert length(exog_idx) ==
-        size(Exogenous_Forecast, 2) "If an exogenous matrix was utilized in the estimation procedure, it must be provided its prediction for the forecast procedure. If no exogenous matrix was utilized, Exogenous_Forecast must be missing"
+)::Union{Matrix{<:AbstractFloat},Vector{<:AbstractFloat}} where {Fl<:AbstractFloat}
+    exog_idx = if typeof(model.output) == Output
+        model.output.components["Exogenous_X"]["Indexes"]
+    else
+        model.output[1].components["Exogenous_X"]["Indexes"]
+    end
+    @assert length(exog_idx) == size(Exogenous_Forecast, 2) "If an exogenous matrix was utilized in the estimation procedure, it must be provided its prediction for the forecast procedure. If no exogenous matrix was utilized, Exogenous_Forecast must be missing"
     @assert size(Exogenous_Forecast, 1) == steps_ahead "Exogenous_Forecast must have the same number of rows as steps_ahead"
 
     Exogenous_X = model.X[:, exog_idx]
@@ -109,11 +119,14 @@ function forecast(
     )
 
     if typeof(model.output) == Output
-        return AbstractFloat.(complete_matrix[(end - steps_ahead + 1):end, :] * model.output.coefs)
+        return AbstractFloat.(
+            complete_matrix[(end - steps_ahead + 1):end, :] * model.output.coefs
+        )
     else
         prediction = Matrix{AbstractFloat}(undef, steps_ahead, length(model.output))
         for i in eachindex(model.output)
-            prediction[:, i] = complete_matrix[(end - steps_ahead + 1):end, :] * model.output[i].coefs
+            prediction[:, i] =
+                complete_matrix[(end - steps_ahead + 1):end, :] * model.output[i].coefs
         end
         return AbstractFloat.(prediction)
     end
@@ -140,16 +153,19 @@ function simulate(
     N_scenarios::Int;
     Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, 0),
     seasonal_innovation_simulation::Int=0,
-)::Union{Vector{Matrix{<:AbstractFloat}}, Matrix{<:AbstractFloat}} where Fl <: AbstractFloat
+)::Union{Vector{Matrix{<:AbstractFloat}},Matrix{<:AbstractFloat}} where {Fl<:AbstractFloat}
     @assert seasonal_innovation_simulation >= 0 "seasonal_innovation_simulation must be a non-negative integer"
     @assert seasonal_innovation_simulation >= 0 "seasonal_innovation_simulation must be a non-negative integer"
 
-    prediction = StateSpaceLearning.forecast(model, steps_ahead; Exogenous_Forecast=Exogenous_Forecast)
+    prediction = StateSpaceLearning.forecast(
+        model, steps_ahead; Exogenous_Forecast=Exogenous_Forecast
+    )
 
     is_univariate = typeof(model.output) == StateSpaceLearning.Output
 
     simulation_X = zeros(steps_ahead, 0)
-    valid_indexes = is_univariate ? model.output.valid_indexes : model.output[1].valid_indexes
+    valid_indexes =
+        is_univariate ? model.output.valid_indexes : model.output[1].valid_indexes
     components_matrix = zeros(length(valid_indexes), 0)
     N_components = 1
 
@@ -168,54 +184,112 @@ function simulate(
 
     if is_univariate
         components_matrix = hcat(components_matrix, model.output.ε[valid_indexes])
-        @assert N_components < length(model.y) // seasonal_innovation_simulation "The parameter `seasonal_innovation_simulation` is too large for the given dataset, please reduce it"
+        @assert N_components < length(model.y)//seasonal_innovation_simulation "The parameter `seasonal_innovation_simulation` is too large for the given dataset, please reduce it"
     else
         for i in eachindex(model.output)
             components_matrix = hcat(components_matrix, model.output[i].ε[valid_indexes])
         end
-        N_mv_components = N_components*length(model.output)
-        @assert N_mv_components < size(model.y, 1) // seasonal_innovation_simulation "The parameter `seasonal_innovation_simulation` is too large for the given dataset, please reduce it"
+        N_mv_components = N_components * length(model.output)
+        @assert N_mv_components < size(model.y, 1)//seasonal_innovation_simulation "The parameter `seasonal_innovation_simulation` is too large for the given dataset, please reduce it"
     end
     simulation_X = hcat(simulation_X, Matrix(1.0 * I, steps_ahead, steps_ahead))
     components_matrix += rand(Normal(0, 1), size(components_matrix)) ./ 1e9 # Make sure matrix is positive definite
 
     MV_dist_vec = Vector{MvNormal}(undef, steps_ahead)
-    o_noises = is_univariate ? zeros(steps_ahead, N_scenarios) : [zeros(steps_ahead, N_scenarios) for _ in 1:length(model.output)]
+    o_noises = if is_univariate
+        zeros(steps_ahead, N_scenarios)
+    else
+        [zeros(steps_ahead, N_scenarios) for _ in 1:length(model.output)]
+    end
 
     if seasonal_innovation_simulation == 0
         ∑ = cov(components_matrix)
         for i in 1:steps_ahead
-            MV_dist_vec[i] = is_univariate ? MvNormal(zeros(N_components), ∑) : MvNormal(zeros(N_mv_components), ∑)
+            MV_dist_vec[i] = if is_univariate
+                MvNormal(zeros(N_components), ∑)
+            else
+                MvNormal(zeros(N_mv_components), ∑)
+            end
         end
 
         if model.outlier
             if is_univariate
-                o_noises = rand(Normal(0, std(model.output.components["o"]["Coefs"])), steps_ahead, N_scenarios)
+                o_noises = rand(
+                    Normal(0, std(model.output.components["o"]["Coefs"])),
+                    steps_ahead,
+                    N_scenarios,
+                )
             else
-                o_noises = [rand(Normal(0, std(model.output[i].components["o"]["Coefs"])), steps_ahead, N_scenarios) for i in eachindex(model.output)]
+                o_noises = [
+                    rand(
+                        Normal(0, std(model.output[i].components["o"]["Coefs"])),
+                        steps_ahead,
+                        N_scenarios,
+                    ) for i in eachindex(model.output)
+                ]
             end
         end
     else
         start_seasonal_term = (size(components_matrix, 1) % seasonal_innovation_simulation)
         for i in 1:steps_ahead
-            ∑ = cov(components_matrix[i + start_seasonal_term:seasonal_innovation_simulation:end, :])
-            MV_dist_vec[i] = is_univariate ? MvNormal(zeros(N_components), ∑) : MvNormal(zeros(N_mv_components), ∑)
+            ∑ = cov(
+                components_matrix[
+                    (i + start_seasonal_term):seasonal_innovation_simulation:end, :,
+                ],
+            )
+            MV_dist_vec[i] = if is_univariate
+                MvNormal(zeros(N_components), ∑)
+            else
+                MvNormal(zeros(N_mv_components), ∑)
+            end
             if is_univariate
-                model.outlier ? o_noises[i, :] = rand(Normal(0, std(model.output.components["o"]["Coefs"][i + start_seasonal_term:seasonal_innovation_simulation:end])), N_scenarios) : nothing
+                if model.outlier
+                    o_noises[i, :] = rand(
+                    Normal(
+                        0,
+                        std(
+                            model.output.components["o"]["Coefs"][(i + start_seasonal_term):seasonal_innovation_simulation:end],
+                        ),
+                    ),
+                    N_scenarios,
+                )
+                else
+                    nothing
+                end
             else
                 for j in eachindex(model.output)
-                    model.outlier ? o_noises[j][i, :] = rand(Normal(0, std(model.output[j].components["o"]["Coefs"][i + start_seasonal_term:seasonal_innovation_simulation:end])), N_scenarios) : nothing
+                    if model.outlier
+                        o_noises[j][i, :] = rand(
+                        Normal(
+                            0,
+                            std(
+                                model.output[j].components["o"]["Coefs"][(i + start_seasonal_term):seasonal_innovation_simulation:end],
+                            ),
+                        ),
+                        N_scenarios,
+                    )
+                    else
+                        nothing
+                    end
                 end
             end
-        end 
+        end
+    end
 
+    simulation = if is_univariate
+        AbstractFloat.(hcat([prediction for _ in 1:N_scenarios]...))
+    else
+        [
+        AbstractFloat.(hcat([prediction[:, i] for _ in 1:N_scenarios]...)) for
+        i in eachindex(model.output)
+    ]
     end
-    
-    simulation = is_univariate ? AbstractFloat.(hcat([prediction for _ in 1:N_scenarios]...)) : [AbstractFloat.(hcat([prediction[:, i] for _ in 1:N_scenarios]...)) for i in eachindex(model.output)]
     if is_univariate
         fill_simulation!(simulation, MV_dist_vec, o_noises, simulation_X)
     else
-        fill_simulation!(simulation, MV_dist_vec, o_noises, simulation_X, length(model_innovations))
+        fill_simulation!(
+            simulation, MV_dist_vec, o_noises, simulation_X, length(model_innovations)
+        )
         simulation = Vector{Matrix{<:AbstractFloat}}(simulation)
     end
 

src/information_criteria.jl

@@ -16,7 +16,7 @@
 """
 function get_information(
     T::Int, K::Int, ε::Vector{Fl}; information_criteria::String="aic"
-)::AbstractFloat where Fl <: AbstractFloat
+)::AbstractFloat where {Fl<:AbstractFloat}
     if information_criteria == "bic"
         return T * log(var(ε)) + K * log(T)
     elseif information_criteria == "aic"