Merge pull request #19 from LAMPSPUC/change_predict_simulate

separate predict and simulate functions
LAMPSPUC · Jul 29, 2024 · 9301d9d · 9301d9d
2 parents 51798c0 + 593e423
commit 9301d9d
Show file tree

Hide file tree

Showing 6 changed files with 1,369 additions and 1,323 deletions.
diff --git a/examples/2_forecast_example.ipynb b/examples/2_forecast_example.ipynb
diff --git a/examples/3_components_example.ipynb b/examples/3_components_example.ipynb
@@ -1,5 +1,14 @@
 {
  "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# UnobservedComponents Example\n",
+    "\n",
+    "In this notebook, we have developed code to show an example of how to obtain and plot the estimated unobserved components of a UnobservedComponentsGAS model."
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},

diff --git a/src/forecast.jl b/src/forecast.jl
@@ -518,6 +518,75 @@ function simulate(gas_model::GASModel, output::Output, dict_hyperparams_and_fitt
     return pred_y
 end
 
+"""
+# simulate(gas_model::GASModel, output::Output, y::Vector{Float64}, steps_ahead::Int64, num_scenarios::Int64
+
+Simulates future values of a time series using the specified GAS model.
+
+## Arguments
+- `gas_model::GASModel`: The GAS model containing parameters and specifications.
+- `output::Output`: The output structure containing the fitted values, residuals, and information criteria.
+- `y::Vector{Float64}`: The vector of observed values for the time series data.
+- `steps_ahead::Int64`: An integer indicating the number of steps ahead to predict.
+- `num_scenarios::Int64`: An integer indicating the number of scenarios to simulate.
+
+## Returns
+- simulated_scenarios: A matrix containing the simulated scenarios, where each column represents a scenario.
+"""
+function simulate(gas_model::GASModel, output::Output, y::Vector{Float64}, steps_ahead::Int64, num_scenarios::Int64)
+
+    new_output = deepcopy(output)
+    if typeof(gas_model.dist) == LogNormalDistribution
+        new_output.fitted_params = convert_to_log_scale(new_output.fitted_params)
+        y = log.(y)
+    end
+
+    dict_hyperparams_and_fitted_components = get_dict_hyperparams_and_fitted_components_with_forecast(gas_model, new_output, steps_ahead, num_scenarios)
+    simulated_scenarios                    = simulate(gas_model, new_output, dict_hyperparams_and_fitted_components, y, steps_ahead, num_scenarios)
+
+    if typeof(gas_model.dist) == LogNormalDistribution
+        simulated_scenarios = convert_forecast_scenarios_to_exp_scale(simulated_scenarios)
+    end
+
+   return simulated_scenarios[end-steps_ahead+1:end, :]
+end
+
+"""
+# simulate(gas_model::GASModel, output::Output, y::Vector{Float64}, X_forecast::Matrix{Fl}, steps_ahead::Int64, num_scenarios::Int64
+
+Simulates future values of a time series using the specified GAS model and explanatory variables for forecasting.
+
+## Arguments
+- `gas_model::GASModel`: The GAS model containing parameters and specifications.
+- `output::Output`: The output structure containing the fitted values, residuals, and information criteria.
+- `y::Vector{Float64}`: The vector of observed values for the time series data.
+- `X_forecast::Matrix{Fl}`: The matrix of explanatory variables for forecasting.
+- `steps_ahead::Int64`: An integer indicating the number of steps ahead to predict.
+- `num_scenarios::Int64`: An integer indicating the number of scenarios to simulate.
+
+## Returns
+- `scenarios`: A matrix containing the simulated scenarios, where each column represents a scenario.
+"""
+function simulate(gas_model::GASModel, output::Output, y::Vector{Float64}, X_forecast::Matrix{Fl}, steps_ahead::Int64, num_scenarios::Int64) where {Ml, Fl}
+
+    new_output = deepcopy(output)
+
+    if typeof(gas_model.dist) == LogNormalDistribution
+        new_output.fitted_params = convert_to_log_scale(new_output.fitted_params)
+        y = log.(y)
+        #X_forecast = log.(X_forecast) #Devo fazer isso?
+    end
+
+    dict_hyperparams_and_fitted_components = get_dict_hyperparams_and_fitted_components_with_forecast(gas_model, new_output, X_forecast, steps_ahead, num_scenarios)
+    simulated_scenarios                    = simulate(gas_model, new_output, dict_hyperparams_and_fitted_components, y, X_forecast, steps_ahead, num_scenarios)
+
+    if typeof(gas_model.dist) == LogNormalDistribution
+        simulated_scenarios = convert_forecast_scenarios_to_exp_scale(simulated_scenarios)
+    end
+
+    return simulated_scenarios[end-steps_ahead+1:end, :]
+end
+
 """
 # get_mean_and_intervals_prediction(pred_y::Matrix{Fl}, steps_ahead::Int64, probabilistic_intervals::Vector{Float64}) where Fl
 
@@ -562,7 +631,7 @@ function get_mean_and_intervals_prediction(pred_y::Matrix{Fl}, steps_ahead::Int6
     end
 
     dict_forec["mean"] = forec
-    dict_forec["scenarios"] = pred_y[end - steps_ahead + 1:end, :]
+    # dict_forec["scenarios"] = pred_y[end - steps_ahead + 1:end, :]
 
     return dict_forec
 end
@@ -590,20 +659,9 @@ Predicts future values of a time series using the specified GAS model.
 """
 function predict(gas_model::GASModel, output::Output, y::Vector{Float64}, steps_ahead::Int64, num_scenarios::Int64; probabilistic_intervals::Vector{Float64} = [0.8, 0.95])
 
-    new_output = deepcopy(output)
-    if typeof(gas_model.dist) == LogNormalDistribution
-        new_output.fitted_params = convert_to_log_scale(new_output.fitted_params)
-        y = log.(y)
-    end
-
-    dict_hyperparams_and_fitted_components = get_dict_hyperparams_and_fitted_components_with_forecast(gas_model, new_output, steps_ahead, num_scenarios)
-    pred_y                                 = simulate(gas_model, new_output, dict_hyperparams_and_fitted_components, y, steps_ahead, num_scenarios)
-
-    if typeof(gas_model.dist) == LogNormalDistribution
-        pred_y = convert_forecast_scenarios_to_exp_scale(pred_y)
-    end
+    simulated_scenarios = simulate(gas_model, output, y, steps_ahead, num_scenarios)
 
-    dict_forec = get_mean_and_intervals_prediction(pred_y, steps_ahead, probabilistic_intervals)
+    dict_forec = get_mean_and_intervals_prediction(simulated_scenarios, steps_ahead, probabilistic_intervals)
 
     return dict_forec
 end
@@ -632,22 +690,9 @@ Predicts future values of a time series using the specified GAS model and explan
 """
 function predict(gas_model::GASModel, output::Output, y::Vector{Float64}, X_forecast::Matrix{Fl}, steps_ahead::Int64, num_scenarios::Int64; probabilistic_intervals::Vector{Float64} = [0.8, 0.95]) where {Ml, Fl}
 
-    new_output = deepcopy(output)
-
-    if typeof(gas_model.dist) == LogNormalDistribution
-        new_output.fitted_params = convert_to_log_scale(new_output.fitted_params)
-        y = log.(y)
-        #X_forecast = log.(X_forecast) #Devo fazer isso?
-    end
-
-    dict_hyperparams_and_fitted_components = get_dict_hyperparams_and_fitted_components_with_forecast(gas_model, new_output, X_forecast, steps_ahead, num_scenarios)
-    pred_y                                 = simulate(gas_model, new_output, dict_hyperparams_and_fitted_components, y, X_forecast, steps_ahead, num_scenarios)
-
-    if typeof(gas_model.dist) == LogNormalDistribution
-        pred_y = convert_forecast_scenarios_to_exp_scale(pred_y)
-    end
+    simulated_scenarios = simulate(gas_model, output, y, X_forecast, steps_ahead, num_scenarios)
 
-    dict_forec = get_mean_and_intervals_prediction(pred_y, steps_ahead, probabilistic_intervals)
+    dict_forec = get_mean_and_intervals_prediction(simulated_scenarios, steps_ahead, probabilistic_intervals)
 
     return dict_forec
 end
diff --git a/test/test_fit_forecast_lognormal.jl b/test/test_fit_forecast_lognormal.jl
@@ -65,6 +65,10 @@
     #forecast_lognormal_X       = UnobservedComponentsGAS.predict(gas_model_lognormal_X, fitted_model_lognormal_X, y, X_lognormal_forec, steps_ahead, num_scenarious)
     forecast_lognormal_2params = UnobservedComponentsGAS.predict(gas_model_lognormal_2params, fitted_model_lognormal_2params, y, steps_ahead, num_scenarious)
 
+    scenarios_lognormal         = UnobservedComponentsGAS.simulate(gas_model_lognormal, fitted_model_lognormal, y, steps_ahead, num_scenarious)
+    # scenarios_lognormal_X       = UnobservedComponentsGAS.simulate(gas_model_lognormal_X, fitted_model_lognormal_X, y, X_lognormal_forec, steps_ahead, num_scenarious)
+    scenarios_lognormal_2params = UnobservedComponentsGAS.simulate(gas_model_lognormal_2params, fitted_model_lognormal_2params, y, steps_ahead, num_scenarious)
+
     @testset " --- Testing create_model functions" begin
         @test(size(parameters_lognormal)         == (T,2))
         @test(size(parameters_lognormal_2params) == (T,2))
@@ -111,21 +115,21 @@
     end
 
     @testset " --- Test forecast function ---" begin
-        @test(isapprox(forecast_lognormal["mean"], vec(mean(forecast_lognormal["scenarios"], dims = 2)); rtol = 1e-3)) 
-        @test(size(forecast_lognormal["scenarios"]) == (steps_ahead, num_scenarious))
+        @test(isapprox(forecast_lognormal["mean"], vec(mean(scenarios_lognormal, dims = 2)); rtol = 1e-3)) 
+        @test(size(scenarios_lognormal) == (steps_ahead, num_scenarious))
 
-        @test(isapprox(forecast_lognormal["intervals"]["80"]["lower"], [quantile(forecast_lognormal["scenarios"][t,:], 0.2/2) for t in 1:steps_ahead]))
-        @test(isapprox(forecast_lognormal["intervals"]["80"]["upper"], [quantile(forecast_lognormal["scenarios"][t,:], 1 - 0.2/2) for t in 1:steps_ahead]))
-        @test(isapprox(forecast_lognormal["intervals"]["95"]["lower"], [quantile(forecast_lognormal["scenarios"][t,:], 0.05/2) for t in 1:steps_ahead]))
-        @test(isapprox(forecast_lognormal["intervals"]["95"]["upper"], [quantile(forecast_lognormal["scenarios"][t,:], 1 - 0.05/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_lognormal["intervals"]["80"]["lower"], [quantile(scenarios_lognormal[t,:], 0.2/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_lognormal["intervals"]["80"]["upper"], [quantile(scenarios_lognormal[t,:], 1 - 0.2/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_lognormal["intervals"]["95"]["lower"], [quantile(scenarios_lognormal[t,:], 0.05/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_lognormal["intervals"]["95"]["upper"], [quantile(scenarios_lognormal[t,:], 1 - 0.05/2) for t in 1:steps_ahead]))
 
-        @test(isapprox(forecast_lognormal_2params["mean"], vec(mean(forecast_lognormal_2params["scenarios"], dims = 2)); rtol = 1e-3)) 
-        @test(size(forecast_lognormal_2params["scenarios"]) == (steps_ahead, num_scenarious))
+        @test(isapprox(forecast_lognormal_2params["mean"], vec(mean(scenarios_lognormal_2params, dims = 2)); rtol = 1e-3)) 
+        @test(size(scenarios_lognormal_2params) == (steps_ahead, num_scenarious))
 
-        @test(isapprox(forecast_lognormal_2params["intervals"]["80"]["lower"], [quantile(forecast_lognormal_2params["scenarios"][t,:], 0.2/2) for t in 1:steps_ahead]))
-        @test(isapprox(forecast_lognormal_2params["intervals"]["80"]["upper"], [quantile(forecast_lognormal_2params["scenarios"][t,:], 1 - 0.2/2) for t in 1:steps_ahead]))
-        @test(isapprox(forecast_lognormal_2params["intervals"]["95"]["lower"], [quantile(forecast_lognormal_2params["scenarios"][t,:], 0.05/2) for t in 1:steps_ahead]))
-        @test(isapprox(forecast_lognormal_2params["intervals"]["95"]["upper"], [quantile(forecast_lognormal_2params["scenarios"][t,:], 1 - 0.05/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_lognormal_2params["intervals"]["80"]["lower"], [quantile(scenarios_lognormal_2params[t,:], 0.2/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_lognormal_2params["intervals"]["80"]["upper"], [quantile(scenarios_lognormal_2params[t,:], 1 - 0.2/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_lognormal_2params["intervals"]["95"]["lower"], [quantile(scenarios_lognormal_2params[t,:], 0.05/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_lognormal_2params["intervals"]["95"]["upper"], [quantile(scenarios_lognormal_2params[t,:], 1 - 0.05/2) for t in 1:steps_ahead]))
 
         # @test(isapprox(forecast_lognormal_X["mean"], vec(mean(forecast_lognormal_X["scenarios"], dims = 2)); rtol = 1e-3)) 
         # @test(size(forecast_lognormal_X["scenarios"]) == (steps_ahead, num_scenarious))

diff --git a/test/test_fit_forecast_normal.jl b/test/test_fit_forecast_normal.jl
@@ -64,6 +64,10 @@
     forecast_normal         = UnobservedComponentsGAS.predict(gas_model_normal, fitted_model_normal, y, steps_ahead, num_scenarious)
     forecast_normal_X       = UnobservedComponentsGAS.predict(gas_model_normal_X, fitted_model_normal_X, y, X_normal_forec, steps_ahead, num_scenarious)
     forecast_normal_2params = UnobservedComponentsGAS.predict(gas_model_normal_2params, fitted_model_normal_2params, y, steps_ahead, num_scenarious)
+
+    scenarios_normal         = UnobservedComponentsGAS.simulate(gas_model_normal, fitted_model_normal, y, steps_ahead, num_scenarious)
+    scenarios_normal_X       = UnobservedComponentsGAS.simulate(gas_model_normal_X, fitted_model_normal_X, y, X_normal_forec, steps_ahead, num_scenarious)
+    scenarios_normal_2params = UnobservedComponentsGAS.simulate(gas_model_normal_2params, fitted_model_normal_2params, y, steps_ahead, num_scenarious)
 
     @testset "create_model" begin
     # Create model with no explanatory series
@@ -116,29 +120,29 @@
     #@info(" --- Test forecast function ---")
     @testset "forecast" begin
 
-        @test(isapprox(forecast_normal["mean"], vec(mean(forecast_normal["scenarios"], dims = 2)); rtol = 1e-3)) 
-        @test(size(forecast_normal["scenarios"]) == (steps_ahead, num_scenarious))
+        @test(isapprox(forecast_normal["mean"], vec(mean(scenarios_normal, dims = 2)); rtol = 1e-3)) 
+        @test(size(scenarios_normal) == (steps_ahead, num_scenarious))
 
-        @test(isapprox(forecast_normal["intervals"]["80"]["lower"], [quantile(forecast_normal["scenarios"][t,:], 0.2/2) for t in 1:steps_ahead]))
-        @test(isapprox(forecast_normal["intervals"]["80"]["upper"], [quantile(forecast_normal["scenarios"][t,:], 1 - 0.2/2) for t in 1:steps_ahead]))
-        @test(isapprox(forecast_normal["intervals"]["95"]["lower"], [quantile(forecast_normal["scenarios"][t,:], 0.05/2) for t in 1:steps_ahead]))
-        @test(isapprox(forecast_normal["intervals"]["95"]["upper"], [quantile(forecast_normal["scenarios"][t,:], 1 - 0.05/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_normal["intervals"]["80"]["lower"], [quantile(scenarios_normal[t,:], 0.2/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_normal["intervals"]["80"]["upper"], [quantile(scenarios_normal[t,:], 1 - 0.2/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_normal["intervals"]["95"]["lower"], [quantile(scenarios_normal[t,:], 0.05/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_normal["intervals"]["95"]["upper"], [quantile(scenarios_normal[t,:], 1 - 0.05/2) for t in 1:steps_ahead]))
 
-        @test(isapprox(forecast_normal_2params["mean"], vec(mean(forecast_normal_2params["scenarios"], dims = 2)); rtol = 1e-3)) 
-        @test(size(forecast_normal_2params["scenarios"]) == (steps_ahead, num_scenarious))
+        @test(isapprox(forecast_normal_2params["mean"], vec(mean(scenarios_normal_2params, dims = 2)); rtol = 1e-3)) 
+        @test(size(scenarios_normal_2params) == (steps_ahead, num_scenarious))
 
-        @test(isapprox(forecast_normal_2params["intervals"]["80"]["lower"], [quantile(forecast_normal_2params["scenarios"][t,:], 0.2/2) for t in 1:steps_ahead]))
-        @test(isapprox(forecast_normal_2params["intervals"]["80"]["upper"], [quantile(forecast_normal_2params["scenarios"][t,:], 1 - 0.2/2) for t in 1:steps_ahead]))
-        @test(isapprox(forecast_normal_2params["intervals"]["95"]["lower"], [quantile(forecast_normal_2params["scenarios"][t,:], 0.05/2) for t in 1:steps_ahead]))
-        @test(isapprox(forecast_normal_2params["intervals"]["95"]["upper"], [quantile(forecast_normal_2params["scenarios"][t,:], 1 - 0.05/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_normal_2params["intervals"]["80"]["lower"], [quantile(scenarios_normal_2params[t,:], 0.2/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_normal_2params["intervals"]["80"]["upper"], [quantile(scenarios_normal_2params[t,:], 1 - 0.2/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_normal_2params["intervals"]["95"]["lower"], [quantile(scenarios_normal_2params[t,:], 0.05/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_normal_2params["intervals"]["95"]["upper"], [quantile(scenarios_normal_2params[t,:], 1 - 0.05/2) for t in 1:steps_ahead]))
 
-        @test(isapprox(forecast_normal_X["mean"], vec(mean(forecast_normal_X["scenarios"], dims = 2)); rtol = 1e-3)) 
-        @test(size(forecast_normal_X["scenarios"]) == (steps_ahead, num_scenarious))
+        @test(isapprox(forecast_normal_X["mean"], vec(mean(scenarios_normal_X, dims = 2)); rtol = 1e-3)) 
+        @test(size(scenarios_normal_X) == (steps_ahead, num_scenarious))
 
-        @test(isapprox(forecast_normal_X["intervals"]["80"]["lower"], [quantile(forecast_normal_X["scenarios"][t,:], 0.2/2) for t in 1:steps_ahead]))
-        @test(isapprox(forecast_normal_X["intervals"]["80"]["upper"], [quantile(forecast_normal_X["scenarios"][t,:], 1 - 0.2/2) for t in 1:steps_ahead]))
-        @test(isapprox(forecast_normal_X["intervals"]["95"]["lower"], [quantile(forecast_normal_X["scenarios"][t,:], 0.05/2) for t in 1:steps_ahead]))
-        @test(isapprox(forecast_normal_X["intervals"]["95"]["upper"], [quantile(forecast_normal_X["scenarios"][t,:], 1 - 0.05/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_normal_X["intervals"]["80"]["lower"], [quantile(scenarios_normal_X[t,:], 0.2/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_normal_X["intervals"]["80"]["upper"], [quantile(scenarios_normal_X[t,:], 1 - 0.2/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_normal_X["intervals"]["95"]["lower"], [quantile(scenarios_normal_X[t,:], 0.05/2) for t in 1:steps_ahead]))
+        @test(isapprox(forecast_normal_X["intervals"]["95"]["upper"], [quantile(scenarios_normal_X[t,:], 1 - 0.05/2) for t in 1:steps_ahead]))
     end
 
     #@info(" --- Test quality of fit and forecast - Normal")