From 758d96716555bd73ed33de71dfa8e2e58a20fe98 Mon Sep 17 00:00:00 2001 From: andreramosfc Date: Thu, 16 Nov 2023 14:11:08 -0300 Subject: [PATCH] =?UTF-8?q?=F0=9F=9A=80=20first=20commit?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .github/workflows/ci.yml | 34 +++ Project.toml | 12 + README.md | Bin 568 -> 191 bytes src/StateSpaceLearning.jl | 61 +++++ src/estimation_procedure/adalasso.jl | 24 ++ src/estimation_procedure/estimation_utils.jl | 36 +++ .../information_criteria.jl | 12 + src/estimation_procedure/lasso.jl | 37 ++++ src/model_utils.jl | 130 +++++++++++ src/structs.jl | 14 ++ src/utils.jl | 17 ++ test/StateSpaceLearning.jl | 41 ++++ test/estimation_procedure/adalasso.jl | 17 ++ test/estimation_procedure/estimation_utils.jl | 42 ++++ .../information_criteria.jl | 14 ++ test/estimation_procedure/lasso.jl | 56 +++++ test/model_utils.jl | 208 ++++++++++++++++++ test/runtests.jl | 9 + test/utils.jl | 49 +++++ 19 files changed, 813 insertions(+) create mode 100644 .github/workflows/ci.yml create mode 100644 Project.toml create mode 100644 src/StateSpaceLearning.jl create mode 100644 src/estimation_procedure/adalasso.jl create mode 100644 src/estimation_procedure/estimation_utils.jl create mode 100644 src/estimation_procedure/information_criteria.jl create mode 100644 src/estimation_procedure/lasso.jl create mode 100644 src/model_utils.jl create mode 100644 src/structs.jl create mode 100644 src/utils.jl create mode 100644 test/StateSpaceLearning.jl create mode 100644 test/estimation_procedure/adalasso.jl create mode 100644 test/estimation_procedure/estimation_utils.jl create mode 100644 test/estimation_procedure/information_criteria.jl create mode 100644 test/estimation_procedure/lasso.jl create mode 100644 test/model_utils.jl create mode 100644 test/runtests.jl create mode 100644 test/utils.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..0b80406 --- /dev/null +++ b/.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()' \ No newline at end of file diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..961cc39 --- /dev/null +++ b/Project.toml @@ -0,0 +1,12 @@ +name = "StateSpaceLearning" +uuid = "971c4b7c-2c4e-4bac-8525-e842df3cde7b" +authors = ["andreramosfc "] +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" diff --git a/README.md b/README.md index 3cb66e70f57c0fbc4b811c0f62848a15245e5499..6ce17cc849f0af6546e4fbfa7896fdf946e06090 100644 GIT binary patch literal 191 zcmY#Z2rfx1NewPYOiuMlO)SdG%uDCuidKvcD$OfNEiNgJ)yOC*DJZtm*H6zZ$tX?I zOU}>L_i^+M2o4B!*2mPWpO{>dnV(mzU!Gr-otBedUaX&-saKhsqo0(RlAfwpT$Zj$ HH8V8499a8{0<$v3EI0b9&Fd$Jc@^`E44#gsnydy{8F8p?`<&l3W5C7@XJ5h zvcjX6MTbHjB|4*xilaEXvLx9Lx;an(L!Trpln!>34wVd&n)3AT%iq=luiTVtU{#(1 zKGm(NmV9tcMPsP*V9!$6k)0eGQCFpy^}TyhbAPbAg5|9CzUy<&y~9!WsnD6QEpM>x zNLc>5aCT|3epEm6v|WqZJdk0=|B#cf#cEjdiMdxZZMmvHIOBAO>E>@NF7)3fzoD;p b-U5RHJRpUntfNd&cBDXXU_etSh#`U-IR?5n diff --git a/src/StateSpaceLearning.jl b/src/StateSpaceLearning.jl new file mode 100644 index 0000000..59561de --- /dev/null +++ b/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 diff --git a/src/estimation_procedure/adalasso.jl b/src/estimation_procedure/adalasso.jl new file mode 100644 index 0000000..73fec27 --- /dev/null +++ b/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 + diff --git a/src/estimation_procedure/estimation_utils.jl b/src/estimation_procedure/estimation_utils.jl new file mode 100644 index 0000000..86cbb9d --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/estimation_procedure/information_criteria.jl b/src/estimation_procedure/information_criteria.jl new file mode 100644 index 0000000..f8be846 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/estimation_procedure/lasso.jl b/src/estimation_procedure/lasso.jl new file mode 100644 index 0000000..8863d43 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/model_utils.jl b/src/model_utils.jl new file mode 100644 index 0000000..959e3f2 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/structs.jl b/src/structs.jl new file mode 100644 index 0000000..2d6755f --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..def4369 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/StateSpaceLearning.jl b/test/StateSpaceLearning.jl new file mode 100644 index 0000000..421672c --- /dev/null +++ b/test/StateSpaceLearning.jl @@ -0,0 +1,41 @@ +@testset "Function: fit_model" begin + y1 = rand(100) + y2 = rand(100); y2[10:20] .= NaN + + output1 = StateSpaceLearning.fit_model(y1) + @test length(output1.ϵ) == 100 + @test length(output1.fitted) == 100 + @test size(output1.X, 1) == 100 + @test size(output1.X, 2) == length(output1.coefs) + + output2 = StateSpaceLearning.fit_model(y2) + @test output2.valid_indexes == setdiff(1:100, 10:20) + @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) - 1 + + @test_throws AssertionError StateSpaceLearning.fit_model(y1; s = 200) + @test_throws AssertionError StateSpaceLearning.fit_model(y1; model_type = "none") + @test_throws AssertionError StateSpaceLearning.fit_model(y1; estimation_procedure = "none") + @test_throws AssertionError StateSpaceLearning.fit_model(y1; α = -0.1) + @test_throws AssertionError StateSpaceLearning.fit_model(y1; α = 1.1) + +end + +@testset "Function: forecast" begin + y1 = rand(100) + y2 = rand(100); y2[10:20] .= NaN + + output1 = StateSpaceLearning.fit_model(y1) + @test length(StateSpaceLearning.forecast(output1, 10)) == 10 + + 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)) +end \ No newline at end of file diff --git a/test/estimation_procedure/adalasso.jl b/test/estimation_procedure/adalasso.jl new file mode 100644 index 0000000..b439398 --- /dev/null +++ b/test/estimation_procedure/adalasso.jl @@ -0,0 +1,17 @@ +@testset "Function: fit_adalasso" begin + Random.seed!(1234) + Exogenous_X = hcat(rand(10, 3), vcat(ones(3), zeros(1), ones(6))) + components_indexes = StateSpaceLearning.get_components_indexes(10, 3, Exogenous_X, true, "Basic Structural", 0) + Estimation_X = StateSpaceLearning.create_X("Basic Structural", 10, 3, Exogenous_X, true, 0) + + estimation_y = Estimation_X*rand(size(Estimation_X, 2)) + rand(10).*5 + + coefs1, ϵ1 = StateSpaceLearning.fit_adalasso(Estimation_X, estimation_y, 0.1, "aic", components_indexes, 0.1, true) + @test length(coefs1) == 43 + @test length(ϵ1) == 10 + + coefs2, ϵ2 = StateSpaceLearning.fit_adalasso(Estimation_X, estimation_y, 0.1, "aic", components_indexes, 10000.0, true) + coefs_lasso, ϵ_lasso = StateSpaceLearning.fit_lasso(Estimation_X, estimation_y, 0.1, "aic", true, components_indexes; intercept = true) + @test all(isapprox.(coefs2, coefs_lasso; atol = 1e-3)) + @test all(isapprox.(ϵ2, ϵ_lasso; atol = 1e-3)) +end \ No newline at end of file diff --git a/test/estimation_procedure/estimation_utils.jl b/test/estimation_procedure/estimation_utils.jl new file mode 100644 index 0000000..815b014 --- /dev/null +++ b/test/estimation_procedure/estimation_utils.jl @@ -0,0 +1,42 @@ +@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) + components_indexes1 = StateSpaceLearning.get_components_indexes(10, 3, Exogenous_X1, true, "Basic Structural", 0) + components_indexes2 = StateSpaceLearning.get_components_indexes(10, 3, Exogenous_X2, true, "Basic Structural", 0) + + Estimation_X1 = StateSpaceLearning.create_X("Basic Structural", 10, 3, Exogenous_X1, true, 0) + outlier_duplicate_columns1 = StateSpaceLearning.get_outlier_duplicate_columns(Estimation_X1, components_indexes1) + @test outlier_duplicate_columns1 == [32] + + Estimation_X2 = StateSpaceLearning.create_X("Basic Structural", 10, 3, Exogenous_X2, true, 0) + outlier_duplicate_columns2 = StateSpaceLearning.get_outlier_duplicate_columns(Estimation_X2, components_indexes2) + @test outlier_duplicate_columns2 == [] +end + +@testset "Function: fit_estimation_procedure" begin + Random.seed!(1234) + Exogenous_X = hcat(rand(10, 3), vcat(ones(3), zeros(1), ones(6))) + components_indexes = StateSpaceLearning.get_components_indexes(10, 3, Exogenous_X, true, "Basic Structural", 0) + Estimation_X = StateSpaceLearning.create_X("Basic Structural", 10, 3, Exogenous_X, true, 0) + estimation_y = Estimation_X*rand(size(Estimation_X, 2)) + rand(10).*5 + + coefs1, ϵ1 = StateSpaceLearning.fit_estimation_procedure("Lasso", Estimation_X, estimation_y, 0.1, "aic", components_indexes, 0.1, true) + @test length(coefs1) == 43 + @test length(ϵ1) == 10 + + coefs2, ϵ2 = StateSpaceLearning.fit_estimation_procedure("AdaLasso", Estimation_X, estimation_y, 0.1, "aic", components_indexes, 0.1, true) + @test length(coefs2) == 43 + @test length(ϵ2) == 10 +end \ No newline at end of file diff --git a/test/estimation_procedure/information_criteria.jl b/test/estimation_procedure/information_criteria.jl new file mode 100644 index 0000000..7152236 --- /dev/null +++ b/test/estimation_procedure/information_criteria.jl @@ -0,0 +1,14 @@ +@testset "Function: get_information" begin + ϵ = [1.1, 2.2, 3.3, 4.4, 5.5] + T = 5 + K = 3 + p = 10 + bic = StateSpaceLearning.get_information(T, K, ϵ; hyperparameter_selection = "bic", p = p) + aic = StateSpaceLearning.get_information(T, K, ϵ; hyperparameter_selection = "aic", p = p) + aicc = StateSpaceLearning.get_information(T, K, ϵ; hyperparameter_selection = "aicc", p = p) + EBIC = StateSpaceLearning.get_information(T, K, ϵ; hyperparameter_selection = "EBIC", p = p) + @test bic == 10.362869194716325 + @test aic == 11.534555457414024 + @test aicc == 35.53455545741402 + @test EBIC == 12.340528691778138 +end \ No newline at end of file diff --git a/test/estimation_procedure/lasso.jl b/test/estimation_procedure/lasso.jl new file mode 100644 index 0000000..36b84ea --- /dev/null +++ b/test/estimation_procedure/lasso.jl @@ -0,0 +1,56 @@ +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, α; hyperparameter_selection="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))) + components_indexes = StateSpaceLearning.get_components_indexes(10, 3, Exogenous_X, true, "Basic Structural", 0) + Estimation_X = StateSpaceLearning.create_X("Basic Structural", 10, 3, Exogenous_X, true, 0) + 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 \ No newline at end of file diff --git a/test/model_utils.jl b/test/model_utils.jl new file mode 100644 index 0000000..0953520 --- /dev/null +++ b/test/model_utils.jl @@ -0,0 +1,208 @@ +@testset "Innovation matrices" begin + @test StateSpaceLearning.ξ_size(10) == 8 + @test StateSpaceLearning.ζ_size(10, 2) == 7 + @test StateSpaceLearning.ζ_size(10, 0) == 9 + @test StateSpaceLearning.ω_size(10, 2) == 8 + @test StateSpaceLearning.o_size(10) == 10 + + X_ξ1 = StateSpaceLearning.create_ξ(5, 0) + X_ξ2 = StateSpaceLearning.create_ξ(5, 2) + + @test X_ξ1 == [0.0 0.0 0.0; + 1.0 0.0 0.0; + 1.0 1.0 0.0; + 1.0 1.0 1.0; + 1.0 1.0 1.0] + + @test X_ξ2 == [0.0 0.0 0.0; + 1.0 0.0 0.0; + 1.0 1.0 0.0; + 1.0 1.0 1.0; + 1.0 1.0 1.0; + 1.0 1.0 1.0; + 1.0 1.0 1.0] + + X_ζ1 = StateSpaceLearning.create_ζ(5, 0, 0) + X_ζ2 = StateSpaceLearning.create_ζ(5, 2, 0) + X_ζ3 = StateSpaceLearning.create_ζ(5, 2, 2) + + @test X_ζ1 == [0.0 0.0 0.0 0.0; + 2.0 0.0 0.0 0.0; + 3.0 2.0 0.0 0.0; + 4.0 3.0 2.0 0.0; + 5.0 4.0 3.0 2.0] + + @test X_ζ2 == [0.0 0.0 0.0 0.0; + 2.0 0.0 0.0 0.0; + 3.0 2.0 0.0 0.0; + 4.0 3.0 2.0 0.0; + 5.0 4.0 3.0 2.0; + 6.0 5.0 4.0 3.0; + 7.0 6.0 5.0 4.0] + + @test X_ζ3 == [0.0 0.0; + 2.0 0.0; + 3.0 2.0; + 4.0 3.0; + 5.0 4.0; + 6.0 5.0; + 7.0 6.0] + + X_ω1 = StateSpaceLearning.create_ω(5, 2, 0) + X_ω2 = StateSpaceLearning.create_ω(5, 2, 2) + + @test X_ω1 == [0.0 0.0 0.0; + -1.0 0.0 0.0; + 1.0 -1.0 0.0; + -1.0 1.0 -1.0; + 1.0 -1.0 1.0] + + @test X_ω2 == [0.0 0.0 0.0; + -1.0 0.0 0.0; + 1.0 -1.0 0.0; + -1.0 1.0 -1.0; + 1.0 -1.0 1.0; + -1.0 1.0 -1.0; + 1.0 -1.0 1.0] + + X_o1 = StateSpaceLearning.create_o_matrix(3, 0) + X_o2 = StateSpaceLearning.create_o_matrix(3, 2) + + @test X_o1 == Matrix(1.0 * I, 3, 3) + @test X_o2 == vcat(Matrix(1.0 * I, 3, 3), zeros(2, 3)) +end + + +@testset "Initial State Matrix" begin + X1 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, "Basic Structural") + X2 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, "Basic Structural") + + @test X1 == [1.0 1.0 1.0 0.0; + 1.0 2.0 0.0 1.0; + 1.0 3.0 1.0 0.0; + 1.0 4.0 0.0 1.0; + 1.0 5.0 1.0 0.0] + + @test X2 == [1.0 1.0 1.0 0.0; + 1.0 2.0 0.0 1.0; + 1.0 3.0 1.0 0.0; + 1.0 4.0 0.0 1.0; + 1.0 5.0 1.0 0.0; + 1.0 6.0 0.0 1.0; + 1.0 7.0 1.0 0.0] + + X3 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, "Local Linear Trend") + X4 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, "Local Linear Trend") + + @test X3 == [1.0 1.0; + 1.0 2.0; + 1.0 3.0; + 1.0 4.0; + 1.0 5.0] + + @test X4 == [1.0 1.0; + 1.0 2.0; + 1.0 3.0; + 1.0 4.0; + 1.0 5.0; + 1.0 6.0; + 1.0 7.0] + X5 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 0, "Local Level") + X6 = StateSpaceLearning.create_initial_states_Matrix(5, 2, 2, "Local Level") + + @test X5 == ones(5, 1) + @test X6 == ones(7, 1) +end + +@testset "Create X matrix" begin + Exogenous_X1 = rand(5, 3) + Exogenous_X2 = zeros(5, 0) + + Exogenous_forecast1 = rand(2, 3) + + param_combination = [ + Any[true, 0, 0], + Any[true, 2, 0], + Any[true, 2, 2], + Any[false, 0, 0], + Any[false, 2, 0], + Any[false, 2, 2] + ] + + size_vec1=[(5, 22), (5, 20), (7, 20), (5, 17), (5, 15), (7, 15), (5, 12), (5, 12), (7, 12), (5, 7), (5, 7), (7, 7), (5, 17), (5, 15), (7, 15), (5, 12), (5, 10), (7, 10)] + size_vec2=[(5, 19), (5, 17), (7, 17), (5, 14), (5, 12), (7, 12), (5, 9), (5, 9), (7, 9), (5, 4), (5, 4), (7, 4), (5, 14), (5, 12), (7, 12), (5, 9), (5, 7), (7, 7)] + count = 1 + for model_type in ["Basic Structural", "Local Level", "Local Linear Trend"] + for param in param_combination + if param[3] != 0 + X1 = StateSpaceLearning.create_X(model_type, 5, 2, Exogenous_X1, param[1], param[2], param[3], Exogenous_forecast1) + else + X1 = StateSpaceLearning.create_X(model_type, 5, 2, Exogenous_X1, param[1], param[2], param[3]) + end + X2 = StateSpaceLearning.create_X(model_type, 5, 2, Exogenous_X2, param[1], param[2], param[3]) + @test size(X1) == size_vec1[count] + @test size(X2) == size_vec2[count] + count += 1 + end + end + +end + +@testset "Function: get_components_indexes" begin + Exogenous_X1 = rand(10, 3) + Exogenous_X2 = zeros(10, 0) + + parameter_combination = [ + ["Basic Structural", true, Exogenous_X1], + ["Local Level", true, Exogenous_X1], + ["Local Linear Trend", true, Exogenous_X1], + ["Basic Structural", false, Exogenous_X1], + ["Basic Structural", true, Exogenous_X2], + ] + + for param in parameter_combination + + components_indexes = StateSpaceLearning.get_components_indexes(10, 3, param[3], param[2], param[1], 0) + + for key in keys(components_indexes) + if param[1] == "Basic Structural" + @test !(key in ["Exogenous_X", "o"]) ? !isempty(components_indexes[key]) : true + elseif param[1] == "Local Level" + @test !(key in ["ω", "γ₁", "ζ", "ν₁", "o", "Exogenous_X"]) ? !isempty(components_indexes[key]) : true + @test !(key in ["o", "Exogenous_X", "μ₁", "ξ", "initial_states"]) ? isempty(components_indexes[key]) : true + elseif param[1] == "Local Linear Trend" + @test !(key in ["ω", "γ₁", "o", "Exogenous_X"]) ? !isempty(components_indexes[key]) : true + @test !(key in ["μ₁", "ξ", "ζ", "ν₁", "o", "Exogenous_X", "initial_states"]) ? isempty(components_indexes[key]) : true + end + @test (param[2] && key == "o") ? !isempty(components_indexes[key]) : true + @test (!param[2] && key == "o") ? isempty(components_indexes[key]) : true + @test key == "Exogenous_X" ? length(components_indexes[key]) == size(param[3], 2) : true + end + end + +end + +@testset "Function: get_variances" begin + Exogenous_X2 = zeros(10, 0) + parameter_combination = [ + ["Basic Structural", true, Exogenous_X2, ["ξ", "ζ", "ω", "ϵ"]], + ["Local Level", true, Exogenous_X2, ["ξ", "ϵ"]], + ["Local Linear Trend", true, Exogenous_X2, ["ξ", "ζ", "ϵ"]] + ] + for param in parameter_combination + components_indexes = StateSpaceLearning.get_components_indexes(10, 3, param[3], param[2], param[1], 0) + variances = StateSpaceLearning.get_variances(rand(100), rand(39), components_indexes) + @test all([key in keys(variances) for key in param[4]]) + end +end + +@testset "Function: forecast_model" begin + Exogenous_X = rand(30, 3) + X = StateSpaceLearning.create_X("Basic Structural", 30, 2, Exogenous_X, true, 0) + coefs = rand(size(X, 2)) + components_indexes = StateSpaceLearning.get_components_indexes(30, 2, Exogenous_X, true, "Basic Structural", 0) + components = StateSpaceLearning.build_components(X, coefs, components_indexes) + + output = StateSpaceLearning.Output("Basic Structural", X, coefs, rand(30), rand(30), components, Dict(), 3, 30, true, collect(1:30), 0) + @test length(StateSpaceLearning.forecast_model(output, 5, rand(5, 3))) == 5 +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..11345cf --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,9 @@ +using StateSpaceLearning, Test, Random, LinearAlgebra, GLMNet, Statistics + +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") +include("StateSpaceLearning.jl") diff --git a/test/utils.jl b/test/utils.jl new file mode 100644 index 0000000..0dc10dd --- /dev/null +++ b/test/utils.jl @@ -0,0 +1,49 @@ +@testset "Function: build_components" begin + Exogenous_X1 = rand(10, 3) + Exogenous_X2 = zeros(10, 0) + + parameter_combination = [ + ["Basic Structural", true, Exogenous_X1, []], + ["Local Level", true, Exogenous_X1, ["ω", "γ₁", "ζ", "ν₁"]], + ["Local Linear Trend", true, Exogenous_X1, ["ω", "γ₁"]], + ["Basic Structural", false, Exogenous_X1, ["o"]], + ["Basic Structural", true, Exogenous_X2, ["Exogenous_X"]], + ] + + for param in parameter_combination + X = StateSpaceLearning.create_X(param[1], 10, 3, param[3], param[2], 0) + components_indexes = StateSpaceLearning.get_components_indexes(10, 3, param[3], param[2], param[1], 0) + coefs = rand(size(X, 2)) + components = StateSpaceLearning.build_components(X, coefs, components_indexes) + + for key in keys(components) + @test "Values" in keys(components[key]) + @test "Coefs" in keys(components[key]) + @test "Indexes" in keys(components[key]) + @test !(key in param[4]) ? !isempty(components[key]["Coefs"]) : isempty(components[key]["Coefs"]) + @test !(key in param[4]) ? !isempty(components[key]["Indexes"]) : isempty(components[key]["Indexes"]) + @test key == "Exogenous_X" ? "Selected" in keys(components[key]) : true + end + end + +end + +@testset "Function: build_complete_variables" begin + + coefs = rand(10) + X = rand(30, 10) + T = 30 + + valid_indexes1 = setdiff(collect(1:30), [11, 12]) + estimation_ϵ1 = rand(length(valid_indexes1)) + ϵ1, fitted1 = StateSpaceLearning.build_complete_variables(estimation_ϵ1, coefs, X, valid_indexes1, T) + @test all(isnan.(ϵ1[11:12])) + @test !all(isnan.(ϵ1[valid_indexes1])) + @test !all(isnan.(fitted1)) + + valid_indexes2 = collect(1:30) + estimation_ϵ2 = rand(length(valid_indexes2)) + ϵ2, fitted2 = StateSpaceLearning.build_complete_variables(estimation_ϵ2, coefs, X, valid_indexes2, T) + @test !all(isnan.(ϵ2[valid_indexes2])) + @test !all(isnan.(fitted2)) +end \ No newline at end of file