compute halfwidth for ARubin in way less likely to produce NaN

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.8.3"
 manifest_format = "2.0"
-project_hash = "303a6ed21ec90fd0f6fc9ac1b7a11781ebdc8c23"
+project_hash = "56b87f7396a7dbd977ff046777573c9c7211aceb"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra"]
@@ -56,6 +56,18 @@ git-tree-sha1 = "0c5f81f47bbbcf4aea7b2959135713459170798b"
 uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b"
 version = "0.1.5"
 
+[[deps.Blosc]]
+deps = ["Blosc_jll"]
+git-tree-sha1 = "310b77648d38c223d947ff3f50f511d08690b8d5"
+uuid = "a74b3585-a348-5f62-a45c-50e91977d574"
+version = "0.7.3"
+
+[[deps.Blosc_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Lz4_jll", "Pkg", "Zlib_jll", "Zstd_jll"]
+git-tree-sha1 = "91d6baa911283650df649d0aea7c28639273ae7b"
+uuid = "0b7ba130-8d10-5ba8-a3d6-c5182647fed9"
+version = "1.21.1+0"
+
 [[deps.CPUSummary]]
 deps = ["CpuId", "IfElse", "Static"]
 git-tree-sha1 = "a7157ab6bcda173f533db4c93fc8a27a48843757"
@@ -153,6 +165,12 @@ deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
 uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 version = "1.6.0"
 
+[[deps.FileIO]]
+deps = ["Pkg", "Requires", "UUIDs"]
+git-tree-sha1 = "7be5f99f7d15578798f338f5433b6c432ea8037b"
+uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+version = "1.16.0"
+
 [[deps.FileWatching]]
 uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 
@@ -168,6 +186,24 @@ git-tree-sha1 = "187198a4ed8ccd7b5d99c41b69c679269ea2b2d4"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
 version = "0.10.32"
 
+[[deps.H5Zblosc]]
+deps = ["Blosc", "HDF5"]
+git-tree-sha1 = "26b22c9039b039e29ec4f4f989946de722e87ab9"
+uuid = "c8ec2601-a99c-407f-b158-e79c03c2f5f7"
+version = "0.1.0"
+
+[[deps.HDF5]]
+deps = ["Compat", "HDF5_jll", "Libdl", "Mmap", "Random", "Requires", "UUIDs"]
+git-tree-sha1 = "b5df7c3cab3a00c33c2e09c6bd23982a75e2fbb2"
+uuid = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+version = "0.16.13"
+
+[[deps.HDF5_jll]]
+deps = ["Artifacts", "JLLWrappers", "LibCURL_jll", "Libdl", "OpenSSL_jll", "Pkg", "Zlib_jll"]
+git-tree-sha1 = "4cc2bb72df6ff40b055295fdef6d92955f9dede8"
+uuid = "0234f1f7-429e-5d53-9886-15a909be8d59"
+version = "1.12.2+2"
+
 [[deps.HostCPUFeatures]]
 deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"]
 git-tree-sha1 = "f64b890b2efa4de81520d2b0fbdc9aadb65bdf53"
@@ -194,6 +230,12 @@ git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
 uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
 version = "0.1.1"
 
+[[deps.JLD]]
+deps = ["Compat", "FileIO", "H5Zblosc", "HDF5", "Printf"]
+git-tree-sha1 = "ec6afa4fd3402e4dd5b15b3e5dd2f7dd52043ce8"
+uuid = "4138dd39-2aa7-5051-a626-17a0bb65d9c8"
+version = "0.13.3"
+
 [[deps.JLLWrappers]]
 deps = ["Preferences"]
 git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
@@ -247,6 +289,12 @@ git-tree-sha1 = "da5317a78e2a9f692730345cf3ff820109f406d3"
 uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
 version = "0.12.141"
 
+[[deps.Lz4_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "5d494bc6e85c4c9b626ee0cab05daa4085486ab1"
+uuid = "5ced341a-0733-55b8-9ab6-a4889d929147"
+version = "1.9.3+0"
+
 [[deps.MacroTools]]
 deps = ["Markdown", "Random"]
 git-tree-sha1 = "42324d08725e200c23d4dfb549e0d5d89dede2d2"
@@ -273,6 +321,9 @@ git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
 uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
 version = "1.0.2"
 
+[[deps.Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
 [[deps.MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
 version = "2022.2.1"
@@ -303,6 +354,12 @@ deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
 version = "0.8.1+0"
 
+[[deps.OpenSSL_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "f6e9dba33f9f2c44e08a020b0caf6903be540004"
+uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
+version = "1.1.19+0"
+
 [[deps.OpenSpecFun_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
@@ -360,6 +417,12 @@ git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
 uuid = "189a3867-3050-52da-a836-e630ba90ab69"
 version = "1.2.2"
 
+[[deps.Requires]]
+deps = ["UUIDs"]
+git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
+uuid = "ae029012-a4dd-5104-9daa-d747884805df"
+version = "1.3.0"
+
 [[deps.Rmath]]
 deps = ["Random", "Rmath_jll"]
 git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
@@ -512,6 +575,12 @@ deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
 version = "1.2.12+3"
 
+[[deps.Zstd_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "e45044cd873ded54b6a5bac0eb5c971392cf1927"
+uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
+version = "1.5.2+0"
+
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

Project.toml

@@ -15,8 +15,8 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 
 [compat]
 Distributions = "0.25"
+LoopVectorization = "0.12.141"
 SnoopPrecompile = "1.0.1"
 SortingAlgorithms = "1.0"
 StableRNGs = "1.0"
 julia = "1.8"
-LoopVectorization = "0.12.141"

src/WRE.jl

@@ -652,7 +652,7 @@ function MakeWREStats!(o::StrBootTest{T}, w::Integer) where T
 					o.YY✻_b[1:i+1,i+1]   = YY✻[i+1][:,b]  # fill uppper triangles, which is all that invsym() looks at
 					o.YPXY✻_b[1:i+1,i+1] = YPXY✻[i+1][:,b]
 				end
-				o.κ = 1/(1 - real(eigvals(invsym(o.YY✻_b) * Symmetric(o.YPXY✻_b))[1]))
+				o.κ = 1/(1 - real(eigvalsNaN(invsym(o.YY✻_b) * Symmetric(o.YPXY✻_b))[1]))
 				!iszero(o.fuller) && (o.κ -= o.fuller / (o._Nobs - o.kX))
 				β̈s[:,b] = (A[b] = invsym(o.κ*o.YPXY✻_b[2:end,2:end] + (1-o.κ)*o.YY✻_b[2:end,2:end])) * (o.κ*o.YPXY✻_b[1,2:end] + (1-o.κ)*o.YY✻_b[1,2:end])
 			end

src/estimators.jl

@@ -577,12 +577,12 @@ function EstimateIV!(o::StrEstimator{T}, parent::StrBootTest{T}, #=_jk::Bool,=#
 
   if o.isDGP
 	  if o.liml
-	    o.κ = 1/(1 - real(eigvals(invsym(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))
 
 			# if parent.jk
 			# 	for g ∈ 1:parent.N✻
-			# 		_κ[g] = 1/(1 - real(eigvals(invsym(_YY) * _YPXY)[1]))
+			# 		_κ[g] = 1/(1 - real(eigvalsNaN(invsym(_YY) * _YPXY)[1]))
 			# 	end
 			# 	!iszero(o.fuller) && (_κ .-= o.fuller ./ (o._Nobsjk .- parent.kX))
 			# end

src/plot-CI.jl

@@ -41,17 +41,17 @@ function plot!(o::StrBootTest{T}) where T
 
 	iszero(nrows(o.gridmin   )) && (o.gridmin    = fill(T(NaN), o.q))
 	iszero(nrows(o.gridmax   )) && (o.gridmax    = fill(T(NaN), o.q))
-	iszero(nrows(o.gridpoints)) && (o.gridpoints = fill(Float32(NaN), o.q))
+	iszero(nrows(o.gridpoints)) && (o.gridpoints = fill(T(NaN), o.q))
 
 	o.gridpoints[isnan.(o.gridpoints) .| isinf.(o.gridpoints)] .= 25
 
   o.boottest!(o)
 
 	Phi = quantile(Normal{T}(zero(T),one(T)), α/2)
   if o.arubin
-		p = o.dist[1] * o.multiplier
-		p = ccdf(Chisq{T}(T(o.dof)), o.sqrt ? p^2 : p)
-		halfwidth = abs.(o.confpeak) * quantile(Normal{T}(zero(T),one(T)), p/2) / Phi
+		# p = o.dist[1] * o.multiplier
+		# p = ccdf(Chisq{T}(T(o.dof)), o.sqrt ? p^2 : p)
+		halfwidth = abs.(o.confpeak * sqrt((o.dist[1] * o.multiplier) / o.dof) / Phi)  # abs.(o.confpeak) * quantile(Normal{T}(zero(T),one(T)), p/2) / Phi
   else
 		halfwidth = T(-1.5) * Phi .* sqrtNaN.(diag(getV(o)))
 		any(isnan.(halfwidth)) && return

src/precompiler.jl

@@ -19,7 +19,6 @@ using SnoopPrecompile, StableRNGs
       wildboottest(T, T[0 0 0 0 1], T[.04]; resp, predexog, predendog, inst, rng, clustid=idcoarse)
       wildboottest(T, T[0 0 0 0 1], T[.04]; resp, predexog, predendog, inst, rng, clustid=idcoarse, arubin=true)
       wildboottest(T, T[0 0 0 1 0], T[.04]; resp, predexog, predendog, inst, rng, clustid=idgranular, liml=true)
-      nothing
     end
   end
 end

src/utilities.jl

@@ -18,6 +18,13 @@ invsym(X::Symmetric) =
 		#=Symmetric=#(fill(eltype(X)(NaN), size(X)))
 	end
 
+eigvalsNaN(X) =
+	try
+		eigvalsNaN(X)
+	catch _
+		fill(eltype(X)(NaN), size(X))
+	end
+
 function invsymsingcheck(X)  # inverse of symmetric matrix, checking for singularity
 	iszero(nrows(X)) && (return (false, #=Symmetric=#(X)))
 	X, ipiv, info = LinearAlgebra.LAPACK.sytrf!('U', Matrix(X))
@@ -26,6 +33,7 @@ function invsymsingcheck(X)  # inverse of symmetric matrix, checking for singula
 	singular, #=Symmetric=#(X)
 end
 
+
 @inline colsum(X::AbstractArray) = iszero(length(X)) ? similar(X, 1, size(X)[2:end]...) : sum(X, dims=1)
 @inline colsum(X::AbstractArray{Bool}) = iszero(length(X)) ? Array{Int}(undef, 1, size(X)[2:end]...) : sum(X, dims=1)  # type-stable
 @inline rowsum(X::AbstractArray) = vec(sum(X, dims=2))

test/unittests.log

@@ -61,7 +61,7 @@ p       = 0.7332
 
 boottest (post_self=.05) (post=-.02) (selfemployed=-.15), reps(9999) weight(webb)
 F(3, 7) = 0.8534
-p       = 0.7442
+p       = 0.7401
 
 
 regress hasinsurance selfemployed post post_self
@@ -236,7 +236,7 @@ CI = [-0.894 -0.3861; 0.3094 1.253]
 
 boottest tenure, noci bootcluster(individual) weight(webb)
 z  = 4.0267
-p  = 0.0000
+p  = 0.0080
 
 boottest tenure, nograph bootcluster(collgrad) cluster(collgrad industry) weight(webb) reps(9999)
 z  = 3.9420