diff --git a/docs/make.jl b/docs/make.jl index 2de0bcd..25c7976 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -15,6 +15,7 @@ makedocs(; "Home" => "index.md" "Executable" => "Executable.md" "Package" => "Package.md" + "KSS pruning" => "Pruning.md" ], ) diff --git a/docs/src/Package.md b/docs/src/Package.md index 1dacf1c..6c4d333 100644 --- a/docs/src/Package.md +++ b/docs/src/Package.md @@ -23,6 +23,10 @@ Depth = 3 ## Main Function +See [KSS leave-one-worker-out pruning](Pruning.md) for the sample-construction +algorithm and its distinction from the subsequent match- or observation-level +bias correction. + ```@docs leave_out_KSS(y,first_id,second_id;controls, do_lincom , Z_lincom , lincom_labels , settings) @@ -90,4 +94,3 @@ region[findall(region.==-1)].=0 ``` - diff --git a/docs/src/Pruning.md b/docs/src/Pruning.md new file mode 100644 index 0000000..1c674df --- /dev/null +++ b/docs/src/Pruning.md @@ -0,0 +1,87 @@ +# KSS leave-one-worker-out pruning + +## What the pruning step guarantees + +Kline, Saggio, and Sølvsten (KSS) estimate variance components on a sample whose +firm mobility network remains connected after removing any one worker's complete +observation history. This is a property of the selected raw-observation sample. +It is not a choice between leaving out a worker, a match, or an observation. + +The input must be sorted by worker and observation order, as required by +`leave_out_KSS`. Pruning preserves original observation indices and returns +dense worker and firm identifiers. + +The implementation follows the released MATLAB routines at commit +[`8b957ff`](https://github.com/rsaggio87/LeaveOutTwoWay/tree/8b957ffeb10b8465a3584fceb0265cccc48379e1): + +- [`connected_set.m`](https://github.com/rsaggio87/LeaveOutTwoWay/blob/8b957ffeb10b8465a3584fceb0265cccc48379e1/codes/connected_set.m) + constructs the initial firm-transition graph. +- [`pruning_unbal_v3.m`](https://github.com/rsaggio87/LeaveOutTwoWay/blob/8b957ffeb10b8465a3584fceb0265cccc48379e1/codes/pruning_unbal_v3.m) + removes articulation workers to a fixed point. +- [`build_adj.m`](https://github.com/rsaggio87/LeaveOutTwoWay/blob/8b957ffeb10b8465a3584fceb0265cccc48379e1/codes/build_adj.m) + rebuilds the mover-induced firm-transition graph. +- [`leave_out_KSS.m`](https://github.com/rsaggio87/LeaveOutTwoWay/blob/8b957ffeb10b8465a3584fceb0265cccc48379e1/codes/leave_out_KSS.m) + applies the final raw-observation singleton filter. + +The Julia implementation is an independent translation of that behavior. + +## Algorithm + +Given sorted raw observations `(worker, firm)`: + +1. Construct the firm-transition graph from consecutive observations within + workers and retain its unique largest component, ranked by number of firms. +2. Identify movers and deduplicate their worker-firm pairs. +3. Construct the mover-only bipartite graph and find articulation workers. +4. Remove every articulation worker's complete raw-observation history. +5. Rebuild the mover-induced firm-transition graph and retain its unique largest + component, again ranked by number of firms. +6. Repeat steps 2--5 until no articulation worker remains. +7. Drop workers represented by only one raw observation and relabel identifiers. + +The routine raises an `ArgumentError` if largest firm components are tied or if +the mover graph disappears. Choosing a tied component by storage order would +make the selected sample arbitrary. + +## Difference from the previous selector + +The previous Julia selector built the full worker-firm bipartite graph and ranked +components by their total number of worker and firm vertices. The released +MATLAB code ranks components of the firm-transition graph by number of firms. +Those objectives can select different surviving samples. + +The graph used to detect articulation workers is not the source of this +difference. Stayers are degree-one leaves. A leaf is never a cut vertex, and +attaching leaves to firm vertices cannot change whether a worker vertex is a cut +vertex. Consequently, the full and mover-only bipartite graphs identify the same +articulation workers on a common sample. They can diverge when the components +remaining after deletion are ranked. + +For example, consider a two-worker, three-firm `K(2, 3)` mover core joined by an +articulation worker to one firm with six two-observation stayers. After removing +the bridge worker, firm ranking retains the three-firm mover core (six +observations), while total-vertex ranking retains the one-firm stayer component +(twelve observations). + +The previous routine also retained a stale worker count across iterations. This +added isolated phantom vertices and wasted memory, but it does not appear to +alter articulation workers or component selection. Graph dimensions are now +recomputed from the current sample in every round. + +Both selectors happen to agree on the public `test.csv` fixture: + +- 71,614 input observations; +- 63,462 observations after the initial firm component; +- 56,044 final observations; +- 28,022 workers, 1,684 firms, 34,436 matches, and 6,414 movers. + +Agreement on this fixture does not establish that the algorithms are equivalent; +the synthetic component-ranking example demonstrates the difference. + +## Pruning versus `leave_out_level` + +`get_leave_one_out_set` always constructs the leave-one-worker-out-connected raw +observation sample described above. The `VCHDFESettings.leave_out_level` value +does not change pruning. Its `"match"` and `"obs"` settings apply only +later in `leave_out_estimation`, when the package computes the KSS bias +correction on the already selected sample. diff --git a/docs/src/index.md b/docs/src/index.md index ac9bf53..ef9151c 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -10,6 +10,6 @@ Please report any bug or problem [here](https://github.com/HighDimensionalEconLa # Contents ```@contents -Pages = ["Executable.md", "Package.md"] +Pages = ["Executable.md", "Package.md", "Pruning.md"] Depth = 3 -``` \ No newline at end of file +``` diff --git a/src/VarianceComponentsHDFE.jl b/src/VarianceComponentsHDFE.jl index c7a7ea0..fe3e18c 100644 --- a/src/VarianceComponentsHDFE.jl +++ b/src/VarianceComponentsHDFE.jl @@ -11,6 +11,7 @@ using ArgParse using FixedEffects #include("init.jl") +include("kss_pruning.jl") include("leave_out_correction.jl") include("parameters_settings.jl") include("laplacians/Laplacians.jl") diff --git a/src/kss_pruning.jl b/src/kss_pruning.jl new file mode 100644 index 0000000..bf9d083 --- /dev/null +++ b/src/kss_pruning.jl @@ -0,0 +1,191 @@ +function _check_pruning_inputs(y, first_id, second_id, obs_id = nothing) + n_observations = length(y) + length(first_id) == n_observations || + throw(ArgumentError("y and first_id must have the same length")) + length(second_id) == n_observations || + throw(ArgumentError("y and second_id must have the same length")) + obs_id === nothing || length(obs_id) == n_observations || + throw(ArgumentError("y and obs_id must have the same length")) + n_observations > 0 || + throw(ArgumentError("KSS pruning requires at least one observation")) + return nothing +end + +function _dense_ids(ids) + levels = unique(sort(ids)) + return Int.(indexin(ids, levels)), levels +end + +function _select_and_relabel(y, first_id, second_id, obs_id, selection) + selected_y = y[selection] + selected_first_id = first_id[selection] + selected_second_id = second_id[selection] + selected_obs_id = obs_id[selection] + + dense_first_id, _ = _dense_ids(selected_first_id) + dense_second_id, _ = _dense_ids(selected_second_id) + return ( + obs_id = selected_obs_id, + y = selected_y, + first_id = dense_first_id, + second_id = dense_second_id, + ) +end + +function _firm_transition_graph(first_id, second_id, firms = unique(sort(second_id))) + firm_id = Int.(indexin(second_id, firms)) + graph = Graph(length(firms)) + for observation = 2:length(first_id) + same_worker = first_id[observation] == first_id[observation - 1] + changed_firm = firm_id[observation] != firm_id[observation - 1] + same_worker && changed_firm && + add_edge!(graph, firm_id[observation - 1], firm_id[observation]) + end + return graph, firms +end + +function _unique_largest_firm_component(graph, stage) + components = connected_components(graph) + isempty(components) && + throw(ArgumentError("KSS pruning found no firm component $stage")) + + component_sizes = length.(components) + largest_size = maximum(component_sizes) + largest_components = findall(==(largest_size), component_sizes) + if length(largest_components) != 1 + throw( + ArgumentError( + "KSS pruning found $(length(largest_components)) tied largest firm " * + "components $stage ($largest_size firms each)", + ), + ) + end + return components[only(largest_components)] +end + +function _mover_workers(first_id, second_id) + first_firm = zeros(Int, maximum(first_id)) + movers = falses(length(first_firm)) + for observation in eachindex(first_id) + worker = first_id[observation] + firm = second_id[observation] + if first_firm[worker] == 0 + first_firm[worker] = firm + elseif first_firm[worker] != firm + movers[worker] = true + end + end + return findall(movers) +end + +function _articulation_workers(first_id, second_id) + mover_workers = _mover_workers(first_id, second_id) + isempty(mover_workers) && throw( + ArgumentError( + "KSS pruning cannot continue because the mover graph has disappeared", + ), + ) + + mover_index = zeros(Int, maximum(first_id)) + for (index, worker) in enumerate(mover_workers) + mover_index[worker] = index + end + mover_rows = findall(worker -> mover_index[worker] > 0, first_id) + pairs = Tuple{Int,Int}[] + for row in mover_rows + push!(pairs, (mover_index[first_id[row]], second_id[row])) + end + unique!(pairs) + + n_movers = length(mover_workers) + graph = Graph(n_movers + maximum(second_id)) + for (worker, firm) in pairs + add_edge!(graph, worker, n_movers + firm) + end + + articulation_workers = filter(vertex -> vertex <= n_movers, articulation(graph)) + return mover_workers[articulation_workers] +end + +function _largest_mover_firm_sample(y, first_id, second_id, obs_id) + mover_workers = _mover_workers(first_id, second_id) + isempty(mover_workers) && throw( + ArgumentError( + "KSS pruning cannot continue because the mover graph has disappeared", + ), + ) + + is_mover = falses(maximum(first_id)) + is_mover[mover_workers] .= true + mover_rows = findall(worker -> is_mover[worker], first_id) + firms = collect(1:maximum(second_id)) + graph, firms = _firm_transition_graph( + first_id[mover_rows], + second_id[mover_rows], + firms, + ) + component = _unique_largest_firm_component(graph, "after articulation deletion") + connected_firms = Set(firms[component]) + selection = findall(firm -> firm in connected_firms, second_id) + return _select_and_relabel(y, first_id, second_id, obs_id, selection) +end + +function _official_find_connected_set(y, first_id, second_id, settings) + _check_pruning_inputs(y, first_id, second_id) + obs_id = collect(1:length(y)) + graph, firms = _firm_transition_graph(first_id, second_id) + component = _unique_largest_firm_component(graph, "in the initial sample") + connected_firms = Set(firms[component]) + selection = findall(firm -> firm in connected_firms, second_id) + return _select_and_relabel(y, first_id, second_id, obs_id, selection) +end + +function _official_prunning_connected_set( + y, + first_idvar, + second_idvar, + obs_id, + settings, +) + _check_pruning_inputs(y, first_idvar, second_idvar, obs_id) + first_id, _ = _dense_ids(first_idvar) + second_id, _ = _dense_ids(second_idvar) + + while true + bad_workers = _articulation_workers(first_id, second_id) + n_bad_workers = length(bad_workers) + settings.print_level > 1 && println( + "Number of $(settings.first_id_display_small) that are articulation points: ", + n_bad_workers, + ) + n_bad_workers == 0 && break + + is_bad_worker = falses(maximum(first_id)) + is_bad_worker[bad_workers] .= true + selection = findall(worker -> !is_bad_worker[worker], first_id) + isempty(selection) && throw( + ArgumentError( + "KSS pruning removed all observations before reaching a fixed point", + ), + ) + + sample = _select_and_relabel(y, first_id, second_id, obs_id, selection) + sample = _largest_mover_firm_sample( + sample.y, + sample.first_id, + sample.second_id, + sample.obs_id, + ) + obs_id = sample.obs_id + y = sample.y + first_id = sample.first_id + second_id = sample.second_id + end + + return ( + obs_id = obs_id, + y = y, + first_id = first_id, + second_id = second_id, + ) +end diff --git a/src/leave_out_correction.jl b/src/leave_out_correction.jl index ef929ff..fc8979b 100644 --- a/src/leave_out_correction.jl +++ b/src/leave_out_correction.jl @@ -96,7 +96,6 @@ end ## Methods and Types Definitions - #1) Finds AKM largest connected set """ $(SIGNATURES) @@ -110,54 +109,7 @@ Returns a tuple of observation belonging to the largest connected set with the c * `settings`: settings based on data type `VCHDFESettings`. Please see the reference provided below. """ function find_connected_set(y, first_idvar, second_idvar, settings) - - #Observation identifier to join later the FE - obs_id = collect(1:size(y,1)) - - seconds = unique(sort(second_idvar)) - second_id = indexin(second_idvar, seconds) - - firsts = unique(sort(first_idvar)) - first_id = indexin(first_idvar, firsts) - - #Save ids just in case - second_id_old=second_idvar; - first_id_old=first_idvar; - - n_seconds=length(seconds) - n_firsts=length(firsts) - second_id=second_id.+n_firsts - - graph_size=length(seconds)+length(firsts) - G=Graph(graph_size) - for i in 1:size(y,1) - add_edge!(G, second_id[i], first_id[i]) - end - - cs=connected_components(G) - #(settings.print_level > 1) && println("Largest connected set has been found") - pos = indexin( [ maximum(size.(cs,1))] , size.(cs,1) )[1] - lcs=cs[pos] - - connected_seconds=lcs[lcs.>n_firsts] - - sel=findall(in(connected_seconds),second_id) - - obs_id = [obs_id[x] for x in sel ] - yvec = [y[x] for x in sel ] - second_id = [second_idvar[x] for x in sel ] - first_id = [first_idvar[x] for x in sel ] - - #Relabel seconds - seconds = unique(sort(second_id)) - second_id = indexin(second_id, seconds) - - #Relabel firsts - first_ids = unique(sort(first_id)) - first_id = indexin(first_id, first_ids) - - return (obs_id = obs_id , y = yvec , first_id = first_id, second_id = second_id ) - + return _official_find_connected_set(y, first_idvar, second_idvar, settings) end #2) Pruning and finding Leave-Out Largest connected set @@ -174,88 +126,13 @@ This function prunes the dataset from articulation points. If the first identifi * `settings`: settings based on data type `VCHDFESettings`. Please see the reference provided below. """ function prunning_connected_set(yvec, first_idvar, second_idvar, obs_id, settings) - - seconds = unique(sort(second_idvar)) - second_id = indexin(second_idvar, seconds) - firsts = unique(sort(first_idvar)) - first_id = indexin(first_idvar, firsts) - - nbadfirsts=1 - while nbadfirsts>0 - - n_seconds=length(seconds) - n_firsts=length(firsts) - second_id=second_id.+n_firsts - - graph_size=n_firsts + n_seconds - lcs_graph=Graph(graph_size) - for i in 1:size(yvec,1) - add_edge!(lcs_graph, second_id[i], first_id[i]) - end - - #get articulation vertex - artic_vertex = articulation(lcs_graph) - - sel=findall(x->x==nothing, indexin(first_id,artic_vertex)) - - bad_firsts = artic_vertex[artic_vertex.<=n_firsts] - nbadfirsts = size(bad_firsts,1) - - settings.print_level > 1 && println("Number of $(settings.first_id_display_small) that are articulation points: ", nbadfirsts) - - #Restrict the sample - yvec = [yvec[x] for x in sel ] - second_id = [second_id[x] for x in sel ] - first_id = [first_id[x] for x in sel ] - obs_id = [obs_id[x] for x in sel ] - - #Relabel seconds - seconds = unique(sort(second_id)) - second_id = indexin(second_id, seconds) - - #Relabel firsts - first_ids = unique(sort(first_id)) - first_id = indexin(first_id, first_ids) - - n_firsts=maximum(first_id) - n_seconds=maximum(second_id); - - second_id = second_id .+ n_firsts - - #Constructing new Graph - G=Graph(n_firsts+n_seconds) - Nprimeprime = size(sel,1) - for i in 1:Nprimeprime - add_edge!(G, second_id[i], first_id[i]) - end - - #Find Largest Connected Set (LCS) - cs=connected_components(G) - - pos = indexin( [ maximum(size.(cs,1))] , size.(cs,1) )[1] - lcs=cs[pos] - - connected_seconds=lcs[lcs.>n_firsts] - - - sel=findall(in(connected_seconds),second_id) - - obs_id = [obs_id[x] for x in sel ] - yvec = [yvec[x] for x in sel ] - second_id = [second_id[x] for x in sel ] - first_id = [first_id[x] for x in sel ] - - #Relabel seconds - seconds = unique(sort(second_id)) - second_id = indexin(second_id, seconds) - - #Relabel firsts - first_ids = unique(sort(first_id)) - first_id = indexin(first_id, first_ids) - - end - - return (obs_id = obs_id , y = yvec , first_id = first_id, second_id = second_id ) + return _official_prunning_connected_set( + yvec, + first_idvar, + second_idvar, + obs_id, + settings, + ) end #3) Drops Single Observations diff --git a/test/runtests.jl b/test/runtests.jl index 106424e..b36bb2c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,7 @@ using VarianceComponentsHDFE using Test +include("test_pruning.jl") #run_full_test = get(ENV, "VCHDFE_FULL_TEST", "1") == "1" ? true : false diff --git a/test/test_pruning.jl b/test/test_pruning.jl new file mode 100644 index 0000000..6579dde --- /dev/null +++ b/test/test_pruning.jl @@ -0,0 +1,196 @@ +using CSV +using DataFrames +using LightGraphs + +const PRUNING_SETTINGS = VCHDFESettings(print_level = 0) + +function count_movers(first_id, second_id) + firms_by_worker = [Set{Int}() for _ = 1:maximum(first_id)] + for (worker, firm) in zip(first_id, second_id) + push!(firms_by_worker[worker], firm) + end + return count(firms -> length(firms) > 1, firms_by_worker) +end + +@testset "Official KSS pruning" begin + @testset "firm-ranked component after articulation deletion" begin + first_id = [ + 1, 1, 1, + 2, 2, 2, + 3, 3, + 4, 4, + 5, 5, + 6, 6, + 7, 7, + 8, 8, + 9, 9, + ] + second_id = [ + 1, 2, 3, + 1, 2, 3, + 3, 4, + 4, 4, + 4, 4, + 4, 4, + 4, 4, + 4, 4, + 4, 4, + ] + y = collect(101.0:120.0) + + result = get_leave_one_out_set( + y, + first_id, + second_id, + PRUNING_SETTINGS, + nothing, + ) + + @test result.obs == collect(1:6) + @test result.y == y[result.obs] + @test result.first_id == [1, 1, 1, 2, 2, 2] + @test result.second_id == [1, 2, 3, 1, 2, 3] + end + + @testset "stayer leaves do not change articulation workers" begin + mover_graph = Graph(7) + for worker = 1:2, firm = 1:3 + add_edge!(mover_graph, worker, 3 + firm) + end + add_edge!(mover_graph, 3, 3 + 3) + add_edge!(mover_graph, 3, 3 + 4) + + full_graph = Graph(13) + for worker = 1:2, firm = 1:3 + add_edge!(full_graph, worker, 9 + firm) + end + add_edge!(full_graph, 3, 9 + 3) + add_edge!(full_graph, 3, 9 + 4) + for worker = 4:9 + add_edge!(full_graph, worker, 9 + 4) + end + + mover_articulation_workers = filter(vertex -> vertex <= 3, articulation(mover_graph)) + full_articulation_workers = filter(vertex -> vertex <= 9, articulation(full_graph)) + @test mover_articulation_workers == full_articulation_workers == [3] + end + + @testset "articulation cascade reaches a later round" begin + first_id = [1, 1, 1, 2, 2] + second_id = [1, 2, 3, 1, 2] + y = zeros(length(first_id)) + error, output = mktemp() do _, stream + error = redirect_stdout(stream) do + try + get_leave_one_out_set( + y, + first_id, + second_id, + VCHDFESettings(print_level = 2), + nothing, + ) + nothing + catch exception + exception + end + end + flush(stream) + seekstart(stream) + return error, read(stream, String) + end + + @test error isa ArgumentError + messages = collect(eachmatch(r"articulation points: 1", output)) + @test length(messages) == 2 + end + + @testset "singletons are counted from raw rows" begin + first_id = [1, 1, 2, 2, 3] + second_id = [1, 2, 1, 2, 1] + y = collect(1.0:5.0) + + result = get_leave_one_out_set( + y, + first_id, + second_id, + PRUNING_SETTINGS, + nothing, + ) + + @test result.obs == collect(1:4) + @test result.first_id == [1, 1, 2, 2] + @test result.second_id == [1, 2, 1, 2] + end + + @testset "ambiguous and unidentified samples fail clearly" begin + tied_first_id = [1, 1, 2, 2, 3, 3, 4, 4] + tied_second_id = [1, 2, 1, 2, 3, 4, 3, 4] + tied_y = zeros(length(tied_first_id)) + @test_throws ArgumentError find_connected_set( + tied_y, + tied_first_id, + tied_second_id, + PRUNING_SETTINGS, + ) + + stayer_first_id = [1, 1, 2, 2] + stayer_second_id = [1, 1, 1, 1] + stayer_y = zeros(length(stayer_first_id)) + error = try + get_leave_one_out_set( + stayer_y, + stayer_first_id, + stayer_second_id, + PRUNING_SETTINGS, + nothing, + ) + nothing + catch exception + exception + end + @test error isa ArgumentError + @test occursin("mover", lowercase(sprint(showerror, error))) + end + + @testset "public KSS fixture" begin + fixture_path = joinpath(pkgdir(VarianceComponentsHDFE), "test.csv") + fixture = DataFrame(CSV.File(fixture_path; header = false)) + first_id = fixture[:, 1] + second_id = fixture[:, 2] + y = fixture[:, 4] + + @test length(y) == 71_614 + + connected = find_connected_set(y, first_id, second_id, PRUNING_SETTINGS) + @test length(connected.obs_id) == 63_462 + + pruned = prunning_connected_set( + connected.y, + connected.first_id, + connected.second_id, + connected.obs_id, + PRUNING_SETTINGS, + ) + final = drop_single_obs( + pruned.y, + pruned.first_id, + pruned.second_id, + pruned.obs_id, + ) + + @test length(final.obs_id) == 56_044 + @test maximum(final.first_id) == 28_022 + @test maximum(final.second_id) == 1_684 + @test length(Set(zip(final.first_id, final.second_id))) == 34_436 + @test count_movers(final.first_id, final.second_id) == 6_414 + + combined = get_leave_one_out_set( + y, + first_id, + second_id, + PRUNING_SETTINGS, + nothing, + ) + @test combined.obs == final.obs_id + end +end