Fixes bug in full ranking (#84)

* Fixed bug in full ranking
Nosferican · Dec 5, 2022 · f8751d1 · f8751d1 · Nosferican · Dec 5, 2022
1 parent 1e0fb4d
commit f8751d1
Show file tree

Hide file tree

Showing 6 changed files with 84 additions and 16 deletions.
diff --git 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
@@ -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,17 +24,17 @@ 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"
 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
@@ -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
@@ -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
@@ -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
@@ -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")