SI

#  Human mobility is well described by closed-form gravity-like models learned
automatically from data

Oriol Cabanas-Tirapu  Department of Chemical Engineering, Universitat Rovira i
Virgili, 43007 Tarragona, Catalonia  Lluís Danús  Department of Chemical
Engineering, Universitat Rovira i Virgili, 43007 Tarragona, Catalonia  Esteban
Moro  Institute for Data, Systems, and Society, Massachusetts Institute of
Technology, Cambridge, MA 02139  Department of Mathematics and GISC,
Universidad Carlos III de Madrid, 28911 Leganés, Spain  Network Science
Institute, Northeastern University, Boston, MA 02115, United States  Marta
Sales-Pardo  Department of Chemical Engineering, Universitat Rovira i Virgili,
43007 Tarragona, Catalonia  Corresponding authors: Marta Sales-Pardo (E-mail:
marta.sales@urv.cat); Roger Guimerà (E-mail<EMAIL_ADDRESS> Roger
Guimerà  Department of Chemical Engineering, Universitat Rovira i Virgili,
43007 Tarragona, Catalonia  ICREA, 08007 Barcelona, Catalonia  Corresponding
authors: Marta Sales-Pardo (E-mail: marta.sales@urv.cat); Roger Guimerà
(E-mail<EMAIL_ADDRESS>
###### Abstract

Modeling of human mobility is critical to address questions in urban planning
and transportation, as well as global challenges in sustainability, public
health, and economic development. However, our understanding and ability to
model mobility flows within and between urban areas are still incomplete. At
one end of the modeling spectrum we have simple so-called gravity models,
which are easy to interpret and provide modestly accurate predictions of
mobility flows. At the other end, we have complex machine learning and deep
learning models, with tens of features and thousands of parameters, which
predict mobility more accurately than gravity models at the cost of not being
interpretable and not providing insight on human behavior. Here, we show that
simple machine-learned, closed-form models of mobility are able to predict
mobility flows more accurately, overall, than either gravity or complex
machine and deep learning models. At the same time, these models are simple
and gravity-like, and can be interpreted in terms similar to standard gravity
models. Furthermore, these models work for different datasets and at different
scales, suggesting that they may capture the fundamental universal features of
human mobility.

## Introduction

Accurate models of population mobility within and between municipalities are
critical to address questions in urban planning and transportation
engineering. Additionally, since municipalities are the main ground on which
societies and cultures develop today, such mobility models are also
instrumental in addressing global challenges in sustainability, public health,
and economic development. Two main factors have driven recent interest in
modeling human mobility patterns [1, 2, 3, 4, 5]. First, accurate models of
human mobility could help identify transportation needs [6], allocate services
and amenities (shopping, health, parks) more efficiently [7], or even
understand and eventually alleviate problems like segregation [8], or epidemic
spreading [9]. But, at the same time, models of human mobility can help
identify the main behavioral components driving people to make large
displacements to, for example, buy a new product, find a new house, or use
physical activity spaces. Better behavioral models can help us implement more
efficient policies to change people’s behavior, rather than urban
environments, in favor of more sustainable attitudes.

Figure 1:  Modeling approaches for mobility flows between test municipalities
in Texas, US. (A) Real mobility flows between municipalities in the test set
in Texas, US (Methods). For each flow, we consider origin $o$ and destination
$d$ features, such as population $m_{o/d}$, aggregate statistics about points
of interest (POI), and the distance between them. (B) Flows predicted by the
deep gravity model [10], which uses a total of 39 features from origin and
destination. (C) Flows predicted by the closed-form, median predictive model
identified by the Bayesian machine scientist (BMS; see text). This model only
uses the population of origin and destination, as well as the distance between
them (Fig. 6B).

Despite these considerations, our understanding of the mobility flows within
and between urban areas is still incomplete. One of the earliest and most
fruitful attempts to model mobility flows between municipalities is the so-
called gravity model [1]. This model assumes that mobility flows depend solely
on the attractiveness or opportunities of the municipalities of origin and
destination (for which population is typically used as a proxy) and the
geographical distance between them, in a fashion that is mathematically
similar to Newton’s law of gravitation. In its different incarnations and
refined versions [11, 4, 12], the gravity model provides a simple
phenomenological description of a very complex phenomenon. Because of this,
while gravity models are not without their limitations, they are very often
used in urban design, transportation, or even commercial applications.
Recently, deep learning algorithms have been proposed, extending the ideas
underlying gravity models; they incorporate many other features besides the
populations of the origin and destination municipalities and their distance
[10]. Although those sophisticated machine learning tools are more accurate at
predicting flows between urban areas, they lack the explanatory power,
analytical tractability, adaptability to different contexts, and connection to
human decision-making of simple gravity models.

Given the reasonable success of simple gravity models in explaining flows in
urban areas, here we investigate the fundamental question of whether we really
need non-interpretable models that are much more complex than the gravity law
to delve deeper into the essence of urban mobility.

Unlike other behavioral models, gravity mobility models are phenomenological.

Because of the lack of precise theoretical underpinnings, their predictive
ability depends on the exact functional specification of the dependency of the
mobility flows on the model features; that is, the mathematical dependency on
origin and destination populations and distance. Here, we leverage recent
developments in Bayesian symbolic regression to obtain closed-form,
interpretable models[13] of mobility from data in a principled and automatic
fashion [14, 15, 16].

In particular, we systematically compare the performance at predicting
mobility flows of simple gravity models, complex machine learning and deep-
learning methods, and closed-form, interpretable models obtained through
Bayesian symbolic regression (Fig. 1). We find that the Bayesian symbolic
regression approach yields simple models that perform better, overall, than
any of the other modeling approaches. Our approach is able to learn accurate
models that, like gravity models, solely take into account the origin and
destination populations and the geographical distance between them.
Importantly, the learned models are gravity-like in their mathematical
dependencies on populations and distance.

Furthermore, exploration of the relationship between the contribution of the
populations of municipalities and their relative distance reveals common
patterns in all the datasets, which suggests a close to universal relationship
between mobility flows and these variables.

## Results

### A Bayesian machine scientist learns closed-form mathematical models from
mobility data

We aim to determine whether it is possible to model mobility flows by means of
closed-form mathematical models that are interpretable like gravity models,
and as predictive as (non-interpretable) machine learning models such as the
deep gravity model [10]. To automatically learn such closed-form models from
data, we use the so-called Bayesian machine scientist (BMS)[14]. Given a
dataset $D$, the BMS samples closed-form mathematical models from the
posterior distribution $p(M|D)$, which gives the probability that a given
model $M$ is the true generating model given the data (Methods). The BMS is
guaranteed to asymptotically identify the true generating model, if one
exists, and makes quasi-optimal predictions for unobserved data [16].

We consider as our main dataset $D$ the set of flows $T_{od}$ between origin
$o$ and destination $d$ municipalities in six states in the USA (New York,
Massachusetts, California, Florida, Washington, and Texas; see Data). The BMS
is fed with $D$, and samples closed-form models from $p(M|D)$ using Markov
chain Monte Carlo[14] (Methods). This sampling yields an ensemble of hundreds
of different closed-form models for the flows $T_{od}$ such as, for example,

$\log
T_{od}=A\left(1+\frac{B\left(\left(m_{d}+C\right)\left(m_{o}+D\right)\right)^{\beta}}{d_{od}}\right)^{\xi}\quad{\rm
or}\quad\log
T_{od}=\log\left[A\left(\frac{B\left(m_{d}m_{o}+Cm_{d}+D\right)}{d_{od}^{\alpha}}+1\right)^{\gamma}\right]\;,$
(1)

where $m_{o/d}$ is the population of the origin/destination municipalities,
$d_{od}$ is the distance between them, and $A$, $B$, $C$, $D$, $\beta$ and
$\xi$ are model parameters. These models are able to make predictions of test
flows (not seen by the BMS during training) that follow real values over
several orders of magnitude (Figs. 1 and 2M-R). In what follows, we analyze in
more depth this ensemble of models and its predictive abilities, vis a vis
gravity models and machine learning models such as the deep gravity model.

### Different models capture flows at different scales

Figure 2: Model predictions of flows between municipalities. Each panel shows,
in logarithmic scale, the scatter plot of predicted flows between
municipalities versus the corresponding real flows, for different states in
the US (columns). Plots show results for test data for different models
(rows): (A-F) the Deep Gravity model, (G-L) a Random Forest regressor, (M-R)
the most plausible model sampled by the Bayesian machine scientist, and (S-X)
a gravity model in its power law version. Supplementary Fig. S1 shows scatter
plots for the full set of models we consider (see Methods for a complete
description of the models and their parameters).

In order to compare the ability of modeling approaches to describe mobility
data, one needs a model selection criterion. In probabilistic terms, selecting
the best model amounts to selecting the most plausible model, that is, the
model that has the highest probability $p(M|D)$ of being the true generating
model given the observed data; or, equivalently, the model with the shortest
description length (Eqs. (2) and (3)). However, this criterion is not always
applicable in practice, because often it is not possible to compute the
description length of a model, as it happens for deep learning and most other
machine learning models.

Alternatively, one can measure performance at certain predictive tasks[17],
which is the approach typically taken in mobility modeling studies, and that
we take here.

Specifically, for each state for which we have data, we split municipalities
in two sets. Flows between municipalities in the first set comprise the
training set, and flows between municipalities in the second set comprise the
test set. By building the training and test sets in this way, all the
information about the municipalities in the test set, their characteristics,
and the distances between them is completely new to the trained algorithm.

We compare the closed-form mobility models identified by the BMS to two
alternative approaches. First, we consider gravity models, in which mobility
flows are directly proportional to the product of masses (that is,
populations) at the origin and destination, and inversely proportional to the
distance between them. These approaches include traditional gravity models
[1], as well as the closely related radiation model [4]. Second, we consider
machine learning approaches that, besides considering population and distance
between municipalities, also consider additional characteristics of
municipalities such as density of shops, entertainment venues, or educational
facilities (Methods). Specifically we consider a random forest regression
model [18] and the deep gravity model [10].

From the ensemble of closed-form models sampled by the BMS, we analyze
(Methods): (i) the most plausible (minimum description length) model found by
the BMS; (ii) the median of the ensemble of models sampled by the BMS (which
is the optimal predictor); and (iii) what we call the median predictive model,
that is, the single model in the ensemble of sampled models whose predictions
are closest to the ensemble median.

In Fig. 2, we show the predicted flows versus the real flows in the test set
(see Supplementary Fig. S1 for results for additional models). Whereas all
models are predictive, we find that different models are differently capable
to describe mobility flows of different orders of magnitude. For instance,
gravity-like models (Fig. 2 S-X and Supplementary Figs. S1 and S2 ) are
typically good at capturing the behavior of large flows, but not of small
flows. Indeed, for small flows (less than around 100 commuters) these
approaches tend to underestimate flows, in some cases by several orders of
magnitude, and even predict flows smaller than 1 person (Supplementary Fig.
S2). This is also the case for the deep gravity model, which again under-
predicts small flows (Figs. 1 and 2 A-F). By contrast, neither the random
forest nor the BMS suffer from this caveat, and both capture the whole range
of flows more consistently and without large systematic deviations (Figs. 1
and 2 G-R).

### Simple closed-form models are overall more accurate than gravity and non-
interpretable machine learning models

Next, we quantify the performance of the models at the task of predicting
unobserved flows. To that end, and considering the qualitative results in the
previous section, we compute several complementary performance metrics. First,
we consider the common part of commuters (CPC; see Methods), which is a usual
choice in the mobility literature[10]. The CPC measures the overlap between
predicted and observed flows, and can take values from 0 to 1; the larger the
CPC, the better the predictions. Despite its popularity, this metric favors
models that predict the larger flows well, but overlooks errors in small flows
(Fig. 3A-F). Since mobility flows typically span several orders of magnitude
(Fig. 2), models with the larger CPC are not necessarily the best models for
the whole range of flows.

To have metrics of performance that cover the whole range of flows, we
consider, in addition and complementary to CPC, the absolute error, the
absolute relative error, and the absolute log-ratio (Fig. 3). For each of
these metrics, and to avoid the disproportionate influence of singular large
errors (especially for non-relative quantities such as the absolute error), we
always show the whole distribution of error values (as a boxplot), and use the
median value to compare models (Fig. 3); the lower the median, the better the
performance of the model. Note that these metrics highlight different aspects
of the prediction. Absolute errors are correlated with the magnitude of the
flow we are trying to predict, so that errors are typically larger for larger
flows. Because of this, and similar to the CPC, average absolute errors are
very sensitive to the errors in predicting large flows but not to errors in
small flows. For the same reason, median values of the absolute error
typically reflect errors in performance for typical flow values, and do not
reflect the ability of a model to predict values in the whole range of flows.

The absolute relative error and the absolute log-ratio do take into account
the effect of the magnitude of the flow, and therefore are more informative of
the global behavior of a model when the range of flows spans several orders of
magnitude (Fig. 3M-X). An issue with the relative error is that while it
penalizes over-prediction, it does not penalize under-prediction; in the
extreme case in which the predicted flow equals zero and the real flow is
larger than zero, the relative error is equal to one. As a result,
distributions for relative errors in gravity-like models and the deep gravity
model, in which small flows are under predicted, are centered around 1 (Fig.
3S-X; Supplementary Fig. S2). By contrast, the absolute log-ratio has the
property that over- and under-prediction are equally penalized (in a
logarithmic scale), that is, predicting the real flow multiplied or divided by
the same factor results in same absolute log-ratio. This metric therefore
captures the ability of a model to predict flows in any range of values.

Figure 3: Model performance at predicting flows between municipalities. For
each one of the model predictions shown in Fig. 2 and Supplementary Fig. S1 we
assess model performance using four different metrics:A-F, Common part of
commuters; G-L Absolute error; M-R, Absolute relative error; S-X, Absolute
log-ratio. Each row corresponds to a different US state as indicated. The
Common Part of Commuters is a global metric, thus we have a single value for
each metric. For the other three metrics, we show the median, 50% confidence
interval (box) and 95% confidence interval (whiskers).Triangles
($\blacktriangleleft$) indicate the best performing model for each metric
(largest CPC or lowest median). See Methods and text for the definition and
discussion of the different metrics. See Supplementary Table S3 for numerical
values. Figure 4: Performance for different flow ranges. The observed flows
between municipalities span six orders of magnitude and are distributed, over
the six states we consider, as shown by the histogram. We measure the
performance of the different models at predicting flows within each bin of
flows: A, CPC; B, absolute error; C, absolute relative error; D, absolute log
ratio.

Using these metrics, we compare the different modeling approaches (Fig. 3).
The first conclusion from this comparison is that gravity models, including
the radiation model, are never optimal; for all states and metrics we
consider, there is always at least one other model that performs better. This
is not surprising, since these models are simple and highly stylized, and they
have already been shown to make less accurate predictions than deep gravity
models[10].

Perhaps more surprisingly, we find that closed-form mathematical models
obtained by the BMS perform, overall, at least as well or even slightly better
than machine learning models, even when the latter are much more complex and
can use many features other than population and distance. For the CPC, machine
learning models are the best performing model for three out of six states (the
random forest is best in two states, and deep gravity in one), whereas closed-
form models identified as most plausible by the BMS are optimal in the three
remaining states. In the case of California, the BMS overestimates a few
extremely large flows (hence the CPC$\simeq 0$), but is able to describe all
other flows well.

Similarly, for the relative error, random forest models are optimal in half of
the States, whereas closed-form BMS models are optimal in the other half. When
models are compared in terms of both absolute error and log-ratio, we find
that closed-form BMS models are optimal in four of the States, whereas random
forest models are optimal in the other two. Remarkably, deep gravity models
are never optimal according to absolute error, relative error or log ratio.

When considering how each model performs for flows in specific ranges (Fig. 4
and Supplementary Figure S3), we find that BMS models, similar to random
forest models, are particularly good at modeling flows in the range that is
the most common in the data.

Taking into account that both random forest and deep gravity models use many
more features for their predictions (39 features in total, in contrast to the
three features used by gravity models and the closed-form models identified by
the BMS, namely, origin and destination population, and origin-destination
distance), we conclude that the symbolic regression approach using the BMS
yields more parsimonious models of human mobility flows between
municipalities. Closed-form models obtained by the BMS also compare well to
the alternatives in terms of the fairness of their predictions[19]
(Supplementary Fig. S4). In particular, we find that the random forest and the
models identified by the BMS are, overall, the most consistent models across
states in terms of the fairness of their predictions.

Figure 5: Model performance at predicting flows at small distances. We
consider the data used by Simini et al.[10] about mobility flows between
census tracts within small geographical (25 km2) regions in NY state. We
evaluate predictions over test data using the same metrics as in Fig. 3.
Different columns correspond to different metrics and triangles
($\blacktriangleleft$) indicate the best performing model according to each
metric. See Methods and text, for a complete description of the metrics, the
models, and the data.

### Closed-form models also describe flows at shorter scales

So far, we have analyzed within-State flows between municipalities, at any
range of distances and of any size. However, it may be that the geographical,
economic and demographic characteristics of smaller areas become relevant when
modeling flows at shorter distances (for example, within neighborhoods or
adjacent towns in large metropolitan areas). To elucidate to what extent
different modeling approaches can accommodate to such short-distance flows, we
adopt the framework used in Ref. [10]: we divide the state of New York in
small tiles of 25 km2, and consider the flows between census tracts within
each tile.[10] We use 50% of the tiles for training the different models, and
then test on the flows between the remaining 50% of the tiles. For this
experiment we find that, regardless of the metric used, closed-form models
identified by the BMS are always more accurate than machine learning and
gravity models (Fig. 5). Our results thus indicate that simple closed-form
models that just consider populations of municipalities/census tracts and the
distances between them provide better descriptions of mobility flows than
complex models that take many more features into account, and have many more
parameters, also for flows at short distances.

### The Bayesian machine scientist finds gravity-like models to describe
mobility flows

Our analysis indicates that the BMS is able to find closed-form mathematical
models that solely consider the populations of the origin and destination and
the geographic distance between them; and that these models provide
predictions of mobility flows that are even more accurate than complex models
such as the deep gravity model. Here, we investigate whether, besides being
predictive, these closed-form models are also interpretable and insightful.

We start by noting, once more, that the BMS samples hundreds of models and
that, overall, they all perform well, which shows that there are many
different models that can describe the data. Such a set of models is sometimes
called a Rashomon set[13]. Then, the relevant question is whether these models
share any common defining properties that could explain why they describe
mobility flows accurately.

To elucidate this question, we analyze two particularly relevant closed-form
models identified by the BMS (Fig. 6): (i) the most plausible model, that is,
the model that has the highest probability $p(M|D)$ given the data (or,
equivalently, the shortest description length ${\mathcal{L}}(M,D)$) among all
those sampled by the BMS (Methods); (ii) the median predictive model, that is,
the closed-form model whose predictions for unobserved data are closest to the
median prediction of the whole ensemble of sampled models (Methods). Formally,
there are marked mathematical differences between both models. In particular,
the most plausible model is an exponential model for the flows, while the
median predictive modesl is a power-law model for the flows. However, the two
models have relevant properties in common. First, both models are gravity-like
models, that is, they depend on a product of the origin and destination
populations (shifted by a certain amount), and inversely on a function of the
geographic distance between origin and destination. This is remarkable because
the BMS has not received any input about the particular shape that models
should take, which suggests that the regularities in the data are well-
described by this general class of models and justifies the historical use of
gravity models.

Second, we find that origin and destination do not necessarily play a
symmetric role, which allows the model to accommodate non-symmetric flows in
contrast to typical gravity models which do not allow for this possibility.
Indeed, an inspection of the parameters shows that in some states such as New
York or Texas, flows are much more symmetric than in others such as Florida
and Massachusetts.

Finally, we also find that the relative contribution of the mass product with
respect to the geographical distance is consistent across models. In both
models, we find a mathematically equivalent dependency on the ratio
$(m_{d}m_{o})/d^{e}$, where $e=1/\beta$ in the most plausible model (Fig. 6A),
and $e=\alpha$ in the median predictive model (Fig. 6B). We find that, for a
given state, $\alpha\approx 1/\beta$ suggesting that this relationship is to a
large extent model-independent (Fig. 6C). Furthermore, we find that the state-
to-state variability is small since all exponents fall within the range
$[1,2]$ (most within the range $[1.5,2]$), which suggests that reasonable
models for mobility flows are gravity-like models with specific constraints in
the relationship between the contributions of the mass product and the
geographical distance.

Figure 6: Closed-form models for mobility flows. We ran the Bayesian machine
scientist (BMS) with a training set of 1000 points and three features: origin
and destination populations and distance between them. We used 5 independent
Markov chains of 12,000 Monte Carlo steps each. (A) Minimum description length
model (Methods) for the logarithm of the data, where $d$ is the inter-
municipality distance, $m_{o}$ is the origin population, and $m_{d}$ is the
destination population. In the table we show the fitting parameters for each
state of the training data. (B) Median predictive model (Methods) for the
logarithm of the data. As before, $d$ is the inter-municipality distance,
$m_{o}$ is the origin population, and $m_{d}$ is the destination population.
In the table we show the fitting parameters for each state of the training
data. (C) Ratio between the distance exponent and the population exponent.
Blue points are obtained from the most plausible model and orange points from
the median predictive model.

## Discussion

Understanding human mobility is critical to address questions in urban
planning and transportation, as well as global challenges in sustainability,
public health, and economic development. Traditionally, mobility flows have
been modeled using simple gravity models, which are conceptually simple and
easy to interpret, but have limited predictive power. Recently, deep learning
models have been proposed as an alternative; whereas these models are
consistently and significantly more predictive than gravity models, they are
not interpretable and provide little insight into human behavior. Here, we
have shown that automated equation discovery approaches lead to parsimonious
models that combine the most desirable aspects of both approaches—the
simplicity of gravity models and a predictive power that is even better,
overall, than deep learning models. Remarkably, the models we identified are
gravity-like in that they are increasing functions of a certain product of
populations of origin and destination, and decreasing functions of the
distance between them. While the ratio between the population and the distance
terms is, in principle, model and data-dependent, in the datasets we explore
the ratio is roughly constant and dataset-independent.

Individual mobility depends critically on the urban environment, personal
preferences, commuting patterns, and accessibility to transportation and
amenities. Thus, modeling mobility at an individual or small spatial scale
might require more complicated models that account for routes, the purpose of
the trip, points of interest, or even the demographic traits of individuals
[19, 10, 20]. However, our results show that by aggregating mobility at a
larger spatial scale, the movements of millions of people can be described by
simple fully explainable models that do not depend on the microscopic
characteristics of the origin, destination, or route taken. This is because
the randomness and variability inherent in individual behaviors tend to cancel
out when looked at collectively, revealing underlying trends and movements
that are driven by the shared needs of large populations, and the structure of
the built environment. Our results show, therefore, that the aggregated flows
in human mobility can be seen as an emergent and universal property of the
complex system of individual movements. More broadly, our results showcase the
potential of using machine scientists[21, 22, 14] to automate the process of
finding similar phenomenological closed-form models from data; and to use
these models to gain insight into the relevant variables and mechanisms to
describe other complex phenomena.

## Methods

### Mobility data between municipalities in the US

We collected weekly flows between United States census tracts [23] for the
first week of January 2019 (2019/01/07 to 2019/01/13) and the first week of
March 2019 (2019/03/04 to 2019/03/9). The data consist of anonymous mobile
data trajectories, regardless of the transportation method. The data set
contains the geographical identifier (GEO ID) for both origin and destination
census tracts as well as their corresponding geographical coordinates, the
estimated number of visitors detected by SafeGraph, and the estimated
population flows inferred from the number of visitors. The datasets are
available at GeoDS.

#### Data processing

In order to obtain the flows between municipalities from the data using census
tract data, we first match municipalities (cities, towns and villages) with
their corresponding census tracts. Then, the total flow between two
municipalities A and B is calculated as the aggregate of flows between the set
of census tracts municipality A is comprised of and the set municipality B is
comprised of.

Specifically, we consider mobility data sets within six states in the US: New
York, Massachusetts, California, Florida, Washington, and Texas. For each
state, our data consist of the origin-destination municipality names, the
flow, the distance, and the origin-destination populations and POI categories.
We only consider municipalities with a non-zero population and pairs of
municipalities with non-zero flow; we do not consider flows within the same
municipality (See Supplementary Table LABEL:t.summary for details).

#### Information about municipalities and census tracts

We obtained shapefiles containing the geographical coordinates of the polygons
delimiting census tracts and municipalities from United States Census Bureau.
We also retrieved population data of each municipality using the GEO ID Data
Commons.

Finally, using a local copy of Open Street Map (OSM) with Overpass API, we
retrieved information about Points of Interest (POI) in each census tract. We
selected 18 categories of OSM elements that represent geographical,
demographic and socio-economic features of the different municipalities.

#### Construction of train and test datasets

To generate the train and test datasets for each state, we divide
municipalities in two approximately equal sets at random. The training fold
comprises flows between municipalities in the first set; the test fold
comprises flows between municipalities in the second set. By doing so, we
ensure that we test the learned models on data never observed before (not only
in terms of flows, but in terms of municipality features as well). Tables 1
and 2 show the characteristics of train and test datasets. To speed up the
training process, we train the models with the same random sample of 1000
points of the training fold except for the Deep Gravity model for which we
have to use a lager number of data points for training.

State | Entries | Municipalities | Flow | Distance (Km) | Population
---|---|---|---|---|---
|  |  | Min | Max | Min | Max | Min | Max
New York | 1000 | 217 | 26 | 122991 | 1.83 | 489.80 | 185 | $8.80\cdot 10^{6}$
Massachusetts | 1000 | 75 | 19 | 66946 | 2.34 | 205.51 | 1029 | 675647
California | 1000 | 260 | 14 | 91267 | 2.12 | 1078.21 | 237 | $1.013\cdot 10^{6}$
Florida | 1000 | 223 | 7 | 78316 | 2.54 | 607.50 | 251 | 949611
Washington | 1000 | 107 | 20 | 71796 | 2.17 | 467.16 | 20 | 228989
Texas | 1000 | 177 | 21 | 192660 | 1.84 | 975.16 | 173 | $2.30\cdot 10^{6}$
Table 1: Train dataset. Number of points and municipalities in the train set for each state, obtained form a random sample of 1000 points of the original train fold. We also detail the lowest and largest flow, distance and municipality population. State | Entries | Municipalities | Flow | Distance (Km) | Population
---|---|---|---|---|---
|  |  | Min | Max | Min | Max | Min | Max
New York | 5952 | 249 | 20 | 35063 | 0.77 | 531.39 | 361 | 211569
Massachusetts | 1180 | 75 | 9 | 55499 | 2.16 | 267.24 | 1517 | 206518
California | 11727 | 319 | 2 | 416696 | 1.09 | 1084.34 | 129 | $3.9\cdot 10^{6}$
Florida | 7092 | 245 | 3 | 79624 | 1.14 | 875.82 | 78 | 442241
Washington | 2083 | 109 | 12 | 54245 | 1.48 | 431.53 | 487 | 737015
Texas | 2782 | 192 | 3 | 112865 | 2.46 | 1196.33 | 106 | $1.43\cdot 10^{6}$
Table 2: Test dataset. Number of points and municipalities in the test set for
each state. We also detail the lowest and largest flow, distance and
municipality population.

#### Economic classification of municipalities

We retrieve the median income per capita of each municipality and state as of
2020 from Data Commons. Then we classify each municipality with the label
_rich_ if the median income of the municipality is above the median income of
the state, or _poor_ if the median income per capita of the municipality is
below the median income of the state.

### Mobility data at small scales from Simini et al. [10].

The code available for the Deep Gravity model provides data at the level of
census tracts for the New York State area. Mobility flows are obtained from
the same dataset. In order to obtain the predictions of the Deep Gravity model
for each individual trip and also the same dataset, we run the program and
save, for each trip, the corresponding tile, the real value, the prediction
and the variables. We store the train and test sets for the comparison with
other methods.

### Bayesian Machine Scientist

The Bayesian machine scientist (BMS) is a Bayesian approach to symbolic
regression that estimates the plausibility of a closed-form mathematical model
$M$ given the observed data $D$ as the posterior probability $p(M|D)$. Without
loss of generality, this posterior can be written [14]

$p(M|D)=\frac{\exp\left[-{\mathcal{L}}(M,D)\right]}{Z}\,,$ (2)

where ${\mathcal{L}}(M,D)$ is the description length[24] of the model (and the
data), and
$Z=\sum_{M^{\prime}}\exp\left[-{\mathcal{L}}(M^{\prime},D)\right]=p(D)$ is the
evidence. Within standard approximations[25, 14], the description length can
be computed as

${\mathcal{L}}(M,D)=\frac{B(M,D)}{2}-\log p(M)\,,$ (3)

where $B(M,D)$ is the Bayesian information criterion[25], which is
straightforward to calculate from the data, and $p(M)$ is a suitable prior
distribution over models[14].

In this work, our goal is to simultaneously model flows from six different
states in the US (that is, six different datasets). To this end, we use a
multi-dataset approach[15] which consists in finding a unique closed-form
model for multiple datasets. For each dataset, we allow model parameters to
take different values[15]. For a single dataset $D=\\{(y_{i},{\bf x}_{i})\\}$,
where $\\{y_{i}\\}$ is the set of observations and $\\{{\mathbb{x}}_{i}\\}$ is
the set of feature values associated to each observation, the description
length of a closed-form model $M$ and the data is given by Eq. (3). In the
case in which our data comprises $K$ independent datasets
$D=\\{D_{k},\leavevmode\nobreak\
k=1\dots,K\\}=\left\\{\\{(y^{1}_{i},{\mathbb{x}}^{1}_{i})\\},\dots,\\{(y^{K}_{i},{\mathbb{x}}^{K}_{i})\\}\right\\}$
that we want to model using a single model $M$, the description length is[15]

${\mathcal{L}}(M,D)=\frac{1}{2}\sum_{k}B(M,D_{k})-\log
p(M)\leavevmode\nobreak\ .$ (4)

The BMS represents closed-form models as labeled trees and uses Markov chain
Monte Carlo (MCMC) to explore the space of closed-form mathematical models by
sampling from the posterior distribution $p(M|D)\propto\exp(-{\mathcal{L}})$.
In this work, we consider models for flows with three independent variables
(populations at origin and destination and geographical distance) and up to
six parameters.

The ensemble of sampled models allows to make predictions on the test set
using three different approaches:

  1. 1.

The most plausible model. This is the model with the shortest description
length that the BMS is able to find.

  2. 2.

The ensemble of models. Ensemble predictions are the optimal predictions since
they correspond to fully integrating over model space. We estimate this
integral by averaging over the predictions made by each one for the models we
sample. Specifically, we perform five independent realizations of 12,000 MCMC
steps. We then collect a model every 100 steps within the last 2,000 steps of
the Markov chain to obtain an ensemble of 100 models.

  3. 3.

The median predictive model. This is the model within the ensemble whose
predictions are closest to the predictions of the ensemble as a whole.

### Benchmark models

#### Gravity model.

We consider the gravity model in its simplest form [1] in which the observed
flow $T_{ij}$ between municipalities $(i,j)$ is a function of the populations
$m_{i}$ and $m_{j}$ of the municipalities and the distance $d_{ij}$ between
them

$T_{ij}=C\frac{m_{i}\,m_{j}}{f(d_{ij})}\,.$ (5)

Here $C$ is a scaling parameter and $f(d)$ is a function of the distance.
Specifically, we consider two possible choices for $f(d)$: i) a power-law
$f_{\rm pow}(d)=d^{\alpha}$; and ii) an exponential-law $f_{\rm
exp}(d)=\exp(\alpha d)$. In both cases, the parameter $\alpha$ is obtained by
fitting the model to the data in the training set. Because flows span several
orders of magnitude, we find that training the model on the logarithm of the
flows, rather than the flows themselves, leads to more predictive models.
Therefore, all results reported here for gravity models correspond to this
approach.

#### Radiation model.

We consider the original formulation of the model [4], in which flow $T_{ij}$
is modeled as the the total outflow of an origin municipality $T_{i}$ times
the probability of going from $i$ to $j$. This probability depends on the
populations of the origin ($m_{i}$) and destination($m_{j}$), as well as the
populations of the municipalities within a radius $d_{ij}$ from the
municipality at the origin:

$T_{ij}=T_{i}\,p_{i\to
j}=T_{i}\frac{m_{i}\,m_{j}}{(m_{i}+s_{ij})\,(m_{i}+m_{j}+s_{ij})},$ (6)

where $s_{ij}=\sum_{k\neq i,j}m_{k}(\forall k\;d_{ik}<d_{ij})$.

Recent works introduce modifications to this model for finite-size systems and
in order to avoid border effects [26, 27, 4, 28, 29](See Supplementary Fig.
S5). However, we find that these models do not outperform the original
formulation .

#### Random Forest.

We implement a Random Forest regressor [30] with 1,000 estimators. As input
data we use a total of 39 features for each origin-destination pair which
include: distance, population at origin and destination, and 36 geographical
and socio-economic features of the origin and destination areas (see Data and
Supplementary Table LABEL:t.poi).

#### Deep Gravity.

Taking as a baseline the algorithm developed in [10], we modified the model to
predict flows between municipalities rather than between small geographic
regions resulting from tessellation. The major difference with the original
model is that municipalities are now the smallest geographic unit, allowing us
to compare with the models evaluated in this study. In all other aspects
(features used, pre-processing of data, an model training) the model remains
the same. The modified version of the code can be consulted and downloaded
here https://github.com/ocabanas/Symbolic_mobility_BMS/DeepGravity.

### Metrics

#### Common Part of Commuters (CPC).

It is a widely used metric to analyze the performance of mobility models that
is defined as

$\mathsf{CPC}=\frac{2\sum_{ij}{\rm
min}(T_{ij},T_{ij}^{*})}{\sum_{ij}T_{ij}+T_{ij}^{*}}\;,$ (7)

where $T_{ij}$ is the predicted value of the flow from $i$ to $j$ and
$T^{*}_{ij}$ is the observed flow. The maximum value of the CPC is $1$ if
there is a complete agreement between the real data and predictions and it
decreases to $0$ if all predictions for any flow are equal to zero. Note that
the CPC is biased toward models that make accurate predictions for large
flows, since smaller flows have marginal contributions to the sums. This is
especially critical in mobility data where flows can span several orders of
magnitude (see Tables 1, 2).

####  Absolute error.

It measures the distance between the real and predicted data

$E_{ij}=|T_{ij}-T_{ij}^{*}|\leavevmode\nobreak\ \leavevmode\nobreak\ .$ (8)

Note that absolute error scales with the size of the flows, so that the
average absolute error is biased towards the errors of average flows. For this
reason we represent the whole distribution.

#### Relative error.

It measures the difference between real and predicted flows relative to the
observed value of the flow:

$\epsilon_{ij}=|\frac{T_{ij}-T_{ij}^{*}}{T_{ij}^{*}}|\leavevmode\nobreak\ .$
(9)

Note that, while over-predicting flows is penalized by the relative error,
under-predicting flows is not, since predicting a zero value for a non-zero
flow results in $\epsilon_{ij}=1$. Because this can bias average relative
errors, we plot the whole distribution.

#### Absolute log-ratio.

It measures the difference in the logarithms of predicted and real flows:

$LR_{ij}=|\log\frac{T_{ij}}{T_{ij}^{*}}|\leavevmode\nobreak\ .$ (10)

Note that for a perfect prediction this metric is equal to zero. Importantly
this metric penalizes both over- and under-predictions equally. For instance,
a prediction of a flow twice as large $T_{ij}=2T^{*}_{ij}$ has
$LR_{ij}=\log\leavevmode\nobreak\ {2}$, and a prediction $T_{ij}=T^{*}_{ij}/2$
has $LR_{ij}=\log\leavevmode\nobreak\ {2}$ as well.

#### Proportional Demographic Parity (PDP).

The goal of this metric is to quantify the fairness of a model when predicting
flows between different demographic or socio-economic groups
$\\{g\in\mathcal{G}\\}$. To do so, it quantifies to what extent the errors of
the predictions for flows across pairs of groups
$\\{f_{ij}\equiv(g_{i},g_{j})\in\mathcal{G}^{2}\equiv\mathcal{G}\times\mathcal{G}\\}$
are equally distributed for all pairs of (different) groups. Consider that
$\bar{l}$ is the median error of all flows, and $\tau$ is a percentile window
around the median ($0\leq\tau\leq 100$). For a pair of flow groups
$(f_{1},f_{2})$, $\mathsf{PDP}$ estimates the difference between their error
distributions as

$\mathsf{PDP}_{f_{1},f_{2}}=\left|P(\overline{l}-\frac{\tau}{2}\leq
l\leq\overline{l}+\frac{\tau}{2}|f_{1})-P(\overline{l}-\frac{\tau}{2}\leq
l\leq\overline{l}+\frac{\tau}{2}|f_{2})\right|$ (11)

where $P(\cdot|f_{i})$ is the probability that a prediction of a flow in flow
group $f_{i}$ has an error $l$ such that $\overline{l}-\frac{\tau}{2}\leq
l\leq\overline{l}+\frac{\tau}{2}$.

To get an overall estimate, then $\mathsf{PDP}$ uses a weighted average[19]

$\mathsf{PDP}=\sum_{f,h\in\mathcal{G}^{2},f\neq
h}w_{f,h}\,\mathsf{PDP}_{f,h}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ ,$ (12)

where the weight
$w_{f,h}=\sum_{k\in\mathcal{G}^{2}}N_{k}/\left(N_{f}+N_{h}\right)$ enhances
the relative contribution of small groups of flows. Note that our approach to
measure $\mathsf{PDP}$ is a generalization of that used in Ref.[19], where,
instead of percentile windows, the authors consider $\tau$ to be a standard
deviation around the mean. However, because error distributions are not
Gaussian in general (see explanation for the different error metrics), we use
a more general definition applicable to any distribution.

In our analysis, we consider two groups of municipalities: above (rich) and
below (poor) the median income per capita (see Data). Therefore we have four
different flow groups $\mathcal{G}^{2}=\\{{\rm poor\to poor},\,{\rm rich\to
poor},\,{\rm poor\to rich},\,{\rm rich\to rich}\\}$.

## References

  * [1] Zipf, G. K.  The p1 p2/d hypothesis: On the intercity movement of persons.  _American Sociological Review_ 11, 677–686 (1946).
  * [2] Erlander, S. & Stewart, N. F.  _The gravity model in transportation analysis: theory and extensions_ , vol. 3 (Vsp, 1990).
  * [3] Guimerà, R., Mossa, S., Turtschi, A. & Amaral, L. A. N.  The worldwide air transportation network: Anomalous centrality, community structure, and cities’ global roles.  _Proc. Natl. Acad. Sci. USA_ 102, 7794–7799, 10.1073/pnas.0407994102 (2005).
  * [4] Simini, F., González, M. C., Maritan, A. & Barabási, A.-L.  A universal model for mobility and migration patterns.  _Nature_ 484, 96–100, 10.1038/nature10856 (2012).
  * [5] Schläpfer, M. _et al._ The universal visitation law of human mobility.  _Nature_ 593, 522–527, 10.1038/s41586-021-03480-9 (2021).
  * [6] Yuan, H. & Li, G.  A survey of traffic prediction: from spatio-temporal data to intelligent transportation.  _Data Science and Engineering_ 6, 63–85 (2021).
  * [7] Haynes, K. E. & Fotheringham, A. S.  _Gravity and Spatial Interaction Models_.  No. 07 in Wholbk (Regional Research Institute, West Virginia University, 1985).
  * [8] Moro, E., Calacci, D., Dong, X. & Pentland, A.  Mobility patterns are associated with experienced income segregation in large us cities.  _Nature Communications_ 12, 4633, 10.1038/s41467-021-24899-8 (2021).
  * [9] Balcan, D. _et al._ Modeling the spatial spread of infectious diseases: The global epidemic and mobility computational model.  _Journal of computational science_ 1, 132–145 (2010).
  * [10] Simini, F., Barlacchi, G., Luca, M. & Pappalardo, L.  A Deep Gravity model for mobility flows generation.  _Nature Communications_ 12, 6576, 10.1038/s41467-021-26752-4 (2021).
  * [11] Pappalardo, L., Rinzivillo, S. & Simini, F.  Human mobility modelling: Exploration and preferential return meet the gravity model.  _Procedia Computer Science_ 83, 934–939, https://doi.org/10.1016/j.procs.2016.04.188 (2016).  The 7th International Conference on Ambient Systems, Networks and Technologies (ANT 2016) / The 6th International Conference on Sustainable Energy Information Technology (SEIT-2016) / Affiliated Workshops.
  * [12] Chen, Y.  The distance-decay function of geographical gravity model: Power law or exponential law?  _Chaos, Solitons & Fractals_ 77, 174–189, https://doi.org/10.1016/j.chaos.2015.05.022 (2015).
  * [13] Rudin, C.  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  _Nat. Mach. Intell._ 1, 206–215 (2019).
  * [14] Guimerà, R. _et al._ A Bayesian machine scientist to aid in the solution of challenging scientific problems.  _Science Advances_ 6, eaav6971, 10.1126/sciadv.aav6971 (2020).  https://www.science.org/doi/pdf/10.1126/sciadv.aav6971.
  * [15] Reichardt, I., Pallarès, J., Sales-Pardo, M. & Guimerà, R.  Bayesian machine scientist to compare data collapses for the Nikuradse dataset.  _Phys. Rev. Lett._ 124, 084503, 10.1103/PhysRevLett.124.084503 (2020).
  * [16] Fajardo-Fontiveros, O. _et al._ Fundamental limits to learning closed-form mathematical models from data.  _Nature Communications_ 14, 1043, 10.1038/s41467-023-36657-z (2023).
  * [17] Vallès-Català, T., Peixoto, T. P., Sales-Pardo, M. & Guimerà, R.  Consistencies and inconsistencies between model selection and link prediction in networks.  _Phys. Rev. E_ 97, 062316, 10.1103/PhysRevE.97.062316 (2018).
  * [18] Ho, T. K.  Random decision forests.  In _Proceedings of 3rd international conference on document analysis and recognition_ , vol. 1, 278–282 (IEEE, 1995).
  * [19] Liu, Z., Huang, L., Fan, C. & Mostafavi, A.  Fairmobi-net: A fairness-aware deep learning model for urban mobility flow generation (2023).  2307.11214.
  * [20] Feng, J. _et al._ DeepMove: Predicting Human Mobility with Attentional Recurrent Networks.  _Proceedings of the 2018 World Wide Web Conference on World Wide Web - WWW ’18_ 1459–1468, 10.1145/3178876.3186058 (2018).
  * [21] Džeroski, S. & Todorovski, L. (eds.) _Computational Discovery of Scientific Knowledge_.  Lecture Notes in Artificial Intelligence (Springer, 2007).
  * [22] Evans, J. & Rzhetsky, A.  Machine science.  _Science_ 329, 399–400 (2010).
  * [23] Kang, Y., Gao, S., Liang, Y., Li, M. & Kruse, J.  Multiscale dynamic human mobility flow dataset in the U.S. during the COVID-19 epidemic.  _Scientific Data_ 1–13 (2020).
  * [24] Grünwald, P. D.  _The Minimum Description Length Principle_ , vol. 1 of _MIT Press Books_ (The MIT Press, 2007).
  * [25] Schwarz, G.  Estimating the Dimension of a Model.  _The Annals of Statistics_ 6, 461 – 464, 10.1214/aos/1176344136 (1978).
  * [26] Masucci, A. P., Serras, J., Johansson, A. & Batty, M.  Gravity versus radiation models: On the importance of scale and heterogeneity in commuting flows.  _Phys. Rev. E_ 88, 022812, 10.1103/PhysRevE.88.022812 (2013).
  * [27] Yang, Y., Herrera, C., Eagle, N. & González, M. C.  Limits of predictability in commuting flows in the absence of data for calibration.  _Scientific Reports_ 4, 5662, 10.1038/srep05662 (2014).
  * [28] Lenormand, M., Huet, S., Gargiulo, F. & Deffuant, G.  A universal model of commuting networks.  _PLOS ONE_ 7, 1–7, 10.1371/journal.pone.0045985 (2012).
  * [29] Lenormand, M., Bassolas, A. & Ramasco, J. J.  Systematic comparison of trip distribution laws and models.  _Journal of Transport Geography_ 51, 158–169, https://doi.org/10.1016/j.jtrangeo.2015.12.008 (2016).
  * [30] Pedregosa, F. _et al._ Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ 12, 2825–2830 (2011).

## Acknowledgements

We thank L. Pappalardo and M. Luca for help with the Deep Gravity algorithm,
and M. Luca for comments and suggestions on the manuscript. This research was
supported by projects PID2019-106811GB-C32 (E.M), PID2019-106811GB-C31 and
PID2022-142600NB-I00 (M.SP. and R.G.), and FPI grant PRE2020-095552 (O.CT.)
from MCIN/AEI/10.13039/501100011033, and by the Government of Catalonia
(2021SGR-633) (M.SP. and R.G.). E.M. acknowledges support from the National
Science Foundation under Grant No. 2218748.

## Author contributions statement

O.C. collected data. O.C. and L.D. wrote code and performed experiments. All
authors designed research, analyzed results, discussed results, and wrote the
paper.

### Data and code availability

All data are available as described in the Methods section. Data used for the
evaluation of the Deep Gravity model are available at Github. The code for the
BMS is available from https://bitbucket.org/rguimera/machine-scientist.

## Competing interests

The authors declare no conflict of interest.


In particle simulations, often the dynamics results from the summation of all
pair-wise forces in the ensemble of particles. Such situations arise in
astrophysics, molecular dynamics, plasma physics, and certain formulations of
fluid dynamics problems, for example. The total field of interest
(gravitational, electrostatic, etc.) at one evaluation point requires adding
the contribution of all source points or particles, and so if both evaluation
points and particles are numbered at $N$, a total of $N^{2}$ operations is
needed. This fact was for a long time an impediment to the wider use of
particle simulations, as the computational effort becomes prohibitive for
large numbers of particles.

The above scenario changed dramatically with the introduction of tree-codes
and the fast multipole method (FMM), which appear in the late 1980s for
accelerated evaluation of $N$-body problems. Tree-codes
Appel1985,BarnesHut1986 are generally perceived to be easier to grasp and
program, and provide a complexity of $\mathcal{O}(N\log N)$. The FMM was
introduced as an algorithm for the rapid evaluation of gravitational or
Coulombic interactions greengard+rokhlin1987 and promises a reduction in
computational complexity to $\mathcal{O}(N)$. It has, since its dissemination,
been adapted for many applications: for fast evaluation of boundary elements
Gaspar1998, for vortex sheet problems with desingularized equations
HamiltonMajda1995, for long-range electrostatics in DNA simulations
FenleyETal1996, and many others. The impact of the FMM has been undeniably
huge, resulting in it being chosen as one of the Top 10 Algorithms of the 20th
Century DongarraSulli2000.

Despite the great volume of work using and adapting the FMM in many
application areas, there remains some lack of insight regarding how the
algorithm can be efficiently used to obtain an accurate representation of the
field of interest. The error of the FMM approximation is estimated by
theoretical bounds, which as could be expected reveal a trade-off between
accuracy and efficiency of the computation. However, there is not much
literature providing measurements of the accuracy of the approximation, in
practice. One may often find such assertions in published works as “only the
first three moments of the expansion were used”, or something to that effect.
But just as often there is no information provided about the actual errors
which are observed. Of course, it is not easy to provide such measures of
observed error, as this would require additional computations using the direct
$\mathcal{O}(N^{2})$ method, for comparison purposes. Nevertheless, it is
important for users of the algorithm to know what to expect in terms of
accuracy and efficiency, depending on the choice of algorithm parameters.

We aim to contribute to this gap in understanding by presenting a methodical
investigation into the errors of the approximation used by the FMM, when the
underlying ‘client’ application is the calculation of the velocity field
induced by $N$ regularized particles of vorticity. This application is rather
more demanding than the Newtonian force calculation, because in the latter
case the gravitational interaction is dominated by the first moment —due to
the fact that all mass is positive. Therefore, keeping only the first three
moments could easily give the desired accuracy. On the other hand, as in
Coulomb electrostatic calculations, the vortex particles can be both positive
and negative, and thus an acceptable accuracy may require that more terms in
the expansion be kept.

For the purposes of this study, a prototype code of the FMM computation of the
velocity induced by $N$ vortex particles was implemented using the
Python††http://www.python.org/ language. The nice features of Python —such as
dynamic typing, extensive numerical libraries, and high programmer
productivity— helped us produce a piece of software which is easy to use and
easy to understand. We are currently using this Python code as a starting
point for a parallel version of the code which, in collaboration with members
of the developer team, will be incorporated to the PETSc library for
scientific computingpetsc-manual. This project will be reported elsewhere, but
preliminary results are being presented in the upcoming Parallel CFD
meetingCruzBarbaKnepley2008a. Our final aim is to contribute to the community
of particle simulations with an open source FMM implementation which is
parallel and portable. For the time being, the Python code is being made
available publicly and we welcome correspondence from interested
researchers††http://www.maths.bris.ac.uk/ aelab/research/pyFMM.html .

Using the Python FMM code, more than 900 calculations were performed, varying
the numerical parameters: $N$, the number of particles, $l$, the number of
levels in the tree, and $p$, the truncation level of the multipole expansion.
We looked not only at the maximum error in the domain, which would be the
conventional approach; we also present results showing how the error varies in
space, revealing some interesting features of the method. Through this
presentation of the results, we believe a clear characterization of the nature
of the FMM approximation is obtained.

The paper is organized as follows. The next section presents an outline of the
vortex particle method, for completeness. In §fmm, we offer an overview of the
FMM, with some details of our implementation. Following, in §errors, we
discuss briefly the sources of errors in the FMM algorithm. And finally,
§results reports the detailed experiments using the FMM for evaluation of the
velocity of $N$ vortex particles; the behavior of the method will be
illustrated for varying parameters, as well as the impact on the efficiency of
the calculation, for different problem sizes.


* Agarwal et al. (2017) Agarwal, Sumit, Itzhak Ben-David, and Vincent Yao, 2017, Systematic mistakes in the mortgage market and lack of financial sophistication, Journal of Financial Economics 123, 42–58.
  * Agarwal et al. (2014) Agarwal, Sumit, Souphala Chomsisengphet, Neale Mahoney, and Johannes Stroebel, 2014, Regulating consumer financial products: Evidence from credit cards, The Quarterly Journal of Economics 130, 111–164.
  * Alan and Loranth (2013) Alan, Sule, and Gyongyi Loranth, 2013, Subprime consumer credit demand: Evidence from a lender’s pricing experiment, The Review of Financial Studies 26, 2353–2374.
  * Allcott et al. (2021) Allcott, Hunt, Joshua Kim, Dmitry Taubinsky, and Jonathan Zinman, 2021, Are high-interest loans predatory? Theory and evidence from payday lending, The Review of Economic Studies 89, 1041–1084.
  * Ally (2017) Ally, 2017, Ally Auto Receivables Trust, Prospectus 2017-3, Ally Auto Assets LLC.
  * Ally (2019) Ally, 2019, Ally Auto Receivables Trust, Prospectus 2019-3, Ally Auto Assets LLC.
  * Andersen et al. (1993) Andersen, Per Kragh, Ørnulf Borgan, Richard D. Gill, and Niels Keiding, 1993, Statistical Models Based on Counting Processes (Springer).
  * Assunção et al. (2013) Assunção, Juliano J., Efraim Benmelech, and Fernando S. S. Silva, 2013, Repossession and the democratization of credit, The Review of Financial Studies 27, 2661–2689.
  * Ausubel (1991) Ausubel, Lawrence M., 1991, The failure of competition in the credit card market, The American Economic Review 81, 50–81.
  * Ayres and Siegelman (1995) Ayres, Ian, and Peter Siegelman, 1995, Race and gender discrimination in bargaining for a new car, The American Economic Review 85, 304–321.
  * Banasik et al. (1999) Banasik, John, Jonathan N. Crook, and L. C. Thomas, 1999, Not if but when will borrowers default, Journal of the Operational Research Society 50, 1185–1190.
  * Bertrand and Morse (2011) Bertrand, Marianne, and Adair Morse, 2011, Information disclosure, cognitive biases, and payday borrowing, The Journal of Finance 66, 1865–1893.
  * Beyersmann et al. (2009) Beyersmann, Jan, Aurélien Latouche, Anika Buchholz, and Martin Schumacher, 2009, Simulating competing risks data in survival analysis, Statistics in Medicine 28, 956–971.
  * Blumenstock et al. (2022) Blumenstock, Gabriel, Stefan Lessmann, and Hsin-Vonn Seow, 2022, Deep learning for survival and competing risk modelling, Journal of the Operational Research Society 73, 26–38.
  * Butler et al. (2022) Butler, Alexander W, Erik J Mayer, and James P Weston, 2022, Racial disparities in the auto loan market, The Review of Financial Studies hhac029.
  * Calem and Mester (1995) Calem, Paul S., and Loretta J. Mester, 1995, Consumer behavior and the stickiness of credit-card interest rates, The American Economic Review 85, 1327–1336.
  * Campbell (2016) Campbell, John Y., 2016, Restoring rational choice: The challenge of consumer financial regulation, American Economic Review 106, 1–30.
  * CarMax (2017) CarMax, 2017, CarMax Auto Owner Trust, Prospectus 2017-2, CarMax Business Services LLC.
  * CarMax (2019) CarMax, 2019, CarMax Auto Owner Trust, Prospectus 2019-4, CarMax Business Services LLC.
  * Cohen (1998) Cohen, Lloyd, 1998, The puzzling phenomenon of interest‐rate discounts on auto loans, The Journal of Legal Studies 27, 483–501.
  * Consumer Financial Protection Bureau (2019) Consumer Financial Protection Bureau, 2019, Borrower risk profiles, url: https://www.consumerfinance.gov/data-research/consumer-credit-trends/auto-loans/borrower-risk-profiles/ (Accessed: 2022-06-15).
  * Crowder (2001) Crowder, Martin J, 2001, Classical Competing Risks (Chapman and Hall/CRC).
  * Dai et al. (2016) Dai, Hongsheng, Marialuisa Restaino, and Huan Wang, 2016, A class of nonparametric bivariate survival function estimators for randomly censored and truncated data, Journal of Nonparametric Statistics 28, 736–751.
  * De Leonardis and Rocci (2008) De Leonardis, Daniele, and Roberto Rocci, 2008, Assessing the default risk by means of a discrete-time survival analysis approach, Appl. Stoch. Model. Bus. Ind. 24, 291–306.
  * Dirick et al. (2017) Dirick, Lore, Gerda Claeskens, and Bart Baesens, 2017, Time to default in credit scoring using survival analysis: A benchmark study, Journal of the Operational Research Society 68, 652–665.
  * Dobbie et al. (2021) Dobbie, Will, Andres Liberman, Daniel Paravisini, and Vikram Pathania, 2021, Measuring bias in consumer lending, The Review of Economic Studies 88, 2799–2832.
  * Edelberg (2006) Edelberg, Wendy, 2006, Risk-based pricing of interest rates for consumer loans, Journal of Monetary Economics 53, 2283–2298.
  * Edelberg (2007) Edelberg, Wendy, 2007, Racial dispersion in consumer credit interest rates, Finance and Economics Discussion Series 2007-28, Board of Governors of the Federal Reserve System (U.S.).
  * Einav et al. (2012) Einav, Liran, Mark Jenkins, and Jonathan Levin, 2012, Contract pricing in consumer credit markets, Econometrica 80, 1387–1432.
  * Fine and Gray (1999) Fine, Jason P, and Robert J Gray, 1999, A proportional hazards model for the subdistribution of a competing risk, Journal of the American Statistical Association 94, 496–509.
  * Frydman and Matuszyk (2022) Frydman, Halina, and Anna Matuszyk, 2022, Random survival forest for competing credit risks, Journal of the Operational Research Society 73, 15–25.
  * Fulford (2015) Fulford, Scott L., 2015, How important is variability in consumer credit limits?, Journal of Monetary Economics 72, 42–63.
  * Geskus (2011) Geskus, Ronald B., 2011, Cause-specific cumulative incidence estimation and the Fine and Gray model under both left truncation and right censoring, Biometrics 67, 39–49.
  * Gross and Souleles (2002) Gross, David B., and Nicholas S. Souleles, 2002, Do liquidity constraints and interest rates matter for consumer behavior? Evidence from credit card data, The Quarterly Journal of Economics 117, 149–185.
  * Grunewald et al. (2020) Grunewald, Andreas, Jonathan A Lanning, David C Low, and Tobias Salz, 2020, Auto dealer loan intermediation: Consumer behavior and competitive effects, Working Paper 28136, National Bureau of Economic Research.
  * Heidhues and Kőszegi (2016) Heidhues, Paul, and Botond Kőszegi, 2016, Naïeté-based discrimination, The Quarterly Journal of Economics 132, 1019–1054.
  * Huang and Wang (1995) Huang, Ying, and Mei-Cheng Wang, 1995, Estimating the occurrence rate for prevalent survival data in competing risks models, Journal of the American Statistical Association 90, 1406–1415.
  * Ishwaran et al. (2014) Ishwaran, Hemant, Thomas A. Gerds, Udaya B. Kogalur, Richard D. Moore, Stephen J. Gange, and Bryan M. Lau, 2014, Random survival forests for competing risks, Biostatistics 15, 757–773.
  * Kalbfleisch and Prentice (2011) Kalbfleisch, John D, and Ross L Prentice, 2011, The Statistical Analysis of Failure Time Data (John Wiley & Sons).
  * Karger (2003) Karger, Howard Jacob, 2003, No deals on wheels: How and why the poor pay more for basic transportation, Journal of Poverty 7, 93–112.
  * Keys et al. (2016) Keys, Benjamin J., Devin G. Pope, and Jaren C. Pope, 2016, Failure to refinance, Journal of Financial Economics 122, 482–499.
  * Klugman et al. (2012) Klugman, Stuart A., Harry H. Panjer, and Gordon E. Willmot, 2012, Loss Models: From Data to Decisions, Fourth Edition (John Wiley & Sons, Inc., Hoboken, New Jersey).
  * Lautier et al. (2021) Lautier, Jackson, Vladimir Pozdnyakov, and Jun Yan, 2021, Estimating a distribution function for discrete data subject to random truncation with an application to structured finance, arXiv 1, 1–26.
  * Lautier et al. (2023) Lautier, Jackson P., Vladimir Pozdnyakov, and Jun Yan, 2023, Pricing time-to-event contingent cash flows: A discrete-time survival analysis approach, Insurance: Mathematics and Economics 110, 53–71.
  * Lee et al. (2018a) Lee, Changhee, William Zame, Jinsung Yoon, and Mihaela van der Schaar, 2018a, Deephit: A deep learning approach to survival analysis with competing risks, Proceedings of the AAAI Conference on Artificial Intelligence 32\.
  * Lee et al. (2018b) Lee, Minjung, Eric J. Feuer, and Jason P. Fine, 2018b, On the analysis of discrete time competing risks data, Biometrics 74, 1468–1481.
  * Lehmann and Casella (1998) Lehmann, E.L., and George Casella, 1998, Theory of Point Estimation, 2nd Edition (Springer).
  * Lim et al. (2014) Lim, Younghee, Trey Bickham, Cassie M. Dinecola, Julia Broussard, Brittany E. Weber, and Alethia Gregory, 2014, Payday loan use and consumer well-being: What consumers and social workers need to know about payday loans, Journal of Poverty 18, 379–398.
  * Livshits (2015) Livshits, Igor, 2015, Recent developments in consumer credit and default literature, Journal of Economic Surveys 29, 594–613.
  * Lusardi and de Bassa Scheresberg (2013) Lusardi, Annamaria, and Carlo de Bassa Scheresberg, 2013, Financial literacy and high-cost borrowing in the United States, Working Paper 18969, National Bureau of Economic Research.
  * Manheim (2023) Manheim, 2023, Used vehicle value index, url: https://publish.manheim.com/en/services/ consulting/used-vehicle-value-index.html (Accessed: 2023-03-14).
  * Melzer (2011) Melzer, Brian T., 2011, The real costs of credit access: Evidence from the payday lending market, The Quarterly Journal of Economics 126, 517–555.
  * Mian and Sufi (2012) Mian, Atif, and Amir Sufi, 2012, The effects of fiscal stimulus: Evidence from the 2009 cash for clunkers program, The Quarterly Journal of Economics 127, 1107–1142.
  * Morse (2011) Morse, Adair, 2011, Payday lenders: Heroes or villains?, Journal of Financial Economics 102, 28–44.
  * Mukhopadhyay (2000) Mukhopadhyay, Nitis, 2000, Probability and Statistical Inference (Marcel Dekker, New York, NY).
  * Obama (2010) Obama, Barack, 2010, Remarks by the President at Signing of Dodd-Frank Wall Street Reform and Consumer Protection Act, Office of the Press Secretary, The White House.
  * Phillips (2013) Phillips, Robert, 2013, Optimizing prices for consumer credit, Journal of Revenue & Pricing Management 12\.
  * Pintilie (2006) Pintilie, Melania, 2006, Competing Risks: A Practical Perspective (John Wiley & Sons).
  * Pollard et al. (2021) Pollard, Jane, Evelyn Blumenberg, and Stephen Brumbaugh, 2021, Driven to debt: Social reproduction and (auto)mobility in los angeles, Annals of the American Association of Geographers 111, 1445–1461.
  * Prentice et al. (1978) Prentice, R. L., J. D. Kalbfleisch, A. V. Peterson, N. Flournoy, V. T. Farewell, and N. E. Breslow, 1978, The analysis of failure times in the presence of competing risks, Biometrics 34, 541–554.
  * Pressman and Scott (2009) Pressman, Steven, and Robert Scott, 2009, Consumer debt and the measurement of poverty and inequality in the US, Review of Social Economy 67, 127–148.
  * R Core Team (2022) R Core Team, 2022, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria.
  * Robb et al. (2015) Robb, Cliff, Patryk Babiarz, Ann Woodyard, and Martin Seay, 2015, Bounded rationality and use of alternative financial services, Journal of Consumer Affairs 49, 407–435.
  * Rosenbaum (2020) Rosenbaum, Eric, 2020, The used car boom is one of the hottest, and trickiest, coronavirus markets for consumers, url: https://www.cnbc.com/2020/10/15/used-car-boom-is-one-of-hottest-coronavirus-markets-for-consumers.html (Accessed: 2023-03-14).
  * Sankaran and Antony (2007) Sankaran, P.G., and Ansa Alphonsa Antony, 2007, Bivariate competing risks models under random left truncation and right censoring, Sankhyā: The Indian Journal of Statistics (2003-2007) 69, 425–447.
  * Santander (2017a) Santander, 2017a, Drive Auto Receivables Trust, Prospectus 2017-1, Santander Drive Auto Receivables LLC.
  * Santander (2017b) Santander, 2017b, Santander Drive Auto Receivables Trust, Prospectus 2017-2, Santander Drive Auto Receivables LLC.
  * Santander (2019a) Santander, 2019a, Drive Auto Receivables Trust, Prospectus 2019-4, Santander Drive Auto Receivables LLC.
  * Santander (2019b) Santander, 2019b, Santander Drive Auto Receivables Trust, Prospectus 2019-3, Santander Drive Auto Receivables LLC.
  * Schmid and Berger (2021) Schmid, Matthias, and Moritz Berger, 2021, Competing risks analysis for discrete time-to-event data, Wiley Interdisciplinary Reviews: Computational Statistics 13, e1529.
  * Securities and Exchange Commission (2014) Securities and Exchange Commission, 2014, 17 CFR Parts 229, 230, 232, 239, 240, 243, and 249 Asset-Backed Securities Disclosure and Registration.
  * Securities and Exchange Commission (2016) Securities and Exchange Commission, 2016, 17 CFR §229.1125 (Item 1125) Schedule AL - Asset-level information.
  * Stango and Zinman (2011) Stango, Victor, and Jonathan Zinman, 2011, Fuzzy math, disclosure regulation, and market outcomes: Evidence from Truth-in-Lending reform, The Review of Financial Studies 24, 506–534.
  * Staten (2015) Staten, Michael, 2015, Risk-based pricing in consumer lending, Journal of Law, Economics & Policy 11, 33–58.
  * Stepanova and Thomas (2002) Stepanova, Maria, and Lyn Thomas, 2002, Survival analysis methods for personal loan data, Operations Research 50, 277–289.
  * Thackham and Ma (2022) Thackham, Mark, and Jun Ma, 2022, On maximum likelihood estimation of competing risks using the cause-specific semi-parametric cox model with time-varying covariates – an application to credit risk, Journal of the Operational Research Society 73, 5–14.
  * Tutz and Schmid (2016) Tutz, Gerhard, and Matthias Schmid, 2016, Modeling Discrete Time-to-Event Data (Springer).
  * U.S. Government Accountability Office (2022) U.S. Government Accountability Office, 2022, Stimulus checks: Direct payments to individuals during the covid-19 pandemic, url: https://www.gao.gov/assets/gao-22-106044.pdf (Accessed: 2023-03-02).
  * Wycinka (2019) Wycinka, Ewa, 2019, Competing risk models of default in the presence of early repayments, Econometrics 23, 99–120.
  * Zhang et al. (2019) Zhang, Nailong, Qingyu Yang, Aidan Kelleher, and Wujun Si, 2019, A new mixture cure model under competing risks to score online consumer loans, Quantitative Finance 19, 1243–1253.
  * Zingales (2015) Zingales, Luigi, 2015, Presidential address: Does finance benefit society?, The Journal of Finance 70, 1327–1363.

Internet Appendix

### F. Proofs: Section II

###### Proof of Proposition 1..

Statement (i) follows from (ii), so it is enough to show (ii). Let
$\Delta+1\leq k\leq\xi$ and observe

$\displaystyle\hat{\lambda}_{\tau,n}^{0i}(k)-\lambda_{\tau}^{0i}(k)$
$\displaystyle=\frac{\frac{1}{n}\sum_{j=1}^{n}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}}{\hat{U}_{\tau,n}(k)}-\frac{f_{*,\tau}^{0i}(k)}{U_{\tau}(k)}$
$\displaystyle=\frac{\\{\sum_{j=1}^{n}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}\\}U_{\tau}(k)-f_{*,\tau}^{0i}(k)\hat{U}_{\tau,n}(k)}{\hat{U}_{\tau,n}(k)U_{\tau}(k)}$
$\displaystyle=\bigg{[}\frac{1}{\hat{U}_{\tau,n}(k)U_{\tau}(k)}\bigg{]}\frac{1}{n}\sum_{j=1}^{n}\\{\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}U_{\tau}(k)-f_{*,\tau}^{0i}(k)\mathbf{1}_{Y_{j}\leq
k\leq\min(X_{j},C_{j})}\\}.$

Define

$H^{0i}_{\tau,k(j)}=\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}U_{\tau}(k)-f_{*,\tau}^{0i}(k)\mathbf{1}_{Y_{j}\leq
k\leq\min(X_{j},C_{j})},$

for $1\leq j\leq n$ and

$\mathbf{A}_{\tau,n}=\text{diag}([\hat{U}_{\tau,n}(\Delta+1)U_{\tau}(\Delta+1)]^{-1},\ldots,[\hat{U}_{\tau,n}(\xi)U_{\tau}(\xi)]^{-1}).$

Then,

$\hat{\bm{\Lambda}}_{\tau,n}^{0i}-\bm{\Lambda}^{0i}_{\tau}=\mathbf{A}_{\tau,n}\frac{1}{n}\sum_{j=1}^{n}\begin{bmatrix}H^{0i}_{\tau,\Delta+1(j)}\\\
\vdots\\\ H^{0i}_{\tau,\xi(j)}\end{bmatrix},$

or, letting
$\mathbf{H}^{0i}_{\tau,(j)}=(H^{0i}_{\tau,\Delta+1(j)},\ldots,H^{0i}_{\tau,\xi(j)})^{\top}$
denote independent and identically distributed random vectors, we have
compactly

$\hat{\bm{\Lambda}}_{\tau,n}^{0i}-\bm{\Lambda}^{0i}_{\tau}=\mathbf{A}_{\tau,n}\frac{1}{n}\sum_{j=1}^{n}\mathbf{H}^{0i}_{\tau,(j)}.$

It is noteworthy the components of $\mathbf{H}^{0i}_{\tau,(j)}$ are
uncorrelated. More specifically,

$\text{Cov}[H^{0i}_{\tau,k(j)},H^{0i}_{\tau,k^{\prime}(j)}]=\begin{cases}U_{\tau}(k)f^{0i}_{*,\tau}(k)[U_{\tau}(k)-f^{0i}_{*,\tau}(k)],&k=k^{\prime}\\\
0,&k\neq k^{\prime}.\end{cases}$ (10)

To see this, first notice the indicator functions $\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}$ and
$\mathbf{1}_{Y_{j}\leq k\leq\min(X_{j},C_{j})}$ are Bernoulli random variables
with probability parameters $f^{0i}_{*,\tau}(k)$ and $U_{\tau}(k)$,
respectively. Hence,

$\displaystyle\mathbf{E}H^{0i}_{\tau,k(j)}$
$\displaystyle=\mathbf{E}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}U_{\tau}(k)-f^{0i}_{*,\tau}(k)\mathbf{E}\mathbf{1}_{Y_{j}\leq
k\leq\min(X_{j},C_{j})}$
$\displaystyle=f^{0i}_{*,\tau}(k)U_{\tau}(k)-f^{0i}_{*,\tau}(k)U_{\tau}(k)$
$\displaystyle=0.$

Therefore,

$\displaystyle\text{Cov}[H^{0i}_{\tau,k(j)},H^{0i}_{\tau,k^{\prime}(j)}]={}$
$\displaystyle\mathbf{E}H^{0i}_{\tau,k(j)}H^{0i}_{\tau,k^{\prime}(j)}$
$\displaystyle={}$ $\displaystyle\mathbf{E}\\{\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}U_{\tau}(k)-f_{*,\tau}^{0i}(k)\mathbf{1}_{Y_{j}\leq
k\leq\min(X_{j},C_{j})}\\}$ $\displaystyle\times\\{\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k^{\prime}}U_{\tau}(k^{\prime})-f_{*,\tau}^{0i}(k^{\prime})\mathbf{1}_{Y_{j}\leq
k^{\prime}\leq\min(X_{j},C_{j})}\\}$ $\displaystyle={}$ $\displaystyle
U_{\tau}(k)U_{\tau}(k^{\prime})\mathbf{E}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k^{\prime}}$
$\displaystyle-
U_{\tau}(k)f^{0i}_{*,\tau}(k^{\prime})\mathbf{E}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}\mathbf{1}_{Y_{j}\leq
k^{\prime}\leq\min(X_{j},C_{j})}$ $\displaystyle-
U_{\tau}(k^{\prime})f^{0i}_{*,\tau}(k)\mathbf{E}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k^{\prime}}\mathbf{1}_{Y_{j}\leq
k\leq\min(X_{j},C_{j})}$
$\displaystyle+f^{0i}_{*,\tau}(k)f^{0i}_{*,\tau}(k^{\prime})\mathbf{E}\mathbf{1}_{Y_{j}\leq
k\leq\min(X_{j},C_{j})}\mathbf{1}_{Y_{j}\leq
k^{\prime}\leq\min(X_{j},C_{j})}.$

We proceed by cases.

Case 1: $k=k^{\prime}$.

Working through each expectation in
$\text{Cov}[H^{0i}_{\tau,k(j)},H^{0i}_{\tau,k^{\prime}(j)}]$, we have

$\displaystyle\mathbf{E}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k^{\prime}}$
$\displaystyle=\mathbf{E}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}$
$\displaystyle=f^{0i}_{*,\tau}(k),$
$\displaystyle\mathbf{E}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}\mathbf{1}_{Y_{j}\leq
k^{\prime}\leq\min(X_{j},C_{j})}$
$\displaystyle=\mathbf{E}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k^{\prime}}\mathbf{1}_{Y_{j}\leq
k\leq\min(X_{j},C_{j})}$ $\displaystyle=\mathbf{E}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}\mathbf{1}_{Y_{j}\leq
k\leq\min(X_{j},C_{j})}$ $\displaystyle=\mathbf{E}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}$
$\displaystyle=f^{0i}_{*,\tau}(k),$

and

$\mathbf{E}\mathbf{1}_{Y_{j}\leq k\leq\min(X_{j},C_{j})}\mathbf{1}_{Y_{j}\leq
k^{\prime}\leq\min(X_{j},C_{j})}=\mathbf{E}\mathbf{1}_{Y_{j}\leq
k\leq\min(X_{j},C_{j})}=U_{\tau}(k).$

Thus,

$\text{Cov}[H^{0i}_{\tau,k(j)},H^{0i}_{\tau,k^{\prime}(j)}]=U_{\tau}(k)f^{0i}_{*,\tau}(k)[U_{\tau}(k)-f^{0i}_{*,\tau}(k)].$

Case 2: $k\neq k^{\prime}$.

Working through each expectation in
$\text{Cov}[H^{0i}_{\tau,k(j)},H^{0i}_{\tau,k^{\prime}(j)}]$, we have

$\mathbf{E}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k^{\prime}}=0,$
$\displaystyle\mathbf{E}$ $\displaystyle\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k}\mathbf{1}_{Y_{j}\leq
k^{\prime}\leq\min(X_{j},C_{j})}$ $\displaystyle=\begin{cases}\Pr(X_{j}\leq
C_{j},Z_{X_{j}}=i,\min(X_{j},C_{j})=k,Y_{j}\leq k^{\prime}),&k>k^{\prime}\\\
0,&k<k^{\prime},\end{cases}$ $\displaystyle\mathbf{E}$
$\displaystyle\mathbf{1}_{X_{j}\leq
C_{j}}\mathbf{1}_{Z_{X_{j}}=i}\mathbf{1}_{\min(X_{j},C_{j})=k^{\prime}}\mathbf{1}_{Y_{j}\leq
k\leq\min(X_{j},C_{j})}$ $\displaystyle=\begin{cases}0,&k>k^{\prime}\\\
\Pr(X_{j}\leq C_{j},Z_{X_{j}}=i,\min(X_{j},C_{j})=k^{\prime},Y_{j}\leq
k),&k<k^{\prime},\end{cases}$

and

$\mathbf{E}\mathbf{1}_{Y_{j}\leq k\leq\min(X_{j},C_{j})}\mathbf{1}_{Y_{j}\leq
k^{\prime}\leq\min(X_{j},C_{j})}=\Pr(Y_{j}\leq
k\leq\min(X_{j},C_{j}),Y_{j}\leq k^{\prime}\leq\min(X_{j},C_{j})).$

Thus,

Cov
$\displaystyle[H^{0i}_{\tau,k(j)},H^{0i}_{\tau,k^{\prime}(j)}]=f^{0i}_{*,\tau}(\min(k,k^{\prime}))\bigg{\\{}$
$\displaystyle-U_{\tau}(\max(k,k^{\prime}))\Pr(X_{j}\leq
C_{j},Z_{X_{j}}=i,\min(X_{j},C_{j})=\max(k,k^{\prime}),Y_{j}\leq\min(k,k^{\prime}))$
$\displaystyle+f^{0i}_{*,\tau}(\max(k,k^{\prime}))\Pr(Y_{j}\leq
k\leq\min(X_{j},C_{j}),Y_{j}\leq k^{\prime}\leq\min(X_{j},C_{j}))\bigg{\\}}.$

However, because of the independence between $Y$ and $(X,Z_{X})$,

$\displaystyle U_{\tau}(\max(k,k^{\prime}))$
$\displaystyle=\Pr(Y_{j}\leq\max(k,k^{\prime})\leq\min(X_{j},C_{j}))$
$\displaystyle=\Pr(Y\leq\max(k,k^{\prime}),X\geq\max(k,k^{\prime}),C\geq\max(k,k^{\prime})\mid
Y\leq X)$ $\displaystyle=\\{\Pr(Y\leq\max(k,k^{\prime})\leq
C)\Pr(X\geq\max(k,k^{\prime}))\\}/\alpha,$ $\displaystyle\Pr$
$\displaystyle(X_{j}\leq
C_{j},Z_{X_{j}}=i,\min(X_{j},C_{j})=\max(k,k^{\prime}),Y_{j}\leq\min(k,k^{\prime}))$
$\displaystyle=\Pr(C\geq\max(k,k^{\prime}),Z_{X}=i,X=\max(k,k^{\prime}),Y\leq\min(k,k^{\prime})\mid
Y\leq X)$
$\displaystyle=\\{\Pr(X=\max(k,k^{\prime}),Z_{X}=i)\Pr(Y\leq\min(k,k^{\prime}),C\geq\max(k,k^{\prime}))\\}/\alpha,$
$\displaystyle f^{0i}_{*,\tau}(\max(k,k^{\prime}))$
$\displaystyle=\Pr(X=\max(k,k^{\prime}),C\geq\max(k,k^{\prime}),Z_{x}=i\mid
Y\leq X)$
$\displaystyle=\\{\Pr(X=\max(k,k^{\prime}),Z_{X}=i)\Pr(Y\leq\max(k,k^{\prime})\leq
C)\\}/\alpha,$

and

$\displaystyle\Pr$ $\displaystyle(Y_{j}\leq k\leq\min(X_{j},C_{j}),Y_{j}\leq
k^{\prime}\leq\min(X_{j},C_{j}))$
$\displaystyle=\Pr(Y\leq\min(k,k^{\prime}),C\geq\max(k,k^{\prime}),X\geq\max(k,k^{\prime})\mid
Y\leq X)$
$\displaystyle=\\{\Pr(Y\leq\min(k,k^{\prime}),C\geq\max(k,k^{\prime}))\Pr(X\geq\max(k,k^{\prime}))\\}/\alpha$

Therefore,

$\displaystyle U_{\tau}$ $\displaystyle(\max(k,k^{\prime}))\Pr(X_{j}\leq
C_{j},Z_{X_{j}}=i,\min(X_{j},C_{j})=\max(k,k^{\prime}),Y_{j}\leq\min(k,k^{\prime}))$
$\displaystyle=f^{0i}_{*,\tau}(\max(k,k^{\prime}))\Pr(Y_{j}\leq
k\leq\min(X_{j},C_{j}),Y_{j}\leq k^{\prime}\leq\min(X_{j},C_{j})),$

and so $\text{Cov}[H^{0i}_{\tau,k(j)},H^{0i}_{\tau,k^{\prime}(j)}]=0$ when
$k\neq k^{\prime}$. This confirms (10). Now define

$\mathbf{D}^{0i}_{\tau}=\text{diag}\big{(}U_{\tau}(\Delta+1)f^{0i}_{*,\tau}(\Delta+1)[U_{\tau}(\Delta+1)-f^{0i}_{*,\tau}(\Delta+1)],\ldots,U_{\tau}(\xi)f^{0i}_{*,\tau}(\xi)[U_{\tau}(\xi)-f^{0i}_{*,\tau}(\xi)]\big{)},$

and

$\bar{\mathbf{H}}_{\tau,n}^{0i}=\frac{1}{n}\sum_{j=1}^{n}\mathbf{H}^{0i}_{\tau,(j)}.$

By the multivariate Central Limit Theorem (Lehmann and Casella, 1998, Theorem
8.21, pg. 61), therefore,

$\sqrt{n}(\bar{\mathbf{H}}_{\tau,n}^{0i}-\bm{0})\overset{\mathcal{L}}{\longrightarrow}N(\bm{0},\mathbf{D}^{0i}_{\tau}),\text{
as }n\rightarrow\infty.$

Next, define
$\mathbf{V}_{\tau}=\text{diag}(U_{\tau}(\Delta+1)^{-2},\ldots,U_{\tau}(\xi)^{-2})$.
By Lemma 1 (Lautier et al., 2023),
$\mathbf{A}_{\tau,n}\overset{\mathcal{P}}{\longrightarrow}\mathbf{V}_{\tau}$,
as $n\rightarrow\infty$. Thus, by multivariate Slutsky’s Theorem (Lehmann and
Casella, 1998, Theorem 5.1.6, pg. 283),

$\sqrt{n}(\mathbf{A}_{\tau,n}\bar{\mathbf{H}}_{\tau,n}^{0i})\overset{\mathcal{L}}{\longrightarrow}N(\bm{0},\mathbf{V}_{\tau}\mathbf{D}^{0i}_{\tau}\mathbf{V}_{\tau}^{\top}),\text{
as }n\rightarrow\infty.$

We may complete the proof by observing
$\mathbf{V}_{\tau}\mathbf{D}^{0i}_{\tau}\mathbf{V}_{\tau}^{\top}=\bm{\Sigma}^{0i}$
and
$\mathbf{A}_{\tau,n}\bar{\mathbf{H}}_{\tau,n}^{0i}=\hat{\bm{\Lambda}}^{0i}_{\tau,n}-\bm{\Lambda}_{\tau}^{0i}$.
∎

###### Proof of Lemma 1.

The classical method dictates first finding a $(1-\theta)$% confidence
interval on a log-scale and then converting back to a standard-scale to ensure
the estimated confidence interval for the hazard rate, which is a probability,
remains in the interval $(0,1)$. By an application of the Delta Method
(Lehmann and Casella, 1998, Theorem 8.12, pg. 58), we have for
$x\in\\{\Delta+1,\ldots,\xi\\}$ and $i=1,2$,

$\sqrt{n}\big{(}\ln\hat{\lambda}_{\tau,n}^{0i}(x)-\ln\lambda^{0i}_{\tau}(x)\big{)}\overset{\mathcal{L}}{\longrightarrow}N\bigg{(}0,\frac{f_{*,\tau}^{0i}(x)\\{U_{\tau}(x)-f_{*,\tau}^{0i}(x)\\}}{U_{\tau}(x)^{3}}\frac{1}{\lambda^{0i}_{\tau}(x)^{2}}\bigg{)}.$

The result follows from (4), the Continuous Mapping Theorem (Mukhopadhyay,
2000, Theorem 5.2.5, pg. 249), the pivotal approach (Mukhopadhyay, 2000,
§9.2.2), and converting back to the standard scale. ∎

### G. Proofs: Section E

###### Proof of Proposition 2.

For a loan with initial balance, $B$, monthly interest rate, $r_{a}$, and
initial term of $\xi$, the monthly payment, $P$, is

$P=B\bigg{[}\frac{1-(1+r_{a})^{-\xi}}{r_{a}}\bigg{]}^{-1}.$

Assume $x\in\\{1,\ldots,\xi\\}$. The balance at month $x$, $B_{x}$ is

$\displaystyle B_{x}$
$\displaystyle=B(1+r_{a})^{x}-P\bigg{[}\frac{(1+r_{a})^{x}-1}{r_{a}}\bigg{]}$
$\displaystyle=B(1+r_{a})^{x}-B\bigg{[}\frac{1-(1+r_{a})^{-\xi}}{r_{a}}\bigg{]}^{-1}\bigg{[}\frac{(1+r_{a})^{x}-1}{r_{a}}\bigg{]}.$
(11)

Thus, $\rho_{a\mid x}$ is the rate such that the expected present value of the
future monthly payments equals $B_{x}$. The payment stream is constant,
however, and so

$\displaystyle B_{x}$ $\displaystyle=P\bigg{[}\frac{1}{(1+\rho_{a\mid
x})}+\cdots+\frac{1}{(1+\rho_{a\mid x})^{\xi-x}}\bigg{]}$
$\displaystyle=B\bigg{[}\frac{1-(1+r_{a})^{-\xi}}{r_{a}}\bigg{]}^{-1}\bigg{[}\frac{1-(1+\rho_{a\mid
x})^{-(\xi-x)}}{\rho_{a\mid x}}\bigg{]}.$

Use (11) and solve for $\rho_{a\mid x}$ to complete the proof. ∎

###### Proof of Lemma 2.

The result follows by Proposition 1, part $(i)$ and the Continuous Mapping
Theorem (Mukhopadhyay, 2000, Theorem 5.2.5, pg. 249). ∎

### H. Large Sample Simulation Study

We present a simulation study in support of Proposition 1 and Lemma 1. Let the
true distribution for the lifetime random variable $X$ and bivariate
distribution of $(X,Z_{X})$ be as in Table 4. The column $p(x)$ denotes the
probability of event type 1 given an event at time $X$. This allows us to
populate the joint distribution for $\Pr(X=x,Z_{X}=i)$ for $i=1,2$. The cause-
specific hazard rates then follow from (3), and we also report the all-cause
hazard rate in the final column. Notice that, for each $x$,

$p(x)=\frac{\lambda^{01}(x)}{\lambda^{01}(x)+\lambda^{02}(x)}.$

Table 4:  Simulation Study Lifetime of Interest Probabilities. The true probabilities of the lifetime random variable, $X$, for the simulation study results of Figure 14. The probabilities $p(x)$ and $\Pr(X=x)$ for $x\in\\{1,\ldots,10\\}$ are selected at onset, and the remaining probabilities in this table may be derived from these quantities. Not summarized here is the truncation random variable, $Y$, which was assumed to be discrete uniform over the integers $\\{1,\ldots,5\\}$.  $p(x)$ | $X$ | $\Pr(X=x)$ | $\Pr(X=x,Z_{x}=1)$ | $\Pr(X=x,Z_{x}=2)$ | $\lambda^{01}(x)$ | $\lambda^{02}(x)$ | $\lambda(x)$
---|---|---|---|---|---|---|---
0.66 | 1 | 0.04 | 0.026 | 0.014 | 0.026 | 0.014 | 0.04
0.20 | 2 | 0.06 | 0.012 | 0.048 | 0.013 | 0.050 | 0.06
0.45 | 3 | 0.10 | 0.045 | 0.055 | 0.050 | 0.061 | 0.11
0.87 | 4 | 0.14 | 0.122 | 0.018 | 0.152 | 0.023 | 0.18
0.20 | 5 | 0.09 | 0.018 | 0.072 | 0.027 | 0.109 | 0.14
0.81 | 6 | 0.06 | 0.049 | 0.011 | 0.085 | 0.020 | 0.11
0.05 | 7 | 0.14 | 0.007 | 0.133 | 0.014 | 0.261 | 0.27
0.78 | 8 | 0.18 | 0.140 | 0.040 | 0.379 | 0.107 | 0.49
0.25 | 9 | 0.07 | 0.018 | 0.053 | 0.092 | 0.276 | 0.37
0.42 | 10 | 0.12 | 0.050 | 0.070 | 0.420 | 0.580 | 1.00

For the truncation random variable, we assume $Y$ is discrete uniform with
sample space $\mathcal{Y}\in\\{1,2,3,4,5\\}$. This results in $\alpha=0.864$.
For the purposes of the simulation, we further assume $\tau=5$. We use the
simulation procedure of Beyersmann et al. (2009) but modified for random
truncation. Specifically,

  1. 1.

Simulate the truncation time, $Y$.

  2. 2.

Set the censoring time to be $Y+\tau$.

  3. 3.

Simulate the event time, $X$.

  4. 4.

Simulate a Bernoulli event with probability $p(x)$ to determine if the event
$X$ was caused by type 1 with probability $p(x)$ or type 2 with probability
$1-p(x)$.

We simulated $n=10{,}000$ lifetimes using the above algorithm. We then tossed
any observations that were truncated (i.e., $Y_{j}>X_{j}$, for
$j=1,\ldots,n$). This left a sample of competing risk events subject to
censoring, which would be the same incomplete data conditions as a trust of
securitized loans. We then used the results of Section II.A to estimate
$\hat{f}^{0i}_{*,\tau,n}(x)$, $\hat{U}_{\tau,n}(x)$, and
$\hat{\lambda}^{0i}_{\tau,n}(x)$ for $i=1,2$ and $x\in\\{1,\ldots,10\\}$ over
$r=1{,}000$ replicates.

To validate the asymptotic results of Proposition 1, we compare the empirical
covariance matrix against the derived asymptotic covariance matrix,
$\bm{\Sigma}^{0i}$, by examining estimates of the confidence intervals using
Lemma 1. Figure 14 presents the results for the cause-specific hazard rate for
cause 01 and 02, respectively. The empirical estimates and 95% confidence
intervals are indistinguishable from the true quantities using Proposition 1
and estimated quantities using Proposition 1 but replacing all quantities with
their respective estimates from Section II.A.

Figure 14:  Simulation Study Results. A comparison of true
$\lambda^{0i}_{\tau}(x)$ and estimated $\hat{\lambda}^{0i}_{\tau,n}(x)$,
including confidence intervals, for the distribution in Table 4 and $i=1,2$.
The “true” values are from Proposition 1 and Lemma 1. The “estimate” values
use the formulas from Proposition 1 and Lemma 1 but replace the true values
with the estimates from Section II.A calculated from the simulated data. The
“empirical” values are empirical confidence interval and mean calculations
directly from the simulated data. All three quantities are indistinguishable
for $n=10{,}000$ and 1,000 replicates, which indicates the asymptotic
properties hold in this instance.

### I. Determination of Loan Outcome

The detail of the loan-level data is extensive, but it remains up to the data
analyst to use the provided fields to determine the outcome of an individual
loan (see Securities and Exchange Commission (2016) for detail on available
field names). To do so, we have aggregated each month of active trust data
into a single source file. This allows us to review each bond’s monthly
outstanding principal balance, monthly payment received from the borrower, and
the portion of each monthly payment applied to principal. Our algorithm to
determine a loan outcome proceeded as follows. For each remaining bond after
the filtering of Section I.A, we extracted three vectors, each of which was
the same length as the number of months a trust was active and paying. The
first vector represented the ordered monthly balance, the second was the
ordered monthly payments, and the third was the ordered monthly amount of
payment applied to principal. We then considered a loan to be repaid if the
sum total principal received was greater than the outstanding loan balance as
of the first month the trust was actively paying. In this case, the timing of
a repayment was set to be the first month with a zero outstanding principal
balance. Note that we do not differentiate between a prepayment or naturally
scheduled loan amortization; i.e., all repayments have been treated as a “non-
default”. If the sum total principal received was less than the first month’s
outstanding loan balance, we then considered a loan outcome to be either
right-censored or defaulted. To make this determination, we searched the
monthly payments received vector for three consecutive zeros (i.e., three
straight months of missed payments). If we found three consecutive missed
payments, we assumed the loan to be defaulted with a time-of-default set to be
the month in which the first of three zeros was observed. If we did not find
three consecutive months of missed payments, the loan was assumed to be a
right-censored observation and assigned an event time as of the last month the
trust was actively paying. For the pseudo-code of this algorithm, see Figure
15.

1:$B\leftarrow\texttt{bond\\_data}$ $\triangleright$ bond_data is a row of the
loan performance data

2:$\texttt{bal\\_vec}\leftarrow\text{each month's sequential outstanding
principal balance}$

3:$\texttt{pmt\\_vec}\leftarrow\text{each month's sequential actual payment}$

4:$\texttt{prc\\_vec}\leftarrow\text{each month's sequential payment applied
to principal}$

5:$\texttt{init\\_bal}\leftarrow\text{current balance as of the first trust
month}$

6:$\texttt{paid\\_princ}\leftarrow\texttt{sum}(\texttt{prc\\_vec})$
$\triangleright$ plus $10 pad to avoid odd tie behavior

7:if $\texttt{paid\\_princ}>=\texttt{init\\_bal}$ then

8: $D=0$

9: $R=1$

10: $C=0$

11: $X\leftarrow$ location of first zero in bal_vec $\triangleright$ loan
repaid

12:else

13: $z\leftarrow$ starting time of three consecutive zero payments in pmt_vec

14: if $z$ empty then

15: $D=0$

16: $R=0$

17: $C=1$

18: $X\leftarrow$ length of pmt_vec $\triangleright$ loan censored

19: else

20: $D=1$

21: $R=0$

22: $C=0$

23: $X\leftarrow z$ $\triangleright$ loan defaults

24: end if

25:end if

Figure 15:  Determination of Loan Outcome. We first extracted three vectors,
each of which was the same length as the number of months a trust was active
and paying. The first vector (bal_vec) represented the ordered monthly
balance, the second (pmt_vec) was the ordered monthly payments, and the third
(prc_vec) was the ordered monthly amount of payment applied to principal. We
then considered a loan to be repaid if the sum total principal received was
greater than the outstanding loan balance as of the first month the trust was
actively paying. In this case, the timing of a repayment was set to be the
first month with a zero outstanding principal balance. If the sum total
principal received was less than the first month’s outstanding loan balance,
we then considered a loan outcome to be either right-censored or defaulted. To
make this determination, we searched the monthly payments received vector for
three consecutive zeros (i.e., three straight months of missed payments). If
we found three consecutive missed payments, we assumed the loan to be
defaulted with a time-of-default set to be the month in which the first of
three zeros was observed. If we did not find three consecutive months of
missed payments, the loan was assumed to be a right-censored observation and
assigned an event time as of the last month the trust was actively paying.


* Dimastrogiovanni  _et al._ (2010) E. Dimastrogiovanni, N. Bartolo, S. Matarrese, and A. Riotto, Adv. Astron. 2010, 752670 (2010), arXiv:1001.4049 [astro-ph.CO] .
  * Bamba  _et al._ (2008) K. Bamba, S. Nojiri, and S. D. Odintsov, Phys. Rev. D 77, 123532 (2008), arXiv:0803.3384 [hep-th] .
  * Garnica  _et al._ (2022) J. C. Garnica, L. G. Gomez, A. A. Navarro, and Y. Rodriguez, Annalen Phys. 534, 2100453 (2022), arXiv:2109.10154 [gr-qc] .
  * Armendariz-Picon (2004) C. Armendariz-Picon, JCAP 0407, 007 (2004), arXiv:astro-ph/0405267 .
  * Koivisto and Mota (2008a) T. Koivisto and D. F. Mota, Astrophys. J. 679, 1 (2008a), arXiv:0707.0279 [astro-ph] .
  * Koivisto and Mota (2008b) T. Koivisto and D. F. Mota, JCAP 06, 018 (2008b), arXiv:0801.3676 [astro-ph] .
  * Mehrabi  _et al._ (2017) A. Mehrabi, A. Maleknejad, and V. Kamali, Astrophys. Space Sci. 362, 53 (2017), arXiv:1510.00838 [astro-ph.CO] .
  * Guarnizo  _et al._ (2020) A. Guarnizo, J. B. Orjuela-Quintana, and C. A. Valenzuela-Toledo, Phys. Rev. D 102, 083507 (2020), arXiv:2007.12964 [gr-qc] .
  * Tasinato (2014a) G. Tasinato, JHEP 1404, 067 (2014a), arXiv:1402.6450 [hep-th] .
  * De Felice  _et al._ (2016) A. De Felice _et al._ , JCAP 06, 048 (2016), arXiv:1603.05806 [gr-qc] .
  * Rodríguez and Navarro (2018) Y. Rodríguez and A. A. Navarro, Phys. Dark Univ. 19, 129 (2018), arXiv:1711.01935 [gr-qc] .
  * Nakamura  _et al._ (2017) S. Nakamura, R. Kase, and S. Tsujikawa, Phys. Rev. D 95, 104001 (2017), arXiv:1702.08610 [gr-qc] .
  * Geng  _et al._ (2021) C.-Q. Geng, Y.-T. Hsu, J.-R. Lu, and L. Yin, Phys. Dark Univ. 32, 100819 (2021), arXiv:2104.06577 [gr-qc] .
  * Heisenberg and Villarrubia-Rojo (2021) L. Heisenberg and H. Villarrubia-Rojo, JCAP 03, 032 (2021), arXiv:2010.00513 [astro-ph.CO] .
  * Zhao (2009) W. Zhao, Int. J. Mod. Phys. D18, 1331 (2009), arXiv:0810.5506 [gr-qc] .
  * Koivisto and Nunes (2013) T. S. Koivisto and N. J. Nunes, Phys. Rev. D88, 123512 (2013), arXiv:1212.2541 [astro-ph.CO] .
  * Ngampitipan and Wongjun (2011) T. Ngampitipan and P. Wongjun, JCAP 1111, 036 (2011), arXiv:1108.0140 [hep-ph] .
  * Yao  _et al._ (2018) Y.-H. Yao, Y.-J. Yan, and X.-H. Meng, Eur. Phys. J. C78, 153 (2018), arXiv:1704.05772 [astro-ph.CO] .
  * Wei and Cai (2006) H. Wei and R.-G. Cai, Phys. Rev. D73, 083002 (2006), arXiv:astro-ph/0603052 .
  * Landim (2016) R. C. G. Landim, Eur. Phys. J. C76, 480 (2016), arXiv:1605.03550 [hep-th] .
  * Koivisto and Mota (2008c) T. Koivisto and D. F. Mota, JCAP 0808, 021 (2008c), arXiv:0805.4229 [astro-ph] .
  * Nakamura  _et al._ (2019) S. Nakamura, R. Kase, and S. Tsujikawa, JCAP 1912, 032 (2019), arXiv:1907.12216 [gr-qc] .
  * Yao and Meng (2021) Y. Yao and X.-H. Meng, Phys. Dark Univ. 33, 100852 (2021), arXiv:2011.09160 [astro-ph.CO] .
  * Allys  _et al._ (2016b) E. Allys, P. Peter, and Y. Rodriguez, JCAP 1602, 004 (2016b), arXiv:1511.03101 [hep-th] .
  * Kimura  _et al._ (2017) R. Kimura, A. Naruko, and D. Yoshida, JCAP 1701, 002 (2017), arXiv:1608.07066 [gr-qc] .
  * Gallego Cadavid  _et al._ (2022) A. Gallego Cadavid, C. M. Nieto, and Y. Rodriguez, Phys. Rev. D 105, 124060 (2022), arXiv:2110.14623 [hep-th] .
  * Gumrukcuoglu and Namba (2019) A. E. Gumrukcuoglu and R. Namba, Phys. Rev. D100, 124064 (2019), arXiv:1907.12292 [hep-th] .
  * Allys  _et al._ (2016c) E. Allys, P. Peter, and Y. Rodriguez, Phys. Rev. D94, 084041 (2016c), arXiv:1609.05870 [hep-th] .
  * Gómez and Rodríguez (2019) L. G. Gómez and Y. Rodríguez, Phys. Rev. D 100, 084048 (2019), arXiv:1907.07961 [gr-qc] .
  * Carroll (1998) S. M. Carroll, Phys. Rev. Lett. 81, 3067 (1998), arXiv:astro-ph/9806099 .
  * Tasinato (2014b) G. Tasinato, Class. Quant. Grav. 31, 225004 (2014b), arXiv:1404.4883 [hep-th] .
  * Silva  _et al._ (2022) H. O. Silva, A. Coates, F. M. Ramazanoğlu, and T. P. Sotiriou, Phys. Rev. D 105, 024046 (2022), arXiv:2110.04594 [gr-qc] .
  * Coates and Ramazanoğlu (2022) A. Coates and F. M. Ramazanoğlu, Phys. Rev. Lett. 129, 151103 (2022), arXiv:2205.07784 [gr-qc] .
  * Coates and Ramazanoğlu (2023) A. Coates and F. M. Ramazanoğlu, (2023), arXiv:2301.02263 [gr-qc] .
  * Clough  _et al._ (2022) K. Clough, T. Helfer, H. Witek, and E. Berti, Phys. Rev. Lett. 129, 151102 (2022), arXiv:2204.10868 [gr-qc] .
  * Mou and Zhang (2022) Z.-G. Mou and H.-Y. Zhang, Phys. Rev. Lett. 129, 151101 (2022), arXiv:2204.11324 [hep-th] .
  * Aoki and Minamitsuji (2022) K. Aoki and M. Minamitsuji, Phys. Rev. D 106, 084022 (2022), arXiv:2206.14320 [gr-qc] .
  * Barausse  _et al._ (2022) E. Barausse, M. Bezares, M. Crisostomi, and G. Lara, JCAP 11, 050 (2022), arXiv:2207.00443 [gr-qc] .
  * Sbisà (2015) F. Sbisà, Eur. J. Phys. 36, 015009 (2015), arXiv:1406.4550 [hep-th] .
  * Gümrükçüoğlu  _et al._ (2016) A. E. Gümrükçüoğlu, S. Mukohyama, and T. P. Sotiriou, Phys. Rev. D 94, 064001 (2016), arXiv:1606.00618 [hep-th] .
  * Langlois and Noui (2016) D. Langlois and K. Noui, JCAP 07, 016 (2016), arXiv:1512.06820 [gr-qc] .
  * Di Valentino  _et al._ (2021) E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, Class. Quant. Grav. 38, 153001 (2021), arXiv:2103.01183 [astro-ph.CO] .