# GEM firm size data
wd = gsub("Figures", "Empirical Data/GEM", dir)
setwd(wd)
energy = fread("energy_v_firm.csv")
# regression
x = log(energy$energy)
y = log(energy$firm_mean)
r = lm(y ~ x)
a = exp(coef(r)[1])
b = coef(r)[2]
e_predict = function(x){
y = (x/a)^(1/b)
return(y)
}
energy_predict = e_predict(mod$firm_size)
mod$energy = energy_predict
# Kohler with energy estimates
####################################################
# energy bounds
wd = gsub("Figures", "Empirical Data/Energy Adaptation", dir)
setwd(wd)
bounds = fread("energy_boundaries.csv")
# kohler
wd = gsub("Figures", "Empirical Data/Kohler", dir)
setwd(wd)
kohler = fread("kohler.csv")
# merge with energy bounds
kohler_energy = merge(bounds, kohler, by = "Adaptation")
# plot
kohler_plot = ggplot() +
geom_point(data = mod, aes(x = energy, y = pay_gini, color = pay_exponent), size = 0.01, alpha = 0.2) +
geom_errorbarh(data = kohler_energy, aes(xmin = `energy_5%`, xmax = `energy_95%`, y = Gini), size = 0.2, height = 0) +
geom_point(data = kohler_energy, aes(x = `energy_50%`, y = Gini), size = 0.8) +
scale_x_log10("Energy Use per Capita (GJ/year)",  breaks = c(2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000)) +
scale_y_continuous("Gini Index", breaks = seq(0,1, 0.1)) +
scale_color_gradientn(expression(beta), colours = rainbow(8), breaks = seq(0, 1, 0.1)) +
coord_cartesian(xlim = seq(2.5, 2000), ylim = c(0.08, 0.8)) +
ggtitle("A.  Ancient Societies") +
theme_bw() +
theme(panel.border = element_rect(color = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(face="bold", size = rel(1), hjust = 0.5),
legend.position = "none",
axis.line = element_line(color = "black"),
axis.title.x=element_text(vjust= 0, size=rel(0.9)),
axis.title.y=element_text(vjust= 1.1, size=rel(0.9)),
axis.text.x = element_text(margin=margin(5,5,0,0,"pt")),
axis.text.y = element_text(margin=margin(3,5,0,3,"pt")),
axis.ticks.length = unit(-0.7, "mm"),
text=element_text(size = text.size, family="Times"))
# milanovic
#####################################
wd = gsub("Figures", "Empirical Data/Pre Industrial", dir)
setwd(wd)
milanovic = fread("milanovic_energy.csv")
wd = gsub("Figures", "Empirical Data/Pre Industrial", dir)
setwd(wd)
milanovic = fread("pre_industrial_ineq_energy.csv")
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Empirical Data/Pre Industrial/gdp_energy.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
kcal_to_gj = 0.000004186798
top_frac = 0.01
frontier_gini_func = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
# energy use of all individuals in population
energy_vec = c(rep(energy_subsistence, n_people-1), energy_top)
gini = ineq(energy_vec)
return(gini)
}
kcal_to_gj = 0.000004186798
top_frac = 0.01
energy_subsistence_low = kcal_to_gj*2000*365
energy_low = exp(seq(log(energy_subsistence_low), log(3000), length.out = 1000))
top_low = mapply(top_share, energy_subsistence_low, energy_low, top_n)
energy_pc = exp(seq(log(energy_subsistence_low), log(3000), length.out = 1000))
top_low = mapply(frontier_gini_func, energy_subsistence_low, energy_pc, top_frac)
library(ggplot2)
library(gridExtra)
library(data.table)
library(hmod)
library(magrittr)
library(scales)
library(here)
library(ineq)
energy_subsistence_high = kcal_to_gj*3000*365
energy_subsistence_high = kcal_to_gj*3000*365
energy_pc = exp(seq(log(energy_subsistence_high), log(3000), length.out = 1000))
top_high = mapply(frontier_gini_func, energy_subsistence_low, energy_pc, top_frac)
library(ggplot2)
library(gridExtra)
library(data.table)
library(hmod)
library(magrittr)
library(scales)
library(here)
library(ineq)
text.size = 10
# inequality frontier GINI
####################################################################################
frontier_gini_func = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
# energy use of all individuals in population
energy_vec = c(rep(energy_subsistence, n_people-1), energy_top)
gini = ineq(energy_vec)
return(gini)
}
top_frac = 0.01
kcal_to_gj = 0.000004186798
# inequality frontier with low subsistence energy
energy_subsistence_low = kcal_to_gj*2000*365
energy_pc = exp(seq(log(energy_subsistence_low), log(3000), length.out = 1000))
frontier_low = mapply(frontier_gini_func, energy_subsistence_low, energy_pc, top_frac)
# inequality frontier with high subsistence energy
energy_subsistence_high = kcal_to_gj*3000*365
energy_pc = exp(seq(log(energy_subsistence_high), log(3000), length.out = 1000))
frontier_high = mapply(frontier_gini_func, energy_subsistence_low, energy_pc, top_frac)
frontier = rbind(
data.table(energy = energy_low, top = top_low),
data.table(energy = energy_high, top = top_high)[order(-top)],
data.table(energy = energy_low, top = top_low)[1,]
)
energy_subsistence_low = kcal_to_gj*2000*365
energy_pc = exp(seq(log(energy_subsistence_low), log(3000), length.out = 1000))
energy_subsistence_low = kcal_to_gj*2000*365
energy_pc = exp(seq(log(energy_subsistence_low), log(3000), length.out = 1000))
frontier_low = mapply(frontier_gini_func, energy_subsistence_low, energy_pc, top_frac)
energy_subsistence_high = kcal_to_gj*3000*365
energy_subsistence_high = kcal_to_gj*3000*365
energy_pc = exp(seq(log(energy_subsistence_high), log(3000), length.out = 1000))
frontier_high = mapply(frontier_gini_func, energy_subsistence_high, energy_pc, top_frac)
energy_subsistence = 3
energy_max = 1000
energy_pc = exp(seq(log(energy_subsistence), log(energy_max), length.out = 1000))
}
frontier_gini_func = function(energy_subsistence, energy_max, top_frac){
# generate energy vector
energy_pc = exp(seq(log(energy_subsistence), log(energy_max), length.out = 1000))
subfunction = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
# energy use of all individuals in population
energy_vec = c(rep(energy_subsistence, n_people-1), energy_top)
gini = ineq(energy_vec)
return(gini)
}
gini_vec = mapply(subfuction, energy_subsistence, energy_pc, top_frac)
return(gini_vec)
}
energy_subsistence_low = kcal_to_gj*2000*365
frontier_low = frontier_gini_func(energy_subsistence_low, 3000, 0.01)
subfunction = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
# energy use of all individuals in population
energy_vec = c(rep(energy_subsistence, n_people-1), energy_top)
gini = ineq(energy_vec)
return(gini)
}
gini_vec = mapply(subfunction, energy_subsistence, energy_pc, top_frac)
return(gini_vec)
frontier_gini_func = function(energy_subsistence, energy_max, top_frac){
# generate energy vector
energy_pc = exp(seq(log(energy_subsistence), log(energy_max), length.out = 1000))
subfunction = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
# energy use of all individuals in population
energy_vec = c(rep(energy_subsistence, n_people-1), energy_top)
gini = ineq(energy_vec)
return(gini)
}
gini_vec = mapply(subfunction, energy_subsistence, energy_pc, top_frac)
return(gini_vec)
}
frontier_gini_func = function(energy_subsistence, energy_max, top_frac){
# generate energy vector
energy_pc = exp(seq(log(energy_subsistence), log(energy_max), length.out = 1000))
subfunction = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
# energy use of all individuals in population
energy_vec = c(rep(energy_subsistence, n_people-1), energy_top)
gini = ineq(energy_vec)
return(gini)
}
gini_vec = mapply(subfunction, energy_subsistence, energy_pc, top_frac)
return(gini_vec)
}
top_frac = 0.01
kcal_to_gj = 0.000004186798
energy_subsistence_low = kcal_to_gj*2000*365
frontier_low = frontier_gini_func(energy_subsistence_low, 3000, 0.01)
plot(frontier_low)
frontier_gini_func = function(energy_subsistence, energy_max, top_frac){
# generate energy vector
energy_pc = exp(seq(log(energy_subsistence), log(energy_max), length.out = 1000))
subfunction = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
# energy use of all individuals in population
energy_vec = c(rep(energy_subsistence, n_people-1), energy_top)
gini = ineq(energy_vec)
return(gini)
}
gini_vec = mapply(subfunction, energy_subsistence, energy_pc, top_frac)
output = data.table(energy_pc, gini = gini_vec)
return(output)
}
top_frac = 0.01
kcal_to_gj = 0.000004186798
energy_subsistence_low = kcal_to_gj*2000*365
frontier_low = frontier_gini_func(energy_subsistence_low, 3000, 0.01)
plot(frontier_low)
plot(frontier_low, log = "x")
energy_subsistence_low = kcal_to_gj*3000*365
kcal_to_gj = 0.000004186798
top_frac = 0.01
energy_max
energy_max = 3000
energy_subsistence_low = kcal_to_gj*3000*365
kcal_to_gj = 0.000004186798
top_frac = 0.01
energy_max = 3000
energy_subsistence_low = kcal_to_gj*3000*365
frontier_low = frontier_gini_func(energy_subsistence_low, energy_max, 0.01)
energy_subsistence_low = kcal_to_gj*3000*365
energy_subsistence_high = kcal_to_gj*4000*365
frontier_low = frontier_gini_func(energy_subsistence_low, energy_max, 0.01)
frontier_low = frontier_gini_func(energy_subsistence_low, energy_max, 0.01)
frontier_high = frontier_gini_func(energy_subsistence_high, energy_max, 0.01)
frontier_low = frontier_gini_func(energy_subsistence_low, energy_max, 0.01)
frontier_high = frontier_gini_func(energy_subsistence_high, energy_max, 0.01)
frontier = rbind(
frontier_low,
frontier_high[order(-gini)],
frontier_low[1,]
)
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
gini_frontier
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
gini_frontier
frontier_top_func = function(energy_subsistence, energy_max, top_frac){
# generate energy vector
energy_pc = exp(seq(log(energy_subsistence), log(energy_max), length.out = 1000))
subfunction = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
top_share = energy_top/energy_bottom
return(top_share)
}
top_share_vec = mapply(subfunction, energy_subsistence, energy_pc, top_frac)
output = data.table(energy_pc, gini = top_share_vec)
return(output)
}
# calculations
kcal_to_gj = 0.000004186798
top_frac = 0.01
energy_max = 3000
energy_subsistence_low = kcal_to_gj*2500*365
energy_subsistence_high = kcal_to_gj*4000*365
frontier_low = frontier_top_func(energy_subsistence_low, energy_max, 0.01)
frontier_high = frontier_top_func(energy_subsistence_high, energy_max, 0.01)
frontier = rbind(
frontier_low,
frontier_high[order(-gini)],
frontier_low[1,]
)
frontier_top_func = function(energy_subsistence, energy_max, top_frac){
# generate energy vector
energy_pc = exp(seq(log(energy_subsistence), log(energy_max), length.out = 1000))
subfunction = function(energy_subsistence, energy_pc, top_frac){
# population size
n_people = 1/top_frac
# total energy consumption of bottom class
energy_bottom = energy_subsistence*(n_people - 1)
# total energy of population
energy_total = energy_pc*n_people
# energy of elite class
energy_top = energy_total - energy_bottom
top_share = energy_top/energy_bottom
return(top_share)
}
top_share_vec = mapply(subfunction, energy_subsistence, energy_pc, top_frac)
output = data.table(energy_pc, top_share = top_share_vec)
return(output)
}
# calculations
kcal_to_gj = 0.000004186798
top_frac = 0.01
energy_max = 3000
energy_subsistence_low = kcal_to_gj*2500*365
energy_subsistence_high = kcal_to_gj*4000*365
# inequality frontier with low subsistence energy
frontier_low = frontier_top_func(energy_subsistence_low, energy_max, 0.01)
# inequality frontier with high subsistence energy
frontier_high = frontier_top_func(energy_subsistence_high, energy_max, 0.01)
frontier = rbind(
frontier_low,
frontier_high[order(-top_share)],
frontier_low[1,]
)
frontier = rbind(
frontier_low,
frontier_high[order(-top_share)],
frontier_low[1,]
)
top_frontier = ggplot() +
geom_polygon(data = frontier, aes(x = energy_pc, y = top_share), fill = "black", alpha = 0.2, col = "black") +
geom_point(data = mod, aes(x = energy, y = power_top_1), size = 0.2, col = "grey30") +
scale_x_log10("Energy Use per Capita (GJ/year)",  breaks = c(2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000)) +
scale_y_continuous("Gini Index", breaks = seq(0,1, 0.1)) +
scale_color_gradientn(expression(beta), colours = rainbow(8), breaks = seq(0, 1, 0.1)) +
coord_cartesian(xlim = seq(2.5, 2000), ylim = c(0, 1)) +
theme_bw() +
theme(panel.border = element_rect(color = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.title = element_text(face="bold", size = rel(1), hjust = 0.5),
legend.title = element_text(hjust = 0.2),
legend.key.height = unit(1.2, "cm"),
legend.justification = c(0, 0.3),
axis.line = element_line(color = "black"),
axis.title.x=element_text(vjust=-0.3, size=rel(0.9)),
axis.title.y=element_text(vjust= 1.1, size=rel(0.9)),
axis.text.x = element_text(margin=margin(5,5,0,0,"pt")),
axis.text.y = element_text(margin=margin(3,5,0,3,"pt")),
axis.ticks.length = unit(-0.7, "mm"),
text=element_text(size = text.size, family="Times"))
top_frontier
View(frontier_low)
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
top_frontier
gA = ggplotGrob(gini_frontier)
# export
####################################################
setwd(dir)
gA = ggplotGrob(gini_frontier)
gB = ggplotGrob(top_frontier)
png("inequality_frontier.png", width = 7.5, height = 4,  units = 'in', res = 600)
grid.arrange(cbind(gA, gB, size = "first"))
dev.off()
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_result.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/energy_firm_power.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_result.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/kuznets.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/kuznets.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/kuznets.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/kuznets.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/kuznets.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/kuznets.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/kuznets.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/kuznets.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/kuznets.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/kuznets.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/kuznets.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/old/beta_function.R')
beta_func(2,10)
beta_func(1000,10^6)
beta_func(100,10^6)
beta_func(2,10)
beta_func(2,10)
beta_func(50,10^6)
beta_func(2,10)
beta_func(60,10^6)
beta_func(70,10^6)
beta_func(5,10)
beta_func(1000,10^6)
beta_func(3,10)
beta_func(1000,10^6)
beta_func(3,10)
beta_func(1000,10^6)
beta_func(3.2,10)
beta_func(1000,10^6)
beta_func(2,5)
beta_func(1000,10^6)
beta_func(3,5)
beta_func(1000,10^6)
beta_func(10,10)
beta_func(3,5)
beta_func(10^6,10^6)
source('~/Desktop/origin_inequality/Supplementary Material/Figures/beta_slave.R')
beta_slave
source('~/Desktop/origin_inequality/Supplementary Material/Figures/beta_slave.R')
beta_slave
source('~/Desktop/origin_inequality/Supplementary Material/Figures/beta_slave.R')
beta_slave
source('~/Desktop/origin_inequality/Supplementary Material/Figures/beta_slave.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/beta_merge.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/mod_empirical_scatter.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
source('~/Desktop/origin_inequality/Supplementary Material/Figures/inequality_frontier.R')
library(data.table)
library(here)
library(magrittr)
library(ggplot2)
library(gridExtra)
dir = here()
wd = gsub("Figures", "Empirical Data/Solomon", dir)
setwd(wd)
s = fread("solomon_data.csv")
n_slaves = s$n_slaves %>% na.omit()
# income
planter_income = s$planter_income %>% na.omit()
slave_expenses = s$slave_expenses %>% na.omit()
# house size
planter_house = s$planter_house %>% na.omit()
slave_house = s$slave_house %>% na.omit()
n = 1000
beta_income = 1:n
beta_house = 1:n
planter_income
i=1
slave_income_pc = slave_expenses/sample(n_slaves,1)
planter_income_pc = sample(planter_income, 1)
planter_slave_income_ratio = planter_income_pc / slave_income_pc
n_slave = sample(n_slaves, 1)
planter_power = n_slave + 1
beta_income[i] = log(planter_slave_income_ratio) / log(planter_power)
slave_house_sample = sample(slave_house, 1)
planter_slave_house_ratio = planter_house / slave_house_sample
beta_house[i] = log(planter_slave_house_ratio) / log(planter_power)
