Created
June 19, 2024 10:41
-
-
Save vankesteren/deeb4576a4d9b593dff6772be50ac7df to your computer and use it in GitHub Desktop.
Revisions
-
vankesteren created this gist
Jun 19, 2024 .There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,134 @@ # 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.