ols-bad-when-sparse.R

# Load necessary libraries
# install.packages(c("igraph", "AER", "dplyr", "tidyr"))
library(igraph)
library(AER)
library(dplyr)
library(tidyr)
library(Matrix)

# Set true parameters
TRUE_BETA <- 0.4 # Endogenous peer effect
TRUE_GAMMA <- 1.0 # Exogenous direct effect
TRUE_DELTA <- 0.5 # Exogenous peer effect
DEGREE <- 4 # Bounded degree for the k-regular graph

# Function to simulate data and estimate models for a given network size n
run_sim_iteration <- function(n, k, beta, gamma, delta) {
  # 1. Generate bounded degree network (k-regular graph)
  # This guarantees max degree is bounded regardless of n
  g <- sample_k_regular(n, k)
  A <- as_adj(g, sparse = TRUE)

  # Row-normalize the adjacency matrix to create G
  rs <- rowSums(A)
  G <- A / rs

  # 2. Generate exogenous data and errors
  X <- matrix(rnorm(n), ncol = 1)
  epsilon <- matrix(rnorm(n), ncol = 1)

  # Calculate spatial lags of X
  GX <- as.matrix(G %*% X)

  # 3. Generate outcome Y via the reduced form
  # Y = (I - beta * G)^(-1) (X*gamma + GX*delta + epsilon)
  I <- diag(n)
  Y <- solve(I - beta * G, as.numeric(X * gamma + GX * delta + epsilon))

  # Calculate spatial lag of Y
  GY <- as.matrix(G %*% Y)

  # 4. Estimation
  # OLS: Biased due to E[GY * epsilon] != 0
  # Note: Suppressing the intercept for strict alignment with the theoretical DGP
  ols_fit <- lm(Y ~ GY + X + GX - 1)
  beta_ols <- coef(ols_fit)["GY"]

  # TSLS: Use G^2X as the instrument for GY
  G2X <- as.matrix(G %*% GX)
  tsls_fit <- ivreg(Y ~ GY + X + GX - 1 | X + GX + G2X - 1)
  beta_tsls <- coef(tsls_fit)["GY"]

  return(c(OLS = unname(beta_ols), TSLS = unname(beta_tsls)))
}

# Define the simulation parameters
n_values <- c(100, 500, 1000, 2000)
simulations_per_n <- 50

# Run the simulation
set.seed(42)
results <- list()

cat("Running simulations...\n")
for (n in n_values) {
  cat("Simulating for n =", n, "...\n")

  # Replicate the data generation and estimation process
  sim_reps <- replicate(
    simulations_per_n,
    run_sim_iteration(n, DEGREE, TRUE_BETA, TRUE_GAMMA, TRUE_DELTA)
  )

  # Store the mean estimates across the iterations
  results[[as.character(n)]] <- data.frame(
    n = n,
    mean_beta_ols = mean(sim_reps["OLS", ]),
    mean_beta_tsls = mean(sim_reps["TSLS", ]),
    var_beta_ols = var(sim_reps["OLS", ]),
    var_beta_tsls = var(sim_reps["TSLS", ])
  )
}

# Combine results into a final table
final_results <- bind_rows(results)

# Display the results
print("Simulation Results (True Beta = 0.4):")
print(final_results, row.names = FALSE)