Merge pull request #11 from LAMPSPUC/format_update

Format update

README.md

@@ -38,8 +38,8 @@ prediction = StateSpaceLearning.forecast(output, 12) #Gets a 12 steps ahead pred
 ## Fit Arguments
 
 * `y::Vector{Fl}`: Vector of data.
-* model_input::Dict: Dictionary containing the model input parameters (default: Dict("level" => true, "stochastic_level" => true, "trend" => true, "stochastic_trend" => true, "seasonal" => true, "stochastic_seasonal" => true, "freq_seasonal" => 12)).
-* estimation_input::Dict: Dictionary containing the estimation input parameters (default: Dict("α" => 0.1, "information_criteria" => "aic", ψ => 0.05, "penalize_exogenous" => true, "penalize_initial_states" => true)).
+* `model_input::Dict`: Dictionary containing the model input parameters (default: Dict("level" => true, "stochastic_level" => true, "trend" => true, "stochastic_trend" => true, "seasonal" => true, "stochastic_seasonal" => true, "freq_seasonal" => 12)).
+* `estimation_input::Dict`: Dictionary containing the estimation input parameters (default: Dict("α" => 0.1, "information_criteria" => "aic", ψ => 0.05, "penalize_exogenous" => true, "penalize_initial_states" => true)).
 * `Exogenous_X::Union{Matrix{Fl}, Missing}`: Exogenous variables matrix (default: missing).
 * `outlier::Bool`: Flag for considering outlier component (default: true).
 * `ζ_ω_threshold::Int64`: ζ_ω_threshold parameter (default: 12).

src/estimation_procedure/default_estimation_procedure.jl

@@ -65,15 +65,15 @@ end
 """
 function get_path_information_criteria(model::GLMNetPath, Estimation_X::Matrix{Tl}, estimation_y::Vector{Fl}, information_criteria::String; intercept::Bool = true)::Tuple{Vector{Float64}, Vector{Float64}} where {Tl, Fl}
     path_size = length(model.lambda)
-    T, p      = size(Estimation_X)
+    T         = size(Estimation_X, 1)
     K         = count(i->i != 0, model.betas; dims = 1)'
 
     method_vec = Vector{Float64}(undef, path_size)
     for i in 1:path_size
         fit = Estimation_X*model.betas[:, i] .+ model.a0[i]
         ϵ   = estimation_y - fit
 
-        method_vec[i] = get_information(T, K[i], ϵ; information_criteria = information_criteria, p = p)
+        method_vec[i] = get_information(T, K[i], ϵ; information_criteria = information_criteria)
     end
 
     best_model_idx = argmin(method_vec)

src/information_criteria.jl

@@ -8,22 +8,19 @@
     - `T::Int64`: Number of observations.
     - `K::Int64`: Number of selected predictors.
     - `ϵ::Vector{Float64}`: Vector of residuals.
-    - `information_criteria::String`: Method for hyperparameter selection (default: "bic").
+    - `information_criteria::String`: Method for hyperparameter selection (default: "aic").
     - `p::Int64`: Number of total predictors (default: 0).
 
     # Returns
     - `Float64`: Information criterion value.
 
 """
-function get_information(T::Int64, K::Int64, ϵ::Vector{Float64}; information_criteria::String = "bic", p::Int64 = 0)::Float64
+function get_information(T::Int64, K::Int64, ϵ::Vector{Float64}; information_criteria::String = "aic")::Float64
     if information_criteria == "bic"
         return T*log(var(ϵ)) + K*log(T)
     elseif information_criteria == "aic"
         return 2*K + T*log(var(ϵ))
     elseif information_criteria == "aicc"
         return 2*K + T*log(var(ϵ))  + ((2*K^2 +2*K)/(T - K - 1))
-    elseif information_criteria == "EBIC"
-        EBIC_comb_term = (K <= 1 || p == K) ? 0 : 2*(sum(log(j) for j in 1:p) - (sum(log(j) for j in 1:K) + sum(log(j) for j in 1:(p-K))))
-        return T*log(var(ϵ)) + K*log(T) + EBIC_comb_term
     end
 end

src/models/unobserved_components.jl

@@ -133,7 +133,7 @@ end
 """
 function create_initial_states_Matrix(T::Int64, freq_seasonal::Int64, steps_ahead::Int64, trend::Bool, seasonal::Bool)::Matrix
 
-    initial_states_matrix = ones(T+steps_ahead)
+    initial_states_matrix = ones(T+steps_ahead, 1)
     trend ? initial_states_matrix = hcat(initial_states_matrix, vcat([0], collect(1:T+steps_ahead-1))) : nothing
 
     if seasonal

test/StateSpaceLearning.jl

@@ -13,13 +13,13 @@
     @test all(isnan.(output2.ϵ[10:20]))
     @test !all(isnan.(output2.fitted[10:20]))
 
-    output3 = StateSpaceLearning.fit_model(y1; stabilize_ζ = 1)
-    @test length(output3.coefs) == length(output1.coefs) - 2
+    output3 = StateSpaceLearning.fit_model(y1; ζ_ω_threshold = 1)
+    @test length(output3.coefs) - 22 == length(output1.coefs)
 
-    @test_throws AssertionError StateSpaceLearning.fit_model(y1; s = 200)
-    @test_throws ErrorException StateSpaceLearning.fit_model(y1; model_type = "none")
-    @test_throws AssertionError StateSpaceLearning.fit_model(y1; α = -0.1)
-    @test_throws AssertionError StateSpaceLearning.fit_model(y1; α = 1.1)
+    @test_throws AssertionError StateSpaceLearning.fit_model(y1; model_input = Dict("stochastic_level" => true, "trend" => true, "stochastic_trend" => true, "seasonal" => true, "stochastic_seasonal" => true, "freq_seasonal" => 1000))
+
+    @test_throws AssertionError StateSpaceLearning.fit_model(y1; estimation_input = Dict("α" => -0.1, "information_criteria" => "aic", "ψ" => 0.05, "penalize_exogenous" => true, "penalize_initial_states" => true))
+    @test_throws AssertionError StateSpaceLearning.fit_model(y1; estimation_input = Dict("α" => 1.1, "information_criteria" => "aic", "ψ" => 0.05, "penalize_exogenous" => true, "penalize_initial_states" => true))
 
 end
 
@@ -33,7 +33,6 @@ end
     output2 = StateSpaceLearning.fit_model(y2; Exogenous_X = rand(100, 3))
     @test length(StateSpaceLearning.forecast(output2, 10; Exogenous_Forecast = rand(10, 3))) == 10
 
-    @test_throws AssertionError StateSpaceLearning.forecast(output1, -1)
     @test_throws AssertionError StateSpaceLearning.forecast(output1, 10; Exogenous_Forecast = rand(5, 3))
     @test_throws AssertionError StateSpaceLearning.forecast(output2, 10)
     @test_throws AssertionError StateSpaceLearning.forecast(output2, 10; Exogenous_Forecast = rand(5, 3))

test/estimation_procedure/default_estimation_procedure.jl

@@ -0,0 +1,112 @@
+Random.seed!(1234)
+Estimation_X = rand(30, 3)
+estimation_y = rand(30)
+α = 0.5
+penalty_factor = ones(3)
+@testset "Function: get_path_information_criteria" begin
+    intercept1 = true
+    intercept2 = false
+
+    model1 = glmnet(Estimation_X, estimation_y, alpha = α, penalty_factor = penalty_factor, intercept = intercept1, dfmax=size(Estimation_X, 2), lambda_min_ratio=0.001)
+    coefs1, ϵ1 = StateSpaceLearning.get_path_information_criteria(model1, Estimation_X, estimation_y, "aic"; intercept = intercept1)
+    @test length(coefs1) == 4
+    @test coefs1[1] != 0
+    @test all(coefs1[2:end] .== 0)
+    @test length(ϵ1) == 30
+
+    model2 = glmnet(Estimation_X, estimation_y, alpha = α, penalty_factor = penalty_factor, intercept = intercept2, dfmax=size(Estimation_X, 2), lambda_min_ratio=0.001)
+    coefs2, ϵ2 = StateSpaceLearning.get_path_information_criteria(model2, Estimation_X, estimation_y, "aic"; intercept = intercept2)
+    @test length(coefs2) == 3
+    @test all(coefs2 .== 0)
+    @test length(ϵ2) == 30
+end
+
+@testset "Function: fit_glmnet" begin
+    coefs, ϵ = StateSpaceLearning.fit_glmnet(Estimation_X, estimation_y, α; information_criteria="aic", penalty_factor=penalty_factor, intercept = true)
+    @test length(coefs) == 4
+    @test length(ϵ) == 30
+end
+
+@testset "Function: fit_lasso" begin
+    Random.seed!(1234)
+    Exogenous_X = hcat(rand(10, 3), vcat(zeros(3), ones(1), zeros(6)))
+    Basic_Structural = Dict("stochastic_level" => true, "trend" => true, "stochastic_trend" => true, "seasonal" => true, "stochastic_seasonal" => true, "freq_seasonal" => 2)
+
+    components_indexes = StateSpaceLearning.get_components_indexes(10, Exogenous_X, Basic_Structural, true, 0)
+
+    Estimation_X = StateSpaceLearning.create_X_unobserved_components(Basic_Structural, Exogenous_X, true, 0, 10)
+    estimation_y = Estimation_X*rand(size(Estimation_X, 2)) + rand(10)
+
+    coefs1, ϵ1 = StateSpaceLearning.fit_lasso(Estimation_X, estimation_y, 0.1, "aic", true, components_indexes; intercept = true)
+    @test length(coefs1) == 43
+    @test length(ϵ1) == 10
+
+    coefs2, ϵ2 = StateSpaceLearning.fit_lasso(Estimation_X, estimation_y, 0.1, "aic", true, components_indexes; intercept = false)
+    @test coefs2[1] == mean(estimation_y)
+    @test length(coefs2) == 43
+    @test length(ϵ2) == 10
+
+    coefs3, ϵ3 = StateSpaceLearning.fit_lasso(Estimation_X, estimation_y, 0.1, "aic", false, components_indexes; intercept = true)
+    @test coefs3[components_indexes["o"][4]] == 0
+    @test all(coefs3[components_indexes["Exogenous_X"]] .!= 0)
+    @test length(coefs3) == 43
+    @test length(ϵ3) == 10
+
+    coefs4, ϵ4 = StateSpaceLearning.fit_lasso(Estimation_X, estimation_y, 0.1, "aic", true, components_indexes; penalty_factor = vcat(ones(1), ones(size(Estimation_X,2) - 2).*Inf), intercept = true)
+    @test all(coefs4[3:end] .== 0)
+    @test length(coefs4) == 43
+    @test length(ϵ4) == 10
+end
+
+@testset "Function: default_estimation_procedure" begin
+    Random.seed!(1234)
+    Exogenous_X = hcat(rand(10, 3), vcat(ones(3), zeros(1), ones(6)))
+    Basic_Structural = Dict("stochastic_level" => true, "trend" => true, "stochastic_trend" => true, "seasonal" => true, "stochastic_seasonal" => true, "freq_seasonal" => 2)
+
+    components_indexes = StateSpaceLearning.get_components_indexes(10, Exogenous_X, Basic_Structural, true, 0)
+
+    Estimation_X = StateSpaceLearning.create_X_unobserved_components(Basic_Structural, Exogenous_X, true, 0, 10)
+
+    estimation_y = Estimation_X*rand(size(Estimation_X, 2)) + rand(10).*5
+
+    estimation_input1 = Dict("α" => 0.1, "information_criteria" => "aic", "ψ" => 0.05, "penalize_exogenous" => true, "penalize_initial_states" => true)
+    coefs1, ϵ1 = StateSpaceLearning.default_estimation_procedure(Estimation_X, estimation_y, components_indexes, estimation_input1)
+    @test length(coefs1) == 43
+    @test length(ϵ1) == 10
+
+    estimation_input2 = Dict("α" => 0.1, "information_criteria" => "aic", "ψ" => 0.05, "penalize_exogenous" => true, "penalize_initial_states" => false)
+    coefs2, ϵ2 = StateSpaceLearning.default_estimation_procedure(Estimation_X, estimation_y, components_indexes, estimation_input2)
+    @test length(coefs2) == 43
+    @test length(ϵ2) == 10
+    @test all(coefs2[components_indexes["initial_states"][2:end] .- 1] .!= 0)
+end
+
+@testset "Function: get_dummy_indexes" begin 
+    Exogenous_X1 = hcat(rand(10, 3), vcat(zeros(3), ones(1), zeros(6)))
+    Exogenous_X2 = hcat(rand(10, 3))
+
+    dummy_indexes1 = StateSpaceLearning.get_dummy_indexes(Exogenous_X1)
+    @test dummy_indexes1 == [4]
+
+    dummy_indexes2 = StateSpaceLearning.get_dummy_indexes(Exogenous_X2)
+    @test dummy_indexes2 == []
+end
+
+@testset "Function: get_outlier_duplicate_columns" begin
+    Random.seed!(1234)
+    Exogenous_X1 = hcat(rand(10, 3), vcat(zeros(3), ones(1), zeros(6)))
+    Exogenous_X2 = rand(10, 3)
+
+    Basic_Structural = Dict("stochastic_level" => true, "trend" => true, "stochastic_trend" => true, "seasonal" => true, "stochastic_seasonal" => true, "freq_seasonal" => 2)
+
+    components_indexes1 = StateSpaceLearning.get_components_indexes(10, Exogenous_X1, Basic_Structural, true, 0)
+    components_indexes2 = StateSpaceLearning.get_components_indexes(10, Exogenous_X2, Basic_Structural, true, 0)
+
+    Estimation_X1 = StateSpaceLearning.create_X_unobserved_components(Basic_Structural, Exogenous_X1, true, 0, 10)
+    outlier_duplicate_columns1 = StateSpaceLearning.get_outlier_duplicate_columns(Estimation_X1, components_indexes1)
+    @test outlier_duplicate_columns1 == [32]
+
+    Estimation_X2 = StateSpaceLearning.create_X_unobserved_components(Basic_Structural, Exogenous_X2, true, 0, 10)
+    outlier_duplicate_columns2 = StateSpaceLearning.get_outlier_duplicate_columns(Estimation_X2, components_indexes2)
+    @test outlier_duplicate_columns2 == []
+end

test/information_criteria.jl

@@ -0,0 +1,11 @@
+@testset "Function: get_information"  begin 
+    ϵ = [1.1, 2.2, 3.3, 4.4, 5.5]
+    T = 5
+    K = 3
+    bic = StateSpaceLearning.get_information(T, K, ϵ; information_criteria = "bic")
+    aic = StateSpaceLearning.get_information(T, K, ϵ; information_criteria = "aic")
+    aicc = StateSpaceLearning.get_information(T, K, ϵ; information_criteria = "aicc")
+    @test round(bic, digits = 5) == 10.36287
+    @test round(aic, digits = 5) == 11.53456
+    @test round(aicc, digits = 5) == 35.53456
+end