diff --git a/.gitignore b/.gitignore index 151d291..a9a71f4 100644 --- a/.gitignore +++ b/.gitignore @@ -15,4 +15,8 @@ longhorizon/ test.jl test_ms.jl new/ -test2.jl \ No newline at end of file +test2.jl + +src/boosting_estimation_procedure.jl +plots/ +paper_tests/m4_test/evaluate_model2.jl \ No newline at end of file diff --git a/Project.toml b/Project.toml index e53f72d..aede96b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StateSpaceLearning" uuid = "971c4b7c-2c4e-4bac-8525-e842df3cde7b" authors = ["andreramosfc "] -version = "2.0.5" +version = "2.0.6" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" diff --git a/README.md b/README.md index e06d4bb..e05603f 100644 --- a/README.md +++ b/README.md @@ -36,6 +36,7 @@ simulation = simulate(model, 12, 1000) # Gets 1000 scenarios path of 12 steps ah * `freq_seasonal::Union{Int,Vector{Int}}`: Seasonal frequency or vector of frequencies (default: 12) * `cycle_period::Union{Union{Int,<:AbstractFloat},Vector{Int},Vector{<:AbstractFloat}}`: Cycle period or vector of periods (default: 0) * `outlier::Bool`: Include outlier component (default: true) +* `ξ_threshold::Int`: Threshold for level innovations (default: 1) * `ζ_threshold::Int`: Threshold for slope innovations (default: 12) * `ω_threshold::Int`: Threshold for seasonal innovations (default: 12) * `ϕ_threshold::Int`: Threshold for cycle innovations (default: 12) diff --git a/src/StateSpaceLearning.jl b/src/StateSpaceLearning.jl index e3dfb01..e09d02f 100644 --- a/src/StateSpaceLearning.jl +++ b/src/StateSpaceLearning.jl @@ -17,6 +17,7 @@ export fit!, forecast, simulate, StructuralModel, + BasicStructuralModel, plot_point_forecast, plot_scenarios, simulate_states diff --git a/src/models/structural_model.jl b/src/models/structural_model.jl index 80cf1c9..813847f 100644 --- a/src/models/structural_model.jl +++ b/src/models/structural_model.jl @@ -65,6 +65,7 @@ mutable struct StructuralModel <: StateSpaceLearningModel freq_seasonal::Union{Int,Vector{Int}} cycle_period::Union{Union{Int,<:AbstractFloat},Vector{Int},Vector{<:AbstractFloat}} outlier::Bool + ξ_threshold::Int ζ_threshold::Int ω_threshold::Int ϕ_threshold::Int @@ -82,6 +83,7 @@ mutable struct StructuralModel <: StateSpaceLearningModel freq_seasonal::Union{Int,Vector{Int}}=12, cycle_period::Union{Union{Int,<:AbstractFloat},Vector{Int},Vector{<:AbstractFloat}}=0, outlier::Bool=true, + ξ_threshold::Int=1, ζ_threshold::Int=12, ω_threshold::Int=12, ϕ_threshold::Int=12, @@ -156,6 +158,7 @@ mutable struct StructuralModel <: StateSpaceLearningModel freq_seasonal, cycle_period, outlier, + ξ_threshold, ζ_threshold, ω_threshold, ϕ_threshold, @@ -180,6 +183,7 @@ mutable struct StructuralModel <: StateSpaceLearningModel freq_seasonal, cycle_period, outlier, + ξ_threshold, ζ_threshold, ω_threshold, ϕ_threshold, @@ -191,6 +195,45 @@ mutable struct StructuralModel <: StateSpaceLearningModel end end +""" +Create a Basic Structural Model. + +# Arguments +- `y::Union{Vector,Matrix}`: Time series to be modeled. +- `level::String="stochastic"`: Level component of the model. +- `slope::String="stochastic"`: Slope component of the model. +- `seasonal::String="stochastic"`: Seasonal component of the model. +- `freq_seasonal::Union{Int,Vector{Int}}=12`: Seasonal period of the model. + +# Returns +- `StructuralModel`: Basic Structural Model. +""" +function BasicStructuralModel( + y::Union{Vector,Matrix}; + level::String="stochastic", + slope::String="stochastic", + seasonal::String="stochastic", + freq_seasonal::Union{Int,Vector{Int}}=12, +) + return StructuralModel( + y; + level=level, + slope=slope, + seasonal=seasonal, + cycle="none", + freq_seasonal=freq_seasonal, + cycle_period=0, + outlier=false, + ξ_threshold=0, + ζ_threshold=0, + ω_threshold=0, + ϕ_threshold=0, + stochastic_start=1, + exog=zeros(length(y), 0), + dynamic_exog_coefs=nothing, + ) +end + """ ξ_size(T::Int, stochastic_start::Int)::Int @@ -204,7 +247,8 @@ end - `Int`: Size of ξ calculated from T. """ -ξ_size(T::Int, stochastic_start::Int)::Int = T - max(2, stochastic_start) +ξ_size(T::Int, ξ_threshold::Int, stochastic_start::Int)::Int = + T - max(2, stochastic_start) + 1 - ξ_threshold """ ζ_size(T::Int, ζ_threshold::Int, stochastic_start::Int)::Int @@ -273,27 +317,28 @@ o_size(T::Int, stochastic_start::Int)::Int = T - max(1, stochastic_start) + 1 (2 * (T - max(2, stochastic_start) + 1) - (max(1, ϕ_threshold) * 2)) """ - create_ξ(T::Int, stochastic_start::Int)::Matrix + create_ξ(T::Int, ξ_threshold::Int, stochastic_start::Int)::Matrix Creates a matrix of innovations ξ based on the input sizes, and the desired steps ahead (this is necessary for the forecast function) # Arguments - `T::Int`: Length of the original time series. + - `ξ_threshold::Int`: Stabilize parameter ξ. - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts. # Returns - `Matrix`: Matrix of innovations ξ constructed based on the input sizes. """ -function create_ξ(T::Int, stochastic_start::Int)::Matrix +function create_ξ(T::Int, ξ_threshold::Int, stochastic_start::Int)::Matrix stochastic_start = max(2, stochastic_start) ξ_matrix = zeros(T, T - stochastic_start + 1) ones_indexes = findall( I -> Tuple(I)[1] - (stochastic_start - 2) > Tuple(I)[2], - CartesianIndices((T, T - stochastic_start)), + CartesianIndices((T, T - stochastic_start + 1)), ) ξ_matrix[ones_indexes] .= 1 - return ξ_matrix[:, 1:(end - 1)] + return ξ_matrix[:, 1:(end - ξ_threshold)] end """ @@ -524,6 +569,7 @@ create_dynamic_exog_coefs_matrix(dynamic_exog_coefs::Vector{<:Tuple}, T::Int,ζ_ function create_dynamic_exog_coefs_matrix( dynamic_exog_coefs::Vector{<:Tuple}, T::Int, + ξ_threshold::Int, ζ_threshold::Int, ω_threshold::Int, ϕ_threshold::Int, @@ -537,7 +583,7 @@ function create_dynamic_exog_coefs_matrix( nothing else state_components_dict["level"] = hcat( - ones(T, 1), create_ξ(T, stochastic_start) + ones(T, 1), create_ξ(T, ξ_threshold, stochastic_start) ) end key_name = "level" @@ -588,6 +634,7 @@ create_forecast_dynamic_exog_coefs_matrix(dynamic_exog_coefs::Vector{<:Tuple}, T - `dynamic_exog_coefs::Vector{<:Tuple}`: Vector of tuples containing the combination components. - `T::Int`: Length of the original time series. - `steps_ahead::Int`: Steps ahead. + - `ξ_threshold::Int`: Stabilize parameter ξ. - `ζ_threshold::Int`: Stabilize parameter ζ. - `ω_threshold::Int`: Stabilize parameter ω. - `ϕ_threshold::Int`: Stabilize parameter ϕ. @@ -600,6 +647,7 @@ function create_forecast_dynamic_exog_coefs_matrix( dynamic_exog_coefs::Vector{<:Tuple}, T::Int, steps_ahead::Int, + ξ_threshold::Int, ζ_threshold::Int, ω_threshold::Int, ϕ_threshold::Int, @@ -613,7 +661,8 @@ function create_forecast_dynamic_exog_coefs_matrix( nothing else state_components_dict["level"] = hcat( - ones(T + steps_ahead, 1), create_ξ(T + steps_ahead, stochastic_start) + ones(T + steps_ahead, 1), + create_ξ(T + steps_ahead, ξ_threshold, stochastic_start), )[ (end - steps_ahead + 1):end, 1:combination[4] ] @@ -680,6 +729,7 @@ create_X( freq_seasonal::Union{Int,Vector{Int}}, cycle_period::Union{Int,Vector{Int}}, outlier::Bool, + ξ_threshold::Int, ζ_threshold::Int, ω_threshold::Int, ϕ_threshold::Int, @@ -701,6 +751,7 @@ create_X( - `freq_seasonal::Union{Int, Vector{Int}}`: Seasonal period. - `cycle_period::Union{Int,Vector{Int}}`: Cycle period. - `outlier::Bool`: Flag for considering outlier component. + - `ξ_threshold::Int`: Stabilize parameter ξ. - `ζ_threshold::Int`: Stabilize parameter ζ. - `ω_threshold::Int`: Stabilize parameter ω. - `ϕ_threshold::Int`: Stabilize parameter ϕ. @@ -722,6 +773,7 @@ function create_X( freq_seasonal::Union{Int,Vector{Int}}, cycle_period::Union{Union{Int,<:AbstractFloat},Vector{Int},Vector{<:AbstractFloat}}, outlier::Bool, + ξ_threshold::Int, ζ_threshold::Int, ω_threshold::Int, ϕ_threshold::Int, @@ -732,7 +784,7 @@ function create_X( T = size(exog, 1) ξ_matrix = if stochastic_level - create_ξ(T, stochastic_start) + create_ξ(T, ξ_threshold, stochastic_start) else zeros(T, 0) end @@ -770,6 +822,7 @@ function create_X( create_dynamic_exog_coefs_matrix( dynamic_exog_coefs, T, + ξ_threshold, ζ_threshold, ω_threshold, ϕ_threshold, @@ -847,7 +900,9 @@ function get_components_indexes(model::StructuralModel)::Dict if model.stochastic_level ξ_indexes = collect( - (FINAL_INDEX + 1):(FINAL_INDEX + ξ_size(T, model.stochastic_start)) + (FINAL_INDEX + 1):(FINAL_INDEX + ξ_size( + T, model.ξ_threshold, model.stochastic_start + )), ) FINAL_INDEX += length(ξ_indexes) else @@ -908,7 +963,6 @@ function get_components_indexes(model::StructuralModel)::Dict FINAL_INDEX += length(exogenous_indexes) dynamic_exog_coefs_indexes = collect((FINAL_INDEX + 1):size(model.X, 2)) - FINAL_INDEX += length(dynamic_exog_coefs_indexes) components_indexes_dict = Dict( "μ1" => μ1_indexes, @@ -940,6 +994,43 @@ function get_components_indexes(model::StructuralModel)::Dict end end + if !isnothing(model.dynamic_exog_coefs) + for i in eachindex(model.dynamic_exog_coefs) + if model.dynamic_exog_coefs[i][2] == "level" + components_indexes_dict["dynamic_exog_coefs_$(i)"] = collect( + (FINAL_INDEX + 1):(FINAL_INDEX + 1 + ξ_size( + T, model.ξ_threshold, model.stochastic_start + )), + ) + FINAL_INDEX += length(components_indexes_dict["dynamic_exog_coefs_$(i)"]) + elseif model.dynamic_exog_coefs[i][2] == "slope" + components_indexes_dict["dynamic_exog_coefs_$(i)"] = collect( + (FINAL_INDEX + 1):(FINAL_INDEX + 1 + ζ_size( + T, model.ζ_threshold, model.stochastic_start + )), + ) + FINAL_INDEX += length(components_indexes_dict["dynamic_exog_coefs_$(i)"]) + elseif model.dynamic_exog_coefs[i][2] == "seasonal" + components_indexes_dict["dynamic_exog_coefs_$(i)"] = collect( + (FINAL_INDEX + 1):(FINAL_INDEX + model.dynamic_exog_coefs[i][3] + ω_size( + T, + model.dynamic_exog_coefs[i][3], + model.ω_threshold, + model.stochastic_start, + )), + ) + FINAL_INDEX += length(components_indexes_dict["dynamic_exog_coefs_$(i)"]) + elseif model.dynamic_exog_coefs[i][2] == "cycle" + components_indexes_dict["dynamic_exog_coefs_$(i)"] = collect( + (FINAL_INDEX + 1):(FINAL_INDEX + 2 + ϕ_size( + T, model.ϕ_threshold, model.stochastic_start + )), + ) + FINAL_INDEX += length(components_indexes_dict["dynamic_exog_coefs_$(i)"]) + end + end + end + return components_indexes_dict end @@ -979,47 +1070,6 @@ 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 @@ -1053,6 +1103,12 @@ function get_model_innovations(model::StructuralModel) push!(model_innovations, "ϕ_$i") end end + + if !isnothing(model.dynamic_exog_coefs) + for i in eachindex(model.dynamic_exog_coefs) + push!(model_innovations, "dynamic_exog_coefs_$i") + end + end return model_innovations end @@ -1083,7 +1139,11 @@ function get_trend_decomposition( end if model.stochastic_level && !isempty(components["ξ"]["Coefs"]) - ξ = vcat(zeros(max(2, model.stochastic_start) - 1), components["ξ"]["Coefs"], 0.0) + ξ = vcat( + zeros(max(2, model.stochastic_start) - 1), + components["ξ"]["Coefs"], + zeros(model.ξ_threshold), + ) @assert length(ξ) == T else ξ = zeros(AbstractFloat, T) @@ -1572,7 +1632,7 @@ function forecast_dynamic_exog_coefs( dynamic_exog_coefs = Vector{Tuple}(undef, length(model.dynamic_exog_coefs)) for i in eachindex(model.dynamic_exog_coefs) if model.dynamic_exog_coefs[i][2] == "level" - n_coefs = 1 + ξ_size(T, model.stochastic_start) + n_coefs = 1 + ξ_size(T, model.ξ_threshold, model.stochastic_start) extra_param = "" elseif model.dynamic_exog_coefs[i][2] == "slope" n_coefs = 1 + ζ_size(T, model.ζ_threshold, model.stochastic_start) @@ -1601,6 +1661,7 @@ function forecast_dynamic_exog_coefs( dynamic_exog_coefs, T, steps_ahead, + model.ξ_threshold, model.ζ_threshold, model.ω_threshold, model.ϕ_threshold, diff --git a/test/models/structural_model.jl b/test/models/structural_model.jl index b2b226d..f6217f9 100644 --- a/test/models/structural_model.jl +++ b/test/models/structural_model.jl @@ -5,11 +5,15 @@ model2 = StateSpaceLearning.StructuralModel(y1; freq_seasonal=[3, 10]) model3 = StateSpaceLearning.StructuralModel(y1; cycle_period=[3, 10.2]) model4 = StateSpaceLearning.StructuralModel(y1; cycle_period=[3, 10.2]) + model5 = StateSpaceLearning.BasicStructuralModel(y1; slope="none") @test typeof(model1) == StateSpaceLearning.StructuralModel @test typeof(model2) == StateSpaceLearning.StructuralModel @test typeof(model3) == StateSpaceLearning.StructuralModel @test typeof(model4) == StateSpaceLearning.StructuralModel + @test typeof(model5) == StateSpaceLearning.StructuralModel + + @test size(model5.X) == (100, 201) @test_throws AssertionError StateSpaceLearning.StructuralModel(y1; stochastic_start=0) @test_throws AssertionError StateSpaceLearning.StructuralModel( @@ -43,7 +47,7 @@ end end @testset "Innovation matrices" begin - @test StateSpaceLearning.ξ_size(10, 1) == 8 + @test StateSpaceLearning.ξ_size(10, 1, 1) == 8 @test StateSpaceLearning.ζ_size(10, 2, 1) == 6 @test StateSpaceLearning.ζ_size(10, 0, 1) == 8 @test StateSpaceLearning.ω_size(10, 2, 0, 1) == 9 @@ -52,7 +56,7 @@ end @test StateSpaceLearning.ϕ_size(10, 0, 1) == 16 @test StateSpaceLearning.ϕ_size(10, 2, 1) == 14 - @test StateSpaceLearning.ξ_size(10, 5) == 5 + @test StateSpaceLearning.ξ_size(10, 1, 5) == 5 @test StateSpaceLearning.ζ_size(10, 2, 5) == 3 @test StateSpaceLearning.ζ_size(10, 0, 5) == 5 @test StateSpaceLearning.ω_size(10, 2, 0, 5) == 6 @@ -60,8 +64,8 @@ end @test StateSpaceLearning.o_size(10, 6) == 5 @test StateSpaceLearning.ϕ_size(10, 0, 5) == 10 - X_ξ1 = StateSpaceLearning.create_ξ(5, 1) - X_ξ2 = StateSpaceLearning.create_ξ(5, 3) + X_ξ1 = StateSpaceLearning.create_ξ(5, 1, 1) + X_ξ2 = StateSpaceLearning.create_ξ(5, 1, 3) @test X_ξ1 == [ 0.0 0.0 0.0 @@ -173,14 +177,14 @@ end dynamic_exog_coefs1 = [(collect(1:5), "level")] X1 = StateSpaceLearning.create_dynamic_exog_coefs_matrix( - dynamic_exog_coefs1, 5, 0, 0, 0, 1 + dynamic_exog_coefs1, 5, 1, 0, 0, 0, 1 ) @test size(X1) == (5, 4) dynamic_exog_coefs2 = [(collect(1:5), "slope")] X2 = StateSpaceLearning.create_dynamic_exog_coefs_matrix( - dynamic_exog_coefs2, 5, 0, 0, 0, 1 + dynamic_exog_coefs2, 5, 1, 0, 0, 0, 1 ) @test size(X2) == (5, 4) @@ -192,7 +196,7 @@ end ] X = StateSpaceLearning.create_dynamic_exog_coefs_matrix( - dynamic_exog_coefs, 5, 0, 0, 0, 1 + dynamic_exog_coefs, 5, 1, 0, 0, 0, 1 ) @test size(X) == (5, 22) @@ -203,21 +207,21 @@ end (collect(6:7), "cycle", 3, 8), ] X_f = StateSpaceLearning.create_forecast_dynamic_exog_coefs_matrix( - dynamic_exog_coefs, 5, 2, 0, 0, 0, 1 + dynamic_exog_coefs, 5, 2, 1, 0, 0, 0, 1 ) @test size(X_f) == (2, 23) dynamic_exog_coefs2 = [(collect(6:7), "level", "", 4)] X_f2 = StateSpaceLearning.create_forecast_dynamic_exog_coefs_matrix( - dynamic_exog_coefs2, 5, 2, 0, 0, 0, 1 + dynamic_exog_coefs2, 5, 2, 1, 0, 0, 0, 1 ) @test size(X_f2) == (2, 4) dynamic_exog_coefs3 = [(collect(6:7), "slope", "", 4)] X_f3 = StateSpaceLearning.create_forecast_dynamic_exog_coefs_matrix( - dynamic_exog_coefs3, 5, 2, 0, 0, 0, 1 + dynamic_exog_coefs3, 5, 2, 1, 0, 0, 0, 1 ) @test size(X_f3) == (2, 4) end @@ -243,6 +247,7 @@ end 3, 3, true, + 1, 0, 0, 0, @@ -271,6 +276,7 @@ end 3, 3, true, + 1, 3, 2, 4,