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..9ef7b71 --- /dev/null +++ b/src/models/locallevelexplanatorytimevarying.jl @@ -0,0 +1,131 @@ +@doc raw""" + LocalLevelExplanatoryTimeVarying(y::Vector{Fl}, X::Matrix{Fl}) where Fl + +A local level model allowing 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*} +``` + +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 = LocalLevelExplanatoryTimeVarying(rand(100), rand(100, 1)) +LocalLevelExplanatoryTimeVarying +``` +""" +mutable struct LocalLevelExplanatoryTimeVarying <: StateSpaceModel + hyperparameters::HyperParameters + system::LinearUnivariateTimeVariant + results::Results + exogenous::Matrix + + function LocalLevelExplanatoryTimeVarying(y::Vector{Fl}, X::Matrix{Fl}; names_of_exogeneous = []) 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) 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, m, m) for _ in 1:num_observations] + + system = LinearUnivariateTimeVariant{Fl}(y, Z, T, R, d, c, H, Q) + + # Define hyperparameters names + 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) + 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) + #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}( + model.hyperparameters.names .=> 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) + 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) + 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] = diagm([get_constrained_value(model, "sigma2_η"); constrained_values_taus]) + 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 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")