Fitted quantiles and SDMForecast (#98)

* added fitted quantiles and SDMforecast

* update default verbose

* remove fitted quantiles

* update

Author: guilhermebodin

.gitignore

@@ -1,4 +1,5 @@
 Manifest.toml
 docs/build/
 .DS_Store
-docs/.DS_Store
+docs/.DS_Store
+paper_scripts/

examples/cpichg.jl

@@ -6,6 +6,6 @@ y = readdlm("./test/data/cpichg.csv")[:]
 gas = Model(1, 1, TDistLocationScale, 0.0, time_varying_params = [1; 2])
 
 # Estimate the model
-@time f = fit!(gas, y)
+f = fit!(gas, y)
 
 estimation_stats = fit_stats(f)

src/MLE.jl

@@ -2,7 +2,7 @@ export fit, fit!, fit_stats
 
 const DEFAULT_INITIAL_PARAM = NaN.*ones(1, 1)
 const DEFAULT_NUM_SEEDS = 3
-const DEFAULT_VERBOSE = 0
+const DEFAULT_VERBOSE = 1
 const VARIANCE_ZERO = 1e-10
 
 struct Fitted{D <: Distribution, T <: AbstractFloat}
@@ -96,6 +96,7 @@ function fit(gas::Model{D, T}, y::Vector{T};
              opt_method::AbstractOptimizationMethod = NelderMead(gas, DEFAULT_NUM_SEEDS),
              verbose::Int = DEFAULT_VERBOSE) where {D, T}
 
+    verbose in [0, 1, 2, 3] || throw(ErrorException, "verbose argument must be in [0, 1, 2, 3]")
     # Number of initial_points and number of params to estimate
     n_initial_points = length(opt_method.initial_points)
     n = length(y)
@@ -114,18 +115,20 @@ function fit(gas::Model{D, T}, y::Vector{T};
 
     for i = 1:n_initial_points
         try 
-            func = TwiceDifferentiable(psi_tilde -> log_lik(psi_tilde, y, gas_fit, initial_params, unknowns, n), opt_method.initial_points[i])
+            func = TwiceDifferentiable(psi_tilde -> log_lik(psi_tilde, y, gas_fit, 
+                                                        initial_params, unknowns, n), 
+                                                        opt_method.initial_points[i])
             opt_result = optimize(func, opt_method, verbose, i)
             update_aux_estimation!(aux_est, func, opt_result)
-            println("Round $i of $n_initial_points - Log-likelihood: $(-opt_result.minimum)")
+            verbose >= 1 && println("Round $i of $n_initial_points - Log-likelihood: $(-opt_result.minimum)")
         catch err
             println(err)
-            println("Round $i diverged")
+            verbose >= 1 && println("Round $i diverged")
         end
     end
 
     if isempty(aux_est.loglikelihood) 
-        println("No initial point converged.")
+        verbose >= 1 && println("No initial point converged.")
         return
     end
 
@@ -135,12 +138,11 @@ function fit(gas::Model{D, T}, y::Vector{T};
     aic = AIC(n_unknowns, best_llk)
     bic = BIC(n, n_unknowns, best_llk)
 
-    if verbose >= 1
+    if verbose >= 2
         println("\nBest optimization result:")
         println(aux_est.opt_result[best_seed])
     end
 
-    println("Finished!")
     return Fitted{D, T}(n, unknowns, aic, bic, best_llk, coefs, num_hessian)
 end
 

src/opt_methods/common_methods.jl

@@ -47,4 +47,4 @@ function optimize(func::Optim.TwiceDifferentiable, opt_method::AbstractOptimizat
                                                                    show_trace = show_trace(verbose) ))
 end
 
-show_trace(verbose::Int) =  (verbose == 2 ? true : false)
+show_trace(verbose::Int) =  (verbose == 3 ? true : false)

src/simulate.jl

@@ -1,5 +1,10 @@
 export simulate, forecast_quantiles
 
+mutable struct Forecast{T <: AbstractFloat}
+    quantiles::Matrix{T}
+    scenarios::Matrix{T}
+end
+
 """
     simulate(series::Vector{T}, gas::Model{D, T}, H::Int, S::Int, kwargs...) where {D, T}
 
@@ -14,17 +19,14 @@ function simulate(series::Vector{T}, gas::Model{D, T}, H::Int, S::Int;
                     initial_params::Matrix{T} = stationary_initial_params(gas)) where {D, T}
     # 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), :]
     # Create scenarios matrix
     scenarios = Matrix{T}(undef, H, 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]
     end
-
     return scenarios
 end
 
@@ -52,8 +54,7 @@ function forecast_quantiles(series::Vector{T}, gas::Model{D, T}, H::Int;
                       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 get_quantiles(quantiles, scenarios)
+    return Forecast(get_quantiles(quantiles, scenarios), scenarios)
 end
 
 function get_quantiles(quantile_probs::Vector{T}, scenarios::Matrix{T}) where T

src/univariate_score_driven_recursion.jl

@@ -1,7 +1,7 @@
 export score_driven_recursion, fitted_mean, simulate_recursion
 
 """
-score_driven_recursion(sd_model::SDM, observations::Vector{T}) where T
+    score_driven_recursion(sd_model::SDM, observations::Vector{T}) where T
 
 start with the stationary params for a 
 """

test/test_estimate.jl

@@ -7,15 +7,15 @@
         @testset "Estimation by passing number of initial_points" begin
             # LBFGS
             gas = Model(1, 1, Beta, 0.0)
-            fit!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
+            fit!(gas, simulation; verbose = 2, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
             test_coefficients_GAS_1_1(gas, ω, A, B)
             # NelderMead
             gas = Model(1, 1, Beta, 0.0)
-            fit!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.NelderMead(gas, 3))
+            fit!(gas, simulation; verbose = 2, opt_method = ScoreDrivenModels.NelderMead(gas, 3))
             test_coefficients_GAS_1_1(gas, ω, A, B)
             # IPNewton
             gas = Model(1, 1, Beta, 0.0)
-            fit!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.IPNewton(gas, 3))
+            fit!(gas, simulation; verbose = 2, opt_method = ScoreDrivenModels.IPNewton(gas, 3))
             test_coefficients_GAS_1_1(gas, ω, A, B)
         end
         @testset "Estimation by passing initial_points" begin

test/test_simulate.jl

@@ -15,9 +15,9 @@
         @test maximum(scenarios) < 1e5
         @test minimum(scenarios) > 0
 
-        quants = forecast_quantiles(y, gas, 50)
-        @test all([quants[i, 1] .< quants[i, 2] .< quants[i, 3] for i in 1:50])
-        @test maximum(quants) < 1e3
-        @test minimum(quants) > 1e-4
+        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
     end
 end

test/utils.jl

@@ -223,7 +223,7 @@ function test_GARCH_1_1(y, seed::Int, optimizer; atol = 1e-4, rtol = 1e-4)
     ub = [1.0; 1.0; 1.0; 1.0]
     lb = [-1.0; 0.0; 0.0; 0.0]
     gas = Model(1, 1, Normal, 1.0, time_varying_params = [2])
-    f = fit!(gas, y; initial_params = ini, verbose = 1, opt_method = optimizer(gas, 10; ub = ub, lb = lb))
+    f = fit!(gas, y; initial_params = ini, verbose = 2, opt_method = optimizer(gas, 10; ub = ub, lb = lb))
     show(stdout, fit_stats(f))
 
     @test gas.ω[1] - -0.00616637237701241 ≈ 0 atol = atol rtol = rtol