Data
For our empirical analysis, we use data from a variety of sources.
Fires: we measure crop burning with data on fires from NASA satellites Moderate Resolution Imaging Spectroradiometer (MODIS) (1 km × 1 km resolution) and Visible Infrared Imaging Radiometer Suite (VIIRS) (375 m × 375 m resolution) from September 2012, the first year in which VIIRS starts operating. These data are generated by detecting thermal signatures that are not impeded in the presence of cloud cover or smoke in the atmosphere. We plot the total number of fires over time in Supplementary Fig. A.1: each year, there are two peaks in the months of April through May and October through November, in correspondence with the crop burning seasons of wheat and rice, respectively49.
Crops: we gather data on average crop area, production and yield in 2010—the latest year in which micro-level crop data are available before the start of our analysis sample (2012–2022). These data are retrieved from Food and Agriculture Organization of United Nations (FAO) maps (https://mapspam.info/index.php/methodology/). They have a resolution of 5 × 5 arcmin, corresponding to approximately the size of a grid cell in our data.
Winds: the data on wind direction and speed are from the European Centre for Medium-Range Weather Forecasts (ECMWF; https://www.ecmwf.int/en/forecasts/datasets) and are precise at a 25 × 25 km resolution. We proxy for the expected wind direction of bureaucrats in a given month–year by computing ten-year rolling averages of wind direction in each month and year by using wind data starting in the 2000s. We plot actual and rolling averages of wind directions in Supplementary Fig. A.1 (right). We can see that rolling averages closely follow the actual wind direction. The two major wind switches happen on average in June and August each year.
Administrative borders: we use the latest 2017 census boundaries for Pakistan, which are available on the Humanitarian Data Exchange website50. For the Indian administrative boundaries, we rely on the latest version generated by Community Created Maps of India51. These maps compile information from various government websites to create district-level polygons. During our analysis period, some districts underwent splits: Charki Dadri was carved out from Bhiwani in Haryana on 1 December 2016, and Malerkotla district was bifurcated from Sangrur in Punjab on 14 May 2021. We account for these splits by modifying our shapefiles accordingly. Our main sample of analyses comprises states in which the problem of crop stubble burning is particularly salient. These states are Punjab in Pakistan, and Punjab, Haryana, Uttar Pradesh, Delhi and Bihar in India. These are the same states considered in other papers focusing on air pollution in South Asia5.
Combining all the data at a 5 km2 grid-cell level: we combine the above datasets into a single gridded dataset of 5 km2. Wind and fire grids do not overlap precisely. To reduce errors in the assignment of fires to wind grids, we split the 25 × 25 km wind grid cells into squares of 5 km2 (around 2.25 km × 2.25 km size). This guarantees that, even if wind and fire grids do not perfectly overlap, data on fires are assigned to the wind grid cell that contains most of the fire grid cell. Splitting wind grids in smaller parts also enables us to address a second measurement concern—some grid cells might overlap across jurisdiction boundaries. By reducing their size, we minimize the amount of misattribution and we are able to assign each grid cell to the appropriate jurisdiction. Notice that any remaining error should bias our findings towards zero, as errors are uncorrelated with the dependent and independent variables. In total, 11% of the grid cells fall between more than one district. These are assigned to the district that contains a larger overlapping area with them (and we show that our results are robust to dropping these from the analysis). For crop data, we assign the pre-period crop production to each 5 km2 grid cell. This enables us to classify what crops are produced in each grid-cell location. Figure 2a illustrates our final gridded dataset for Hardoi District in Uttar Pradesh, where each square refers to a 5 km2 grid cell. Within the sample area, we are able to observe 149,830 grid cells in each month and year from September 2012 to August 2022, for a total of 17,979,600 observations. Supplementary Table A.1 presents descriptive statistics on these data.
Air pollution: we use data on air pollution from the Copernicus Atmosphere Monitoring Service (CAMS; https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global-reanalysis-eac4?tab=overview), measuring the concentration in the air of particulate matter with a diameter of 2.5 μm or less (PM2.5). The data are available at a resolution of approximately 80 km. In Supplementary Fig. A.2, we map the total concentration of PM2.5 in the air in 2014. The north of the subcontinent is covered by a dense cloud of air pollution (darker red), unlike any other parts of the globe.
PM2.5 data from CAMS are valuable as they offer high-resolution information on air quality in all the regions of our sample. However, these data are estimates obtained by combining satellite observations, ground-based measurements and atmospheric models. To validate these estimates, we gather data that are purely ground based: measurement of air pollution from air monitors52. The database contains data from 2,110 monitors operated by the Central Pollution Control Board (CPCB) in the period 2015–2019. These data cover 76 CAMS grids out of the 2,156 in our entire study area, and not in all months, making it too sparse for use as a regressor. We instead use air monitors to validate CAMS estimates in Supplementary Fig. E.2, which shows a positive correlation (0.726) between the total pollution recorded each month by the monitoring stations and CAMS estimates.
Bureaucratic punishment: for India, penalty data come from the Socioeconomic High-resolution Rural-Urban Geographic Platform for India (SHRUG) database (https://www.devdatalab.org/judicial-data), which contains all judicial case records from the Indian eCourts platform. From this database, we extracted all judicial cases in which there was a violation of the Air (Prevention and Control of Pollution) Act of 1981, which includes crop burning offences from 2010–2018, the years during which the data overlap with our sample. Measurement error in these data should bias our results downwards. In the case of Pakistan, we obtained new data on police records of criminal activity from the Central Police Office in Lahore, called First Information Reports (FIRs), that relate specifically to crop burning offences for the years 2017–2022.
Bureaucratic control of fires
Treatment definition
As discussed above, the treatment turns on when smoke from a potential fire in a grid cell would pollute the home district of the administrator. For each month–year, we proxy for administrators’ expectations on actual wind directions by using 10 month–year rolling average of wind direction. In practice, these beliefs are formed through long-term local knowledge, as well as through weather data that are made available as part of bureaucrats’ routine work53.
Once we have this measure of wind direction, we use it to calculate the area of the jurisdiction that is downwind, as well as upwind, from each grid cell in each month. We do this by sketching a 180° line that passes through the grid-cell centroid and is perpendicular to the wind direction vector such that the district is split into a downwind area and an upwind area54. Using this definition, a grid cell is treated in a month–year if the area of the jurisdiction downwind from the grid cell is larger than the area upwind. Figure 2b provides a visual illustration of this treatment definition. Figure 2c illustrates the actual representation of treatment status variation in different months.
Estimation
$${\rm{Fire}}{{\rm{s}}}_{it}=\mathop{\sum }\limits_{r=-s}^{z}\zeta {D}_{ir}+{\alpha }_{i}+{\gamma }_{j}\times {\theta }_{t}+\delta {{\bf{X}}}_{it}+{\nu }_{it}$$
(1)
where our dependent variable Firesit is the number of fires in a grid cell and month. The coefficient of interest ζ is the effect of treatment in each period before and after treatment, relative to the switching month–year (the omitted category). For example, if changes in wind direction convert grid-cell i from control to treatment in May, relative time will take value 0 in May, +1 in June, −1 in April, and so on. We restrict attention to observations in a six-month window around treatment, where the vast majority of our sample lies (Supplementary Fig. A.3) and drop treated observations in the first period in the sample, September 2012, as their relative time cannot be established.
The high dimensionality of our data allows for the inclusion of fine-grained fixed effects, although as we show in the Supplementary Information, the treatment effect size is quite stable when these are included. First, αi are 5 km2 grid-cell fixed effects that absorb time-invariant factors for over 149,000 grid cells in our data, such as crop suitability or economic characteristics of small locations. Second, γj × θt are jurisdiction by month–year fixed effects that account for time-varying shocks at the district level, such as changes in administrative personnel. X is a vector of controls that includes average wind speed and direction. Standard errors are clustered at the grid cell and month–year by district levels to account for temporal dependence of observations as well as spatial correlation in wind.
In addition to the event study design, we also run a difference in differences strategy. For a 5 km2 grid-cell i in jurisdiction (district) j and month–year t, we estimate:
$${{\rm{Fires}}}_{it}=\beta {{\rm{Downwind}}}_{it}+{\alpha }_{i}+{\gamma }_{j}\times {\theta }_{t}+\delta {{\bf{X}}}_{it}+{{\epsilon }}_{it}$$
(2)
The coefficient of interest β captures the overall effect of treatment on the number of fires in that grid cell. Downwind represents the treatment variable: a grid switches to treatment when, depending on wind direction, it would expose most of the home district to smoke. Fixed effects, controls and clustering are as in the event study specification described above.
Action and inaction in the control of fires
We build on the above specification to separate the effect of action and inaction towards the control of fires. For each grid cell, we calculate the distance from its respective borders and code whether these borders are upwind or downwind in a given month. This yields a grid-border–grid-cell pair database with spatial variation in the distance of each grid cell to the border and monthly variation in whether pollution exits or enters the border corresponding to each grid cell.
The key advantage of this strategy is that we can separately estimate the effort bureaucrats exert to prevent fires in areas close to upwind borders from the effort not exerted to prevent fires in areas close to downwind borders. Calling a border segment b, we estimate, for each segment b corresponding to grid-cell i in month–year t:
$$\begin{array}{l}{{\rm{F}}{\rm{i}}{\rm{r}}{\rm{e}}{\rm{s}}}_{ibt}=\mathop{\sum }\limits_{q=1}^{5}\kappa {{\rm{D}}{\rm{i}}{\rm{s}}{\rm{t}}{\rm{a}}{\rm{n}}{\rm{c}}{\rm{e}}}_{it}+\mathop{\sum }\limits_{q=1}^{5}\lambda {\rm{D}}{\rm{i}}{\rm{s}}{\rm{t}}{\rm{a}}{\rm{n}}{\rm{c}}{\rm{e}}\\ \,\,\,\,\,\,\,\,\,\,\,\times \,{{\rm{D}}{\rm{o}}{\rm{w}}{\rm{n}}{\rm{w}}{\rm{i}}{\rm{n}}{\rm{d}}{\rm{N}}{\rm{e}}{\rm{i}}{\rm{g}}{\rm{h}}{\rm{b}}{\rm{o}}{\rm{u}}{\rm{r}}}_{ibt}+{\mu }_{i}+{\eta }_{b}\times {\theta }_{t}+{\nu }_{ibt}\end{array}$$
(3)
where q are quintiles of the distance of each grid cell to the bordering district segment b, and we control for grid cell and for border times month–year fixed effects. We use the fifth quintile (the farthest) as the absorbed category. We are interested in two coefficients: κ captures the effect of proximity to the upwind border, an area where bureaucrats have the largest incentives to curb fires, as those would affect their own district the most. Instead, λ captures the differential effect of what happens as we approach the downwind border, which might be the same location as an upwind area in other periods. Our expectation is that bureaucrats have little incentive to act on fires when those are located close to downwind borders, as smoke would exit their district. The logic of this estimation is further illustrated in Fig. 3a, where κ coefficients are represented in blue and λ in red.
Bureaucratic action and farmer deterrence
We code a new district-level variable called PostPunishment, taking value one if a farmer was punished in a district in the previous 12 months. To stay consistent with our calendar month analysis, we start the 12-month duration from the filing month itself if the penalty was filed before the 15th day of the month. Otherwise, it starts the following month. We use this variable in a triple difference-in-differences event study design by interacting it with the down > up dummy referenced as Dir in equation (1). We show the event study results in Fig. 3c. We also plot these results in 1-month windows and find consistent, but noisier, results (Supplementary Fig. D.1). Finally, we show that these results are robust using a simple triple difference-in-differences setup in Supplementary Table D.2.
Mortality impacts of crop burning
Linking fires to downstream recipients of pollution
We measure infant and child mortality using geolocalized individual-level health information from surveys collected by the Demographic and Health Surveys (DHS) program (https://dhsprogram.com/data/) in India (2015–2016 and 2019–2021) and Pakistan (2016–2017). The DHS data enable us to match each pregnancy to a specific small location, called a DHS cluster, across India and Pakistan. In total, we have 542,150 pregnancies across 58,103 DHS clusters for which we observed child, mother and family characteristics.
For each birth, we calculate the pollution exposure from crop burning fires as the total number of particles coming from fires started in our analysis sample (Punjab in India and Pakistan, Haryana, Delhi, Uttar Pradesh and Bihar) to the home DHS cluster location of the birth, for the expected length of the pregnancy. This process enables us to link pollution particles from each upwind fire to pregnancies, thereby generating a precise measure of the expected air pollution from crop burning in our sample for each pregnant mother. Besides the number of particles received in the home DHS cluster location of the pregnancy, we also calculate the total PM2.5 pollution experienced in that location using satellite data from the CAMS.
We proceed to collapse the data at the pregnancy level, calculating the total local air pollution, as well as total particles from crop burning fires received by each child in utero. Children born in the same location and from a similar family will experience different levels of in utero pollution from PM2.5 depending on the timing of the pregnancy.
Estimation strategy
Our quantity of interest is the impact of air pollution (PM2.5 concentration) caused by crop burning on child and infant mortality. This quantity is generally confounded by characteristics of the location, period of birth and of the family, which might be correlated with both air pollution and worse health outcomes. We can address these concerns in our setting using time variation in exposure to pollution. We consider two models. First, we use a reduced form regression of the impact of total crop burning fire particles on mortality outcomes. Second, we capture the impact that fires have on overall levels of air pollution using an instrumental variable strategy. For this, we first estimate the impact of total fire particles on PM2.5 concentration experienced by each child in utero (first stage). We then estimate the impact of the predicted air pollution from fires on children’s mortality (second stage). For each child (c) born in district d during month–year t, we estimate:
$$(\text{first stage})\,{\rm{PM}}{2.5}_{ct}=\pi {{\rm{FireParticles}}}_{ct}+{\xi }_{d}+{\tau }_{t}+{\phi }_{b}+\sigma {X}_{t,dhs}+{\upsilon }_{ct}$$
(4)
$$(\text{second stage})\,{{\rm{M}}{\rm{o}}{\rm{r}}{\rm{t}}{\rm{a}}{\rm{l}}{\rm{i}}{\rm{t}}{\rm{y}}}_{ct}=\rho \hat{{\rm{P}}{\rm{M}}{2.5}_{ct}}+{\xi }_{d}+{\tau }_{t}+{\phi }_{b}+\sigma {X}_{t,dhs}+{\omega }_{ct}$$
(5)
Both the first and second stage control for district ξ, month–year from birth τ, and birth order ϕ fixed effects. All regressions also control for the total number of fires during the pregnancy within a range of 10 km from the DHS cell, X—while most of the particles affecting a pregnancy will be carried from different and farther locations, fires originating around the DHS cluster are likely to have a larger impact on exposure, and locations with more fires might be different on average than other places. Our coefficient of interest, ρ, identifies the effect of exposure to fire-related air pollution in utero on child and infant mortality. We provide evidence consistent with the validity of the IV identifying assumptions in Supplementary Table E.1 (c), which shows a strong first stage of fire particles on PM2.5 concentration, and Supplementary Fig. E.3, which shows that their relationship is monotonic. The extent to which a pregnancy is affected by pollution particles from a specific fire is likely as if random, because this is determined by a host of upwind meteorological factors interacted with the timing of the pregnancy. Finally, the case for exclusion restriction in our setting is strong because it is unlikely that particles from fires have an impact on mortality from channels other than air pollution.
Linking bureaucratic incentives with mortality outcomes
In the main section, we have shown that (1) bureaucratic strategic behaviour impacts the incidence of fires, and that (2) fires have a substantial impact on child and infant mortality. We then connect these two findings to assess how much of the increase in mortality determined by fires can be attributed to jurisdictional incentives. We do a back-of-the-envelope calculation to understand what would happen if bureaucrats had the same incentives and capacity across the entire district area, and exerted control in the downwind parts of their districts the same way they do in the upwind parts. Here we describe the details of our calculation.
We have shown in Supplementary Table B.1 that treatment reduces fires by 10–13% on average in the upwind part of the district, or by approximately 5–7.5% in the full district. We have also estimated that a 1-log increase in particles from fires over the 9 months of the pregnancy raises child mortality by up to 36 deaths in 1,000 births (Supplementary Table E.1). Now, consider what would happen if exposure to fires in utero dropped by the extent that bureaucrats can affect fires in the district. If bureaucrats, assuming proper incentives and capacity, reduced crop burning in the downwind parts of their district to match the upwind parts, child mortality would reduce by between 1.8 to 2.7 fewer deaths in 1,000, which corresponds to a 4.4 to 6.6% reduction over the average child mortality in the sample. Of course, a problem with bureaucratic action is the lack of capacity, but this scenario illustrates one way in which effects could be obtained by creating slack in bureaucratic capacity through the prioritization of environmental management.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.