A Mathematical Model of Universal Basic Income and Its Numerical Simulations

In this paper, an elementary mathematical model describing the introduction of a universal basic income in a closed market society is constructed. The model is formulated in terms of a system of nonlinear ordinary differential equations, each of which gives account of how the number of individuals in a certain income class changes in time. Societies ruled by different fiscal systems (with no taxes, with taxation and redistribution, with a welfare system) are considered and the effect of the presence of a basic income in the various cases is analysed by means of numerical simulations. The main findings are that basic income effectively acts as a tool of poverty alleviation: indeed, in its presence the portion of individuals in the poorest classes and economic inequality diminish. Of course, the issue of a universal basic income in the real world is more complex and involves a variety of aspects. The goal here is simply to show how mathematical models can help in forecasting scenarios resulting from one or the other policy.


Introduction
A lively debate on the appropriateness of introducing some form of universal basic income (UBI) has been ongoing in several countries for a long time (for concepts, methodologies, practice and much more see for example the recent papers [1][2][3] and books [4,5] which in turn also provide several further references). Additionally, economic difficulties which grow stronger in times of crisis spark the interest for the subject. The problem is in fact a complex one which involves not only political and economic, but also sociological and psychological aspects. Engaging in the discussion to endorse one or the other position towards is certainly not the goal of this paper. Rather, the goal is to illustrate, through a technical approach, how numerical simulations of mathematical models can help predict long-term effects of different policies, in particular related to UBI. This can then help to explore different possible scenarios emerging in correspondence of different political programs. Of course, a serious and attentive contribution (not to be found here) cannot ignore and should definitely take into account several related socio-economic facets and elements.
The investigation and the simulations here described have roots in a mathematical "micro-to-macro" model proposed by this author in [6] and later generalised and studied in other works (see e.g., [7,8]). The model describes the evolution in time, in the case of constant global income, of the income distribution over a population, resulting out of a whole of economic interactions which take place between individuals of the population. Individuals are grouped in a finite number of classes distinguished by their income and exchange money according to frequency and rules first suggested in [7]. A fiscal system with progressive taxation together with a redistribution process, integrated in paper [8] with a welfare policy, was considered in the mentioned papers. The mathematical framework developed for the purpose is provided by a set of nonlinear ordinary differential equations, each of which gives account of how the number of individuals in a certain income class changes in time.
The analysis of the dynamics, based to a great extent on a large number of simulations performed in correspondence to different choices of initial conditions and different values of the model parameters allows the establishment of the existence, for fixed parameter values and for fixed global income, of a unique asymptotic (long-run) stationary distribution. Such distribution has been intended as a dynamic equilibrium, because, albeit the number of individuals in each class remains constant from a certain time on, single individuals can move from one class to another. It has also been shown in [7] that suitable choices of the parameters allow the recovery of distributions which are well in agreement with empirical data. Among other features, economic inequality, typically measured by means of the Gini coefficient, has been the object of particular attention.
The mentioned framework turns out to be suitable also to investigate the effect of the introduction of a universal basic income. Indeed, if a basic income provides a guaranteed financial grant to everybody, then everybody has a minimal income which is different from zero. This can be expressed in the context of the framework by assuming that the lowest income class is empty, or the lowest income classes are empty. Specifically, distributions of the population can be considered for which the lowest income class/classes are empty. The interest is then to examine the differences between long-run income profiles emerging in the case in which no basic income allowance is provided and in the case in which a basic income allowance is provided, under the condition of a same global income (i.e., assuming that the "richness" of the population is the same), in correspondence to different government schemes: This is what we do in this paper. Actually, in carrying out this project, we introduce some modifications in the original model which let it become even more appropriate to represent real situations. Basically, the intervals characterising the income classes are here assumed to be of growing length in passing from the lowest income to the highest income one. This requires that also some terms in the expressions representing taxes and welfare be modified.
The main results from the numerical simulations are the following ones: (1) when basic income is introduced in a specific model, the portion of individuals in the poorest classes and economic inequality diminish; (2) if two societies are compared which are ruled by the same fiscal system, but have a different total income, the portion of individuals belonging to a certain class at the bottom of the income ladder is smaller in the "richer" society; (3) the effect of reduction of economic inequality is more pronounced in poorer societies and, in order, is more pronounced when no taxation and redistribution exist, when taxation with a unique tax rate exist, when a progressive taxation system exist, when taxation, redistribution and welfare exist.
The rest of the paper is organised as follows. In the next section, we introduce the modified model. The results of some representative simulations (among the many executed), obtained by numerically solving, as detailed below, systems of differential equations are presented and discussed in Section 3. These include the representation of the income distribution as well as the representation of the "richness" of the various classes in the long run in various cases, together with a comparison of the Gini coefficients corresponding to these distributions. Section 4 contains a summarising discussion together with a critical analysis. Details on the numerical solution methods are contained in Appendix A. To give a short preview here, the numerical solutions were found with two different methods, essentially leading to the same results. First, a numerical solver from the software Wolfram Mathematica 12.3.0.0 [9] was employed to solve a suitably reduced equation system arising from a differential-algebraic system of equations. Then, a second method was used that is in fact a combination of Euler and Romberg methods and that does not require the reduction of the equations.

The Model
We introduce here a modified version of a model first proposed in [6] and then further developed in [7,8]. For the convenience of the reader we explicitly point out since the beginning that the novelties in the present version consist of a different subdivision of income ranges compared to the previous model and in the formalisation of the allocation of a basic income through the assumption of emptiness of the poorest class/classes (see below).
We start by considering n + 1 numbers 0 = ρ 0 < ρ 1 < . . . < ρ n , with ρ n intended to be large enough to represent an upper limit (that cannot be exceeded) for the maximal conceivable income of an individual of a given society. In particular, differently from what was done in the mentioned previous works, we will take ρ k+1 − ρ k = α (ρ k − ρ k−1 ) with α > 1, for k = 1, . . ., n − 1. In this way, n intervals can be considered, whose length increases with k. These intervals identify n income classes. We emphasise that this provides a more realistic description. Indeed, in this way, one can, by suitably choosing α, describe societies in which the richest individuals (those belonging to the n-th income class) have an average income which is order of magnitude larger that of the poorest individuals (those in the first class). We suppose for simplicity that the fraction of individuals whose income falls between ρ k−1 and ρ k (more precisely, in [ρ k−1 , ρ k )) is the same in all points of the interval. We denote by x k (t), with x k : R → [0, +∞), this fraction at time t, and denote r k = (ρ k+1 − ρ k )/2 the average income of the k-th class.
The overall dynamics, which determines at each instant the income distribution, results from a whole of economic exchanges between individuals. These exchanges represent both earnings (due for example to the provision of some good or service) and losses (due to the payment in exchange for some good or service). They induce "small" movements of the individuals undergoing the exchanges from a class to a closer one. Beside such earnings and losses, which occur as direct interactions of pairs of individuals, others can be thought to take place, representing the benefit from the tax revenue redistribution (provided through education, health care and similar supports) and the payment of taxes accompanying any gain respectively. Indeed, we postulate the existence of a fiscal system based on a progressive taxation. Specifically, denoting by τ k ∈ [0, 1) the tax rate for the k-th class, we assume where τ min = τ 1 and τ max = τ n respectively denote the minimum and maximum tax rates. The way how taxation (together with redistribution) is taken into account by the algorithm described below is as follows. Every time − an individual of the i-th class pays an amount S (S (r k+1 − r k ) for all k = 1, . . ., n) to one of the j-th class and − the latter one pays to the government a tax on her earnings according to her income class (her tax rate being τ j ) and − the government redistributes the tax revenue to the population. This process is equivalently described by the following actions: − the i-individual pays an amount S (1 − τ j ) to the j-individual and − the i-individual pays an amount S τ j which is divided among all individuals.
Further details can be found in [6], where this mechanism was first introduced. We notice here that due to constraints inherent to the model, individuals of the n-th class are excluded by the redistribution and cannot in general receive money which would let them exit from the highest income class. Analogously, individuals of the poorest class cannot pay money to end up with a negative income.
What we want to emphasise here is that this elementary model treats society as a collection of individuals who interact directly (through money exchange) or indirectly (through taxation and redistribution). No central government is effectively present in the algorithm.
A system of n nonlinear ordinary differential equations, one for each income class in which the population is divided, describes the variation in time of the fraction of individuals in these classes: The meaning and the form of the C k ij 's and T k [ij] 's appearing on the r.h.s. of (1) are detailed next. For the convenience of the reader, we repeat here a description necessarily similar to one already given in [6][7][8]. A warning is however in order that some terms below are here different than the corresponding ones in those papers.

−
Each C k ij ∈ [0, +∞) expresses the probability density that an individual of the i-th class will belong to the k-th class after a direct interaction with an individual of the j-th class. The identity ∑ n k=1 C k ij = 1 has hence to hold true for any fixed i and j; − Each T k [ij] : R n → R expresses the variation density in the k-th class due to an interaction between an individual of the i-th class with an individual of the j-th class. The functions T k [ij] are continuous and need to satisfy ∑ n k=1 T k [ij] (x) = 0 for any fixed i, j and x ∈ R n .
To specify the form of the C k ij 's and T k [ij] 's we still need some more definitions: − Denote by p i,j (for i, j = 1, 2, . . ., n) the probability that in an interaction between an individual of the i-th class with an individual of the j-th class, the one who pays is the former one. Specifically, we assume with the exception of the terms p h,h = r h /2r n for h = 2, . . ., n − 1, p i,1 = r 1 /2r n for i = 2, . . ., n, p n,j = r j /2r n for j = 1, . . ., n − 1, p 1,j = 0 for j = 1, . . ., n and p i,n = 0 for i = 1, . . ., n. − Denote by S the minimal amount of money that individuals may exchange; − Denote by w h ∈ R (for h = 1, 2, . . ., n) n coefficients playing the role of weights, suitable to incorporate welfare in the redistribution terms. Specifically, we assume here (this expression is different from that one adopted in [8]) (2) one easily checks that if β = 1/2, then w h takes the same value, 1/2, for each h = 1, . . ., n, whereas if β > 1/2, then w 1 = β, w n = 1 − β and w h is decreasing as a function of h.
Referring to [6] for technical details on the range of the indices in the formulas below (essentially, all expressions are to be thought as present only for meaningful values of the indices), we can now write the expressions of the C k ij 's and T k [ij] 's as follows: the only nonzero elements C k ij 's are those of the form describing the advancement from a class to the subsequent one, due to the benefit of tax revenue redistribution and describing the retrocession from a class to the preceding one, due to the payment of some tax. The symbol δ i,j is nothing but the Kronecker delta.

Properties of the Model
We point out here that independently of some differences between the expressions of certain terms in the present version of the model and those appearing in versions discussed in previous works, the following properties hold true.
The existence for any x 0 = (x 01 , . . . , x 0n ) with x 0k ≥ 0 for all k = 1, 2, . . ., n and ∑ n k=1 x 0k = 1, of a unique solution x(t) = (x 1 (t), . . ., x n (t)) of (1), satisfying x(0) = x 0 has been proved in [6]. In particular, it has been proved there that such a solution is defined for all t ≥ 0 and satisfies x k (t) ≥ 0 for k = 1, 2, . . ., n, as well as Additionally, it has been proved in [6] that the function µ : {x ∈ R n : ∑ n k=1 x k = 1} → R, defined as µ(x) = ∑ n k=1 r k x k and which represents the global income is a first integral for (1): Lastly, an extensive number of numerical simulations, also carried out with different parameter values, provides evidence that for fixed model parameters, in correspondence to any conceivable value of the global income µ, a unique stationary solution of (1) exists to which all solutions x(t) = (x 1 (t), . . ., x n (t)) satisfying x(0) = x 0 with µ(x 0 ) = µ tend as t → ∞.

Allocation of a Basic Income
The way a basic income system is incorporated into the model is as follows. We suppose that the bottom class or classes in the income ladder are empty. The construction of the model with parameters chosen as above implies that if there are not individuals belonging to the poorest income class [respectively, to the poorest m income classes at the bottom of the income ladder] at a time that we consider to be the initial time, t = 0, then the same also holds true for all t ≥ 0. This follows from the technical constraint included in the model and recalled above, expressed by the requirement that p 1,j = 0 for j = 1, . . ., n. According to this constraint, individuals in the poorest class are supposed to never pay. By contrast, they could end up in a class with "negative" income, something that we do not consider here. Yet, it is possible that individuals of other classes, characterised by higher income than the poorest one, worsen in time their status and end up in this last one. In short, the mentioned constraint provides a kind of barrier to the transition to poverty beyond the poorest class.
According to the model, if a society ruled by a certain fiscal system is divided into n income classes, all of which, from the first one to the n-th one, are effectively populated, the income classes effectively populated of the same society when basic income provision is introduced are those going for example from the second one to the n-th one, or those going from the third one to the n-th one [or those going from the (m + 1)-th one to the n-th one, if the poorest m income classes are to be empty].

Results from Numerical Simulations
To find numerical solutions of the equation system (1), we obtain insight on its dynamics and make comparisons between the emerging scenarios in the presence of different fiscal systems, and in particular, in the presence of a basic income allocation, we need to fix the values of the parameters entering in the model.
We want to be able to acquire a graphical representation of the income distribution, somehow resembling, albeit within a discrete approach, one of a continuous distribution. We hence choose here a rather high number of classes, n = 15. We then choose ρ 1 = 1 and α = 1.3, which yields (rounding off numbers) Furthermore, whereas for the cases (I) and (I I) we assume β = 0.50 (i.e., no additional welfare), for the case (I I I) we consider both the possibility that β = 0.50 (no additional welfare) and β = 0.67 (additional welfare). We also fix S = 1.
We consider, in correspondence of each of the cases (I), (I I), (I I I) above (and in the case (I I I) in correspondence of the two versions, with β = 0.50 and β = 0.67), relative to societies ruled by different fiscal systems, the time evolution of system (1) for two different values of the total income, here µ = 16 and µ = 24 (to exemplify a "poor" and a "rich" population), in three cases: -(i) when no basic income allocation is present, -(ii) with basic income allocation and the class with average income r 1 empty, -(iii) with basic income allocation and the classes with average income r 1 , r 2 empty.
Summarising, the cases here considered are those in the following list, where we also introduce a suitable label for each of them, to be then employed in the captions of the panels of some figures: As recalled in Section 2.1, for any choice of model parameters, and any conceivable value of the total income µ, a unique stationary income distribution exists to which solutions tend in the long run.
The "asymptotic" distributions relative to the mentioned cases are displayed in the Figures 1 and 2. More specifically, the three histograms in each panel of Figure 1 display the fraction of individuals in each of the 15, 14, 13 non-empty income classes relative to the model quoted in the caption of the panel, in the three cases (i), (ii), (iii). A similar representation of the long-run distributions, in which account is taken also of the width of the income intervals characterising the income classes, is given in Figure 2.
One may easily notice that when basic income is introduced in a specific model, the fraction of individuals belonging to the poorest classes diminishes. Even if in some panels in Figure 1 the first column in the second and third histogram is higher than that one in the first histogram, one should remember that the fraction of individuals in the poorest class of (iii) (i.e., the class with average income r 3 ) must be compared with-and it is always smaller than-the fraction of individuals belonging to the two poorest classes of (ii) (the classes with average income r 2 an r 3 ) and, in turn, this latter one must be compared with-and it is always smaller than-the fraction of individuals belonging to the three poorest classes of (i) (the classes with average income r 1 , r 2 an r 3 ).
This leads one to think that for each of the examined cases economic inequality is smaller provided basic income allocation is present.
Additionally, a comparison between stationary distributions emerging in the same taxation model (namely for fixed values of τ min , τ max and β), but in correspondence to different total income values is interesting. It can be achieved by looking at the histograms contained in two panels on the same row, one in the right column and the other in the left column, e.g., in Figure 1a,b or in Figure 1c,d, and so on. By looking at the left and right panels in Figure 1 one may observe that, if µ = 24 (larger total income), the fraction of individuals belonging to a class at the bottom of the income ladder is systematically smaller than that one of the same income class, of the same model, with (smaller total income) µ = 16. Here and in the Figure 2 also notice that the histograms on the left and on the right panel column are scaled differently.
To obtain an explicit measure of the income inequality characterising the various cases, we calculate the Gini coefficient. Its definition, next recalled, is based on the Lorenz curve. This curve plots on the y axis the fraction of the cumulative total income of the lowest-income fraction of the population, represented on the x axis. In other words, if the population and the total income are normalised to 1, the poorest cumulative fraction x of the population possesses the cumulative fraction y(x) of the total money. The line at 45 degrees represents perfect income equality. The Lorenz curve satisfies y(0) = 0, y(1) = 1 and y(x) ≤ x for all x ∈ [0, 1]. The Gini coefficient is defined as the ratio A/B of the area A that lies between the line of equality and the Lorenz curve over the total area B of the triangle which lies under the line of equality. Consequently, it takes values in [0, 1]. The larger the Gini coefficient is, the greater the income inequality is.
Each panel of Figure 3 shows the Lorenz curves relative to the income distribution of the model quoted in its caption in the three cases (i), (ii), (iii). All graphs in the eight panels confirm that the Gini coefficient for a society, and equivalently also the economic inequality, diminishes provided a basic income is introduced. This can be seen also quantitatively: the approximate values of the Gini coefficients in the cases here considered are reported in Table 1.  (h) [20-45 , 67 , 24], cases (i), (ii), (iii) Figure 3. Each of the panels (a-h) displays the Lorenz curves relative to the income distribution of the model quoted in the relative caption, in the three cases (i), (ii), (iii). It is evident from all panels that the Gini coefficient, and correspondingly economic inequality, diminishes when a basic income is introduced.

Conclusions
In this paper, a toy model is proposed which provides a stylised description of the introduction of a universal basic income in a closed society. The consequent effects of such a policy on the income distribution and the resulting level of economic inequality are investigated in the presence of different taxation schemes. The model is expressed in terms of differential equations, whose dynamics is explored by means of numerical simulations. These suggest that a basic income allocation acts as a tool of poverty reduction and leads to a decrease of economic inequality. As already anticipated in the introduction, we are well conscious that the issue of universal basic income in the real world is much more complex and cannot be faced without economic, social and political analysis. For example, the problem concerning the sustainability of the costs of such a program, an evaluation as to whether UBI could become a disincentive to work and commitment, to mention at least a couple of aspects, should be taken into account.
The purpose of this paper is simply to show that mathematical models can also help predict scenarios arising from one policy or another.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A
We give here some details on the numerical method employed to find the solutions of the equations (numbered as (1) in the text) x j , k = 1 , . . . , n .
We can write more concisely this differential equation system as with obvious definition of the vector field F = (F 1 , . . ., F n ) and with x = (x 1 , x 2 , . . ., x n ) .
If the two vectors v 1 = (1, 1, . . ., 1) ∈ R n , v 2 = (r 1 , r 2 , . . ., r n ) ∈ R n are introduced, the constraints (3) and (4) can be expressed as The Cauchy problem in correspondence to a fixed value of the total income (and for a population normalised to 1) amounts to search the solution x(t) of (A2) which satisfies for a given vector x 0 = (x 0,1 , x 0,2 , . . ., x 0,n ), where x 0 · v 1 = 1 and x 0 · v 2 = µ. The Equation (A1) (or, equivalently (A2)) together with the constrains (A3) and (A4) form a semi-explicit DAE (differential algebraic equation). By directly applying a standard numerical solver to the n Equation (A1) together with some initial condition x(0) = x 0 , one would not obtain the conservation of the population and of the total income; namely the constraints (A3) and (A4) would not be satisfied. In view of that, the numerical results described in the paper have been obtained by first performing a reduction of the number of system equations in the following way.
The constraints (A3) and (A4) allow the explicit expression of two of the components of a vector x satisfying them in terms of the others. For instance, given the vector x * = (x 1 , x 2 , . . ., x n−2 ) , one can obtain from (A3) and (A4) (and this was done here) x n−1 = ξ n−1 (x * ) , x n = ξ n (x * ) .