Merge pull request #279 from LAMPSPUC/auto_arima

iagochavarry · web-flow · commit b1723e74f3c6 · 2021-06-15T16:26:01.000-03:00
Auto arima
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "StateSpaceModels"
 uuid = "99342f36-827c-5390-97c9-d7f9ee765c78"
 authors = ["raphaelsaavedra <raphael.saavedra93@gmail.com>, guilhermebodin <guilherme.b.moraes@gmail.com>, mariohsouto"]
-version = "0.5.16"
+version = "0.5.17"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
@@ -12,6 +12,7 @@ OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+SeasonalTrendLoess = "42fb36cb-998a-4034-bf40-4eee476c43a1"
 ShiftedArrays = "1277b4bf-5013-50f5-be3d-901d8477a67a"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
@@ -23,6 +24,7 @@ Optim = "0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 1.2"
 OrderedCollections = "1"
 Polynomials = "1"
 RecipesBase = "1"
+SeasonalTrendLoess = "0.1"
 ShiftedArrays = "1"
 StatsBase = "0.27, 0.28, 0.29, 0.30, 0.31, 0.32, 0.33"
 julia = "1"
diff --git a/docs/src/manual.md b/docs/src/manual.md
@@ -76,6 +76,7 @@ In any case if the documentation explains the procedures and indicates if there
 
 ```
 auto_ets
+auto_arima
 ```
 
 ### Implementing a custom model
diff --git a/src/StateSpaceModels.jl b/src/StateSpaceModels.jl
@@ -13,6 +13,7 @@ using Optim
 using OrderedCollections
 using RecipesBase
 using StatsBase
+using SeasonalTrendLoess
 
 abstract type StateSpaceModel end
 
@@ -21,7 +22,8 @@ include("datasets.jl")
 include("hyperparameters.jl")
 include("systems.jl")
 include("kalman_filter_and_smoother.jl")
-
+include("statistical_tests/kpss.jl")
+include("statistical_tests/canova_hensen.jl")
 include("filters/univariate_kalman_filter.jl")
 include("filters/multivariate_kalman_filter.jl")
 include("filters/scalar_kalman_filter.jl")
@@ -81,6 +83,7 @@ export UnivariateKalmanFilter
 export UnobservedComponents
 
 # Exported functions
+export auto_arima
 export auto_ets
 export cross_validation
 export constrain_box!
diff --git a/src/fit.jl b/src/fit.jl
@@ -30,6 +30,7 @@ function fit!(
     filter::KalmanFilter=default_filter(model),
     optimizer::Optimizer=default_optimizer(model),
 )
+    isfitted(model) && return model
     @assert has_fit_methods(typeof(model))
     initial_unconstrained_hyperparameter = handle_optim_initial_hyperparameters(model)
     # TODO Should there be a try catch?
@@ -144,26 +145,30 @@ function fill_results!(model::StateSpaceModel, llk::Fl, std_err::Vector{Fl}) whe
     model.results.model_name = model_name(model)
     model.results.coef_table = coef_table
     model.results.llk = llk
-    model.results.aic = AIC(num_hyperparameters, llk)
-    model.results.aicc = AICc(n_obs, num_hyperparameters, llk)
-    model.results.bic = BIC(n_obs, num_hyperparameters, llk)
+    model.results.aic = aic(num_hyperparameters, llk)
+    model.results.aicc = aicc(n_obs, num_hyperparameters, llk)
+    model.results.bic = bic(n_obs, num_hyperparameters, llk)
     model.results.num_observations = n_obs
     model.results.num_hyperparameters = num_hyperparameters
     return model
 end
 
-function AIC(k::Int, llk::Fl) where Fl
+function aic(k::Int, llk::Fl) where Fl
     return convert(Fl, 2 * k - 2 * llk)
 end
 
-function AICc(n::Int, k::Int, llk::Fl) where Fl
-    return convert(Fl, AIC(k, llk) + (2 * k * (k + 1) / (n - k - 1)))
+function aicc(n::Int, k::Int, llk::Fl) where Fl
+    return convert(Fl, aic(k, llk) + (2 * k * (k + 1) / (n - k - 1)))
 end
 
-function BIC(n::Int, k::Int, llk::Fl) where Fl
+function bic(n::Int, k::Int, llk::Fl) where Fl
     return convert(Fl, log(n) * k - 2 * llk)
 end
 
+aic(model::StateSpaceModel) = model.results.aic
+aicc(model::StateSpaceModel) = model.results.aicc
+bic(model::StateSpaceModel) = model.results.bic
+
 function build_coef_table(model::StateSpaceModel, std_err::Vector{Fl}) where Fl
     all_coef = get_constrained_values(model)
     all_std_err = handle_std_err(model, std_err)
diff --git a/src/forecast.jl b/src/forecast.jl
@@ -39,7 +39,7 @@ function forecast(
     # Associate with the model hyperparameters
     forecasting_model.hyperparameters = model_hyperparameters
     # Perform the kalman filter
-    fo = kalman_filter(forecasting_model)
+    fo = kalman_filter(forecasting_model; filter = filter)
     # fill forecast matrices
     expected_value = Vector{Vector{Fl}}(undef, steps_ahead)
     covariance = Vector{Matrix{Fl}}(undef, steps_ahead)
@@ -82,7 +82,7 @@ function forecast(
     # Associate with the model hyperparameters
     forecasting_model.hyperparameters = model_hyperparameters
     # Perform the kalman filter
-    fo = kalman_filter(forecasting_model)
+    fo = kalman_filter(forecasting_model; filter = filter)
     # fill forecast matrices
     expected_value = Vector{Vector{Fl}}(undef, steps_ahead)
     covariance = Vector{Matrix{Fl}}(undef, steps_ahead)
@@ -119,7 +119,7 @@ function simulate_scenarios(
 )
     # Query the type of model elements
     Fl = typeof_model_elements(model)
-    fo = kalman_filter(model)
+    fo = kalman_filter(model; filter = filter)
     last_state = fo.a[end]
     num_series = size(model.system.y, 2)
 
diff --git a/src/models/basicstructural.jl b/src/models/basicstructural.jl
@@ -17,7 +17,7 @@ component. It is defined by:
 # Example
 ```jldoctest
 julia> model = BasicStructural(rand(100), 12)
-BasicStructural model
+BasicStructural
 ```
 
 # References
diff --git a/src/models/basicstructural_explanatory.jl b/src/models/basicstructural_explanatory.jl
@@ -17,7 +17,7 @@ It is defined by:
 # Example
 ```jldoctest
 julia> model = BasicStructuralExplanatory(rand(100), 12, rand(100, 2))
-BasicStructuralExplanatory model
+BasicStructuralExplanatory
 ```
 
 # References
diff --git a/src/models/basicstructural_multivariate.jl b/src/models/basicstructural_multivariate.jl
@@ -18,7 +18,7 @@ It is defined by:
 # Example
 ```jldoctest
 julia> model = MultivariateBasicStructural(rand(100, 2), 12)
-MultivariateBasicStructural model
+MultivariateBasicStructural
 ```
 
 # References
diff --git a/src/models/dar.jl b/src/models/dar.jl
@@ -22,7 +22,7 @@ A Dynamic Autorregressive model is defined by:
 # Example
 ```jldoctest
 julia> model = DAR(randn(100), 2)
-DAR model
+DAR
 ```
 """
 mutable struct DAR <: StateSpaceModel
diff --git a/src/models/exponential_smoothing.jl b/src/models/exponential_smoothing.jl
@@ -325,6 +325,6 @@ function auto_ets(y::Vector{Fl}; seasonal::Int = 0) where Fl
     fit!(m3)
     push!(models, m3)
     push!(models_aic, m3.results.aic)
-    best_aic_idx = findmin(models_aic)[2] # index of the best BIC
+    best_aic_idx = findmin(models_aic)[2] # index of the best bic
     return models[best_aic_idx]
 end
diff --git a/src/models/linear_regression.jl b/src/models/linear_regression.jl
@@ -23,7 +23,7 @@ The linear regression state-space model is defined by:
 # Example
 ```jldoctest
 julia> model = LinearRegression(rand(100, 2), rand(100))
-LinearRegression model
+LinearRegression
 ```
 """
 mutable struct LinearRegression <: StateSpaceModel
diff --git a/src/models/locallevel.jl b/src/models/locallevel.jl
@@ -14,7 +14,7 @@ The local level model is defined by:
 # Example
 ```jldoctest
 julia> model = LocalLevel(rand(100))
-LocalLevel model
+LocalLevel
 ```
 
 See more on [Nile river annual flow](@ref)
diff --git a/src/models/locallevelcycle.jl b/src/models/locallevelcycle.jl
@@ -16,7 +16,7 @@ The local level model with a cycle component is defined by:
 # Example
 ```jldoctest
 julia> model = LocalLevelCycle(rand(100))
-LocalLevelCycle model
+LocalLevelCycle
 ```
 
 See more on TODO RJ_TEMPERATURE
diff --git a/src/models/locallevelexplanatory.jl b/src/models/locallevelexplanatory.jl
@@ -17,7 +17,7 @@ A local level model with explanatory variables is defined by:
 # Example
 ```jldoctest
 julia> model = LocalLevelExplanatory(rand(100), rand(100, 1))
-LocalLevelExplanatory model
+LocalLevelExplanatory
 ```
 """
 mutable struct LocalLevelExplanatory <: StateSpaceModel
diff --git a/src/models/locallineartrend.jl b/src/models/locallineartrend.jl
@@ -13,7 +13,7 @@ The linear trend model is defined by:
 # Example
 ```jldoctest
 julia> model = LocalLinearTrend(rand(100))
-LocalLinearTrend model
+LocalLinearTrend
 ```
 
 See more on [Finland road traffic fatalities](@ref)
diff --git a/src/models/sarima.jl b/src/models/sarima.jl
diff --git a/src/prints.jl b/src/prints.jl
diff --git a/src/statistical_tests/canova_hensen.jl b/src/statistical_tests/canova_hensen.jl
diff --git a/src/statistical_tests/kpss.jl b/src/statistical_tests/kpss.jl
diff --git a/src/systems.jl b/src/systems.jl
diff --git a/test/models/locallevel.jl b/test/models/locallevel.jl
diff --git a/test/models/sarima.jl b/test/models/sarima.jl
