In what level does inheritance affect the capital distribution of a population?

Find out with the help of Python

Project summary

Footnote

Person stats model

Creating a python module with functions that statistically calculate the following variables of a person

classPerson:def__init__(self,gender,social_class,age,income,capital,marital_status):globalid_countself.id=id_countself.social_class=social_classself.gender=genderself.age=ageself.income=incomeself.husband_wife=0self.children_id= []self.mother_id=0self.father_id=0self.children_status=0self.marital_status=marital_statusself.capital=capitalid_count+=1defget_married(self,id):self.marital_status=1self.husband_wife=iddefassign_child(self,id):self.children_status=1self.children_id.append(id)defassign_parents(self,father_id,mother_id):self.father_id=father_idself.mother_id=mother_iddefset_social_class(self,new_class):self.social_class=new_class#def update_income(self):defupdate_capital(self,costs):self.capital+=self.incomeself.capital-=costsdefincrement_age(self,increment):self.age+=incrementdefpass_capital(self,cap,siblings):self.capital+= (cap/siblings)def__repr__(self):return"<Person ID:%s, Gender:%s, Social Class:%s, Age:%s, Income:%s, Capital:%s, Husband/Wife ID:%s Children: %s Father: %s Mother: %s>"% (self.id,self.gender,self.social_class,self.age,self.income,self.capital,self.husband_wife,self.children_id,self.father_id,self.mother_id)

First, to create a person we need to make it a he or she. So we will be taking a random binary probability of thesex being 0 for males and 1 for females.

importnumpyasnpnp.random.seed(1)gender=np.random.randint(2)gender

1

Next we will define that person'sage based on his/her gender and according to the UN's data (http://data.un.org/Data.aspx?d=POP&f=tableCode%3A22)

The age data follows this dataframe:

importpandasaspdages=pd.read_csv("data/ages.csv")ages_df=pd.DataFrame(ages)ages_df

.dataframe tbody tr th {    vertical-align: top;}.dataframe thead th {    text-align: right;}

The Age column contains data in bins of 5, but in theget_age_percents() function, we convert the bins to bins of 5 and then calculate theprobability density function of the data manually.

defget_age_percents():ages=pd.read_csv("data/ages.csv")ages_df=pd.DataFrame(ages)starting_columns= ["Total","0 - 4","5 - 9","10 - 14","15 - 19","20 - 24","25 - 29","30 - 34","35 - 39","40 - 44","45 - 49","50 - 54","55 - 59","60 - 64","65 - 69","70 - 74","75 - 79","80 - 84","85 - 89","90 - 94","95 - 99","100 +"]final_columns= ["Total","0 - 9","10 - 19","20 - 29","30 - 39","40 - 49","50 - 59","60 - 69","70 - 79","80 - 89","90 - 99","100 +"]# ages_1 = ages_df["Value"].loc[(ages_df["Age"].isin(starting_columns[5:7])) & (ages_df["Sex"] == "Male")]# print ages_1ages_m= []ages_f= []t=0foriinrange(1,11):ages_m .append(sum(ages_df["Value"].loc[(ages_df["Age"].isin(starting_columns[i+t:i+t+2]))& (ages_df["Sex"]=="Male")]))ages_f .append(sum(ages_df["Value"].loc[(ages_df["Age"].isin(starting_columns[i+t:i+t+2]))& (ages_df["Sex"]=="Female")]))t=t+1total_m=sum(ages_df["Value"].loc[(ages_df["Age"]=="Total")& (ages_df["Sex"]=="Male")])hundred_m=sum(ages_df["Value"].loc[(ages_df["Age"]=="100 +")& (ages_df["Sex"]=="Male")])total_f=sum(ages_df["Value"].loc[(ages_df["Age"]=="Total")& (ages_df["Sex"]=="Female")])hundred_f=sum(ages_df["Value"].loc[(ages_df["Age"]=="100 +")& (ages_df["Sex"]=="Female")])ages_m .append(hundred_m)ages_f .append(hundred_f)#Get percents of agesimportnumpyasnpages_m=np.array(ages_m)ages_f=np.array(ages_f)ages_m_percent= (ages_m/total_m)*100ages_f_percent= (ages_f/total_f)*100ages_m_percent= ["%.2f"%itemforiteminages_m_percent ]ages_f_percent= ["%.2f"%itemforiteminages_f_percent ]returnages_m_percent,ages_f_percent

When we have the p.d.f of the ages, we will use aget_age function to get a random age based on the gender of the individual like this:

importagesage=ages.get_age(gender)age

8

Then we have to find the Social Class of the individual which follows the Thompson, Hickey 2005 model I described in the table above.

Concretely:

importscipy.statsasstsdefget_class():xk=np.arange(5)pk= (0.2,0.32,0.32,0.15,0.01)custm=sts.rv_discrete(name='custm',values=(xk,pk))returncustm.rvs()+1social_class=get_class()social_class

1

I used R and some data collected from IPUMS.org, to achieve my goal.Then I used some R code to divide the data instances to the classes based on the proposed model[1], after I analysed the income to the gender and age of the people as shown in Figure 1.

Figure 1:Age to Income Scatterplot (Blue:Male, Red: Female)

With the following R code I divided the data I gathered from IPUMS.org to social classes based on the model proposed by Thompson and Hickey[1].

df<- read.csv(file="cps_00005.csv",header=TRUE,sep=",")#Remove nan, d$INCTOT == 99999999, d$INCTOT == 99999998,df<- na.omit(df)df<-df[!(df$FTOTVAL==99999999|df$FTOTVAL==99999998),]df<-df[!(df$INCTOT==99999999|df$INCTOT==99999998),]df<-df[order(df$FTOTVAL),]#We get the 5 classes. According to William Thompson & Joseph Hickey, 2005we havelc= (floor(.2*nrow(df)))wc= (floor(.52*nrow(df)))lmc= (floor(.84*nrow(df)))umc= (floor(.99*nrow(df)))uc= (floor(.01*nrow(df)))lower_class=df[1:lc,]working_class=df[(lc+1):wc,]lower_middle_class=df[(wc+1):lmc,]upper_middle_class=df[(lmc+1):umc,]upper_class=df[(umc+1):nrow(df),]#print(nrow(lower_class)/nrow(df)*100)  -> 20%#print(nrow(working_class)/nrow(df)*100) -> 32%#print(nrow(lower_middle_class)/nrow(df)*100) -> 32%#print(nrow(upper_middle_class)/nrow(df)*100) -> 15%#print(nrow(upper_class)/nrow(df)*100) -> 1%print(fitdist(lower_class$INCTOT,distr="unif",method="mme"))print(fitdist(working_class$INCTOT,distr="norm",method="mme"))print(fitdist(lower_middle_class$INCTOT,distr="norm",method="mme"))print(fitdist(upper_middle_class$INCTOT,distr="norm",method="mme"))print(fitdist(upper_class$INCTOT,distr="gamma",method="mme"))

Figure 2:Lower Class Distribution fit

Figure 3:Working Class Distribution fit

Figure 4:Lower Middle Class Distribution fit

Figure 5:Upper Middle Class Distribution fit

Figure 6:Upper Class Distribution fit

Python

**Lower Class**: Uniform**Working Class**: Normal**Lower Middle Class**: Normal**Upper Middle Class**: Normal**Upper Class**: Gamma

defget_income(social_class):#Got classes from https://en.wikipedia.org/wiki/Household_income_in_the_United_States  Social_Class#lower_class = 1 (14-20%) , working_class =  2 (20 -52%), lower_middle_class = 3 (52-84%), upper_middle_class= 4 (84-99%), upper_class = 5 (99-100%)if(social_class==1):inc=np.random.uniform(-3744.26,15749.78)elif(social_class==2):inc=np.random.normal(15879.603,9528.182)elif(social_class==3):inc=np.random.normal(27564.15,19854.43)elif(social_class==4):inc=np.random.normal(55957.25,48711.24)else:inc=np.random.gamma(1/0.6346247,1/0.000002982905)returninc# ro.r.source("income.R")# ro.r('''source("income.R")''')# inc = ro.globalenv['income']# x = inc(18,2,3)# return xincome=get_income(social_class)if(age<15):income=0income

0

Capital

importstatscapital=stats.get_capital(age,social_class,income)capital

0.0

fromfamilyimport*marital_status=family.get_marital_status(gender,age)#1 for yes, 0 for nomarital_status

0

person=Person(gender,social_class,age,income,capital,marital_status)printperson

<Person ID:1, Gender:1, Social Class:1, Age:8, Income:0, Capital:0.0, Husband/Wife ID:0 Children: [] Father: 0 Mother: 0>

Generation 0

defmodel(N):population= []while(len(population)<N ):# Create a person, if it aged more than 18, run a probability of it being married, else just give his parents# If he/she is create a partner for him/her and give them a child# If it's a male the social class has chances of being lower, females on the other hand, tend to marry upgender=np.random.randint(2)social_class=get_class()age=ages.get_age(gender)marital_status=0partner=0if(age<15):inc=0else:inc=get_income(social_class)cap=get_capital(age,social_class,inc)person=Person(gender,social_class,age,inc,cap,marital_status)#print ( "Created person-id %d" % person.id)marital_status=get_marital_status(gender,age)if(marital_status==1):#Create another human being with different gender, if male then social class is same or lower,#if female marry only to same classpartner=create_partner(person)#Give parents to partnerpartner_mother,partner_father=create_parents(partner)if(partner_mother!=0):partner.mother_id=partner_mother.idpopulation.append(partner_mother)if(partner_father!=0):partner.father_id=partner_father.idpopulation.append(partner_father)#Marry Thempopulation.append(partner)person.get_married(partner.id)partner.get_married(person.id)#Check for and give them Childrenif(person.gender==1):m=personf=partnerelse:m=partnerf=person#m = mother, f = father of childif(get_children_status(m,f)):no_of_children=get_number_of_children(m,f)if(no_of_childrenisNone):no_of_children=1for_inrange(no_of_children):child=create_child(person,partner,1)if (child!=0):partner.assign_child(child.id)person.assign_child(child.id)child.assign_parents(f.id,m.id)population.append(child)mother,father=create_parents(person)population.append(person)if(mother!=0):person.mother_id=mother.idpopulation.append(mother)if(father!=0):person.father_id=father.idpopulation.append(father)returnpopulation

Results

importnumpyasnpnp.random.seed(1)population1=model(1000)population2=model(2000)population3=model(5000)

defplot_capital_to_class(population,N):# Check the capital to social class distributioncapital= []total_cap=0class_size= [0]*5class_cap= [0]*5tot_inc=0class_inc= [0]*5forpinpopulation:capital.append(p.capital)total_cap+=p.capitaltot_inc+=p.incomeclass_size[p.social_class-1]+=1class_cap[p.social_class-1]+=p.capitalclass_inc[p.social_class-1]+=p.incomeclass_inc=np.array(class_inc)print"Income:"print  (class_inc/tot_inc)*100class_cap=np.array(class_cap)print"Capital:"arr= (class_cap/total_cap)*100printarrclass_size=np.array(class_size)print"Population:"pop= (class_size/float(N))*100printpopprint"Capital:Population ratio"print (arr/pop)importmatplotlib.pyplotaspltind=np.arange(5)plt.bar(ind,arr,align='center',width=0.5)plt.xticks(ind,("Lower Class","Working Class","Lower Middle Class","Upper Middle Class","Upper class") )plt.ylabel("Capital")plt.title("Capital distribution per social class (Simulation)")plt.tight_layout()plt.ylim(0.0,60.0)plt.show()returncapital

cap1=plot_capital_to_class(population1,1000)cap2=plot_capital_to_class(population2,2000)cap3=plot_capital_to_class(population3,5000)

Income:[ 3.25754056 17.64907108 27.34955917 36.03464325 15.70918595]Capital:[ 2.94931656 14.77192988 27.60128711 35.05683699 19.62062946]Population:[16.9 33.3 31.1 18.   1.4]Capital:Population ratio[ 0.17451577  0.4436015   0.88750119  1.94760205 14.01473533]

Income:[ 3.60685866 17.54332042 32.38809354 26.961163   19.50056438]Capital:[ 2.92336124 14.31667082 31.78975533 29.25399279 21.71621981]Population:[17.95 32.5  33.4  14.65  1.55]Capital:Population ratio[ 0.16286135  0.44051295  0.95178908  1.99685958 14.01046439]

Income:[ 4.81334394 17.38036787 33.43219905 28.23697775 16.13711139]Capital:[ 3.92188299 14.30665683 32.61451713 29.31294689 19.84399616]Population:[21.52 29.48 33.76 14.1   1.16]Capital:Population ratio[ 0.18224363  0.48530044  0.96606982  2.0789324  17.10689324]

defgini_index(capital):capital=capital.flatten()ifnp.amin(capital)<0:# Values cannot be negative:capital-=np.amin(capital)# Values cannot be 0:capital+=0.0000001# Values must be sorted:capital=np.sort(capital)# Index per array element:index=np.arange(1,capital.shape[0]+1)# Number of array elements:n=capital.shape[0]# Gini coefficient:return ((np.sum((2*index-n-1)*capital))/ (n*np.sum(capital)))printgini_index(np.array(cap1))printgini_index(np.array(cap2))printgini_index(np.array(cap3))

0.196175876285178350.182059035083937060.12080484386107655

Simulation

defmarriage_query(person,partner):if(partner.marital_status==0):if(partner.gender!=person.gender):if(person.gender==0):#if(partner.age in xrange((person.age/2) + 7, person.age + 8)):if(partner.social_class<=person.social_class):person_marital_status=get_marital_status(person.gender,person.age)partner_marital_status=get_marital_status(partner.gender,partner.age)if((person_marital_status)& (partner_marital_status)):return1else:return0#else:#    return 0else:#if(partner.age in xrange((person.age) , person.age + 21)):person_marital_status=get_marital_status(person.gender,person.age)partner_marital_status=get_marital_status(partner.gender,partner.age)if((person_marital_status)& (partner_marital_status)):return1#    else:#return 0else:return0else:return0# This function updates the people's classes after each iterationdefupdate_class(population):# Sort ascending based on capitalpopulation.sort(key=lambdax:x.capital)# Get classesT=len(population)l=int(0.2*T)w=int(0.32*T)lm=int(0.32*T)um=int(0.15*T)u=int(0.01*T)# Set new social classfori,pinenumerate(population):if(i<=l):p.social_class=1elif(i<=l+w):p.social_class=2elif(i<=l+w+lm):p.social_class=3elif(i<=l+w+lm+um):p.social_class=4else:p.social_class=5returnpopulationdefrun_simulation(population,years):foryearinrange(years):printlen(population)forpinpopulation:if(check_death(p)==0):p.increment_age(1)# TODO:#Create a costs variable from data and initialize that too#Then adjust income and capital according to the already given data?! or make them interact?!#ADJUST Incomeif(p.age==15):p.income=get_income(p.social_class)else:#p.update_income()p.income+=p.income*0.01costs=0p.update_capital(0)#If not married, run probability of him getting get_marriedif(p.marital_status==0):#Find a potential marriage partnerpartners= [partnerforpartnerinpopulationif(marriage_query(p,partner))]if(partners):prtnr=random.choice(partners)prtnr.get_married(p)p.get_married(prtnr)#If married  run probability of giving birthif((p.marital_status==1)):forpersoninpopulation:if(person.id==p.husband_wife):if(person.gender==0):husband=personwife=pelse:wife=personhusband=pchild_status=family.get_children_status(wife,husband)if(child_status==1):child=family.create_child(husband,wife,0)husband.assign_child(child.id)wife.assign_child(child.id)child.assign_parents(husband.id,wife.id)population.append(child)breakelse:#If that person has children, then pass its capital to his/her kidschildren= [childforchildinpopulationif (child.idinp.children_id) ]no_of_children=len(children)forcinchildren:c.pass_capital(p.capital,no_of_children)population= [personforpersoninpopulationifp.id!=person.id]population=update_class(population)cap=plot_capital_to_class(population,len(population))printgini_index(np.array(cap))

run_simulation(population1,20)

10071462188122772650304334043723402843794688498252635543582060936358662168567104Income:[ 1.96802942 20.26639224 34.34630328 30.3449897  13.07428536]Capital:[-2.49912499  9.9415439  33.48622741 41.4580972  17.61325647]Population:[20.01363327 31.99727335 31.99727335 14.99659168  0.99522836]Capital:Population ratio[-0.12487113  0.31069972  1.04653378  2.7645013  17.69770359]

0.18507320299921423

Can we do better?