From cf580f0ebfec246ad84633b2d193e2d3ce79d2a4 Mon Sep 17 00:00:00 2001 From: Paul Mainwood Date: Sat, 26 Jun 2021 10:01:00 +0100 Subject: [PATCH 1/6] Add locallevelexplanatorytimevarying with one var. --- src/StateSpaceModels.jl | 2 + src/fit.jl | 2 +- .../locallevelexplanatorytimevarying.jl | 123 ++++++++++++++++++ 3 files changed, 126 insertions(+), 1 deletion(-) create mode 100644 src/models/locallevelexplanatorytimevarying.jl diff --git a/src/StateSpaceModels.jl b/src/StateSpaceModels.jl index 646ce2d..fa92fd0 100644 --- a/src/StateSpaceModels.jl +++ b/src/StateSpaceModels.jl @@ -51,6 +51,7 @@ include("models/unobserved_components.jl") include("models/exponential_smoothing.jl") include("models/naive_models.jl") include("models/dar.jl") +include("models/locallevelexplanatorytimevarying.jl") include("visualization/forecast.jl") include("visualization/components.jl") @@ -71,6 +72,7 @@ export LinearUnivariateTimeVariant export LocalLevel export LocalLevelCycle export LocalLevelExplanatory +export LocalLevelExplanatoryTimeVarying export LocalLinearTrend export MultivariateBasicStructural export Naive diff --git a/src/fit.jl b/src/fit.jl index af982b9..050fe38 100644 --- a/src/fit.jl +++ b/src/fit.jl @@ -46,7 +46,7 @@ function fit!( opt_hyperparameters = opt.minimizer update_model_hyperparameters!(model, opt_hyperparameters) # TODO - # I leaned that this is not a good way to compute the covariance matrix of paarameters + # I leaned that this is not a good way to compute the covariance matrix of parameters # we should investigate other methods numerical_hessian = Optim.hessian!(func, opt_hyperparameters) std_err = diag(pinv(numerical_hessian)) diff --git a/src/models/locallevelexplanatorytimevarying.jl b/src/models/locallevelexplanatorytimevarying.jl new file mode 100644 index 0000000..58424a8 --- /dev/null +++ b/src/models/locallevelexplanatorytimevarying.jl @@ -0,0 +1,123 @@ +@doc raw""" + LocalLevelExplanatoryTimeVarying(y::Vector{Fl}, X::Matrix{Fl}) where Fl + +A local level model with time-varying regressors is defined by: +```math +\begin{gather*} + \begin{aligned} + y_{t} &= \mu_{t} + X_{1,t} \cdot \beta_{1,t} + \dots + X_{n,t} \cdot \beta_{n,t} + \varepsilon_{t} \quad &\varepsilon_{t} \sim \mathcal{N}(0, \sigma^2_{\varepsilon})\\ + \mu_{t+1} &= \mu_{t} + \xi_{t} &\xi_{t} \sim \mathcal{N}(0, \sigma^2_{\xi})\\ + \beta_{1,t+1} &= \beta_{1,t} + \tau_{1, t} &\tau_{1, t} \sim \mathcal{N}(0, \sigma^2_{\tau_{1}})\\ + \dots &= \dots\\ + \beta_{n,t+1} &= \beta_{n,t} + \tau_{n, t} &\tau_{n, t} \sim \mathcal{N}(0, \sigma^2_{\tau_{n}})\\\\ + \end{aligned} +\end{gather*} +``` + +# Example +```jldoctest +julia> model = tvr(rand(100), rand(100, 1)) +tvr +``` +""" +mutable struct LocalLevelExplanatoryTimeVarying <: StateSpaceModel + hyperparameters::HyperParameters + system::LinearUnivariateTimeVariant + results::Results + exogenous::Matrix + + function LocalLevelExplanatoryTimeVarying(y::Vector{Fl}, X::Matrix{Fl}) where Fl + + @assert length(y) == size(X, 1) + + num_observations = size(X, 1) + num_exogenous = size(X, 2) + m = num_exogenous + 1 + + # Define system matrices + Z = [vcat(ones(Fl, 1), X[t, :]) for t in 1:num_observations] + T = [Matrix{Fl}(I, m, m) for _ in 1:num_observations] + R = [Matrix{Fl}(I, m, m)[:, 1:2] for _ in 1:num_observations] + d = [zero(Fl) for _ in 1:num_observations] + c = [zeros(m) for _ in 1:num_observations] + H = [one(Fl) for _ in 1:num_observations] + Q = [Matrix{Fl}(I, 2, 2) for _ in 1:num_observations] + + system = LinearUnivariateTimeVariant{Fl}(y, Z, T, R, d, c, H, Q) + + # Define hyperparameters names + names = [["sigma2_ε", "sigma2_η"];["β_$i" for i in 1:num_exogenous]; ["τ_1"]] + hyperparameters = HyperParameters{Fl}(names) + + return new(hyperparameters, system, Results{Fl}(), X) + end +end + +# Obligatory methods +function default_filter(model::LocalLevelExplanatoryTimeVarying) + Fl = typeof_model_elements(model) + a1 = zeros(Fl, num_states(model)) + P1 = Fl(1e6) .* Matrix{Fl}(I, num_states(model), num_states(model)) + steadystate_tol = Fl(1e-5) + return UnivariateKalmanFilter(a1, P1, num_states(model), steadystate_tol) +end + +function get_beta_name(model::LocalLevelExplanatoryTimeVarying, i::Int) + return model.hyperparameters.names[i + 2] +end + +function initial_hyperparameters!(model::LocalLevelExplanatoryTimeVarying) + #Placeholder initial hyperparamters for tau + Fl = typeof_model_elements(model) + observed_variance = var(model.system.y[findall(!isnan, model.system.y)]) + initial_hyperparameters = Dict{String,Fl}( + "sigma2_ε" => observed_variance, "sigma2_η" => observed_variance, "τ_1" => observed_variance + ) + # This is an heuristic for a good approximation + 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 constrain_hyperparameters!(model::LocalLevelExplanatoryTimeVarying) + 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, "τ_1") + return model +end + +function unconstrain_hyperparameters!(model::LocalLevelExplanatoryTimeVarying) + 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, "τ_1") + return model +end + +function fill_model_system!(model::LocalLevelExplanatoryTimeVarying) + 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][end] = get_constrained_value(model, "τ_1") + end + return model +end + +function fill_H_in_time(model::LocalLevelExplanatoryTimeVarying, H::Fl) where Fl + return fill_system_matrice_with_value_in_time(model.system.H, H) +end + +function reinstantiate(::LocalLevelExplanatoryTimeVarying, y::Vector{Fl}, X::Matrix{Fl}) where Fl + return LocalLevelExplanatoryTimeVarying(y, X) +end + +has_exogenous(::LocalLevelExplanatoryTimeVarying) = true From e14625dc2fe39586d99ee15f6ec21b8ebd2daac8 Mon Sep 17 00:00:00 2001 From: Paul Mainwood Date: Tue, 29 Jun 2021 07:08:26 +0100 Subject: [PATCH 2/6] Convert lletv to all variable. --- .../locallevelexplanatorytimevarying.jl | 41 ++++++++++--------- 1 file changed, 22 insertions(+), 19 deletions(-) diff --git a/src/models/locallevelexplanatorytimevarying.jl b/src/models/locallevelexplanatorytimevarying.jl index 58424a8..abbd8ee 100644 --- a/src/models/locallevelexplanatorytimevarying.jl +++ b/src/models/locallevelexplanatorytimevarying.jl @@ -37,16 +37,16 @@ mutable struct LocalLevelExplanatoryTimeVarying <: StateSpaceModel # Define system matrices Z = [vcat(ones(Fl, 1), X[t, :]) for t in 1:num_observations] T = [Matrix{Fl}(I, m, m) for _ in 1:num_observations] - R = [Matrix{Fl}(I, m, m)[:, 1:2] for _ in 1:num_observations] + R = [Matrix{Fl}(I, m, m) for _ in 1:num_observations] d = [zero(Fl) for _ in 1:num_observations] c = [zeros(m) for _ in 1:num_observations] H = [one(Fl) for _ in 1:num_observations] - Q = [Matrix{Fl}(I, 2, 2) for _ in 1:num_observations] + Q = [Matrix{Fl}(I, m, m) for _ in 1:num_observations] system = LinearUnivariateTimeVariant{Fl}(y, Z, T, R, d, c, H, Q) # Define hyperparameters names - names = [["sigma2_ε", "sigma2_η"];["β_$i" for i in 1:num_exogenous]; ["τ_1"]] + names = [["sigma2_ε", "sigma2_η"];["β_$i" for i in 1:num_exogenous]; ["tau2_$i" for i in 1:num_exogenous]] hyperparameters = HyperParameters{Fl}(names) return new(hyperparameters, system, Results{Fl}(), X) @@ -67,11 +67,11 @@ function get_beta_name(model::LocalLevelExplanatoryTimeVarying, i::Int) end function initial_hyperparameters!(model::LocalLevelExplanatoryTimeVarying) - #Placeholder initial hyperparamters for tau + #Fill all with observed variance - betas redone later, taus are pretty bad Fl = typeof_model_elements(model) observed_variance = var(model.system.y[findall(!isnan, model.system.y)]) initial_hyperparameters = Dict{String,Fl}( - "sigma2_ε" => observed_variance, "sigma2_η" => observed_variance, "τ_1" => observed_variance + model.hyperparameters.names .=> observed_variance ) # This is an heuristic for a good approximation initial_exogenous = model.exogenous \ model.system.y @@ -83,31 +83,34 @@ function initial_hyperparameters!(model::LocalLevelExplanatoryTimeVarying) end function constrain_hyperparameters!(model::LocalLevelExplanatoryTimeVarying) - 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, "τ_1") + betas = [get_beta_name(model, i) for i in axes(model.exogenous, 2)] + constrain_identity!.(Ref(model), betas) + + non_betas = setdiff(model.hyperparameters.names, betas) + constrain_variance!.(Ref(model), non_betas) return model end function unconstrain_hyperparameters!(model::LocalLevelExplanatoryTimeVarying) - 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, "τ_1") + betas = [get_beta_name(model, i) for i in axes(model.exogenous, 2)] + unconstrain_identity!.(Ref(model), betas) + + non_betas = setdiff(model.hyperparameters.names, betas) + unconstrain_variance!.(Ref(model), non_betas) return model end function fill_model_system!(model::LocalLevelExplanatoryTimeVarying) H = get_constrained_value(model, "sigma2_ε") fill_H_in_time(model, H) + + numtaus = Int((length(model.hyperparameters.names) - 2) / 2) + taus = model.hyperparameters.names[Int(end - numtaus + 1): end] + + constrained_values_taus = get_constrained_value.(Ref(model), taus) + for t in 1:length(model.system.Q) - model.system.Q[t][1] = get_constrained_value(model, "sigma2_η") - model.system.Q[t][end] = get_constrained_value(model, "τ_1") + model.system.Q[t] = diagm([get_constrained_value(model, "sigma2_η"); constrained_values_taus]) end return model end From 09da9f75ab8d398a8574559946a53f7a7ae85b7d Mon Sep 17 00:00:00 2001 From: Paul Mainwood Date: Tue, 29 Jun 2021 07:47:27 +0100 Subject: [PATCH 3/6] Update docstrings --- src/models/locallevelexplanatorytimevarying.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/models/locallevelexplanatorytimevarying.jl b/src/models/locallevelexplanatorytimevarying.jl index abbd8ee..d36ac52 100644 --- a/src/models/locallevelexplanatorytimevarying.jl +++ b/src/models/locallevelexplanatorytimevarying.jl @@ -1,7 +1,7 @@ @doc raw""" LocalLevelExplanatoryTimeVarying(y::Vector{Fl}, X::Matrix{Fl}) where Fl -A local level model with time-varying regressors is defined by: +A local level model allowing time-varying regressors is defined by: ```math \begin{gather*} \begin{aligned} @@ -14,9 +14,10 @@ A local level model with time-varying regressors is defined by: \end{gather*} ``` +This model can quickly overfit, so it is useful to fix the variance of one or more of the tau hyperparamters at zero - with all fixed to zero, this model becomes the LocalLevelExplanatory # Example ```jldoctest -julia> model = tvr(rand(100), rand(100, 1)) +julia> model = LocalLevelExplanatoryTimeVarying(rand(100), rand(100, 1)) tvr ``` """ From 65f704a42cde0583788c2bb467193fd67d862939 Mon Sep 17 00:00:00 2001 From: Paul Mainwood Date: Tue, 29 Jun 2021 10:35:33 +0100 Subject: [PATCH 4/6] Add tests - just same as lle --- .../locallevelexplanatorytimevarying.jl | 46 +++++++++++++++++++ test/runtests.jl | 1 + 2 files changed, 47 insertions(+) create mode 100644 test/models/locallevelexplanatorytimevarying.jl diff --git a/test/models/locallevelexplanatorytimevarying.jl b/test/models/locallevelexplanatorytimevarying.jl new file mode 100644 index 0000000..56417be --- /dev/null +++ b/test/models/locallevelexplanatorytimevarying.jl @@ -0,0 +1,46 @@ +@testset "Local Level With Explanatory Variables Time Varying Model" begin + @test has_fit_methods(LocalLevelExplanatoryTimeVarying) + y = [ + 0.8297606696482265 + 0.5151262659173474 + 0.6578708499069135 + 0.15479984562134974 + 0.41836893137914033 + 0.46381576757801835 + 0.18373491281112986 + 0.09712212189897196 + 0.9696168042329114 + 0.8623172534114631 + 0.4312662075760265 + 0.17938380947382027 + 0.6265684534572271 + 0.6395165239895426 + 0.370025441453405 + ] + X = [ + 0.221409 0.646901 + 0.943735 0.163734 + 0.19864 0.620721 + 0.462932 0.876136 + 0.88467 0.673716 + 0.387496 0.415591 + 0.29221 0.989594 + 0.487926 0.721127 + 0.0321253 0.318292 + 0.742561 0.582234 + 0.664695 0.749821 + 0.791058 0.423212 + 0.0652515 0.955714 + 0.842513 0.788755 + 0.198105 0.79092 + ] + model = LocalLevelExplanatoryTimeVarying(y, X) + fit!(model) + + # forecasting + # For a fixed forecasting explanatory the variance must not decrease + forec = forecast(model, ones(5, 2)) + @test monotone_forecast_variance(forec) + kf = kalman_filter(model) + ks = kalman_smoother(model) +end diff --git a/test/runtests.jl b/test/runtests.jl index 8ce3c74..2e1a168 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,6 +18,7 @@ include("models/locallevel.jl") include("models/locallineartrend.jl") include("models/locallevelcycle.jl") include("models/locallevelexplanatory.jl") +include("models/locallevelexplanatorytimevarying.jl") include("models/basicstructural.jl") include("models/basicstructural_explanatory.jl") include("models/basicstructural_multivariate.jl") From 9c00531e9ed30dec859c05a967937890ee433bff Mon Sep 17 00:00:00 2001 From: Paul Mainwood Date: Tue, 29 Jun 2021 10:43:58 +0100 Subject: [PATCH 5/6] Fix docstring --- src/models/locallevelexplanatorytimevarying.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/models/locallevelexplanatorytimevarying.jl b/src/models/locallevelexplanatorytimevarying.jl index d36ac52..f9946a3 100644 --- a/src/models/locallevelexplanatorytimevarying.jl +++ b/src/models/locallevelexplanatorytimevarying.jl @@ -18,7 +18,7 @@ This model can quickly overfit, so it is useful to fix the variance of one or mo # Example ```jldoctest julia> model = LocalLevelExplanatoryTimeVarying(rand(100), rand(100, 1)) -tvr +LocalLevelExplanatoryTimeVarying ``` """ mutable struct LocalLevelExplanatoryTimeVarying <: StateSpaceModel From 5acb6c571dfbfdb3369e65478cf3713345cad0d8 Mon Sep 17 00:00:00 2001 From: Paul Mainwood Date: Mon, 6 Jun 2022 19:05:35 +0100 Subject: [PATCH 6/6] Update locallevelexplanatorytimevarying.jl --- src/models/locallevelexplanatorytimevarying.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/models/locallevelexplanatorytimevarying.jl b/src/models/locallevelexplanatorytimevarying.jl index f9946a3..9ef7b71 100644 --- a/src/models/locallevelexplanatorytimevarying.jl +++ b/src/models/locallevelexplanatorytimevarying.jl @@ -27,7 +27,7 @@ mutable struct LocalLevelExplanatoryTimeVarying <: StateSpaceModel results::Results exogenous::Matrix - function LocalLevelExplanatoryTimeVarying(y::Vector{Fl}, X::Matrix{Fl}) where Fl + function LocalLevelExplanatoryTimeVarying(y::Vector{Fl}, X::Matrix{Fl}; names_of_exogeneous = []) where Fl @assert length(y) == size(X, 1) @@ -47,7 +47,11 @@ mutable struct LocalLevelExplanatoryTimeVarying <: StateSpaceModel system = LinearUnivariateTimeVariant{Fl}(y, Z, T, R, d, c, H, Q) # Define hyperparameters names - names = [["sigma2_ε", "sigma2_η"];["β_$i" for i in 1:num_exogenous]; ["tau2_$i" for i in 1:num_exogenous]] + if isempty(names_of_exogeneous) + names = [["sigma2_ε", "sigma2_η"];["β_$i" for i in 1:num_exogenous]; ["tau2_$i" for i in 1:num_exogenous]] + else + names = [["sigma2_ε", "sigma2_η"];[name * "_β" for name in names_of_exogeneous]; [name * "_τ2" for name in names_of_exogeneous]] + end hyperparameters = HyperParameters{Fl}(names) return new(hyperparameters, system, Results{Fl}(), X)