Measuring Commuting and Economic Activity Inside Pt 1

Measuring Commuting and Economic Activity inside Cities with Cell Phone Records Abstract We show how to use commuting flows to infer the spatial distribution of income within a city. A simple workplace choice model predicts a gravity equation for commuting flows whose destination fixed effects correspond to wages. We implement this method with cell phone transaction data from Dhaka and Colombo. Model-predicted income predicts separate income data, at the workplace and residential level, and by skill group. Unlike machine learning approaches, our method does not require training data, yet achieves comparable predictive power. We show that hartals (transportation strikes) in Dhaka reduce commuting more for high model-predicted wage and high skill commuters. Measures of urban economic activity at fine temporal and spatial scales are important yet rare. Such data is necessary to understand how cities respond to localized shocks such as changes in transportation infrastructure or floods, and to help governments target scarce public resources. These issues are especially salient in large cities in developing countries, which are growing fast yet are least covered by conventional data sources. At the same time, comprehensive new data sources on urban behavior, especially individual mobility and commuting, are becoming available worldwide. Indeed, academic research publications using call detail records (CDR) or cell phone data covered at least 19 out of 62 countries in Africa and Asia as of 2020. Meanwhile, less than 10% of the urban population in sub-saharan African countries is covered by a census of firms with wage data.1 In this paper, we provide a theory-based method to predict the spatial distribution of urban economic activity from commuting choices. The revealed-preference logic of our approach is simple. A core function of cities is to connect workers and jobs (Duranton and Puga 2015). While many factors enter into workplace choice decisions, areas with high wages should disproportionately attract workers, keeping distance and home locations fixed. We propose inverting this reasoning to infer the relative wage at a location based on its “attractiveness” as a commuting destination. We formalize this intuition by building on recent urban economics models of commuting choices. In these models, work location decisions aggregate up to a gravity equation on commuting flows, and destination fixed effects are proportional to log wages. This property holds for a general class of models developed to evaluate urban policies and transport infrastructure (Redding and Turner 2015, Redding and Rossi-Hansberg 2017). We also show how to apply a similar method when wages and the commuting elasticity differ flexibly by skill group. We implement our approach using call detail record (CDR) data from two large metropolises: Colombo, Sri Lanka and Dhaka, Bangladesh. CDR data is a prototypical example of “big data” available in developing countries (Björkegren 2018), and it contains phone user location for every transaction. We construct individual home and work locations by observing a user’s 1 Authors’ calculations (Appendix A). location at different times of the day over time. Commuting flows constructed this way have very fine geographic resolution and correlate strongly with commuting flows from a transportation survey from Dhaka. We use this data to estimate the model gravity equation. Next, we assess how well this simple measure captures real differences in income. First, we show that model workplace income is a robust predictor of workplace commuter income data from a large transportation survey in Dhaka. Second, in both cities, model-predicted residential income is a robust predictor of a census residential income proxy. To set the right benchmark for the predictive performance, we also implement a supervised learning approach (elastic net regularization) using hundreds of features from cell-phone data and geographic measures (Zou and Hastie 2005, Blumenstock et al. 2015). The model prediction, despite being computed without any training data, explains between 70% and 90% as much income variation compared to the supervised learning approach. Hence, the destination fixed effects act as a “summary statistic” for most information in cell phone data. We also estimate and validate the model extension with multiple worker skills. The ideal application of our income-prediction method and of the high-frequency commuting data is to trace out heterogeneous impacts of urban events and policies, such as transportation shocks, lockdowns or floods. We study hartals, a type of strike intended to disrupt transportation and economic activity in Bangladesh. We find larger reductions in commuting on hartal days for commuters with higher model-predicted income, high skill, and long commute. We build on a growing literature using quantitative urban models. While papers in this literature often use gravity equations to estimate structural parameters as part of a broader exercise (Ahlfeldt et al. 2015, Monte et al. 2018, Owens et al. 2020, Tsivanidis 2019, Severen 2019, Heblich et al. 2018, Dingel and Tintelnot 2020), our focus is to use gravity equations to construct proxies of the spatial distribution of wages. Another related literature uses machine learning to empirically predict wealth and consumption at individual or regional level (Blumenstock et al. 2015, Jean et al. 2016, Glaeser et al. 2017). A key feature of our approach is that it does not require training data, relying instead on a simple and general theory of commuting behavior. We believe that one promising path for “big data” in economics is using revealed preference techniques to infer attributes of choice options, such as workplace wages in our paper or spatial aspects of consumption behavior (Athey et al. 2018, Davis et al. 2018, Agarwal et al. 2018). 1 Cell-Phone Data and Commuting Flows Cell phone transaction data. We use call detail record (CDR) data from large operators in Sri Lanka and Bangladesh to compute detailed commuting matrices. CDR data includes an observation for each transaction, such as outgoing or incoming voice call and text messages, or GPRS internet connections. Each observation has a timestamp, the anonymized participant user identifiers, and their cell tower locations. Towers are unevenly distributed in space; they are denser in urban and developed areas. We focus on the greater metropolitan areas around the capital cities of Colombo and Dhaka. The data covers a little over a year in Sri Lanka and four months in Bangladesh in the early 2010’s.2 We construct commuting trips by assigning “home” and “work” locations for each user. Home (work) locations are identified as the most frequent towers with a transaction between 9pm to 5am of the next day (10am to 3pm) during weekdays excluding hartal days. For robustness, we also construct daily commuting trips.3 We then aggregate over users to obtain an origin-destination (OD) matrix of commuting flows between every pair of cell towers. Google Maps travel time. As a proxy for travel costs, we obtain estimated typical driving travel times between pairs of cell towers using the Google Maps API. In each city we obtain Google data for 90,000 randomly selected pairs of towers, and interpolate to pairs with nearby origin and nearby destination. We use predicted time without traffic congestion. Using predicted time with traffic congestion in Colombo, where such data was available, yields virtually identical model-predicted wages (Table G.4). 2 In Bangladesh, the data only covers outgoing voice calls. Our sample covers the Western Province in Sri Lanka, and the Dhaka, Narayanganj, and Gazipur Districts in Bangladesh. 3 To construct daily commuting trips, on a given day, we define a user’s origin as the location of the first transaction between 5am to 10am, and the user’s destination as the location of the last transaction between 10am and 3pm. If transaction data is missing in either time interval, commuting behavior is not observed for that user-day (Table G.1). Household transportation survey. We use individual survey data from the 2009 Dhaka Urban Transport Network Development Study or DHUTS (JICA 2010). The survey covers 16,394 randomly selected households in the Dhaka City Corporation (DCC), Dhaka’s urban core. Home and work locations are at the level of 90 “survey areas.” Our main analysis sample covers 12,510 commuters who live and work within the DCC, with positive income from work, excluding students, homemakers, and the unemployed. Population Censuses. We use census data from 2011 in Bangladesh and 2012 in Sri Lanka, the closest years to our cell phone data. Since the census does not report income in either country, we obtain the first principal component of houshold assets (house building materials, toilet facilities, water and electricity connection) at the finest geographic unit available.4 The residential income proxy at the cell tower level is the average across overlapping census units, weighted by overlap area with the tower. Representativeness of Commuters in Cell Phone Data. In Dhaka, commuting flows derived from cell phone data are strongly related to those from the DHUTS commuting survey, including when controlling for log travel time, origin and destination survey area fixed effects (Table G.2). This is consistent with previous research validating cell-phone-based commuting flows (Calabrese et al. 2011, Wang et al. 2012, Iqbal et al. 2014). The decay of commuting flows with travel time is virtually identical between the two data sources (Figure G.2, Panel A). Residential population density from cell phone data is correlated with census population density at the level of 1,866 and 1,201 cell phone towers in the two cities, with R2 = 0.61 in Dhaka and R2 = 0.49 in Colombo (Table G.3). The slope is 1.16 for both cities, hence cell phone data slightly over-represents population in denser areas. This type of bias does not automatically affect our results as our approach uses workplace choice shares conditional on a residential location. In section 2.4 we further allow choices to differ by worker skill. 4 In the study areas, there are 2,381 Grama Niladhari (GN) in Sri Lanka, and 3,704 mauza in Dhaka. 2 Model: Commuting Flows, Gravity, and Wages Is it possible to infer the spatial distribution of wages from commuting flows? The interaction between wages and commuting costs to determine urban structure is fundamental in classical urban economics models (Alonso 1960, Mills 1967, Muth 1968). Here, we explore this insight using a new generation of models inspired from the trade literature, designed to better match spatially disaggregated urban data (Ahlfeldt et al. 2015). 2.1 Workplace Choice Model Space is partitioned into a finite set of locations, which may serve as both residential and work locations. We discuss how to map model locations to the data in section 2.3. There is a unit mass of workers, and each worker ω sequentially decides her residential location (or origin) i, and then her work location (or destination) j. Conditional on residing in location i, the utility of worker ω if she chooses destination j is: Uijω = Wj Zijω Dijτ (1) Wj is the wage per effective unit of labor supply at location j (all firms at location j offer the same wage), Dij is the travel time between i and j, and Zijω is an idiosyncratic utility shock that is i.i.d. following the Fréchet distribution, with scale parameter normazlied to one and shape parameter e.5 Each worker supplies one unit of labor and earns income Wj . Each worker observes the shocks Zijω and chooses the work location j where Uijω is maximized. The probability that a worker commutes to j conditional on residing in i is given by6 πij|i = (Wj /Dijτ )e ∑s (Ws /Disτ )e (2) follow Ahlfeldt et al. (2014) and assume that Zijω are preference shocks. Tsivanidis (2019) alternatively assumes that Zijω are productivity shocks and derives expected income as “commuter market access”. Appendix C shows that our results are robust to alternate assumptions, and estimates a model where Zijω and Dij partly affect productivity and utility. 6 Assuming joint home and work location choice leads to the same conditional choice probabilities π ij|i as residential terms (amenities, rent) cancel out (Ahlfeldt et al. 2015). However, if workers choose their workplace first and then the home location (e.g. migrants), destination fixed effects would not exclusively capture destination wages. 5 We Taking logs, and denoting log quantities by lowercase letters: log(πij|i ) = ew j − eτdij − log ∑ exp (ews − eτdis ) ! (3) s 2.2 Estimating the Gravity Equation We estimate equation (3) using the empirical Poisson pseudo-maximum likelihood (PPML) method with two-way fixed effects: log( E[πij|i ]) = ψj − βdij + µi (4) where µi and ψj are origin and destination fixed effects. We use PPML, rather than OLS, to deal with zero commuting flows (Silva and Tenreyro 2006, Dingel and Tintelnot 2020).7 This equation allows us to measure the attractiveness of each destination j after accounting for the fact that conditional commuting probabilities also depend on log commute duration. Importantly, in our model, ψj is proportional to the (relative) log wage at j with a factor of e, the Fréchet dispersion parameter. Our main goal is to recover the ψj ’s from observed commuting choices. For this purpose, it is not necessary to model explicitly how wages are determined in equilibrium. The mapping between commuting choices and wages holds in any general equilibrium model that micro-founds the gravity equation for commuting flows with a discrete commuting choice model.8 To obtain relative wage levels, we further need to know e, the Fréchet parameter, which governs the variance of idiosyncratic preference shocks. In section 3, we will estimate e in Dhaka using data on average wages by location in the city. When such data is not available, e may also be estimated from the overall variance of wages in the entire city, as in Ahlfeldt et al. (2015). 7 Log travel time as a measure of commuting costs offers a good fit (Figure G.2). model does not include workplace amenities. If these differ considerably across space, the gravity destination fixed effects will capture the combined effect of wages and amenities. Our empirical results in section 3 address empirically the extent to which our measure is correlated with wages. 8 Our 2.3 Mapping Model Locations to Geographic Areas A key advantage of the model is that locations can be mapped directly to two-dimensional urban data. However, exactly how we map model locations to locations in the data may not be innocuous for inferring wages. The difficulty of defining appropriate geographic units is related to the “modifiable areal unit problem” in geography (Fotheringham and Wong 1991), and to the problem of defining urban areas in economics (Rozenfeld et al. 2011, Baragwanath et al. 2019). We show how to solve this issue in a theory-consistent way. Assume that each Voronoi cell j consists of Nj sub-locations and workers draw independent Fréchet shocks at the levels of these sub-locations. Hence, the model implies a gravity equation at the level of sub-locations. However, a gravity equation continues to hold at the level of the larger Voronoi cells under two simplifying assumptions: that all sub-locations within j offer the same log wage w Rj , and that all sub-location in j are located in the same place in space (and hence they are equally distant to any given origin i). Using standard Fréchet properties, the underlying log wage w Rj is expressed as:9
ψjR = ew Rj = ψj − log Nj , (5) where ψj is the destination fixed effect from the gravity equation between Voronoi cells. In other words, we need to adjust the destination fixed effects by subtracting a term (log( Nj )) that is increasing in the number of sub-locations in each cell. Implementing this adjustment requires taking a stand on Nj , i.e. the number of independent Fréchet shocks drawn per Voronoi cell j. Theory does not offer strong guidance on how to select Nj . Here, we assume that Nj is proportional to the geographic area of Voronoi cell j, that is, Fréchet shocks are drawn for each vertex on a square lattice. We also show robustness to assuming Nj = 1 and to aggregating our raw data at the level of square grid cells (Tables G.6, G.9, and Figure G.5). 9 When wages differ across sub-locations within j, ψjR corresponds to a CES aggregate of these wages. 2.4 Estimating Wages Separately by Skill Group Skill inequality is salient in cities like Dhaka and Colombo, and our baseline model does not capture this. We now extend our method to estimate wages separately by skill group. We consider two types of models. In our first model, workers are either low-skill (L) or high-skill (H) and face separate log wage profiles w jL and w jH and elasticities of commuting cost with respect to travel time τ L and τ H . All agents have a common Fréchet shape parameter e. As before, our goal is to estimate composite wage terms ψsj = ewsj for all j up to a constant and distance terms βs = eτ s for each skill s ∈ { L, H }. Equation (2) expresses commute probabilities πijs |i separately by skill s. The average commuting probability (that we observe in our cell phone data) is πij|i = (1 − λiH )πijL|i + λiH πijH|i where λiH is the share of high skill at residential location i (measured from census data). Because πij|i is not log-linear in the parameters of interest, we estimate the model using maximum likelihood (Appendix B). In our second model, each residential location i has a “representative agent” with preferences given by a weighted mean, with weights that depend on λiH , the share of high skill at i. Agents from i earn log wages w jH λiH + w jL (1 − λiH ) at j and have commuting elasticity τ H λiH + τ L (1 − λiH ). Plugging into equation (4), our estimating equation becomes: log(E[πij|i ]) = ψjL (1 − λiH ) + ψjH λiH − β L (1 − λiH )dij − β H λiH dij + µi (6) This equation is intuitive: each destination has two levels of “attractiveness,” ψjL = ew jL and ψjH = ew jH . Loosely speaking, ψjH are identified from commuting probability patterns from origins with large shares of high-skilled in the census. We estimate (6) using PPML with origin and destination fixed effects, log duration interacted with skill shares, and destination-specific linear slopes in λiH . The two models are identical when residences are completely segregated by skills (λiH = 0 or 1). In the intermediate case, we show by simulation that the two data generating process are closely aligned when using the same profiles of w jL and w jH (Appendix B). Because of the simpler estimation procedure, we use the second model as our benchmark. Appendix B has results from the first model. 2.5 Estimation Results: Gravity and Wages We estimate gravity equation (4) using cell phone commuting flows and Google Maps travel times. Our goal is to recover the destination fixed effects, which in the model are proportional to workplace log wages. The estimation sample is non-holiday weekday commuting trips between pairs of towers excluding nearby and very distant towers.10 Table 1 reports the results, based on commuting flows between 1,859 locations in Dhaka (columns 1-3) and between 1,201 locations in Colombo (columns 4-6). The gravity equation is estimated with commuting flows constructed from assigned home and work locations for 1.5 and 1 million commuters in the two cities (columns 1,3,4, and 6), and using the commuting flows identified at the daily level (columns 2 and 5). Commuting probability decreases strongly with travel time. Interestingly, although the average commuting trip is 25% longer on average in Sri Lanka, the coefficients are similar (-2.44 in Column 1 and -2.19 in Column 4). This is a substantive finding, as the two cities differ in terms of economic development, population, and urban structure (mono- vs polycentric). At the same time, these elasticities are substantially lower in magnitudes than in the United States (between 4 to 8, Monte et al. (2018) and Owens et al. (2020)), potentially reflecting differences in the disutility of commuting time. Turning to our model with skill heterogeneity, we estimate equation (6) using PPML, data on aggregate bilateral commuting flows, and residential-level skill shares from the population census. High-skilled is defined as literate in Dhaka, and having secondary education in Colombo. The travel time elasticity is smaller in absolute value for the high skill group, implying that low skill workers commute closer to home. Tsivanidis (2019) finds a similar result using commuting flows data by skill. Our method only requires aggregate commuting flows and the residential skill shares. Figure 1 displays smoothed estimated wages in Dhaka and Colombo using choropleth maps. Our estimated, area adjusted measure ψjR is proportional to log wages, with factor e 10 In Dhaka, we exclude 31 days with transportation strikes (hartals). Tower pairs closer than 3 minutes are excluded as they may capture calls randomly connecting to different towers (“tower-bouncing”) rather than real commuting. Destination fixed effects estimated including nearby and same tower pairs are virtually identical (Table G.4). (the Fréchet shape parameter). Estimated wages are higher near city centers and alongside some (but not all) major road corridors. Moreover, secondary centers are visible, especially in Dhaka. The next sections will compare these results with independent income proxies. Destination fixed effects using different estimation methods are highly correlated: estimating on disjoint samples (Table D.1), using daily commuting flows instead of home and work assignment, when we use travel times with congestion in Colombo, and when we include neighboring and same tower pairs in the samples (Table G.4). Using OLS instead of PPML leads to a flatter profile of destination fixed effects due to many zero commuting flows (57% of all possible tower pairs in Bangladesh and 15% in Sri Lanka). We obtain virtually identical results when allowing travel time to enter non-parametrically (as deciles of the travel time distribution). Low- and high-skill destination fixed effects are positively correlated with the destination fixed effects from the benchmark model, with respective correlation coefficients 0.71 and 0.80 in Dhaka and 0.43 and 0.67 in Colombo. However, the model also points to independent variation by skill: low- and high-skill destination fixed effects are weakly correlated in Dhaka and negatively correlated in Colombo, with correlation coefficients 0.17 and -0.37, respectively. 3 Validation: Model-Predicted and Survey Workplace Income in Dhaka Our first validation exercise compares income from the model and survey income from the DHUTS survey. We compute average income at the workplace level in each survey area in the DCC, the finest geographic location available in the DHUTS survey.11 Our procedure predicts that the area adjusted destination fixed effects ψ̂jR corresponds to the log labor income multiplied by e, the Fréchet shape parameter of worker unobserved preferences. Panel (A) of Table 2 shows the scatter plot between ψ̂jR against survey income, at the level of 88 survey areas. We expect a slope of e−1 and obtain a regression slope of 0.12 (with standard deviation 0.03), implying ê = 8.3. This estimate is slightly higher than previous estimates in Berlin (6.83; Ahlfeldt et al. 2015) and London in the 19th century (5.25; 11 Given that government jobs are typically paid less yet include large non-monetary benefits (such as job tenure) and are centrally located, our baseline estimation sample excludes government workers. Including them does not substantially change our results (Table G.7). Heblich et al. 2018), suggesting that idiosyncratic shocks are less important in our context. Panel (A) shows that model-predicted wages are significantly correlated with survey wages. In Panel (B), we repeat this exercise separately by skill, using the extended model described in Section 2.4. We find that model income for one skill group predicts survey income for the same skill, with a much weaker relationship across skills. Turning to predictive power, the bivariate relationship in Panel (A) has R2 = 0.25 and Root-Mean Squared Error (RMSE) = 0.22. A perfect fit would imply R2 = 1 and RMSE= 0 if wages are measured without noise in the DHUTS data. Our results might indicate that our model has poor predictive power, or that wages are measured with noise. We believe the latter concern is warranted. For example, 25% of survey areas have fewer than 50 observations, and the survey records total monthly income, which may include non-labor income. To set the right benchmark given the potential for measurement error in this validation data, we compare the model prediction with a supervised-learning approach (elastic net regularization) using 498 features from cell-phone data and geographic data. We view this exercise as the state of the art method to do prediction when training data is available (Blumenstock et al. 2015). To implement this, we randomly select half of all survey areas as “training data," and predict survey income in the other half as “test data," using either OLS or elastic net regularization. See Appendix E for details. Panel (C) of Table 2 reports the average test and training R2 and RMSE statistics over 100 random splits. Model-predicted income alone predicts 22% of the variation in the test data (column 1). The area of the tower voronoi cell, an intuitive predictor of economic activity from cell phone data, has test R2 = 0.09 (column 2).12 Including all features from cell-phone data raises test R2 to 0.24, a slight improvement over just using the model-predicted income (Column 3). The same pattern holds for RMSE. This result indicates that the model-predicted income (one statistic computed from cell phone data) summarizes nearly all information about predicting workplace income in this context, despite the parsimonious model and functional form specification for travel cost. Another important factor for R2 is geographic extent. In Section 5, we show that including 12 Cell phone operators tend to locate more towers in locations with high activity (Figure G.1). peri-urban areas drastically increases R2 .13 Severen (2019) and Tsivanidis (2019) perform validation exercises for commuting-choice models similar to ours, estimated using commuting data from survey or administrative data. Using fine divisions (census tracts), Severen (2019) finds that model wages estimated using commuting flows barely predict tract-level wages in Los Angeles. Tsivanidis (2019) finds that model-predicted wages across 19 urban areas in Bogotá predict survey wages, with an R2 = 0.36. These results suggest that the geographic aggregation level may also matter. However, in our setting, aggregating up does not change predictive power significantly (Figure G.5). Table G.5 shows that the model income is statistically significantly correlated with survey income after controlling for employment density, distance to the central business district, and log model residential income. Table G.6 shows that our results are robust to several alternate gravity equation specifications. Table G.7 uses an individual-level specification and shows that our main result is robust to controlling for workplace sorting along observable worker characteristics and other controls. 4 Validation: Model-Predicted Income and Residential Income Proxies We next use a residential income proxy constructed from population census data to validate the model prediction at the residential location level. Model-predicted residential income at tower i is defined as ∑ ψ̂jR Vij /Vi (7) j where j indexes workplace towers, ψ̂jR is the area adjusted destination fixed effect at j, Vi is total residential population at i, and Vij is the commuting volume from i to j. We focus on model fit and not the magnitude of the slope in this exercise, because the income proxy is not measured in the same units as actual income. Table 3 shows the results in Dhaka and Colombo. Model residential income is a strong predictor of the income proxy at the cell tower level (panel A). The R2 = 0.55 (Dhaka) and a machine learning application, Blumenstock et al. (2015) finds R2 = 0.41 when restricting to 37 urban DHS clusters in Rwanda, compared to 0.62 when including both rural and urban clusters. Jean et al. (2016) do not report results separately by urban areas. Within entire countries, R2 ranges between 0.37 and 0.55 for DHS-cluster level predicted consumption. 13 In 0.77 (Colombo) is high, partly because of the coverage of suburban areas. We next benchmark the predictive power to a supervised learning method (panel (B) Table 3) in Dhaka.14 The procedure is similar to Panel (B) in Table 2. Test R2 is 0.55 when using model-predicted income alone (column 1). The test R2 when using the cell phone tower Voronoi cell area alone is 0.71 (column 2), and the supervised-learning method using all features increases it to 0.73. Model-predicted residential income alone achieves about 75% of the predictive power of using all the cell-phone data metrics. 5 Practical Guidance: Spatial Coverage and Resolution What factors affect the predictive performance of our method? We explore this issue using model-predicted residential income and census income proxy data to guide future applications of this method. Our approach performs better when the analysis includes both urban core and peri-urban areas. In both cities, the explanatory power of model residential income increases as we include areas further away from the city center. Adjusted R2 goes from 0.2 at the 10 km cutoff to 0.5 and 0.7 at 30 km cutoff in Dhaka and Colombo, respectively (Figure G.3).15 Model performance is not sensitive to the level of spatial aggregation. We combine cell phone towers into square grid cells, and estimate the gravity equation at this level. In both cities, regressing the census income proxy on model-predicted residential income yields stable adjusted R2 between 0.6 and 0.9 for grid cells of between 2 and 10km wide (Figure G.5). (Gravity distance slopes are also stable.) Overall, the method is best suited to applications that cover a large urban area, and is likely not sensitive to the spatial aggregation unit. 14 Unfortunately, we do not have access to features of cell phone data to implement the supervised learning approach in Colombo. 15 Related, average residential income from the (DHUTS) survey data – only available in Dhaka’s urban core – is difficult to predict, both using model residential income, and using other measures such as distance to CBD or residential density.