Skip to content

Commit

Permalink
Merge pull request #65 from s-broda/subsets
Browse files Browse the repository at this point in the history
  • Loading branch information
@s-broda
s-broda committed Jun 29, 2020
2 parents a3ec3b4 + b111ad4 commit b7cd0dc
Show file tree
Hide file tree
Showing 4 changed files with 179 additions and 24 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
@@ -1,7 +1,7 @@
name = "ARCHModels"
uuid = "6d3278bc-c23a-5105-85e5-0d57d2bf684f"
authors = ["Simon Broda <simon.broda@uzh.ch>"]
version = "1.1.0"
version = "1.2.0"

[deps]
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Expand Down
47 changes: 47 additions & 0 deletions src/EGARCH.jl
View file
Expand Up @@ -29,6 +29,7 @@ Parameters: -0.1 0.1 0.9 0.04
EGARCH{o, p, q}(coefs::Vector{T}) where {o, p, q, T} = EGARCH{o, p, q, T}(coefs)

@inline nparams(::Type{<:EGARCH{o, p, q}}) where {o, p, q} = o+p+q+1
@inline nparams(::Type{<:EGARCH{o, p, q}}, subset) where {o, p, q} = isempty(subset) ? 1 : sum(subset) + 1

@inline presample(::Type{<:EGARCH{o, p, q}}) where {o, p, q} = max(o, p, q)

Expand Down Expand Up @@ -70,6 +71,19 @@ function startingvals(spec::Type{<:EGARCH{o, p, q}}, data::Array{T}) where {o, p
return x0
end

function startingvals(TT::Type{<:EGARCH}, data::Array{T} , subset::Tuple) where {T}
o, p, q = subsettuple(TT, subsetmask(TT, subset)) # defend against (p, q) instead of (o, p, q)
x0 = zeros(T, o+p+q+1)
x0[2:o+1] .= 0.04/o
x0[o+2:o+p+1] .= 0.9/p
x0[o+p+2:end] .= o>0 ? 0.01/q : 0.05/q
x0[1] = var(data)*(one(T)-sum(x0[2:o+1])/2-sum(x0[o+2:end]))
mask = subsetmask(TT, subset)
x0long = zeros(T, length(mask))
x0long[mask] .= x0
return x0long
end

function constraints(::Type{<:EGARCH{o, p,q}}, ::Type{T}) where {o, p, q, T}
lower = zeros(T, o+p+q+1)
upper = zeros(T, o+p+q+1)
Expand All @@ -89,3 +103,36 @@ function coefnames(::Type{<:EGARCH{o, p, q}}) where {o, p, q}
names[o+p+2:o+p+q+1] .= (i -> "α"*subscript(i)).([1:q...])
return names
end

@inline function subsetmask(VS_large::Union{Type{EGARCH{o, p, q}}, Type{EGARCH{o, p, q, T}}}, subs) where {o, p, q, T}
ind = falses(nparams(VS_large))
subset = zeros(Int, 3)
subset[4-length(subs):end] .= subs
ind[1] = true
os = subset[1]
ps = subset[2]
qs = subset[3]
@assert os <= o
@assert ps <= p
@assert qs <= q
ind[2:2+os-1] .= true
ind[2+o:2+o+ps-1] .= true
ind[2+o+p:2+o+p+qs-1] .= true
ind
end

@inline function subsettuple(VS_large::Union{Type{EGARCH{o, p, q}}, Type{EGARCH{o, p, q, T}}}, subsetmask) where {o, p, q, T}
os = 0
ps = 0
qs = 0
@inbounds @simd ivdep for i = 2 : o + 1
os += subsetmask[i]
end
@inbounds @simd ivdep for i = o + 2 : o + p + 1
ps += subsetmask[i]
end
@inbounds @simd ivdep for i = o + p + 2 : o + p + q + 1
qs += subsetmask[i]
end
(os, ps, qs)
end
48 changes: 48 additions & 0 deletions src/TGARCH.jl
View file
Expand Up @@ -71,9 +71,11 @@ Parameters: 1.0 0.3 0.4
const ARCH = GARCH{0}

@inline nparams(::Type{<:TGARCH{o, p, q}}) where {o, p, q} = o+p+q+1
@inline nparams(::Type{<:TGARCH{o, p, q}}, subset) where {o, p, q} = isempty(subset) ? 1 : sum(subset) + 1

@inline presample(::Type{<:TGARCH{o, p, q}}) where {o, p, q} = max(o, p, q)


Base.@propagate_inbounds @inline function update!(
ht, lht, zt, at, ::Type{<:TGARCH{o, p, q}}, garchcoefs
) where {o, p, q}
Expand Down Expand Up @@ -112,6 +114,19 @@ function startingvals(::Type{<:TGARCH{o,p,q}}, data::Array{T}) where {o, p, q, T
return x0
end

function startingvals(TT::Type{<:TGARCH}, data::Array{T} , subset::Tuple) where {T}
o, p, q = subsettuple(TT, subsetmask(TT, subset)) # defend against (p, q) instead of (o, p, q)
x0 = zeros(T, o+p+q+1)
x0[2:o+1] .= 0.04/o
x0[o+2:o+p+1] .= 0.9/p
x0[o+p+2:end] .= o>0 ? 0.01/q : 0.05/q
x0[1] = var(data)*(one(T)-sum(x0[2:o+1])/2-sum(x0[o+2:end]))
mask = subsetmask(TT, subset)
x0long = zeros(T, length(mask))
x0long[mask] .= x0
return x0long
end

function constraints(::Type{<:TGARCH{o,p,q}}, ::Type{T}) where {o,p, q, T}
lower = zeros(T, o+p+q+1)
upper = ones(T, o+p+q+1)
Expand All @@ -128,3 +143,36 @@ function coefnames(::Type{<:TGARCH{o,p,q}}) where {o,p, q}
names[o+p+2:o+p+q+1] .= (i -> "α"*subscript(i)).([1:q...])
return names
end

@inline function subsetmask(VS_large::Union{Type{TGARCH{o, p, q}}, Type{TGARCH{o, p, q, T}}}, subs) where {o, p, q, T}
ind = falses(nparams(VS_large))
subset = zeros(Int, 3)
subset[4-length(subs):end] .= subs
ind[1] = true
os = subset[1]
ps = subset[2]
qs = subset[3]
@assert os <= o
@assert ps <= p
@assert qs <= q
ind[2:2+os-1] .= true
ind[2+o:2+o+ps-1] .= true
ind[2+o+p:2+o+p+qs-1] .= true
ind
end

@inline function subsettuple(VS_large::Union{Type{TGARCH{o, p, q}}, Type{TGARCH{o, p, q, T}}}, subsetmask) where {o, p, q, T}
os = 0
ps = 0
qs = 0
@inbounds @simd ivdep for i = 2 : o + 1
os += subsetmask[i]
end
@inbounds @simd ivdep for i = o + 2 : o + p + 1
ps += subsetmask[i]
end
@inbounds @simd ivdep for i = o + p + 2 : o + p + q + 1
qs += subsetmask[i]
end
(os, ps, qs)
end
106 changes: 83 additions & 23 deletions src/univariatearchmodel.jl
View file
Expand Up @@ -41,6 +41,23 @@ mutable struct UnivariateARCHModel{T<:AbstractFloat,
end
end

mutable struct UnivariateSubsetARCHModel{T<:AbstractFloat,
VS<:UnivariateVolatilitySpec,
SD<:StandardizedDistribution{T},
MS<:MeanSpec{T},
N
} <: ARCHModel
spec::VS
data::Vector{T}
dist::SD
meanspec::MS
fitted::Bool
subset::NTuple{N, Int}
function UnivariateSubsetARCHModel{T, VS, SD, MS, N}(spec, data, dist, meanspec, fitted, subset) where {T, VS, SD, MS, N}
new(spec, data, dist, meanspec, fitted, subset)
end
end

"""
UnivariateARCHModel(spec::UnivariateVolatilitySpec, data::Vector; dist=StdNormal(),
meanspec=NoIntercept(), fitted=false
Expand Down Expand Up @@ -75,14 +92,38 @@ function UnivariateARCHModel(spec::VS,
UnivariateARCHModel{T, VS, SD, MS}(spec, data, dist, meanspec, fitted)
end

function UnivariateSubsetARCHModel(spec::VS,
data::Vector{T};
dist::SD=StdNormal{T}(),
meanspec::MS=NoIntercept{T}(),
fitted::Bool=false,
subset::NTuple{N, Int}
) where {T<:AbstractFloat,
VS<:UnivariateVolatilitySpec,
SD<:StandardizedDistribution,
MS<:MeanSpec,
N
}
UnivariateSubsetARCHModel{T, VS, SD, MS, N}(spec, data, dist, meanspec, fitted, subset)
end
loglikelihood(am::UnivariateARCHModel) = loglik(typeof(am.spec), typeof(am.dist),
am.meanspec, am.data,
vcat(am.spec.coefs, am.dist.coefs,
am.meanspec.coefs
)
)

loglikelihood(am::UnivariateSubsetARCHModel) = loglik(typeof(am.spec), typeof(am.dist),
am.meanspec, am.data,
vcat(am.spec.coefs, am.dist.coefs,
am.meanspec.coefs
),
subsetmask(typeof(am.spec), am.subset)
)


dof(am::UnivariateARCHModel) = nparams(typeof(am.spec)) + nparams(typeof(am.dist)) + nparams(typeof(am.meanspec))
dof(am::UnivariateSubsetARCHModel) = nparams(typeof(am.spec), am.subset) + nparams(typeof(am.dist)) + nparams(typeof(am.meanspec))
coef(am::UnivariateARCHModel)=vcat(am.spec.coefs, am.dist.coefs, am.meanspec.coefs)
coefnames(am::UnivariateARCHModel) = vcat(coefnames(typeof(am.spec)),
coefnames(typeof(am.dist)),
Expand Down Expand Up @@ -232,19 +273,17 @@ end
#but grows them by length(data); hence it should be called with an empty one-
#dimensional array of the right type.
@inline function loglik!(ht::AbstractVector{T2}, lht::AbstractVector{T2},
zt::AbstractVector{T2}, at::AbstractVector{T2}, ::Type{VS}, ::Type{SD}, meanspec::MS,
data::Vector{T1}, coefs::AbstractVector{T3}
zt::AbstractVector{T2}, at::AbstractVector{T2}, vs::Type{VS}, ::Type{SD}, meanspec::MS,
data::Vector{T1}, coefs::AbstractVector{T3}, subsetmask=trues(nparams(vs))
) where {VS<:UnivariateVolatilitySpec, SD<:StandardizedDistribution,
MS<:MeanSpec, T1<:AbstractFloat, T2, T3
}
garchcoefs, distcoefs, meancoefs = splitcoefs(coefs, VS, SD, meanspec)
#the below 6 lines can be removed when using Fminbox
lowergarch, uppergarch = constraints(VS, T1)
lowerdist, upperdist = constraints(SD, T1)
lowergarch, uppergarch = constraints(VS, T1)
lowerdist, upperdist = constraints(SD, T1)
lowermean, uppermean = constraints(MS, T1)
lower = vcat(lowergarch, lowerdist, lowermean)
upper = vcat(uppergarch, upperdist, uppermean)
all(lower.<coefs.<upper) || return T2(-Inf)
all(lowerdist.<distcoefs.<upperdist) && all(lowermean.<meancoefs.<uppermean) && all(lowergarch[subsetmask].<garchcoefs[subsetmask].<uppergarch[subsetmask]) || return T2(-Inf)
garchcoefs .*= subsetmask
T = length(data)
r1 = presample(VS)
r2 = presample(meanspec)
Expand Down Expand Up @@ -280,7 +319,7 @@ end
end#function

function loglik(spec::Type{VS}, dist::Type{SD}, meanspec::MS,
data::Vector{<:AbstractFloat}, coefs::AbstractVector{T2}
data::Vector{<:AbstractFloat}, coefs::AbstractVector{T2}, subsetmask=trues(nparams(spec))
) where {VS<:UnivariateVolatilitySpec, SD<:StandardizedDistribution,
MS<:MeanSpec, T2
}
Expand All @@ -290,7 +329,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)
loglik!(ht, lht, zt, at, spec, dist, meanspec, data, coefs, subsetmask)

end

Expand Down Expand Up @@ -325,13 +364,6 @@ function _fit!(garchcoefs::Vector{T}, distcoefs::Vector{T},
}
obj = x -> -loglik(VS, SD, meanspec, data, x)
coefs = vcat(garchcoefs, distcoefs, meancoefs)
#for fminbox:
# lowergarch, uppergarch = constraints(VS, T)
# lowerdist, upperdist = constraints(SD, T)
# lowermean, uppermean = constraints(MS, T)
# lower = vcat(lowergarch, lowerdist, lowermean)
# upper = vcat(uppergarch, upperdist, uppermean)
# res = optimize(obj, lower, upper, coefs, Fminbox(algorithm); autodiff=autodiff, kwargs...)
res = optimize(obj, coefs, algorithm; autodiff=autodiff, kwargs...)
coefs .= Optim.minimizer(res)
ng = nparams(VS)
Expand Down Expand Up @@ -396,6 +428,34 @@ function fit(::Type{VS}, data::Vector{T}; dist::Type{SD}=StdNormal{T},
return UnivariateARCHModel(VS(coefs), data; dist=SD(distcoefs), meanspec=ms, fitted=true)
end

function fitsubset(::Type{VS}, data::Vector{T}, maxlags::Int, subset::Tuple; dist::Type{SD}=StdNormal{T},
meanspec::Union{MS, Type{MS}}=Intercept{T}(T[0]), algorithm=BFGS(),
autodiff=:forward, kwargs...
) where {VS<:UnivariateVolatilitySpec, SD<:StandardizedDistribution,
MS<:MeanSpec, T<:AbstractFloat
}
#can't use dispatch for this b/c meanspec is a kwarg
meanspec isa Type ? ms = meanspec(zeros(T, nparams(meanspec))) : ms = deepcopy(meanspec)
VS_large = VS{ntuple(i->maxlags, length(subset))...}
ng = nparams(VS_large)
ns = nparams(SD)
nm = nparams(typeof(ms))
mask = subsetmask(VS_large, subset)
garchcoefs = startingvals(VS_large, data, subset)
distcoefs = startingvals(SD, data)
meancoefs = startingvals(ms, data)

obj = x -> -loglik(VS_large, SD, ms, data, x, mask)
coefs = vcat(garchcoefs, distcoefs, meancoefs)
res = optimize(obj, coefs, algorithm; autodiff=autodiff, kwargs...)
coefs .= Optim.minimizer(res)
garchcoefs .= coefs[1:ng]
distcoefs .= coefs[ng+1:ng+ns]
meancoefs .= coefs[ng+ns+1:ng+ns+nm]
ms.coefs .= meancoefs
return UnivariateSubsetARCHModel(VS_large(garchcoefs), data; dist=SD(distcoefs), meanspec=ms, fitted=true, subset=subset)
end

function fit!(am::UnivariateARCHModel; algorithm=BFGS(), autodiff=:forward, kwargs...)
am.spec.coefs.=startingvals(typeof(am.spec), am.data)
am.dist.coefs.=startingvals(typeof(am.dist), am.data)
Expand Down Expand Up @@ -475,14 +535,15 @@ function selectmodel(::Type{VS}, data::Vector{T};
mylock=Threads.ReentrantLock()
ndims = max(my_unwrap_unionall(VS)-1, 0) # e.g., two (p and q) for GARCH{p, q, T}
ndims2 = max(my_unwrap_unionall(MS)-1, 0 )# e.g., two (p and q) for ARMA{p, q, T}
res = Array{UnivariateARCHModel, ndims+ndims2}(undef, ntuple(i->maxlags - minlags + 1, ndims+ndims2))
res = Array{UnivariateSubsetARCHModel, ndims+ndims2}(undef, ntuple(i->maxlags - minlags + 1, ndims+ndims2))
Threads.@threads for ind in collect(CartesianIndices(size(res)))
VSi = VS{ind.I[1:ndims] .+ minlags .-1...}
MSi = (ndims2==0 ? meanspec : meanspec{ind.I[ndims+1:end] .+ minlags .- 1...})
res[ind] = fit(VSi, data; dist=dist, meanspec=MSi,
tup = (ind.I[1:ndims] .+ minlags .-1)
MSi = (ndims2==0 ? deepcopy(meanspec) : meanspec{ind.I[ndims+1:end] .+ minlags .- 1...})
res[ind] = fitsubset(VS, data, maxlags, tup; dist=dist, meanspec=MSi,
algorithm=algorithm, autodiff=autodiff, kwargs...)
if show_trace
lock(mylock)
VSi = VS{tup...}
Core.print(modname(VSi))
ndims2>0 && Core.print("-", modname(MSi))
Core.println(" model has ",
Expand All @@ -494,10 +555,9 @@ function selectmodel(::Type{VS}, data::Vector{T};
end
crits = criterion.(res)
_, ind = findmin(crits)
return res[ind]
return fit(VS{res[ind].subset...}, data; dist=dist, meanspec=res[ind].meanspec, algorithm=algorithm, autodiff=autodiff, kwargs...)
end


function coeftable(am::UnivariateARCHModel)
cc = coef(am)
se = stderror(am)
Expand Down

2 comments on commit b7cd0dc

@s-broda
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/17176

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v1.2.0 -m "<description of version>" b7cd0dce6bb50bec957910ed83fe3c5f2118315f
git push origin v1.2.0

Please sign in to comment.