From b88b875daefbcfbfc7ee57fa0d8e4663da052765 Mon Sep 17 00:00:00 2001
From: andre_ramos <andre_ramos+1-at-806357755127>
Date: Tue, 3 Dec 2024 13:18:47 -0300
Subject: [PATCH] Add Seasonal Simulation and Multivariate Innovation
 Simulation Features

---
 .github/dependabot.yml                        |   7 +
 .gitignore                                    |   7 +-
 Project.toml                                  |   2 +-
 README.md                                     |   2 +-
 docs/src/manual.md                            |   2 +-
 paper_tests/m4_test/evaluate_model.jl         |   6 +-
 paper_tests/m4_test/m4_test.jl                |  12 +-
 paper_tests/m4_test/metrics.jl                |   2 +-
 paper_tests/m4_test/prepare_data.jl           |   4 +-
 .../simulation_test/evaluate_models.jl        |  26 +-
 paper_tests/simulation_test/simulation.jl     |   2 +-
 .../simulation_test/simulation_generator.jl   |   4 +-
 src/estimation_procedure.jl                   | 167 +++++++++---
 src/fit_forecast.jl                           | 144 ++++++----
 src/information_criteria.jl                   |  12 +-
 src/models/structural_model.jl                |  96 +++++--
 src/structs.jl                                |   2 +-
 src/utils.jl                                  | 255 ++++++++++++++++--
 test/estimation_procedure.jl                  |  10 +-
 test/fit_forecast.jl                          |  30 +++
 test/utils.jl                                 |  23 +-
 21 files changed, 633 insertions(+), 182 deletions(-)
 create mode 100644 .github/dependabot.yml

diff --git a/.github/dependabot.yml b/.github/dependabot.yml
new file mode 100644
index 0000000..ff6499d
--- /dev/null
+++ b/.github/dependabot.yml
@@ -0,0 +1,7 @@
+# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "weekly"
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
index 25394c9..e15812f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -5,6 +5,11 @@ local_test/
 paper_tests/m4_test/metrics_results/
 paper_tests/m4_test/results_ARIMA/
 paper_tests/m4_test/results_SSL/
+paper_tests/m4_test/init_SSL/
+
 paper_tests/simulation_test/results_simulation_raw/
+paper_tests/mv_probabilistic_forecast/
 
-docs/build
\ No newline at end of file
+docs/build
+longhorizon/
+test.jl
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index d3864c0..992ffe7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "StateSpaceLearning"
 uuid = "971c4b7c-2c4e-4bac-8525-e842df3cde7b"
 authors = ["andreramosfc <andreramosfdc@gmail.com>"]
-version = "1.0.2"
+version = "1.1.0"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
diff --git a/README.md b/README.md
index 9f86d92..06b902e 100644
--- a/README.md
+++ b/README.md
@@ -153,7 +153,7 @@ using Plots
 airp = CSV.File(StateSpaceLearning.AIR_PASSENGERS) |> DataFrame
 log_air_passengers = log.(airp.passengers)
 
-airpassengers = Float64.(airp.passengers)
+airpassengers = AbstractFloat.(airp.passengers)
 log_air_passengers[60:72] .= NaN
 
 model = StructuralModel(log_air_passengers)
diff --git a/docs/src/manual.md b/docs/src/manual.md
index 4fd8d60..4855bbc 100644
--- a/docs/src/manual.md
+++ b/docs/src/manual.md
@@ -149,7 +149,7 @@ using Plots
 airp = CSV.File(StateSpaceLearning.AIR_PASSENGERS) |> DataFrame
 log_air_passengers = log.(airp.passengers)
 
-airpassengers = Float64.(airp.passengers)
+airpassengers = AbstractFloat.(airp.passengers)
 log_air_passengers[60:72] .= NaN
 
 model = StructuralModel(log_air_passengers)
diff --git a/paper_tests/m4_test/evaluate_model.jl b/paper_tests/m4_test/evaluate_model.jl
index faace47..1f1055a 100644
--- a/paper_tests/m4_test/evaluate_model.jl
+++ b/paper_tests/m4_test/evaluate_model.jl
@@ -3,9 +3,9 @@ function evaluate_SSL(
     results_df::DataFrame,
     input::Dict,
     outlier::Bool,
-    α::Float64,
-    H::Int64,
-    sample_size::Int64,
+    α::AbstractFloat,
+    H::Int,
+    sample_size::Int,
     information_criteria::String,
 )
     normalized_y = input["normalized_train"]
diff --git a/paper_tests/m4_test/m4_test.jl b/paper_tests/m4_test/m4_test.jl
index c7a580f..890ed14 100644
--- a/paper_tests/m4_test/m4_test.jl
+++ b/paper_tests/m4_test/m4_test.jl
@@ -33,9 +33,9 @@ function run_config(
     results_table::DataFrame,
     outlier::Bool,
     information_criteria::String,
-    α::Float64,
+    α::AbstractFloat,
     save_init::Bool,
-    sample_size::Int64,
+    sample_size::Int,
 )
     NAIVE_sMAPE = 14.427 #M4 Paper
     NAIVE_MASE = 1.063   #M4 Paper
@@ -98,9 +98,9 @@ end
 # Main script
 function main()
     results_table = DataFrame()
-    for outlier in [true, false]
-        for information_criteria in ["aic", "bic"]
-            for α in vcat([0.0], collect(0.1:0.2:0.9), [1.0])
+    for outlier in [true]
+        for information_criteria in ["aic"]
+            for α in [0.1]
                 @info "Running SSL with outlier = $outlier, information_criteria = $information_criteria, α = $α"
                 results_table = run_config(
                     results_table, outlier, information_criteria, α, false, 60
@@ -140,4 +140,4 @@ create_dirs()
 
 main()
 
-run_config(DataFrame(), false, "aic", 0.1, true, 2794)#max sample size
+#run_config(DataFrame(), false, "aic", 0.1, true, 2794)#max sample size
diff --git a/paper_tests/m4_test/metrics.jl b/paper_tests/m4_test/metrics.jl
index ed5f8ae..c7714c6 100644
--- a/paper_tests/m4_test/metrics.jl
+++ b/paper_tests/m4_test/metrics.jl
@@ -5,7 +5,7 @@ function sMAPE(y_test::Vector, prediction::Vector)
     return (200 / H) * sum(abs(y_test[i] - prediction[i]) / (denominator[i]) for i in 1:H)
 end
 
-function MASE(y_train::Vector, y_test::Vector, prediction::Vector; m::Int64=12)
+function MASE(y_train::Vector, y_test::Vector, prediction::Vector; m::Int=12)
     T = length(y_train)
     H = length(y_test)
 
diff --git a/paper_tests/m4_test/prepare_data.jl b/paper_tests/m4_test/prepare_data.jl
index 5a81cc0..9d69a57 100644
--- a/paper_tests/m4_test/prepare_data.jl
+++ b/paper_tests/m4_test/prepare_data.jl
@@ -1,8 +1,8 @@
-function normalize(y::Vector, max_y::Float64, min_y::Float64)
+function normalize(y::Vector, max_y::AbstractFloat, min_y::AbstractFloat)
     return (y .- min_y) ./ (max_y - min_y)
 end
 
-function de_normalize(y::Vector, max_y::Float64, min_y::Float64)
+function de_normalize(y::Vector, max_y::AbstractFloat, min_y::AbstractFloat)
     return (y .* (max_y - min_y)) .+ min_y
 end
 
diff --git a/paper_tests/simulation_test/evaluate_models.jl b/paper_tests/simulation_test/evaluate_models.jl
index 4700204..8843c7c 100644
--- a/paper_tests/simulation_test/evaluate_models.jl
+++ b/paper_tests/simulation_test/evaluate_models.jl
@@ -8,13 +8,13 @@ function get_confusion_matrix(selected, true_features, false_features)
 end
 
 function get_SSL_results(
-    y_train::Vector{Float64},
+    y_train::Vector{Fl},
     true_features::Vector{Int64},
     false_features::Vector{Int64},
-    X_train::Matrix{Float64},
+    X_train::Matrix{Fl},
     inf_criteria::String,
-    true_β::Vector{Float64},
-)
+    true_β::Vector{Fl},
+) where Fl <: AbstractFloat
     series_result = nothing
 
     model = StateSpaceLearning.StructuralModel(
@@ -72,13 +72,13 @@ function get_SSL_results(
 end
 
 function get_SS_res_results(
-    y_train::Vector{Float64},
+    y_train::Vector{Fl},
     true_features::Vector{Int64},
     false_features::Vector{Int64},
-    X_train::Matrix{Float64},
+    X_train::Matrix{Fl},
     inf_criteria::String,
-    true_β::Vector{Float64},
-)
+    true_β::Vector{Fl},
+) where Fl <: AbstractFloat
     py"""
     import math
     import statsmodels.api as sm
@@ -141,7 +141,7 @@ function get_SS_res_results(
     return series_result, converged
 end
 
-function get_exogenous_ss_inf_criteria(y_train::Vector{Float64}, X_train::Matrix{Float64})
+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
@@ -160,13 +160,13 @@ function get_exogenous_ss_inf_criteria(y_train::Vector{Float64}, X_train::Matrix
 end
 
 function get_forward_ss(
-    y_train::Vector{Float64},
+    y_train::Vector{Fl},
     true_features::Vector{Int64},
     false_features::Vector{Int64},
-    X_train::Matrix{Float64},
+    X_train::Matrix{Fl},
     inf_criteria::String,
-    true_β::Vector{Float64},
-)
+    true_β::Vector{Fl},
+) where Fl <: AbstractFloat
     best_inf_crit = Inf
     current_inf_crit = 0
     coefs = nothing
diff --git a/paper_tests/simulation_test/simulation.jl b/paper_tests/simulation_test/simulation.jl
index b55ebc1..de289dd 100644
--- a/paper_tests/simulation_test/simulation.jl
+++ b/paper_tests/simulation_test/simulation.jl
@@ -327,7 +327,7 @@ CSV.write("paper_tests/simulation_test/results_metrics/metrics_confusion_matrix.
 
 df_mse_bias = DataFrame()
 
-function convert_to_sci_notation(num::Float64)
+function convert_to_sci_notation(num::AbstractFloat)
     # Get the exponent part of the number in scientific notation
     exp_part = floor(log10(abs(num)))
 
diff --git a/paper_tests/simulation_test/simulation_generator.jl b/paper_tests/simulation_test/simulation_generator.jl
index 3b04bc8..acb4a34 100644
--- a/paper_tests/simulation_test/simulation_generator.jl
+++ b/paper_tests/simulation_test/simulation_generator.jl
@@ -23,7 +23,7 @@ function assert_invertibility(q::Vector{Fl}) where {Fl}
     return all(abs.(roots_of_inverse_polinomial(poly)) .< 1)
 end
 
-function generate_sarima_exog(T::Int64, M::Int64)
+function generate_sarima_exog(T::Int, M::Int)
     X = zeros(T, M)
     s = 12
     for j in 1:M
@@ -87,7 +87,7 @@ function generate_sarima_exog(T::Int64, M::Int64)
     return X
 end
 
-function generate_subset(T::Int64, M::Int64, K::Int64)
+function generate_subset(T::Int, M::Int, K::Int)
     s = 12
 
     μ1 = 1.0
diff --git a/src/estimation_procedure.jl b/src/estimation_procedure.jl
index 367215a..a1a1bc6 100644
--- a/src/estimation_procedure.jl
+++ b/src/estimation_procedure.jl
@@ -24,12 +24,12 @@ function get_dummy_indexes(Exogenous_X::Matrix{Fl}) where {Fl}
 end
 
 """
-    get_outlier_duplicate_columns(Estimation_X::Matrix{Tl}, components_indexes::Dict{String, Vector{Int}}) where{Tl}
+    get_outlier_duplicate_columns(Estimation_X::Matrix{Fl}, components_indexes::Dict{String, Vector{Int}}) where{Fl}
 
     Identifies and returns the indexes of outlier columns that are duplicates of dummy variables in the exogenous matrix.
 
     # Arguments
-    - `Estimation_X::Matrix{Tl}`: Matrix used for estimation.
+    - `Estimation_X::Matrix{Fl}`: Matrix used for estimation.
     - `components_indexes::Dict{String, Vector{Int}}`: Dictionary containing indexes for different components.
 
     # Returns
@@ -37,8 +37,8 @@ end
 
 """
 function get_outlier_duplicate_columns(
-    Estimation_X::Matrix{Tl}, components_indexes::Dict{String,Vector{Int}}
-) where {Tl}
+    Estimation_X::Matrix{Fl}, components_indexes::Dict{String,Vector{Int}}
+) where {Fl}
     if !haskey(components_indexes, "o")
         return []
     else
@@ -52,8 +52,13 @@ function get_outlier_duplicate_columns(
 end
 
 """
-    get_path_information_criteria(model::GLMNetPath, Lasso_X::Matrix{Tl}, Lasso_y::Vector{Fl},
-        information_criteria::String; intercept::Bool = true)::Tuple{Vector{Float64}, Vector{Float64}} where {Tl, Fl}
+    get_path_information_criteria(
+    model::GLMNetPath,
+    Lasso_X::Matrix{Tl},
+    Lasso_y::Vector{Fl},
+    information_criteria::String;
+    intercept::Bool=true,
+)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}} where {Fl <: AbstractFloat, Tl <: AbstractFloat}
 
     Calculates the information criteria along the regularization path of a GLMNet model and returns coefficients and residuals of the best model based on the selected information criteria.
 
@@ -65,7 +70,7 @@ end
     - `intercept::Bool`: Flag for intercept inclusion in the model (default: true).
 
     # Returns
-    - `Tuple{Vector{Float64}, Vector{Float64}}`: Tuple containing coefficients and residuals of the best model.
+    - `Tuple{Vector{AbstractFloat}, Vector{AbstractFloat}}`: Tuple containing coefficients and residuals of the best model.
 
 """
 function get_path_information_criteria(
@@ -74,12 +79,12 @@ function get_path_information_criteria(
     Lasso_y::Vector{Fl},
     information_criteria::String;
     intercept::Bool=true,
-)::Tuple{Vector{Float64},Vector{Float64}} where {Tl,Fl}
+)::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)'
 
-    method_vec = Vector{Float64}(undef, path_size)
+    method_vec = Vector{AbstractFloat}(undef, path_size)
     for i in 1:path_size
         fit = Lasso_X * model.betas[:, i] .+ model.a0[i]
         ε = Lasso_y - fit
@@ -101,33 +106,37 @@ function get_path_information_criteria(
 end
 
 """
-    fit_glmnet(Lasso_X::Matrix{Tl}, Lasso_y::Vector{Fl}, α::Float64;
-               information_criteria::String = "aic",
-               penalty_factor::Vector{Float64}=ones(size(Lasso_X,2) - 1),
-               intercept::Bool = intercept)::Tuple{Vector{Float64}, Vector{Float64}} where {Tl, Fl}
+    fit_glmnet(
+    Lasso_X::Matrix{Tl},
+    Lasso_y::Vector{Fl},
+    α::AbstractFloat;
+    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}
 
     Fits a GLMNet model to the provided data and returns coefficients and residuals based on selected criteria.
 
     # Arguments
     - `Lasso_X::Matrix{Tl}`: Matrix of predictors for estimation.
     - `Lasso_y::Vector{Fl}`: Vector of response values for estimation.
-    - `α::Float64`: Elastic net control factor between ridge (α=0) and lasso (α=1) (default: 0.1).
+    - `α::AbstractFloat`: Elastic net control factor between ridge (α=0) and lasso (α=1) (default: 0.1).
     - `information_criteria::String`: Information Criteria method for hyperparameter selection (default: aic).
-    - `penalty_factor::Vector{Float64}`: Penalty factors for each predictor (default: ones(size(Lasso_X, 2) - 1)).
+    - `penalty_factor::Vector{Pl}`: Penalty factors for each predictor (default: ones(size(Lasso_X, 2) - 1)).
     - `intercept::Bool`: Flag for intercept inclusion in the model (default: true).
 
     # Returns
-    - `Tuple{Vector{Float64}, Vector{Float64}}`: Tuple containing coefficients and residuals of the best model.
+    - `Tuple{Vector{AbstractFloat}, Vector{AbstractFloat}}`: Tuple containing coefficients and residuals of the best model.
 
 """
 function fit_glmnet(
     Lasso_X::Matrix{Tl},
     Lasso_y::Vector{Fl},
-    α::Float64;
+    α::AbstractFloat;
     information_criteria::String="aic",
-    penalty_factor::Vector{Float64}=ones(size(Lasso_X, 2) - 1),
+    penalty_factor::Vector{Pl}=ones(size(Lasso_X, 2) - 1),
     intercept::Bool=intercept,
-)::Tuple{Vector{Float64},Vector{Float64}} where {Tl,Fl}
+)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}} where {Fl <: AbstractFloat, Tl <: AbstractFloat, Pl <: AbstractFloat}
     model = glmnet(
         Lasso_X,
         Lasso_y;
@@ -143,36 +152,43 @@ function fit_glmnet(
 end
 
 """
-    fit_lasso(Estimation_X::Matrix{Tl}, estimation_y::Vector{Fl}, α::Float64, information_criteria::String,
-              penalize_exogenous::Bool, components_indexes::Dict{String, Vector{Int}}, penalty_factor::Vector{Float64};
-              rm_average::Bool = false)::Tuple{Vector{Float64}, Vector{Float64}} where {Tl, Fl}
+    fit_lasso(
+    Estimation_X::Matrix{Tl},
+    estimation_y::Vector{Fl},
+    α::AbstractFloat,
+    information_criteria::String,
+    penalize_exogenous::Bool,
+    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}
 
     Fits a Lasso regression model to the provided data and returns coefficients and residuals based on selected criteria.
 
     # Arguments
-    - `Estimation_X::Matrix{Tl}`: Matrix of predictors for estimation.
+    - `Estimation_X::Matrix{Fl}`: Matrix of predictors for estimation.
     - `estimation_y::Vector{Fl}`: Vector of response values for estimation.
-    - `α::Float64`: Elastic net control factor between ridge (α=0) and lasso (α=1) (default: 0.1).
+    - `α::AbstractFloat`: Elastic net control factor between ridge (α=0) and lasso (α=1) (default: 0.1).
     - `information_criteria::String`: Information Criteria method for hyperparameter selection (default: aic).
     - `penalize_exogenous::Bool`: Flag for selecting exogenous variables. When false the penalty factor for these variables will be set to 0.
     - `components_indexes::Dict{String, Vector{Int}}`: Dictionary containing indexes for different components.
-    - `penalty_factor::Vector{Float64}`: Penalty factors for each predictor.
+    - `penalty_factor::Vector{Fl}`: Penalty factors for each predictor.
     - `rm_average::Bool`: Flag to consider if the intercept will be calculated is the average of the time series (default: false).
 
     # Returns
-    - `Tuple{Vector{Float64}, Vector{Float64}}`: Tuple containing coefficients and residuals of the fitted Lasso model.
+    - `Tuple{Vector{AbstractFloat}, Vector{AbstractFloat}}`: Tuple containing coefficients and residuals of the fitted Lasso model.
 
 """
 function fit_lasso(
     Estimation_X::Matrix{Tl},
     estimation_y::Vector{Fl},
-    α::Float64,
+    α::AbstractFloat,
     information_criteria::String,
     penalize_exogenous::Bool,
     components_indexes::Dict{String,Vector{Int}},
-    penalty_factor::Vector{Float64};
+    penalty_factor::Vector{Pl};
     rm_average::Bool=false,
-)::Tuple{Vector{Float64},Vector{Float64}} where {Tl,Fl}
+)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}} where {Fl <: AbstractFloat, Tl <: AbstractFloat, Pl <: AbstractFloat}
     outlier_duplicate_columns = get_outlier_duplicate_columns(
         Estimation_X, components_indexes
     )
@@ -226,39 +242,43 @@ function fit_lasso(
 end
 
 """
-    estimation_procedure(Estimation_X::Matrix{Tl}, estimation_y::Vector{Fl}, 
-                        components_indexes::Dict{String, Vector{Int}}, α::Float64, 
-                        information_criteria::String, 
-                        ϵ::Float64, 
-                        penalize_exogenous::Bool, 
-                        penalize_initial_states::Bool)::Tuple{Vector{Float64}, Vector{Float64}} where {Tl, Fl}
+    estimation_procedure(
+    Estimation_X::Matrix{Tl},
+    estimation_y::Vector{Fl},
+    components_indexes::Dict{String,Vector{Int}},
+    α::AbstractFloat,
+    information_criteria::String,
+    ϵ::AbstractFloat,
+    penalize_exogenous::Bool,
+    penalize_initial_states::Bool,
+)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}} where {Fl <: AbstractFloat, Tl <: AbstractFloat}
 
     Fits an Adaptive Lasso (AdaLasso) regression model to the provided data and returns coefficients and residuals.
 
     # Arguments
-    - `Estimation_X::Matrix{Tl}`: Matrix of predictors for estimation.
+    - `Estimation_X::Matrix{Fl}`: Matrix of predictors for estimation.
     - `estimation_y::Vector{Fl}`: Vector of response values for estimation.
     - `components_indexes::Dict{String, Vector{Int}}`: Dictionary containing indexes for different components.
-    - `α::Float64`: Elastic net control factor between ridge (α=0) and lasso (α=1) (default: 0.1).
+    - `α::AbstractFloat`: Elastic net control factor between ridge (α=0) and lasso (α=1) (default: 0.1).
     - `information_criteria::String`: Information Criteria method for hyperparameter selection (default: aic).
-    - `ϵ::Float64`: Non negative value to handle 0 coefs on the first lasso step (default: 0.05).
+    - `ϵ::AbstractFloat`: Non negative value to handle 0 coefs on the first lasso step (default: 0.05).
     - `penalize_exogenous::Bool`: Flag for selecting exogenous variables. When false the penalty factor for these variables will be set to 0.
     - `penalize_initial_states::Bool`: Flag for selecting initial states. When false the penalty factor for these variables will be set to 0.
 
     # Returns
-    - `Tuple{Vector{Float64}, Vector{Float64}}`: Tuple containing coefficients and residuals of the fitted AdaLasso model.
+    - `Tuple{Vector{AbstractFloat}, Vector{AbstractFloat}}`: Tuple containing coefficients and residuals of the fitted AdaLasso model.
 
 """
 function estimation_procedure(
     Estimation_X::Matrix{Tl},
     estimation_y::Vector{Fl},
     components_indexes::Dict{String,Vector{Int}},
-    α::Float64,
+    α::AbstractFloat,
     information_criteria::String,
-    ϵ::Float64,
+    ϵ::AbstractFloat,
     penalize_exogenous::Bool,
     penalize_initial_states::Bool,
-)::Tuple{Vector{Float64},Vector{Float64}} where {Tl,Fl}
+)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}} where {Fl <: AbstractFloat, Tl <: AbstractFloat}
     @assert 0 <= α <= 1 "α must be in [0, 1]"
     @assert ϵ > 0 "ϵ must be positive"
 
@@ -279,7 +299,7 @@ function estimation_procedure(
         )
     else
         penalty_factor = ones(size(Estimation_X, 2))
-        penalty_factor[components_indexes["initial_states"][2:end]] .= 0
+        penalty_factor[components_indexes["initial_states"][1:end]] .= 0
         coefs, _ = fit_lasso(
             Estimation_X,
             estimation_y,
@@ -323,7 +343,7 @@ function estimation_procedure(
         end
     else
         if !penalize_initial_states
-            ts_penalty_factor[components_indexes["initial_states"][2:end]] .= 0
+            ts_penalty_factor[components_indexes["initial_states"][1:end]] .= 0
         else
             nothing
         end
@@ -340,3 +360,62 @@ function estimation_procedure(
         rm_average=false,
     )
 end
+
+"""
+    estimation_procedure(
+    Estimation_X::Matrix{Tl},
+    estimation_y::Matrix{Fl},
+    components_indexes::Dict{String,Vector{Int}},
+    α::AbstractFloat,
+    information_criteria::String,
+    ϵ::AbstractFloat,
+    penalize_exogenous::Bool,
+    penalize_initial_states::Bool,
+)::Tuple{Vector{Vector{AbstractFloat}},Vector{Vector{AbstractFloat}}} where {Fl <: AbstractFloat, Tl <: AbstractFloat}
+
+    Fits an Adaptive Lasso (AdaLasso) regression model to the provided data and returns coefficients and residuals.
+
+    # Arguments
+    - `Estimation_X::Matrix{Fl}`: Matrix of predictors for estimation.
+    - `estimation_y::Matrix{Fl}`: Matrix of response values for estimation.
+    - `components_indexes::Dict{String, Vector{Int}}`: Dictionary containing indexes for different components.
+    - `α::AbstractFloat`: Elastic net control factor between ridge (α=0) and lasso (α=1) (default: 0.1).
+    - `information_criteria::String`: Information Criteria method for hyperparameter selection (default: aic).
+    - `ϵ::AbstractFloat`: Non negative value to handle 0 coefs on the first lasso step (default: 0.05).
+    - `penalize_exogenous::Bool`: Flag for selecting exogenous variables. When false the penalty factor for these variables will be set to 0.
+    - `penalize_initial_states::Bool`: Flag for selecting initial states. When false the penalty factor for these variables will be set to 0.
+
+    # Returns
+    - `Tuple{Vector{AbstractFloat}, Vector{AbstractFloat}}`: Tuple containing coefficients and residuals of the fitted AdaLasso model.
+
+"""
+function estimation_procedure(
+    Estimation_X::Matrix{Tl},
+    estimation_y::Matrix{Fl},
+    components_indexes::Dict{String,Vector{Int}},
+    α::AbstractFloat,
+    information_criteria::String,
+    ϵ::AbstractFloat,
+    penalize_exogenous::Bool,
+    penalize_initial_states::Bool,
+)::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)
+
+    for i in 1:N_series
+        coef_i, ε_i = estimation_procedure(
+            Estimation_X,
+            estimation_y[:, i],
+            components_indexes,
+            α,
+            information_criteria,
+            ϵ,
+            penalize_exogenous,
+            penalize_initial_states,
+        )
+        push!(coefs_vec, coef_i)
+        push!(ε_vec, ε_i)
+    end
+    return coefs_vec, ε_vec
+end
\ No newline at end of file
diff --git a/src/fit_forecast.jl b/src/fit_forecast.jl
index 47e2ea5..9900150 100644
--- a/src/fit_forecast.jl
+++ b/src/fit_forecast.jl
@@ -1,9 +1,9 @@
 
 """
 function fit!(model::StateSpaceLearningModel,
-                    α::Float64 = 0.1, 
+                    α::AbstractFloat = 0.1, 
                     information_criteria::String = "aic", 
-                    ϵ::Float64 = 0.05, 
+                    ϵ::AbstractFloat = 0.05, 
                     penalize_exogenous::Bool = true, 
                     penalize_initial_states::Bool = true,
                     )
@@ -12,17 +12,17 @@ function fit!(model::StateSpaceLearningModel,
 
     # Arguments
     model::StateSpaceLearningModel: Model to be fitted.
-    α::Float64: Elastic net mixing parameter (default: 0.1).
+    α::AbstractFloat: Elastic net mixing parameter (default: 0.1).
     information_criteria::String: Method for hyperparameter selection (default: "aic").
-    ϵ::Float64: Non negative value to handle 0 coefs on the first lasso step (default: 0.05).
+    ϵ::AbstractFloat: Non negative value to handle 0 coefs on the first lasso step (default: 0.05).
     penalize_exogenous::Bool: If true, penalize exogenous variables (default: true).
     penalize_initial_states::Bool: If true, penalize initial states (default: true).
 """
 function fit!(
     model::StateSpaceLearningModel;
-    α::Float64=0.1,
+    α::AbstractFloat=0.1,
     information_criteria::String="aic",
-    ϵ::Float64=0.05,
+    ϵ::AbstractFloat=0.05,
     penalize_exogenous::Bool=true,
     penalize_initial_states::Bool=true,
 )
@@ -50,16 +50,25 @@ function fit!(
 
     residuals_variances = get_variances(model, estimation_ε, coefs, components_indexes)
 
+    T = typeof(model.y) <:Vector ? length(model.y) : size(model.y, 1)
+
     ε, fitted = get_fit_and_residuals(
-        estimation_ε, coefs, model.X, valid_indexes, length(model.y)
+        estimation_ε, coefs, model.X, valid_indexes, T
     )
 
-    output = Output(coefs, ε, fitted, residuals_variances, valid_indexes, components)
+    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]))
+        end
+    end
     return model.output = output
 end
 
 """
-    forecast(model::StateSpaceLearningModel, steps_ahead::Int; Exogenous_Forecast::Union{Matrix{Fl}, Missing}=missing)::Vector{Float64} where Fl
+    forecast(model::StateSpaceLearningModel, steps_ahead::Int; Exogenous_Forecast::Union{Matrix{Fl}, Missing}=missing)::Vector{AbstractFloat} where Fl
 
     Returns the forecast for a given number of steps ahead using the provided StateSpaceLearning output and exogenous forecast data.
 
@@ -69,19 +78,21 @@ end
     - `Exogenous_Forecast::Matrix{Fl}`: Exogenous variables forecast (default: zeros(steps_ahead, 0))
 
     # Returns
-    - `Vector{Float64}`: Vector containing forecasted values.
+    - `Union{Matrix{AbstractFloat}, Vector{AbstractFloat}}`: Matrix or vector of matrices containing forecasted values.
 
 """
 function forecast(
     model::StateSpaceLearningModel,
     steps_ahead::Int;
     Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, 0),
-)::Vector{Float64} where {Fl}
-    @assert length(model.output.components["Exogenous_X"]["Indexes"]) ==
+)::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"
     @assert size(Exogenous_Forecast, 1) == steps_ahead "Exogenous_Forecast must have the same number of rows as steps_ahead"
 
-    Exogenous_X = model.X[:, model.output.components["Exogenous_X"]["Indexes"]]
+    Exogenous_X = model.X[:, exog_idx]
     complete_matrix = create_X(
         model.level,
         model.stochastic_level,
@@ -97,12 +108,20 @@ function forecast(
         Exogenous_Forecast,
     )
 
-    return complete_matrix[(end - steps_ahead + 1):end, :] * model.output.coefs
+    if typeof(model.output) == Output
+        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
+        end
+        return AbstractFloat.(prediction)
+    end
 end
 
 """
 simulate(model::StateSpaceLearningModel, steps_ahead::Int, N_scenarios::Int;
-                                 Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, 0))::Matrix{Float64} where Fl
+                                 Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, 0))::Matrix{AbstractFloat} where Fl
 
     Generate simulations for a given number of steps ahead using the provided StateSpaceLearning output and exogenous forecast data.
 
@@ -113,60 +132,91 @@ simulate(model::StateSpaceLearningModel, steps_ahead::Int, N_scenarios::Int;
     - `Exogenous_Forecast::Matrix{Fl}`: Exogenous variables forecast (default: zeros(steps_ahead, 0))
 
     # Returns
-    - `Matrix{Float64}`: Matrix containing simulated values.
+    - `Union{Vector{Matrix{AbstractFloat}}, Matrix{AbstractFloat}}`: Matrix or vector of matrices containing simulated values.
 """
 function simulate(
     model::StateSpaceLearningModel,
     steps_ahead::Int,
     N_scenarios::Int;
     Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, 0),
-)::Matrix{Float64} where {Fl}
-    prediction = forecast(model, steps_ahead; Exogenous_Forecast=Exogenous_Forecast)
+    seasonal_innovation_simulation::Int=0,
+)::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)
+
+    is_univariate = typeof(model.output) == StateSpaceLearning.Output
 
     simulation_X = zeros(steps_ahead, 0)
-    components_matrix = zeros(length(model.output.valid_indexes), 0)
+    valid_indexes = is_univariate ? model.output.valid_indexes : model.output[1].valid_indexes
+    components_matrix = zeros(length(valid_indexes), 0)
     N_components = 1
 
-    model_innovations = get_model_innovations(model)
+    model_innovations = StateSpaceLearning.get_model_innovations(model)
     for innovation in model_innovations
-        if innovation in keys(model.output.residuals_variances)
-            simulation_X = hcat(
-                simulation_X,
-                get_innovation_simulation_X(model, innovation, steps_ahead)[
-                    (end - steps_ahead):(end - 1), (end - steps_ahead + 1):end
-                ],
-            )
-            comp = fill_innovation_coefs(model, innovation)
-            components_matrix = hcat(components_matrix, comp[model.output.valid_indexes])
-            N_components += 1
-        end
+        simulation_X = hcat(
+            simulation_X,
+            StateSpaceLearning.get_innovation_simulation_X(model, innovation, steps_ahead)[
+                (end - steps_ahead):(end - 1), (end - steps_ahead + 1):end
+            ],
+        )
+        comp = StateSpaceLearning.fill_innovation_coefs(model, innovation, valid_indexes)
+        components_matrix = hcat(components_matrix, comp)
+        N_components += 1
     end
 
-    components_matrix = hcat(components_matrix, model.output.ε[model.output.valid_indexes])
-    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
-
-    ∑ = cov(components_matrix)
-    MV_dist = MvNormal(zeros(N_components), ∑)
-    o_noises = if model.outlier
-        rand(Normal(0, std(model.output.components["o"]["Coefs"])), steps_ahead, N_scenarios)
+    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"
     else
-        zeros(steps_ahead, N_scenarios)
+        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"
     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
 
-    simulation = hcat([prediction for _ in 1:N_scenarios]...)
-    for s in 1:N_scenarios
-        sim_coefs = ones(size(simulation_X, 2)) .* NaN
+    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)]
 
+    if seasonal_innovation_simulation == 0
+        ∑ = cov(components_matrix)
         for i in 1:steps_ahead
-            rand_inovs = rand(MV_dist)
+            MV_dist_vec[i] = is_univariate ? MvNormal(zeros(N_components), ∑) : MvNormal(zeros(N_mv_components), ∑)
+        end
 
-            for comp in eachindex(rand_inovs)
-                sim_coefs[i + (comp - 1) * steps_ahead] = rand_inovs[comp]
+        if model.outlier
+            if is_univariate
+                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)]
             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), ∑)
+            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
+            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
+                end
+            end
+        end 
 
-        simulation[:, s] += (simulation_X * sim_coefs + o_noises[:, s])
+    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))
+        simulation = Vector{Matrix{<:AbstractFloat}}(simulation)
     end
 
     return simulation
diff --git a/src/information_criteria.jl b/src/information_criteria.jl
index e205778..fc0103e 100644
--- a/src/information_criteria.jl
+++ b/src/information_criteria.jl
@@ -1,22 +1,22 @@
 """
-    get_information(T::Int, K::Int, ε::Vector{Float64};
-                    information_criteria::String = "bic")::Float64
+    get_information(T::Int, K::Int, ε::Vector{Fl};
+                    information_criteria::String = "bic")::AbstractFloat where Fl <: AbstractFloat
 
     Calculates information criterion value based on the provided parameters and residuals.
 
     # Arguments
     - `T::Int`: Number of observations.
     - `K::Int`: Number of selected predictors.
-    - `ε::Vector{Float64}`: Vector of residuals.
+    - `ε::Vector{Fl}`: Vector of residuals.
     - `information_criteria::String`: Method for hyperparameter selection (default: "aic").
 
     # Returns
-    - `Float64`: Information criterion value.
+    - `AbstractFloat`: Information criterion value.
 
 """
 function get_information(
-    T::Int, K::Int, ε::Vector{Float64}; information_criteria::String="aic"
-)::Float64
+    T::Int, K::Int, ε::Vector{Fl}; information_criteria::String="aic"
+)::AbstractFloat where Fl <: AbstractFloat
     if information_criteria == "bic"
         return T * log(var(ε)) + K * log(T)
     elseif information_criteria == "aic"
diff --git a/src/models/structural_model.jl b/src/models/structural_model.jl
index 981c9d6..392aa0a 100644
--- a/src/models/structural_model.jl
+++ b/src/models/structural_model.jl
@@ -1,5 +1,5 @@
 mutable struct StructuralModel <: StateSpaceLearningModel
-    y::Vector
+    y::Union{Vector, Matrix}
     X::Matrix
     level::Bool
     stochastic_level::Bool
@@ -11,10 +11,10 @@ mutable struct StructuralModel <: StateSpaceLearningModel
     outlier::Bool
     ζ_ω_threshold::Int
     n_exogenous::Int
-    output::Union{Output,Nothing}
+    output::Union{Vector{Output}, Output, Nothing}
 
     function StructuralModel(
-        y::Vector{Fl};
+        y::Union{Vector{Fl}, Matrix{Fl}};
         level::Bool=true,
         stochastic_level::Bool=true,
         trend::Bool=true,
@@ -24,12 +24,15 @@ mutable struct StructuralModel <: StateSpaceLearningModel
         freq_seasonal::Int=12,
         outlier::Bool=true,
         ζ_ω_threshold::Int=12,
-        Exogenous_X::Matrix{Fl}=zeros(length(y), 0),
+        Exogenous_X::Matrix{Fl}=typeof(y) <:Vector ? zeros(length(y), 0) : zeros(size(y, 1), 0),
     ) where {Fl}
         n_exogenous = size(Exogenous_X, 2)
         @assert !has_intercept(Exogenous_X) "Exogenous matrix must not have an intercept column"
-        @assert seasonal ? length(y) > freq_seasonal : true "Time series must be longer than the seasonal period"
-
+        if typeof(y) <:Vector
+            @assert seasonal ? length(y) > freq_seasonal : true "Time series must be longer than the seasonal period"
+        else
+            @assert seasonal ? size(y, 1) > freq_seasonal : true "Time series must be longer than the seasonal period"
+        end
         X = create_X(
             level,
             stochastic_level,
@@ -119,7 +122,7 @@ end
 
 """
 function create_ξ(T::Int, steps_ahead::Int)::Matrix
-    ξ_matrix = Matrix{Float64}(undef, T + steps_ahead, T - 1)
+    ξ_matrix = Matrix{AbstractFloat}(undef, T + steps_ahead, T - 1)
     for t in 1:(T + steps_ahead)
         ξ_matrix[t, :] = t < T ? vcat(ones(t - 1), zeros(T - t)) : ones(T - 1)
     end
@@ -142,7 +145,7 @@ create_ζ(T::Int, steps_ahead::Int, ζ_ω_threshold::Int)::Matrix
 
 """
 function create_ζ(T::Int, steps_ahead::Int, ζ_ω_threshold::Int)::Matrix
-    ζ_matrix = Matrix{Float64}(undef, T + steps_ahead, T - 2)
+    ζ_matrix = Matrix{AbstractFloat}(undef, T + steps_ahead, T - 2)
     for t in 1:(T + steps_ahead)
         ζ_matrix[t, :] = if t < T
             vcat(collect((t - 2):-1:1), zeros(T - 2 - length(collect((t - 2):-1:1))))
@@ -230,9 +233,20 @@ function create_initial_states_Matrix(
 end
 
 """
-create_X(level::Bool, stochastic_level::Bool, trend::Bool, stochastic_trend::Bool,
-                  seasonal::Bool, stochastic_seasonal::Bool, freq_seasonal::Int, outlier::Bool, ζ_ω_threshold::Int, Exogenous_X::Matrix{Fl},
-                  steps_ahead::Int=0, Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, size(Exogenous_X, 2))) where Fl
+create_X(
+    level::Bool,
+    stochastic_level::Bool,
+    trend::Bool,
+    stochastic_trend::Bool,
+    seasonal::Bool,
+    stochastic_seasonal::Bool,
+    freq_seasonal::Int,
+    outlier::Bool,
+    ζ_ω_threshold::Int,
+    Exogenous_X::Matrix{Fl},
+    steps_ahead::Int=0,
+    Exogenous_Forecast::Matrix{Tl}=zeros(steps_ahead, size(Exogenous_X, 2)),
+) where {Fl <: AbstractFloat, Tl <: AbstractFloat}
 
     Creates the StateSpaceLearning matrix X based on the model type and input parameters.
 
@@ -265,8 +279,8 @@ function create_X(
     ζ_ω_threshold::Int,
     Exogenous_X::Matrix{Fl},
     steps_ahead::Int=0,
-    Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, size(Exogenous_X, 2)),
-) where {Fl}
+    Exogenous_Forecast::Matrix{Tl}=zeros(steps_ahead, size(Exogenous_X, 2)),
+) where {Fl <: AbstractFloat, Tl <: AbstractFloat}
     T = size(Exogenous_X, 1)
 
     ξ_matrix = stochastic_level ? create_ξ(T, steps_ahead) : zeros(T + steps_ahead, 0)
@@ -296,7 +310,7 @@ function create_X(
 end
 
 """
-function get_components_indexes(model::StructuralModel)::Dict where Fl
+function get_components_indexes(model::StructuralModel)::Dict
 
     Generates indexes dict for different components based on the model type and input parameters.
 
@@ -308,7 +322,7 @@ function get_components_indexes(model::StructuralModel)::Dict where Fl
 
 """
 function get_components_indexes(model::StructuralModel)::Dict
-    T = length(model.y)
+    T = typeof(model.y) <:Vector ? length(model.y) : size(model.y, 1)
 
     FINAL_INDEX = 0
 
@@ -387,7 +401,12 @@ function get_components_indexes(model::StructuralModel)::Dict
 end
 
 """
-    function get_variances(model::StructuralModel, ε::Vector{Fl}, coefs::Vector{Fl}, components_indexes::Dict{String, Vector{Int}})::Dict where Fl
+get_variances(
+    model::StructuralModel,
+    ε::Vector{Fl},
+    coefs::Vector{Tl},
+    components_indexes::Dict{String,Vector{Int}},
+)::Dict where {Fl, Tl}
 
     Calculates variances for each innovation component and for the residuals.
 
@@ -404,9 +423,9 @@ end
 function get_variances(
     model::StructuralModel,
     ε::Vector{Fl},
-    coefs::Vector{Fl},
+    coefs::Vector{Tl},
     components_indexes::Dict{String,Vector{Int}},
-)::Dict where {Fl}
+)::Dict where {Fl, Tl}
     model_innovations = get_model_innovations(model)
 
     variances = Dict()
@@ -417,6 +436,47 @@ function get_variances(
     return variances
 end
 
+"""
+get_variances(
+    model::StructuralModel,
+    ε::Vector{Vector{Fl}},
+    coefs::Vector{Vector{Tl}},
+    components_indexes::Dict{String,Vector{Int}},
+)::Vector{Dict} where {Fl, Tl}
+
+    Calculates variances for each innovation component and for the residuals.
+
+    # Arguments
+    - `model::StructuralModel`: StructuralModel object.
+    - `ε::Vector{Vector{Fl}}`: Vector of residuals.
+    - `coefs::Vector{Vector{Fl}}`: Vector of coefficients.
+    - `components_indexes::Dict{String, Vector{Int}}`: Dictionary containing indexes for different components.
+
+    # Returns
+    - `Vector{Dict}`: Dictionary containing variances for each innovation component.
+
+"""
+function get_variances(
+    model::StructuralModel,
+    ε::Vector{Vector{Fl}},
+    coefs::Vector{Vector{Tl}},
+    components_indexes::Dict{String,Vector{Int}},
+)::Vector{Dict} where {Fl, Tl}
+    model_innovations = get_model_innovations(model)
+
+    variances_vec = Dict[]
+
+    for i in eachindex(coefs)
+        variances = Dict()
+        for component in model_innovations
+            variances[component] = var(coefs[i][components_indexes[component]])
+        end
+        variances["ε"] = var(ε[i])
+        push!(variances_vec, variances)
+    end
+    return variances_vec
+end
+
 """
     get_model_innovations(model::StructuralModel)::Vector
 
diff --git a/src/structs.jl b/src/structs.jl
index 641b729..706cb40 100644
--- a/src/structs.jl
+++ b/src/structs.jl
@@ -18,4 +18,4 @@ mutable struct Output
     residuals_variances::Dict
     valid_indexes::Vector{Int}
     components::Dict
-end
+end
\ No newline at end of file
diff --git a/src/utils.jl b/src/utils.jl
index c559913..953259f 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,11 +1,13 @@
 """
-    build_components(X::Matrix{Tl}, coefs::Vector{Float64}, components_indexes::Dict{String, Vector{Int}}) -> Dict where Tl
-    
+    build_components(
+    X::Matrix{Tl}, coefs::Vector{Fl}, components_indexes::Dict{String,Vector{Int}}
+)::Dict where {Fl <: AbstractFloat, Tl <: AbstractFloat}
+
     Constructs components dict containing values, indexes and coefficients for each component.
 
     # Arguments
     - X::Matrix{Tl}: Input matrix.
-    - coefs::Vector{Float64}: Coefficients.
+    - coefs::Vector{Fl}: Coefficients.
     - components_indexes::Dict{String, Vector{Int}}: Dictionary mapping component names to their indexes.
 
     # Returns
@@ -16,8 +18,8 @@
 
 """
 function build_components(
-    X::Matrix{Tl}, coefs::Vector{Float64}, components_indexes::Dict{String,Vector{Int}}
-)::Dict where {Tl}
+    X::Matrix{Tl}, coefs::Vector{Fl}, components_indexes::Dict{String,Vector{Int}}
+)::Dict where {Fl <: AbstractFloat, Tl <: AbstractFloat}
     components = Dict()
     for key in keys(components_indexes)
         components[key] = Dict()
@@ -35,36 +37,127 @@ function build_components(
 end
 
 """
-get_fit_and_residuals(estimation_ε::Vector{Float64}, coefs::Vector{Float64}, X::Matrix{Tl}, valid_indexes::Vector{Int}, T::Int) -> Tuple{Vector{Float64}, Vector{Float64}} where Tl
+    build_components(
+    X::Matrix{Tl}, coefs::Vector{Vector{Fl}}, components_indexes::Dict{String,Vector{Int}}
+)::Vector{Dict} where {Fl <: AbstractFloat, Tl <: AbstractFloat}
+
+    Constructs components dict containing values, indexes and coefficients for each component.
+
+    # Arguments
+    - X::Matrix{Fl}: Input matrix.
+    - coefs::Vector{Vector{Fl}}: Coefficients for each time series.
+    - components_indexes::Dict{String, Vector{Int}}: Dictionary mapping component names to their indexes.
+
+    # Returns
+    - components::Vector{Dict}: Dictionary containing components, each represented by a dictionary with keys:
+        - "Coefs": Coefficients for the component.
+        - "Indexes": Indexes associated with the component.
+        - "Values": Values computed from `X` and component coefficients.
+
+"""
+function build_components(
+    X::Matrix{Tl}, coefs::Vector{Vector{Fl}}, components_indexes::Dict{String,Vector{Int}}
+)::Vector{Dict} where {Fl <: AbstractFloat, Tl <: AbstractFloat}
+    components_vec = Dict[]
+    for coef_el in coefs
+        components = Dict()
+        for key in keys(components_indexes)
+            components[key] = Dict()
+            components[key]["Coefs"] = coef_el[components_indexes[key]]
+            components[key]["Indexes"] = components_indexes[key]
+            components[key]["Values"] =
+                X[:, components_indexes[key]] * coef_el[components_indexes[key]]
+        end
+        if haskey(components, "Exogenous_X")
+            components["Exogenous_X"]["Selected"] = findall(
+                i -> i != 0, components["Exogenous_X"]["Coefs"]
+            )
+        end
+        push!(components_vec, components)
+    end
+    return components_vec
+end
+
+"""
+get_fit_and_residuals(
+    estimation_ε::Vector{Fl},
+    coefs::Vector{Pl},
+    X::Matrix{Tl},
+    valid_indexes::Vector{Int},
+    T::Int,
+)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}} where {Fl <: AbstractFloat, Pl <: AbstractFloat, Tl <: AbstractFloat}
 
     Builds complete residuals and fit in sample. Residuals will contain nan values for non valid indexes. Fit in Sample will be a vector of fitted values computed from input data and coefficients (non valid indexes will also be calculated via interpolation).
 
     # Arguments
-    - `estimation_ε::Vector{Float64}`: Vector of estimation errors.
-    - `coefs::Vector{Float64}`: Coefficients.
+    - `estimation_ε::Vector{Fl}`: Vector of estimation errors.
+    - `coefs::Vector{Pl}`: Coefficients.
     - `X::Matrix{Tl}`: Input matrix.
     - `valid_indexes::Vector{Int}`: Valid indexes.
     - `T::Int`: Length of the original time series.
 
     # Returns
     - Tuple containing:
-        - `ε::Vector{Float64}`: Vector containing NaN values filled with estimation errors at valid indexes.
-        - `fitted::Vector{Float64}`: Vector of fitted values computed from input data and coefficients.
+        - `ε::Vector{AbstractFloat}`: Vector containing NaN values filled with estimation errors at valid indexes.
+        - `fitted::Vector{AbstractFloat}`: Vector of fitted values computed from input data and coefficients.
 
 """
 function get_fit_and_residuals(
-    estimation_ε::Vector{Float64},
-    coefs::Vector{Float64},
+    estimation_ε::Vector{Fl},
+    coefs::Vector{Pl},
     X::Matrix{Tl},
     valid_indexes::Vector{Int},
     T::Int,
-)::Tuple{Vector{Float64},Vector{Float64}} where {Tl}
+)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}} where {Fl <: AbstractFloat, Pl <: AbstractFloat, Tl <: AbstractFloat}
     ε = fill(NaN, T)
     ε[valid_indexes] = estimation_ε
     fitted = X * coefs
     return ε, fitted
 end
 
+"""
+get_fit_and_residuals(
+    estimation_ε::Vector{Vector{Fl}},
+    coefs::Vector{Vector{Pl}},
+    X::Matrix{Tl},
+    valid_indexes::Vector{Int},
+    T::Int,
+)::Tuple{Vector{Vector{AbstractFloat}},Vector{Vector{AbstractFloat}}} where {Fl <: AbstractFloat, Pl <: AbstractFloat, Tl <: AbstractFloat}
+
+    Builds complete residuals and fit in sample. Residuals will contain nan values for non valid indexes. Fit in Sample will be a vector of fitted values computed from input data and coefficients (non valid indexes will also be calculated via interpolation).
+
+    # Arguments
+    - `estimation_ε::Vector{Vector{Fl}}`: Vector of estimation errors.
+    - `coefs::Vector{Vector{Pl}}`: Coefficients.
+    - `X::Matrix{Tl}`: Input matrix.
+    - `valid_indexes::Vector{Int}`: Valid indexes.
+    - `T::Int`: Length of the original time series.
+
+    # Returns
+    - Tuple containing:
+        - `ε::Vector{AbstractFloat}`: Vector containing NaN values filled with estimation errors at valid indexes.
+        - `fitted::Vector{AbstractFloat}`: Vector of fitted values computed from input data and coefficients.
+
+"""
+function get_fit_and_residuals(
+    estimation_ε::Vector{Vector{Fl}},
+    coefs::Vector{Vector{Pl}},
+    X::Matrix{Tl},
+    valid_indexes::Vector{Int},
+    T::Int,
+)::Tuple{Vector{Vector{AbstractFloat}},Vector{Vector{AbstractFloat}}} where {Fl <: AbstractFloat, Pl <: AbstractFloat, Tl <: AbstractFloat}
+    ε_vec      = Vector{AbstractFloat}[]
+    fitted_vec = Vector{AbstractFloat}[]
+
+    for i in eachindex(coefs)
+        ε = fill(NaN, T)
+        ε[valid_indexes] = estimation_ε[i]
+        fitted = X * coefs[i]
+        push!(ε_vec, ε)
+        push!(fitted_vec, fitted)
+    end
+    return ε_vec, fitted_vec
+end
 """
 o_size(T::Int)::Int
 
@@ -97,23 +190,25 @@ function create_o_matrix(T::Int, steps_ahead::Int)::Matrix
 end
 
 """
-handle_missing_values(X::Matrix{Tl}, y::Vector{Fl}) -> Tuple{Vector{Fl}, Matrix{Tl}} where {Tl, Fl}
+    handle_missing_values(
+    X::Matrix{Tl}, y::Vector{Fl}
+)::Tuple{Vector{Fl},Matrix{Fl},Vector{Int}} where {Fl <: AbstractFloat, Tl <: AbstractFloat}
 
     Removes missing values from input data and returns the time series and matrix without missing values.
 
     # Arguments
-    - `X::Matrix{Tl}`: Input matrix.
+    - `X::Matrix{Fl}`: Input matrix.
     - `y::Vector{Fl}`: Time series.
 
     # Returns
     - Tuple containing:
         - `y::Vector{Fl}`: Time series without missing values.
-        - `X::Matrix{Tl}`: Input matrix without missing values.
+        - `X::Matrix{Fl}`: Input matrix without missing values.
         - `valid_indexes::Vector{Int}`: Vector containing valid indexes of the time series.
 """
 function handle_missing_values(
     X::Matrix{Tl}, y::Vector{Fl}
-)::Tuple{Vector{Fl},Matrix{Tl},Vector{Int}} where {Tl,Fl}
+)::Tuple{Vector{Fl},Matrix{Fl},Vector{Int}} where {Fl <: AbstractFloat, Tl <: AbstractFloat}
     invalid_indexes = unique(
         vcat([i[1] for i in findall(i -> any(isnan, i), X)], findall(i -> isnan(i), y))
     )
@@ -123,22 +218,58 @@ function handle_missing_values(
 end
 
 """
-has_intercept(X::Matrix{Tl})::Bool where Tl
+handle_missing_values(
+    X::Matrix{Tl}, y::Matrix{Fl}
+)::Tuple{Matrix{Fl},Matrix{Fl},Vector{Int}} where {Fl <: AbstractFloat, Tl <: AbstractFloat}
 
-    Checks if the input matrix has a constant column (intercept).
+    Removes missing values from input data and returns the time series and matrix without missing values.
 
     # Arguments
     - `X::Matrix{Tl}`: Input matrix.
+    - `y::Matrix{Fl}`: Time series.
+
+    # Returns
+    - Tuple containing:
+        - `y::Vector{Fl}`: Time series without missing values.
+        - `X::Matrix{Fl}`: Input matrix without missing values.
+        - `valid_indexes::Vector{Int}`: Vector containing valid indexes of the time series.
+"""
+function handle_missing_values(
+    X::Matrix{Tl}, y::Matrix{Fl}
+)::Tuple{Matrix{Fl},Matrix{Fl},Vector{Int}} where {Fl <: AbstractFloat, Tl <: AbstractFloat}
+    invalid_cartesian_indexes = unique(
+        vcat([i[1] for i in findall(i -> any(isnan, i), X)], findall(i -> isnan(i), y))
+    )
+
+    invalid_indexes = Int64[]
+    for i in invalid_cartesian_indexes
+        if !(i[1] in invalid_indexes)
+            push!(invalid_indexes, i[1])
+        end
+    end
+
+    valid_indexes = setdiff(1:size(y, 1), invalid_indexes)
+
+    return y[valid_indexes, :], X[valid_indexes, :], valid_indexes
+end
+
+"""
+has_intercept(X::Matrix{Fl})::Bool where Fl <: AbstractFloat
+
+    Checks if the input matrix has a constant column (intercept).
+
+    # Arguments
+    - `X::Matrix{Fl}`: Input matrix.
 
     # Returns
     - `Bool`: True if the input matrix has a constant column, false otherwise.
 """
-function has_intercept(X::Matrix{Tl})::Bool where {Tl}
+function has_intercept(X::Matrix{Fl})::Bool where Fl <: AbstractFloat
     return any([all(X[:, i] .== 1) for i in 1:size(X, 2)])
 end
 
 """
-fill_innovation_coefs(model::StructuralModel, T::Int, component::String)::Vector{Float64}
+fill_innovation_coefs(model::StructuralModel, T::Int, component::String, valid_indexes::Vector{Int64})::Vector{AbstractFloat}
 
     Build the innovation coefficients for a given component with same length as the original time series and coefficients attributed to the first observation they are associated with.
 
@@ -146,15 +277,87 @@ fill_innovation_coefs(model::StructuralModel, T::Int, component::String)::Vector
     - `model::StructuralModel`: Structural model.
     - `T::Int`: Length of the original time series.
     - `component::String`: Component name.
+    - `valid_indexes::Vector{Int64}`: Valid Indexes in the time series
 
     # Returns
-    - `Vector{Float64}`: Vector containing innovation coefficients for the given component.
+    - `Union{Vector{AbstractFloat}, Matrix{AbstractFloat}}`: Vector or matrix containing innovation coefficients for the given component.
 """
-function fill_innovation_coefs(model::StructuralModel, component::String)::Vector{Float64}
+function fill_innovation_coefs(model::StructuralModel, component::String, valid_indexes::Vector{Int64})::Union{Vector, Matrix}
     T = length(model.y)
-    inov_comp = zeros(T)
-    for (i, idx) in enumerate(model.output.components[component]["Indexes"])
-        inov_comp[findfirst(i -> i != 0, model.X[:, idx])] = model.output.components[component]["Coefs"][i]
+    if typeof(model.output) == Output 
+        inov_comp = zeros(T)
+        for (i, idx) in enumerate(model.output.components[component]["Indexes"])
+            inov_comp[findfirst(i -> i != 0, model.X[:, idx])] = model.output.components[component]["Coefs"][i]
+        end
+        inov_comp = inov_comp[valid_indexes]
+    else
+        inov_comp = zeros(T, length(model.output))
+        for j in eachindex(model.output)
+            for (i, idx) in enumerate(model.output[j].components[component]["Indexes"])
+                inov_comp[findfirst(i -> i != 0, model.X[:, idx]), j] = model.output[j].components[component]["Coefs"][i]
+            end
+        end
+        inov_comp = inov_comp[valid_indexes, :]
     end
     return inov_comp
 end
+
+"""
+fill_simulation!(simulation::Matrix{Tl}, MV_dist_vec::Vector{MvNormal}, o_noises::Matrix{Fl}, simulation_X::Matrix{Pl}) where {Fl <: AbstractFloat, Pl <: AbstractFloat, Tl <: AbstractFloat}
+
+    Fill the simulation matrix with the generated values
+
+    # Arguments
+    - `simulation::Matrix{Tl}`: Matrix to be filled with simulated values.
+    - `MV_dist_vec::Vector{MvNormal}`: Vector of MvNormal distributions.
+    - `o_noises::Matrix{Fl}`: Matrix of outliers.
+    - `simulation_X::Matrix{Pl}`: Matrix of simulation coefficients.
+"""
+function fill_simulation!(simulation::Matrix{Tl}, MV_dist_vec::Vector{MvNormal}, o_noises::Matrix{Fl}, simulation_X::Matrix{Pl}) where {Fl <: AbstractFloat, Pl <: AbstractFloat, Tl <: AbstractFloat}
+    steps_ahead, N_scenarios = size(simulation)
+    for s in 1:N_scenarios
+        sim_coefs = ones(size(simulation_X, 2)) .* NaN
+
+        for i in 1:steps_ahead
+            rand_inovs = rand(MV_dist_vec[i])
+
+            for comp in eachindex(rand_inovs)
+                sim_coefs[i + (comp - 1) * steps_ahead] = rand_inovs[comp]
+            end
+        end
+
+        simulation[:, s] += (simulation_X * sim_coefs + o_noises[:, s])
+    end
+end
+
+"""
+fill_simulation!(simulation::Vector{Matrix{Tl}}, MV_dist_vec::Vector{MvNormal}, o_noises::Vector{Matrix{Fl}}, simulation_X::Matrix{Pl}, N_innovations::Int) where {Fl <: AbstractFloat, Pl <: AbstractFloat, Tl <: AbstractFloat}
+
+    Fill the simulation matrix with the generated values
+
+    # Arguments
+    - `simulation::Vector{Matrix{Tl}}`: Vector of matrices to be filled with simulated values.
+    - `MV_dist_vec::Vector{MvNormal}`: Vector of MvNormal distributions.
+    - `o_noises::Vector{Matrix{Fl}}`: Vector of matrices of outliers.
+    - `simulation_X::Matrix{Pl}`: Matrix of simulation coefficients.
+    - `N_innovations::Int`: Number of innovations.
+"""
+function fill_simulation!(simulation::Vector{Matrix{Tl}}, MV_dist_vec::Vector{MvNormal}, o_noises::Vector{Matrix{Fl}}, simulation_X::Matrix{Pl}, N_innovations::Int) where {Fl <: AbstractFloat, Pl <: AbstractFloat, Tl <: AbstractFloat}
+    steps_ahead, N_scenarios = size(simulation[1])
+    for j in eachindex(simulation)
+       
+        for s in 1:N_scenarios
+            sim_coefs = ones(size(simulation_X, 2)) .* NaN
+
+            for i in 1:steps_ahead
+                rand_inovs = rand(MV_dist_vec[i])[j:N_innovations:end]
+
+                for comp in eachindex(rand_inovs)
+                    sim_coefs[i + (comp - 1) * steps_ahead] = rand_inovs[comp]
+                end
+            end
+
+            simulation[j][:, s] += (simulation_X * sim_coefs + o_noises[j][:, s])
+        end
+    end
+end
\ No newline at end of file
diff --git a/test/estimation_procedure.jl b/test/estimation_procedure.jl
index d5cb75e..7f122b5 100644
--- a/test/estimation_procedure.jl
+++ b/test/estimation_procedure.jl
@@ -95,26 +95,26 @@ end
     X2 = Basic_Structural_w_level.X
 
     coefs0, ε0 = StateSpaceLearning.fit_lasso(
-        X, y, 0.1, "aic", true, components_indexes, ones(size(X, 2) - 1); rm_average=true
+        X, AbstractFloat.(y), 0.1, "aic", true, components_indexes, ones(size(X, 2) - 1); rm_average=true
     )
     @test length(coefs0) == 43
     @test length(ε0) == 10
 
     coefs1, ε1 = StateSpaceLearning.fit_lasso(
-        X2, y, 0.1, "aic", true, components_indexes2, ones(size(X2, 2)); rm_average=false
+        X2, AbstractFloat.(y), 0.1, "aic", true, components_indexes2, ones(size(X2, 2)); rm_average=false
     )
     @test length(coefs1) == 42
     @test length(ε1) == 10
 
     coefs2, ε2 = StateSpaceLearning.fit_lasso(
-        X, y, 0.1, "aic", true, components_indexes, ones(size(X, 2) - 1); rm_average=true
+        X, AbstractFloat.(y), 0.1, "aic", true, components_indexes, ones(size(X, 2) - 1); rm_average=true
     )
     @test coefs2[1] == mean(y)
     @test length(coefs2) == 43
     @test length(ε2) == 10
 
     coefs3, ε3 = StateSpaceLearning.fit_lasso(
-        X, y, 0.1, "aic", false, components_indexes, ones(size(X, 2) - 1); rm_average=true
+        X, AbstractFloat.(y), 0.1, "aic", false, components_indexes, ones(size(X, 2) - 1); rm_average=true
     )
     @test coefs3[components_indexes["o"][4]] == 0
     @test all(coefs3[components_indexes["Exogenous_X"]] .!= 0)
@@ -123,7 +123,7 @@ end
 
     coefs4, ε4 = StateSpaceLearning.fit_lasso(
         X,
-        y,
+        AbstractFloat.(y),
         0.1,
         "aic",
         true,
diff --git a/test/fit_forecast.jl b/test/fit_forecast.jl
index bfc5d0e..d24570b 100644
--- a/test/fit_forecast.jl
+++ b/test/fit_forecast.jl
@@ -22,6 +22,19 @@
     @test length(model2.output.valid_indexes) == 89
     @test length(model2.output.residuals_variances) == 4
     @test length(keys(model2.output.components)) == 9
+
+    model3 = StateSpaceLearning.StructuralModel(rand(100, 3))
+    StateSpaceLearning.fit!(model3)
+
+    @test length(model3.output) == 3
+    for i in eachindex(model3.output)
+        @test length(model3.output[i].ε) == 100
+        @test length(model3.output[i].fitted) == 100
+        @test length(model3.output[i].coefs) == 375
+        @test length(model3.output[i].valid_indexes) == 100
+        @test length(model3.output[i].residuals_variances) == 4
+        @test length(keys(model3.output[i].components)) == 9
+    end
 end
 
 @testset "Function: forecast" begin
@@ -226,9 +239,26 @@ end
     StateSpaceLearning.fit!(model1)
     @test size(StateSpaceLearning.simulate(model1, 10, 100)) == (10, 100)
 
+    @test size(StateSpaceLearning.simulate(model1, 10, 100; seasonal_innovation_simulation = 10)) == (10, 100)
+
     model2 = StateSpaceLearning.StructuralModel(y2; Exogenous_X=rand(100, 3))
     StateSpaceLearning.fit!(model2)
     @test size(
         StateSpaceLearning.simulate(model2, 10, 100; Exogenous_Forecast=rand(10, 3))
     ) == (10, 100)
+
+    model3 = StateSpaceLearning.StructuralModel(rand(100, 3))
+    StateSpaceLearning.fit!(model3)
+    simulations  = StateSpaceLearning.simulate(model3, 10, 100)
+
+    # test assert error
+    @test_throws AssertionError StateSpaceLearning.simulate(model3, 10, 100; seasonal_innovation_simulation = 10)
+    simulations2 = StateSpaceLearning.simulate(model3, 10, 100; seasonal_innovation_simulation = 3)
+
+    @test length(simulations) == 3
+    @test length(simulations2) == 3
+    for i in eachindex(model3.output)
+        @test size(simulations[i]) == (10, 100)
+        @test size(simulations2[i]) == (10, 100)
+    end
 end
diff --git a/test/utils.jl b/test/utils.jl
index 09466b4..74e1df4 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -110,11 +110,28 @@ end
     StateSpaceLearning.fit!(model)
     components = ["ξ", "ζ", "ω"]
 
-    inov_comp1 = StateSpaceLearning.fill_innovation_coefs(model, components[1])
-    inov_comp2 = StateSpaceLearning.fill_innovation_coefs(model, components[2])
-    inov_comp3 = StateSpaceLearning.fill_innovation_coefs(model, components[3])
+    valid_indexes = model.output.valid_indexes
+
+    inov_comp1 = StateSpaceLearning.fill_innovation_coefs(model, components[1], valid_indexes)
+    inov_comp2 = StateSpaceLearning.fill_innovation_coefs(model, components[2], valid_indexes)
+    inov_comp3 = StateSpaceLearning.fill_innovation_coefs(model, components[3], valid_indexes)
 
     @test length(inov_comp1) == 100
     @test length(inov_comp2) == 100
     @test length(inov_comp3) == 100
+
+    model = StateSpaceLearning.StructuralModel(rand(100, 3))
+    StateSpaceLearning.fit!(model)
+    components = ["ξ", "ζ", "ω"]
+
+    valid_indexes = model.output[1].valid_indexes
+
+    inov_comp1 = StateSpaceLearning.fill_innovation_coefs(model, components[1], valid_indexes)
+    inov_comp2 = StateSpaceLearning.fill_innovation_coefs(model, components[2], valid_indexes)
+    inov_comp3 = StateSpaceLearning.fill_innovation_coefs(model, components[3], valid_indexes)
+
+    @test size(inov_comp1) == (100, 3)
+    @test size(inov_comp2) == (100, 3)
+    @test size(inov_comp3) == (100, 3)
+    
 end