diff --git a/src/StateSpaceLearning.jl b/src/StateSpaceLearning.jl
index d19d07e..baee639 100644
--- a/src/StateSpaceLearning.jl
+++ b/src/StateSpaceLearning.jl
@@ -1,6 +1,6 @@
 module StateSpaceLearning
 
-using LinearAlgebra, Statistics, GLMNet, Distributions, SparseArrays
+using LinearAlgebra, Statistics, GLMNet, Distributions, SparseArrays, Random
 
 abstract type StateSpaceLearningModel end
 
diff --git a/src/fit_forecast.jl b/src/fit_forecast.jl
index 5836d7f..c7689c8 100644
--- a/src/fit_forecast.jl
+++ b/src/fit_forecast.jl
@@ -63,8 +63,12 @@ function fit!(
 
     ε, fitted = get_fit_and_residuals(estimation_ε, coefs, model.X, valid_indexes, T)
 
+    components_ts = get_components_ts(model, components)
+
     if typeof(model.y) <: Vector
-        output = Output(coefs, ε, fitted, residuals_variances, valid_indexes, components)
+        output = Output(
+            coefs, ε, fitted, residuals_variances, valid_indexes, components, components_ts
+        )
     else
         output = Output[]
         for i in eachindex(coefs)
@@ -77,6 +81,7 @@ function fit!(
                     residuals_variances[i],
                     valid_indexes,
                     components[i],
+                    components_ts[i],
                 ),
             )
         end
diff --git a/src/models/structural_model.jl b/src/models/structural_model.jl
index 1f79683..f8ae63c 100644
--- a/src/models/structural_model.jl
+++ b/src/models/structural_model.jl
@@ -956,4 +956,281 @@ function get_innovation_simulation_X(
     end
 end
 
+"""
+    get_μ(model::StructuralModel, components::Dict, ν::Vector{AbstractFloat})::Tuple{Vector{AbstractFloat}, Vector{AbstractFloat}}
+
+    Returns the level component and associated innovation vectors.
+
+    # Arguments
+    - `model::StructuralModel`: StructuralModel object.
+    - `components::Dict` : Components dict.
+    - `ν::Vector{AbstractFloat}`: Time-series of the slope component.
+
+    # Returns
+    - `Tuple{Vector{AbstractFloat}, Vector{AbstractFloat}}`: Tuple of vectors containing the time series state and innovations.
+
+"""
+function get_μ(
+    model::StructuralModel, components::Dict, ν::Vector{AbstractFloat}
+)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}}
+    T = size(model.y, 1)
+    μ = Vector{AbstractFloat}(undef, T)
+
+    if model.level
+        μ[1] = components["μ1"]["Coefs"][1]
+    else
+        μ[1] = 0.0
+    end
+
+    if model.stochastic_level
+        ξ = vcat(0.0, components["ξ"]["Coefs"], 0.0)
+    else
+        ξ = zeros(AbstractFloat, T)
+    end
+
+    for t in 2:T
+        μ[t] = μ[t - 1] + ν[t - 1] + ξ[t]
+    end
+
+    return μ, ξ
+end
+
+"""
+    get_ν(model::StructuralModel, components::Dict)::Tuple{Vector{AbstractFloat}, Vector{AbstractFloat}} 
+
+    Returns the slope component and associated innovation vectors.
+
+    # Arguments
+    - `model::StructuralModel`: StructuralModel object.
+    - `components::Dict`: Components dict..
+
+    # Returns
+    - `Tuple{Vector{AbstractFloat}, Vector{AbstractFloat}}`: Tuple of vectors containing the time series state and innovations.
+
+"""
+function get_ν(
+    model::StructuralModel, components::Dict
+)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}}
+    T = size(model.y, 1)
+    ν = Vector{AbstractFloat}(undef, T)
+
+    ζ_ω_threshold = model.ζ_ω_threshold
+
+    if model.trend
+        ν[1] = components["ν1"]["Coefs"][1]
+    else
+        ν[1] = 0.0
+    end
+
+    if model.stochastic_trend
+        ζ = vcat(0.0, 0.0, components["ζ"]["Coefs"], zeros(ζ_ω_threshold))
+    else
+        ζ = zeros(AbstractFloat, T)
+    end
+
+    for t in 2:T
+        ν[t] = ν[t - 1] + ζ[t]
+    end
+
+    return ν, ζ
+end
+
+"""
+    get_γ(model::StructuralModel, components::Dict, s::Int)::Tuple{Vector{AbstractFloat}, Vector{AbstractFloat}} 
+
+    Returns the seasonality component and associated innovation vectors.
+
+    # Arguments
+    - `model::StructuralModel`: StructuralModel object.
+    - `components::Dict`: Components dict.
+    - `s::Int`: Seasonal frequency.	
+
+    # Returns
+    - `Tuple{Vector{AbstractFloat}, Vector{AbstractFloat}}`: Tuple of vectors containing the time series state and innovations.
+
+"""
+function get_γ(
+    model::StructuralModel, components::Dict, s::Int
+)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}}
+    T = size(model.y, 1)
+    γ = Vector{AbstractFloat}(undef, T)
+
+    ζ_ω_threshold = model.ζ_ω_threshold
+
+    if model.seasonal
+        γ[1:s] = components["γ1_$(s)"]["Coefs"]
+    else
+        γ[1:s] = zeros(AbstractFloat, s)
+    end
+
+    if model.stochastic_seasonal
+        ω = vcat(zeros(s - 1), components["ω_$(s)"]["Coefs"], zeros(ζ_ω_threshold))
+    else
+        ω = zeros(AbstractFloat, T)
+    end
+
+    for t in (s + 1):T
+        γ[t] = -sum(γ[t - j] for j in 1:(s - 1)) + ω[t]
+    end
+
+    return γ, ω
+end
+
+"""
+    get_components_ts(model::StructuralModel, components::Dict)::Dict  
+
+    Returns a dictionary with the time series state and innovations for each component.
+
+    # Arguments
+    - `model::StructuralModel`: StructuralModel object.
+    - `components::Dict`: Components dict.
+
+    # Returns
+    - `Dict`: Dictionary of time-series states and innovations.
+
+"""
+function get_components_ts(model::StructuralModel, components::Dict)::Dict
+    freq_seasonal = model.freq_seasonal
+    components_ts_dict = Dict()
+
+    ν, ζ = get_ν(model, components)
+    μ, ξ = get_μ(model, components, ν)
+
+    components_ts_dict["ν"] = ν
+    components_ts_dict["μ"] = μ
+    components_ts_dict["ζ"] = ζ
+    components_ts_dict["ξ"] = ξ
+
+    for s in freq_seasonal
+        γ, ω = get_γ(model, components, s)
+
+        components_ts_dict["γ_$s"] = γ
+        components_ts_dict["ω_$s"] = ω
+    end
+
+    return components_ts_dict
+end
+
+"""
+    get_components_ts(model::StructuralModel, components::Vector{Dict})::Vector{Dict}  
+
+    Returns a vector of dictionaries with the time series state and innovations for each component, of each dependent time-series.
+
+    # Arguments
+    - `model::StructuralModel`: StructuralModel object.
+    - `components::Vector{Dict}`: Vector of components dict.
+
+    # Returns
+    - `Vector{Dict}`: Vector of dictionaries with the time-series states and innovations.
+
+"""
+function get_components_ts(
+    model::StructuralModel, components::Vector{Dict}
+)::Vector{Dict} where {Al}
+    components_ts = []
+    freq_seasonal = model.freq_seasonal
+
+    for component in components
+        components_ts_dict = Dict()
+
+        ν, ζ = get_ν(model, component)
+        μ, ξ = get_μ(model, component, ν)
+
+        components_ts_dict["ν"] = ν
+        components_ts_dict["μ"] = μ
+        components_ts_dict["ζ"] = ζ
+        components_ts_dict["ξ"] = ξ
+
+        for s in freq_seasonal
+            γ, ω = get_γ(model, component, s)
+
+            components_ts_dict["γ_$s"] = γ
+            components_ts_dict["ω_$s"] = ω
+        end
+
+        push!(components_ts, components_ts_dict)
+    end
+
+    return components_ts
+end
+
+"""
+    simulate_states(model::StructuralModel, steps_ahead::Int, N_scenarios::Int)::Vector{Dict}  
+
+    Returns a vector of dictionaries with the scenarios of each component, for each dependent time-series.
+
+    # Arguments
+    - `model::StructuralModel`: StructuralModel object.
+    - `steps_ahead::Int`: Number of steps ahead for forecasting.
+    - `N_scenarios::Int`: Number of scenarios.
+
+    # Returns
+    - `Vector{Dict}`: Vector of dictionaries with the scenarios of the components of each dependent time-series.
+
+"""
+function simulate_states(
+    model::StructuralModel, steps_ahead::Int, N_scenarios::Int
+)::Vector{Dict}
+    Random.seed!(1234)
+
+    freq_seasonal = model.freq_seasonal
+    T = size(model.y, 1)
+    idxs = rand(1:T, steps_ahead, N_scenarios)
+
+    components_ts_vec = []
+    scenarios_ts = []
+    if typeof(model.output) == Vector{StateSpaceLearning.Output}
+        for i in eachindex(model.output)
+            push!(components_ts_vec, get_components_ts(model, model.output[i].components))
+        end
+    else
+        push!(components_ts_vec, get_components_ts(model, model.output.components))
+    end
+
+    for components_ts in components_ts_vec
+        scenarios_ts_dict = Dict()
+
+        scenarios_ts_dict["μ"] = Matrix{AbstractFloat}(undef, steps_ahead, N_scenarios)
+        scenarios_ts_dict["ν"] = Matrix{AbstractFloat}(undef, steps_ahead, N_scenarios)
+
+        for f in freq_seasonal
+            scenarios_ts_dict["γ_$f"] = Matrix{AbstractFloat}(
+                undef, steps_ahead, N_scenarios
+            )
+        end
+
+        for s in 1:N_scenarios
+            ν_vec = vcat(components_ts["ν"], zeros(steps_ahead))
+            μ_vec = vcat(components_ts["μ"], zeros(steps_ahead))
+            γ_matrix = vcat(
+                hcat([components_ts["γ_$s"] for s in freq_seasonal]...),
+                zeros(steps_ahead, length(freq_seasonal)),
+            )
+
+            for t in 1:steps_ahead
+                ν_vec[T + t] = ν_vec[T + t - 1] + components_ts["ζ"][idxs[t, s]]
+                μ_vec[T + t] =
+                    μ_vec[T + t - 1] + ν_vec[T + t - 1] + components_ts["ξ"][idxs[t, s]]
+
+                for (i, f) in enumerate(freq_seasonal)
+                    γ_matrix[T + t, i] =
+                        -sum(γ_matrix[T + t - j, i] for j in 1:(f - 1)) +
+                        components_ts["ω_$f"][idxs[t, s]]
+                end
+            end
+
+            scenarios_ts_dict["μ"][:, s] = μ_vec[(T + 1):(T + steps_ahead)]
+            scenarios_ts_dict["ν"][:, s] = ν_vec[(T + 1):(T + steps_ahead)]
+
+            for (i, f) in enumerate(freq_seasonal)
+                scenarios_ts_dict["γ_$f"][:, s] = γ_matrix[(T + 1):(T + steps_ahead), i]
+            end
+        end
+
+        push!(scenarios_ts, scenarios_ts_dict)
+    end
+
+    return scenarios_ts
+end
+
 isfitted(model::StructuralModel) = isnothing(model.output) ? false : true
diff --git a/src/structs.jl b/src/structs.jl
index 641b729..18384dd 100644
--- a/src/structs.jl
+++ b/src/structs.jl
@@ -18,4 +18,5 @@ mutable struct Output
     residuals_variances::Dict
     valid_indexes::Vector{Int}
     components::Dict
+    components_ts::Dict
 end