Add EstVcov

Author: junyuan-chen

Project.toml

@@ -18,6 +18,7 @@ StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
+Vcov = "ec2bfdc2-55df-4fc9-b9ae-4958c2cf2486"
 
 [compat]
 Combinatorics = "1"
@@ -32,6 +33,7 @@ StatsFuns = "0.9"
 StatsModels = "0.6"
 Tables = "1"
 TypedTables = "1.2"
+Vcov = "0.4"
 julia = "1.3"
 
 [extras]

src/InteractionWeightedDIDs.jl

@@ -3,24 +3,27 @@ module InteractionWeightedDIDs
 using Base: Callable
 using Combinatorics: combinations
 using Dictionaries: Dictionary
-using FixedEffectModels: Vcov
 using FixedEffectModels: FixedEffectTerm, Combination,
     fe, has_fe, parse_fixedeffect, omitsintercept, hasintercept,
-    basecol, isnested, tss, Fstat
+    basecol, isnested, tss, Fstat, nunique
 using FixedEffects
-using LinearAlgebra: Symmetric, cholesky!
+using LinearAlgebra: Factorization, Symmetric, cholesky!
 using Printf
 using Reexport
 using SplitApplyCombine: group, groupfind, groupreduce
-using StatsBase: AbstractWeights, UnitWeights, NoQuote
+using StatsBase: AbstractWeights, UnitWeights, NoQuote, vcov
 using StatsModels: termvars, FullRank
 using Tables: columntable, getcolumn, rowtable
 using TypedTables: Table
+using Vcov
 @reexport using DiffinDiffsBase
 
 import Base: show
 import DiffinDiffsBase: required, default, transformed, combinedargs, _getsubcolumns
 
+# Handle naming conflicts
+const getvcov = vcov
+
 export Vcov,
        fe
 
@@ -30,6 +33,7 @@ export CheckVcov,
        MakeYXCols,
        MakeTreatCols,
        SolveLeastSquares,
+       EstVcov,
 
        RegressionBasedDID,
        Reg

src/procedures.jl

@@ -309,3 +309,49 @@ may be shared across multiple specifications.
 const SolveLeastSquares = StatsStep{:SolveLeastSquares, typeof(solveleastsquares)}
 
 required(::SolveLeastSquares) = (:tr, :yterm, :xterms, :yxcols, :treatcols)
+
+function _vcov(data, esample::BitArray,
+        vcov::Union{Vcov.SimpleCovariance,Vcov.RobustCovariance}, fes::Vector{FixedEffect})
+    dof_absorb = 0
+    for fe in fes
+        dof_absorb += nunique(fe)
+    end
+    return vcov, dof_absorb
+end
+
+function _vcov(data, esample::BitArray, vcov::Vcov.ClusterCovariance,
+        fes::Vector{FixedEffect})
+    cludata = _getsubcolumns(data, vcov.clusters, esample)
+    concrete_vcov = Vcov.materialize(cludata, vcov)
+    dof_absorb = 0
+    for fe in fes
+        any(c->isnested(fe, c.refs), concrete_vcov.clusters) && (dof_absorb += 1)
+    end
+    return concrete_vcov, dof_absorb
+end
+
+"""
+    estvcov(args...)
+
+Estimate variance-covariance matrix.
+See also [`EstVcov`](@ref).
+"""
+function estvcov(data, esample::BitArray, vcov::Vcov.CovarianceEstimator,
+        X::Matrix, crossx::Factorization, residuals::Vector, fes::Vector{FixedEffect})
+    concrete_vcov, dof_absorb = _vcov(data, esample, vcov, fes)
+    dof_residuals = max(1, sum(esample) - size(X,2) - dof_absorb)
+    vcov_data = Vcov.VcovData(X, crossx, residuals, dof_residuals)
+    vcov_mat = getvcov(vcov_data, concrete_vcov)
+    return (vcov_mat=vcov_mat, dof_residuals=dof_residuals), true
+end
+
+"""
+    EstVcov <: StatsStep
+
+Call [`InteractionWeightedDIDs.estvcov`](@ref) to estimate variance-covariance matrix.
+The returned object named `vcov_mat` and `dof_residuals`
+may be shared across multiple specifications.
+"""
+const EstVcov = StatsStep{:EstVcov, typeof(estvcov)}
+
+required(::EstVcov) = (:data, :esample, :vcov, :X, :crossx, :residuals, :fes)

test/procedures.jl

@@ -201,6 +201,9 @@ end
     nt = (tr=tr, yterm=term(:oop_spend), xterms=(term(1),), yxcols=yxcols0, treatcols=tcols0)
     ret, share = solveleastsquares(nt...)
     # Compare estimates with Stata
+    # gen col0 = wave_hosp==10 & wave==10
+    # gen col1 = wave_hosp==10 & wave==11
+    # reg oop_spend col0 col1
     @test ret.coef[1] ≈ 2862.4141 atol=1e-4
     @test ret.coef[2] ≈ 490.44869 atol=1e-4
     @test ret.coef[3] ≈ 3353.6565 atol=1e-4
@@ -224,3 +227,79 @@ end
     res = SolveLeastSquares()(nt)
     @test res.coef == ret.coef
 end
+
+@testset "EstVcov" begin
+    hrs = exampledata("hrs")
+    nobs = size(hrs, 1)
+    col0 = convert(Vector{Float64}, (hrs.wave_hosp.==10).&(hrs.wave.==10))
+    col1 = convert(Vector{Float64}, (hrs.wave_hosp.==10).&(hrs.wave.==11))
+    y = convert(Vector{Float64}, hrs.oop_spend)
+    X = hcat(col0, col1, ones(nobs, 1))
+    crossx = cholesky!(Symmetric(X'X))
+    coef = crossx \ (X'y)
+    residuals = y - X * coef
+    nt = (data=hrs, esample=trues(nobs), vcov=Vcov.simple(), X=X, crossx=crossx,
+        residuals=residuals, fes=FixedEffect[])
+    ret, share = estvcov(nt...)
+    # Compare estimates with Stata
+    # reg oop_spend col0 col1
+    # mat list e(V)
+    @test ret.vcov_mat[1,1] ≈ 388844.2 atol=0.1
+    @test ret.vcov_mat[2,2] ≈ 388844.2 atol=0.1
+    @test ret.vcov_mat[3,3] ≈ 20334.169 atol=1e-3
+    @test ret.vcov_mat[2,1] ≈ 20334.169 atol=1e-3
+    @test ret.vcov_mat[3,1] ≈ -20334.169 atol=1e-3
+    @test ret.dof_residuals == nobs - 3
+
+    nt = merge(nt, (vcov=Vcov.robust(), fes=FixedEffect[]))
+    ret, share = estvcov(nt...)
+    # Compare estimates with Stata
+    # reg oop_spend col0 col1, r
+    # mat list e(V)
+    @test ret.vcov_mat[1,1] ≈ 815817.44 atol=0.1
+    @test ret.vcov_mat[2,2] ≈ 254993.93 atol=1e-2
+    @test ret.vcov_mat[3,3] ≈ 19436.209 atol=1e-3
+    @test ret.vcov_mat[2,1] ≈ 19436.209 atol=1e-3
+    @test ret.vcov_mat[3,1] ≈ -19436.209 atol=1e-3
+    @test ret.dof_residuals == nobs - 3
+
+    nt = merge(nt, (vcov=Vcov.cluster(:hhidpn),))
+    ret, share = estvcov(nt...)
+    # Compare estimates with Stata
+    # reghdfe oop_spend col0 col1, noa clu(hhidpn)
+    # mat list e(V)
+    @test ret.vcov_mat[1,1] ≈ 744005.01 atol=0.1
+    @test ret.vcov_mat[2,2] ≈ 242011.45 atol=1e-2
+    @test ret.vcov_mat[3,3] ≈ 28067.783 atol=1e-3
+    @test ret.vcov_mat[2,1] ≈ 94113.386 atol=1e-2
+    @test ret.vcov_mat[3,1] ≈ 12640.559 atol=1e-2
+    @test ret.dof_residuals == nobs - 3
+
+    fes = FixedEffect[FixedEffect(hrs.hhidpn)]
+    wt = uweights(nobs)
+    feM = AbstractFixedEffectSolver{Float64}(fes, wt, Val{:cpu}, Threads.nthreads())
+    X = hcat(col0, col1)
+    _feresiduals!(Combination(y, X), feM, 1e-8, 10000)
+    crossx = cholesky!(Symmetric(X'X))
+    coef = crossx \ (X'y)
+    residuals = y - X * coef
+    nt = merge(nt, (X=X, crossx=crossx, residuals=residuals, fes=fes, vcov=Vcov.robust()))
+    ret, share = estvcov(nt...)
+    # Compare estimates with Stata
+    # reghdfe oop_spend col0 col1, a(hhidpn) vce(robust)
+    # mat list e(V)
+    @test ret.vcov_mat[1,1] ≈ 654959.97 atol=0.1
+    @test ret.vcov_mat[2,2] ≈ 503679.27 atol=0.1
+    @test ret.vcov_mat[2,1] ≈ 192866.2 atol=0.1
+    @test ret.dof_residuals == nobs - nunique(fes[1]) - 2
+
+    nt = merge(nt, (vcov=Vcov.cluster(:hhidpn),))
+    ret, share = estvcov(nt...)
+    # Compare estimates with Stata
+    # reghdfe oop_spend col0 col1, a(hhidpn) clu(hhidpn)
+    # mat list e(V)
+    @test ret.vcov_mat[1,1] ≈ 606384.66 atol=0.1
+    @test ret.vcov_mat[2,2] ≈ 404399.89 atol=0.1
+    @test ret.vcov_mat[2,1] ≈ 106497.43 atol=0.1
+    @test ret.dof_residuals == nobs - 3
+end

test/runtests.jl

@@ -4,9 +4,11 @@ using InteractionWeightedDIDs
 using DataFrames
 using Dictionaries
 using DiffinDiffsBase: required, default, transformed, combinedargs
+using FixedEffectModels: Combination, nunique
 using FixedEffects
 using InteractionWeightedDIDs: checkvcov!, checkfes!, makefesolver,
-    _feresiduals!, makeyxcols, maketreatcols, solveleastsquares
+    _feresiduals!, makeyxcols, maketreatcols, solveleastsquares, estvcov
+using LinearAlgebra
 using StatsBase: Weights, uweights
 
 import Base: ==