# Temporally Consistent Present Population from Mobile Network Signaling Data for Official Statistics

Milena Suarez Castillo, Francois Sémécurbe, Cezary Ziemlicki, Haixuan

Xavier Tao, Tom Seimandi

► To cite this version:

Milena Suarez Castillo, Francois Sémécurbe, Cezary Ziemlicki, Haixuan Xavier Tao, Tom Seimandi. Temporally Consistent Present Population from Mobile Network Signaling Data for Official Statistics. Journal of Official Statistics, 2023, 39 (4), pp.535-570. 10.2478/jos-2023-0025. hal-05024940

HAL Id: hal-05024940

https://hal.science/hal-05024940v1

Submitted on 8 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Temporally Consistent Present Population from Mobile Network Signaling Data for Official Statistics

Milena Suarez Castillo<fnref n="1" />, Francois Sémécurbe<fnref n="1" />, Cezary Ziemlicki<fnref n="2" />,
Haixuan Xavier Tao<fnref n="1" />, and Tom Seimandi<fnref n="1" />

Mobile network data records are promising for measuring temporal changes in present populations. This promise has been boosted since high-frequency passively-collected signaling data became available. Its temporal event rate is considerably higher than that of Call Detail Records – on which most of the previous literature is based. Yet, we show it remains a challenge to produce statistics consistent over time, robust to changes in the “measuring instruments” and conveying spatial uncertainty to the end user. In this article, we propose a methodology to estimate – consistently over several months – hourly population presence over France based on signaling data spatially merged with fine-grained official population counts. We draw particular attention to consistency at several spatial scales and over time and to spatial mapping reflecting spatial accuracy. We compare the results with external references and discuss the challenges which remain. We argue data fusion approaches between fine-grained official statistics data sets and mobile network data, spatially merged to preserve privacy, are promising for future methodologies.

**Key words:** Big data; High frequency statistics; dynamic population mapping; spatial accuracy.

## 1. Introduction

Mobile network data has the potential to significantly enhance population statistics by increasing their levels of spatiotemporal details and their timeliness. Several international initiatives aimed at incorporating this new data source into the production of official statistics, such as the United Nations Big Data dedicated Task Team (U.N. Global Working Group 2019), the European Statistical System working groups (ESSnet Big Data 2021), in addition to national initiatives (Statistics Netherlands 2020; Coudin et al. 2021). Yet, despite the new interest in timely measuring dynamic population during the Covid crisis, producing statistics from mobile network data has remained a great challenge for

<footnote n="1">¹INSEE, SSP Lab, 88 avenue Verdier, Montrouge, Île-de-France, 92120, France. Email: milena.suarez-castillo@sante.gouv.fr, francois semecurbe@agriculture.gouv.fr, tao.xavier@outlook.com, and tom.seimandi@insee.fr</footnote>

<footnote n="2">²Orange Labs R&D Châtillon, Chatillon, Île-de-France, France. Email: cezary.ziemlicki@orange.com

**Acknowledgments:** Since 2016, INSEE, Eurostat and Orange Labs have been collaborating on exploring the usefulness of Mobile Phone Data for official statistics. In this context, we have been able to benefit from the work of a French national collaborative research project (ANR Cancan†) which collected three months of mobile signaling data in 2019. The raw data was deleted after twelve months. The approach described in this document only required exchanges of anonymous aggregates between INSEE and ORANGE. Processing of individual data was performed by each data owner. However, methods and algorithms starting from the raw data were developed jointly and transparently allowing both parties to evaluate and validate the outputs. We are grateful to Elise Coudin, Mathilde Poulhes, Etienne Come, Angelo Furno, Latifa Oukhellou and Zbigniew Smoreda for insightful comments, and to European Statistical System collegues for great discussions, in particular Fabio Ricciato, David Salgado and May Offermans. We also thank the audience at the BDES 2020 and at the NTTS 2021 conference. We also thank Romain Avouac and Vincent Attia who first explored the data, and the TIC survey team.</footnote>
official statistics. Lack of access to real data, privacy protection and models to efficiently cooperate with private data holders have been long discussed challenges. The need to invest in a transparent methodology may be one of the distinctive features of official statistics relative to private producers who already disseminate statistical products worldwide.

Historically, a large body of research provided insights on human mobility by relying on Call Detail Records (CDR), data generated by the user when actively communicating through their cell phone (texting, calling...) and collected for billing purposes (Blondel et al. 2015). The reliance on user activity and the low time event rate has proven an obstacle for inferring from CDR present populations fitted for official statistics purposes (Sakarovitch et al. 2018; Vanhoof et al. 2018). Today in France and generally in Europe, private actors producing population statistics from their networks instead use re-purposed signaling data as base material when available. This new generation of mobile network data is characterized by a much greater spatiotemporal event rate. To ensure the centralized network knows about the mobile's state and location to reach it efficiently, passive communications are pervasive with the device when switched on, and are increasingly collected by Mobile Network Operators (MNOs) for network optimization and monitoring purposes. For statistical purposes, both types of data are however similar in the unitary information they contain: records convey *non-continuous* information on the proximate *radio cell* ($\approx$ antennas) mobile *devices* are wirelessly communicating with. Thus, three main dimensions of uncertainty should be addressed to infer from MNOs records information on population presence: temporal, spatial, and population coverage uncertainties (Ricciato et al. 2020).

To be disseminated as official statistics, dynamic population counts should be steadily comparable over time, that is, as insensitive as possible to network-related and user behavioral effects or MNOs client turnover. They should also be consistent with other approaches to estimate population counts when other high-quality data sets are available. Indeed, if many case studies providing dynamic maps of mobile phone users exist, focusing only on mobile phone users with-out any extrapolation to the total population is of limited interest for official population statistics (Panczak et al. 2020). Mobile network data represent an opportunity for developing countries where alternative data sets are scarce and high-quality census data are not always available. In these countries, static population estimates can be vital for development and infrastructure planning (Salat et al. 2020). In developed countries, dynamic population counts, within typical days or over the year, are not part of official statistics production today. Yet population is a highly temporally dynamic variable, with considerable shifts in its distribution occurring in daily cycles. Today the majority of published indices are monthly or quarterly and correspond to the night-time population distribution. Hence, all applied sciences and policy support that require spatially detailed information on population distribution must rely on a fractional and static representation of reality (Batista e Silva et al. 2020). This is a major issue for planning of infrastructure and transport for instance, as well as assessment of human exposure to natural, environmental, epidemiological, and technological hazards (Freire and Aubrecht 2012; Panczak et al. 2020).

The literature on addressing spatial uncertainty is insightful. Location at radio-cell level is imprecise and depends on the network's local density. Voronoi tessellations are
commonly used models for radio-cell coverage as they require limited information (the cell tower physical coordinates) and are easy to implement by partitioning space. Yet, they are imprecise and do not correspond to the way the network functions (Sakarovitch et al. 2018; Ricciato et al. 2017). Probabilistic approaches accounting for overlapping coverage of cells are more realistic and thus preferable in a context where the mapping choice entails large discrepancies in outputs (Tennekes and Gootzen 2022; Ricciato et al. 2020; Ricciato and Coluccia 2021; Salgado et al. 2021). As for temporal uncertainty, interpolation techniques were explored for mobility analysis when targeted time granularity is high or when working with sparse Call Detail Records data (Hoteit et al. 2016; Chen et al. 2019; Bonnetain et al. 2021). Typically, CDR records a few events per device per day. In turn, signaling data ensures a very large detection, even at night – but the literature is still in its infancy due to lack of access. The availability of signaling data may allow for simpler strategies for the use case of dynamic population measurement and a reassessment of this dimension. Finally, only a specific fraction of the population is observed, and thus scaling is needed. To link mobile network data to other data sources, such as census data, home detection algorithms aim at identifying where the user lives, to rescale mobile-phone counts by residency locations, for example, as in Fekih et al. (2021). In the absence of a home detection step, cruder alternatives exist, using counts of mobile-phone users at night or census counts as rescaling factors. Home detection can be considered a data augmentation technique and rescaling can be stratified across more device characteristics to improve representativeness. Data fusion approaches – ensuring coherency across several sources with distinct strengths and weaknesses – can enhance mobile network data using additional sources. In particular, combining other data sources also allows to estimate daytime and nighttime populations. For instance, Batista e Silva et al. (2020) combine official statistics (residents, employees, students, tourist, and other population counts and estimated flows between region) with geospatial data to produce and validate a European data set of population grids taking into account intraday and monthly population variations, called “ENACT” (Schiavina et al. 2020).

In this article, we focus on estimating hourly population counts at a fine spatial scale over several months over France's metropolitan territory. Our approach is based on spatially merging passively-collected cellular network signaling data with fine-grained official population counts by place of residence. The structure of this article is as follows. In Section 2, we define the theoretical framework and draw particular attention to various dimensions of uncertainty observed in the data. In Section 3, we describe our estimation method, driven by consistency constraints at several spatial scales and over time. We require that the dynamic population counts stratified by living places are locally consistent with the official number of residents. This requires estimating a residency for each device and allows to break down dynamic population counts by place of residence. We focus on building dynamic flows for French residents only (more precisely residents of metropolitan France) and exclude from the onset the existence of inbound and outbound trips, absent reliable sources on daily population flows in and out of the country. We design the spatial mapping of population counts to reflect spatial accuracy by building an adapted quadtree grid – independent of administrative borders. In Section 4, we compare the results with the ENACT dynamic population counts. Finally, in Section 5, we discuss challenges that remain and argue that data fusion approaches between fine-grained official
statistics data sets and mobile network data, spatially merged to preserve privacy, are promising for future methodologies. We review how this work relates to recent initiatives in the European Statistical System to build a reference methodological framework for mobile network data integration into official statistics production (Salgado et al. 2021; Ricciato et al. 2020; Statistics Netherlands 2020).

## 2. Statistics of Interest, Theoretical Framework, and Challenges Deserving Particular Attention

### 2.1. Data

Three months of raw signaling data from Orange clients were collected from March 16th to June 15th, 2019. These data included all Orange client device interactions (active, such as text or calls, and passive, such as hand-overs between antennas) with the Orange metropolitan France 2G, 3G, and 4G networks. All personal information was removed and device identifiers were pseudonymized before storing the data, which were erased twelve months after collection. Data are collected through probes positioned on the Orange networks for monitoring performance. Events are information exchanges between the devices and radio cells, base station antennas which communicate wirelessly with mobile devices, hereafter, cell or radio cell. Users of the 4G networks generate on average about a thousand events per day, while users of the 3G and 2G networks generate respectively about 50 and 20 events per day.

Orange provides two types of radio cell location information. First, weekly extractions of a cell registry covering the data collection period were performed. The cell registry exists for network maintenance purposes and sees regular entries and exits following network life. The extracted data were limited to cell tower coordinates and cell technology type. Second, Orange Fluxvision provided a cell-specific coverage map, which modelizes the network coverage as of February 2019. The latter is static and is obtained from a radio-propagation model taking into account network specificities, local topography and land use, and device diversities through simulations. Applications for these data include helping to provide emergency call locations.

The official source which serves as a reference for localized population counts is fiscal data from 2016 (Filosofi for Fichier Localisé Social et Fiscal), obtained from merging fiscal income tax files (at the individual and household levels) and local residence tax files. It includes households with at least one income declaration and liable for the local residence tax (all owners, tenants, or occupants on a free basis of a piece of real estate). The situation is evaluated at the address on the first of January. Filosofi excludes the homeless population and people living in collective housing (retirement homes, young worker homes, prison, university residences, etc.). The native individual data is geolocalized at the tax address and can be aggregated over our regular grid of interest covering metropolitan France. As this data records tax addresses, it tends to be of lower quality to determine residency for young adults when they are fiscally attached to the home address of their parents. Yet, it is the most granular source of resident population localization available for France (Lamarche and Lollivier 2021). We note that the forthcoming gridded populations from the EU 2021 censuses could be used in future development as an alternative.
Finally, the TIC survey (Enquête sur les technologies de l'information et de la communication) is a survey conducted by Insee every year, with the objective to collect information on the usage of information and communication technology in households. In 2021, 31,000 different households were interrogated across metropolitan and overseas France (excluding Mayotte). A question was specifically added asking for the MNO of the surveyed individuals. This allows us to get a sense of how representative Orange subscribers are of French residents in 2021, which is close to the study period.

We provide the definitions of the populations we use to construct our estimates in Table 1.

## 2.2. Notations

Let us introduce the notations summarized in Table 2. The device set $\{d \in D\}$ is defined as all devices (as identified by their International Mobile Subscriber Identity) appearing at least 30 distinct days over Orange 2G, 3G, and 4G networks within the three-month time window $T$, identified as an Orange client (based on the Mobile Network Code of their IMSI) and as a mobile phone (based on their Type Allocation Code and an external register of known mobile phones TAC). The hourly time grid of interest $\{t \in T\}$ contains each hour in the three consecutive months. $P$ is the total population of interest, which is assumed constant over $T$. $D_t$ denotes devices observed during $t$ and $D_t$ devices observed at date $l$. The grid $\{i \in I\}$ on which we want to estimate present populations is made up of tiles, and $\{j \in J\}$ is the set of cells in the MNO network.

We would be in an ideal setting for statistical purposes if the number of observed devices were constantly equal to the size of the target population, that is $D_t \Leftrightarrow P$. In this case, a count of all devices in location $i$ at hour $t$ would provide our target statistics $u_{i,t}$, the

<table><caption>Table 1. Population definitions by data source.</caption><thead><tr><th></th><th>Source</th><th>Period (location)</th><th>Population</th></tr></thead><tbody><tr><td>Residents <i>P</i></td><td>Filosofi (tax data)</td><td>01-01-2016</td><td>Households with an income declaration and liable for residency taxation</td></tr><tr><td>Devices <i>D</i></td><td>Orange Probes</td><td>16/03 – 15/06/2019</td><td>Users of a mobile phone, with an Orange contract, interacting at least 30 distinct days with the MNO network</td></tr></tbody></table>

<table><caption>Table 2. Notations.</caption><tr><td><i>d</i> &isin; <i>D</i></td><td>Devices in the scope, detected on the network over the period of interest</td></tr><tr><td><i>t</i> &isin; <i>T</i></td><td>Hourly time grid of interest (several weeks)</td></tr><tr><td><i>P</i></td><td>The target population set</td></tr><tr><td><i>D</i><sub>t</sub></td><td>Set of devices in the scope, detected on the network during <i>t</i></td></tr><tr><td><i>D</i><sub>l</sub></td><td>Set of devices in the scope, detected on the network at date <i>l</i></td></tr><tr><td><i>i</i> &isin; <i>I</i></td><td>Tiles that cover the territory of interest</td></tr><tr><td><i>j</i> &isin; <i>J</i></td><td>Cells of the MNO network</td></tr><tr><td><i>u</i><sub><i>i</i>,<i>t</i></sub></td><td>Population count in location <i>i</i> at time <i>t</i> (target statistics)</td></tr></table>
present population for this location and hour. In fact, several dimensions of uncertainty in the data must be dealt with to compute valid estimates for present populations. We detail these dimensions now.

## 2.3. *Temporal Uncertainty: Only Counting Active Mobile Devices Leads to Unreasonable Variations in Aggregates*

The first uncertainty observed in the data is the high temporal variation in counts of active devices. This variation confounds many mechanisms unrelated to population variation and is actually a major stylized fact of mobile network data. The presence of a given device may be very sporadic in the collected data, both within and across days. Although we improve time event rate by relying on both passive and active data, this remains a crucial issue when estimating population variation.

Within a given day, the behaviors of mobile phone users and network coverage are the main reasons for the observed variation: users may choose to shut down their phones, run out of battery, or out of signal. Network coverage is uneven across space – and can not be ignored in a country such as France where low-density areas still host a large fraction of the population. Across days, we may expect a large role for client turnover, the telecom market in France being competitive and changing operators being increasingly easy for customers. For instance, over the first semester of 2018, four million mobile phone numbers were kept while changing MNO. Telco market dynamics are also at play: the sim cards number quarterly growth was more than 200,000 in the second quarter of 2019, over our period of interest (for all MNOs, excluding M2M). We also have to take into account failure in data collection of the richer passive data – as probes may punctually malfunction, without a too high cost for the MNO as it does not affect the network communications but is rather a monitoring tool. Investing in their continuous reliability is not a priority. In addition to these mechanisms entailing undesirable user disappearance, users may disappear due to outbound trips – but for now, they are indistinguishable from the former, although of interest for present population estimation. Characterizing places where outbound trips originate (airports, train stations, borders) could be considered to discriminate absence from the territory from other phenomena.

Figure 1 illustrates the issue on the data set. The first panel represents the ratio of observed devices $\frac{|D_i|}{|D|}$ for each date $i$. Over a long period, aggregate variations in the percentage of users observed on a given date vary considerably: by several percentage points regularly, occasionally by more than ten percentage points. Underlying causes of the aggregate variations remain ultimately speculative (reported data collection failure, outbound trips e.g., on bank holidays, etc.). While we can detect relatively easily unreasonable variations at the aggregate level, we may suspect that the same issues go undetected at local levels. Within a day, the hourly detection rate varies between 70 and almost 90% (second panel of Figure 1), a rate considerably higher than for CDR data. For instance, Orange CDR data from 2007 recorded only a few events per device per day and a percentage of observed devices which goes from less than five at nighttime to about 50% in the late afternoon (Galiana et al. 2020a). Variations seem highly driven by devices disconnecting at nighttime. Although the detection rate is arguably rather high, only counting active mobile devices would lead to unreasonable temporal variations in aggregates.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Proportion of observed devices by date and hour of the day.</figcaption></figure>

Note: Left: % among all Orange devices. We distinguish all users (dotted line) and users present for at least 30 days. We represent dates with reported data collection issues (in red) and within-week bank holidays (in blue). Right: % among devices appearing from March 16th to March 31st, 2019. Scope: Orange metropolitan France network, Orange-client devices identified as mobile phones.

A simple method to deal with this source of inconsistency is to account for all devices across the entire time grid by adopting a device-centric view and interpolating device trajectories when unobserved – building a panel of device trajectories. We detail this method in Section 3.

## 2.4. Spatial Uncertainty: Mapping Presence Over the Network in Space

Recorded events are located at the level of radio cells, not on the grid of interest *I* on which we estimate present populations. Spatial information on devices over time must be mapped onto the grid of interest. In general, radio cells overlap which means that a mobile phone in tile *i* at time *t* can potentially be detected in several cells, as pictured in Figure 2.

Let us define coverage probability matrix *A*, such that $A_{ji}$ represents the probability of being detected at cell *j* while being in tile *i*:

$$A_{ji} = \mathbb{P}\{\text{device detected in cell } j \mid \text{device in tile } i\}.$$

Orange Fluxvision provided an instantiated cell-specific coverage map matrix *A*, which modelizes the network and devices population as of February 2019. Note that without this

<figure><img src="image_3.png" /><figcaption>Fig. 2. Illustration of the relationship between tiles and radio cells.</figcaption></figure>

Note: The grid of interest is made up of tiles represented by dotted squares. In this illustration two radio cells A and B are pictured. In general a cell covers more than one single tile, and cells overlap, meaning that a single tile is covered by more than one cell.
information coming from the MNO itself one can use a Voronoï tessellation of the cell plan to model radio-cell coverage or design a more complex model like in Tennekes and Gootzen (2022). Figure 3 illustrates how the coverage probability matrix encodes the device connexion information from tiles to cells, and the spatial mapping task which aim at locating devices connected to cells into tiles. Note: Left: 4 tiles are represented by rectangles and 3 cells by ellipses. A mobile phone in tile *i* has probability $A_{ji}$ to be detected in cell *j*. For instance, the mobile phone in tile 2 has probability $A_{12}$ to be detected in cell 1 and probability $A_{22}$ to be detected in cell 2, with $A_{12} + A_{22} = 1$. Right: The phone detected in cell 2 is mapped to the tile grid using spatial mapping *Q*. It has probability $Q_{22}$ to be located in tile 2 and probability $Q_{32}$ to be located in tile 3. Using a uniform prior, we have $Q_{22} = A_{22}/(A_{22} + A_{23})$ and $Q_{32} = 1 - Q_{22}$.

The dimensionality of matrix *A* is high when considering the entire French territory (approximately 55 million tiles with sides of 100 meters for several hundred thousand cells). If we take this matrix as a good approximation of the ground truth, radio cells are massively overlapping with each other over the territory and the signal is quite diluted, in the sense that the cell-tiles mapping is highly non-exclusive. On average over France in the coverage map, there are 33 cells with positive coverage per 100-meter tile. However, most of these links are really low: if we restrict to links with significant coverage (say $A_{ji} > 0.1$) there are only 2.4 cells per tile and 16% of tiles without any linked cell. The first panel of Figure 4 illustrates the distribution of areas covered (resp. significantly covered) by cells. The median cell covers 25 tiles significantly (1706 tiles with non-zero coverage). The second panel of Figure 4 illustrates the number of cells per tile. The median tile is covered by 30 distinct cells, but significantly by two cells.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Coverage probability matrix and spatial mapping.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Distributions of cells coverage areas and of cells per tiles.
Note: Here, a tile *i* is said covered (significantly covered) by a cell *j* if $A_{ji} > 0$ ($A_{ji} > 0.1$).</figcaption></figure>
In practice, the location accuracy we can expect when mapping in space an event located at a given cell relies on the extent of the cell-covered area. If taken as the ground truth, the estimation A encodes a relatively low network precision which is highly heterogeneous over space. The method to map cells to tiles based on A as precisely as possible while managing dimensionality is detailed in Section 3.

## 2.5. Coverage Uncertainty: Mapping Devices to the Population

The third dimension of uncertainty in the data is linked to population coverage. By definition, we do not observe the target population (all French residents) but a selected subset through their device(s). It is a subset to the extent that we have a reliable method to exclude devices which are not used by human beings (M2M, IoT), and that the remaining Orange devices identified as mobile phones indeed belong to French residents who carry a single Orange device along. We maintain this assumption throughout this work, although some devices may not be carried out by French residents in the fiscal sense. Present population estimates relying on device counts are consistent under the assumption that the mobility and presence patterns we observe for the selected subset of the target population can be extrapolated to the target population.

However, active Orange devices are not exactly representative of the target population. Orange Fluxvision estimated that the national market share of Orange was around 37% (29% for individual consumers) in 2020 (Insee 2020). The TIC survey from 2021 indicates that biases first include an over-representation of the elderly (from 55 to 84 years old) amongst Orange subscribers in comparison with the general population, and an under-representation of younger individuals (in particular from 15 to 24) as shown in Figure 5. The youngest individuals interrogated in the TIC survey are 15 years old: children younger than 15 are also naturally under-represented in our signaling data, as many do not possess a mobile phone. Individuals belonging to households with lower incomes are also less likely to be Orange mobile users. Finally, non-French citizens living in France are also under-represented amongst Orange subscribers.

A milder assumption can be formulated when some characteristics of the devices are observed: mobility and presence patterns of unobserved residents can be extrapolated from mobility and presence patterns of persons carrying an Orange mobile device when they share these characteristics. High-frequency signaling data fittingly allows observing one such characteristic – the home environment.

A related question is then whether the places of residence of Orange clients are representative of places of residence in France. Locally, we may derive a ratio between residents and MNO-detected residents. We expect that the lower this ratio is, the better the local representativity, as we observe more devices per resident. This ratio informs on differential local representativity but also divulgates Orange market shares and is thus not reported on a map. In practice, it is highly heterogeneous over space. The local representativity tends to deteriorate in poor neighborhoods in some urban areas. As an illustration, Figure 6 represents how this ratio distribution evolves by municipality median disposable income. At the municipality level, there is no clear correlation between disposable income and local representativity as captured by the ratio between residents and MNO-detected residents, except at the lower hand of disposable incomes. If the
<figure><img src="image_6.png" /><figcaption>Fig. 5. Probability of being an Orange subscriber as a function of age, income quintile and nationality on the TIC survey data, modeled with logistic regression.</figcaption></figure>

Note: The odd ratios minus 1 as percentages are displayed. Reference modalities are “from 45 to 54” for the age variable, “2” for the income quintile and “French” for the nationality.

median ratio is about 0.33 detected resident per resident, the first decile (D1) of the ratio is 0.16 while the ninth decile (D9) is 0.58. A large heterogeneity may also exist at lower spatial scales.

On top of users’ socio-economic characteristics which may differ from one MNO to another, mobile network data better captures the behavior of active users who benefit from good network coverage. In this case, raw mobile network activity data under-represents the population from less dense areas where the network is less developed.

### 3. Measuring Present Population

We now present the methodology we use to compute hourly present population estimates over metropolitan France using three months of signaling data from the Orange network and the geography of French residents from fiscal data. Formally, we aim at computing estimates $\hat{u}_{i,t}$
<figure><img src="image_7.png" /><figcaption>Fig. 6. Detected residents per actual residents and Municipality-level Disposable Income.
Note: Distribution (D1, Median, D9) of the municipality-level ratio by municipality median disposable income range.</figcaption></figure>

giving the distribution of France metropolitan residents at hours *t* in *T* over several months and over tiles *i* ∈ *I*. As an intermediate output, we aim at estimating $\hat{u}_{i,t,r}$ giving the distribution of France metropolitan residents in location *r* present in tile *i* at hour *t*. Thus, $\hat{u}_{i,t} = \sum_{r} \hat{u}_{i,t,r}$.

**Minimal consistency constraints.** For consistency, we require that the present population estimate total matches an external official source: $\forall t \in T, \sum_{i} \hat{u}_{i,t} = P$. We will further require that we have as many residents of *r* contributing to $\hat{u}$ as there are residents of *r* in the official source to balance our estimates across residencies: $\forall t \in T, \sum_{i} \hat{u}_{i,t,r} = P_r$. This requires estimating a residency for each device. It allows to break down population presence by place of residence. We note that we exclude from the onset the existence of inbound and outbound trips, due to the absence of reliable sources on daily population flows in and out of the country.

As we do not restrict the analysis to a sub-period or part of the territory, we favored a simple method to deal with the high dimensionality of the data.

<table><caption>Table 3. Additional notations.</caption><tr><td>$j_{d,t} \in J$</td><td>Estimated presence cell for device d during time interval t</td></tr><tr><td>$P_r$</td><td>Official residents in place r</td></tr><tr><td>$D_r$</td><td>Devices with estimated residency in r</td></tr><tr><td>$\mathbf{m}_{r,t} \in \mathbb{N}^{[J]}$</td><td>Devices count in presence cells at time t with residency r</td></tr><tr><td>$\mathbf{u}^0 \in \mathbb{N}^{[I]}$</td><td>Population count from official source over tiles</td></tr><tr><td>$\hat{\mathbf{u}}_{i,t} \in \mathbb{R}^{[I]}$</td><td>Estimated population count at time t over tiles</td></tr><tr><td>$\hat{\mathbf{u}}_{i,t} \in \mathbb{E}^{[I]}$</td><td>Estimated population count at time t who are resident in r over tiles</td></tr></table>
### 3.1. Overview of the Method

Our population estimation relies on several modules, the critical ones being the construction of a device presence panel and residency-based weighting.

**Device presence panel.** This first module's role is to bypass the temporally sporadic presence of users over the network, by interpolating the trajectory of each device before any aggregation. In practice, we define and estimate for each hour and each device a presence cell $j_{d,t}$ – whether or not the device was observed during t. The presence cell $j_{d,t}$ represents the cell where the device d would mostly connect during the time interval t if it was active.

**Residency characterization.** This second module's role is to estimate the residency r of each device d at an adapted geographical level to be defined. We denote $D_r$ the set of devices with residency in place r, which can be compared to the set of residents in the official source, $P_r$.

We can then define presence over cells of devices residing in r, denoted $\mathbf{m}_{r,t} \in \mathbb{R}^{|J|}$ with

$$m_{j,r,t} = \sum_{d \in D_r} 1\{j_{d,t} = j\}.$$

$m_{j,r,t}$ is the count of devices which are resident in r and are considered present in cell $j \in J$.

**Residency-based weighting.** This third module's role is to extrapolate the number of devices to an estimate of the present population. If we assume that the sample $D_r$ has been randomly chosen among $P_r$ with sampling rate $|D_r|/|P_r| = \frac{1}{w_r}$, a valid estimation of the expected presence of residents of r over the network cells is $w_r \times \mathbf{m}_{r,t}$. We extrapolate the presence patterns of $D_r$ to $P_r$. It amounts to applying a rescaling factor $|P_r|$ to the density of resident devices. The pseudo-weight $w_r = |P_r|/|D_r|$ is the ratio of residents from the official sources to the resident devices in place r. Instead of counting for 1 person, each device in the scope will participate in the aggregate with weight $w_d = w_{r(d)}$. Of course, this approach is valid if $D_r$ is indeed close to a random sample drawn from $P_r$. If we could add additional inferred characteristics to the devices, we could improve this stage by stratifying weights beyond residency.

**Spatial Mapping.** This fourth module's role is to transform an active device at the cell level into an active device at the tile level. We do it by defining a linear spatial mapping $Q: \mathbb{R}^{|J|} \rightarrow \mathbb{R}^{|J|}$ which distributes a vector of presence over the network cells in the tiles with $\sum_i Q_{ij} = 1$, by specifying

$$Q_{ij} = \mathbb{P}\{\text{device mapped to tile } i \mid \text{device connected to cell } j\}.$$

Then, $\mathbf{Q}\mathbf{m}_{r,t} \in \mathbb{R}^{|J|}$ represents the presence over tiles of devices which reside in r, at time interval t.

**Presence Estimation.** We estimate the presence over tiles of the residents of place r denoted $\hat{\mathbf{u}}_{r,t}$ by attributing to residents of r the presence distribution of devices who are resident in r. That is,

$$\hat{\mathbf{u}}_{t,r} = w_r \mathbf{Q} \mathbf{m}_{r,t} \quad (1)$$

We finally estimate the total population presence over tiles with

$$\hat{\mathbf{u}}_t = \sum_j \hat{\mathbf{u}}_{t,j}$$
These definitions enforce our minimal consistency constraints. Note that $\hat{u}_t$ can be written as a reweighted projection of device-level trajectories:

$$ \hat{\mathbf{u}}_{i,t} = \sum_{j \in J} Q_{ij} \sum_{d \in D} w_d 1\{j_{d,t} = j\} \quad (2) $$

## 3.2. Implementation

The raw signaling data contain about 20 billion events per date, totaling 130 Tb of parquet files. We handled them using the big data framework Spark on an HDFS infrastructure at the MNO office. The Spark cluster was configured to stop any job lasting more than 24 hours. Given that the cluster was shared with other projects, the estimated resources available for our project were 300 CPUs and 1.2 Tb of RAM. We queried the data through PySpark, the Python API for Spark. We did not use any algorithm that could not be simply implemented using the PySpark API functions on the raw signaling data, which limited device-level algorithms. Even with simple algorithms, the whole-network longitudinal analysis entailing device-level sorting over a long period was challenging given the available resources.

### 3.2.1. Device-level Simplifications to Manage Dimensionality

Given these constraints, the following simplifications drove the choice of the method presented in Subsection 3.1 and were adopted in the implementation:

1. **For device presence, the location information is restricted to one cell per hour.** It seems a reasonable simplification for official statistics low-frequency purposes (at most, presence per hour). It stabilizes by design the oscillation phenomenon by which a motionless device may switch cells for network-related reasons (Katsikouli et al. 2019). However, if it is probably a good approximation for motionless devices, it is not for non-stationary devices visiting a high number of distant cells during an hour. Therefore, we oversimplify the presence of non-stationary devices by attributing them to one cell on their path.

2. **The location information was kept at the cell level for all device-level calculations.** We were only allowed to export aggregates to the national statistical institute premises, as opposed to device-level data. Hence, the spatial mapping was performed on aggregates only. This avoids bottleneck operations on $|D| \times |I|$-sized matrices during calculations (such as operations involving A). On one hand, this leaves the spatial mapping from cells to tiles easily adjustable downstream and leaves room for comparison between several choices of spatial mapping, which has been shown to matter a lot (Ricciato et al. 2020). On the other hand, as the spatial links between cells are never considered at the device level, some loss of spatial information is likely. In particular, the location information was kept at the cell level for device-level residency characterization. To bridge the tiles resident counts and the cells resident devices so as to build weights with counts in the same geography, we can map residents to cells or resident devices to tiles. Instead of affecting to each device a home tile i, we distribute each resident in a home cell j, the cell they would have probabilistically connected to if they were Orange users, using the probability distribution encoded in A. This choice avoids resorting to a home-cell spatial mapping $\mathbb{R}^{|J|} \to \mathbb{R}^{|I|}$ which requires dealing with
an additional probabilistic error. This approach is not straightforward in a multi-operator scenario where there is in general a specific radio cell network for each MNO.

While these simplifications are acceptable given our constraints, they are not necessarily recommended for future standard methodologies. Figure 7 summarizes the different steps of the method now described in more depth.

### 3.2.2. Device Presence Panel

We first filter devices identified as mobile phones (to filter M2M). Note that about 20% of daily unique identifiers are filtered as they can not be identified as mobile phones based on their TAC. We expect a large part of these excluded devices to be carried out by machines rather than persons (M2M and IoT devices, such as cameras, vehicles, alarms, sensors, etc.). We retain devices that are present at least 30 distinct days out of the three months so as to ensure relative stability of the scope (e.g., to filter movements due to client turnover as they are irrelevant to inform on total counts). Here, a user is present on a given day if there is at least one radio-cell event associated with its IMSI on this day. Figure 18 in the Appendix (Section 7) gives the distribution of the lengths of observations (in days) over the period. We keep 87% of devices when filtering those present at least 30 days out of the three months. 67% of the devices are observed for 80 days or more, and 9% are observed for ten days or less.

We attribute to each device appearing on date *l* a presence cell $j_{d,t}$ for each hour interval from 5 a.m. on date *l* to 5 a.m. the next day *l* + 1 as follows. When the device appears during an hour interval *t*, the presence cell at *t* is taken as the cell recording the most events during this hour. When the device does not appear during an hour interval *t*, we consider as candidates all cells with an event from the device on hour *t'* from 0 a.m. on date *l* to 5 a.m. the next day. We extended the search window to 29 hours for helping interpolation in the early morning with nighttime observations. We rank the candidate cells for hour *t* first by the time interval distance between hour *t* and the hour of the observed event *t'* and

<figure><img src="image_8.png" /><figcaption>Fig. 7. Overview of the method implementation.</figcaption></figure>

Note: Red boxes represent source data sets, dark blue boxes represent produced data sets and the key variables. Arrows represent transforms, which pertain to a given module referenced in light blue boxes.
second by the number of signalling events generated by the device at $t'$. We retain the cell ranked first. This cell selection based on the number of events is also a simplification as generating a higher number of events in a certain cell for a given hour does not automatically indicate the device spent more time in that cell than in others. One could choose the presence cell using more complex rules based on the distribution of signaling events on a sub-hour level.

Figure 8 presents the mean interpolation error, defined as the absolute difference between the reference hour *t* and the hour of actual observation, used to estimate the *presence cell*. It is on average lower than one hour but varies following the daily user behaviors and technology. It is particularly high when the presence cell is a 2G or 3G cell.

### 3.2.3. Residency Characterization

We can draw a picture of the typical day (including weekends and holidays) of (over-15) persons in France using time use surveys. On average in 2010, 8:30 hours are spent sleeping, 1:02 hours for washing/health care, 2:13 hours for eating, 4:04 hours for leisure (including 2:06 hours spent watching tv), 3:10 hours for domestic work, 2:51 hours working or studying and 0:24 minutes of work-home commute (Ricroch and Roumier 2011). This average includes unemployed, retirees, and housewives and may be heterogeneous across population types. However, for the “average” person, the vast majority of the time is spent at home.

For the characterization of residency, we therefore choose the cell with the most time spent rather than favoring a heuristic using time spent at night, although we run it as an alternative. Observation at night suffers from undersampling (Figure 1). In addition, we do not want to assume or constrain where the population is at night but rather deduce it from the data. For instance, we note that 1.8 million employees work at nighttime (8 p.m. to 5 a.m.) more than half of their working hours a given month (Létroublon and Daniel 2018).

To estimate the residency of each device, we thus want to identify recurring points of presence. A device is present at cell *j* in hour *t* as soon as an event is recorded at this cell during the time interval *t*. For this task, a device is therefore counted as present in all cells which have detected it at some point – even for a single event. This is therefore distinct from the presence cell as defined to track the longitudinal presence of devices. Here, all

<figure><img src="image_9.png" /><figcaption>Fig. 8. Active devices and Interpolation in Time.</figcaption></figure>

Note: The right panel presents the proportion of active devices among all devices in dotted lines, in total and by cell technology. The left panel presents the mean interpolation error defined as the absolute difference between the reference hour *t* and the hour of actual observation.
events are used and unobserved periods are not inferred. Then, the max-presence cell is the cell where the device has been recorded present the most. For simplicity in computation and scalability, we kept the analysis at the cell level but note that this approximation could probably be improved by considering several cells and their geography. It turns out that we find evidence of at least one strong “anchor point” for most of the devices in the scope. Figure 9 represents the distribution of the number of distinct hours of presence in the max-presence cell when the latter is defined over two weeks. 75% of the devices in the scope are observed at least 27% of the hours over the period in the same max-presence cell. Overall, signaling data prove very promising for pinpointing anchor points. Note that for this step, we did not interpolate device trajectories over an unobserved time period. If we define residency as the place with the most time spent, surely longitudinal signaling data offer many perspectives.

This step requires us to use all events and to sort the longitudinal data by device pseud identifier to be able to rank cells. To derive the max-presence cell over three months, we run a max-presence cell algorithm by two-weeks windows and kept per device only ten max-presence cell candidates. Filtering the least likely candidate cells allowed to keep the computational burden manageable. Finally, we define the home cell as the max-presence cell over the pooled max-cell candidates. At this stage, this step is highly stylized from a methodological point of view but benefits from the richness of the data over a long period.

Note that the international definition of population recommended in official statistics is based on the concept of ‘usual residence’, which is the place of daily rest assessed over a period of at least 12 months. We only had three months of data for this work, but the same methodology could essentially be applied with a longer observation window to produce comparable population estimates, consistent with the international definition of population.

<figure><img src="image_10.png" /><figcaption>Fig. 9. *Distinct hours of presence in the max-presence cell.*</figcaption></figure>

Note: The max-presence cell is defined as the cell recording the highest number of distinct hours of presence over a two-week time period. One event within the hour is enough to consider presence in this cell at hour *t*. Scope: Orange metropolitan France network, Orange-client devices identified as mobile phones, from March 16th to March 31st, 2019 (384 hours).
### 3.2.4. Residency-Based Weighting

In this step, we map French residents over home cells using realistic information on the coverage of each tile of 100 meters by Orange cells as provided by Orange (matrix A). Let us denote $\mathbf{u}^0$ the resident counts over grid I estimated from 2016 fiscal data. If all French residents were Orange clients, active and at home, we would expect to observe over the Orange network cells the following counts:

$$ \mathbf{Au}^0 $$

We define places of residency $r$ as home cells or groups of contiguous home cells gathering at least 20 detected resident devices. This allows us to adapt our definition of residency to the MNOs' varying local market shares, having enough devices per residency place while keeping the place of residence relatively precise. We start from all home cells, search for the closest home cells for home cells with less than 20 detected resident devices, group both, and iterate until all groups of contiguous home cells reach the condition. This ensures a minimal size for $D_r$ while keeping a high level of disaggregation in $r$. In turn, $P_r = \sum_{j \in r} [\mathbf{Au}^0]_j$ is the expected number of residents in the home-cell group $r$. Figure 10 presents at level $r$ (group of contiguous home cells) and at the municipality level the number of detected residents against actual residents, hinting at weight heterogeneity.

We define weights with the ratio of actual residents $P_r$ divided by the network-detected residents $D_r$ at the level of contiguous groups of home cells $r$ which contain at least twenty detected resident devices. We end by trimming weights at their 2nd and 98th percentiles – weights fall in [0.07, 53.5].

Ideally, $D_r$ and $P_r$ should be consistently defined and both tend to match the concept of “usual residence” to avoid bias.

### 3.2.5. Spatial Mapping

We use the coverage probability matrix A provided by Orange Fluxvision that modelizes the network as of February 2019. $A_{ji}$ represents the probability of being detected at cell j while being in tile i:

$$ A_{ji} = \mathbb{P}\{\text{device detected in cell } j \mid \text{device in tile } i\}. $$

<figure><img src="image_11.png" /><figcaption>Fig. 10. Detected residents and actual residents, by aggregation level.</figcaption></figure>
In particular $\sum_j A_{ji} = 1$. From a vector of presence over all tiles $\mathbf{u}_t$, we expect to observe on network cells $\mathbb{E}[\mathbf{m}_t] = A\mathbf{u}_t$ translating the presence of devices (we here assume $D \Leftrightarrow P$ for clarity of exposition). The estimate $\hat{\mathbf{u}}_t$ can be written in general as $\hat{\mathbf{u}}_t = g(A, \mathbf{m}_t)$ where g is a chosen spatial mapping (Ricciato et al. 2020). In this work, we focus on a linear estimator $\hat{\mathbf{u}}_t = Q\mathbf{m}_t$ (see Figure 3). Although any spatial mapping could be used, for our empirical results we follow Tennekes and Gootzen (2022) who suggest deducing $Q$ from A using Bayes' rule by introducing a prior that reflects where the population is most likely located (e.g., based on land-use). In the results presented here, we use a simpler uniform prior. Specifically,

$$Q_{ij} = \frac{A_{ji}}{\sum_{j'} A_{ji'}}$$

In addition, we propose a general framework to evaluate the location estimation accuracy of cellular network events. This evaluation combined with a quadtree algorithm enables us to build an adaptive spatial grid featuring small tiles for high accuracy areas and large tiles for low accuracy areas. The spatial accuracy is embedded within the dissemination grid.

**Estimating accuracy locally.** The accuracy of the linear estimator $Q$ can be approached locally by defining the probability to localize in $i$ a device that is in $i_0$ and connects to the network probabilistically through A.

$$N_{i,i_0} = \mathbb{P}\{\text{device mapped to tile } i | \text{device in tile } i_0\}$$

Formally, $N = QA$. A good estimator $Q$ should lead to a high $N_{i_0,i_0}$ probability (correct mapping), or at least a high probability of tiles $i$ in the neighborhood of $i_0$. With previous notations, if we take the example of localizing a single device $d$ which is in $i_0$, that is $\mathbf{u}_t = 1_{i_0}$, $N1_{i_0}$ can be interpreted as $\mathbb{E}[\hat{\mathbf{u}}_t | \text{device in tile } i_0]$. $N$ encodes the spatial error by integrating the uncertainty from A and $Q$.

**Embedding accuracy within dissemination.** We build a quadtree that directly embeds the calculated spatial accuracy by gathering tiles until the probability of correct location in the macro tile $I_0$ (group of tiles) is higher than a threshold: $N_{i_0,i_0} > s$. We derive present population estimates within this reduced spatial grid, which visually provide a clear idea of the achievable accuracy (Figure 11). In what follows, the tile grid $i \in I$ should be understood as this quadtree-derived grid for $s = 5\%$, referred to as the dissemination grid. Table 11 describes the link between tile size and resident density in further detail.

### 3.2.6. Presence Estimation (Aggregation)

In practice, the set of devices present a given day varies (Figure 1). We denote this set $D_l \subset D$ for each date $l$. We only interpolate device-level trajectories within days. To respect our consistency constraint, we finally define $w_{d,l} = w_{r(d)} \times \frac{|D_r|}{|D_r \cap D_l|}$. If all detected residents of $r$ are here on date $l$, their weights are set to $w_r$. If some are missing, their mass is transferred to the remaining residents' devices.

Precisely, our final estimator writes, for the chosen spatial mapping $Q$:

$$\hat{u}_{it} = \sum_{j} Q_{ij} \sum_{d \in D_l} w_{d,l} 1\{j_{d,t} = j\}$$

To illustrate the difference with counts of active devices, we reproduce Figure 8 by reweighting each device with $w_{d,l}$ in Figure 12. The mean interpolation error increases –
<figure><img src="image_12.png" /><figcaption>Fig. 11. Reduced grid with a threshold s = 1%. The larger the tiles, the less accurate the spatial mapping is locally.</figcaption></figure>

Note: The grid was based on the Orange matrix A and a uniform prior.

<figure><img src="image_13.png" /><figcaption>Fig. 12. Estimated present population and interpolation in time.</figcaption></figure>

Note: The right panel presents the proportion of estimated present residents among all residents, in total (by assumption, all residents are represented) and by cell technology. The left panel presents the weighted mean interpolation error defined as the absolute difference between the reference hour *t* and the hour of actual observation. Weights are *w*_{d,t}.

showing that devices with high interpolation error have been upweighted. This suggests that this weighting scheme corrects for the under-representation of users with limited network access or limited network interactions. By construction, the total number of present residents is constant. Estimated presence over the 3G network is now comparable to estimated presence over the 4G network: $\sum_{d \in D_t} w_{d,l} 1\{j_{d,t} = j \& j \in 3G\} \approx \sum_{d \in D_t} w_{d,l} 1\{j_{d,t} = j \& j \in 4G\}$.

### 3.2.7. Code Availability

The code implementing these modules is available on github at https://github.com/In-seeFrLab/presentpop.

## 4. Results and Comparison with External Sources

The clear advantages of present population estimates derived from mobile network data are their timeliness, their granularity in time, and (relatively) in space. We first provide a rapid overview of the hourly and weekly fine-grained patterns which can be uncovered. This dynamic nature and this spatial extent are rarely achievable with other sources. We
then compare some present population estimate snapshots to other external, more static,
sources of population density.

4.1. *Daily and Weekly Cycle, Local and National Variations*

Figure 13 illustrates the daily variations of present population density at the national level and in Paris. To represent the 24-hour *within-day* variations, we subtract from the hourly density its 24 hours mean. The 24-hour variation features the daily pendulum movement of suburban commuters: while the present population tends to be higher in the periphery of urban areas at night (≥ 100 persons per square kilometer compared to average), at 9 a.m. the population in these peripheries decreases for the benefit of urban centers until to the evening. At a finer scale in the Paris surroundings, the variation of present population density discriminates places mainly characterized by their economic, leisure, and touristic activities from more residential areas, and shows the attractiveness of a multi-polarized center.

<figure><img src="image_14.png" /><figcaption>Fig. 13. Variation of present population density (persons per square kilometer) at daytime and nighttime at national level and in Paris.</figcaption></figure>

Note: To represent the 24-hour *within-day* variations, we subtract from the hourly density its 24 hours mean (Wednesday, March the 20th).
Figure 14 illustrates the weekly variations of present population density at the national level and in Paris. To represent the seven days *within-week* variations, we subtract from the present population density at 3 to 4 a.m. its seven days average. The variation of nighttime present population density within the week discriminates the nights from Friday to Saturday and from Saturday to Sunday, where some locations in coastal areas and the mountains fill up. In Paris, the present population is overall lower during the nights from Friday to Saturday and Saturday to Sunday. However, we observe some nighttime excess in the present population in some places on these nights, probably reflecting nightlife activity or touristic overnight stays. The animated and colored version of Figures 13 and 14 is available at https://github.com/InseeFrLab/presentpop.

## 4.2. Comparison with External Sources

Measuring the quality of present population statistics is difficult as there is no source of truth. But, we can assess how comparable our estimates are to other high-quality

<figure><img src="image_15.png" /><figcaption>Fig. 14. Variation of present population density (persons per square kilometer), Tuesday or Saturday at nighttime, at national level and in Paris.</figcaption></figure>

Note: To represent the seven days *within-week* variations, we subtract from the present population density at 3 to 4 a.m. its seven days average (Monday 18th to the Sunday 24th of March 2019).
population measures. We here choose two points of comparison: residents geolocalized at their tax address and day and nighttime population density estimation from Batista e Silva et al. (2020).

We have considered comparisons with estimates of other external actors like Meta’s high-resolution population density maps and WorldPop population distributions. Indeed, the general idea of the methodology behind these estimations is – just like in our case – to calibrate population counts originating from new data sources (in the case of Meta and WorldPop, mostly satellite imagery) with spatially detailed population census data. However, we believe that a comparison with each of these two data sets is not entirely relevant. Indeed, their estimates are static and do not vary as a function of time, which is precisely the aspect where mobile phone data make a difference. The main interest of these static estimates lies in countries in which fine-grained census or fiscal data does not exist at a national level. This is not the case for the methodology outlined in this article, designed to estimate dynamic present populations. We still include a static comparison with the resident populations $u^0$ computed according to geolocalized fiscal data from 2016, which should be more accurate than Meta or WorldPop estimates in a country like France. Note that we also used this data source to build our weights.

The second comparison source is the ENACT database which uses data fusion to estimate the population at daytime and nighttime at the European level on a 1km grid. Batista e Silva et al. (2020) use a top-down approach disaggregating NUTS3 population counts based on the assumed places of activities of 16 different population groups, using land use information. In contrast to our estimation, foreign tourists’ presence is estimated and contributes to the population density.

We compare population densities in our dissemination grid. Figure 15 presents maps over France and Figure 16 a focus on the Paris area. At the country scale, the present population estimate respects the distribution of the population found in the other sources. However, the present population estimate shows a dilution of the population mass in space. Around dense urban areas, we observe a halo of presence absent from other sources. The first reason for this is probably the lack of accuracy of the cell-level localization. Another reason could be that by definition, the external sources considered here locate the population in buildings. Adopting a land-use prior in the spatial mapping task (Avouac et al. 2019) could therefore partially close the gap between our estimate and external sources and in any case improve its quality. However, this may create bias, in particular during daytime and weekends. We here chose a static spatial mapping, that is, independent of $t$. Another way forward could be to experiment with other estimators $g$, where $\hat{u}_t = g(A, m_t)$. On simulated data, Ricciato and Coluccia (2021) show that some estimators (called data-first, ML/EM) may be able to reduce the dilution effect compared to our linear spatial mapping with a uniform prior (Simple Bayes-rule estimator in the terminology of the article).

At the level of the Paris area, the structure of the present population distribution differs strongly during the day from during the night. The present population at 3 a.m. on a weekday (f) tends to offer a smoothed but quite similar version to the high-resolution image of the resident population (h). We report a contextual map of the Paris region in the Appendix (Section 7). For instance, the large parks in the east and the west (Boulogne and Vincennes) display a non-null density at nighttime according to our present population estimates, most likely due to a lack of accuracy in the spatial mapping. The population
<figure><img src="image_16.png" /><figcaption>Fig. 15. Population Densities in the Dissemination Grid (a-d). Present Population: March 2019 (week day). ENACT: March 2011. Resident Population: 2016.</figcaption></figure>

variation from nighttime to daytime is similar if we consider either present populations as estimated from mobile network data from (f) to (e) or from disaggregation of NUTS3 official sources counts as obtained in ENACT data, from (j) to (i). For instance, the core center near the Seine river and the Défense neighborhood attract population during the day, as predicted from activity-related presence with the ENACT methodology. Figure 20 in the Appendix (Section 7) presents the differences between the present population statistics and both external sources over France. It makes the tendency of mobile network data estimates to create halos around dense areas clearer. It shows that at night in dense areas such as Paris, except for particular places such as parks, the error resembles a white noise (no tendencies over space to either overestimate or under-estimate the population compared to the resident population). During the day, it tends to offer even more contrast between places' densities than the ENACT estimates.

In addition, Table 4 reports comparisons along three metrics: correlation, rank correlation and allocation accuracy. Allocation accuracy can be interpreted as the percentage of population density allocated in the same tiles in both sources and is defined as:

$$AA(\rho^1, \rho^0) = 1 - \sum_i \frac{\frac{1}{2} \times |\rho_i^0 - \rho_i^1|}{\sum_i \rho_i^0}.$$
<figure><img src="image_17.png" /><figcaption>Fig. 16. Population Densities in the Dissemination Grid (e-h). Present Population: March 2019 (week day). ENACT: March 2011, daytime (original grid 1km² in (i-j)). Resident Population: 2016.</figcaption></figure>

The present population at nighttime is as close to the resident population as is the ENACT estimation during the night (slightly closer according to correlation and allocation accuracy – and farther according to rank correlation). Both present population measures get more distant from the resident population during the day, although the night/day difference is more pronounced in the MNO-derived estimation. Finally, the MNO-derived presence estimation and the ENACT presence estimation are closer to each other during both day and night than they are to the resident population.

Overall, at nighttime, MNO and ENACT densities fall in the same metrics range when compared with the resident population and are aligned with each other. If we had used a prior based on land uses for spatial mapping *Q*, we would probably be even closer to the ENACT estimation – which by definition follows land use.

Table 5 reports correlations with external sources by urban area type. Consistency between sources is highest in urban areas of intermediate size (between 200,000 and
<table><caption>Table 4. Population density comparisons. Note: All densities are computed in our dissemination grid, using proportional area estimation for ENACT data and direct calculations for resident density from tax files. ENACT: March 2011. Present population: five working days in March 2019. Resident population: 2016. Daytime (resp. nighttime) estimates are taken at 3 p.m. (resp. 3 a.m.).</caption><thead><tr><th></th><th>Allocation accuracy</th><th>Correlation</th><th>Rank correlation</th></tr></thead><tbody><tr><td><strong>Residents</strong></td><td></td><td></td><td></td></tr><tr><td>Present at daytime</td><td>0.57</td><td>0.56</td><td>0.73</td></tr><tr><td>Present at nighttime</td><td>0.74</td><td>0.83</td><td>0.75</td></tr><tr><td>ENACT – Day</td><td>0.66</td><td>0.77</td><td>0.78</td></tr><tr><td>ENACT – Night</td><td>0.71</td><td>0.78</td><td>0.81</td></tr><tr><td><strong>ENACT – Day</strong></td><td></td><td></td><td></td></tr><tr><td>Present at daytime</td><td>0.74</td><td>0.75</td><td>0.89</td></tr><tr><td><strong>ENACT – Night</strong></td><td></td><td></td><td></td></tr><tr><td>Present at nighttime</td><td>0.79</td><td>0.87</td><td>0.88</td></tr></tbody></table>

<table><caption>Table 5. Population density correlations by urban area type. Note: All densities are computed in our dissemination grid, using proportional area estimation for ENACT data and direct calculations for resident density from tax files. ENACT: March 2011. Present Population: March 2019 (five working days). Resident Population: 2016. Daytime (resp. nighttime) estimates are taken at 3 p.m. (resp. 3 a.m.).</caption><thead><tr><th></th><th>Paris urban area</th><th>Other urban areas ≥ 700,000 inhabitants</th><th>Urban areas 200,000–700,000 inhabitants</th></tr></thead><tbody><tr><td><strong>Residents</strong></td><td></td><td></td><td></td></tr><tr><td>Present at daytime</td><td>0.42</td><td>0.62</td><td>0.74</td></tr><tr><td>Present at nighttime</td><td>0.77</td><td>0.81</td><td>0.86</td></tr><tr><td>ENACT – Day</td><td>0.69</td><td>0.75</td><td>0.77</td></tr><tr><td>ENACT – Night</td><td>0.70</td><td>0.77</td><td>0.79</td></tr><tr><td><strong>ENACT – Day</strong></td><td></td><td></td><td></td></tr><tr><td>Present at daytime</td><td>0.70</td><td>0.76</td><td>0.76</td></tr><tr><td><strong>ENACT – Night</strong></td><td></td><td></td><td></td></tr><tr><td>Present at nighttime</td><td>0.81</td><td>0.84</td><td>0.83</td></tr><tr><td></td><td>Urban areas 50,000–200,000 inhabitants</td><td>Urban areas ≤ 50,000 inhabitants</td><td>Municipalities outside of cities attraction</td></tr><tr><td><strong>Residents</strong></td><td></td><td></td><td></td></tr><tr><td>Present at daytime</td><td>0.76</td><td>0.69</td><td>0.58</td></tr><tr><td>Present at nighttime</td><td>0.81</td><td>0.73</td><td>0.62</td></tr><tr><td>ENACT – Day</td><td>0.75</td><td>0.71</td><td>0.69</td></tr><tr><td>ENACT – Night</td><td>0.79</td><td>0.73</td><td>0.74</td></tr><tr><td><strong>ENACT – Day</strong></td><td></td><td></td><td></td></tr><tr><td>Present at daytime</td><td>0.70</td><td>0.65</td><td>0.52</td></tr><tr><td><strong>ENACT – Night</strong></td><td></td><td></td><td></td></tr><tr><td>Present at nighttime</td><td>0.76</td><td>0.70</td><td>0.59</td></tr></tbody></table>
700,000 inhabitants) and at night. The Paris urban area stands out as where the correlation between present population at daytime and resident population is the lowest (0.42) and is the most in contrast with the correlation between present population at nighttime and resident population (0.77). The correlation between MNO and ENACT densities is particularly low for the smaller urban areas and municipalities outside of cities attraction. This reflects that our methodology probably suffers in areas with bad network coverage where the spatial mapping is imprecise.

# 5. Discussion

## 5.1. Other Works on the Present Population Estimation for Official Statistics

These experimental present population estimates were built with knowledge and inspiration from a number of existing works (Salgado et al. 2021; Statistics Netherlands 2020; Ricciato et al. 2020; Ricciato and Coluccia 2021). However, we found that off-the-shelf solutions were never fully applicable to our case, and resorted to a number of simplifications which we discuss here.

Salgado et al. (2021) propose an ambitious bayesian general framework for producing statistics from mobile network data, based on a Hidden Markov Model (HMM) modelization for device-level trajectories. In contrast, we do not rely in this work on any inference framework to quantify the uncertainties in our final estimates. Of course, this is a downside, but the level of computational complexity entailed by resorting to this modelization seemed prohibitive in our context. Static spatial mapping on cell-level aggregates, as opposed to the dynamical spatial mapping at the device level delivered by the HMM model, was chosen to avoid the computational burden. One significant advantage of our method is that it gives the ability to vary the spatial mapping after the computationally intensive aggregation step, as various spatial mappings have been shown to provide quite different results. On a similar signaling data set restricted to a large urban area, Bonnetain et al. (2019) resort to a hidden Markov chain modelization for map-matching device trajectories on the transportation network, but simplify the cell spatial coverage information to a Voronoï tesselation. The computational complexity of resorting to simple temporal interpolation has nothing to do with setting up a HMM estimation in a high-dimensional states space (up to 55 million tiles), emissions probabilities (connecting these millions of tiles to hundreds of thousands of cells), and devices (about twenty millions three-months trajectories). Although the problem is in theory parallelizable at the device level, the single device problem can be quickly high dimensional in space and time.

One of the strengths of the HMM model is to probabilistically recover trajectories when the device is unobserved from future and past observations. Given Figure 1, it is a guarantee against network and behavioral effects which seems highly desirable. We see this figure as urging for longitudinal views to derive sensible statistics. We therefore resorted to a simple interpolation method. Only a few works mention interpolation as a key feature for deriving reliable present population statistics, whereas we tend to consider interpolation as essential for sustainable and comparable-in-time statistics. Ricciato et al. (2020) point it has a promising line for future research. Interpolation techniques were explored mostly for mobility analysis. The issues of data time sparsity and sensitivity to user behavior are generic, and we show that they apply as well when working with
signaling data on a present population estimation use case. Up to now, existing literature derived snapshots of “dynamic population” (e.g., within a given day) due to the lack of access to longitudinal data for research. Throughout the lockdown period imposed in France during the Covid-19 crisis, some methods used by MNO showed sensitivity to changes in behavior (increased usage of the phone, change in the timing of usage). Roughly, the increased presence over the network translated into more detected devices and therefore to unexpected large mechanical increases in present population estimates. The French National Statistical Institute computed estimates based on reprocessed multi-MNO data (Galiana et al. 2020b; Coudin et al. 2021). In addition, interpolation has a positive effect on representativeness if the extent of network detection is correlated with socio-economic background: less active users (or users having access to a less performant mobile phone or local network) contribute to aggregates more equally after interpolation relatively to more active users. Network usage has indeed been shown related to individual characteristics, such as gender (Jahani et al. 2017).

CBS (Statistics Netherlands) has recently published a report on its methodology (Statistics Netherlands 2020) which shares similarities with our implementation. Working on signaling data of one MNO in the Netherlands, they integrate a device presence estimation with a residency detection module. They perform a home-cell detection step, where the computation barrier appears to be important as well – and based on a similar heuristic. Their calibration step is based on rescaling the estimated number of active devices to the number of local residents independently of their places of presence. Implicitly, the minimal consistency constraints are therefore the same as the ones we impose. One difference is that Statistics Netherlands (2020) does not interpolate device-level trajectories. Only the calibration step ensures the consistency constraints.

Ricciato et al. (2020) stress the practical importance of the geolocation step. In this work as here, the geolocation step is performed after cell-level aggregation. Ricciato and Coluccia (2021) propose several classes of estimators based on matrix A and cell-level aggregates which could be considered in place of our bayesian spatial mapping, for instance, to deliver confidence intervals and improve spatial accuracy. In this article, the adopted linear estimator corresponds to the "Simple Bayes-rule estimator" described in Ricciato and Coluccia (2021) applied to reweighted device counts. Our framework, which introduces device-level weights and a home estimation, could easily be integrated with the other estimators introduced in Ricciato and Coluccia (2021) by replacing the simple count vector c with a reweighed version.

## 5.2. Discussion on Device-Level Weights

We argue that residency detection is not only useful for statistical filtering, but also to balance the estimates across residencies to correct for unbalanced representativeness (as illustrated in Figures 6 and 10). Imposing a constraint of equality at a local level between “usual” residents detected from mobile-phone data and actual residents as estimated from official sources appears milder than an alternative that would consist of equalizing resident population $P_r$ to population present at night: $u_{r,t_0} = P_r$. A significant part of the population spends nights regularly outside of its main residency or works at night. These atypical location behaviors are better captured by MNO data than by traditional sources. Figure 14
shows how the present population at night can vary greatly during the week – and it is also even more the case for holidays and bank holidays. Note that if the pseudo-weights are based only on residency location, they could be based on as many characteristics as we can accurately recover at the device level and for which we have an external population-level estimate (see the discussion on pseudo-weights in Beresewicz et al. (2018)). This framework could be promising for future combinations of sources aiming at facing selectivity issues to derive representative statistics. To support this statement, we report in Table 6 how the consistency between our estimate and ENACT estimate drops when we do not use individual residency-based weights but simply multiply the device density by an hourly constant to match the total population size. This consistency decreases noticeably both in terms of allocation accuracy and correlation, at night and daytime. Rank correlation comparisons are similar, suggesting that rank ordering across tiles is rather close when using a constant weighting scheme.

Using device-level weights based on individual characteristics may lift part of the concern of using a single MNO data set. If mobility patterns of the MNO customers are similar to that of the resident populations which share these same characteristics, the estimation would be unbiased. Of course, the larger and the more representative the mobile phone users used in the analysis, the better will be the quality of the estimates. Therefore, multiple MNOs data sets could be fruitfully employed with a similar methodology. Device-level weights are obtained from merging MNO data and reliable and localized estimates of the resident populations. The latter is more and more commonly produced for a large number of countries. For instance in Europe, population density have been disseminated on a common European 1 km² reference grid in 2011 as a prototype (GEOSTAT-1) and a legal act was adopted for Census 2021 for a harmonised publication of key census topics on an EU-wide 1 km² grid by the European Statistical System.

### 5.3. Quantifying Uncertainty

This work does not provide an end-to-end variance estimator for our final estimates but analyzes uncertainties at various stages while building the estimates, which will

<table><caption>Table 6. Comparisons between the ENACT estimate and both our present population estimate as well as globally reweighted device densities. Note: All densities are computed in our dissemination grid, using proportional area estimation for ENACT data and direct calculations for resident density from tax files. ENACT: March 2011. Present population and rescaled device densities: Five working days in March 2019. Resident Population: 2016. Rescaled devices estimates are obtained when omitting reweighting by residency and rescaling device densities hour by hour with a single scalar weight. Daytime (resp. nighttime) estimates are taken at 3 p.m. (resp. 3 a.m.).</caption><thead><tr><th></th><th>Allocation accuracy</th><th>Correlation</th><th>Rank correlation</th></tr></thead><tbody><tr><td><strong>ENACT – Day</strong></td><td></td><td></td><td></td></tr><tr><td>Present at daytime</td><td>0.74</td><td>0.75</td><td>0.89</td></tr><tr><td>Rescaled devices at daytime</td><td>0.68</td><td>0.66</td><td>0.89</td></tr><tr><td><strong>ENACT – Night</strong></td><td></td><td></td><td></td></tr><tr><td>Present at nighttime</td><td>0.79</td><td>0.87</td><td>0.88</td></tr><tr><td>Rescaled devices at nighttime</td><td>0.73</td><td>0.78</td><td>0.88</td></tr></tbody></table>
hopefully help future theoretical development. We estimate the spatial mapping accuracy locally through the probability to localize in *i* a device that is in *i*₀, *N*ᵢᵢ₀, where *N* = *QA*. With previous notations, if we take the example of localizing a single device *d* which is in *i*₀, that is **u**ᵢ = 1ᵢ₀, N1ᵢ₀ can be interpreted as E[|**u**ᵢ|device in tile i₀]. Thus, ||N1ᵢ₀ - 1ᵢ₀|| provides a local evaluation of the spatial accuracy: it is the root mean square error of the spatialisation task of a device which is in i₀. We provide in Figure 17 maps of this spatial uncertainty based on this framework, with *Q* a spatial mapping base on a uniform or a resident density prior. Dense areas have a small spatial error compared to rural areas, as could be expected from the network density. In addition, using a prior helps in reducing this spatial uncertainty. As for the temporal uncertainties, with the large time event rate, we find that interpolation errors are limited for hourly estimates of present population (Figure 12). Finally, part of population coverage uncertainties may be captured through the weights wᵣ = Pᵣ/Dᵣ, as this ratio informs on differential local representativity. Although informative, a map of wᵣ also divulgates Orange market shares and is thus not reported.

<figure><img src="image_18.png" /><figcaption>Fig. 17. RMSE of the spatial mapping, with a linear spatial mapping with an uniform or resident density prior.</figcaption></figure>
### 5.4. Other Limitations

In addition to the limitations discussed above, some issues were left aside in this work.

First, we hypothesize that a person carries a single device and we did not attempt to filter out potential additional devices. The issue is probably less salient than when relying on multiple MNO data, but even in a single MNO framework, identifying from the data devices owned by the same person can be a complex task. In our work, without additional filtering, the hourly population is biased toward the mobility of persons carrying multiple devices with an Orange contract which seemingly represent multiple persons.

Second, our methodology should be adapted next to border areas – as we only observe the Orange network on the French territory. We here detail the consequences of daily worker flows, but note that similar issue could be described for touristic outflows and location where devices are last seen over the French territory (such as airport). Take the example of daily outflows of cross-border workers leaving the French territory in the morning and coming back at night. Due to the interpolation method (device presence panel step), the trajectories of devices carried by these workers will stop and remain at the border in the meantime. Thus, present population will be overestimated in tiles next to the border during the day. As for cross-border inflows, the presence of foreign residents carrying an Orange device with a French Orange contract will be interpolated to cells next to the border at night. In addition, their attributed residence will most likely be their workplace cells. This complex subject should be examined in detail in future work.

## 6. Conclusion

We derived a prototype for an experimental population statistic for France. Daily and weekly cycles, both at the local and national levels offer unprecedented insight into population dynamics in France. However, to develop a full-fledged methodology, access is of utmost importance. A legal basis for processing MNO data for official statistics under due privacy protection, as well as the cooperation of MNOs, is of primary importance. It is all the more challenging today that this work tends to demonstrate that some form of access to longitudinal individual data would be needed for deriving reliable population estimates, even aside from an interest in mobility analysis (interpolation to avoid being plagued with activity bias, home detection or device-level characterization to ensure representativeness, deduplication...). Data management from the MNO side seems decisive for the quality of the final statistics (network topology modelization and cell register management, filtering of IoT and M2M, reports on data collection failure...). A combination of sources from the MNO and the NSI sides at a fine-grained level (e.g., build weights for representativeness) seems very promising. Data fusion approaches, merging high-quality population and mobile network data sets could be considered and studied in the future to make the most of both sources. Nonetheless, the combination of sources requires cooperation, privacy-preserving transfer of information, and a transparent sharing of computations – which has been at this stage possible only within research projects, if at all, in European countries.
7. Appendix

<figure><img src="image_19.png" /><figcaption>Fig. 18. Distinct dates of observation over the period.</figcaption></figure>

<table><caption>Table 7. Dissemination grid structure and average resident density. In our dissemination grid, presented in the section "Embedding accuracy within dissemination", and illustrated in Figure 11, there are 4,062 tiles of 200 meters square. The average resident density in these tiles is 15,810 persons per square kilometer (source: Filosofi).</caption><thead><tr><td>Tile side (meters)</td><td>Resident density (persons/km²)</td><td>Number of tiles</td></tr></thead><tbody><tr><td>100</td><td>22970</td><td>2816</td></tr><tr><td>200</td><td>15810</td><td>4062</td></tr><tr><td>400</td><td>5720</td><td>8796</td></tr><tr><td>800</td><td>1670</td><td>12721</td></tr><tr><td>1600</td><td>280</td><td>20743</td></tr><tr><td>3200</td><td>60</td><td>28339</td></tr><tr><td>6400</td><td>30</td><td>5024</td></tr><tr><td>12800</td><td>10</td><td>2</td></tr></tbody></table>
<figure><img src="image_20.png" /><figcaption>Fig. 19. Paris area context, with Département administrative boundaries.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 20. Differences between present population at 3 a.m. and resident population (left) and between present population at 3 p.m. and ENACT day time population (right). The first panel shows urban attraction areas with more than 700,000 inhabitants.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. 21. Differences between present population at 3 a.m. and resident population (left) and between present population at 3 p.m. and ENACT day time population (right).</figcaption></figure>

<table><caption>Table 8. Notations used throughout the document.</caption><tr><td>d &isin; D</td><td>Devices in the scope, detected on the network over the period of interest</td></tr><tr><td>t &isin; T</td><td>Hourly time grid of interest (several weeks)</td></tr><tr><td>P</td><td>The target population set</td></tr><tr><td>D<sub>t</sub></td><td>Set of devices in the scope, detected on the network during t</td></tr><tr><td>D<sub>l</sub></td><td>Set of devices in the scope, detected on the network at date l</td></tr><tr><td>i &isin; I</td><td>Tiles that cover the territory of interest</td></tr><tr><td>j &isin; J</td><td>Cells of the MNO network</td></tr><tr><td>j<sub>d,t</sub> &isin; J</td><td>Estimated presence cell for device d during time interval t</td></tr><tr><td>P<sub>r</sub></td><td>Official residents in place r</td></tr><tr><td>D<sub>r</sub></td><td>Devices with estimated residency in r</td></tr><tr><td>u<sub>t</sub> in R<sup>l</sup></td><td>Population counts at time t over tiles (target statistics)</td></tr><tr><td>m<sub>r,t</sub> &isin; N<sup>l</sup></td><td>Devices count in presence cells at time t with residency r</td></tr><tr><td>u<sub>0</sub> &isin; N<sup>l</sup></td><td>Population count from official source over tiles</td></tr><tr><td>u&#770;<sub>r</sub> &isin; R<sup>l</sup></td><td>Estimated population count at time t over tiles</td></tr><tr><td>u&#770;<sub>r,t</sub> &isin; R<sup>l</sup></td><td>Estimated population count at time t who are resident in r over titles</td></tr><tr><td>w<sub>r</sub></td><td>Pseudo-weight of devices residing in place r</td></tr><tr><td>w<sub>r,l</sub></td><td>Final weight of devices in place r on date l</td></tr><tr><td>A</td><td>Coverage probability matrix</td></tr><tr><td>Q: R<sup>l</sup> &isin; R<sup>l</sup></td><td>Linear spatial mapping from cells to tiles</td></tr><tr><td>N<sub>i,t<sub>0</sub></sub></td><td>Probability that a device in tile i<sub>0</sub> is mapped in tile i</td></tr></table>

## 8. References

Avouac, R., B. Sakarovitch, and Z. Smoreda. 2019. “A Bayesian Approach to Improve the Estimation of Population using Mobile Phone Data.” Research Workshop on Digital Demography in the Era of Big Data, Seville, Spain, 6–7 June 2019. Available at: https://iussp.org/en/digital - demography - era - big - data (accessed March 2023).

Batista e Silva, F., S. Freire, M. Schiavina, K. Rosina, M.A. Marìn-Herrera, L. Ziemba, M. Craglia, E. Koomen, and C. Lavalle. 2020. “Uncovering Temporal Changes in Europe’s Population Density Patterns using a Data Fusion Approach.” *Nature communications* 11(1): 1–11. DOI: https://doi.org/10.1038/s41467-020-18344-5.
Beresewicz, M., R. Lehtonen, F. Reis, L. Di Consiglio, and M. Karlberg. 2018. *An Overview of Methods for Treating Selectivity in Big Data Sources*. Technical report. Eurostat Statistical Working Article. Available at: https://doi.org/10.2785/312232.

Blondel, V.D, A. Decuyper, and G. Krings. 2015. "A Survey of Results on Mobile Phone Data Sets Analysis." *EPJ data science* 4(1): 10. DOI: https://doi.org/10.1140/epjds/s13688-015-0046-0.

Bonnetain, L., A. Furno, N.-E. el Faouzi, M. Fiore, R. Stanica, Z. Smoreda, and C. Ziemlicki. 2021. "TRANSIT: Fine-Grained Human Mobility Trajectory Inference at Scale with Mobile Network Signaling Data." *Transportation Research Part C: Emerging Technologies* 130: 103257. DOI: https://doi.org/10.1016/j.trc.2021.103257.

Bonnetain, L., A. Furno, J. Krug, and N.-E. el Faouzi. 2019. "Can We Map-Match Individual Cellular Network Signaling Trajectories in Urban Environments? Data-Driven Study." *Transportation Research Record* 2673(7): 74–88. DOI: https://doi.org/10.1177/0361198119847472.

Chen, G., A. Carneiro Viana, M. Fiore, and C. Sarraute. 2019. "Complete Trajectory Reconstruction from Sparse Mobile Phone Data." *EPJ Data Science* 8(1): 30. DOI: https://doi.org/10.1140/epjds/s13688-019-0206-8.

Coudin, E., M. Poulhes, and M. Suarez Castillo. 2021. "The French Official Statistics Strategy: Combining Signaling Data from Various Mobile Network Operators for Documenting COVID-19 Crisis Effects on Population Movements and Economic Outlook." *Data & Policy* 3. DOI: https://doi.org/10.1017/dap.2021.1.

ESSnet Big Data, II. 2021. *Work Package I – Mobile Networks Data Deliverables*. Available at: https://ec.europa.eu/eurostat/cros/content/wpi-milestones-and-deliverables_en. (accessed March 2023).

Fekih, M., T. Bellemans, Z. Smoreda, P. Bonnel, A. Furno, and S. Galland. 2021. "A Data-Driven Approach for Origin – Destination Matrix Construction from Cellular Network Signalling Data: a Case Study of Lyon Region (France)." *Transportation* 48(4): 1671–1702. DOI: https://doi.org/10.1007/s11116-020-10108-w.

Freire, S., and C. Aubrecht. 2012. "Integrating population dynamics into mapping human exposure to seismic hazard." *Natural Hazards and Earth System Sciences* 12(11): 3533–3543. DOI: https://doi.org/10.5194/nhess-12-3533-2012.

Galiana, L., B. Sakarovitch, F. Sémécurbe, and Z. Smoreda. 2020. *Residential Segregation, Daytime Segregation and Spatial Frictions: an Analysis from Mobile Phone Data, Insee Working article G2020-12*. Available at: https://www.insee.fr/en/-statistiques/4925202.

Galiana, L., M. Suarez Castillo, F. Sémécurbe, E. Coudin, and M.-P. de Bellefon. 2020b. *Retour Partiel des Mouvements de Population avec le Déconfinement, Insee Analyses* 54. Available at: https://www.insee.fr/fr/statistiques/4635407.

Hoteit, S., G. Chen, A. Viana, and M. Fiore. 2016. "Filling the Gaps: On the Completion of Sparse Call Detail Records for Mobility Analysis." In *Proceedings of the eleventh ACM workshop on challenged networks: 45–50*. DOI: https://doi.org/10.1145/2979683.2979685.

Insee. 2020. *Population présente sur le territoire avant et après le début du confinement – Premiers résultats (communiqué de presse)*. Available at: https://www.insee.fr/fr/information/4477356 (accessed March 2023).
Jahani, E., P. Sundsoy, J. Bjelland, L. Bengtsson, A. Pentland, and Y.-A. de Montjoye. 2017. "Improving Official Statistics in Emerging Markets using Machine Learning and Mobile Phone Data." *EPJ Data Science* 6(1): 3. DOI: https://doi.org/10.1140/epjds/s13688-017-0099-3.

Katsikouli, P., M. Fiore, A. Furno, and R. Stanica. 2019. "Characterizing and Removing Oscillations in Mobile Phone Location Data." In 2019 IEEE 20th International Symposium on "A World of Wireless, Mobile and Multimedia Networks" (WoWMoM): 1-9. IEEE. DOI: https://doi.org/10.1109/WoWMoM.2019.8793034.

Lamarche, P., and S. Lollivier. 2021. "Fidéli, l'intégration des sources fiscales dans les données sociales *Courrier des statistiques* 6. Available at: https://www.insee.fr/fr/information/5398683 (accessed March 2023).

Léetroublon, C., and C. Daniel. 2018. *Le Travail en Horaires Atypiques: Quels Salariés pour quelle Organisation du Temps de Travail? Dares Analyses 2018-30*. Available at: https://dares.travail-emploi.gouv.fr/publications/le-travail-en-horaires-atypiques. (accessed March 2023).

Panczak, R., E. Charles-Edwards, and J. Corcoran. 2020. "Estimating Temporary Populations: a Systematic Review of the Empirical Literature." *Humanities and Social Sciences Communications* 6(1): 1-10. DOI: https://doi.org/10.1057/s41599-020-0455-y.

Ricciato, F., and A. 2021. "On the Estimation of Spatial Density from Mobile Network Operator Data." *IEEE Transactions on Mobile Computing*, 22(6): 1. DOI: https://doi.org/10.1109/TMC.2021.3134561.

Ricciato, F., G. Lanzieri, A. Wirthmann, and G. Seynaeve. 2020. "Towards a Methodological Framework For Estimating Present Population Density from Mobile Network Operator Data." *Pervasive and Mobile Computing* 101263. DOI: https://doi.org/10.1016/j.pmcj.2020.101263.

Ricciato, F., P. Widhalm, F. Pantisano, and M. Craglia. 2017. "Beyond the "Single-Operator, CDR-only Paradigm: An Interoperable Framework for Mobile Phone Network Data Analyses and Population Density Estimation." *Pervasive and Mobile Computing* 35: 65-82. DOI: https://doi.org/10.1016/j.pmcj.2016.04.009.

Ricroch, L., and B. Roumier. 2011. *Depuis 11 Ans, Moins de Taches Ménagères, Plus d'Internet. Insee premièreure 1377*. Available at: https://www.insee.fr/fr/statistiques/1281050. (accessed March 2023).

Sakarovitch, B., M.-P. de Bellefon, P. Givord, and M. Vanhoof. 2018. "Estimating the Residential Population from Mobile Phone Data, an Initial Exploration." *Economie et Statistique* 505(1): 109-132. DOI: https://doi.org/10.24187/ecostat.2018.505d.1968.

Salat, H., Z. Smoreda, and M. Schläpfer. 2020. "A Method to Estimate Population Densities and Electricity Consumption from Mobile Phone Data in Developing Countries." *PLOS ONE* 15(6): 1-11. DOI: https://doi.org/10.1371/journal.pone.0235224.

Salgado, D., L. Sanguiao, B. Oancea, S. Barragàn, and M. Necula. 2021. "An End-to-End Statistical Process with Mobile Network Data for Official Statistics." *EPJ Data Science* 10(1): 1-46. DOI: https://doi.org/10.1140/epjds/s13688-021-00275-w.

Schiavina, M., S. Freire, K. Rosina, L. Ziemba, M. Marin Herrera, M. Craglia, C. Lavalle, T. Kemper, and F. Batista. 2020. *ENACT-POP R2020A-ENACT 2011 Population Grid*.
European Commission, Joint Research Centre (JRC). DOI: https://doi.org/10.2905/-BE02937C-5A08-4732-A24A-03E0A48BDCDA.

Statistics Netherlands, C.B.S. 2020. *Estimating Hourly Population Flows in the Netherlands*. Avalialble at: https://www.cbs.nl/en-gb/longread/diversen/2020/pilot-study-mobile-phone-meta-data-records-introduction-to-the-research-method (accessed March 2023).

Tennekes, M., and Y. Gootzen. 2022. “Bayesian location estimation of mobile devices using a signal strength model.” *Journal of Spatial Information Science* 25: 29–66. DOI: https://doi.org/10.5311/JOSIS.2022.25.166.

U.N. Global Working Group. 2019. *Handbook on the Use of Mobile Phone Data for Official Statistics*. Technical report. United Nations. Available at: https://unstats.un.org/bigdata/task-teams/mobile-phone/MPD%20Handbook%2020191004.pdf.

Vanhoof, M., F. Reis, T. Ploetz, and Z. Smoreda. 2018. “Assessing the Quality of Home Detection from Mobile Phone Data for Official Statistics.” *Journal of Official Statistics* 34(4): 935–960. DOI: https://doi.org/10.2478/jos-2018-0046.

Received March 2022

Revised August 2022

Accepted March 2023