Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
119,687 characters · 21 sections · 74 citation commands
Time is limited on the road to asymptopia
\noindentKeywords\ Broken Ergodicity, Simulated Method of Moments, Validation, Calibration, Agent-based Models
\noindentJEL Classification\ \href{https://www.aeaweb.org/econlit/jelCodes.php?view=jel#C}{ C61} $\cdot$ \href{https://www.aeaweb.org/econlit/jelCodes.php?view=jel#D}{ D01 $\cdot$ D81 $\cdot$ D9}
Financial crashes, global pandemics or disruptive innovations all qualify as formative events during which economic systems reach a new order. On an aggregate level such events can trigger the transition from an old equilibrium to a new equilibrium. A rapid sequence of such formative events can prevent the economy, financial markets and observables from ever reaching their long-time equilibria and thus keeping them steadily out-of-equilibrium. In an effort to use a new methodology to understand such complex emergent behaviors and answer urging policy questionsBouchaud2008,FarmerFoley2009,LuxWesterhoff2009 a community of researchers adopted financial agent-based models (FABMs) that deliberately allow for more complexity than traditional models.TesfatsionJudd2006,HommesLeBaron2018 The scientific challenge is to work with FABMs that are capable to endogenously reproduce known statistical regularities and that remain statistically and computationally tractable for estimation exercises.\footnote{The most important stylized facts of financial time series include phenomena like excess volatility, fat-tailed return distributions, absence of autocorrelations in the raw returns but long memory in absolute returns and volatility clustering. Regarding the study of financial markets, there exists a large strand of literature which provides models that explain the stylized facts well. Surveys about the statistical properties of financial markets can be found in Cont2001,LuxAusloos2002,Hommes2006,LeBaron2006,Lux2009.} The ever increasing amounts of data and computational resources make FABMs a promising candidate to ultimately join the set of models which drive informed policy decisions.WesterhoffFranke2018
The development of a consensual validation and estimation protocol for (F)ABMs is important for scientific standards of reproducibility. However, the estimation and validation of FABMs remains challenging and despite a growing body of research there exists no consensus in the literature, yet.\footnote{RichiardiEtAl2006 is one attempt in this direction, DelliGattiEtAl2018 is a helpful textbook. See FagioloEtAl2019,LuxZwinkels2018 for recent surveys.} Having a standardized protocol after all would address main methodological critiques raised against agent-based modeling.LeombruniRichiardi2005,RichiardiEtAl2006,FagioloRoventini2017
Here we follow a direction outlined in GaffeoEtAl2007 of a \textquote{descriptive output validation}, where we focus on simulated moment estimators which have been studied already intensively for FABMsAmilon2008,GilliWinker2003,WinkerEtAl2007,Franke2009,Schmitt2021,Schwartz2022. Apart from the simulated method of moments (SMM) used here, different estimation approaches have recently emerged in the literature including likelihood-based methods, information-theoretic similarity-based measures and more recently a number of Bayesian estimations of FABMs appeared, too.KukackaBarunik2017,Lamperti2018,Barde2016,GrazziniEtAl2017,Platt2021,BertschingerMozzhorin2021,DelliGattiGrazzini2020,Lux2018,Lux2021,Lux2021a
Reliable estimation rests on sufficient convergence speed given finite-sample lengths and limited computational resources. Surprisingly, the development of better tools and protocols in the validation of (F)ABM has often overlooked the ergodicity issue so far. Our contribution enters exactly at this issue. Already the seminal contribution by Cont2001 emphasised ergodicity in the process of estimation of parameters in models of financial market: \blockquote[emphasis added]{One needs an ergodic property which ensures that the time average of a quantity converges to its expectation. Ergodicity is typically satisfied by IID observations but it is not obvious -- in fact it may be very hard to prove or disprove -- for processes with complicated dependence properties such as the ones observed in asset returns.} Admittedly, ergodicity is an elusive property and the fact that no ready-to-use econometric test for ergodicity exists stands in the way of a comprehensive treatment. It may come as less of a surprise then that only few publications identify ergodicity as an important assumption at all.\footnote{To convey a rough idea, the four volumes of the \citetitle{TesfatsionJudd2006} only casually mention 'ergodic'; esp. in the volumes most relevant, i.e.\xspace Vol. 2 on \citefield{TesfatsionJudd2006}{subtitle} TesfatsionJudd2006 and Vol. 4 on \citefield{HommesLeBaron2018}{subtitle} HommesLeBaron2018, 'ergodic' is mentioned in Ch. 22 & 32 and resp. in Ch. 8 & 14, however with no relation to the validation of ABMs.} Regardless, ergodicity often lurks in the background as an implicit regularity condition.\autocites[19]{SoltykChan2021}[p. 25, footnote 36]{FrankeWesterhoff2016}[p. 91, footnote 8]{LeBaron2021} As mentioned by Cont2001 a proof of ergodicity of the relevant observables is not trivial, because it depends on the convergence speed and the respective time scales are often large. The danger for modellers is to run into unnoticed trouble during the process of estimation and validation since knowledge about the estimator's convergence behaviour is instrumental.
In this paper, we aim to contribute to a growing body of research on the validation and estimation of FABMs. In particular, we add useful steps towards a coherent study of the effects of broken ergodicity. We do not propose any new estimation technique, but rather enhance the understanding and applicability of existing ones with regard to the ergodicity issue. Our main contributions are the following. First, there exists a nucleus of publications studying ergodicity in the estimation of (F)ABMs especially by Jakob Grazzini and Matteo Richiardi, which started with the study of models much simpler than the FABMs we study hereGrazzini2011,GrazziniRichiardi2013,GrazziniRichiardi2015. We analyse broken ergodicity in the SMM estimation of two well established FABMs, the AlfaranoEtAl2008 (ALW model\xspace) and the FrankeWesterhoff2012 (FW model\xspace), which reproduce a much larger set of stylized facts of financial markets. We perform Monte Carlo simulations and use the simulated method of moments (SMM) approach for the estimation of both ABMs.\footnote{Standard references for the SMM are McFadden1989,PakesPollard1989,LeeIngram1991,DuffieSingleton1993. A textbook treatment is in DavidsonMacKinnon2003.}
Second, we use the terminology of time averages and ensemble averages which allows to express many challenges in the estimation process as depending on taking different types of averages. For instance, we explicitly analyse the convergence behaviour over time vis-a-vis the convergence behaviour over replications. Doing so, we show how to reduce the uncertainty coming from broken ergodicity by analysing the estimation for different choices of ensemble size and simulation time length under limited computational resources, for which most estimators remain in a pre-asymptotic regime.
A brief survey of related work puts our contribution in context of the existing literature. Our work is closely related to a series of publications initiated by Grazzini and Richiardi who in Grazzini2012,GrazziniRichiardi2015 introduced a useful terminology and also use the SMM approach. However, they analyse the ergodic properties of much simpler models than those FABMs that we study here. As stated by GrazziniRichiardi2015, ergodicity assures consistent estimates even in the adjustment phase, i.e.\xspace when the system has not yet fully reached the limiting distribution, because the properties are invariant across different realizations of the stochastic process. This assumes that certain regularities are \textquote{stable enough across different replications of the model} to be identified. In general, Grazzini2012,GrazziniRichiardi2015 use non-parametric statistical tests to test for ergodicity (and stationarity). Grazzini2012 applies a Wald-Wolfowitz runs test to check whether the observed and simulated data come from the same population. This can also be done using a Kolmogorov-Smirnov test.\footnote{For ergodic Markov processes, the tests by DomowitzElGamal1993,DomowitzElGamal2001 are applicable.} These tests can offer a first starting point. However, DosiEtAl2018 find conflicting results using these econometric tests and thus conclude that they are not reliable to detect (non-)ergodicity. Another approach based on robust statistics is given by WinkerEtAl2007 who study situations when a typical distribution or time series is better represented by the median instead of the mean. The divergence of the median from the mean is a typical symptom of broken ergodicity as the mean or ensemble average is prone to outliers and because of that tracks the exceptional individual whereas the median is resistant to extremes and approximates the time-averaged behaviour and thus the behavior of typical individuals.
Furthermore, sufficient asymptotic convergence might be beyond reach for all practical purposes. BouchaudFarmer2021 recently termed these situations quasi non-ergodicity, i.e.\xspace convergence times are astronomically-long such that for most practical purposes the estimation takes place in a pre-asymptotic regime where ergodicity is notoriously hard to achieve. Similarly DessertaineEtAl2022 draw another analogy from physics and refer to non-ergodic effects of estimation in pre-asymptopia as violations of the adiabatic assumption, \textquote{in the sense that the time needed for the system to reach equilibrium is much shorter than the time over which the environment changes, so out-of-equilibrium effects can be neglected}. The model data-generating process (mDGP) of a (F)ABM is usually not analytically solvable let alone available in closed-form, which makes it even harder to choose the number of simulation runs and length without any prior knowledge about the speed of convergence. Our study helps to reduce this uncertainty.
The remainder of the paper is organised as follows. In Section (ref) we provide a basic understanding of ergodicity for the purpose of moment selection. It turns out that the simulated method of moments like any estimation approach closely relies on ergodicity somewhere along the process. Section (ref) shows where exactly ergodicity enters. In Section (ref) we investigate how (broken) ergodicity influences the validation of FABMs. In doing so, we run a bunch of numerical experiments on two financial market models and present the results of our simulations. Section (ref) contains a summary and concludes the paper.
This section starts with a brief discussion of ergodicity for the purpose to guide the selection of proper moments. In short, proper moments must be ergodic to be invariants of the model and not artefacts of particular simulated trajectories. With the help of an example of (geometric) Brownian motion we show that it is possible and necessary to derive ergodic observables from non-ergodic processes. These insights about ergodicity then guide the process of moment selection in Subsec. (ref). The moments enter the method of simulated moments which is briefly described in Subsec. (ref).
Ergodic theory can be understood as a mathematical theory that studies the asymptotic behaviour of averages. Results from this mathematical field become relevant in an economics context whenever an economic observable is modelled as a stochastic process, In our context of financial markets, the asset price as our economic observable of interest is modelled as a random variable in a static context or as a stochastic process if we are interested in the time evolution as well. Figure (ref) provides an intuitive understanding of the core aspects of ergodicity and shows five sample time series.
From a modelling standpoint we think of the five time series as five realizations of a stochastic process $\left\{ X_t \right\}_{t\in\ensuremath{\mathbb{T}}\xspace}$ with five different initial conditions $x_1, \ldots, x_5$. \footnote{For simplicity we will omit the initial conditions in the notation if they are not relevant.} For FABMs, the raw observables will be time series of asset prices or returns initialized by the series' first observations. We think of the observed data as being generated by some unknown real-world data-generating stochastic process (rwDGP), which remains a purely mental construction, of course. The FABM is the model data-generating process (mDGP), from which we are able to sample. Validation and estimation in econometrics can then be understood as a test of the match between the mDGP and the rwDGP with respect to some metric. In SMM, this metric consists of a vector of moments. The selection of moments is discussed in Subsec. (ref).
With the help of ergodicity and Fig. (ref), we distinguish between two types of the physical identity of an average, ensemble average (blue shade) and time average (red shade). Ensemble and time averages perform different aggregations on different instances of the observables. A time average computes an average along a single time series (or a single realization of a stochastic process) which is is depicted in Fig. (ref) by the red shade tracking the time series $X_t(x_3)$. An ensemble average computes an average over different realizations of the random variable (or stochastic process) at a single point in time which is depicted in Fig. (ref) by the blue shade covering realizations at time $t^*$. In econometric terminology time averages correspond to longitudinal averages and ensemble averages to cross-sectional averages. Intuitively, a time average gives an aggregate statement about a single system over time, an ensemble average gives an aggregate statement about many systems at an instance in time. \footnote{A key insight that has recently become better understood is to distinguish not only between the physical identity of an average (ensemble/time) but to further abstract and distinguish the physical identity from the mathematical identity or mathematical functional form of an average (arithmetic/geometric/harmonic/etc.\xspace) Kirstein2019.}
Ergodicity is then a mathematical property of an observable that describes the conditions under which both averages exist and coincide. A common definition of ergodicity establishes some type of convergence of longitudinal or time series averages on the LHS towards the cross-sectional average on the RHS:
where $f$ is an arbitrary function of the observables. This handy notation allows to cover the case of analysing the raw observables as well as some transformation of them, in the former case $f$ is simply the identity mapping, $f(x) = x$. In our context, the function $f$ potentially represents some summary statistic of the data, i.e.\xspace a moment function which we discuss below in more detail. The LHS of Eq. ((ref)) yields a time average in the long-time limit. It is important to note that the RHS is also the result of a limit, it yields the ensemble average in the large-ensemble limit. The latter is rarely made explicit and instead the RHS is often presented just as the expectation value, i.e.\xspace \ensuremath{\left\langle X \right\rangle}\xspace = \ensuremath{\operatorname{E}\left[ X \right]}\xspace.
In this section we briefly discuss the possibility to extract ergodic observables from non-ergodic processes. This corresponds to extracting moment functions from some rwDGP of which we don't need to know whether it is ergodic or not as long as we select an ergodic moment function. The goal is to identify invariants of a model which is a key goal in science.Petersen1996 Invariants are adequate quantities when describing and comparing models of real-world phenomena. Identifying these invariants and package them in moment functions will be instrumental for the selection of proper moments discussed in the subsequent section.
Ergodicity as a property of mathematical models characterises such possible invariants in the specific sense of Eq. ((ref)), i.e.\xspace an ergodic observable converges (almost always) towards a unique invariant limit distribution.\footnote{In a sense, an ergodic observable ensures that more data is better for estimation and validation purposes, because more observations are likely to lead to better convergence.} Most models of real-world phenomena are per se\xspace non-ergodic processes, because some dynamic aspect (like growth) is relevant.\footnote{Let us emphasise that ergodicity is not a property of reality but of the mathematical model.} For illustration purposes let us look at the standard model of asset price evolution which is that of (geometric) Brownian motion (GBM).
As such the model of GBM is clearly non-ergodic. Yet it is possible to derive ergodic observables from this non-ergodic process. Quite generally, a stochastic process is a model of an uncertain economic observable that can vary in time and over realizations. A simple stochastic process is the Wiener process $W_t$. In finance the Wiener process is used as a model of the price evolution, i.e.\xspace the repeated measurements of the price of an asset. \footnote{The Wiener process is sometimes also called Brownian motion because it was initially conceived to model the real-world phenomenon of Brownian motion, i.e.\xspace the jittery motion of pollen suspended in fluids.} The Wiener process $W_t$ is non-ergodic, as its time average behaves as $\ensuremath{\overline{W_t}}\xspace \sim \ensuremath{\mathcal{N}\left(0,\nicefrac{t}{3}\right)}\xspace$ for $t\to\infty$, and thus does not converge, because the distribution has a time-dependent variance, i.e.\xspace the distribution is broadening over time. Therefore, the time-dependent time average cannot almost always be equal to its time-independent ensemble average $\ensuremath{\left\langle W_t \right\rangle}\xspace = 0$, if so then only occasionally (by chance).
Again, we are interested in identifying ergodic observables or summary statistics of a process. In the unlikely case that the raw observable is ergodic we are lucky. More generally some transformations might be necessary to arrive at an ergodic observable. Typical derived observables are increments, squared increments (sometimes referred to as squared displacements), correlations, multiples, returns or growth rates. If we analyse the squared displacement from the origin of a Wiener process as our observable of interest, $D_T^2 = \left(W_T - W_0 \right)^2$, then this observable is ergodic, $\ensuremath{\left\langle D_T^2 \right\rangle}\xspace = \ensuremath{\overline{D_T^2}}\xspace$.MetzlerEtAl2014a The increments of a Brownian motion converge towards a unique invariant time-independent distribution, and are in this sense an ergodic observable. The increments of a geometric Brownian motion, however, are not converging but remain level-dependent. For GBM it is the exponential growth rate that converges to an ergodic invariant distribution instead. This little example is meant to show that it is possible to derive ergodic observables from non-ergodic processes. In fact, we use the ergodic observables to characterize non-ergodic processes and thus in a similar way like in this little example use ergodic moments to characterize the mDGP of a (F)ABM.
It is the ergodic property that makes a summary statistic a suitable moment for the validation of (F)ABMs. However, the high degrees of freedom in (F)ABMs often prevent the existence of a closed-form solution, which would simplify the direct analysis of the stochastic data-generating process and to analytically derive the limit distributions of observables. Therefore, we have to rely on numerical simulations to study the convergence behavior and ergodic properties of moment functions. With this in mind, we now turn to the process of moment selection.
The term moment is central in all moment-based estimations as the generalized method of moments (GMM) and the simulated method of moments (SMM) which we use in this paper. However, the word moment has different meanings in different fields. Let us therefore briefly clarify how we use it here. Whenever we refer to moments, we use this term in the narrow sense of their ordinary meaning in stochastics as properties of a probability distribution \ensuremath{\mathcal{P}}\xspace in Eq. ((ref)). In the context of SMM estimation of (F)ABMs, the literature refers to moments in a wider sense as some summary statistic. To avoid confusion we will thus refer to them as moment functions in this paper.
In stochastics a moment \ensuremath{m}\xspace of order $\ensuremath{q}\xspace = 0,1,2,3, \ldots ~; \ensuremath{q}\xspace\in\ensuremath{\mathbb{N}}\xspace$ of a distribution is computed as the mean of the \ensuremath{q}\xspace-th power of its realizations, $\ensuremath{m}\xspace_\ensuremath{q}\xspace = \int_{-\infty}^{\infty} x^\ensuremath{q}\xspace \ensuremath{P}\xspace(x) \ensuremath{\dif x}\xspace = \ensuremath{\operatorname{E}\left[ X^\ensuremath{q}\xspace \right]}\xspace.$ It is often convenient and therefore desirable to work with probability density functions such as \ensuremath{P}\xspace, but the density of a measure doesn't always exist, thus we can express moments using only the probability measure \ensuremath{\mathcal{P}}\xspace instead,
As the mDGP for most FABMs is not readily available in closed form, the moment functions can't simply be calculated analytically. Thus we need to simulate the data first and then compute the moments numerically. Inference from a given sequence of moments to a distribution is a classical problem in mathematics known as the moment problem. The estimation task of a (F)ABM is an inverse form of the classical moment problem. Generally, moment problems are only solvable if long sequences of moments are available, which is usually not the case for distributions of observables of financial markets because they turn out to be heavy-tailed. \footnote{A classic text on the moment problem is Akhiezer1965, more recent research reveals the solvability of the Hamburger moment problem depends on the non-negative extendability of the Hankel matrix, see Bolotnikov1996,Bolotnikov2008,ChenHu1998,DyukarevEtAl2008.} Higher-order moments of heavy-tailed distributions (HTDs) quickly vanish and in extreme cases even the first moment doesn't exist.Rachev2003 It seems therefore unlikely that solutions to classical moment problems are of immediate help for our estimation task. However, it provides an additional motivation to extend the selection of moment functions in SMM beyond the ordinary moments. \footnote{SoltykChan2021 use solutions to classical moment problems for modelling time-varying higher-order conditional moments of financial time series.}
Empirical analyses of economic and financial data at least since the 1960s\footnote{Mandelbrot1960,Mandelbrot1963,Fama1963 are prominent references that find heavy-tailed price increments (returns).} have repeatedly confirmed that many observables are not very well described by the family of Gaussian distributions but instead follow HTDs, often with a power law tail. Since the increasing magnitude of an extreme event raised to the \ensuremath{q}\xspace-th power, $x^\ensuremath{q}\xspace$, can't be offset by its decreasing probability, $\ensuremath{\dif \ensuremath{\mathcal{P}}\xspace}\xspace(x)$, the integral in Eq. ((ref)) diverges and higher-order moments cease to exist for HTDs. For instance, the calculation of the second central moment, the variance, is numerically always possible but does not converge for larger samples. Associated with the (existence of the) variance is a typical scale of dispersion. If the distribution has no second moment, then such the variations are so extreme that they lack a typical scale. Obviously, non-existing moments would not be good candidates for moment functions. An invariant of such volatile behavior might be the scaling exponent of the increments or returns. In fact, it is now understood that power laws are invariants of scaling relationships which again are imprints of complex systems and their otherwise out-of-equilibrium dynamics.
Moments \textquote[{Greene2018}]{provide a natural source of information about the parameters, other functions of the data may also be useful}. A natural motivation to move beyond moments to moment functions is to not only use static properties of distributions but also dynamic properties of the underlying processes, e.g.\xspace correlations are not a property of a distribution but of observations at different times, locations or both. We refer to these 'other functions of the data' as moment functions. Our notation proves convenient here, as it is possible to interpret the function $f$ in Eq. ((ref)) simply as some moment function $f(x) = m(x)$. Eventually, our vector of moments \ensuremath{ \ensuremath{\mathbold{m}} }\xspace can contain ordinary moments in the sense of Eq. ((ref)) as well as moment functions.
\paragraph{A hierarchy of convergences in ergodic theorems} In order to emphasise how strong a statement about an ergodic property actually is, it is worthwhile to make explicit a hierarchy of convergences contained in ergodic theorems like Eq. ((ref)). For an ergodic theorem to hold:
Let us comment on this hierarchy. On the first two points, most raw time series do not converge at all to a single value and are thus not stationary in this sense, e.g.\xspace all growth processes, security prices, business cycle dynamics or exchange rates. Some transformation like first or higher differencing is necessary to impose stationarity and thus convergence.DelliGattiEtAl2018 However, such transformations do not necessarily impose ergodicity, as they operate only on a single time series. On the third point, the expectation value is a particular average, i.e.\xspace the probability-weighted arithmetic mean of all realizations of a random variable. Thus its mathematical identity is an arithmetic mean and its physical identity is an ensemble average. Note that in general both time and ensemble averages do not coincide and hence non-ergodicity or broken ergodicity is our default condition. In fact, ergodicity only holds for very special cases or carefully chosen transformations of the observables.
This convergence hierarchy also applies to moment functions. Already \textcites[805]{Franke2009}[918]{RugeMurcia2012} note that many empirical moment functions are computed as time averages of some function of the (time series) data, $\nicefrac1T \sum_{t=1}^{\ensuremath{T}\xspace} \ensuremath{m}\xspace \left( x_t \right).$ We only take the obvious next step to demand the ergodicity of moment functions in the selection process such that they converge towards a unique value or distribution and capture invariants of the mDGP. So far it is only tacitly assumed in the literature that moment functions are ergodic but rarely explicitly analyzed, which implies that the calculation for (almost) all time series samples converges (almost) always towards a particular and identical value (namely the expectation value). For FABMs that offer no closed-form analytical solution, the convergence of moment functions can only be assessed numerically which will be done in Subsec. (ref). Thus, the hierarchy of convergences is directly informative in the context of the validation of ABMs. Firstly, the time averages of the moment functions need to converge at all and additionally to some value, $\nicefrac1T \sum_{t=1}^{\ensuremath{T}\xspace} \ensuremath{m}\xspace \left( x_t \right) \to \ensuremath{\overline{\ensuremath{m}\xspace}}\xspace$. Corresponding to the third aspect in the hierarchy, this convergence needs to go towards the expected value we observed empirically, $\ensuremath{\overline{\ensuremath{m}\xspace}}\xspace = \nicefrac1T \sum_{t=1}^{\ensuremath{T}\xspace} \ensuremath{m}\xspace \left( x_t \right) \to \ensuremath{\left\langle \ensuremath{m}\xspace \right\rangle}\xspace$. Eventually, proper moment functions have to rely be ergodic invariants of the mDGP.
Ergodic theory implies a further highly relevant aspect to the validation of FABMs which is its asymptotic nature. Roughly speaking, ergodic theorems live in the land of asymptopia, while reality and the estimation of (F)ABMs take place in pre-asymptopia with finite samples. For all practical purposes real-world data will always be limited even if the available computational resources are constantly expanding. If there are only finite observations in the time series of the FABMs in pre-asymptopia then the time scales of the convergences will become crucially important sooner or later. Any uncertainty about the LHS in Eq. ((ref)) vanishes only in the large-time limit.\footnote{ But there is also some uncertainty about the true expectation value on the RHS in Eq. ((ref)), which only vanishes in the large-ensemble limit.} In pre-asymptopia all quantities are computed under the constraint of finite computational resources. Thus we operate with samples finite in time and ensemble size, $\lim_{t\to\ensuremath{T}\xspace}, \ensuremath{T}\xspace \ll \infty$ and $\lim_{n\to\ensuremath{N}\xspace}, \ensuremath{N}\xspace \ll \infty$. Given any limited budget of computational resources, the uncertainties in computing moment functions of our FABM depend on such limits and behave differently and might never vanish sufficiently fast.
So far we have discussed (i) ergodicity in general, (ii) the importance of the ergodic property, (iii) assumptions hidden in the ergodic theorem about specific convergences that reappear in the selection of suitable moments of FABMs and (iv) the fact that validation of FABMs always takes place in pre-asymptopia where ergodicity is generically broken due to slow convergences. In this section we focus on the different dimensions where convergences take place when performing Monte Carlo simulations and how to visualize them. This section prepares for a better understanding of the effects of broken ergodicity along the different dimensions as discussed in Sec. (ref). It will turn out to be instructive for the choice of how to split a limited number of total observations over (i) the simulation length, (ii) the number of simulation runs and (iii) the number of Monte Carlo repetitions to improve validation efforts.
We model observables as realizations of a random variable $X$.\footnote{Throughout this paper we denote random variables by capital letters and their realizations by the respective lower case.} Repeated observations or realizations $x$ of the same random variable at regular intervals from time 0 to time $T$ form a time series of our observable, $\left\{ x_0, x_1, x_2, \ldots , x_T \right\}$, also denoted by $\left\{x_t\right\}_0^T$. Within this stochastic model a specific observed time series can therefore be interpreted as a realization of finite length of a stochastic process $\left\{ x_t \right\}_{t\in\ensuremath{\mathbb{T}}\xspace}$, where \ensuremath{\mathbb{T}}\xspace denotes the time domain of the time index. \footnote{The flexible notation of the time domain as \ensuremath{\mathbb{T}}\xspace allows to easily adopt different time domains. For a continuous time model the time parameter $t$ of the stochastic process is $\ensuremath{\mathbb{T}}\xspace = \ensuremath{\mathbb{R}}\xspace$ or $\left\{ x_t \right\}_{t\in\ensuremath{\mathbb{R}}\xspace}$. For a model in discrete time $\ensuremath{\mathbb{T}}\xspace = \ensuremath{\mathbb{N}}\xspace$ which yields the stochastic process $\left\{ x_t \right\}_{t\in\ensuremath{\mathbb{N}}\xspace}$.} A stochastic process, $\left\{ x_t \right\}_{t\in\ensuremath{\mathbb{T}}\xspace}$, is then a family of all infinite time series.
Our notation establishes the following relation between the real world (which is observable) and our model. Nature or the financial markets are thought of as a rwDGP which is unknown and of which we often observe only exactly one unique time series. On the other hand we have an economic model of a financial market -- a FABM in our case -- which is our model data-generating stochastic process (mDGP), from which we can generate an ensemble of (at least in principle) arbitrarily many $N$ time series of (finite) length $T$. At least two facts distinguish two time series generated by the same stochastic process. First, the two time series (or realizations of the stochastic process) differ in their random initial conditions or random seeds.DelliGattiEtAl2018 The dependence on some initial conditions is indicated in Fig. (ref) by the argument in the function, e.g.\xspace $X_t\left(x_1\right)$ signifies a process $X_t$ that started in $x_1$. Second, randomness in the noise realizations -- there are simply different realizations of the random variable that appear over time or are generated and used in the simulation.
Throughout the paper, the distinction between the time dimension and the ensemble dimension will be crucial. We denote time by $T$ and denote the size of the ensemble by $N$. A Monte Carlo simulation experiment is then fully determined by the
Figure (ref) contains a visualization of the three dimensions of simulations of (F)ABMs. Commonly, this third dimension of the Monte Carlo runs is understood as belonging to the ensemble dimension as only the random seeds might be different, but additional numerical effects appear that make a distinction between the two dimensions possible which are discussed in more detail in Sec. (ref).\footnote{Related two-dimensional visualisations in the context of tests for ergodicity can be found in GueriniMoneta2017. GrazziniRichiardi2015 refer to the three-dimensional mental model of a cube and refer to it as the \textquote{replications} dimension. Similar reasoning about correct alignment of research question and statistical identification led to the idiographic paradigm in behavioral psychology and physiology and similar visualisations MolenaarCampbell2009,NeumannEtAl2022.} The convergence behavior in the numerical estimation crucially depends on the ensemble size \ensuremath{N}\xspace, the time length \ensuremath{T}\xspace and the number of Monte Carlo simulations \ensuremath{M}\xspace in a non-intuitive and non-linear way. In principle different convergence behaviours in these three dimensions have not been investigated before. Our analysis is thus contributing to the joint community efforts of improving the validation of (F)ABMs.
Let us now turn to the simulated method of moment estimation. This estimation method is studied for different kind of models such as dynamic stochastic general equilibrium (DSGE) modelsRugeMurcia2012,RugeMurcia2013 as well as (F)ABMsFranke2009,FagioloRoventini2017,FagioloEtAl2019. In general SMM is more robust to misspecification than for example maximum likelihood methodsRugeMurcia2007 and performs better for large-scale modelsPlatt2020. SMM belongs to the broader class of simulated minimum distance methods whose goal is to identify the unknown parameters of the model, which generate the least distance between the simulated moment functions and the empirical observed moment functions.
Let us now introduce a simulated method of moment estimator. Therefore, we consider a FABM with unknown parameter vector $\ensuremath{ \ensuremath{\mathbold{\theta}} }\xspace \in \Theta$, which we ultimately want to estimate from the parameter space $\Theta$. In many empirical studies it is the primary goal to fit parameter values of the mDGP e.g.\xspace by matching empirical and simulated data moment functions. In our context, the goal is to assess the quality of the estimation approach itself, i.e.\xspace how well does the estimation approach identify some known 'true' parameter values \ensuremath{\ensuremath{ \ensuremath{\mathbold{\theta}} }\xspace_0}\xspace. If no empirical data is used for the evaluation of the estimation approach but data simulated from an ABM -- like in our study -- then it is a common practice to use known parameter values from the literature as the benchmark of the estimation that match the moment functions well and refer to them as 'true' values.\footnote{They are sometimes also referred to as 'pseudo-true' values as they are derived from synthetic/simulated data. See also Subsubsec. (ref).} Since SMM belongs to the broader class of simulated minimum distance methods, a distance function $\ensuremath{ \ensuremath{\mathbold{d}} }\xspace$ measures the difference between empirical and simulated data given some parameters \ensuremath{ \ensuremath{\mathbold{\theta}} }\xspace.\footnote{To be very precise, the empirical data is produced given the true parameters \ensuremath{\ensuremath{ \ensuremath{\mathbold{\theta}} }\xspace_0}\xspace, which we want the estimation to recover reliably.} Thus, we compute the distance between the two moment function vectors for the simulated and empirical moment functions $\ensuremath{ \ensuremath{\mathbold{d}} }\xspace = \ensuremath{\ensuremath{ \ensuremath{\mathbold{m}} }\xspace^{\text{emp}}}\xspace\left(x_t\right) - \ensuremath{\ensuremath{ \ensuremath{\mathbold{m}} }\xspace^{\text{sim}}}\xspace\left(x_t|\ensuremath{ \ensuremath{\mathbold{\theta}} }\xspace\right)$ given some parameters \ensuremath{ \ensuremath{\mathbold{\theta}} }\xspace. In the following, we explicitly consider two different types of averaging over simulated observables. First, simulated moment functions can be computed as time averages over one realization of observables $x_t$ of time length $\ensuremath{T}\xspace$:
Or second, they can also be captured as averages over an ensemble of size \ensuremath{N}\xspace of time averages of length \ensuremath{T}\xspace:
The objective function or criterion function for a given a set of model parameters \ensuremath{ \ensuremath{\mathbold{\theta}} }\xspace aggregates then the distances
for a given weighting matrix \ensuremath{\ensuremath{\mathbf{W}}}\xspace which is positive semi-definite. The estimator \ensuremath{\hat{\ensuremath{ \ensuremath{\mathbold{\theta}} }\xspace}}\xspace yields that vector of model parameters for which the objective function is minimized:
Under standard regularity conditions, the distance function \ensuremath{ \ensuremath{\mathbold{d}} }\xspace is assumed to be stationary and ergodic resulting in an asymptotically consistent estimator.LeeIngram1991,DuffieSingleton1993 For most (F)ABMs, however, a corresponding SMM estimator might not have such properties. In fact, we will show that in pre-asymptopia the following inequality holds due to non-commutativity of the limits in ensemble size, time and also of the number of MC runs:
Put simply, under broken ergodicity and/or in pre-asymptopia the order matters in which the limits are taken.Grazzini2011 We thus study the effect of an efficient allocation of a limited budget of observations and how to get an estimator with improved properties.
The efficiency of the SMM estimator is affected by the design of the weighting matrix \ensuremath{\ensuremath{\mathbf{W}}}\xspace. The optimal choice is given by a weighting matrix with the smallest asymptotic covariance for the estimator. One popular choice would be the use of the Newey-West estimator. In our numerical study, we are interested in the uncertainty that can be solely associated with broken ergodicity. Therefore, we will consider the inverse of the long-run covariance matrix of the true data as the optimal weighting matrix.
To summarize this section, we have discussed how the ergodic property plays a crucial role in moment selection. As we will see next, convergences along the three different dimensions might behave differently due to broken ergodicity in pre-asymptopia. How exactly shows the evaluation of the simulation experiments in the next section.
\FloatBarrier
To conduct our study of (broken) ergodicity in the convergence of moment functions, we take two established financial agent-based models: AlfaranoEtAl2008 (ALW model\xspace) and FrankeWesterhoff2012 (FW model\xspace). They are both based on a herding mechanism that has its roots in Kirman1993. Both models are commonly used for estimation exercises in the community.ChenLux2018,Lux2018,KukackaKristoufek2020,BertschingerMozzhorin2021 They replicate most of the stylized facts mentioned above and require comparatively little computational resources despite their prototype-nature.
In this subsection we briefly discuss the key mechanisms underlying the ALW model\xspace and the FW model\xspace which are the basis of our numerical experiments hereafter.
The first model by AlfaranoEtAl2008 incorporates a behavioural herding mechanism based on Kirman1993 with precedent analysis in AlfaranoEtAl2005. \footnote{For a simple variant of the herding model by Kirman1993 with local interaction between the agents (or equivalently extensive transition rates) a closed-form exists for which analytical solutions of the time-variation of moments and some related quantities of interest exist can be computed (see AlfaranoEtAl2005 and esp. AlfaranoEtAl2008).}
The ALW model assumes two types of financial speculators, fundamentalist and chartist traders. Fundamentalists' excess demand is given by $\ensuremath{D_\text{f}}\xspace = \ensuremath{N_\text{f}}\xspace \ensuremath{V_\text{f}}\xspace \left(p^{\star}_{t}-p_t\right)$ where \ensuremath{V_f}\xspace is the average demand of \ensuremath{N_f}\xspace fundamental speculators. The fundamental value $p^{\star}_{t}$ is assumed to follow a random walk $p^{\star}_{t}=p^{\star}_{t-1}+\ensuremath{\sigma_\text{f}}\xspace \cdot \epsilon_{f,t}$ with $\epsilon_{f}\sim \ensuremath{\mathcal{N}\left(0,1\right)}\xspace$. The excess demand of chartist traders is given by $\ensuremath{D_\text{c}}\xspace = \ensuremath{N_\text{c}}\xspace \ensuremath{V_\text{c}}\xspace x_t$. Chartists are in one of two opinion states, either optimistic or pessimistic. A sentiment index $x_t$ is defined as $x_t = \nicefrac{2n_t}{\ensuremath{N_\text{c}}\xspace}-1$ with $n_t$ optimistic traders at time $t$ of a total number of \ensuremath{N_c}\xspace chartists. They are assumed to change their sentiment based on the (extensive) transition rates $\pi^+ = a + bn$ and $\pi^- = a + b(N-n)$, where parameter $a$ indicates idiosyncratic switches and $b$ measures the herding intensity. The resulting sentiment dynamics is approximated by the following Langevin equation with drift component $A(x) = -2ax$ and diffusion term $D(x) = 2b(1-x^2_t) + \nicefrac{4a}{N}$ which gives
Eq. ((ref)) can be discretized with $\Delta t = 1$:
with $\epsilon\sim \ensuremath{\mathcal{N}\left(0,1\right)}\xspace$. The distribution of the sentiment index is known to be bimodal for $a<b$ and unimodal for $a>b$.
Price dynamics are governed by a standard Walrasian adjustment mechanism, depending on total excess demand of both trading groups:
where $\beta$ is the assumed price adjustment speed. With instantaneous market clearing, i.e.\xspace $\beta\to\infty$, and setting $\nicefrac{\ensuremath{N_\text{c}}\xspace\ensuremath{V_\text{c}}\xspace}{\ensuremath{N_\text{f}}\xspace \ensuremath{V_\text{f}}\xspace}=1$, we get the following evolution of returns:
Finally, the parameter vector to be estimated for the ALW model\xspace contains three items $\ensuremath{\ensuremath{ \ensuremath{\mathbold{\theta}} }\xspace^{\text{ALW}}}\xspace = \left( a,b,\ensuremath{\sigma_\text{f}}\xspace \right)^\top$. In accordance with the literature we use for all our simulation experiments the following true model parameters $\ensuremath{\ensuremath{ \ensuremath{\mathbold{\theta}} }\xspace_0}\xspace = \left( 0.3, 1.4, 30 \right)^\top$. \footnote{ See GhonghadzeLux2016,ChenLux2018. The (pseudo-)true parameter setting for the ALW model\xspace that is used in this paper represents the bimodal case of the underlying sentiment index given that $b > a$. Since ChenLux2018 have shown that the bimodal case gives the best fit to empirical data, we will only focus on this model scenario. Note that for better readability parameters of the ALW model\xspace are always multiplied by $10^3$ throughout the paper. }
FrankeWesterhoff2012 propose an entire model zoo for which they run a model contest. Here, we will only consider their top performing DCA-HPM model version (discrete choice with herding, predisposition and price misalignment). FrankeWesterhoff2012 apply a market maker model that considers as well two types of speculators: chartist and fundamentalists, whose fractions are denoted by \ensuremath{n_c}\xspace and \ensuremath{n_f}\xspace, respectively. The evolution of the log prices is then determined by
with $\mu$ reflecting the speed of price adjustment. Excess demand of chartists \ensuremath{D_c}\xspace and fundamentalists \ensuremath{D_f}\xspace is given by
with $\ensuremath{\epsilon_\text{f}}\xspace \sim \ensuremath{\mathcal{N}\left(0,\ensuremath{\sigma_\text{f}}\xspace^2\right)}\xspace$ and $\ensuremath{\epsilon_\text{c}}\xspace \sim \ensuremath{\mathcal{N}\left(0,\ensuremath{\sigma_\text{c}}\xspace^2\right)}\xspace$. Parameters $\chi$ and $\phi$ are always positive and indicate the strength of reaction. The switching between both trading strategies is governed by a fitness function measuring the attractiveness $a_t$ of fundamentalism over chartism:
where $\alpha_0$ represents a constant idiosyncratic predisposition for one of the two trading strategies, $\alpha_n$ relates to herding intensity and $\alpha_p$ accounts for price misalignments from the fundamental value $p^\star$. Note that while $\alpha_n$ and $\alpha_p$ are always strictly positive, parameter $\alpha_0$ might be negative as well. The current market shares of fundamentalists \ensuremath{n_{f,t}}\xspace and chartists \ensuremath{n_{c,t}}\xspace at time $t$ are then updated according to the following discrete choice approach:
As in FrankeWesterhoff2012 we set the intensity of choice parameter to $\beta=1$, the speed of price adjustment parameter to $\mu=0.01$ and the fundamental price to $p^*=0$. Finally, the parameter vector to be estimated for the FW model\xspace contains a total of seven parameters $\ensuremath{\ensuremath{ \ensuremath{\mathbold{\theta}} }\xspace^{\text{FW}}}\xspace = \left( \phi,\chi,\alpha_0,\alpha_n,\alpha_p,\ensuremath{\sigma_\text{f}}\xspace,\ensuremath{\sigma_\text{c}}\xspace \right)^\top$. In accordance with FrankeWesterhoff2012 we use for all our simulation experiments the following true model parameters $\ensuremath{\ensuremath{ \ensuremath{\mathbold{\theta}} }\xspace_0}\xspace = \left(0.12,1.5,-0.336,1.839,19.671,0.708,2.147\right)^\top$.
To sum up, the ALW model\xspace by AlfaranoEtAl2008 is a small-scale model which can be reduced to only three estimation parameters $\ensuremath{\ensuremath{ \ensuremath{\mathbold{\theta}} }\xspace^{\text{ALW}}}\xspace = \left(a,b,\sigma_f\right)^\top$. The FW model\xspace by FrankeWesterhoff2012 is more complex with a total of seven parameters $\ensuremath{\ensuremath{ \ensuremath{\mathbold{\theta}} }\xspace^{\text{FW}}}\xspace= \left( \phi,\chi,\alpha_0,\alpha_n,\alpha_p,\ensuremath{\sigma_\text{f}}\xspace,\ensuremath{\sigma_\text{c}}\xspace \right)^\top$. Sample simulation runs for both model dynamics can be found in appendix (ref).
The following subsection builds on Sec. (ref) and explains the motivation behind the choice of our set of moment functions. As explained in Subsec. (ref), we focus on two agent-based asset pricing models which replicate many of the stylized facts of financial markets. The vector of moment functions \ensuremath{ \ensuremath{\mathbold{m}} }\xspace should contain a set of reasonable summary statistics that (partly) capture these stylized facts expressed as observables and thus measurable statistical quantities. A necessary condition for the choice of moment functions is given by the order condition, i.e.\xspace for estimations with more than one model parameter the number of moment functions \ensuremath{K}\xspace needs to be greater or equal than the number of model parameters or the cardinality of the parameter vector \ensuremath{ \ensuremath{\mathbold{\theta}} }\xspace. Thus, the order condition provides a lower bound for the number of moment conditions. Theoretically, the (full) rank condition is a sufficient condition for identification assuring that only moment functions without linear dependence are included. While the order condition is easy to meet, the rank condition is barely testable for most (F)ABMs given their non-linearity and the lack of analytical closed-form expressions.
Besides these restrictions the number and choice of moment conditions is not strictly limited which may render it arbitrary. However, this is not a weakness per se\xspace, because SMM is designed to capture complex patterns which likely escape any single moment condition. Thus, if empirical data show complex patterns that are hard to squeeze into a single metric, it can be necessary to wrap the patterns in more than a single moment condition. The stylized facts of financial markets show such complex patterns like bubbles and crashes, excess volatility, heavy-tailed return distributions, absence of autocorrelation in raw returns or slow decay in volatility. Especially the statistical pattern of slow non-linear decay in autocorrelation over many lags requires more than one moment condition. This justifies the comparatively large size of our vector of moment functions with $\ensuremath{K}\xspace = 18$ listed in Table (ref). For the SMM approach this is a common size of the moment vector, see e.g.\xspace ChenLux2018 who analyse the ALW model\xspace with up to 15 moment functions or RugeMurcia2012 who uses 16 moment functions for the estimation of a macroeconomic DSGE model.
As the first two moment functions we use the mean of absolute returns and the variance of the return series. As mentioned in Sec. (ref), ordinary moments in the sense of Eq. ((ref)) for HTDs do only exist for orders lower than the tail index. For example, if a return series has a tail index of $\alpha = 3$, only the first and second moment exist. This means that sample estimates of the third or higher moments will never converge with increasing sample size. Hence, such moments might not be suitable as moment functions for estimation purposes. We yet consider the kurtosis of the return series as one moment function and additionally include the Hill tail index estimates at \SI{2.5}{\percent} and \SI{5.0}{%}. The autocorrelation of raw returns at lag 1 checks for the absence of serial correlation in returns. Finally, in order to capture the slow decay of volatility we include autocorrelations of absolute and squared returns for the lags $1,5,10,25,50$ and $100$.
Table (ref) offers an extensive yet non-comprehensive list of possible summary statistics considered for the estimation of financial market models. We abstain from modelling moment functions as further derived processes, like a GARCH$(1,1)$, firstly because of the additional computational cost associated with the estimation of the GARCH parameters. Secondly and more importantly, they introduce additional uncertainty leading to possible biases in the estimates.
\afterpage{
}
In order to apply any moment estimator, it is important that simulated moments (asymptotically) equal their true values. This means that true values exist and convergence happens with an increasing number of observations (either over time \ensuremath{T}\xspace and/or the ensemble \ensuremath{N}\xspace). Otherwise, meaningful inference is not possible. Hence, the very first step of our analysis is to numerically check if the moment functions do exist, i.e.\xspace if they quickly converge to identical values for some finite number of observations. Since the mDGP for most (F)ABMs is not available in closed form, the ensemble average \ensuremath{\left\langle X \right\rangle}\xspace is not readily at hand. Living in a pre-asymptotic numerical simulation, sufficient convergence to the theoretical average needs to take place rather quickly for finite observations. Recall that we talked about the conditions for the existence of higher moments in the previous section. To find the moments' true values, we run simulations over different time lengths \ensuremath{T}\xspace with $\ensuremath{T}\xspace_1=\num{10 000}$, $\ensuremath{T}\xspace_2=\num{100 000}$, $\ensuremath{T}\xspace_3=\num{1 000 000}$ and $\ensuremath{T}\xspace_4=\num{10000000}$. We run each scenario over $\ensuremath{M}\xspace = \num{5000}$ replications. Results are presented in Tab. (ref). For each moment function, we report mean, variance and $p$-value of a Kolmogorov-Smirnov test statistic.
We observe convergence behaviour for all moment functions for both the ALW model\xspace and FW model\xspace for an increased number of observations. In fact, the variance decreases consistently for all considered moment functions and models. Given that the variance is low for all considered moment functions, we can assume that the long-run estimates are informative for the estimation. Therefore, we will consider the estimates of the longest time length $T_4$ to compute the inverse of the long-run covariance matrix as the optimal weighting matrix in our SMM estimator. High $p$-values for the KS-test also suggest normally distributed samples of the moment estimates. Given that the true values exit and converge, we can move on to the next step which is to test if the chosen moment functions are suited for identification of the model parameters.
\FloatBarrier
In this subsection we are concerned with the pre-asymptotic properties of our moment functions and their identification power. We want to check if the mapping between model parameters and model output in terms of moment functions is unique. Only a one-to-one mapping works as a sufficient condition (i.e.\xspace full rank condition) for identification. For the one parameter case, the identification issue is quite obvious. A sufficient condition for identification would be a strictly monotone relationship between model parameter and moment function. For higher dimensional models, this is less trivial. In fact, as stated earlier, the rank condition is barely directly testable for ABMs.
Nevertheless, running sensitivity analyses allows us to study how strongly moment functions react towards changes in model parameters.\footnote{Since not all moment functions are equally informative, another way to identify relevant moment functions is in terms of statistical efficiency. GallantTauchen1996 suggest using scores of an auxiliary model as moment functions in a generalized method of moment (GMM) estimator in order to reduce loss of efficiency due to uninformative moments. Hence, by searching for the statistically most informative moments, such an approach tries to reach the efficiency of maximum-likelihood. Auxiliary models are considered to be especially attractive when the number of moment functions is rather limited.} Therefore, we break down the model's dimensionality and study how sensitive our moment functions are to variations of the model parameters.
For the ALW model\xspace, we sample 151 equidistant points over a large range of parameter variations: $a \in [0.06, 0.66]$, $b \in [0.28, 3.28]$ and $\sigma_f \in [6.0, 66.0]$. We plot the relative responsiveness of moments $m_1$, $m_2$, $m_3$, $m_4$, $m_6$, $m_7$ and $m_8$ for the ALW model\xspace in Fig. (ref).\footnote{Moment functions $m_9$ to $m_18$ for the autocorrelations of absolute and squared returns at higher lags show qualitatively the same functional behaviour as $m_7$ and $m_8$, respectively.} We find non-linear relationships for all moment-parameter pairs. Moment function $m_6$ (autocorrelation of raw returns at lag 1) shows highly non-smooth behaviour towards changes in model parameters $a$ and $b$ of the ALW model\xspace. Regarding the strength of reactions, moment functions react most sensitive to parameter $\sigma_f$. This is in line with earlier findings by ChenLux2018. In fact, the moment functions are least sensitive towards varying values of parameter $a$. \afterpage{
}
For the FW model\xspace, we also run simulation for 151 equidistant points over a broad range of parameter variations: $\phi \in [0.02, 1.52]$, $\chi \in [0.75, 2.25]$, $\alpha_0 \in [-0.4, 2.0]$, $\alpha_n \in [0.839, 2.339]$, $\alpha_p \in [11.671, 23.671]$, $\ensuremath{\sigma_\text{f}}\xspace \in [0.58, 0.88]$ and $\ensuremath{\sigma_\text{c}}\xspace \in [1.447, 2.947]$. Results are shown in Fig. (ref). We find again highly non-linear functional relationships for almost all moment-parameter combinations. We further see parameter $\alpha_0$ as probably hard to identify correctly given the wild and untamed nature of the moment function responses. Additionally, moment functions react comparatively insensitive towards parameter $\chi$ suggesting that it might trigger identification problems, too.
To conclude, sensitivity analyses reveal that most functional relationships are non-linear which might result in (small) biases of the estimates. As mentioned in GrazziniEtAl2012,GrazziniRichiardi2015, if a moment function happens to be non-linear, there will be a small bias with the direction of the bias depending on moment functions' derivatives. Thereby concave (convex) moment functions lead to an upward (downward) bias. Such biases can be reduced in different ways: (i) given the analytical expression of the moment function, however, this is unknown for most ABMs, (ii) through monotonic transformations and (iii) through increased number of observations. The latter brings us to the focus of our study.
Recall that biases decrease for increasing number of observations since the simulated moment functions converge to their true theoretical value. Hence, the next step is to take a closer look at the convergence speed of the moment functions. The goal here is twofold. First, we want to check if an ensemble \ensuremath{N}\xspace of simulation runs converges to the exact same value as the corresponding long-run realization, i.e.\xspace how strong the effect of broken ergodicity is. Second, we aim at evaluating the speed of convergence.
For this experiment, we run simulations over $\ensuremath{M}\xspace = \num{5000}$ repetitions for two different scenarios: (i) $\ensuremath{T}\xspace = \num{400000}$, $\ensuremath{N}\xspace = 1$ and (ii) $\ensuremath{T}\xspace = \num{40000}$, $\ensuremath{N}\xspace = 10$. In Fig. (ref) we plot again the convergence behaviour for the moment functions $m_1$, $m_2$, $m_3$, $m_4$, $m_6$, $m_7$ and $m_8$ defined as relative deviations from their true values, i.e.\xspace their long-run estimates reported in Subsubsec. (ref). The top and lower panels of Fig. (ref) and Fig. (ref) show results for scenario (i) and (ii), respectively. The left (right) panels refer to the selected moment function including $m_6$ (excluding $m_6$). We find for both the ALW model\xspace and FW model\xspace that moment functions $m_1$, $m_2$ and $m_4$ closely converge towards their true values, while $m_3$, $m_7$ and $m_8$ show more or less pronounced and persistent biases. Regarding the speed of convergence, we observe long transition phases for scenario (i) and comparatively fast adjustments for scenario (ii). Yet, deviations are more severe here, too. Looking at the right panels including moment function $m_6$ which is the autocorrelation of raw returns at lag 1, we find very volatile and probably non-convergent behaviour for both models. This suggests that moment $m_6$ might not be suited for identification of the given model outcomes.
\FloatBarrier
Next, we study the properties of the SMM's objective function with explicit regard to the impact of broken ergodicity. Let us start with investigating the objective function's response surface. The aim is to check for discontinuities, (non-)smoothness and flat valleys, which create problems during the optimization as they impede finding global optima during the estimation. Such plots provide further help in gaining an understanding of the mapping between model parameters and the objective function.
\FloatBarrier
For this analysis, we will focus only on the most problematic model parameters in terms of identification power. For the ALW model\xspace, we have identified both herding parameters $a$ and $b$ as troublesome. We run simulation scenarios summarized in Tab. (ref) with a constant observation budget of \num{8000000} throughout all experiments, i.e.\xspace one short time series realization $\ensuremath{\ensuremath{T}\xspace_{\text{short}}}\xspace = \num{40 000}$, a long time series with $\ensuremath{\ensuremath{T}\xspace_{\text{long}}}\xspace = \num{400 000}$ and over an ensemble with $\ensuremath{N}\xspace = 10$ and $\ensuremath{T}\xspace = \num{40 000}$. The parameter grid consists of 41 equidistant points between $\left[0.06,0.66\right]$ for parameter $a$ and between $\left[0.28,3.1\right]$ for parameter $b$. The moment conditions are based on all 18 moment functions of Tab. (ref). Figure (ref) shows the corresponding results for three different sample simulation runs using different random seeds. The solid lines mark the true value of the model parameters $\left(a,b\right)=\left(0.0003,0.0014\right)$. The intersection of the two lines build the theoretical optimum of the objective function. Note that we plot the response surface of the inverse of the objective function to better visualize the global optimum. Ideally, the highest peak on the response surface coincides with the intersection of the true parameter values.
The response surface of the objective function is clearly non-smooth and contains multiple local optima. Depending on the choice of random seeds, the fitness of the objective function differs greatly. While parameter $b$ is correctly identified in all three simulation samples, estimates of parameter $a$ happen to be far (left plot) to slightly off (middle plot) compared with the true value of $a=0.0003$. One might be even lucky enough to identify both parameters correctly (right plot).
Increasing the number of observations from $\ensuremath{T}\xspace = \num{40 000}$ to $\ensuremath{T}\xspace = \num{400 000}$, we find a smoother surface as well as less variations as expected. However, even for this higher sample size, local optima arise (see the middle plot) depending again on the choice (or luck) of random seeds. Comparing this with an ensemble of $\ensuremath{N}\xspace = 10$ simulation runs for $\ensuremath{T}\xspace = \num{40 000}$, we qualitatively observe the same patterns.
During a pre-exploratory experiment for the FW model\xspace, we identified the most problematic parameter pairs ($\chi, \alpha_0$) and ($\chi$,$\alpha_p$). The latter pair is visualized in Fig. (ref).\footnote{Graphical results for the objective function's response surface of the FW model\xspace are almost identical for parameter pairs ($\chi$,$\alpha_0$) and ($\chi$,$\alpha_p$).} In the following, we set the (pseudo-)true values $\left(\chi,\alpha_p\right)=\left(1.5,19.671\right)$ and grid points result from 41 equidistant parameter variations over $\left[0.3,2.7\right]$ for parameter $\chi$ and over $\left[3.9342,3.9342\right]$ for parameter $\alpha_p$. Let us first note that we find qualitatively the same results as for the ALW model\xspace. Yet, the contour plots look rather different for the FW model\xspace. We observe that the surface is less elevated with peaks at much lower levels than we see for the ALW model\xspace. Especially for the \ensuremath{\ensuremath{T}\xspace_{short}}\xspace-setting, we find very flat hilltops for parameter $\chi$ and even flat plateaus for both parameters. Increasing the number of observations helps to smoothen the surface. However, identification issues leading to biases in the estimates might be more pronounced here which is most probably due to the higher dimensionality and complexity of the \textsc{FW} model\xspace.
While these contour plots of the objective function reveal a partially non-smooth surface, they allow no assessment of the estimator's further distributional properties. To do so, we need to include the third dimension of our simulation cube shown in Fig. (ref) which is referred to the Monte Carlo repetitions. Hence, we run the simulations repeatedly over $\ensuremath{M}\xspace = 200$ and present graphical results in Figure (ref). The left panel shows boxplots for the three different simulation settings. First, we notice that for the ALW model\xspace the estimator is consistent and overall able to detect the true value of parameter $a$ on average. Unsurprisingly, variations are comparatively larger for the \ensuremath{\ensuremath{T}\xspace_{short}}\xspace-setting. Regarding the other two settings, variations differ only slightly with almost unbiased estimates for the \ensuremath{\ensuremath{T}\xspace_{long}}\xspace-setting and a small downward bias in the mix-setting. The paired histogram plot on the right panel in Fig. (ref) confirms this observation. While for the \ensuremath{\ensuremath{T}\xspace_{long}}\xspace-simulation the histogram is quite symmetric, the histogram for the mix-simulation is skewed towards lower values of $a$.
For the FW model\xspace results differ slightly. We find small biases of the estimates for the \ensuremath{\ensuremath{T}\xspace_{short}}\xspace and \ensuremath{\ensuremath{T}\xspace_{long}}\xspace-cases while the \ensuremath{\ensuremath{T}\xspace \times \ensuremath{N}\xspace }\xspace-setting is able to perfectly identify the true value of parameter $\chi$ on average. This suggest that the ensemble setting might lead to more consistent and efficient estimates. However, we are cautious to not overinterpret our results here given that we have included only two varying parameters while keeping the rest fixed.
This brings us directly to the next experiment where we include all model parameters for estimation. In the following, we will run a small estimation exercise over a Sobol-sequenced sample of $\num{2000}$ parameter combinations. The included parameter space is given by $\pm25\%$ of the true values which results in the following parameter ranges: $a \in \left[0.225,0.375\right]$, $b \in \left[1.05,1.75\right]$, $\sigma_f \in \left[22.5,37.5\right]$ for the ALW model\xspace and $\phi \in \left[0.09,0.15\right]$, $\chi \in \left[1.125,1.875 \right]$, $\alpha_0 \in \left[-0.252,-0.42\right]$, $\alpha_n \in \left[1.379,2.299\right]$, $\alpha_p \in \left[14.753,24.589\right]$, $\ensuremath{\sigma_\text{f}}\xspace \in \left[0.531,0.885\right]$, $\ensuremath{\sigma_\text{c}}\xspace \in \left[1.61,2.684\right]$ for the FW model\xspace.\footnote{The defined parameter ranges provide broad parameter variations while keeping execution times of the simulation runs feasible.} The number of Monte Carlo runs is set again to $\ensuremath{M}\xspace = 200$. Since we are interested in the degree of uncertainty that solely comes from broken ergodicity, we abstain from running a full estimation exercise given that the choice of the optimization algorithm would highly influence and distort our results.\footnote{WinkerMaringer2009 provide a study of the joint convergence of an estimator and a heuristic optimization algorithm. They compute the necessary number of Monte Carlo repetitions of the optimization routine to derive robust estimates.} We report graphical results in Fig. (ref) and the corresponding mean estimates together with their standard deviations and root-mean squared errors in Tab. (ref). Looking at Fig. (ref), the left vertical axis shows the value of the objective function \ensuremath{J}\xspace. Settings (a) and (b) result in the lowest values of \ensuremath{J}\xspace. For the ALW model\xspace, it seems that both (a) and (b) are compatible in terms of minimized objective function. Mean sample estimates in Tab. (ref) confirm these findings. Decreasing the length of \ensuremath{T}\xspace while increasing the ensemble size \ensuremath{N}\xspace raises the level of uncertainty in the estimates, pushing mean estimates slowly away from their true values. Accordingly, standard deviations and root-mean squared errors tend to increase, too.
Regarding the FW model\xspace and its comparatively higher dimensionality, results are qualitatively similar yet the loss in efficiency for increasing ensemble sizes \ensuremath{N}\xspace is more pronounced here. This loss in efficiency becomes particularly obvious for setting (d) with $\ensuremath{T}\xspace = \num{10 000}$ and $\ensuremath{N}\xspace = 40$. In fact, the number of observations per realized time series is too short here resulting in uninformative samples.
Summing up, we find different values of the objective function \ensuremath{J}\xspace for different allocations of the total budget of observations. While the differences might be negligible for longer time lengths \ensuremath{T}\xspace, they significantly deteriorate estimates for shorter time lengths due to broken ergodicity. Recall the ergodic theorem in Eq. ((ref)) which is effective here.
However, investigating only the finite sample properties of the parameter estimates neglects the computational resources needed to run these simulations. Therefore, we add the execution times as black stars, belonging to the right vertical axis, to Fig. (ref). Setting (a) is by far the most computationally demanding simulation. The three \ensuremath{\ensuremath{T}\xspace \times \ensuremath{N}\xspace }\xspace-settings (b), (c) and (d) have similar execution times which reduce the benchmark setting (a) by a factor of three for both the ALW model\xspace and FW model\xspace. The reason is that single simulations over \ensuremath{\ensuremath{T}\xspace_{long}}\xspace are memory intense and inefficient. In fact, repeated simulations over different random seeds are so called embarrassingly parallel tasks and can, thus, be easily distributed over multiple processing units which is done here.\footnote{All numerical simulations in this paper have been performed using Julia 1.0.5 on a standard desktop computer with an eight core 3.00 Ghz AMD Ryzen 7 1700 CPU and 32GB of memory. We made sure to run all simulation settings most efficiently to keep the different scenarios comparable.} We finally run a fifth scenario (e) for which we take the best performing \ensuremath{\ensuremath{T}\xspace \times \ensuremath{N}\xspace }\xspace-setting (b) with $\ensuremath{T}\xspace = \num{40 000}$ but double the amount of ensemble size $\ensuremath{N}\xspace = 20$. This setting is able to easily outperform all other settings, even scenario (a) while taking only two third of setting (a)'s computational time. We find this to be a consistent pattern for both FABMs considered here.
With this analysis, we are not only aiming for an adequate mixture of \ensuremath{\ensuremath{T}\xspace \times \ensuremath{N}\xspace }\xspace from the perspective of (broken) ergodicity but also considering the associated computational costs. We conclude that for a limited budget of computational resources, using a mix of \ensuremath{\ensuremath{T}\xspace \times \ensuremath{N}\xspace }\xspace is favourable over a single long simulation run with the same number of observations.
\FloatBarrier
How robust are these numerical results especially regarding variations of the moment function set? In the following, we try to answer this question by running sensitivity analyses for varying sizes \ensuremath{K}\xspace of the moment vector \ensuremath{ \ensuremath{\mathbold{m}} }\xspace. For this, we randomly sample 17 different moment sets per number of moment vector size and evaluate them for the numerical estimation data of setting (a). Fitness is defined as the euclidean distance between the (pseudo-)true and estimated parameter vector. We decide to investigate two different scenarios here. For the first scenario we sample from the complete range of 18 moment functions as listed in Tab. (ref). For the second scenario we exclude $m_6$ known as the autocorrelation of raw returns at lag 1. The reason is that we have seen for example in Fig. (ref) problematic to non-convergent behaviour for moment function $m_6$. We suspect improved results in terms of fitness given that $m_6$ offers very limited identification power. We present graphical results in Fig. (ref). The two figures show boxplots for the ALW model\xspace on the left and the FW model\xspace on the right. Unsurprisingly, results improve for an increasing number of included moment functions. More interestingly and rather unexpectedly, the improvement does not happen linearly. While variations in the fitness are quite large for smaller moment vector sizes, we observe a sudden drop for moment vector sizes of nine (ALW model\xspace) and ten (FW model\xspace), respectively. After that drop results are pretty robust for further increasing number of moments. This is especially true for the FW model\xspace while we see slight improvements for the ALW model\xspace. Comparing both scenarios (with and without $m_6$, graphically represented by dark red and orange) we find consistently better fitness for moment sample sets excluding the autocorrelation of raw returns at lag 1. Additionally, estimates related to this scenario tend to produce less outliers (see the grey and black dots).
This suggest that re-evaluating our estimation results from Tab. (ref) for an adjusted moment vector should lead to significantly better results. Therefore, we take our previous results as presented in Tab. (ref) and recompute them for a reduced moment vector excluding $m_6$. We relabel these results $J_{18}$ since they have been actually estimated for a complete set of 18 moments. Now we can compare them with the $J_{17}$ results that we estimated based on 17 moment functions. We report the differences between $J_{18}$ and $J_{17}$ values for all five settings in Tab. (ref). As we anticipated, we find clearly improved \ensuremath{J}\xspace-values as a common and consistent pattern. Furthermore, one-tailed $t$-tests reveal statistical significance for all results except for setting (d).
\FloatBarrier
(F)ABMs are often conceived as black boxes lacking a sound and well-behaved mathematical model. This has caused considerable resistance in the (economics) community to accept computer simulations as proper research methods. The development of estimation and validation tools that are particularly suited for the special properties of most (F)ABMs is a very active field of research. This paper aims to contribute to this line of research by explicitly considering the uncertainty coming from broken ergodicity.
We systematically study the properties and convergence behaviour of a SMM estimator and its individual moment functions. We have seen that assuring a one-to-one mapping between model output in terms of moment functions and model parameters is not trivial. Most moment functions are indeed non-linear and non-monotone. This can lead to biased estimations. More important is however the relative responsiveness of moment functions towards changes in the model parameters, i.e.\xspace how strongly they react. Parameters for which we find weak responsiveness of moment functions tend to be the most troublesome during estimation. We have further shown that not all moment functions which are a considered to be a popular choice for the estimation of univariate asset pricing models are actually suited. As it turns out, the convergence of the autocorrelation of raw returns at lag 1 is highly distorted making it an unfavourable candidate statistic. We show that estimation results improve significantly when leaving out this summary statistic. For all other moment functions we find robust results for the moment matching fitness. Therefore, our analysis confirms our choice of moment functions for (univariate) asset price models and thereby offers a response towards the arbitrariness critique of SMM.
Estimation methods like SMM are perfect tools for perfect models. Living in an ever-changing world, it should be no design fault to create models that possess non-ergodic properties. Even if regularity conditions are not fully met, models and tools may still be useful. Being aware of broken ergodicity is the key here. We further suggest researchers to not blindly follow previous studies based on other models. In fact, every (F)ABM behaves differently as we have seen here. Therefore, we run a bunch of Monte Carlo analyses to learn about the FABMs' sensitivity and responsiveness towards changes in model parameters and simulation settings. We have shown that for a given computational time budget an adequate mixture of ensemble size \ensuremath{N}\xspace and time length \ensuremath{T}\xspace is better suited than the same number of observations in one realization. We can conclude that due to careful selection of moment functions and consideration of ergodicity, we can improve the objective function's potential to correctly identify model parameters for a limited budget of computational resources.
So far, the energy demand of simulation models seems to play a minor role. Access to institute-owned high performance computing clusters and the recent global increase in cloud computing may have made the necessity for energy-efficient simulations obsolete. Yet, the opposite is true given the looming threat of climate change and a likely exponential increase in demand for simulations like FABMs especially for policy purposes in the near future. A holistic approach to fight global warming will include many different aspects of our lives. This also includes an efficient and mindful use of resources when it comes to research. Our work provides steps into this direction.
We have limited our focus in this paper on the use of an SMM estimator. The scale of our work can be extended to including other estimation methods as well. We specifically aim to apply our approach to likelihood-based methods. Since computing the likelihood function is expensive in terms of computational costs, applying our insights might improve cost efficiency there too.
\printbibliography