Merge pull request #266 from LAMPSPUC/basic_structural_ex

Add BasicStructuralExplanatory

Project.toml

@@ -1,7 +1,7 @@
 name = "StateSpaceModels"
 uuid = "99342f36-827c-5390-97c9-d7f9ee765c78"
 authors = ["raphaelsaavedra <[email protected]>, guilhermebodin <[email protected]>, mariohsouto"]
-version = "0.5.12"
+version = "0.5.13"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

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

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

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

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")