From 7a219e3d7ca87b5aee83aaa8f24978a994194c71 Mon Sep 17 00:00:00 2001 From: andre_ramos Date: Wed, 9 Oct 2024 09:43:53 -0400 Subject: [PATCH 1/2] add explanatory variables to unobserved_components --- src/StateSpaceModels.jl | 2 + src/models/unobserved_components.jl | 250 ++++++++++++++++++++------- test/models/unobserved_components.jl | 52 ++++++ 3 files changed, 240 insertions(+), 64 deletions(-) diff --git a/src/StateSpaceModels.jl b/src/StateSpaceModels.jl index 3ef82a1..8d09196 100644 --- a/src/StateSpaceModels.jl +++ b/src/StateSpaceModels.jl @@ -50,6 +50,7 @@ include("models/basicstructural_multivariate.jl") include("models/sarima.jl") include("models/linear_regression.jl") include("models/unobserved_components.jl") + include("models/exponential_smoothing.jl") include("models/naive_models.jl") include("models/dar.jl") @@ -87,6 +88,7 @@ export SparseUnivariateKalmanFilter export StateSpaceModel export UnivariateKalmanFilter export UnobservedComponents +export UnobservedComponentsExplanatory export VehicleTracking # Exported functions diff --git a/src/models/unobserved_components.jl b/src/models/unobserved_components.jl index 172fe31..b9462ac 100644 --- a/src/models/unobserved_components.jl +++ b/src/models/unobserved_components.jl @@ -109,6 +109,7 @@ end @doc raw""" UnobservedComponents( y::Vector{Fl}; + X::Matrix{Fl} = zeros(Fl, length(y), 0), trend::String = "local level", seasonal::String = "no" cycle::String = "no" @@ -247,34 +248,59 @@ mutable struct UnobservedComponents <: StateSpaceModel has_cycle::Bool stochastic_cycle::Bool damped_cycle::Bool - system::LinearUnivariateTimeInvariant + system::Union{LinearUnivariateTimeVariant, LinearUnivariateTimeInvariant} results::Results + exogenous::Matrix function UnobservedComponents(y::Vector{Fl}; + X::Matrix{Fl} = zeros(Fl, length(y), 0), trend::String = "local level", seasonal::String = "no", # for example "deterministic 3" or "stochastic 12" cycle::String = "no" # for example "deterministic damped" or "stochastic" ) where Fl + @assert length(y) == size(X, 1) + num_observations = size(X, 1) + num_exogenous = size(X, 2) + + if cycle != "no" && num_exogenous > 0 + @warn "UnobservedComponents is currently unstable when cycle and exogenous variables are included." + end + (has_irregular, has_trend, stochastic_trend, has_slope, stochastic_slope) = parse_trend(trend) (has_seasonal, stochastic_seasonal, seasonal_freq) = parse_seasonal(seasonal) (has_cycle, stochastic_cycle, damped_cycle) = parse_cycle(cycle) # Define system matrices - Z = build_Z(Fl, has_trend, has_slope, has_seasonal, seasonal_freq, has_cycle) - T = build_T(Fl, has_trend, has_slope, has_seasonal, seasonal_freq, has_cycle) - R = build_R(Fl, has_trend, has_slope, has_seasonal, seasonal_freq, has_cycle, - stochastic_trend, stochastic_slope, stochastic_seasonal, stochastic_cycle) - d = zero(Fl) - c = build_c(T) - H = zero(Fl) - Q = build_Q(R) - - system = LinearUnivariateTimeInvariant{Fl}(y, Z, T, R, d, c, H, Q) + if num_exogenous > 0 + Z = [vcat(build_Z(Fl, has_trend, has_slope, has_seasonal, seasonal_freq, has_cycle), X[t, :]) for t in 1:num_observations] + T = build_T(Fl, has_trend, has_slope, has_seasonal, seasonal_freq, has_cycle; num_observations = num_observations, num_exogenous = num_exogenous) + R = build_R(Fl, has_trend, has_slope, has_seasonal, seasonal_freq, has_cycle, + stochastic_trend, stochastic_slope, stochastic_seasonal, stochastic_cycle; num_observations = num_observations, num_exogenous = num_exogenous) + d = [zero(Fl) for _ in 1:num_observations] + c = [build_c(T[1]) for _ in 1:num_observations] + H = [zero(Fl) for _ in 1:num_observations] + Q = [build_Q(R[1]) for _ in 1:num_observations] + + system = LinearUnivariateTimeVariant{Fl}(y, Z, T, R, d, c, H, Q) + else + Z = build_Z(Fl, has_trend, has_slope, has_seasonal, seasonal_freq, has_cycle) + T = build_T(Fl, has_trend, has_slope, has_seasonal, seasonal_freq, has_cycle) + R = build_R(Fl, has_trend, has_slope, has_seasonal, seasonal_freq, has_cycle, + stochastic_trend, stochastic_slope, stochastic_seasonal, stochastic_cycle) + d = zero(Fl) + c = build_c(T) + H = zero(Fl) + Q = build_Q(R) + + system = LinearUnivariateTimeInvariant{Fl}(y, Z, T, R, d, c, H, Q) + + end # Define hyperparameters names names = build_names(has_irregular, stochastic_trend, stochastic_slope, - stochastic_seasonal, has_cycle, stochastic_cycle, damped_cycle) + stochastic_seasonal, has_cycle, stochastic_cycle, damped_cycle; num_exogenous = num_exogenous) + hyperparameters = HyperParameters{Fl}(names) return new(hyperparameters, trend, seasonal, cycle, @@ -282,12 +308,14 @@ mutable struct UnobservedComponents <: StateSpaceModel has_slope, stochastic_slope, has_seasonal, stochastic_seasonal, seasonal_freq, has_cycle, stochastic_cycle, damped_cycle, - system, Results{Fl}()) + system, Results{Fl}(), X) end end +has_exogenous(model::UnobservedComponents) = size(model.exogenous, 2) > 0 + num_components(model::UnobservedComponents) = (model.has_trend + model.has_slope + model.has_seasonal + - model.has_cycle) + model.has_cycle + size(model.exogenous, 2)) function dict_components(model::UnobservedComponents) dict_components = OrderedDict{String, Int}() i = 1 @@ -307,6 +335,11 @@ function dict_components(model::UnobservedComponents) dict_components["Cycle"] = i i += model.has_cycle + 1 end + for j in 1:size(model.exogenous, 2) + dict_components["Exogenous $j"] = i + i += 1 + end + return dict_components end @@ -336,12 +369,14 @@ function build_T(Fl::DataType, has_slope::Bool, has_seasonal::Bool, seasonal_freq::Int, - has_cycle::Bool) + has_cycle::Bool; + num_observations::Int = 0, + num_exogenous::Int = 0) # Caalculate how maany states for each component num_trend_states = has_trend + has_slope num_seasonal_states = has_seasonal ? seasonal_freq - 1 : 0 num_cycle_states = 2 * has_cycle - num_states = num_trend_states + num_seasonal_states + num_cycle_states + num_states = num_trend_states + num_seasonal_states + num_cycle_states + num_exogenous T = zeros(Fl, num_states, num_states) if has_trend T[1, 1] = one(Fl) @@ -361,7 +396,13 @@ function build_T(Fl::DataType, if has_cycle # Do nothing T should be filled with the cycle expression end - return T + if num_exogenous > 0 + T[end-num_exogenous+1:end, end-num_exogenous+1:end] = Matrix{Fl}(I, num_exogenous, num_exogenous) + return [T for _ in 1:num_observations] + else + return T + end + end function build_R(Fl::DataType, has_trend::Bool, @@ -372,20 +413,26 @@ function build_R(Fl::DataType, stochastic_trend::Bool, stochastic_slope::Bool, stochastic_seasonal::Bool, - stochastic_cycle::Bool) + stochastic_cycle::Bool; + num_observations::Int = 0, + num_exogenous::Int = 0) # Assign some model properties num_trend_states = has_trend + has_slope num_stochastic_trend = stochastic_trend + stochastic_slope num_seasonal_states = has_seasonal ? seasonal_freq - 1 : 0 num_cycle_states = 2 * has_cycle - num_states = num_trend_states + num_seasonal_states + num_cycle_states + num_states = num_trend_states + num_seasonal_states + num_cycle_states + num_exogenous num_stochastic_states = num_stochastic_trend + stochastic_seasonal + 2 * stochastic_cycle R = zeros(num_states, num_stochastic_states) for i in 1:num_stochastic_states R[i, i] = one(Fl) end - return R + if num_exogenous > 0 + return [R for _ in 1:num_observations] + else + return R + end end function build_c(T::Matrix{Fl}) where Fl m = size(T, 1) @@ -401,7 +448,8 @@ function build_names(has_irregular::Bool, stochastic_seasonal::Bool, has_cycle::Bool, stochastic_cycle::Bool, - damped_cycle::Bool) + damped_cycle::Bool; + num_exogenous::Int = 0) names = String[] if has_irregular push!(names, "sigma2_irregular") @@ -424,7 +472,7 @@ function build_names(has_irregular::Bool, push!(names, "ρ_cycle") end end - return names + return vcat(names, ["β_$i" for i in 1:num_exogenous]) end has_sigma2(str::String) = occursin("sigma2", str) @@ -470,11 +518,23 @@ function initial_hyperparameters!(model::UnobservedComponents) initial_hyperparameters["ρ_cycle"] = Fl(0.7) end end + 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::UnobservedComponents, i::Int) + num_non_exog_components = num_components(model) - size(model.exogenous, 2) + return model.hyperparameters.names[i + num_non_exog_components + 1] +end + function constrain_hyperparameters!(model::UnobservedComponents) + for i in axes(model.exogenous, 2) + constrain_identity!(model, get_beta_name(model, i)) + end Fl = typeof_model_elements(model) for variance_name in sigma2_names(model) constrain_variance!(model, variance_name) @@ -490,6 +550,9 @@ function constrain_hyperparameters!(model::UnobservedComponents) end function unconstrain_hyperparameters!(model::UnobservedComponents) + for i in axes(model.exogenous, 2) + unconstrain_identity!(model, get_beta_name(model, i)) + end Fl = typeof_model_elements(model) for variance_name in sigma2_names(model) unconstrain_variance!(model, variance_name) @@ -507,59 +570,118 @@ end function fill_model_system!(model::UnobservedComponents) # TODO (performance) maybe cache this indices # This creates some extra allocations - diag_Q_idx = diag_indices(size(model.system.Q, 1)) - idx = 1 + if has_exogenous(model) + diag_Q_idx = diag_indices(size(model.system.Q[1], 1)) + + for t in 1:length(model.system.Q) + idx = 1 + if model.has_irregular + H = get_constrained_value(model, "sigma2_irregular") + fill_H_in_time(model, H) + end + if model.stochastic_trend + model.system.Q[t][diag_Q_idx[idx]] = get_constrained_value(model, "sigma2_trend") + idx += 1 + end + if model.stochastic_slope + model.system.Q[t][diag_Q_idx[idx]] = get_constrained_value(model, "sigma2_slope") + idx += 1 + end + if model.stochastic_seasonal + model.system.Q[t][diag_Q_idx[idx]] = get_constrained_value(model, "sigma2_seasonal") + idx += 1 + end + if model.stochastic_cycle + model.system.Q[t][diag_Q_idx[idx]] = get_constrained_value(model, "sigma2_cycle") + idx += 1 + model.system.Q[t][diag_Q_idx[idx]] = get_constrained_value(model, "sigma2_cycle") + idx += 1 + end + + if model.has_cycle + num_trend_states = model.has_trend + model.has_slope + num_seasonal_states = model.has_seasonal ? model.seasonal_freq - 1 : 0 + num_trend_seasonal_states = num_trend_states + num_seasonal_states + cycle_rows = num_trend_seasonal_states+1:num_trend_seasonal_states+2 + cycle_cols = cycle_rows + if model.damped_cycle + ρ = get_constrained_value(model, "ρ_cycle") + model.system.T[t][cycle_rows[1], cycle_rows[1]] = ρ * cos(get_constrained_value(model, "λ_cycle")) + model.system.T[t][cycle_rows[2], cycle_rows[1]] = ρ * -sin(get_constrained_value(model, "λ_cycle")) + model.system.T[t][cycle_rows[1], cycle_rows[2]] = ρ * sin(get_constrained_value(model, "λ_cycle")) + model.system.T[t][cycle_rows[2], cycle_rows[2]] = ρ * cos(get_constrained_value(model, "λ_cycle")) + else + model.system.T[t][cycle_rows[1], cycle_rows[1]] = cos(get_constrained_value(model, "λ_cycle")) + model.system.T[t][cycle_rows[2], cycle_rows[1]] = -sin(get_constrained_value(model, "λ_cycle")) + model.system.T[t][cycle_rows[1], cycle_rows[2]] = sin(get_constrained_value(model, "λ_cycle")) + model.system.T[t][cycle_rows[2], cycle_rows[2]] = cos(get_constrained_value(model, "λ_cycle")) + end + end + end + else + diag_Q_idx = diag_indices(size(model.system.Q, 1)) + idx = 1 - if model.has_irregular - model.system.H = get_constrained_value(model, "sigma2_irregular") - end - if model.stochastic_trend - model.system.Q[diag_Q_idx[idx]] = get_constrained_value(model, "sigma2_trend") - idx += 1 - end - if model.stochastic_slope - model.system.Q[diag_Q_idx[idx]] = get_constrained_value(model, "sigma2_slope") - idx += 1 - end - if model.stochastic_seasonal - model.system.Q[diag_Q_idx[idx]] = get_constrained_value(model, "sigma2_seasonal") - idx += 1 - end - if model.stochastic_cycle - model.system.Q[diag_Q_idx[idx]] = get_constrained_value(model, "sigma2_cycle") - idx += 1 - model.system.Q[diag_Q_idx[idx]] = get_constrained_value(model, "sigma2_cycle") - idx += 1 - end + if model.has_irregular + model.system.H = get_constrained_value(model, "sigma2_irregular") + end + if model.stochastic_trend + model.system.Q[diag_Q_idx[idx]] = get_constrained_value(model, "sigma2_trend") + idx += 1 + end + if model.stochastic_slope + model.system.Q[diag_Q_idx[idx]] = get_constrained_value(model, "sigma2_slope") + idx += 1 + end + if model.stochastic_seasonal + model.system.Q[diag_Q_idx[idx]] = get_constrained_value(model, "sigma2_seasonal") + idx += 1 + end + if model.stochastic_cycle + model.system.Q[diag_Q_idx[idx]] = get_constrained_value(model, "sigma2_cycle") + idx += 1 + model.system.Q[diag_Q_idx[idx]] = get_constrained_value(model, "sigma2_cycle") + idx += 1 + end - if model.has_cycle - num_trend_states = model.has_trend + model.has_slope - num_seasonal_states = model.has_seasonal ? model.seasonal_freq - 1 : 0 - num_trend_seasonal_states = num_trend_states + num_seasonal_states - cycle_rows = num_trend_seasonal_states+1:num_trend_seasonal_states+2 - cycle_cols = cycle_rows - if model.damped_cycle - ρ = get_constrained_value(model, "ρ_cycle") - model.system.T[cycle_rows[1], cycle_rows[1]] = ρ * cos(get_constrained_value(model, "λ_cycle")) - model.system.T[cycle_rows[2], cycle_rows[1]] = ρ * -sin(get_constrained_value(model, "λ_cycle")) - model.system.T[cycle_rows[1], cycle_rows[2]] = ρ * sin(get_constrained_value(model, "λ_cycle")) - model.system.T[cycle_rows[2], cycle_rows[2]] = ρ * cos(get_constrained_value(model, "λ_cycle")) - else - model.system.T[cycle_rows[1], cycle_rows[1]] = cos(get_constrained_value(model, "λ_cycle")) - model.system.T[cycle_rows[2], cycle_rows[1]] = -sin(get_constrained_value(model, "λ_cycle")) - model.system.T[cycle_rows[1], cycle_rows[2]] = sin(get_constrained_value(model, "λ_cycle")) - model.system.T[cycle_rows[2], cycle_rows[2]] = cos(get_constrained_value(model, "λ_cycle")) + if model.has_cycle + num_trend_states = model.has_trend + model.has_slope + num_seasonal_states = model.has_seasonal ? model.seasonal_freq - 1 : 0 + num_trend_seasonal_states = num_trend_states + num_seasonal_states + cycle_rows = num_trend_seasonal_states+1:num_trend_seasonal_states+2 + cycle_cols = cycle_rows + if model.damped_cycle + ρ = get_constrained_value(model, "ρ_cycle") + model.system.T[cycle_rows[1], cycle_rows[1]] = ρ * cos(get_constrained_value(model, "λ_cycle")) + model.system.T[cycle_rows[2], cycle_rows[1]] = ρ * -sin(get_constrained_value(model, "λ_cycle")) + model.system.T[cycle_rows[1], cycle_rows[2]] = ρ * sin(get_constrained_value(model, "λ_cycle")) + model.system.T[cycle_rows[2], cycle_rows[2]] = ρ * cos(get_constrained_value(model, "λ_cycle")) + else + model.system.T[cycle_rows[1], cycle_rows[1]] = cos(get_constrained_value(model, "λ_cycle")) + model.system.T[cycle_rows[2], cycle_rows[1]] = -sin(get_constrained_value(model, "λ_cycle")) + model.system.T[cycle_rows[1], cycle_rows[2]] = sin(get_constrained_value(model, "λ_cycle")) + model.system.T[cycle_rows[2], cycle_rows[2]] = cos(get_constrained_value(model, "λ_cycle")) + end end end return model end +function fill_H_in_time(model::UnobservedComponents, H::Fl) where Fl + return fill_system_matrice_with_value_in_time(model.system.H, H) +end + function reinstantiate(model::UnobservedComponents, y::Vector{Fl}) where Fl - return UnobservedComponents(y; + return UnobservedComponents(y; trend = model.trend, seasonal = model.seasonal, cycle = model.cycle) end -has_exogenous(::UnobservedComponents) = false +function reinstantiate(model::UnobservedComponents, y::Vector{Fl}, X::Matrix{Fl}) where Fl + return UnobservedComponents(y; X = X, + trend = model.trend, + seasonal = model.seasonal, + cycle = model.cycle) +end \ No newline at end of file diff --git a/test/models/unobserved_components.jl b/test/models/unobserved_components.jl index 0ce640e..2de8150 100644 --- a/test/models/unobserved_components.jl +++ b/test/models/unobserved_components.jl @@ -45,4 +45,56 @@ rj_temp = CSV.File(StateSpaceModels.RJ_TEMPERATURE) |> DataFrame model = UnobservedComponents(rj_temp.Values; trend = "smooth trend", cycle = "stochastic") fit!(model) +end + +@testset "UnobservedComponentsExplanatory" begin + Random.seed!(123) + @test has_fit_methods(UnobservedComponents) + + nile = CSV.File(StateSpaceModels.NILE) |> DataFrame + steps_ahead = 10 + nile_y = nile.flow + X = [0.906 0.247 0.065; 0.443 0.211 0.778; 0.745 0.32 0.081; 0.512 0.011 0.95; 0.253 0.99 0.136; 0.334 0.285 0.345; 0.427 0.078 0.449; 0.867 0.468 0.119; 0.099 0.476 0.04; 0.125 0.157 0.979; 0.692 0.328 0.558; 0.136 0.279 0.91; 0.032 0.022 0.346; 0.35 0.321 0.748; 0.93 0.852 0.115; 0.959 0.023 0.291; 0.774 0.631 0.059; 0.183 0.62 0.563; 0.297 0.434 0.124; 0.15 0.788 0.941; 0.893 0.851 0.298; 0.354 0.623 0.814; 0.131 0.355 0.719; 0.941 0.847 0.823; 0.057 0.719 0.176; 0.245 0.064 0.287; 0.805 0.366 0.027; 0.337 0.28 0.191; 0.823 0.893 0.4; 0.45 0.054 0.463; 0.891 0.661 0.365; 0.711 0.379 0.718; 0.36 0.792 0.0; 0.259 0.398 0.783; 0.39 0.832 0.179; 0.461 0.828 0.061; 0.934 0.575 0.829; 0.753 0.903 0.832; 0.729 0.28 0.71; 0.162 0.505 0.684; 0.637 0.266 0.309; 0.991 0.824 0.73; 0.383 0.371 0.828; 0.618 0.542 0.516; 0.484 0.503 0.101; 0.599 0.443 0.559; 0.453 0.087 0.414; 0.324 0.545 0.931; 0.51 0.596 0.049; 0.656 0.677 0.855; 0.869 0.206 0.392; 0.373 0.78 0.09; 0.692 0.564 0.374; 0.746 0.839 0.649; 0.096 0.342 0.578; 0.459 0.999 0.078; 0.608 0.776 0.99; 0.836 0.719 0.063; 0.283 0.652 0.737; 0.631 0.752 0.883; 0.691 0.081 0.348; 0.821 0.364 0.711; 0.564 0.881 0.045; 0.03 0.6 0.734; 0.091 0.975 0.586; 0.766 0.965 0.567; 0.246 0.14 0.612; 0.096 0.371 0.774; 0.092 0.797 0.055; 0.437 0.096 0.641; 0.231 0.895 0.384; 0.647 0.634 0.325; 0.681 0.612 0.507; 0.122 0.304 0.524; 0.648 0.364 0.023; 0.587 0.294 0.481; 0.377 0.61 0.275; 0.24 0.267 0.63; 0.222 0.24 0.283; 0.509 0.195 0.741; 0.218 0.944 0.804; 0.834 0.72 0.432; 0.563 0.016 0.229; 0.42 0.571 0.976; 0.62 0.31 0.238; 0.292 0.786 0.647; 0.145 0.938 0.018; 0.864 0.918 0.967; 0.926 0.09 0.712; 0.068 0.688 0.407; 0.047 0.739 0.862; 0.426 0.323 0.057; 0.972 0.092 0.697; 0.928 0.761 0.841; 0.077 0.757 0.522; 0.538 0.98 0.388; 0.357 0.843 0.108; 0.471 0.2 0.491; 0.97 0.026 0.486; 0.612 0.999 0.249; 0.539 0.877 0.886; 0.156 0.892 0.617; 0.159 0.57 0.768; 0.485 0.081 0.468; 0.924 0.506 0.68; 0.219 0.068 0.711; 0.457 0.461 0.514; 0.87 0.21 0.412; 0.907 0.305 0.581; 0.688 0.441 0.311] + β = [ -1.937, 4.389, -4.870] + X_train = X[1:end - steps_ahead, :] + X_test = X[end - steps_ahead + 1:end, :] + nile_y += X_train*β + + model = UnobservedComponents(nile_y; X = X_train) + fit!(model) + @test loglike(model) ≈ -620.58519 atol = 1e-5 rtol = 1e-5 + filt = kalman_filter(model) + smoother = kalman_smoother(model) + forec = forecast(model, X_test) + @test monotone_forecast_variance(forec) + + finland_fatalities = CSV.File(StateSpaceModels.VEHICLE_FATALITIES) |> DataFrame + log_finland_fatalities = log.(finland_fatalities.ff) + X = [0.906 0.484 0.926; 0.443 0.599 0.068; 0.745 0.453 0.047; 0.512 0.324 0.426; 0.253 0.51 0.972; 0.334 0.656 0.928; 0.427 0.869 0.077; 0.867 0.373 0.538; 0.099 0.692 0.357; 0.125 0.746 0.471; 0.692 0.096 0.97; 0.136 0.459 0.612; 0.032 0.608 0.539; 0.35 0.836 0.156; 0.93 0.283 0.159; 0.959 0.631 0.485; 0.774 0.691 0.924; 0.183 0.821 0.219; 0.297 0.564 0.457; 0.15 0.03 0.87; 0.893 0.091 0.907; 0.354 0.766 0.688; 0.131 0.246 0.247; 0.941 0.096 0.211; 0.057 0.092 0.32; 0.245 0.437 0.011; 0.805 0.231 0.99; 0.337 0.647 0.285; 0.823 0.681 0.078; 0.45 0.122 0.468; 0.891 0.648 0.476; 0.711 0.587 0.157; 0.36 0.377 0.328; 0.259 0.24 0.279; 0.39 0.222 0.022; 0.461 0.509 0.321; 0.934 0.218 0.852; 0.753 0.834 0.023; 0.729 0.563 0.631; 0.162 0.42 0.62; 0.637 0.62 0.581; 0.991 0.292 0.311; 0.383 0.145 0.121; 0.618 0.864 0.204] + β = [3.028, 0.086, 0.512] + X_train = X[1:end - steps_ahead, :] + X_test = X[end - steps_ahead + 1:end, :] + log_finland_fatalities += X_train*β + model = UnobservedComponents(log_finland_fatalities;X = X_train, trend = "local linear trend") + fit!(model) + @test loglike(model) ≈ 21.6918 atol = 1e-5 rtol = 1e-5 + forec = forecast(model, X_test) + @test monotone_forecast_variance(forec) + + air_passengers = CSV.File(StateSpaceModels.AIR_PASSENGERS) |> DataFrame + log_air_passengers = log.(air_passengers.passengers) + X = [0.121 0.575 0.614; 0.137 0.619 0.881; 0.05 0.45 0.098; 0.386 0.515 0.536; 0.846 0.388 0.115; 0.835 0.731 0.818; 0.16 0.04 0.365; 0.904 0.794 0.9; 0.991 0.497 0.354; 0.74 0.396 0.492; 0.229 0.855 0.15; 0.969 0.628 0.544; 0.54 0.964 0.391; 0.546 0.514 0.763; 0.783 0.316 0.799; 0.813 0.152 0.763; 0.83 0.845 0.824; 0.656 0.651 0.374; 0.492 0.655 0.382; 0.1 0.528 0.67; 0.878 0.98 0.573; 0.959 0.956 0.487; 0.634 0.279 0.057; 0.064 0.634 0.291; 0.018 0.443 0.494; 0.16 0.445 0.85; 0.605 0.241 0.269; 0.732 0.819 0.088; 0.439 0.201 0.683; 0.265 0.068 0.738; 0.966 0.321 0.001; 0.783 0.563 0.935; 0.66 0.175 0.698; 0.357 0.975 0.065; 0.455 0.295 0.256; 0.519 0.656 0.04; 0.17 0.448 0.445; 0.072 0.972 0.129; 0.428 0.378 0.492; 0.907 0.067 0.947; 0.991 0.802 0.413; 0.048 0.756 0.986; 0.695 0.641 0.311; 0.319 0.008 0.154; 0.084 0.398 0.522; 0.379 0.21 0.627; 0.298 0.343 0.537; 0.968 0.681 0.117; 0.32 0.469 0.884; 0.315 0.767 0.392; 0.869 0.427 0.456; 0.41 0.392 0.237; 0.204 0.874 0.104; 0.25 0.512 0.376; 0.417 0.078 0.434; 0.806 0.733 0.628; 0.307 0.78 0.457; 0.731 0.628 0.034; 0.553 0.467 0.155; 0.6 0.245 0.304; 0.139 0.825 0.547; 0.052 0.618 0.801; 0.269 0.607 0.652; 0.295 0.861 0.409; 0.602 0.169 0.891; 0.907 0.277 0.555; 0.127 0.504 0.294; 0.865 0.724 0.235; 0.228 0.013 0.193; 0.584 0.698 0.457; 0.157 0.868 0.705; 0.302 0.673 0.641; 0.486 0.359 0.784; 0.915 0.477 0.329; 0.741 0.193 0.023; 0.386 0.936 0.01; 0.92 0.03 0.42; 0.037 0.838 0.788; 0.923 0.745 0.225; 0.199 0.28 0.761; 0.702 0.372 0.745; 0.284 0.81 0.861; 0.301 0.52 0.452; 0.712 0.953 0.624; 0.459 0.574 0.769; 0.097 0.257 0.332; 0.594 0.526 0.014; 0.533 0.906 0.049; 0.94 0.149 0.629; 0.496 0.736 0.392; 0.597 0.313 0.123; 0.129 0.683 0.218; 0.407 0.317 0.398; 0.363 0.509 0.369; 0.487 0.695 0.98; 0.047 0.134 0.24; 0.273 0.727 0.069; 0.457 0.259 0.338; 0.825 0.509 0.401; 0.371 0.092 0.267; 0.186 0.198 0.363; 0.88 0.116 0.889; 0.387 0.041 0.812; 0.535 0.034 0.194; 0.923 0.59 0.38; 0.955 0.073 0.299; 0.765 0.169 0.8; 0.39 0.871 0.788; 0.907 0.919 0.067; 0.46 0.252 0.988; 0.615 0.363 0.861; 0.769 0.098 0.741; 0.046 0.096 0.104; 0.391 0.033 0.121; 0.66 0.142 0.455; 0.589 0.937 0.518; 0.13 0.059 0.895; 0.215 0.644 0.998; 0.938 0.51 0.378; 0.124 0.498 0.954; 0.304 0.966 0.948; 0.411 0.998 0.863; 0.214 0.724 0.539; 0.093 0.907 0.409; 0.33 0.325 0.438; 0.944 0.574 0.093; 0.24 0.36 0.569; 0.533 0.42 0.907; 0.252 0.951 0.448; 0.588 0.154 0.344; 0.019 0.961 0.015; 0.355 0.221 0.353; 0.812 0.145 0.933; 0.118 0.127 0.227; 0.181 0.352 0.502; 0.019 0.45 0.533; 0.757 0.76 0.485; 0.614 0.861 0.739; 0.829 0.892 0.79; 0.602 0.16 0.32; 0.814 0.653 0.382; 0.626 0.185 0.525; 0.133 0.571 0.423; 0.077 0.166 0.477; 0.601 0.862 0.135; 0.243 0.772 0.579; 0.335 0.329 0.651; 0.359 0.462 0.588; 0.229 0.04 0.932; 0.152 0.161 0.84; 0.949 0.695 0.795; 0.92 0.528 0.794; 0.895 0.911 0.64; 0.552 0.543 0.445] + β = [-0.408, 2.838, 4.718] + X_train = X[1:end - steps_ahead, :] + X_test = X[end - steps_ahead + 1:end, :] + log_air_passengers += X_train*β + model = UnobservedComponents(log_air_passengers;X = X_train, trend = "local linear trend", seasonal = "stochastic 12") + fit!(model) + @test loglike(model) ≈ 219.978 atol = 1e-5 rtol = 1e-5 + forec = forecast(model, X_test) + @test monotone_forecast_variance(forec) + smoother = kalman_smoother(model) + alpha = get_smoothed_state(smoother) + @test size(alpha) == (144, 16) + end \ No newline at end of file From 5c4ea4b2a13d013e717c445e2cb1d6b1657bdf8f Mon Sep 17 00:00:00 2001 From: andre_ramos Date: Wed, 9 Oct 2024 16:35:57 -0400 Subject: [PATCH 2/2] add simulation functions for unobserved components with explanatory --- src/StateSpaceModels.jl | 1 - src/forecast.jl | 8 +-- src/models/unobserved_components.jl | 82 ++++++++++++++++++++++++++++ test/models/unobserved_components.jl | 13 +++++ 4 files changed, 99 insertions(+), 5 deletions(-) diff --git a/src/StateSpaceModels.jl b/src/StateSpaceModels.jl index 8d09196..2324994 100644 --- a/src/StateSpaceModels.jl +++ b/src/StateSpaceModels.jl @@ -88,7 +88,6 @@ export SparseUnivariateKalmanFilter export StateSpaceModel export UnivariateKalmanFilter export UnobservedComponents -export UnobservedComponentsExplanatory export VehicleTracking # Exported functions diff --git a/src/forecast.jl b/src/forecast.jl index 28296cc..bb1c299 100644 --- a/src/forecast.jl +++ b/src/forecast.jl @@ -89,13 +89,13 @@ function forecast( for i in 1:steps_ahead if isunivariate(model) expected_value[i] = [ - dot(model.system.Z[end - steps_ahead + i - 1], fo.a[end - steps_ahead + i - 1]) + - model.system.d[end - steps_ahead + i - 1], + dot(forecasting_model.system.Z[end - steps_ahead + i], fo.a[end - steps_ahead + i - 1]) + + forecasting_model.system.d[end - steps_ahead + i], ] else expected_value[i] = - model.system.Z[end - steps_ahead + i - 1] * fo.a[end - steps_ahead + i - 1] + - model.system.d[end - steps_ahead + i - 1] + forecasting_model.system.Z[end - steps_ahead + i] * fo.a[end - steps_ahead + i - 1] + + forecasting_model.system.d[end - steps_ahead + i] end covariance[i] = fo.F[end - steps_ahead + i] end diff --git a/src/models/unobserved_components.jl b/src/models/unobserved_components.jl index b9462ac..2137738 100644 --- a/src/models/unobserved_components.jl +++ b/src/models/unobserved_components.jl @@ -684,4 +684,86 @@ function reinstantiate(model::UnobservedComponents, y::Vector{Fl}, X::Matrix{Fl} trend = model.trend, seasonal = model.seasonal, cycle = model.cycle) +end + + +# UnobservedComponents with explanatory requires a custom simulation + +function simulate_scenarios( + model::UnobservedComponents, + 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_scenarios( + model::UnobservedComponents, + steps_ahead::Int, + n_scenarios::Int, + new_exogenous::Array{Fl, 3}; + filter::KalmanFilter=default_filter(model), +) where Fl + @assert steps_ahead == size(new_exogenous, 1) "new_exogenous must have the same dimension of steps_ahead" + @assert n_scenarios == size(new_exogenous, 3) "new_exogenous must have the same number of scenarios of n_scenarios" + # 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[:, :, s]) + end + return scenarios +end + +function simulate( + model::UnobservedComponents, + 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_decomposition(sys.Q[1]) + standard_ε = randn(n) + standard_η = randn(n + 1, size(sys.Q[1], 1)) + + num_exogenous = size(model.exogenous, 2) + @assert num_exogenous == size(new_exogenous, 2) "You must have the same number of exogenous variables of the model." + + # The first state of the simulation is the update of a_0 + alpha[1, :] .= initial_state + sys.Z[1][end-num_exogenous+1: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][end-num_exogenous+1: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/unobserved_components.jl b/test/models/unobserved_components.jl index 2de8150..43ae819 100644 --- a/test/models/unobserved_components.jl +++ b/test/models/unobserved_components.jl @@ -97,4 +97,17 @@ end alpha = get_smoothed_state(smoother) @test size(alpha) == (144, 16) + @test_throws AssertionError simulate_scenarios(model, 10, 1000, ones(5, 2)) + model_sim = deepcopy(model) + scenarios = simulate_scenarios(model_sim, 10, 100000, X_test) + test_scenarios_adequacy_with_forecast(forec, scenarios) + #build X_test as 3D array, where all the scenarios are the same + X_test_3d = ones(10, 3, 500) + for i in 1:500 + X_test_3d[:, :, i] = X_test + end + model_sim2 = deepcopy(model) + scenarios = simulate_scenarios(model_sim2, 10, 500, X_test_3d) + test_scenarios_adequacy_with_forecast(forec, scenarios) + end \ No newline at end of file