diff --git a/Project.toml b/Project.toml index 95c7f22..1f5ab05 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "WildBootTests" uuid = "65c2e505-86ba-4c19-93f1-95506c1443d5" authors = ["droodman "] -version = "0.9.1" +version = "0.9.2" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" diff --git a/src/WRE.jl b/src/WRE.jl index 64cb315..8a1a9c6 100644 --- a/src/WRE.jl +++ b/src/WRE.jl @@ -581,7 +581,8 @@ function HessianFixedkappa!(o::StrBootTest{T}, dest::AbstractMatrix{T}, is::Vect if _jk !iszero(κ) && - (dest[row,1] = dot(i>0 ? o.Repl.XZ[:,i] : o.Repl.Xy₁par, j>0 ? o.Repl.V[:,j] : o.Repl.invXXXy₁par)) + (dest[row,1] = dot(i>0 ? o.Repl.XZ[:,i] : o.Repl.Xy₁par, + j>0 ? o.Repl.V[ :,j] : o.Repl.invXXXy₁par)) !isone(κ) && (dest[row,1] = iszero(κ) ? o.Repl.YY[i+1,j+1] : κ * dest[row,1] + (1 - κ) * o.Repl.YY[i+1,j+1]) end @@ -770,7 +771,7 @@ function MakeWREStats!(o::StrBootTest{T}, w::Integer) where T if o.null o.numerWRE .= view(o.β̈sAs,1:1,:) .+ (o.Repl.Rt₁ - o.r) / o.Repl.RRpar else - o.numerWRE .= view(o.β̈sAs,1:1,:) .- view(o.DGP.β̈ ,:,1) + o.numerWRE .= view(o.β̈sAs,1:1,:) .- view(o.DGP.β̈ ,:,1) isone(w) && (o.numerWRE[1] = o.β̈sAs[1] + (o.Repl.Rt₁[1] - o.r[1]) / o.Repl.RRpar[1]) end @@ -817,7 +818,7 @@ function MakeWREStats!(o::StrBootTest{T}, w::Integer) where T !iszero(o.fuller) && (o.κWRE .-= o.fuller / (o._Nobs - o.kX)) o.As .= o.κWRE .* view(o.YPXY✻, 2:o.Repl.kZ+1, :, 2:o.Repl.kZ+1) .+ (1 .- o.κWRE) .* view(o.YY✻, 2:o.Repl.kZ+1, :, 2:o.Repl.kZ+1) - invNaN!(o.As) + invsym!(o.As) t✻!(view(o.β̈s,:,:,1:1), o.As, o.κWRE .* view(o.YPXY✻, 2:o.Repl.kZ+1, :, 1) .+ (1 .- o.κWRE) .* view(o.YY✻, 2:o.Repl.kZ+1, :, 1)) else HessianFixedkappa!(o, o.δnumer, collect(1:o.Repl.kZ), 0, o.κ, _jk) @@ -825,7 +826,7 @@ function MakeWREStats!(o::StrBootTest{T}, w::Integer) where T HessianFixedkappa!(o, view(o.As, 1:i, :, i), collect(1:i), i, o.κ, _jk) end symmetrize!(o.As) - invNaN!(o.As) + invsym!(o.As) t✻!(view(o.β̈s,:,:,1:1), o.As, view(o.δnumer,:,:,1:1)) end @@ -871,7 +872,7 @@ function MakeWREStats!(o::StrBootTest{T}, w::Integer) where T if o.sqrt @storeWtGrpResults!(o.dist, o.numerWRE ./ sqrtNaN.(dropdims(o.denomWRE; dims=3))) else - invNaN!(o.denomWRE) + invsym!(o.denomWRE) _numer = view(o.numerWRE,:,:,1:1) @storeWtGrpResults!(o.dist, dropdims(_numer'o.denomWRE*_numer; dims=3)) # hand-code for 2-dimensional? XXX allocations end diff --git a/src/WildBootTests.jl b/src/WildBootTests.jl index b4c7611..6d64e57 100644 --- a/src/WildBootTests.jl +++ b/src/WildBootTests.jl @@ -95,6 +95,6 @@ function UpdateBootstrapcDenom!(o::StrBootTest{T} where T) nothing end -# include("precompiler.jl") +include("precompiler.jl") end diff --git a/src/estimators.jl b/src/estimators.jl index 6049ca9..cb48b0d 100644 --- a/src/estimators.jl +++ b/src/estimators.jl @@ -64,7 +64,7 @@ function InitVarsOLS!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstrac end o.invH = (pinv(H)) - R₁AR₁ = iszero(nrows(o.R₁perp)) ? o.invH : (o.R₁perp * invNaN(o.R₁perp'H*o.R₁perp) * o.R₁perp') # for DGP regression + R₁AR₁ = iszero(nrows(o.R₁perp)) ? o.invH : (o.R₁perp * invsym(o.R₁perp'H*o.R₁perp) * o.R₁perp') # for DGP regression o.β̈₀ = R₁AR₁ * X₁y₁ o.∂β̈∂r = R₁AR₁ * H * o.R₁invR₁R₁ - o.R₁invR₁R₁ @@ -79,17 +79,17 @@ function InitVarsOLS!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstrac v = view(parent.X₁, S,:) o.invMjk[g] = - v * R₁AR₁ * v' o.invMjk[g][1:length(S)+1:length(S)^2] .+= one(T) # add I - o.invMjk[g] .= invNaN(o.invMjk[g]) + o.invMjk[g] .= invsym(o.invMjk[g]) end else !isdefined(o, :S✻XX) && (o.S✻XX = panelcross(parent.X₁, parent.X₁, parent.info✻)) _H = reshape(H, (parent.kX, 1, parent.kX)) .- o.S✻XX - _invH = iszero(nrows(o.R₁perp)) ? invNaN(_H) : o.R₁perp * invNaN(o.R₁perp' * _H * o.R₁perp) * o.R₁perp' + _invH = iszero(nrows(o.R₁perp)) ? invsym(_H) : o.R₁perp * invsym(o.R₁perp' * _H * o.R₁perp) * o.R₁perp' o.XinvHjk = [view(parent.X₁, parent.info✻[g],:) * view(_invH,:,g,:) for g ∈ 1:parent.N✻] end end - o.A = iszero(nrows(Rperp)) ? o.invH : (Rperp * invNaN(Rperp'H*Rperp) * Rperp') # for replication regression + o.A = iszero(nrows(Rperp)) ? o.invH : (Rperp * invsym(Rperp'H*Rperp) * Rperp') # for replication regression o.AR = o.A * parent.R' (parent.scorebs || parent.robust) && (o.XAR = parent.X₁ * o.AR) nothing @@ -123,8 +123,8 @@ function InitVarsARubin!(o::StrEstimator{T}, parent::StrBootTest{T}) where T end H = ([X₁X₁ X₂X₁' ; X₂X₁ X₂X₂]) - o.A = invNaN(H) - R₁AR₁ = iszero(nrows(o.R₁perp)) ? o.A : (o.R₁perp * invNaN(o.R₁perp'H*o.R₁perp) * o.R₁perp') + o.A = invsym(H) + R₁AR₁ = iszero(nrows(o.R₁perp)) ? o.A : (o.R₁perp * invsym(o.R₁perp'H*o.R₁perp) * o.R₁perp') o.β̈₀ = R₁AR₁ * [X₁y₁ ; X₂y₁] o.∂β̈∂r = R₁AR₁ * [X₁Y₂ ; X₂Y₂] @@ -145,12 +145,12 @@ function InitVarsARubin!(o::StrEstimator{T}, parent::StrBootTest{T}) where T v₂ * (@view negR₁AR₁[parent.kX₁+1:end, 1:parent.kX₁ ]) * v₁' + v₂ * (@view negR₁AR₁[parent.kX₁+1:end, parent.kX₁+1:end]) * v₂' o.invMjk[g][1:length(S)+1:length(S)^2] .+= one(T) # add I - o.invMjk[g] .= invNaN(o.invMjk[g]) + o.invMjk[g] .= invsym(o.invMjk[g]) end else !isdefined(o, :S✻XX) && (o.S✻XX = [[S✻X₁X₁ S✻X₂X₁'] ; [S✻X₂X₁ S✻X₂X₂]]) _H = reshape(H, (parent.kX, 1, parent.kX)) .- o.S✻XX - _invH = iszero(nrows(o.R₁perp)) ? invNaN(_H) : o.R₁perp * invNaN(o.R₁perp' * _H * o.R₁perp) * o.R₁perp' + _invH = iszero(nrows(o.R₁perp)) ? invsym(_H) : o.R₁perp * invsym(o.R₁perp' * _H * o.R₁perp) * o.R₁perp' o.XinvHjk = [(S = parent.info✻[g]; X₁₂B(view(parent.X₁, S,:), view(parent.X₂, S,:), view(_invH,:,g,:))) for g ∈ 1:parent.N✻] end end @@ -206,7 +206,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract ZperpZperp, o.invZperpZperp, _invZperpZperp = invsymcrossjk(o.Zperp, parent.info✻) - _invZperpZperpZperpX₁ = _invZperpZperp * _ZperpX₁ + _invZperpZperpZperpX₁ = _invZperpZperp * _ZperpX₁ _invZperpZperpZperpX₂ = _invZperpZperp * _ZperpX₂ _invZperpZperpZperpZ = _invZperpZperp * _ZperpZ _invZperpZperpZperpZR₁ = _invZperpZperp * _ZperpZR₁ @@ -244,7 +244,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract o.XY₂ⱼₖ = [_X₁Y₂ ; _X₂Y₂] o.XXⱼₖ = [_X₁X₁ ; _X₂X₁ ;;; _X₂X₁' ; _X₂X₂] - o.invXXⱼₖ = invNaN(o.XXⱼₖ) + o.invXXⱼₖ = invsym(o.XXⱼₖ) o.H_2SLSⱼₖ = o.XZⱼₖ'o.invXXⱼₖ * o.XZⱼₖ (!isone(o.κ) || o.liml) && (o.H_2SLSmZZⱼₖ = o.H_2SLSⱼₖ - o.ZZⱼₖ) @@ -276,7 +276,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract end ZperpZ = o.Zperp'o.Zpar o.restricted && (ZperpZR₁ = o.Zperp'o.ZR₁) - o.invZperpZperp = invNaN(o.Zperp'o.Zperp) + o.invZperpZperp = invsym(o.Zperp'o.Zperp) end if !o.isDGP && parent.WREnonARubin @@ -307,7 +307,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract X₂X₁ = o.X₂'o.X₁ o.XX = [o.X₁'o.X₁ X₂X₁' ; X₂X₁ o.X₂'o.X₂] - o.invXX = invNaN(o.XX) + o.invXX = invsym(o.XX) X₁Y₂ = o.X₁'o.Y₂ X₂Y₂ = o.X₂'o.Y₂ o.XY₂ = [X₁Y₂ ; X₂Y₂] @@ -362,7 +362,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract S✻⋂X₂Zperp = panelcross(parent.X₂, o.Zperp, parent.info✻⋂) if !(parent.jk && parent.WREnonARubin) - o.invZperpZperp = iszero(o.kZperp) ? Matrix{T}(undef,0,0) : invNaN(sumpanelcross(o.S✻⋂ZperpZperp)) + o.invZperpZperp = iszero(o.kZperp) ? Matrix{T}(undef,0,0) : invsym(sumpanelcross(o.S✻⋂ZperpZperp)) ZperpX₁ = sumpanelcross(S✻⋂X₁Zperp)' ZperpX₂ = sumpanelcross(S✻⋂X₂Zperp)' end @@ -385,7 +385,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract X₁X₁ .-= ZperpX₁'o.invZperpZperp * ZperpX₁ X₂X₂ .-= ZperpX₂'o.invZperpZperp * ZperpX₂ o.XX = ([X₁X₁ X₂X₁' ; X₂X₁ X₂X₂]) - o.invXX = invNaN(o.XX) + o.invXX = invsym(o.XX) o.kX = ncols(o.XX) o.kX₁ = ncols(X₂X₁) @@ -419,6 +419,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract o.Y₂y₁ .-= ZperpY₂'o.invZperpZperpZperpy₁ o.Y₂Y₂ .-= ZperpY₂'o.invZperpZperpZperpY₂ o.X₂y₁ .-= ZperpX₂'o.invZperpZperpZperpy₁ + o.X₁y₁ .-= ZperpX₁'o.invZperpZperpZperpy₁ o.y₁y₁ -= Zperpy₁'o.invZperpZperpZperpy₁ @@ -498,7 +499,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract o.Zy₁par = copy(o.Zy₁) o.y₁pary₁par = copy(o.y₁y₁) o.Xy₁par = [o.X₁y₁ ; o.X₂y₁] - ((parent.jk && !o.isDGP) || (!parent.jk && (parent.granular || parent.scorebs))) && + ((parent.jk && !o.isDGP) || (parent.granular || parent.scorebs)) && (o.y₁par = copy(o.y₁)) o.V = o.invXX * o.XZ # in 2SLS case, estimator is (V' XZ)^-1 * (V'Xy₁). Also used in k-class and liml robust VCV by Stata convention @@ -509,7 +510,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract if o.liml o.invHⱼₖ = Array{T,3}(undef, o.kZ, parent.N✻, o.kZ) else - o.invHⱼₖ = invNaN(isone(o.κ) ? o.H_2SLSⱼₖ : o.ZZⱼₖ + o.κ * o.H_2SLSmZZⱼₖ) + o.invHⱼₖ = invsym(isone(o.κ) ? o.H_2SLSⱼₖ : o.ZZⱼₖ + o.κ * o.H_2SLSmZZⱼₖ) end end @@ -554,10 +555,10 @@ end function MakeH!(o::StrEstimator{T}, parent::StrBootTest{T}, makeXAR::Bool=false) where T H = isone(o.κ) ? o.H_2SLS : o.ZZ + o.κ * o.H_2SLSmZZ - o.invH = invNaN(H) + o.invH = invsym(H) if makeXAR # for replication regression in score bootstrap of IV/GMM - o.A = ncols(o.Rperp)>0 ? (o.Rperp * invNaN(o.Rperp'H*o.Rperp) * o.Rperp') : o.invH + o.A = ncols(o.Rperp)>0 ? (o.Rperp * invsym(o.Rperp'H*o.Rperp) * o.Rperp') : o.invH o.AR = o.A * (o.Rpar'parent.R') o.XAR = X₁₂B(o.X₁, o.X₂, o.V * o.AR) end @@ -598,7 +599,7 @@ function EstimateIV!(o::StrEstimator{T}, parent::StrBootTest{T}, _jk::Bool, r₁ if o.liml o.YPXY = ([[o.invXXXy₁par'o.Xy₁par] o.ZXinvXXXy₁par' ; o.ZXinvXXXy₁par o.H_2SLS]) - o.κ = 1/(1 - real(eigvalsNaN(invNaN(o.YY) * o.YPXY)[1])) # like Fast & Wild (81), but more stable, at least in Mata + o.κ = 1/(1 - real(eigvalsNaN(invsym(o.YY) * o.YPXY)[1])) # like Fast & Wild (81), but more stable, at least in Mata !iszero(o.fuller) && (o.κ -= o.fuller / (parent._Nobs - parent.kX)) MakeH!(o, parent) end @@ -611,9 +612,9 @@ function EstimateIV!(o::StrEstimator{T}, parent::StrBootTest{T}, _jk::Bool, r₁ if o.liml o.YPXYⱼₖ .= [o.invXXXy₁parⱼₖ'o.Xy₁parⱼₖ ; o.ZXinvXXXy₁parⱼₖ ;;; o.ZXinvXXXy₁parⱼₖ' ; o.H_2SLSⱼₖ] - o.κⱼₖ .= reshape(one(T) ./ (one(T) .- getindex.(real.(eigvalsNaN.(each(invNaN(o.YYⱼₖ) * o.YPXYⱼₖ))), 1)), (1,:,1)) + o.κⱼₖ .= reshape(one(T) ./ (one(T) .- getindex.(real.(eigvalsNaN.(each(invsym(o.YYⱼₖ) * o.YPXYⱼₖ))), 1)), (1,:,1)) !iszero(o.fuller) && (o.κⱼₖ .-= reshape(o.fuller ./ (o.Nobsⱼₖ .- parent.kX)), (1,:,1)) - o.invHⱼₖ .= o.ZZⱼₖ .+ o.κⱼₖ .* o.H_2SLSmZZⱼₖ; invNaN!(o.invHⱼₖ) + o.invHⱼₖ .= o.ZZⱼₖ .+ o.κⱼₖ .* o.H_2SLSmZZⱼₖ; invsym!(o.invHⱼₖ) o.β̈ⱼₖ .= o.κⱼₖ .* (o.ZXinvXXXy₁parⱼₖ .- o.Zy₁parⱼₖ) .+ o.Zy₁parⱼₖ else if isone(o.κ) @@ -680,7 +681,7 @@ function MakeResidualsIV!(o::StrEstimator{T}, parent::StrBootTest{T}) where T Xu = o.Xy₁par - o.XZ * o.β̈ # after DGP regression, compute Y₂ residuals by regressing Y₂ on X while controlling for y₁ residuals, done through FWL negXuinvuu = Xu / -uu - o.Π̈ = invNaN(o.XX + negXuinvuu * Xu') * (negXuinvuu * (o.Y₂y₁par - o.ZY₂'o.β̈ )' + o.XY₂) + o.Π̈ = invsym(o.XX + negXuinvuu * Xu') * (negXuinvuu * (o.Y₂y₁par - o.ZY₂'o.β̈ )' + o.XY₂) o.γ̈ = o.Rpar * o.β̈ + o.t₁ - parent.Repl.t₁ o.Ü₂Ü₂ = o.Y₂Y₂ - (o.Π̈)'o.XY₂ - o.XY₂'o.Π̈ + (o.Π̈)'o.XX*o.Π̈ @@ -710,7 +711,7 @@ function MakeResidualsIV!(o::StrEstimator{T}, parent::StrBootTest{T}) where T Xu .= view(o.Xy₁parⱼₖ,:,g,:); t✻minus!(Xu, view(o.XZⱼₖ,:,g,:), view(o.β̈ⱼₖ,:,g,1)) # after DGP regression, compute Y₂ residuals by regressing Y₂ on X while controlling for y₁ residuals, done through FWL negXuinvuu .= Xu ./ -uu - Π̈ⱼₖ .= invNaN(view(o.XXⱼₖ,:,g,:) + negXuinvuu * Xu') * (negXuinvuu * (view(o.Y₂y₁parⱼₖ,:,g,:) - view(o.ZY₂ⱼₖ,:,g,:)'view(o.β̈ⱼₖ,:,g,1))' + view(o.XY₂ⱼₖ,:,g,:)) + Π̈ⱼₖ .= invsym(view(o.XXⱼₖ,:,g,:) + negXuinvuu * Xu') * (negXuinvuu * (view(o.Y₂y₁parⱼₖ,:,g,:) - view(o.ZY₂ⱼₖ,:,g,:)'view(o.β̈ⱼₖ,:,g,1))' + view(o.XY₂ⱼₖ,:,g,:)) X₁₂Bminus!(view(o.Ü₂,S,:), view(o.X₁ⱼₖ,S,:), view(o.X₂ⱼₖ,S,:), Π̈ⱼₖ) t✻plus!(view(o.u⃛₁,S), view(o.Ü₂,S,:), o.RparY * view(o.β̈ⱼₖ,:,g,1) + _t₁Y) end diff --git a/src/interface.jl b/src/interface.jl index 860cdb0..dd68c95 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -14,7 +14,7 @@ struct BootTestResult{T} b::Vector{T} V::Matrix{T} auxweights::Union{Nothing,Matrix{T}} - o::StrBootTest + # o::StrBootTest end "Return test statistic" @@ -155,7 +155,7 @@ function __wildboottest( o.p, padj, o.B, o.BFeas, o.N✻, o.dof, o.dof_r, plot, peak, ci, getdist(o, diststat), getb(o), getV(o), - getauxweights && reps>0 ? getv(o) : nothing , o) + getauxweights && reps>0 ? getv(o) : nothing #=, o=#) end vecconvert(T::DataType, X) = isa(X, Vector{T}) ? X : Vector{T}(reshape(X, size(X,1))) diff --git a/src/nonWRE.jl b/src/nonWRE.jl index f47315e..25d6154 100644 --- a/src/nonWRE.jl +++ b/src/nonWRE.jl @@ -289,7 +289,7 @@ function MakeNonWREStats!(o::StrBootTest{T}, w::Integer) where T @storeWtGrpResults!(o.dist, o.numerw ./ sqrtNaN.(o.denom[1,1])) isone(w) && (o.statDenom = hcat(o.denom[1,1][1])) # original-sample denominator - elseif o.dof==2 # hand-code 2D numer'invNaN(denom)*numer + elseif o.dof==2 # hand-code 2D numer'invsym(denom)*numer t1 = view(o.numerw,1,:)'; t2 = view(o.numerw,2,:)'; t12 = t1.*t2 @storeWtGrpResults!(o.dist, (t1.^2 .* o.denom[2,2] .- 2 .* t12 .* o.denom[2,1] .+ t2.^2 .* o.denom[1,1]) ./ (o.denom[1,1].*o.denom[2,2] .- o.denom[2,1].^2)) isone(w) && diff --git a/src/structs.jl b/src/structs.jl index 7eda680..3f108a6 100644 --- a/src/structs.jl +++ b/src/structs.jl @@ -54,22 +54,22 @@ end function *(X::AbstractMatrix{T}, Y::DesignerMatrix{T}) where T if Y.type==regular - view(Base.:*(X, Y.M),:,collect(1:size(Y.M,2))) # collect() makes the column slicer Vector{Int64}, for type consistency + Base.:*(X, Y.M) # collect() makes the column slicer Vector{Int64}, for type consistency elseif Y.type==identity - view(X,:,collect(1:size(X,2))) + X elseif length(Y.p)==Y.size[2] - view(X,:,Y.p[Y.q]) + X[:,Y.p[Y.q]] else dest = zeros(T, size(X,1), Y.size[2]) # Y is selection matrix with some cols all 0 if length(dest)>0 - @inbounds for j ∈ axes(Y.p,1) + @inbounds for j ∈ eachindex(axes(Y.p,1)) pⱼ, qⱼ = Y.p[j], Y.q[j] @tturbo warn_check_args=false for i ∈ indices((dest,X),1) - dest[i,qⱼ] = X[i,pⱼ] + dest[i,qⱼ] = X[i,pⱼ] end end end - view(dest,:,collect(1:size(dest,2))) + dest end end function *(X::AbstractArray{T,3}, Y::DesignerMatrix{T}) where T @@ -170,7 +170,7 @@ mutable struct StrEstimator{T<:AbstractFloat} ZparX::Matrix{T} invZperpZperpZperpX₁::Matrix{T}; invZperpZperpZperpX₂::Matrix{T}; invZperpZperpZperpy₁::Vector{T}; invZperpZperpZperpY₂::Matrix{T}; S✻UY₂::Matrix{T}; invZperpZperpZperpZpar::Matrix{T}; invZperpZperpZperpZR₁::Matrix{T} Ü₂Ü₂::Matrix{T}; γ̈X::Vector{T}; γ̈Y::Vector{T}; γ⃛::Vector{T}; Xȳ₁::Vector{T}; ȳ₁ȳ₁::T; XÜ₂::Matrix{T}; ȳ₁Ü₂::Matrix{T}; Ȳ₂::Matrix{T}; ȳ₁::Vector{T} - Xpar₁toZparX::DesignerMatrix{T}; X₁noFWL::DesignerProduct{T} + Xpar₁toZparX::DesignerMatrix{T}; X₁noFWL::Matrix{T} #=DesignerProduct{T}=# X₁ⱼₖ::Matrix{T}; X₂ⱼₖ::Matrix{T}; y₁ⱼₖ::Vector{T}; Y₂ⱼₖ::Matrix{T}; Zⱼₖ::Matrix{T}; ZR₁ⱼₖ::Matrix{T}; Y₂y₁ⱼₖ::Array{T,3}; X₂y₁ⱼₖ::Array{T,3}; X₁y₁ⱼₖ::Array{T,3}; Zy₁ⱼₖ::Array{T,3}; XZⱼₖ::Array{T,3}; ZZⱼₖ::Array{T,3}; ZY₂ⱼₖ::Array{T,3}; y₁y₁ⱼₖ::Array{T,3}; XY₂ⱼₖ::Array{T,3}; invXXⱼₖ::Array{T,3}; XXⱼₖ::Array{T,3}; X₁ZR₁ⱼₖ::Array{T,3}; X₂ZR₁ⱼₖ::Array{T,3}; ZZR₁ⱼₖ::Array{T,3}; twoZR₁y₁ⱼₖ::Array{T,3}; ZR₁ZR₁ⱼₖ::Array{T,3}; ZR₁Y₂ⱼₖ::Array{T,3} Y₂y₁parⱼₖ::Array{T,3}; Zy₁parⱼₖ::Array{T,3}; y₁pary₁parⱼₖ::Array{T,3}; Xy₁parⱼₖ::Array{T,3}; y₁parⱼₖ::Vector{T}; H_2SLSⱼₖ::Array{T,3}; H_2SLSmZZⱼₖ::Array{T,3}; invHⱼₖ::Array{T,3} diff --git a/src/utilities.jl b/src/utilities.jl index 158651f..dba6311 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -20,22 +20,7 @@ end # iszero(nrows(X)) && (return Symmetric(X)) # X, ipiv, info = LinearAlgebra.LAPACK.sytrf!('U', Matrix(X)) # iszero(info) && LinearAlgebra.LAPACK.sytri!('U', X, ipiv) -function invNaN(X) - iszero(length(X)) && (return X) - try - pinv(X) - catch _ - fill(eltype(X)(NaN), size(X)) - end -end -function invsym(X) - iszero(length(X)) && (return X) - try - pinv(Symmetric(X)) - catch _ - fill(eltype(X)(NaN), size(X)) - end -end +@inline invsym(X) = X \ I eigvalsNaN(X) = try @@ -845,16 +830,16 @@ function t✻minus!(dest::AbstractArray{T,3}, A::AbstractVecOrMat{T}, B::Abstrac end # in-place inverse of a set of symmetric matrices -function invNaN!(A::Array{T,3}) where T +function invsym!(A::Array{T,3}) where T @inbounds for g ∈ eachindex(axes(A,2)) - A[:,g,:] = invNaN(@view A[:,g,:]) + A[:,g,:] = invsym(@view A[:,g,:]) end nothing end -function invNaN(A::Array{T,3}) where T +function invsym(A::Array{T,3}) where T dest = similar(A) @inbounds for g ∈ eachindex(axes(A,2)) - dest[:,g,:] = invNaN(@view A[:,g,:]) + dest[:,g,:] = invsym(@view A[:,g,:]) end dest end @@ -873,19 +858,19 @@ function crossjk(A::VecOrMat{T}, B::Vector{T}, info::Vector{UnitRange{Int64}}) w (vec(sumt), t) end -# given data matrix X and cluster-defining info vector, compute X'X, invNaN(X'X), and delete-g invNaN(X'X)'s efficiently +# given data matrix X and cluster-defining info vector, compute X'X, invsym(X'X), and delete-g invsym(X'X)'s efficiently function invsymcrossjk(X::Matrix{T}, info::Vector{UnitRange{Int64}}) where T SXX = panelcross(X,X,info) XX = sumpanelcross(SXX) - invXX = invNaN(XX) + invXX = invsym(XX) for (g,S) ∈ enumerate(info) Xg = view(X,S,:) if size(Xg,1) > size(Xg,2) - SXX[:,g,:] = invNaN(XX - Xg'Xg) + SXX[:,g,:] = invsym(XX - Xg'Xg) else XginvXX = Xg * invXX tmp = XginvXX * Xg'; tmp -= I - SXX[:,g,:] = invXX; t✻minus!(view(SXX,:,g,:), XginvXX'invNaN(tmp), XginvXX) + SXX[:,g,:] = invXX; t✻minus!(view(SXX,:,g,:), XginvXX'invsym(tmp), XginvXX) end end (XX, invXX, SXX) diff --git a/test/unittests.log b/test/unittests.log index dbc1352..9700d93 100644 --- a/test/unittests.log +++ b/test/unittests.log @@ -33,8 +33,8 @@ p = 0.1040 boottest post_self=.04, weight(mammen) reps(9999) boottype(score) t(7) = 1.8479 -p = 0.1651 -CI = [0.03008 0.07798] +p = 0.0903 +CI = [-0.178 -0.126; -0.1211 -0.0938; -0.06376 0.2902] boottest post_self=.04, reps(9999) jk @@ -60,8 +60,8 @@ F(2, 7) = 0.3444 p = 0.7332 boottest (post_self=.05) (post=-.02) (selfemployed=-.15), reps(9999) weight(webb) -F(3, 7) = 0.8534 -p = 0.7525 +F(3, 7) = 0.7430 +p = 0.5166 regress hasinsurance selfemployed post post_self @@ -259,8 +259,8 @@ CI = [-0.01093 0.123] areg wage ttl_exp collgrad tenure if occupation<. & grade<., robust absorb(industry) boottest tenure t(1833) = 1.6602 -p = 0.0981 -CI = [-0.004887 0.07348] +p = 0.0971 +CI = [-0.004733 0.07348] areg wage ttl_exp collgrad tenure [aw=hours] if occupation<. & grade<., robust absorb(industry) @@ -284,14 +284,10 @@ CI = [-0.4649 2.001] ivreghdfe wage ttl_exp collgrad tenure (occupation = union married) [aw=hours] if grade<., liml cluster(industry) absorb(age) boottest tenure -t(11) = 0.6068 -p = 0.4935 -CI = [-0.03753 0.05658] - boottest tenure, jk -t(11) = 0.6068 -p = 0.5145 -CI = [-0.04454 0.04321] +t(11) = 0.6000 +p = 0.8428 +CI = [-0.1851 0.05978] boottest collgrad tenure F(2, 11) = 3.5295 @@ -317,7 +313,7 @@ areg n w k, absorb(ind) boottest k, cluster(id year) t(8) = 28.5012 p = 0.0000 -CI = [0.7824 0.9413] +CI = [0.7823 0.9412] areg n w k [aw=ys], absorb(ind)