Local level explanatory time varying

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

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

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 

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

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