A visit generation process for human mobility random graphs with location-specific latent-variables: from land use to travel demand

This research introduces a mathematical framework to comprehending human mobility patterns, integrating mathematical modeling and economic analysis. The study focuses on latent-variable networks, investigating the dynamics of human mobility using stochastic models. By examining actual origin-destination data, the research reveals scaling relations and uncovers the economic implications of mobility patterns, such as the income elasticity of travel demand. The mathematical analysis commences with the development of a stochastic model based on inhomogeneous random graphs to construct a visitation model with multipurpose drivers for travel demand. A directed multigraph with weighted edges is considered, incorporating trip costs and labels to represent factors like distance traveled and travel time. The study gains insights into the structural properties and dynamic correlations of human mobility networks, to derive analytical and computational solutions for key network metrics, including scale-free behavior of the strength and degree distribution, together with the estimation of assortativity and clustering coefficient. Additionally, the model's validity is assessed through a real-world case study of the New York metropolitan area. The analysis of this data exposes clear scaling relations in commuting patterns, confirming theoretical predictions and validating the efficacy of the mathematical model. The model further explains a series of scaling behaviors in origin-destination flows among areas of a region, successfully reproducing statistical regularities observed in real-world cases using extensive human mobility datasets. In particular, the model's application to estimating income elasticity of travel demand bears significant implications for urban and transport economics.


Introduction
How people move from one place to another represents a crucial key to depicting human relations and social interactions in complex societies.Data-driven structured studies on human mobility are valuable to identify the factors that drive the movement of individuals and goods, pointing out how they affects economic outcomes in field as labor market dynamics, economic growth, transportation planning and consumer behavior [1,2].Human mobility encompasses a wide range of spatial and temporal scales, from daily commuting within a city to long-term migration across countries and its dynamics are influenced by factors like economic development, technological advancements, political stability and natural features of the territory.Human mobility and transport dynamics has been studied through different mathematical approaches ranging from individual path-based models to collective population based ones starting from single individual's behavior up to examining collective mobility trends [3][4][5].The study presented in this work aims to provide a stochastic model for complex network of human mobility based on an origin-destination structure which represents the backbone of mobility visitation flows with individuals moving from origin locations to destinations.Then, inferring latent-variables in a specific case study, the graph statistical patterns observed in origin-destination network flow can be recovered to have the same expected properties computed trough the latent-variable estimations.Specifically, the model is grounded in the class of inhomogeneous random graphs [6,7] also known as latent-variable networks [8][9][10][11][12].Moreover, the mobility network is a directed multigraph with weighted edges, since each travel is labeled with weights (or trip jumps) that represent trip costs such as distance traveled, travel time or emissions per trip.The use of a graph representation of human mobility helps to investigate the dynamic behavior of the system starting from the properties of single constituents, and to estimate how much one part of the network influences another by using graph metrics like degree distribution, centrality, assortativity and transitivity of trip mobility network.Following the approach of dynamical formation of latent-variable graphs [13][14][15], the network has been, then, formalized by using a master equation for the probability density function which encodes the distribution of the number of visits in destination locations weighted by the costs of trips to get there.In particular, the evolution of visit distribution over time is described trough an integro-differential equation with an explicit asymptotic solution which highlights the relation of visiting generation patterns to intrinsic and environmental features of locations as specified by latent-variable properties.In stochastic theory, the analytical interpretation achieved is related to the literature of continuous-time Markov processes with the use of a Kolmogorov-Feller partial differential equation [16,17].The real-valued solution of the visit ditribution of the trip mobility network has been as compound distribution of traveling packet of normal like conditional probabilities for independent visitation process for each type of destinations.Such visit generation process of the temporal random graph can be is then proofed to be equivalently formalized in terms of a mixture of compound Poisson counting processes commonly used in financial and actuarial science literature [18,19], for example in models of option pricing, risk and insurance analysis, high-frequency trading and market microstructure.In the case of interest, the process is driven by hidden graph structures represented by latent variables associated to each location, that will be identified with an attractiveness attribute and a productiveness one.The first variable encodes an intrinsic capability of a location to attract visits due to some specific properties of the destination, as extensively discussed in the paper.The second variable encodes, complementary, the ability to produce new travelers, and it depends on some intrinsic feature of the origin area.Moreover the way how each area produces or attracts new trips is achieved according to arrival and departure rates that will be express by some function of the latent variables.In summary, the mobility graph process is fully characterized by the latent-variable statistical characteristics, and different configurations of the latent-variable patterns will determine the mobility network scaling properties with a particular attention on scale-free behavior of the degree distribution and the spectral trend of graph correlations.In a economic modeling framework, those latent variables allows the characterization of travel demand for different areas in the region, and they can be inferred through the analysis of important demographic, geographic and economic indicators.As discussed in literature [20][21][22], many different factors are involved for a place to be considered as attractive, such as trip purposes, job or leisure opportunities, infrastructure facilities, geographical characteristics, urban zoning planing.Generally, all those physical and human factors can be captured by land use and travel behavior analysis which are at the basis of two primary approaches in transportation economics and engineering literature i.e. trip-based models and activity-based models.In such literature, in fact, travel is considered to be a derived demand, that it is generated in response to people satisfying personal needs and desires [23,24].Moreover, the costs and efforts to realize the necessary trips are then influenced and determined by geographical, social and economic factors such as the distance from means of transport, house affordability, economic status, social equality, the proximity to nature and many others.The research direction is then focused on the use of a real mobility network in order to validate the model presented.In this work, New York metropolitan area has been used as case study by using detailed origin-destination tables where many important information such as the number of people that travel between different locations, along with the information on costs of the visits as well as demographic and economic properties of the travelers.The analysis of such data reveals clear scaling relations in commuting patterns of movement flow in agreement with empirical studies and consistent with theoretical predictions as in [25][26][27][28].
The first contribution of the paper is to provide a theoretical network framework able to describe human mobility flow based on the existence of latent variables which characterize the travel demand features in urban and transport dynamics.Specifically, analytical, numerical and computational solutions are provided for the strength distribution, assortativity and clustering coefficient.The second contribution consists on the analysis of real origin destination network, so revealing the existence of scale-free behaviors in the frequency of visits of locations and scaling relations between visits and trip costs.Consequently, a more pragmatic original contribution of the study is to determine what are the latent variable statistical features able to reproduce the observed patterns of the trip mobility network in terms of distribution and correlations.The study emphasizes a reciprocal interplay between human mobility and urban structure, whi is an important topic widely discussed in the literature [21,22,[29][30][31][32].Finally, as economical application, I will show the relation between the attractiveness scaling exponents and the value of income elasticity passing through the information about allocation and utilization of land resources for various economic and social activities.The estimated income elasticity of demand is often used to predict future changes in consumption in response to changes in income.[33][34][35][36].
The paper is organized in the following way: the section 2 is devoted to the definition and mathematical modeling of the latent variable mobility network and, then, analytical estimates of scaling laws and graph correlations of visitation patterns.The section 3 will be devoted to network measurements in the case study of New York metropolitan area with the analysis of origin destination tables obtained from Safegraph dataset [37].Visit distribution and network correlations will be measured.In Section 4, the network model is implemented and validated on the basis of latent variable specification.Finally, the latent variable formalism will be applied to an economic problem of the estimation of income elasticity of travel demand.An appendix section and supplementary material are provided to enhance the discussion on the stochastic intepretation of the visit process, further statistical analysis of the data, and a detailed description of data and interpretation of the latent variables in travel demand modeling.Figure 1: Example of a tessellation of a large urban area (i.e.New York City) in (Census) blocks (a) as in [38].In particular, a block is characterized by an intrinsic attractiveness x as a latent variable depending on many features of location (b).

Model
A trip mobility network can be build upon an origin-destination rationale and it will be represented as a directed graph where the nodes are administrative units of an urban or regional area (a city, a county, a state etc...) and the directed edges represent number of visits from an origin to a destination area.Let us observe, that one latent-variable is the attractiveness of an area and it can be thought as an ancestral property of the block which captures ideally all the variables that can determine it, such as the number residential apartments, job positions, retail stores, geographical features, geopolitical conditions, social and economical factors, transportation systems and facilities including also the average distance of a block with attractiveness from the rest of the urban area (i.e. the physical average accessibility).Essentially those latent-variables, from an economic perspective, can be seen as drivers of travel demand and trip behavior in general.Each node is characterized by an intrinsic attractiveness x representing the travel demand for a destination area.Such attribute is related to various properties of the location such as job opportunities, number of retail stores, geographic features, infrastructure, facilities, school districts etc.Similarly, each node is also characterized by another property y that characterizes the location as an origin of trips.That attribute captures the possible users which can depart from the area, so capturing the potential of the are to produce trips, see the Supplementary Materials for a detailed discussion.The model presented describes the occurrence dynamics of trips that occur between origins with a given population y and destination blocks with a given attractiveness label x.In this paper, I will, primarily, focus on the case of exchangeable origin-destination blocks, so that the trip generating process can be represented by two independent process: the trip production and the trip attraction process.The first accounts for the number of trips (departures) originating from a block and the latter accounts for the number of trips (visits) ending in each destination block.
First, the visitation model is defined as a time-varying network made up of three ingredients: 1) the location structure of the graph in terms of latent variables, 2) the visit arrival process to destinations, and 3) the effect of trips as sizes assigned to each visit.Second, properties of dynamic mobility network will be presented in terms of the its strength degree distribution and degree-correlations, where the strength of a location is intended as number of visits weighted by their correspondent trip size.

Definition
A geographical region of interest, R, is specified as the portion of territory for which we are interested in generating the flows.Over the region of interest, a set of geographical tiles called tessellation, T, made up of locations l i so that T = {l i : i = 1 . . .n} so that the locations are non-overlapping, l i ∩ l j = ∅, ∀i ̸ = j, and the union of all locations completely covers the region of interest, ∪ n i=1 l i = R.The tessellation for real geographical regions can be obtained in many ways according to the scope.In the case of interest, location tiles are the census areas defined by national authorities for administrative and demographic purposes.In particular e census blocks are the smallest geographic unit used by the United States Census Bureau for tabulation of data collected from all houses.An example of tessellation is given in Fig. 1 for Ney York city, where tiles are the census block division of the urban area.We denote a graph by G = (T, E), where T is the set of n locations (nodes) and E is the set of trips (edges).The graph is directed (trip direction) and it allows self-edges and multiple edges (many travelers from and to the same origin and destination).Furthermore, the network is a labeled graph since the nodes will be specified by intrinsic location attributes.Finally, the trip mobility network will be defined on top of three fundamental assumptions which will be explained in the next paragraph.The first assumption defines the latent variable backbone of the mobility graph.

A. 1 (Latent variable graph)
Each location ℓ i is labeled with a node type by the latent variable x i , which represents the attractiveness of destination where a trip ends.It will be interpreted as the driver for the travel demand to each destination.Let (x i ) i∈[n] be a sequence of latent variables with values in a node-type space Ω x ⊆ R such that the empirical distribution of (x i ) i∈[n] approximates a probability measure µ x as n → ∞ [6,7].Consequently, the set of attractiveness variables {x 1 , . . ., x n } assocated to a location are considered realizations of independent and identically distributed latent random variables with an empirical distribution that converges almost surely to the cumulative distribution function F (x).
Let us observe that, in the paper, the latent-variable probability measure is defined on (R, B(R)) where B(R) is the Borel σ-algebra generated from the real line and the probability will be assumed absolutelycontinuous respect to the Lebesgue measure, so that the probability density function can be defined as ρ(x) = F ′ (x).Such condition will allow for various mathematical and statistical manipulations in the paper for practical applications.Similarly, the locations, other than being destinations, are, at the same time, labeled as origins.The correspondent intrinsic feature variable will represent the resident active population which stays in the area from which a trip originates1 .The latent variable in this case is named as productiveness of the location.Consequently, the set of productiveness variables {y 1 , . . ., y n } in nodetype space Ω y ⊆ R consists of realizations of independent and identically distributed latent random variables with probabilty measure µ y where a probability density function is defined ϕ(y), which can be different from the attractiveness distribution.From a dynamical perspective, the model is presented as random graph process that is a stochastic process that describes a random graph evolving in time [39,40].The processes studied here will have a fixed vertex set, and they will start without any edges and grow by adding edges according to linking rule, without deleting any, since the total number of visits up to a certain time is studied.Once a regional tessellation has been defined in terms of a graph as in A.1, the graph processes for directed multiedges can be defined [41] as random graph G t evolving so that ∀t a new edge is added ) is chosen randomly with replacement with probability proportional to a kernel function K(x, y) defined as Ω x × Ω y → [0, ∞).Such kernel is rules out the chances that a location of attractiveness x will be a destination of a visit whose trip originated in a location of feature y.Knowing that µ y (dy) represents the measure assigned to the infinitesimal interval of productiveness for origin locations.

A. 2 (Trip generation)
The trip arrival process is defined through the infinitesimal arrival intensity dV x = K(x, y)µ y (dy) which drives the dynamics of new travels landing in destinations of attractiveness x conditional to trips that departed from locations of infinitesimal productiveness y.
Such rate can be interpreted as the propensity (or intensity) for a traveler to move towards a destination of a given attractiveness x.So, the graph evolves when a new visit from an origin to a destination is completed (a new link in the network) according to the attraction and the production rates in the degreespace of the mobility graph.Similarly, as dual problem, a trip departure process can be defined by the infinitesimal departure intensity dV y = K(x, y)µ y (dy).
The evolving mobility graph {G (t) } is fully described by means of a time-varying adjacency matrix A (t) which represents the origin-destination table of the mobility problem at time t and the sum along the columns represent the in-degree of the nodes or equivalently the number k of visits received up to time t.Finally, an extension of the model is considered where typical lengths of single trips is associated to the corresponded visits.In the case of weighted network where each links is associated to a weight as a random variable independent from visitation process and the graph is expressed in terms of a weighted directed multigraph adjacency matrix Ã(t) = C (t) • A (t) , which indicates a element-wise matrix multiplication and C is considered either a coefficient matrix or a random one according to the real-world phenomena under study alongside the data that can be used.Specifically:

A. 3 (Trip weights)
All the arrivals of new trips are weighted in terms of some trip features.The weights r are independent and identically distributed random variables, with a specified probability density function ϱ x (r).
Those weights have the meaning of trip 'size' as for example the distance traveled from an origin to the selected destination, or the emission impact of each trip, or any type of travel cost or visit benefit.The weight probability ϱ x (r) may depend or may not on the attractiveness of the destination node according to the particular type of weight considered in the specific situation under investigation.Let us notice that the in-strength (i.e.weighted in-degree) of a destination node κ is defined as the sum of weights of all the visits received.In the case that weights are all equal to 1 the strength κ is equivalent to the degree k.At this point, finally, the model can be defined:

Definition 1
The trip mobility network is a temporal inhomogeneous random graph model that describes a visit generation process satisfying the assumptions A.1, A.2, A.3.
The weighted mobility multigraph process G t can be studied in terms of the adjacency matrix Ã(t) representation of the graph.A general framework which describes ensembles of dynamic networks can be assessed by making a Markov assumption on the evolution of the network by studying the probability of realization of a member of configuration ensemble graph over time [42,43].Specifically, the temporal evolution of finding the trip mobility network in the configuration A at time t after L steps (trips) can be written as: where Λ is any transition rule built on latent variables, which describes the creation, at each time t l , of a new trip from an origin towards a destination.From a modeling perspective, the complexity of knowing the configurational distribution of the origin-destination mobility network can be reduced trough the study of the degree (or strength) distribution and its higher order statistics.We actually model the degree distribution, that in the in-degree case results to be the visiting distribution defined as: where δ is the delta function, and P(A, t) is the probability to find our trip mobility network in the configuration A at time t, where each node has a strength degree κ i .If the in-degree of a node is simply the number of arrivals of new travelers, the in-strength of a location node is defined as the number of arrivals of new travelers who have faced a cost.Such weighted version of arrivals is named visits where each trip has an intensity.

Visit distribution
In the case of this work, the derivation of the visit distribution is assessed on the basis of the latent variable framework.The model is developed upon the asymptotic regime assumption of an infinite network where the number of locations in a tessellation is extremely large, where each location can generate and attract an unlimited number of edges (trips).Consequently, the continuous mean-field approach is used so that single locations can be studied as uncorrelated nodes in the same class of locations with the same attractiveness [14,15,44].So rather than studying the single location l i , one analyzes those locations which belong to the same class of attractiveness level, i.e. l x seen as a continuous random variable with probability density ρ(x).Let us call conditional visit distribution p(κ, t|x) the strength distribution conditional to destinations of the attractiveness type x, and each class evolves independently one from any another, so that a superposition of co-evolving conditional degree distributions is possible.The overall visit distribution can be derived according to the following proposition:

Proposition 1
According the visitation model's assumptions as in the Definition 1, the visit distribution of the mobility network is fully characterized by the attraction rate ν x that defines the transition probability, per unit of time, that a destination of attractiveness x increases its number of visits by one from any destination.The attraction rate is defined as the mean intensity of the trip arrival process ν x = Ωy dV x .
• The evolution of the conditional visit distribution can be described by a master equation for destinations with attractiveness x as an integro-differential Kolmogorov-Feller equation: with the initial condition p x (κ, 0) = δ(κ).In the asymptotic regime, the conditional probability p(κ, t|x) can be written as: , where where the correction factor ϵ x → 1 in the asymptotic limit t → ∞, and ⟨r⟩ x and ⟨r 2 ⟩ x are the first and the second moment of the trip-weight distribution ϱ x (r).
• Finally, the temporal asymptotic expression of the visit distribution is a mixture distribution as: where x 0,i = x 0,i (κ, t) are the zeros of the expression z(x) = κ − ν x ⟨r⟩ x t, and P (κ, t) the in-strength distribution of the overall trip mobility network.
Proof.Let us notice that the attraction rate νx is transition rate of new visits per unit of time interval, as the chance of having a new arrival in a destination of type x originated from any location, so that the attraction rate is the mean intensity as in [6,7,13,15]: The master equation for the evolution of the conditional probabilities px(κ, n) for locations with attractiveness x, where the step size is the correspondent ∆κ = r which represents the weight of each link i.e. the decision heuristic variable to move in the selected destination.It can be written as: If we also assume that px is slow, so that it changes only slightly during this time step τ = ∆t and redefine νx := νx/τ , then we can write the continuum master equation like: and ϱx(r) is the distribution of the trip distances covered to reach destination blocks of attractiveness x.At this point let us apply the Laplace transformation in the variable κ so L{p(k, t)} ≡ p(s, t) and the master equation transforms as: where the convolution product has been used.The solution can be written as: with the initial condition p(s, 0) = L{δ(k)} = e 0 = 1 so that c = 1.One can notice that the characteristic function is equivalent to the laplace transform as expressed before.It is possible to write the Laplace transform of the jump distribution ϱx(r) in terms of its moments as: If one assumes that ϱ(r) is a peaked distribution, following the central limit theorem rationale, one can assume it is described by the first two (finite) moments so that the solution in Eq.( 6) can be approximately: its inverse Laplace transform can be considered [45] the case asymptotic case of s ≪ ⟨r⟩x ⟨r 2 ⟩x of the truncated normal distribution as stated in the eq.( 2) For the second point, the visit (in-strength) distribution of the mobiliy network is expressed as a compound probability distribution that results from assuming that a random variable κ is distributed according to some parametrized distribution with the latent parameter x distributed according to some attractiveness distribution [46,47].So, the (unconditional) visiting in-strength distribution results from marginalizing the conditional distribution p(k, t|x) of the non-negative real-valued random variable x.So, the probability density function of the visiting distribution is given by the following the mixture density: Here, p(κ|x, t) is the distribution of κ when we know x at time t, in which the relation between κ and x can be seen as deterministic, i.e. κ = F (x, t) = Et[κ|x] = νx⟨r⟩xt, defining the distribution by its expected degree value through the momentgenerating function from the Laplace transform above.So, κ can only be single value, whose distribution is represented by dirac-delta functions δ(κ − F (x, t)).Consequently, the empirical visiting density probability function can be written as: July 29, 2023 which consists in approximating the visiting probability density by means of a Dirac mixture [48,49], where ρ(x) is the attractiveness probability density.Such procedure is equivalent to a change of variable respect to the deterministic one-to-one function in the static model as in [12].At this point we use the property: |∂z(x)/∂x| where x (i) 0 are the m-roots of z(x) = 0 where in the transport model z(x) = κ − νx⟨r⟩xt where z is a continuously differentiable function with z ′ nowhere zero.So: which represent a general formula for the tail behavior of the degree distribution for a mobility network with a generic attraction rate νx.
Let us observe that if the trip weight distribution ϱ x (r) is a dirac delta, then the strength distribution is equivalent to the degree distribution.Another particular case is when the trip weight distribution is identical over the attractiveness variable so that ϱ x (r) = ϱ(r) then the strength is proportional to the degree κ = ⟨r⟩k.In the next remark, the same results can be expressed in terms of stochastic process terminology rather than in terms of an integral-differential master equation.Let us observe that the type of processes described in Proposition 1 can be reinterpreted in terms of a mixture of compound Poisson processes as described in the Appendix which represents to be very popular in financial mathematics to model stock prices, insurance claims, and other financial phenomena [50,51].The visitation model of the mobility network can also be interpreted in combinatorial terms by using urn processes for solving balls in bins problems and finding the occupancy distributions as sketched in SM4.Such approach is used in world trade literature as in [52][53][54].Despite different approaches, the one introduced in the paper has the advantage to provided direct, though asymptotic, solutions to the scaling relations in the mobility networks.
As a particular choice of the occupation probability p, for every couple of origin-destination locations, trips can be realized according to a kernel function K(x, y) [12,[55][56][57].As a crucial case in the context of scale-free networks, one can consider the attraction rate to be proportional to some power of the destinations' attractiveness:
In the particular case that the attractiveness distribution is ρ(x) ∼ ρ 0 x −η , the visiting in-degree distribution has the following asymptotic tail distribution which shows the typical scale-free structure of an inverse power law distribution for the visiting degree of the mobility network.
The analytical results in the previous remark is confirmed by numerical integration of the compound distribution eq.( 3) by using the truncated normal conditional probability eq.( 2).Moreover, a graph process is performed via monte carlo (MC) simulation of the network where occupation probability is expressed via a separable linking function K(x, y) = g(x)h(y) for the evolution of the adjacency matrix.Such kernel gives arise to an attraction rate as ν x = ν 0 g(x) , where ν 0 is a normalization constant as shown in details in the supplementary materials.So by choosing g(x) = x α0 we are in the case as specified in Remark 1.At this point, it is possible to compare the three approaches and confirm the consistency of results obtained.In Fig. 2 it can be observed how the three approaches provide the same vist distribution for a particular choice of the parameters as discussed in the caption.As expressed in different research works [11,[57][58][59][60], different combinations of the attraction rate and attractiveness distribution can generate the same trip-visit distribution.In particular a scale-free behavior can be recovered trough exponential drivers where the attraction rate is in the form ν x = ν 0 e γx and the attractiveness distribution is ρ(x) ∼ ρ 0 e −λx .Also in this case the the asymptotic trip-visit distribution is scale free, in particular, via the transformation variable x 0 = x 0 (κ, t) = 1 γ log κ ν0t , the distribution becomes P (κ, t) ∼ t λ γ κ −(1+ λ γ ) .In particular, there is Figure 2: Visit distribution computed in the case of constant trip weights, and the attractiveness distribution is ρ(x) ∼ ρ 0 x −2 and kernel function K(x, y) ∝ x 1.5 h(y) so that νx = ν 0 x 1. 5 .In (a), the probability density P (κ) has been estimate with three different approaches: (1) it is evaluated through the numerical integration of the compound probability as in eq.( 3).( 2) It is evaluated through montecarlo (MC) graph simulation of sequential adjacency matrices with N = 300 locations and with a simulation time of t = 10 5 time steps.The three approaches provide the same scale-free behavior of the visit distribution as P (κ, t) ∼ t 2 3 κ − 5 3 as expected by the analytical asymptotic estimate eq.( 12).In (b) the compound distribution approach is calculated at three different time snapshots.a model selection issue, since different choices of attractiveness features can generate the same effect on the strength-distribution, one could clarify the ambiguity investigating higher order characterization of the degree distribution of the mobility network.

Visit correlations
A more detailed characterization concerns the exploration of the connectivity correlations in the origindestination correspondences of trip mobility network.Higher order statistics of a network in the degree space, can be obtained trough by the conditional probability P (κ (1) , κ (2) , . . ., κ (k) |κ ′ , t) that a node with strength κ ′ connects to nodes with strength κ (1) , k (2) , . . ., κ (k) at time t.The simplest of these degree correlations is the two-point correlation being described by the conditional probability P (κ|κ ′ , t) as the probability that a trip departing from an origin location of out-strength (departure) κ ′ reaches a destination node of instrength (visit) κ.The correlations between degrees of the nearest-neighbouring vertices are described by the probability distribution: However, the empirical evaluation of such conditional probability in real networks is cumbersome, so the weighted degree-degree correlations are commonly accounted by average-nearest-neighbor's strength function k nn (κ, t) which makes use of a smoothed conditional probability [61] often used as a measure of degree homophily of the nodes.In the latent variable framework, as shown in Fig. 3a, the conditional assortativity k nn (x) measures how much a location with attractiveness x tend to be a destination of an origin location of population y as defined in [14,62,63].In a similar way, the three-point correlations can be studied in terms of the clustering coefficient spectrum c(κ, t) which indicates the probability that two neighbors of strength-κ node are neighbors themselves.In the case of weighted and directed networks there many different ways to define the cluster coefficient [64,65].At the latent variable level, the conditional clustering coefficient of a destination with attractiveness x can be interpreted as the probability that two randomly chosen locations with trips towards a destination with attractiveness x are neighbors.Consequently, the Markovian property at the latent variable level [62,66,67] allows to calculate analytical expressions for the assortativity k nn (κ), quantifying two vertices correlations, and clustering coefficient spectrum c(κ), as a measure of three vertices correlation.A very important result is that the degree correlations of trip-visit distributions are completely determined by the attraction (and production) rate and by the origin-destination conditional probability χ(y|x).

Proposition 2
In the visiting mobility network under the latent variable assumption, the origin-destination correlation is defined as the conditional probability that a visit in the destination of attractiveness x has originated from a location of population y, and it is written as: As consequence, the following estimates of the two-point and three-point correlations hold: • the average out-strength of origin neighbors of destinations with in-strength κ, can be written as: If destinations and origins are independent χ(y|x) = χ(y) and ⟨k nn ⟩(κ, t) = const.
• the clustering coefficient for destinations of in-strength κ is: Proof.The conditional origin-destination probability is the conditional probability that a destination block of attractiveness x is connected to an origin block of population y is: , y) ∂y where Vx is the primitive function of the attraction rate νx.Io order to write the explicit expression of the origindestination correlation it is necessary to know the pairing rule, for example trough the connection kernel K(x, y) so that νx = ν 0 Ωy K(x, y)ϕ(y)dy, or imagining a generic primitive function for the attraction rate.So the conditional origindestination probability is an important indicator for correlations between origins and destinations in the degree-distribution of the visiting mobility network.The conditional average-nearest-neighbor's in-degree for destinations of attractiveness x can be written for directed multigraph in a continuous limit as in [12,13]: where in the mobility model the conditional expected strength is E[κ|y] ∝ νyt.Since, for the markovian degree property, the degree two point correlation can be fully determined by the conditional probability P (κ ′ |κ), the average degree of neighbors of an in-degree κ destination is known to be calculated as [14,63]: which is an in-out (origin-destination) assortativity measure, which is independent of κ for χ(y|x) = χ(y) as in the case of multiplicative separable linking function K(x, y).Let us notice that since the network is directed, other than weighted, one could define, similarly, other three average nearest neighbor's degree functions: destination-origin, origin-origin and destinationdestination.
The clustering coefficient of a destination with attractiveness x can be interpreted as the probability that two randomly chosen edges from x are origin-neighbors.The clustering of a destination of degree one or zero is defined as zero.In the space of latent variables, consider a destination i of attractiveness x i and population y i , which is connected with with probability p(y j , y k |x i ) through trips originated from two other locations j and k which have attractiveness x j and x k and population y j and y k respectively.Since the network is markovian at the latent variable level, p(y j , y k |x i ) = p(y j |x i )p(y k |x i ).Thus, similarly to the definition in [14,68] together with modifications [64,69] for the directed and weighted case, the local origins-destination clustering for locations of attractiveness x i can be written as: where p (x j , y j ), (x k , y k ) is the probability that the two origin nodes are connected one to the other in both directions.Now, in the asymptotic continuous regime the clustering coefficient can be rewritten: where M = ( K(x, y ′ )ϕ(y ′ )dy ′ K(x, y ′′ )ϕ(y ′′ )dy ′′ ) −1 , so it is possible to write: knowing that νy = ν 0 Ωx K(x, y)ρ(x)dx and the definition of χ(y|x).Let us notice that for independent origins and destinations then c(x) = c 0 = const.Moreover in the case of clustering coefficient in a multigraph one can calculate the number of triangles repeated κ times, which in a markovian graph can be approximated on average as ct(x) = tc(x), since at each time step a possible link is considered as a bernoullian trail, so that the observed number of links in t trials follows a binomial distribution and so the expected value is tK(x, y).Consequently, since the clustering coefficient has values in [0, 1], we normalize the adjacency matrix respect to t, so that the average local clustering coefficient of a node with strength κ, denoted by c(κ), [14,62,63,70] is given by: which represents the local in-clustering spectrum for destination locations and it is independent of κ for χ(y|x) = χ(x) as in the case of multiplicative separable linking function K(x, y), where P (κ, t) represents the in-strength distribution.Let us notice that in the case of the clustering coefficient the adjacency matrix and the latent variables are needed to be normalized in order to transform a multigraph in a weighted graph and from there a clustering coefficient no larger than 1 is guaranteed.
The Markovian nature of this class of networks implies that all higher-order correlations can be expressed as a function of the attraction and production rates ν x , ν y and the conditional origin-destination probability χ(y|x), allowing an exact treatment of mobility models at the mean-field level.Under the hypothesis that origins and destinations are independent, that is χ(y|x) = χ(y), then the average-nearest-neighbor's strength function and the clustering coefficient are constant along κ as shown, as an example, in the simulation shown in Fig. 4a and Fig. 4b.Consequently, for neutral networks the two and three point correlations can be obtained with three different approaches that will provide the same estimate by using the input approach of latent variables ('latent estimate'), by using the output approach trough the adjacency matrix ('expected value') and, finally, the algorithm computation of assortativity and clustering coefficients for directed and weighted networks ('simulation approach').

Remark 2
Under the hypothesis that visit production process and visit attraction process are independent the average in-strength of nearest neighbor function is constant, and in the asymptotic limit: As regard with the clustering coefficient spectrum, under the same hypothesis: Proof.The equations for the average in-strength of nearest neighbors and the in-clustering coefficient can be derived directly from Proposition 2 after some algebraic manipulations considering the production rate νy and the conditional probability χ(y|x) are independent from x as, in the case when the attraction and production rates are recovered from a linking function that is multiplicative separable, i.e.K(x, y) = g(x)h(y).In fact, when origins and destinations are independent, that is χ(y|x) = χ(y), then the average-nearest-neighbor's strength function knn(κ, t) = 1 + t νyχ(y)dy which is a constant over the strength degrees κ and where E[g(x)] and E[h(y)] are the expectation values of the function g(x) and h(y) under the circumstances they have finite values, which occurs even for fat-tail distributions in a finite set for the latent variables [13].In the case of infinite moments then the equation above represents a unreliable estimation but still shows the neutral assortativity of the graph.As regard with the clustering coefficient, the result in (2.3) is straightforward by using the multiplicative separable linking function as above, with the only difference that it has been normalized in order to provide a weighted matrix with link weights not larger that one.Another approaches for the estimates of the assortativity and clustering coefficient can be derived in terms of the strength degrees of the nodes as provided in [14,61,62] with proper modifications for directed weighted graph, see [13,71,72].Their expected values for knn(κ) and c(κ) are analytically known in literature for neutral networks, i.e. no degree correlations, as derived in terms of degrees the expected average in-strength of nearest neighbors can be written: where in the absence of correlations P (κout|κ in ) = koutP (κout)/⟨κ in ⟩ has been used.In the case of clustering coefficient, after normalization of the multigraph, the in-clustering coefficient can be written as in [62,73]: That in the case of uncorrelated networks: Let us notice the the clustering coefficient is normalized respect to time t since c(κ) ∈ [0, 1] and so it results to be the generalization for directed multigraphs without correlations as for simple graphs in [14,57,68].Simulations for such results are shown in Fig. 5 where several computational simulations of a graph for different time length t is presented alongside the prediction results for uncorrelated networks for assortativity and clustering coefficient.It is worth noticing that the analytical prediction are asymptotically valid so that no isolated nodes or leafs exists since local neighborhood clustering is typically not defined if a node has one or no neighbor and such situation influences the estimation of the global clustering in sparse networks [74].In the present work, the clustering coefficient algorithm removes all the local clustering of all the nodes with less than 2 neighbors, so the global clustering coefficient is over-estimated2 .

Results
In the present section a network analysis of the main graph measures and topology will be conveyed in the particular case study of New York metropolitan area by using Safegraph Mobility Dataset [75] for the year 2019.

Data
Origin-destination (OD) data represent movement flows through geographic space, from an origin (O) to a destination (D).OD datasets represent information on trips between two geographic areas often represented by the geographical centroids of the areas.Typically encoded with a square symmetric matrix, OD flow data contain numerical data on the aggregate quantity of individuals travelling from one geographic area to another over a specific time period.Mostly used in transportation planning, OD flows are an invaluable source of data for understanding spatial and temporal patterns of urban mobility and dynamics [2,[76][77][78].Visit flows can be in practice estimated in various ways in real world data.In particular, mobile phone location data are provided by SafeGraph trough dynamic population Origin-Destination flow matrices with hourly temporal resolution and aggregated by census block groups (CBG) in the USA as discussed in [37].In the daily CBG to CBG visitor flows metric, each row contains an origin CBG and a destination CBG, as well as the number of mobile phone-based visitor flows from the origin CBG to the destination CBG.Every day, the number of unique mobile phone users who live in the origin CBG and visits to the destination CBG are recorded.As regarding with visit production model, the population in each block is the key information to obtain from data in order to define the variable y and its respective distribution.However the population data is susceptible to the way data are collected and sampled by the provider.In fact the demographic sampling depends on many factors as the geographical boundaries which define a block, a tract, or any  administrative tessellation.Moreover in the statistical sampling methodology the individual measurements in each block go through a few transformations and aggregations which impacts the final measurement [79].
From SafeGraph data it is possible to build a matrix of trip flows between locations in a day for arbitrary large region R of the US.In particular, a county is considered with a tessellation at the resolution of census block group level, then it is possible to reconstruct the adjacency matrix of directed trips from an origin location towards a destination as stops of devices as described in [75].Fig. 6 plots a sequence of visitation counts in different census block group areas during different time windows of a day for New York city.
The data has been re-organized as shown in Table 1 where the (N + 1) × (N + 1) matrix A indicates the global origin-destination table visiting flow between the N blocks of the region R plus one external node which represents the resto of the world available in the data outseide the region of interes R. In particular, a ij is the number of visits registered as stop in the destination j in R that originated in the location i in R. For locations outside the selected region, the entry w i counts the visits in the destination i of the region R that originated from a location outside the region R.The in-degree for the destination location i is the sum along the columns of the row i from the matrix A from which the empirical visit distribution is evaluated.In addition, the trip-visit distribution, aka the in-strength distribution, is evaluated by associating a weight to each visit.As a typical choice, the weight is taken to be the distance between the origin block and the destination one in kilometers as estimate of the distance from home travelled by devices.Such information   is recovered by the census bureau geographical data using Census Block Group geometries with longitude and latitude coordinates of the block centroid [75,81], calculated as the haversine distance between the visitor's home geohash-7 and the destination location geohash-7 for each visit.A more detailed estimate would be the effective distance traveled by each visitor in the trip between its main location to the selected destination.Such information is not reported in SafeGraph at the moment neither in other data sources consistent with the data structure in the study.However, the radial movement approximation is motivated by the fact that travelers typically seek the shortest route [27,78].

Mobility network analysis
Let us start with the estimate of the distribution of visits among different destinations, which namely represent the in-strength distribution so that the in-strength of destinations is κ i = j C ij A ij where A ij is the entry of the adjacency matrix indicating the number of arrivals in destination i of a trip originated in the location j.Whereas C ij is the visit "size" of the traveler who departed from origin j and has arrived to destination i.Such value is taken from weight matrix C that represents, in this particular case, the distances between origin-destination pairs.Such cost matrices are directly recovered from census data included in the dataset used.In Fig. 7(a) the empirical complementary cumulative distribution function is plotted for the case of New York metropolitan area in a typical day of November 2019.The inspection of in-strength distribution shows that the visit distribution has a scale-free asymptotic behavior as P (κ) ∼ κ −µ with power law coefficient of µ ≈ 1.8.Let us now discuss the degree correlations such as assortativity and clustering just in the case of the origin-destination matrix.As already discussed, the average out-strength of neighbors of destinations of in-strength κ measures the tendency of having a directed trips from an origin location to a destination, defined as in [13,82], and from here the assortativity spectrum can be built.As plotted in Fig. 7(b), the average nearest neighbor in-strength function is flat, and this shows a neural assortativity behavior with a mean value of ⟨k nn (κ)⟩ κ ≈ 15.4.Such estimate is in agreement with the analytical prediction of the expected average in-strength of the nearest neighbor for uncorrelated networks out ⟩/⟨κ⟩ = 15.7 as proposed in Remark 2 .Under the same conditions, the average in-clustering coefficient can be computed accordingly to the definition [69,83] adapted to the in-clustering coefficient defined in the present work.The clustering spectrum for the data is shown in Fig. 7(c) where the clustering spectrum is flat as for uncorrelated-graphs, and the global clustering coefficient is given by ⟨c in (κ)⟩ ≈ 0.02 consistently with the analytical prediction of E[c (u) ] for uncorrelated networks as reported in Remark 2. This shows that the O-D SafeGraph mobility network is consistent with the hypothesis of uncorrelated graph with scale-free visit distribution at a macroscopic scale as also discussed in [2].The absence of degree-correlations allows to consider the origin-destination conditional probability to be χ(y|x) = χ(y).This means that a destination receives a visit from a randomly chosen origin location, without any particular choice of the origin but the number of resident populations in it.Consequently, the correlations between origins and destinations are entirely due to trip costs represented in the weight matrix is a well-mixed locations at level of attractiveness property as specified in the model assumptions.
The scale free behavior of the visit distribution together with the neutral tendency of degree correlations, allows us to narrow the type of kernel function to be considered in the model.Despite that, we do not have enough information to uniquely determine the attraction rate and the latent variable distribution.We will face such issue in the nex section when we will discuss possible proxies of attractiveness latent-variable.

Network scaling topology
A very important analysis of the mobility network, is to study how the trip costs affect the topology of network.The number of trip arrivals at the i-th destination can be written as the in-degree k i = j A ij , meanwhile the visit strength of the i-th destination can be written as κ i = j C ij A ij where C ij is the entry of the Origin-Destination distance matrix C. The in-degree and in-strength distribution have been plotted in Fig. 8a and Fig. 8b respectively, with a clear scale-free asymptotic behavior but with different power law coefficients.Such evidence suggests a very interesting aspect of visiting patterns where trip cost weights have a significant effect of on the mobility network structure.In particular, the relation between strength and degree of a location node can be written the average strength of destinations with degree k changes as: where the exponent δ represents the rescaling factor and δ = 0 occurs in the absence of correlations between the weight of links and the degree of nodes [82] so that the strength of a node is simply proportional to its degree and the two quantities provide therefore the same information on the system.The action of some correlation in the weight can bring cases where δ ̸ = 0.In such situation such relation induces a change in the scaling of the degree distribution (i.e.visit distribution) P (k) ∼ k −µ0 and the strength distribution (i.e.trip-visit distribution) P (κ) ∼ κ −µ according to the relation: and in the case of δ = −1, P (κ) is a Delta distribution of constant strength.In the case of the data under study, there is a clear linear relation between strength and degree as in Fig. 8c, consequently the strength distribution shows a scale free coefficient P (κ) ∼ κ −µ different from the one in the degree distribution P (k) ∼ k −µ0 consistently to the trasformation in eq.( 21).The value of the rescaling exponent δ has been estimated through a regression analysis which also allows to perform the significance level for a linear relation between degree and distance-strength.As regarding with the data used in this study, the weight of a links correspond to the physical distance between the origin and the destination of the trip which defines the trip-visit distribution as indicated in Fig. 3.3.As reported in the caption, there is a neat change of the slope in the scale-free distribution, the visit-distribution based on the degree fig.8a shows a slower power law coefficient respect to trip-visit distribution based on node strengths Fig. 8b.This is due to the linear relation between the degree and its weighted version trough the origin-destination distances.Let us notice that the weights have an impact on the degree which is significant in determining a change in the scaling relation of the distribution.Such result of δ > 0 suggests that the strength of nodes grows faster than their degree, in other words, the trip distances associated to highly visited locations have higher values than those expected if the trip distances were assigned at random.Such tendency denotes a strong correlation between the trip weight and the topological properties in the mobility network, where the higher the number of visits in a location, the more traffic the location can handle.As a conclusion, a general scheme arises where the number of visitors to any destinations decreases as the inverse power law of the product of their visiting frequency and travel distance, as already suggested by [27].On the contrary, random weights would have produced a δ close to zero, which occurs in the case when link weights are independent from the network topology, so that the strength distribution would carry no information than the degree distribution.For a more detailed estimate of the rescaling exponent δ, a regression analysis is reported in table Tab.2, for different types of trip weights.As a additional study to the linear regression statistics, performing a residual analysis makes it possible to test the assumption of a linear regression model such as the errors are independent and normally distributed, as shown in the Supporting Materials (SM3).

Economic applications: from land-use to travel demand
In this section, the investigation on possible properties of latent variables will be useful to select which kernel function is suitable to meet the trip mobility network characteristics.Later, the latent variable will be proved to be a crucial key that justify the income elasticity of a multi-purpose travel demand for any transportation mode.

The latent variable interpretation
The model dynamic is driven by the presence of latent-variables that are related to intrinsic properties of the areas but they cannot be not directly given (as for deterministic measurements) but they can only be inferred through statistical indicators (probabilistic measurement or proxy).At this point, it is possible to formulate the observed mobility network in terms of the latent-variable model so that the scale-free distribution shown in the real-world trip mobility network can be stated in terms of the latent variable statistical attributes.In particular, the attractiveness variable x and the productiveness variable y have been considered as hidden variables in the visitation model.The first summarizes the notion of the ability of a destination to attract visitors, the latter describes the ability of an origin location to produce and generate travelers with specific characteristics.However, latent variables can be interpreted more like proxies or indexes rather than proper measurements of observed phenomena.Under such perspective, there are several studies which try to estimate the number of visits by a combination of urban features, such as job opportunities, retail shops, business activities, infrastructural capacity, geographical positions etc...The standard approach to categorize urban areas classifies regions by their physical features and land use which refers to the way in which land is utilized, developed, and transformed for different purposes such as residential, commercial, industrial, and agricultural purposes.Urban morphology is seen as the result  of dynamic interactions between multiple factors, such as transportation efficiency, population size, and local land use.In particular, regional movement patterns, and consequently travelers' distribution, can be explained from land use, since purposes of people's trips are strongly correlated with the land use of the trip's origin and destination [22,32,77,84].For example, the latent-variable framework can provide an interesting interpretation of the effect of trip distances on the visit distribution, which can be formalized in terms of the attractiveness latent variable model.It is out of scope of this work to detect the best combination of factors which define the attractiveness of locations, however, how supported by some studies [21,24,85], let us take the non-residential land-use of census blocks to be the primary cause of travel demand and so it could be consider a proxy of destination attractiveness x in a multi-purpose travel model.In such perspective, the data from the New York Open Data [86] has been used where the land use zones is reported as the square feet occupied by building, parks and areas with a given use of destination (except residential), see the Supplementary Materials for more details on data used and correspondent interpretations.Land-use can be seen as a "ceteris paribus" candidate for attractiveness of locations and in Fig. 9a, the plot shows the probability density function of square feet of land-use lots and it shows a scale-free behavior with a power law exponent of 2 as plotted in Fig. 9a.The same analysis is performed after aggregating tax lots into census block groups, the land-use areas for non-residential purposes keep the same asymptotic fat-tail distribution with an inverse power law probability density function ρ(x) ∼ x −η with η ≈ 2 as plotted in Fig. 9b.
At this point it is possible to investigate the relation between mobility and land-use data (as attractiveness indicator).First, it is possible to analyze the relation between the number of arrivals k in each destination with the non-residential land-use of that area x.Similarly the relation of land-use versus the visits κ as weighted arrivals is analyzed as well.So, let us investigate a log-log linear regression analysis of such variables.The linear fit analysis of the relation k x ∼ x α0 reveals an estimated value α 0 ≈ 0.842 as shown in Fig. 10a and linear fit analysis of the relation κ x ∼ x α reveals an estimated value α ≈ 1.322 as shown in Fig. 10b.In table 2 such estimations for α 0 and α are reported for different types of trip weights.By knowing that κ x = E[κ|x] ∝ ⟨r⟩ x ν x ∼ x θ+α0 , it is possible estimate the scaling of the trip size goes as ⟨r⟩ x ∼ x θ , as confirmed by a direct linear fit analysis where ⟨r⟩ x ∝ x θ with θ ≈ 0.48.In such circumstances, by using the latent variable framework one can recover the visit distribution as P (κ) ∼ κ −(1+ η−1 α 0 +θ ) , which is the same scale-free distribution directly observed during the analysis of the origin-destination network.It is worth noticing the relation between the scaling exponents θ and δ.It can be easily verified that the strength distribution in eq (21) has µ = µ0+δ 1+δ = 1 + η−1 α , where µ 0 = 1 + η−1 α0 and α = α 0 + θ.Solving, we find that, the following relation holds: which can be checked by comparing the values reported in the Table 2 under their relative error margins.In conclusion, the attraction rate of a location is higher than another destination in the sense that the travelers are so motivated to travel a longer distance to get there.Such effect reveals the action of human travel demand on mobility, as formulated in the present paper, in terms of latent variable network model.The scaling relation between land use versus both travel distances and visits are power law like so that the attraction rate must be of a power law type and the distribution of land use that is the attractiveness latent variable has a pareto-type distribution.So no other attraction rate is compatible with the observation.Morever since degree correlation are absent, it is plausible that χ(y|x) = χ(y) so that the kernel function K(x, y) is multiplicative separable function.

Income elasticity
The visitation model approach can provide economical interpretations of some empirical urban scaling evidences [2], where scaling laws are also present in economical values of each location respect to its attractiveness.A crucial attention will be focused on the income elasticity of visit demand.Let us, first, define the "benefit" that travelers received by visiting a location, and in particular, the variable I x indicates the income level associated to the visitors who have traveled towards a given location with attractiveness x for a job purposes.In this way the benefit can be determined trough the strength-by-income variable κ i which converts a visit into potential economic output through a conversion factor i 0 that is the traveler's income per unit of visit time3 .On the other side, let us define the "cost" Q x faced by travelers to reach a location  as the strength-by-distance variable κ q taking into consideration that commuting costs are proportional to distance or time traveled, and the proportionality conversion factor c 0 is the cost of transportation per unit of quantity traveled 4 .Finally, the relation between the two variables can be written in terms of the attractiveness variable as: where θ Q is the scaling exponent derived from regression slope by income in Table 2 and the exponent θ from regression slope by amount of travel given by distance traveled in the same analysis.The relation between three variables is represented graphically in Fig. 11a.At this point it is possible to write the income elasticity of travel demand which has the meaning of how sensitive the demand for traveling a certain distance is to changes in income levels, the direct relation between the income reward and the quantity of travel demanded for visiting locations can be derived from eq.( 23) as: destinations are, again, randomly independent i 0 can be considered a mean-field constant correction factor.Further it can change over the observation time during which data have been collected.Considering i 0 to be location independent it can be considered a proportionality constant. 4Quantity can be measures in distance or travel time.A lower c 0 correspond to a better transportation system.Let us observed that in the case of cities the transportation costs and distance are well approximated by a linear relation.However for larger areas, empirical studies [87,88] show that transport cost is an increasing and concave function of distance so that c(r) = c 0 r v with 0 < v < 1, since travelers switch to faster transport modes for longer trips.
where the exponent ε is the income elasticity as estimated in Fig. 11, and in the case under study ε < 1, indicating that distance traveled is a necessity good, i.e., as income increases, people spend proportionally less on traveling when the income levels of travelers increase for any transportation mode 5 .Such prediction of an income elasticity of about ε ≈ 0.65 can be compared to results presented in literature reviews of empirical evidences and meta-analysis studies [35,[89][90][91] and also discussed with a theoretical interpretations [36,88], where average income elasticity for aggregated travel demand is estimated to be in a range values compatible with the elasticity ε estimated here 6 .It's important to notice that the income elasticity of travel demand can be influenced by many factors, such as the availability of transportation options, the price of transportation, and individual preferences.In conclusion, the in eq.( 23) shed lights on an intrinsic relation between the income elasticity of travel demand and the attractiveness scaling of urban areas.So urban planning policies, market conditions, and social dynamics emerges but also explains the interplay between attractiveness and economic through a mobility network flow which drives the travel behavior.From a network analysis perspective, it would be convenient to compare the number of trip arrivals at destinations against the travelers' income.Using eq.( 23) and eq.( 24), and visit elasticity of travel demand ε ′ can be defined through: where the income scales as respect to the number of trip arrivals.The visit elasticity of travel demand is a measure of how sensitive the number of trips to a certain destination is to changes in key travel attributes, such as fare level, service quality, journey time components, income and car ownership, and price of competing modes.The previous scaling relations of income elasticity, it is possible to range from the microeconomics perspective of transport economics to the macroeconomic outputs such as employment and the economic growth.For example, during a period of an economic growth (recession) where incomes are rising (falling) the distribution and the magnitudes of attractiveness can change and then modify the mobility pattern which, in its turn, has an impact on the global economic performance as well.

Conclusions and perspectives
In conclusion, the study presents a data-driven model for human mobility network based on an origindestination structure, which serves as the foundation for understanding mobility visitation flows.The model utilizes latent variables associated with each location, representing attractiveness and productiveness, to capture the intrinsic characteristics of destinations and origins.The contributions of this research are twofold.Firstly, it provides a theoretical framework that describes and reproduces a visit generation stochastic process through the use of intergral-differential equation for the evolution of the visit probability density and degree correlations in a trip mobility network.Consequently, analytical, numerical, and computational solutions are provided for important network characteristics, such as the strength distribution, assortativity, and clustering coefficient.A collateral impact of the visit generation model in travel dynamics is the introduction of a mathematical formalism of compound renewal processes commonly used in financial and actuarial science literature.Such stochastic models in a network perspective are useful for capturing the randomness and dependence between the arrival times and event sizes, providing a flexible framework for modeling and analyzing various mobility dynamics and financial and actuarial risks as well.The second contribution of this research has to do with the empirical analysis of real world phenomena.By analyzing origin-destination 5 For any mode and any purpose condition, undistinguished transportation mode is considered here, so all the possible modes are combined and only the distance necessary to reach the destination is taken into account. 6Let us notice that income elasticity shows a large variability in the empirical evidences since distance traveled is not homogeneous across different sources of income, type of jobs and age [34].Moreover, travel demand is reported in different units (individual or aggregate distance km/day, travel time, fuel consumption) and in different behaviors (commuting vs noncommuting, essential vs non-essential) or travel purposes (business, job, shopping, leisure).Travel can also vary in terms of different traveling modes according to transportation infrastructure.Moreover the estimations reported can even change over the years.In particular, the Q-I plane shows the projected regression analysis an elasticity of ε ≃ 0.76 consistently with the theoretical prediction from eq.( 24) and from eq.( 24) where the income Ix ∼ k 1+δ I /α 0 x .This indicates that a 10 per cent increase in income leads to about 6.5 percent increase in distance traveled or, conversely, a 10 percent increase in the demand of travel distance requires an increase of the income by about 15.4 percent.
data, the study reveals the presence of scale-free behaviors in visit frequencies and identifies correlations between visits and trip costs.The research also explores the statistical characteristics that latent variables should have to reproduce the observed patterns in the trip mobility network.Hence, the model permits to disentangle the effect of attractiveness (as land-use), population, trip costs and economic features of travelers on the visit dynamics in a mobility network.Clear scaling laws emerges between latent variables and travel demand.Finally, the model points out the effect of income on travel behavior depends strictly from the latent-variables that, therefore, can be considered as decision variable from a economic policy viewpoint.The possibility that human mobility belongs to the class of scale-free networks has impacted on the economic, engineering and mathematical communities in the multidisciplinary field of sustainable urban transportation, smart cities and world trade webs.As future outlooks, from a modeling side, the mathematical formalism discussed in the paper could also be extended to more general mobility graphs by considering more granular interactions on a time interval much shorter than the one used in the data explored here.For example, one can investigate more on the distribution of inter-arrival times of new trips to be non-Poissonian, and the the trip size are events with intensity that does not have finite moments (Levy-jumps).Another forthcoming study will be the study of emission mobility reductions for sustainable city planning.The approach can help the decision for the optimal allocation of economically attractive urban areas by transportation mode and land-use polices, and at the same time, minimizing the emission of pollutants and the impact of mobility trips among origin and destinations.From the data-side, moreover, a larger scale panel analysis is required for other regions to increase the robustness and the validity of the model.Moreover, the study shed lights on possible applications in fields like urban, transportation and environmental economics.These findings offer valuable insights to policymakers and urban planners, aiding in the comprehension and prediction of mobility patterns for informed decision-making and sustainable development.

Declaration of Competing Interest
The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Origin Destination mobility data has been retreived by SafeGraph for academic research program [75], and they cannot made public available.Social, demographic and geographical data are public accessible at [38,81,86,92].Codes for the main network analysis are povide in the Github repository [93].
(a) destination assortativity for the node x (b) destination clustering for the node (x, y)

Figure 3 :
Figure3: The two point and three point correlations of the trip mobility network can be calculated in terms of the hidden variables x and y.In particular the conditional average-nearest-neighbor's strength at the latent variable level (a) is the in-out strength (origin-destination) assortativity coefficient, and the cluster coefficient at the latent variable level (b) is the "In" clustering (or destination clustering) coefficient for weighted and directed networks as[64,65], i.e. a triangle such that there are two trips coming into of the destination node (x ← y ′ , x ← y ′′ , x ′ ← y ′′ ∨ x ′′ ← y ′ ).

Figure 4 :
Figure 4: The estimate of assortativity and transitivity of a latent variable network with the same structure as in Fig.2 with the additional specification of ϕ(y) ∼ ϕ 0 y −2 .The average-nearest-neighbor's strength function (a) which show a neutral assortativity in the network, the dashed line represents the assortativity mean value.The local clustering coefficient for different values of κ in (b), so that the transitivity is constant, the dashed line represent the average global cluster coefficient.
(a) mean in-assortativity for different simulation times (b) mean in-clustering coefficient for different simulation times

Figure 5 :
Figure 5: Two points and three point network correlations using three different approaches, in the case of uncorrelated graph as reported in Remark 2. The overall average strength of nearest neighbors ⟨knn(κ)⟩κ is replicated for each t so to obtain a mean global value ⟨knn⟩ over a ensemble of S = 50 replications as in (a).Similarly one can obtain the global mean in-clustering coefficient ⟨c⟩ in (b).

Figure 6 :
Figure 6: Number of cumulative visits in New York city [80] for different time windows.

Figure 7 :
Figure 7: Visit distribution and correlations of the Safegraph Neighbor Patterns mobility network for New York city in November 2019.The complementary cumulative in-strength distribution function is plotted in (a) of the trip-visit probability density function with a clear scale-free behavior.As for the trip-visit correlations, the network shows a neutral assortativity (b) and a neutral clustering (c) computed on the normalized version of the OD adjacency matrix.

Figure 8 :
Figure 8: Visit patterns for in-strength defined in terms of distance traveled in trips.The linear fitting of degree-strength scatter plot provide a significant scaling coefficient as κ ∼ k 1+0.53 , so that the visit distribution asymptotic behavior is P (k) ∼ k −2.2 meanwhile the trip-visit distribution by distance weights is P (κ) ∼ κ −1.8 .
July 29, 2023 (a) Tax lot square feet for non-residential land use (b) Tax lot square feet aggregated in census block groups

Figure 9 :
Figure 9: Log-binning procedure for the probability density function of square feet for land-use zones as hypothetical index for destination attractiveness.It shows an inverse power law fat-tail distribution so that the asymptotic behavior the probability density function can be written as ρ(x) ∼ x −2 .In the inset the whole distribution is plotted where the data is fully represented even at a low spatial scale.In (b) the complementary cumulative density function of the land-use square feet aggregated for census block from tax lot data.
(a) Scatter log-log plot of the destination with land-use x versus the number of arrivals kx to reach that destination.The slope is α0 ≈ 0.84.(b) Scatter log-log plot of the destination with land-use x versus the number of visits κx (i.e.number of arrivals weighted by distance traveled) to reach that destination.The slope is α ≈ 1.32.

Figure 10 :
Figure 10: Regression fit between the land-use of census block group and (a) arrivals and (b) visits.
(a) Three dimensional scatter plot for the variable {x, Ix, Qx}.(b) Scatter plot for the evaluation of the income elasticity of travel demand

Figure 11 :
Figure 11: Using the dataset, for New York in 2019, in (a) a 3D scatter plot shows the relation among attractiveness x as land-use, the income level of visitors and travel demand as distance traveled for each visit.The projection planes show three 2D scatter plots for pairs of the previous variables.In particular, the Q-I plane shows the projected regression analysis an elasticity of ε ≃ 0.76 consistently with the theoretical prediction from eq.(24) and from eq.(24) where the income Ix ∼ k

Table 1 :
Data matrix format.The vector D represents the array of locations as destination units in the tessellation region.Similarly O represents the same array locations but as origins of trips inside the tessellation region.The array W D is the set of destinations located outside the region and W O is the set of all the origin locations outside the region.So A 0 is the open O-D table, and A is the close O-D table.

Table 2 :
In-degree vs in-strength regression analysis for different trip-weights.The slope of linear fit reported is the coefficient 1 + δ.