Add parameters forecast (#100)

* Remove Chi and Chisq, add parameters to forecasting

* add forecast

* tests working

* address sugestions

Author: guilhermebodin

docs/src/manual.md

@@ -93,8 +93,6 @@ this is handled internally.
 | :------- |:---------------:|:--------------------:|
 |[`Beta`](@ref)| ✓ | ✓  |
 |[`BetaLocationScale`](@ref)| ✓ |  x  |
-|[`Chi`](@ref)| ✓ | ✓  |
-|[`Chisq`](@ref)| ✓ | ✓  |
 |[`Exponential`](@ref)| ✓ | ✓  |
 |[`Gamma`](@ref)| ✓ | ✓  |
 |[`LogitNormal`](@ref)| ✓ | ✓  |
@@ -110,8 +108,6 @@ this is handled internally.
 ```@docs
 ScoreDrivenModels.Beta
 ScoreDrivenModels.BetaLocationScale
-ScoreDrivenModels.Chi
-ScoreDrivenModels.Chisq
 ScoreDrivenModels.Exponential
 ScoreDrivenModels.Gamma
 ScoreDrivenModels.LogitNormal

src/ScoreDrivenModels.jl

@@ -33,8 +33,6 @@ include("distributions/non_native_dists.jl")
 include("distributions/common_interface.jl")
 include("distributions/beta.jl")
 include("distributions/betalocationscale.jl")
-include("distributions/chi.jl")
-include("distributions/chisq.jl")
 include("distributions/exponential.jl")
 include("distributions/gamma.jl")
 include("distributions/logitnormal.jl")

src/distributions/chi.jl

@@ -2,8 +2,6 @@
 const DISTS = [
     Beta,
     BetaLocationScale,
-    Chi,
-    Chisq,
     Exponential,
     Gamma,
     LogitNormal,

src/distributions/chisq.jl

@@ -1,8 +1,12 @@
-export simulate, forecast_quantiles
+export simulate, forecast
 
 mutable struct Forecast{T <: AbstractFloat}
-    quantiles::Matrix{T}
-    scenarios::Matrix{T}
+    observation_quantiles::Matrix{T}
+    observation_forecast::Vector{T}
+    observation_scenarios::Matrix{T}
+    parameter_quantiles::Array{T, 3}
+    parameter_forecast::Matrix{T}
+    parameter_scenarios::Array{T, 3}
 end
 
 """
@@ -20,19 +24,29 @@ function simulate(series::Vector{T}, gas::Model{D, T}, H::Int, S::Int;
     # Filter params estimated on the time series
     params = score_driven_recursion(gas, series; initial_params = initial_params)
     biggest_lag = number_of_lags(gas)
-    params_simulation = params[(end - biggest_lag):(end - 1), :]
+    n_params = num_params(D)
+    params_simulation = params[(end - biggest_lag + 1):end, :]
     # Create scenarios matrix
-    scenarios = Matrix{T}(undef, H, S)
+    observation_scenarios = Matrix{T}(undef, H, S)
+    parameter_scenarios = Array{T, 3}(undef, H, n_params, S)
     for scenario in 1:S
-        sim, param = simulate_recursion(gas, H + biggest_lag + 1; initial_params = params_simulation)
-        scenarios[:, scenario] = sim[biggest_lag + 2:end]
+        # Notice that we know the parameter time_varying parameters for T + 1 
+        # So the last initial_params is already a part of the future simulation
+        # And we must take the 
+        sim, param = simulate_recursion(gas, H + biggest_lag; initial_params = params_simulation)
+        observation_scenarios[:, scenario] = sim[biggest_lag+1:end]
+        # The first param is known
+        parameter_scenarios[1, :, scenario] = params_simulation[end, :]
+        # The last param (param[end, :]) is actually in a time step bigger 
+        # than H
+        parameter_scenarios[2:end, :, scenario] = param[biggest_lag+1:end-1, :]
     end
-    return scenarios
+    return observation_scenarios, parameter_scenarios
 end
 
 
 """
-    forecast_quantiles(series::Vector{T}, gas::Model{D, T}, H::Int; kwargs...) where {D, T}
+    forecast(series::Vector{T}, gas::Model{D, T}, H::Int; kwargs...) where {D, T}
 
 Forecast quantiles for future values of a time series by updating the GAS recursion `H` times and 
 using Monte Carlo method as in Blasques, Francisco, Siem Jan Koopman,
@@ -49,17 +63,40 @@ By default this method uses the `stationary_initial_params` method to perform th
 score driven recursion. If you estimated the model with a different set of `initial_params`
 use them here to maintain the coherence of your estimation.
 """
-function forecast_quantiles(series::Vector{T}, gas::Model{D, T}, H::Int;
-                      initial_params::Matrix{T} = stationary_initial_params(gas),
-                      quantiles::Vector{T} = T.([0.025, 0.5, 0.975]), S::Int = 10_000) where {D, T}
+function forecast(series::Vector{T}, gas::Model{D, T}, H::Int;
+                    initial_params::Matrix{T} = stationary_initial_params(gas),
+                    quantiles::Vector{T} = T.([0.025, 0.5, 0.975]), S::Int = 10_000) where {D, T}
 
-    scenarios = simulate(series, gas, H, S; initial_params = initial_params)
-    return Forecast(get_quantiles(quantiles, scenarios), scenarios)
+    observation_scenarios, parameter_scenarios = simulate(series, gas, H, S; 
+                                                    initial_params = initial_params)
+    return Forecast(
+        get_quantiles(quantiles, observation_scenarios), 
+        mean(observation_scenarios, dims = 2)[:],
+        observation_scenarios,
+        get_quantiles(quantiles, parameter_scenarios), 
+        mean(parameter_scenarios, dims = 3)[:, :, 1],
+        parameter_scenarios
+        )
 end
 
-function get_quantiles(quantile_probs::Vector{T}, scenarios::Matrix{T}) where T
-    @assert all((quantile_probs .< 1.0) .& (quantile_probs .> 0.0))
-    unique!(sort!(quantile_probs))
-    quantiles = mapslices(x -> quantile(x, quantile_probs), scenarios; dims = 2)
+function get_quantiles(quantiles_probs::Vector{T}, scenarios::Matrix{T}) where T
+    @assert all((quantiles_probs .< 1.0) .& (quantiles_probs .> 0.0))
+    unique!(sort!(quantiles_probs))
+    quantiles = mapslices(x -> quantile(x, quantiles_probs), scenarios; dims = 2)
     return quantiles
 end
+
+function get_quantiles(quantiles_probs::Vector{T}, scenarios::Array{T, 3}) where T
+    @assert all((quantiles_probs .< 1.0) .& (quantiles_probs .> 0.0))
+    unique!(sort!(quantiles_probs))
+    quantiles_per_parameter = Array{T, 3}(undef, 
+                                size(scenarios, 1),
+                                size(scenarios, 2),
+                                length(quantiles_probs)
+                                )
+    for j in 1:size(scenarios, 2)
+        quantiles_per_parameter[:, j, :] = mapslices(x -> quantile(x, quantiles_probs), 
+                                                        scenarios[:, j, :]; dims = 2)
+    end
+    return quantiles_per_parameter
+end

src/distributions/common_interface.jl

@@ -11,13 +11,23 @@
         gas.B[1]  = convert(Matrix{Float64}, Diagonal(3*v))
         gas.B[12] = convert(Matrix{Float64}, Diagonal(-3*v))
 
-        scenarios = simulate(y, gas, 50, 1000)
-        @test maximum(scenarios) < 1e5
-        @test minimum(scenarios) > 0
+        scenarios_observations, scenarios_params = simulate(y, gas, 50, 1000)
+        @test maximum(scenarios_observations) < 1e5
+        @test minimum(scenarios_observations) > 0
 
-        forecast = forecast_quantiles(y, gas, 50)
-        @test all([forecast.quantiles[i, 1] .< forecast.quantiles[i, 2] .< forecast.quantiles[i, 3] for i in 1:50])
-        @test maximum(forecast.quantiles) < 1e3
-        @test minimum(forecast.quantiles) > 1e-4
+        forec = forecast(y, gas, 50)
+        @test all([forec.observation_quantiles[i, 1] .< 
+                   forec.observation_quantiles[i, 2] .< 
+                   forec.observation_quantiles[i, 3] for i in 1:50])
+        @test all([forec.parameter_quantiles[i, 1, 1] <=
+                   forec.parameter_quantiles[i, 1, 2] <= 
+                   forec.parameter_quantiles[i, 1, 3] for i in 1:50])
+        @test all([forec.parameter_quantiles[i, 2, 1] <=
+                   forec.parameter_quantiles[i, 2, 2] <= 
+                   forec.parameter_quantiles[i, 2, 3] for i in 1:50])
+        @test maximum(forec.observation_quantiles) < 1e3
+        @test minimum(forec.observation_quantiles) > 1e-4
+
+        @test minimum(forec.parameter_scenarios[:, 2, :]) > 0.0
     end
 end