From 56fcae7a412f4b482edd3083e4b7d338f9c04808 Mon Sep 17 00:00:00 2001 From: marinadietze Date: Mon, 31 Mar 2025 11:12:09 -0300 Subject: [PATCH 1/6] exporting components time-series (in-sample) --- src/fit_forecast.jl | 5 +- src/models/structural_model.jl | 100 +++++++++++++++++++++++++++++++++ src/structs.jl | 1 + 3 files changed, 105 insertions(+), 1 deletion(-) diff --git a/src/fit_forecast.jl b/src/fit_forecast.jl index 5836d7f..4982b54 100644 --- a/src/fit_forecast.jl +++ b/src/fit_forecast.jl @@ -63,8 +63,10 @@ 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 +79,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..33d7288 100644 --- a/src/models/structural_model.jl +++ b/src/models/structural_model.jl @@ -956,4 +956,104 @@ function get_innovation_simulation_X( end end +function get_μ(model::StructuralModel, components::Dict, ν::Vector{Float64}) + + T = size(model.y, 1) + + μ = Vector{Float64}(undef, T) + μ[1] = components["μ1"]["Coefs"][1] + ξ = vcat(0.0, components["ξ"]["Coefs"], 0.0) + + for t in 2:T + μ[t] = μ[t - 1] + ν[t - 1] + ξ[t] + end + + return μ, ξ +end + +function get_ν(model::StructuralModel, components::Dict) + + T = size(model.y, 1) + ζ_ω_threshold = model.ζ_ω_threshold + + ν = Vector{Float64}(undef, T) + ν[1] = components["ν1"]["Coefs"][1] + ζ = vcat(0.0, 0.0, components["ζ"]["Coefs"], zeros(ζ_ω_threshold)) + + for t in 2:T + ν[t] = ν[t - 1] + ζ[t] + end + + return ν, ζ +end + +function get_γ(model::StructuralModel, components::Dict, s::Int64) + + T = size(model.y, 1) + ζ_ω_threshold = model.ζ_ω_threshold + + γ = Vector{Float64}(undef, T) + γ[1:s] = components["γ1_$(s)"]["Coefs"] + ω = vcat(zeros(s - 1), components["ω_$(s)"]["Coefs"], zeros(ζ_ω_threshold)) + + for t in s + 1:T + γ[t] = -sum(γ[t - j] for j in 1:s - 1) + ω[t] + end + + return γ, ω +end + +function get_components_ts(model::StructuralModel, components::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 + +function get_components_ts(model::StructuralModel, components::Vector{Dict}) + + 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 + 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 From 26c82191820940221264b38e187c486ec79f27fb Mon Sep 17 00:00:00 2001 From: marinadietze Date: Mon, 31 Mar 2025 11:20:54 -0300 Subject: [PATCH 2/6] fixed some types --- src/models/structural_model.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/models/structural_model.jl b/src/models/structural_model.jl index 33d7288..7218e21 100644 --- a/src/models/structural_model.jl +++ b/src/models/structural_model.jl @@ -956,11 +956,11 @@ function get_innovation_simulation_X( end end -function get_μ(model::StructuralModel, components::Dict, ν::Vector{Float64}) +function get_μ(model::StructuralModel, components::Dict, ν::Vector{AbstractFloat}) T = size(model.y, 1) - μ = Vector{Float64}(undef, T) + μ = Vector{AbstractFloat}(undef, T) μ[1] = components["μ1"]["Coefs"][1] ξ = vcat(0.0, components["ξ"]["Coefs"], 0.0) @@ -976,7 +976,7 @@ function get_ν(model::StructuralModel, components::Dict) T = size(model.y, 1) ζ_ω_threshold = model.ζ_ω_threshold - ν = Vector{Float64}(undef, T) + ν = Vector{AbstractFloat}(undef, T) ν[1] = components["ν1"]["Coefs"][1] ζ = vcat(0.0, 0.0, components["ζ"]["Coefs"], zeros(ζ_ω_threshold)) @@ -987,12 +987,12 @@ function get_ν(model::StructuralModel, components::Dict) return ν, ζ end -function get_γ(model::StructuralModel, components::Dict, s::Int64) +function get_γ(model::StructuralModel, components::Dict, s::Int) T = size(model.y, 1) ζ_ω_threshold = model.ζ_ω_threshold - γ = Vector{Float64}(undef, T) + γ = Vector{AbstractFloat}(undef, T) γ[1:s] = components["γ1_$(s)"]["Coefs"] ω = vcat(zeros(s - 1), components["ω_$(s)"]["Coefs"], zeros(ζ_ω_threshold)) From 928bddfba34fe24fd5711f44ce22defdcf235660 Mon Sep 17 00:00:00 2001 From: marinadietze Date: Mon, 31 Mar 2025 11:37:32 -0300 Subject: [PATCH 3/6] added documentation --- src/models/structural_model.jl | 77 +++++++++++++++++++++++++++++++--- 1 file changed, 72 insertions(+), 5 deletions(-) diff --git a/src/models/structural_model.jl b/src/models/structural_model.jl index 7218e21..b235a2c 100644 --- a/src/models/structural_model.jl +++ b/src/models/structural_model.jl @@ -956,7 +956,21 @@ function get_innovation_simulation_X( end end -function get_μ(model::StructuralModel, components::Dict, ν::Vector{AbstractFloat}) +""" + 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) @@ -971,7 +985,20 @@ function get_μ(model::StructuralModel, components::Dict, ν::Vector{AbstractFlo return μ, ξ end -function get_ν(model::StructuralModel, components::Dict) +""" + 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) ζ_ω_threshold = model.ζ_ω_threshold @@ -987,7 +1014,21 @@ function get_ν(model::StructuralModel, components::Dict) return ν, ζ end -function get_γ(model::StructuralModel, components::Dict, s::Int) +""" + 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) ζ_ω_threshold = model.ζ_ω_threshold @@ -1003,7 +1044,20 @@ function get_γ(model::StructuralModel, components::Dict, s::Int) return γ, ω end -function get_components_ts(model::StructuralModel, components::Dict) +""" + 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() @@ -1026,7 +1080,20 @@ function get_components_ts(model::StructuralModel, components::Dict) return components_ts_dict end -function get_components_ts(model::StructuralModel, components::Vector{Dict}) +""" + 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} components_ts = [] freq_seasonal = model.freq_seasonal From 735ad9f4107c5569eebfb83203d5000d908ef3ae Mon Sep 17 00:00:00 2001 From: dietzemarina Date: Thu, 3 Apr 2025 17:21:19 -0300 Subject: [PATCH 4/6] generalize the presence of components --- src/models/structural_model.jl | 77 +++++++++++++++++++++++----------- 1 file changed, 53 insertions(+), 24 deletions(-) diff --git a/src/models/structural_model.jl b/src/models/structural_model.jl index b235a2c..7850ef4 100644 --- a/src/models/structural_model.jl +++ b/src/models/structural_model.jl @@ -970,13 +970,23 @@ end - `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}} - +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 - μ = Vector{AbstractFloat}(undef, T) - μ[1] = components["μ1"]["Coefs"][1] - ξ = vcat(0.0, components["ξ"]["Coefs"], 0.0) + 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] @@ -998,14 +1008,25 @@ end - `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}} +function get_ν( + model::StructuralModel, components::Dict +)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}} + T = size(model.y, 1) + ν = Vector{AbstractFloat}(undef, T) - T = size(model.y, 1) ζ_ω_threshold = model.ζ_ω_threshold - ν = Vector{AbstractFloat}(undef, T) - ν[1] = components["ν1"]["Coefs"][1] - ζ = vcat(0.0, 0.0, components["ζ"]["Coefs"], zeros(ζ_ω_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] @@ -1028,17 +1049,28 @@ end - `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}} +function get_γ( + model::StructuralModel, components::Dict, s::Int +)::Tuple{Vector{AbstractFloat},Vector{AbstractFloat}} + T = size(model.y, 1) + γ = Vector{AbstractFloat}(undef, T) - T = size(model.y, 1) ζ_ω_threshold = model.ζ_ω_threshold - γ = Vector{AbstractFloat}(undef, T) - γ[1:s] = components["γ1_$(s)"]["Coefs"] - ω = vcat(zeros(s - 1), components["ω_$(s)"]["Coefs"], zeros(ζ_ω_threshold)) - - for t in s + 1:T - γ[t] = -sum(γ[t - j] for j in 1:s - 1) + ω[t] + 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 γ, ω @@ -1058,15 +1090,14 @@ end """ function get_components_ts(model::StructuralModel, components::Dict)::Dict - - freq_seasonal = model.freq_seasonal + 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["ζ"] = ζ components_ts_dict["ξ"] = ξ @@ -1094,19 +1125,17 @@ end """ function get_components_ts(model::StructuralModel, components::Vector{Dict})::Vector{Dict} - 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["ζ"] = ζ components_ts_dict["ξ"] = ξ From 1709afb5b7270b1a3885f0454195591860ab81bb Mon Sep 17 00:00:00 2001 From: dietzemarina Date: Thu, 3 Apr 2025 17:25:08 -0300 Subject: [PATCH 5/6] format fit_forecast file --- src/fit_forecast.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/fit_forecast.jl b/src/fit_forecast.jl index 4982b54..c7689c8 100644 --- a/src/fit_forecast.jl +++ b/src/fit_forecast.jl @@ -66,7 +66,9 @@ function fit!( components_ts = get_components_ts(model, components) if typeof(model.y) <: Vector - output = Output(coefs, ε, fitted, residuals_variances, valid_indexes, components, components_ts) + output = Output( + coefs, ε, fitted, residuals_variances, valid_indexes, components, components_ts + ) else output = Output[] for i in eachindex(coefs) From f4c08b8fd259165c38380afa6f3b72d4e45abf80 Mon Sep 17 00:00:00 2001 From: dietzemarina Date: Thu, 3 Apr 2025 19:12:43 -0300 Subject: [PATCH 6/6] exporting components scenarios --- src/StateSpaceLearning.jl | 2 +- src/models/structural_model.jl | 83 +++++++++++++++++++++++++++++++++- 2 files changed, 83 insertions(+), 2 deletions(-) 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/models/structural_model.jl b/src/models/structural_model.jl index 7850ef4..f8ae63c 100644 --- a/src/models/structural_model.jl +++ b/src/models/structural_model.jl @@ -1124,7 +1124,9 @@ end - `Vector{Dict}`: Vector of dictionaries with the time-series states and innovations. """ -function get_components_ts(model::StructuralModel, components::Vector{Dict})::Vector{Dict} +function get_components_ts( + model::StructuralModel, components::Vector{Dict} +)::Vector{Dict} where {Al} components_ts = [] freq_seasonal = model.freq_seasonal @@ -1152,4 +1154,83 @@ function get_components_ts(model::StructuralModel, components::Vector{Dict})::Ve 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