---
title: "The gender wage gap in Krabbistan, 1995–2015"
subtitle: "A synthetic three-census study inspired by Blau and Kahn (2017)"
author: "Krabbs"
date: 2026-08-11
format:
  html:
    embed-resources: true
    toc: true
    code-fold: true
    code-summary: "Show R code"
    theme: cosmo
execute:
  echo: true
  warning: false
  message: false
  freeze: false
---

## Overview

Krabbistan is a fictional small nation whose labor market is calibrated to be broadly reminiscent of the United States. This report creates census microdata for 1995, 2005, and 2015 and studies the gender wage gap using annual own earnings divided by annual hours.

The analysis follows the organizing logic of Blau and Kahn (2017):

1. raw mean log-wage gaps overall and by parenthood × marriage;
2. male-reference Oaxaca–Blinder decompositions of the mean gap, paralleling their Table 4;
3. decompositions of changes from 1995 into changing characteristics, changing male prices, and changing unexplained gaps, paralleling their Table 5; and
4. RIF-regression decompositions at the 10th, 50th, and 90th unconditional percentiles, as a rough analogue to their Table 6.

The data are entirely synthetic and the results are pedagogical, not estimates for any real country. The executable analysis is in R because this rendering environment does not contain a licensed Stata runtime. A complete Stata translation is supplied alongside the report.

```{r}
#| label: setup
set.seed(42042017)

library(dplyr)
library(tidyr)
library(ggplot2)
library(knitr)
library(scales)

dir.create("data", showWarnings = FALSE)
dir.create("results", showWarnings = FALSE)
dir.create("figures", showWarnings = FALSE)

years <- c(1995L, 2005L, 2015L)
n_per_wave <- 20000L
beta_groups <- list(
  Education = paste0("educ_", 2:5),
  Experience = c("experience", "experience_sq"),
  Industry = paste0("industry_", 2:8),
  Occupation = paste0("occupation_", 2:8)
)

theme_set(theme_minimal(base_size = 12.5))
gap_colors <- c(
  Education = "#56B4E9", Experience = "#E69F00",
  Industry = "#009E73", Occupation = "#CC79A7",
  Unexplained = "#777777"
)
```

## Simulating the census

Each wave contains 20,000 residents ages 18–64. Education rises over time, women overtake men educationally, occupational and industrial segregation decline, women’s labor-force participation rises, and the direct wage penalty attached to being female shrinks. Parenthood and marriage retain additional female wage penalties, although both also diminish.

```{r}
#| label: simulation-functions
softmax_draw <- function(scores) {
  scores <- scores - apply(scores, 1, max)
  probabilities <- exp(scores)
  probabilities <- probabilities / rowSums(probabilities)
  cumulative <- t(apply(probabilities, 1, cumsum))
  u <- runif(nrow(scores))
  1L + rowSums(u > cumulative)
}

simulate_wave <- function(year, n = n_per_wave) {
  t <- (year - 1995) / 10
  sex <- rbinom(n, 1, 0.505) # 0 male, 1 female

  # Triangular working-age distribution with an early male age advantage.
  age <- round(runif(n, 18, 65) + (1 - sex) * (3.5 - 1.2 * t))
  age <- pmin(64L, pmax(18L, age))

  educ_latent <- rnorm(n) + 0.16 * t + sex * (-0.06 + 0.03 * t) +
    0.010 * (age - 40)
  education <- as.integer(cut(
    educ_latent,
    breaks = c(-Inf, -0.90, -0.20, 0.50, 1.18, Inf),
    labels = FALSE
  ))
  education_years <- c(10, 12, 14, 16, 18)[education]
  college <- education >= 4

  p_married <- plogis(
    -0.15 + 0.105 * (age - 25) - 0.0031 * (age - 41)^2 +
      0.17 * college - 0.10 * sex - 0.10 * t
  )
  married <- rbinom(n, 1, p_married)
  p_former <- plogis(-2.2 + 0.075 * (age - 30) + 0.12 * sex)
  formerly <- married == 0 & runif(n) < p_former
  marital_status <- ifelse(married == 1, 2L, ifelse(formerly, 3L, 1L))

  p_parent <- plogis(
    -1.15 + 0.115 * (age - 24) - 0.0030 * (age - 39)^2 +
      0.72 * married - 0.13 * (education - 3) - 0.05 * t
  )
  parent <- rbinom(n, 1, p_parent)

  edc <- education - 3
  occ_scores <- cbind(
    -0.55 + 0.48 * edc - sex * (0.72 - 0.20 * t),
    -0.75 + 0.70 * edc + sex * (0.12 + 0.03 * t),
    0.05 + 0.18 * edc + sex * (0.45 - 0.09 * t),
    -0.15 + 0.05 * edc + sex * 0.10,
    -0.05 - 0.18 * edc + sex * (0.53 - 0.07 * t),
    -0.35 - 0.08 * edc - sex * (1.15 - 0.16 * t),
    -0.25 - 0.16 * edc - sex * (0.76 - 0.11 * t),
    -0.45 - 0.24 * edc - sex * (0.28 - 0.04 * t)
  )
  occupation <- softmax_draw(occ_scores)

  professional <- occupation %in% 1:3
  manual <- occupation %in% 6:8
  ind_scores <- cbind(
    -0.15 + 0.35 * edc + 0.25 * professional - sex * (0.25 - 0.06 * t),
    -0.10 + 0.22 * edc + sex * (0.65 - 0.08 * t),
    0.05 - 0.12 * edc + 0.38 * manual - sex * (0.50 - 0.08 * t),
    -0.55 - 0.08 * edc + 0.70 * (occupation == 6) - sex * (1.20 - 0.16 * t),
    0.05 - 0.05 * edc + 0.18 * (occupation %in% 3:5),
    -0.35 + 0.42 * edc + 0.20 * professional - sex * (0.18 - 0.04 * t),
    -0.50 + 0.20 * edc + sex * 0.02,
    -0.20 - 0.10 * edc + 0.20 * (occupation == 5) + sex * 0.18
  )
  industry <- softmax_draw(ind_scores)

  potential_experience <- pmax(age - education_years - 6, 0)
  career_break <- sex * parent * (4.0 - 0.75 * t) +
    sex * married * (0.8 - 0.15 * t)
  effective_experience <- pmax(potential_experience - career_break, 0)

  p_lfp <- plogis(
    1.50 - 0.72 * sex + 0.22 * sex * t - 0.48 * sex * parent -
      0.14 * sex * married + 0.16 * (education - 3) -
      0.0035 * (age - 42)^2
  )
  lfp <- rbinom(n, 1, p_lfp)

  educ_premium <- c(0.00, 0.10, 0.22, 0.48, 0.72)[education]
  occ_premium <- c(0.45, 0.38, 0.15, 0.05, -0.15, 0.08, -0.05, -0.15)[occupation]
  ind_premium <- c(0.25, 0.05, 0.12, 0.15, -0.05, 0.30, 0.05, -0.10)[industry]
  base_penalty <- c(`1995` = 0.100, `2005` = 0.090, `2015` = 0.080)[as.character(year)]
  parent_penalty <- c(`1995` = 0.040, `2005` = 0.040, `2015` = 0.040)[as.character(year)]
  married_penalty <- c(`1995` = 0.020, `2005` = 0.020, `2015` = 0.015)[as.character(year)]

  log_hourly_wage <- 2.42 + educ_premium +
    0.045 * effective_experience - 0.00055 * effective_experience^2 +
    occ_premium + ind_premium -
    sex * (base_penalty + parent_penalty * parent + married_penalty * married) +
    rnorm(n, 0, 0.32)
  hourly_wage <- pmin(250, pmax(4, exp(log_hourly_wage)))

  log_hours <- log(1930) - 0.10 * sex - 0.13 * sex * parent -
    0.04 * sex * married + 0.035 * t * sex + rnorm(n, 0, 0.20)
  annual_hours <- ifelse(lfp == 1, round(pmin(3200, pmax(250, exp(log_hours)))), 0)
  own_earnings <- ifelse(lfp == 1, round(hourly_wage * annual_hours), 0)

  data.frame(
    year = as.integer(year), sex = as.integer(sex),
    own_earnings, age = as.integer(age),
    marital_status = as.integer(marital_status), married = as.integer(married),
    parent = as.integer(parent), education = as.integer(education),
    education_years = as.integer(education_years),
    annual_hours = as.integer(annual_hours), lfp = as.integer(lfp),
    industry = as.integer(industry), occupation = as.integer(occupation)
  )
}
```

```{r}
#| label: simulate-and-save
census <- bind_rows(lapply(years, simulate_wave))

# R-native and portable CSV files are generated directly. Stata-format files
# are written with the recommended foreign package included with R.
saveRDS(census, "data/krabbistan_census_1995_2015.rds")
write.csv(census, "data/krabbistan_census_1995_2015.csv", row.names = FALSE)
foreign::write.dta(census, "data/krabbistan_census_1995_2015.dta", version = 12)
for (yr in years) {
  wave <- census[census$year == yr, ]
  saveRDS(wave, sprintf("data/krabbistan_census_%s.rds", yr))
  foreign::write.dta(wave, sprintf("data/krabbistan_census_%s.dta", yr), version = 12)
}

coverage <- census |>
  group_by(year) |>
  summarise(
    observations = n(), female_share = mean(sex),
    labor_force_participation = mean(lfp),
    mean_age = mean(age), .groups = "drop"
  )
kable(coverage, digits = 3, caption = "Synthetic census coverage")
```

The public-use files retain all residents, including zero earnings and zero hours for nonparticipants. The wage analysis restricts attention to ages 25–64, positive earnings, labor-force participation, and at least 500 annual hours. Hourly wages are calculated—not simulated as a released variable—as own earnings divided by annual hours.

```{r}
#| label: analysis-sample
work <- census |>
  filter(age >= 25, age <= 64, lfp == 1, own_earnings > 0, annual_hours >= 500) |>
  mutate(
    hourly_wage = own_earnings / annual_hours,
    log_wage = log(hourly_wage),
    experience = pmax(age - education_years - 6, 0),
    experience_sq = experience^2
  )

for (j in 2:5) work[[paste0("educ_", j)]] <- as.numeric(work$education == j)
for (j in 2:8) {
  work[[paste0("industry_", j)]] <- as.numeric(work$industry == j)
  work[[paste0("occupation_", j)]] <- as.numeric(work$occupation == j)
}
saveRDS(work, "data/krabbistan_wage_analysis.rds")
foreign::write.dta(work, "data/krabbistan_wage_analysis.dta", version = 12)
```

## Raw wage gaps

```{r}
#| label: raw-gap-estimates
raw_sex <- work |>
  group_by(year, sex) |>
  summarise(mean_log_wage = mean(log_wage), N = n(), .groups = "drop") |>
  pivot_wider(names_from = sex, values_from = c(mean_log_wage, N), names_prefix = "sex_")

raw_overall <- raw_sex |>
  transmute(
    year,
    male_mean_log_wage = mean_log_wage_sex_0,
    female_mean_log_wage = mean_log_wage_sex_1,
    male_minus_female_log_gap = mean_log_wage_sex_0 - mean_log_wage_sex_1,
    female_male_wage_ratio = exp(mean_log_wage_sex_1 - mean_log_wage_sex_0),
    N = N_sex_0 + N_sex_1
  )

family_sex <- work |>
  group_by(year, parent, married, sex) |>
  summarise(mean_log_wage = mean(log_wage), N = n(), .groups = "drop") |>
  pivot_wider(names_from = sex, values_from = c(mean_log_wage, N), names_prefix = "sex_")

raw_family <- family_sex |>
  mutate(
    family_group = paste0(
      ifelse(parent == 1, "Parent", "Not parent"), " × ",
      ifelse(married == 1, "married", "not married")
    ),
    male_minus_female_log_gap = mean_log_wage_sex_0 - mean_log_wage_sex_1,
    female_male_wage_ratio = exp(mean_log_wage_sex_1 - mean_log_wage_sex_0),
    N = N_sex_0 + N_sex_1
  ) |>
  select(year, parent, married, family_group, male_minus_female_log_gap,
         female_male_wage_ratio, N)

write.csv(raw_overall, "results/raw_gap_overall.csv", row.names = FALSE)
write.csv(raw_family, "results/raw_gap_family.csv", row.names = FALSE)

kable(
  raw_overall |>
    mutate(
      male_minus_female_log_gap = round(male_minus_female_log_gap, 3),
      female_male_wage_ratio = percent(female_male_wage_ratio, accuracy = 0.1)
    ),
  caption = "Raw mean wage gaps"
)
```

```{r}
#| label: raw-gap-plot
#| fig-width: 7.5
#| fig-height: 4.5
p_raw <- ggplot(raw_overall, aes(year, male_minus_female_log_gap)) +
  geom_line(linewidth = 1.1, color = "#0072B2") +
  geom_point(size = 3, color = "#0072B2") +
  scale_x_continuous(breaks = years) +
  labs(
    title = "Krabbistan's raw gender wage gap narrows",
    x = "Census year", y = "Male − female mean log wage gap"
  )
p_raw
ggsave("figures/raw_gap_over_time.png", p_raw, width = 7.5, height = 4.5, dpi = 180)
```

```{r}
#| label: family-table
family_display <- raw_family |>
  transmute(
    year, `Family category` = family_group,
    `Log gap` = round(male_minus_female_log_gap, 3),
    `Female/male wage ratio` = percent(female_male_wage_ratio, accuracy = 0.1),
    N
  )
kable(family_display, caption = "Raw wage gaps by parenthood × marriage")
```

```{r}
#| label: family-plot
#| fig-width: 8.5
#| fig-height: 5
p_family <- ggplot(
  raw_family,
  aes(year, male_minus_female_log_gap, color = family_group)
) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.6) +
  scale_x_continuous(breaks = years) +
  scale_color_brewer(palette = "Dark2") +
  labs(
    title = "Raw gaps by parenthood and marriage",
    x = "Census year", y = "Male − female mean log wage gap", color = NULL
  ) +
  theme(legend.position = "bottom")
p_family
ggsave("figures/raw_gap_family_groups.png", p_family, width = 8.5, height = 5, dpi = 180)
```

Parents generally have larger gaps because family status affects women’s wage offers, labor-force attachment, and hours more strongly. Married parents have the largest gap in the first two waves; by 2015, compositional selection among working married mothers changes the ordering. The persistent parenthood gradient is more robust than the marriage gradient.

## Mean Oaxaca–Blinder decomposition

For year $t$, let $\Delta\bar X_t=\bar X_{mt}-\bar X_{ft}$ and let $\widehat\beta_{mt}$ be the male wage coefficients. The male-reference decomposition is

$$
G_t = \Delta\bar X_t'\widehat\beta_{mt} + U_t,
$$

where $G_t$ is the male–female mean log-wage gap. Group contributions sum the relevant dummy or polynomial terms. The experience proxy is potential experience, $\max(\text{age}-\text{schooling}-6,0)$, and its square.

```{r}
#| label: ob-functions
design_matrix <- function(frame, specification) {
  groups <- c("Education", "Experience")
  if (specification == "Full") groups <- c(groups, "Industry", "Occupation")
  variables <- unlist(beta_groups[groups], use.names = FALSE)
  X <- cbind(`_cons` = 1, as.matrix(frame[, variables, drop = FALSE]))
  list(X = X, variables = c("_cons", variables), groups = groups)
}

fit_ols <- function(y, X) as.numeric(qr.solve(X, y))

table4_rows <- list()
ob_store <- list()
row_id <- 1L

for (yr in years) {
  wave <- work[work$year == yr, ]
  men <- wave[wave$sex == 0, ]
  women <- wave[wave$sex == 1, ]
  gap <- mean(men$log_wage) - mean(women$log_wage)

  for (specification in c("Human capital", "Full")) {
    spec_key <- ifelse(specification == "Full", "Full", "HC")
    dm <- design_matrix(men, spec_key)
    df <- design_matrix(women, spec_key)
    beta <- fit_ols(men$log_wage, dm$X)
    names(beta) <- dm$variables
    dx <- colMeans(dm$X) - colMeans(df$X)
    contributions <- sapply(dm$groups, function(group) {
      variables <- beta_groups[[group]]
      sum(dx[variables] * beta[variables])
    })
    explained <- sum(contributions)
    unexplained <- gap - explained
    ob_store[[paste(yr, specification, sep = "::")]] <- list(
      beta = beta, dx = dx, gap = gap,
      explained = explained, unexplained = unexplained
    )
    values <- c(
      contributions,
      `Total explained` = explained,
      Unexplained = unexplained,
      `Total gap` = gap
    )
    table4_rows[[row_id]] <- data.frame(
      year = yr, specification, component = names(values),
      log_points = as.numeric(values),
      percent_of_gap = 100 * as.numeric(values) / gap
    )
    row_id <- row_id + 1L
  }
}

table4 <- bind_rows(table4_rows)
write.csv(table4, "results/table4_ob_decomposition.csv", row.names = FALSE)
```

### Table 4 analogue

```{r}
#| label: table4
table4_display <- table4 |>
  mutate(entry = sprintf("%.3f (%.1f%%)", log_points, percent_of_gap)) |>
  select(specification, component, year, entry) |>
  pivot_wider(names_from = year, values_from = entry)
kable(
  table4_display,
  caption = "Male-reference decomposition: log points (percent of total gap)"
)
```

```{r}
#| label: table4-plot
#| fig-width: 8.5
#| fig-height: 5
table4_plot_data <- table4 |>
  filter(
    specification == "Full",
    component %in% c("Education", "Experience", "Industry", "Occupation", "Unexplained")
  ) |>
  mutate(component = factor(
    component,
    levels = c("Education", "Experience", "Industry", "Occupation", "Unexplained")
  ))

p_table4 <- ggplot(table4_plot_data, aes(factor(year), log_points, fill = component)) +
  geom_col(width = 0.68) +
  scale_fill_manual(values = gap_colors) +
  labs(
    title = "Full Oaxaca–Blinder decomposition of the mean gap",
    x = "Census year", y = "Log points", fill = NULL
  ) +
  theme(legend.position = "bottom")
p_table4
ggsave("figures/table4_stacked_components.png", p_table4, width = 8.5, height = 5, dpi = 180)
```

Education contributes little because women’s schooling catches up and then surpasses men’s. Industry and especially occupation remain important sorting margins. The unexplained component includes the simulated direct gender penalty, unmeasured career interruptions, family penalties, and sampling noise; it must not be interpreted mechanically as discrimination alone.

## Change decomposition relative to 1995

Blau and Kahn’s Table 5 emphasizes that the change in the gap depends on both changing gender differences in characteristics and changing prices. Two exact paths are reported. With 1995 coefficients as the base,

$$
G_t-G_{95} = (\Delta\bar X_t-\Delta\bar X_{95})'\widehat\beta_{m,95}
+\Delta\bar X_t'(\widehat\beta_{mt}-\widehat\beta_{m,95})
+(U_t-U_{95}).
$$

The alternative evaluates mean changes at current coefficients and coefficient changes at the 1995 characteristic gap.

```{r}
#| label: table5-calculate
table5_rows <- list()
row_id <- 1L

for (specification in c("Human capital", "Full")) {
  base <- ob_store[[paste(1995, specification, sep = "::")]]
  active_groups <- if (specification == "Full") names(beta_groups) else c("Education", "Experience")

  for (yr in c(2005L, 2015L)) {
    current <- ob_store[[paste(yr, specification, sep = "::")]]
    for (path in c("1995 coefficients / current gaps", "Current coefficients / 1995 gaps")) {
      group_rows <- lapply(active_groups, function(group) {
        variables <- beta_groups[[group]]
        dx0 <- base$dx[variables]
        dxt <- current$dx[variables]
        b0 <- base$beta[variables]
        bt <- current$beta[variables]
        if (startsWith(path, "1995")) {
          means <- sum((dxt - dx0) * b0)
          coefficients <- sum(dxt * (bt - b0))
        } else {
          means <- sum((dxt - dx0) * bt)
          coefficients <- sum(dx0 * (bt - b0))
        }
        data.frame(
          year = yr, specification, path, component = group,
          changing_means = means, changing_coefficients = coefficients
        )
      })
      group_rows <- bind_rows(group_rows)
      extras <- data.frame(
        year = yr, specification, path,
        component = c("All covariates", "Unexplained gap", "Total gap change"),
        changing_means = c(
          sum(group_rows$changing_means),
          current$unexplained - base$unexplained,
          current$gap - base$gap
        ),
        changing_coefficients = c(sum(group_rows$changing_coefficients), NA, NA)
      )
      table5_rows[[row_id]] <- bind_rows(group_rows, extras)
      row_id <- row_id + 1L
    }
  }
}

table5 <- bind_rows(table5_rows)
write.csv(table5, "results/table5_changes_from_1995.csv", row.names = FALSE)
```

### Table 5 analogue

```{r}
#| label: table5
table5_display <- table5 |>
  filter(specification == "Full") |>
  mutate(
    `Changing means` = ifelse(is.na(changing_means), "", sprintf("%.3f", changing_means)),
    `Changing coefficients` = ifelse(
      is.na(changing_coefficients), "", sprintf("%.3f", changing_coefficients)
    )
  ) |>
  select(year, path, component, `Changing means`, `Changing coefficients`)
kable(table5_display, caption = "Contributions to the change in the gap relative to 1995: full specification")
```

```{r}
#| label: table5-plot
#| fig-width: 8.8
#| fig-height: 5.2
means_plot <- table5 |>
  filter(
    specification == "Full",
    path == "1995 coefficients / current gaps",
    component %in% names(beta_groups)
  ) |>
  transmute(year, component = paste("Means:", component), contribution = changing_means)

other_plot <- table5 |>
  filter(
    specification == "Full",
    path == "1995 coefficients / current gaps",
    component %in% c("All covariates", "Unexplained gap")
  ) |>
  transmute(
    year,
    component = ifelse(component == "All covariates", "All coefficient changes", "Unexplained-gap change"),
    contribution = ifelse(component == "All coefficient changes", changing_coefficients, changing_means)
  )

table5_plot <- bind_rows(means_plot, other_plot)
p_table5 <- ggplot(table5_plot, aes(factor(year), contribution, fill = component)) +
  geom_col(width = 0.68) +
  geom_hline(yintercept = 0, linewidth = 0.4) +
  labs(
    title = "Why the gap changed relative to 1995",
    subtitle = "1995-coefficient/current-characteristic-gap path",
    x = "Ending census year", y = "Contribution to change in log gap", fill = NULL
  ) +
  theme(legend.position = "bottom")
p_table5
ggsave("figures/table5_stacked_changes.png", p_table5, width = 8.8, height = 5.2, dpi = 180)
```

Negative bars narrow the gap. The two paths allocate interactions between changing characteristics and changing coefficients differently, but both sum to the same observed gap change once the unexplained component is included.

## Distributional decomposition: Table 6 analogue

The original paper uses the Chernozhukov–Fernández-Val–Melly distribution-regression framework. Here, a deliberately rougher recentered-influence-function (RIF) approximation is used. For percentile $p$, each sex-specific wage is transformed to

$$
\operatorname{RIF}(Y;q_p)=q_p+\frac{p-\mathbf{1}\{Y\le q_p\}}{f_Y(q_p)},
$$

where $f_Y(q_p)$ is estimated with a Gaussian kernel. Separate male and female RIF regressions then yield an Oaxaca-style covariate effect and wage-coefficient effect at each unconditional percentile.

```{r}
#| label: table6-calculate
rif_values <- function(y, p) {
  q <- as.numeric(quantile(y, p, names = FALSE, type = 7))
  h <- max(1.06 * sd(y) * length(y)^(-1 / 5), 1e-4)
  density <- mean(dnorm((y - q) / h) / h)
  list(rif = q + (p - as.numeric(y <= q)) / density, quantile = q)
}

table6_rows <- list()
row_id <- 1L
for (yr in years) {
  wave <- work[work$year == yr, ]
  men <- wave[wave$sex == 0, ]
  women <- wave[wave$sex == 1, ]
  for (p in c(0.10, 0.50, 0.90)) {
    rm <- rif_values(men$log_wage, p)
    rf <- rif_values(women$log_wage, p)
    for (specification in c("Human capital", "Full")) {
      key <- ifelse(specification == "Full", "Full", "HC")
      dm <- design_matrix(men, key)
      df <- design_matrix(women, key)
      bm <- fit_ols(rm$rif, dm$X)
      bf <- fit_ols(rf$rif, df$X)
      covariate <- sum((colMeans(dm$X)[-1] - colMeans(df$X)[-1]) * bm[-1])
      coefficient <- sum(colMeans(df$X) * (bm - bf))
      table6_rows[[row_id]] <- data.frame(
        year = yr, percentile = as.integer(100 * p), specification,
        covariate_effect = covariate,
        coefficient_effect = coefficient,
        sum_effects = covariate + coefficient,
        raw_quantile_gap = rm$quantile - rf$quantile
      )
      row_id <- row_id + 1L
    }
  }
}

table6 <- bind_rows(table6_rows)
write.csv(table6, "results/table6_rif_quantile_decomposition.csv", row.names = FALSE)
```

```{r}
#| label: table6
table6_display <- table6 |>
  mutate(across(
    c(covariate_effect, coefficient_effect, sum_effects, raw_quantile_gap),
    ~ round(.x, 3)
  ))
kable(table6_display, caption = "RIF-regression decomposition by unconditional wage percentile")
```

```{r}
#| label: table6-plot
#| fig-width: 11
#| fig-height: 5.5
table6_plot <- table6 |>
  select(year, percentile, specification, covariate_effect, coefficient_effect) |>
  pivot_longer(
    c(covariate_effect, coefficient_effect),
    names_to = "effect", values_to = "log_points"
  ) |>
  mutate(
    effect = recode(
      effect,
      covariate_effect = "Covariates",
      coefficient_effect = "Wage coefficients"
    ),
    cell = paste0(year, "\np", percentile)
  )

p_table6 <- ggplot(table6_plot, aes(cell, log_points, fill = effect)) +
  geom_col(width = 0.72) +
  facet_wrap(~ specification, nrow = 1) +
  scale_fill_manual(values = c(Covariates = "#0072B2", `Wage coefficients` = "#D55E00")) +
  labs(
    title = "RIF approximation to Blau–Kahn Table 6",
    x = "Year and unconditional percentile", y = "Log-point contribution", fill = NULL
  ) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "bottom")
p_table6
ggsave("figures/table6_rif_stacked.png", p_table6, width = 11, height = 5.5, dpi = 180)
```

The percentile results separate compositional sorting from differences in the wage structure. Because RIF regression is a local linear approximation, `sum_effects` need not equal the raw sample quantile gap exactly. Standard errors are omitted; a substantive application should bootstrap the entire quantile, density, and regression procedure.

## Where the design is feasible in IPUMS International

The [IPUMS International sample catalogue](https://international.ipums.org/international-action/samples) and harmonized-variable documentation were audited on 11 August 2026 to identify real-data settings for the same exercise. A *direct match* must have (i) at least two census microdata samples whose endpoints are at least ten years apart, (ii) the demographic, earnings, education, employment, industry, and occupation variables below in each qualifying wave, and (iii) at least 100,000 person records in each wave. The record threshold is an analytical rule used here, not an IPUMS standard; it leaves ample scope for employed-sex and family-status cells after sample restrictions.

The hours requirement is deliberately flexible. A series passes if it supports at least one of three wage-sample constructions:

1. **Computed hourly wage:** earnings divided by exact usual or actual hours (`HRSWORK1`, `HRSUSUAL1`, `HRSMAIN`, or `HRSACTUAL1`), after aligning reference periods;
2. **Reported hourly wage:** a wage rate reported directly in hourly units, without dividing earnings by hours; or
3. **Full-time earnings sample:** wage earnings analyzed only among workers identifiable as full-time, using `HRSFULL` or a source variable with a documented full-time threshold.

Route 3 estimates an earnings gap among full-time workers, not a literal hourly-wage gap. It is therefore a useful fallback but not estimand-equivalent to routes 1 or 2.

```{r}
#| label: ipums-variable-crosswalk
ipums_crosswalk <- tibble::tribble(
  ~Krabbistan_measure, ~IPUMS_variable, ~Interpretation,
  "Sex", "SEX", "Reported sex",
  "Own earnings", "INCWAGE or INCEARN", "Wage-and-salary income, or broader earned income when the former is unavailable",
  "Age", "AGE", "Age in years",
  "Marital status", "MARST", "Harmonized marital status",
  "Parenthood", "NCHILD/MOMLOC/POPLOC; CHBORN/CHSURV/LASTBYR", "Co-resident parent for both sexes; fertility-history proxy for women",
  "Education", "EDATTAIN", "Internationally harmonized attainment",
  "Wage denominator/sample", "HRSWORK1, HRSUSUAL1, HRSMAIN, HRSACTUAL1, or HRSFULL", "Exact hours, a direct hourly rate, or a documented full-time indicator",
  "Labor-force status", "EMPSTAT", "Employment status",
  "Industry", "INDGEN", "Coarse general industry",
  "Occupation", "OCCISCO", "ISCO-harmonized occupation"
)

kable(
  ipums_crosswalk,
  col.names = c("Krabbistan measure", "IPUMS field", "Interpretation"),
  caption = "Crosswalk used in the feasibility screen"
)
```

### Reconstructing parenthood

`NCHILD` is not essential. The preferred replacement is a **co-resident-parent indicator** constructed by reversing the IPUMS parent pointers: a person is a parent if another household member's `MOMLOC` or `POPLOC` equals that person's `PERNUM`. This has the same conceptual limitation as `NCHILD`---children outside the household are missed---but it works symmetrically for women and men. Where pointers are unavailable, household relationship and spouse links can often identify children of the head or spouse, though they miss more complex families.

Fertility histories provide a second route. `CHBORN > 0`, `CHSURV > 0`, or a valid `LASTBYR` identifies women who have ever had a birth. These variables capture nonresident children but usually apply only to women, sometimes only within a restricted age range. They therefore cannot by themselves define the same parent category for men and women. The recommended design uses co-resident parenthood as the common main measure and women's fertility history as a sensitivity check.

Removing `NCHILD` as a hard screen and admitting `INCEARN` (total own labor income) materially expands the usable set. **Brazil, Jamaica, Mexico, the United States, and Venezuela** support repeated-census analyses with a common parent construct; **Canada** is conditionally usable with a source-level relationship reconstruction or an asymmetric fertility-history sensitivity analysis. **Puerto Rico** is a technical match but a U.S. territory. Brazil, Mexico, and Venezuela use broader labor earnings rather than wage-and-salary income alone, so their estimand includes self-employment income.

```{r}
#| label: ipums-feasibility-table
ipums_feasibility <- tibble::tribble(
  ~jurisdiction, ~assessment, ~samples, ~minimum_records, ~parent_route, ~wage_status,
  "Brazil", "Feasible; broader earnings", "1991, 2000, 2010 censuses", "17,045,712", "Co-resident links; CHBORN sensitivity", "INCEARN / exact hours; includes business and farm income",
  "Canada", "Conditional", "1981, 1991, 2001 censuses", "486,875", "CHBORN for women; source relationship reconstruction needed for a sex-symmetric measure", "INCWAGE / exact hours or full-time restriction",
  "Jamaica", "Feasible", "1982, 1991, 2001 censuses", "205,179", "Co-resident links; CHBORN sensitivity", "INCWAGE / exact hours",
  "Mexico", "Feasible; broader earnings", "1990, 2000, 2010 censuses", "8,118,242", "Co-resident links; CHBORN sensitivity", "INCEARN / exact hours; includes self-employment income",
  "United States", "Feasible", "1980, 1990, 2000 censuses", "11,343,120", "Co-resident links; CHBORN available through 1990", "INCWAGE / exact hours",
  "Venezuela", "Feasible; broader earnings", "1981 and 2001 censuses", "1,441,266", "Co-resident links; CHBORN available in 2001", "INCEARN / exact hours; 20-year endpoint comparison",
  "Puerto Rico", "Technical match; territory", "1990 and 2000 censuses", "177,655", "Co-resident links; CHBORN in 1990", "INCWAGE / exact hours"
)

write.csv(
  ipums_feasibility,
  "results/ipums_international_feasibility.csv",
  row.names = FALSE
)

kable(
  ipums_feasibility,
  col.names = c(
    "Country or territory", "Assessment", "Potential waves",
    "Minimum IPUMS records", "Parent construction", "Wage construction/deviation"
  ),
  caption = "Repeated-census settings that can support the analysis after flexible parent construction"
)
```

### Twenty-country longlist and the binding constraint

The broader audit produces a longlist of **21 sovereign countries**, plus Puerto Rico. All have large IPUMS samples and some substantial portion of the desired demographic/labor design. It is not honest, however, to label all twenty-one feasible wage-gap settings: once parenthood is reconstructed, the binding limitation is **repeated individual earnings**, not children. The table makes those deviations explicit instead of silently changing the estimand to household income or occupation-imputed wages.

```{r}
#| label: ipums-longlist-table
ipums_longlist <- tibble::tribble(
  ~country, ~candidate_censuses, ~minimum_records, ~parent_option, ~status_or_blocker,
  "Brazil", "1991/2000/2010", "17,045,712", "Resident child links; fertility history", "Feasible with broader earned income",
  "Canada", "1981/1991/2001", "486,875", "Fertility history; source relationship reconstruction", "Conditional: common male/female parent proxy needs source work",
  "Jamaica", "1982/1991/2001", "205,179", "Resident child links; fertility history", "Feasible",
  "Mexico", "1990/2000/2010", "8,118,242", "Resident child links; fertility history", "Feasible with broader earned income",
  "United States", "1980/1990/2000", "11,343,120", "Resident child links; fertility history", "Feasible",
  "Venezuela", "1981/2001", "1,441,266", "Resident child links; fertility history", "Feasible with broader earned income",
  "Indonesia", "1980/1990 censuses", "912,544", "Resident child links; fertility history", "Censuses lack earnings; 1976/1995 surveys are usable alternatives",
  "Panama", "1990/2000/2010", "232,737", "Resident child links; fertility history", "Earnings present; qualifying waves lack hours/full-time route",
  "Israel", "1983/1995/2008", "403,474", "Resident links through 1995; fertility history", "Earnings present; no repeated hours/full-time route",
  "Trinidad and Tobago", "1980/1990/2000/2011", "105,464", "Resident child links; fertility history", "Individual earnings only in 2000",
  "Uruguay", "1975/1985/1996/2011", "279,994", "Resident child links; fertility history", "Individual wage income only in the 2006 household survey",
  "Switzerland", "1970/1980/1990/2000", "312,538", "Resident child links; fertility history in 2000", "No individual earnings measure",
  "Ecuador", "1990/2001/2010", "966,234", "Resident child links; fertility history", "No individual earnings measure",
  "Spain", "2001/2011", "2,039,274", "Resident child links; earlier fertility history", "No individual earnings measure",
  "France", "1990/1999/2006/2011", "2,360,854", "Resident child links", "No individual earnings measure",
  "Greece", "1991/2001/2011", "969,407", "Resident child links; fertility history", "No individual earnings measure",
  "Italy", "2001/2011", "2,968,065", "Resident child links", "No census individual earnings measure",
  "Mauritius", "1990/2000/2011", "106,710", "Resident child links; fertility history", "No individual earnings measure",
  "Nicaragua", "1995/2005", "435,728", "Resident child links; fertility history", "No individual earnings measure",
  "El Salvador", "1992/2007", "510,760", "Resident child links; fertility history", "No individual earnings measure",
  "South Africa", "1996/2001/2011", "3,621,164", "Resident child links; fertility history", "Excellent labor covariates, but IPUMS provides household—not own—income in these waves",
  "Puerto Rico", "1990/2000", "177,655", "Resident child links; fertility history", "Feasible, but a U.S. territory"
)

write.csv(ipums_longlist, "results/ipums_twenty_country_longlist.csv", row.names = FALSE)

kable(
  ipums_longlist,
  col.names = c(
    "Country or territory", "Candidate census years", "Minimum records",
    "Parent option", "Assessment or binding blocker"
  ),
  caption = "Expanded IPUMS International longlist: 21 sovereign countries plus Puerto Rico"
)
```

This expanded screen also clarifies why India is not in the census longlist: its income-bearing IPUMS files are employment surveys, not census microdata. Mexico has both survey and census material, but the census sequence itself is usable through `INCEARN`.

South Africa is indeed a strong candidate for the **employment** and child-penalty parts of the project. IPUMS supplies the 1996, 2001, and 2011 censuses, millions of observations, parent links and fertility histories, schooling, labor-force status, and industry/occupation (through 2007 in the harmonized series). The available income fields for these census samples are household income or household-head income, not the worker's own earnings. Household income is mechanically contaminated by partners' earnings and cannot serve as the dependent variable in an individual Oaxaca–Blinder gender-pay decomposition. South Africa should therefore remain on the longlist, with another microdata source required for the wage portion.

### Child Penalty Atlas source audit

Kleven, Landais, and Leite-Mariante's *Child Penalty Atlas* covers **134 countries**. Its core outcome is employment, so inclusion in the Atlas does not imply that its cited file contains own earnings, hours, education, occupation, and industry. The following audit starts from Appendix Table A.1 and asks whether the cited source could support the richer gender-pay analysis here. “Accessible” includes free registration or a research-use application; it does not mean unrestricted anonymous download.

```{r}
#| label: atlas-country-audit
atlas_source_text <- "country|source|years
Afghanistan|DHS|2015
Albania|DHS|2008-2017
Algeria|MICS|2012-2019
Angola|Census|2014
Argentina|IPUMS|1970-2001
Armenia|IPUMS|2001-2011
Australia|Panel Data|2001-2019
Austria|Panel Data|1981-2017
Bangladesh|IPUMS|1991-2011
Belarus|IPUMS|1999-2009
Belgium|LIS|1985-2017
Benin|IPUMS|1979-2013
Bolivia|IPUMS|1976-2001
Botswana|IPUMS|1991-2011
Brazil|IPUMS|1991-2010
Bulgaria|SILC|2007-2020
Burkina Faso|IPUMS|1996-2006
Burundi|DHS|2010-2016
Cambodia|IPUMS|1998-2008
Cameroon|IPUMS|1976-2005
Canada|IPUMS|2011
Chad|DHS|1996-2014
Chile|LIS|1990-2017
China|Mini Census|2005
Colombia|IPUMS|1973-2005
Congo, Dem. Rep.|DHS|2007-2013
Congo, Rep.|DHS|2005-2011
Costa Rica|IPUMS|1973-2011
Cote d'Ivoire|LIS|2002-2015
Croatia|SILC|2010-2020
Cuba|IPUMS|2002-2012
Cyprus|SILC|2009-2020
Czech Republic|LIS|1992-2016
Denmark|Panel Data|1980-2017
Dominican Rep.|IPUMS|1981-2010
Ecuador|IPUMS|1982-2010
Egypt|IPUMS|1996-2006
El Salvador|IPUMS|1992-2007
Estonia|SILC|2009-2020
Ethiopia|DHS|2000-2016
Fiji|IPUMS|1976-2014
Finland|SILC|2009-2020
France|LFS|1990-2020
Gabon|DHS|2000-2012
Gambia|LFS|2010-2015
Georgia|LFS|2020-2021
Germany|LIS|1989-2005
Ghana|IPUMS|2000-2010
Greece|IPUMS|1981-2001
Guatemala|IPUMS|1964-2002
Guinea|IPUMS|1983-2014
Guyana|LFS|2017-2021
Haiti|IPUMS|1971-2003
Honduras|IPUMS|1974-2001
Hungary|IPUMS|1990-2011
Iceland|SILC|2004-2018
India|DHS|2005-2015
Indonesia|IPUMS|1971-2010
Iran|IPUMS|2006
Iraq|IPUMS|1997
Ireland|LIS|1994-2018
Israel|LIS|1986-2013
Italy|LIS|1986-2020
Jamaica|IPUMS|1982-2001
Japan|Panel Data|2004-2020
Jordan|IPUMS|2004
Kenya|IPUMS|1989-2009
Kyrgyz Republic|DHS|1997-2012
Laos|IPUMS|2005
Latvia|SILC|2009-2020
Lesotho|IPUMS|1996-2006
Liberia|IPUMS|2008
Lithuania|LIS|2009-2018
Luxembourg|LIS|1985-2013
Madagascar|DHS|1992-2008
Malawi|IPUMS|1987-2008
Malaysia|IPUMS|1991-2000
Maldives|DHS|2009-2016
Mali|IPUMS|1987-2009
Mauritius|IPUMS|1990-2011
Mexico|IPUMS|1970-2015
Moldova|DHS|2005
Mongolia|IPUMS|2000
Morocco|IPUMS|1982-2004
Mozambique|IPUMS|1997-2007
Myanmar|IPUMS|2014
Namibia|DHS|1992-2013
Nepal|IPUMS|2001-2011
Netherlands|LIS|1990-2018
Nicaragua|IPUMS|1995-2005
Niger|DHS|1992-2012
Nigeria|DHS|1990-2018
Norway|Panel Data|1993-2017
Pakistan|LFS|2010-2021
Panama|IPUMS|1960-2010
Papua New Guinea|IPUMS|1980-2000
Paraguay|IPUMS|1962-2002
Peru|IPUMS|1993-2007
Philippines|IPUMS|1990
Poland|LIS|1992-2020
Portugal|IPUMS|1981-2011
Puerto Rico|IPUMS|1990-2010
Romania|IPUMS|1992-2011
Russia|IPUMS|2002-2010
Rwanda|IPUMS|2002-2012
Senegal|IPUMS|1988-2002
Serbia|SILC|2013-2020
Sierra Leone|IPUMS|2004
Slovakia|LIS|1992-2018
Slovenia|LIS|1997-2012
South Africa|IPUMS|1996-2011
South Korea|LIS|2006-2016
South Sudan|IPUMS|2008
Spain|IPUMS|1991-2001
Sudan|IPUMS|2008
Suriname|IPUMS|2012
Sweden|Panel Data|1997-2017
Switzerland|Panel Data|1981-2020
Taiwan|LIS|1981-2016
Tanzania|IPUMS|1988-2012
Thailand|IPUMS|1990-2000
Timor-Leste|DHS|2009-2016
Togo|IPUMS|2010
Trinidad & Tobago|IPUMS|1970-2011
Tunisia|Census|2004
Turkey|IPUMS|1985-2000
Uganda|IPUMS|1991-2014
United Kingdom|APS|2012-2020
United States|CPS/ACS|1968-2020
Uruguay|IPUMS|1963-2011
Venezuela|IPUMS|1971-2001
Vietnam|IPUMS|1989-2009
Zambia|IPUMS|1990-2010
Zimbabwe|IPUMS|2012"

atlas_audit <- read.delim(
  textConnection(atlas_source_text), sep = "|", quote = "",
  stringsAsFactors = FALSE, check.names = FALSE
)
stopifnot(nrow(atlas_audit) == 134)

ipums_pay_yes <- c("Brazil", "Canada", "Jamaica", "Mexico", "Panama", "Puerto Rico", "Venezuela")
ipums_household_only <- "South Africa"

atlas_audit <- atlas_audit |>
  mutate(
    public_access = case_when(
      source == "IPUMS" ~ "Yes: free registration and approved extract",
      source == "DHS" ~ "Yes: free registration and project request",
      source == "MICS" ~ "Yes: free registration/download",
      source == "SILC" ~ "Restricted: Eurostat scientific-use application",
      source == "LIS" ~ "Restricted: LIS membership/remote execution",
      source == "Panel Data" ~ "Restricted; country-specific application",
      source == "LFS" ~ "Varies by national statistical office; verify access",
      source == "APS" ~ "Yes: UK Data Service registration/licence",
      source == "CPS/ACS" ~ "Yes: public-use microdata",
      source == "Mini Census" ~ "Restricted; public microdata not confirmed",
      source == "Census" ~ "Not confirmed from the paper"
    ),
    pay_gap_data = case_when(
      country %in% ipums_pay_yes ~ "Yes/partial: own earnings plus hours/full-time route; verify wave-specific covariates",
      country %in% ipums_household_only ~ "No: household income only; own earnings absent in cited IPUMS censuses",
      source == "IPUMS" ~ "Not confirmed: Atlas employment file lacks a verified own-earnings + hours/full-time combination",
      source %in% c("DHS", "MICS") ~ "No: lacks monetary own earnings and hours; industry/occupation detail is also insufficient",
      source == "SILC" ~ "Yes: earnings, hours/full-time, education, occupation, industry, and household structure; access restricted",
      source == "LIS" ~ "Partial/likely: labor earnings and demographics exist; hours and industry/occupation vary by country-wave",
      source == "Panel Data" ~ "Yes for annual-earnings gaps; exact hours/occupation/industry vary and administrative access is restricted",
      source == "LFS" ~ "Unconfirmed: hours, occupation, industry, and education exist, but earnings and public access vary",
      source == "APS" ~ "Yes: earnings, paid hours, education, occupation, industry, and family variables",
      source == "CPS/ACS" ~ "Yes: earnings, hours, education, occupation, industry, and household parent links",
      source == "Mini Census" ~ "Unconfirmed: access and individual earnings documentation not established",
      source == "Census" ~ "Unconfirmed: Atlas documents employment, not a public own-earnings/hourly-wage file"
    )
  )

write.csv(atlas_audit, "results/child_penalty_atlas_data_audit.csv", row.names = FALSE)

atlas_summary <- atlas_audit |>
  mutate(category = case_when(
    grepl("^Yes", pay_gap_data) ~ "Yes",
    grepl("^(Partial|Unconfirmed|Not confirmed)", pay_gap_data) ~ "Partial or unconfirmed",
    TRUE ~ "No"
  )) |>
  count(category, name = "countries")

kable(
  atlas_summary,
  col.names = c("Assessment", "Number of countries"),
  caption = "Summary of the 134-country Atlas audit"
)

kable(
  atlas_audit,
  col.names = c("Country", "Atlas data source", "Atlas years", "Public/research access", "Gender-pay-gap readiness"),
  caption = "Country-level review of the data sources in Kleven et al. (2025), Appendix Table A.1"
)
```

The Atlas table is a source audit, not a guarantee of variable access. In particular, “Panel Data” often means confidential administrative registers; LIS provides harmonized output through a remote system rather than distributing national microdata; EU-SILC scientific-use files require approval; and national labor-force-survey earnings modules vary. The machine-readable audit is available in `results/child_penalty_atlas_data_audit.csv`.

Access documentation: [IPUMS International](https://international.ipums.org/international-action/faq), [DHS Program](https://dhsprogram.com/data/Access-Instructions.cfm), [UNICEF MICS](https://mics.unicef.org/surveys), [Luxembourg Income Study](https://www.lisdatacenter.org/data-access/), [Eurostat microdata access](https://ec.europa.eu/eurostat/web/microdata/overview), [UK Annual Population Survey](https://beta.ukdataservice.ac.uk/datacatalogue/series/series?id=200002), and [U.S. CPS](https://www.census.gov/programs-surveys/cps/data.html).

Five comparability cautions apply even to feasible settings. First, co-resident parenthood misses nonresident and older children, while fertility-history parenthood is generally women-only. Second, `INCWAGE` and `INCEARN` are different concepts; the latter includes business and farm earnings. Third, income can be reported weekly, monthly, or annually and is subject to sample-specific universes and top-codes; numerator and hours must be placed on a common period within each sample. Fourth, `HRSFULL` definitions vary across countries and waves. Fifth, broad industry and occupation harmonization can still contain source-classification breaks. A real analysis should read every sample-specific comparability tab and apply IPUMS person weights.

Sources: [sample catalogue](https://international.ipums.org/international-action/samples) and country sample-detail pages; variable documentation for [INCWAGE](https://international.ipums.org/international-action/variables/INCWAGE), [INCEARN](https://international.ipums.org/international-action/variables/INCEARN), [HRSWORK1](https://international.ipums.org/international-action/variables/HRSWORK1), [HRSUSUAL1](https://international.ipums.org/international-action/variables/HRSUSUAL1), [HRSMAIN](https://international.ipums.org/international-action/variables/HRSMAIN), [HRSACTUAL1](https://international.ipums.org/international-action/variables/HRSACTUAL1), [HRSFULL](https://international.ipums.org/international-action/variables/HRSFULL), [NCHILD](https://international.ipums.org/international-action/variables/NCHILD), [MOMLOC](https://international.ipums.org/international-action/variables/MOMLOC), [POPLOC](https://international.ipums.org/international-action/variables/POPLOC), [CHBORN](https://international.ipums.org/international-action/variables/CHBORN), [CHSURV](https://international.ipums.org/international-action/variables/CHSURV), [LASTBYR](https://international.ipums.org/international-action/variables/LASTBYR), [INDGEN](https://international.ipums.org/international-action/variables/INDGEN), and [OCCISCO](https://international.ipums.org/international-action/variables/OCCISCO).

## Conclusions

- The synthetic mean log-wage gap narrows over the three censuses.
- Parenthood remains associated with larger gaps, while selection into employment makes the marriage gradient vary across waves.
- Women’s educational gains remove education as an explanation and can make its contribution negative.
- Occupational and industrial sorting remain meaningful explained components even as segregation declines.
- Both changing characteristics and a declining unexplained component contribute to convergence from 1995.
- Distributional decompositions reveal that a single mean gap can conceal different composition and coefficient effects across the wage distribution.

## Reproduction and files

The report source is `krabbistan-gender-gap.qmd`. Running

```bash
quarto render krabbistan-gender-gap.qmd
```

regenerates the census files, all CSV result tables, figures, and this self-contained HTML. The optional `code/krabbistan_gender_gap.do` contains a standalone Stata translation of the simulation and mean/distributional decomposition workflow.

Downloads:

- [Combined census, Stata format](data/krabbistan_census_1995_2015.dta)
- [Combined census, CSV format](data/krabbistan_census_1995_2015.csv)
- [1995 Stata file](data/krabbistan_census_1995.dta), [2005 Stata file](data/krabbistan_census_2005.dta), and [2015 Stata file](data/krabbistan_census_2015.dta)
- [R/Quarto source](krabbistan-gender-gap.qmd)
- [Optional Stata translation](code/krabbistan_gender_gap.do)
- [Machine-readable result tables](results/)

## Reference

Blau, Francine D., and Lawrence M. Kahn. 2017. “The Gender Wage Gap: Extent, Trends, and Explanations.” *Journal of Economic Literature* 55 (3): 789–865. [https://doi.org/10.1257/jel.20160995](https://doi.org/10.1257/jel.20160995)

Kleven, Henrik, Camille Landais, and Gabriel Leite-Mariante. 2025. “The Child Penalty Atlas.” *Review of Economic Studies* 92 (5): 3174–3207. [https://doi.org/10.1093/restud/rdae104](https://doi.org/10.1093/restud/rdae104). Country sources are transcribed from Appendix Table A.1.