2

## Dispatch → Meta → Gameplay

## 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[0]

y = xn[1]

x1 = x - np.matmul(np.linalg.inv(j(x,y)),fv(x,y))[0]

y1 = y - np.matmul(np.linalg.inv(j(x,y)),fv(x,y))[1]

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[0],xn[1],', which took',iteration,'iterations to produce. The 2-norm of the residual is',np.linalg.norm(fv(xn[0],xn[1])),'.')