Maximize mean loglikelihood (#114)

* Update dependencies

* Maximize the mean llk

* Remove Beta distribution

Author: guilhermebodin

Project.toml

@@ -11,10 +11,10 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-julia = "1"
 Distributions = "0.23"
-Optim = "0.20"
+Optim = "0.20, 0.21, 0.22, 1.2"
 SpecialFunctions = "0.8"
+julia = "1"
 
 [extras]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

src/MLE.jl

@@ -120,7 +120,7 @@ function fit(gas::Model{D, T}, y::Vector{T};
                                                         opt_method.initial_points[i])
             opt_result = optimize(func, opt_method, verbose, i)
             update_aux_estimation!(aux_est, func, opt_result)
-            verbose >= 1 && println("Round $i of $n_initial_points - Log-likelihood: $(-opt_result.minimum)")
+            verbose >= 1 && println("Round $i of $n_initial_points - Log-likelihood: $(-opt_result.minimum * length(y))")
         catch err
             println(err)
             verbose >= 1 && println("Round $i diverged")
@@ -133,6 +133,8 @@ function fit(gas::Model{D, T}, y::Vector{T};
     end
 
     best_llk, best_seed = findmax(aux_est.loglikelihood)
+    # We optimize the average loglikelihood inside log_lik
+    best_llk = best_llk * length(y)
     num_hessian = aux_est.numerical_hessian[best_seed]
     coefs = aux_est.psi[best_seed]
     aic = AIC(n_unknowns, best_llk)

src/ScoreDrivenModels.jl

@@ -31,7 +31,6 @@ include("opt_methods/NelderMead.jl")
 # Distributions
 include("distributions/non_native_dists.jl")
 include("distributions/common_interface.jl")
-include("distributions/beta.jl")
 include("distributions/betalocationscale.jl")
 include("distributions/exponential.jl")
 include("distributions/gamma.jl")

src/distributions/beta.jl

@@ -27,10 +27,6 @@ function score!(score_til::Matrix{T}, y::T, ::Type{BetaLocationScale}, param::Ma
     return
 end
 
-function fisher_information!(aux::AuxiliaryLinAlg{T}, ::Type{BetaLocationScale}, param::Matrix{T}, t::Int) where T
-    return error("Fisher information not implemented for BetaLocationScale distribution.")
-end
-
 function log_likelihood(::Type{BetaLocationScale}, y::Vector{T}, param::Matrix{T}, n::Int) where T
     loglik = 0.0
     for t in 1:n

src/distributions/betalocationscale.jl

@@ -1,6 +1,5 @@
 # Currently supported distributions
 const DISTS = [
-    Beta,
     BetaLocationScale,
     Exponential,
     Gamma,
@@ -14,10 +13,7 @@ const DISTS = [
     Weibull
 ]
 
-export Beta,
-    BetaLocationScale,
-    Chi,
-    Chisq,
+export BetaLocationScale,
     Exponential,
     Gamma,
     LogitNormal,

src/distributions/common_interface.jl

@@ -168,5 +168,5 @@ function log_lik(psitilde::Vector{T}, y::Vector{T}, gas::Model{D, T},
         params = score_driven_recursion(gas, y; initial_params = initial_params)
     end
 
-    return log_likelihood(D, y, params, n)
+    return log_likelihood(D, y, params, n) / length(y)
 end

src/model.jl

@@ -1,7 +1,6 @@
 @testset "quantile and pearson residuals" begin
     n = 10^3
     lags = 1
-    normality_quantile_and_pearson_residuals(Beta, n, lags)
     normality_quantile_and_pearson_residuals(Normal, n, lags)
     normality_quantile_and_pearson_residuals(LogNormal, n, lags)
     normality_quantile_and_pearson_residuals(Gamma, n, lags)

test/test_diagnostics.jl

@@ -1,36 +1,36 @@
 @testset "Estimate" begin
-    @testset "Beta" begin
+    @testset "Normal" begin
         ω = [0.1, 0.1]
         A = [0.5 0; 0 0.5]
         B = [0.5 0; 0 0.5]
-        simulation = simulate_GAS_1_1(Beta, 0.0, ω, A, B, 1)
+        simulation = simulate_GAS_1_1(Normal, 0.0, ω, A, B, 1)
         @testset "Estimation by passing number of initial_points" begin
             # LBFGS
-            gas = Model(1, 1, Beta, 0.0)
+            gas = Model(1, 1, Normal, 0.0)
             fit!(gas, simulation; verbose = 2, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
             test_coefficients_GAS_1_1(gas, ω, A, B)
             # NelderMead
-            gas = Model(1, 1, Beta, 0.0)
+            gas = Model(1, 1, Normal, 0.0)
             fit!(gas, simulation; verbose = 2, opt_method = ScoreDrivenModels.NelderMead(gas, 3))
             test_coefficients_GAS_1_1(gas, ω, A, B)
             # IPNewton
-            gas = Model(1, 1, Beta, 0.0)
+            gas = Model(1, 1, Normal, 0.0)
             fit!(gas, simulation; verbose = 2, opt_method = ScoreDrivenModels.IPNewton(gas, 3))
             test_coefficients_GAS_1_1(gas, ω, A, B)
         end
         @testset "Estimation by passing initial_points" begin
             # LBFGS
-            gas = Model(1, 1, Beta, 0.0)
+            gas = Model(1, 1, Normal, 0.0)
             initial_points = [[0.1, 0.1, 0.5, 0.5, 0.5, 0.5]]
             fit!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.LBFGS(gas, initial_points))
             test_coefficients_GAS_1_1(gas, ω, A, B)
             # NelderMead
-            gas = Model(1, 1, Beta, 0.0)
+            gas = Model(1, 1, Normal, 0.0)
             initial_points = [[0.1, 0.1, 0.5, 0.5, 0.5, 0.5]]
             fit!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.NelderMead(gas, initial_points))
             test_coefficients_GAS_1_1(gas, ω, A, B)
             # IPNewton
-            gas = Model(1, 1, Beta, 0.0)
+            gas = Model(1, 1, Normal, 0.0)
             initial_points = [[0.1, 0.1, 0.5, 0.5, 0.5, 0.5]]
             fit!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.IPNewton(gas, initial_points))
             test_coefficients_GAS_1_1(gas, ω, A, B)