Mathematical models describing the effects of different tax evasion behaviors

Microscopic models describing a whole of economic interactions in a closed society are considered. The presence of a tax system combined with a redistribution process is taken into account, as well as the occurrence of tax evasion. In particular, the existence is postulated, in relation to the level of evasion, of different individual taxpayer behaviors. The effects of the mentioned different behaviors on shape and features of the emerging income distribution profile are investigated qualitatively and quantitatively. Numerical solutions show that the Gini inequality index of the total population increases when the evasion level is higher, but does not depend significantly on the evasion spread. For fixed spread, the relative difference between the average incomes of the worst evaders and honest taxpayers increases approximately as a quadratic function of the evasion level.


Introduction
The phenomenon of tax evasion represents a serious problem for several countries. Having as a main consequence a reduction of the tax revenue, it negatively affects a Faculty of Science and Technology, Free University of Bozen-Bolzano, Piazza Università 5, 39100 Bolzano, Italy correct functioning of the public sector: it hurts the supply of education, health care and services in general. In turn, this contributes, together with other factors, to an increase of economic inequality.
In this paper, we try and look at some aspects of the problem through a mathematicalmodelling approach. Specifically, we discuss a kinetic-type model for economic exchanges in a closed society in the presence of taxation and redistribution, within which occurrence of tax evasion to various extents is assumed: we admit the possibility that different citizens pay different percentages of the taxes they should pay.
We emphasize that the illegal practice under consideration involves in fact a large number of interacting agents. These include clearly the evaders themselves, but in addition also all other citizens who, by the situation, are deprived of the access to the benefits deriving from the revenue redistribution. In view of this and of the various levels at which it can take place, tax evasion can be thought of as an example of complex system. With this denomination, systems composed by a high number of heterogeneous units are meant, whose collective and macroscopic behavior is not derivable from the simple summation of the single units properties, but inherently depends on their nonlinear interplay. An interaction-based approach seems to be quite a natural one in the study of economic questions, but it was only during the last decades, especially thanks to the new opportunities offered by the increased computer power, that it started to be pursued. Arguments in favour and related work can be found e.g. in Aoki (2002), Arthur et al. (1997), Gallegati and Kirman (2012), Kirman (2006Kirman ( , 2010, Landini et al. (2015) and Tesfatsion and Judd (2006). The technical tools more frequently employed in connection with this approach are agent-based computational algorithms and simulations, possibly combined with complex networks theory. Also, starting in the mid-1990s an interdisciplinary research field called econophysics has been growing, which explores the dynamical behaviour of economic and financial markets by means of methods taken from statistical mechanics and gas kinetic theory, see e.g. Chakrabarti (2009), Chatterjee et al. (2005), Heinsalu and Patriarca (2014) and Yakovenko and Barkley Rosser (2009).
Adopting here a mathematics-supported complex system perspective, we aim at deriving and explaining the emergence of a population's aggregate feature like the income distribution, as a result of the whole of economic exchanges and interactions which take place between the individuals. More specifically, our focus here is on the effects that the heterogeneity of taxpayer behaviors has on the income distribution profiles of the categories of individuals evading to different degrees. To this end we build on our previous work, incorporating suitable additional facets into the models discussed in Bertotti (2010) and Bertotti and Modanese (2011, 2014a. We also emphasize that our goal here is mainly a methodological one: more than perfectly representing real world features, we aim at constructing a tool endowed with explorative ability. Our model in fact has this character. Indeed, the possibility of finding numerical solutions which evolve from differently chosen "initial conditions" and in the presence of differently tuned parameters amounts to the ability to forecast the emergence of different scenarios. In turn, this can give insights as to which policies could be adopted to favour or prevent desired and undesired trends. The tax evasion process is the object of a large amount of literature. Works among those which involve agent based modeling and simulations include e.g. Bloomquist (2006), Crokidakis (2014), Hokamp and Seibold (2014), and Zaklan et al. (2009), to name but a few. Specific aspects therein investigated concern the occurrence of behavioral changes of agent types with diverse moral attitude, due to imitation or also to some audit procedure. In particular, in Hokamp and Seibold (2014) and Zaklan et al. (2009) a variant of the Ising model originally developed within the theory of magnetism is considered. In that context, each spin represents a citizen, which can be either in the tax compliant state +1 or in the tax evader state −1. Citizens undergo transitions from +1 to −1 caused by imitation of their nearest neighbours, and from −1 to +1 induced by tax audits. Indeed, the consequence of an audit on an evader is assumed to be the fact that she/he will remain honest for a certain number of steps. This approach is helpful for the analysis of evasion phenomena in relation to local interaction and external controls, but not for studying the effect of evasion on the income distribution as we do here.
The paper is organised as follows. In Sect. 2 we describe the model and some of its analytical properties. In Sect. 3 we report and discuss results from a set of numerical solutions. Section 4 contains our conclusions.

The model
This section is devoted to a short description of the model proposed. More details on the primary mechanism underlying it can be found in our papers (Bertotti 2010;Bertotti and Modanese 2011, 2014a. We point out however that in Bertotti (2010), Bertotti andModanese (2011, 2012) the phenomenon of tax evasion was not taken into account, whereas in Bertotti and Modanese (2014a, b) all individuals were assumed to engage to the same degree in a kind of evasion different from that dealt with here. 1 The main novelty here is given by the assumption of the existence of different degrees of evasion.
The totality of individuals (supposed to remain constant in time, in fact during the time period under consideration) is divided into a number n of classes, each one characterized by its average income, the average incomes being the positive numbers r 1 < r 2 < · · · < r n , and, in turn, every income class is divided into a number m of sectors characterized by possible evasion behaviors.
We denote by x α j the fraction of individuals belonging to the j-th income class and to the α-th evasion behavior sector. We will say for brevity that x α j is the fraction of individuals of type ( j, α), also called ( j, α)-individuals, the number of different groups being n × m.
We suppose here that the evasion behavior of each individual remains constant in time. In contrast, individuals may move through different income classes. The model is formulated by a system of n × m ordinary differential equations describing the variation in time of x α j for j = 1, . . . , n and α = 1, . . . , m. Such a variation is the result of direct economic interactions, in which pairs of individuals exchange some money. And is affected as well by the payment of taxes (in some cases, the due taxes and in other cases, partial quotes of them) and by the revenue redistribution, represented in the real world by healthcare, education and services in general. More specifically, the dynamic process is as follows: a whole of interactions between pairs of individuals occur simultaneously: for any h and k in {1, . . . , n}, for any β and γ in {1, . . . , m} individuals belonging to the h-th income class and the β-th evasion sector meet individuals of the k-th income class and the γ -th evasion sector and some money exchange between such pairs takes place. How many of these economic exchanges do we have? If at the considered time the fraction of (h, β)-individuals is x β h and the fraction of (k, γ )-individuals is x γ k , the number of encounters of these two categories of individuals is the product x β h x γ k . And any single encounter contributes, albeit to a very small extent, to a change of the fraction of individuals in some income classes. The differential equations contain several parameters, which express for example transition probabilities, the probability that in an encounter between two individuals of different classes the one or the other is paying, the tax rates relative to the different income classes and the percentages of evasion.
Similarly as in Bertotti and Modanese (2014a, b) we first define for h, k = 1, . . . , n certain coefficients p h,k , aimed to specify the probability that in an encounter between an individual of the h-th income class and one of the k-th class, the one who pays is the former. We take with the exception of the terms p j, j = r j /2r n for j = 2, . . . , n − 1, p h,1 = r 1 /2r n for h = 2, . . . , n, p n,k = r k /2r n for k = 1, . . . , n − 1, p 1,k = 0 for k = 1, . . . , n and p h,n = 0 for h = 1, . . . , n. The introduction of these coefficients implies that, with only reference here to the income class (and independently of the evasion sector), for each h − k pair there is -a probability denoted by p h,k ∈ [0, 1] that the h-individual will transfer some money to the k-th one, -a probability p k,h ∈ [0, 1] that the k-individual will transfer some money to the h-th one, -a probability 1 − p h,k − p k,h ∈ [0, 1] that the two do not exchange money. Correspondingly, we are assuming that for any h ∈ {1, . . . , n} and any k ∈ {1, . . . , n} the frequency of payment of individuals of the h-th income class to individuals of the k-th income class is a fixed one. We could call this a compartmental representative agent behavior [see in this connection (Bischi 1998;Tramontana and Gallegati 2010)].
Then, we introduce S (with S << (r i+1 − r i ) for all i = 1, . . . , n), the amount of money that in each direct transaction one individual is supposed to pay to another. The individual who receives the money is expected to pay a part of this as a tax to the government. If this individual belongs to the k-th income class, we may for sake of simplicity assume that he should pay an amount corresponding to S τ k , τ k being the tax rate of his income class.
At this point we notice that in a tax compliance case the effect of a direct interaction with an individual of the h-th income class paying S to one of the k-th class and this paying the due tax would be equal to that of the first individual paying an amount S (1 − τ k ) to the k-th income class one and paying as well a quantity S τ k to the government or equivalently, due to the redistribution, to the community of individuals. 2 Being interested in a tax evasion case study, we now introduce, beside the tax rates τ k , some other parameters. For any α = 1, . . . , m let denote the percentage of the due taxes payed by individuals characterized by an evasion behavior index α. Then, we define for any k = 1, . . . , n and α = 1, . . . , m, The quantity θ k,α in (2) expresses the fraction a (k, α)-individual actually pays as a tax; this percentage depends both on the income class represented by the index k and on the evasion index α.
Example 1 As an example to illustrate the situation, take m = 3 and consider three evasion behaviors, described by In such a case one would have -individuals in the first sector paying all they should and not evading at all, -individuals in the second sector paying half of what they should, -individuals in the third sector paying one quarter of what they should.
We also emphasize that the coefficients p h,k enter in the formulae (4)-(6) in such a way that their effect can be also interpreted as weighting the amount of money exchanged. In other words, the situation is the same one would have assuming the frequency of payment of individuals independent on the income class, but with the amount of money paid in each transaction by individuals of the h-th income class to individuals of the k-th income class equal to p h,k S instead of S. The specific choice of p h,k adopted here is suggested by the phenomenological observation that typically poor people pay and earn less than rich people.
Remark 1 It may be helpful stressing here the fact that, even if the equations (1) describe migrations of aggregate fractions of individuals, in fact a probabilistic microinteraction modelling underlies the dynamical process expressed by these equations.
The right hand sides of (1) contain quadratic (and other nonlinear) terms, exactly because they give account of a large number of pairwise interactions (and, through the taxation and redistribution process, also of interactions involving more individuals). For example, the origin of the coefficients C ( j,α) (h,β);(k,γ ) (which refer to direct money exchange) is the following. The interaction between an (h, α)-individual and a (k, β)individual with the (h, α)-individual paying, produces the variation of the fraction of individuals in some groups (in general, four). Indeed, the (h, α)-individual becomes a little bit poorer, inducing a partial migration from the h-th income class to the (h−1)-th one and the (k, β)-individual) becomes a little bit richer, inducing a partial migration from the k-th income class to the (k + 1)-th one. The mentioned variation in the four groups is described through the coefficients: The coefficients in the formulae (4)-(6), namely those appearing in the equation (1) and which refer to the ( j, α) group, are then obtained from these, observing that each probability C ( j,α) (h,β);(k,γ ) can be written as a sum where the "absence-of-variation term" a ( j,α) (h,β);(k,γ ) = 1 only if j = h and α = β, independently of (k, γ ), and the b ( j,α) (h,β); (k,γ ) are as in (8) (see Bertotti (2010) for a similar discussion in a simpler case).
We also observe that the structure of the C ( j,α) (h,β);(k,γ ) and the T ( j,α) [(h,β);(k,γ )] (x) in (4)-(7) is determined by the conservation requirements of the global mechanism. The stochastic character they enjoy is due to the presence of the coefficients p h,k which can be defined with some degrees of freedom. It is in view of the coefficients p h,k that we call this a probabilistic micro-interaction modelling. We do not know who exactly is going to interact with whom: we only know at a probabilistic level how often individuals in a group interact with individuals of another group.
Of course, assuming that each individual of the h-th income class has the same probability p h,k of paying in an encounter with an individual of the k-th class corresponds to attribute the same behavior (intended as attitude to pay) to all pairs of individuals of the same two specific classes. This reminds a mean-field approach, see e.g. Aoki (2002) and Aoki and Yoshikawa (2007). But, in fact, the underlying modelling just described provides our approach with a different characterization. 3 By suitably adapting proofs, which are relative to a previous, less general, version of the model and can be found in Bertotti (2010), one may check that the following properties, amounting to conservation in time of both the number of individuals and the global income as well, hold true. x α j (t) = 1 for all t ≥ 0.
As a consequence of this property, the expressions of the T ( j,α) [(h,β);(k,γ )] (x)'s in (7) somehow simplify and the right hand sides of system (1) turn out to be in fact polynomials of degree 3.

Property 2
The scalar function μ(x) = n j=1 r j m α=1 x α j remains constant along each solution of system (1).
In addition, the running of several numerical solutions provides evidence of the following fact.
Property 3 If the parameters of the model (r 1 , . . . , r n , S, the τ k 's, the θ ev (α)'s) and also the fraction of individuals with different behavior 4 are fixed, if μ ∈ [r 1 , r n ] is fixed, 3 In the Boltzmann approach to statistical mechanics, to which our approach is inspired, the variables are described by means of a probability distribution function. The Maxwell-Boltzmann statistics gives the expected number of particles in a given volume of the phase space. Space and velocity are continuous variables. The discretized Boltzmann approach is somewhat more manageable, because it groups together particles with the same velocity or, in the socio-economic version, individuals in the same income class. Integrals are then replaced by sums and the number of admissible velocities or classes can be increased as needed to ensure the precision required in any specific case. In comparison to the discretized Boltzmann, the mean-field approach entails a loss of information, since the evolution equation for each individual implies that it interacts with an average value of all the others, and the whole evolution process is self-consistent. In the discretized Boltzmann approach and in our model the interactions are microscopic and occur between all individuals. Finally, if one compares a Boltzmann approach with an agent-based simulation, it is fair to say that the Boltzmann agents have properties typical of deterministic particles, while the behavior of simulated agents can be much more flexible. However, the Boltzmann approach offers the advantage of a closed mathematical formulation, independent from the software. 4 The fraction of individuals with a specific evasion behavior is assumed here to be the same in each income class. x α 0 j = μ tend asymptotically to a same stationary distribution as t → +∞.

The evidence from numerical solutions
Our specific interest here is to analyze the income evolution in the long time limit of groups of individuals characterized by different behaviors. In other words, whereas our previous investigation in Bertotti and Modanese (2014a, b) was especially addressed to detect the effects of (a different kind of) tax evasion on the population as a whole, here we also focus on the evasion effects on the different mentioned groups. Finding analytical solutions of the nonlinear differential equations (1) is of course hopeless. However, numerical solutions provide sufficient information. In order to obtain them, one has to first fix several parameters. We take here n = 9, m = 3, S = 1, r j = 10 j for j = 1, . . . , n and assume tax rates increasing from a minimal one, τ min , to a maximal one, τ max , according to a progressive taxation system, as given by Still, the values of τ min , τ max and θ ev (α) for α = 1, . . . , m remain to be chosen. Each time we explore aggregate effects in the asymptotic stationary income distribution of the population and compare situations of tax compliance and tax evasion, we find that evasion leads to an increment in the number of individuals in the poorest and in the richest classes, at the detriment of the middle classes. Figure 1 illustrates the typical situation. In Bertotti and Modanese (2014a, b) we have shown that tax evasion has in general the effect to increase economic inequality, as measured for instance by the Gini index. We have also studied situations in which evasion grows in response to an increase of the tax rates, finding an optimal "compromise" characterized by minimum inequality. These results were obtained in conditions of homogeneous evasion rates.
Assume now, to fix ideas, that the evasion behaviors are as in (3) in Example 1 and that each of them is present in one third of the population. 5 Accordingly, in each income class one third of the individuals pays all due taxes, one third pays half of them and one third pays one quarter of them. Not only can we obtain information on the aggregate shape of the asymptotic income distribution. We also get a more detailed picture of the effects of tax evasion within different behavioural sectors. Specifically, Fig. 1 Collective effect of tax evasion. The asymptotic income distribution on the left refers to a case with tax evasion, the one in the middle refers to a tax compliance case with the same initial conditions and the panel on the right displays the difference of the fraction of individuals in the various classes in the first and second case. Note that the figures are scaled differently Fig. 2 In the two rows the asymptotic income distribution referring to two different initial data are plotted. The two panels on the left display the distributions by income classes. The histograms in the central panels represent the sequence for each income class of three evasion behavior sectors with θ ev (1) = 1, 1/2, 1/4. The panels on the right provide an alternative bi-dimensional representation, in which the sectors can also be distinguished it turns out that in low income classes the numbers of "honest individuals", "halfevaders" and "three-quarters evaders", counted in this order, are decreasing from the largest to the smallest, while the situation is reversed in the high income classes. A graphic illustration of this is provided by the Fig. 2.
In order to understand which variations of the model parameters can best represent some real situations, and what results we can expect, we can compare situations in which evasion is widespread with situations in which it is confined only to a part of the population, the total evasion level being the same. By total evasion level we denote the total fraction of tax payments which is eluded, summed over all sectors of the population. With three sectors, for instance, we can consider the two following situations, which have both total evasion level 1/6, but exhibit respectively widespread and concentrated evasion: 1. one sector is honest, with θ = θ ev = 1, i.e. 100% of taxes are paid, while the other two sectors are slightly dishonest, with θ = 0.75, i.e. 75% of taxes are paid; 2. two sectors are honest and the third is quite dishonest, with θ = 0.5.
We wonder whether the Gini index is significantly different in the two cases, or in other words, if evasion has a different impact on inequality when it is widespread or confined to a part of the population. Furthermore, we might wonder whether the partial Gini indices of the behavioral sectors are substantially different from that of the A gradual spread of evasion in three behavioral sectors is assumed. For instance, when the total evasion level is 10% (the first sector pays 100% of the due taxes, the second pays 90% and the third 80%), the average income of the third sector is 6.8% larger than the income of the first sector. The applied tax rates grow linearly with income between τ min = 10% and τ max = 45%. The dependence d(η) can be approx. fitted as d(η) 0.42η 2 + 0.62η total population, for instance in the sense that among evaders there is more inequality than among honest taxpayers. The numerical solutions show, however, that in all these cases the Gini index does not exhibit any significant variations, depending only on the total evasion level. The introduction of behavioral sectors into the model allows to observe another important effect, namely the appearance of a clear difference between the average incomes of honest taxpayers and tax evaders. It is already apparent from the income histograms in Fig. 2 that, as expected, tax evaders tend to get richer than honest taxpayers: the bars representing the lowest-income classes clearly show a larger share of honest people, while the opposite happens for the highest-income classes, where more evaders are found. It is interesting to give a quantitative estimate of this difference and to study its dependence on the evasion level.
To this end, let us consider the relative difference d between the average income of the worst evaders and the average income of the honest taxpayers: d = (μ m −μ 1 )/μ 1 . This relative difference is quite large, typically of the order of 10 to 20% in the histograms of Figs. 1 and 2. In order to evaluate its dependence on the evasion level, we need to choose a fixed and reasonable "evasion spread" pattern.
The dependence is manifestly non-linear, showing that the phenomenon is complex and its interpretation not simple. In fact, an increase in the evasion level affects the income distribution in at least two ways: (a) through direct interactions, because evaders gain systematically more from any interaction; (b) through indirect interactions (tax redistribution), because when tax evaders tend to outnumber honest taxpayers in the higher-income classes, which should pay higher tax rates, the total tax collection is diminished. It is not obvious, however, that this diminution should further increase the difference d(η) defined above, since in the present version of the model redistribution is uniform. In the version with heterogeneous redistribution Bertotti and Modanese (2015) we could take into account more subtle real effects: for instance, supposing that welfare provisions are means-tested, we could assign them based on the tax paid, instead than on the real income. This is known to give a further unjust and detestable advantage to tax evaders. We plan to address these issues in future work.

Concluding remarks and further perspectives
In this paper, a kinetic-type model describing economic interactions, taxation and redistribution in a closed society is discussed. The focus is on the effects produced by tax evasion by individuals who under-declare their income in different measure.
The model suggests the following considerations. From the sound point of view of individuals who care for society, tax compliance plays an important role towards the overcoming of economic inequality. From the point of view of selfish individuals, the probability of improving their own economic status is higher when evading. The result is not surprising. Certainly, it has to be noticed that no audit actions or punishment are taken here into consideration. Incorporating them in the model and investigating possible impacts would be an interesting point to explore in future work.
Another, related, point which deserves further work is the necessity of inclusion into the model of possible changes of the behavioral taxpayer attitudes. Of course, in real life, the propensity of individuals to be compliant or non compliant does not remain always constant in time: in particular, it can be influenced by the behavior of others and by specific experiences (such as fines) with fiscal agencies. Combining into the model the treatment of the money distributional aspect with behavioral and psychological factors remains a major challenge to be faced in the future. As well, the possibility of bankruptcy and production should be taken into account in an effort of an improved realism.
In this connection and in conclusion, we like to emphasize that our goal here has been mainly a methodological one. The interest of this model (possibly bound to be further enhanced) lies in our opinion in its possible explorative use: simulations corresponding to different conceivable parameters allow to understand and forecast the emergence of different scenarios and could possibly suggest policies addressed to favour desired trends or prevent undesired ones.