diff --git a/Project.toml b/Project.toml index e850918..b4ea079 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StateSpaceModels" uuid = "99342f36-827c-5390-97c9-d7f9ee765c78" authors = ["raphaelsaavedra , guilhermebodin , mariohsouto"] -version = "0.5.12" +version = "0.5.13" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" diff --git a/src/StateSpaceModels.jl b/src/StateSpaceModels.jl index e8b3c50..b6b5343 100644 --- a/src/StateSpaceModels.jl +++ b/src/StateSpaceModels.jl @@ -41,6 +41,7 @@ include("models/locallevelcycle.jl") include("models/locallevelexplanatory.jl") include("models/locallineartrend.jl") include("models/basicstructural.jl") +include("models/basicstructural_explanatory.jl") include("models/basicstructural_multivariate.jl") include("models/sarima.jl") include("models/linear_regression.jl") @@ -55,6 +56,7 @@ include("visualization/diagnostics.jl") # Exported types and structs export BasicStructural +export BasicStructuralExplanatory export DAR export ExponentialSmoothing export LinearMultivariateTimeInvariant diff --git a/src/models/basicstructural_explanatory.jl b/src/models/basicstructural_explanatory.jl new file mode 100644 index 0000000..7777e0c --- /dev/null +++ b/src/models/basicstructural_explanatory.jl @@ -0,0 +1,203 @@ +@doc raw""" + BasicStructuralExplanatory(y::Vector{Fl}, s::Int, X::Matrix{Fl}) where Fl + +It is defined by: +```math +\begin{gather*} + \begin{aligned} + y_{t} &= \mu_{t} + \gamma_{t} + \beta_{t, i}X_{t, i} \varepsilon_{t} \quad &\varepsilon_{t} \sim \mathcal{N}(0, \sigma^2_{\varepsilon})\\ + \mu_{t+1} &= \mu_{t} + \nu_{t} + \xi_{t} \quad &\xi_{t} \sim \mathcal{N}(0, \sigma^2_{\xi})\\ + \nu_{t+1} &= \nu_{t} + \zeta_{t} \quad &\zeta_{t} \sim \mathcal{N}(0, \sigma^2_{\zeta})\\ + \gamma_{t+1} &= -\sum_{j=1}^{s-1} \gamma_{t+1-j} + \omega_{t} \quad & \omega_{t} \sim \mathcal{N}(0, \sigma^2_{\omega})\\ + \beta_{t+1} &= \beta_{t} + \end{aligned} +\end{gather*} +``` + +# Example +```jldoctest +julia> model = BasicStructuralExplanatory(rand(100), 12, rand(100, 2)) +BasicStructuralExplanatory model +``` + +# References + * Durbin, James, & Siem Jan Koopman. (2012). "Time Series Analysis by State Space Methods: Second Edition." Oxford University Press. +""" +mutable struct BasicStructuralExplanatory <: StateSpaceModel + hyperparameters::HyperParameters + system::LinearUnivariateTimeVariant + seasonality::Int + results::Results + exogenous::Matrix + + function BasicStructuralExplanatory(y::Vector{Fl}, s::Int, X::Matrix{Fl}) where Fl + + @assert length(y) == size(X, 1) + + num_observations = size(X, 1) + num_exogenous = size(X, 2) + + Z = [vcat([1; 0; 1; zeros(Fl, s - 2)], X[t, :]) for t in 1:num_observations] + T = [[ + 1 1 zeros(Fl, 1, s - 1) zeros(Fl, 1, num_exogenous) + 0 1 zeros(Fl, 1, s - 1) zeros(Fl, 1, num_exogenous) + 0 0 -ones(Fl, 1, s - 1) zeros(Fl, 1, num_exogenous) + zeros(Fl, s - 2, 2) Matrix{Fl}(I, s - 2, s - 2) zeros(Fl, s - 2) zeros(Fl, 10, num_exogenous) + zeros(Fl, num_exogenous, 13) Matrix{Fl}(I, num_exogenous, num_exogenous) + ] for _ in 1:num_observations] + R = [[ + Matrix{Fl}(I, 3, 3) + zeros(Fl, s - 2, 3) + zeros(num_exogenous, 3) + ] for _ in 1:num_observations] + d = [zero(Fl) for _ in 1:num_observations] + c = [zeros(Fl, size(T[1], 1)) for _ in 1:num_observations] + H = [one(Fl) for _ in 1:num_observations] + Q = [zeros(Fl, 3, 3) for _ in 1:num_observations] + + system = LinearUnivariateTimeVariant{Fl}(y, Z, T, R, d, c, H, Q) + + names = [["sigma2_ε", "sigma2_ξ", "sigma2_ζ", "sigma2_ω"]; ["β_$i" for i in 1:num_exogenous]] + + hyperparameters = HyperParameters{Fl}(names) + + return new(hyperparameters, system, s, Results{Fl}(), X) + end +end + +function default_filter(model::BasicStructuralExplanatory) + Fl = typeof_model_elements(model) + steadystate_tol = Fl(1e-5) + a1 = zeros(Fl, num_states(model)) + P1 = zeros(Fl, num_states(model), num_states(model)) + P1[1:13, 1:13] = Fl(1e6) .* Matrix{Fl}(I, 13, 13) + return UnivariateKalmanFilter(a1, P1, 13, steadystate_tol) +end + +function initial_hyperparameters!(model::BasicStructuralExplanatory) + Fl = typeof_model_elements(model) + initial_hyperparameters = Dict{String,Fl}( + "sigma2_ε" => one(Fl), + "sigma2_ξ" => one(Fl), + "sigma2_ζ" => one(Fl), + "sigma2_ω" => one(Fl), + ) + initial_exogenous = model.exogenous \ model.system.y + for i in axes(model.exogenous, 2) + initial_hyperparameters[get_beta_name(model, i)] = initial_exogenous[i] + end + set_initial_hyperparameters!(model, initial_hyperparameters) + return model +end + +function get_beta_name(model::BasicStructuralExplanatory, i::Int) + return model.hyperparameters.names[i + 4] +end + +function constrain_hyperparameters!(model::BasicStructuralExplanatory) + for i in axes(model.exogenous, 2) + constrain_identity!(model, get_beta_name(model, i)) + end + constrain_variance!(model, "sigma2_ε") + constrain_variance!(model, "sigma2_ξ") + constrain_variance!(model, "sigma2_ζ") + constrain_variance!(model, "sigma2_ω") + return model +end + +function unconstrain_hyperparameters!(model::BasicStructuralExplanatory) + for i in axes(model.exogenous, 2) + unconstrain_identity!(model, get_beta_name(model, i)) + end + unconstrain_variance!(model, "sigma2_ε") + unconstrain_variance!(model, "sigma2_ξ") + unconstrain_variance!(model, "sigma2_ζ") + unconstrain_variance!(model, "sigma2_ω") + return model +end + +function fill_model_system!(model::BasicStructuralExplanatory) + H = get_constrained_value(model, "sigma2_ε") + fill_H_in_time(model, H) + for t in 1:length(model.system.Q) + model.system.Q[t][1] = get_constrained_value(model, "sigma2_ξ") + model.system.Q[t][5] = get_constrained_value(model, "sigma2_ζ") + model.system.Q[t][end] = get_constrained_value(model, "sigma2_ω") + end + return model +end + +function fill_model_filter!(filter::KalmanFilter, model::BasicStructuralExplanatory) + for i in axes(model.exogenous, 2) + filter.kalman_state.a[i + 13] = get_constrained_value(model, get_beta_name(model, i)) + end + return filter +end + + +function fill_H_in_time(model::BasicStructuralExplanatory, H::Fl) where Fl + return fill_system_matrice_with_value_in_time(model.system.H, H) +end + +function reinstantiate(model::BasicStructuralExplanatory, y::Vector{Fl}, X::Matrix{Fl}) where Fl + return BasicStructuralExplanatory(y, model.seasonality, X) +end + +has_exogenous(::BasicStructuralExplanatory) = true + +# BasicStructuralExplanatory requires a custom simulation + +function simulate_scenarios( + model::BasicStructuralExplanatory, + steps_ahead::Int, + n_scenarios::Int, + new_exogenous::Matrix{Fl}; + filter::KalmanFilter=default_filter(model), +) where Fl + @assert steps_ahead == size(new_exogenous, 1) "new_exogenous must have the same dimension as steps_ahead" + # Query the type of model elements + fo = kalman_filter(model) + last_state = fo.a[end] + num_series = size(model.system.y, 2) + + scenarios = Array{Fl,3}(undef, steps_ahead, num_series, n_scenarios) + for s in 1:n_scenarios + scenarios[:, :, s] = simulate(model, last_state, steps_ahead, new_exogenous) + end + return scenarios +end + +function simulate( + model::BasicStructuralExplanatory, + initial_state::Vector{Fl}, + n::Int, + new_exogenous::Matrix{Fl}; + return_simulated_states::Bool=false, +) where Fl + sys = model.system + m = size(sys.T[1], 1) + y = Vector{Fl}(undef, n) + alpha = Matrix{Fl}(undef, n + 1, m) + # Sampling errors + chol_H = sqrt(sys.H[1]) + chol_Q = cholesky(sys.Q[1]) + standard_ε = randn(n) + standard_η = randn(n + 1, size(sys.Q[1], 1)) + + # The first state of the simulation is the update of a_0 + alpha[1, :] .= initial_state + sys.Z[1][14:end] .= new_exogenous[1, :] + y[1] = dot(sys.Z[1], initial_state) + sys.d[1] + chol_H * standard_ε[1] + alpha[2, :] = sys.T[1] * initial_state + sys.c[1] + sys.R[1] * chol_Q.L * standard_η[1, :] + # Simulate scenarios + for t in 2:n + sys.Z[t][14:end] .= new_exogenous[t, :] + y[t] = dot(sys.Z[t], alpha[t, :]) + sys.d[t] + chol_H * standard_ε[t] + alpha[t + 1, :] = sys.T[t] * alpha[t, :] + sys.c[t] + sys.R[t] * chol_Q.L * standard_η[t, :] + end + + if return_simulated_states + return y, alpha[1:n, :] + end + return y +end \ No newline at end of file diff --git a/test/models/basicstructural_explanatory.jl b/test/models/basicstructural_explanatory.jl new file mode 100644 index 0000000..1daa861 --- /dev/null +++ b/test/models/basicstructural_explanatory.jl @@ -0,0 +1,21 @@ +@testset "Basic Structural With Explanatory Model" begin + @test has_fit_methods(BasicStructuralExplanatory) + y = CSV.File(StateSpaceModels.AIR_PASSENGERS) |> DataFrame + logap = log.(y.passengers) + X = rand(length(logap), 2) + model = BasicStructuralExplanatory(logap, 12, X) + fit!(model) + model.results + # forecasting + # For a fixed forecasting explanatory the variance must not decrease + forec = forecast(model, ones(10, 2)) + @test monotone_forecast_variance(forec) + kf = kalman_filter(model); + ks = kalman_smoother(model); + a = get_predictive_state(kf) + @test a[1, 14] ≈ a[end, 14] atol=1e-3 + @test a[1, 15] ≈ a[end, 15] atol=1e-3 + @test_throws AssertionError simulate_scenarios(model, 10, 1000, ones(5, 2)) + scenarios = simulate_scenarios(model, 10, 1000, ones(10, 2)) + test_scenarios_adequacy_with_forecast(forec, scenarios) +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index b0737ec..41ee4f8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,6 +19,7 @@ include("models/locallineartrend.jl") include("models/locallevelcycle.jl") include("models/locallevelexplanatory.jl") include("models/basicstructural.jl") +include("models/basicstructural_explanatory.jl") include("models/basicstructural_multivariate.jl") include("models/sarima.jl") include("models/unobserved_components.jl")