model name change and add more examples

Author: andre_ramos

README.md

@@ -17,17 +17,15 @@ y = randn(100)
 output = StateSpaceLearning.fit_model(y)
 
 #Main output options 
-model_type = output.model_input # State Space Equivalent model utilized in the estimation (default = Basic Structural).
+model_input         = output.model_input # Model inputs that were utilized to build the regression matrix.
+Create_X            = output.Create_X # The function utilized to build the regression matrix.
 X                   = output.X # High Dimension Regression utilized in the estimation.
 coefs               = output.coefs # High Dimension Regression coefficients estimated in the estimation.
 ϵ                   = output.ϵ # Residuals of the model.
 fitted              = output.fitted # Fit in Sample of the model.
 components          = output.components # Dictionary containing information about each component of the model, each component has the keys: "Values" (The value of the component in each timestamp) , "Coefs" (The coefficients estimated for each element of the component) and "Indexes" (The indexes of the elements of the component in the high dimension regression "X").
 residuals_variances = output.residuals_variances # Dictionary containing the estimated variances for the innovations components (that is the information that can be utilized to initialize the state space model).
-T                   = output.T # The length of the original time series.
-outlier             = output.outlier # Boolean indicating the presence of outlier component (default = false).
 valid_indexes       = output.valid_indexes # Vector containing valid indexes of the time series (non valid indexes represent NaN values in the time series).
-ζ_ω_threshold         = output.ζ_ω_threshold # ζ_ω_threshold parameter (default = 0). A non 0 value for this parameter might be important in terms of forecast for some time series to lead to more stable predictions (we recommend ζ_ω_threshold = 11 for monthly series).
 
 #Forecast
 prediction = StateSpaceLearning.forecast(output, 12) #Gets a 12 steps ahead prediction
@@ -37,11 +35,9 @@ prediction = StateSpaceLearning.forecast(output, 12) #Gets a 12 steps ahead pred
 ## Fit Arguments
 
 * `y::Vector{Fl}`: Vector of data.
-* `model_input::Dict`: Dictionary containing the model input parameters (default: Dict("level" => true, "stochastic_level" => true, "trend" => true, "stochastic_trend" => true, "seasonal" => true, "stochastic_seasonal" => true, "freq_seasonal" => 12)).
+* `model_input::Dict`: Dictionary containing the model input parameters (default: Dict("level" => true, "stochastic_level" => true, "trend" => true, "stochastic_trend" => true, "seasonal" => true, "stochastic_seasonal" => true, "freq_seasonal" => 12, "outlier" => true, "ζ_ω_threshold" => 12)).
 * `estimation_input::Dict`: Dictionary containing the estimation input parameters (default: Dict("α" => 0.1, "information_criteria" => "aic", ψ => 0.05, "penalize_exogenous" => true, "penalize_initial_states" => true)).
 * `Exogenous_X::Union{Matrix{Fl}, Missing}`: Exogenous variables matrix (default: missing).
-* `outlier::Bool`: Flag for considering outlier component (default: true).
-* `ζ_ω_threshold::Int64`: ζ_ω_threshold parameter (default: 12).
 
 ## Features
 
@@ -51,9 +47,8 @@ Current features include:
 * Forecasting
 * Completion of missing values
 * Predefined models, including:
-  * Basic Structural"
-  * Local Linear Trend
-  * Local Level
+* Outlier detection
+* Outlier robust models
 
 ## Quick Examples
 
@@ -70,23 +65,49 @@ log_air_passengers = log.(airp.passengers)
 steps_ahead = 30
 
 output = StateSpaceLearning.fit_model(log_air_passengers)
-prediction_raw = StateSpaceLearning.forecast(output, steps_ahead)
-prediction = exp.(prediction_raw)
+prediction_log = StateSpaceLearning.forecast(output, steps_ahead)
+prediction = exp.(prediction_log)
 
 plot(airp.passengers, w=2 , color = "Black", lab = "Historical", legend = :outerbottom)
 plot!(vcat(ones(output.T).*NaN, prediction), lab = "Forcast", w=2, color = "blue")
 
 ```
 ![quick_example_airp](./docs/assets/quick_example_airp.PNG)
 
+### Component Extraction
+Quick example on how to perform component extraction in time series utilizing StateSpaceLearning.
+
+```julia
+using CSV
+using DataFrames
+using Plots
+
+airp = CSV.File(StateSpaceLearning.AIR_PASSENGERS) |> DataFrame
+log_air_passengers = log.(airp.passengers)
+
+output = StateSpaceLearning.fit_model(log_air_passengers)
+
+level = output.components["μ1"]["Values"] + output.components["ξ"]["Values"]
+slope = output.components["ν1"]["Values"] + output.components["ζ"]["Values"]
+seasonal = output.components["γ1"]["Values"] + output.components["ω"]["Values"]
+trend = level + slope
+
+plot(trend, w=2 , color = "Black", lab = "Trend Component", legend = :outerbottom)
+plot(seasonal, w=2 , color = "Black", lab = "Seasonal Component", legend = :outerbottom)
+
+```
+
+| ![quick_example_trend](./docs/assets/trend.svg) | ![quick_example_seas](./docs/assets/seasonal.svg)|
+|:------------------------------:|:-----------------------------:|
+
+
 ### Best Subset Selection
 Quick example on how to perform best subset selection in time series utilizing StateSpaceLearning.
 
 ```julia
 using StateSpaceLearning
 using CSV
 using DataFrames
-using Plots
 using Random
 
 Random.seed!(2024)
@@ -98,18 +119,15 @@ X = rand(length(log_air_passengers), 10) # Create 10 exogenous features
 
 y = log_air_passengers + X[:, 1:3]*β # add to the log_air_passengers series a contribution from only 3 exogenous features.
 
-plot(y)
-
-output = StateSpaceLearning.fit_model(y; Exogenous_X = X, estimation_input = Dict("α" => 1.0, "information_criteria" => "bic", "ϵ" => 0.05, 
-                                                   "penalize_exogenous" => true, "penalize_initial_states" => true))
+output = StateSpaceLearning.fit_model(y; Exogenous_X = X, estimation_input = Dict("α" => 1.0, "information_criteria" => "bic", "ϵ" => 0.05, "penalize_exogenous" => true, "penalize_initial_states" => true))
 
 Selected_exogenous = output.components["Exogenous_X"]["Selected"]
 
 ```
 
 In this example, the selected exogenous features were 1, 2, 3, as expected.
 
-### Completion of missing values
+### Missing values imputation
 Quick example of completion of missing values for the air passengers time-series (artificial NaN values are added to the original time-series).
 
 ```julia
@@ -135,6 +153,30 @@ plot!(fitted_completed_missing_values, lab = "Fit in Sample completed values", w
 ```
 ![quick_example_completion_airp](./docs/assets/quick_example_completion_airp.PNG)
 
+### Outlier Detection
+Quick example of outlier detection for an altered air passengers time-series (artificial NaN values are added to the original time-series).
+
+```julia
+using CSV
+using DataFrames
+using Plots
+
+airp = CSV.File(StateSpaceLearning.AIR_PASSENGERS) |> DataFrame
+log_air_passengers = log.(airp.passengers)
+
+log_air_passengers[60] = 10
+log_air_passengers[30] = 1
+log_air_passengers[100] = 2
+
+output = StateSpaceLearning.fit_model(log_air_passengers)
+detected_outliers = findall(i -> i != 0, output.components["o"]["Coefs"])
+
+plot(log_air_passengers, w=2 , color = "Black", lab = "Historical", legend = :outerbottom)
+scatter!([detected_outliers], log_air_passengers[detected_outliers], lab = "Detected Outliers")
+
+```
+![quick_example_completion_airp](./docs/assets/outlier.svg)
+
 ### StateSpaceModels initialization
 Quick example on how to use StateSpaceLearning to initialize  StateSpaceModels
 
@@ -168,7 +210,6 @@ To reproduce M4 paper results you can clone the repository and run the following
 ```shell
 julia paper_tests/m4_test/m4_test.jl
 python paper_tests/m4_test/m4_test.py
-1
 ```
 
 The results for SSL model in terms of MASE and sMAPE for all 48000 series will be stored in folder "paper_tests/m4_test/results_SSL". The average results of MASE, sMAPE and OWA will be saved in file "paper_tests/m4_test/metric_results/SSL_METRICS_RESULTS.csv".

docs/assets/outlier.svg

@@ -14,11 +14,12 @@ model_input = Dict()
 ```
 
 ### create_X
-The create_X function constructs the matrices in the State Space Learning format. It must accept the following inputs: (model_input::Dict, Exogenous_X::Matrix{Fl}, outlier::Bool, ζ_ω_threshold::Int64, T::Int64, steps_ahead::Int64=0, Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, size(Exogenous_X, 2))). This function may not use parameters such as outlier, ζ_ω_threshold, or Exogenous_X. It must return a matrix.
+The create_X function constructs the matrices in the State Space Learning format. It must accept the following inputs: (model_input::Dict, Exogenous_X::Matrix{Fl}, steps_ahead::Int64=0, Exogenous_Forecast::Matrix{Fl}). It must return a matrix.
 
 ```julia
-function create_X_LocalLevel(model_input::Dict, Exogenous_X::Matrix{Fl}, outlier::Bool, ζ_ω_threshold::Int64, T::Int64,
+function create_X_LocalLevel(model_input::Dict, Exogenous_X::Matrix{Fl},
                   steps_ahead::Int64=0, Exogenous_Forecast::Matrix{Fl}=zeros(steps_ahead, size(Exogenous_X, 2))) where Fl
+    T = size(Exogenous_X, 1)
     initial_states_matrix = ones(T+steps_ahead, 1)
     ξ_matrix = Matrix{Float64}(undef, T+steps_ahead, T - 1)
     for t in 1:T+steps_ahead
@@ -30,10 +31,11 @@ end
 ```
 
 ### get_components_indexes
-The get_components_indexes function outputs a dictionary indicating the indexes of each model component, including a set of indexes for all initial states. For the Local Level Model, the only components are the initial state μ1 and its innovations ξ. The function must accept the following inputs: (T::Int64, Exogenous_X::Matrix{Fl}, model_input::Dict, outlier::Bool, ζ_ω_threshold::Int64). This function may not use parameters such as outlier, ζ_ω_threshold, or Exogenous_X. It must return a dictionary.
+The get_components_indexes function outputs a dictionary indicating the indexes of each model component, including a set of indexes for all initial states. For the Local Level Model, the only components are the initial state μ1 and its innovations ξ. The function must accept the following inputs: (Exogenous_X::Matrix{Fl}, model_input::Dict). It must return a dictionary.
 
 ```julia
-function get_components_indexes_LocalLevel(T::Int64, Exogenous_X::Matrix{Fl}, model_input::Dict, outlier::Bool, ζ_ω_threshold::Int64)::Dict where Fl
+function get_components_indexes_LocalLevel(Exogenous_X::Matrix{Fl}, model_input::Dict)::Dict where Fl
+    T = size(Exogenous_X, 1)
     μ1_indexes = [1]
     initial_states_indexes = [1]
     ξ_indexes = collect(2:T)

docs/assets/seasonal.svg

@@ -4,7 +4,6 @@
 |:-----------------:|:-----------------:|:-----------------:|
 | [![ci](https://github.com/LAMPSPUC/StateSpaceLearning.jl/actions/workflows/ci.yml/badge.svg)](https://github.com/LAMPSPUC/StateSpaceLearning.jl/actions/workflows/ci.yml) | [![codecov](https://codecov.io/gh/LAMPSPUC/StateSpaceLearning.jl/graph/badge.svg?token=VDpuXvPSI2)](https://codecov.io/gh/LAMPSPUC/StateSpaceLearning.jl) | [![](https://img.shields.io/badge/docs-latest-blue.svg)](https://lampspuc.github.io/StateSpaceLearning.jl/latest/)
 
-
 StateSpaceLearning.jl is a package for modeling and forecasting time series in a high-dimension regression framework.
 
 ## Quickstart
@@ -18,17 +17,15 @@ y = randn(100)
 output = StateSpaceLearning.fit_model(y)
 
 #Main output options 
-model_type = output.model_input # State Space Equivalent model utilized in the estimation (default = Basic Structural).
+model_input         = output.model_input # Model inputs that were utilized to build the regression matrix.
+Create_X            = output.Create_X # The function utilized to build the regression matrix.
 X                   = output.X # High Dimension Regression utilized in the estimation.
 coefs               = output.coefs # High Dimension Regression coefficients estimated in the estimation.
 ϵ                   = output.ϵ # Residuals of the model.
 fitted              = output.fitted # Fit in Sample of the model.
 components          = output.components # Dictionary containing information about each component of the model, each component has the keys: "Values" (The value of the component in each timestamp) , "Coefs" (The coefficients estimated for each element of the component) and "Indexes" (The indexes of the elements of the component in the high dimension regression "X").
 residuals_variances = output.residuals_variances # Dictionary containing the estimated variances for the innovations components (that is the information that can be utilized to initialize the state space model).
-T                   = output.T # The length of the original time series.
-outlier             = output.outlier # Boolean indicating the presence of outlier component (default = false).
 valid_indexes       = output.valid_indexes # Vector containing valid indexes of the time series (non valid indexes represent NaN values in the time series).
-ζ_ω_threshold         = output.ζ_ω_threshold # ζ_ω_threshold parameter (default = 0). A non 0 value for this parameter might be important in terms of forecast for some time series to lead to more stable predictions (we recommend ζ_ω_threshold = 11 for monthly series).
 
 #Forecast
 prediction = StateSpaceLearning.forecast(output, 12) #Gets a 12 steps ahead prediction
@@ -38,11 +35,9 @@ prediction = StateSpaceLearning.forecast(output, 12) #Gets a 12 steps ahead pred
 ## Fit Arguments
 
 * `y::Vector{Fl}`: Vector of data.
-* `model_input::Dict`: Dictionary containing the model input parameters (default: Dict("level" => true, "stochastic_level" => true, "trend" => true, "stochastic_trend" => true, "seasonal" => true, "stochastic_seasonal" => true, "freq_seasonal" => 12)).
+* `model_input::Dict`: Dictionary containing the model input parameters (default: Dict("level" => true, "stochastic_level" => true, "trend" => true, "stochastic_trend" => true, "seasonal" => true, "stochastic_seasonal" => true, "freq_seasonal" => 12, "outlier" => true, "ζ_ω_threshold" => 12)).
 * `estimation_input::Dict`: Dictionary containing the estimation input parameters (default: Dict("α" => 0.1, "information_criteria" => "aic", ψ => 0.05, "penalize_exogenous" => true, "penalize_initial_states" => true)).
 * `Exogenous_X::Union{Matrix{Fl}, Missing}`: Exogenous variables matrix (default: missing).
-* `outlier::Bool`: Flag for considering outlier component (default: true).
-* `ζ_ω_threshold::Int64`: ζ_ω_threshold parameter (default: 12).
 
 ## Features
 
@@ -52,9 +47,8 @@ Current features include:
 * Forecasting
 * Completion of missing values
 * Predefined models, including:
-  * Basic Structural"
-  * Local Linear Trend
-  * Local Level
+* Outlier detection
+* Outlier robust models
 
 ## Quick Examples
 
@@ -71,23 +65,49 @@ log_air_passengers = log.(airp.passengers)
 steps_ahead = 30
 
 output = StateSpaceLearning.fit_model(log_air_passengers)
-prediction_raw = StateSpaceLearning.forecast(output, steps_ahead)
-prediction = exp.(prediction_raw)
+prediction_log = StateSpaceLearning.forecast(output, steps_ahead)
+prediction = exp.(prediction_log)
 
 plot(airp.passengers, w=2 , color = "Black", lab = "Historical", legend = :outerbottom)
 plot!(vcat(ones(output.T).*NaN, prediction), lab = "Forcast", w=2, color = "blue")
 
 ```
 ![quick_example_airp](./docs/assets/quick_example_airp.PNG)
 
+### Component Extraction
+Quick example on how to perform component extraction in time series utilizing StateSpaceLearning.
+
+```julia
+using CSV
+using DataFrames
+using Plots
+
+airp = CSV.File(StateSpaceLearning.AIR_PASSENGERS) |> DataFrame
+log_air_passengers = log.(airp.passengers)
+
+output = StateSpaceLearning.fit_model(log_air_passengers)
+
+level = output.components["μ1"]["Values"] + output.components["ξ"]["Values"]
+slope = output.components["ν1"]["Values"] + output.components["ζ"]["Values"]
+seasonal = output.components["γ1"]["Values"] + output.components["ω"]["Values"]
+trend = level + slope
+
+plot(trend, w=2 , color = "Black", lab = "Trend Component", legend = :outerbottom)
+plot(seasonal, w=2 , color = "Black", lab = "Seasonal Component", legend = :outerbottom)
+
+```
+
+| ![quick_example_trend](./docs/assets/trend.svg) | ![quick_example_seas](./docs/assets/seasonal.svg)|
+|:------------------------------:|:-----------------------------:|
+
+
 ### Best Subset Selection
 Quick example on how to perform best subset selection in time series utilizing StateSpaceLearning.
 
 ```julia
 using StateSpaceLearning
 using CSV
 using DataFrames
-using Plots
 using Random
 
 Random.seed!(2024)
@@ -99,18 +119,15 @@ X = rand(length(log_air_passengers), 10) # Create 10 exogenous features
 
 y = log_air_passengers + X[:, 1:3]*β # add to the log_air_passengers series a contribution from only 3 exogenous features.
 
-plot(y)
-
-output = StateSpaceLearning.fit_model(y; Exogenous_X = X, estimation_input = Dict("α" => 1.0, "information_criteria" => "bic", "ϵ" => 0.05, 
-                                                   "penalize_exogenous" => true, "penalize_initial_states" => true))
+output = StateSpaceLearning.fit_model(y; Exogenous_X = X, estimation_input = Dict("α" => 1.0, "information_criteria" => "bic", "ϵ" => 0.05, "penalize_exogenous" => true, "penalize_initial_states" => true))
 
 Selected_exogenous = output.components["Exogenous_X"]["Selected"]
 
 ```
 
 In this example, the selected exogenous features were 1, 2, 3, as expected.
 
-### Completion of missing values
+### Missing values imputation
 Quick example of completion of missing values for the air passengers time-series (artificial NaN values are added to the original time-series).
 
 ```julia
@@ -136,6 +153,30 @@ plot!(fitted_completed_missing_values, lab = "Fit in Sample completed values", w
 ```
 ![quick_example_completion_airp](./docs/assets/quick_example_completion_airp.PNG)
 
+### Outlier Detection
+Quick example of outlier detection for an altered air passengers time-series (artificial NaN values are added to the original time-series).
+
+```julia
+using CSV
+using DataFrames
+using Plots
+
+airp = CSV.File(StateSpaceLearning.AIR_PASSENGERS) |> DataFrame
+log_air_passengers = log.(airp.passengers)
+
+log_air_passengers[60] = 10
+log_air_passengers[30] = 1
+log_air_passengers[100] = 2
+
+output = StateSpaceLearning.fit_model(log_air_passengers)
+detected_outliers = findall(i -> i != 0, output.components["o"]["Coefs"])
+
+plot(log_air_passengers, w=2 , color = "Black", lab = "Historical", legend = :outerbottom)
+scatter!([detected_outliers], log_air_passengers[detected_outliers], lab = "Detected Outliers")
+
+```
+![quick_example_completion_airp](./docs/assets/outlier.svg)
+
 ### StateSpaceModels initialization
 Quick example on how to use StateSpaceLearning to initialize  StateSpaceModels
 
@@ -169,7 +210,6 @@ To reproduce M4 paper results you can clone the repository and run the following
 ```shell
 julia paper_tests/m4_test/m4_test.jl
 python paper_tests/m4_test/m4_test.py
-1
 ```
 
 The results for SSL model in terms of MASE and sMAPE for all 48000 series will be stored in folder "paper_tests/m4_test/results_SSL". The average results of MASE, sMAPE and OWA will be saved in file "paper_tests/m4_test/metric_results/SSL_METRICS_RESULTS.csv".