Pearson res (#128)

* Add backtest and time_limit_sec

* Add backtest

* Add JarqueBera on pearson residuals

Author: guilhermebodin

Project.toml

@@ -1,10 +1,11 @@
 name = "ScoreDrivenModels"
 uuid = "4a87933e-d659-11e9-0e65-7f40dedd4a3a"
 authors = ["guilhermebodin <[email protected]>, raphaelsaavedra <[email protected]>"]
-version = "0.1.6"
+version = "0.1.7"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -15,6 +16,7 @@ Distributions = "0.23"
 Optim = "0.20, 0.21, 0.22, 1.2"
 SpecialFunctions = "0.8"
 julia = "1"
+HypothesisTests = "0.9"
 
 [extras]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

src/MLE.jl

@@ -13,6 +13,7 @@ struct Fitted{D <: Distribution, T <: AbstractFloat}
     llk::T
     coefs::Vector{T}
     numerical_hessian::Matrix{T}
+    pearson_residuals::Vector{T}
 end
 
 struct CoefsStats{T <: AbstractFloat}
@@ -29,13 +30,16 @@ struct EstimationStats{D <: Distribution, T <: AbstractFloat}
     aic::T
     bic::T
     np::T
+    jarquebera_p_value::T
     coefs_stats::CoefsStats{T}
 end
 
 function results(f::Fitted{D, T}) where {D, T}
     estim_results = eval_coefs_stats(f)
     np = length(f.unknowns)
-    return EstimationStats{D, T}(f.num_obs, f.llk, f.aic, f.bic, np, estim_results)
+    jarquebera_p_value = pvalue(JarqueBeraTest(f.pearson_residuals))
+    return EstimationStats{D, T}(f.num_obs, f.llk, f.aic, f.bic, np, 
+                                 jarquebera_p_value, estim_results)
 end
 
 function eval_coefs_stats(f::Fitted{D, T}) where {D, T}
@@ -91,13 +95,18 @@ function update_aux_estimation!(aux_est::AuxEstimation{T}, func::Optim.TwiceDiff
     return
 end
 
-function fit(gas::Model{D, T}, y::Vector{T};
+function fit!(gas::Model{D, T}, y::Vector{T};
              initial_params::Matrix{T} = DEFAULT_INITIAL_PARAM,
              opt_method::AbstractOptimizationMethod = NelderMead(gas, DEFAULT_NUM_SEEDS),
              verbose::Int = DEFAULT_VERBOSE,
              throw_errors::Bool = false,
              time_limit_sec::Int = 10^8) where {D, T}
 
+    unknowns = find_unknowns(gas)
+    # Check if the model has no unknowns
+    n_unknowns = length(unknowns)
+    check_model_estimated(n_unknowns) && return gas
+
     verbose in [0, 1, 2, 3] || throw(ErrorException, "verbose argument must be in [0, 1, 2, 3]")
     # Number of initial_points and number of params to estimate
     n_initial_points = length(opt_method.initial_points)
@@ -150,21 +159,12 @@ function fit(gas::Model{D, T}, y::Vector{T};
         println(aux_est.opt_result[best_seed])
     end
 
-    return Fitted{D, T}(n, unknowns, aic, bic, best_llk, coefs, num_hessian)
-end
+    fill_psitilde!(gas, coefs, unknowns)
 
-function fit!(gas::Model{D, T}, y::Vector{T};
-              initial_params::Matrix{T} = DEFAULT_INITIAL_PARAM,
-              opt_method::AbstractOptimizationMethod = NelderMead(gas, DEFAULT_NUM_SEEDS),
-              verbose::Int = DEFAULT_VERBOSE,
-              throw_errors::Bool = false) where {D, T}
+    # Calculate pearson residuals
+    pearson_res = isnan(initial_params[1]) ? 
+                pearson_residuals(y, gas) :
+                pearson_residuals(y, gas; initial_params = initial_params)
 
-    unknowns = find_unknowns(gas)
-    # Check if the model has no unknowns
-    n_unknowns = length(unknowns)
-    check_model_estimated(n_unknowns) && return gas
-    
-    f = fit(gas, y; initial_params = initial_params, opt_method = opt_method, verbose = verbose, throw_errors = throw_errors)
-    fill_psitilde!(gas, f.coefs, unknowns)
-    return f
+    return Fitted{D, T}(n, unknowns, aic, bic, best_llk, coefs, num_hessian, pearson_res)
 end

src/ScoreDrivenModels.jl

@@ -1,6 +1,6 @@
 module ScoreDrivenModels
 
-using Distributions, Optim, SpecialFunctions
+using Distributions, Optim, SpecialFunctions, HypothesisTests
 using LinearAlgebra, Printf
 
 import Base: length, deepcopy, show

src/prints.jl

@@ -6,6 +6,7 @@ function Base.show(io::IO, est::EstimationStats{D, T}) where {D, T}
     println(io, "Log-likelihood:               ", @sprintf("%.4f", est.loglikelihood))
     println(io, "AIC:                          ", @sprintf("%.4f", est.aic))
     println(io, "BIC:                          ", @sprintf("%.4f", est.bic))
+    println(io, "Jarque Bera p value:          ", @sprintf("%.4f", est.jarquebera_p_value))
     print_coefs_stats(est.coefs_stats)
     return nothing
 end