Forecast with confidence intervals (#62)

* forecast working

* fix 0.95 ci

* take out push loadpath

Author: guilhermebodin

examples/inflow_lognormal.jl

@@ -46,19 +46,27 @@ inflow = [
 
 # Convert data to vector
 y = Vector{Float64}(vec(inflow'))
+y_train = y[1:400]
+y_test = y[401:end]
 
 # Specify GAS model: here we use lag 1 for trend characterization and lag 12 for seasonality characterization
-gas = GAS.Model([1, 12], [1, 12], LogNormal, 0.0)
+gas = GAS.Model([1, 2, 11, 12], [1, 2, 11, 12], LogNormal, 0.0; time_varying_params = [1])
 
 # Define initial_params with
-initial_params = dynamic_initial_params(y, gas)
+initial_params = dynamic_initial_params(y_train, gas)
 
 # Estimate the model via MLE
-fit!(gas, y; initial_params = initial_params)
+fit!(gas, y_train; initial_params = initial_params, opt_method = NelderMead(gas, 100))
 
 # Obtain in-sample estimates for the inflow
-y_gas = fitted_mean(gas, y; initial_params = initial_params)
+y_fitted = fitted_mean(gas, y_train; initial_params = initial_params)
 
 # Compare observations and in-sample estimates
-plot(y, label = "historical inflow")
-plot!(y_gas, label = "in-sample estimates")
+plot(y_train, label = "In-sample inflow")
+plot!(y_fitted, label = "in-sample estimates")
+
+# Forecasts with 95% confidence interval
+forec = GAS.forecast(y_train, gas, 80; initial_params = initial_params, ci = [0.95])
+
+plot(y_test, label = "Out-of-sample inflow")
+plot!(forec, label = "Forecast", color = "Steel Blue")

src/simulate.jl

@@ -17,8 +17,8 @@ function simulate(serie::Vector{T}, gas::Model{D, T}, N::Int, S::Int;
     # Create scenarios matrix
     scenarios = Matrix{T}(undef, N, S)
     for scenario in 1:S
-        sim, param = simulate_recursion(gas, N + biggest_lag; initial_params = params_simulation)
-        scenarios[:, scenario] = sim[biggest_lag + 1:end]
+        sim, param = simulate_recursion(gas, N + biggest_lag + 1; initial_params = params_simulation)
+        scenarios[:, scenario] = sim[biggest_lag + 2:end]
     end
 
     return scenarios
@@ -29,20 +29,45 @@ end
     forecast(serie::Vector{T}, gas::Model{D, T}, N::Int; kwargs...) where {D, T}
 
 Forecast future values of a time series by updating the GAS recursion `N` times and 
-taking the mean of the distribution at each time.
+taking the mean of the distribution at each time. You can pass the desired confidence interval 
+as a `Vector{T}`. The forecast will be the first column and the confidence intervals are the remaining
+columns.
+
+The forecast is build using Monte Carlo method as in Blasques, Francisco, Siem Jan Koopman, 
+Katarzyna Lasak and Andre Lucas (2016): "In-Sample Confidence Bounds and Out-of-Sample Forecast 
+Bands for Time-Varying Parameters in Observation Driven Models", International Journal of Forecasting, 32(3), 875-887. 
+
+By default 1000 scenarios are used but one can change by switching the `S` keyword argument.
 """
 function forecast(serie::Vector{T}, gas::Model{D, T}, N::Int;
-                    initial_params::Matrix{T} = stationary_initial_params(gas)) where {D, T}
-    # Filter params estimated on the time series
-    params = score_driven_recursion(gas, serie; initial_params = initial_params)
+                    initial_params::Matrix{T} = stationary_initial_params(gas),
+                    ci::Vector{T} = T.([0.95]), S::Int = 1000) where {D, T}
 
-    biggest_lag = number_of_lags(gas)
-
-    params_simulation = params[(end - biggest_lag):(end - 1), :]
-    # Create scenarios matrix
-    forec = Vector{T}(undef, N)
-    sim, param = simulate_recursion(gas, N + biggest_lag; initial_params = params_simulation, update = mean)
-    forec = sim[biggest_lag + 1:end]
+    scenarios = simulate(serie, gas, N, S; initial_params = initial_params)
+    
+    quantiles = get_quantiles(ci)
 
+    forec = Matrix{T}(undef, N, length(quantiles) + 1)
+    
+    forec[:, 1] = mean(scenarios, dims = 2)
+    for t in 1:N
+        for (i, q) in enumerate(quantiles)
+            forec[t, i + 1] = quantile(scenarios[t, :], q)
+        end
+    end
+    
     return forec
+end
+
+function get_quantiles(ci::Vector{T}) where T
+    @assert all((ci .< 1.0) .& (ci .> 0.0))
+    quantiles = zeros(2 * length(ci))
+    i = 1
+    for v in ci
+        quantiles[i] = 0.5 + v/2
+        i += 1
+        quantiles[i] = 0.5 - v/2
+        i += 1
+    end
+    return quantiles
 end

src/univariate_score_driven_recursion.jl

@@ -134,6 +134,10 @@ return the fitted mean.. #TODO
 function fitted_mean(gas::Model{D, T}, observations::Vector{T}; 
                      initial_params::Matrix{T} = stationary_initial_params(gas)) where {D, T}
     params_fitted = score_driven_recursion(gas, observations; initial_params = initial_params)
+
+    # Discard last step of the update
+    params_fitted = params_fitted[1:end-1, :]
+
     n = size(params_fitted, 1)
     fitted_mean = Vector{T}(undef, n)