Fix compatibility issues

.github/workflows/CI-latest.yml

@@ -23,12 +23,12 @@ jobs:
         arch:
           - 'x64'
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v3
+      - uses: actions/cache@v4
         env:
           cache-name: cache-artifacts
         with:

.github/workflows/CI-stable.yml

@@ -17,8 +17,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.3'
-          - '1.6'
+          - '1.9'
           - '1'
         os:
           - 'ubuntu-latest'
@@ -27,12 +26,12 @@ jobs:
         arch:
           - 'x64'
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v3
+      - uses: actions/cache@v4
         env:
           cache-name: cache-artifacts
         with:
@@ -45,6 +44,7 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v3
+      - uses: codecov/codecov-action@v4
         with:
           file: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}

Project.toml

@@ -9,7 +9,6 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FixedEffectModels = "9d5cd8c9-2029-5cab-9928-427838db53e3"
 GroupedArrays = "6407cd72-fade-4a84-8a1e-56e431fc1533"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
@@ -19,26 +18,23 @@ Vcov = "ec2bfdc2-55df-4fc9-b9ae-4958c2cf2486"
 [compat]
 BSplines = "0.3"
 CSV = "0.9, 0.10"
-CovarianceMatrices = "0.10"
 DataFrames = "1"
 FFTW = "1"
-FixedEffectModels = "1.6"
+FixedEffectModels = "1.11"
 GroupedArrays = "0.3"
-Requires = "1"
-ShiftedArrays = "1"
+ShiftedArrays = "1, 2"
 StatsAPI = "1.2"
-StatsBase = "0.33"
+StatsBase = "0.33, 0.34"
 StatsFuns = "0.9, 1"
 Tables = "1"
-Vcov = "0.6"
-julia = "1.3"
+Vcov = "0.8"
+julia = "1.6"
 
 [extras]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
-CovarianceMatrices = "60f91f6f-d783-54cb-84f9-544141854719"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 ShiftedArrays = "1277b4bf-5013-50f5-be3d-901d8477a67a"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["CSV", "CovarianceMatrices", "DataFrames", "ShiftedArrays", "Test"]
+test = ["CSV", "DataFrames", "ShiftedArrays", "Test"]

src/LocalProjections.jl

@@ -2,23 +2,21 @@ module LocalProjections
 
 using BSplines: BSplineBasis, basismatrix
 using FFTW: fft!
-using FixedEffectModels: AbstractFixedEffectSolver, Combination, FixedEffect,
+using FixedEffectModels: AbstractFixedEffectSolver, FixedEffect,
     solve_residuals!, isnested, nunique
 using GroupedArrays: GroupedArray
 using LinearAlgebra: I, Symmetric, cholesky!, svd!, ldiv!, inv!, mul!
-using Requires
 using StatsAPI: RegressionModel, StatisticalModel
 using StatsBase: CovarianceEstimator, CoefTable, TestStat, PValue,
     AbstractWeights, Weights, UnitWeights, uweights
 using StatsFuns: tdistccdf, tdistinvccdf, chisqccdf
 using Tables
-using Vcov: ClusterCovariance, VcovData, robust, cluster, names, nclusters,
+using Vcov: ClusterCovariance, S_hat, robust, cluster, names, nclusters,
     ranktest!, pinvertible
 
 import Base: ==, show, size, length, vec, view
 import StatsAPI: coef, vcov, stderror, confint, coeftable, modelmatrix, residuals, dof_residual
 import Tables: getcolumn
-import Vcov: S_hat, dof_tstat
 
 # Reexport objects from StatsAPI
 export coef, vcov, stderror, confint, coeftable, modelmatrix, residuals, dof_residual
@@ -75,17 +73,4 @@ include("slp.jl")
 include("vce.jl")
 include("irf.jl")
 
-function __init__()
-    @require CovarianceMatrices = "60f91f6f-d783-54cb-84f9-544141854719" begin
-        function vcov(m::OLS, vce::CovarianceMatrices.RobustVariance)
-            dof = size(m.X,1) - dof_residual(m)
-            return CovarianceMatrices.sandwich(vce, m.invXX, m.score, dof=dof)
-        end
-        function vcov(m::Ridge, vce::CovarianceMatrices.RobustVariance)
-            dof = size(m.C,1) - dof_residual(m)
-            return CovarianceMatrices.sandwich(vce, m.invCCP, m.score, dof=dof)
-        end
-    end
-end
-
 end # module

src/lp.jl

@@ -388,9 +388,14 @@ function _firststage(nendo, niv, ys, xs, ws, sts, fes, clus, pw, nlag::Int, horz
         nX = size(X,2) + nendo - niv
         # Only allow heteroskedasticity-robust or cluster-robust VCE for rank test
         vcerank = vce isa ClusterCovariance ? cluster(names(vce), (CLU...,)) : robust()
-        r_kp = ranktest!(Endores, Zres, Pip, vcerank, nX, doffe)
-        F_kp = r_kp / size(Zres, 2)
-        p_kp = chisqccdf(size(Zres, 2) - size(Endores, 2) + 1, r_kp)
+        try
+            # ranktest! sometimes gives LAPACK errors even though coef estimates seem correct
+            r_kp = ranktest!(Endores, Zres, Pip, vcerank, nX, doffe)
+            F_kp = r_kp / size(Zres, 2)
+            p_kp = chisqccdf(size(Zres, 2) - size(Endores, 2) + 1, r_kp)
+        catch
+            r_kp, F_kp, p_kp = NaN, NaN, NaN
+        end
     else
         F_kp, p_kp = nothing, nothing
     end
@@ -453,14 +458,18 @@ function _est(::LeastSquaresLP, data, xnames, ys, xs, ws, sts, fes, clus, pw,
     if !firststagebyhorz && !any(x->x isa Cum, xs)
         dt = LPData(ys, xs, ws, sts, fes, clus, pw, nlag, minhorz, subset, groups, TF)
     end
-    F_kps = firststagebyhorz ? Vector{Float64}(undef, nhorz) : nothing
-    p_kps = firststagebyhorz ? Vector{Float64}(undef, nhorz) : nothing
+    F_kps = firststagebyhorz && testweakiv ? Vector{Float64}(undef, nhorz) : nothing
+    p_kps = firststagebyhorz && testweakiv ? Vector{Float64}(undef, nhorz) : nothing
     for h in minhorz:minhorz+nhorz-1
         i = h - minhorz + 1
         # Handle cases where all data need to be regenerated for each horizon
         if firststagebyhorz
-            fitted, F_kps[i], p_kps[i] = _firststage(nendo, niv, yfs, xfs,
+            fitted, F_kpsi, p_kpsi = _firststage(nendo, niv, yfs, xfs,
                 ws, sts, fes, clus, pw, nlag, h, subset, groups, testweakiv, vce; TF=TF)
+            if testweakiv
+                F_kps[i] = F_kpsi
+                p_kps[i] = p_kpsi
+            end
             xs[ix_iv] .= fitted
             dt = LPData(ys, xs, ws, sts, fes, clus, pw, nlag, h, subset, groups, TF)
         elseif any(x->x isa Cum, xs)
@@ -555,8 +564,8 @@ function lp(estimator, data, ynames;
         F_kp, p_kp = F_kps, p_kps
     end
 
-    return LocalProjectionResult(B, V, T, estimator, er, vce, ynames, xnames, wnames,
-        stnames, fenames,
+    return LocalProjectionResult(B, V, T, estimator, er, vce,
+        ynames, xnames, wnames, stnames, fenames,
         Dict{VarName,Int}(n=>i for (i,n) in enumerate(ynames)),
         Dict{VarName,Int}(n=>i for (i,n) in enumerate(xnames)),
         Dict{VarName,Int}(n=>i for (i,n) in enumerate(wnames)),

src/utils.jl

@@ -277,13 +277,30 @@ function _group(col::AbstractVector)
     return inds
 end
 
+# A work around for the type restrictions in FixedEffects.jl
+function _solve_residuals!(Y::AbstractVecOrMat, X::AbstractMatrix,
+        feM::AbstractFixedEffectSolver; progress_bar=false, kwargs...)
+    iterations = Int[]
+    convergeds = Bool[]
+    for j in axes(Y, 2)
+        _, iteration, converged = solve_residuals!(view(Y,:,j), feM; kwargs...)
+        push!(iterations, iteration)
+        push!(convergeds, converged)
+    end
+    for j in axes(X, 2)
+        _, iteration, converged = solve_residuals!(view(X,:,j), feM; kwargs...)
+        push!(iterations, iteration)
+        push!(convergeds, converged)
+    end
+    return iterations, convergeds
+end
+
 # Residualize columns in Y and X with weights W for fixed effects FE
 function _feresiduals!(Y::AbstractVecOrMat, X::AbstractMatrix, FE::Vector{FixedEffect},
         W::AbstractWeights; nfethreads::Int=Threads.nthreads(),
         fetol::Real=1e-8, femaxiter::Int=10000)
     feM = AbstractFixedEffectSolver{Float64}(FE, W, Val{:cpu}, nfethreads)
-    M = Combination(Y, X)
-    _, iters, convs = solve_residuals!(M, feM;
+    iters, convs = _solve_residuals!(Y, X, feM;
         tol=fetol, maxiter=femaxiter, progress_bar=false)
     iter = maximum(iters)
     conv = all(convs)

test/lp.jl

@@ -230,13 +230,6 @@ end
     @test stderror(f1)[9] ≈ 0.05917848282306491 atol=1e-8
     @test stderror(f1)[17] ≈ 0.04441907325684153 atol=1e-8
 
-    # Try Newey-West standard errors from CovarianceMatrices.jl
-    r1 = lp(df, Cum(:y), xnames=Cum(:g), wnames=(:newsy, :y, :g), iv=Cum(:g)=>:newsy,
-        nlag=4, nhorz=17, addylag=false, firststagebyhorz=true, vce=NeweyWest94())
-    f1 = irf(r1, Cum(:y), Cum(:g))
-    @test stderror(f1)[9] ≈ 0.06207907472277425 atol=1e-8
-    @test stderror(f1)[17] ≈ 0.043723073565539526 atol=1e-8
-
     # Construct lags for comparing results with FixedEffectModels.jl
     # For rank test, heteroskedasticity-robust VCE is used even with vce=HARVCE(EWC())
     r1_0 = lp(df, Cum(:y), xnames=Cum(:g), wnames=(:newsy, :y, :g), iv=Cum(:g)=>:newsy,
@@ -312,8 +305,8 @@ end
     # wwii
     r2 = lp(df, Cum(:y), xnames=(Cum(:g,:rec), Cum(:g,:exp), :rec), wnames=(:newsy, :y, :g),
         iv=(Cum(:g,:rec), Cum(:g,:exp))=>(:recnewsy, :expnewsy, :recg, :expg),
-        states=(:rec, :exp), nlag=4, nhorz=16, minhorz=1, addylag=false, firststagebyhorz=true,
-        subset=df.wwii.==0)
+        states=(:rec, :exp), nlag=4, nhorz=16, minhorz=1, addylag=false,
+        firststagebyhorz=true, subset=df.wwii.==0)
     f2rec = irf(r2, Cum(:y), Cum(:g,:rec))
     @test coef(f2rec)[1] ≈ .345617820757 atol=1e-9
     @test coef(f2rec)[8] ≈ 1.350903038984 atol=1e-8
@@ -323,10 +316,16 @@ end
     @test coef(f2exp)[8] ≈ .2171636455168 atol=1e-9
     @test coef(f2exp)[16] ≈ .2140385336746 atol=1e-9
 
+    r21 = lp(df, Cum(:y), xnames=(Cum(:g,:rec), Cum(:g,:exp), :rec), wnames=(:newsy, :y, :g),
+        iv=(Cum(:g,:rec), Cum(:g,:exp))=>(:recnewsy, :expnewsy, :recg, :expg),
+        states=(:rec, :exp), nlag=4, nhorz=16, minhorz=1, addylag=false,
+        firststagebyhorz=true, subset=df.wwii.==0, testweakiv=false)
+    @test r21.B ≈ r2.B
+
     r3 = lp(df, Cum(:y), xnames=(Cum(:g,:rec), Cum(:g,:exp), :rec), wnames=(:newsy, :y, :g),
         iv=(Cum(:g,:rec), Cum(:g,:exp))=>(:recnewsy, :expnewsy, :recg, :expg),
-        states=(:rec, :exp), nlag=4, nhorz=16, minhorz=1, addylag=false, firststagebyhorz=true,
-        subset=df.wwii.==0, panelid=:gid, panelweight=:wt)
+        states=(:rec, :exp), nlag=4, nhorz=16, minhorz=1, addylag=false,
+        firststagebyhorz=true, subset=df.wwii.==0, panelid=:gid)#, panelweight=:wt)
     @test r3.B ≈ r2.B[(1:28).!=4,:,:]
 
     @test_logs (:warn, "panelweight is ignored when panelid is nothing")

test/runtests.jl

@@ -2,20 +2,17 @@ using Test
 using LocalProjections
 
 using CSV
-using CovarianceMatrices
 using DataFrames
 using FixedEffectModels
 using LinearAlgebra: I
 using LocalProjections: kron_fastl, kron_fastr, getscore, _geto,
     OLS, VarName, _makeYX, _firststage, _lp, _toname,
     Ridge, _basismatrix, _makeYSr, _makeP
-using ShiftedArrays
+using ShiftedArrays: lag
 using Tables: getcolumn
 
 exampledata(name::Union{Symbol,String}) = CSV.read(datafile(name), DataFrame)
 
-const NeweyWest94 = CovarianceMatrices.ParzenKernel{CovarianceMatrices.NeweyWest}
-
 const tests = [
     "utils",
     "data",

test/slp.jl

@@ -295,14 +295,6 @@ end
     @test stderror(f1)[9] ≈ 0.0858939264679245 atol=1e-8
     @test stderror(f1)[17] ≈ 0.1270379208490072 atol=1e-8
 
-    # Try Newey-West standard errors from CovarianceMatrices.jl
-    r1 = lp(est, df, Cum(:y), xnames=Cum(:g), wnames=(:newsy, :y, :g), iv=Cum(:g)=>:newsy,
-        nlag=4, nhorz=17, addylag=false, firststagebyhorz=true, subset=df.wwii.==0,
-        vce=NeweyWest94())
-    f1 = irf(r1, Cum(:y), Cum(:g))
-    @test stderror(f1)[9] ≈ 0.08694010697048188 atol=1e-8
-    @test stderror(f1)[17] ≈ 0.13627053333187472 atol=1e-8
-
     r2 = lp(est, df, Cum(:y), xnames=Cum(:g), wnames=(:newsy, :y, :g),
         iv=Cum(:g)=>(:newsy, :g), nlag=4, nhorz=16, minhorz=1, addylag=false,
         firststagebyhorz=true, subset=df.wwii.==0, vce=HARVCE(EWC()))