return some estimation results (#44)

Author: guilhermebodin

src/MLE.jl

@@ -1,42 +1,61 @@
 export estimate!
 
 const DEFAULT_INITIAL_PARAM = NaN.*ones(1, 1)
+const DEFAULT_NUM_SEEDS = 3
+
+struct EstimationResults{T <: AbstractFloat}
+    aic::T
+    bic::T
+    llk::T
+    minimizer::Vector{T}
+    numerical_hessian::Matrix{T}
+end
 
 function log_lik(psitilde::Vector{T}, y::Vector{T}, sdm::SDM{D, T}, 
                  initial_params::Vector{Vector{T}}, unknowns::Unknowns_SDM, n::Int) where {D, T}
     return error("log_lik not defined for a model of type ", typeof(sdm))
 end
 
+function AIC(n_unknowns::Int, log_lik::T) where T
+    return T(2 * n_unknowns - 2 * log_lik)
+end
+
+function BIC(n::Int, n_unknowns::Int, log_lik::T) where T
+    return T(log(n) * n_unknowns - 2 * log_lik)
+end
+
 function estimate!(sdm::SDM{D, T}, y::Vector{T};
                    initial_params::Matrix{T} = DEFAULT_INITIAL_PARAM,
-                   opt_method::AbstractOptimizationMethod = LBFGS(sdm, 3),
+                   opt_method::AbstractOptimizationMethod = LBFGS(sdm, DEFAULT_NUM_SEEDS),
                    verbose::Int = 0) where {D, T}
 
     # Number of seed and number of params to estimate
     nseeds = length(opt_method.seeds)
     n = length(y)
 
     unknowns = find_unknowns(sdm)
-    len_unknowns = length(unknowns)
+    n_unknowns = length(unknowns)
 
     # Check if the model has no unknowns
-    check_model_estimated(len_unknowns) && return sdm
+    check_model_estimated(n_unknowns) && return sdm
 
     # optimize for each seed
-    psi = Vector{Vector{Float64}}(undef, 0)
-    loglikelihood = Vector{Float64}(undef, 0)
+    psi = Vector{Vector{T}}(undef, 0)
+    numerical_hessian = Vector{Matrix{T}}(undef, 0)
+    loglikelihood = Vector{T}(undef, 0)
     optseeds = Vector{Optim.OptimizationResults}(undef, 0)
 
     for i = 1:nseeds
         try 
-            optseed = optimize(psi_tilde -> log_lik(psi_tilde, y, sdm, initial_params, unknowns, n), 
-                                                    opt_method.seeds[i],
-                                                    opt_method.method, Optim.Options(f_tol = opt_method.f_tol, 
-                                                                                     g_tol = opt_method.g_tol, 
-                                                                                     iterations = opt_method.iterations,
-                                                                                     show_trace = (verbose == 2 ? true : false) ))
+            func = TwiceDifferentiable(psi_tilde -> log_lik(psi_tilde, y, sdm, initial_params, unknowns, n), opt_method.seeds[i])
+            optseed = optimize(func, opt_method.seeds[i],
+                                     opt_method.method, Optim.Options(f_tol = opt_method.f_tol, 
+                                                                      g_tol = opt_method.g_tol, 
+                                                                      iterations = opt_method.iterations,
+                                                                      show_trace = (verbose == 2 ? true : false) ))
             push!(loglikelihood, -optseed.minimum)
             push!(psi, optseed.minimizer)
+            push!(numerical_hessian, Optim.hessian!(func, optseed.minimizer))
             push!(optseeds, optseed)
             println("seed $i of $nseeds - $(-optseed.minimum)")
         catch err
@@ -50,15 +69,20 @@ function estimate!(sdm::SDM{D, T}, y::Vector{T};
         return
     end
 
-    best_llk = argmax(loglikelihood)
-    bestpsi = psi[best_llk]
+    best_llk, best_seed = findmax(loglikelihood)
+    num_hessian = numerical_hessian[best_seed]
+    best_psi = psi[best_seed]
+    aic = AIC(n_unknowns, best_llk)
+    bic = BIC(n, n_unknowns, best_llk)
+
     if verbose >= 1
         println("\nBest seed optimization result:")
-        println(optseeds[best_llk])
+        println(optseeds[best_seed])
     end
 
     # return the estimated 
-    fill_psitilde!(sdm, bestpsi, unknowns)
+    fill_psitilde!(sdm, best_psi, unknowns)
 
     println("Finished!")
+    return EstimationResults{T}(aic, bic, best_llk, best_psi, num_hessian)
 end

src/gas/gas.jl

@@ -133,7 +133,7 @@ function log_lik(psitilde::Vector{T}, y::Vector{T}, gas::GAS{D, T},
     if isnan(initial_params[1]) # Means default stationary initialization
         params = score_driven_recursion(gas, y)
     else
-        params = score_driven_recursion(gas, y, initial_params)
+        params = score_driven_recursion(gas, y; initial_params = initial_params)
     end
 
     return log_likelihood(D, y, params, n)