From 8971a91feeabae0cafdd39199a85a89fd61c536e Mon Sep 17 00:00:00 2001 From: droodman Date: Thu, 19 Jan 2023 19:29:56 -0500 Subject: [PATCH] compute halfwidth for ARubin in way less likely to produce NaN --- Manifest.toml | 71 +++++++++++++++++++++++++++++++++++++++++++++- Project.toml | 2 +- src/WRE.jl | 2 +- src/estimators.jl | 4 +-- src/plot-CI.jl | 8 +++--- src/precompiler.jl | 1 - src/utilities.jl | 8 ++++++ test/unittests.log | 4 +-- 8 files changed, 88 insertions(+), 12 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 11c0e32..1ff5436 100644 --- a/Manifest.toml +++ b/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" diff --git a/Project.toml b/Project.toml index 581b542..a5fe63c 100644 --- a/Project.toml +++ b/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" diff --git a/src/WRE.jl b/src/WRE.jl index 473807d..f6c422a 100644 --- a/src/WRE.jl +++ b/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 diff --git a/src/estimators.jl b/src/estimators.jl index 9a27cd6..2fff4a0 100644 --- a/src/estimators.jl +++ b/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 diff --git a/src/plot-CI.jl b/src/plot-CI.jl index 4ffbcb3..6f8da2b 100644 --- a/src/plot-CI.jl +++ b/src/plot-CI.jl @@ -41,7 +41,7 @@ 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 @@ -49,9 +49,9 @@ function plot!(o::StrBootTest{T}) where T 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 diff --git a/src/precompiler.jl b/src/precompiler.jl index 08b8f99..6ee37ad 100644 --- a/src/precompiler.jl +++ b/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 diff --git a/src/utilities.jl b/src/utilities.jl index 3473da1..84d1e03 100644 --- a/src/utilities.jl +++ b/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)) diff --git a/test/unittests.log b/test/unittests.log index 20be9da..7ac76eb 100644 --- a/test/unittests.log +++ b/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