Bounded Rationality and Animal Spirits: A Fluctuation-Response Approach to Slutsky Matrices

The Slutsky equation, central in consumer choice theory, is derived from the usual hypotheses underlying most standard models in Economics, such as full rationality, homogeneity, and absence of interactions. We present a statistical physics framework that allows us to relax such assumptions. We first derive a general fluctuation-response formula that relates the Slutsky matrix to spontaneous fluctuations of consumption rather than to response to changing prices and budget. We then show that, within our hypotheses, the symmetry of the Slutsky matrix remains valid even when agents are only boundedly rational but non-interacting. We then propose a model where agents are influenced by the choice of others, leading to a phase transition beyond which consumption is dominated by herding (or `"fashion") effects. In this case, the individual Slutsky matrix is no longer symmetric, even for fully rational agents. The vicinity of the transition features a peak in asymmetry.


INTRODUCTION
Economic theory still appears to be wedded to the strict Homo economicus paradigm, despite numerous well established facts that contradict its very essence. In most models, economic agents are assumed to be fully rational, isolated, and homogeneous optimizers, capable of maximizing utility efficiently and instantaneously, at odds with many behavioral experiments and common sense. Even though boundedly rational [1][2][3] or event irrational agents could behave collectively as a unique representative rational agent, this is very much a leap of faith, specially in the presence of strong interactions that can generate trends, fads and fashions, herding or mass panics, to name a few -see e.g. [4][5][6][7]. Empirical data also reveal very strong, Pareto heterogeneities in wealth, firm sizes, market capitalisation, etc. [8][9][10][11] that may considerably affect aggregate quantities: this is Gabaix's "granularity hypothesis" [12]. 1 Other effects, such as habit formation and hysteresis can also play a major role in the behaviour of consumers, firms and investors alike [14][15][16][17].
The Slutsky matrix [18], central in consumer choice theory, stands as a textbook example of a mathematically sound economic construction that most likely lacks realism. It describes the response of consumption of good i to a change of price of good j. From a theoretical standpoint, the Slutsky equation has been extensively investigated, see e.g. [19][20][21], focusing notably on its geometrical interpretation [22]. Under the classical assumptions mentioned above, very clear conclusions can be drawn: the Slutsky matrix is symmetric, negative semi-definite, * jerome.garnier-brun@polytechnique.edu 1 Heterogenous agents are progressively (but only partly) accounted for in recent macroeconomic models, see e.g. [13]. and satisfies "homogeneity" -meaning the price vector of the products considered is an eigenvector of the matrix associated to a zero eigenvalue.
Comparatively very little has been done from the empirical perspective. Looking at private consumption in the pre-and post-WWII periods in the US [23] and the Netherlands [24], statistical tests have ruled in favor of negative semi-definiteness but remained rather inconclusive regarding symmetry and homogeneity. Other studies have supported the symmetry condition [25,26], or found that it is satisfied for single households and violated for couples, but in a way that can be accounted for [27]. Despite these conflicting conclusions, the Slutsky matrix and its associated conditions still stand at the core of consumer choice theory, or as nicely put by Barten [24]: "Such a situation is not uncommon in economics [...]. Many statements are generally accepted as being a true description of behaviour without a thorough questioning of their empirical validity." For further criticism of Slutsky theory and a particular focus on how income distribution changes are not taken into account, see [28]. In response to the inconclusiveness of empirical tests, some authors have shown that there is in fact no reason for symmetry to hold when taking into account the fact that consumers are many [29], non-cooperative households [30], or bounded rationality in the form of limited attention [31]. More recently, there have also been some attempts at taking into account the somewhat unrealistic nature of the assumptions behind the classical Slutsky equation, with the idea of providing a measure of some loosely defined rationality based on the nearness between the theoretical prediction and estimates constructed from observed demand functions [32,33].
In this paper, we introduce a rigorous mathematical framework inspired by statistical physics that is able to accommodate bounded rationality, population heterogeneity, and interactions between agents, with in particu-lar herding tendencies. We investigate how the inclusion of such behavioral features impacts the classical properties of the Slutsky matrix. One of our central results is a novel representation of the Slutsky matrix in terms of the consumption fluctuations for fixed prices and budgets, instead of the traditional response to a change of prices and budgets. Using this language, the symmetry property of the Slutsky matrix for a unique, boundedly rational agent becomes quite transparent. However, we show that including interactions between agents can and does affect this property -although all eigenvalues still retain a non-positive real part.
The model we propose to account for interactions between consumer choices turns out to be very interesting in its own right. We show that beyond a certain interaction strength, the choice of consumers condenses onto a limited set of products, that become desirable just because others find it desirable, in a self-fulfilling fashion. We determine the boundaries of this condensed regime in parameter space. It turns out that the Slutsky matrix is maximally non-symmetric close to such boundaries.
The paper is organized as follows. In Section I, we briefly introduce classical consumer choice theory and the Slutsky equation. Its generalization to bounded rationality and interacting agents is given in Section II, where the fluctuation-response formulation of the problem is also derived. Two possible definitions of the aggregate Slutsky matrix are also discussed. In Section III, we propose an interacting model and study the condensation of consumer choices. The consequences of this collective effect on the Slutsky matrix and its properties is then presented in Section IV. The equivalence of the thermodynamic ensembles used and possible out-of-equilibrium effects are discussed in Section V. We finally summarize our findings and discuss possible future work in Section VI.

I. CONSUMER CHOICE THEORY
Consumer choice theory is based on the idea that, for a given bundle of M goods with prices p ∈ R M + , agents choose a basket x to maximize their utility function u(x), while subject to the constraint p · x = w, where w is the consumption budget.
The utility function must then satisfy some elementary conditions. In particular, it is taken to be an increasing function of the quantity of goods, i.e. more is always better and there is no satiation. While the utility must therefore have a positive first derivative, its rate of increase should be a decreasing function of the absolute quantity of goods. In other words, utility is postulated to be concave, i.e. the marginal utility gain diminishes as goods are accumulated.
Assuming that the agent is fully rational, the optimal basket that maximizes the agent's satisfaction x * is then simply given by Conventionally one defines the Marshallian demand x(p, w) which corresponds to that obtained from solving problem (1), and the Hicksian demand h(p, u) defined as the demand that minimizes the expenditure e(p, u) for a fixed utility level u [34]. Setting h(p, u) = x * (p, e(p, u)) and differentiating, one obtains the Slutsky equation describing the change in consumption of good i following a change in the price of good j Changes δx in the optimal basket's composition in response to a price change δp can thus be separated into two contributions: the substitution effect (first term in the RHS of Eq. (2)), describing how consumption is impacted by changes in relative prices of goods, and the income effect (second term in the RHS of Eq. (2)), expressing the impact of changes in purchasing power. The substitution effect is often described in terms of the Slutsky matrix S, with entries therefore defined as Provided that the utility function is sufficiently regular, S can then be shown to be symmetric, negative semidefinite, and equal to the Hessian of the expenditure function. In practice, the Slutsky matrix cannot be observed directly, but can be estimated as the other two terms in Eq. (3) should be accessible empirically. In the following, we will give an alternative theoretical expression for the Slutsky matrix in terms of consumption fluctuations.

II. BOUNDED RATIONALITY
As described in the introduction, taking agents to be perfect optimizers seems to be an unrealistic assumption in most contexts. There are several ways to relax such an assumption. One is that agents have a limited attention and cannot process all the information accessible to them, see e.g. [31] and refs. therein. As argued by Gabaix, this may effectively lead to perceived prices that differ from real prices. This in turn affects the symmetry of the Slutsky matrix.
Another traditional line of thought in the literature on choice theory is to replace the deterministic utility optimization prescription by a stochastic choice rule, such that the stationary choice distribution is given by a Boltzmann-Gibbs measure (see below) with an intensity of choice parameter that allows one to interpolate between full rationality and completely random choicesfor a classical review, see [35]; see also [6] for more recent developments. This is the prescription we will adopt in the rest of the paper; as we shall see this allows us to describe a rather wide range of phenomena while ensuring mathematical tractability of the model.

A. A Single Agent
Formally, considering first a single agent, we postulate that the probability density for selecting the basket of goods x is given by where u(x) is the utility function and the parameter β is known as the intensity of choice (or the inverse temperature in physics parlance). Taking β → ∞, one recovers the fully rational case previously introduced and equation (1), while β = 0 means the agent picks any basket satisfying the constraint with equal probability. The normalization factor Z is known in statistical physics as the partition function, which, given the hard budget constraint enforced here, writes where + means that we integrate over non-negative baskets x ≥ 0. For finite values of β, the basket x must now be statistically described, since different realizations of the system will lead to different outcomes. A suited definition of the Slutsky matrix must therefore be considered.
Here, we propose to replace all x i by their averages x i , with angular brackets referring to an average over the distribution given in Eq. (4), to wit Note that in the limit β → ∞, one recovers the standard Slutsky matrix since x i → x * i . Our set up of the problem allows us to draw several analogies with the statistical mechanics foundation of thermodynamics. For example, the strict application of the budget constraint δ(p·x−w) is reminiscent of the socalled canonical ensemble, where the conservation of the number of particles in the system is strictly enforced. We will see later that one can also work in the analogue of the grand-canonical ensemble where the budget constraint is only enforced on average. This eases some analytical calculations while being equivalent to the canonical ensemble in some limits (for example when the number of goods is large). One can argue that in some cases, allowing the budget to fluctuate (due to loans for example) can be realistic as well.

B. A Fluctuation-Response Relation
More interestingly, statistical mechanics also provides relations between the response of certain quantities to external perturbations to spontaneous fluctuations of these quantities in the absence of perturbations. These relations can be established using the derivatives of the partition function, assuming an equilibrium state has indeed been reached. In particular, the Slutsky matrix can indeed be expressed in terms of other correlations, as was first mentioned in one of the present authors (JPB) contribution to [36]. In Appendix A 1, we derive the following "fluctuation-response" formula in the single agent case: with Γ = ∂ w log Z and x i x j c := x i x j − x i x j . Equation (7) is manifestly symmetric in i,j. It shows that even with bounded rationality, the Slutsky matrix is still symmetric, for any value of β, not only in the rational limit β → ∞. Hence, symmetry of the Slutsky matrix may not be used as a proof of the rationality of economic agents, as was argued in [36] for example. Our fluctuation formula Eq. (7) is also interesting from an econometric standpoint, as it provides a way to measure the Slutsky matrix without varying prices. Measuring response quantities from equilibrium correlations is in fact commonly used in statistical mechanics, through what is referred to as fluctuation-dissipation relations [37], or, in a restricted context, to Einstein's relation relating mobility to diffusion for Brownian particles.
One can go one step further and eliminate all derivatives from Eq. (7), provided the utility function is known. One finds (see Appendix A 3) and both equations being valid for an arbitrary choice of k.

C. Many Agents
In practice, agents make correlated choices and we must adapt our formalism to treat interactions. We thus consider N agents, indexed by α = 1, . . . , N , and denote U ({x α }) the aggregate utility of all the agents. Interactions mean that U cannot in general be written as a sum of individual utility functions u α . The probability that the N agents choose a bundle of goods {x α } is proportional to exp(βU )/Z N , where the aggregate partition function Z N writes Here, we integrate over the M × N degrees of freedom, i.e. the quantities x α i of good i consumed by agent α, while enforcing that all agents respect their own specific budget w α . 2 Analytically computing this partition function is usually very difficult due to the product of Dirac δ distributions in the integrand. Depending on the form of the utility, one might need to slightly relax the budget constraints. As mentioned above, a way to do this is to move to the so-called grand-canonical partition function Z N , defined as where the µ α , known in physics as "chemical potentials", are fixed by enforcing that budgets are satisfied on average, i.e. p · x α = w α , ∀α. In general, the two partition functions are not equivalent; however, many quantities calculated from the two ensembles become identical in the large M limit and/or in the rational limit β → ∞, see Section V B for a detailed discussion.
Consistent with the single agent definition, we take the individual Slutsky matrix of agent α to be given by where we have assumed that prices are the same for all agents. Note that we have also taken a uniform systemwide rationality parameter β, although a generalization to different β α 's is possible and would be an interesting extension of our work. As in the single agent case, the partition function allows one to derive a fluctuation-response expression for the Slutsky matrix in terms of correlations. Equation (12) may be rewritten (see Appendix A 2) in what will be referred to as its "thermodynamic" form with Γ γ = ∂ w γ log Z N . This is the central theoretical result of this paper, which, to the best of our knowledge, is new. Although its derivation makes explicit use of the Boltzmann-Gibbs distribution, such an expression holds more generally in the near-rational limit β → ∞ where small deviations around the rational solution are always Gaussian, see Section II E below. Importantly, unlike in the single-agent case, Eq. (13) does not allow one to infer anything general about the 2 It should be noted that with this partition function we implicitly assume that all agents work at improving the aggregate utility of the system, which may not be the case in reality. For a complete discussion addressing the possible differences when considering purely individualistic agents, see Section V A.
symmetry of the Slutsky matrices. In the case where interactions between agents are negligible, i.e. when U is the sum of individual utility functions, correlations between agents are zero whenever γ = α and we recover the single agent expression, as expected. Assuming the utility function is known, one can again go further and express derivatives with respect to budgets w γ as a function of some correlations. For our problem we find a generalization of the formula obtained for a single agent in the canonical ensemble (see Appendix A 3) as well as both again valid for any good k. From these expressions and given a utility function U , all terms in Eq. (13) can be computed in principle, at least numerically. Bear in mind, however, that these relations are only valid if the system has reached equilibrium, which might take a very long time, for example near a phase transition point.

D. Aggregate Slutsky Matrices
There are a priori two possible definitions of the Slutsky matrix at the aggregate level. One is simply to take the average over all agents of the individual Slutsky matrices, to wit However, empirical measurements often rely on estimates at the aggregate level. In that case, a better suited definition uses aggregate consumption: with overlines indicating averaging over agents e.g.
and similarly for the average budget w. In that case, as shown in Appendix A 4, the thermodynamic expression becomes with the possibly heterogeneous factor Clearly, if consumption is proportional to wages and if all wages scale with the average wage, i.e. w γ = κ γ w, then x γ j = κ γ x j , and there is no contribution from the last term. In this case, we find an expression very close to the single-agent case. More generally, S ij has no reason to be symmetric, except when all agents are identical. In such a case, even in the presence of interactions and for bounded rationality (β < +∞), S ij is always symmetric, whereas S ij is not, as we in the next sections.

E. Near-Rational Limit β → ∞
In order to simplify the problem and get some intuition, we place ourselves in the near-rational case where β → ∞. This corresponds to the low temperature case in a physical system, where one can expect that all relevant configurations {x α } are close to the optimal configuration {x α * } that maximizes the global utility function U subject to budget constraints. We thus write δx α i := x α i −x α * i and Taylor-expand the utility function to second order, resulting in with H the (M × N ) × (M × N ) Hessian of the system evaluated at the maximum of the utility function, for a given set of budget constraints w α . Here we only consider deviations {δx α } that all satisfy the budget constraints, so that the first derivative terms U are zero as we expand around a maximum along all directions but one. Calculating the partition function and correlations then simply amounts to computing Gaussian integrals, although the budget constraints must still be handled with care. 3 By taking the Fourier representation of the Dirac δ (see Appendix A 5), we finally find, to leading order in β −1 : with the N × N matrix G defined as and where the second term in the RHS of Eq. (20) is the result of the constraint being applied.
Let us first illustrate this formula in the N = 1 case. In the canonical ensemble and large β regime, one has where λ is the Lagrange parameter enforcing that the budget constraint is satisfied. Hence the first term Eq. (7) remains finite since Γ diverges as β when x i x j c tends to zero as β −1 . The second term, on the other hand, tends to zero at least as β −1 . It then follows that for a single and near-rational agent, This is the classic expression for the Slutsky matrix, which our framework therefore allows to recover in the corresponding limit, with corrections in β −1 that can be computed. Since the Hessian is both symmetric and negative semi-definite at a utility maximum, we recover the classic properties of the Slutsky matrix. "Homogeneity" is also easily recovered by checking that multiplying the matrix by p indeed gives a zero eigenvalue. When N > 1, we will need to specify the utility function U to make the final result more explicit. Keeping U fully general and taking the limit β → ∞ only allows one to simplify the general expression Eq. (13) to where we have used Eq. (A20) and expanded ∂ x U to first order in δx. A similar expression can be derived for the aggregate Slutsky matrix S ij as well. Remember that C is of order β −1 (see Eq. (20)) so S α ij is well behaved in the large β limit. Although the final expression is not transparent, it is clear that S α ij has no reason to be always symmetric, except when agents' choices are uncorrelated, in which case C αγ ij = 0 whenever α = γ and C αα ij is symmetric by construction. We will now turn to an explicit model with herding, where the asymmetry of Slutsky matrix can be made apparent.

III. ANIMAL SPIRITS
As argued in the introduction, both anecdotal evidence about fads and fashion and more serious scientific studies point to the fact that agents' choices can be strongly influenced by the choice of others (see e.g. [4,6,[38][39][40][41]), an effect also known as "keeping up with the Joneses" [15,42], see also [7]. We now propose a family of models which account for interactions between boundedly rational agents. Interestingly, these models lead, in some regions of parameters, to "concentration" (or "condensation") of choices, much in the spirit of the model proposed by Borghesi and Bouchaud [43] (see also [44] 0.16 for a recent extension). As we shall see, close to the concentration transition, the non-symmetric contribution to the Slutsky matrix reaches a maximum.

A. Interactions and Herding
In order to study imitation or fashion effects, we first consider agents with log-utilities and take the aggregate utility to simply be the sum of all agent-specific utilities, where a α i describes the preference of agent α for good i. In the following, we assume that agents are homogeneous (a α i = a i , ∀α) but interacting, by which we mean that the preference for good i depends on how much good i is consumed by other agents. Mathematically, we posit that the preference for goods increases with the k-th power of their average consumption: where we remind x i is the average consumption of good i (over all agents), and c and k are non-negative parameters that describe the strength and nature of the interactions. 4 The non-interacting case c = 0 (or equivalently k = 0, up to some rescaling of the a i ) may be treated exactly in the canonical ensemble, i.e. strictly enforcing the budget constraint. As detailed in Appendix C 1, the equilibrium configurations are given by which will henceforth be referred to as the non-condensed or uniform solution. This solution matches results from numerical experiments, as illustrated in the two top panels of Figure 1(a) in the fully rational case, for example. For details about how numerical simulations have been performed, see Appendix B. The agent specific Slutsky matrices S α ij can then also be written explicitly from the original definition, verifying both symmetry and negative semi-definiteness for any β. In this non-interacting case, budget heterogeneities simply affect the magnitude of any given agent Slutsky matrix entries and are thus inconsequential for the properties of interest.
Taking k > 0 and increasing c, we expect the system to progressively depart from this solution and concentrate on some product(s) as c → ∞, as observed in the numerical simulations shown in Fig. 1(a), bottom panel. To evaluate how this concentration occurs, we start by taking w α = w, ∀α, that is a system of interacting but identical agents, both in wealth and preferences. As is often the case in statistical mechanics, the partition function cannot be computed exactly for interacting systems for general N , but can be more and more accurately approximated in the large N limit.
In our case, it is convenient to relax the budget constraint and to place ourselves in the grand-canonical ensemble where the budget constraint is only enforced on average. The procedure, detailed in Appendix C 2, allows us to rewrite the grand-canonical partition function as an integral over the mean consumption vector x where f is usually called the "free energy density", here given by 4 For an alternative specification, see Appendix D.
with µ the "chemical potential" used to enforce the constraint on average, taken to be identical for all agents and Γ(·) is the Gamma function. Importantly, this grandcanonical description allows one to have a free energy density f (x) that is a sum over entirely decoupled goods. The only coupling is through the value of µ. Given that N → ∞, the shape of the free energy then completely determines the state of the system. Indeed, for large N the partition function can be estimated using Laplace's method, such that the values x * i = x i that minimize the free energy are overwhelmingly more probable than any other values. Setting c = 0 and solving the set of equations ∂f ∂xi = 0, one can for example check that the previously obtained solution is recovered.
As in statistical physics, we expect to identify the phase transition from the uniform to the condensed (concentrated) phase where herding dominates [6,43]. This occurs when the single minimum in free energy associated to the non-condensed solution becomes a maximum while one or several new minima appear. Such a change of topology occurs for some value c crit , as illustrated in Fig. 1(b) in the case where there are only two products. For given values of c, the depth of the minima in the concentrated phase will depend on the a i and p i . The most favorable configuration (in this case having more of the least expensive of the two products, as here a i = 1, ∀i) is associated to a lower free energy.

B. Concentration for β → ∞
In order to find precisely when and how concentration occurs, we first study the case of fully rational agents, β → ∞. Carefully rescaling the free energy density, we find that an extremum is reached for the configurations x * satisfying with the value of the chemical potential being such that As previously described, the critical value c crit = c ∞ where this transition occurs in the fully rational case can be found by looking at the Hessian of the free energy evaluated at the non-condensed solution. Doing so (see Appendix C 2), and choosing for the sake of simplicity a i = p i = 1 for all agents, one finds As long as the right hand side is strictly positive, there will therefore be a value of c above which concentration occurs when β → ∞, while if this is not the case concentration never occurs, regardless of the strength of interactions. The resulting theoretical phase diagram is shown in Fig. 1(c) for k = 2, a value that we shall keep fixed henceforth. This phase diagram is in perfect agreement with our numerical simulations, as shown in Fig. 1(a). Note that this procedure can be repeated for different values of a i and p i , leading qualitatively similar results. Since the Hessian of the free energy is diagonal in our case, the critical value c ∞ would then correspond to the first change of sign of a diagonal element of the matrix.

C. Finite β Effects
Placing ourselves in the regions where condensation does occur in the fully rational limit, we now set out to understand how bounded rationality might alter the phase transition. In the general case, analytical expressions are difficult to obtain. However, numerically finding where the Hessian (which is still diagonal in i, j) loses stability for the non-condensed solution yields a semianalytical critical line in (c, β) space. It should be noted that this condition, explicitly given in Appendix C 2, is independent of the value of µ, suggesting that the loca-tion of the transition is identical in the canonical and grand-canonical descriptions.
Our analytical result can be compared to numerical simulations, for which the transition to fashion dominated consumption is identified by looking at the rescaled Herfindahl index This index takes the values 0 and 1 in the uniform (x i = w/(M p i ), ∀i) and fully concentrated (x i = w/p i for one product and = 0 for all others) cases respectively. As shown in Fig. 2(a), the numerical phase diagram and the theoretical critical line match very well, despite the fact that the semi-analytical calculation is based on the grand-canonical ensemble, whereas numerical simulations strictly enforce the budget constraints for all agents. Furthermore, using the theoretical values for c crit (β) to rescale the evolution of the average basket as a function of c/c crit (as plotted in Fig. 2(c)), we observe that the evolution of the mean basket and related quantities appears to be largely independent of β.

IV. CONSEQUENCES ON THE SLUTSKY MATRIX
Using a simple interacting model, we have shown that introducing herding in the problem leads to a concentration transition and radical changes in the way agents allocate their budget among the M available goods. We now set out to evaluate the impact of such a transition on the Slutsky matrix and its properties, in particular its negative semi-definiteness and symmetry.
Whereas the grand-canonical theory allowed us to calculate x α i = x * i for the entire range of β and c selfconsistently, the absence of an explicit expressions prevents us from directly computing S α ij . Instead, because we have found that results are largely independent of β provided c is rescaled as c/c crit , we may gain insight from the Gaussian approximation of the Slutsky matrix introduced in Section II C and valid for β → ∞. Equation (23) can now be made explicit and writes (see Appendix C 3) find that C has the following form where ϕ ij and ψ ij are both symmetric and O(1) in N , and, of course, depend on x * . Such correlations, which can be checked to match numerical simulations very well, finally yield an identical Slutsky matrix for all agents, Together with Eq. (13), this is a key result of our work. Thus far, we have considered the case where all agents and goods are identical. In the present context, however, it is more illustrative to break the symmetry between products by introducing heterogeneous p i and a i , leading to clear preferences between products. An example with M = 4 goods is given in Fig. 3, showing a very good agreement between the fully rational theory and numerical simulations with finite β and N sufficiently far from the transition region. 5 Upon inspection of opposing entries, it quickly appears that the symmetry of the matrix is not satisfied near the critical value of c where condensation first occurs. Note that all matrix entries quickly become very small once the system has entered the concentrated phase. This can be understood for large c by making the ansatz x * i = w/p i − 1 c j =i y j for the dominant product and x * j = y j /c for others. Plugging such a guess in Eq. 29 indeed solves the equations, and predicts S ij ∼ 1/c to leading order in c −1 .
The first property of the matrix that interests us is its spectrum, and in particular the non-positivity of its eigenvalues. As shown in Fig 4(a), the fully rational theory provides a very good description of the matrix eigenvalues, which remain non-positive for the entire range of c. The leading eigenvalue actually peaks close to the transition. Consistent with the decay of the matrix entries themselves, the magnitude of the eigenvalues also vanish as c increases beyond the transition. As a result, our theory and numerical experiments show that the Slutsky matrix does remain negative semi-definite for all value of c and β. The main consequence of the herding 5 Close to the transition, agents flip-flop between different basket compositions ( Fig. 2(b)), which leads to anomalously strong fluctuations and corrections to the Laplace saddle point method used in all our analytical calculations. This is a well known effect in statistical physics, which leads to interesting phenomena in their own right, but that we do not explore further here.
transition on the spectrum is the decay in the magnitude of eigenvalues, as confirmed by looking at the trace of the matrix shown in Fig. 4(b). The other essential property of the Slutsky matrix is its symmetry. Due to the decay in the magnitude of the matrix entries once the system has entered the herding phase, this property becomes difficult to measure as c is increased beyond c crit . Indeed, as both S ij and S ji become very small, any finite numerical error ε affecting either entries will result in a very large relative asymmetry, at which point most common measures of asymmetry will fail. To minimize the impact of such numerical errors on the conclusions of our study, we propose an asymmetry measure χ defined as which should be equal to zero for symmetric matrices, and diverge in the anti-symmetric case. Employing this metric with the β → ∞ theoretical Slutsky matrix, we find that the matrix is very close to (but not exactly) symmetric far from the condensation transition. In the vicinity of c crit however, strong interactions give rise to a significant value of χ, contradicting the conventional lore even in the fully rational case. This theoretical result, compared with numerical results, is shown in Fig. 4(c), while also visible in Fig. 3. Note that, as expected given our choice of identical agents with identical budgets, the previously introduced aggregate Slutsky matrix S remains symmetric even as the concentration of choice occurs. The theoretical results for this alternate definition are shown by the dashed lines in Fig. 4, where it is clear that the asymmetry measure χ is always zero. Interestingly, the eigenvalues of the matrix are very similar for the individual (S ij ) and aggregate (S ij ) definitions.
Regardless of the metric, Fig. 4 also illustrates the discrepancy between the equilibrium theory we have devised and the numerical measurements in the transition region. As expected, the system indeed takes a very long time to reach the Boltzmann-Gibbs distribution when the transition occurs. These non-equilibrium effects are also reflected in the difference between the numerical measurements obtained with finite differences and those calculated using the fluctuation-response relations and the associated "thermodynamic" expression of the Slutsky matrices. We expect that such effects will also be present in real empirical data, specially if herding effects bring the system close to a transition point, as seemed to be the case in the Salganik et al. experiment [40,43].

V. DISCUSSION
Before concluding our study, let us examine two subtle points that the above analysis has treated in a somewhat cavalier way.

A. Global vs. Individual Utilities
From the very start of our statistical mechanics description of the multi-agent problem, we have emphasized that we are assuming that individual agents change their basket of goods according to the change of the global utility of the population rather than of their own utility. In other words, agents also take into account the change of utility of others when they update their choices. The main motivation behind this specification is that the dynamics will then spontaneously reach a Boltzmann-Gibbs equilibrium measure. Indeed, following an agent-based framework where individuals set out to improve their own utility function, including the herding component c(x) k would bring us in the realm of non-Hamiltonian dynamics, for which general analytical tools are still lacking. Natural questions are now (i) how may the current model still be interpreted at the agent level and (ii) how is the phenomenology of the system impacted if we abandon the Hamiltonian description and take agents to maximize their own utilities.
Regarding the former question, we start by writing the change in global utility ∆U if the randomly selected agent γ in the Monte Carlo dynamics accepts the proposed basket of goods x γ + ∆x. Assuming N 1, we have where we remind the reader that the overline notation refers to arithmetic averages over the agents. This equation shows that the dynamics can in fact be interpreted in terms of agents only concerned with their own utility, albeit with modified values of the interaction parameter c (c → c i = c[1 + klog x i ]). The individual utility thus becomes configuration dependent, but only through aggregate quantities (the average of the logarithm of x i ). This logarithm correction does lead to a supplementary penalty for products that have lost the favor of the crowd, but is not expected to radically change the concentration phenomenon reported in the above sections. Consistent with this interpretation, simulating the system with a decision rule based on the individual utility but with a constant c -i.e. placing ourselves in a non-Hamiltonian setting -indeed results in a largely unchanged phenomenology, with the same herding transition as was observed in the Hamiltonian case. As expected, the absence of this logarithmic penalty for small x i pushes the transition to slightly higher values of c and leads to more volatile Monte Carlo trajectories, with the appearance of more frequent switches in the vicinity of the transition. Although the absence of a solid theoretical framework to describe the steady state in that case prevents us from drawing definitive conclusions at this stage, we conjecture that most of the results obtained at equilibrium regarding the Slutsky matrix continue to hold, with only minor quantitative modifications.
Alternatively, one could also purposefully write a global utility for which the detailed-balance condition matches the maximization of an agent-specific utility. To do so, the interaction term must be entirely symmetric in the sense that the change in global utility is identical regardless of the randomly selected agent (which was not the case with the previously studied model). For instance, one could take with J i a symmetric interaction matrix. The mean-field approximation of this model J αγ i = a i J/N is studied in the limit β → ∞ in Appendix D. For ρ > 1 2 , concentration will occur for sufficiently large values of J, and we therefore expect our results for the Slutsky matrix properties to remain largely unchanged.
In the case where J i is not symmetric, the problem may no longer be treated with the standard techniques of equilibrium statistical mechanics. The system can then still be simulated with a decision rule maximizing the individual utility of the agent, but we can no longer assume that the system reaches an equilibrium given by the Boltzmann-Gibbs distribution.

B. Equivalence of Ensembles
In order to study analytically the equilibrium properties of our interacting set of agents and to determine when the herding transition occurs, we had to relax the budget constraint and place ourselves in the "grand-canonical" ensemble. As previously mentioned, the equivalence of results between the canonical and grand-canonical ensembles is not guaranteed a priori, although in the present case the analytical (grand-canonical) results appear to match (canonical) simulations extremely well.
To formally assess the possible differences between the two ensembles, the budget fluctuations in the grand-canonical ensemble must be studied. The ensemble equivalence corresponds to cases where the variance of the realized budget are vanishingly small. Hence, we set out to compute in the grand-canonical ensemble. This computation may be performed by introducing small heterogeneous perturbations to the chemical potential, µ → µ + δµ α and differentiating log Z N with respect to δµ α . The calculation, detailed in Appendix C 4, finally yields The evolution of this quantity for different values of β as a function of c in the previously discussed case a i = p i = 1 is shown in Fig. 5. As expected, this quantity vanishes for β → ∞, as well as for c → ∞. In the c < c crit region, where fluctuations are expected to be the strongest and thus where the difference between the ensembles is expected to be the largest, the non-condensed solution (dotted lines on Fig. 5) gives σ 2 ∼ O(M −1 ). As such, the grand-canonical and canonical ensembles will become strictly equivalent only in the limit M → ∞. This being said, the Hessian of the free energy being independent of the chemical potential suggests that despite this apparent absence of strict equivalence between the ensembles for finite M , the onset of the transition appears to be largely unaffected. In any case, although precise measures of the fluctuations of the system in the grand-canonical ensemble should be inaccurate for smaller values of β, this difference will not change the key results presented in this work. Besides, one should expect that some amount of budget fluctuations are indeed present in the real world!

VI. CONCLUSION
Let us summarize what we have achieved in this paper. By introducing a rationality parameter (or "intensity of choice") β to account for the fact that agents are not strict utility maximizers, we have first reformulated the Slutsky equation within a general "fluctuation-response" framework, which allows one to express the Slutsky matrix in terms of consumption fluctuations only, without having to measure changes of consumptions when prices are slightly modified.
We have then shown that irrationality does not necessarily result in the breakdown of the symmetry of the Slutsky matrix. As a result, the hypothetical symmetry of empirically measured Slutsky matrices cannot be used as a general argument against bounded rationality [36].
When accounting for herding within large assemblies of agents, we found that symmetry is no longer guaranteed in general. Introducing a simple model of utility with interactions, we have indeed shown using the powerful methods offered by statistical mechanics that a concentration transition may occur, at which point strong selection of goods occurs. At the transition, the individual Slutsky matrix becomes markedly asymmetric, although the Slutsky matrix constructed using aggregate consumption can still remain symmetric when all agents are identical. Hence our result is not necessarily incompatible with existing work on non-unitary households [27]. From simulations, we also found that asymmetry is further amplified by out-of-equilibrium effects near the critical point. In line with standard consumer choice theory and most empirical studies, our model preserves the negative semi-definiteness of the Slutsky matrix regardless of interactions. Nonetheless and although not studied here, it should be noted that introducing some interactions in between products (representing redundancy or complementarity for example) in a similar framework appears to lead to positive eigenvalues [45].
Central to our theoretical framework, the assumption that the system may be written in Hamiltonian form, meaning agents are not strictly individualistic when interactions are included, does not appear to be necessary for our results to hold. This being said, further investigations in particular regarding imitation-lead switches in an out-of-equilibrium model would be a natural path to explore, although we expect doing so may only lead to stronger asymmetry of the Slutsky matrix. Considering different forms of utility and their impact on the symmetry of the Slutsky matrix is another potential aspect to expand upon. In particular, adding disorder in the interactions (for example by taking J i to be a random matrix in Eq. (36)) would likely lead to very interesting effects.
Of course, further empirical studies on the properties of the Slutsky matrix would be of great interest, in particular to contrast our model of interacting, bounded rational agents with the recent sparsity-based approach of Gabaix [31]. To this end, we believe that fluctuation-response relations such as those presented in this work could be a very valuable econometric tool. Given the difficulty of conducting repeatable experiments in socioeconomic systems, the ability to estimate derivative quantities without actually requiring prices to change appears quite promising.
Appendix A: General "Thermodynamic" Relations

Single Agent Slutksy Matrix
Assuming the equilibrium distribution is reached, Eq. (6) can be rearranged to be rewritten through correlation functions. Indeed, from the quotient rule the first term becomes Now, using the fact that then for a test function O(x) Therefore bringing everything together we have or in a more compact form, with xixj c = xixj − xi xj and Γ = ∂w log Z,

Many Agent Slutsky Matrix
The agent-specific Slutsky matrix is defined as Once again assuming equilibrium is reached, we wish to rewrite this expression through correlations to retrieve some information on the symmetry of the matrix. As before, we look at the first term in the above equation, with hereafter the shorthand notation Dx := + α,i dx α i . Now, by the generalized chain rule we have which may be plugged back into Eq. (A7). Noticing that with Γγ = ∂ ∂w γ log ZN , then Bringing everything together, we finally have When there are no interactions, γ = α, x α i x γ j c = 0, as well as ∂ ∂w γ x α i = 0. As a result, we are simply left with which is symmetric and identical to the single-agent system, as expected.

Fluctuation-Response Relations
We start from the most general expression for the partition function and set-out to find a readily measurable expression for Γα = ∂ ∂w α log ZN . Starting from the computation of we employ the generalized chain rule, Now, regardless of the assumption on their distribution, the agents' individual budgets are not related, and as such ∂w γ ∂w α = δαγ. Reinjecting in the previous expression and integrating by parts, we obtain which finally gives The other terms in the thermodynamic Slutsky matrix expression can be simplified in similar ways. Indeed starting from and reusing Eq. (A15) and subsequent iterations by parts, it is straightforward to show that from which we can also derive

Aggregate Slutsky Matrix
We now take with w the arithmetic mean over the agents' budgets. Now, the above can simply be rewritten as The fluctuation-dissipation relation for the first term still holds as before, while for the derivative with respect to the mean budget we have, in general, with κ η := ∂w η ∂w .
As a result, we may write Bringing both expressions together and changing indices appropriately, we find In the limit β → ∞, we have shown that the second term in the sum vanishes. In this case, we therefore recover a symmetric Slutsky matrix whenever x γ j = κ γ xj. This is for example the case when the budgets are proportional to the mean i.e. w γ = κ γ w and all agents have identical preferences. Now, in a simple model where w γ = κ γ w, ∀γ, and all agents have identical preferences, then x γ j ≡ κ γ xj and the second contribution to Sij identically vanishes, so which is always symmetric when all agents are identical.
Interesting things may happen when agents are more heterogeneous, or in a model where a change of average wealth does not affect all agents in an affine way. Imagine agents are characterized by a label κ α ∈ R such that where w0 is a certain fixed wealth scale, w a varying parameter and F a possibly non-linear function. We must impose that For example when F (x) = x one finds the above linear model: More generally, one has ∂w ∂w = κF κw w0 and ∂w α ∂w = κ α F κ α w w0 , so in general ∂w α ∂w = w α w and the resulting Slutsky matrix has no longer obvious reasons to be symmetric.

Gaussian Fluctuations
We take x α − x α = δx α and set out to calculate x α i x γ j c = δx α i δx α j . At this stage, for simplicity, we combine the N × M degrees of freedom in the single column vector v = [x 1 1 , . . . , x 1 M , . . . , x N 1 , . . . , x N M ] such that the correlations of interest are given by with the resized price vectorp ∈ R M ×N ,p k =μ α pi, k = M (α − 1) + i, and where the product of budget constraints has been rewritten using the Fourier representation of the Dirac δ. In this Gaussian approximation, the partition function can be calculated exactly with two consecutive integration by parts, Going back to Eqs. (A28), completing the square in the exponent gives the change of variable δu = δv + i β H −1p . The integral then reads The first term is straightforward, and simply gives − 1 β (H −1 ) k 1 k 2 as without the constraint. The second term requires more care, as we have As a result, Bringing everything together, we finally have which, after replacing with the correct indices, is the expression given in the main text.
Appendix B: Numerical Methods

General Idea
To measure equilibrium properties of the system, we choose to follow a Monte Carlo approach, although a stochastic gradient descent of the global utility (i.e. Langevin dynamics) should be equivalent. The system is first initialized by randomly allocating the N baskets following a uniform distribution. We then employ a Metropolis-Hastings algorithm (see next subsection), as a modification ∆x to a randomly selected agent's basket is proposed and accepted with a rate verifying detailed balance [46]. To satisfy the non-negativity of the baskets and to correctly sample the small x α i region if high concentration occurs, the move is actually constructed in terms of logarithms. For a randomly selected agent γ, we take log(x γ + ∆x) = log x γ + ξ, with ξi ∼ N (0, 1), resulting in a log-normally distributed multiplicative noise. The budget constraint is then enforced by a simple rescaling of the resulting basket of goods. Here, it is essential to note that the proposal distribution is not symmetric. As a result, the asymmetry must be accounted for in the acceptance rate, giving in our specific case with ∆U the change in the global utility 6 caused by the move (see next subsection). From each run, average quantities such as x α are measured by taking arithmetic means in algorithmic time once the system is equilibrated.

Metropolis-Hastings Acceptance Rate
In its most general form, the Metropolis-Hastings acceptance rate is given by where π(x) is the probability density function we wish to sample, and q(y | x) is the conditional probability of proposing the state y given the current state x. In our case, and as usual in statistical mechanics, we simply have π(y)/π(x) = e β∆U . The ratio of conditional probabilities requires more care. Given the noise yi = xie ξ , with ξ a Gaussian noise, the probability density of which will be noted ρ, we have and very similarly Given the symmetry of ρ for zero-mean noise as is the case here, and the fact that all dimensions are statistically independent, yielding the expression in the previous subsection.

Computing the Slutsky Matrix
To compute the entries of the Slutsky matrix, one can use the fluctuation-response relations derived here, which do not require derivatives. However, out-of-equilibrium effects can become significant, in which case these relations do not hold. For this reason, it is essential to compute as a check the original expression given in Eq. (12). This is not entirely straightforward, as taking discrete derivatives after time averaging may induce some biases. Instead, we rely on so-called pathwise derivative estimates, often used in mathematical finance to compute risk sensitivities [47,48]. In a nutshell, we generate perturbed trajectories based on the original simulation that will use the same random numbers (both for proposed moves and acceptance), and measure finite differences at every step.

Appendix C: Interacting Model
We consider the global utility

Non-interacting Limit -Canonical Ensemble
We start by solving the much simplified non-interacting limit, setting c = 0. When there are no interactions in general, the canonical partition function entirely decouples and can be written as a product, Now, plugging in the desired utility function and rewriting the Dirac delta with its integral representation, we have This integral now decouples in a product over i, and performing the integral on x α i with the appropriate change of variable, where we have used the known inverse Laplace transform of a power law. Taking the product over all agents then gives the result in the text for the partition function.
To find the average x α i , we proceed directly, by computing i.e. exactly as before except for the agent and product considered, for which we must now calculate where we have used the property Γ(x + 1) = xΓ(x). The integration over µ is then identical with only the change in the exponent. At this stage, we can plug back everything together, all γ = α, j = i terms will cancel out from the ZN at the denominator, and we are left with Thanks to this explicit form, the agent-specific Slutsky matrix entries read which, as predicted by thermodynamic relations, is clearly symmetric, and where the heterogeneous budgets only act a a prefactor altering the magnitude, such that we can write S α ij = w α Kij. The distribution of S α ij is then simply a rescaling of the original distribution of w α in the large N limit.

Finite Interactions -Grand-Canonical Ensemble
To tackle finite interactions, we need to place ourselves in the grand-canonical ensemble, where calculations are possible. Enforcing the value of xi to be the arithmetic mean over all agents with a Dirac distribution expressed in integral form, we have We may begin by integrating over x α i . Given the positive support and logarithmic utility, the result is simply a Gamma function for the N identical agents, The entire integrand can thus be rewritten as an exponential with N as a prefactor, a form which naturally points to the saddle-point approximation as N → ∞. Indeed, the integral over λi above may first be be rewritten as approximated as where λ * is a minimum of g, and H * is the Hessian with respect to λ evaluated at that point. At the saddle, we have and the Hessian is entirely diagonal and given by As a result, the partition function can finally be expressed in the desired form, with F (x) = N f (x), and where from the above equations we have which is completely decoupled in between degrees of freedom. Keeping only O(1) terms, we find the expression given in the main text that we will use onward.
In the thermodynamic limit, the state of the system can once again be determined through the saddle-point approximation, as the system will reach a minimum of f with overwhelming likelihood. This amounts to solving with ψ the digamma function. To determine the nature of stationary points, we will also require where ψ (1) is the first polygamma function. Remarkably, the second derivative is not dependent on the chemical potential µ. Importantly, one can evaluate the partition function by steepest descent to show that, consistent with intuition, when N → ∞, straight from Eq. (C16) as all derivatives other that those with respect to pi are zero at the saddle by construction.
a. c = 0 solution Taking c = 0, we easily recover the unconcentrated solution b. β → ∞ limit In order to better understand the influence of the different terms on the equilibrium solution of the system, we start by studying the zero temperature, or fully rational, limit of the system. Knowing the asymptotics of the digamma function ψ(z) ∼ log z − 1 2z , we obtain i.e. a set of M + 1 equations for the M variables x * i and the correct value of µ for the given budget w. To obtain the solution as a function of c, one can numerically solve the system of equations.
Knowing what the uniform solution is, we may then identify where the phase transition occurs, i.e. what the value of c∞ is, by studying the stability of the free energy at that point, as it should go from being a minimum to a maximum at the transition. Now in this limit, we have Setting to zero at the saddle, and solving for c∞ finally yields In the case where all products are equivalent, we recover the expression in the main text.
c. Finite β In the general case, the equations for x * i are much harder to manipulate due to the presence of the highly nonlinear logarithm and digamma functions. However, the method outlined to obtain the critical value of c remains valid: the c = 0 solution is to be plugged into the second derivative that is set to zero. Solving numerically (Eq. (C18)) then yields the critical line in (c, β) shown in the main text and in Fig. 2.

Slutsky Matrix in the β → ∞ Limit
We wish to employ the previously devised Gaussian approximation for the Slutsky matrix at β → ∞ for our interacting model. To do so means evaluating the "thermodynamic" expression of the Slutsky matrix that requires computing the cross-agent term x γ j ∂ ∂w γ x α i . For our specific choice of utility, ∂U ∂x γ j = aj x γ j [1 + c(xj) k ] + kcaj(xj) k−1 log xj, with overlines indicating arithmetic averages over the N agents. As we are interested in the N → ∞ behavior, we can first notice that in the vicinity of the maximum by virtue of the central limit theorem. Similarly, To leading order, we therefore have As a result, developing to the second order in δx γ so plugging back this expression in the fluctuation-response relation gives Given the Gaussian fluctuations scale as x α i x γ j c ∼ β −1 , this term will clearly have a non-negligible contribution to the Slutsky matrix at β → ∞.
Using the fluctuation-response relation to compute the first term in the Slutsky matrix, we can finally write, choosing = j There now remains to compute the pairwise correlations, which requires us to invert the Hessian of the utility function. In our specific case Due to the homogeneity in between agents, this M × N matrix has the very simple structure with A, B ∈ R M ×M , both having the additional simplification of being symmetric matrices. We take the ansatz that the inverse has an identical structure, obtained from the multiplication of these simple block matrices. As a result, setting HH −1 = I amounts to Going back to Eq. (C32) (evaluated at the maximum x * ), we have and The diagonal matrix can be inverted (almost) exactly in the large N limit, yielding and similarly at the leading order finally giving Using the inverse of the Hessian, we can then compute the correlations of interest. We first need which can clearly be split in diagonal and off-diagonal parts, therefore the matrix has a similar structure to the Hessian, but reduced in dimension. We have G = aI + bE with I and E respectively and identity matrix and a matrix full of ones, both of dimension N × N . Taking the same ansatz as before but with scalars, we have , a = p Fp, b = p Dp (C43) Bringing everything together, we finally find, Ultimately, in the N → ∞, β → ∞ limits, the Slutsky matrix is thus given by The evolution of this analytical expression for some a, p with c, compared with numerical simulations, is shown in Figure 3 of the main text.