diff --git a/.gitignore b/.gitignore index 3f0ffe5..c8ce933 100644 --- a/.gitignore +++ b/.gitignore @@ -7,4 +7,7 @@ test_m3.jl test_dynamic_comb.jl inovations_ssmodels.jl test_component_expression.jl -results_variables_expressions.csv \ No newline at end of file +test_gas.jl + +# Files with csv extension +*.csv \ No newline at end of file diff --git a/src/UnobservedComponentsGAS.jl b/src/UnobservedComponentsGAS.jl index 4fd5383..d8be480 100644 --- a/src/UnobservedComponentsGAS.jl +++ b/src/UnobservedComponentsGAS.jl @@ -23,6 +23,7 @@ module UnobservedComponentsGAS include("distributions/normal.jl") include("distributions/t_location_scale.jl") include("distributions/log_normal.jl") + include("distributions/gamma.jl") include("initialization.jl") include("fit.jl") include("utils.jl") @@ -32,17 +33,22 @@ module UnobservedComponentsGAS include("residuals_diagnostics.jl") const DICT_CODE = Dict(1 => "Normal", - 2 => "tLocationScale" ) + 2 => "tLocationScale", + 3 => "Gamma") const DICT_SCORE = Dict("Normal" => score_normal, - "tLocationScale" => score_tlocationscale) + "tLocationScale" => score_tlocationscale, + "Gamma" => score_gamma) const DICT_FISHER_INFORMATION = Dict("Normal" => fisher_information_normal, - "tLocationScale" => fisher_information_tlocationscale) + "tLocationScale" => fisher_information_tlocationscale, + "Gamma" => fisher_information_gamma) const DICT_LOGPDF = Dict("Normal" => logpdf_normal, - "tLocationScale" => logpdf_tlocationscale) + "tLocationScale" => logpdf_tlocationscale, + "Gamma" => logpdf_gamma) const DICT_CDF = Dict("Normal" => cdf_normal, - "tLocationScale" => cdf_tlocationscale) + "tLocationScale" => cdf_tlocationscale, + "Gamma" => cdf_gamma) end diff --git a/src/components_dynamics.jl b/src/components_dynamics.jl index 579cc68..0340e0b 100644 --- a/src/components_dynamics.jl +++ b/src/components_dynamics.jl @@ -233,7 +233,7 @@ function add_AR!(model::Ml, s::Vector{Fl}, T::Int64, ar::Union{Int64, Vector{Int @variable(model, ϕ[1:max_order, idx_params]) @variable(model, κ_AR[idx_params]) - @constraint(model, [i in idx_params], -2 ≤ κ_AR[i] ≤ 2) + @constraint(model, [i in idx_params], 0 ≤ κ_AR[i] ≤ 2) # Revisar essas restrições com o Alves !! for i in unique_orders @@ -273,8 +273,8 @@ function add_random_walk_slope!(model::Ml, s::Vector{Fl}, T::Int64, random_walk_ @constraint(model, [t = 2:T, j in idx_params], b[t, j] == b[t - 1, j] + κ_b[j] * s[j][t]) @constraint(model, [t = 2:T, j in idx_params], RWS[t, j] == RWS[t - 1, j] + b[t - 1, j] + κ_RWS[j] * s[j][t]) - @constraint(model, [j in idx_params], -2 ≤ κ_RWS[j] ≤ 2) - @constraint(model, [j in idx_params], -2 ≤ κ_b[j] ≤ 2) + @constraint(model, [j in idx_params], 0 ≤ κ_RWS[j] ≤ 0) + @constraint(model, [j in idx_params], 0 ≤ κ_b[j] ≤ 0) end """ @@ -299,7 +299,7 @@ function add_random_walk!(model::Ml, s::Vector{Fl}, T::Int64, random_walk::Dict{ @variable(model, κ_RW[idx_params]) @constraint(model, [t = 2:T, j in idx_params], RW[t, j] == RW[t-1, j] + κ_RW[j] * s[j][t]) - @constraint(model, [j in idx_params], -2 ≤ κ_RW[j] ≤ 2) + @constraint(model, [j in idx_params], 0 ≤ κ_RW[j] ≤ 2) end """ @@ -325,7 +325,7 @@ function add_ar1!(model::Ml, s::Vector{Fl}, T::Int64, ar1::Dict{Int64, Bool}) w @variable(model, κ_AR1_LEVEL[idx_params]) @variable(model, ω_AR1_LEVEL[idx_params]) - @constraint(model, [i in idx_params], -2 ≤ κ_AR1_LEVEL[i] ≤ 2) + @constraint(model, [i in idx_params], 0 ≤ κ_AR1_LEVEL[i] ≤ 2) @constraint(model, [i in idx_params], -0.95 <= ϕ_AR1_LEVEL[i] <= 0.95) @constraint(model, [t = 2:T, j in idx_params], AR1_LEVEL[t, j] == ω_AR1_LEVEL[j] + ϕ_AR1_LEVEL[j] * AR1_LEVEL[t-1, j] + κ_AR1_LEVEL[j] * s[j][t]) @@ -435,7 +435,7 @@ function add_trigonometric_seasonality!(model::Ml, s::Vector{Fl}, T::Int64, seas @variable(model, κ_S[idx_params_stochastic]) - @constraint(model, [i in idx_params_stochastic], -2 ≤ κ_S[i] ≤ 2) + @constraint(model, [i in idx_params_stochastic], 0 ≤ κ_S[i] ≤ 2) #JuMP.fix.(model[:κ_S][idx_params_deterministic], 1e-4) @variable(model, γ_sto[1:unique_num_harmonic, 1:T, idx_params_stochastic]) diff --git a/src/distributions/gamma.jl b/src/distributions/gamma.jl new file mode 100644 index 0000000..283a0fd --- /dev/null +++ b/src/distributions/gamma.jl @@ -0,0 +1,203 @@ +" +Defines a Gamma distribution with parameters α λ. +From a shape (α) and ratio (β) parametrization, we obtain our parametrization making λ = α/β +" +mutable struct GammaDistribution <: ScoreDrivenDistribution + λ::Union{Missing, Float64} + α::Union{Missing, Float64} +end + +" +Outer constructor for the Normal distribution. +" +function GammaDistribution() + return GammaDistribution(missing, missing) +end + +" +Gamma Function Γ(x) +" +function Γ(x) + return SpecialFunctions.gamma(x) +end + + +"Auxiliar function Ψ1(α)" +function Ψ1(α) + return SpecialFunctions.digamma(α) +end + + +"Auxiliar function Ψ2(α)" +function Ψ2(α) + return SpecialFunctions.trigamma(α) +end + + +" +Evaluate the score of a Normal distribution with mean μ and variance σ², in observation y. +Colocar link aqui +" +function score_gamma(λ, α, y) + + α <= 0 ? α = 1e-2 : nothing + λ <= 0 ? λ = 1e-4 : nothing + + ∇_α = log(y) - y/λ + log(α) - Ψ1(α) - log(λ) + 1 + ∇_λ = (α/λ)*((y/λ)-1) + + return [∇_λ; ∇_α] +end + +" +Evaluate the fisher information of a Normal distribution with mean μ and variance σ². +Colocar link aqui +" +function fisher_information_gamma(λ, α) + + α <= 0 ? α = 1e-2 : nothing + λ <= 0 ? λ = 1e-4 : nothing + + I_λ = α/(λ^2) + I_α = Ψ2(α) - 1/α + + return [I_λ 0; 0 I_α] +end + +" +Evaluate the log pdf of a Normal distribution with mean μ and variance σ², in observation y. +" +function logpdf_gamma(λ, α, y) + + return logpdf_gamma([λ, α], y) +end + +" +Evaluate the log pdf of a Normal distribution with mean μ and variance σ², in observation y. + param[1] = α + param[2] = λ +" +function logpdf_gamma(param, y) + + param[2] > 0 ? α = param[2] : α = 1e-2 + param[1] > 0 ? λ = param[1] : λ = 1e-2 + + return Distributions.logpdf(Distributions.Gamma(α, λ/α), y) +end + +" +Evaluate the cdf of a Gamma distribution with α,λ, in observation y. +" +function cdf_gamma(param::Vector{Float64}, y::Fl) where Fl + + param[2] > 0 ? α = param[2] : α = 1e-2 + param[1] > 0 ? λ = param[1] : λ = 1e-2 + + return Distributions.cdf(Distributions.Gamma(α, λ/α), y) +end + +" +Returns the code of the Normal distribution. Is the key of DICT_CODE. +" +function get_dist_code(dist::GammaDistribution) + return 3 +end + +" +Returns the number of parameters of the Normal distribution. +" +function get_num_params(dist::GammaDistribution) + return 2 +end + +" +Simulates a value from a given Normal distribution. + param[1] = λ + param[2] = α +" +function sample_dist(param::Vector{Float64}, dist::GammaDistribution) + + "A Gamma do pacote Distributions é parametrizada com shape α e scale θ" + "Como θ = 1/β e β = α/λ, segue-se que θ = λ/α" + param[2] > 0 ? α = param[2] : α = 1e-2 + param[1] > 0 ? λ = param[1] : λ = 1e-2 + + return rand(Distributions.Gamma(α, λ/α)) +end + +" +Indicates which parameters of the Normal distribution must be positive. +" +function check_positive_constrainst(dist::GammaDistribution) + return [true, true] +end + + +function get_initial_α(y::Vector{Float64}) + + T = length(y) + model = JuMP.Model(Ipopt.Optimizer) + set_silent(model) + @variable(model, α >= 1e-2) # Ensure α is positive + @variable(model, λ[1:T] .>= 1e-4) + register(model, :Γ, 1, Γ; autodiff = true) + @NLobjective(model, Max, sum(-log(Γ(α)) - α*log(1/α) - α*log(λ[i]) +(α-1)*log(y[i]) - (α/λ[i])*y[i] for i in 1:T)) + optimize!(model) + if has_values(model) + return JuMP.value.(α) + else + return fit_mle(Gamma, y).α + end +end + +" +Returns a dictionary with the initial values of the parameters of Normal distribution that will be used in the model initialization. +" +function get_initial_params(y::Vector{Fl}, time_varying_params::Vector{Bool}, dist::GammaDistribution, seasonality::Dict{Int64, Union{Bool, Int64}}) where Fl + + T = length(y) + dist_code = get_dist_code(dist) + seasonal_period = get_num_harmonic_and_seasonal_period(seasonality)[2] + + initial_params = Dict() + fitted_distribution = fit_mle(Gamma, y) + + # param[2] = λ = média + if time_varying_params[1] + initial_params[1] = y + else + initial_params[1] = fitted_distribution.α*fitted_distribution.θ + end + + # param[1] = α + if time_varying_params[2] + initial_params[2] = get_seasonal_var(y, maximum(seasonal_period), dist)#(scaled_score.(y, ones(T) * var(diff(y)) , y, 0.5, dist_code, 2)).^2 + else + initial_params[2] = get_initial_α(y)#mean(y)^2/var((y)) + end + + return initial_params +end + + +function get_seasonal_var(y::Vector{Fl}, seasonal_period::Int64, dist::GammaDistribution) where Fl + + num_periods = ceil(Int, length(y) / seasonal_period) + seasonal_variances = zeros(Fl, length(y)) + + for i in 1:seasonal_period + month_data = y[i:seasonal_period:end] + num_observations = length(month_data) + if num_observations > 0 + g = Distributions.fit(Gamma, month_data) + α = g.α + θ = g.θ + variance = α*(θ^2) + for j in 1:num_observations + seasonal_variances[i + (j - 1) * seasonal_period] = variance + end + end + end + return seasonal_variances +end + \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 25d73e4..a77437a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,6 +8,7 @@ using JSON3, CSV, DataFrames, StateSpaceModels include("test_fit_forecast_normal.jl") include("test_fit_forecast_lognormal.jl") include("test_fit_forecast_t.jl") +include("test_fit_forecast_gamma.jl") include("test_components_dynamics.jl") include("test_optimization.jl") include("test_distributions.jl") diff --git a/test/test_fit_forecast_gamma.jl b/test/test_fit_forecast_gamma.jl new file mode 100644 index 0000000..7af2009 --- /dev/null +++ b/test/test_fit_forecast_gamma.jl @@ -0,0 +1,199 @@ +@testset "Fit & Forecast Gamma" begin + + time_series = CSV.read(joinpath(@__DIR__, "data/timeseries_lognormal_rws_d1.csv"), DataFrame) + T = size(time_series, 1) + y = time_series[:,2] + X = [2*y y/2 rand(T)] + + steps_ahead = 12 + num_scenarious = 500 + X_gamma_forec = hcat(y[end-steps_ahead+1:end].+5*rand(steps_ahead), + y[end-steps_ahead+1:end].+10*rand(steps_ahead), + rand(steps_ahead)) + + function test_initial_values_components(initial_values, rw, rws, ar, seasonality) + ismissing(seasonality) ? s = false : s = true + ar == false ? ar_bool = false : ar_bool = true + dict_has_component = Dict("rws" => rws, "rw" => rw, "ar" => ar_bool, "seasonality" => s) + + dict_tests = Dict() + for component in keys(initial_values) + if !occursin("param", component) && component != "explanatories" + dict_tests[component] = !all(iszero.(initial_values[component]["values"])) + elseif occursin("param", component) && component != "explanatories" + dict_tests[component] = !all(iszero.(initial_values[component])) + end + end + + return all([dict_tests["rw"] == rw, dict_tests["rws"] == rws, dict_tests["ar"] == ar_bool, dict_tests["seasonality"] == s]) + end + + function convert_dict_keys_to_string(dict)::Dict{String, Any} + new_dict = Dict{String, Any}() + for (key, value) in dict + new_dict[string(key)] = value + end + return new_dict + end + + rw = false + rws = true + ar = false + seasonality = 12 + + dist_gamma = UnobservedComponentsGAS.GammaDistribution() + + gas_model_gamma = UnobservedComponentsGAS.GASModel(dist_gamma, [true, false], 1.0, "random walk slope", "deterministic 12", missing) + gas_model_gamma_2params = UnobservedComponentsGAS.GASModel(dist_gamma, [true, true], 1.0, ["random walk slope", "random walk"], + ["deterministic 12", "deterministic 12"], [missing, missing]) + + gas_model_gamma_X = deepcopy(gas_model_gamma) + gas_model_gamma_X_2params = deepcopy(gas_model_gamma_2params) + + model_gamma, parameters_gamma, initial_values_gamma = UnobservedComponentsGAS.create_model(gas_model_gamma, y, missing) + model_gamma_2params, parameters_gamma_2params, initial_values_gamma_2params = UnobservedComponentsGAS.create_model(gas_model_gamma_2params, y, missing) + + model_gamma_X, parameters_gamma_X, initial_values_gamma_X = UnobservedComponentsGAS.create_model(gas_model_gamma_X, y, X, missing); + model_gamma_X_2params, parameters_gamma_X_2params, initial_values_gamma_X_2params = UnobservedComponentsGAS.create_model(gas_model_gamma_X_2params, y, X, missing); + + fitted_model_gamma = UnobservedComponentsGAS.fit(gas_model_gamma, y; tol = 5e-2) + fitted_model_gamma_2params = UnobservedComponentsGAS.fit(gas_model_gamma_2params, y; tol = 5e-2) + fitted_model_gamma_X = UnobservedComponentsGAS.fit(gas_model_gamma_X, y, X; tol = 5e-2) + fitted_model_gamma_X_2params = UnobservedComponentsGAS.fit(gas_model_gamma_X_2params, y, X; tol = 5e-2) + + forecast_gamma = UnobservedComponentsGAS.predict(gas_model_gamma, fitted_model_gamma, y, steps_ahead, num_scenarious) + forecast_gamma_X = UnobservedComponentsGAS.predict(gas_model_gamma_X, fitted_model_gamma_X, y, X_gamma_forec, steps_ahead, num_scenarious) + forecast_gamma_2params = UnobservedComponentsGAS.predict(gas_model_gamma_2params, fitted_model_gamma_2params, y, steps_ahead, num_scenarious) + + scenarios_gamma = UnobservedComponentsGAS.simulate(gas_model_gamma, fitted_model_gamma, y, steps_ahead, num_scenarious) + scenarios_gamma_X = UnobservedComponentsGAS.simulate(gas_model_gamma_X, fitted_model_gamma_X, y, X_gamma_forec, steps_ahead, num_scenarious) + scenarios_gamma_2params = UnobservedComponentsGAS.simulate(gas_model_gamma_2params, fitted_model_gamma_2params, y, steps_ahead, num_scenarious) + + @testset "create_model" begin + # Create model with no explanatory series + + @test(size(parameters_gamma) == (T,2)) + @test(size(parameters_gamma_2params) == (T,2)) + @test(typeof(model_gamma) == JuMP.Model) + @test(typeof(model_gamma_2params) == JuMP.Model) + @test(test_initial_values_components(initial_values_gamma, rw, rws, ar, seasonality)) + @test(test_initial_values_components(initial_values_gamma_2params, rw, rws, ar, seasonality)) + + @test(size(parameters_gamma_X) == (T,2)) + @test(size(parameters_gamma_X_2params) == (T,2)) + @test(typeof(model_gamma_X) == JuMP.Model) + @test(typeof(model_gamma_X_2params) == JuMP.Model) + @test(test_initial_values_components(initial_values_gamma_X, rw, rws, ar, seasonality)) + @test(test_initial_values_components(initial_values_gamma_X_2params, rw, rws, ar, seasonality)) + end + + #@info(" --- Testing fit functions") + @testset "fit" begin + # "Test if termination_status is correct" + possible_status = ["LOCALLY_SOLVED", "TIME_LIMIT"] + @test(fitted_model_gamma.model_status in possible_status) + @test(fitted_model_gamma_2params.model_status in possible_status) + @test(fitted_model_gamma_X.model_status in possible_status) + # "Test if selected_variables is missing " + @test(ismissing(fitted_model_gamma.selected_variables)) + @test(ismissing(fitted_model_gamma_2params.selected_variables)) + # "Test if fitted_params has the right keys -> order may be a problem" + @test(all(keys(fitted_model_gamma.fitted_params) .== ["param_2","param_1"])) + @test(all(keys(fitted_model_gamma_2params.fitted_params) .== ["param_2","param_1"])) + + # "Test if all time varying and fixed params are time varying and fixed" + @test(!all(y->y==fitted_model_gamma.fitted_params["param_1"][1],fitted_model_gamma.fitted_params["param_1"])) + @test(all(y->y==fitted_model_gamma.fitted_params["param_2"][1],fitted_model_gamma.fitted_params["param_2"])) + @test(!all(y->y==fitted_model_gamma_2params.fitted_params["param_1"][1],fitted_model_gamma_2params.fitted_params["param_1"])) + @test(!all(y->y==fitted_model_gamma_2params.fitted_params["param_2"][1],fitted_model_gamma_2params.fitted_params["param_2"])) + @test(!all(y->y==fitted_model_gamma_X.fitted_params["param_1"][1],fitted_model_gamma_X.fitted_params["param_1"])) + @test(all(y->y==fitted_model_gamma_X.fitted_params["param_2"][1],fitted_model_gamma_X.fitted_params["param_2"])) + @test(!all(y->y==fitted_model_gamma_X_2params.fitted_params["param_1"][1],fitted_model_gamma_X_2params.fitted_params["param_1"])) + @test(!all(y->y==fitted_model_gamma_X_2params.fitted_params["param_2"][1],fitted_model_gamma_X_2params.fitted_params["param_2"])) + + # "Test if all residuals are being generated" + residuals_types = ["q_residuals", "std_residuals", "cs_residuals"] + @test(all(keys(fitted_model_gamma.residuals) .== residuals_types)) + @test(all(keys(fitted_model_gamma_2params.residuals) .== residuals_types)) + end + + #@info(" --- Test forecast function ---") + @testset "forecast" begin + + @test(isapprox(forecast_gamma["mean"], vec(mean(scenarios_gamma, dims = 2)); rtol = 1e-3)) + @test(size(scenarios_gamma) == (steps_ahead, num_scenarious)) + + @test(isapprox(forecast_gamma["intervals"]["80"]["lower"], [quantile(scenarios_gamma[t,:], 0.2/2) for t in 1:steps_ahead])) + @test(isapprox(forecast_gamma["intervals"]["80"]["upper"], [quantile(scenarios_gamma[t,:], 1 - 0.2/2) for t in 1:steps_ahead])) + @test(isapprox(forecast_gamma["intervals"]["95"]["lower"], [quantile(scenarios_gamma[t,:], 0.05/2) for t in 1:steps_ahead])) + @test(isapprox(forecast_gamma["intervals"]["95"]["upper"], [quantile(scenarios_gamma[t,:], 1 - 0.05/2) for t in 1:steps_ahead])) + + @test(isapprox(forecast_gamma_2params["mean"], vec(mean(scenarios_gamma_2params, dims = 2)); rtol = 1e-3)) + @test(size(scenarios_gamma_2params) == (steps_ahead, num_scenarious)) + + @test(isapprox(forecast_gamma_2params["intervals"]["80"]["lower"], [quantile(scenarios_gamma_2params[t,:], 0.2/2) for t in 1:steps_ahead])) + @test(isapprox(forecast_gamma_2params["intervals"]["80"]["upper"], [quantile(scenarios_gamma_2params[t,:], 1 - 0.2/2) for t in 1:steps_ahead])) + @test(isapprox(forecast_gamma_2params["intervals"]["95"]["lower"], [quantile(scenarios_gamma_2params[t,:], 0.05/2) for t in 1:steps_ahead])) + @test(isapprox(forecast_gamma_2params["intervals"]["95"]["upper"], [quantile(scenarios_gamma_2params[t,:], 1 - 0.05/2) for t in 1:steps_ahead])) + + @test(isapprox(forecast_gamma_X["mean"], vec(mean(scenarios_gamma_X, dims = 2)); rtol = 1e-3)) + @test(size(scenarios_gamma_X) == (steps_ahead, num_scenarious)) + + @test(isapprox(forecast_gamma_X["intervals"]["80"]["lower"], [quantile(scenarios_gamma_X[t,:], 0.2/2) for t in 1:steps_ahead])) + @test(isapprox(forecast_gamma_X["intervals"]["80"]["upper"], [quantile(scenarios_gamma_X[t,:], 1 - 0.2/2) for t in 1:steps_ahead])) + @test(isapprox(forecast_gamma_X["intervals"]["95"]["lower"], [quantile(scenarios_gamma_X[t,:], 0.05/2) for t in 1:steps_ahead])) + @test(isapprox(forecast_gamma_X["intervals"]["95"]["upper"], [quantile(scenarios_gamma_X[t,:], 1 - 0.05/2) for t in 1:steps_ahead])) + end + + #@info(" --- Test quality of fit and forecast - Gamma") + @testset "quality of fit - Gamma" begin + y = time_series[1:end-steps_ahead,5] + y_test = time_series[end-steps_ahead+1:end, 5] + + gas_model = UnobservedComponentsGAS.GASModel(UnobservedComponentsGAS.GammaDistribution(), [true, false], + 1.0, "random walk slope", "deterministic 12", missing) + fitted_model = UnobservedComponentsGAS.fit(gas_model, y) + forec = UnobservedComponentsGAS.predict(gas_model, fitted_model, y, steps_ahead, num_scenarious) + + @test(isapprox(fitted_model.fit_in_sample[2:end], y[2:end]; rtol = 1e0)) + @test(isapprox(forec["mean"], y_test; rtol = 1e0)) + end + + #@info(" --- Test quality of fit - Gamma with 2 params") + @testset "quality of fit - Gamma with 2 params" begin + y = time_series[1:end-steps_ahead,5] + y_test = time_series[end-steps_ahead+1:end, 5] + gas_model = UnobservedComponentsGAS.GASModel(UnobservedComponentsGAS.GammaDistribution(), [true, true], + 0.0, ["random walk slope", "random walk"], ["deterministic 12", "deterministic 12"], [missing, missing]) + fitted_model = UnobservedComponentsGAS.fit(gas_model, y) + forec = UnobservedComponentsGAS.predict(gas_model, fitted_model, y, steps_ahead, num_scenarious) + + @test(isapprox(fitted_model.fit_in_sample[2:end], y[2:end]; rtol = 1e0)) + @test(isapprox(forec["mean"], y_test; rtol = 1e0)) + end + + # #@info(" --- Test quality of fit - Gamma with robust") + @testset "quality of fit - Gamma with robust" begin + y = time_series[1:end-steps_ahead,5] + y_test = time_series[end-steps_ahead+1:end, 5] + gas_model = UnobservedComponentsGAS.GASModel(UnobservedComponentsGAS.GammaDistribution(), [true, false], + 1.0, "random walk slope", "deterministic 12", 1) + fitted_model = UnobservedComponentsGAS.fit(gas_model, y; α = 0.0, robust = true) + forec = UnobservedComponentsGAS.predict(gas_model, fitted_model, y, steps_ahead, num_scenarious) + + @test(isapprox(fitted_model.fit_in_sample[2:end], y[2:end]; rtol = 1e0)) + @test(isapprox(forec["mean"], y_test; rtol = 1e0)) + end + + @testset "AR(1) level" begin + y = time_series[1:end-steps_ahead,5] + y_test = time_series[end-steps_ahead+1:end, 5] + gas_model = UnobservedComponentsGAS.GASModel(UnobservedComponentsGAS.GammaDistribution(), [true, false], + 0.5, ["ar(1)", ""], ["deterministic 12", ""], [missing, missing]) + fitted_model = UnobservedComponentsGAS.fit(gas_model, y) + + @test(all(fitted_model.components["param_1"]["level"]["value"] .!= zeros(length(y)))) + end + + +end \ No newline at end of file