Add visualization and fix backtests (#129)

Author: guilhermebodin

Project.toml

@@ -9,14 +9,17 @@ HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
 Distributions = "0.23"
+HypothesisTests = "0.10"
 Optim = "0.20, 0.21, 0.22, 1.2"
+RecipesBase = "1"
 SpecialFunctions = "0.8"
 julia = "1"
-HypothesisTests = "0.9"
 
 [extras]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

src/ScoreDrivenModels.jl

@@ -1,6 +1,6 @@
 module ScoreDrivenModels
 
-using Distributions, Optim, SpecialFunctions, HypothesisTests
+using Distributions, Optim, SpecialFunctions, HypothesisTests, RecipesBase, StatsBase
 using LinearAlgebra, Printf
 
 import Base: length, deepcopy, show
@@ -45,4 +45,7 @@ include("distributions/tdist.jl")
 include("distributions/tdistlocationscale.jl")
 include("distributions/weibull.jl")
 
+include("visualization/forecast.jl")
+include("visualization/residuals.jl")
+
 end

src/backtest.jl

@@ -35,11 +35,11 @@ end
 TODO
 """
 function backtest(gas::Model{<:Distribution, T}, y::Vector{T}, steps_ahead::Int, start_idx::Int;
-                  S::Int = 1000,
-                  initial_params::Matrix{T} = stationary_initial_params(gas),
-                  opt_method = NelderMead(gas, DEFAULT_NUM_SEEDS)) where T
+                  S::Int = 10_000,
+                  initial_params = stationary_initial_params(gas),
+                  opt_method = NelderMead(gas, 3)) where T
     num_mle = length(y) - start_idx - steps_ahead
-    backtest = Backtest(num_mle, steps_ahead)
+    b = Backtest(num_mle, steps_ahead)
     for i in 1:num_mle
         println("Backtest: step $i of $num_mle")
         gas_to_fit = deepcopy(gas)
@@ -49,8 +49,8 @@ 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)
-        backtest.abs_errors[i, :] = abs_errors
-        backtest.crps_scores[i, :] = crps_scores
+        b.abs_errors[i, :] = abs_errors
+        b.crps_scores[i, :] = crps_scores
     end
-    return backtest
+    return b
 end

src/simulate.jl

@@ -36,10 +36,7 @@ function simulate(series::Vector{T}, gas::Model{D, T}, H::Int, S::Int;
         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, :]
+        parameter_scenarios[:, :, scenario] = param[biggest_lag+1:end, :]
     end
     return observation_scenarios, parameter_scenarios
 end
@@ -69,12 +66,18 @@ function forecast(series::Vector{T}, gas::Model{D, T}, H::Int;
 
     observation_scenarios, parameter_scenarios = simulate(series, gas, H, S; 
                                                     initial_params = initial_params)
+    parameters_forecast = mean(parameter_scenarios, dims = 3)[:, :, 1]
+    observations_forecast = Vector{T}(undef, H)
+    for i in 1:H
+        dist = update_dist(D, parameters_forecast, i)
+        observations_forecast[i] = mean(dist)
+    end
     return Forecast(
         get_quantiles(quantiles, observation_scenarios), 
-        mean(observation_scenarios, dims = 2)[:],
+        observations_forecast,
         observation_scenarios,
         get_quantiles(quantiles, parameter_scenarios), 
-        mean(parameter_scenarios, dims = 3)[:, :, 1],
+        parameters_forecast,
         parameter_scenarios
         )
 end

src/visualization/forecast.jl

@@ -0,0 +1 @@
+#TODO

src/visualization/residuals.jl

@@ -0,0 +1,58 @@
+RecipesBase.@recipe function f(fi::ScoreDrivenModels.Fitted)
+    layout := (2, 2)
+    acf = autocor(fi.pearson_residuals)[2:end]
+    @series begin
+        seriestype := :path
+        label := ""
+        seriescolor := "black"
+        subplot := 1
+        fi.pearson_residuals
+    end
+    @series begin
+        seriestype := :bar
+        label := ""
+        seriescolor := "black"
+        subplot := 2
+        acf
+    end
+    ub = ones(length(acf)) * 1.96 / sqrt(fi.num_obs)
+    lb = ones(length(acf)) * -1.96 / sqrt(fi.num_obs)
+    @series begin
+        seriestype := :path
+        linestyle := :dash
+        seriescolor := "red"
+        label := ""
+        subplot := 2
+        ub
+    end
+    @series begin
+        seriestype := :path
+        linestyle := :dash
+        label := ""
+        seriescolor := "red"
+        subplot := 2
+        lb
+    end
+    @series begin
+        seriestype := :histogram
+        label := ""
+        seriescolor := "black"
+        subplot := 3
+        fi.pearson_residuals
+    end
+    qqpair = qqbuild(Normal(), fi.pearson_residuals)
+    @series begin
+        seriestype := :scatter
+        label := ""
+        seriescolor := "black"
+        subplot := 4
+        qqpair.qx, qqpair.qy
+    end
+    @series begin
+        seriestype := :path
+        seriescolor := "red"
+        label := ""
+        subplot := 4
+        collect(-3:0.01:3), collect(-3:0.01:3)
+    end
+end

test/runtests.jl

@@ -1,4 +1,4 @@
-using ScoreDrivenModels, Distributions, LinearAlgebra
+using ScoreDrivenModels, Distributions, LinearAlgebra, RecipesBase
 
 using Test, Random, HypothesisTests, DelimitedFiles
 
@@ -12,3 +12,4 @@ include("test_diagnostics.jl")
 include("test_estimate.jl")
 include("test_simulate.jl")
 include("test_backtest.jl")
+include("test_visualization.jl")

test/test_visualization.jl

@@ -0,0 +1,11 @@
+@testset "Visualization Residuals" begin
+    ω = [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)
+    # LBFGS
+    gas = Model(1, 1, Normal, 0.0)
+    f = fit!(gas, simulation; verbose = 2, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
+    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), f)
+    @test length(rec) == 7
+end