Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
102,109 characters · 18 sections · 18 citation commands
Aggregation Bias in Proxy Measurement: Nighttime Lights and Local Economic Activity
JEL Classification Numbers: C23, C52, R12, R15
Keywords: aggregation bias, attenuation bias, spatial aggregation, proxy measurement, out-of-sample validation, Black Marble VIIRS
{\scriptsize Acknowledgements: We thank Piotr Wójcik and the participants at presentations given at ERSA 2023, SEA 2024, SIE 2024, AISRE 2025, IBEO 2026, SEW 2026 for useful comments and suggestions. The usual disclaimer applies. The authors have been supported by the Italian Ministry of University and Research (MIUR), in the framework of PRIN project 2017FKHBA8 001, by the project MAPPE (GRINS, PNRR), and by the University of Pisa, in the framework of the PRA Project PRA_2022_86.}
Many empirical questions require local measures of economic activity at a spatial resolution at which official statistics are unavailable, incomplete, or not comparable across countries. A common response is to use a high-resolution signal, satellite imagery, remotely sensed data, or other spatially detailed measurements, as a proxy for the missing economic variable. This strategy creates a distinctive econometric problem. The signal is typically generated by the underlying economic outcome, but the empirical objective is to use the signal to predict that outcome. Moreover, the signal is observed on a fine spatial grid and then aggregated to administrative units whose size, shape, and internal heterogeneity differ across places. Proxy measurement is therefore not only a question of predictive fit; it is also a question of reverse regression, measurement error, and spatial aggregation.
Nighttime lights (NTL) provide a leading example of this problem. Since the pioneering contribution of nordhaus2006geography, an expanding literature has used satellite-recorded luminosity to augment official income measures, study growth at sub- and supranational scales henderson2012measuring, audit GDP statistics where institutions are weak martinez2022much, measure poverty and wealth where household data are scarce jean2016combining,abbes2024deepwealth, and construct high-resolution economic maps. The appeal is clear: NTL are global, spatially detailed, comparable across borders, and repeatedly observed. They appear to provide exactly the local and frequent information that official statistics often lack.
However, “using lights as a proxy” is not a single empirical operation. Lights have been used to predict levels of GDP or income, growth rates, rankings, poverty, wealth, and local well-being henderson2012measuring,galimberti2020forecasting,jean2016combining,abbes2024deepwealth,huber2024economic, and the light--activity relationship varies with the spatial support (country to grid cell) and with development, sectoral composition, electrification, informality, geography, and time galimberti2020forecasting,bluhm2022can,gibson2024luminosity,lehnert2023proxying,huber2024economic. Empirically, estimated GDP--NTL elasticities vary widely across studies, differ across countries, become unstable at large geographies, and change across aggregation levels and NTL sources galimberti2020forecasting,bluhm2022can,gibson2024luminosity; in developing settings, where large shares of territory may be unlit, NTL may also carry non-classical measurement error huber2024economic. The relevant question is therefore not whether lights are a universally valid proxy, but what they proxy well, at which scale, in which development context, and with how much uncertainty.
This paper argues that part of this instability is the outcome of two biases that are usually discussed separately. First, although economic activity generates lights, empirical proxy regressions typically place lights on the right-hand side. The resulting coefficient is therefore a predictive reverse-regression elasticity, not the structural elasticity of light with respect to activity.\footnote{We use “reverse regression” descriptively, to signal that the direction of the fitted regression (activity on lights) is the reverse of the data-generating direction (lights emitted by activity). This is distinct from the classical reverse-regression device of wald1940 and durbin1954, in which the regression of the regressor on the outcome is used to bound an errors-in-variables coefficient. Here the reverse direction is the object of interest, and the attenuation it induces is a feature to be characterised, not a bound to be exploited.} When lights are an imperfect signal of activity, this predictive elasticity is attenuated by the noise component of the luminosity signal. Second, observed units are spatial aggregates of more elementary locations. If the elementary light--activity relationship is nonlinear, the elasticity estimated on aggregates differs systematically from the elementary elasticity. We characterise these two forces in a single probability-limit decomposition. The decomposition shows that reverse-regression attenuation is strongest at fine scales, whereas spatial aggregation becomes more important at coarser scales and contracts elasticities toward one. Unit elasticity is therefore the only aggregation-invariant benchmark.
This result links the NTL literature to classic problems in spatial statistics and econometrics: the change-of-support and modifiable areal unit problems, under which estimates depend on the spatial support at which variables are observed or aggregated cressie1996change; and the long-recognised distinction between measurement error, grouping, and aggregation: suitably constructed grouping estimators can address errors-in-variables bias, while grouping need not itself introduce additional bias wald1940,durbin1954; by contrast, aggregation over heterogeneous units can cause aggregate coefficients to differ from the underlying micro-level parameters theil1954,stoker1993. In the NTL setting these mechanisms coexist, so the theoretical framework is not specific to lights: any finely observed signal used to measure an economic variable on heterogeneous administrative units inherits the same tension between signal noise, regression direction, and spatial support.
The empirical analysis combines high-resolution Black Marble VIIRS data with official local GDP or income data for Brazil, Italy, the United States, Indonesia, and Kenya over 2012--2019, the common period over which consistently processed nightlight data and comparable local GDP or income series are available. Throughout the paper, “local economic activity” denotes the real outcome available at the relevant spatial support. This is local GDP for Brazil, the United States, Kenya, and Indonesia, and taxable personal income for Italian municipalities. We therefore interpret the Italian municipal estimates as income-density estimates, while using NUTS2 and NUTS3 GDP data to assess whether the income-based elasticity is informative about GDP-based activity at coarser spatial supports. The sample includes high-income settings where local accounts are relatively strong, an upper-middle-income economy with large territorial disparities, a lower-middle-income archipelagic economy observed at two nested administrative scales, and a lower-income agricultural economy where lights are most likely to miss part of production.
The results support a conditional, rather than universal, use of NTL. In Italy and the United States, and approximately in Brazil, the density elasticity of economic activity with respect to NTL is close to one at the finest scales, so proportional changes in lights track proportional changes in local GDP or income density and level prediction is feasible once the light-to-activity conversion is locally anchored. Indonesia and Kenya differ: their elasticities are substantially below one and their level conversions much higher. Indonesia, observed at two nested scales (kabupaten/kota and provinces), is especially informative---the estimated elasticity rises toward one under aggregation, as the mechanism predicts when the finest-scale elasticity is below one. A near-unit elasticity is therefore an empirical property of some contexts, not a mechanical feature of the data.
The out-of-sample exercises show both the promise and the limits of the proxy. Cross-sectional, cross-scale, and temporal exercises show that NTL improve prediction over models with only area and time effects, but errors widen in the lowest luminosity ranges and transferring coefficients across scales requires re-anchoring the intercept. We translate the benchmark-country evidence into a pooled finest-level proxying rule for Brazil, Italy, and the United States: a transparent first approximation for developed and upper-middle-income economies with near-unit elasticity, not to be applied mechanically to lower-income, agricultural, informal, or weakly illuminated settings without further calibration.
In summary, the paper makes three contributions. First, it provides an analytical framework for proxy regressions based on aggregated high-resolution signals. The main theorem decomposes the probability limit of the predictive elasticity into the elementary elasticity, a reverse-regression attenuation component, and a spatial aggregation component driven by unit size and within-unit dispersion. Second, it uses Monte Carlo experiments, in which the elementary relationship is known and aggregation is controlled, to assess the finite-sample implications of the theorem. The simulations show when attenuation or aggregation dominates, why density specifications with area controls reduce but do not eliminate aggregation bias, and why slopes and intercepts have different transferability properties. Third, it evaluates these implications empirically using a common NTL source, a common econometric specification, and official local economic accounts for five heterogeneous countries.
The paper is organised as follows. Section (ref) develops the methodological framework and reports the Monte Carlo evidence. Section (ref) implements the framework empirically: it describes the data, estimates the predictive elasticity and level conversion across countries and scales, evaluates out-of-sample transferability, and reports the pooled finest-level calibration. Section (ref) concludes. The Appendix contains the proof, while the Online Appendix contains Monte Carlo tables, robustness checks, additional validation results, and supplementary descriptive material.
This section develops the econometric framework that underpins the empirical analysis. The problem is not only that nighttime lights are an imperfect signal of economic activity. It is that the signal is generated at a fine spatial support, observed with noise, aggregated to heterogeneous administrative units, and then used in the reverse direction to predict the economic variable that generated it. The object of interest is therefore a predictive elasticity: the slope of the best linear projection of local economic activity on observed luminosity at a given spatial support. This coefficient is useful for proxy construction and out-of-sample validation, but it should not be interpreted as a structural light-generation parameter or as a causal effect of lights on economic activity.
The framework proceeds in three steps. First, we define an elementary-level light--activity relationship. Second, we derive the aggregate relationship that arises when elementary locations are grouped into larger spatial units. This yields a theorem on the probability-limit decomposition of the predictive elasticity into reverse-regression attenuation and spatial aggregation bias. Third, we use Monte Carlo experiments to assess the finite-sample relevance of the theorem and to motivate the density-with-area specification estimated on the real data in Section (ref).
The framework separates two possible sources of bias that are otherwise confounded in regressions of economic activity on nighttime lights. The first is a reverse-regression bias: lights are generated by activity with an idiosyncratic component, but the empirical specification uses lights to predict activity. The second is an aggregation bias: sums of nonlinear elementary relationships need not preserve the elementary elasticity once locations are combined into heterogeneous administrative units.
\paragraph{Elementary relationship and aggregation}
Consider $K$ elementary locations grouped into $P$ spatial units $m_p$, with $\abs{m_p}$ denoting the number of elementary locations contained in unit $p$. Nighttime light emissions follow a log-linear function of output with an independent shock:
where $\mu$ is the true elasticity, $\phi$ is the scale parameter, and $\eps_i$ an idiosyncratic shock. At the elementary level, write the inverse constant-elasticity relationship as
where $y_i$ and $ntl_i$ denote output and nighttime lights in elementary location $i$, $\mu>0$ is the elementary elasticity of output with respect to lights, $\phi$ is a common scale parameter, and $\eps_i$ is an idiosyncratic output-to-light shock. We assume that the shocks $\eps_i$ are i.i.d. with $\EE{\eps_i}=0$ and $\mathrm{Var}(\eps_i)=\sigma_\eps^2$, and that $\eps_i$ and $y_i$ are independent. Equivalently, with $\nu_i\equiv\log y_i$, the physical direction is
This ordering matters. The maintained exogeneity condition is imposed in the light-generation equation, not in the reverse predictive regression. Hence, even when lights are generated from activity with an exogenous idiosyncratic component, regressing activity on observed lights produces a reverse-regression attenuation term whenever lights contain such a component.
Let $Y_p=\sum_{i\in m_p}y_i$, $NTL_p=\sum_{i\in m_p}ntl_i$, and $s_{ip}=ntl_i/NTL_p$ be the within-unit light share. Since $ntl_i=s_{ip}NTL_p$,
Aggregation therefore affects the slope of a linear regression only through $\Lambda_p$, a term that depends on within-unit light shares and on the output-to-light shocks. For a cross section with $P$ aggregate units, let $\hat\beta_P$ denote the OLS slope from regressing $\log Y_p$ on $\log NTL_p$ across units $p=1,\ldots,P$ at the chosen aggregation scale. Under the usual large-cross-section regularity conditions, as $P\to\infty$,
Throughout this subsection, $\mathrm{Cov}(\cdot,\cdot)$ and $\mathrm{Var}(\cdot)$ refer to the corresponding cross-unit population moments at the chosen aggregation scale. Eq. (ref) is the central object: the overall bias is the projection of the aggregate residual $\Lambda_p$ on aggregate lights.
\paragraph{Decomposing attenuation and aggregation} To compute and interpret this projection, we decompose both objects entering the covariance: the residual $\Lambda_p$ and the regressor $\log NTL_p$. This separates the two sources of bias: the aggregation channel, which comes from summing nonlinear elementary relationships, and the reverse-regression channel, which comes from the fact that lights are generated by activity but used as a predictor of activity.
First, decompose the aggregate residual. Define
Since
we have the exact algebraic split
The aggregation component depends only on within-unit light shares and captures the mechanical effect of aggregation. The noise component captures the contribution of output-to-light shocks to the aggregate residual.
Second, the regressor splits, to first order, into a noise-free activity index and an aggregate shock,
where $A_p\equiv\log\!\left(\sum_{i\in m_p}\exp\{(\nu_i-\phi)/\mu\}\right)$ is the aggregate light that would be observed in the absence of shocks, $\tilde\eps_p\equiv\sum_{i\in m_p}s_{ip}^{(0)}\eps_i$ is the first-order aggregate shock entering the regressor, and $s_{ip}^{(0)}=\exp\{(\nu_i-\phi)/\mu\}\big/\sum_{j\in m_p}\exp\{(\nu_j-\phi)/\mu\}$ is the noise-free within-unit light share; the expansion is derived in the Appendix. Under the maintained independence between activity and shocks, $A_p$ is orthogonal to $\tilde\eps_p$ to first order.
Combining the exact residual split with the local regressor split, the numerator in (ref) can be read, to first order, as four terms:
The first and fourth terms are the two leading components. To see their approximate magnitudes, define the squared coefficient of variation of elementary-location lights within unit $p$ as
A second-order expansion of $\Lambda_p^{\mathrm{agg}}$ around equal within-unit shares gives
so that, normalised by $\mathrm{Var}(\log NTL_p)$, the aggregation--activity term is approximately $(1-\mu)\left[\delta_P-\tfrac{\mu}{2}\rho_P\right]$, with the finite-$P$ projection coefficients
Here $\delta_P$ measures the unit-size channel: larger aggregates mechanically contain more elementary locations and therefore more total lights. The coefficient $\rho_P$ measures how within-unit light dispersion covaries with aggregate lights.
The fourth term is the reverse-regression channel. With the weighted aggregate shock entering the residual, $\bar\eps_{p,w}\equiv\sum_{i\in m_p}w_{ip}\eps_i$, define the cross-sectional attenuation coefficient at the chosen aggregation scale as
The numerator pairs the shock aggregator entering the residual, $\bar\eps_{p,w}$, with that entering aggregate lights, $\tilde\eps_p$; the denominator is, to first order, $\mathrm{Var}(\log NTL_p)\approx\mu^{-2}\left[\mathrm{Var}(\mu A_p)+\mathrm{Var}(\tilde\eps_p)\right]$. At the elementary level $\bar\eps_{p,w}=\tilde\eps_p=\eps_i$ and $\kappa_P=\sigma_\eps^2/(\sigma_\nu^2+\sigma_\eps^2)$ (with $\sigma_\nu^2\equiv\mathrm{Var}(\nu_i)$), the classical noise-to-total-variance ratio. Since shocks enter the residual through $\bar\eps_{p,w}$ but the regressor through $-\tilde\eps_p/\mu$, the normalised noise--noise term in (ref) equals $-\mu\kappa_P$ to first order (see the Appendix): the two powers of $\mu$ from normalising $\mathrm{Var}(\log NTL_p)$ cancel the $1/\mu$ of the regressor expansion, so attenuation is proportional to $\mu$, exactly the errors-in-variables slope $\mu(1-\kappa_P)$.
The remaining two terms in (ref), aggregation--noise and noise--activity, are cross terms. Under the maintained independence between activity and output-to-light shocks, they vanish to first order, up to higher-order interactions between within-unit dispersion and shocks. Hence:
The reasoning is based on a local approximation: it combines an exact decomposition of $\Lambda_p$ with a first-order expansion of $\log NTL_p$ and a second-order expansion of $\Lambda_p^{\mathrm{agg}}$.
The following result formalises this approximation.
Before proceeding, we make some remarks on Theorem (ref).
\paragraph{Implications for the empirical analysis}
Equation (ref) has three implications for the empirical analysis. First, in the noiseless environment (i.e.\ no attenuation bias), $\mu=1$ is the unique elasticity that is invariant to arbitrary aggregation patterns: when $\mu=1$, both the unit-size term and the within-unit-dispersion term disappear. For $\mu\neq1$, aggregation generally changes the estimated elasticity, with the sign governed by the bracket $\delta-\mu\rho/2$. In the leading case in which the unit-size channel is positive and dominates within-unit dispersion, aggregation pushes estimates upward when $\mu<1$ and downward when $\mu>1$, contracting them toward unity.
Second, idiosyncratic output-to-light shocks attenuate the reverse predictive slope, but this attenuation weakens as aggregation averages out the shock component. Expressing variables as densities, or including area controls in levels, targets the first-order unit-size channel $\delta$; what remains is the reverse-regression attenuation and the smaller within-unit-dispersion component. The Monte Carlo exercises below verify these mechanisms in controlled environments.
Third, under the maintained assumption that reverse-regression attenuation is negligible at the relevant scale ($\kappa_\infty\approx0$), Eq. (ref) can be inverted to recover an aggregation-corrected estimate of the elementary elasticity $\mu$ from the observed predictive slope $\hat\beta_P$ and empirical estimates of $\delta_P$ and $\rho_P$ computed from the observed NTL distribution. We use this inversion only as a diagnostic, because it relies on the strong assumption that the unobserved output-to-light shock is small.
In each simulation run, we generate data for $K = 100{,}000$ elementary locations, indexed by $i$ as in the analytical benchmark. Elementary output is drawn from a log-normal distribution, $y_i=\exp(\nu_i)$ with $\nu_i\sim\mathcal{N}(1,1)$. Nighttime light emissions follow Eq. (ref), with an independent shock $\eps_i\sim\mathcal{N}(0,\sigma_\eps^2)$ and scale parameter $\phi=11.5$. Elementary locations are then aggregated into $P\in\{100,1{,}000,10{,}000,20{,}000,50{,}000,100{,}000\}$ spatial units with heterogeneous elementary-location counts. For each aggregation level we estimate the regression covered by Theorem (ref):
where $Y_{p}$ and $NTL_{p}$ are total output and total lights of spatial unit $p$. Sections (ref)--(ref) study this specification, whose bias is fully characterised by Equation (ref); Section (ref) then explains why the correction implied by the theorem is not directly operational, introduces unit areas, and evaluates the empirical specifications used in Section (ref).
The first Monte Carlo exercise verifies the finite-sample accuracy of Eq. (ref). For each combination of $\mu\in\{0.7,1,1.3\}$ and $\sigma_\eps\in\{0.1,0.3,0.5,1\}$, and for each aggregation ratio $K/P$ with $P\in\{100,1{,}000,10{,}000,20{,}000,50{,}000,100{,}000\}$, we regress $\log Y_p$ on $\log NTL_p$ and compare the simulated bias $\hat\beta_P-\mu$ with the approximation $-\mu\kappa_P+(1-\mu)\left[\delta_P-\tfrac{\mu}{2}\rho_P\right]$. The projection coefficients $\delta_P$ and $\rho_P$ are computed from the simulated lights using Eq. (ref), while $\kappa_P$ is computed from Eq. (ref) using the true shocks $\eps_i$ and the true $\mu$, so that the only discrepancy left between the two sides is the remainder $R(\chi,\sigma_\eps)$ of Theorem (ref).
Figure (ref) summarises the comparison; the full set of estimates is reported in Section (ref) of the Online Appendix (Table (ref)). Three findings stand out. First, at the elementary scale ($K/P=1$) the approximation is essentially exact for every $(\mu,\sigma_\eps)$: with no aggregation, $\delta_P=\rho_P=0$ and the bias reduces to pure attenuation, $-\mu\kappa_P$. Second, for moderate noise ($\sigma_\eps\le0.5$) and moderate aggregation ($K/P\le10$), the remainder is uniformly small (below $0.04$ in absolute value, and typically below $0.02$), so Eq. (ref) tracks both the attenuation and the aggregation components of the bias remarkably well. Third, the approximation deteriorates exactly where the theorem says it should: for $\sigma_\eps=1$ the remainder grows with the $O(\sigma_\eps^2)$ term, and for $\mu=0.7$ at coarse scales ($K/P\ge100$) the within-unit light shares become highly unequal --- the variance of elementary log lights is $\sigma_\nu^2/\mu^2$, so low $\mu$ inflates within-unit dispersion --- and the local expansion in $\chi$ is no longer accurate, with $\rho_P$ taking large values. In this region the approximation still gets the sign and the qualitative contraction toward one right, but no longer the magnitude. The working range of Eq. (ref) therefore covers the empirically relevant configurations in which the calibrated light signal is not dominated by measurement noise and administrative aggregation is not at the extreme coarse limit of the simulation design.
Figures (ref) and (ref) decompose the approximation into the two components of Theorem (ref), computed separately in the same experiments. The decomposition shows that the two channels operate at opposite ends of the aggregation range. The attenuation component $-\mu\kappa_P$ is maximal at the elementary scale and vanishes monotonically as $K/P$ grows: within-unit averaging of the output-to-light shocks raises the signal content of aggregate lights, with a profile that is common across $\mu$ up to the scale factor. The aggregation component $(1-\mu)\left[\delta_P-\tfrac{\mu}{2}\rho_P\right]$ behaves in the opposite way: it is zero at the elementary scale, grows in magnitude with $K/P$, and is identically zero at every scale when $\mu=1$ --- the knife-edge property that makes the unit elasticity invariant to aggregation. Its sign is positive for $\mu<1$ and negative for $\mu>1$, generating the contraction toward one.
The crossing of the two profiles also explains the non-monotone total bias visible in Figure (ref) for high noise levels: for $\mu=0.7$ and $\sigma_\eps=1$ the bias moves from $-0.35$ at the elementary scale, where it is pure attenuation, to $+0.21$ at the coarsest scale, where the aggregation channel dominates.
We set $\sigma_\eps=0.1$ (i.e.\ $\sigma^2_\eps=0.01$), so that attenuation is small relative to aggregation. At the elementary scale, $\kappa_P=\sigma^2_\eps/(\sigma^2_\nu+\sigma^2_\eps)\approx0.01$, implying an attenuation bias of about one percent of $\mu$ in Eq. (ref).
Figure (ref) in the Online Appendix reports the Monte Carlo estimates of Model (ref) across aggregation levels. It confirms the analytical benchmark: aggregation leaves the elasticity unchanged at $\mu=1$, pushes estimates upward when $\mu<1$, and pushes them downward when $\mu>1$. Thus aggregation contracts the estimated elasticity toward one rather than merely adding noise. The bottom row shows the implication for level recovery: the scale parameter $\alpha$ is well recovered at the unit-elastic benchmark at every scale, but when $\mu\neq1$ the size-driven distortion of $\hat\beta$ is absorbed by the intercept, which drifts away from $\phi$ as aggregation coarsens. Level prediction from lights therefore requires an anchored light-to-output conversion whenever the elasticity is away from one. Full estimates and confidence bands are reported in Section (ref) of the Online Appendix (Table (ref)).
Figure (ref) in the Online Appendix varies the idiosyncratic-shock standard deviation $\sigma_\eps$ to put attenuation and aggregation at work jointly in Model (ref). At the finest scale, the slope is $\hat\beta \approx \lambda\mu$, with $\lambda \equiv \sigma_\nu^2/(\sigma_\nu^2+\sigma_\eps^2)$; noisier light therefore implies stronger reverse-regression attenuation. As aggregation increases, the slope is pulled toward one from whatever attenuated starting point: the attenuation component dies out while the aggregation component takes over, exactly the crossing of the two profiles documented in Figures (ref) and (ref). Section (ref) in the Online Appendix reports the underlying numbers (Table (ref)).
Theorem (ref) characterises the bias of Model (ref), but the correction it implies is not directly operational. Undoing the bias requires the attenuation coefficient $\kappa_P$, which depends on the unobserved output-to-light shocks $\eps_i$, and the projection coefficients $\delta_P$ and $\rho_P$. Without auxiliary information, the theorem describes the bias but cannot remove it. Section (ref) shows what can be achieved when partial auxiliary information is available (grid cells as elementary units and $\kappa_\infty\approx0$ as a maintained assumption), together with the failures of that route at coarse scales; here we take the complementary, fully operational route: specifications that mitigate the observable aggregation channels by construction. In particular, the empirical specification is written in densities and includes an area control. Under the log-log specification a totals regression with an area control and a density regression with an area control are equivalent; we adopt the density form for comparability across aggregation scales.
We therefore introduce unit land areas $S_p$ into the simulations. Area is the observable counterpart of the unit-size channel: if elementary locations had equal area, $S_p$ would be proportional to $\abs{m_p}$, and an area control would target the $\delta_P$ term of Theorem (ref). To allow for imperfect proportionality, elementary areas are drawn separately. In the baseline, $S_p$ is positively correlated with the elementary-location count, as in the natural case where larger administrative units contain more land and more light-emitting locations. We also consider independent and negatively correlated areas. The negative-correlation case is a useful stress test for settings with compact, dense urban units and large, sparsely populated rural units. Two alternatives to Model (ref) are considered:
which expresses variables as spatial densities and adds an explicit area control, and
which further allows a luminosity shift $\ell_0$ estimated from the data. These are the Monte Carlo counterparts of the Baseline Model and of the Nonlinear Model of Section (ref).
\paragraph{Bias in the estimates}
Figure (ref) reports the in-sample evidence for the area-control specification, with variables in totals or densities and under the three area--size correlation regimes. The first row shows that the density specification with an area control keeps $\hat\beta$ close to the true elasticity at every aggregation level, removing the contraction toward one that afflicts the totals regression. The second row reports recovery of the scale parameter, and the third row reports the residual role of unit size after density normalisation. The area control is therefore substantive: when area tracks the aggregation support, it prevents the NTL coefficient from absorbing variation due to observational-unit size. Adding noise to the picture, under aggregation the naive $\hat\beta/\mu$ is pulled toward the fixed point $1/\mu$, while the density-and-area ratio stays close to the reliability ratio $\lambda$, with residual drift only at the highest noise levels. Densities and the area control therefore largely neutralise the first-order aggregation channel, but not the attenuation channel; attenuation fades only as aggregation averages the shocks out (Figure (ref) in the Online Appendix).
\paragraph{Out-of-sample prediction} Finally, the specifications are compared on the criterion that matters for using lights as a proxy: out-of-sample prediction. Within each replication, units are split into five folds. Each specification is estimated on four folds and used to predict the fifth. Performance is measured by the root mean squared out-of-fold error, which is comparable across specifications because predicting $\log Y_p$ or $\log(Y_p/S_p)$ is equivalent once area is observed. To isolate the pure aggregation channel, the independent and negative area--size regimes are evaluated at low noise ($\sigma_\eps=0.1$).
Table (ref) reports the results. Three facts stand out. First, the density-with-area specification weakly dominates the totals regression in every regime. It coincides with the totals regression wherever area carries no size information, and predicts systematically better wherever it does. The gains are largest at coarse scales and high noise under positive area--size correlation: for example, at $\mu=0.7$, $\sigma_\eps=0.5$, and $K/P=1000$, the RMSE falls from $0.135$ to $0.082$. The same logic applies under negative correlation: at $\mu=0.7$ and $K/P=10$, the RMSE falls from $0.183$ to $0.150$. The area control exploits the area--size relation in either direction.
Second, density normalisation alone is not sufficient. The density specification without the area control is dominated at fine scales whenever normalisation injects area noise into the regression. Removing the unit-size channel requires the control, not just the transformation. Conversely, a totals specification with the area control is not reported because it is observationally equivalent to the density specification with the area control: the two share the same regressor span and their out-of-fold predictions differ exactly by the observed $\log S_p$, so the errors coincide. What matters econometrically is controlling for area; densities remain preferable for interpretation and comparability.
Third, the Nonlinear Model performs exactly like the Baseline Model in this DGP\@. Since there is no luminosity floor, the estimated shift $\ell_0$ correctly collapses toward zero, so allowing for it costs nothing out of sample. Its value emerges on real data at low luminosity, as documented in Section (ref). At $\mu=1$ all specifications are equivalent, the knife-edge case in which aggregation is harmless. These results provide the controlled-environment justification for adopting the density-with-area specification as the Baseline Model of Section (ref), with the Nonlinear Model as a flexible extension.
This section takes the aggregation framework to the data. We first describe the harmonised light and economic panel data and construct empirical counterparts of the aggregation terms in Theorem (ref) (Section (ref)). We then estimate the predictive light--activity relationship at each available spatial support (Section (ref)). The final two subsections assess transferability: out-of-sample validation tests whether estimated relationships predict held-out units, scales, and years (Section (ref)), and a pooled finest-level calibration evaluates how far a common proxying rule can be pushed in benchmark economies (Section (ref)).
We combine a common high-resolution signal, nighttime lights (NTL), with official local measures of GDP or personal income for Brazil, Italy, the United States, Indonesia, and Kenya. The sample period is 2012--2019, the common window over which consistently processed VIIRS data and the required local economic series are available, and it stops before the COVID-19 disruption. The countries are deliberately heterogeneous: the United States and Italy are high-income benchmarks, Brazil is an upper-middle-income economy with large territorial disparities, Indonesia is a lower-middle-income archipelagic economy, and Kenya is a lower-income agricultural economy where lights are more likely to miss non-electrified, agricultural, or informal production.
NTL are sourced from NASA's Black Marble project,\footnote{\url{https://blackmarble.gsfc.nasa.gov/}.} based on VIIRS observations from the Suomi NPP satellite at approximately 500-meter resolution. Black Marble corrects for atmospheric effects, terrain, vegetation, snow cover, lunar illumination, stray light, saturation, and calibration issues roman2018nasa. These corrections are important for the proxy problem because measurement noise contributes directly to reverse-regression attenuation. Raw NTL maps are reported in Section (ref) of the Online Appendix.
For every country, year, and spatial support, we construct total NTL by summing grid-cell radiance within administrative boundaries. Brazil is observed at municipality, microregion, mesoregion, and federal-state levels; Italy at municipality, local labour area, NUTS3, and NUTS2 levels; the United States at county, commuting-zone, and state levels; Kenya at county level; and Indonesia at kabupaten/kota and province levels. For compact notation, we use MUN for municipality, MICRO for microregion, MESO for mesoregion, STATE for federal state or U.S. state, LLA for local labour area, CZ for commuting zone, COUNTY for county, KAB/KOTA for kabupaten/kota, and PROV for province. NUTS2 and NUTS3 denote European Union statistical regions at levels 2 and 3. Figure (ref) in the Online Appendix reports NTL density and economic density at the lowest available administrative level for each country.
Local economic activity is measured from official sources. Brazilian municipal GDP comes from IBGE and is aggregated to microregions, mesoregions, and federal states.\footnote{\url{https://www.ibge.gov.br/estatisticas/economicas/contas-nacionais/9088-produto-interno-bruto-dos-municipios.html}.} U.S. county- and state-level GDP comes from the BEA regional accounts.\footnote{\url{https://apps.bea.gov/iTable/}.} Italian municipal personal income is derived from Agenzia delle Entrate tax declarations, while NUTS2 and NUTS3 GDP comes from ARDECO.\footnote{\url{https://www1.finanze.gov.it/finanze/pagina_dichiarazioni/public/dichiarazioni.php}. For ARDECO, see \url{https://knowledge4policy.ec.europa.eu/territorial/ardeco-database_en\#database}.} Kenya's Gross County Product comes from the Kenya National Bureau of Statistics,\footnote{\url{https://www.knbs.or.ke/all-reports/}.} and Indonesian real GDP at province and kabupaten/kota levels comes from INDO-DAPOER.\footnote{Indonesia Database for Policy and Economic Research, World Bank DataBank: \url{https://databank.worldbank.org/source/indonesia-database-for-policy-and-economic-research}.} The effective time coverage and unit counts reported below refer to the merged analytical panel after sample restrictions.\footnote{Alaska, Hawaii, and U.S. territories are excluded from county-level maps and regressions because their extreme geography and low population density would distort county-level estimates. State-level regressions retain aggregate state observations, including the District of Columbia where available, as coarse aggregation benchmarks.}
Monetary variables are first expressed in 2015 prices using CPI or GDP deflators and are then converted into international dollars using World Bank PPP conversion factors. This harmonisation supports cross-country comparison of levels, while the main elasticity estimates rely on within-country variation. Where available, we also apply within-country cost-of-living adjustments: NUTS2 PPP data for Italy from cannari2009differenze, and state-level PPP data for the United States from the BEA\@. No comparable intra-country PPP adjustment is available for Brazil, Kenya, or Indonesia. The baseline regressions use density measures: NTL, GDP, income, and area are all harmonised at the same administrative support, with economic activity and NTL divided by square kilometres.
As defined in Section (ref), “local economic activity” is the real outcome available at the relevant spatial support: GDP for Brazil, the United States, Kenya, and Indonesia, and taxable personal income for Italian municipalities. Only for Italy does the finest-level variable differ from GDP; we therefore read the Italian municipal estimates as income-density estimates and use the NUTS2 and NUTS3 GDP data to check whether they carry over to GDP-based activity at coarser supports (Section (ref)).
Table (ref) reports descriptive statistics for the analytical samples. The table also provides a first comparison of economic density across countries after PPP conversion.\footnote{As a consistency check, we aggregate local monetary values to the country-year level and compare them with World Bank GDP in current international dollars (NY.GDP.MKTP.PP.CD). The ratios are close to one for Brazil and Indonesia, about 0.94 for the contiguous United States, consistent with the exclusion of Alaska, Hawaii, and U.S. territories, and about 0.91 for Kenya, where Gross County Product does not exactly match national GDP\@. For Italy, the municipal variable is personal income rather than GDP; aggregating municipalities gives about 42--45 percent of GDP in current international dollars, i.e.\ roughly one half.} Within-country dispersion is large, especially at fine scales, and cross-country comparisons should be interpreted with caution where subnational cost-of-living adjustments are unavailable. Administrative labels are also not directly comparable across countries: municipalities in Brazil and Italy, counties in the United States and Kenya, and kabupaten/kota in Indonesia are different spatial and institutional objects.
We next construct empirical diagnostics for the aggregation terms in Theorem (ref). In particular, Table (ref) reports three empirical counterparts of the theory. For this diagnostic only, the elementary units are lit Black Marble grid cells. The variables below continue to aggregate total radiance within administrative boundaries; the lit-cell restriction is used only to compute within-unit light shares and coefficients of variation in a way that corresponds to the positive-light elementary locations of Section (ref). Including fully dark cells would make $\mathrm{CV}_p^2$ dominated by large, sparsely lit units and would turn the diagnostic into a measure of darkness rather than of dispersion among light-emitting locations. The cell-level variance of log radiance, $\sigma_n^2$, is a proxy for the variance of elementary (log) local economic activity under the light-generation equation when measurement error is negligible. The projection coefficient $\delta_P$ captures the unit-size channel: luminous aggregate units tend to contain more lit elementary locations. The projection coefficient $\rho_P$ captures the within-unit dispersion channel: conditional on luminosity, units with more unequal light shares can have different aggregate elasticities.
Two regularities are important for the estimates below. First, $\delta_P$ is positive and sizeable in all five countries, already at the finest available administrative levels. The unit-size channel is therefore empirically relevant even before moving to coarse regions. Second, $\rho_P$ is heterogeneous. It is moderate in Brazil, Italy, the United States, and Kenya, but large in Indonesia, where archipelagic units combine bright urban cores with extensive dim areas; it can also turn negative at coarse scales, as in Brazilian mesoregions, where large Amazonian units combine low total light with high dispersion. Because both coefficients enter the aggregation component multiplied by $(1-\mu)$, they generate little bias when the elementary elasticity is close to one, but they contract estimated elasticities toward unity when $\mu$ differs from one. Section (ref) uses these diagnostics to interpret the country and scale patterns of $\hat\beta$.
We now estimate the predictive NTL--activity relationship across all available countries and administrative levels. The baseline specification is
where $\text{EconomicActivity\_km2}$ is GDP density per square kilometre for Brazil, the United States, Indonesia, and Kenya, and personal-income density per square kilometre for Italy. $\text{NTL\_km2}$ is nighttime lights density per square kilometre; $\alpha_t$ is a time-varying intercept, $\beta$ is the NTL elasticity, and $\gamma$ captures the residual role of observational-unit area.
The three parameters of this Baseline Model (BM) map directly into the proxy problem. The intercept $\alpha_t$ anchors levels. The elasticity $\beta$ governs proportional variation: when it is close to one, NTL density approximates relative differences and growth rates in economic density. The area coefficient $\gamma$ implements the aggregation correction motivated by Section (ref): if unit size is correlated with luminosity and economic density, omitting area can contaminate the NTL elasticity. The nonlinear specification introduced in Section (ref) relaxes the constant-elasticity restriction and lets the NTL elasticity vary with luminosity.
Figure (ref) reports $\hat\beta$ for all countries and administrative levels. Confidence intervals use standard errors clustered by spatial unit; Section (ref) of the Online Appendix reports the corresponding BM estimates and higher-geography clustering checks where available.
The estimates separate the benchmark economies from the lower-income cases. Italy and the United States are close to unit elasticity at the finest available scales; Brazil is slightly below one at the municipal level and moves above one at coarser levels. In these more luminous settings, proportional variation in NTL density is therefore a useful approximation to proportional variation in GDP or income density. This is a statement about economic magnitudes, not about failing to reject $\beta=1$: with large panels, small deviations can be statistically precise.
Kenya and Indonesia instead have elasticities below one, especially Kenyan counties and Indonesian kabupaten/kota. This is consistent with low luminosity, agriculture, informality, and weaker electrification, although data quality may also contribute: subnational accounts in these settings are noisier and sometimes partly imputed. Indonesia is particularly informative because it is observed at two scales. Moving from kabupaten/kota to provinces, the elasticity rises toward one, matching the aggregation benchmark for $\mu<1$; Section (ref) in the Online Appendix finds the same pattern under random aggregation of kabupaten/kota. A near-unit estimate at a coarse scale can therefore mask a lower fine-scale elasticity.
That the near-unit estimates in Brazil, Italy, and the United States are already present at the finest levels, where the density specification removes the first-order aggregation channel (Eq. (ref)), indicates they are not merely an aggregation artefact, which would appear mainly at coarse scales, but a property of the local NTL--activity relationship, with residual deviations measuring the approximation error of a one-for-one proxy. The inference is robust to clustering by unit, higher-geography clustering, and a wild cluster bootstrap (Section (ref) in the Online Appendix). The equivalence tests in Table (ref) sharpen the magnitude claim: under conservative clustering, Italian municipalities and U.S. counties are equivalent to one within $\pm 10\%$, Brazilian municipalities within $\pm 20\%$, while Kenya and Indonesia reject equivalence at either tolerance. For Italy, where the finest-level variable is income rather than GDP, we re-estimate the model at NUTS3 and NUTS2 using GDP\@. The elasticities are very similar, especially at NUTS3, supporting the use of personal income at the municipal level while keeping the GDP interpretation for coarser Italian scales.
Finally, high in-sample fit should not be overinterpreted: both luminosity and economic activity density are strongly related to population density, so fit can be high even when the proxy is imperfect. Section (ref) therefore assesses proxy quality out of sample.
The near-unit fine-scale estimates also receive indirect support from the aggregation formula. Under negligible measurement error ($\kappa_\infty\approx0$), Theorem (ref) can be inverted to recover an implied elementary, cell-level elasticity $\mu$ from the naive totals slope and the empirical projection coefficients in Table (ref). This is a diagnostic exercise rather than an estimator we rely on for calibration, because it imposes the strong assumption that reverse-regression attenuation is negligible. The inversion (Table (ref) in Section (ref) of the Online Appendix) delivers $\hat\mu\approx1.00$ for Italian municipalities, $0.91$ for U.S. counties, $0.76$ for Brazilian municipalities, and $0.70$ for Indonesian kabupaten/kota. The ordering matches the BM estimates. Kenya is informative through failure: its naive slope lies below $\delta_P$, a configuration no non-negative $\mu$ can generate under $\kappa_\infty\approx0$, signalling that measurement error is non-negligible where lights miss much agricultural and informal activity.
The time-varying intercept $\alpha_t$ is the level conversion between NTL density and economic density. It absorbs common changes in the baseline light--activity relationship over time and is the key parameter when NTL are used to predict levels rather than proportional differences.
Figure (ref) reports the estimated intercepts. At the finest available levels, Brazil, Italy, and the United States have broadly comparable conversions, whereas Indonesia and especially Kenya require much higher intercepts. In logs, these are large multiplicative differences: conditional on NTL density and area, a unit of observed luminosity corresponds to substantially more measured GDP density in the lower-luminosity economies. This is consistent with a larger agricultural, informal, rural, or otherwise low-emission component of activity henderson2018global,jean2016combining,abbes2024deepwealth. One qualification: at the finest scale Italy is measured with personal income, not GDP (municipal income is about 42--45% of GDP; Section (ref)). The like-for-like NUTS3 and NUTS2 checks suggest the mismatch shifts the Italian intercept by less than the raw income-to-GDP ratio would imply; still, the Brazil--Italy--United States level comparison should be read as approximate, and the pooled calibration (Section (ref)) handles it through country-year fixed effects.
Within countries, $\alpha_t$ is fairly stable over time. Brazil shows a mild decline, especially between 2012 and 2017, while Italy and the United States show small increases at finer levels. This stability indicates that annual recalibration is not the main difficulty; cross-country level conversion is.
Figure (ref) in the Online Appendix reports the area coefficient $\gamma$, which captures the residual association between observational-unit size and economic density after conditioning on NTL density and time effects. In the density specification, the first-order unit-size term is removed by construction, so $\gamma$ should be close to zero unless area remains correlated with residual aggregation effects.
This is what we find for Brazil, Italy, and the United States, where $\hat\gamma$ is generally small. Kenya and Indonesia instead display larger and statistically different area coefficients, suggesting that unit size and within-unit heterogeneity still matter after density normalisation. The area correction is therefore marginal in the benchmark economies but material in lower-luminosity and geographically uneven settings.
Taken together, $\beta$, $\alpha_t$, and $\gamma$ point to two calibration problems in poorer and lower-luminosity economies. Levels require a country-specific conversion factor, and growth rates require an elasticity that may be well below one. This complements gibson2024luminosity: NTL remain informative in developing settings, but they need explicit calibration for both levels and proportional changes.
We finally allow the NTL elasticity to vary with luminosity using the shifted-log Nonlinear Model (NLM):
where the luminosity shift $\ell_0$ governs curvature and satisfies $\text{NTL\_km2}_{it}+\ell_0>0$ throughout the sample. The implied local elasticity is $\beta\,\text{NTL\_km2}_{it}/(\text{NTL\_km2}_{it}+\ell_0)$. When $\ell_0>0$, elasticity is lower in the darkest units and rises toward $\beta$ as luminosity increases; when $\ell_0<0$, the opposite pattern can arise over the admissible range.
Figure (ref) in the Online Appendix reports the implied local elasticity $\varepsilon(x)=\hat\beta\,x/(x+\hat\ell_0)$, with estimates in Section (ref) of the Online Appendix. In Brazil, Italy, and the United States, elasticity is close to one around typical luminosity levels but falls in the darkest units. Kenya and Indonesian kabupaten/kota remain below one over most of their support. The nonlinear model therefore refines the linear evidence by locating the low-luminosity regime in which the constant-elasticity approximation becomes fragile.
For a proxy, the relevant test is not in-sample fit but prediction for units, scales, or years not used in estimation. We compare the Baseline Model (BM), the Nonlinear Model (NLM), and a deliberately No-Information Model (NIM), which omits NTL density and keeps only the time-varying intercept and $\log \text{Km2}$. The gain over NIM measures the information in lights beyond unit size and time effects, while the comparison between BM and NLM tests whether low-luminosity curvature improves prediction.
The validation separates two uses of the proxy. For growth rates or relative comparisons, the key object is the slope: with $\beta$ close to one, light density moves nearly proportionally with economic density. For levels, the slope is not enough: the intercept must anchor the light-to-activity conversion. The exercises below therefore ask how well slopes transfer across units, scales, and years, and how much accuracy is lost when the level correction is unavailable or observed only at an aggregate scale. This design is feasible because the model uses common time-varying intercepts rather than unit fixed effects, so held-out units can still be predicted.
We consider three complementary exercises.
(A) Cross-sectional (leave-regions-out). We hold out groups of spatial units, estimate the model on the remaining units, and predict all years of the held-out units. At the coarsest levels, where the number of units is small, this becomes leave-one-unit-out. This tests transferability across units.
(B) Cross-scale (downscaling). We estimate the model at a coarse administrative level and use the coefficients, together with fine-level NTL, to predict fine-level income or GDP density. This tests whether aggregate estimates extrapolate to finer spatial resolutions.
(C) Temporal (leave-years-out). We hold out one year and predict it using slopes estimated on the remaining years and the average estimated intercept $\bar{\hat\alpha}$. The error captures any drift in the light-to-activity conversion.
In every exercise, the forecast error is
where the three forecasting functions are \[ f_{\mathrm{BM}}(x)=\log x, \qquad f_{\mathrm{NLM}}(x)=\log(x+\ell_0), \qquad f_{\mathrm{NIM}}(x)=0, \] and $(-\ell,m)$ denotes coefficients of model $m$ estimated without held-out fold $\ell$. Predictive performance is summarised by the out-of-sample $R^2$,
by root mean square error (RMSE), and by the distribution of $e_{it}$. The increment over NIM, $\Delta R^2_{\text{oos}}$, is the information content of NTL beyond area and time effects; the gain from NLM is measured relative to BM.
Table (ref) reports the main cross-sectional validation exercise by country and administrative level.
Figure (ref) shows the distribution of temporal forecast errors conditional on NTL density at the finest available level: municipalities for Brazil and Italy, kabupaten/kota for Indonesia, and counties for Kenya and the United States. The complete cross-scale and temporal validation tables are reported in Section (ref) of the Online Appendix. The results answer the three validation questions. First, NTL add substantial predictive content over NIM. Across countries and scales, the NTL models explain much more out-of-sample variation than area and time effects alone, so predictive power is not merely a mechanical artefact of unit size. The nonlinear correction, however, changes average fit only modestly; its role is mainly in the tails of the error distribution.
Second, coarse-to-fine prediction supports downscaling, but only for the slope structure. Once the systematic level shift between scales is removed, coefficients estimated at coarse levels recover the cross-sectional pattern of finer-level density well. What does not transfer automatically is the intercept: the naive level prediction becomes biased as the distance between estimation and target scales grows. Aggregate estimates therefore extrapolate to finer resolutions in relative terms, but level prediction requires re-anchoring with some local information.
Third, the error distribution shows where the proxy degrades. Forecast errors widen in the lowest NTL deciles and are tighter elsewhere (Figure (ref)). This is why Table (ref) reports $L_4$ in addition to RMSE: the nonlinear correction is most visible in tail errors rather than average fit. These patterns mirror the curvature documented in Section (ref): the proxy is least reliable in the darkest, least dense units, where the constant-elasticity approximation begins to fail. Overall, NTL are informative out of sample, but accuracy degrades predictably at low luminosity and downscaling requires re-anchoring the level.
The validation exercises show that NTL contain out-of-sample information and identify where the proxy becomes fragile. We now ask a more operational question: can the benchmark-country evidence be summarised by a pooled finest-level calibration usable when local economic data are unavailable? This is not a new identification design, but a calibration exercise that turns the preceding estimates into a proxying rule and measures the cost of different level anchors.
We restrict the benchmark pool to Brazil, Italy, and the United States, where finest-level elasticities are close enough to unity to justify a common slope (Table (ref)). The pooled sample uses municipalities for Brazil and Italy and counties for the United States. Italy contributes income density rather than GDP density. Including it in a GDP-based pool is therefore an empirical decision, not an assumption: it rests on the finding that the income- and GDP-based elasticities are near-identical where both are observed (Italian NUTS3 and NUTS2), so that the Italian slope is informative about GDP even though the finest-level variable is income. Country-year intercepts then absorb the residual level difference between income and GDP. Kenya and Indonesia are excluded because their elasticities and level conversions are not comparable to this benchmark group.
We compare three level anchors. The first uses country-year intercepts estimated on the pooled local sample. The second imposes a single common intercept. The third keeps the pooled slopes but recovers the intercept from aggregate country-level totals, either by country-year or by country. Predicted local income or GDP density is then evaluated in log errors and by deciles of observed economic density.
The main lesson is that slopes transfer better than levels. With country-year intercepts, the pooled specification has mean errors essentially equal to zero by construction, and the nonlinear term modestly improves the error distribution (Table (ref)). A balanced-sample robustness check in Section (ref) of the Online Appendix gives essentially unchanged results, so the calibration is not driven by the larger number of Brazilian and Italian municipalities.
A common intercept is less damaging than expected in this benchmark group, but intercepts recovered only from aggregate country-level totals perform substantially worse. Thus $\beta$, $\gamma$, and $\ell_0$ are reasonably transferable within the benchmark economies, whereas an $\alpha$ calibrated at the country scale does not automatically anchor local density. Aggregate information can provide a pilot estimate, but quantitative fine-scale prediction still requires some local-level anchoring. The pooled calibration should therefore be used primarily to recover relative economic density across local units. When the objective is level prediction in a new country, the intercept must be re-anchored using aggregate national or regional accounts and, where possible, some local observations. Without this re-anchoring, the systematic level shifts documented in the validation exercises are likely to translate into biased local predictions.
This paper has studied nighttime lights as an application of a more general econometric problem: when can a high-resolution signal, observed on a fine spatial grid and aggregated to administrative units, be used to predict an economic variable that is not observed locally? The answer depends on two distinct margins. The first is reverse-regression attenuation: lights are generated by economic activity, but the empirical objective is often to predict economic activity from observed lights. The second is spatial aggregation: the relationship estimated at an administrative support need not coincide with the elementary relationship operating at finer locations.
The analytical framework clarifies how these two margins interact. The probability limit of the predictive elasticity differs from the elementary elasticity because of a noise component in the light signal and an aggregation component driven by unit size and within-unit dispersion. A central implication is that aggregation tends to contract elasticities toward unity. Unit elasticity is therefore not merely a convenient empirical regularity; it is the only benchmark under which the elasticity is invariant to aggregation. The Monte Carlo exercises confirm this decomposition. They also show that density specifications with area controls reduce the first-order unit-size channel, but do not remove all aggregation effects when the underlying relationship is nonlinear or when the light signal is noisy.
The empirical results are consistent with this conditional view of proxy validity. Using Black Marble VIIRS nighttime lights and official local GDP or income data, we find that Brazil, Italy, and the United States form a relatively stable benchmark group: at the finest available spatial supports, predictive elasticities are economically close to one and remain broadly stable across scales and over time once densities and area are used. Out-of-sample validation shows that lights add substantial predictive content beyond area and time effects alone, although prediction errors remain larger among low-luminosity units. In this group, pooled calibration is informative about relative local economic density, provided that levels are appropriately anchored.
The results for Indonesia and Kenya show the limits of mechanical transfer. Their predictive elasticities are substantially below one at fine scales, implying that differences in luminosity translate less than proportionally into differences in measured economic activity. Aggregation may partly conceal this pattern by moving estimated elasticities toward unity at coarser scales. Level conversion is even more fragile: for a given amount of light and area, lower-income, more agricultural, more informal, or weakly illuminated contexts may contain more economic activity than a benchmark-country calibration would imply.
The practical implication is that nighttime lights should be treated as a calibrated proxy rather than as a universal measure of local output. Applications should estimate the elasticity at the relevant spatial support, work with densities and unit area when possible, test whether the predictive elasticity is close enough to one for the intended use, and validate predictions out of sample. When the slope is transferable but levels are not, lights are better suited to recovering relative economic density than to producing unanchored level estimates. When the elasticity itself differs across contexts, country- or context-specific calibration is necessary.
Future work should extend this framework to longer panels and a wider set of low- and middle-income countries, where official local economic data are scarce and the value of a reliable proxy is highest. The same logic can also be applied to other high-resolution spatial signals used in economics, such as remotely sensed land use, building footprints, mobile-phone activity, or transaction-based indicators. In all these cases, the key question is not only whether the signal correlates with the target variable, but whether its predictive relationship survives aggregation, noise, and transfer across contexts.