1 Introduction

What follows is a revision of my thesis research project, Do Poor Countries Catch Up to Rich Countries? Structural Change in the World Economy, 1816-1916. It is a Markov analysis of country-level social mobility with inputs upper class, middle class, and lower class derived from frequency distributions. Histograms were plotted using national primary energy consumption weighted by population with a density overlay. Raw data were taken from the National Materiel Capabilities data set (version 6) hosted by the Correlates of War Project. Primary energy consumption was used as a proxy for industrial output or level of industrialization. Minima in the density overlay were interpreted as class boundaries. Bandwidth was methodologically determined. This approach produced three classes between 1860 and 1880 and four classes between 1881 and 1916. The fourth class in the second era was the result of non-industrial countries developing industrial capabilities and breaking into those markets. It can be seen as a bifurcation of the lower class and an adaptation of its role as the source of cheap labor in the evolving world economy. As the lower class became increasingly industrialized, primary energy consumption was no longer a distinguishing feature of class, and after 1916, this approach to operationalizing class structure was unreliable.

The Markov analysis assessed social mobility, answering the question: Do poor countries catch up to rich countries? The question was examined with a balanced panel of countries composed of five-year snapshots between 1860 and 1916. The findings indicate a persistent class structure with upper class being an absorbing state in both eras. Structural change was greater in the second era, which is to be expected in moments of hegemonic transition. This can be seen in transitions from lower class 1 to lower class 2 (the non-industrial lower class to the industrializing lower class), and in transitions from middle class to lower class 2, with lower class 2 being an absorbing state. Thus, structural change in the second era consisted of both upward and downward movement in the economic hierarchy.

An evaluation of Markovian memory rejected independence in favor of first-order dependence, indicating strong state dependence–a country’s position at one time point strongly predicts its position at the next. The second order test failed to reject first order sufficiency, not because higher order dependence was absent, but because the high degree of persistence rendered the test uninformative. Over the long interval, Pearson and likelihood-ratio chi-square tests strongly rejected independence, implying substantial long-run persistence in countries’ hierarchical positions. The Markov chain failed the test of approximate ergodicity. That is to say, the long-term distribution of countries in the upper, middle, and lower classes was not random. A country’s starting position had a lasting effect for the period of examination.

These findings support the contention that economic development, while reducing poverty and improving living conditions, does not change the relative differences of countries in the global economic hierarchy. In this sense, findings undermine the central claim of economic convergence, namely that poor countries are catching up to rich countries.

# Preferences
options(scipen = 999,
        digits = 2)
# Package Dependencies
library(dplyr)
library(tidyr)
library(ggplot2)
library(MASS)
library(parallel)
library(readxl)
library(knitr)

# Data Import and Recoding
df <- read_excel("Class Assignments.xlsx") %>% 
  dplyr::select(stateabb, year, class) %>%
  dplyr::mutate(class = dplyr::recode(class,
                                      "core"          = "upper_class",
                                      "semiperiphery" = "middle_class",
                                      "periphery_1"   = "lower_class_1",
                                      "periphery_2"   = "lower_class_2"))

2 Transition Sequence Construction

build_transitions <- function(data,
                                       snapshot_years,
                                       step,
                                       full_levels) {
  df_snap <- data %>%
    dplyr::filter(year %in% snapshot_years) %>%
    dplyr::mutate(class = factor(class, levels = full_levels))
  df_balanced <- df_snap %>%
    arrange(stateabb, year) %>%
    group_by(stateabb) %>%
    mutate(n_obs = n()) %>%
    ungroup()
  message("Balanced-panel countries for snapshots {",
          paste(snapshot_years, collapse = ", "),
          "}: ", paste("n =",
          length(unique(df_balanced$stateabb))))
  df_trans <- df_balanced %>%
    dplyr::arrange(stateabb, year) %>%
    dplyr::group_by(stateabb) %>%
    dplyr::mutate(
      class_prev = dplyr::lag(class),
      year_prev  = dplyr::lag(year)) %>%
    dplyr::ungroup() %>%
    dplyr::filter(
      !is.na(class_prev),
      (year - year_prev) %in% step) %>%
    dplyr::mutate(
      step_id = factor(paste0(year_prev, "→", year)))
  df_trans
}

# ERA 1: 1860-1880
# 3-class system: upper_class, middle_class, lower_class_1
# Snapshots: 1860, 1865, 1870, 1875, 1880 (5-year steps)

# ---- Subset data ----
era1_years <- 1860:1880
df1 <- df %>% dplyr::filter(year %in% era1_years)

levels_era1 <- c("upper_class", "middle_class", "lower_class_1")

snap1 <- c(1860, 1865, 1870, 1875, 1880)
step1 <- 5

# ---- Build transitions ----
trans1 <- build_transitions(
  data           = df1,
  snapshot_years = snap1,
  step           = step1,
  full_levels    = levels_era1)

# ERA 2: (1881-1916)
# 3-class system: upper_class, middle_class, lower_class_2, lower_class_1
# Snapshots: 1881, 1886, 1891, 1896, 1901, 1906, 1911, 1916 (5-year steps)

# ---- Subset data ----
era2_years <- 1881:1916
df2 <- df %>% dplyr::filter(year %in% era2_years)

levels_era2 <- c("upper_class", "middle_class", "lower_class_2", "lower_class_1")

snap2 <- c(1881, 1886, 1891, 1896, 1901, 1906, 1911, 1916)
step2 <- 5

# ---- Build transitions ----
trans2 <- build_transitions(
  data           = df2,
  snapshot_years = snap2,
  step           = step2,
  full_levels    = levels_era2)

# Dataset Overview
era_summary <- tibble::tibble(
  Era = c("ERA 1 (1860-1880)", "ERA 2 (1881-1916)"),
  Snapshots = c(
  paste(snap1, collapse = ", "),
  paste(snap2, collapse = ", ")),
  Transitions = c(
  nrow(trans1),
  nrow(trans2)))

cat('Balanced-panel snapshot coverage and transition counts')
era_summary_ordered <- era_summary %>% 
  dplyr::select(Era, Transitions, Snapshots)

kable(era_summary_ordered,
      align = c("l", "l", "r"))

Balanced-panel snapshot coverage and transition counts

Era	Transitions	Snapshots
ERA 1 (1860-1880)	138	1860, 1865, 1870, 1875, 1880
ERA 2 (1881-1916)	274	1881, 1886, 1891, 1896, 1901, 1906, 1911, 1916

3 Sample Attrition Diagnostics

attrition_diagnostic <- function(data, 
                                 snapshot_years, 
                                 balanced_df, 
                                 full_levels) {
  kept_states <- unique(balanced_df$stateabb)
  year0 <- min(snapshot_years)
  all_year0 <- data %>%
    dplyr::filter(year == year0) %>%
    dplyr::mutate(class = factor(class, levels = full_levels)) %>%
    dplyr::count(class, name = "n") %>%
    dplyr::mutate(group = "All")
  kept_year0 <- data %>%
    dplyr::filter(year == year0, stateabb %in% kept_states) %>%
    dplyr::mutate(class = factor(class, levels = full_levels)) %>%
    dplyr::count(class, name = "n") %>%
    dplyr::mutate(group = "Balanced")
  df_comp <- dplyr::bind_rows(all_year0, kept_year0)
  tab <- xtabs(n ~ group + class, data = df_comp)
  chisq <- suppressWarnings(chisq.test(tab))
  list(counts = df_comp,
       table = tab,
       chisq = chisq)
}

interpret_attrition <- function(chisq_obj, alpha = 0.05) {
  p <- chisq_obj$p.value
  if (p < alpha) {
    paste0("Reject equality of baseline class distributions (p < ", alpha,
           "). Evidence of selective attrition.")
  } else {
    paste0("Fail to reject equality of baseline class distributions (p ≥ ", alpha,
           "). No evidence of selective attrition.")
  }
}

print_chisq_compact <- function(chisq_obj) {
  cat("X^2 =", round(chisq_obj$statistic, 3),
      " df =", chisq_obj$parameter,
      " p =", signif(chisq_obj$p.value, 4), "\n")
}

3.1 ERA 1 1860-1880

# ---- Attrition diagnostic (ERA 1) ----
df1_snap <- df1 %>% dplyr::filter(year %in% snap1)
df1_balanced <- df1_snap %>%
  dplyr::group_by(stateabb) %>%
  dplyr::filter(dplyr::n() == length(snap1)) %>%
  dplyr::ungroup()

attr1 <- attrition_diagnostic(df1, snap1, df1_balanced, levels_era1)

cat("\n----------------------------------------\n")
cat("ERA 1 Attrition diagnostic\n")
cat("----------------------------------------\n")

cat("\nTest statistics:\n")
print_chisq_compact(attr1$chisq)

cat("\nInterpretation:\n")
cat(interpret_attrition(attr1$chisq), "\n")

cat("\nObserved counts:\n")
print(attr1$counts)


----------------------------------------
ERA 1 Attrition diagnostic
----------------------------------------

Test statistics:
X^2 = 0.097  df = 2  p = 0.95 

Interpretation:
Fail to reject equality of baseline class distributions (p ≥ 0.05). No evidence of selective attrition. 

Observed counts:
# A tibble: 6 × 3
  class             n group   
  <fct>         <int> <chr>   
1 upper_class       8 All     
2 middle_class     16 All     
3 lower_class_1    13 All     
4 upper_class       6 Balanced
5 middle_class     14 Balanced
6 lower_class_1    12 Balanced

3.2 ERA 2 1881-1916

# ---- Attrition diagnostic (ERA 2) ----
df2_snap <- df2 %>% dplyr::filter(year %in% snap2)
df2_balanced <- df2_snap %>%
  dplyr::group_by(stateabb) %>%
  dplyr::filter(dplyr::n() == length(snap2)) %>%
  dplyr::ungroup()

attr2 <- attrition_diagnostic(df2, snap2, df2_balanced, levels_era2)

cat("\n----------------------------------------\n")
cat("ERA 2 Attrition diagnostic\n")
cat("----------------------------------------\n")

cat("\nTest statistics:\n")
print_chisq_compact(attr2$chisq)

cat("\nInterpretation:\n")
cat(interpret_attrition(attr2$chisq), "\n")

cat("\nObserved counts:\n")
print(attr2$counts)


----------------------------------------
ERA 2 Attrition diagnostic
----------------------------------------

Test statistics:
X^2 = 0.5  df = 3  p = 0.92 

Interpretation:
Fail to reject equality of baseline class distributions (p ≥ 0.05). No evidence of selective attrition. 

Observed counts:
# A tibble: 8 × 3
  class             n group   
  <fct>         <int> <chr>   
1 upper_class       7 All     
2 middle_class     16 All     
3 lower_class_2     1 All     
4 lower_class_1    13 All     
5 upper_class       7 Balanced
6 middle_class     16 Balanced
7 lower_class_2     1 Balanced
8 lower_class_1     9 Balanced

4 5-Year Step Transition Matrices

# ---- Heatmap Function (Probabilities Only) ----
plot_prob_heatmap <- function(P, title, subtitle = NULL) {
  df_plot <- as.data.frame(as.table(P))
  colnames(df_plot) <- c("origin", "dest", "prob")
  df_plot$origin <- factor(df_plot$origin, levels = rev(rownames(P)))
  df_plot$dest   <- factor(df_plot$dest,   levels = colnames(P))
  ggplot2::ggplot(df_plot, ggplot2::aes(x = dest, y = origin, fill = prob)) +
    ggplot2::geom_tile(color = "grey70", linewidth = 0.3) +
    ggplot2::geom_text(ggplot2::aes(label = sprintf("%.2f", prob)), size = 4) +
    ggplot2::scale_fill_gradient(low = "white", high = "steelblue") +
    ggplot2::scale_x_discrete(position = "top") +
    ggplot2::labs(
      title = title,
      subtitle = subtitle,
      x = "Destination class",
      y = "Origin class",
      fill = "Prob"
    ) +
    ggplot2::theme_minimal(base_size = 13) +
    ggplot2::theme(
      axis.text.x.top = ggplot2::element_text(margin = ggplot2::margin(b = 10)),
      axis.title.x.top = ggplot2::element_text(margin = ggplot2::margin(b = 20)),
      axis.title.y = ggplot2::element_text(margin = ggplot2::margin(r = 20)),
      axis.ticks = ggplot2::element_line(color = "black", linewidth = 0.4),
      panel.grid = ggplot2::element_blank()
    ) +
    ggplot2::coord_fixed(ratio = 0.8)
}

# ERA 1 COMPUTATIONS (5-Year)
# Counts
C1 <- xtabs(~ class_prev + class, data = trans1)
C1 <- C1[levels_era1, levels_era1, drop = FALSE]
C1_df <- as.data.frame.matrix(C1) 

# Probabilities
P1 <- prop.table(C1, margin = 1)
P1_df <- as.data.frame.matrix(round(P1, 3))

# Plot
p1_probs <- plot_prob_heatmap(P1, "ERA 1 Probabilities (1860-1880)", "5-year steps")

# ERA 2 COMPUTATIONS (5-Year)
# Counts
C2 <- xtabs(~ class_prev + class, data = trans2)
C2 <- C2[levels_era2, levels_era2, drop = FALSE]
C2_df <- as.data.frame.matrix(C2)

# Probabilities
P2 <- prop.table(C2, margin = 1)
P2_df <- as.data.frame.matrix(round(P2, 3))

# Plot
p2_probs <- plot_prob_heatmap(P2, "ERA 2 Probabilities (1881-1916)", "5-year steps")

4.1 ERA 1 1860-1880

# Print Count Matrix
cat("Transition Counts")
knitr::kable(C1_df, caption = "Era 1 Counts")

# Print Probability Matrix
#cat("Transition Probabilities")
#knitr::kable(P1_df, caption = "Era 1 Probabilities")

# Plot Probability Heatmap
print(p1_probs)

Transition Counts

Era 1 Counts
	upper_class	middle_class	lower_class_1
upper_class	28	1	0
middle_class	1	57	0
lower_class_1	0	1	50

4.2 ERA 2 1881-1916

# Print Count Matrix
cat("<strong>Transition Counts</strong>")
knitr::kable(C2_df, caption = "Era 2 Counts")

# Print Probability Matrix
#cat("<br><strong>Transition Probabilities</strong>")
#knitr::kable(P2_df, caption = "Era 2 Probabilities")

# Plot Probability Heatmap
print(p2_probs)

<strong>Transition Counts</strong>

Era 2 Counts
	upper_class	middle_class	lower_class_2	lower_class_1
upper_class	36	3	0	0
middle_class	3	100	8	0
lower_class_2	0	3	41	0
lower_class_1	0	0	6	74

5 5-Year Step Time Homogeneity Tests

markov_stationarity_test_dezzani <- function(transitions_df) {

  # Build C×C×T counts: step × origin × dest
  counts <- transitions_df %>%
    dplyr::count(step_id, class_prev, class, name = "n") %>%
    tidyr::complete(step_id, class_prev, class, fill = list(n = 0)) %>%
    dplyr::mutate(
      step_id    = factor(step_id),
      class_prev = factor(class_prev),
      class      = factor(class)
    )

  df_counts <- as.data.frame(counts)
  colnames(df_counts) <- c("step", "origin", "dest", "n")

  # Dezzani-style stationarity / time-homogeneity null:
  # dest ⟂ step | origin  <=>  n ~ origin*dest + origin*step
  fit_H0 <- MASS::loglm(
    n ~ origin*dest + origin*step,
    data = df_counts
  )

  # Saturated model (allows transition matrix to vary by step):
  fit_H1 <- MASS::loglm(
    n ~ origin*dest*step,
    data = df_counts
  )

  # Likelihood-ratio test
  G2      <- fit_H0$deviance - fit_H1$deviance
  df_diff <- fit_H0$df - fit_H1$df
  p_val   <- 1 - pchisq(G2, df_diff)

  list(
    counts = df_counts,
    G2 = G2,
    df = df_diff,
    p_value = p_val,
    fit_H0 = fit_H0,
    fit_H1 = fit_H1
  )
}

interpret_stationarity_dezzani <- function(p, alpha = 0.05) {
  if (is.na(p)) {
    return("p-value is NA (often due to sparse/degenerate cells); Stationarity cannot be assessed reliably with this table.")
  }
  if (p < alpha) {
    paste0("Reject stationarity/time-homogeneity (p < ", alpha,
           "): The 5-year transition probabilities vary across steps.")
  } else {
    paste0("Fail to reject stationarity/time-homogeneity (p ≥ ", alpha,
           "): No evidence that the 5-year transition probabilities vary across steps.")
  }
}

5.1 ERA 1 1860-1880

# ERA 1: 1860-1880 ---- Homogeneity test
test1_homog <- markov_stationarity_test_dezzani(trans1)

cat("\nERA 1 Stationarity / Time-Homogeneity (Dezzani-style)\n")
cat("G^2 =", round(test1_homog$G2, 3),
    " df =", test1_homog$df,
    " p =", signif(test1_homog$p_value, 4), "\n")
cat("Interpretation:", interpret_stationarity_dezzani(test1_homog$p_value), "\n")


ERA 1 Stationarity / Time-Homogeneity (Dezzani-style)
G^2 = 8.4  df = 18  p = 0.97 
Interpretation: Fail to reject stationarity/time-homogeneity (p ≥ 0.05): No evidence that the 5-year transition probabilities vary across steps.

5.2 ERA 2 1881-1916

# ERA 2: 1881-1916 ---- Homogeneity test
test1_homog <- markov_stationarity_test_dezzani(trans1)

cat("\nERA 1 Stationarity / Time-Homogeneity (Dezzani-style)\n")
cat("G^2 =", round(test1_homog$G2, 3),
    " df =", test1_homog$df,
    " p =", signif(test1_homog$p_value, 4), "\n")
cat("Interpretation:", interpret_stationarity_dezzani(test1_homog$p_value), "\n")


ERA 1 Stationarity / Time-Homogeneity (Dezzani-style)
G^2 = 8.4  df = 18  p = 0.97 
Interpretation: Fail to reject stationarity/time-homogeneity (p ≥ 0.05): No evidence that the 5-year transition probabilities vary across steps.

6 Long Interval Transition Matrices

# Function for computing long-interval transitions
compute_long_trans <- function(data, start_year, end_year, full_levels) {
  df_start <- data %>%
    dplyr::filter(year == start_year) %>%
    dplyr::select(stateabb, class_start = class)
  
  df_end <- data %>%
    dplyr::filter(year == end_year) %>%
    dplyr::select(stateabb, class_end = class)
  
  df_long <- dplyr::inner_join(df_start, df_end, by = "stateabb") %>%
    dplyr::mutate(
      class_start = factor(class_start, levels = full_levels),
      class_end   = factor(class_end,   levels = full_levels)) %>%
    dplyr::filter(!is.na(class_start), !is.na(class_end))
  
  tab <- with(df_long, table(class_start, class_end))
  list(df_long = df_long, tab = tab)
}

long_interval_probs <- function(tab) {
  pt <- prop.table(tab, margin = 1)
  pt[is.nan(pt)] <- 0 
  return(pt)
}

# ERA 1 COMPUTATION
long1 <- compute_long_trans(df1, 1860, 1880, levels_era1)
long1$tab <- long1$tab[levels_era1, levels_era1, drop = FALSE]

# Create Data Frames for HTML printing
long1_counts_df <- as.data.frame.matrix(long1$tab)
P_long1 <- long_interval_probs(long1$tab)
long1_probs_df  <- as.data.frame.matrix(round(P_long1, 3))

# Generate Plot
p_long1_probs <- plot_prob_heatmap(P_long1, "ERA 1 Long Interval (1860-1880)", "20-year transition")

# ERA 2 COMPUTATION
# Note: Using 1916 to prevent blank matrix
long2 <- compute_long_trans(df2, 1881, 1916, levels_era2)
long2$tab <- long2$tab[levels_era2, levels_era2, drop = FALSE]

# Create Data Frames for HTML printing
long2_counts_df <- as.data.frame.matrix(long2$tab)
P_long2 <- long_interval_probs(long2$tab)
long2_probs_df  <- as.data.frame.matrix(round(P_long2, 3))

# Generate Plot
p_long2_probs <- plot_prob_heatmap(P_long2, "ERA 2 Long Interval (1881-1916)", "35-year transition")

6.1 ERA 1 1860-1880

# Print Counts
cat("Transition Counts")
knitr::kable(long1_counts_df, caption = "Era 1 Long Interval Counts")

# Print Probs
#cat("<Transition Probabilities")
#knitr::kable(long1_probs_df, caption = "Era 1 Long Interval Probabilities")

# Print Map
print(p_long1_probs)

Transition Counts

Era 1 Long Interval Counts
	upper_class	middle_class	lower_class_1
upper_class	6	0	0
middle_class	0	14	0
lower_class_1	0	1	12

6.2 ERA 2 1881-1916

# Print Counts
cat("Transition Counts")
knitr::kable(long2_counts_df, caption = "Era 2 Long Interval Counts")

# Print Probs
#cat("Transition Probabilities")
#knitr::kable(long2_probs_df, caption = "Era 2 Long Interval Probabilities")

# Print Map
print(p_long2_probs)

Transition Counts

Era 2 Long Interval Counts
	upper_class	middle_class	lower_class_2	lower_class_1
upper_class	7	0	0	0
middle_class	0	11	5	0
lower_class_2	0	0	1	0
lower_class_1	0	0	5	4

7 Long-Interval Transition Diagnostics

# Bowker symmetry test
bowker_symmetry_safe <- function(tab) {
  tab <- as.matrix(tab)
  storage.mode(tab) <- "numeric"
  if (nrow(tab) != ncol(tab)) stop("tab must be square")
  k <- nrow(tab)

  X2 <- 0
  df <- 0
  for (i in 1:(k - 1)) {
    for (j in (i + 1):k) {
      nij <- tab[i, j]; nji <- tab[j, i]
      s <- nij + nji
      if (s > 0) {
        X2 <- X2 + (nij - nji)^2 / s
        df <- df + 1
      }
    }
  }
  p <- if (df > 0) 1 - pchisq(X2, df) else NA_real_
  list(X2 = X2, df = df, p = p)
}

# Bowker X^2, Pearson X^2 & Likelihood Ratio G^2 long-interval tests
db_long_interval_tests <- function(tab) {
  pear <- suppressWarnings(stats::chisq.test(tab, correct = FALSE))

  df_counts <- as.data.frame(tab)
  colnames(df_counts) <- c("origin", "dest", "n")
  fit_ind <- MASS::loglm(n ~ origin + dest, data = df_counts)
  LR_G2   <- fit_ind$deviance
  LR_df   <- fit_ind$df
  LR_p    <- 1 - pchisq(LR_G2, LR_df)

  bow <- bowker_symmetry_safe(tab)

  list(
    Pearson = list(X2 = unname(pear$statistic), df = unname(pear$parameter), p = pear$p.value),
    LR      = list(G2 = LR_G2,                  df = LR_df,                  p = LR_p),
    Bowker  = list(X2 = bow$X2,                 df = bow$df,                 p = bow$p)
  )
}

interpret_long_association <- function(p, alpha = 0.05) {
  if (is.na(p)) return("Not estimable.")
  if (p < alpha) "Reject independence. Strong long-run persistence." else "Fail to reject independence."
}

interpret_bowker <- function(p, df, alpha = 0.05) {
  if (is.na(p)) return("Not estimable.")
  if (p < alpha) {
    paste0("Reject symmetry. Directional/asymmetric movement (df=", df, ").")
  } else {
    paste0("Fail to reject symmetry (df=", df, ").")
  }
}

p_fmt <- function(p) {
  if (is.na(p)) return("NA")
  if (p < 0.01) return("<0.01")
  formatC(p, format = "f", digits = 2)
}

# Output function (P_long removed because unused)
print_long_interval_block <- function(label, start_year, end_year, long_obj, alpha = 0.05) {
  tests <- db_long_interval_tests(long_obj$tab)
  cat("\n----------------------------------------\n")
  cat(label, "Long-Interval Diagnostics (", start_year, " -> ", end_year, ")\n", sep = " ")
  cat("----------------------------------------\n")
  cat("\nTest statistics:\n")
  cat("Pearson Chi-square:          X^2 =", round(tests$Pearson$X2, 2),
      " df =", tests$Pearson$df,
      " p =", p_fmt(tests$Pearson$p), "\n")
  cat("Likelihood Ratio Chi-square: G^2 =", round(tests$LR$G2, 2),
      " df =", tests$LR$df,
      " p =", p_fmt(tests$LR$p), "\n")
  cat("Bowker Symmetry Chi-square:  X^2 =", round(tests$Bowker$X2, 2),
      " df =", tests$Bowker$df,
      " p =", p_fmt(tests$Bowker$p), "\n")
  cat("\nInterpretation:\n")
  cat("Pearson:", interpret_long_association(tests$Pearson$p, alpha = alpha), "\n")
  cat("Likelihood Ratio:", interpret_long_association(tests$LR$p, alpha = alpha), "\n")
  cat("Bowker:", interpret_bowker(tests$Bowker$p, df = tests$Bowker$df, alpha = alpha), "\n")

  invisible(list(tests = tests))
}

7.1 ERA 1 1860-1880

print_long_interval_block("ERA 1", 1860, 1880, long1, alpha = 0.05)


----------------------------------------
ERA 1 Long-Interval Diagnostics ( 1860  ->  1880 )
----------------------------------------

Test statistics:
Pearson Chi-square:          X^2 = 61  df = 4  p = <0.01 
Likelihood Ratio Chi-square: G^2 = 61  df = 4  p = <0.01 
Bowker Symmetry Chi-square:  X^2 = 1  df = 1  p = 0.32 

Interpretation:
Pearson: Reject independence. Strong long-run persistence. 
Likelihood Ratio: Reject independence. Strong long-run persistence. 
Bowker: Fail to reject symmetry (df=1).

7.2 ERA 2 1881-1916

print_long_interval_block("ERA 2", 1881, 1916, long2, alpha = 0.05)


----------------------------------------
ERA 2 Long-Interval Diagnostics ( 1881  ->  1916 )
----------------------------------------

Test statistics:
Pearson Chi-square:          X^2 = 53  df = 9  p = <0.01 
Likelihood Ratio Chi-square: G^2 = 55  df = 9  p = <0.01 
Bowker Symmetry Chi-square:  X^2 = 10  df = 2  p = <0.01 

Interpretation:
Pearson: Reject independence. Strong long-run persistence. 
Likelihood Ratio: Reject independence. Strong long-run persistence. 
Bowker: Reject symmetry. Directional/asymmetric movement (df=2).

8 Higher-Order Dependency Structure

build_markov_pairs <- function(data,
                               snapshot_years,
                               step,
                               full_levels) {
  df_snap <- data %>%
    dplyr::filter(year %in% snapshot_years) %>%
    dplyr::mutate(class = factor(class, levels = full_levels))
  df_pairs <- df_snap %>%
    dplyr::arrange(stateabb, year) %>%
    dplyr::group_by(stateabb) %>%
    dplyr::mutate(class_lag1 = dplyr::lag(class, 1),
                  year_lag1  = dplyr::lag(year, 1)) %>%
    dplyr::ungroup() %>%
    dplyr::filter(!is.na(class_lag1),
                  (year - year_lag1) %in% step) %>%
    dplyr::rename(class_curr = class)
  df_pairs
}

build_markov_triples <- function(data,
                                 snapshot_years,
                                 step,
                                 full_levels) {
  df_snap <- data %>%
    dplyr::filter(year %in% snapshot_years) %>%
    dplyr::mutate(class = factor(class, levels = full_levels))
  df_with_counts <- df_snap %>%
    dplyr::arrange(stateabb, year) %>%
    dplyr::group_by(stateabb) %>%
    dplyr::mutate(n_obs = dplyr::n()) %>%
    dplyr::ungroup()
  df_triple <- df_with_counts %>%
    dplyr::arrange(stateabb, year) %>%
    dplyr::group_by(stateabb) %>%
    dplyr::mutate(class_lag1 = dplyr::lag(class, 1),
                  year_lag1  = dplyr::lag(year,  1),
                  class_lag2 = dplyr::lag(class, 2),
                  year_lag2  = dplyr::lag(year,  2)) %>%
    dplyr::ungroup() %>%
    dplyr::filter(!is.na(class_lag1),
                  !is.na(class_lag2),
                  (year - year_lag1) %in% step,
                  (year_lag1 - year_lag2) %in% step) %>%
    dplyr::rename(class_curr = class)
}

# ERA 1: 1860-1880 ---- Build Markov pairs & triples
pairs1 <- build_markov_pairs(
  data           = df1,
  snapshot_years = snap1,
  step           = step1,
  full_levels    = levels_era1
)

triples1 <- build_markov_triples(
  data           = df1,
  snapshot_years = snap1,
  step           = step1,
  full_levels    = levels_era1
)

# ERA 2: 1881-1916 ---- Build Markov pairs & triples
pairs2 <- build_markov_pairs(
  data           = df2,
  snapshot_years = snap2,
  step           = step2,
  full_levels    = levels_era2
)

triples2 <- build_markov_triples(
  data           = df2,
  snapshot_years = snap2,
  step           = step2,
  full_levels    = levels_era2
)

9 Evaluation of Markovian Memory

test_first_order_markov_pairs <- function(pairs_df, full_levels) {
  df_counts <- pairs_df %>%
    dplyr::count(class_lag1, class_curr, name = "n") %>%
    dplyr::mutate(class_lag1 = factor(class_lag1, levels = full_levels),
                  class_curr = factor(class_curr, levels = full_levels)) %>%
    tidyr::complete(class_lag1, class_curr, fill = list(n = 0))
  colnames(df_counts) <- c("M", "N", "n")
  fit_H0 <- MASS::loglm(n ~ M + N, data = df_counts)   # 0th order (indep)
  fit_H1 <- MASS::loglm(n ~ M * N, data = df_counts)   # 1st order
  G2      <- fit_H0$deviance - fit_H1$deviance
  df_diff <- fit_H0$df - fit_H1$df
  p_val   <- 1 - pchisq(G2, df_diff)
  list(counts  = df_counts,
       G2      = G2,
       df      = df_diff,
       p_value = p_val,
       fit_H0  = fit_H0,
       fit_H1  = fit_H1)
}

test_second_order_markov_triples <- function(triples_df, full_levels) {
  df_counts <- triples_df %>%
    dplyr::count(class_lag2, class_lag1, class_curr, name = "n") %>%
    dplyr::mutate(class_lag2 = factor(class_lag2, levels = full_levels),
                  class_lag1 = factor(class_lag1, levels = full_levels),
                  class_curr = factor(class_curr, levels = full_levels)) %>%
    tidyr::complete(class_lag2, class_lag1, class_curr, fill = list(n = 0))
  colnames(df_counts) <- c("L", "M", "N", "n")
  fit_H0 <- MASS::loglm(n ~ L + M + N + L:M + M:N, data = df_counts)  # 1st-order Markov
  fit_H1 <- MASS::loglm(n ~ L * M * N, data = df_counts)              # 2nd-order (saturated)
  G2      <- fit_H0$deviance - fit_H1$deviance
  df_diff <- fit_H0$df - fit_H1$df
  p_val   <- 1 - pchisq(G2, df_diff)
  list(counts  = df_counts,
       G2      = G2,
       df      = df_diff,
       p_value = p_val,
       fit_H0  = fit_H0,
       fit_H1  = fit_H1)
}

interpret_markov_order <- function(p, comparison, alpha = 0.05) {

  if (comparison == "0vs1") {
    if (p < alpha) {
      return(paste0("Reject independence / zero-order (p < ", alpha,
                    "). Evidence of first-order state dependence."))
    } else {
      return(paste0("Fail to reject independence / zero-order (p ≥ ", alpha,
                    "). Transitions consistent with independence."))
    }
  }

  if (comparison == "1vs2") {
    if (p < alpha) {
      return(paste0("Reject first-order sufficiency (p < ", alpha,
                    "). Evidence of second-order dependence."))
    } else {
      return(paste0("Fail to reject first-order sufficiency (p ≥ ", alpha,
                    "). First-order Markov assumption appears sufficient."))
    }
  }

  stop("Unknown comparison type: use '0vs1' or '1vs2'")
}

9.0.1 ERA 1 1860-1880

# ERA 1: 1860-1880 ---- Zero vs First-order
test1_markov01 <- test_first_order_markov_pairs(pairs1, levels_era1)

cat("\nERA 1: Zero vs First-order\n")
cat("G^2 =", round(test1_markov01$G2, 3),
    " df =", test1_markov01$df,
    " p =", signif(test1_markov01$p_value, 4), "\n")
cat("Interpretation:",
    interpret_markov_order(test1_markov01$p_value, "0vs1"), "\n")

# ERA 1: 1860-1880 ---- First vs Second-order
test1_markov12 <- test_second_order_markov_triples(triples1, levels_era1)

cat("\nERA 1: First vs Second-order\n")
cat("G^2 =", round(test1_markov12$G2, 3),
    " df =", test1_markov12$df,
    " p =", signif(test1_markov12$p_value, 4), "\n")
cat("Interpretation:",
    interpret_markov_order(test1_markov12$p_value, "1vs2"), "\n")


ERA 1: Zero vs First-order
G^2 = 264  df = 4  p = 0 
Interpretation: Reject independence / zero-order (p < 0.05). Evidence of first-order state dependence. 

ERA 1: First vs Second-order
G^2 = 0.1  df = 12  p = 1 
Interpretation: Fail to reject first-order sufficiency (p ≥ 0.05). First-order Markov assumption appears sufficient.

9.0.2 ERA 2 1881-1916

# ERA 2: 1881-1916 ---- Zero vs First-order
test2_markov01 <- test_first_order_markov_pairs(pairs2, levels_era2)

cat("\nERA 2: Zero vs First-order\n")
cat("G^2 =", round(test2_markov01$G2, 3),
    " df =", test2_markov01$df,
    " p =", signif(test2_markov01$p_value, 4), "\n")
cat("Interpretation:",
    interpret_markov_order(test2_markov01$p_value, "0vs1"), "\n")

# ERA 2: 1881-1916 ---- First vs Second-order
test2_markov12 <- test_second_order_markov_triples(triples2, levels_era2)

cat("\nERA 2: First vs Second-order\n")
cat("G^2 =", round(test2_markov12$G2, 3),
    " df =", test2_markov12$df,
    " p =", signif(test2_markov12$p_value, 4), "\n")
cat("Interpretation:",
    interpret_markov_order(test2_markov12$p_value, "1vs2"), "\n")


ERA 2: Zero vs First-order
G^2 = 553  df = 9  p = 0 
Interpretation: Reject independence / zero-order (p < 0.05). Evidence of first-order state dependence. 

ERA 2: First vs Second-order
G^2 = 21  df = 36  p = 0.98 
Interpretation: Fail to reject first-order sufficiency (p ≥ 0.05). First-order Markov assumption appears sufficient.

10 Ergodicity Diagnostics

check_ergodicity <- function(P, max_power = 10, tol = 1e-12) {
  A <- (P > tol) * 1
  A_power <- A
  reachable_all <- FALSE
  for (k in 1:max_power) {
    if (all(A_power > 0)) {
      reachable_all <- TRUE
      break
    }
    A_power <- (A_power %*% A > 0) * 1
  }
  P_power <- P
  positive_diag_all <- FALSE
  for (k in 1:max_power) {
    if (all(diag(P_power) > tol)) {
      positive_diag_all <- TRUE
      break
    }
    P_power <- P_power %*% P
  }
  list(
    irreducible_approx = reachable_all,
    aperiodic_approx   = positive_diag_all,
    max_power_checked  = max_power
  )
}

interpret_ergodicity <- function(erg) {
  if (isTRUE(erg$irreducible_approx) && isTRUE(erg$aperiodic_approx)) {
    return("Chain is ergodic (approximately irreducible and aperiodic).")
  }
  if (isFALSE(erg$irreducible_approx) && isTRUE(erg$aperiodic_approx)) {
    return("Chain is not ergodic (fails approximate irreducibility).")
  }
  if (isTRUE(erg$irreducible_approx) && isFALSE(erg$aperiodic_approx)) {
    return("Chain is not ergodic (fails approximate aperiodicity).")
  }
  return("Ergodicity status unclear under approximation criteria.")
}

10.0.1 ERA 1 1860-1880

# -------- TESTS ----------

# ERA 1: 1860-1880 ---- Ergodicity
erg1 <- check_ergodicity(P1)
cat("\nERA 1 Ergodicity:\n")
print(erg1)
cat("Interpretation:",
    interpret_ergodicity(erg1), "\n")


ERA 1 Ergodicity:
$irreducible_approx
[1] FALSE

$aperiodic_approx
[1] TRUE

$max_power_checked
[1] 10

Interpretation: Chain is not ergodic (fails approximate irreducibility).

10.0.2 ERA 2 1881-1916

# ERA 2: 1881-1916 ---- Ergodicity
erg2 <- check_ergodicity(P2)
cat("\nERA 2 Ergodicity:\n")
print(erg2)
cat("Interpretation:",
    interpret_ergodicity(erg2), "\n")


ERA 2 Ergodicity:
$irreducible_approx
[1] FALSE

$aperiodic_approx
[1] TRUE

$max_power_checked
[1] 10

Interpretation: Chain is not ergodic (fails approximate irreducibility).

11 Comparison of Structural Change Across Eras

# Equality of Transition Matrices Between Eras
compare_eras_markov <- function(trans1_df, trans2_df,
                                era1_label = "Era1",
                                era2_label = "Era2",
                                common_levels = c("upper_class","middle_class","lower_class_1")) {

  prep <- function(df, era) {
    df %>%
      dplyr::mutate(
        class_prev = factor(class_prev, levels = common_levels),
        class      = factor(class,      levels = common_levels)
      ) %>%
      dplyr::filter(!is.na(class_prev), !is.na(class)) %>%
      dplyr::count(class_prev, class, name = "n") %>%
      dplyr::mutate(era = era)
  }
  df_counts <- dplyr::bind_rows(
    prep(trans1_df, era1_label),
    prep(trans2_df, era2_label)
  )
  colnames(df_counts) <- c("origin", "dest", "n", "era")
  fit_common <- MASS::loglm(
    n ~ era + origin + dest + origin:dest,
    data = df_counts
  )
  fit_diff <- MASS::loglm(
    n ~ era * origin * dest,
    data = df_counts
  )
  G2      <- fit_common$deviance - fit_diff$deviance
  df_diff <- fit_common$df - fit_diff$df
  p_val   <- 1 - pchisq(G2, df_diff)
  list(
    counts = df_counts,
    G2 = G2,
    df = df_diff,
    p_value = p_val,
    fit_common = fit_common,
    fit_diff = fit_diff
  )
}

res <- compare_eras_markov(
  trans1_df = trans1,          # Era 1 transitions
  trans2_df = trans2,          # Era 2 transitions
  era1_label = "1860-1880",
  era2_label = "1881-1916",
  common_levels = c("upper_class","middle_class","lower_class_1")
)

interpret_matrix_equality <- function(p, alpha = 0.05) {
  if (p < alpha) {
    paste0("Reject equality of transition matrices (p < ", alpha,
           "). Evidence of structural change between eras.")
  } else {
    paste0("Fail to reject equality of transition matrices (p ≥ ", alpha,
           "). Transition dynamics are statistically indistinguishable across eras.")
  }
}

# -------- TEST: Equality of transition matrices across eras ----------

res <- compare_eras_markov(
  trans1_df = trans1,
  trans2_df = trans2,
  era1_label = "1860-1880",
  era2_label = "1881-1916",
  common_levels = c("upper_class","middle_class","lower_class_1")
)

cat("\nEquality of transition matrices across eras:\n")
cat("G^2 =", round(res$G2, 3),
    " df =", res$df,
    " p =", signif(res$p_value, 4), "\n")
cat("Interpretation:",
    interpret_matrix_equality(res$p_value), "\n")


Equality of transition matrices across eras:
G^2 = 1.9  df = 8  p = 0.98 
Interpretation: Fail to reject equality of transition matrices (p ≥ 0.05). Transition dynamics are statistically indistinguishable across eras.

12 Appendix: Methodological Notes

# Balanced-panel of snapshot transitions

# Purpose:
#   Given a set of snapshot years {t_1, ..., t_T} and step length k,
#   construct transitions (X_t -> X_{t+k}) for a *balanced panel* of
#   countries (stateabb).
#
#   In Markov-chain terms, we observe a finite-state stochastic process
#   X_t(u) taking values in a fixed state space S = {upper_class, semi, lower_class_1*}
#   at discrete times t ∈ {t_1, ..., t_T}. A balanced panel holds S (the set
#   of *units* u) fixed across snapshots so that transitions are between
#   *states* rather than driven by changing sample composition
#   (see Dezzani 2002 on state-level Markov models for the world-system).
#
# Inputs:
#   data          : data frame with columns stateabb, year, class.
#   snapshot_years: vector of years (e.g., 1860, 1865, 1704, 1875, 1880).
#   step          : step size (e.g., 5 years).
#   full_levels   : complete factor levels for class in this era,
#                   enforcing a fixed state space S across models.
#
# Steps (data-analytic):
#   1. Filter to snapshot years t_s ∈ {t_1, ..., t_T}.
#   2. Within these years, keep only countries that appear in *all*
#      snapshot years → ensures a balanced panel:
#        U* = {u : u observed at every t_s}.
#      This isolates mobility *within* a fixed set of units, which is
#      critical for interpretable Markov mobility (Dezzani 2002; Babones
#      & Dezzani 2008).
#   3. For each u ∈ U*, sort by year and compute lagged pairs
#         X_t(u) and X_{t+k}(u).
#      These give observed transitions between states i → j.
#   4. Keep only pairs where (year - year_prev) == step, so that we
#      construct *k-step* chains (Anderson & Goodman 1957).
#   5. Define step_id = "t → t+k" to index the snapshot step s in
#      the 3-way table N_{sij}.
#
# Math:
#   For each country u and snapshot s (t_s, t_{s+1} = t_s + k):
#     X_t(u)      = class at time t_s
#     X_{t+k}(u)  = class at time t_{s+1}
#
#   These transitions produce 3D counts N_{sij} used later in:
#     - homogeneity tests across steps (N_{sij} as a 3-way table),
#     - pooled transition matrix estimation N_{ij} = Σ_s N_{sij}.
#
# References:
#   Dezzani (2002) "Measuring Transition and Mobility in the
#   Hierarchical World-Economy", Journal of Regional Science.
#   Babones & Dezzani (2008) "Mobility in the Modern World-Economy".
#   Anderson & Goodman (1957) "Statistical Inference about Markov Chains".


# Balanced-panel attrition diagnostic


# Purpose:
#   Compare class distributions for:
#     - all countries in a year vs
#     - countries retained in the balanced panel for the era.
#
#   This checks whether panel restrictions introduce compositional
#   bias in the Markov estimates (Dezzani 2002; Babones & Dezzani
#   2008).
#
# INTERPRETATION:
#   A chi-square comparison of class distributions for all
#   countries versus balanced-panel countries in the first
#   snapshot year provides a diagnostic for potential attrition
#   bias (Dezzani 2002; Babones & Dezzani 2008).


# 5-year pooled transition matrices


# Purpose:
#   Pool transitions across steps to estimate:
#     - transition probabilities \hat{P} = (P_ij),
#     - stationary distribution \hat{π}.
#
# Data:
#   From all transitions (across steps), we construct a 2D table:
#     N_{ij} = Σ_s N_{sij},
#   where i = origin, j = destination.
#
# Estimation:
#   Define:
#     N_{ij}  = total transitions from i to j,
#     N_{i+}  = Σ_j N_{ij} (row totals).
#   Then MLE for one-step Markov transition probabilities is:
#     \hat{P}_{ij} = N_{ij} / N_{i+}.                (Anderson & Goodman 1957)
#
# Stationary distribution:
#   The stationary distribution π satisfies:
#     π \hat{P} = π,     Σ_i π_i = 1.
#
#   Numerically, we obtain π as the dominant left eigenvector of \hat{P}:
#     Let v be a right eigenvector of \hat{P}^T:
#       \hat{P}^T v = λ v.
#     If λ = 1 is the largest eigenvalue (for an irreducible, aperiodic chain),
#     then v (normalized to sum to 1) is the stationary distribution:
#       π̂ = v / Σ_i v_i.
#
# References:
#   Anderson & Goodman (1957) on ML estimation and asymptotic theory
#   for discrete-state Markov chains.
#   Bickenbach & Bode (2001) for applications of such matrices to
#   regional/income convergence.


# 5-year step homogeneity tests


# Purpose:
#   Test whether transition probabilities P_ij are identical across
#   snapshot steps s (time-homogeneous Markov chain within an era).
#
#   Data structure:
#     - We construct a 3D contingency table N_{sij}:
#         S = step (snapshot interval),
#         O = origin class,
#         D = destination class.
#     - In log-linear form, this is a standard three-way table S × O × D
#       (Bishop, Fienberg & Holland 1975/2007; Goodman 1971).
#
# Model (log-linear representation):
#   Let μ_{sij} = E[N_{sij}] denote expected counts under a given model.
#
#   H0: Homogeneous Markov model (OD association invariant in s)
#     log μ_{sij} = λ + λ_s^(S) + λ_i^(O) + λ_j^(D) + λ_{ij}^{(OD)}
#     → includes:
#         - main effects for step, origin, destination,
#         - origin-destination association λ_{ij}^{(OD)} *not* varying
#           with s, corresponding to a time-homogeneous chain.
#
#   H1: Fully heterogeneous model
#     log μ_{sij} = λ + all interactions among (S, O, D)
#                 = λ + λ_s + λ_i + λ_j
#                   + λ_{si} + λ_{sj} + λ_{ij} + λ_{sij}.
#     → OD association allowed to differ by step via λ_{sij}.
#
# Test statistic:
#   G^2 = deviance(H0) - deviance(H1),
#   df  = df(H0) - df(H1),
#   p-value from χ^2_df.
#
#   Here, deviance is the standard log-likelihood ratio statistic:
#     G^2(model) = 2 Σ_{sij} n_{sij} log( n_{sij} / μ̂_{sij} ),
#   which is 2×(sample size) times a Kullback-Leibler divergence between
#   observed and fitted cell probabilities (Kullback & Leibler 1951;
#   Kullback, Kupperman & Ku 1962). Under regularity conditions, G^2
#   is asymptotically χ^2 (Anderson & Goodman 1957; Bishop et al. 1975).
#
# Relation to Markov-chain diagnostics:
#   - This is essentially the log-linear version of testing whether
#     transition matrices differ across time slices, analogous to the
#     homogeneity tests in Bickenbach & Bode’s (2001) discussion of
#     Markov chains for income distributions.


# Long-interval transition diagnostics (D&B-style)


# Purpose:
#   Assess long-run persistence and directionality by testing the
#   association between countries’ origin class at the first snapshot
#   year t0 and destination class at the last snapshot year tT.
#
# Data structure:
#   - Construct a 2D contingency table N_ij:
#       O = origin class at t0,
#       D = destination class at tT,
#     where N_ij counts countries observed at BOTH endpoints.
#   - Note: this is an endpoint-matched sample (not a fully balanced
#     panel across all intermediate snapshots).
#
# Models / null hypotheses:
#   1) Independence (association / persistence):
#        H0: O ⟂ D
#      Tested with:
#        - Pearson X^2 on N_ij, and
#        - Likelihood-ratio G^2 from the log-linear independence model:
#            log μ_ij = λ + λ_i^(O) + λ_j^(D)
#
#   2) Symmetry (directionality):
#        H0: N_ij = N_ji for all i ≠ j
#      Tested with Bowker’s symmetry chi-square (generalized McNemar).
#
# Interpretation:
#   - Reject independence (Pearson or G^2): origin predicts destination
#     over the long interval (long-run persistence / stratification).
#   - Reject symmetry (Bowker): net directional movement between at least
#     one class pair (asymmetric i→j vs j→i flows).
#
# Relation to other diagnostics:
#   - Complements step-homogeneity tests: it does NOT test time-homogeneity
#     across intermediate 5-year steps; it summarizes start–end dependence
#     over the full era (Dezzani 2002; Babones & Dezzani 2008).


# Markov pairs X_t, X_{t+1} (0 → 1 order)


# Purpose:
#   Construct ordered pairs (X_t, X_{t+1}) for estimating and testing
#   *zero-order vs first-order* dependence in a discrete-state process.
#
#   Specifically, the pair table supports the hypothesis test:
#
#     H0 (0th order / independence):   X_{t+1} ⟂ X_t
#     H1 (1st order / Markov):         X_{t+1} depends on X_t
#
#   In a Markov-chain setting, H1 corresponds to the standard
#   first-order dependence assumption used to estimate a transition
#   matrix P with entries P_ij = Pr(X_{t+1}=j | X_t=i), as developed in
#   Anderson & Goodman (1957, Annals of Mathematical Statistics).
#
#   In contingency-table terms, we represent transitions as a 2-way table
#   N_{ij} of origin (i = X_t) by destination (j = X_{t+1}) counts and fit
#   hierarchical log-linear models:
#
#     - Independence (0th order):
#         log μ_{ij} = λ + λ_i^(O) + λ_j^(D)
#         (no origin-destination association)
#
#     - Association (1st order):
#         log μ_{ij} = λ + λ_i^(O) + λ_j^(D) + λ_{ij}^{(OD)}
#
#   The likelihood-ratio deviance G^2 compares these models:
#
#     G^2 = 2 Σ_{ij} n_{ij} log(n_{ij} / μ̂_{ij})
#
#   which is 2N times a Kullback-Leibler divergence between empirical
#   and fitted cell probabilities (Kullback & Leibler 1951; Kullback,
#   Kupperman & Ku 1962). Under standard regularity conditions, G^2
#   is asymptotically χ² (Anderson & Goodman 1957; Bishop, Fienberg &
#   Holland 1975/2007; Goodman 1971).
#
# Inputs:
#   data          : data frame with columns stateabb, year, class.
#   snapshot_years: ordered vector of observation years.
#   step          : spacing constraint between consecutive snapshots
#                   (e.g., 7; 8; or a vector like c(7,7,6)).
#   full_levels   : factor levels for class (fixed state space S).
#
# Output:
#   data frame with columns:
#     stateabb, year, class_prev, class_curr
#   restricted to balanced-panel countries and exact step spacing.


# Markov triples X_{t-1}, X_t, X_{t+1} (1 → 2 order)


# Purpose:
#   Construct triples (X_{t-1}, X_t, X_{t+1}) for testing the
#   first-order Markov property:
#
#     H0: X_{t+1} ⟂ X_{t-1} | X_t
#
#   Following Anderson & Goodman (1957, Annals of Mathematical
#   Statistics) and the log-linear representation of Markov chains
#   in Bishop, Fienberg & Holland (1975/2007, Discrete Multivariate
#   Analysis) and Goodman (1971, Technometrics).
#
# Inputs:
#   data          : df with stateabb, year, class
#   snapshot_years: ordered vector of observation years
#   step          : spacing between snapshots (e.g., 7, 8)
#   full_levels   : factor levels for class
#
# Output:
#   data frame with columns:
#     stateabb, year, class_lag2, class_lag1, class_curr
#   restricted to balanced-panel countries and exact step spacing.
#
# References:
#   Anderson, T.W. & Goodman, L.A. (1957). "Statistical Inference about
#     Markov Chains." Annals of Mathematical Statistics.
#   Bishop, Y.M.M., Fienberg, S.E., & Holland, P.W. (1975/2007).
#     Discrete Multivariate Analysis. MIT Press.
#   Goodman, L.A. (1971). "The Analysis of Multidimensional Contingency
#     Tables." Technometrics.
#   Kullback, S. & Leibler, R.A. (1951). "On Information and Sufficiency."
#     Annals of Mathematical Statistics.
#   Kullback, S., Kupperman, M., & Ku, H. (1962). "Tests for Contingency
#     Tables and Markov Chains." Technometrics.


# 0 vs 1 order tests (X_{t-1}, X_t)


# Purpose:
#   Test the first-order Markov property:
#
#     H0: X_{t+1} ⟂ X_{t-1} | X_t
#
#   using a three-way contingency table:
#     L = X_{t-1}, M = X_t, N = X_{t+1}.
#
# Log-linear models:
#
#   Let μ_{ijk} = E[N_{ijk}] = expected counts for L=i, M=j, N=k.
#
#   H0 (first-order Markov):
#     log μ_{ijk} = λ
#                 + λ_i^L + λ_j^M + λ_k^N
#                 + λ_{ij}^{LM} + λ_{jk}^{MN}
#
#     → includes L-M and M-N interactions, but no L-N or L-M-N
#       interactions, corresponding to X_{t+1} ⟂ X_{t-1} | X_t.
#
#   H1 (saturated three-way):
#     log μ_{ijk} = λ + all main effects + all interactions
#                 = L * M * N in loglm() formula notation.
#
# LR test:
#   G^2 = dev(H0) - dev(H1),
#   df  = df(H0) - df(H1),
#   p   = 1 - F_χ²(G^2; df).
#
# References:
#   Anderson & Goodman (1957) "Statistical Inference about Markov
#   Chains", Annals of Mathematical Statistics.
#   Bishop, Fienberg & Holland (1975/2007) Discrete Multivariate
#   Analysis, MIT Press.
#   Goodman (1971) "The Analysis of Multidimensional Contingency
#   Tables", Technometrics.


# 1 vs 2 order tests (X_{t-2}, X_{t-1}, X_t)


# Purpose:
#   Test whether X_{t} depends on X_{t-2k} given X_{t-k}:
#
#       H0 (1st-order): N ⟂ L | M
#       H1 (2nd-order): N depends on both M and L
#
# Models (log-linear):
#
#   H0: log μ_LMN = L + M + N + L:M + M:N
#   H1: log μ_LMN = L * M * N
#
# LR statistic:
#   G^2 = dev(H0) - dev(H1)
#   df  = df(H0) - df(H1)
#
# Citations:
#   Anderson & Goodman (1957)
#   Bishop, Fienberg & Holland (1975/2007)
#   Goodman (1971)


# Ergodicity diagnostic (irreducibility & aperiodicity)


# Purpose:
#   Provide diagnostics for:
#     - irreducibility (all states communicate),
#     - aperiodicity (no deterministic cycles),
#   using powers of the transition matrix P.
#
# Approach:
#   - Construct adjacency matrix A with A_ij = 1 if P_ij > 0.
#   - Compute reachability via powers of A.
#   - For aperiodicity, check whether some power of P has strictly
#     positive diagonal entries for all states.
#
# Interpretation:
#   If some k ≤ max_power has all entries of A^k > 0, the chain
#   is regular (thus irreducible and aperiodic).
#
# References:
#   Anderson & Goodman (1957) on long-run properties of Markov chains.
#   Standard Markov chain theory for ergodicity criteria.


# Equality of transition matrices between eras


# Purpose:
#   Compare transition structures between Era 1 and Era 2 on the
#   shared class set (upper_class, middle_class, lower_class_1), using a
#   4-way contingency table with factors:
#
#     Era (E), Origin (O), Destination (D)
#
# Log-linear setup:
#   H0: common Markov structure across eras:
#       log μ_{eij} = λ + λ_e^E + λ_i^O + λ_j^D + λ_{ij}^{OD}
#
#   H1: era-specific OD association:
#       log μ_{eij} = λ + all main effects + all interactions
#                   = E * O * D
#
# LR test:
#   G^2 = dev(H0) - dev(H1),
#   df  = df(H0) - df(H1).
#
# References:
#   Bishop, Fienberg & Holland (1975/2007), Goodman (1971),
#   Bickenbach & Bode (2001) on comparing Markov structures
#   across groups/periods.

13 Bibliography

Anderson, T. W., & Goodman, L. A. (1957). Statistical inference about Markov chains. Annals of Mathematical Statistics.

Babones, S. (2005). The country-level income structure of the world-economy. Journal of World-Systems Research.

Babones, S., & Dezzani, R. (2008). Mobility in the modern world-economy. Conference paper.

Dezzani, R. (2002). Measuring transition and mobility in the hierarchical world-economy. Journal of Regional Science.

Bickenbach, F., & Bode, E. (2001). Markov or not Markov - This should be a question. Empirical Economics.

Bishop, Y. M. M., Fienberg, S. E., & Holland, P. W. (1975/2007). Discrete Multivariate Analysis. MIT Press.

Goodman, L. A. (1971). The analysis of multidimensional contingency tables. Technometrics.

Kullback, S., & Leibler, R. A. (1951). On information and sufficiency. Annals of Mathematical Statistics.

Kullback, S., Kupperman, M., & Ku, H. (1962). Tests for contingency tables and Markov chains. Technometrics.

A Markov Analysis of Country-Level Social Mobility, 1860-1916

Jared Walker

ChatGPT