Several fixes (#103)

* also catch LapackException

* make sure startingvals are feasible for ARMA

* Don't return early from loglik! when not optimizing

* bump version number

* skip a doctest

* add a test

* typo

Project.toml

@@ -1,7 +1,7 @@
 name = "ARCHModels"
 uuid = "6d3278bc-c23a-5105-85e5-0d57d2bf684f"
 authors = ["Simon Broda <simon.broda@alumni.ethz.ch> and contributors."]
-version = "2.2"
+version = "2.2.1"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"

src/general.jl

@@ -44,7 +44,7 @@ function vcov(am::ARCHModel)
     Ji = try
         inv(J)
     catch e
-        if e isa LinearAlgebra.SingularException
+        if e in [LinearAlgebra.SingularException, LinearAlgebra.LAPACKException(1)]
             @warn "Fisher information is singular; vcov matrix is inaccurate."
             pinv(J)
         else

src/meanspecs.jl

@@ -171,14 +171,17 @@ function constraints(::Type{<:ARMA{p, q}}, ::Type{T})  where {T<:AbstractFloat,
     return lower, upper
 end
 
-function startingvals(::ARMA{p, q, T}, data::Vector{T})  where {p, q, T<:AbstractFloat}
+function startingvals(mod::ARMA{p, q, T}, data::Vector{T})  where {p, q, T<:AbstractFloat}
     N = length(data)
     X = Matrix{T}(undef, N-p, p+1)
     X[:, 1] .= T(1)
     for i = 1:p
         X[:, i+1] .= data[p-i+1:N-i]
     end
     phi = X \ data[p+1:end]
+    lower, upper = constraints(ARMA{p, q}, T)
+    phi[2:end] .= max.(phi[2:end], lower[2:p+1]*.99)
+    phi[2:end] .= min.(phi[2:end], upper[2:p+1]*.99)
     return T[phi..., zeros(T, q)...]
 end
 

src/univariatearchmodel.jl

@@ -287,21 +287,22 @@ end
 #dimensional array of the right type.
 @inline function loglik!(ht::AbstractVector{T2}, lht::AbstractVector{T2},
                          zt::AbstractVector{T2}, at::AbstractVector{T2}, vs::Type{VS}, ::Type{SD}, meanspec::MS,
-                         data::Vector{T1}, coefs::AbstractVector{T3}, subsetmask=trues(nparams(vs))
+                         data::Vector{T1}, coefs::AbstractVector{T3}, subsetmask=trues(nparams(vs)), returnearly=false
                          ) where {VS<:UnivariateVolatilitySpec, SD<:StandardizedDistribution,
                                   MS<:MeanSpec, T1<:AbstractFloat, T2, T3
                                   }
     garchcoefs, distcoefs, meancoefs = splitcoefs(coefs, VS, SD, meanspec)
 	lowergarch, uppergarch = constraints(VS, T1)
 	lowerdist, upperdist = constraints(SD, T1)
     lowermean, uppermean = constraints(MS, T1)
-    all(lowerdist.<distcoefs.<upperdist) && all(lowermean.<meancoefs.<uppermean) && all(lowergarch[subsetmask].<garchcoefs[subsetmask].<uppergarch[subsetmask]) || return T2(-Inf)
+	all_inbounds = all(lowerdist.<distcoefs.<upperdist) && all(lowermean.<meancoefs.<uppermean) && all(lowergarch[subsetmask].<garchcoefs[subsetmask].<uppergarch[subsetmask])
+    returnearly && !all_inbounds && return T2(-Inf)
 	garchcoefs .*= subsetmask
     T = length(data)
 	r1 = presample(VS)
 	r2 = presample(meanspec)
     r = max(r1, r2)
-    T > r || error("Sample too small.")
+    T - r > 0 || error("Sample too small.")
 	ki = kernelinvariants(SD, distcoefs)
     @inbounds begin
         h0 = var(data) # could be moved outside
@@ -329,10 +330,11 @@ end
         end#for
     end#inbounds
     LL += T*logconst(SD, distcoefs)
+	return all_inbounds ? LL : T2(-Inf)
 end#function
 
 function loglik(spec::Type{VS}, dist::Type{SD}, meanspec::MS,
-                   data::Vector{<:AbstractFloat}, coefs::AbstractVector{T2}, subsetmask=trues(nparams(spec))
+                   data::Vector{<:AbstractFloat}, coefs::AbstractVector{T2}, subsetmask=trues(nparams(spec)), returnearly=false
                    ) where {VS<:UnivariateVolatilitySpec, SD<:StandardizedDistribution,
                             MS<:MeanSpec, T2
                             }
@@ -342,7 +344,7 @@ function loglik(spec::Type{VS}, dist::Type{SD}, meanspec::MS,
     lht = CircularBuffer{T2}(r)
     zt = CircularBuffer{T2}(r)
 	at = CircularBuffer{T2}(r)
-    loglik!(ht, lht, zt, at, spec, dist, meanspec, data, coefs, subsetmask)
+    loglik!(ht, lht, zt, at, spec, dist, meanspec, data, coefs, subsetmask, returnearly)
 
 end
 
@@ -375,7 +377,7 @@ function _fit!(garchcoefs::Vector{T}, distcoefs::Vector{T},
               ) where {VS<:UnivariateVolatilitySpec, SD<:StandardizedDistribution,
                        MS<:MeanSpec, T<:AbstractFloat
                        }
-    obj = x -> -loglik(VS, SD, meanspec, data, x)
+    obj = x -> -loglik(VS, SD, meanspec, data, x, trues(length(garchcoefs)), true)
     coefs = vcat(garchcoefs, distcoefs, meancoefs)
     res = optimize(obj, coefs, algorithm; autodiff=autodiff, kwargs...)
     coefs .= Optim.minimizer(res)
@@ -458,7 +460,7 @@ function fitsubset(::Type{VS}, data::Vector{T}, maxlags::Int, subset::Tuple; dis
 	distcoefs = startingvals(SD, data)
     meancoefs = startingvals(ms, data)
 
-	obj = x -> -loglik(VS_large, SD, ms, data, x, mask)
+	obj = x -> -loglik(VS_large, SD, ms, data, x, mask, true)
     coefs = vcat(garchcoefs, distcoefs, meancoefs)
     res = optimize(obj, coefs, algorithm; autodiff=autodiff, kwargs...)
     coefs .= Optim.minimizer(res)
@@ -511,7 +513,7 @@ minimizes the [BIC](https://en.wikipedia.org/wiki/Bayesian_information_criterion
 - `algorithm=BFGS(), autodiff=:forward, kwargs...`: passed on to the optimizer.
 
 # Example
-```jldoctest
+```
 julia> selectmodel(EGARCH, BG96)
 
 EGARCH{1, 1, 2} model with Gaussian errors, T=1974.

test/runtests.jl

@@ -454,3 +454,9 @@ end
     str = sprint(io -> show(io, ccc.spec))
     @test startswith(str, "DCC{0, 0")
 end
+@testset "fixes" begin
+    X = [-49.78749999996362, 2951.7375000000347, 1496.437499999923, 973.8375, 2440.662500000128, 2578.062500000019, 1064.42500000032, 3378.0625000002415, -1971.5000000001048, 4373.899999999894]
+    am = fit(GARCH{2, 2}, X; meanspec = ARMA{2, 2});
+    @test length(volatilities(am)) == 10
+    @test loglikelihood(am) ≈ -86.01774
+end