From f8751d15229d72252365c974dd1084e6fd2931f0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Bayo=C3=A1n=20Santiago=20Calder=C3=B3n=2C=20Ph?= =?UTF-8?q?D?= Date: Sun, 4 Dec 2022 21:34:27 -0500 Subject: [PATCH] Fixes bug in full ranking (#84) * Fixed bug in full ranking --- .github/workflows/ci.yml | 4 ++-- Project.toml | 13 +++++++------ data/_72.csv | 41 ++++++++++++++++++++++++++++++++++++++++ src/solvers.jl | 16 ++++++++-------- test/linear_models.jl | 25 ++++++++++++++++++++++++ test/runtests.jl | 1 + 6 files changed, 84 insertions(+), 16 deletions(-) create mode 100644 data/_72.csv diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a31d16d..a2735a8 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -35,9 +35,9 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - julia-version: ['1.7'] + julia-version: ['1.8'] julia-arch: [x64] - os: [ubuntu-20.04] + os: [ubuntu-22.04] steps: - uses: actions/checkout@v2 - uses: julia-actions/setup-julia@v1 diff --git a/Project.toml b/Project.toml index a68c5e3..5165754 100644 --- a/Project.toml +++ b/Project.toml @@ -3,7 +3,7 @@ uuid = "4d6a76a9-bfbc-5492-8924-cf6ed7875f06" license = "ISC" description = "Econometrics in Julia" maintainer = "José Bayoán Santiago Calderón" -version = "0.2.9" +version = "0.2.10" [deps] CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597" @@ -24,9 +24,9 @@ TableOperations = "ab02a1b2-a7df-11e8-156e-fb1833f50b87" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" [compat] -CategoricalArrays = "0.8, 0.9, 0.10" -Distributions = "0.23, 0.24, 0.25" -FillArrays = "0.9, 0.10, 0.11" +CategoricalArrays = "0.10" +Distributions = "0.25" +FillArrays = "0.13" ForwardDiff = "0.10" JLLWrappers = "1" MLJModelInterface = "1" @@ -34,7 +34,7 @@ Optim = "1" Parameters = "0.12" StatsAPI = "1" StatsBase = "0.33" -StatsFuns = "0.9, 1" +StatsFuns = "1" StatsModels = "0.6" TableOperations = "1" Tables = "1" @@ -49,6 +49,7 @@ RDatasets = "ce6b1742-4840-55fa-b093-852dadbb1d8b" StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Weave = "44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9" +WooldridgeDatasets = "b46f53ff-21ea-4df8-83e7-e19d2972755a" [targets] -test = ["CSV", "Documenter", "MLJBase", "MLJModelInterface", "RDatasets", "StatsAPI", "Test", "Weave"] +test = ["CSV", "Documenter", "MLJBase", "MLJModelInterface", "RDatasets", "StatsAPI", "Test", "Weave", "WooldridgeDatasets"] diff --git a/data/_72.csv b/data/_72.csv new file mode 100644 index 0000000..a4982c8 --- /dev/null +++ b/data/_72.csv @@ -0,0 +1,41 @@ +V,d,p,c1,c2,r2,p̃,j +1.2027419218197932,27.703656833029616,0.0,false,false,false,0.0,1 +1.075147895949145,26.275771622919024,0.0,false,false,false,0.0,1 +0.9313659184439649,30.378364038034068,0.0,false,false,true,0.0,1 +1.0761276167040825,29.045430133131287,0.0,false,false,true,0.0,1 +1.157924142455271,26.6897777510207,0.0,false,false,true,0.0,1 +1.2602159294494486,23.564996779476353,0.0,false,false,true,0.0,1 +1.0691480410425707,24.302754839676016,0.0,false,true,false,0.0,1 +0.9747923308125769,27.93962459498397,0.0,false,false,false,0.0,1 +1.0864642002730331,29.43234468759983,0.0,false,false,true,0.0,1 +0.8280195291069027,36.454816370348006,0.0,true,false,false,0.0,1 +0.9938150611635577,25.882930650932604,0.0,false,false,false,0.0,1 +1.1902971518664769,28.003360206792166,0.0,false,false,false,0.0,1 +0.8440652135748151,27.186667855146023,0.0,false,false,true,0.0,1 +1.1942135229240978,24.793572848996735,0.0,false,false,false,0.0,1 +0.9749415395991344,32.231070524734776,0.0,false,false,false,0.0,1 +1.0630157511408558,29.954964998040836,0.0,false,false,true,0.0,1 +1.2088620878197598,28.88111761520257,0.0,false,false,false,0.0,1 +1.1554818612875875,26.25362099456051,0.0,false,false,true,0.0,1 +1.0683169195107736,26.66594108388161,0.0,false,false,false,0.0,1 +1.1825776372250463,28.501424506184375,0.0,false,false,false,0.0,1 +0.9216968403274366,21.43410899036141,5.000056366090773,false,false,false,5.125,2 +0.9038509857862627,21.49017975132197,5.25003596784075,false,false,false,5.37498733358276,2 +0.3543068600742174,25.487400244731294,6.25,false,false,true,6.189546924592175,2 +0.7443615843260774,22.498270803257157,5.5,false,false,true,5.43255861288492,2 +0.9204453025709346,20.470734056384202,4.75,false,false,true,4.735834317058153,2 +0.8185203035625981,21.275447440483994,4.75,false,false,true,4.64249459776033,2 +0.663460159111698,25.179971661184624,6.5,false,true,false,6.56955705314014,2 +0.33270708470409316,26.91964183437916,6.75,false,false,false,6.77710246638078,2 +0.723848557740205,25.53805792453679,7.0,false,false,true,6.8758529838688025,2 +0.6844470209911019,21.98950570308865,5.5,true,false,false,5.50813576231531,2 +0.6778781468474899,24.248051662239188,6.25,false,false,false,6.291796915708559,2 +0.7395466413433043,24.37288297318075,6.000006846676847,false,false,false,6.1249742106215415,2 +0.6944880579705579,23.729820384574147,6.4999977785530465,false,false,true,6.375047150072365,2 +0.9026005539496791,20.93161210697483,4.75,false,false,false,4.838826807533771,2 +0.6972541480874805,25.660159077369066,7.0,false,false,false,7.0883822269378065,2 +0.4954034749087639,24.646645084956198,6.0,false,false,true,5.880590893610366,2 +0.9095561199444187,22.33350651705032,5.5,false,false,false,5.527812278473579,2 +0.8492111796989497,23.63584813451155,6.249953503069051,false,false,true,6.125016675359465,2 +0.9170291009885081,22.16336314179665,5.75,false,false,false,5.763628709164413,2 +0.6726866287469144,25.692005034400466,6.500134661603563,false,false,false,6.625,2 diff --git a/src/solvers.jl b/src/solvers.jl index 0fc75f0..cdecc40 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -19,7 +19,7 @@ F = bunchkaufman(Hermitian(Z̃' * Diagonal(wts) * Z̃), true, check = false) bkr = count(x -> abs(x) ≥ √eps(), F.D) if bkr < size(F, 2) - lin_ind = sort!(invperm(F.p)[1:bkr]) + lin_ind = sort!(F.p[sortperm(F.p)[1:bkr]]) count(x -> lin_ind > size(X, 2) + 1) ≥ size(z, 2) || throw(ArgumentError("Insufficient number of instruments.")) Z̃ = Z̃[:,lin_ind] @@ -35,7 +35,7 @@ F = bunchkaufman(Hermitian(X̃' * Diagonal(wts) * X̃), true, check = false) bkr = count(x -> abs(x) ≥ √eps(), F.D) if bkr < size(F, 2) - lin_ind = sort!(invperm(F.p)[1:bkr]) + lin_ind = sort!(F.p[sortperm(F.p)[1:bkr]]) X̃ = convert(Matrix{Float64}, X̃[:,lin_ind]) F = bunchkaufman(Hermitian(X̃' * Diagonal(wts) * X̃), true) else @@ -63,7 +63,7 @@ end F = bunchkaufman(Hermitian(Z̃' * Diagonal(wts) * Z̃), true, check = false) bkr = count(x -> abs(x) ≥ √eps(), F.D) if bkr < size(F, 2) - lin_ind = sort!(invperm(F.p)[1:bkr]) + lin_ind = sort!(F.p[sortperm(F.p)[1:bkr]]) count(x -> lin_ind > size(X, 2) + 1) ≥ size(z, 2) || throw(ArgumentError("Insufficient number of instruments.")) Z̃ = Z̃[:,lin_ind] @@ -79,7 +79,7 @@ end F = bunchkaufman(Hermitian(X̃' * Diagonal(wts) * X̃), true, check = false) bkr = count(x -> abs(x) ≥ √eps(), F.D) if bkr < size(F, 2) - lin_ind = sort!(invperm(F.p)[1:bkr]) + lin_ind = sort!(F.p[sortperm(F.p)[1:bkr]]) X̃ = convert(Matrix{Float64}, X̃[:,lin_ind]) F = bunchkaufman(Hermitian(X̃' * Diagonal(wts) * X̃), true) else @@ -112,7 +112,7 @@ end F = bunchkaufman(Hermitian(Z̃' * Diagonal(w) * Z̃), true, check = false) bkr = count(x -> abs(x) ≥ √eps(), F.D) if bkr < size(F, 2) - lin_ind = sort!(invperm(F.p)[1:bkr]) + lin_ind = sort!(F.p[sortperm(F.p)[1:bkr]]) F₀ = bunchkaufman(Hermitian(X' * Diagonal(w) * X), true, check = false) num_of_linear_independent_exogenous_vars = count(x -> abs(x) ≥ √eps(), F₀.D) (length(lin_ind) - num_of_linear_independent_exogenous_vars) ≥ size(z, 2) || @@ -130,7 +130,7 @@ end F = bunchkaufman(Hermitian(X̃' * Diagonal(w) * X̃), true, check = false) bkr = count(x -> abs(x) ≥ √eps(), F.D) if bkr < size(F, 2) - lin_ind = sort!(invperm(F.p)[1:bkr]) + lin_ind = sort!(F.p[sortperm(F.p)[1:bkr]]) X̃ = convert(Matrix{Float64}, X̃[:,lin_ind]) F = bunchkaufman(Hermitian(X̃' * Diagonal(w) * X̃), true) else @@ -179,7 +179,7 @@ end F = qr(X, ColumnNorm()) qrr = count(x -> abs(x) ≥ √eps(), diag(F.R)) if qrr < size(F, 2) - lin_ind = sort!(invperm(F.p)[1:qrr]) + lin_ind = sort!(F.p[sortperm(F.p)[1:bkr]]) X = convert(Matrix{Float64}, X[:,lin_ind]) F = qr(X, ColumnNorm()) else @@ -234,7 +234,7 @@ end F = qr(X, ColumnNorm()) qrr = count(x -> abs(x) ≥ √eps(), diag(F.R)) if qrr < size(F, 2) - lin_ind = sort!(invperm(F.p)[1:qrr]) + lin_ind = sort!(F.p[sortperm(F.p)[1:bkr]]) X = convert(Matrix{Float64}, X[:,lin_ind]) else lin_ind = collect(1:size(F, 2)) diff --git a/test/linear_models.jl b/test/linear_models.jl index 89955cf..7fee7ce 100644 --- a/test/linear_models.jl +++ b/test/linear_models.jl @@ -1,6 +1,16 @@ # Linear Models @testset "Linear Models" begin @testset "Balanced Panel Data" begin + # Full ranking instruments #72 + data = CSV.read( + joinpath(pkgdir(Econometrics), "data", "_72.csv"), + DataFrame + ) + model = fit( + EconometricModel, + @formula(V ~ 1 + (d + p ~ (c1 + c2 + r2 + p̃) * j)), + data) + @test isa(model, EconometricModel) # Balanced Panel Data data = dataset("Ecdat", "Crime") # Pooling @@ -76,6 +86,21 @@ ) @test sprint(show, model) == "One-way Random Effect Model\nLongitudinal dataset: County, Year\nBalanced dataset with 90 panels of length 7\nindividual error component: 0.0162\nidiosyncratic error component: 0.0071\nρ: 0.8399\nNumber of observations: 630\nNull Loglikelihood: 2225.79\nLoglikelihood: 2226.01\nR-squared: 0.0007\nWald: 0.15 ∼ F(3, 626) ⟹ Pr > F = 0.9291\nFormula: CRMRTE ~ PrbConv + AvgSen + PrbPris\nVariance Covariance Estimator: OIM\n──────────────────────────────────────────────────────────────────────────────────────\n PE SE t-value Pr > |t| 2.50% 97.50%\n──────────────────────────────────────────────────────────────────────────────────────\n(Intercept) 0.0309293 0.00266706 11.5968 <1e-27 0.0256918 0.0361668\nPrbConv -3.96634e-5 0.000198471 -0.199845 0.8417 -0.000429413 0.000350086\nAvgSen 8.2428e-5 0.000127818 0.644887 0.5192 -0.000168575 0.000333431\nPrbPris -0.000123327 0.00404272 -0.030506 0.9757 -0.00806226 0.00781561\n──────────────────────────────────────────────────────────────────────────────────────" + wagepan = DataFrame(wooldridge("wagepan")) + # From #83 + model = fit( + RandomEffectsEstimator, + @formula(lwage ~ educ + black + hisp + exper + exper^2 + married + union + year), + wagepan, + panel = :nr, + time = :year, + contrasts = Dict(:year => DummyCoding()) + ) + @test coef(model) ≈ [ + 0.023586377, 0.091876276, -0.139376726, 0.021731732, 0.105754520, + -0.004723943, 0.063986022, 0.106134429, 0.040462003, 0.030921157, + 0.020280640, 0.043118708, 0.057815458, 0.091947584, 0.134928917 + ] rtol = 1e-3 # IV Pooling model = fit(EconometricModel, @formula(CRMRTE ~ PrbConv + (AvgSen ~ PrbPris)), data) @test sprint(show, model) == diff --git a/test/runtests.jl b/test/runtests.jl index fffe640..59304db 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using Econometrics: ConvergenceException using CSV, RDatasets, CategoricalArrays using MLJBase: MLJBase using MLJModelInterface: MLJModelInterface +using WooldridgeDatasets: wooldridge for file in ["exceptions", "linear_models", "mlogit", "ologit", "formula_display", "mlj", "docs"] include("$file.jl")