Merge pull request #272 from LAMPSPUC/cross_validation

Change from backtest to cross validation

Author: guilhermebodin

docs/src/examples.md

@@ -193,10 +193,10 @@ smoother_output = kalman_smoother(model)
 plot(df.date, get_smoothed_state(smoother_output)[:, 2], label = "slope")
 ```
 
-## Backtest the forecasts of a model
+## Cross validation of the forecasts of a model
 
 Often times users would like to compare the forecasting skill of different models. The function 
-[`backtest`](@ref) makes it easy to make a rolling window scheme of estimations and forecasts 
+[`cross_validation`](@ref) makes it easy to make a rolling window scheme of estimations and forecasts 
 that allow users to track each model forecasting skill per lead time. A simple plot recipe is
 implemented to help users to interpret the results easily.
 
@@ -210,6 +210,6 @@ using CSV, DataFrames
 air_passengers = CSV.File(StateSpaceModels.AIR_PASSENGERS) |> DataFrame
 log_air_passengers = log.(air_passengers.passengers)
 model = BasicStructural(log_air_passengers, 12)
-b = backtest(model, 24, 50)
+b = cross_validation(model, 24, 50)
 plot(b, "Basic structural model")
 ```

docs/src/manual.md

@@ -154,12 +154,12 @@ isfitted
 
 StateSpaceModels.jl has functions to make forecasts of the predictive densities multiple steps ahead and to
 simulate scenarios based on those forecasts. The package also has a functions to benchmark the model forecasts 
-using backtest techniques.
+using cross_validation techniques.
 
 ```@docs
 forecast
 simulate_scenarios
-backtest
+cross_validation
 ```
 
 ## Visualization

src/StateSpaceModels.jl

@@ -33,7 +33,7 @@ include("optimizers.jl")
 include("fit.jl")
 include("prints.jl")
 include("forecast.jl")
-include("backtest.jl")
+include("cross_validation.jl")
 
 include("models/common.jl")
 include("models/locallevel.jl")
@@ -52,7 +52,7 @@ include("models/dar.jl")
 
 include("visualization/forecast.jl")
 include("visualization/components.jl")
-include("visualization/backtest.jl")
+include("visualization/cross_validation.jl")
 include("visualization/diagnostics.jl")
 
 # Exported types and structs
@@ -82,7 +82,7 @@ export UnobservedComponents
 
 # Exported functions
 export auto_ets
-export backtest
+export cross_validation
 export constrain_box!
 export constrain_identity!
 export constrain_variance!

src/backtest.jl renamed to src/cross_validation.jl

@@ -1,9 +1,9 @@
-struct Backtest{Fl <: AbstractFloat}
+struct CrossValidation{Fl <: AbstractFloat}
     abs_errors::Matrix{Fl}
     mae::Vector{Fl}
     crps_scores::Matrix{Fl}
     mean_crps::Vector{Fl}
-    function Backtest{Fl}(n::Int, steps_ahead::Int) where Fl
+    function CrossValidation{Fl}(n::Int, steps_ahead::Int) where Fl
         abs_errors = Matrix{Fl}(undef, steps_ahead, n)
         crps_scores = Matrix{Fl}(undef, steps_ahead, n)
         mae = Vector{Fl}(undef, steps_ahead)
@@ -34,27 +34,27 @@ function evaluate_crps(y::Vector{Fl}, scenarios::Matrix{Fl}) where {Fl}
 end
 
 """
-    backtest(model::StateSpaceModel, steps_ahead::Int, start_idx::Int;
+    cross_validation(model::StateSpaceModel, steps_ahead::Int, start_idx::Int;
              n_scenarios::Int = 10_000,
              filter::KalmanFilter=default_filter(model),
              optimizer::Optimizer=default_optimizer(model)) where Fl
 
 Makes rolling window estimating and forecasting to benchmark the forecasting skill of the model
 in for different time periods and different lead times. The function returns a struct with the MAE
-and mean CRPS per lead time. See more on [Backtest the forecasts of a model](@ref)
+and mean CRPS per lead time. See more on [CrossValidation the forecasts of a model](@ref)
 
 # References
  * DTU course "31761 - Renewables in electricity markets" available on youtube https://www.youtube.com/watch?v=Ffo8XilZAZw&t=556s
 """
-function backtest(model::StateSpaceModel, steps_ahead::Int, start_idx::Int;
+function cross_validation(model::StateSpaceModel, steps_ahead::Int, start_idx::Int;
                   n_scenarios::Int = 10_000,
                   filter::KalmanFilter=default_filter(model),
                   optimizer::Optimizer=default_optimizer(model))
     Fl = typeof_model_elements(model)
     num_mle = length(model.system.y) - start_idx - steps_ahead
-    b = Backtest{Fl}(num_mle, steps_ahead)
+    cv = CrossValidation{Fl}(num_mle, steps_ahead)
     for i in 1:num_mle
-        println("Backtest: step $i of $num_mle")
+        println("CrossValidation: step $i of $num_mle")
         y_to_fit = model.system.y[1:start_idx - 1 + i]
         y_to_verify = model.system.y[start_idx + i:start_idx - 1 + i + steps_ahead]
         model_to_fit = reinstantiate(model, y_to_fit)
@@ -64,12 +64,12 @@ function backtest(model::StateSpaceModel, steps_ahead::Int, start_idx::Int;
         expected_value_vector = forecast_expected_value(forec)[:]
         abs_errors = evaluate_abs_error(y_to_verify, expected_value_vector)
         crps_scores = evaluate_crps(y_to_verify, scenarios[:, 1, :])
-        b.abs_errors[:, i] = abs_errors
-        b.crps_scores[:, i] = crps_scores
+        cv.abs_errors[:, i] = abs_errors
+        cv.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, :])
+        cv.mae[i] = mean(cv.abs_errors[i, :])
+        cv.mean_crps[i] = mean(cv.crps_scores[i, :])
     end
-    return b
+    return cv
 end

src/models/naive_models.jl

@@ -12,7 +12,7 @@ typeof_model_elements(model::NaiveModel) = eltype(model.y)
 
 function assert_zero_missing_values(model::NaiveModel)
     for i in 1:length(model.y)
-        if isnan(y[i])
+        if isnan(model.y[i])
             return error("model $(typeof(model)) does not support missing values.)")
         end
     end
@@ -44,6 +44,7 @@ mutable struct Naive <: NaiveModel
 end
 
 function fit!(model::Naive)
+    assert_zero_missing_values(model)
     residuals = model.y[2:end] - model.y[1:end-1]
     model.residuals = residuals
     model.sigma2 = var(residuals)
@@ -111,6 +112,7 @@ mutable struct SeasonalNaive <: NaiveModel
 end
 
 function fit!(model::SeasonalNaive)
+    assert_zero_missing_values(model)
     residuals = model.y[model.seasonal+1:end] - model.y[1:end-model.seasonal]
     model.residuals = residuals
     model.sigma2 = var(residuals)
@@ -185,6 +187,7 @@ mutable struct ExperimentalSeasonalNaive <: NaiveModel
 end
 
 function fit!(model::ExperimentalSeasonalNaive)
+    assert_zero_missing_values(model)
     residuals = model.y[model.seasonal+1:end] - model.y[1:end-model.seasonal]
     model.residuals = residuals
     model.sigma2 = var(residuals)
@@ -226,13 +229,13 @@ function reinstantiate(model::ExperimentalSeasonalNaive, y::Vector{<:Real})
     return ExperimentalSeasonalNaive(y, model.seasonal; S = model.S)
 end
 
-function backtest(model::NaiveModel, steps_ahead::Int, start_idx::Int;
+function cross_validation(model::NaiveModel, steps_ahead::Int, start_idx::Int;
                   n_scenarios::Int = 10_000)
     Fl = typeof_model_elements(model)
     num_fits = length(model.y) - start_idx - steps_ahead
-    b = Backtest{Fl}(num_fits, steps_ahead)
+    b = CrossValidation{Fl}(num_fits, steps_ahead)
     for i in 1:num_fits
-        println("Backtest: step $i of $num_fits")
+        println("CrossValidation: step $i of $num_fits")
         y_to_fit = model.y[1:start_idx - 1 + i]
         y_to_verify = model.y[start_idx + i:start_idx - 1 + i + steps_ahead]
         model_to_fit = reinstantiate(model, y_to_fit)

src/visualization/backtest.jl renamed to src/visualization/cross_validation.jl

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

test/models/naive_models.jl

@@ -3,11 +3,13 @@
 
     model = Naive(nile.flow)
     fit!(model)
+    StateSpaceModels.get_standard_residuals(model)
     forec = forecast(model, 10; bootstrap = true)
     forec = forecast(model, 10)
     scenarios = simulate_scenarios(model, 10, 1_000)
     @test monotone_forecast_variance(forec)
     test_scenarios_adequacy_with_forecast(forec, scenarios; rtol = 1e-1)
+    cross_validation(model, 10, 70; n_scenarios=100)
 
     air_passengers = CSV.File(StateSpaceModels.AIR_PASSENGERS) |> DataFrame
     log_air_passengers = log.(air_passengers.passengers)
@@ -18,10 +20,19 @@
     forec = forecast(model, 60)
     scenarios = simulate_scenarios(model, 60, 1_000)
     @test monotone_forecast_variance(forec)
+    StateSpaceModels.reinstantiate(model, model.y)
 
     # Just see if it runs
     model = ExperimentalSeasonalNaive(log_air_passengers, 12)
     fit!(model)
     forec = forecast(model, 60)
     scenarios = simulate_scenarios(model, 60, 1_000)
+    StateSpaceModels.reinstantiate(model, model.y)
+
+    using StateSpaceModels
+    using Test
+    y = randn(100)
+    y[10] = NaN
+    model = Naive(y)
+    @test_throws ErrorException fit!(model)
 end

test/runtests.jl

@@ -31,5 +31,5 @@ include("models/dar.jl")
 # Visualization
 include("visualization/forecast.jl")
 include("visualization/components.jl")
-include("visualization/backtest.jl")
+include("visualization/cross_validation.jl")
 include("visualization/diagnostics.jl")

test/visualization/components.jl

@@ -20,4 +20,24 @@
     @test length(rec) == 3
     rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), model, ks)
     @test length(rec) == 3
+
+    model = ExponentialSmoothing(log_air_passengers; trend = true, seasonal = 12)
+    fit!(model)
+    kf = kalman_filter(model)
+    ks = kalman_smoother(model)
+
+    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), model, kf)
+    @test length(rec) == 4
+    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), model, ks)
+    @test length(rec) == 4
+
+    finland_fatalities = CSV.File(StateSpaceModels.VEHICLE_FATALITIES) |> DataFrame
+    log_finland_fatalities = log.(finland_fatalities.ff)
+    model = ExponentialSmoothing(log_finland_fatalities; trend = true)
+    fit!(model)
+
+    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), model, kf)
+    @test length(rec) == 3
+    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), model, ks)
+    @test length(rec) == 3
 end

test/visualization/backtest.jl renamed to test/visualization/cross_validation.jl

@@ -3,7 +3,7 @@
     log_air_passengers = log.(air_passengers.passengers)
     model = BasicStructural(log_air_passengers, 12)
     # forecasting
-    b = backtest(model, 24, 110)
+    b = cross_validation(model, 24, 110)
     rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), b, "str")
     @test length(rec) == 2
 end