Estimation statistics and prints (#57)

* add fit_stats method

* add estimation results prints

* change NelderMead to the default algorithm

Author: guilhermebodin

Project.toml

@@ -7,6 +7,7 @@ version = "0.1.0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]

examples/garch.jl

@@ -51,10 +51,14 @@ gas = GAS(1, 1, Normal, 1.0, time_varying_params = [2])
 initial_point = [0.0; 0.5; 0.25; 0.75]
 
 # Estimate the model
-f = fit(gas, y; initial_params = initial_params, 
+f_ip_newton = fit(gas, y; initial_params = initial_params, 
                   opt_method = IPNewton(gas, [initial_point]; ub = ub, lb = lb))
 
+estimation_stats_ip_newton = fit_stats(f_ip_newton)
+
 # Another way to estimate the model is to rely on luck and give many random initial_points
 # for a non constrained optimization 
-f = fit(gas, y; initial_params = initial_params, 
-                  opt_method = NelderMead(gas, 100))
+f_nelder_mead = fit(gas, y; initial_params = initial_params, 
+                  opt_method = NelderMead(gas, 100))
+
+estimation_stats_nelder_mead = fit_stats(f_nelder_mead)

src/MLE.jl

@@ -1,17 +1,60 @@
-export fit, fit!
+export fit, fit!, fit_stats
 
 const DEFAULT_INITIAL_PARAM = NaN.*ones(1, 1)
 const DEFAULT_NUM_SEEDS = 3
 const DEFAULT_VERBOSE = 0
+const VARIANCE_ZERO = 1e-10
 
 struct FittedSDM{T <: AbstractFloat}
+    unknowns::UnknownsSDM
     aic::T
     bic::T
     llk::T
     coefs::Vector{T}
     numerical_hessian::Matrix{T}
 end
 
+struct CoefsStatsSDM{T <: AbstractFloat}
+    unknowns::UnknownsSDM
+    coefs::Vector{T}
+    std_errors::Vector{T}
+    t_stat::Vector{T}
+    p_values::Vector{T}
+end
+
+struct EstimationStatsSDM{T <: AbstractFloat}
+    loglikelihood::T
+    aic::T
+    bic::T
+    np::T
+    coefs_stats::CoefsStatsSDM{T}
+end
+
+function fit_stats(f::FittedSDM{T}) where T
+    estim_results = eval_coefs_stats(f)
+    np = length(f.unknowns)
+    return EstimationStatsSDM{T}(f.llk, f.aic, f.bic, np, estim_results)
+end
+
+function eval_coefs_stats(f::FittedSDM{T}) where T
+    np = length(f.unknowns)
+    inv_H = inv(f.numerical_hessian)
+    vars = diag(inv_H)
+    for i in eachindex(vars)
+        if vars[i] <= VARIANCE_ZERO 
+            vars[i] = VARIANCE_ZERO
+        end
+    end
+    std_errors = sqrt.(vars)
+    t_stats = f.coefs ./ std_errors
+
+    # Calculate p-values of the the t statistics
+    t_dist = TDist(np)
+    p_values = 1 .- 2*abs.(cdf.(t_dist, t_stats) .- 0.5)
+
+    return CoefsStatsSDM{T}(f.unknowns, f.coefs, std_errors, t_stats, p_values)
+end
+
 mutable struct AuxEstimation{T <: AbstractFloat}
     psi::Vector{Vector{T}}
     numerical_hessian::Vector{Matrix{T}}
@@ -48,7 +91,7 @@ end
 
 function fit(sdm::SDM{D, T}, y::Vector{T};
              initial_params::Matrix{T} = DEFAULT_INITIAL_PARAM,
-             opt_method::AbstractOptimizationMethod = LBFGS(sdm, DEFAULT_NUM_SEEDS),
+             opt_method::AbstractOptimizationMethod = NelderMead(sdm, DEFAULT_NUM_SEEDS),
              verbose::Int = DEFAULT_VERBOSE) where {D, T}
 
     # Number of seed and number of params to estimate
@@ -96,7 +139,7 @@ function fit(sdm::SDM{D, T}, y::Vector{T};
     end
 
     println("Finished!")
-    return FittedSDM{T}(aic, bic, best_llk, coefs, num_hessian)
+    return FittedSDM{T}(unknowns, aic, bic, best_llk, coefs, num_hessian)
 end
 
 function fit!(sdm::SDM{D, T}, y::Vector{T};

src/ScoreDrivenModels.jl

@@ -1,9 +1,9 @@
 module ScoreDrivenModels
 
 using Distributions, Optim, SpecialFunctions
-using LinearAlgebra
+using LinearAlgebra, Printf
 
-import Base.length, Base.deepcopy
+import Base: length, deepcopy, show
 
 const SCALINGS  = [0.0, 1/2, 1.0]
 const BIG_NUM = 1e8
@@ -15,6 +15,7 @@ include("utils.jl")
 include("link_functions.jl")
 include("score.jl")
 include("MLE.jl")
+include("prints.jl")
 
 # Optimization methods
 include("opt_methods/common_methods.jl")

src/prints.jl

@@ -0,0 +1,62 @@
+function Base.show(io::IO, est::EstimationStatsSDM)
+    println("--------------------------------------------------------")
+    println("log-likelihood: ", @sprintf("%.4f", est.loglikelihood))
+    println("            np: ", Int(est.np))
+    println("           AIC: ", @sprintf("%.4f", est.aic))
+    println("           BIC: ", @sprintf("%.4f", est.bic))
+    print_coefs_stats(est.coefs_stats)
+    return nothing
+end
+
+function print_coefs_stats(coefs_stats)
+    println("--------------------------------------------------------")
+    println("Parameter      Estimate   Std.Error     t stat   p-value")
+    offset = 1
+    for i in coefs_stats.unknowns.ω
+        p_c, p_std, p_t_stat, p_p_val = print_coefs_sta(coefs_stats, offset)
+        p = build_print(p_c, p_std, p_t_stat, p_p_val, "ω_$i")
+        println(p)
+        offset += 1
+    end
+    for k in sort(collect(keys(coefs_stats.unknowns.A)))
+        for i in coefs_stats.unknowns.A[k]
+            p_c, p_std, p_t_stat, p_p_val = print_coefs_sta(coefs_stats, offset)
+            ind = round(Int, sqrt(i))
+            p = build_print(p_c, p_std, p_t_stat, p_p_val, "A_$(k)_$ind$ind")
+            println(p)
+            offset += 1
+        end
+    end
+    for k in sort(collect(keys(coefs_stats.unknowns.B)))
+        for i in coefs_stats.unknowns.B[k]
+            p_c, p_std, p_t_stat, p_p_val = print_coefs_sta(coefs_stats, offset)
+            ind = round(Int, sqrt(i))
+            p = build_print(p_c, p_std, p_t_stat, p_p_val, "B_$(k)_$ind$ind")
+            println(p)
+            offset += 1
+        end
+    end
+    return nothing
+end
+
+function print_coefs_sta(sta, offset::Int)
+    p_c      = @sprintf("%.4f", sta.coefs[offset])
+    p_std    = @sprintf("%.4f", sta.std_errors[offset])
+    p_t_stat = @sprintf("%.4f", sta.t_stat[offset])
+    p_p_val  = @sprintf("%.4f", sta.p_values[offset])
+    return p_c, p_std, p_t_stat, p_p_val
+end
+
+function build_print(p_c, p_std, p_t_stat, p_p_val, param)
+    p = ""
+    p *= param
+    p *= " "^max(0, 23 - length(p_c) - length(param))
+    p *= p_c
+    p *= " "^max(0, 12 - length(p_std))
+    p *= p_std
+    p *= " "^max(0, 11 - length(p_t_stat))
+    p *= p_t_stat
+    p *= " "^max(0, 10 - length(p_p_val))
+    p *= p_p_val
+    return p
+end