Viz backtest (#130)

* Add backtest recipes

* update main file and tests

* Update scripts

Author: guilhermebodin

Project.toml

@@ -1,7 +1,7 @@
 name = "ScoreDrivenModels"
 uuid = "4a87933e-d659-11e9-0e65-7f40dedd4a3a"
 authors = ["guilhermebodin <[email protected]>, raphaelsaavedra <[email protected]>"]
-version = "0.1.7"
+version = "0.1.8"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

examples/ane.jl

@@ -9,8 +9,8 @@ y_test = y[401:460]
 # Set RNG seed to guarantee consistent results
 Random.seed!(123)
 
-# Specify GAS model: a lognormal model with time-varying μ, constant σ, and lags 1 and 12
-gas = Model([1, 2, 11, 12], [1, 2, 11, 12], LogNormal, 0.0; time_varying_params=[1])
+# Specify GAS model: a lognormal model with time-varying μ, constant σ, and lags 4 and 12
+gas = Model(4, 12, LogNormal, 1.0; time_varying_params=[1])
 
 # Obtain initial parameters to start the GAS recursion
 initial_params = dynamic_initial_params(y_train, gas)
@@ -20,6 +20,7 @@ f = ScoreDrivenModels.fit!(gas, y_train; initial_params=initial_params)
 
 # Print estimation statistics
 results(f)
+plot(f)
 
 # Simulate 1000 future scenarios and obtain the 5% and 95% quantiles in each time period
 forec = forecast(y_train, gas, 60; S=1000, initial_params=initial_params)

examples/cpichg.jl

@@ -28,6 +28,7 @@ f = fit!(gas, y; opt_method=LBFGS(gas, 5))
 
 # Print estimation statistics
 results(f)
+plot(f)
 
 # Perform forecast via simulations for 12 time periods ahead
 forec = forecast(y, gas, 12)
@@ -37,3 +38,11 @@ forec.parameter_forecast
 
 # Similarly, we can access the simulated observation scenarios
 forec.observation_scenarios
+
+gas_t = Model(1, 1, TDistLocationScale, 0.0; time_varying_params = [1])
+steps_ahead = 50
+first_idx = 150
+b_t = backtest(gas_t, y, steps_ahead, first_idx)
+plot(b_t, "GAS(1, 1) Student t")
+using Plots
+plot!(b_t, "GAS(1, 1) Student t")

src/ScoreDrivenModels.jl

@@ -47,5 +47,6 @@ include("distributions/weibull.jl")
 
 include("visualization/forecast.jl")
 include("visualization/residuals.jl")
+include("visualization/backtest.jl")
 
 end

src/backtest.jl

@@ -2,11 +2,15 @@ export backtest
 
 struct Backtest
     abs_errors::Matrix{Float64}
+    mae::Vector{Float64}
     crps_scores::Matrix{Float64}
+    mean_crps::Vector{Float64}
     function Backtest(n::Int, steps_ahead::Int)
-        abs_errors = Matrix{Float64}(undef, n, steps_ahead)
-        crps_scores = Matrix{Float64}(undef, n, steps_ahead)
-        return new(abs_errors, crps_scores)
+        abs_errors = Matrix{Float64}(undef, steps_ahead, n)
+        crps_scores = Matrix{Float64}(undef, steps_ahead, n)
+        mae = Vector{Float64}(undef, steps_ahead)
+        mean_crps = Vector{Float64}(undef, steps_ahead)
+        return new(abs_errors, mae, crps_scores, mean_crps)
     end
 end
 
@@ -36,8 +40,8 @@ TODO
 """
 function backtest(gas::Model{<:Distribution, T}, y::Vector{T}, steps_ahead::Int, start_idx::Int;
                   S::Int = 10_000,
-                  initial_params = stationary_initial_params(gas),
-                  opt_method = NelderMead(gas, 3)) where T
+                  initial_params::Matrix{T} = stationary_initial_params(gas),
+                  opt_method = NelderMead(gas, DEFAULT_NUM_SEEDS)) where T
     num_mle = length(y) - start_idx - steps_ahead
     b = Backtest(num_mle, steps_ahead)
     for i in 1:num_mle
@@ -49,8 +53,12 @@ function backtest(gas::Model{<:Distribution, T}, y::Vector{T}, steps_ahead::Int,
         forec = forecast(y_to_fit, gas_to_fit, steps_ahead; S=S, initial_params=initial_params)
         abs_errors = evaluate_abs_error(y_to_verify, forec.observation_forecast)
         crps_scores = evaluate_crps(y_to_verify, forec.observation_scenarios)
-        b.abs_errors[i, :] = abs_errors
-        b.crps_scores[i, :] = crps_scores
+        b.abs_errors[:, i] = abs_errors
+        b.crps_scores[:, i] = crps_scores
+    end
+    for i in 1:steps_ahead
+        b.mae[i] = mean(b.abs_errors[i, :])
+        b.mean_crps[i] = mean(b.crps_scores[i, :])
     end
     return b
 end

src/visualization/backtest.jl

@@ -0,0 +1,15 @@
+RecipesBase.@recipe function f(b::ScoreDrivenModels.Backtest, name::String)
+    xguide := "lead times"
+    @series begin
+        seriestype := :path
+        label := "MAE " * name
+        marker := :circle
+        b.mae
+    end
+    @series begin
+        seriestype := :path
+        label := "Mean CRPS " * name
+        marker := :circle
+        b.mean_crps
+    end
+end

src/visualization/residuals.jl

@@ -1,11 +1,12 @@
-RecipesBase.@recipe function f(fi::ScoreDrivenModels.Fitted)
+RecipesBase.@recipe function f(fi::Fitted)
     layout := (2, 2)
     acf = autocor(fi.pearson_residuals)[2:end]
     @series begin
         seriestype := :path
         label := ""
         seriescolor := "black"
         subplot := 1
+        marker := :circle
         fi.pearson_residuals
     end
     @series begin

test/test_visualization.jl

@@ -1,4 +1,4 @@
-@testset "Visualization Residuals" begin
+@testset "Visualizations" begin
     ω = [0.1, 0.1]
     A = [0.5 0; 0 0.5]
     B = [0.5 0; 0 0.5]
@@ -8,4 +8,12 @@
     f = fit!(gas, simulation; verbose = 2, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
     rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), f)
     @test length(rec) == 7
+    ω = [0.1, 0.1]
+    A = [0.5 0; 0 0.5]
+    B = [0.5 0; 0 0.5]
+    simulation = simulate_GAS_1_1(Normal, 0.0, ω, A, B, 1)
+    gas = ScoreDrivenModels.Model(1, 1, Normal, 0.0)
+    bac = backtest(gas, simulation, 10, 4985)
+    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), bac, "name")
+    @test length(rec) == 2
 end