s-broda / ARCHModels.jl
Showing 3 of 4 files from the diff.

@@ -41,6 +41,23 @@
Loading
41 41
    end
42 42
end
43 43
44 +
mutable struct UnivariateSubsetARCHModel{T<:AbstractFloat,
45 +
                 				   VS<:UnivariateVolatilitySpec,
46 +
                 		  	  	   SD<:StandardizedDistribution{T},
47 +
                 				   MS<:MeanSpec{T},
48 +
								   N
49 +
                 				   } <: ARCHModel
50 +
	spec::VS
51 +
    data::Vector{T}
52 +
    dist::SD
53 +
    meanspec::MS
54 +
	fitted::Bool
55 +
	subset::NTuple{N, Int}
56 +
    function UnivariateSubsetARCHModel{T, VS, SD, MS, N}(spec, data, dist, meanspec, fitted, subset) where {T, VS, SD, MS, N}
57 +
        new(spec, data, dist, meanspec, fitted, subset)
58 +
    end
59 +
end
60 +
44 61
"""
45 62
    UnivariateARCHModel(spec::UnivariateVolatilitySpec, data::Vector; dist=StdNormal(),
46 63
	          			meanspec=NoIntercept(), fitted=false
@@ -75,6 +92,20 @@
Loading
75 92
    UnivariateARCHModel{T, VS, SD, MS}(spec, data, dist, meanspec, fitted)
76 93
end
77 94
95 +
function UnivariateSubsetARCHModel(spec::VS,
96 +
          		 			 data::Vector{T};
97 +
          					 dist::SD=StdNormal{T}(),
98 +
          				 	 meanspec::MS=NoIntercept{T}(),
99 +
		  			 		 fitted::Bool=false,
100 +
							 subset::NTuple{N, Int}
101 +
							 ) where {T<:AbstractFloat,
102 +
                    			 	  VS<:UnivariateVolatilitySpec,
103 +
                   					  SD<:StandardizedDistribution,
104 +
                   					  MS<:MeanSpec,
105 +
									  N
106 +
                   			 		  }
107 +
    UnivariateSubsetARCHModel{T, VS, SD, MS, N}(spec, data, dist, meanspec, fitted, subset)
108 +
end
78 109
loglikelihood(am::UnivariateARCHModel) = loglik(typeof(am.spec), typeof(am.dist),
79 110
                                      am.meanspec, am.data,
80 111
                                      vcat(am.spec.coefs, am.dist.coefs,
@@ -82,7 +113,17 @@
Loading
82 113
                                           )
83 114
                                      )
84 115
116 +
loglikelihood(am::UnivariateSubsetARCHModel) = loglik(typeof(am.spec), typeof(am.dist),
117 +
                                      am.meanspec, am.data,
118 +
                                      vcat(am.spec.coefs, am.dist.coefs,
119 +
                                           am.meanspec.coefs
120 +
                                           ),
121 +
									  subsetmask(typeof(am.spec), am.subset)
122 +
									  )
123 +
124 +
85 125
dof(am::UnivariateARCHModel) = nparams(typeof(am.spec)) + nparams(typeof(am.dist)) + nparams(typeof(am.meanspec))
126 +
dof(am::UnivariateSubsetARCHModel) = nparams(typeof(am.spec), am.subset) + nparams(typeof(am.dist)) + nparams(typeof(am.meanspec))
86 127
coef(am::UnivariateARCHModel)=vcat(am.spec.coefs, am.dist.coefs, am.meanspec.coefs)
87 128
coefnames(am::UnivariateARCHModel) = vcat(coefnames(typeof(am.spec)),
88 129
                                coefnames(typeof(am.dist)),
@@ -232,19 +273,17 @@
Loading
232 273
#but grows them by length(data); hence it should be called with an empty one-
233 274
#dimensional array of the right type.
234 275
@inline function loglik!(ht::AbstractVector{T2}, lht::AbstractVector{T2},
235 -
                         zt::AbstractVector{T2}, at::AbstractVector{T2}, ::Type{VS}, ::Type{SD}, meanspec::MS,
236 -
                         data::Vector{T1}, coefs::AbstractVector{T3}
276 +
                         zt::AbstractVector{T2}, at::AbstractVector{T2}, vs::Type{VS}, ::Type{SD}, meanspec::MS,
277 +
                         data::Vector{T1}, coefs::AbstractVector{T3}, subsetmask=trues(nparams(vs))
237 278
                         ) where {VS<:UnivariateVolatilitySpec, SD<:StandardizedDistribution,
238 279
                                  MS<:MeanSpec, T1<:AbstractFloat, T2, T3
239 280
                                  }
240 281
    garchcoefs, distcoefs, meancoefs = splitcoefs(coefs, VS, SD, meanspec)
241 -
    #the below 6 lines can be removed when using Fminbox
242 -
    lowergarch, uppergarch = constraints(VS, T1)
243 -
    lowerdist, upperdist = constraints(SD, T1)
282 +
	lowergarch, uppergarch = constraints(VS, T1)
283 +
	lowerdist, upperdist = constraints(SD, T1)
244 284
    lowermean, uppermean = constraints(MS, T1)
245 -
    lower = vcat(lowergarch, lowerdist, lowermean)
246 -
    upper = vcat(uppergarch, upperdist, uppermean)
247 -
    all(lower.<coefs.<upper) || return T2(-Inf)
285 +
    all(lowerdist.<distcoefs.<upperdist) && all(lowermean.<meancoefs.<uppermean) && all(lowergarch[subsetmask].<garchcoefs[subsetmask].<uppergarch[subsetmask]) || return T2(-Inf)
286 +
	garchcoefs .*= subsetmask
248 287
    T = length(data)
249 288
	r1 = presample(VS)
250 289
	r2 = presample(meanspec)
@@ -280,7 +319,7 @@
Loading
280 319
end#function
281 320
282 321
function loglik(spec::Type{VS}, dist::Type{SD}, meanspec::MS,
283 -
                   data::Vector{<:AbstractFloat}, coefs::AbstractVector{T2}
322 +
                   data::Vector{<:AbstractFloat}, coefs::AbstractVector{T2}, subsetmask=trues(nparams(spec))
284 323
                   ) where {VS<:UnivariateVolatilitySpec, SD<:StandardizedDistribution,
285 324
                            MS<:MeanSpec, T2
286 325
                            }
@@ -290,7 +329,7 @@
Loading
290 329
    lht = CircularBuffer{T2}(r)
291 330
    zt = CircularBuffer{T2}(r)
292 331
	at = CircularBuffer{T2}(r)
293 -
    loglik!(ht, lht, zt, at, spec, dist, meanspec, data, coefs)
332 +
    loglik!(ht, lht, zt, at, spec, dist, meanspec, data, coefs, subsetmask)
294 333
295 334
end
296 335
@@ -325,13 +364,6 @@
Loading
325 364
                       }
326 365
    obj = x -> -loglik(VS, SD, meanspec, data, x)
327 366
    coefs = vcat(garchcoefs, distcoefs, meancoefs)
328 -
    #for fminbox:
329 -
    # lowergarch, uppergarch = constraints(VS, T)
330 -
    # lowerdist, upperdist = constraints(SD, T)
331 -
    # lowermean, uppermean = constraints(MS, T)
332 -
    # lower = vcat(lowergarch, lowerdist, lowermean)
333 -
    # upper = vcat(uppergarch, upperdist, uppermean)
334 -
    # res = optimize(obj, lower, upper, coefs, Fminbox(algorithm); autodiff=autodiff, kwargs...)
335 367
    res = optimize(obj, coefs, algorithm; autodiff=autodiff, kwargs...)
336 368
    coefs .= Optim.minimizer(res)
337 369
    ng = nparams(VS)
@@ -396,6 +428,34 @@
Loading
396 428
	return UnivariateARCHModel(VS(coefs), data; dist=SD(distcoefs), meanspec=ms, fitted=true)
397 429
end
398 430
431 +
function fitsubset(::Type{VS}, data::Vector{T}, maxlags::Int, subset::Tuple; dist::Type{SD}=StdNormal{T},
432 +
             meanspec::Union{MS, Type{MS}}=Intercept{T}(T[0]), algorithm=BFGS(),
433 +
             autodiff=:forward, kwargs...
434 +
             ) where {VS<:UnivariateVolatilitySpec, SD<:StandardizedDistribution,
435 +
                      MS<:MeanSpec, T<:AbstractFloat
436 +
                      }
437 +
	#can't use dispatch for this b/c meanspec is a kwarg
438 +
	meanspec isa Type ? ms = meanspec(zeros(T, nparams(meanspec))) : ms = deepcopy(meanspec)
439 +
	VS_large = VS{ntuple(i->maxlags, length(subset))...}
440 +
	ng = nparams(VS_large)
441 +
	ns = nparams(SD)
442 +
	nm = nparams(typeof(ms))
443 +
	mask = subsetmask(VS_large, subset)
444 +
	garchcoefs = startingvals(VS_large, data, subset)
445 +
	distcoefs = startingvals(SD, data)
446 +
    meancoefs = startingvals(ms, data)
447 +
448 +
	obj = x -> -loglik(VS_large, SD, ms, data, x, mask)
449 +
    coefs = vcat(garchcoefs, distcoefs, meancoefs)
450 +
    res = optimize(obj, coefs, algorithm; autodiff=autodiff, kwargs...)
451 +
    coefs .= Optim.minimizer(res)
452 +
    garchcoefs .= coefs[1:ng]
453 +
	distcoefs .= coefs[ng+1:ng+ns]
454 +
    meancoefs .= coefs[ng+ns+1:ng+ns+nm]
455 +
	ms.coefs .= meancoefs
456 +
    return UnivariateSubsetARCHModel(VS_large(garchcoefs), data; dist=SD(distcoefs), meanspec=ms, fitted=true, subset=subset)
457 +
end
458 +
399 459
function fit!(am::UnivariateARCHModel; algorithm=BFGS(), autodiff=:forward, kwargs...)
400 460
    am.spec.coefs.=startingvals(typeof(am.spec), am.data)
401 461
    am.dist.coefs.=startingvals(typeof(am.dist), am.data)
@@ -475,14 +535,15 @@
Loading
475 535
	mylock=Threads.ReentrantLock()
476 536
    ndims = max(my_unwrap_unionall(VS)-1, 0) # e.g., two (p and q) for GARCH{p, q, T}
477 537
	ndims2 = max(my_unwrap_unionall(MS)-1, 0 )# e.g., two (p and q) for ARMA{p, q, T}
478 -
    res = Array{UnivariateARCHModel, ndims+ndims2}(undef, ntuple(i->maxlags - minlags + 1, ndims+ndims2))
538 +
    res = Array{UnivariateSubsetARCHModel, ndims+ndims2}(undef, ntuple(i->maxlags - minlags + 1, ndims+ndims2))
479 539
    Threads.@threads for ind in collect(CartesianIndices(size(res)))
480 -
		VSi = VS{ind.I[1:ndims] .+ minlags .-1...}
481 -
		MSi = (ndims2==0 ? meanspec : meanspec{ind.I[ndims+1:end] .+ minlags .- 1...})
482 -
		res[ind] = fit(VSi, data; dist=dist, meanspec=MSi,
540 +
		tup = (ind.I[1:ndims] .+ minlags .-1)
541 +
		MSi = (ndims2==0 ? deepcopy(meanspec) : meanspec{ind.I[ndims+1:end] .+ minlags .- 1...})
542 +
		res[ind] = fitsubset(VS, data, maxlags, tup; dist=dist, meanspec=MSi,
483 543
                       algorithm=algorithm, autodiff=autodiff, kwargs...)
484 544
        if show_trace
485 545
            lock(mylock)
546 +
			VSi = VS{tup...}
486 547
            Core.print(modname(VSi))
487 548
			ndims2>0 && Core.print("-", modname(MSi))
488 549
			Core.println(" model has ",
@@ -494,10 +555,9 @@
Loading
494 555
    end
495 556
    crits = criterion.(res)
496 557
    _, ind = findmin(crits)
497 -
    return res[ind]
558 +
	return fit(VS{res[ind].subset...}, data; dist=dist, meanspec=res[ind].meanspec, algorithm=algorithm, autodiff=autodiff, kwargs...)
498 559
end
499 560
500 -
501 561
function coeftable(am::UnivariateARCHModel)
502 562
    cc = coef(am)
503 563
    se = stderror(am)

@@ -71,9 +71,11 @@
Loading
71 71
const ARCH = GARCH{0}
72 72
73 73
@inline nparams(::Type{<:TGARCH{o, p, q}}) where {o, p, q} = o+p+q+1
74 +
@inline nparams(::Type{<:TGARCH{o, p, q}}, subset) where {o, p, q} = isempty(subset) ? 1 : sum(subset) + 1
74 75
75 76
@inline presample(::Type{<:TGARCH{o, p, q}}) where {o, p, q} = max(o, p, q)
76 77
78 +
77 79
Base.@propagate_inbounds @inline function update!(
78 80
        ht, lht, zt, at, ::Type{<:TGARCH{o, p, q}}, garchcoefs
79 81
        ) where {o, p, q}
@@ -112,6 +114,19 @@
Loading
112 114
    return x0
113 115
end
114 116
117 +
function startingvals(TT::Type{<:TGARCH}, data::Array{T} , subset::Tuple) where {T}
118 +
	o, p, q = subsettuple(TT, subsetmask(TT, subset)) # defend against (p, q) instead of (o, p, q)
119 +
	x0 = zeros(T, o+p+q+1)
120 +
    x0[2:o+1] .= 0.04/o
121 +
    x0[o+2:o+p+1] .= 0.9/p
122 +
    x0[o+p+2:end] .= o>0 ? 0.01/q : 0.05/q
123 +
    x0[1] = var(data)*(one(T)-sum(x0[2:o+1])/2-sum(x0[o+2:end]))
124 +
	mask = subsetmask(TT, subset)
125 +
	x0long = zeros(T, length(mask))
126 +
	x0long[mask] .= x0
127 +
    return x0long
128 +
end
129 +
115 130
function constraints(::Type{<:TGARCH{o,p,q}}, ::Type{T}) where {o,p, q, T}
116 131
    lower = zeros(T, o+p+q+1)
117 132
    upper = ones(T, o+p+q+1)
@@ -128,3 +143,36 @@
Loading
128 143
    names[o+p+2:o+p+q+1] .= (i -> "α"*subscript(i)).([1:q...])
129 144
    return names
130 145
end
146 +
147 +
@inline function subsetmask(VS_large::Union{Type{TGARCH{o, p, q}}, Type{TGARCH{o, p, q, T}}}, subs) where {o, p, q, T}
148 +
	ind = falses(nparams(VS_large))
149 +
	subset = zeros(Int, 3)
150 +
	subset[4-length(subs):end] .= subs
151 +
	ind[1] = true
152 +
	os = subset[1]
153 +
	ps = subset[2]
154 +
	qs = subset[3]
155 +
	@assert os <= o
156 +
	@assert ps <= p
157 +
	@assert qs <= q
158 +
	ind[2:2+os-1] .= true
159 +
	ind[2+o:2+o+ps-1] .= true
160 +
	ind[2+o+p:2+o+p+qs-1] .= true
161 +
	ind
162 +
end
163 +
164 +
@inline function subsettuple(VS_large::Union{Type{TGARCH{o, p, q}}, Type{TGARCH{o, p, q, T}}}, subsetmask) where {o, p, q, T}
165 +
	os = 0
166 +
	ps = 0
167 +
	qs = 0
168 +
	@inbounds @simd ivdep for i = 2 : o + 1
169 +
		os += subsetmask[i]
170 +
	end
171 +
	@inbounds @simd ivdep for i = o + 2 : o + p + 1
172 +
		ps += subsetmask[i]
173 +
	end
174 +
	@inbounds @simd ivdep for i = o + p + 2 : o + p + q + 1
175 +
		qs += subsetmask[i]
176 +
	end
177 +
	(os, ps, qs)
178 +
end

@@ -29,6 +29,7 @@
Loading
29 29
EGARCH{o, p, q}(coefs::Vector{T}) where {o, p, q, T}  = EGARCH{o, p, q, T}(coefs)
30 30
31 31
@inline nparams(::Type{<:EGARCH{o, p, q}}) where {o, p, q} = o+p+q+1
32 +
@inline nparams(::Type{<:EGARCH{o, p, q}}, subset) where {o, p, q} = isempty(subset) ? 1 : sum(subset) + 1
32 33
33 34
@inline presample(::Type{<:EGARCH{o, p, q}}) where {o, p, q} = max(o, p, q)
34 35
@@ -70,6 +71,19 @@
Loading
70 71
    return x0
71 72
end
72 73
74 +
function startingvals(TT::Type{<:EGARCH}, data::Array{T} , subset::Tuple) where {T}
75 +
	o, p, q = subsettuple(TT, subsetmask(TT, subset)) # defend against (p, q) instead of (o, p, q)
76 +
	x0 = zeros(T, o+p+q+1)
77 +
    x0[2:o+1] .= 0.04/o
78 +
    x0[o+2:o+p+1] .= 0.9/p
79 +
    x0[o+p+2:end] .= o>0 ? 0.01/q : 0.05/q
80 +
    x0[1] = var(data)*(one(T)-sum(x0[2:o+1])/2-sum(x0[o+2:end]))
81 +
	mask = subsetmask(TT, subset)
82 +
	x0long = zeros(T, length(mask))
83 +
	x0long[mask] .= x0
84 +
    return x0long
85 +
end
86 +
73 87
function constraints(::Type{<:EGARCH{o, p,q}}, ::Type{T}) where {o, p, q, T}
74 88
    lower = zeros(T, o+p+q+1)
75 89
    upper = zeros(T, o+p+q+1)
@@ -89,3 +103,36 @@
Loading
89 103
    names[o+p+2:o+p+q+1] .= (i -> "α"*subscript(i)).([1:q...])
90 104
    return names
91 105
end
106 +
107 +
@inline function subsetmask(VS_large::Union{Type{EGARCH{o, p, q}}, Type{EGARCH{o, p, q, T}}}, subs) where {o, p, q, T}
108 +
	ind = falses(nparams(VS_large))
109 +
	subset = zeros(Int, 3)
110 +
	subset[4-length(subs):end] .= subs
111 +
	ind[1] = true
112 +
	os = subset[1]
113 +
	ps = subset[2]
114 +
	qs = subset[3]
115 +
	@assert os <= o
116 +
	@assert ps <= p
117 +
	@assert qs <= q
118 +
	ind[2:2+os-1] .= true
119 +
	ind[2+o:2+o+ps-1] .= true
120 +
	ind[2+o+p:2+o+p+qs-1] .= true
121 +
	ind
122 +
end
123 +
124 +
@inline function subsettuple(VS_large::Union{Type{EGARCH{o, p, q}}, Type{EGARCH{o, p, q, T}}}, subsetmask) where {o, p, q, T}
125 +
	os = 0
126 +
	ps = 0
127 +
	qs = 0
128 +
	@inbounds @simd ivdep for i = 2 : o + 1
129 +
		os += subsetmask[i]
130 +
	end
131 +
	@inbounds @simd ivdep for i = o + 2 : o + p + 1
132 +
		ps += subsetmask[i]
133 +
	end
134 +
	@inbounds @simd ivdep for i = o + p + 2 : o + p + q + 1
135 +
		qs += subsetmask[i]
136 +
	end
137 +
	(os, ps, qs)
138 +
end
Files Coverage
src 99.41%
Project Totals (12 files) 99.41%
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading