2

## A treatise on the mathematical modeling of income inequality

Income inequality is one of the most confusing gameplay statistics. We know the average income, and we know the average incomes of the bottom 10% and the top 10%, but we do not know about anything in the middle. This is important, as the median income is much more important than the average in determining the quality of life in your nation. First, it should be noted that this is a statistical examination based on a mathematical model. This is not a guide on how to raise or lower income inequality in your country.

The Model

Here, we will be using an income inequality model described by two parameters, χ and α. The model is I = χ/(1-P)^α, where I is income, P is proportion of the population are below you in income level, χ is the minimum income for a full-time worker, and α is the inequality parameter. It is closer to 1 with higher inequality, and closer to zero with lower inequality. Using this equation, we can calculate many of the parameters shown in NationStates, and some not shown, but that are important, using χ and α, by applying or integrating the formula shown above. It should be noted that this formula assumes a minimum income level, which is simply the constant χ. These derived formulas are shown in appendix A.

This model comes from a Stanford paper on modelling inequality (see references), the description of which is as follows. "It is well known that the upper tail of the income distribution follows a power law. One way of thinking about this is to note that income inequality is fractal in nature, as we document more carefully below. In particular, the following questions all have essentially the same answer: What fraction of the income going to the top 10 percent of earners accrues to the top 1 percent? What fraction of the income going to the top 1 percent of earners accrues to the top 0.1 percent? What fraction of the income going to the top 0.1 percent of earners accrues to the top 0.01 percent? The answer to each of these questions - which turns out to be around 40 percent in the United States today - is a simple function of the parameter that characterizes the power law." (Jones and Kim, 1786)

Applications

As formulas for the average income (which we will call μ) and the average income of the rich (top 10%) divided by the average income of the poor (bottom 10%) (we will call this massive number ι) can be derived, and these formulas are the ones provided by NationStates, we can use them to derive the values of χ and α, and thus the rest of the formulas in appendix A can be expressed in terms of ι and μ. Unfortunately, the only way to compute these numbers is using numerical methods (such as newton's method), as the system of equations generated cannot be inverted using ordinary mathematical functions. The code required to do this is in appendix C.

Results

As an example, Rho Ophiuchi has an average income of 88,084, and has ι, the "inequality score" equal to 14.2647. Using newton's method, we can calculate that χ = 28,383.9 and α = 0.677764. In contrast, the Municipalities of Antarctica has an average income of 199,607, and an "inequality score" equal to 1.20169. Using newton's method, we can calculate that χ = 177,408 and α = 0.111211. Found in appendix B are the other parameters, calculated form the derived χ and α values.

Notice how both the average income of the poor and the average income of the rich calculated using this model are higher than the displayed value in NationStates, when calculated from the income inequality and mean income. This is due to much less complex modelling being used in the internals of NationStates itself, and the values displayed here are likely more realistic.

References

Jones, Charles, and Jihee Kim. "A Schumpeterian Model of TopIncome Inequality." Journal of Political Economy, 2018, vol. 126, no. 5. https://web.stanford.edu/~chadj/inequality.pdf.

Appendix A

 Value Type Formula using χ and α Formula using ι and μ Average income χ/(1-α) μ Average income of poor (bottom 10%) (1-0.9^(1-α))∙10χ/(1-α) N/A Average income of rich (top 10%) 10χ∙0.1^(1-α)/(1-α) N/A Lowest income for a full-time job χ N/A 10th percentile income χ/0.9^α N/A Median income χ/0.5^α N/A 90th percentile income χ/0.1^α N/A 99th percentile income χ/0.01^α N/A Income of rich divided by income of poor (0.1^(1-α))/(1-0.9^(1-α)) ι

Appendix B

 Value Type Rho Ophiuchi: (ι, μ) = (14.2647, 88,084) The Municipalities of Antarctica: (ι, μ) = (1.20169, 199,607) α 0.677764 0.111211 χ 28,383.9 177,408 Average income 84,084.2 199,606 Average income of poor (bottom 10%) 29,403.5 178,433 Average income of rich (top 10%) 419,432 257,861 Lowest income for a full-time job 28,383.9 177,408 10th percentile income 30,484.9 179,499 Median income 45,404.5 191,624 90th percentile income 135,156 229,184 99th percentile income 643,576 296,090 Income of rich divided by income of poor 14.2647 1.20169

Appendix C

i = 14.2647#this is inequality
m = 88084.0#this is average income

import numpy as np

def f1x(x,y):
f1x = 0
return f1x
def f1y(x,y):
f1y = (((0.1**(1-y)) + (0.1**(1-y) * y * np.log(10))) / (1-(9/10)**(1-y))) + (((0.09**(1-y) * y * np.log(10 / 9))) / (1-(9/10)**(1-y)) ** 2)
return f1y
def f1(x,y):
f1 = (0.1**(1-y))/(1-0.9**(1-y)) - i
return f1
def f2(x,y):
f2 = x/(1-y) - m
return f2
def f2x(x,y):
f2x = 1/(1-y)
return f2x
def f2y(x,y):
f2y = x/(1-y)**2
return f2y
def j(x,y):#jacobian
j = np.array([[f1x(x,y),f1y(x,y)],[f2x(x,y),f2y(x,y)]])
return j
def fv(x,y):#vector function
fv2 = np.array([f1(x,y),f2(x,y)])
return fv2
def newton_i(xn,j,fv):# iteration of multivariable newton method
x = xn
y = xn
x1 = x - np.matmul(np.linalg.inv(j(x,y)),fv(x,y))
y1 = y - np.matmul(np.linalg.inv(j(x,y)),fv(x,y))
xn1 = np.array([x1,y1])
return xn1
xn = np.array([0.5,0.5])#starting point
x = [xn]# all x
xnm1 = [0,2]#previous estimate of x
iteration = 0
while iteration < 500:
xnm1 = xn
xn = newton_i(xn,j,fv)
iteration += 1
x.append(xn)
error = np.linalg.norm(xn-xnm1)
print('The solution is',xn,xn,', which took',iteration,'iterations to produce. The 2-norm of the residual is',np.linalg.norm(fv(xn,xn)),'.')