Garch 1 1 (#47)

* Able to estimate GARCH(1,1)

* more opt_methods

* remove unecessary method

* update GARCH

Author: guilhermebodin

Project.toml

@@ -16,9 +16,10 @@ SpecialFunctions = "~0.8"
 julia = "~1"
 
 [extras]
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Random", "HypothesisTests"]
+test = ["DelimitedFiles", "Test", "Random", "HypothesisTests"]

docs/src/manual.md

@@ -1,8 +1,24 @@
 # Manual
 
+## Model Specification
+
+## Recursion
+
 ## Links
 
-Links are reparametrizations utilized in the estimation. For example, if we want to estimate a parameter ``f`` which is by definition strictly positive, then an obvious way to estimate ``f`` via numerical optimization is to model ``\tilde{f} = \ln{f}``. We refer to this procedure as **linking**. After obtaining the optimal value of ``\tilde{f}``, we can then **unlink** it to obtain ``f`` by computing ``f = e^{\tilde{f}}``.
+Links are reparametrizations utilized to ensure certain parameter is within its original domain, i.e. in a distribution one would like to ensure that the time varying parameter ``f \in \mathbb{R}^+``. The way to do this is to model ``\tilde{f} = \ln{f}``. More generally one can stablish that ``\tilde{f} = h(f)``. We refer to this procedure as **linking**. When the parameter is linked the GAS recursion happens in the domain of ``\tilde{f}`` and then one can recover the orginal parameter by ``f = \left(h\right)^-1(\tilde f)``. We refer to this procedure as **unlinking**. The new GAS recursion becomes.
+
+```math
+\begin{equation*}\left\{\begin{array}{ccl}
+    f_{t} &=& h^{-1}(\tilde f_t), \\
+    \tilde f_{t+1} &=& \omega + \sum_{i=1}^p A_{i}\tilde s_{t-i+1} + \sum_{j=1}^q B_{j}\tilde f_{t-j+1}
+    \end{array}
+    \right.
+\end{equation*}
+```
+
+Notice that the change in parametrization changes the dynamics of the model. The GAS(1,1) for a Normal distribution with inverse scaling ``d = 1`` is equivalent to the GARCH(1, 1) model, but only on the original parameter, if you work with a different parametrization the model is no longer equivalent.
+
 
 ### Types of links
 
@@ -11,6 +27,7 @@ The abstract type `Link` subsumes any type of link that can be expressed.
 ```@docs
 ScoreDrivenModels.IdentityLink
 ScoreDrivenModels.LogLink
+ScoreDrivenModels.LogitLink
 ```
 
 ### Link functions

examples/garch.jl

@@ -0,0 +1,65 @@
+using ScoreDrivenModels, Statistics, DelimitedFiles
+
+const SDM = ScoreDrivenModels
+
+# The GARCH(1, 1) is equivaent to a GAS(1, 1) Normal model with inverse scaling. The models are only equivalent under 
+# no reparametrizations, here this means that the you should fit the model overwriting the default link!, unlink! and
+# jacobian_link! methods for the Normal distribution to use only the IdentityLink.
+
+# Overwrite link interface to work with Identity links.
+function SDM.link!(param_tilde::Matrix{T}, ::Type{Normal}, param::Matrix{T}, t::Int) where T 
+    param_tilde[t, 1] = link(IdentityLink, param[t, 1])
+    param_tilde[t, 2] = link(IdentityLink, param[t, 2])
+    return
+end
+function SDM.unlink!(param::Matrix{T}, ::Type{Normal}, param_tilde::Matrix{T}, t::Int) where T 
+    param[t, 1] = unlink(IdentityLink, param_tilde[t, 1])
+    param[t, 2] = unlink(IdentityLink, param_tilde[t, 2])
+    return
+end
+function SDM.jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{Normal}, param::Matrix{T}, t::Int) where T 
+    aux.jac[1] = jacobian_link(IdentityLink, param[t, 1])
+    aux.jac[2] = jacobian_link(IdentityLink, param[t, 2])
+    return
+end
+
+# Data from [Bollerslev and Ghysels (JBES 1996)](https://doi.org/10.2307/1392425). Took from ARCHModels.jl
+y = readdlm("./test/data/BG96.csv")[:]
+
+# Evaluate some approximate initial params
+initial_params = [mean(y) var(y)]
+
+# There are 2 possibilities to estimate the model.
+# The first one is to understand that in the GARCH estimation process usually one 
+# defines that the parameters ω, α_0, α_1 and β_1 are constrained. The problem is 
+# that the equivalent GAS model relies that 
+# ω_1 = ω
+# ω_2 = α_0
+# A_1 = α_1
+# B_1 - A_1 = β_1
+# Here there is no easy way of adding a constraint B_1 - A_1 >= 0. An alternative is to define 
+# bounds assuring the upper bound of A_1 is smaller than the lower bound for B_1.
+
+# Define upper bounds and lower bounds for parameters ω_1, ω_2, A_1, B_1
+ub = [1.0; 1.0; 0.5; 1.0]
+lb = [-1.0; 0.0; 0.0; 0.5]
+
+# Define a GAS model where only \sigma^2 is varying over time
+gas = GAS(1, 1, Normal, 1.0, time_varying_params = [2])
+
+# Give an initial_point in the interior of the bounds.
+initial_point = [0.0; 0.5; 0.25; 0.75]
+
+# Estimate the model
+res = estimate!(gas, y; 
+                                  initial_params = initial_params, 
+                                  opt_method = IPNewton(gas, [initial_point]; ub = ub, lb = lb))
+
+
+# Another way to estimate the model is to rely on luck and give many random initial_points
+# for a non constrained optimization 
+
+gas = GAS(1, 1, Normal, 1.0, time_varying_params = [2])
+res = estimate!(gas, y; 
+                                  initial_params = initial_params, 
+                                  opt_method = NelderMead(gas, 100))

src/MLE.jl

@@ -7,13 +7,24 @@ struct EstimationResults{T <: AbstractFloat}
     aic::T
     bic::T
     llk::T
-    minimizer::Vector{T}
+    coefs::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))
+mutable struct AuxEstimation{T <: AbstractFloat}
+    psi::Vector{Vector{T}}
+    numerical_hessian::Vector{Matrix{T}}
+    loglikelihood::Vector{T}
+    opt_result::Vector{Optim.OptimizationResults}
+
+    function AuxEstimation{T}() where T
+        return new(
+            Vector{Vector{T}}(undef, 0), #psi
+            Vector{Matrix{T}}(undef, 0), 
+            Vector{T}(undef, 0), # loglikelihood
+            Vector{Optim.OptimizationResults}(undef, 0) # opt_result
+            )
+    end
 end
 
 function AIC(n_unknowns::Int, log_lik::T) where T
@@ -24,13 +35,23 @@ function BIC(n::Int, n_unknowns::Int, log_lik::T) where T
     return T(log(n) * n_unknowns - 2 * log_lik)
 end
 
+function update_aux_estimation!(aux_est::AuxEstimation{T}, func::Optim.TwiceDifferentiable,
+                                opt_result::Optim.OptimizationResults) where T
+                                
+    push!(aux_est.loglikelihood, -opt_result.minimum)
+    push!(aux_est.psi, opt_result.minimizer)
+    push!(aux_est.numerical_hessian, Optim.hessian!(func, opt_result.minimizer))
+    push!(aux_est.opt_result, opt_result)
+    return
+end
+
 function estimate!(sdm::SDM{D, T}, y::Vector{T};
                    initial_params::Matrix{T} = DEFAULT_INITIAL_PARAM,
                    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_seeds = length(opt_method.seeds)
     n = length(y)
 
     unknowns = find_unknowns(sdm)
@@ -40,49 +61,39 @@ function estimate!(sdm::SDM{D, T}, y::Vector{T};
     check_model_estimated(n_unknowns) && return sdm
 
     # optimize for each seed
-    psi = Vector{Vector{T}}(undef, 0)
-    numerical_hessian = Vector{Matrix{T}}(undef, 0)
-    loglikelihood = Vector{T}(undef, 0)
-    optseeds = Vector{Optim.OptimizationResults}(undef, 0)
+    aux_est = AuxEstimation{T}()
 
-    for i = 1:nseeds
+    for i = 1:n_seeds
         try 
             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)")
+            opt_result = optimize(func, opt_method, verbose, i)
+            update_aux_estimation!(aux_est, func, opt_result)
+            println("seed $i of $n_seeds - $(-opt_result.minimum)")
         catch err
             println(err)
             println("seed $i diverged")
         end
     end
 
-    if isempty(loglikelihood) 
+    if isempty(aux_est.loglikelihood) 
         println("No seed converged.")
         return
     end
 
-    best_llk, best_seed = findmax(loglikelihood)
-    num_hessian = numerical_hessian[best_seed]
-    best_psi = psi[best_seed]
+    best_llk, best_seed = findmax(aux_est.loglikelihood)
+    num_hessian = aux_est.numerical_hessian[best_seed]
+    coefs = aux_est.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_seed])
+        println(aux_est.opt_result[best_seed])
     end
 
     # return the estimated 
-    fill_psitilde!(sdm, best_psi, unknowns)
+    fill_psitilde!(sdm, coefs, unknowns)
 
     println("Finished!")
-    return EstimationResults{T}(aic, bic, best_llk, best_psi, num_hessian)
+    return EstimationResults{T}(aic, bic, best_llk, coefs, num_hessian)
 end

src/ScoreDrivenModels.jl

@@ -15,7 +15,12 @@ include("utils.jl")
 include("link_functions.jl")
 include("score.jl")
 include("MLE.jl")
-include("opt_methods.jl")
+
+# Optimization methods
+include("opt_methods/common_methods.jl")
+include("opt_methods/LBFGS.jl")
+include("opt_methods/IPNewton.jl")
+include("opt_methods/NelderMead.jl")
 
 # GAS
 include("gas/gas.jl")

src/abstracts.jl

@@ -1,6 +1,6 @@
 abstract type ScoreDrivenModel{D, T} end
 
-abstract type Unknowns_SDM end
+abstract type UnknownsSDM end
 
 """
     AbstractOptimizationMethod

src/distributions/common_interface.jl

@@ -8,6 +8,8 @@ const DISTS = [
     Weibull
 ]
 
+export Normal, Beta, Poisson, LogNormal, Gamma, Weibull
+
 function score!(score_til::Matrix{T}, y::T, ::Type{<:Distribution}, param::Matrix{T}, t::Int) where T
     return error("score! not implemented for $D distribution")
 end

src/gas/gas.jl

@@ -64,13 +64,13 @@ find_unknowns
 dim_unknowns
 length
 """
-mutable struct Unknowns_GAS <: Unknowns_SDM
+mutable struct UnknownsGAS <: UnknownsSDM
     ω::Vector{Int}
     A::Dict{Int, Vector{Int}}
     B::Dict{Int, Vector{Int}}
 end
 
-function fill_psitilde!(gas::GAS, psitilde::Vector{T}, unknowns::Unknowns_GAS) where T
+function fill_psitilde!(gas::GAS, psitilde::Vector{T}, unknowns::UnknownsGAS) where T
     offset = 0
     # fill ω
     for i in unknowns.ω
@@ -106,14 +106,14 @@ function find_unknowns(gas::GAS)
     for (k, v) in gas.B
         unknowns_B[k] = find_unknowns(v)
     end
-    return Unknowns_GAS(unknowns_ω, unknowns_A, unknowns_B)
+    return UnknownsGAS(unknowns_ω, unknowns_A, unknowns_B)
 end
 
 function dim_unknowns(gas::GAS)
     return length(find_unknowns(gas))
 end
 
-function length(unknowns::Unknowns_GAS)
+function length(unknowns::UnknownsGAS)
     len = length(values(unknowns.ω))
     for (k, v) in unknowns.A
         len += length(v)
@@ -125,7 +125,7 @@ function length(unknowns::Unknowns_GAS)
 end
 
 function log_lik(psitilde::Vector{T}, y::Vector{T}, gas::GAS{D, T}, 
-                 initial_params::Matrix{T}, unknowns::Unknowns_GAS, n::Int) where {D, T}
+                 initial_params::Matrix{T}, unknowns::UnknownsGAS, n::Int) where {D, T}
 
     # Use the unkowns vectors to fill the right positions
     fill_psitilde!(gas, psitilde, unknowns)

src/link_functions.jl

@@ -1,3 +1,6 @@
+export IdentityLink, LogLink, LogitLink
+export link, unlink, jacobian_link
+
 abstract type Link end
 
 """