diff --git a/src/VARHAC.jl b/src/VARHAC.jl index 2c33a70..8be2b9e 100644 --- a/src/VARHAC.jl +++ b/src/VARHAC.jl @@ -1,14 +1,242 @@ -# VARHAC(; maxlag=12, lagstrategy=2, selectionstrategy=:aic) = VARHAC(maxlag, lagstrategy, selectionstrategy) +function VARHAC(; maxlag=12, lagstrategy=2, selectionstrategy=:aic) + @assert maxlag > 0 "maxlag must be greater than 0" + @assert lagstrategy in 1:3 "lagstrategy must be 1, 2, or 3" + @assert selectionstrategy in (:aic, :bic, :fixed) "selectionstrategy must be :aic, :bic, or :fixed" + imodel = selectionstrategy == :fixed ? 3 : selectionstrategy == :aic ? 1 : 2 + VARHAC(maxlag, lagstrategy, imodel) +end + +function avar(k::VARHAC, X::AbstractMatrix{T}; kwargs...) where T<:AbstractFloat + varhac(X, k.maxlag, k.ilag, k.lagstrategy) +end + +function avar(k::VARHAC, X::AbstractMatrix{T}; kwargs...) where T + varhac(float.(X), k.maxlag, k.ilag, k.lagstrategy) +end + +function varhac(dat, imax, ilag, imodel) + + ## VARHAC ESTIMATOR FOR THE COVARIANCE MATRIX (SPECTRAL DENSITY AT FREQUENCY ZERO) + + ## INPUT: + + ## dat: input matrix where each row corresponds to a different observation + ## imax: maximum lag order considered + ## 0 = no correction for serial correlation will be made + ## ilag: 1 = all elements enter with the same lag + ## 2 = own element can enter with a different lag + ## 3 = only the own lag enters + ## imodel: 1 = AIC is used + ## 2 = BIC is used + ## 3 = a fixed lag order equal to imax is used + + ## OUTPUT: + + ## ccc: VARHAC estimator + + nt = size(dat, 1) + kdim = size(dat, 2) + + ex = dat[imax:nt-1, :] + dat2 = dat[imax:nt-1, :] + ex1 = dat[imax:nt-1, :] + ex2 = dat[imax:nt-1, :] + dep = dat[imax+1:nt, 1] + + ddd = dat[imax+1:nt, :] + minres = ddd + aic = log.(sum(minres .^ 2, dims=1)) + minorder = zeros(kdim, 2) + + if imax > 0 + minpar = zeros(kdim, kdim * imax) + ## ALL ELEMENTS ENTER WITH THE SAME LAG (ILAG = 1) + if ilag == 1 + for k = 1:kdim + dep = ddd[:, k] + ex = dat[imax:nt-1, :] + + for iorder = 1:imax + if iorder == 1 + ex = dat[imax:nt-1, :] + else + ex = [ex dat[imax+1-iorder:nt-iorder, :]] + end + + if imodel != 3 + ## Run the VAR + b = (dep \ ex)' + resid = dep - ex[:, :] * b + ## Compute the model selection criterion + npar = kdim * iorder + if imodel == 1 + aicnew = log.(sum(resid .^ 2, dims=1)) .+ 2 * npar / (nt - imax) + else + aicnew = log.(sum(resid .^ 2, dims=1)) .+ log(nt - imax) * npar / (nt - imax) + end + if aicnew[1] < aic[k] + aic[k] = aicnew[1] + minpar[k, 1:kdim*iorder] = b' + minres[:, k] = resid + minorder[k] = iorder + end + end + + if imodel == 3 + b = (dep \ ex)' + resid = dep - ex * b + minpar = zeros(kdim, 2 * iorder) + minpar[k, :] = b' + minres[:, k] = resid + minorder[k] = imax + end + end + end + ## ONLY THE OWN LAG ENTERS (ILAG = 3 + elseif ilag == 3 + + for k = 1:kdim + dep = ddd[:, k] + + for iorder = 1:imax + if iorder == 1 + ex = dat[imax:nt-1, k] + else + ex = [ex dat[imax+1-iorder:nt-iorder, k]] + end + + if imodel != 3 + ## Run the VAR + b = ex \ dep + resid = dep - ex[:, :] * b + ## Compute the model selection criterion + npar = iorder + if imodel == 1 + aicnew = log.(sum(resid .^ 2, dims=1)) .+ 2 * npar / (nt - imax) + else + aicnew = log.(sum(resid .^ 2, dims=1)) .+ log(nt - imax) * npar / (nt - imax) + end + + if aicnew[1] < aic[k] + aic[k] = aicnew[1] + for i = 1:iorder + minpar[k, k+(i-1)*kdim] = b[i] + end + minres[:, k] = resid + minorder[k] = iorder + end + end + end + + if imodel == 3 + b = (dep \ ex)' + resid = dep - ex[:, :] * b + for i = 1:imax + minpar[k, k+(i-1)*kdim] = b[i] + end + minres[:, k] = resid + minorder[k] = imax + end + end + ## OWN ELEMENT CAN ENTER WITH A DIFFERENT LAG (ILAG = 2) + else + begin + for k = 1:kdim + dep = ddd[:, k] + if k == 1 + dat2 = dat[:, 2:kdim] + elseif k == kdim + dat2 = dat[:, 1:kdim-1] + else + dat2 = [dat[:, 1:k-1] dat[:, k+1:kdim]] + end + + for iorder = 0:imax + if iorder == 1 + ex1 = dat[imax:nt-1, k] + elseif iorder > 1 + ex1 = [ex1 dat[imax+1-iorder:nt-iorder, k]] + end + + for iorder2 = 0:imax + if iorder2 == 1 + ex2 = dat2[imax:nt-1, :] + elseif iorder2 > 1 + ex2 = [ex2 dat2[imax+1-iorder2:nt-iorder2, :]] + end + + if iorder + iorder2 == 0 + break + elseif iorder == 0 + ex = ex2 + elseif iorder2 == 0 + ex = ex1 + else + ex = [ex1 ex2] + end + ## Run the VAR + b = ex \ dep + #println(size(b)) + #println(size(ex)) + resid = dep .- ex[:, :] * b + ## Compute the model selection criterion + + npar = iorder + iorder2 * (kdim - 1) + + if imodel == 1 + aicnew = log.(sum(resid .^ 2, dims=1)) .+ 2 * npar / (nt - imax) + else + aicnew = log.(sum(resid .^ 2, dims=1)) .+ log(nt - imax) * npar / (nt - imax) + end + #println(aic[k]) + if aicnew[1] < aic[k] + aic[k] = aicnew[1] + minpar[k, :] = zeros(1, kdim * imax) + for i = 1:iorder + minpar[k, k+(i-1)*kdim] = b[i] + end + + for i = 1:iorder2 + if k == 1 + minpar[k, 2+(i-1)*kdim:kdim+(i-1)*kdim] = + b[iorder+1+(i-1)*(kdim-1):iorder+i*(kdim-1)]' + elseif k == kdim + minpar[k, 1+(i-1)*kdim:kdim-1+(i-1)*kdim] = + b[iorder+1+(i-1)*(kdim-1):iorder+i*(kdim-1)]' + else + minpar[k, 1+(i-1)*kdim:k-1+(i-1)*kdim] = + b[iorder+1+(i-1)*(kdim-1):iorder+k-1+(i-1)*(kdim-1)]' + minpar[k, k+1+(i-1)*kdim:kdim+(i-1)*kdim] = + b[iorder+k+(i-1)*(kdim-1):iorder+i*(kdim-1)]' + end + end + minres[:, k] = resid + minorder[k, :] = [iorder iorder2] + end + end + end + end + end + end + end + + ## COMPUTE THE VARHAC ESTIMATOR + covar = minres'minres #/ (nt - imax) + + bbb = diagm(0 => 1:kdim) + + if imax > 0 + for iorder = 1:imax + bbb = bbb - minpar[:, (iorder-1)*kdim+1:iorder*kdim] + end + end + + inv(bbb) * covar * inv(bbb') + +end -# function VARHAC(maxlag::Int, lagstrategy::Int, selectionstrategy::Symbol) -# @assert maxlag > 0 "maxlag must be greater than 0" -# @assert lagstrategy in 1:3 "lagstrategy must be 1, 2, or 3" -# @assert selectionstrategy in (:aic, :bic, :fixed) "selectionstrategy must be :aic, :bic, or :fixed" -# return VARHAC(maxlag, lagstrategy, selectionstrategy) -# end # function avar(k::VARHAC, X::Union{Matrix{F},Vector{F}}; kwargs...) where {F<:AbstractFloat} -# #avarhac(k, X; maxlag=k.maxlag, lagstrategy=k.lagstrategy, selectionstrategy=k.selectionstrategy) # if k.lagstrategy == 3 # D = fitvar(X, maxlag) # end @@ -126,228 +354,3 @@ # end - - - - - -# function varhac(dat, imax, ilag, imodel) - -# ## VARHAC ESTIMATOR FOR THE COVARIANCE MATRIX (SPECTRAL DENSITY AT FREQUENCY ZERO) - -# ## INPUT: - -# ## dat: input matrix where each row corresponds to a different observation -# ## imax: maximum lag order considered -# ## 0 = no correction for serial correlation will be made -# ## ilag: 1 = all elements enter with the same lag -# ## 2 = own element can enter with a different lag -# ## 3 = only the own lag enters -# ## imodel: 1 = AIC is used -# ## 2 = BIC is used -# ## 3 = a fixed lag order equal to imax is used - -# ## OUTPUT: - -# ## ccc: VARHAC estimator - -# nt = size(dat, 1) -# kdim = size(dat, 2) - -# ex = dat[imax:nt-1, :] -# dat2 = dat[imax:nt-1, :] -# ex1 = dat[imax:nt-1, :] -# ex2 = dat[imax:nt-1, :] -# dep = dat[imax+1:nt, 1] - -# ddd = dat[imax+1:nt, :] -# minres = ddd -# aic = log.(sum(minres .^ 2, dims=1)) -# minorder = zeros(kdim, 2) - -# if imax > 0 -# minpar = zeros(kdim, kdim * imax) -# ## ALL ELEMENTS ENTER WITH THE SAME LAG (ILAG = 1) -# if ilag == 1 -# for k = 1:kdim -# dep = ddd[:, k] -# ex = dat[imax:nt-1, :] - -# for iorder = 1:imax -# if iorder == 1 -# ex = dat[imax:nt-1, :] -# else -# ex = [ex dat[imax+1-iorder:nt-iorder, :]] -# end - -# if imodel != 3 -# ## Run the VAR -# b = (dep \ ex)' -# resid = dep - ex[:, :] * b -# ## Compute the model selection criterion -# npar = kdim * iorder -# if imodel == 1 -# aicnew = log.(sum(resid .^ 2, dims=1)) .+ 2 * npar / (nt - imax) -# else -# aicnew = log.(sum(resid .^ 2, dims=1)) .+ log(nt - imax) * npar / (nt - imax) -# end -# if aicnew[1] < aic[k] -# aic[k] = aicnew[1] -# minpar[k, 1:kdim*iorder] = b' -# minres[:, k] = resid -# minorder[k] = iorder -# end -# end - -# if imodel == 3 -# b = (dep \ ex)' -# resid = dep - ex * b -# minpar = zeros(kdim, 2 * iorder) -# minpar[k, :] = b' -# minres[:, k] = resid -# minorder[k] = imax -# end -# end -# end -# ## ONLY THE OWN LAG ENTERS (ILAG = 3 -# elseif ilag == 3 - -# for k = 1:kdim -# dep = ddd[:, k] - -# for iorder = 1:imax -# if iorder == 1 -# ex = dat[imax:nt-1, k] -# else -# ex = [ex dat[imax+1-iorder:nt-iorder, k]] -# end - -# if imodel != 3 -# ## Run the VAR -# b = ex \ dep -# resid = dep - ex[:, :] * b -# ## Compute the model selection criterion -# npar = iorder -# if imodel == 1 -# aicnew = log.(sum(resid .^ 2, dims=1)) .+ 2 * npar / (nt - imax) -# else -# aicnew = log.(sum(resid .^ 2, dims=1)) .+ log(nt - imax) * npar / (nt - imax) -# end - -# if aicnew[1] < aic[k] -# aic[k] = aicnew[1] -# for i = 1:iorder -# minpar[k, k+(i-1)*kdim] = b[i] -# end -# minres[:, k] = resid -# minorder[k] = iorder -# end -# end -# end - -# if imodel == 3 -# b = (dep \ ex)' -# resid = dep - ex[:, :] * b -# for i = 1:imax -# minpar[k, k+(i-1)*kdim] = b[i] -# end -# minres[:, k] = resid -# minorder[k] = imax -# end -# end -# ## OWN ELEMENT CAN ENTER WITH A DIFFERENT LAG (ILAG = 2) -# else -# begin -# for k = 1:kdim -# dep = ddd[:, k] -# if k == 1 -# dat2 = dat[:, 2:kdim] -# elseif k == kdim -# dat2 = dat[:, 1:kdim-1] -# else -# dat2 = [dat[:, 1:k-1] dat[:, k+1:kdim]] -# end - -# for iorder = 0:imax -# if iorder == 1 -# ex1 = dat[imax:nt-1, k] -# elseif iorder > 1 -# ex1 = [ex1 dat[imax+1-iorder:nt-iorder, k]] -# end - -# for iorder2 = 0:imax -# if iorder2 == 1 -# ex2 = dat2[imax:nt-1, :] -# elseif iorder2 > 1 -# ex2 = [ex2 dat2[imax+1-iorder2:nt-iorder2, :]] -# end - -# if iorder + iorder2 == 0 -# break -# elseif iorder == 0 -# ex = ex2 -# elseif iorder2 == 0 -# ex = ex1 -# else -# ex = [ex1 ex2] -# end -# ## Run the VAR -# b = ex \ dep -# #println(size(b)) -# #println(size(ex)) -# resid = dep .- ex[:, :] * b -# ## Compute the model selection criterion - -# npar = iorder + iorder2 * (kdim - 1) - -# if imodel == 1 -# aicnew = log.(sum(resid .^ 2, dims=1)) .+ 2 * npar / (nt - imax) -# else -# aicnew = log.(sum(resid .^ 2, dims=1)) .+ log(nt - imax) * npar / (nt - imax) -# end -# #println(aic[k]) -# if aicnew[1] < aic[k] -# aic[k] = aicnew[1] -# minpar[k, :] = zeros(1, kdim * imax) -# for i = 1:iorder -# minpar[k, k+(i-1)*kdim] = b[i] -# end - -# for i = 1:iorder2 -# if k == 1 -# minpar[k, 2+(i-1)*kdim:kdim+(i-1)*kdim] = -# b[iorder+1+(i-1)*(kdim-1):iorder+i*(kdim-1)]' -# elseif k == kdim -# minpar[k, 1+(i-1)*kdim:kdim-1+(i-1)*kdim] = -# b[iorder+1+(i-1)*(kdim-1):iorder+i*(kdim-1)]' -# else -# minpar[k, 1+(i-1)*kdim:k-1+(i-1)*kdim] = -# b[iorder+1+(i-1)*(kdim-1):iorder+k-1+(i-1)*(kdim-1)]' -# minpar[k, k+1+(i-1)*kdim:kdim+(i-1)*kdim] = -# b[iorder+k+(i-1)*(kdim-1):iorder+i*(kdim-1)]' -# end -# end -# minres[:, k] = resid -# minorder[k, :] = [iorder iorder2] -# end -# end -# end -# end -# end -# end -# end - -# ## COMPUTE THE VARHAC ESTIMATOR -# covar = minres'minres / (nt - imax) - -# bbb = diagm(0 => 1:kdim) - -# if imax > 0 -# for iorder = 1:imax -# bbb = bbb - minpar[:, (iorder-1)*kdim+1:iorder*kdim] -# end -# end - -# inv(bbb) * covar * inv(bbb') - -# end diff --git a/src/api.jl b/src/api.jl index 2fb78d6..4a0131c 100644 --- a/src/api.jl +++ b/src/api.jl @@ -1,8 +1,46 @@ invpseudohessian(x) = (t = typeof(x); error("Not defined for type $t", )) + +""" +### Description +The `momentmatrix` function returns the matrix of moment conditions for the estimation problem defined by `x`. Moment conditions are crucial for various estimation procedures, such as Generalized Method of Moments (GMM) and Maximum Likelihood Estimation (MLE). + +For MLE, the moment matrix corresponds to the inverse of the Hessian of the (pseudo-)log-likelihood function, evaluated at the data. For GMM, it represents the matrix of moment conditions evaluated at the observed data. + +This function can be extended to support user-defined types, allowing flexibility for different estimation methods. + +### Usage + +```julia +momentmatrix(x) -> Matrix +momentmatrix(x, t) -> Matrix +momentmatrix!(x, y) -> Matrix +``` + +- `momentmatrix(x)`: Returns the moment matrix for the estimation problem `x`. +- `momentmatrix(x, t)`: Returns the moment matrix for the estimation problem `x` when the parameter is equal to t. +- `momentmatrix!(x, y)`: In-place version, updating the matrix `y` with the moment conditions evaluated for `x`. + +The matrix returned is typically of size `(obs x m)`, where `obs` refers to the number of observations, and `m` refers to the number of moments. Users can define their own moment matrices for custom types by overloading this function. +""" momentmatrix(x) = (t = typeof(x); error("Not defined for type $t")) momentmatrix!(x, y) = (t = typeof(x); error("Not defined")) resid(x) = (t = typeof(x); error("Not defined for type $t")) +""" + bread(x) +Return the bread matrix for the estimation problem `x`. + +Note: This function is not defined for all types and must be extended for specific types. +""" bread(x) = (t = typeof(x); error("Not defined for type $t")) leverage(x) = (t = typeof(x); error("Not defined for type $t")) residualadjustment(k::AVarEstimator, x::Any) = (t = typeof(x); error("Not defined for type $t")) +function StatsBase.vcov(𝒦::AVarEstimator, e) + gᵢ= momentmatrix(e) + ## Bread mut return a k×m + B = bread(e) + Ω = aVar(𝒦, gᵢ) + B*Ω*B'/size(gᵢ,1) +end + +StatsBase.stderror(𝒦::AVarEstimator, e) = sqrt.(diag(vcov(𝒦, e))) diff --git a/src/types.jl b/src/types.jl index a9e1303..a84162c 100644 --- a/src/types.jl +++ b/src/types.jl @@ -7,7 +7,12 @@ Abstraction abstract type AVarEstimator end abstract type HAC{G} <: AVarEstimator end -abstract type VARHAC{G} <: AVarEstimator end + +struct VARHAC{T<:Int} <: AVarEstimator + maxlag::T + ilag::T + lagstrategy::T +end abstract type CrossSectionEstimator <: AVarEstimator end abstract type HR <: AVarEstimator end @@ -134,11 +139,7 @@ end const QuadraticSpectral = QuadraticSpectralKernel -# mutable struct VARHAC <: TimeSeries -# maxlag::Int -# lagstrategy::Int -# selectionstrategy::Symbol -# end + ## ------------------------------------------------------------------------------------------ # Define HAC constructor diff --git a/test/runtests.jl b/test/runtests.jl index bc678c5..08bf73d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -440,10 +440,6 @@ df.w = rand(StableRNG(123), size(df,1)) V2 = vcov(HC1(), lmm) idxr = isempty.(map(i->findall(all(isnan.(i))), eachrow(V1))) idxc = isempty.(map(i->findall(all(isnan.(i))), eachcol(V1))) - @show V1 - @show V2 - @show idxr - @show idxc end @testset "LM CR" begin @@ -681,3 +677,30 @@ end end ## TODO: Add more tests for the GLM interface for CR & DK + +@testset "Smoothed HAC" begin + X = randn(StableRNG(12322), 3700000, 3) + k = BartlettSmoother(3) + Σₛ = a𝕍ar(k, X; demean = true) + Σₕ = a𝕍ar(Parzen(3), X; demean = true) + @test Σₛ ≈ Σₕ rtol = 1e-3 + k = TruncatedSmoother(3) + Σₛ = a𝕍ar(k, X; demean = true) + Σₕ = a𝕍ar(Bartlett(3), X; demean = true) + @test Σₛ ≈ Σₕ rtol = 1e-3 +end + +@testset "Promotion" begin + X = randn(StableRNG(123), 100, 3) + Z = mapreduce(x->x.<=0, hcat, eachcol(X)) + aVar(HC0(), Z) + aVar(Bartlett(2), Z) + aVar(Bartlett{NeweyWest}(), Z) +end + +@testset "Stability" begin + X = randn(StableRNG(123), Float32, 100, 3) + @test eltype(aVar(HC0(), X)) == Float32 + @test eltype(aVar(Parzen(2), X)) == Float32 + @test eltype(aVar(Parzen{NeweyWest}(), X)) == Float32 +end