key: cord-0147071-1rau2arh
authors: Halperin, Igor
title: Phases of MANES: Multi-Asset Non-Equilibrium Skew Model of a Strongly Non-Linear Market with Phase Transitions
date: 2022-03-14
journal: nan
DOI: nan
sha: 235f31618c98114d89ab3cf694383d1ea4d4bd95
doc_id: 147071
cord_uid: 1rau2arh

This paper presents an analytically tractable and practically-oriented model of non-linear dynamics of a multi-asset market in the limit of a large number of assets. The asset price dynamics are driven by money flows into the market from external investors, and their price impact. This leads to a model of a market as an ensemble of interacting non-linear oscillators with the Langevin dynamics. In a homogeneous portfolio approximation, the mean field treatment of the resulting Langevin dynamics produces the McKean-Vlasov equation as a dynamic equation for market returns. Due to the strong non-linearity of the McKean-Vlasov equation, the resulting dynamics give rise to ergodicity breaking and first- or second-order phase transitions under variations of model parameters. Using a tractable potential of the Non-Equilibrium Skew (NES) model previously suggested by the author for a single-stock case, the new Multi-Asset NES (MANES) model enables an analytically tractable framework for a multi-asset market. The equilibrium expected market log-return is obtained as a self-consistent mean field of the McKean-Vlasov equation, and derived in closed form in terms of parameters that are inferred from market prices of S&P 500 index options. The model is able to accurately fit the market data for either a benign or distressed market environments, while using only a single volatility parameter.

When financial practitioners and academics talk about the behavior of the market, they normally refer to the price behavior of stock market indexes such as the S&P 500 or the Dow Jones index. For both these market indexes (and other similar ones), returns are defined as weighted averages of returns of their constituent stocks, the differences being in the chosen universe of stocks and the weighting scheme. In particular, the S&P 500 index weighs all stocks by their total market capitalization, and is therefore heavily influenced by mega-stocks.

Given that a market index such as the S&P 500 (the SPX index) is a weighted average of individual stock prices, one might be tempted to assume that the dynamics of market returns would be similar to dynamics of a single 'representative' stock in the market. Such interchangeability of modeling returns for the whole market versus modeling returns for single stocks is commonly assumed by both academics and practitioners alike. In particular, in both the theory and practice of derivatives pricing, the same models such as stochastic volatility models are used for either single stocks or market indexes, albeit with different parameters.

However, the dynamics of the mean return of all stocks (i.e. the market return) can only be simply expressed as the mean of dynamics of individual returns if these dynamics are linear. In a more general case, we may think of a market as an ensemble of individual stocks whose individual dynamics are generally non-linear due to market friction effects as discussed in more details below. Furthermore, these non-linear dynamics for individual stocks are not independent of each other (again, specific market mechanisms producing such co-dependencies will be presented below). Therefore, in general the market dynamics should be viewed as statistical mechanics of an interacting ensemble of self-interacting nonlinear 'particles' representing individual stocks. Quantities such as market returns would be computed in such a framework as ensemble averages.

For non-linear interacting systems, ensemble averages do not in general reduce to some sort of expectations within a single particle dynamics. Therefore, a relation between dynamic properties of a markets index such as SPX and properties of a 'typical' or 'representative' stock is in general unknown. However, in statistical physics there exists an approach that essentially constructs such effective singleparticle dynamics starting with an initial multi-particle interacting system. This approach, or rather a family of approaches, are known in physics as mean field approximations, or MFA for short. Mean field approximations used in statistical physics are known to become exact in the limit when the number of particles in the systems goes to infinity, → ∞, known as the thermodynamic limit. As the S&P 500 index is composed of = 500 stocks, such value of might be already sufficiently large to make the MFA qualitatively or even quantitatively accurate.

In this paper, I present a simple and tractable model that builds the dynamics of a market index as the dynamics of the mean log-return of all stocks in the market. In other words, the market index is identified with the mean field in an ensemble of stocks that compose the market index. The way the MFA is constructed and used in this work is similar to how it is used in one of its most canonical applications to the classical Ising model, where the mean field is identified with the mean magnetization of Ising spins, and computed as a solution of a certain self-consistency equation (see e.g. [16] or [18] ).

Similarly, in the model presented below the MFA gives rise to the mean field (i.e. the equilibrium expected market log-return) as a solution of a particular MFA self-consistency equation. Differently from the Ising model that deals with binary spins without self-interactions, the model developed below deals with continuous real-valued non-linear (i.e. self-interacting) oscillators as building blocks representing individual stocks. Nonetheless, the mean field approximation applied to the model of this paper similarly produces a self-consistency equation whose solution gives the predicted equilibrium market log-return. Moreover, many other details of the model developed in this paper will also bring strong analogies with the Ising model, including the phase structure of the model that admits both first-and second-order phase transitions in different regimes of parameters, and the values of critical exponents.

In this work, the market index portfolio is considered as an ensemble of stocks with individual logreturns ( = 1, . . . , ). They can be viewed as 'particles' with 'positions' . The key input to modeling ensembles of interacting particles in statistical physics is a potential function ( 1 , . . . , ). The choice of the potential defines all further properties of a model, including both static and dynamic properties. In general, a potential ( 1 , . . . , ) can be decomposed into a sum of self-interaction potentials and interaction potentials. For example, if only pairwise interactions are allowed, the decomposition takes the following form

where ( ) is a self-interaction potential for stock , ( − ) is a pairwise interaction potential, and is a coupling constant that regulates the strength of interactions in the system.

In this paper, I motivate the choices for both self-interaction and interaction potentials using the analysis of money flows and their price impact in a multi-asset market. This approach follows the previous work by the author [9, 10] 2 where non-linear models of a single stock price dynamics were obtained starting with a similar analysis of money flows and their impact. 3 As was shown in [9] , a combined effect of money flows from external investors and their price impact gives rise to a non-linear (cubic) drift ( ) in the diffusion law for the stock price. This translates into stochastic Langevin dynamics where diffusion in the price space is described as a Brownian motion of a particle that is additionally subject to an external non-linear potential ( ). The latter is defined to satisfy the relation ( ) = − / , therefore a cubic drift ( ) translates into a quartic polynomial as a model of the potential ( ).

In [10] , a closely related model called the Non-Equilibrium Skew (NES) model was presented for the log-return space. Unlike more traditional approaches in physics where a potential is an input and stationary or non-stationary (transition) probability distributions are outputs, the NES model starts with a parameterized model for a stationary distribution, chosen to be a square of a simple two-component Gaussian mixture. The potential in the NES model is given by a negative of a logarithm of this Gaussian mixture. This provides a flexible and highly analytically tractable self-interaction potential ( ) with five parameters, which can produce either single-well or double-well potentials that have, respectively, either one or two local minima. The outputs of the model are transition probabilities for a pre-asymptotic regime, before settling to an equilibrium steady state (which is in fact an input to the model as mentioned above) in the long run. As was shown in [10] , these pre-asymptotic, non-equilibrium corrections to the asymptotic steady-state return distribution impact estimated moments of the return distribution such as variance, skewness and kurtosis -which explains the name of NES model. This paper presents a multi-asset extension of the NES model, to be referred to as the MANES model. 4 In this framework, the NES potential serves as a self-interaction potential ( ) in Eq.(1), while the interaction potential ( − ) is found to be quadratic ( ) = 1 2 2 . The MANES model can be viewed as a new statistical mechanics model with a highly tractable potential that can describe different dynamics depending on the model parameters. In particular, the self-interaction potential in the MANES model can be of a double-well form, depending on model parameters. In the latter case, the model behavior is similar to the Desai-Zwanzig model of interacting double-well anharmonic oscillators [3] .

As will be shown in detail below, the approach developed in this paper offers a number of insights. First, it links the dynamics of market returns with the dynamics of an equivalent single stock that arises in 2 See also [11] for a non-technical presentation. 3 The critical role of money flows and their impact on asset returns were highlighted in recent important papers by Gabaix and Koijen [6] and Bouchaud [1] . 4 While the NES model developed in [10] is a single-stock model, its performance was explored in [10] using options on the S&P 500 indexes rather than single stocks. In this paper, the single-stock model proposed in [10] will be properly used as a building block of a multi-asset MANES model. the mean field approximation to the multi-particle dynamics of a market made of stocks. Second, the MFA approximation applied to a multi-particle interacting system with the potential (1) produces a nonlinear extension of the classical Fokker-Planck equation called the McKean-Vlasov equation [2] . Because the McKean-Vlasov equation is non-linear, it leads to ergodicity breaking and phase transitions [2] . As the mean field approximation employed in the McKean-Vlasov equation is accurate in the thermodynamical limit → ∞, this suggests that the dynamics of the market index can be successfully modeled as the mean field dynamics for a large ensemble of particles/stocks. This approach is able to produce a large set of dynamics scenarios. In addition to a benign market regime of small fluctuations around some stationary or time-varying deterministic level (trend), the model also admits regimes of large fluctuations involving both first-and second-order phase transitions that can be realized in different parameter regimes.

Furthermore, I show how the model can be used in practice by calibrating it to market quotes on options on market indexes such as S&P 500 (SPX) options. Using a homogeneous approximation, the mean field potential arising with the McKean-Vlasov equation describing an ensemble of non-linear oscillators with the NES self-interaction potentials can be represented as an effective single-stock NES potential with parameters that are modified ('renormalized') by interactions. This provides the aforementioned missing link between the dynamics of the market index and the dynamics of a 'representative' stock that mimics the index. When parameters of this effective single-stock NES potential are inferred from the market options data, they are used to compute the equilibrium mean field of the MANES model. The latter is interpreted as a model-based, option-implied prediction of the equilibrium log-return of the market, and can be used as a signal driving asset allocation decisions. Furthermore, other moments of the index return distribution inferred from market option quotes can also be used as predictive signals that can be used for investment decisions.

The paper is organized as follows. Sect. 2 gives the derivation of non-linear stochastic dynamics of a multi-asset market that is driven by money flows and their impact. These dynamics are then reformulated as multi-particle Langevin and Fokker-Planck dynamics in Sect. 3. Sect. 3.2 introduces the NES potential as a tractable approximation to a non-linear self-interaction potential obtained with the approach of Sect. 2. Sect. 3.4 provides a derivation of the McKean-Vlasov equation -a non-linear version of the Fokker-Planck equation for a multi-particle system that arises within a mean field approximation. Sect. 4 derives the self-consistency equations resulting from using the NES potential within the McKean-Vlasov equation, and then explores the phase structure of the model including both first-and second-order phase transitions in different parameter regimes, and computes critical exponents. The next Sect. 5 derives a closed-form relation for the equilibrium expected log-return of the market in terms of parameters obtained by calibration to index options, and considers examples of calibration to the market data on the SPX options. The final Section 6 concludes.

Let x with components be a vector of asset values with = 1, . . . , in an investment universe of assets at the beginning of the period [ , + Δ ], where is the current time and Δ is a time step size. Let us start with a discrete-time asset price dynamics in the presence of outside investors described by the following equations:

where • stands for the element-wise product. Here the first equation in (2) defines the change of asset values in the time step [ , + Δ ] as a composition of two changes to their time-values x . First, at the beginning of the interval, investors adjust positions in each stock by buying (or selling, depending on the sign of ) the amount Δ of the stock (so that the action vector a is defined as an instantaneous rate of change with the dimension of inverse time). Therefore, immediately after that the asset value is deterministically changed to x + := x + a • x Δ . After that, the new portfolio x + = x + a • x Δ grows at rate r +Δ . The latter is given by the second of Eqs.(2) that defines the vector of equity returns as a combination of a risk-free rate , predictors z with weights w, a multivariate noise ∼ N (·|0, ), and a vector-valued market impact factor f (a ) with components ( ). A particular form for function f (a ) will be presented in the next section, while in this section I proceed with a general form of this function. 5 Equations (2) can also be written for increments Δx = x +Δ − x :

where terms ∼ (Δ ) 2 are omitted assuming that the time step Δ is small enough to justify a continuoustime limit Δ → 0. In the strict limit Δ = → 0, Eq.(3) transforms into the following stochastic differential equation (SDE):

where W is a standard -dimensional Brownian motion, and is a × volatility matrix such that = . Next consider a vector y of period-log-returns with = 1, . . . , , defined as follows

Using Itô's lemma, we obtain the SDE for the log-return vector y :

where x − stands for the derivative of the time-lagged asset price x − with respect to the calendar time . Note that without control a and the price impact f (a ), Eq.(4) becomes a standard multivariate lognormal process with the growth rate given by + wz , while Eq.(6) becomes a multivariate normal process with a time-dependent drift due to the term x − /x − . On the other hand, Eq.(2) or its continuous-time limits (4) or (6) describe a controlled system where investors buy or sell stocks at rates at each time step. If investors are rational or boundedrational, their decisions would depend on the previous performance of assets. The problem of finding optimal control a from the viewpoint of market investors could be further formalized by specifying their utility (reward) functions, and then solving corresponding Bellman or Hamilton-Jacobi-Bellman (HJB) equations. Instead of following this route, in this paper I employ a more phenomenological way, and choose a simple functional form of the optimal investment rate a as a function of state variables based on general arguments.

To this end, our specification of investment rates a should encode some simple stylized facts about retail investors. In particular, they usually buy stocks when prices go up (i.e. log-returns y are positive), and sell when prices go down and log-returns are negative. This can be formalized by specifying a parametric function a = a(y ) where a possible dependence on other state variables can be encoded in parameters of this function. Assuming for simplicity that their actions are perfectly asymmetric with respect to such scenarios, it produces asymmetry relations a(−y ) = −a(y ) that should be imposed on all admissible functional specifications a(y ). In other words, a(y ) should be an odd function.

Based on this, the following specification of the investment rate for stock as a simple function of log-return will be used in this work:

where , ≥ 0 are parameters capturing, respectively, a linear and non-linear dependence of investors' allocation decisions on the log-return of asset , and parameter ≥ 0 determines how it depends on the average performance of other assets (i.e. performance of the market). As will be shown below, parameters control non-linear effects, while parameters control interactions in the system. Note that the functional specification (7) can be viewed as a low-order Taylor expansion of a general function, where a constant term and a coefficient in front of 2 are set to zero in order to produce an odd function.

The definition of the investor flow rate given by Eq. (7) both refines the previous similar specifications suggested in [9] and [10] for a single-stock (1D) case, and extends the approach suggested in this previous work to the current case of a multi-asset market with stocks. Furthermore, the model developed below assumes that is large, which seems to be a good assumption if for example we consider the case = 500, by the number of constituents in the S&P 500 index.

To complete the model specification, we need to define the model of market impact f (a t ). First, note that on the grounds of dimensional analysis, as both a and f have dimensions of rates, it means that in addition to the ratios a , the market impact can only depend on some dimensionless variables such as e.g. a ratio of the asset price to average trading volume over a lookback time window [ − , ]. In what follows, the market impact f (a ) will be modeled as a parameterized function of a with constant parameters, thus neglecting such possible additional dependencies. 6 We need a functional specification of the price impact function ( ) that captures the most important effects of impact of market flows on market prices. First, we want to ensure consistency with the presence of momentum in stock prices. As typically a good recent performance of a stock leads to an increased demand, in a short run this typically further increases returns of this stock. Therefore, at least until cumulative flows are small enough, one should expect a positive co-dependence between market flows and asset returns.

However, such positive co-dependence will only persist until the stock becomes "saturated" or "crowded". "Crowding" in a stock occurs when many market participants simultaneously hold large positions in this stock. During periods of market downturns or high volatility, when many market participants simultaneous unwind or reduce their positions in a crowded stock, it creates a further downward price pressure on the stock, producing diminishing or even negative returns on positions in the stock.

To capture both effects discussed above, we need to produce scenarios where money flows initially produce a positive impact on asset returns, but switch to a negative impact once cumulative inflows exceed some threshold value. Such a model would be consistent with the 'dumb money' effect of [5] that predicts that an initial flow into a stock should increase expected returns, but a continuous buildup of inflows into the stock (a 'crowding') leads to diminishing long-term returns.

To capture such saturation effects, a proper impact function should depend on previous inflows. Let be the cumulative inflow rate in the stock over the last periods excluding the current period:

Now consider a simple model that produces an increasing impact for small values of until the sum of and¯does not exceed some fixed valueˆ, and a decreasing impact for larger values:

where > 0 is a parameter. This piece-linear function vanishes at = 0 and = 2 ˆ−¯ , and reaches its maximum equal to ★ = ˆ−¯ at =ˆ−¯. For large values of | |, this functions asymptotically behaves as ( ) ∼ − | |.

While the functional form (9) produces the desired behavior, it amounts to a non-differentiable function. We can consider a soft relaxation of this function that produces a differentiable approximation:

where the function ( ) is defined as follows:

Note that this function provides a soft differentiable relaxation of function ( ) = | |, so that ∞ ( ) = | |. Furthermore, for finite values of , the function ( ) is convex with a minimum at = 0 and (0) = 0. For small values of , the Taylor expansion of ( ) produces

while for large values of | |, we have ( ) → | |. Therefore, the impact function (10) produces the following behavior:

For small values of the argument, we obtain

and for | | 1, we have ( ) = − | |. Both (9) and its relaxed form (10) describe a concave impact function that vanishes at = 0, reaches its peak at = , and then decreases and eventually becomes negative for larger values ≥ 2 . 7 Note for what follows the asymptotic behavior of the combination + ( ) that enters the SDE (6): 7 A similar profile could be reached with an inverted parabola function such as ( ) = ( − ) with parameters , > 0, however I choose (10) over such specification as I want to have an asymptotically linear, rather that quadratic, behavior of ( ).

Using the cubic law (7) for the flow rate, this can also be written as a function of = :

where I retained only a linear dependence on the 'mean field' 1

, and parameters are defined in terms of previously defined model parameters as follows:

3

The functional form of dependence on log-returns found in Eq.(16) will be used in the next section to introduce a nonlinear potential that will be then used for analysis of the joint dynamics of all stocks in the market.

Focusing on an idealized market portfolio made of identical stocks, let us consider the case when all model parameters , , , w · , ℎ etc. are the same for all stocks, so that we write them as , , etc. 8 Furthermore, while parameters defined in Eqs. (17) are generally time-dependent, here I neglect such potential time dependence, and treat them as constants.

The stochastic differential equation (6) that describes the continuous-time stochastic dynamics of any individual stock log-return = in this setting can be represented as an (overdamped) Langevin equation [13] . In the Langevin approach, the diffusion drift is obtained as a negative gradient of a potential function (y ):

where ℎ is a volatility parameter 9 , is a standard Brownian motion, and (y ) is a potential that is decomposed into a sum of self-interaction and pair interaction terms:

Here ( ) stands for a single-stock self-interaction potential, is an interaction constant, and ( − ) is an interaction potential. The Langevin equation (18) describes diffusion of identical particles placed in a common potential ( ), and in addition interacting via a pairwise interaction potential . Using Eq.(19), we can also write Eq.(18) as follows:

In our setting, the choices for potentials , are informed by Eqs. (6) and (16) adapted to a homogeneous portfolio setting (with the coupling constant obtained similarly from Eq. (17)). In particular, interactions are induced by the third term in the control law (7) . The interaction potential ( − ) is therefore given by the Curie-Weiss quadratic potential:

While the interaction potential in Eq.(21) is quadratic, the self-interaction potential ( ) as suggested by Eq.(16) should have higher-order non-linearities:

where is an additional parameter given by the sum of all a -independent terms in Eq. (6) . Therefore, for small values of log-returns, the self-potential (22) behaves as a quartic potential. It also scales as 4 asymptotically for | | → ∞, though with different coefficients from those appearing in its smallexpansion given by the first line in Eq.(22).

Instead of literally using Eq.(22) as the specification of a self-interaction potential ( ) , in this paper I will use a somewhat different potential that, while retaining the important non-linearity of Eq.(22), also offers better analytical tractability, as well as fixes some potential issues with Eq.(22). These issues are related to the asymptotic behavior of the potential ( ) at → ±∞, as will be discussed next.

First, note that the expression (22) implies that in order to have a confining potential that grows as → ∞, we need to have parameter > 1. 10 On the other hand, even if we take > 1, the potential (22) is unbounded from below for large negative values → −∞. Such an unbounded behavior could describe a corporate default or bankruptcy event when the stock price drops to zero, which would be equivalent to the strict limit → −∞. This behavior of the model is obtained as a direct consequence of the specification of the cubic flow rate (7) which produces a strong negative selloff rate for large negative values of log-returns which, together with the price impact function ( ) ∝ − | |, can drive the stock price all the way to zero, or equivalently to the infinite negative log-return → −∞.

Note however that for all practical purposes, instead of being identified with price drops to the strict zero level, corporate defaults or bankruptcies may be associated with events when the stock price drops to a very small but a non-zero value (e.g. a few cents given the previous price of $10, or $1 given the previous price of $50, etc.). In terms of log-returns , such events would correspond to sudden drops to a large negative value. The infinite value of obtained when the stock price hits the exact zero is therefore solely due to singularity of the transformation = log ( / − ) at = 0. If a small but non-vanishing stock price is used as the default boundary, the range of practically admissible values of becomes finite.

The last observation implies that we can replace a potential that may become unbounded as → −∞ with a confining (i.e. growing at → ±∞) potential where corporate defaults would correspond to events of reaching a certain arbitrarily negative but finite threshold ★ < 0. If the probability of reaching this barrier is small enough to ensure consistency with the market data, further details of model behavior for yet smaller values of , which could otherwise differentiate between a confining and unbounded potential, would be immaterial for all practical purposes. This is illustrated in Fig. 1 that shows the original potential ( ) according to Eq.(22) alongside a tractable confining potential referred to as the MANES potential, to be presented in details in the next section. Depending on the parameters, the MANES potential can have one or more minima. In Fig. 1 , parameters are such that the MANES potential is a double well potential with a second minimum at a negative value of , which is not shown in Fig. 1 . The important point to emphasize here is that as long as barriers for both potentials are sufficiently high and have similar heights, both the dynamics of small fluctuations around the potential minimum at 0 and probabilities of transitions from 1 to 2 will be similar in both models. When these probabilities are small, details of behavior in the left tail between the two potentials (an unbounded original potential vs a confining MANES potential) would have a negligible impact on observable consequences of the model. Therefore, instead of literally using Eq.(22) as a specification of the self-interaction potential ( ), this work will assume a confining potential that may be a single-well or a multi-well potential, depending on the model parameters. Confining potentials are easier to work with than non-confining ones, as they give rise to discrete spectra and stationary states in the long run → ∞, which would not exist for non-confining potentials. For a multi-well potential with two or more local minima, the model dynamics would be similar to the Desai-Zwanzig model [3] which chooses a symmetric quartic double well potential as a model of ( ), and its generalized version in [8] that considers more complex multi-well potentials. The next section provides details of the MANES potential that will be used later in this paper for the analysis of dynamics and phase structure of the model.

In this paper, the dynamics implied by the Langevin equation (20) are explored using a simple and highly tractable non-linear self-interaction potential ( ) given by the logarithm of a two-component Gaussian Mixture (GM):

where ( | , 2 ) is a Gaussian density, is a constant that will be defined below, and 0 is another constant chosen such that the minimum value of ( ) is zero. This potential was introduced in [10] within the context of a model for a single stock called the Non-Equilibrium Skew (NES) model. The choice (23) is motivated by both the ease of the analytical treatment and the flexibility of the parametric family described by Eq.(23) which, depending on parameters, describes either a single-well anharmonic potential or a double-well potential. As market flows and their price impact produce a nonlinear potential, the parametric family of potentials presented by Eq.(23) make it possible for the model to capture these effects. As the present work can be considered a multi-asset (MA) generalization of the NES model, in what follows I will refer to the form in Eq.(23) as the MANES potential.

As explained in more detail in [10] , the MANES potential (23) can also be represented as

where Ψ 0 ( ) is the ground state wave function (WF) of a quantum mechanical system corresponding to the classical stochastic dynamics with potential ( ) [10] . 11 In what follows, the WF (24) will be occasionally referred to as the MANES WF. The stationary state of the original classical stochastic dynamics is then given by its square Ψ 2 0 ( ). The constant introduced in (23) is in fact a normalization constant that can be obtained from the requirement that the ground state WF Ψ 0 is squared-normalized, i.e. ∫ Ψ 2 0 ( ) = 1. Note that while Ψ 0 ( ) is proportional to a two-component Gaussian mixture, its square is proportional to a three-component Gaussian mixture:

where the additional third Gaussian component has the following mean and variance:

The normalization condition thus fixes the value of the constant as follows:

The three-component Gaussian mixture density Ψ 2 0 ( ) given by Eq.(25) can be written more compactly using Gaussian mixture weights

This produces the following three-component Gaussian mixture model for the stationary distribution ( ) = Ψ 2 0 ( ):

Examples of trial ground state WFs Ψ 0 and resulting potentials ( ) are shown in Fig. 2 . Figure 2 : The ground state wave function Ψ 0 ( ), the stationary distribution Ψ 2 0 ( ) and the Langevin potential (23) as a function of the log-return , for a few values of the asymmetry parameter , and with different values of parameters 2 , with fixed values 1 = 0.2, 1 = 0.4, 2 = −0.4, = 1. Graphs on the left describe a healthy stock where the right minimum is a global minimum. Graphs on the right correspond to a severely distressed stock, when the left minimum becomes a global minimum. Graphs in the middle column describe intermediate scenarios. In particular, when = 0.5 and 1 = 2 , the resulting potential shown in the blue line is symmetric. For a multi-asset setting, this potential gives rise to a spontaneous breaking of the Z 2 symmetry → − .

While this expression produces a non-linear behavior for small positive or negative values of , its limiting behavior at → ±∞ is rather simple and coincides with a harmonic (quadratic) potential:

The fact that the limiting behavior of the potential coincides with a harmonic potential as → ±∞ means that in this asymptotic regime the model behavior is described by a harmonic oscillator, and thus is fully analytically tractable.

For small values 1, the NES potential (23) can provide a good approximation to a quartic potential that we found for small values 1 in the previous analysis, see Eq.(22). On the other hand, the asymptotic behavior of the NES potential (23) which coincides with a harmonic potential as | | → ∞ is different from the asymptotically quartic potential implied by Eq.(22). However, if the 'physics' of the system is determined by a region of small values of , replacing an asymptotically quartic potential by an asymptotically harmonic potential should produce a negligible impact on observable consequences of the model while significantly simplifying the analysis.

As Gaussian mixtures are known to be universal approximations for an arbitrary non-negative functions given enough components, this implies that an arbitrary potential that asymptotically coincides with a harmonic oscillator potential can be represented as a negative logarithm of a Gaussian mixture. We can refer to such class of potentials as Log-Gaussian Mixture (LGM) potentials. In this paper, I only consider a two-component LGM potential. 12 While with general parameters 1 ≠ 2 , 1 ≠ 2 and ≠ 1/2 the potential (24) is non-symmetric under reflections → − , it becomes symmetric, with (− ) = ( ), for the special choice 1 = − 2 = , 1 = 2 = and = 1/2, see Fig. 2 . It is useful for what follows to consider a potential which is only slightly asymmetric. This can be done by considering the following specification of model parameters in (24):

with small parameters , , 1. Using these parameters in Eq.(24), we obtain, to the linear order in the asymmetry,

is a symmetric ground state wave function, and parameter 0 is a linear function of , , :

Eq.(32) shows that a slightly asymmetric potential corresponding to model parameters in Eq.(31) can be approximated for small values of by adding a linear term to the potential ( ) ( ) obtained in the symmetric limit, where the coefficient 0 can be directly computed from the original model parameters.

As the term 0 in Eq.(32) can be interpreted as the contribution of an external field 0 to the potential energy of the oscillator , this implies that the dynamics in a slightly asymmetric potentials ( ) can be approximated by the dynamics in a symmetric potential ( ) ( ) with an additional fictitious external field 0 . This observation will be used below.

Now, after we specified the particular self-interaction MANES potential ( ) given by Eq.(23), we proceed with an equivalent probabilistic approach to the dynamics described by the Langevin equation (20). Such probabilistic method is provided by a corresponding Fokker-Planck equation (FPE). This approach will be presented in the next two sections. Note that in both of them the explicit form of the potential ( ) is not used, and thus all equations in this section and the next Sect. 3.4 are general and valid for an arbitrary confining potential ( ).

The FPE corresponding to the Langevin equation (20) is a linear partial differential equation for the joint probability (y, ) = ( 1 , . . . , , ) of a state y = [ 1 , . . . , ] describing stocks at time , given an initial position y 0 at time = 0:

For applications, the most interesting probability distributions for a system of identical non-linear oscillators are one-particle density (1) and pair-density function (2) defined as follows: 

Note that while a single-particle FPE is a partial differential equation (PDE) for a one-particle density, Eq.(37) is an integro-differential equation that relates two different densities (1) and (2) . Due to multibody interactions whose strength is controlled by the coupling constant , the densities (1) and (2) are coupled in Eq.(37) which in fact represents the first equation in an infinite hierarchy of the BBGKY type, see e.g. [16] .

To proceed, I follow the traditional approach in the physics literature (see e.g. [14] , [8] ), and rely on the mean field approximation (MFA) where the probability density factorizes into a product of single-particle densities:

Therefore, with the MFA, dynamics amounts to a system of independent and identical particles, such that the coordinate of any of them is given by the average 1

of all original (interacting) particles in the system. This approximation becomes exact in the thermodynamic limit → ∞.

Note that the mean field 1 =1 of a homogeneous system of identical non-linear oscillators can be viewed as a reasonable proxy to the market returns, which are usually proxied by returns of the S&P 500 index. Of course, the real stock market is quite heterogeneous, and furthermore different firms are weighted in the S&P 500 index by their total capitalization. Therefore, identifying the mean field 1 =1 for a homogeneous system with market returns is not expected to provide a good approximation for the price dynamics of any individual stock. However, more importantly, the present mean field approach emphasizes the difference between a single-stock dynamics and the group dynamics of a market made of many stocks, while still retaining a link between them and keeping the whole approach practical.

Plugging the MFA ansatz (38) into Eq.(37), we obtain

This non-linear integro-differential equation holds for an arbitrary interaction potential ( − ). For a particular case of the quadratic Curie-Weiss potential (21), Eq.(39) produces the following equation:

This equation is non-linear due to the dependence of the coefficient − in the last term on the density ( , ). The nonlinear Fokker-Planck equation (40) is known in the literature as the McKean-Vlasov equation [15, 2, 4] . Critically important is the fact the unlike the initial linear FPE equation (35) for a finite -particle system, the the McKean-Vlasov equation that describes the thermodynamic limit → ∞ is nonlinear. As a result, it produces far richer dynamics including phase transitions [2] . Note that if we formally treat as an independent parameter, the McKean-Vlasov equation (40) can be viewed as a linear FPE with the following 'effective' potential:

The stationary density for a given value of = is therefore obtained as a Boltzmann distribution with the 'inverse temperature' = 2/ℎ 2 :

This solution should satisfy the self-consistency condition for = = ∫ ( , ):

Once the value of is found by solving the self-consistency equation (43), it is substituted back to Eq.(42) to find the stationary density. The number of equilibrium states in the system is therefore given by the number of solutions to Eq.(43). The solution of this equation for the MANES potential will be presented below in Sect. 4.

As was noted above, the McKean-Vlasov equation (40) can be viewed as a single-particle linear FPE with the effective potential (41) that I repeat here for convenience:

provided is found from self-consistency equations that would be introduced below. However, prior to this, we can note that for the specific choice of the NES potential (23), the effective potential ( ) will have exactly the same functional form as the single-stock potential ( ), albeit with different parameters. Using the physics nomenclature, one can say that parameters of the original single-particle system are renormalized by interactions which are represented by the second term in Eq.(44).

To find the new renormalized parameters, we write Eq.(44) using Eq.(23) as follows:

where¯1,¯2,¯1,¯2,¯are new 'renormalized' parameters for a single-particle NES potential, andˆ( ) stands for terms that depend on but not on . They can be easily computed using the well-known relation expressing a product of two Gaussian densities as a rescaled third Gaussian density: 

Using these relations, we obtain the renormalized parameters and the functionˆ( ):

These formulae show that additional quadratic and linear terms in Eq.(44) can be re-absorbed into rescaled parameters of the Gaussian mixture, where the Gaussian means become linear functions of the mean field , while the mixing coefficient depends on non-linearly. This implies that within the mean field approximation, the dynamics of market log-returns can be modeled as the dynamics of a fictitious single stock, with initial single-stock model's parameters being 'dressed' by interactions according to Eq.(48).

For a confining potential ( ), the original finite dimensional Langevin equation (18) or (20) produces ergodic and reversible dynamics under the Gibbs measure

where is a normalization factor and ( 1 , . . . , ) is the potential (19), as a consequence of the linearity of the corresponding FPE, and uniqueness of the stationary state of this FPE. On the other hand, the McKean-Vlasov with a confining but non-convex self-interaction potential ( ) can lead to violations of ergodicity and phase transitions [2] .

To explore the phase structure in our setting with the NES potential (23), first note that using Eq.(24), the integral entering the self-consistency condition (43) can be expressed in terms of an integral that involves Ψ 2 0 :

Analysis of the phase structure of the model under variations of parameters is based on Eq.(50) where is viewed as an order parameter. This is similar to the classical Ising model where the mean magnetization is analogously used as an order parameter, and a second order phase transition occurs at a critical temperature at which the self-consistency equation = tanh( / ) of the Ising model bifurcates and produces a solution ≠ 0, in addition to the 'trivial' solution = 0 which describes a high-temperature regime → ∞.

A similar pattern of bifurcations and phase transitions can also be found, for certain types of the potential ( ), for the self-consistency equation (50) that deals with continuous random variables rather than binary spins of the Ising model. In particular, in the Desai-Zwangiz model, Eq.(50) is solved with the quartic double well potential ( ) = − 2 + 4 with , > 0. However, some general properties of the resulting stochastic system are determined by general properties of the potential, such as the number of local minima, and the asymptotic behavior at | | → ∞, rather than its particular functional form. In this work I will use the Log-GM NES potential (24) which produces a tractable double well potential for certain values of parameters, however the general analysis below holds for an arbitrary confining potential ( ).

If the potential ( ) is symmetric, (− ) = ( ), then = 0 is always a solution to the selfconsistency equation. This is similar to how the state with zero magnetization is the only stable solution for the Ising model in the high-temperature limit (which corresponds to the limit ℎ → ∞ in our conventions). However, for lower temperatures, the Ising spins becomes aligned with the mean magnetization = ±1, producing a bifurcation of the state equation in the parameter space. The bifurcation provides a mean field approximation description of the second order phase transition in the Ising model. Similarly in the current setting, for certain shapes of the potential ( ) and sufficiently low volatilities ℎ, the selfconsistency equation (50) can produce multiple solutions below some critical volatility ℎ . The analysis of the solution of this equation as a function of the 'temperature' parameter ℎ should be performed at different values of parameter that drives interactions in the system, as well as other parameters of the WF (24). The equilibrium second order phase transition describing the order-disorder transitions for the Desai-Zwanzig model was established in [2] . As will be shown next, a similar second-order phase transition can also occur for the MANES model when the potential ( ) is symmetric.

While for many choice of interesting non-convex potentials the analysis requires numerical integration that needs some extra care in the presence of multiple solutions (see [8] ), in our case the integral involved in the self-consistency relation (50) can be easily computed analytically for the WF Ψ 0 defined in Eq.(24), with its square given by Eq.(25). Using Eq.(46), we obtain

The partition function ( ) that enters this expression can be computed using Eqs. (44) and (45):

whereˆ( ) and parameters¯,¯1,¯2 etc. are defined in Eq.(48). For a symmetric NES potential with 1 = − 2 = , 1 = 2 = and = 1/2, this expression is further simplified:

whereˆ( ) stands for the symmetric version ofˆ( ) in Eqs. (48):

The derivative log / with the symmetric potential is therefore as follows:

Substituting this expression into (51), we obtain the MANES self-consistency equation for a symmetric potential:

= sinh 2

which looks similar to the self-consistency equation for the Ising model, see Fig. 3 . Note that a symmetric potential assumed in Eq.(56) may not necessarily match the market data, and actual self-interaction potentials ( ) implied by market prices are typically not symmetric (see Sect. 5). Nevertheless, analysis of symmetric potentials is of interest because it connects with the theory of second-order phase transitions. As shown below in Sect. 4.3, for a symmetric potential, the mean field vanishes if the volatility ℎ exceeds a certain critical value ℎ , and becomes non-zero, = ± 0 for some value 0 , for yet lower values ℎ < ℎ , with a continuous change from = 0 for ℎ > ℎ to non-vanishing values of for ℎ < ℎ . This describes a continuous (second-order) phase transition, similar to the one obtained for the Ising model with the vanishing magnetic field field ( = 0).

On the other hand, a non-symmetric self-interaction potential ( ) can be approximated by a symmetric potential ( ) ( ) with a fictitious 'magnetic field' 0 , see Eq.(32). Driven by the analogy with the Ising model, we could expect that this setting would produce scenarios for a first-order phase transition describing a decay of a metastable state, rather than a second-order phase transition. As we will see next, this is indeed the case in the present model.

When the potential ( ) is not symmetric, the self-consistent mean field should be computed using the general formula (51). Such analysis can be further simplified using a linear approximation to an asymmetric potential by adding a fictitious external field 0 to a symmetric potential ( ) ( ), see Eqs.(32) and (34). This is equivalent to replacing → + 1 0 and ( ) → ( ) ( ) in Eq.(44) that defines the effective potential ( ). Therefore, the generalization of Eq.(53) to the case of a slightly asymmetric potential ( ) can be obtained by the same replacement → + 1 0 :

Using this expression, we can obtain a generalization of Eq.(56) for the case 0 ≠ 0:

As implied by this equation, when 0 ≠ 0, the mean field does not vanish for any value of ℎ, and therefore there is no second-order phase transition when 0 ≠ 0. Instead, in this case we obtain a first-order phase transition under variations of the mean field 0 . When 0 is very small but still nonvanishing, its sign determines whether the negative or positive solution with = − 0 or = 0 will have the lowest energy, and thus will be the true ground state, where ± 0 are two degenerate mean field solutions obtained in the strict limit 0 = 0. The discrete Z 2 symmetry ↔ − is explicitly broken when 0 ≠ 0. If an initial state with e.g. = + 0 is observed at time = 0 and the 'magnetic field' 0 > 0, then it is the left-well state = − 0 that would be the true ground state, while = + 0 will be a metastable state. Vice versa, for 0 < 0, the state = 0 would be the true ground state, while = − 0 could be a metastable state released at = 0. A decay of such a metastable state is described as a first-order phase transition, which amounts to a sudden discontinuous jump from = − 0 to = 0 (or vice versa, depending on the sign of 0 and the initial state). An example of an asymmetric potential leading to such a first-order phase transition will be shown in the next section. The model thus suggests an interplay between the first-and second-order phase transitions in different regimes of parameters (ℎ, 0 ) which is similar to the phase transitions pattern in the Ising model.

Analysis of the phase structure of the model is performed using the traditional approach similar to the mean field analysis of the Ising model. We solve the self-consistency equation to compute a function = (ℎ), while keeping other model parameters fixed. The same exercise can be repeated for different values of another important model parameter that controls interactions between individual particles (assets).

The phase diagram obtained with this method for a symmetric potential assumed in Eq.(56) is shown in Fig. 4 . The order parameter obtained as a function of a continuously varying volatility parameter ℎ bifurcates at the critical value ℎ . The transition from a zero to a non-zero mean field at the bifurcation point ℎ = ℎ is continuous, as it should be for a second order phase transition. To identify stable versus unstable solutions of the self-consistency equation for either a symmetric or asymmetric potential ( ) (see, respectively, Eqs.(56) and (58)), we need to compute the Gibbs free energy F which is given by the following expression [8] :

For the stationary density (42) and the quadratic Curie-Weiss interaction potential, the free energy F = F ( ) can be computed in a more explicit form:

The gradient of the free energy with respect to is therefore as follows:

The stationary points of the free energy are obtained by setting this expression to zero, which again produces the self-consistency equation (51). While the relations (59) and (61) are general and apply for the generalized Desai-Zwanzig model with an arbitrary self-interaction potential, in this work that uses the log-GM NES potential (23), the free energy can be computed in closed form for either symmetric or asymmetric potentials ( ) using, respectively Eq.(53) or Eq.(52). To have even more compact formulae, it is convenient to use the partition function (57) corresponding to a weakly asymmetric potential whose asymmetry is controlled by a fictitious external field 0 . This produces the following expression:

where the ellipses stand for omitted constant terms that do not depend on , and parameter is defined as follows:

Clearly, the free energy F ( ) is not symmetric, i.e. F ( ) ≠ F (− ) for 0 ≠ 0 as F ( ) ∼ 0 + ( 2 ) for small values of . More generally, as can be seen from Eq.(60), the partition function ( ) and the free energy F ( ) are symmetric as long as Ψ 0 ( ) is symmetric in the -space, i.e. Ψ 0 ( )

For examples of shapes of the free energy leading to scenarios with phase transitions, see Fig. 5 . Note that for the left graph obtained with a symmetric potential with 1 = − 2 , 1 = 2 , = 0.5, the two minima at = ± with 0.4 are degenerate, and the point = 0 is a local maximum and is therefore unstable. This setting corresponds to the second-order phase transition for ℎ < ℎ , and spontaneous breaking of the Z 2 symmetry ↔ − of the free energy F ( ) for 0 = 0. On the other hand, on the right graph, the free energy is non-symmetric for a non-symmetric potential obtained with 

Critical behavior in the present model is defined in a similar way to the Ising model. For the latter, the critical behavior and bifurcation diagram is usually considered with the temperature being the control parameter and average magnetization being the order parameter, while the external magnetic field is used as an additional control parameter. Similarly, in the framework considered here, the volatility parameter ℎ serves in a similar way to the temperature parameter in the Ising model, while the coupling constant is used as an additional degree of freedom similar to the magnetic field in the Ising model. To investigate the critical behavior of the model in the vicinity of the phase transition at ℎ = ℎ , the free energy (62) is expanded into a Taylor series around = 0 to the fourth order in and second order in 0 :

where constant terms and higher order terms in 0 , are omitted, and parameters 1 , 2 , 4 are defined in terms of parameter introduced in Eq.(63) and other model parameters as follows:

The free energy F is written in Eq.(64) as a function of the noise variance ℎ 2 rather than of the mean field because in this section we want to explore its dependence on ℎ 2 . Note that Eq.(63) implies that 0 ≤ ≤ 1, therefore the coefficient in from of 4 is always positive, ensuring stability of any approximate solution that would be based on the small-expansion (64). On the other hand, one can see that the coefficient 2 can change the sign depending on the value of ℎ 2 . A bifurcation point ℎ 2 = ℎ 2 corresponds to the value of ℎ at which the coefficient in front of 2 in Eq.(64) vanishes, and then becomes negative for yet smaller values ℎ 2 < ℎ 2 . This produces the following relation for the critical volatility parameter ℎ :

To produce a real-valued parameter ℎ , the expression under the square root should be positive, producing a constraint on parameter combinations that may lead to bifurcation scenarios:

The coefficient 2 in the expansion (64) can therefore be written in a more suggestive form:

In the vicinity of the bifurcation point ℎ = ℎ , the solution for in the limit 0 → 0 can be well approximated by a solution to Eq.(64), which reads

This produces the following expression for the mean field in the vicinity of the critical point ℎ = ℎ :

which again looks very similar to the relation ∼ ( − ) 1/2 with the same critical exponent = 1/2 arising for the Ising model. 13 . As in the Ising model, the order parameter is continuous across the critical point ℎ 2 = ℎ 2 , indicating that we deal here with a second-order (continuous) phase transition. Substituting Eq.(69) back into Eq.(64) gives an approximate expression for the free energy for values of ℎ that are slightly below the critical value ℎ :

While this expression approximates the free energy F (ℎ) for 0 = 0 for values of ℎ that are approaching ℎ from below, for values of ℎ that are above ℎ , the value of F (ℎ) in Eq.(64) will be zero. Further following the analogy with the Ising model, we next define the 'specific heat' to be proportional to the second derivative of the free energy with respect to noise variance parameter ℎ 2 :

Because the expression in Eq.(71) arises only for ℎ 2 < ℎ 2 but vanishes for ℎ 2 > ℎ 2 , it is clear that the first derivative F (ℎ 2 )/ ℎ 2 is continuous at the critical point ℎ 2 = ℎ 2 , but its second derivative

has a finite jump at ℎ 2 = ℎ 2 , translating into a finite jump of the specific heat at the critical point:

This is again the behavior characterizing a second-order (continuous) phase transition which is similar to the second-order phase transition of the Ising model. Due to a finite jump, we obtain a vanishing value for the critical exponent entering the formula for the specific heat ∝ |ℎ − ℎ | , i.e. = 0.

To consider other properties of the model such as pairwise correlations, in this section we depart from the approximation of a homogeneous market portfolio that was used above, and consider a heterogeneous market where different stocks may have different parameters , . The partition function relevant for this setting is given by the following expression:

Here control parameters are introduced to facilitate calculations of various expectations and correlation functions, and are similar in their meaning to an external magnetic field used in the Ising model and other models of phase transitions for similar purposes.

With the Curie-Weiss quadratic interaction potential and MANES self-interaction potential, the partition function can be computed analytically in the limit → ∞. To this end, first note that the interaction term in the potential can be written as follows: 13 Here stands for the temperature, not the time interval as in the previous formulas.

where in the last step I replaced ( − 1)/ → 1 assuming that 1, and introduced matrix J with matrix coefficients = 2 ℎ 2 1 − . This can also be written in the matrix form as follows:

where I is a unit × matrix, and 1 is a vector of ones of size . The interaction term can now be represented using integration over auxiliary variables with = 1, . . . , using the Hubbard-Stratonovich transformation:

which holds for any real symmetric and invertible matrix J, whose inverse matrix has matrix elements −1 . The inverse of matrix J defined in Eq.(76) is computed using the Sherman-Morrison formula: 14

Using Eq.(77), the partition function (74) can be written as follows:

The last expression contains a product of one-dimensional integrals, which can be evaluated in close form by noting that terms proportional to and 2 in the exponential can be combined with the potential ( ) into a new effective potential similarly to Eqs.(44), (45):

Using Eq.(45) with replaced by , we obtain

where functionˆ( ) is defined as in Eq.(48), and

and parameters¯1,¯2,¯1,¯2,¯are defined as in Eq.(48) for each stock (I write them here as¯1( ) etc. to emphasizes their dependence on parameters ).

14 The Sherman-Morrison formula

holds for a non-singular matrix A and column vectors b, c such that the combination A + bc is non-singular.

Substituting (81) into Eq.(79) and changing for convenience the integration variables from to according to Eq.(80), the latter can be written as follows:

where

is the effective Hamiltonian for variables . The multi-dimensional integral with respect to can be well approximated in the large-limit → ∞ by a saddle point solution, i.e. a solution of the variational equation H ( , )/ = 0. The saddle point equation is therefore

The solution =¯of this equation hence satisfies the following equation

The partition function (83) with the saddle point approximation then reads

where¯is a solution to Eq.(86), and ( ) = H (¯, ) is the free energy. The local mean field = is defined as a partial derivative of :

This equation can be inverted to express variables¯in terms of local mean fields :

Using Eqs.(88) and (89), we can now write the saddle point equation (85) in terms of local mean fields :

The last equation is a general mean field self-consistency equation for the MANES model that defines the local mean field (i.e. the expectation of the log-return for the -th stock) in terms of the expected log-returns for other stocks. Other versions of the self-consistency equation can be obtained from Eq.(90) if we make further assumptions. For example, if we use it with a symmetric single-stock NES potential ( ), we obtain

For a homogeneous version of this self-consistency equation, we set → , and also remove indices from model parameters → , → . Furthermore, we have in this case → = −1 + 1 which is well approximated by = + 1 in the limit of large . This produces the self-consistency equation for the homogeneous portfolio:

This coincides with Eq.(58) provided we identify the external control field introduced in Eq.(74) with the fictitious field 0 introduced in Eq.(34) to mimic slightly asymmetric potentials.

Formulas for the local mean field approximation developed above enable computing susceptibilities for both the homogeneous and heterogeneous market portfolio settings. Starting with a homogeneous setting, the susceptibility defined in a similar way to the Ising model:

While in the Ising model such expression computes the sensitivity of the mean magnetization to changes of an external magnetic field, in the present setting, is the sensitivity of the expected market log-return to the amount of asymmetry in the the single-stock self-interaction potential.

To compute the susceptibility , assume a small but non-vanishing value of , and expand the right hand side of Eq.(92) to the first order in . After re-grouping terms, this gives

where the critical volatility parameter ℎ is defined in Eq.(66). This produces the following result for the susceptibility :

Therefore, at the critical volatility value ℎ 2 = ℎ 2 , diverges as (ℎ 2 − ℎ 2 ) −1 , similar to the Ising model behavior.

A similar analysis can be performed without using the homogeneous portfolio setting but rather working with Eq.(90) defined for the local mean fields (i.e. expected values) . Again, for small external fields → 0, we also expect the local fields to be small. In this regime, one can retain only the leading linear term in in the expansion of the second term in (91), and write

log Ω(¯)

where := ℎ 2 ℎ 2 + 2 ℎ 2 2 + 2 4 2

Plugging Eq.(96) back into Eq.(90), the latter can be written as a linear system of equations

where G is a matrix with matrix elements

The solution of (98) is

where A • B stands for a direct product of vectors A and B. This produces the following result for local susceptibilities

The inverse of matrix G can be found using the Sherman-Morrison formula:

In the limit of large , this can be well approximated by a simpler expression

The covariance of log-returns for two assets and is defined as follows:

This can be computed from the the partition function as follows:

where is the local mean field defined in Eq.(88). Using Eqs.(101) and (103), we obtain

The first relation here shows that in the mean field approach, we obtain ∝ implying that the covariance matrix C has rank one. This is similar to the relation ∝ that arises in a one-factor model with a common 'market' factor for all stocks, with being the regression coefficient of stock 's return on the market return.

A link with the calculations performed for the homogeneous setting is provided by considering a special case of the partition function (74) for a homogeneous external field which is the same for all stocks, i.e.

→ . For this case, the first derivative of log reads log = 2

Differentiating once more, one obtains

where Eq.(107) is used at the last step. Re-arranging the last equation here, we obtain the following relation between the "magnetic susceptibility" and the average covariance:

where is computed in Eq.(95). This relation shows that the susceptibility is driven by the fluctuations in the system. Such relations are known in statistical physics as fluctuation-response formulae. For a homogeneous setting with = , the mean covariance entering Eq.(109) can be computed as follows:¯= 1

Using here Eq.(97) with replaced by , we finally obtain = 1 2

where is the susceptibility for the homogeneous market computed in Eq.(95). We have therefore verified that covariances for a heterogeneous market defined in Eq.(106) reduce for a homogeneous market portfolio to the average covariance¯which is proportional to the susceptibility . In addition to establishing a correspondence with the homogeneous market setting, the local mean field approach of this section leads to the following observation. While the susceptibility diverges as (ℎ 2 − ℎ 2 ) −1 as ℎ → ℎ , this divergence originates in off-diagonal elements with ≠ , and arises due to the factor 1 − in the denominator in the first relation in Eq.(106). Therefore, while off-diagonal covariances are proportional to 1/ and are hence parametrically small in the limit of large , their values are increased as the volatility parameter ℎ approaches its critical value ℎ . This implies that, for certain combinations of parameters that produce scenarios with ℎ 2 ℎ 2 , covariances and/or correlations between different stocks obtained in this modeling framework can be made comparable with the average correlation between stocks in the real market, which is around 0.4 for stocks in the S&P 500 universe.

As was shown in Sect. 3.5, with the mean field approximation, the dynamics of the market index can be represented as a single-stock dynamics with renormalized parameters¯and¯,¯(with = 1, 2) given by Eq.(48) which I repeat here for relations involving¯,¯:

Therefore, the resulting effective single-stock model has six parameters¯1,¯2,¯1,¯2, , ℎ. In [10] , this model specification was further reduced to five parameters by setting¯1 = −¯2 =¯. In this paper, such constraint will not be imposed. Calibration of the model to market prices of S&P 500 options (SPX options) or other index options therefore produces six parameters¯1,¯2,¯1,¯2,¯, ℎ. On the other hand, the full set of model parameters in the original multi-asset version of the model involves seven parameters 1 , 2 , 1 , 2 , , ℎ, , plus one more unknown value of the expected log-return , which thus effectively serves as the eighth parameter. Clearly, as eight parameters cannot be uniquely recovered from six parameters¯1,¯2,¯1,¯2, , ℎ that could be found by calibration to SPX options, we need additional constraints to fix their values. The next few sub-sections develop such constraints on model parameters.

To estimate the coupling constant that controls interactions in the system, one can try to fix it by fitting the single-stock volatility and pair-wise return correlation obtained in our homogeneous portfolio setting to the average stock vol and correlations obtained in the real market. This can be readily done using the homogeneous version of Eqs.(106). It gives the mean single-stock volatilitȳ

where Eq.(97) was used at the last step. Note that the factor 1/ √ in the right-hand side of this relation arises because¯is defined in annualized terms. The mean pairwise correlation is obtained from Eqs.(106) as the ratio of the first row to the second row:

Note that¯∼ 1/ , which implies that correlations die off in the strict thermodynamic limit → ∞.

The mean correlation can also be estimated differently by computing the variance of the mean field in terms of and the mean single-stock volatility :

Neglecting the first term in the right hand side, this relation gives an approximation for in terms of the ratio of annualized variance of the market index' log-return 2 to the variance¯2 of the representative stock in the portfolio:¯ This suggests that should be very close to one, approaching it from below, and respectively ℎ 2 /(2¯2 ). Next, we can invert the second relation in Eq.(112) to obtain

This relation shows that should be such that ℎ 2¯2 < 1 in order to keep the single-stock volatility real-valued and finite, should satisfy the following constraint:

Using the second relation in (117), the constraint (119) can also be re-stated as the following constraint on the volatility of the 'representative' stock:

where the second approximate form follows as long as 1 as suggested by Eq. (117) . The second observation with Eq.(118) is that the single-stock variance 2 is higher than the variance of the mean log-return¯2 in both states = 1, 2. Equivalently, it means that the variance of the mean¯2 is smaller than the individual variance 2 -which is as expected from the central limit theorem.

If we use the values¯= 0.3 − 0.4 and = 500, then Eqs.(117) imply that should approach one from below, and ∼ ℎ 2 /(2¯2 ). On the other hand, using Eq.(97), we can write the expression for as follows:

Given that should be slightly below one, this implies that ℎ 2 should be slightly above the critical value ℎ 2 . This means that the model is in a high-temperature phase, yet close to the bifurcation point ℎ 2 = ℎ 2 .

In this regime, the non-vanishing mean field is solely due to the explicit breaking of Z 2 symmetry due to asymmetry of the potential, which can be mimicked by adding the fictitious field 0 , see Eq.(34). This behavior of the model appears to be quite reasonable. Indeed, a benign market regime is typically associated with a steadily growing market which, while exhibiting some volatility, has a positive trend. While it is not easy to estimate this trend accurately, a market with a positive market trend is certainly very different from a market with a negative trend. In a world with a symmetric potential and a second-order phase transition, a non-vanishing expected market return can only arise due to a spontaneous breaking of the Z 2 symmetry, which implies that both choices are 'physically' equivalent, i.e. equally preferable, in contradiction with the common sense. This suggests that the scenario in which the Z 2 symmetry is broken explicitly due to asymmetry of the potential, while the sign of the the mean field is fixed by the potential, appears more plausible in the present context. The next section shows how the mean field can be computed for this scenario using the MANES model calibrated to option prices.

As was discussed above, the equilibrium expected market log-return is given by the mean field parameter in the McKean-Vlasov equation (40), and it should satisfy the self-consistency equation (50). Here I present an explicit solution of this equation, assuming that the model is calibrated to option prices, so that renormalized parameters¯1,¯2,¯1,¯2, , ℎ are known.

Using Eq.(118) in Eq.(112), the relations between the 'bare' and 'renormalized' parameters, resp. and¯can be written as follows:

According to Eq.(26), the mean of the third Gaussian component in Eq.(25) is as follows:

The self-consistency condition (50) is equivalent to the constraint on the expected value of obtained in the model with the effective potential with the interaction-dressed parameters defined in Eqs.(122) and (123). Therefore, the self-consistency condition (50) now takes the following simple form:

where weights¯are defined as in Eq.(28) using the 'interaction-dressed' parameters¯1,¯2,¯1,¯2, ℎ,¯.

Using here Eq.(123) and re-grouping terms, we obtain

This relation provides a closed-form solution of the self-consistency condition (50) in terms of renormalized parameters that are directly calibrated to option prices. Recall that Eqs.(48) imply that renormalized weights¯depend on 'bare' weights and the mean field . Therefore, have we used the bare weights instead of renormalized weights¯and expressed¯1,¯2 in terms of 1 , 2 and according to Eq.(122), Eq.(125) would amount to a non-linear equation similar to the self-consistency equation (58). Instead, by working with the MANES effective potential (45) and renormalized parameters (122), the original self-consistency condition (43) of the McKean-Vlasov equation (40) is resolved here analytically rather than numerically. Interestingly, when expressed in terms of renormalized parameters in Eq.(125), the expected market log-return does not explicitly depend on the coupling constant , and is found in terms of parameters directly calibrated to option prices.

On the other hand, the 'bare' model parameters 1 , 2 etc. do depend on the value of . In particular, after the mean field is computed from Eq.(125), the 'bare' model parameters 1 , 2 can be obtained from Eq.(122), volatilities 1 , 2 are obtained using Eq.(118), and the bare mixing coefficient can be computed by inverting the formula for¯in Eq.(48):

In this section, I present three sets of examples of calibration to market quotes on European options on SPX (the S&P 500 index). I use the same set of option quotes that were used in [10] to illustrate the working of the NES model. In all three sets of experiments, I calibrate to market quotes on 10 call options and 10 put options. The strikes are chosen among available market quotes to cover the the range of option deltas between 0.02 and 0.5, in absolute terms, so that the model is calibrated to both ATM strikes and deep OTM strikes. Details of the loss function are described in [10] . Optimization of the loss function is done using the shgo algorithm available in the Python scientific computing package scipy.

In all experiments presented below, parameters¯1,¯2,¯1,¯2,¯, ℎ are found by calibration to SPX options. Inferred parameters are different from those reported in [10] , because here I do not enforce the constraint¯2 = −¯1 that was used in [10] . In addition, I display estimated parameters and . The coupling constant is roughly estimated according to the following formula that is obtained by combining Eqs. (117) and (120):

where Δ 2 is a margin that controls the strength of the inequality (120). For numerical examples, the value Δ 2 = 0.05 will be used in the examples below. Furthermore, the equilibrium log-return (mean field) is computed according to Eq.(125). In the first example, I consider SPX 1M options on 07/12/2021 with maturity on 08/09/2021. The results of calibration are presented in Table 1 and Figs. 6 and 7. Potentials shown in these figures and in the examples to follow are effective potentials (45), and are computed as a single-stock NES potential with parameters¯1,¯2,¯1,¯2,¯, ℎ inferred from calibration to index options. Note the difference in implied potentials for puts and calls, as well as different values of inferred model parameters. For both potentials, the current value of the log-return 0 , as shown by the vertical red lines, is located near the bottom of the potential well. Furthermore, the equilibrium log-return is also located near the bottom of the potential for for potentials. This suggests that the price dynamics on this date correspond to an equilibrium regime of small fluctuations around a stable minimum.

In the second example, I look at longer option tenors, and consider SPX 1Y options on 11/06/2020 with maturity on 09/21/2021. The results of calibration are presented in Table 2 and Figs. 8 and 9. Again, we can note the difference in implied effective potentials for puts and calls, as well as different values of inferred model parameters. Also as in the previous example, for both potentials, the current value of the log-return 0 and equilibrium expected log-return are located near the bottom of the potential wells, suggesting an equilibrium regime of small price fluctuations on this date. Finally, the last example considers 6M options with the expiry on 09/18/2020 on 03/16/2020, at the peak of Covid-19 crisis where the SPX index had the largest drop. This example thus illustrates the model behavior for a severely distressed market. The results of calibration are presented in Table 3 and Figs. 10 and 11. As in the previous examples, again note the difference in implied effective potentials for puts and calls, as well as different values of inferred model parameters.

Unlike the previous examples, for the present case of a distressed market, the initial log-return 0 on 03/16/2020 is located far away from the global minima for both effective potentials, indicating that this initial state is strongly non-equilibrium, consistently with a prevailing market sentiment on that date.

Interestingly, while both implied effective potentials suggest that the current value 0 is far from the global minimum of the potential and hence describes a non-equilibrium scenario, they differ in thē character of a subsequent relaxation mechanism for this initial state. The implied potential for puts in Fig. 10 is a single-well potential, therefore the initial state 0 is unstable. If the potential itself remains constant through time, this initial state would eventually relax into its global minimum. This would happen even in the limit of zero volatility (zero noise). When the noise is present with ℎ > 0, it produces an uncertainty for the time needed to reach the minimum of the potential, as well as small fluctuations around this minimum. Differently from the puts, the call options seem to suggest a different type of relaxation dynamics, as the implied potential for the calls in Fig. 11 is a double well potential quite similar to the one shown in the left column in Fig. 2 . The initial location is in the vicinity of the local minimum corresponding to the left well, while the right well corresponds to the global minimum. As they are separated by the barrier, the potential in Fig. 11 describes a scenario of metastability. The relaxation to the global minimum (the true ground state) proceeds via instanton transitions as discussed in [9, 10] . Instanton transitions are only possibly when the volatility parameter is non-zero, ℎ > 0, however small it can be in practice. This suggests that double well potentials implying metastable dynamics may occur when markets are in distress or during periods of a high market uncertainty, e.g. during general crises such as the 2008 crisis or the Covid-19 crisis of 2020, or during general elections. 15 

This paper proposed a tractable non-linear model of interacting and non-equilibrium market, formulated as statistical mechanics of interacting non-linear oscillators where individual stocks' log-returns are viewed as coordinates of 'particles' describing these stocks. Both self-interactions of oscillators and their pairwise interactions are explained in terms of money flows into the market from external investors. In particular, correlations between log-returns of individual stocks originate from the dependence of money flow into stock on the average previous performance of all other stocks.

The main idea of this paper was to approximate the return of a market index in a heterogeneous market by a mean field of a homogeneous market made of replicas of the same 'representative' stock with a self-interaction potential ( ). This enables modeling the dynamics of the market log-return using the mean field approximation that produces the McKean-Vlasov equation as the equation governing the dynamics of the system. As the McKean-Vlasov equation is a non-linear equation corresponding to the thermodynamic limit → ∞, it produces far richer dynamics than the original multi-particle Fokker-Planck equation (35), giving rise to ergodicity breaking and first-and second-order phase transitions in different parameter regimes. The resulting dynamics of the mean field resembles the Desai-Zwanzig model of interacting non-linear oscillators [3, 2] and its generalized version in [8] .

Furthermore, while the exploration of the phase structure in the generalized Desai-Zwanzig model requires dedicated numerical methods to identify stable and unstable solutions [8] , the analysis in the present paper is considerably simplified due to its reliance on the Non-Equilibrium Skew (NES) potential (24). The NES model with the NES potential (24) was introduced in [10] as a flexible and highly tractable non-linear model of single-stock dynamics that is capable of describing either a benign or stressed market regime, and tracks the pre-asymptotic dynamic behavior of transition probabilities. In particular, non-equilibrium, pre-asymptotic corrections to the moments of a stationary (asymptotic) distribution of log-returns are explicitly controlled in the NES model in terms of the model parameters entering the NES potential (24).

In the present paper, the NES potential from [10] is used to model a multi-asset market. In Sect. 3.2, the NES potential is introduced as a tractable approximation to a non-linear self-interaction potential ( ) that originates in money flows and their impact effect on market returns. In addition, the quadratic Curie-Weiss interaction potential (21) is obtained as an approximate interaction potential following the same lines of analysis. This results in a multi-asset extension of the previous single-stock NES model, referred to as the Multi-Asset NES (MANES) model in this paper.

As was shown in this paper, the new multi-asset MANES model is as tractable as the previous singlestock NES model. This is made possible due to the fact that with the NES self-interaction potential ( ), the new effective potential ( ) that incorporates interactions in the system can be expressed in terms of a single-stock NES potential ( ) with renormalized parameters, see Eqs.(45) and (48). This suggests that for the purpose of calibration to market prices of index options, the multi-asset MANES model is computationally equivalent to a single-stock NES model applied to the 'representative' stock. Also due to this property of the model, the self-consistency equation of the mean field approximation that is normally expressed as a non-linear equation (43) is resolved here analytically in terms of renormalized parameters that are directly calibrated to market prices of index options. Furthermore, estimates made in Sect. 5.2 suggest that the model operates in a 'high-temperature' phase close to the criticality point, where the volatility parameter ℎ is slightly higher than the critical value ℎ = ℎ that corresponds to the bifurcation point of the order parameter .

Just as the single-stock NES model, the MANES model demonstrates that a single volatility parameter is sufficient to accurately match available market prices of index options. These results stand in stark contrast to the most of other option pricing models such as local, stochastic, or rough volatility models that need more complex specifications of noise to fit the market data. Alongside a single volatility parameter ℎ, the effective single-stock potential ( ) (45) has parameters¯1,¯2,¯1,¯2,¯that can be calibrated to market prices of index options. By fitting these parameters to the market data, we produce implied potentials which replace implied volatility smiles as a way to fit market data. If desired, the calibrated model parameters¯1,¯2,¯1,¯2,¯can be approximately mapped onto the parameters that enter Eq.(22), and thus could be interpreted as as implied money flow and market impact parameters. Clearly, the success of this approach for matching market prices of vanilla index options does not preclude one from using more sophisticated models of noise such as stochastic or rough volatility models -but only if needed beyond the need to explain market prices of vanilla options.

The MANES model can be used for several applications or practical interest. In particular, optionimplied moments of future returns can be used as predictors for actual future returns, volatilities and skewness, and employed for portfolio trading, see e.g. [7, 17] . While such analyses typically use riskneutral moments implied by option prices, the MANES model enables extracting both risk-neutral and real-measure moments, and thus enriches the set of predictors for such tasks. The model can also use other market data for model calibration. In particular, in addition to using market prices of options, we could incorporate open interest data in the model calibration. The MANES model could also be jointly calibrated to the equity and credit markets data, by adding a fit to credit indices such as CDX as proxies to probabilities of large market drops. Such applications and extensions are left here for a future research.

The Inelastic Market Hypothesis: A Microstructural Interpretation

Critical Dynamics and Fluctuations for a Mean-Field Model of Cooperative Behavior

Statistical Mehanics of a Nonlinear Stochastic Model

Nonlinear Fokker-Planck Equations

Dumb Money: Mutual Fund Flows and the Cross-Section of Stock Returns

Search of the Origins of Financial Fluctuations: The Inelastic Markets Hypothesis

How Useful are Implied Distributions? Evidence from

Dynamics of the Desai-Zwanzig Model in Multi-Well and Random Energy Landscapes

Quantum Equilibrium-Disequilibrium: Asset Price Dynamics, Symmetry Breaking, and Defaults as Dissipative Instantons

Non-Equilibrium Skewness, Market Crises, and Option Pricing: Non-Linear Langevin Model of Markets with Supersymmetry

The Inverted World of Classical Quantitative Finance: a Non-Equilibrium and Non-Perturbative Finance Perspective

Quantum Mechanics

Sur la Théorie du Mouvement Brownien

Mean-Field Treatment of the Many-Body Fokker-Planck Equation

A Class of Markov Processes Associated with Nonlinear Parabolic Equations

Statistical Mechanics, Harper & Row

What Does Risk-Neutral Skewness Tell Us About Future Stock Returns?

Quantum Field Theory and Critical Phenomena, Fourth Edition