Example of a probabilistic social simulation

probsocsim.R

# Simple probabilistic simulation script

# Causal graph: 
# NetIncome -> + CulturalActivities
# NetIncome -> + SportsActivities
# NetIncome -> - Debts
# NetIncome -> + SocialComparison
# SportsActivities -> + Health
# SportsActivities -> + Partnership
# CulturalActivities -> + SocialComparison
# CulturalActivities -> + Partnership
# Partnership -> + Welfare
# SocialComparison -> - Welfare
# Health -> + Welfare
# Debts -> - Welfare

library(pbapply)
library(parallel)


N <- 200000

NetIncome <- exp(rnorm(N, 7, 2))
CulturalActivities <- rpois(N, 4 + .07 * log(NetIncome))
SportsActivities <- rpois(N, 10 + .14 * log(NetIncome))
Debt <- exp(rnorm(N, 9 - .4 * log(NetIncome), min(0.1, 2 - 0.1 * log(NetIncome))))
SocialComparison <- pmax(0, pmin(10, rnorm(N, 4 + .5 * log(NetIncome) + 0.1*log(CulturalActivities), 1)))
Health <- pmax(0, pmin(10, rnorm(N, 6.5 + 0.2*log(SportsActivities), 2.5)))
Partnership <- rbinom(N, 1, 1/(1 + exp(-0.1*log(SportsActivities)-0.1*log(CulturalActivities))))
Welfare_raw <- rnorm(N, 0.12*Partnership - 0.03*SocialComparison + 0.4*Health - (Debt > 1000)*.33*log(Debt), 1.5)
Welfare <- 1 / (1 + exp(-Welfare_raw)) * 10


df_old <- data.frame(
  net_income = NetIncome,
  cultural_activity = CulturalActivities,
  sports_activity = SportsActivities, 
  debt = Debt,
  social_comparison = SocialComparison,
  health = Health,
  partnership = Partnership,
  welfare = Welfare
)

# now, we give everyone a basic minimum income and run the sim again
NetIncome[NetIncome<150] <- 150
CulturalActivities <- rpois(N, 4 + .07 * log(NetIncome))
SportsActivities <- rpois(N, 10 + .14 * log(NetIncome))
Debt <- exp(rnorm(N, 9 - .4 * log(NetIncome), min(0.1, 2 - 0.1 * log(NetIncome))))
SocialComparison <- pmax(0, pmin(10, rnorm(N, 4 + .5 * log(NetIncome) + 0.1*log(CulturalActivities), 1)))
Health <- pmax(0, pmin(10, rnorm(N, 6.5 + 0.2*log(SportsActivities), 2.5)))
Partnership <- rbinom(N, 1, 1/(1 + exp(-0.1*log(SportsActivities)-0.1*log(CulturalActivities))))
Welfare_raw <- rnorm(N, 0.12*Partnership - 0.03*SocialComparison + 0.4*Health - (Debt > 1000)*.33*log(Debt), 1.5)
Welfare <- 1 / (1 + exp(-Welfare_raw)) * 10

df_new <- data.frame(
  net_income = NetIncome,
  cultural_activity = CulturalActivities,
  sports_activity = SportsActivities, 
  debt = Debt,
  social_comparison = SocialComparison,
  health = Health,
  partnership = Partnership,
  welfare = Welfare
)

# let's look at the new distribution of welfare
plot(density(df_new$welfare), lty = 2)
lines(density(df_old$welfare))

# welfare improvement (ATE)
hist(df_new$welfare - df_old$welfare, breaks = "FD")
mean(df_new$welfare - df_old$welfare)

# for the ones at the bottom (ATT)
hist(df_new$welfare[df_old$net_income<150] - df_old$welfare[df_old$net_income<150], breaks = "FD")
mean(df_new$welfare[df_old$net_income<150] - df_old$welfare[df_old$net_income<150])

# for inference: just do this 1000 times
clus <- makeCluster(10)
effect <- pbsapply(1:1000, function(i) {
  N <- 200000
  NetIncome <- exp(rnorm(N, 7, 2))
  CulturalActivities <- rpois(N, 4 + .07 * log(NetIncome))
  SportsActivities <- rpois(N, 10 + .14 * log(NetIncome))
  Debt <- exp(rnorm(N, 9 - .4 * log(NetIncome), min(0.1, 2 - 0.1 * log(NetIncome))))
  SocialComparison <- pmax(0, pmin(10, rnorm(N, 4 + .5 * log(NetIncome) + 0.1*log(CulturalActivities), 1)))
  Health <- pmax(0, pmin(10, rnorm(N, 6.5 + 0.2*log(SportsActivities), 2.5)))
  Partnership <- rbinom(N, 1, 1/(1 + exp(-0.1*log(SportsActivities)-0.1*log(CulturalActivities))))
  Welfare_raw <- rnorm(N, 0.12*Partnership - 0.03*SocialComparison + 0.4*Health - (Debt > 1000)*.33*log(Debt), 1.5)
  Welfare <- 1 / (1 + exp(-Welfare_raw)) * 10
  
  
  df_old <- data.frame(
    net_income = NetIncome,
    cultural_activity = CulturalActivities,
    sports_activity = SportsActivities, 
    debt = Debt,
    social_comparison = SocialComparison,
    health = Health,
    partnership = Partnership,
    welfare = Welfare
  )
  
  # now, we give everyone a basic minimum income
  NetIncome[NetIncome<150] <- 150
  CulturalActivities <- rpois(N, 4 + .07 * log(NetIncome))
  SportsActivities <- rpois(N, 10 + .14 * log(NetIncome))
  Debt <- exp(rnorm(N, 9 - .4 * log(NetIncome), min(0.1, 2 - 0.1 * log(NetIncome))))
  SocialComparison <- pmax(0, pmin(10, rnorm(N, 4 + .5 * log(NetIncome) + 0.1*log(CulturalActivities), 1)))
  Health <- pmax(0, pmin(10, rnorm(N, 6.5 + 0.2*log(SportsActivities), 2.5)))
  Partnership <- rbinom(N, 1, 1/(1 + exp(-0.1*log(SportsActivities)-0.1*log(CulturalActivities))))
  Welfare_raw <- rnorm(N, 0.12*Partnership - 0.03*SocialComparison + 0.4*Health - (Debt > 1000)*.33*log(Debt), 1.5)
  Welfare <- 1 / (1 + exp(-Welfare_raw)) * 10
  
  df_new <- data.frame(
    net_income = NetIncome,
    cultural_activity = CulturalActivities,
    sports_activity = SportsActivities, 
    debt = Debt,
    social_comparison = SocialComparison,
    health = Health,
    partnership = Partnership,
    welfare = Welfare
  )
  
  c(ATE = mean(df_new$welfare - df_old$welfare), ATT = mean(df_new$welfare[df_old$net_income<150] - df_old$welfare[df_old$net_income<150]))
}, cl = clus)

# Histogram of ATE
hist(effect["ATE",], breaks = "FD")
hist(effect["ATT",], breaks = "FD")

# NB: this does not include parameter & model uncertainty, which it should.