diff --git a/Project.toml b/Project.toml index 7aaf3dd..80f3404 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StateSpaceLearning" uuid = "971c4b7c-2c4e-4bac-8525-e842df3cde7b" authors = ["andreramosfc "] -version = "1.4.3" +version = "2.0.0" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" diff --git a/README.md b/README.md index 8388d9b..e06d4bb 100644 --- a/README.md +++ b/README.md @@ -20,35 +20,41 @@ model = StructuralModel(y) fit!(model) # Point Forecast -prediction = StateSpaceLearning.forecast(model, 12) #Gets a 12 steps ahead prediction +prediction = forecast(model, 12) # Gets a 12 steps ahead prediction # Scenarios Path Simulation -simulation = StateSpaceLearning.simulate(model, 12, 1000) #Gets 1000 scenarios path of 12 steps ahead predictions +simulation = simulate(model, 12, 1000) # Gets 1000 scenarios path of 12 steps ahead predictions ``` ## StructuralModel Arguments -* `y::Vector`: Vector of data. -* `level::Bool`: Boolean where to consider intercept in the model (default: true) -* `stochastic_level::Bool`: Boolean where to consider stochastic level component in the model (default: true) -* `trend::Bool`: Boolean where to consider trend component in the model (default: true) -* `stochastic_trend::Bool`: Boolean where to consider stochastic trend component in the model (default: true) -* `seasonal::Bool`: Boolean where to consider seasonal component in the model (default: true) -* `stochastic_seasonal::Bool`: Boolean where to consider stochastic seasonal component in the model (default: true) -* `freq_seasonal::Int`: Seasonal frequency to be considered in the model (default: 12) -* `outlier::Bool`: Boolean where to consider outlier component in the model (default: true) -* `ζ_ω_threshold::Int`: Argument to stabilize `stochastic trend` and `stochastic seasonal` components (default: 12) +* `y::Union{Vector,Matrix}`: Time series data +* `level::String`: Level component type: "stochastic", "deterministic", or "none" (default: "stochastic") +* `slope::String`: Slope component type: "stochastic", "deterministic", or "none" (default: "stochastic") +* `seasonal::String`: Seasonal component type: "stochastic", "deterministic", or "none" (default: "stochastic") +* `cycle::String`: Cycle component type: "stochastic", "deterministic", or "none" (default: "none") +* `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 slope innovations (default: 12) +* `ω_threshold::Int`: Threshold for seasonal innovations (default: 12) +* `ϕ_threshold::Int`: Threshold for cycle innovations (default: 12) +* `stochastic_start::Int`: Starting point for stochastic components (default: 1) +* `exog::Matrix`: Matrix of exogenous variables (default: zeros(length(y), 0)) +* `dynamic_exog_coefs::Union{Vector{<:Tuple}, Nothing}`: Dynamic exogenous coefficients (default: nothing) ## Features Current features include: -* Estimation -* Components decomposition -* Forecasting -* Completion of missing values -* Predefined models, including: -* Outlier detection -* Outlier robust models +* Model estimation using elastic net based regularization +* Automatic component decomposition (trend, seasonal, cycle) +* Point forecasting and scenario simulation +* Missing value imputation +* Outlier detection and robust modeling +* Multiple seasonal frequencies support +* Deterministic and stochastic component options +* Dynamic exogenous variable handling +* Best subset selection for exogenous variables ## Quick Examples @@ -67,7 +73,7 @@ steps_ahead = 30 model = StructuralModel(log_air_passengers) fit!(model) -prediction_log = StateSpaceLearning.forecast(model, steps_ahead) # arguments are the output of the fitted model and number of steps ahead the user wants to forecast +prediction_log = forecast(model, steps_ahead) # arguments are the output of the fitted model and number of steps ahead the user wants to forecast prediction = exp.(prediction_log) plot_point_forecast(airp.passengers, prediction) @@ -76,7 +82,7 @@ plot_point_forecast(airp.passengers, prediction) ```julia N_scenarios = 1000 -simulation = StateSpaceLearning.simulate(model, steps_ahead, N_scenarios) # arguments are the output of the fitted model, number of steps ahead the user wants to forecast and number of scenario paths +simulation = simulate(model, steps_ahead, N_scenarios) # arguments are the output of the fitted model, number of steps ahead the user wants to forecast and number of scenario paths plot_scenarios(airp.passengers, exp.(simulation)) @@ -97,22 +103,20 @@ log_air_passengers = log.(airp.passengers) model = StructuralModel(log_air_passengers) fit!(model) -level = model.output.components["μ1"]["Values"] + model.output.components["ξ"]["Values"] -slope = model.output.components["ν1"]["Values"] + model.output.components["ζ"]["Values"] -seasonal = model.output.components["γ1_12"]["Values"] + model.output.components["ω_12"]["Values"] -trend = level + slope - -plot(trend, w=2 , color = "Black", lab = "Trend Component", legend = :outerbottom) -plot(seasonal, w=2 , color = "Black", lab = "Seasonal Component", legend = :outerbottom) +# Access decomposed components directly +trend = model.output.decomposition["trend"] +seasonal = model.output.decomposition["seasonal_12"] +plot(trend, w=2, color = "Black", lab = "Trend Component", legend = :outerbottom) +plot(seasonal, w=2, color = "Black", lab = "Seasonal Component", legend = :outerbottom) ``` | ![quick_example_trend](./docs/src/assets/trend.svg) | ![quick_example_seas](./docs/src/assets/seasonal.svg)| |:------------------------------:|:-----------------------------:| -### Best Subset Selection -Quick example on how to perform best subset selection in time series utilizing StateSpaceLearning. +### Best Subset Selection and Dynamic Coefficients +Example of performing best subset selection and using dynamic coefficients: ```julia using StateSpaceLearning @@ -122,22 +126,33 @@ using Random Random.seed!(2024) +# Load data airp = CSV.File(StateSpaceLearning.AIR_PASSENGERS) |> DataFrame log_air_passengers = log.(airp.passengers) -X = rand(length(log_air_passengers), 10) # Create 10 exogenous features -β = rand(3) - -y = log_air_passengers + X[:, 1:3]*β # add to the log_air_passengers series a contribution from only 3 exogenous features. - -model = StructuralModel(y; Exogenous_X = X) -fit!(model; α = 1.0, information_criteria = "bic", ϵ = 0.05, penalize_exogenous = true, penalize_initial_states = true) - -Selected_exogenous = model.output.components["Exogenous_X"]["Selected"] +# Create exogenous features +X = rand(length(log_air_passengers), 10) +β = rand(3) +y = log_air_passengers + X[:, 1:3]*β + +# Create model with exogenous variables +model = StructuralModel(y; + exog = X +) + +# Fit model with elastic net regularization +fit!(model; + α = 1.0, # 1.0 for Lasso, 0.0 for Ridge + information_criteria = "bic", + ϵ = 0.05, + penalize_exogenous = true, + penalize_initial_states = true +) + +# Get selected features +selected_exog = model.output.components["exog"]["Selected"] ``` -In this example, the selected exogenous features were 1, 2, 3, as expected. - ### Missing values imputation Quick example of completion of missing values for the air passengers time-series (artificial NaN values are added to the original time-series). @@ -209,7 +224,7 @@ fit!(model) residuals_variances = model.output.residuals_variances ss_model = BasicStructural(log_air_passengers, 12) -set_initial_hyperparameters!(ss_model, Dict("sigma2_ε" => residuals_variances["ε"], +StateSpaceModels.set_initial_hyperparameters!(ss_model, Dict("sigma2_ε" => residuals_variances["ε"], "sigma2_ξ" =>residuals_variances["ξ"], "sigma2_ζ" =>residuals_variances["ζ"], "sigma2_ω" =>residuals_variances["ω_12"])) diff --git a/docs/src/assets/dynamic_exog.png b/docs/src/assets/dynamic_exog.png new file mode 100644 index 0000000..bbab8cd Binary files /dev/null and b/docs/src/assets/dynamic_exog.png differ diff --git a/docs/src/examples.md b/docs/src/examples.md index d0f2d3f..57ed198 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -46,7 +46,7 @@ y = solars[!, "y1"] T = length(y) steps_ahead = 48 -model = StructuralModel(y; freq_seasonal=24, trend=false, level=false) +model = StructuralModel(y; freq_seasonal=24, slope="none", level="none", outlier=false) fit!(model; penalize_initial_states=false) simulation = StateSpaceLearning.simulate(model, steps_ahead, 100) #Gets a 12 steps ahead prediction plot_scenarios(y, simulation) @@ -63,7 +63,7 @@ y = solars[!, "y1"] T = length(y) steps_ahead = 48 -model = StructuralModel(y; freq_seasonal=24, trend=false, level=false) +model = StructuralModel(y; freq_seasonal=24, slope="none", level="none", outlier = false) fit!(model; penalize_initial_states=false) simulation = StateSpaceLearning.simulate(model, steps_ahead, 100; seasonal_innovation_simulation=24) #Gets a 12 steps ahead prediction plot_scenarios(y, simulation) @@ -114,4 +114,35 @@ plot_point_forecast(y, prediction) ``` ![two_seas](assets/two_seas.png) -Note that the model was able to capture both seasonalities in this case. \ No newline at end of file +Note that the model was able to capture both seasonalities in this case. + +## Dynamic Exogenous Coefficients + +Dynamic exogenous coefficients allow the effect of exogenous variables to vary over time with specific patterns (e.g., level, slope, seasonal or cyclical). This is configured through the `dynamic_exog_coefs` parameter in the StructuralModel constructor. + +The `dynamic_exog_coefs` parameter accepts a vector of tuples, where each tuple contains: +- First element: A vector of an exogenous variable +- Second element: The name of the component that the exogenous variable will be associated with +- Third element (optional): For the seasonal component, the freq_seasonal parameter and for cycle component, the cycle_period parameter. + +For example: +```julia +# Make X1's effect vary annually and X2's effect vary semi-annually +airp = CSV.File(StateSpaceLearning.AIR_PASSENGERS) |> DataFrame +y = log.(airp.passengers) +X = vcat(collect(1:90), collect(90.5:-0.5:64)) + (rand(144) .* 10) +y += X * -0.03 + +dynamic_coefs = [(X, "level")] + +model = StructuralModel(y; + dynamic_exog_coefs=dynamic_coefs +) + +fit!(model) + +prediction = forecast(model, 30; dynamic_exog_coefs_forecasts = [collect(63.5:-0.5:49)]) + +plot_point_forecast(y, prediction) +``` +![dynamic_exog](assets/dynamic_exog.png) diff --git a/docs/src/features.md b/docs/src/features.md index 372cb0a..edd1e2c 100644 --- a/docs/src/features.md +++ b/docs/src/features.md @@ -104,10 +104,10 @@ X = rand(length(log_air_passengers), 10) # Create 10 exogenous features y = log_air_passengers + X[:, 1:3]*β # add to the log_air_passengers series a contribution from only 3 exogenous features. -model = StructuralModel(y; Exogenous_X = X) +model = StructuralModel(y; exog = X) fit!(model; α = 1.0, information_criteria = "bic", ϵ = 0.05, penalize_exogenous = true, penalize_initial_states = true) -Selected_exogenous = model.output.components["Exogenous_X"]["Selected"] +Selected_exogenous = model.output.components["exog"]["Selected"] ``` diff --git a/docs/src/manual.md b/docs/src/manual.md index e8d20e9..ce4f096 100644 --- a/docs/src/manual.md +++ b/docs/src/manual.md @@ -23,7 +23,14 @@ simulation = StateSpaceLearning.simulate(model, 12, 1000) #Gets 1000 scenarios p ``` ## Models -The package currently supports the implementation of the StructuralModel. If you have suggestions for additional models to include, we encourage you to contribute by opening an issue or submitting a pull request. +The package currently supports the implementation of the StructuralModel, which includes capabilities for handling dynamic exogenous coefficients. If you have suggestions for additional models to include, we encourage you to contribute by opening an issue or submitting a pull request. + +``` + +When using dynamic coefficients: +- The model will create time-varying coefficients for each specified exogenous variable +- Each coefficient will follow the specified cyclical pattern +- When forecasting, you must provide future values for exogenous variables using the `Exogenous_Forecast` parameter ```@docs StateSpaceLearning.StructuralModel @@ -39,7 +46,11 @@ StateSpaceLearning.fit! ## Forecasting and Simulating -The package has functions to make point forecasts multiple steps ahead and to simulate scenarios based on those forecasts. These functions are implemented both for the univariate and to the multivariate cases. +The package has functions to make point forecasts multiple steps ahead and to simulate scenarios based on those forecasts. These functions are implemented for both univariate and multivariate cases, with support for exogenous variables and dynamic coefficients. + +When using models with exogenous variables: +- For standard exogenous variables, provide future values using the `Exogenous_Forecast` parameter +- For dynamic coefficients, use the same `Exogenous_Forecast` parameter with values for each exogenous variable ```@docs StateSpaceLearning.forecast diff --git a/paper_tests/m4_test/evaluate_model.jl b/paper_tests/m4_test/evaluate_model.jl index 852da58..987cab7 100644 --- a/paper_tests/m4_test/evaluate_model.jl +++ b/paper_tests/m4_test/evaluate_model.jl @@ -19,15 +19,13 @@ function evaluate_SSL( model = StateSpaceLearning.StructuralModel( normalized_y; - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=true, + level="stochastic", + slope="stochastic", + seasonal="stochastic", freq_seasonal=12, outlier=outlier, - ζ_ω_threshold=12, + ζ_threshold=12, + ω_threshold=12, ) StateSpaceLearning.fit!( model; diff --git a/paper_tests/m4_test/m4_test.jl b/paper_tests/m4_test/m4_test.jl index 7b18025..9c9e723 100644 --- a/paper_tests/m4_test/m4_test.jl +++ b/paper_tests/m4_test/m4_test.jl @@ -50,7 +50,7 @@ function run_config( save_init ? CSV.write(init_filepath, initialization_df) : nothing # Initialize empty CSV for i in 1:48000 - if i in [10001, 20001, 30001, 40001] # Clear DataFrame to save memory + if i % 1000 == 1 # Clear DataFrame to save memory results_df = DataFrame() initialization_df = DataFrame() end @@ -66,7 +66,8 @@ function run_config( information_criteria, ) - if i in [10000, 20000, 30000, 40000, 48000] + if i % 1000 == 0 + @info "Saving results for $i series" !save_init ? append_results(filepath, results_df) : nothing save_init ? append_results(init_filepath, initialization_df) : nothing end @@ -100,9 +101,9 @@ end # Main script function main() results_table = DataFrame() - for outlier in [true] - for information_criteria in ["aic"] - for α in [0.1] + for outlier in [true, false] + for information_criteria in ["aic", "bic"] + for α in [0.0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0] @info "Running SSL with outlier = $outlier, information_criteria = $information_criteria, α = $α" results_table = run_config( results_table, outlier, information_criteria, α, false, 60 diff --git a/paper_tests/m4_test/prepare_data.jl b/paper_tests/m4_test/prepare_data.jl index 9d69a57..4ea1be0 100644 --- a/paper_tests/m4_test/prepare_data.jl +++ b/paper_tests/m4_test/prepare_data.jl @@ -2,7 +2,7 @@ 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::AbstractFloat, min_y::AbstractFloat) +function de_normalize(y, 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 c1d029d..72eddf8 100644 --- a/paper_tests/simulation_test/evaluate_models.jl +++ b/paper_tests/simulation_test/evaluate_models.jl @@ -28,7 +28,7 @@ function get_SSL_results( freq_seasonal=12, outlier=false, ζ_ω_threshold=12, - Exogenous_X=X_train, + exog=X_train, ) t = @elapsed StateSpaceLearning.fit!( model; @@ -39,13 +39,13 @@ function get_SSL_results( penalize_initial_states=true, ) - selected = model.output.components["Exogenous_X"]["Selected"] + selected = model.output.components["exog"]["Selected"] true_positives, false_positives, false_negatives, true_negatives = get_confusion_matrix( selected, true_features, false_features ) - mse = mse_func(model.output.components["Exogenous_X"]["Coefs"], true_β) - bias = bias_func(model.output.components["Exogenous_X"]["Coefs"], true_β) + mse = mse_func(model.output.components["exog"]["Coefs"], true_β) + bias = bias_func(model.output.components["exog"]["Coefs"], true_β) series_result = DataFrame( [ diff --git a/src/StateSpaceLearning.jl b/src/StateSpaceLearning.jl index d19d07e..e3dfb01 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 @@ -13,6 +13,12 @@ include("datasets.jl") include("fit_forecast.jl") include("plots.jl") -export fit!, forecast, simulate, StructuralModel, plot_point_forecast, plot_scenarios +export fit!, + forecast, + simulate, + StructuralModel, + plot_point_forecast, + plot_scenarios, + simulate_states end # module StateSpaceLearning diff --git a/src/estimation_procedure.jl b/src/estimation_procedure.jl index d3548da..e0bad0f 100644 --- a/src/estimation_procedure.jl +++ b/src/estimation_procedure.jl @@ -1,22 +1,22 @@ """ - get_dummy_indexes(Exogenous_X::Matrix{Fl}) where {Fl} + get_dummy_indexes(exog::Matrix{Fl}) where {Fl} Identifies and returns the indexes of columns in the exogenous matrix that contain dummy variables. # Arguments - - `Exogenous_X::Matrix{Fl}`: Exogenous variables matrix. + - `exog::Matrix{Fl}`: Exogenous variables matrix. # Returns - `Vector{Int}`: Vector containing the indexes of columns with dummy variables. """ -function get_dummy_indexes(Exogenous_X::Matrix{Fl}) where {Fl} - T, p = size(Exogenous_X) +function get_dummy_indexes(exog::Matrix{Fl}) where {Fl} + T, p = size(exog) dummy_indexes = [] for i in 1:p - if count(iszero.(Exogenous_X[:, i])) == T - 1 - push!(dummy_indexes, findfirst(i -> i != 0.0, Exogenous_X[:, i])) + if count(iszero.(exog[:, i])) == T - 1 + push!(dummy_indexes, findfirst(i -> i != 0.0, exog[:, i])) end end @@ -43,7 +43,7 @@ function get_outlier_duplicate_columns( return [] else o_indexes = components_indexes["o"] - exogenous_indexes = components_indexes["Exogenous_X"] + exogenous_indexes = components_indexes["exog"] dummy_indexes = get_dummy_indexes(Estimation_X[:, exogenous_indexes]) @@ -203,14 +203,14 @@ function fit_lasso( hasintercept = has_intercept(Estimation_X) if hasintercept if !penalize_exogenous - penalty_factor[components_indexes["Exogenous_X"] .- 1] .= 0 + penalty_factor[components_indexes["exog"] .- 1] .= 0 else nothing end Lasso_X = Estimation_X[:, 2:end] else if !penalize_exogenous - penalty_factor[components_indexes["Exogenous_X"]] .= 0 + penalty_factor[components_indexes["exog"]] .= 0 else nothing end @@ -257,6 +257,8 @@ end ϵ::AbstractFloat, penalize_exogenous::Bool, penalize_initial_states::Bool, + innovations_names::Vector{String}, + initial_states_names::Vector{String}, )::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. @@ -270,6 +272,9 @@ end - `ϵ::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. + - `innovations_names::Vector{String}`: Vector of strings containing the names of the innovations. + - `initial_states_names::Vector{String}`: Vector of strings containing the names of the initial states. + # Returns - `Tuple{Vector{AbstractFloat}, Vector{AbstractFloat}}`: Tuple containing coefficients and residuals of the fitted AdaLasso model. @@ -284,6 +289,7 @@ function estimation_procedure( ϵ::AbstractFloat, penalize_exogenous::Bool, penalize_initial_states::Bool, + innovations_names::Vector{String}, )::Tuple{ Vector{AbstractFloat},Vector{AbstractFloat} } where {Fl<:AbstractFloat,Tl<:AbstractFloat} @@ -292,11 +298,15 @@ function estimation_procedure( hasintercept = has_intercept(Estimation_X) + # all zero columns in X + all_zero_idx = findall(i -> all(iszero, Estimation_X[:, i]), 1:size(Estimation_X, 2)) + if hasintercept penalty_factor = ones(size(Estimation_X, 2) - 1) if length(penalty_factor) != length(components_indexes["initial_states"][2:end]) penalty_factor[components_indexes["initial_states"][2:end] .- 1] .= 0 end + penalty_factor[all_zero_idx .- 1] .= Inf coefs, _ = fit_lasso( Estimation_X, estimation_y, @@ -312,6 +322,7 @@ function estimation_procedure( if length(penalty_factor) != length(components_indexes["initial_states"]) penalty_factor[components_indexes["initial_states"][1:end]] .= 0 end + penalty_factor[all_zero_idx] .= Inf coefs, _ = fit_lasso( Estimation_X, estimation_y, @@ -330,10 +341,7 @@ function estimation_procedure( for key in keys(components_indexes) if key != "initial_states" && key != "μ1" component = components_indexes[key] - if key != "Exogenous_X" && - key != "o" && - !(key in ["ν1"]) && - !(occursin("γ", key)) + if key in innovations_names κ = count(i -> i != 0, coefs[component]) < 1 ? 0 : std(coefs[component]) if hasintercept ts_penalty_factor[component .- 1] .= (1 / (κ + ϵ)) @@ -356,12 +364,14 @@ function estimation_procedure( else nothing end + ts_penalty_factor[all_zero_idx .- 1] .= Inf else if !penalize_initial_states ts_penalty_factor[components_indexes["initial_states"][1:end]] .= 0 else nothing end + ts_penalty_factor[all_zero_idx] .= Inf end return fit_lasso( @@ -375,64 +385,3 @@ 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 diff --git a/src/fit_forecast.jl b/src/fit_forecast.jl index 5836d7f..cfed4a1 100644 --- a/src/fit_forecast.jl +++ b/src/fit_forecast.jl @@ -44,6 +44,8 @@ function fit!( components_indexes = get_components_indexes(model) + innovations_names = get_model_innovations(model) + coefs, estimation_ε = estimation_procedure( Estimation_X, estimation_y, @@ -53,6 +55,7 @@ function fit!( ϵ, penalize_exogenous, penalize_initial_states, + innovations_names, ) components = build_components(Estimation_X, coefs, components_indexes) @@ -63,292 +66,9 @@ function fit!( ε, fitted = get_fit_and_residuals(estimation_ε, coefs, model.X, valid_indexes, T) - 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 - -@doc raw""" -Returns the forecast for a given number of steps ahead using the provided StateSpaceLearning output and exogenous forecast data. - -forecast(model::StateSpaceLearningModel, steps\_ahead::Int; Exogenous\_Forecast::Union{Matrix{Fl}, Missing}=missing)::Vector{AbstractFloat} where Fl - -# Arguments -- `model::StateSpaceLearningModel`: Model obtained from fitting. -- `steps_ahead::Int`: Number of steps ahead for forecasting. -- `Exogenous_Forecast::Matrix{Fl}`: Exogenous variables forecast (default: zeros(steps_ahead, 0)) - -# Returns -- `Union{Matrix{AbstractFloat}, Vector{AbstractFloat}}`: Matrix or vector of matrices containing forecasted values. - -# Example -```julia -y = rand(100) -model = StructuralModel(y) -fit!(model) -steps_ahead = 12 -point_prediction = forecast(model, steps_ahead) -``` -""" -function forecast( - model::StateSpaceLearningModel, - steps_ahead::Int; - Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, 0), -)::Union{Matrix{<:AbstractFloat},Vector{<:AbstractFloat}} where {Fl<:AbstractFloat} - @assert isfitted(model) "Model must be fitted before simulation" - exog_idx = if typeof(model.output) == Output - model.output.components["Exogenous_X"]["Indexes"] - else - model.output[1].components["Exogenous_X"]["Indexes"] - end - @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[:, exog_idx] - complete_matrix = create_X(model, Exogenous_X, steps_ahead, Exogenous_Forecast) - - 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 - -@doc raw""" -Generate simulations for a given number of steps ahead using the provided StateSpaceLearning output and exogenous forecast data. - -simulate(model::StateSpaceLearningModel, steps\_ahead::Int, N\_scenarios::Int; - Exogenous\_Forecast::Matrix{Fl}=zeros(steps_ahead, 0))::Matrix{AbstractFloat} where Fl - -# Arguments -- `model::StateSpaceLearningModel`: Model obtained from fitting. -- `steps_ahead::Int`: Number of steps ahead for simulation. -- `N_scenarios::Int`: Number of scenarios to simulate (default: 1000). -- `Exogenous_Forecast::Matrix{Fl}`: Exogenous variables forecast (default: zeros(steps_ahead, 0)) - -# Returns -- `Union{Vector{Matrix{AbstractFloat}}, Matrix{AbstractFloat}}`: Matrix or vector of matrices containing simulated values. - -# Example (Univariate Case) -```julia -y = rand(100) -model = StructuralModel(y) -fit!(model) -steps_ahead = 12 -N_scenarios = 1000 -simulation = simulate(model, steps_ahead, N_scenarios) -``` + decomposition = get_model_decomposition(model, components) -# Example (Multivariate Case) -```julia -y = rand(100, 3) -model = StructuralModel(y) -fit!(model) -steps_ahead = 12 -N_scenarios = 1000 -simulation = simulate(model, steps_ahead, N_scenarios) -``` -""" -function simulate( - model::StateSpaceLearningModel, - steps_ahead::Int, - N_scenarios::Int; - Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, 0), - 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" - @assert isfitted(model) "Model must be fitted before simulation" - - prediction = forecast(model, steps_ahead; Exogenous_Forecast=Exogenous_Forecast) - - is_univariate = typeof(model.output) == Output - - simulation_X = zeros(steps_ahead, 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) - for innovation in model_innovations - 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, valid_indexes) - components_matrix = hcat(components_matrix, comp) - N_components += 1 - end - - 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 - 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 - - MV_dist_vec = Vector{MvNormal}(undef, steps_ahead) - o_noises = if is_univariate - zeros(steps_ahead, N_scenarios) - else - [zeros(steps_ahead, N_scenarios) for _ in 1:length(model.output)] - end - - if seasonal_innovation_simulation == 0 - ∑ = if is_univariate - Diagonal([var(components_matrix[:, i]) for i in 1:N_components]) - else - Diagonal([var(components_matrix[:, i]) for i in 1:N_mv_components]) - end - for i in 1:steps_ahead - MV_dist_vec[i] = if is_univariate - MvNormal(zeros(N_components), ∑) - else - MvNormal(zeros(N_mv_components), ∑) - end - end - - 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:min(seasonal_innovation_simulation, steps_ahead) - ∑ = if is_univariate - Diagonal([ - var( - components_matrix[ - (i + start_seasonal_term):seasonal_innovation_simulation:end, - j, - ], - ) for j in 1:N_components - ]) - else - Diagonal([ - var( - components_matrix[ - (i + start_seasonal_term):seasonal_innovation_simulation:end, - j, - ], - ) for j in 1:N_mv_components - ]) - end - - MV_dist_vec[i] = if is_univariate - MvNormal(zeros(N_components), ∑) - else - MvNormal(zeros(N_mv_components), ∑) - end - if is_univariate - if model.outlier - o_noises[i, :] = rand( - Normal( - 0, - std( - model.output.components["o"]["Coefs"][(i + start_seasonal_term):seasonal_innovation_simulation:end], - ), - ), - N_scenarios, - ) - else - nothing - end - else - for j in eachindex(model.output) - if 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, - ) - else - nothing - end - end - end - end - for i in (seasonal_innovation_simulation + 1):steps_ahead - MV_dist_vec[i] = MV_dist_vec[i - seasonal_innovation_simulation] - if model.outlier - if is_univariate - o_noises[i, :] = o_noises[i - seasonal_innovation_simulation, :] - else - for j in eachindex(model.output) - o_noises[j][i, :] = o_noises[j][ - i - seasonal_innovation_simulation, :, - ] - end - end - end - end - end - - simulation = if is_univariate - AbstractFloat.(hcat([prediction for _ in 1:N_scenarios]...)) - else - [ - AbstractFloat.(hcat([prediction[:, i] for _ in 1:N_scenarios]...)) for - i in eachindex(model.output) - ] - end - 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 + return model.output = Output( + coefs, ε, fitted, residuals_variances, valid_indexes, components, decomposition + ) end diff --git a/src/models/structural_model.jl b/src/models/structural_model.jl index 1f79683..c973b36 100644 --- a/src/models/structural_model.jl +++ b/src/models/structural_model.jl @@ -10,14 +10,15 @@ Instantiates a Structural State Space Learning model. seasonal::Bool=true, stochastic_seasonal::Bool=true, freq_seasonal::Union{Int, Vector{Int}}=12, + cycle_period::Union{Union{Int,<:AbstractFloat},Vector{Int},Vector{<:AbstractFloat}}=0, + stochastic_cycle::Bool=false, outlier::Bool=true, - ζ_ω_threshold::Int=12, + ζ_threshold::Int=12, + ω_threshold::Int=12, + ϕ_threshold::Int=12, stochastic_start::Int=1, - Exogenous_X::Matrix=if typeof(y) <: Vector - zeros(length(y), 0) - else - zeros(size(y, 1), 0) - end, + exog::Matrix=zeros(length(y), 0), + dynamic_exog_coefs::Union{Vector{<:Tuple}, Nothing}=nothing ) A Structural State Space Learning model that can have level, stochastic_level, trend, stochastic_trend, seasonal, stochastic_seasonal, outlier and Exogenous components. Each component should be specified by Booleans. @@ -55,118 +56,136 @@ mutable struct StructuralModel <: StateSpaceLearningModel X::Matrix level::Bool stochastic_level::Bool - trend::Bool - stochastic_trend::Bool + slope::Bool + stochastic_slope::Bool seasonal::Bool stochastic_seasonal::Bool + cycle::Bool + stochastic_cycle::Bool freq_seasonal::Union{Int,Vector{Int}} cycle_period::Union{Union{Int,<:AbstractFloat},Vector{Int},Vector{<:AbstractFloat}} - cycle_matrix::Vector{Matrix} - stochastic_cycle::Bool outlier::Bool - ζ_ω_threshold::Int + ζ_threshold::Int + ω_threshold::Int + ϕ_threshold::Int stochastic_start::Int n_exogenous::Int + dynamic_exog_coefs::Union{Vector{<:Tuple},Nothing} output::Union{Vector{Output},Output,Nothing} function StructuralModel( y::Union{Vector,Matrix}; - level::Bool=true, - stochastic_level::Bool=true, - trend::Bool=true, - stochastic_trend::Bool=true, - seasonal::Bool=true, - stochastic_seasonal::Bool=true, + level::String="stochastic", + slope::String="stochastic", + seasonal::String="stochastic", + cycle::String="none", freq_seasonal::Union{Int,Vector{Int}}=12, cycle_period::Union{Union{Int,<:AbstractFloat},Vector{Int},Vector{<:AbstractFloat}}=0, - dumping_cycle::Float64=1.0, - stochastic_cycle::Bool=false, outlier::Bool=true, - ζ_ω_threshold::Int=12, + ζ_threshold::Int=12, + ω_threshold::Int=12, + ϕ_threshold::Int=12, stochastic_start::Int=1, - Exogenous_X::Matrix=if typeof(y) <: Vector - zeros(length(y), 0) - else - zeros(size(y, 1), 0) - end, + exog::Matrix=zeros(length(y), 0), + dynamic_exog_coefs::Union{Vector{<:Tuple},Nothing}=nothing, ) - n_exogenous = size(Exogenous_X, 2) - @assert !has_intercept(Exogenous_X) "Exogenous matrix must not have an intercept column" - if typeof(y) <: Vector - @assert seasonal ? length(y) > minimum(freq_seasonal) : true "Time series must be longer than the seasonal period" - else - @assert seasonal ? size(y, 1) > minimum(freq_seasonal) : true "Time series must be longer than the seasonal period" - end + n_exogenous = size(exog, 2) + + @assert !has_intercept(exog) "Exogenous matrix must not have an intercept column" @assert 1 <= stochastic_start < length(y) "stochastic_start must be greater than or equal to 1 and smaller than the length of the time series" - @assert 0 < dumping_cycle <= 1 "dumping_cycle must be greater than 0 and less than or equal to 1" - if cycle_period != 0 && !isempty(cycle_period) - if typeof(cycle_period) <: Vector - cycle_matrix = Vector{Matrix}(undef, length(cycle_period)) - for i in eachindex(cycle_period) - A = dumping_cycle * cos(2 * pi / cycle_period[i]) - B = dumping_cycle * sin(2 * pi / cycle_period[i]) - cycle_matrix[i] = [A B; -B A] - end - else - cycle_matrix = Vector{Matrix}(undef, 1) - A = dumping_cycle * cos(2 * pi / cycle_period) - B = dumping_cycle * sin(2 * pi / cycle_period) - cycle_matrix[1] = [A B; -B A] - end - else - cycle_matrix = Vector{Matrix}(undef, 0) - end + @assert level in ["deterministic", "stochastic", "none"] "level must be either deterministic, stochastic or no" + @assert slope in ["deterministic", "stochastic", "none"] "slope must be either deterministic, stochastic or no" + @assert seasonal in ["deterministic", "stochastic", "none"] "seasonal must be either deterministic, stochastic or no" + @assert cycle in ["deterministic", "stochastic", "none"] "cycle must be either deterministic, stochastic or no" + @assert seasonal != "none" ? length(y) > minimum(freq_seasonal) : true "Time series must be longer than the seasonal period if seasonal is added" if typeof(freq_seasonal) <: Vector - @assert all(freq_seasonal .> 0) "Seasonal period must be greater than 0" + (@assert all(freq_seasonal .> 0) "Seasonal period must be greater than 0") + else + (@assert freq_seasonal > 0 "Seasonal period must be greater than 0") end if typeof(cycle_period) <: Vector - @assert all(cycle_period .>= 0) "Cycle period must be greater than or equal to 0" + (@assert all(cycle_period .>= 0) "Cycle period must be greater than or equal to 0") + else + (@assert cycle_period >= 0 "Cycle period must be greater than or equal to 0") end if cycle_period == 0 - @assert !stochastic_cycle "stochastic_cycle must be false if cycle_period is 0" + (@assert cycle == "none" "stochastic_cycle and cycle must be false if cycle_period is 0") + else + nothing + end + + if !isnothing(dynamic_exog_coefs) + @assert all( + typeof(dynamic_exog_coefs[i][1]) <: Vector for + i in eachindex(dynamic_exog_coefs) + ) "The first element of each tuple in dynamic_exog_coefs must be a vector" + @assert all( + typeof(dynamic_exog_coefs[i][2]) <: String for + i in eachindex(dynamic_exog_coefs) + ) "The second element of each tuple in dynamic_exog_coefs must be a string" + @assert all([ + length(dynamic_exog_coefs[i][1]) .== length(y) for + i in eachindex(dynamic_exog_coefs) + ]) "The exogenous features that will be combined with state space components must have the same length as the time series" + @assert all( + dynamic_exog_coefs[i][2] in ["level", "slope", "seasonal", "cycle"] for + i in eachindex(dynamic_exog_coefs) + ) "The second element of each tuple in dynamic_exog_coefs must be a string that is either level, slope, seasonal or cycle" + for i in eachindex(dynamic_exog_coefs) + if dynamic_exog_coefs[i][2] == "seasonal" || + dynamic_exog_coefs[i][2] == "cycle" + @assert length(dynamic_exog_coefs[i]) == 3 "The tuple in dynamic_exog_coefs must have 3 elements if the second element is seasonal or cycle" + @assert typeof(dynamic_exog_coefs[i][3]) <: Int "The third element of each tuple in dynamic_exog_coefs must be an integer if the second element is seasonal or cycle" + @assert dynamic_exog_coefs[i][3] > 1 "The third element of each tuple in dynamic_exog_coefs must be greater than 1 if the second element is seasonal or cycle" + end + end end X = create_X( - level, - stochastic_level, - trend, - stochastic_trend, - seasonal, - stochastic_seasonal, + level in ["stochastic", "deterministic"], + level == "stochastic", + slope in ["stochastic", "deterministic"], + slope == "stochastic", + seasonal in ["stochastic", "deterministic"], + seasonal == "stochastic", + cycle in ["stochastic", "deterministic"], + cycle == "stochastic", freq_seasonal, - cycle_matrix, - stochastic_cycle, + cycle_period, outlier, - ζ_ω_threshold, + ζ_threshold, + ω_threshold, + ϕ_threshold, stochastic_start, - Exogenous_X, + exog, + dynamic_exog_coefs, ) - # convert y format into vector or matrix of AbstractFloat - if typeof(y) <: Vector - y = convert(Vector{AbstractFloat}, y) - else - y = convert(Matrix{AbstractFloat}, y) - end + # convert y format into vector of AbstractFloat + y = convert(Vector{AbstractFloat}, y) + return new( y, X, - level, - stochastic_level, - trend, - stochastic_trend, - seasonal, - stochastic_seasonal, + level in ["stochastic", "deterministic"], + level == "stochastic", + slope in ["stochastic", "deterministic"], + slope == "stochastic", + seasonal in ["stochastic", "deterministic"], + seasonal == "stochastic", + cycle in ["stochastic", "deterministic"], + cycle == "stochastic", freq_seasonal, cycle_period, - cycle_matrix, - stochastic_cycle, outlier, - ζ_ω_threshold, + ζ_threshold, + ω_threshold, + ϕ_threshold, stochastic_start, n_exogenous, + dynamic_exog_coefs, nothing, ) end @@ -188,38 +207,39 @@ end ξ_size(T::Int, stochastic_start::Int)::Int = T - max(2, stochastic_start) """ -ζ_size(T::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Int +ζ_size(T::Int, ζ_threshold::Int, stochastic_start::Int)::Int Calculates the size of ζ innovation matrix based on the input T. # Arguments - `T::Int`: Length of the original time series. - - `ζ_ω_threshold::Int`: Stabilize parameter ζ. + - `ζ_threshold::Int`: Stabilize parameter ζ. - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts. # Returns - `Int`: Size of ζ calculated from T. """ -ζ_size(T::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Int = - max(0, T - ζ_ω_threshold - max(2, stochastic_start)) +ζ_size(T::Int, ζ_threshold::Int, stochastic_start::Int)::Int = + max(0, T - ζ_threshold - max(2, stochastic_start)) """ -ω_size(T::Int, s::Int, stochastic_start::Int)::Int +ω_size(T::Int, s::Int, ω_threshold::Int, stochastic_start::Int)::Int Calculates the size of ω innovation matrix based on the input T. # Arguments - `T::Int`: Length of the original time series. - `s::Int`: Seasonal period. + - `ω_threshold::Int`: Stabilize parameter ω. - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts. # Returns - `Int`: Size of ω calculated from T. """ -ω_size(T::Int, s::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Int = - max(0, T - ζ_ω_threshold - s + 1 - max(0, max(2, stochastic_start) - s)) +ω_size(T::Int, s::Int, ω_threshold::Int, stochastic_start::Int)::Int = + max(0, T - ω_threshold - s + 1 - max(0, max(2, stochastic_start) - s)) """ o_size(T::Int, stochastic_start::Int)::Int @@ -237,74 +257,64 @@ o_size(T::Int, stochastic_start::Int)::Int o_size(T::Int, stochastic_start::Int)::Int = T - max(1, stochastic_start) + 1 """ - ϕ_size(T::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Int + ϕ_size(T::Int, ϕ_threshold::Int, stochastic_start::Int)::Int Calculates the size of ϕ innovation matrix based on the input T. # Arguments - `T::Int`: Length of the original time series. - - `ζ_ω_threshold::Int`: Stabilize parameter ζ. + - `ϕ_threshold::Int`: Stabilize parameter ϕ. - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts. # Returns - `Int`: Size of ϕ calculated from T. """ -function ϕ_size(T::Int, ζ_ω_threshold::Int, stochastic_start::Int) - ζ_ω_threshold = ζ_ω_threshold == 0 ? 1 : ζ_ω_threshold - if stochastic_start == 1 - return (2 * (T - max(2, stochastic_start) + 1) - (ζ_ω_threshold * 2)) - 2 - else - return (2 * (T - max(2, stochastic_start) + 1) - (ζ_ω_threshold * 2)) - end -end +ϕ_size(T::Int, ϕ_threshold::Int, stochastic_start::Int)::Int = + (2 * (T - max(2, stochastic_start) + 1) - (max(1, ϕ_threshold) * 2)) """ - create_ξ(T::Int, steps_ahead::Int, stochastic_start::Int)::Matrix + create_ξ(T::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. - - `steps_ahead::Int`: Number of steps ahead (for estimation purposes this should be set at 0). - `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, steps_ahead::Int, stochastic_start::Int)::Matrix +function create_ξ(T::Int, stochastic_start::Int)::Matrix stochastic_start = max(2, stochastic_start) - ξ_matrix = zeros(T + steps_ahead, T - stochastic_start + 1) + ξ_matrix = zeros(T, T - stochastic_start + 1) ones_indexes = findall( I -> Tuple(I)[1] - (stochastic_start - 2) > Tuple(I)[2], - CartesianIndices((T + steps_ahead, T - stochastic_start)), + CartesianIndices((T, T - stochastic_start)), ) ξ_matrix[ones_indexes] .= 1 return ξ_matrix[:, 1:(end - 1)] end """ -create_ζ(T::Int, steps_ahead::Int, ζ_ω_threshold::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. - - `steps_ahead::Int`: Number of steps ahead (for estimation purposes this should be set at 0). - - `ζ_ω_threshold::Int`: Stabilize parameter ζ. + - `ζ_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, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int -)::Matrix +function create_ζ(T::Int, ζ_threshold::Int, stochastic_start::Int)::Matrix stochastic_start = max(2, stochastic_start) - ζ_matrix = zeros(T + steps_ahead, T - stochastic_start) + ζ_matrix = zeros(T, T - stochastic_start) - for t in 2:(T + steps_ahead) + for t in 2:T if t < T len = t - stochastic_start ζ_matrix[t, 1:len] .= len:-1:1 @@ -312,33 +322,32 @@ function create_ζ( ζ_matrix[t, :] .= (t - stochastic_start):-1:(t - T + 1) end end - return ζ_matrix[:, 1:(end - ζ_ω_threshold)] + return ζ_matrix[:, 1:(end - ζ_threshold)] end """ -create_ω(T::Int, freq_seasonal::Int, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Matrix +create_ω(T::Int, freq_seasonal::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. - `freq_seasonal::Int`: Seasonal period. - - `steps_ahead::Int`: Number of steps ahead (for estimation purposes this should be set at 0). - - `ζ_ω_threshold::Int`: Stabilize parameter ζ. + - `ω_threshold::Int`: Stabilize parameter ω. # Returns - `Matrix`: Matrix of innovations ω constructed based on the input sizes. """ function create_ω( - T::Int, freq_seasonal::Int, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int + T::Int, freq_seasonal::Int, ω_threshold::Int, stochastic_start::Int )::Matrix stochastic_start = max(2, stochastic_start) - ω_matrix_size = T - freq_seasonal + 1 + ω_matrix_size = max(0, T - freq_seasonal + 1) stochastic_start_diff = max(0, stochastic_start - freq_seasonal) - ω_matrix = zeros(T + steps_ahead, ω_matrix_size - stochastic_start_diff) - for t in (freq_seasonal + 1):(T + steps_ahead) + ω_matrix = zeros(T, ω_matrix_size - stochastic_start_diff) + for t in (freq_seasonal + 1):T ωₜ_coefs = zeros(ω_matrix_size) Mₜ = Int(floor((t - 1) / freq_seasonal)) lag₁ = [t - j * freq_seasonal for j in 0:(Mₜ - 1)] @@ -349,104 +358,98 @@ function create_ω( -1 ω_matrix[t, :] = ωₜ_coefs[(1 + stochastic_start_diff):end] end - return ω_matrix[:, 1:(end - ζ_ω_threshold)] + return ω_matrix[:, 1:(end - ω_threshold)] end """ -create_o_matrix(T::Int, steps_ahead::Int, stochastic_start::Int)::Matrix +create_o_matrix(T::Int, stochastic_start::Int)::Matrix Creates a matrix of outliers 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. - - `steps_ahead::Int`: Number of steps ahead (for estimation purposes this should be set at 0). - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts. # Returns - `Matrix`: Matrix of outliers constructed based on the input sizes. """ -function create_o_matrix(T::Int, steps_ahead::Int, stochastic_start::Int)::Matrix +function create_o_matrix(T::Int, stochastic_start::Int)::Matrix stochastic_start = max(1, stochastic_start) rows = stochastic_start:T cols = 1:(T - stochastic_start + 1) values = ones(length(rows)) o_matrix = sparse(rows, cols, values, T, T - stochastic_start + 1) - return vcat(o_matrix, zeros(steps_ahead, length(cols))) + return o_matrix end """ -create_ϕ(X_cycle::Matrix, T::Int, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int)::Matrix +create_ϕ(c_period::Union{Int, Fl}, 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 - - `X_cycle::Matrix`: deterministic Cycle matrix. + - `c_period::Union{Int, Fl}`: Cycle period. - `T::Int`: Length of the original time series. - - `steps_ahead::Int64`: Number of steps ahead (for estimation purposes this should be set at 0). - - `ζ_ω_threshold::Int`: Stabilize parameter ζ. + - `ϕ_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_ϕ( - c_matrix::Matrix, T::Int, steps_ahead::Int, ζ_ω_threshold::Int, stochastic_start::Int -)::Matrix - num_cols = 2 * (T - stochastic_start + 1) - X = Matrix{Float64}(undef, T + steps_ahead, num_cols) - - for (idx, t) in enumerate(stochastic_start:T) - X[:, 2 * (idx - 1) + 1] = vcat( - zeros(t - 1), c_matrix[1:(T - t + 1 + steps_ahead), 1] - ) - X[:, 2 * (idx - 1) + 2] = vcat( - zeros(t - 1), c_matrix[1:(T - t + 1 + steps_ahead), 2] - ) + c_period::Union{Int,Fl}, T::Int, ϕ_threshold::Int, stochastic_start::Int +)::Matrix where {Fl<:AbstractFloat} + X = Matrix{Float64}(undef, T, 0) + λ = 2 * pi * (1:T) / c_period + + for t in max(2, stochastic_start):(T - max(1, ϕ_threshold)) # one of last two columns might be full of zeros + X_t = hcat(cos.(λ), sin.(λ)) + X_t[1:(t - 1), :] .= 0 + X = hcat(X, X_t) end - ζ_ω_threshold = ζ_ω_threshold == 0 ? 1 : ζ_ω_threshold - if stochastic_start == 1 - return X[:, 3:(end - (ζ_ω_threshold * 2))] - else - return X[:, 1:(end - (ζ_ω_threshold * 2))] - end + return round.(X, digits=5) end """ - create_deterministic_cycle_matrix(cycle_matrix::Vector{Matrix}, T::Int, steps_ahead::Int)::Vector{Matrix} +create_deterministic_seasonal(T::Int, s::Int)::Matrix - Creates a deterministic cycle matrix based on the input parameters. + Creates a matrix of deterministic seasonal components based on the input sizes. # Arguments - - `cycle_matrix::Vector{Matrix}`: Vector of cycle matrices. - `T::Int`: Length of the original time series. - - `steps_ahead::Int`: Number of steps ahead. + - `s::Int`: Seasonal period. +""" +function create_deterministic_seasonal(T::Int, s::Int)::Matrix + γ1_matrix = zeros(T, s) + for t in 1:T + γ1_matrix[t, t % s == 0 ? s : t % s] = 1.0 + end + return γ1_matrix +end - # Returns - - `Vector{Matrix}`: Deterministic cycle matrix constructed based on the input parameters. -""" -function create_deterministic_cycle_matrix( - cycle_matrix::Vector{Matrix}, T::Int, steps_ahead::Int -)::Vector{Matrix} - deterministic_cycle_matrix = Vector{Matrix}(undef, length(cycle_matrix)) - for (idx, c_matrix) in enumerate(cycle_matrix) - X_cycle = Matrix{Float64}(undef, T + steps_ahead, 2) - cycle_matrix_term = c_matrix^0 - X_cycle[1, :] = cycle_matrix_term[1, :] - for t in 2:(T + steps_ahead) - cycle_matrix_term *= c_matrix - X_cycle[t, :] = cycle_matrix_term[1, :] - end - deterministic_cycle_matrix[idx] = X_cycle - end - return deterministic_cycle_matrix +""" +create_deterministic_cycle(T::Int, c_period::Union{Int, Fl})::Matrix where {Fl<:AbstractFloat} + + Creates a matrix of deterministic cycle components based on the input sizes. + + # Arguments + - `T::Int`: Length of the original time series. + - `c_period::Int`: Cycle period. +""" +function create_deterministic_cycle( + T::Int, c_period::Union{Int,Fl} +)::Matrix where {Fl<:AbstractFloat} + λ = 2 * pi * (1:T) / c_period + cycle1_matrix = hcat(cos.(λ), sin.(λ)) + return cycle1_matrix end """ create_initial_states_Matrix( - T::Int, freq_seasonal::Union{Int, Vector{Int}}, steps_ahead::Int, level::Bool, trend::Bool, seasonal::Bool + T::Int, freq_seasonal::Union{Int, Vector{Int}}, level::Bool, trend::Bool, seasonal::Bool, cycle::Bool, cycle_period::Union{Int,Vector{Int}} )::Matrix Creates an initial states matrix based on the input parameters. @@ -454,10 +457,11 @@ end # Arguments - `T::Int`: Length of the original time series. - `freq_seasonal::Union{Int, Vector{Int}}`: Seasonal period. - - `steps_ahead::Int`: Number of steps ahead. - `level::Bool`: Flag for considering level component. - `trend::Bool`: Flag for considering trend component. - `seasonal::Bool`: Flag for considering seasonal component. + - `cycle::Bool`: Flag for considering cycle component. + - `cycle_period::Union{Int,Vector{Int}}`: Cycle period. # Returns - `Matrix`: Initial states matrix constructed based on the input parameters. @@ -466,78 +470,242 @@ end function create_initial_states_Matrix( T::Int, freq_seasonal::Union{Int,Vector{Int}}, - steps_ahead::Int, level::Bool, trend::Bool, seasonal::Bool, - deterministic_cycle_matrix::Vector{Matrix}, + cycle::Bool, + cycle_period::Union{Union{Int,<:AbstractFloat},Vector{Int},Vector{<:AbstractFloat}}, )::Matrix - initial_states_matrix = zeros(T + steps_ahead, 0) + initial_states_matrix = zeros(T, 0) if level - initial_states_matrix = hcat(initial_states_matrix, ones(T + steps_ahead, 1)) + initial_states_matrix = hcat(initial_states_matrix, ones(T, 1)) else nothing end if trend - initial_states_matrix = hcat( - initial_states_matrix, vcat([0], collect(1:(T + steps_ahead - 1))) - ) + initial_states_matrix = hcat(initial_states_matrix, vcat([0], collect(1:(T - 1)))) else nothing end if seasonal for s in freq_seasonal - γ1_matrix = zeros(T + steps_ahead, s) - for t in 1:(T + steps_ahead) - γ1_matrix[t, t % s == 0 ? s : t % s] = 1.0 - end + γ1_matrix = create_deterministic_seasonal(T, s) initial_states_matrix = hcat(initial_states_matrix, γ1_matrix) end end - if !isempty(deterministic_cycle_matrix) - for c_matrix in deterministic_cycle_matrix - initial_states_matrix = hcat(initial_states_matrix, c_matrix) + if cycle + for c_period in cycle_period + cycle1_matrix = create_deterministic_cycle(T, c_period) + initial_states_matrix = hcat(initial_states_matrix, cycle1_matrix) end end return initial_states_matrix end +""" +create_dynamic_exog_coefs_matrix(dynamic_exog_coefs::Vector{<:Tuple}, T::Int,ζ_threshold::Int, ω_threshold::Int, ϕ_threshold::Int, stochastic_start::Int)::Matrix + + Creates a matrix of combination components based on the input parameters. + + # Arguments + - `dynamic_exog_coefs::Vector{<:Tuple}`: Vector of tuples containing the combination components. + - `T::Int`: Length of the original time series. + - `ζ_threshold::Int`: Stabilize parameter ζ. + - `ω_threshold::Int`: Stabilize parameter ω. + - `ϕ_threshold::Int`: Stabilize parameter ϕ. + - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts. + + # Returns + - `Matrix`: Matrix of combination components constructed based on the input parameters. +""" +function create_dynamic_exog_coefs_matrix( + dynamic_exog_coefs::Vector{<:Tuple}, + T::Int, + ζ_threshold::Int, + ω_threshold::Int, + ϕ_threshold::Int, + stochastic_start::Int, +)::Matrix + state_components_dict = Dict{String,Matrix}() + dynamic_exog_coefs_matrix = zeros(T, 0) + for combination in dynamic_exog_coefs + if combination[2] == "level" + if haskey(state_components_dict, "level") + nothing + else + state_components_dict["level"] = hcat( + ones(T, 1), create_ξ(T, stochastic_start) + ) + end + key_name = "level" + elseif combination[2] == "slope" + if haskey(state_components_dict, "slope") + nothing + else + state_components_dict["slope"] = hcat( + vcat([0], collect(1:(T - 1))), + create_ζ(T, ζ_threshold, stochastic_start), + ) + end + key_name = "slope" + elseif combination[2] == "seasonal" + if haskey(state_components_dict, "seasonal_$(combination[3])") + nothing + else + state_components_dict["seasonal_$(combination[3])"] = hcat( + create_deterministic_seasonal(T, combination[3]), + create_ω(T, combination[3], ω_threshold, stochastic_start), + ) + end + key_name = "seasonal_$(combination[3])" + elseif combination[2] == "cycle" + if haskey(state_components_dict, "cycle_$(combination[3])") + nothing + else + state_components_dict["cycle_$(combination[3])"] = hcat( + create_deterministic_cycle(T, combination[3]), + create_ϕ(combination[3], T, ϕ_threshold, stochastic_start), + ) + end + key_name = "cycle_$(combination[3])" + end + dynamic_exog_coefs_matrix = hcat( + dynamic_exog_coefs_matrix, combination[1] .* state_components_dict[key_name] + ) + end + return dynamic_exog_coefs_matrix +end + +""" +create_forecast_dynamic_exog_coefs_matrix(dynamic_exog_coefs::Vector{<:Tuple}, T::Int, steps_ahead::Int, ζ_threshold::Int, ω_threshold::Int, ϕ_threshold::Int, stochastic_start::Int)::Matrix + + Creates a matrix of combination components based on the input parameters. + + # Arguments + - `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 ϕ. + - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts. + + # Returns + - `Matrix`: Matrix of combination components constructed based on the input parameters. +""" +function create_forecast_dynamic_exog_coefs_matrix( + dynamic_exog_coefs::Vector{<:Tuple}, + T::Int, + steps_ahead::Int, + ζ_threshold::Int, + ω_threshold::Int, + ϕ_threshold::Int, + stochastic_start::Int, +)::Matrix + state_components_dict = Dict{String,Matrix}() + dynamic_exog_coefs_matrix = zeros(steps_ahead, 0) + for combination in dynamic_exog_coefs + if combination[2] == "level" + if haskey(state_components_dict, "level") + nothing + else + state_components_dict["level"] = hcat( + ones(T + steps_ahead, 1), create_ξ(T + steps_ahead, stochastic_start) + )[ + (end - steps_ahead + 1):end, 1:combination[4] + ] + end + key_name = "level" + elseif combination[2] == "slope" + if haskey(state_components_dict, "slope") + nothing + else + state_components_dict["slope"] = hcat( + vcat([0], collect(1:(T + steps_ahead - 1))), + create_ζ(T + steps_ahead, ζ_threshold, stochastic_start), + )[ + (end - steps_ahead + 1):end, 1:combination[4] + ] + end + key_name = "slope" + elseif combination[2] == "seasonal" + if haskey(state_components_dict, "seasonal_$(combination[3])") + nothing + else + state_components_dict["seasonal_$(combination[3])"] = hcat( + create_deterministic_seasonal(T + steps_ahead, combination[3]), + create_ω( + T + steps_ahead, combination[3], ω_threshold, stochastic_start + ), + )[ + (end - steps_ahead + 1):end, 1:combination[4] + ] + end + key_name = "seasonal_$(combination[3])" + elseif combination[2] == "cycle" + if haskey(state_components_dict, "cycle_$(combination[3])") + nothing + else + state_components_dict["cycle_$(combination[3])"] = hcat( + create_deterministic_cycle(T + steps_ahead, combination[3]), + create_ϕ( + combination[3], T + steps_ahead, ϕ_threshold, stochastic_start + ), + )[ + (end - steps_ahead + 1):end, 1:combination[4] + ] + end + key_name = "cycle_$(combination[3])" + end + dynamic_exog_coefs_matrix = hcat( + dynamic_exog_coefs_matrix, combination[1] .* state_components_dict[key_name] + ) + end + return dynamic_exog_coefs_matrix +end + """ create_X( level::Bool, stochastic_level::Bool, - trend::Bool, - stochastic_trend::Bool, + slope::Bool, + stochastic_slope::Bool, seasonal::Bool, stochastic_seasonal::Bool, - freq_seasonal::Union{Int, Vector{Int}}, + cycle::Bool, + stochastic_cycle::Bool, + freq_seasonal::Union{Int,Vector{Int}}, + cycle_period::Union{Int,Vector{Int}}, outlier::Bool, - ζ_ω_threshold::Int, - stochastic_start::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} + ζ_threshold::Int, + ω_threshold::Int, + ϕ_threshold::Int, + stochastic_start::Int, + exog::Matrix{Fl}, +) where {Fl<:AbstractFloat} Creates the StateSpaceLearning matrix X based on the model type and input parameters. # Arguments - `level::Bool`: Flag for considering level component. - `stochastic_level::Bool`: Flag for considering stochastic level component. - - `trend::Bool`: Flag for considering trend component. - - `stochastic_trend::Bool`: Flag for considering stochastic trend component. + - `slope::Bool`: Flag for considering slope component. + - `stochastic_slope::Bool`: Flag for considering stochastic slope component. - `seasonal::Bool`: Flag for considering seasonal component. - `stochastic_seasonal::Bool`: Flag for considering stochastic seasonal component. + - `cycle::Bool`: Flag for considering cycle component. + - `stochastic_cycle::Bool`: Flag for considering stochastic cycle component. - `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 ϕ. - `stochastic_start::Int`: parameter to set at which time stamp the stochastic component starts. - - `Exogenous_X::Matrix{Fl}`: Exogenous variables matrix. - - `steps_ahead::Int`: Number of steps ahead (default: 0). - - `Exogenous_Forecast::Matrix{Fl}`: Exogenous variables forecast matrix (default: zeros). + - `exog::Matrix{Fl}`: Exogenous variables matrix. # Returns - `Matrix`: StateSpaceLearning matrix X constructed based on the input parameters. @@ -545,64 +713,72 @@ create_X( function create_X( level::Bool, stochastic_level::Bool, - trend::Bool, - stochastic_trend::Bool, + slope::Bool, + stochastic_slope::Bool, seasonal::Bool, stochastic_seasonal::Bool, - freq_seasonal::Union{Int,Vector{Int}}, - cycle_matrix::Vector{Matrix}, + cycle::Bool, stochastic_cycle::Bool, + 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, stochastic_start::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} - T = size(Exogenous_X, 1) + exog::Matrix{Fl}, + dynamic_exog_coefs::Union{Vector{<:Tuple},Nothing}, +) where {Fl<:AbstractFloat} + T = size(exog, 1) ξ_matrix = if stochastic_level - create_ξ(T, steps_ahead, stochastic_start) + create_ξ(T, stochastic_start) else - zeros(T + steps_ahead, 0) + zeros(T, 0) end - ζ_matrix = if stochastic_trend - create_ζ(T, steps_ahead, ζ_ω_threshold, stochastic_start) + ζ_matrix = if stochastic_slope + create_ζ(T, ζ_threshold, stochastic_start) else - zeros(T + steps_ahead, 0) + zeros(T, 0) end - ω_matrix = zeros(T + steps_ahead, 0) + ω_matrix = zeros(T, 0) if stochastic_seasonal for s in freq_seasonal - ω_matrix = hcat( - ω_matrix, create_ω(T, s, steps_ahead, ζ_ω_threshold, stochastic_start) - ) + ω_matrix = hcat(ω_matrix, create_ω(T, s, ω_threshold, stochastic_start)) end end - deterministic_cycle_matrix = create_deterministic_cycle_matrix( - cycle_matrix, T, steps_ahead - ) - ϕ_matrix = zeros(T + steps_ahead, 0) + ϕ_matrix = zeros(T, 0) if stochastic_cycle - for c_matrix in deterministic_cycle_matrix - ϕ_matrix = hcat( - ϕ_matrix, - create_ϕ(c_matrix, T, steps_ahead, ζ_ω_threshold, stochastic_start), - ) + for c_period in cycle_period + ϕ_matrix = hcat(ϕ_matrix, create_ϕ(c_period, T, ϕ_threshold, stochastic_start)) end end o_matrix = if outlier - create_o_matrix(T, steps_ahead, stochastic_start) + create_o_matrix(T, stochastic_start) else - zeros(T + steps_ahead, 0) + zeros(T, 0) end initial_states_matrix = create_initial_states_Matrix( - T, freq_seasonal, steps_ahead, level, trend, seasonal, deterministic_cycle_matrix + T, freq_seasonal, level, slope, seasonal, cycle, cycle_period ) + + dynamic_exog_coefs_matrix = if !isnothing(dynamic_exog_coefs) + create_dynamic_exog_coefs_matrix( + dynamic_exog_coefs, + T, + ζ_threshold, + ω_threshold, + ϕ_threshold, + stochastic_start, + ) + else + zeros(T, 0) + end + return hcat( initial_states_matrix, ξ_matrix, @@ -610,51 +786,8 @@ function create_X( ω_matrix, ϕ_matrix, o_matrix, - vcat(Exogenous_X, Exogenous_Forecast), - ) -end - -""" -create_X( - model::StructuralModel, - 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 and input parameters. - - # Arguments - - `model::StructuralModel`: StructuralModel object. - - `Exogenous_X::Matrix{Fl}`: Exogenous variables matrix. - - `steps_ahead::Int`: Number of steps ahead (default: 0). - - `Exogenous_Forecast::Matrix{Fl}`: Exogenous variables forecast matrix (default: zeros). - - # Returns - - `Matrix`: StateSpaceLearning matrix X constructed based on the input parameters. -""" -function create_X( - model::StructuralModel, - Exogenous_X::Matrix{Fl}, - steps_ahead::Int=0, - Exogenous_Forecast::Matrix{Tl}=zeros(steps_ahead, size(Exogenous_X, 2)), -) where {Fl<:AbstractFloat,Tl<:AbstractFloat} - return create_X( - model.level, - model.stochastic_level, - model.trend, - model.stochastic_trend, - model.seasonal, - model.stochastic_seasonal, - model.freq_seasonal, - model.cycle_matrix, - model.stochastic_cycle, - model.outlier, - model.ζ_ω_threshold, - model.stochastic_start, - Exogenous_X, - steps_ahead, - Exogenous_Forecast, + exog, + dynamic_exog_coefs_matrix, ) end @@ -671,7 +804,7 @@ function get_components_indexes(model::StructuralModel)::Dict """ function get_components_indexes(model::StructuralModel)::Dict - T = typeof(model.y) <: Vector ? length(model.y) : size(model.y, 1) + T = length(model.y) FINAL_INDEX = 0 @@ -684,7 +817,7 @@ function get_components_indexes(model::StructuralModel)::Dict initial_states_indexes = Int[] end - if model.trend + if model.slope ν1_indexes = [FINAL_INDEX + 1] initial_states_indexes = vcat(initial_states_indexes, ν1_indexes) FINAL_INDEX += length(ν1_indexes) @@ -703,8 +836,8 @@ function get_components_indexes(model::StructuralModel)::Dict end c_indexes = Vector{Int}[] - if !isempty(model.cycle_matrix) - for _ in eachindex(model.cycle_matrix) + if model.cycle + for _ in eachindex(model.cycle_period) c_i_indexes = collect((FINAL_INDEX + 1):(FINAL_INDEX + 2)) initial_states_indexes = vcat(initial_states_indexes, c_i_indexes) FINAL_INDEX += length(c_i_indexes) @@ -721,10 +854,10 @@ function get_components_indexes(model::StructuralModel)::Dict ξ_indexes = Int[] end - if model.stochastic_trend + if model.stochastic_slope ζ_indexes = collect( (FINAL_INDEX + 1):(FINAL_INDEX + ζ_size( - T, model.ζ_ω_threshold, model.stochastic_start + T, model.ζ_threshold, model.stochastic_start )), ) FINAL_INDEX += length(ζ_indexes) @@ -737,7 +870,7 @@ function get_components_indexes(model::StructuralModel)::Dict for s in model.freq_seasonal ω_s_indexes = collect( (FINAL_INDEX + 1):(FINAL_INDEX + ω_size( - T, s, model.ζ_ω_threshold, model.stochastic_start + T, s, model.ω_threshold, model.stochastic_start )), ) FINAL_INDEX += length(ω_s_indexes) @@ -749,10 +882,10 @@ function get_components_indexes(model::StructuralModel)::Dict ϕ_indexes = Vector{Int}[] if model.stochastic_cycle - for _ in eachindex(model.cycle_matrix) + for _ in eachindex(model.cycle_period) ϕ_i_indexes = collect( (FINAL_INDEX + 1):(FINAL_INDEX + ϕ_size( - T, model.ζ_ω_threshold, model.stochastic_start + T, model.ϕ_threshold, model.stochastic_start )), ) FINAL_INDEX += length(ϕ_i_indexes) @@ -773,14 +906,17 @@ function get_components_indexes(model::StructuralModel)::Dict exogenous_indexes = collect((FINAL_INDEX + 1):(FINAL_INDEX + model.n_exogenous)) + dynamic_exog_coefs_indexes = collect((FINAL_INDEX + 1):size(model.X, 2)) + components_indexes_dict = Dict( "μ1" => μ1_indexes, "ν1" => ν1_indexes, "ξ" => ξ_indexes, "ζ" => ζ_indexes, "o" => o_indexes, - "Exogenous_X" => exogenous_indexes, + "exog" => exogenous_indexes, "initial_states" => initial_states_indexes, + "dynamic_exog_coefs" => dynamic_exog_coefs_indexes, ) for (i, s) in enumerate(model.freq_seasonal) @@ -793,11 +929,11 @@ function get_components_indexes(model::StructuralModel)::Dict end end - if !isempty(model.cycle_matrix) + if model.cycle for i in eachindex(model.cycle_period) - components_indexes_dict["c1_$i"] = c_indexes[i] + components_indexes_dict["c1_$(model.cycle_period[i])"] = c_indexes[i] if model.stochastic_cycle - components_indexes_dict["ϕ_$i"] = ϕ_indexes[i] + components_indexes_dict["ϕ_$(model.cycle_period[i])"] = ϕ_indexes[i] end end end @@ -900,7 +1036,7 @@ function get_model_innovations(model::StructuralModel) push!(model_innovations, "ξ") end - if model.stochastic_trend + if model.stochastic_slope push!(model_innovations, "ζ") end @@ -911,7 +1047,7 @@ function get_model_innovations(model::StructuralModel) end if model.stochastic_cycle - for i in eachindex(model.cycle_period) + for i in model.cycle_period push!(model_innovations, "ϕ_$i") end end @@ -919,41 +1055,644 @@ function get_model_innovations(model::StructuralModel) end """ - get_innovation_functions(model::StructuralModel, innovation::String)::Function + get_trend_decomposition(model::StructuralModel, components::Dict, slope::Vector{AbstractFloat})::Vector{AbstractFloat} + + Returns the level component and associated innovation vectors. + + # Arguments + - `model::StructuralModel`: StructuralModel object. + - `components::Dict` : Components dict. + - `slope::Vector{AbstractFloat}`: Time-series of the slope component. + + # Returns + - `Vector{AbstractFloat}`: Time-series of the level component. + +""" +function get_trend_decomposition( + model::StructuralModel, components::Dict, slope::Vector{AbstractFloat} +)::Vector{AbstractFloat} + T = size(model.y, 1) + trend = Vector{AbstractFloat}(undef, T) + + if model.level + trend[1] = components["μ1"]["Coefs"][1] + else + trend[1] = 0.0 + end + + if model.stochastic_level + ξ = vcat(zeros(max(2, model.stochastic_start) - 1), components["ξ"]["Coefs"], 0.0) + @assert length(ξ) == T + else + ξ = zeros(AbstractFloat, T) + end + + for t in 2:T + trend[t] = trend[t - 1] + slope[t] + ξ[t] + end + + return trend +end + +""" + get_slope_decomposition(model::StructuralModel, components::Dict)::Vector{AbstractFloat} + + Returns the slope component and associated innovation vectors. + + # Arguments + - `model::StructuralModel`: StructuralModel object. + - `components::Dict`: Components dict.. + + # Returns + - `Vector{AbstractFloat}`: Time-series of the slope component. + +""" +function get_slope_decomposition( + model::StructuralModel, components::Dict +)::Vector{AbstractFloat} + T = size(model.y, 1) + slope = Vector{AbstractFloat}(undef, T) + + if model.slope + slope[1] = components["ν1"]["Coefs"][1] + else + slope[1] = 0.0 + end + + if model.stochastic_slope + ζ = vcat( + zeros(max(2, model.stochastic_start)), + components["ζ"]["Coefs"], + zeros(model.ζ_threshold), + ) + @assert length(ζ) == T + else + ζ = zeros(AbstractFloat, T) + end + + for t in 2:T + slope[t] = slope[t - 1] + ζ[t] + end + + return slope +end + +""" + get_seasonal_decomposition(model::StructuralModel, components::Dict, s::Int)::Vector{AbstractFloat} - Returns the innovation function based on the input innovation string. + Returns the seasonality component and associated innovation vectors. # Arguments - `model::StructuralModel`: StructuralModel object. - - `innovation::String`: Innovation string. - - steps_ahead::Int: Number of steps ahead. + - `components::Dict`: Components dict. + - `s::Int`: Seasonal frequency. # Returns + - `Vector{AbstractFloat}`: Time-series of the seasonality component. """ -function get_innovation_simulation_X( - model::StructuralModel, innovation::String, steps_ahead::Int -) - if innovation == "ξ" - return create_ξ(length(model.y) + steps_ahead + 1, 0, model.stochastic_start) - elseif innovation == "ζ" - return create_ζ(length(model.y) + steps_ahead + 1, 0, 1, model.stochastic_start) - elseif occursin("ω_", innovation) - s = parse(Int, split(innovation, "_")[2]) - return create_ω(length(model.y) + steps_ahead + 1, s, 0, 1, model.stochastic_start) - elseif occursin("ϕ_", innovation) - i = parse(Int, split(innovation, "_")[2]) - deterministic_cycle_matrix = create_deterministic_cycle_matrix( - model.cycle_matrix, length(model.y), steps_ahead + 1 +function get_seasonal_decomposition( + model::StructuralModel, components::Dict, s::Int +)::Vector{AbstractFloat} + T = size(model.y, 1) + seasonal = Vector{AbstractFloat}(undef, T) + + if model.seasonal + seasonal[1:s] = components["γ1_$(s)"]["Coefs"] + else + seasonal[1:s] = zeros(AbstractFloat, s) + end + + if model.stochastic_seasonal + ω = vcat( + zeros(s - 1 + max(0, max(2, model.stochastic_start) - s)), + components["ω_$(s)"]["Coefs"], + zeros(model.ω_threshold), + ) + @assert length(ω) == T + else + ω = zeros(AbstractFloat, T) + end + + for t in (s + 1):T + seasonal[t] = seasonal[t - s] + ω[t] - ω[t - 1] + end + + return seasonal +end + +""" + get_cycle_decomposition(model::StructuralModel, components::Dict, cycle_period::Union{AbstractFloat, Int})::Tuple{Vector{AbstractFloat}, Vector{AbstractFloat}} + + Returns the cycle component and associated innovation vectors. + + # Arguments + - `model::StructuralModel`: StructuralModel object. + - `components::Dict`: Components dict. + - `cycle_period::Union{AbstractFloat, Int}`: Cycle period. + + # Returns + - `Tuple{Vector{AbstractFloat}, Vector{AbstractFloat}}`: Tuple containing the cycle component and the cycle component hat. + +""" +function get_cycle_decomposition( + model::StructuralModel, components::Dict, cycle_period::Union{AbstractFloat,Int} +)::Vector{AbstractFloat} + T = size(model.y, 1) + cycle = Vector{AbstractFloat}(undef, T) + + if cycle_period != 0 + λ = 2 * pi * (1:T) / cycle_period + c1 = components["c1_$(cycle_period)"]["Coefs"] + + cycle[1] = (dot(c1, [cos(λ[1]), sin(λ[1])])) + + if model.stochastic_cycle + ϕ_cos = vcat( + zeros(max(2, model.stochastic_start) - 1), + components["ϕ_$(cycle_period)"]["Coefs"][1:2:end], + zeros(max(1, model.ϕ_threshold)), + ) + ϕ_sin = vcat( + zeros(max(2, model.stochastic_start) - 1), + components["ϕ_$(cycle_period)"]["Coefs"][2:2:end], + zeros(max(1, model.ϕ_threshold)), + ) + @assert length(ϕ_cos) == T + @assert length(ϕ_sin) == T + else + ϕ_cos = zeros(AbstractFloat, T) + ϕ_sin = zeros(AbstractFloat, T) + end + + for t in 2:T + ϕ_indexes = + max(2, model.stochastic_start):min(t, (T - max(1, model.ϕ_threshold))) + cycle[t] = + dot(c1, [cos(λ[t]), sin(λ[t])]) + sum( + ϕ_cos[i] * cos(λ[t]) + ϕ_sin[i] * sin(λ[t]) for + i in eachindex(ϕ_indexes) + ) + end + + else + cycle = zeros(AbstractFloat, T) + end + + return cycle +end + +""" + get_model_decomposition(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_model_decomposition(model::StructuralModel, components::Dict)::Dict + freq_seasonal = model.freq_seasonal + cycle_period = model.cycle_period + model_decomposition = Dict() + + if model.slope + slope = get_slope_decomposition(model, components) + model_decomposition["slope"] = slope + end + + if model.level || model.slope + slope = model.slope ? slope : convert(Vector{AbstractFloat}, zeros(length(model.y))) + trend = get_trend_decomposition(model, components, slope) + model_decomposition["trend"] = trend + end + + if model.seasonal + for s in freq_seasonal + seasonal = get_seasonal_decomposition(model, components, s) + model_decomposition["seasonal_$s"] = seasonal + end + end + + if model.cycle + for i in cycle_period + cycle, cycle_hat = get_cycle_decomposition(model, components, i) + model_decomposition["cycle_$i"] = cycle + model_decomposition["cycle_hat_$i"] = cycle_hat + end + end + return model_decomposition +end + +""" + simulate_states( + model::StructuralModel, steps_ahead::Int, punctual::Bool, seasonal_innovation_simulation::Int + )::Vector{AbstractFloat} + + Simulates the states of the model. + + # Arguments + - `model::StructuralModel`: StructuralModel object. + - `steps_ahead::Int`: Steps ahead. + - `punctual::Bool`: Flag for considering punctual forecast. + - `seasonal_innovation_simulation::Int`: Flag for considering seasonal innovation simulation. + + # Returns + - `Vector{AbstractFloat}`: Vector of states. +""" +function simulate_states( + model::StructuralModel, + steps_ahead::Int, + punctual::Bool, + seasonal_innovation_simulation::Int, +)::Vector{AbstractFloat} + T = length(model.y) + + prediction = AbstractFloat[] + + if model.slope + slope = deepcopy(model.output.decomposition["slope"]) + start_idx = max(2, model.stochastic_start) + 1 + final_idx = T - model.ζ_threshold + if model.stochastic_slope && !punctual + if seasonal_innovation_simulation != 0 + ζ_values = vcat( + zeros(start_idx - 1), + model.output.components["ζ"]["Coefs"], + zeros(model.ζ_threshold), + ) + else + ζ_values = model.output.components["ζ"]["Coefs"] + end + else + ζ_values = zeros(T) + end + stochastic_slope_set = get_stochastic_values( + ζ_values, steps_ahead, T, start_idx, final_idx, seasonal_innovation_simulation + ) + else + slope = zeros(T) + end + + if model.level || model.slope + trend = deepcopy(model.output.decomposition["trend"]) + start_idx = max(2, model.stochastic_start) + final_idx = T - 1 + if model.stochastic_level && !punctual + if seasonal_innovation_simulation != 0 + ξ_values = vcat( + zeros(start_idx - 1), model.output.components["ξ"]["Coefs"], zeros(1) + ) + else + ξ_values = model.output.components["ξ"]["Coefs"] + end + else + ξ_values = zeros(T) + end + stochastic_level_set = get_stochastic_values( + ξ_values, steps_ahead, T, start_idx, final_idx, seasonal_innovation_simulation + ) + end + + if model.seasonal + seasonals = [ + deepcopy(model.output.decomposition["seasonal_$s"]) for s in model.freq_seasonal + ] + start_idx = [ + model.freq_seasonal[i] - 1 + + max(0, max(2, model.stochastic_start) - model.freq_seasonal[i]) for + i in eachindex(model.freq_seasonal) + ] + final_idx = [T - model.ω_threshold for _ in eachindex(model.freq_seasonal)] + if model.ω_threshold == 0 + final_ω = [ + model.output.components["ω_$(s)"]["Coefs"][end] for s in model.freq_seasonal + ] + else + final_ω = [0.0 for _ in model.freq_seasonal] + end + if model.stochastic_seasonal && !punctual + if seasonal_innovation_simulation != 0 + ω_values = [ + vcat( + zeros(s - 1 + max(0, max(2, model.stochastic_start) - s)), + model.output.components["ω_$(s)"]["Coefs"], + zeros(model.ω_threshold), + ) for s in model.freq_seasonal + ] + else + ω_values = [ + model.output.components["ω_$(s)"]["Coefs"] for s in model.freq_seasonal + ] + end + else + ω_values = [zeros(T) for _ in model.freq_seasonal] + end + stochastic_seasonals_set = [ + vcat( + final_ω[i], + get_stochastic_values( + ω_values[i], + steps_ahead, + T, + start_idx[i], + final_idx[i], + seasonal_innovation_simulation, + ), + ) for i in eachindex(model.freq_seasonal) + ] + end + + if model.cycle + start_idx = [max(2, model.stochastic_start) for _ in eachindex(model.cycle_period)] + final_idx = [T - max(1, model.ϕ_threshold) for _ in eachindex(model.cycle_period)] + if model.stochastic_cycle && !punctual + if seasonal_innovation_simulation != 0 + ϕ_cos_values = [ + vcat( + zeros(max(2, model.stochastic_start) - 1), + model.output.components["ϕ_$(i)"]["Coefs"][1:2:end], + zeros(max(1, model.ϕ_threshold)), + ) for i in model.cycle_period + ] + ϕ_sin_values = [ + vcat( + zeros(max(2, model.stochastic_start) - 1), + model.output.components["ϕ_$(i)"]["Coefs"][2:2:end], + zeros(max(1, model.ϕ_threshold)), + ) for i in model.cycle_period + ] + else + ϕ_cos_values = [ + model.output.components["ϕ_$(i)"]["Coefs"] for i in model.cycle_period + ] + ϕ_sin_values = [ + model.output.components["ϕ_$(i)"]["Coefs"][2:2:end] for + i in model.cycle_period + ] + end + else + ϕ_cos_values = [zeros(T) for _ in model.cycle_period] + ϕ_sin_values = [zeros(T) for _ in model.cycle_period] + end + stochastic_cycles_cos_set = [ + get_stochastic_values( + ϕ_cos_values[i], + steps_ahead, + T, + start_idx[i], + final_idx[i], + seasonal_innovation_simulation, + ) for i in eachindex(model.cycle_period) + ] + stochastic_cycles_sin_set = [ + get_stochastic_values( + ϕ_sin_values[i], + steps_ahead, + T, + start_idx[i], + final_idx[i], + seasonal_innovation_simulation, + ) for i in eachindex(model.cycle_period) + ] + end + + if model.outlier && !punctual + start_idx = 1 + final_idx = T + outlier_values = model.output.components["o"]["Coefs"] + #stochastic_outliers_set = get_stochastic_values(outlier_values, steps_ahead, T, 1, T, seasonal_innovation_simulation) + stochastic_outliers_set = rand(outlier_values, steps_ahead) + end + + if !punctual + #stochastic_residuals_set = get_stochastic_values(model.output.ε, steps_ahead, T, 1, T, seasonal_innovation_simulation) + stochastic_residuals_set = rand(model.output.ε, steps_ahead) + end + + for t in (T + 1):(T + steps_ahead) + slope_t = model.slope ? slope[end] + stochastic_slope_set[t - T] : 0.0 + + trend_t = if (model.level || model.slope) + trend[end] + slope[end] + stochastic_level_set[t - T] + else + 0.0 + end + + if model.seasonal + seasonals_t = [ + seasonals[i][t - model.freq_seasonal[i]] + + stochastic_seasonals_set[i][t - T + 1] - stochastic_seasonals_set[i][t - T] + for i in eachindex(model.freq_seasonal) + ] + else + seasonals_t = zeros(AbstractFloat, length(model.freq_seasonal)) + end + + if model.cycle_period != 0 + cycles_t = zeros(AbstractFloat, length(model.cycle_period)) + for i in eachindex(model.cycle_period) + ϕ_cos = model.output.components["ϕ_$(model.cycle_period[i])"]["Coefs"][1:2:end] + ϕ_sin = model.output.components["ϕ_$(model.cycle_period[i])"]["Coefs"][2:2:end] + λ = 2 * pi * (1:(T + steps_ahead)) / model.cycle_period[i] + + cycle_t = + dot( + model.output.components["c1_$(model.cycle_period[i])"]["Coefs"], + [cos(λ[t]), sin(λ[t])], + ) + + sum( + ϕ_cos[j] * cos(λ[t]) + ϕ_sin[j] * sin(λ[t]) for + j in eachindex(ϕ_cos) + ) + + sum( + stochastic_cycles_cos_set[i][j] * cos(λ[t]) + + stochastic_cycles_sin_set[i][j] * sin(λ[t]) for + j in eachindex(stochastic_cycles_cos_set[i][1:(t - T)]) + ) + cycles_t[i] = cycle_t + end + else + cycles_t = zeros(AbstractFloat, length(model.cycle_period)) + end + + outlier_t = (model.outlier && !punctual) ? stochastic_outliers_set[t - T] : 0.0 + residuals_t = !punctual ? stochastic_residuals_set[t - T] : 0.0 + + push!( + prediction, trend_t + sum(seasonals_t) + sum(cycles_t) + outlier_t + residuals_t ) - return create_ϕ( - deterministic_cycle_matrix[i], - length(model.y) + steps_ahead + 1, - 0, - model.ζ_ω_threshold, + model.slope ? push!(slope, slope_t) : nothing + model.level ? push!(trend, trend_t) : nothing + if model.seasonal + for i in eachindex(model.freq_seasonal) + seasonals[i] = vcat(seasonals[i], seasonals_t[i]) + end + end + end + + return prediction +end + +""" + forecast_dynamic_exog_coefs(model::StructuralModel, steps_ahead::Int, dynamic_exog_coefs_forecasts::Vector{<:Vector})::Vector{AbstractFloat} + + Returns the prediction of the combination components terms. + + # Arguments + - `model::StructuralModel`: StructuralModel object. + - `steps_ahead::Int`: Number of steps ahead for forecasting. + - `dynamic_exog_coefs_forecasts::Vector{<:Vector}`: Vector of vectors of combination components forecasts. + + # Returns + - `Vector{AbstractFloat}`: Vector of combination components forecasts. +""" +function forecast_dynamic_exog_coefs( + model::StructuralModel, steps_ahead::Int, dynamic_exog_coefs_forecasts::Vector{<:Vector} +)::Vector{AbstractFloat} + if !isempty(dynamic_exog_coefs_forecasts) + T = length(model.y) + 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) + extra_param = "" + elseif model.dynamic_exog_coefs[i][2] == "slope" + n_coefs = 1 + ζ_size(T, model.ζ_threshold, model.stochastic_start) + extra_param = "" + elseif model.dynamic_exog_coefs[i][2] == "seasonal" + n_coefs = + model.dynamic_exog_coefs[i][3] + ω_size( + T, + model.dynamic_exog_coefs[i][3], + model.ω_threshold, + model.stochastic_start, + ) + extra_param = model.dynamic_exog_coefs[i][3] + elseif model.dynamic_exog_coefs[i][2] == "cycle" + n_coefs = 2 + ϕ_size(T, model.ϕ_threshold, model.stochastic_start) + extra_param = model.dynamic_exog_coefs[i][3] + end + dynamic_exog_coefs[i] = ( + dynamic_exog_coefs_forecasts[i], + model.dynamic_exog_coefs[i][2], + extra_param, + n_coefs, + ) + end + dynamic_exog_coefs_forecasts_matrix = create_forecast_dynamic_exog_coefs_matrix( + dynamic_exog_coefs, + T, + steps_ahead, + model.ζ_threshold, + model.ω_threshold, + model.ϕ_threshold, model.stochastic_start, ) + dynamic_exog_coefs_prediction = + dynamic_exog_coefs_forecasts_matrix * + model.output.components["dynamic_exog_coefs"]["Coefs"] + else + dynamic_exog_coefs_prediction = zeros(steps_ahead) + end + return dynamic_exog_coefs_prediction +end + +""" + forecast(model::StructuralModel, steps_ahead::Int; Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, 0), dynamic_exog_coefs_forecasts::Vector{<:Vector}=Vector{Vector}(undef, 0))::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. + - `Exogenous_Forecast::Matrix{Fl}`: Matrix of forecasts of exogenous variables. + - `dynamic_exog_coefs_forecasts::Vector{<:Vector}`: Vector of vectors of combination components forecasts. + +""" +function forecast( + model::StructuralModel, + steps_ahead::Int; + Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, 0), + dynamic_exog_coefs_forecasts::Vector{<:Vector}=Vector{Vector}(undef, 0), +)::Vector{AbstractFloat} where {Fl<:AbstractFloat} + states_prediction = simulate_states(model, steps_ahead, true, 0) + + @assert size(Exogenous_Forecast, 1) == steps_ahead + @assert all( + length(dynamic_exog_coefs_forecasts[i]) == steps_ahead for + i in eachindex(dynamic_exog_coefs_forecasts) + ) + if !isnothing(model.dynamic_exog_coefs) + (@assert length(dynamic_exog_coefs_forecasts) == length(model.dynamic_exog_coefs)) + else + nothing + end + if (dynamic_exog_coefs_forecasts == Vector{Vector}(undef, 0)) + (@assert isnothing(model.dynamic_exog_coefs)) + else + nothing + end + @assert size(Exogenous_Forecast, 2) == model.n_exogenous + + dynamic_exog_coefs_prediction = forecast_dynamic_exog_coefs( + model, steps_ahead, dynamic_exog_coefs_forecasts + ) + + prediction = + states_prediction + + (Exogenous_Forecast * model.output.components["exog"]["Coefs"]) + + dynamic_exog_coefs_prediction + + return prediction +end + +""" + simulate(model::StructuralModel, steps_ahead::Int, N_scenarios::Int; Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, 0), dynamic_exog_coefs_forecasts::Vector{<:Vector}=Vector{Vector}(undef, 0), seasonal_innovation_simulation::Int=0, seed::Int=1234)::Matrix{AbstractFloat} + + Returns a matrix of scenarios of the states of the model. + + # Arguments + - `model::StructuralModel`: StructuralModel object. + - `steps_ahead::Int`: Number of steps ahead for forecasting. + - `N_scenarios::Int`: Number of scenarios to simulate. + - `Exogenous_Forecast::Matrix{Fl}`: Matrix of forecasts of exogenous variables. + - `dynamic_exog_coefs_forecasts::Vector{<:Vector}`: Vector of vectors of combination components forecasts. + - `seasonal_innovation_simulation::Int`: Number of seasonal innovation simulation. + - `seed::Int`: Seed for the random number generator. + + # Returns + - `Matrix{AbstractFloat}`: Matrix of scenarios of the states of the model. +""" +function simulate( + model::StructuralModel, + steps_ahead::Int, + N_scenarios::Int; + Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, 0), + dynamic_exog_coefs_forecasts::Vector{<:Vector}=Vector{Vector}(undef, 0), + seasonal_innovation_simulation::Int=0, + seed::Int=1234, +)::Matrix{AbstractFloat} where {Fl<:AbstractFloat} + scenarios = Matrix{AbstractFloat}(undef, steps_ahead, N_scenarios) + Random.seed!(seed) + for s in 1:N_scenarios + scenarios[:, s] = simulate_states( + model, steps_ahead, false, seasonal_innovation_simulation + ) end + + dynamic_exog_coefs_prediction = forecast_dynamic_exog_coefs( + model, steps_ahead, dynamic_exog_coefs_forecasts + ) + scenarios .+= + (Exogenous_Forecast * model.output.components["exog"]["Coefs"]) + + dynamic_exog_coefs_prediction + + return scenarios end isfitted(model::StructuralModel) = isnothing(model.output) ? false : true diff --git a/src/structs.jl b/src/structs.jl index 641b729..48a3e9b 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 + decomposition::Dict end diff --git a/src/utils.jl b/src/utils.jl index 93a170a..0511530 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -28,56 +28,12 @@ function build_components( components[key]["Values"] = X[:, components_indexes[key]] * coefs[components_indexes[key]] end - if haskey(components, "Exogenous_X") - components["Exogenous_X"]["Selected"] = findall( - i -> i != 0, components["Exogenous_X"]["Coefs"] - ) + if haskey(components, "exog") + components["exog"]["Selected"] = findall(i -> i != 0, components["exog"]["Coefs"]) end return components end -""" - 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}, @@ -117,52 +73,6 @@ function get_fit_and_residuals( 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 - """ handle_missing_values( X::Matrix{Tl}, y::Vector{Fl} @@ -191,42 +101,6 @@ function handle_missing_values( return y[valid_indexes], X[valid_indexes, :], valid_indexes end -""" -handle_missing_values( - X::Matrix{Tl}, y::Matrix{Fl} -)::Tuple{Matrix{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. - - `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 = Int[] - 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 @@ -243,107 +117,47 @@ function has_intercept(X::Matrix{Fl})::Bool where {Fl<:AbstractFloat} end """ -fill_innovation_coefs(model::StructuralModel, T::Int, component::String, valid_indexes::Vector{Int})::Vector{AbstractFloat} +get_stochastic_values(estimated_stochastic::Vector{Fl}, steps_ahead::Int, T::Int, start_idx::Int, final_idx::Int, seasonal_innovation_simulation::Int)::Vector{AbstractFloat} where {Fl<: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. + Generates stochastic seasonal values for a given time series. # Arguments - - `model::StructuralModel`: Structural model. - - `T::Int`: Length of the original time series. - - `component::String`: Component name. - - `valid_indexes::Vector{Int}`: Valid Indexes in the time series + - `estimated_stochastic::Vector{Fl}`: Vector of estimated stochastic terms. + - `steps_ahead::Int`: Number of steps ahead to generate. + - `T::Int`: Length of the time series. + - `start_idx::Int`: Starting index of the time series. + - `final_idx::Int`: Final index of the time series. + - `seasonal_innovation_simulation::Int`: Seasonal innovation simulation. # Returns - - `Union{Vector{AbstractFloat}, Matrix{AbstractFloat}}`: Vector or matrix containing innovation coefficients for the given component. + - `Vector{AbstractFloat}`: Vector of stochastic seasonal values. """ -function fill_innovation_coefs( - model::StructuralModel, component::String, valid_indexes::Vector{Int} -)::Union{Vector,Matrix} - T = length(model.y) - 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] +function get_stochastic_values( + estimated_stochastic::Vector{Fl}, + steps_ahead::Int, + T::Int, + start_idx::Int, + final_idx::Int, + seasonal_innovation_simulation::Int, +)::Vector{AbstractFloat} where {Fl<:AbstractFloat} + if seasonal_innovation_simulation != 0 + stochastic_term = Vector{AbstractFloat}(undef, steps_ahead) + for t in 1:steps_ahead + + # Generate potential seasonal indices + seasonal_indices = + ((T + t) % seasonal_innovation_simulation):seasonal_innovation_simulation:T + + # Filter indices to be within the valid range + valid_indices = filter(idx -> start_idx <= idx <= final_idx, seasonal_indices) + + # Sample with randomness and sign flip + stochastic_term[t] = rand(estimated_stochastic[valid_indices]) * rand([1, -1]) 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]) + stochastic_term = + rand(estimated_stochastic, steps_ahead) .* rand([1, -1], steps_ahead) 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 + return stochastic_term end diff --git a/test/estimation_procedure.jl b/test/estimation_procedure.jl index 842aebe..5084f47 100644 --- a/test/estimation_procedure.jl +++ b/test/estimation_procedure.jl @@ -57,33 +57,29 @@ end @testset "Function: fit_lasso" begin Random.seed!(1234) y = rand(10) - Exogenous_X = hcat(rand(10, 3), vcat(zeros(3), ones(1), zeros(6))) + exog = hcat(rand(10, 3), vcat(zeros(3), ones(1), zeros(6))) Basic_Structural = StateSpaceLearning.StructuralModel( y; - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=true, + level="stochastic", + slope="stochastic", + seasonal="stochastic", freq_seasonal=2, outlier=true, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X, + ζ_threshold=0, + ω_threshold=0, + exog=exog, ) Basic_Structural_w_level = StateSpaceLearning.StructuralModel( y; - level=false, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=true, + level="none", + slope="stochastic", + seasonal="stochastic", freq_seasonal=2, outlier=true, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X, + ζ_threshold=0, + ω_threshold=0, + exog=exog, ) components_indexes = StateSpaceLearning.get_components_indexes(Basic_Structural) @@ -107,6 +103,19 @@ end @test length(coefs0) == 43 @test length(ε0) == 10 + coefs0, ε0 = StateSpaceLearning.fit_lasso( + X, + AbstractFloat.(y), + 0.1, + "aic", + false, + components_indexes, + ones(size(X, 2) - 1); + rm_average=true, + ) + @test length(coefs0) == 43 + @test length(ε0) == 10 + coefs1, ε1 = StateSpaceLearning.fit_lasso( X2, AbstractFloat.(y), @@ -117,7 +126,7 @@ end ones(size(X2, 2)); rm_average=false, ) - @test length(coefs1) == 42 + @test length(coefs1) == 34 @test length(ε1) == 10 coefs2, ε2 = StateSpaceLearning.fit_lasso( @@ -145,7 +154,7 @@ end rm_average=true, ) @test coefs3[components_indexes["o"][4]] == 0 - @test all(coefs3[components_indexes["Exogenous_X"]] .!= 0) + @test all(coefs3[components_indexes["exog"]] .!= 0) @test length(coefs3) == 43 @test length(ε3) == 10 @@ -167,33 +176,29 @@ end @testset "Function: estimation_procedure" begin Random.seed!(1234) y = rand(10) - Exogenous_X = hcat(rand(10, 3), vcat(zeros(3), ones(1), zeros(6))) + exog = hcat(rand(10, 3), vcat(zeros(3), ones(1), zeros(6))) Basic_Structural = StateSpaceLearning.StructuralModel( y; - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=true, + level="stochastic", + slope="stochastic", + seasonal="stochastic", freq_seasonal=2, outlier=true, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X, + ζ_threshold=0, + ω_threshold=0, + exog=exog, ) Basic_Structural_w_level = StateSpaceLearning.StructuralModel( y; - level=false, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=true, + level="none", + slope="stochastic", + seasonal="stochastic", freq_seasonal=2, outlier=true, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X, + ζ_threshold=0, + ω_threshold=0, + exog=exog, ) components_indexes = StateSpaceLearning.get_components_indexes(Basic_Structural) @@ -204,20 +209,22 @@ end X = Basic_Structural.X X2 = Basic_Structural_w_level.X + innovations_names = StateSpaceLearning.get_model_innovations(Basic_Structural) + coefs1, ε1 = StateSpaceLearning.estimation_procedure( - X, y, components_indexes, 0.1, "aic", 0.05, true, true + X, y, components_indexes, 0.1, "aic", 0.05, true, true, innovations_names ) @test length(coefs1) == 43 @test length(ε1) == 10 coefs1, ε1 = StateSpaceLearning.estimation_procedure( - X2, y, components_indexes2, 0.1, "aic", 0.05, true, true + X2, y, components_indexes2, 0.1, "aic", 0.05, true, true, innovations_names ) - @test length(coefs1) == 42 + @test length(coefs1) == 34 @test length(ε1) == 10 coefs2, ε2 = StateSpaceLearning.estimation_procedure( - X, y, components_indexes, 0.1, "aic", 0.05, true, false + X, y, components_indexes, 0.1, "aic", 0.05, true, false, innovations_names ) @test length(coefs2) == 43 @test length(ε2) == 10 @@ -225,47 +232,43 @@ end end @testset "Function: get_dummy_indexes" begin - Exogenous_X1 = hcat(rand(10, 3), vcat(zeros(3), ones(1), zeros(6))) - Exogenous_X2 = hcat(rand(10, 3)) + exog1 = hcat(rand(10, 3), vcat(zeros(3), ones(1), zeros(6))) + exog2 = hcat(rand(10, 3)) - dummy_indexes1 = StateSpaceLearning.get_dummy_indexes(Exogenous_X1) + dummy_indexes1 = StateSpaceLearning.get_dummy_indexes(exog1) @test dummy_indexes1 == [4] - dummy_indexes2 = StateSpaceLearning.get_dummy_indexes(Exogenous_X2) + dummy_indexes2 = StateSpaceLearning.get_dummy_indexes(exog2) @test dummy_indexes2 == [] end @testset "Function: get_outlier_duplicate_columns" begin Random.seed!(1234) y = rand(10) - Exogenous_X = hcat(rand(10, 3), vcat(zeros(3), ones(1), zeros(6))) - Exogenous_X2 = rand(10, 3) + exog = hcat(rand(10, 3), vcat(zeros(3), ones(1), zeros(6))) + exog2 = rand(10, 3) Basic_Structural = StateSpaceLearning.StructuralModel( y; - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=true, + level="stochastic", + slope="stochastic", + seasonal="stochastic", freq_seasonal=2, outlier=true, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X, + ζ_threshold=0, + ω_threshold=0, + exog=exog, ) Basic_Structural_w_out = StateSpaceLearning.StructuralModel( y; - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=true, + level="stochastic", + slope="stochastic", + seasonal="stochastic", freq_seasonal=2, outlier=true, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X2, + ζ_threshold=0, + ω_threshold=0, + exog=exog2, ) components_indexes = StateSpaceLearning.get_components_indexes(Basic_Structural) diff --git a/test/fit_forecast.jl b/test/fit_forecast.jl index 06827c9..4a8993e 100644 --- a/test/fit_forecast.jl +++ b/test/fit_forecast.jl @@ -1,4 +1,4 @@ -@testset "Function: fit_model" begin +@testset "Function: fit!" begin y1 = rand(100) y2 = rand(100) y2[10:20] .= NaN @@ -11,7 +11,8 @@ @test length(model1.output.coefs) == 375 @test length(model1.output.valid_indexes) == 100 @test length(model1.output.residuals_variances) == 4 - @test length(keys(model1.output.components)) == 9 + @test length(keys(model1.output.components)) == 10 + @test length(keys(model1.output.decomposition)) == 3 model2 = StateSpaceLearning.StructuralModel(y2) StateSpaceLearning.fit!(model2) @@ -21,256 +22,6 @@ @test length(model2.output.coefs) == 375 @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 - y1 = rand(100) - y2 = rand(100) - y2[10:20] .= NaN - - model1 = StateSpaceLearning.StructuralModel(y1) - StateSpaceLearning.fit!(model1) - @test length(StateSpaceLearning.forecast(model1, 10)) == 10 - - model2 = StateSpaceLearning.StructuralModel(y2; Exogenous_X=rand(100, 3)) - StateSpaceLearning.fit!(model2) - @test length(StateSpaceLearning.forecast(model2, 10; Exogenous_Forecast=rand(10, 3))) == - 10 - - @test_throws AssertionError StateSpaceLearning.forecast( - model1, 10; Exogenous_Forecast=rand(5, 3) - ) - @test_throws AssertionError StateSpaceLearning.forecast(model2, 10) - @test_throws AssertionError StateSpaceLearning.forecast( - model2, 10; Exogenous_Forecast=rand(5, 3) - ) - - y3 = [ - 4.718, - 4.77, - 4.882, - 4.859, - 4.795, - 4.905, - 4.997, - 4.997, - 4.912, - 4.779, - 4.644, - 4.77, - 4.744, - 4.836, - 4.948, - 4.905, - 4.828, - 5.003, - 5.135, - 5.135, - 5.062, - 4.89, - 4.736, - 4.941, - 4.976, - 5.01, - 5.181, - 5.093, - 5.147, - 5.181, - 5.293, - 5.293, - 5.214, - 5.087, - 4.983, - 5.111, - 5.141, - 5.192, - 5.262, - 5.198, - 5.209, - 5.384, - 5.438, - 5.488, - 5.342, - 5.252, - 5.147, - 5.267, - 5.278, - 5.278, - 5.463, - 5.459, - 5.433, - 5.493, - 5.575, - 5.605, - 5.468, - 5.351, - 5.192, - 5.303, - 5.318, - 5.236, - 5.459, - 5.424, - 5.455, - 5.575, - 5.71, - 5.68, - 5.556, - 5.433, - 5.313, - 5.433, - 5.488, - 5.451, - 5.587, - 5.594, - 5.598, - 5.752, - 5.897, - 5.849, - 5.743, - 5.613, - 5.468, - 5.627, - 5.648, - 5.624, - 5.758, - 5.746, - 5.762, - 5.924, - 6.023, - 6.003, - 5.872, - 5.723, - 5.602, - 5.723, - 5.752, - 5.707, - 5.874, - 5.852, - 5.872, - 6.045, - 6.142, - 6.146, - 6.001, - 5.849, - 5.72, - 5.817, - 5.828, - 5.762, - 5.891, - 5.852, - 5.894, - 6.075, - 6.196, - 6.224, - 6.001, - 5.883, - 5.736, - 5.82, - 5.886, - 5.834, - 6.006, - 5.981, - 6.04, - 6.156, - 6.306, - 6.326, - 6.137, - 6.008, - 5.891, - 6.003, - 6.033, - 5.968, - 6.037, - 6.133, - 6.156, - 6.282, - 6.432, - 6.406, - 6.23, - 6.133, - 5.966, - 6.068, - ] - model3 = StateSpaceLearning.StructuralModel(y3) - StateSpaceLearning.fit!(model3) - forecast3 = trunc.(StateSpaceLearning.forecast(model3, 18); digits=3) - @assert forecast3 == [ - 6.11, - 6.082, - 6.221, - 6.19, - 6.197, - 6.328, - 6.447, - 6.44, - 6.285, - 6.163, - 6.026, - 6.142, - 6.166, - 6.138, - 6.278, - 6.246, - 6.253, - 6.384, - ] - - model4 = StateSpaceLearning.StructuralModel(y3; freq_seasonal=[12, 36]) - StateSpaceLearning.fit!(model4) - forecast4 = trunc.(StateSpaceLearning.forecast(model4, 18); digits=3) - - @test length(forecast4) == 18 -end - -@testset "Function: simulate" begin - y1 = rand(100) - y2 = rand(100) - y2[10:20] .= NaN - - model1 = StateSpaceLearning.StructuralModel(y1) - 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 + @test length(keys(model2.output.components)) == 10 + @test length(keys(model2.output.decomposition)) == 3 end diff --git a/test/models/structural_model.jl b/test/models/structural_model.jl index f692c54..becdbe4 100644 --- a/test/models/structural_model.jl +++ b/test/models/structural_model.jl @@ -4,9 +4,7 @@ model1 = StateSpaceLearning.StructuralModel(y1) 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], dumping_cycle=0.5 - ) + model4 = StateSpaceLearning.StructuralModel(y1; cycle_period=[3, 10.2]) @test typeof(model1) == StateSpaceLearning.StructuralModel @test typeof(model2) == StateSpaceLearning.StructuralModel @@ -20,76 +18,28 @@ @test_throws AssertionError StateSpaceLearning.StructuralModel(y1; freq_seasonal=1000) - @test_throws AssertionError StateSpaceLearning.StructuralModel( - y1; cycle_period=[3, 10.2], dumping_cycle=-1.0 - ) - @test_throws AssertionError StateSpaceLearning.StructuralModel( - y1; cycle_period=[3, 10.2], dumping_cycle=2.0 - ) - exog_error = ones(100, 3) - @test_throws AssertionError StateSpaceLearning.StructuralModel( - y1; Exogenous_X=exog_error - ) + @test_throws AssertionError StateSpaceLearning.StructuralModel(y1; exog=exog_error) end -@testset "create_deterministic_cycle_matrix" begin - cycle_matrix = Vector{Matrix}(undef, 1) - A = 1 * cos(2 * pi / 12) - B = 1 * sin(2 * pi / 12) - cycle_matrix[1] = [A B; -B A] - det_matrix1 = StateSpaceLearning.create_deterministic_cycle_matrix(cycle_matrix, 5, 0) - exp_det_matrix1 = [ - 1.0 0.0 - 0.866025 0.5 - 0.5 0.866025 - 2.77556e-16 1.0 - -0.5 0.866025 - ] +@testset "create deterministic matrices" begin + X1 = StateSpaceLearning.create_deterministic_cycle(100, 12) + @test size(X1) == (100, 2) - n, p = size(det_matrix1[1]) - for i in 1:n - for j in 1:p - @test isapprox(det_matrix1[1][i, j], exp_det_matrix1[i, j]; atol=1e-6) - end - end + X2 = StateSpaceLearning.create_deterministic_seasonal(100, 12) + @test size(X2) == (100, 12) - cycle_period = [3, 12.2] - cycle_matrix = Vector{Matrix}(undef, length(cycle_period)) - for i in eachindex(cycle_period) - A = 1 * cos(2 * pi / cycle_period[i]) - B = 1 * sin(2 * pi / cycle_period[i]) - cycle_matrix[i] = [A B; -B A] - end + X3 = StateSpaceLearning.create_initial_states_Matrix( + 100, [12, 20], true, true, true, true, [3, 10.2] + ) - #### - - det_matrix2 = StateSpaceLearning.create_deterministic_cycle_matrix(cycle_matrix, 5, 0) - exp_det_matrix2 = [ - [ - 1.0 0.0 - -0.5 0.866025 - -0.5 -0.866025 - 1.0 -6.10623e-16 - -0.5 0.866025 - ], - [ - 1.0 0.0 - 0.870285 0.492548 - 0.514793 0.857315 - 0.0257479 0.999668 - -0.469977 0.882679 - ], - ] + @test size(X3) == (100, 38) - n, p = 5, 2 - for h in eachindex(det_matrix2) - for i in 1:n - for j in 1:p - @test isapprox(det_matrix2[h][i, j], exp_det_matrix2[h][i, j]; atol=1e-6) - end - end - end + X4 = StateSpaceLearning.create_initial_states_Matrix( + 100, 12, true, true, true, false, 3 + ) + + @test size(X4) == (100, 14) end @testset "Innovation matrices" begin @@ -99,8 +49,8 @@ end @test StateSpaceLearning.ω_size(10, 2, 0, 1) == 9 @test StateSpaceLearning.ω_size(10, 2, 2, 1) == 7 @test StateSpaceLearning.o_size(10, 1) == 10 - @test StateSpaceLearning.ϕ_size(10, 0, 1) == 14 - @test StateSpaceLearning.ϕ_size(10, 2, 1) == 12 + @test StateSpaceLearning.ϕ_size(10, 0, 1) == 16 + @test StateSpaceLearning.ϕ_size(10, 2, 1) == 14 @test StateSpaceLearning.ξ_size(10, 5) == 5 @test StateSpaceLearning.ζ_size(10, 2, 5) == 3 @@ -110,8 +60,8 @@ end @test StateSpaceLearning.o_size(10, 6) == 5 @test StateSpaceLearning.ϕ_size(10, 0, 5) == 10 - X_ξ1 = StateSpaceLearning.create_ξ(5, 0, 1) - X_ξ2 = StateSpaceLearning.create_ξ(5, 2, 1) + X_ξ1 = StateSpaceLearning.create_ξ(5, 1) + X_ξ2 = StateSpaceLearning.create_ξ(5, 3) @test X_ξ1 == [ 0.0 0.0 0.0 @@ -122,19 +72,6 @@ end ] @test X_ξ2 == [ - 0.0 0.0 0.0 - 1.0 0.0 0.0 - 1.0 1.0 0.0 - 1.0 1.0 1.0 - 1.0 1.0 1.0 - 1.0 1.0 1.0 - 1.0 1.0 1.0 - ] - - X_ξ3 = StateSpaceLearning.create_ξ(5, 0, 3) - X_ξ4 = StateSpaceLearning.create_ξ(5, 2, 3) - - @test X_ξ3 == [ 0.0 0.0 0.0 0.0 1.0 0.0 @@ -142,19 +79,9 @@ end 1.0 1.0 ] - @test X_ξ4 == [ - 0.0 0.0 - 0.0 0.0 - 1.0 0.0 - 1.0 1.0 - 1.0 1.0 - 1.0 1.0 - 1.0 1.0 - ] - - X_ζ1 = StateSpaceLearning.create_ζ(5, 0, 0, 1) - X_ζ2 = StateSpaceLearning.create_ζ(5, 2, 0, 1) - X_ζ3 = StateSpaceLearning.create_ζ(5, 2, 2, 1) + X_ζ1 = StateSpaceLearning.create_ζ(5, 0, 1) + X_ζ2 = StateSpaceLearning.create_ζ(6, 2, 1) + X_ζ3 = StateSpaceLearning.create_ζ(5, 2, 3) @test X_ζ1 == [ 0.0 0.0 0.0 @@ -165,48 +92,19 @@ end ] @test X_ζ2 == [ - 0.0 0.0 0.0 - 0.0 0.0 0.0 - 1.0 0.0 0.0 - 2.0 1.0 0.0 - 3.0 2.0 1.0 - 4.0 3.0 2.0 - 5.0 4.0 3.0 + 0.0 0.0 + 0.0 0.0 + 1.0 0.0 + 2.0 1.0 + 3.0 2.0 + 4.0 3.0 ] - @test X_ζ3 == reshape( - [ - 0.0 - 0.0 - 1.0 - 2.0 - 3.0 - 4.0 - 5.0 - ], - 7, - 1, - ) - - X_ζ4 = StateSpaceLearning.create_ζ(6, 2, 2, 3) - - @test X_ζ4 == reshape( - [ - 0.0 - 0.0 - 0.0 - 1.0 - 2.0 - 3.0 - 4.0 - 5.0 - ], - 8, - 1, - ) + @test X_ζ3 == zeros(5, 0) - X_ω1 = StateSpaceLearning.create_ω(5, 2, 0, 0, 1) - X_ω2 = StateSpaceLearning.create_ω(5, 2, 2, 0, 1) + X_ω1 = StateSpaceLearning.create_ω(5, 2, 0, 1) + X_ω2 = StateSpaceLearning.create_ω(5, 2, 2, 1) + X_ω3 = StateSpaceLearning.create_ω(5, 2, 0, 3) @test X_ω1 == [ 0.0 0.0 0.0 0.0 @@ -217,17 +115,13 @@ end ] @test X_ω2 == [ - 0.0 0.0 0.0 0.0 - 0.0 0.0 0.0 0.0 - -1.0 1.0 0.0 0.0 - 0.0 -1.0 1.0 0.0 - -1.0 1.0 -1.0 1.0 - 0.0 -1.0 1.0 -1.0 - -1.0 1.0 -1.0 1.0 + 0.0 0.0 + 0.0 0.0 + -1.0 1.0 + 0.0 -1.0 + -1.0 1.0 ] - X_ω3 = StateSpaceLearning.create_ω(5, 2, 0, 0, 3) - @test X_ω3 == [ 0.0 0.0 0.0 0.0 0.0 0.0 @@ -236,341 +130,196 @@ end 1.0 -1.0 1.0 ] - X_o1 = StateSpaceLearning.create_o_matrix(3, 0, 1) - X_o2 = StateSpaceLearning.create_o_matrix(3, 2, 1) - X_o3 = StateSpaceLearning.create_o_matrix(3, 0, 2) + X_o1 = StateSpaceLearning.create_o_matrix(3, 1) + X_o3 = StateSpaceLearning.create_o_matrix(3, 2) @test X_o1 == Matrix(1.0 * I, 3, 3) - @test X_o2 == vcat(Matrix(1.0 * I, 3, 3), zeros(2, 3)) - @test X_o3 == Matrix(1.0 * I, 3, 3)[:, 2:3] - - cycle_matrix = Vector{Matrix}(undef, 1) - A = 1 * cos(2 * pi / 12) - B = 1 * sin(2 * pi / 12) - cycle_matrix[1] = [A B; -B A] - det_matrix1 = StateSpaceLearning.create_deterministic_cycle_matrix(cycle_matrix, 5, 0) - X_ϕ1 = StateSpaceLearning.create_ϕ(det_matrix1[1], 5, 0, 0, 1) - - @test all( - isapprox.( - X_ϕ1, - [ - 0.0 0.0 0.0 0.0 0.0 0.0 - 1.0 0.0 0.0 0.0 0.0 0.0 - 0.866025 0.5 1.0 0.0 0.0 0.0 - 0.5 0.866025 0.866025 0.5 1.0 0.0 - 2.77556e-16 1.0 0.5 0.866025 0.866025 0.5 - ], - atol=1e-6, - ), - ) -end + @test X_o3 == [ + 0.0 0.0 + 1.0 0.0 + 0.0 1.0 + ] -@testset "Initial State Matrix" begin - X1 = StateSpaceLearning.create_initial_states_Matrix( - 5, 2, 0, true, true, true, Vector{Matrix}(undef, 0) - ) - X2 = StateSpaceLearning.create_initial_states_Matrix( - 5, 2, 2, true, true, true, Vector{Matrix}(undef, 0) - ) + X_ϕ1 = StateSpaceLearning.create_ϕ(3, 5, 0, 1) + X_ϕ2 = StateSpaceLearning.create_ϕ(3, 5, 3, 1) + X_ϕ3 = StateSpaceLearning.create_ϕ(3, 5, 0, 2) - @test X1 == [ - 1.0 0.0 1.0 0.0 - 1.0 1.0 0.0 1.0 - 1.0 2.0 1.0 0.0 - 1.0 3.0 0.0 1.0 - 1.0 4.0 1.0 0.0 + @test X_ϕ1 == [ + 0.0 0.0 0.0 0.0 0.0 0.0 + -0.5 -0.86603 0.0 0.0 0.0 0.0 + 1.0 -0.0 1.0 -0.0 0.0 0.0 + -0.5 0.86603 -0.5 0.86603 -0.5 0.86603 + -0.5 -0.86603 -0.5 -0.86603 -0.5 -0.86603 ] - @test X2 == [ - 1.0 0.0 1.0 0.0 - 1.0 1.0 0.0 1.0 - 1.0 2.0 1.0 0.0 - 1.0 3.0 0.0 1.0 - 1.0 4.0 1.0 0.0 - 1.0 5.0 0.0 1.0 - 1.0 6.0 1.0 0.0 + @test X_ϕ2 == [ + 0.0 0.0 + -0.5 -0.86603 + 1.0 -0.0 + -0.5 0.86603 + -0.5 -0.86603 ] - X3 = StateSpaceLearning.create_initial_states_Matrix( - 5, 2, 0, true, true, false, Vector{Matrix}(undef, 0) - ) - X4 = StateSpaceLearning.create_initial_states_Matrix( - 5, 2, 2, true, true, false, Vector{Matrix}(undef, 0) - ) - - @test X3 == [ - 1.0 0.0 - 1.0 1.0 - 1.0 2.0 - 1.0 3.0 - 1.0 4.0 + @test X_ϕ3 == [ + 0.0 0.0 0.0 0.0 0.0 0.0 + -0.5 -0.86603 0.0 0.0 0.0 0.0 + 1.0 -0.0 1.0 -0.0 0.0 0.0 + -0.5 0.86603 -0.5 0.86603 -0.5 0.86603 + -0.5 -0.86603 -0.5 -0.86603 -0.5 -0.86603 ] +end - @test X4 == [ - 1.0 0.0 - 1.0 1.0 - 1.0 2.0 - 1.0 3.0 - 1.0 4.0 - 1.0 5.0 - 1.0 6.0 - ] +@testset "dynamic_exog_coefs" begin + dynamic_exog_coefs1 = [(collect(1:5), "level")] - X5 = StateSpaceLearning.create_initial_states_Matrix( - 5, 2, 0, true, false, false, Vector{Matrix}(undef, 0) - ) - X6 = StateSpaceLearning.create_initial_states_Matrix( - 5, 2, 2, true, false, false, Vector{Matrix}(undef, 0) + X1 = StateSpaceLearning.create_dynamic_exog_coefs_matrix( + dynamic_exog_coefs1, 5, 0, 0, 0, 1 ) + @test size(X1) == (5, 4) - @test X5 == ones(5, 1) - @test X6 == ones(7, 1) + dynamic_exog_coefs2 = [(collect(1:5), "slope")] - X7 = StateSpaceLearning.create_initial_states_Matrix( - 5, 2, 0, false, true, false, Vector{Matrix}(undef, 0) + X2 = StateSpaceLearning.create_dynamic_exog_coefs_matrix( + dynamic_exog_coefs2, 5, 0, 0, 0, 1 ) - X8 = StateSpaceLearning.create_initial_states_Matrix( - 5, 2, 2, false, true, false, Vector{Matrix}(undef, 0) - ) - - @test X7 == [0.0; 1.0; 2.0; 3.0; 4.0][:, :] - @test X8 == [0.0; 1.0; 2.0; 3.0; 4.0; 5.0; 6.0][:, :] + @test size(X2) == (5, 4) - cycle_matrix = Vector{Matrix}(undef, 1) - A = 1 * cos(2 * pi / 12) - B = 1 * sin(2 * pi / 12) - cycle_matrix[1] = [A B; -B A] - det_matrix1 = StateSpaceLearning.create_deterministic_cycle_matrix(cycle_matrix, 5, 0) - det_matrix2 = StateSpaceLearning.create_deterministic_cycle_matrix(cycle_matrix, 5, 2) + dynamic_exog_coefs = [ + (collect(1:5), "level"), + (collect(1:5), "slope"), + (collect(1:5), "seasonal", 2), + (collect(1:5), "cycle", 3), + ] - X9 = StateSpaceLearning.create_initial_states_Matrix( - 5, 2, 0, true, true, true, det_matrix1 + X = StateSpaceLearning.create_dynamic_exog_coefs_matrix( + dynamic_exog_coefs, 5, 0, 0, 0, 1 ) - X10 = StateSpaceLearning.create_initial_states_Matrix( - 5, 2, 2, true, true, true, det_matrix2 + @test size(X) == (5, 22) + + dynamic_exog_coefs = [ + (collect(6:7), "level", "", 4), + (collect(6:7), "slope", "", 4), + (collect(6:7), "seasonal", 2, 7), + (collect(6:7), "cycle", 3, 8), + ] + X_f = StateSpaceLearning.create_forecast_dynamic_exog_coefs_matrix( + dynamic_exog_coefs, 5, 2, 0, 0, 0, 1 ) + @test size(X_f) == (2, 23) - @test all( - isapprox.( - X9, - [ - 1.0 0.0 1.0 0.0 1.0 0.0 - 1.0 1.0 0.0 1.0 0.866025 0.5 - 1.0 2.0 1.0 0.0 0.5 0.866025 - 1.0 3.0 0.0 1.0 2.77556e-16 1.0 - 1.0 4.0 1.0 0.0 -0.5 0.866025 - ], - atol=1e-6, - ), + 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 ) + @test size(X_f2) == (2, 4) - @test all( - isapprox.( - X10, - [ - 1.0 0.0 1.0 0.0 1.0 0.0 - 1.0 1.0 0.0 1.0 0.866025 0.5 - 1.0 2.0 1.0 0.0 0.5 0.866025 - 1.0 3.0 0.0 1.0 2.77556e-16 1.0 - 1.0 4.0 1.0 0.0 -0.5 0.866025 - 1.0 5.0 0.0 1.0 -0.866025 0.5 - 1.0 6.0 1.0 0.0 -1.0 5.55112e-16 - ], - atol=1e-6, - ), + 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 ) + @test size(X_f3) == (2, 4) end @testset "Create X matrix" begin - Exogenous_X1 = rand(5, 3) - Exogenous_X2 = zeros(5, 0) - - Exogenous_forecast1 = rand(2, 3) - - cycle_matrix = Vector{Matrix}(undef, 0) - stochastic_cycle = false - stochastic_start = 1 - - param_combination = [ - Any[true, 0, 0], - Any[true, 2, 0], - Any[true, 2, 2], - Any[false, 0, 0], - Any[false, 2, 0], - Any[false, 2, 2], + exog1 = rand(5, 3) + dynamic_exog_coefs = [ + (collect(1:5), "level"), + (collect(1:5), "slope"), + (collect(1:5), "seasonal", 2), + (collect(1:5), "cycle", 3), ] - size_vec1 = [ - (5, 22), - (5, 18), - (7, 18), - (5, 17), - (5, 13), - (7, 13), - (5, 12), - (5, 12), - (7, 12), - (5, 7), - (5, 7), - (7, 7), - (5, 16), - (5, 14), - (7, 14), - (5, 11), - (5, 9), - (7, 9), - (5, 13), - (5, 11), - (7, 11), - (5, 8), - (5, 6), - (7, 6), - ] - size_vec2 = [ - (5, 19), - (5, 15), - (7, 15), - (5, 14), - (5, 10), - (7, 10), - (5, 9), - (5, 9), - (7, 9), - (5, 4), - (5, 4), - (7, 4), - (5, 13), - (5, 11), - (7, 11), - (5, 8), - (5, 6), - (7, 6), - (5, 10), - (5, 8), - (7, 8), - (5, 5), - (5, 3), - (7, 3), - ] - counter = 1 - for args in [ - [true, true, true, true, true, true, 2, cycle_matrix, stochastic_cycle, true, 12], - [ - true, - true, - false, - false, - false, - false, - 2, - cycle_matrix, - stochastic_cycle, - true, - 12, - ], - [true, true, true, true, false, false, 2, cycle_matrix, stochastic_cycle, true, 12], - [ - true, - false, - true, - true, - false, - false, - 2, - cycle_matrix, - stochastic_cycle, - true, - 12, - ], + X1 = StateSpaceLearning.create_X( + true, + true, + true, + true, + true, + true, + true, + true, + 3, + 3, + true, + 0, + 0, + 0, + 1, + exog1, + dynamic_exog_coefs, + ) + + exog2 = zeros(10, 3) + dynamic_exog_coefs2 = [ + (collect(1:10), "level"), + (collect(1:10), "slope"), + (collect(1:10), "seasonal", 2), + (collect(1:10), "cycle", 3), ] - args = [x in [0, 1] ? Bool(x) : x for x in args] - for param in param_combination - args[end - 1] = param[1] - args[end] = param[2] - if param[3] != 0 - X1 = StateSpaceLearning.create_X( - args..., stochastic_start, Exogenous_X1, param[3], Exogenous_forecast1 - ) - else - X1 = StateSpaceLearning.create_X( - args..., stochastic_start, Exogenous_X1, param[3] - ) - end - X2 = StateSpaceLearning.create_X( - args..., stochastic_start, Exogenous_X2, param[3] - ) - @test size(X1) == size_vec1[counter] - @test size(X2) == size_vec2[counter] - counter += 1 - end - end + + X2 = StateSpaceLearning.create_X( + true, + true, + true, + true, + true, + true, + true, + true, + 3, + 3, + true, + 3, + 2, + 4, + 1, + exog2, + dynamic_exog_coefs2, + ) + + @test size(X1) == (5, 52) + @test size(X2) == (10, 85) end @testset "Function: get_components" begin - Exogenous_X1 = rand(10, 3) - Exogenous_X2 = zeros(10, 0) + exog1 = rand(10, 3) + exog2 = zeros(10, 0) Basic_Structural = StateSpaceLearning.StructuralModel( - rand(10); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=true, - freq_seasonal=2, - outlier=true, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X1, + rand(10); freq_seasonal=2, outlier=true, ζ_threshold=0, ω_threshold=0, exog=exog1 ) Local_Level = StateSpaceLearning.StructuralModel( rand(10); - level=true, - stochastic_level=true, - trend=false, - stochastic_trend=false, - seasonal=false, - stochastic_seasonal=false, + slope="none", + seasonal="none", + cycle="none", freq_seasonal=2, outlier=true, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X1, + ζ_threshold=0, + exog=exog1, ) Local_Linear_Trend1 = StateSpaceLearning.StructuralModel( - rand(10); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=false, - stochastic_seasonal=false, - freq_seasonal=2, - outlier=false, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X1, + rand(10); seasonal="none", cycle="none", outlier=false, ζ_threshold=0, exog=exog1 ) Local_Linear_Trend2 = StateSpaceLearning.StructuralModel( rand(10); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=false, - stochastic_seasonal=false, + seasonal="none", + cycle="none", freq_seasonal=2, outlier=false, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X2, + ζ_threshold=0, + exog=exog2, ) models = [Basic_Structural, Local_Level, Local_Linear_Trend1, Local_Linear_Trend2] empty_keys_vec = [ - [], ["ν1", "ζ", "γ₁", "ω_2"], ["γ₁", "ω_2", "o"], ["γ₁", "ω_2", "o", "Exogenous_X"] + [], ["ν1", "ζ", "γ₁", "ω_2"], ["γ₁", "ω_2", "o"], ["γ₁", "ω_2", "o", "exog"] ] - exogs = [Exogenous_X1, Exogenous_X1, Exogenous_X1, Exogenous_X2] + exogs = [exog1, exog1, exog1, exog2] for idx in eachindex(models) model = models[idx] @@ -579,7 +328,7 @@ end for key in keys(components_indexes) @test (key in empty_keys_vec[idx]) ? isempty(components_indexes[key]) : true - @test if key == "Exogenous_X" + @test if key == "exog" length(components_indexes[key]) == size(exogs[idx], 2) else true @@ -588,75 +337,42 @@ end end end -@testset "Function: get_variances" begin - Exogenous_X2 = zeros(10, 0) +@testset "Functions: get_variances and get_model_innovations" begin + exog2 = zeros(10, 0) Basic_Structural = StateSpaceLearning.StructuralModel( - rand(10); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=true, - freq_seasonal=2, - outlier=true, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X2, + rand(10); freq_seasonal=2, outlier=true, ζ_threshold=0, exog=exog2 ) Basic_Structural2 = StateSpaceLearning.StructuralModel( - rand(10); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=true, - freq_seasonal=[2, 5], - outlier=true, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X2, + rand(10); freq_seasonal=[2, 5], outlier=true, ζ_threshold=0, exog=exog2 ) Local_Level = StateSpaceLearning.StructuralModel( rand(10); - level=true, - stochastic_level=true, - trend=false, - stochastic_trend=false, - seasonal=false, - stochastic_seasonal=false, + slope="none", + seasonal="none", freq_seasonal=2, outlier=true, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X2, + ζ_threshold=0, + exog=exog2, ) Local_Linear_Trend = StateSpaceLearning.StructuralModel( rand(10); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=false, - stochastic_seasonal=false, + seasonal="none", + cycle="none", freq_seasonal=2, outlier=false, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X2, + ζ_threshold=0, + exog=exog2, ) Local_Linear_Trend_cycle = StateSpaceLearning.StructuralModel( rand(10); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=false, - stochastic_seasonal=false, + seasonal="none", + cycle="stochastic", freq_seasonal=2, cycle_period=3, - stochastic_cycle=true, outlier=false, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X2, + ζ_threshold=0, + exog=exog2, ) models = [ @@ -672,7 +388,11 @@ end ["ξ", "ζ", "ω_2", "ω_5", "ε"], ["ξ", "ε"], ["ξ", "ζ", "ε"], - ["ξ", "ζ", "ε", "ϕ_1"], + ["ξ", "ζ", "ε", "ϕ_3"], + ] + + models_innovations = [ + ["ξ", "ζ", "ω_2"], ["ξ", "ζ", "ω_2", "ω_5"], ["ξ"], ["ξ", "ζ"], ["ξ", "ζ", "ϕ_3"] ] for idx in eachindex(models) @@ -682,207 +402,451 @@ end model, rand(100), rand(100), components_indexes ) @test all([key in keys(variances) for key in params_vec[idx]]) + model_innovations = StateSpaceLearning.get_model_innovations(model) + @test model_innovations == models_innovations[idx] end end -@testset "Function: get_model_innovations" begin - model1 = StateSpaceLearning.StructuralModel( - rand(10); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=true, - freq_seasonal=2, - outlier=true, - ζ_ω_threshold=0, - Exogenous_X=zeros(10, 0), - ) - model2 = StateSpaceLearning.StructuralModel( - rand(10); - level=true, - stochastic_level=false, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=true, - freq_seasonal=2, - outlier=true, - ζ_ω_threshold=0, - Exogenous_X=zeros(10, 0), - ) - model3 = StateSpaceLearning.StructuralModel( - rand(10); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=false, - seasonal=true, - stochastic_seasonal=true, - freq_seasonal=2, - outlier=true, - ζ_ω_threshold=0, - Exogenous_X=zeros(10, 0), +@testset "Decomposion functions" begin + Random.seed!(123) + model = StateSpaceLearning.StructuralModel( + vcat(collect(1:5), collect(5:-1:1)); + level="deterministic", + seasonal="none", + cycle="none", + outlier=false, + slope="stochastic", + ζ_threshold=0, ) - model4 = StateSpaceLearning.StructuralModel( - rand(10); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=false, - freq_seasonal=2, - outlier=true, - ζ_ω_threshold=0, - Exogenous_X=zeros(10, 0), + StateSpaceLearning.fit!(model) + slope = StateSpaceLearning.get_slope_decomposition(model, model.output.components) + @test length(slope) == 10 + + model = StateSpaceLearning.StructuralModel( + vcat(rand(5) .+ 5, rand(5) .- 5) + vcat(collect(1:5), collect(5:-1:1)); + seasonal="none", + cycle="none", + outlier=false, + slope="stochastic", + ζ_threshold=0, ) - model5 = StateSpaceLearning.StructuralModel( - rand(10); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=false, - seasonal=true, - stochastic_seasonal=true, - freq_seasonal=[2, 5], - outlier=true, - ζ_ω_threshold=0, - Exogenous_X=zeros(10, 0), + StateSpaceLearning.fit!(model) + trend = StateSpaceLearning.get_trend_decomposition( + model, model.output.components, slope ) - model6 = StateSpaceLearning.StructuralModel( + @test length(trend) == 10 + + model = StateSpaceLearning.StructuralModel( rand(10); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=false, - seasonal=true, - stochastic_seasonal=true, - freq_seasonal=[2, 5], + cycle="stochastic", cycle_period=3, - stochastic_cycle=true, - outlier=true, - ζ_ω_threshold=0, - Exogenous_X=zeros(10, 0), + outlier=false, + slope="stochastic", + ζ_threshold=0, + freq_seasonal=3, + ω_threshold=0, + ϕ_threshold=0, + ) + StateSpaceLearning.fit!(model) + @test all( + isapprox.( + StateSpaceLearning.get_seasonal_decomposition( + model, model.output.components, 3 + ), + [ + -0.011114430313782316 + 0.016993897901375513 + 0.0 + -0.06137711460224293 + -0.028221145277986203 + 0.045215043179361716 + -0.0027753371516487588 + -0.08682292272858037 + -0.036923166352424444 + 0.13243351808078532 + ], + atol=1e-6, + ), ) - models = [model1, model2, model3, model4] + @test all( + isapprox.( + StateSpaceLearning.get_cycle_decomposition(model, model.output.components, 3), + [ + 0.0 + 0.0 + 1.6111635465327112e-18 + -0.005696779520104198 + -0.030765036256794644 + 0.08224069806471179 + 0.030974382058151183 + -0.15828117942894446 + 0.037361403241860734 + 0.17009485813575637 + ], + atol=1e-6, + ), + ) - keys_vec = [ - ["ξ", "ζ", "ω_2"], - ["ζ", "ω_2"], - ["ξ", "ω_2"], - ["ξ", "ζ"], - ["ξ", "ω_2", "ω_5", "ϕ_1"], - ] + model_decomposition = StateSpaceLearning.get_model_decomposition( + model, model.output.components + ) + @test sort(collect(keys(model_decomposition))) == + sort(["cycle_3", "cycle_hat_3", "seasonal_3", "slope", "trend"]) - for idx in eachindex(models) - model = models[idx] - model_innovations = StateSpaceLearning.get_model_innovations(model) - @test model_innovations == keys_vec[idx] - end + model = StateSpaceLearning.StructuralModel( + vcat(rand(5) .+ 5, rand(5) .- 5) + vcat(collect(1:5), collect(5:-1:1)); + level="deterministic", + seasonal="none", + cycle="none", + outlier=false, + slope="stochastic", + ζ_threshold=0, + ) + StateSpaceLearning.fit!(model) + trend = StateSpaceLearning.get_trend_decomposition( + model, model.output.components, slope + ) + @test length(trend) == 10 end -@testset "Function: get_innovation_simulation_X" begin - innovation1 = "ξ" - innovation2 = "ζ" - innovation3 = "ω_2" - innovation4 = "ϕ_1" +@testset "Function: simulate_states" begin + model = StateSpaceLearning.StructuralModel(rand(100)) + StateSpaceLearning.fit!(model) + @test length(StateSpaceLearning.simulate_states(model, 10, true, 12)) == 10 + @test length(StateSpaceLearning.simulate_states(model, 8, false, 12)) == 8 + @test length(StateSpaceLearning.simulate_states(model, 10, false, 0)) == 10 model = StateSpaceLearning.StructuralModel( - rand(3); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=true, - freq_seasonal=2, - outlier=true, - cycle_period=3, - stochastic_cycle=true, - ζ_ω_threshold=0, - Exogenous_X=zeros(10, 0), + rand(100); seasonal="none", cycle="stochastic", cycle_period=3, outlier=false ) - steps_ahead = 2 + StateSpaceLearning.fit!(model) + @test length(StateSpaceLearning.simulate_states(model, 10, true, 12)) == 10 + @test length(StateSpaceLearning.simulate_states(model, 8, false, 12)) == 8 + @test length(StateSpaceLearning.simulate_states(model, 10, false, 0)) == 10 +end - X1 = StateSpaceLearning.get_innovation_simulation_X(model, innovation1, steps_ahead) - @assert X1 == [ - 0.0 0.0 0.0 0.0 - 1.0 0.0 0.0 0.0 - 1.0 1.0 0.0 0.0 - 1.0 1.0 1.0 0.0 - 1.0 1.0 1.0 1.0 - 1.0 1.0 1.0 1.0 +@testset "Function: forecast_dynamic_exog_coefs" begin + Random.seed!(1234) + model = StateSpaceLearning.StructuralModel( + rand(100); seasonal="none", cycle="stochastic", cycle_period=3, outlier=false + ) + StateSpaceLearning.fit!(model) + @test StateSpaceLearning.forecast_dynamic_exog_coefs( + model, 10, Vector{Vector}(undef, 0) + ) == zeros(10) + @test StateSpaceLearning.forecast_dynamic_exog_coefs( + model, 8, Vector{Vector}(undef, 0) + ) == zeros(8) + + dynamic_exog_coefs = [ + (collect(1:100), "level"), + (collect(1:100), "slope"), + (collect(1:100), "seasonal", 2), + (collect(1:100), "cycle", 3), ] - - X2 = StateSpaceLearning.get_innovation_simulation_X(model, innovation2, steps_ahead) - @assert X2 == [ - 0.0 0.0 0.0 - 0.0 0.0 0.0 - 1.0 0.0 0.0 - 2.0 1.0 0.0 - 3.0 2.0 1.0 - 4.0 3.0 2.0 + forecast_dynamic_exog_coefs = [ + collect(101:110), collect(101:110), collect(101:110), collect(101:110) ] + model2 = StateSpaceLearning.StructuralModel( + rand(100); dynamic_exog_coefs=dynamic_exog_coefs + ) + StateSpaceLearning.fit!(model2) + @test StateSpaceLearning.forecast_dynamic_exog_coefs( + model2, 10, forecast_dynamic_exog_coefs + ) != zeros(10) + @test length( + StateSpaceLearning.forecast_dynamic_exog_coefs( + model2, 10, forecast_dynamic_exog_coefs + ), + ) == 10 +end - X3 = StateSpaceLearning.get_innovation_simulation_X(model, innovation3, steps_ahead) - @assert X3 == [ - 0.0 0.0 0.0 0.0 - 0.0 0.0 0.0 0.0 - -1.0 1.0 0.0 0.0 - 0.0 -1.0 1.0 0.0 - -1.0 1.0 -1.0 1.0 - 0.0 -1.0 1.0 -1.0 +@testset "Function: forecast" begin + y1 = rand(100) + y2 = rand(100) + y2[10:20] .= NaN + + model1 = StateSpaceLearning.StructuralModel(y1) + StateSpaceLearning.fit!(model1) + @test length(StateSpaceLearning.forecast(model1, 10)) == 10 + + model2 = StateSpaceLearning.StructuralModel(y2; exog=rand(100, 3)) + StateSpaceLearning.fit!(model2) + @test length(StateSpaceLearning.forecast(model2, 10; Exogenous_Forecast=rand(10, 3))) == + 10 + + @test_throws AssertionError StateSpaceLearning.forecast( + model1, 10; Exogenous_Forecast=rand(5, 3) + ) + @test_throws AssertionError StateSpaceLearning.forecast(model2, 10) + @test_throws AssertionError StateSpaceLearning.forecast( + model2, 10; Exogenous_Forecast=rand(5, 3) + ) + + y3 = [ + 4.718, + 4.77, + 4.882, + 4.859, + 4.795, + 4.905, + 4.997, + 4.997, + 4.912, + 4.779, + 4.644, + 4.77, + 4.744, + 4.836, + 4.948, + 4.905, + 4.828, + 5.003, + 5.135, + 5.135, + 5.062, + 4.89, + 4.736, + 4.941, + 4.976, + 5.01, + 5.181, + 5.093, + 5.147, + 5.181, + 5.293, + 5.293, + 5.214, + 5.087, + 4.983, + 5.111, + 5.141, + 5.192, + 5.262, + 5.198, + 5.209, + 5.384, + 5.438, + 5.488, + 5.342, + 5.252, + 5.147, + 5.267, + 5.278, + 5.278, + 5.463, + 5.459, + 5.433, + 5.493, + 5.575, + 5.605, + 5.468, + 5.351, + 5.192, + 5.303, + 5.318, + 5.236, + 5.459, + 5.424, + 5.455, + 5.575, + 5.71, + 5.68, + 5.556, + 5.433, + 5.313, + 5.433, + 5.488, + 5.451, + 5.587, + 5.594, + 5.598, + 5.752, + 5.897, + 5.849, + 5.743, + 5.613, + 5.468, + 5.627, + 5.648, + 5.624, + 5.758, + 5.746, + 5.762, + 5.924, + 6.023, + 6.003, + 5.872, + 5.723, + 5.602, + 5.723, + 5.752, + 5.707, + 5.874, + 5.852, + 5.872, + 6.045, + 6.142, + 6.146, + 6.001, + 5.849, + 5.72, + 5.817, + 5.828, + 5.762, + 5.891, + 5.852, + 5.894, + 6.075, + 6.196, + 6.224, + 6.001, + 5.883, + 5.736, + 5.82, + 5.886, + 5.834, + 6.006, + 5.981, + 6.04, + 6.156, + 6.306, + 6.326, + 6.137, + 6.008, + 5.891, + 6.003, + 6.033, + 5.968, + 6.037, + 6.133, + 6.156, + 6.282, + 6.432, + 6.406, + 6.23, + 6.133, + 5.966, + 6.068, + ] + model3 = StateSpaceLearning.StructuralModel(y3) + StateSpaceLearning.fit!(model3) + forecast3 = trunc.(StateSpaceLearning.forecast(model3, 18); digits=3) + @assert forecast3 == [ + 6.11, + 6.082, + 6.221, + 6.19, + 6.197, + 6.328, + 6.447, + 6.44, + 6.285, + 6.163, + 6.026, + 6.142, + 6.166, + 6.138, + 6.278, + 6.246, + 6.253, + 6.384, ] - X4 = StateSpaceLearning.get_innovation_simulation_X(model, innovation4, steps_ahead) - @assert all( - isapprox.( - X4, - [ - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - -0.5 0.866025 1.0 0.0 0.0 0.0 0.0 0.0 - -0.5 -0.866025 -0.5 0.866025 1.0 0.0 0.0 0.0 - 1.0 -6.10623e-16 -0.5 -0.866025 -0.5 0.866025 1.0 0.0 - -0.5 0.866025 1.0 -6.10623e-16 -0.5 -0.866025 -0.5 0.866025 - ], - atol=1e-6, + model4 = StateSpaceLearning.StructuralModel(y3; freq_seasonal=[12, 36]) + StateSpaceLearning.fit!(model4) + forecast4 = trunc.(StateSpaceLearning.forecast(model4, 18); digits=3) + + @test length(forecast4) == 18 + + exog = rand(length(y3), 3) + model5 = StateSpaceLearning.StructuralModel(y3; exog=exog) + StateSpaceLearning.fit!(model5) + exog_forecast = rand(18, 3) + forecast5 = + trunc.( + StateSpaceLearning.forecast(model5, 18; Exogenous_Forecast=exog_forecast); + digits=3, + ) + @test length(forecast5) == 18 + + dynamic_exog_coefs = [ + (collect(1:length(y3)), "level"), + (collect(1:length(y3)), "slope"), + (collect(1:length(y3)), "seasonal", 2), + (collect(1:length(y3)), "cycle", 3), + ] + forecast_dynamic_exog_coefs = [ + collect((length(y3) + 1):(length(y3) + 10)), + collect((length(y3) + 1):(length(y3) + 10)), + collect((length(y3) + 1):(length(y3) + 10)), + collect((length(y3) + 1):(length(y3) + 10)), + ] + model6 = StateSpaceLearning.StructuralModel(y3; dynamic_exog_coefs=dynamic_exog_coefs) + StateSpaceLearning.fit!(model6) + @test StateSpaceLearning.forecast_dynamic_exog_coefs( + model6, 10, forecast_dynamic_exog_coefs + ) != zeros(10) + @test length( + StateSpaceLearning.forecast_dynamic_exog_coefs( + model6, 10, forecast_dynamic_exog_coefs ), - ) + ) == 10 +end - model2 = StateSpaceLearning.StructuralModel( - rand(4); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=true, - freq_seasonal=2, - outlier=true, - cycle_period=3, - stochastic_cycle=true, - ζ_ω_threshold=0, - Exogenous_X=zeros(10, 0), - stochastic_start=3, - ) - X5 = StateSpaceLearning.get_innovation_simulation_X(model2, innovation4, steps_ahead) - @assert all( - isapprox.( - X5, - [ - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 - -0.5 0.866025 1.0 0.0 0.0 0.0 0.0 0.0 - -0.5 -0.866025 -0.5 0.866025 1.0 0.0 0.0 0.0 - 1.0 -6.10623e-16 -0.5 -0.866025 -0.5 0.866025 1.0 0.0 - -0.5 0.866025 1.0 -6.10623e-16 -0.5 -0.866025 -0.5 0.866025 - ], - atol=1e-6, +@testset "Function: simulate" begin + y1 = rand(100) + y2 = rand(100) + y2[10:20] .= NaN + + model1 = StateSpaceLearning.StructuralModel(y1) + 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; exog=rand(100, 3)) + StateSpaceLearning.fit!(model2) + @test size( + StateSpaceLearning.simulate(model2, 10, 100; Exogenous_Forecast=rand(10, 3)) + ) == (10, 100) + + exog = rand(length(y2), 3) + model5 = StateSpaceLearning.StructuralModel(y2; exog=exog) + StateSpaceLearning.fit!(model5) + exog_forecast = rand(10, 3) + @test size( + StateSpaceLearning.simulate(model5, 10, 100; Exogenous_Forecast=exog_forecast) + ) == (10, 100) + + dynamic_exog_coefs = [ + (collect(1:length(y2)), "level"), + (collect(1:length(y2)), "slope"), + (collect(1:length(y2)), "seasonal", 2), + (collect(1:length(y2)), "cycle", 3), + ] + forecast_dynamic_exog_coefs = [ + collect((length(y2) + 1):(length(y2) + 10)), + collect((length(y2) + 1):(length(y2) + 10)), + collect((length(y2) + 1):(length(y2) + 10)), + collect((length(y2) + 1):(length(y2) + 10)), + ] + model6 = StateSpaceLearning.StructuralModel(y2; dynamic_exog_coefs=dynamic_exog_coefs) + StateSpaceLearning.fit!(model6) + @test size( + StateSpaceLearning.simulate( + model6, 10, 100; dynamic_exog_coefs_forecasts=forecast_dynamic_exog_coefs ), - ) + ) == (10, 100) +end + +@testset "Basics" begin + model = StateSpaceLearning.StructuralModel(rand(100)) + @test StateSpaceLearning.isfitted(model) == false + StateSpaceLearning.fit!(model) + + @test StateSpaceLearning.isfitted(model) == true end diff --git a/test/utils.jl b/test/utils.jl index 42ccfc8..6041a32 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,58 +1,42 @@ @testset "Function: build_components" begin - Exogenous_X1 = rand(10, 3) - Exogenous_X2 = zeros(10, 0) + exog1 = rand(10, 3) + exog2 = zeros(10, 0) Basic_Structural = StateSpaceLearning.StructuralModel( rand(10); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=true, - stochastic_seasonal=true, + level="stochastic", + slope="stochastic", + seasonal="stochastic", freq_seasonal=2, outlier=true, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X1, + ζ_threshold=0, + ω_threshold=0, + exog=exog1, ) Local_Level = StateSpaceLearning.StructuralModel( rand(10); - level=true, - stochastic_level=true, - trend=false, - stochastic_trend=false, - seasonal=false, - stochastic_seasonal=false, - freq_seasonal=2, + level="stochastic", + slope="none", + seasonal="none", outlier=true, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X1, + exog=exog1, ) Local_Linear_Trend1 = StateSpaceLearning.StructuralModel( rand(10); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=false, - stochastic_seasonal=false, - freq_seasonal=2, + level="stochastic", + slope="stochastic", + seasonal="none", outlier=false, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X1, + exog=exog1, ) Local_Linear_Trend2 = StateSpaceLearning.StructuralModel( rand(10); - level=true, - stochastic_level=true, - trend=true, - stochastic_trend=true, - seasonal=false, - stochastic_seasonal=false, - freq_seasonal=2, + level="stochastic", + slope="stochastic", + seasonal="none", outlier=false, - ζ_ω_threshold=0, - Exogenous_X=Exogenous_X2, + exog=exog2, + ζ_threshold=0, ) models = [Basic_Structural, Local_Level, Local_Linear_Trend1, Local_Linear_Trend2] @@ -69,7 +53,7 @@ @test "Values" in keys(components[key]) @test "Coefs" in keys(components[key]) @test "Indexes" in keys(components[key]) - @test key == "Exogenous_X" ? "Selected" in keys(components[key]) : true + @test key == "exog" ? "Selected" in keys(components[key]) : true end end end @@ -105,44 +89,72 @@ end @test StateSpaceLearning.has_intercept(X) end -@testset "Function: fill_innovation_coefs" begin - model = StateSpaceLearning.StructuralModel(rand(100)) - StateSpaceLearning.fit!(model) - components = ["ξ", "ζ", "ω_12"] - - valid_indexes = model.output.valid_indexes +@testset "Function: handle_missing_values" begin + y = rand(10) + X = rand(10, 3) + y[1] = NaN + X[1, :] .= NaN + y_treated, X_treated, valid_indexes = StateSpaceLearning.handle_missing_values(X, y) + @test y_treated == y[2:end] + @test X_treated == X[2:end, :] + @test valid_indexes == 2:10 +end - inov_comp1 = StateSpaceLearning.fill_innovation_coefs( - model, components[1], valid_indexes +@testset "Function: get_stochastic_values" begin + Random.seed!(1234) + + estimated_stochastic = rand(10) + steps_ahead = 5 + T = 10 + start_idx = 1 + final_idx = 10 + seasonal_innovation_simulation1 = 0 + seasonal_innovation_simulation2 = 2 + + st_values1 = StateSpaceLearning.get_stochastic_values( + estimated_stochastic, + steps_ahead, + T, + start_idx, + final_idx, + seasonal_innovation_simulation1, ) - inov_comp2 = StateSpaceLearning.fill_innovation_coefs( - model, components[2], valid_indexes + st_values2 = StateSpaceLearning.get_stochastic_values( + estimated_stochastic, + steps_ahead, + T, + start_idx, + final_idx, + seasonal_innovation_simulation2, ) - 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 = ["ξ", "ζ", "ω_12"] - 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 + @test length(st_values1) == steps_ahead + @test length(st_values2) == steps_ahead + + @test all( + isapprox.( + st_values1, + [ + 0.6395615996802734 + -0.8396219340580711 + 0.6395615996802734 + -0.5798621201341324 + 0.967142768915383 + ], + atol=1e-6, + ), ) - inov_comp3 = StateSpaceLearning.fill_innovation_coefs( - model, components[3], valid_indexes + @test all( + isapprox.( + st_values2, + [ + 0.520354993723718 + -0.014908849285099945 + -0.13102565622085904 + -0.6395615996802734 + -0.520354993723718 + ], + atol=1e-6, + ), ) - - @test size(inov_comp1) == (100, 3) - @test size(inov_comp2) == (100, 3) - @test size(inov_comp3) == (100, 3) end