Merge pull request #1 from LAMPSPUC/first_commit

🚀 first commit

.github/workflows/ci.yml

@@ -0,0 +1,34 @@
+# This is a basic workflow to help you get started with Actions
+
+name: CI
+
+# Controls when the workflow will run
+on:
+  # Triggers the workflow on push or pull request events but only for the "main" branch
+  push:
+    branches: [ "main" ]
+  pull_request:
+    branches: [ "main" ]
+
+  # Allows you to run this workflow manually from the Actions tab
+  workflow_dispatch:
+
+# A workflow run is made up of one or more jobs that can run sequentially or in parallel
+jobs:
+  test:
+    runs-on: ubuntu-latest
+
+    steps:
+      - name: Checkout repository
+        uses: actions/checkout@v2
+
+      - name: Set up Julia 1.8
+        uses: julia-actions/setup-julia@v1
+        with:
+          julia-version: 1.8
+          project: .
+
+      - name: Install dependencies and run tests
+        run: |
+          julia --color=yes -e 'import Pkg; Pkg.build()'
+          julia --color=yes -e 'import Pkg; Pkg.activate("."); Pkg.instantiate(); Pkg.test()'

Project.toml

@@ -0,0 +1,12 @@
+name = "StateSpaceLearning"
+uuid = "971c4b7c-2c4e-4bac-8525-e842df3cde7b"
+authors = ["andreramosfc <[email protected]>"]
+version = "0.1.0"
+
+[deps]
+Coverage = "a2441757-f6aa-5fb2-8edb-039e3f45d037"
+GLMNet = "8d5ece8b-de18-5317-b113-243142960cc6"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/StateSpaceLearning.jl

@@ -0,0 +1,61 @@
+module StateSpaceLearning
+
+const AVAILABLE_MODELS = ["Basic Structural", "Local Linear Trend", "Local Level"]
+const AVAILABLE_ESTIMATION_PROCEDURES = ["Lasso", "AdaLasso"]
+const AVAILABLE_HYPERPARAMETER_SELECTION = ["aic", "bic", "aicc", "EBIC"]
+
+using LinearAlgebra, Statistics, GLMNet
+
+include("structs.jl")
+include("model_utils.jl")
+include("estimation_procedure/information_criteria.jl")
+include("estimation_procedure/lasso.jl")
+include("estimation_procedure/adalasso.jl")
+include("estimation_procedure/estimation_utils.jl")
+include("utils.jl")
+
+export fit_model, forecast
+
+function fit_model(y::Vector{Fl}; model_type::String="Basic Structural", Exogenous_X::Union{Matrix{Fl}, Missing}=missing,
+                    estimation_procedure::String="AdaLasso", s::Int64=12, outlier::Bool=false, stabilize_ζ::Int64=0,
+                    α::Float64=0.1, hyperparameter_selection::String="aic", adalasso_coef::Float64=0.1, select_exogenous::Bool=true)::Output where Fl
+
+    T = length(y)
+    @assert T > s "Time series must be longer than the seasonal period"
+    @assert (model_type in AVAILABLE_MODELS) "Unavailable Model"
+    @assert (estimation_procedure in AVAILABLE_ESTIMATION_PROCEDURES) "Unavailable estimation procedure"
+    @assert (hyperparameter_selection in AVAILABLE_HYPERPARAMETER_SELECTION) "Unavailable hyperparameter selection method"
+    @assert 0 < α <= 1 "α must be in (0, 1], Lasso.jl cannot handle α = 0"
+
+    valid_indexes = findall(i -> !isnan(i), y)
+    estimation_y  = y[valid_indexes]
+
+    Exogenous_X = ismissing(Exogenous_X) ? zeros(T, 0) : Exogenous_X
+
+    X = create_X(model_type, T, s, Exogenous_X, outlier, stabilize_ζ)
+    Estimation_X = X[valid_indexes, :]
+
+    components_indexes  = get_components_indexes(T, s, Exogenous_X, outlier, model_type, stabilize_ζ)
+
+    coefs, estimation_ϵ = fit_estimation_procedure(estimation_procedure, Estimation_X, estimation_y, α, hyperparameter_selection, components_indexes, adalasso_coef, select_exogenous)
+
+    components          = build_components(Estimation_X, coefs, components_indexes)
+
+    residuals_variances = get_variances(estimation_ϵ, coefs, components_indexes)
+
+    ϵ, fitted = build_complete_variables(estimation_ϵ, coefs, X, valid_indexes, T)
+
+    return Output(model_type, X, coefs, ϵ, fitted, components, residuals_variances, s, T, outlier, valid_indexes, stabilize_ζ)
+end
+
+function forecast(output::Output, steps_ahead::Int64; Exogenous_Forecast::Union{Matrix{Fl}, Missing}=missing)::Vector{Float64} where Fl
+    @assert steps_ahead > 0 "steps_ahead must be a positive integer"
+    Exogenous_Forecast = ismissing(Exogenous_Forecast) ? zeros(steps_ahead, 0) : Exogenous_Forecast
+
+    @assert length(output.components["Exogenous_X"]["Indexes"]) == 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"
+
+    return forecast_model(output, steps_ahead, Exogenous_Forecast)
+end
+
+end # module StateSpaceLearning

src/estimation_procedure/adalasso.jl

@@ -0,0 +1,24 @@
+function fit_adalasso(Estimation_X::Matrix{Tl}, estimation_y::Vector{Fl}, α::Float64, 
+                        hyperparameter_selection::String, 
+                        components_indexes::Dict{String, Vector{Int64}},
+                        adalasso_coef::Float64, select_exogenous::Bool)::Tuple{Vector{Float64}, Vector{Float64}} where {Tl, Fl}
+
+    penalty_factor = ones(size(Estimation_X, 2) - 1); penalty_factor[components_indexes["initial_states"][2:end] .- 1] .= 0
+    coefs, _  = fit_lasso(Estimation_X, estimation_y, α, hyperparameter_selection, select_exogenous, components_indexes; penalty_factor = penalty_factor, intercept = false)
+
+    #AdaLasso per component
+    penalty_factor = zeros(size(Estimation_X, 2) - 1)
+    for key in keys(components_indexes)
+        if key != "initial_states" && key != "μ₁"
+            component = components_indexes[key]
+            if key != "Exogenous_X" && key != "o" && !(key in ["ν₁", "γ₁"])
+                κ = count(i -> i == 0, coefs[component]) < 1 ? 0 : std(coefs[component])
+                penalty_factor[component .- 1] .= (1 / (κ + adalasso_coef))
+            else
+                penalty_factor[component .- 1]  = (1 ./ (abs.(coefs[component]) .+ adalasso_coef))
+            end
+        end
+    end 
+    return fit_lasso(Estimation_X, estimation_y, α, hyperparameter_selection, select_exogenous, components_indexes; penalty_factor=penalty_factor)
+end
+

src/estimation_procedure/estimation_utils.jl

@@ -0,0 +1,36 @@
+function get_dummy_indexes(Exogenous_X::Matrix{Fl}) where{Fl}
+
+    T, p = size(Exogenous_X)
+    dummy_indexes = []
+
+    for i in 1:p
+        if count(iszero.(Exogenous_X[:, i])) == T - 1
+            push!(dummy_indexes, findfirst(i -> i != 0.0, Exogenous_X[:, i]))
+        end
+    end
+
+    return dummy_indexes
+end
+
+function get_outlier_duplicate_columns(Estimation_X::Matrix{Tl}, components_indexes::Dict{String, Vector{Int64}}) where{Tl}
+    o_indexes = components_indexes["o"]
+    exogenous_indexes = components_indexes["Exogenous_X"]
+
+    dummy_indexes = get_dummy_indexes(Estimation_X[:, exogenous_indexes])
+
+    return o_indexes[dummy_indexes] .- 1
+end
+
+function fit_estimation_procedure(estimation_procedure::String, Estimation_X::Matrix{Tl}, estimation_y::Vector{Fl}, α::Float64, 
+                                    hyperparameter_selection::String, components_indexes::Dict{String, Vector{Int64}}, 
+                                    adalasso_coef::Float64, select_exogenous::Bool)::Tuple{Vector{Float64}, Vector{Float64}} where {Tl, Fl}
+
+    estimation_arguments     = Dict("Lasso" => (Estimation_X, estimation_y, α, hyperparameter_selection, select_exogenous, components_indexes),
+                                "AdaLasso" => (Estimation_X, estimation_y, α, hyperparameter_selection, 
+                                               components_indexes, adalasso_coef, select_exogenous))
+
+    available_estimation = Dict("Lasso" => fit_lasso,
+                                "AdaLasso" => fit_adalasso)
+
+    return available_estimation[estimation_procedure](estimation_arguments[estimation_procedure]...)
+end

src/estimation_procedure/information_criteria.jl

@@ -0,0 +1,12 @@
+function get_information(T::Int64, K::Int64, ϵ::Vector{Float64}; hyperparameter_selection::String = "bic", p::Int64 = 0)::Float64
+    if hyperparameter_selection == "bic"
+        return T*log(var(ϵ)) + K*log(T)
+    elseif hyperparameter_selection == "aic"
+        return 2*K + T*log(var(ϵ))
+    elseif hyperparameter_selection == "aicc"
+        return 2*K + T*log(var(ϵ))  + ((2*K^2 +2*K)/(T - K - 1))
+    elseif hyperparameter_selection == "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/estimation_procedure/lasso.jl

@@ -0,0 +1,37 @@
+function get_path_information_criteria(model::GLMNetPath, Estimation_X::Matrix{Tl}, estimation_y::Vector{Fl}, hyperparameter_selection::String; intercept::Bool = true)::Tuple{Vector{Float64}, Vector{Float64}} where {Tl, Fl}
+    path_size = length(model.lambda)
+    T, p      = size(Estimation_X)
+    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], ϵ; hyperparameter_selection = hyperparameter_selection, p = p)
+    end
+
+    best_model_idx = argmin(method_vec)
+    coefs = intercept ? vcat(model.a0[best_model_idx], model.betas[:, best_model_idx]) : model.betas[:, best_model_idx]
+    fit   = intercept ? hcat(ones(T), Estimation_X)*coefs : Estimation_X*coefs
+    ϵ   = estimation_y - fit
+    return coefs, ϵ
+end
+
+function fit_glmnet(Estimation_X::Matrix{Tl}, estimation_y::Vector{Fl}, α::Float64; hyperparameter_selection::String = "aic", penalty_factor::Vector{Float64}=ones(size(Estimation_X,2) - 1), intercept::Bool = intercept)::Tuple{Vector{Float64}, Vector{Float64}} where {Tl, Fl}
+    model = glmnet(Estimation_X, estimation_y, alpha = α, penalty_factor = penalty_factor, intercept = intercept, dfmax=size(Estimation_X, 2), lambda_min_ratio=0.001)
+    return get_path_information_criteria(model, Estimation_X, estimation_y, hyperparameter_selection; intercept = intercept)
+end
+
+function fit_lasso(Estimation_X::Matrix{Tl}, estimation_y::Vector{Fl}, α::Float64, hyperparameter_selection::String, select_exogenous::Bool, components_indexes::Dict{String, Vector{Int64}}; penalty_factor::Vector{Float64}=ones(size(Estimation_X,2) - 1), intercept::Bool = true)::Tuple{Vector{Float64}, Vector{Float64}} where {Tl, Fl}
+
+    outlier_duplicate_columns = get_outlier_duplicate_columns(Estimation_X, components_indexes)
+    penalty_factor[outlier_duplicate_columns] .= Inf
+
+    !select_exogenous ? penalty_factor[components_indexes["Exogenous_X"] .- 1] .= 0 : nothing
+    mean_y = mean(estimation_y); Lasso_y = intercept ? estimation_y : estimation_y .- mean_y
+
+    coefs, ϵ =  fit_glmnet(Estimation_X[:, 2:end], Lasso_y, α; hyperparameter_selection=hyperparameter_selection, penalty_factor=penalty_factor, intercept = intercept)
+    return !intercept ? (vcat(mean_y, coefs), ϵ) : (coefs, ϵ)
+
+end

src/model_utils.jl

@@ -0,0 +1,130 @@
+ξ_size(T::Int64)::Int64 = T-2
+
+ζ_size(T::Int64, stabilize_ζ::Int64)::Int64 = T-stabilize_ζ-1
+
+ω_size(T::Int64, s::Int64)::Int64 = T-s
+
+o_size(T::Int64)::Int64 = T
+
+function create_ξ(T::Int64, steps_ahead::Int64)::Matrix
+    ξ_matrix = Matrix{Float64}(undef, T+steps_ahead, ξ_size(T))
+    for t in 1:T+steps_ahead
+        ξ_matrix[t, :] = t < T ? vcat(ones(t-1), zeros(T-t-1)) : ξ_matrix[T, :] = ones(ξ_size(T))
+    end
+
+    return ξ_matrix
+end
+
+function create_ζ(T::Int64, steps_ahead::Int64, stabilize_ζ::Int64)::Matrix
+    ζ_matrix = Matrix{Float64}(undef, T+steps_ahead, ζ_size(T, stabilize_ζ) + stabilize_ζ)
+    for t in 1:T+steps_ahead
+        ζ_matrix[t, :] = t < T ? vcat(collect(t:-1:2), zeros(T-t)) : collect(t:-1:t-ζ_size(T, stabilize_ζ)-stabilize_ζ+1)
+    end
+    return ζ_matrix[:, 1:end-stabilize_ζ]
+end
+
+function create_ω(T::Int64, s::Int64, steps_ahead::Int64)::Matrix
+    ω_matrix_size = ω_size(T, s)
+    ω_matrix = Matrix{Float64}(undef, T+steps_ahead, ω_matrix_size)
+    for t in 1:T+steps_ahead
+        ωₜ_coefs = zeros(ω_matrix_size)
+        Mₜ = Int64(floor((t-1)/s))
+        lag₁ = [t - j*s for j in 1:Mₜ]
+        lag₂ = [t - j*s - 1 for j in 0:Mₜ]
+        ωₜ_coefs[lag₁[lag₁.<=ω_matrix_size]] .= 1
+        ωₜ_coefs[lag₂[0 .< lag₂.<=ω_matrix_size]] .= -1
+        ω_matrix[t, :] = ωₜ_coefs
+    end
+    return ω_matrix
+end
+
+function create_o_matrix(T::Int64, steps_ahead::Int64)::Matrix
+    return vcat(Matrix(1.0 * I, T, T), zeros(steps_ahead, T))
+end
+
+function create_initial_states_Matrix(T::Int64, s::Int64, steps_ahead::Int64, model_type::String)::Matrix
+    μ₀_coefs = ones(T+steps_ahead)
+    ν₀_coefs = collect(1:T+steps_ahead)
+
+    γ₀_coefs = zeros(T+steps_ahead, s)
+    for t in 1:T+steps_ahead
+        γ₀_coefs[t, t % s == 0 ? s : t % s] = 1.0
+    end
+
+    if model_type == "Basic Structural"
+        return hcat(μ₀_coefs, ν₀_coefs, γ₀_coefs)
+    elseif model_type == "Local Linear Trend"
+        return hcat(μ₀_coefs, ν₀_coefs)
+    elseif model_type == "Local Level"
+        return hcat(μ₀_coefs)
+    end
+
+end
+
+function create_X(model_type::String, T::Int64, s::Int64, Exogenous_X::Matrix{Fl}, outlier::Bool, stabilize_ζ::Int64,
+                  steps_ahead::Int64=0, Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, size(Exogenous_X, 2))) where Fl
+
+    ξ_matrix = create_ξ(T, steps_ahead)
+    ζ_matrix = create_ζ(T, steps_ahead, stabilize_ζ)
+    ω_matrix = create_ω(T, s, steps_ahead)
+    o_matrix = outlier ? create_o_matrix(T, steps_ahead) : zeros(T+steps_ahead, 0)
+
+    initial_states_matrix = create_initial_states_Matrix(T, s, steps_ahead, model_type)
+    if model_type == "Basic Structural"
+        return hcat(initial_states_matrix, ξ_matrix, ζ_matrix, ω_matrix, o_matrix, vcat(Exogenous_X, Exogenous_Forecast))
+    elseif model_type == "Local Level"
+        return hcat(initial_states_matrix, ξ_matrix, o_matrix, vcat(Exogenous_X, Exogenous_Forecast))
+    elseif model_type == "Local Linear Trend"
+        return hcat(initial_states_matrix, ξ_matrix, ζ_matrix, o_matrix, vcat(Exogenous_X, Exogenous_Forecast))
+    end
+
+end
+
+function get_components_indexes(T::Int64, s::Int64, Exogenous_X::Matrix{Fl}, outlier::Bool, model_type::String, stabilize_ζ::Int64)::Dict where Fl
+    μ₁_indexes = [1]
+    ν₁_indexes = model_type in ["Local Linear Trend", "Basic Structural"] ? [2] : Int64[]
+    γ₁_indexes = model_type == "Basic Structural" ? collect(3:2+s) : Int64[]
+
+    initial_states_indexes = μ₁_indexes
+    !isempty(ν₁_indexes) ?  initial_states_indexes = vcat(initial_states_indexes, ν₁_indexes) : nothing
+    !isempty(γ₁_indexes) ?  initial_states_indexes = vcat(initial_states_indexes, γ₁_indexes) : nothing 
+
+    FINAL_INDEX = initial_states_indexes[end]
+
+    ξ_indexes   = collect(FINAL_INDEX + 1:FINAL_INDEX + ξ_size(T))
+    FINAL_INDEX = ξ_indexes[end]
+
+    ζ_indexes   = !isempty(ν₁_indexes) ? collect(FINAL_INDEX + 1:FINAL_INDEX + ζ_size(T, stabilize_ζ)) : Int64[]
+    FINAL_INDEX = !isempty(ν₁_indexes) ? ζ_indexes[end] : FINAL_INDEX
+
+    ω_indexes   = !isempty(γ₁_indexes) ? collect(FINAL_INDEX + 1:FINAL_INDEX + ω_size(T, s)) : Int64[]
+    FINAL_INDEX = !isempty(γ₁_indexes) ? ω_indexes[end] : FINAL_INDEX
+
+    o_indexes   = outlier ? collect(FINAL_INDEX + 1:FINAL_INDEX + o_size(T)) : Int64[]
+    FINAL_INDEX = outlier ? o_indexes[end] : FINAL_INDEX
+
+    exogenous_indexes = collect(FINAL_INDEX + 1:FINAL_INDEX + size(Exogenous_X, 2))
+
+    return Dict("μ₁" => μ₁_indexes, "ν₁" => ν₁_indexes, "γ₁" => γ₁_indexes, 
+                "ξ" => ξ_indexes, "ζ" => ζ_indexes, "ω" => ω_indexes, "o" => o_indexes, 
+                "Exogenous_X" => exogenous_indexes, "initial_states" => initial_states_indexes)
+end
+
+function get_variances(ϵ::Vector{Fl}, coefs::Vector{Fl}, components_indexes::Dict{String, Vector{Int64}})::Dict where Fl
+
+    variances = Dict()
+    for component in ["ξ", "ζ", "ω"]
+        !isempty(components_indexes[component]) ? variances[component] = var(coefs[components_indexes[component]]) : nothing
+    end
+    variances["ϵ"] = var(ϵ)
+    return variances
+end
+
+function forecast_model(output::Output, steps_ahead::Int64, Exogenous_Forecast::Matrix{Fl})::Vector{Float64} where Fl
+    @assert output.model_type in AVAILABLE_MODELS "Unavailable Model"
+    Exogenous_X = output.X[:, output.components["Exogenous_X"]["Indexes"]]
+    complete_matrix = create_X(output.model_type, output.T, output.s, Exogenous_X, 
+                                output.outlier, output.stabilize_ζ, steps_ahead, Exogenous_Forecast)
+
+    return complete_matrix[end-steps_ahead+1:end, :]*output.coefs
+end

src/structs.jl

@@ -0,0 +1,14 @@
+mutable struct Output 
+    model_type::String
+    X::Matrix
+    coefs::Vector
+    ϵ::Vector
+    fitted::Vector
+    components::Dict
+    residuals_variances::Dict
+    s::Int64
+    T::Int64
+    outlier::Bool
+    valid_indexes::Vector{Int64}
+    stabilize_ζ::Int64
+end

src/utils.jl

@@ -0,0 +1,17 @@
+function build_components(X::Matrix{Tl}, coefs::Vector{Float64}, components_indexes::Dict{String, Vector{Int64}})::Dict where Tl
+    components = Dict()
+    for key in keys(components_indexes)
+        components[key] = Dict()
+        components[key]["Coefs"]   = coefs[components_indexes[key]]
+        components[key]["Indexes"] = components_indexes[key]
+        components[key]["Values"]  = X[:, components_indexes[key]]*coefs[components_indexes[key]]
+    end
+    components["Exogenous_X"]["Selected"] = findall(i -> i != 0, components["Exogenous_X"]["Coefs"])
+    return components
+end
+
+function build_complete_variables(estimation_ϵ::Vector{Float64}, coefs::Vector{Float64}, X::Matrix{Tl}, valid_indexes::Vector{Int64}, T::Int64)::Tuple{Vector{Float64}, Vector{Float64}} where Tl
+    ϵ      = fill(NaN, T); ϵ[valid_indexes] = estimation_ϵ
+    fitted = X*coefs
+    return ϵ, fitted
+end