Particulate matter and COVID-19 excess deaths: Decomposing long-term exposure and short-term effects.

We investigate the time-varying effect of particulate matter (PM) on COVID-19 deaths in Italian municipalities. We find that the lagged moving averages of PM2.5 and PM10 are significantly related to higher excess deceases during the first wave of the disease, after controlling, among other factors, for time-varying mobility, regional and municipality fixed effects, the nonlinear contagion trend, and lockdown effects. Our findings are confirmed after accounting for potential endogeneity, heterogeneous pandemic dynamics, and spatial correlation through pooled and fixed-effect instrumental variable estimates using municipal and provincial data. In addition, we decompose the overall PM effect and find that both pre-COVID long-term exposure and short-term variation during the pandemic matter. In terms of magnitude, we observe that a 1 μg/m3 increase in PM2.5 can lead to up to 20% more deaths in Italian municipalities, which is equivalent to a 5.9% increase in mortality rate.

JEL numbers

I18

Q53

J18

H12

Introduction

The impact of the COVID-19 pandemic in relation to contagion and deaths in the first half of 2020 is markedly heterogeneous from a geographical perspective. Several studies have tried to identify the causal factors of this puzzling outcome. Epidemiological literature has identified the frequency of human physical interactions as a leading causal factor of contagions. However, even after controlling for these factors, a significant part of the observed variability of COVID-19-related outcomes remains unexplained. This paper aims to shed light on this issue by investigating the role of particulate matter (PM) in the pandemic's high mortality rate.

The theoretical background for our research hypothesis can be summarized into two main literature strands. The first strand deems long-term exposure to PM as a contributing factor to COVID-19-related deaths. The research hypothesis relies on the maintained assumption that PM inhalation induces inflammation and oxidative stress, thereby reducing lung efficiency and contributing to respiratory and pulmonary diseases (see Pope and Dockery, 2006). Over the years, several empirical papers have estimated the relationship between long-term exposure to PM and total mortality, and mortality from cardiovascular and respiratory diseases (Kim et al., 2015; Pelucchi et al., 2009; Pinault et al., 2017; Faustini et al., 2011; Anderson, 2020; Ciencewicki and Jaspers, 2007; Sedlmaier et al., 2009) as well as the effect of fine PM as a factor in cardiovascular and respiratory morbidity and mortality (McGuinn et al., 2017; Jeong et al., 2017; Yin et al., 2017; Cakmak et al., 2018).

Based on this literature strand, several researchers have tried to ascertain if the effect of COVID-19 on respiratory and pulmonary diseases can be enhanced by PM exposure. Wu et al. (2020) found that an increase in exposure to PM2.5 is associated with increased COVID-19 fatality in the US, Cole et al. (2020) find similar results for the Netherlands.1 Focusing on Italy, Cartenì et al. (2020) find that the number of days in 2019 in which the national PM10 exceeded the 50 μg/m3 daily limit is positively correlated with the number of certified daily cases. Perone (2020) finds that the case fatality rate is affected by ozone and nitrogen dioxide beyond PM, while Becchetti et al. (2022) used mortality data at province level to confirm this relationship controlling for several concurring factors. Coker et al. (2020) found support for the relationship in Northern Italy, estimating that a 1 μg/m3 increase in PM2.5 is associated with a 9% increase in COVID-19-related mortality. Many other studies have been conducted using different geographical samples to confirm the positive relationship between long-term exposure to PM and increase in COVID-19-related deaths (Ogen, 2020; Magazzino et al., 2020; Yongjian et al., 2020; Travaglio et al., 2021; Setti et al., 2020; Comunian et al., 2020).

The second strand of literature developed a theoretical hypothesis on the relevance of short term effects and on the role of PM as a carrier of the SARS-CoV-2 virus. Preliminary findings in this direction are found by Setti et al. (2020). The authors show the presence of SARS-CoV-2 viral RNA by detecting highly specific RtDR gene on eight filters in two parallel PCR analyses on 34 PM10 samples of outdoor/airborne PM10 in Bergamo province. However, they state that it is impossible to assess the viral charge of the carried virus outside the human body. Following this hypothesis, Delnevo et al. (2020) find that the lagged PM Granger causes adverse COVID-19 outcomes in several Italian provinces located in the Emilia-Romagna region. Similarly, Zoran et al. (2020) find a correlation between daily average ground levels of particulate matter concentrations and new cases and deaths in Milan. Isphording and Pestel (2020) conduct the same analysis for German regions. In the same vein, Austin et al. (2020) find that contemporary variation in PM significantly affects COVID-19 contagions and deaths in US counties, while Becchetti et al. (2020) use daily atmospheric data from the European regional level provided by the Copernicus Atmosphere Monitoring Service (CAMS) and find that PM2.5 and PM10 concentration positively affects confirmed cases and deaths. They estimated that the effect peaks at the 6th to the 8th day lag for confirmed cases and the 13th day lag for deaths. In the opposite direction, Bontempi (2020) find, after assessing data from Lombardia and Piemonte, that it is impossible to conclude that COVID-19 diffusion also occurs through the air using PM10 as a carrier.

Our contribution is original with respect to the existing literature in several aspects. First, we disentangle the long-term exposure and the short-term effects to test the two aforementioned research hypotheses on the effect of PM on COVID-19 related outcomes in Italy simultaneously. Second, we control for heterogeneous pandemic dynamics and spatial correlation providing empirical evidence at both the municipality and province levels. Third, we control for the differential introduction of lockdown measures adopted by the Italian government. Fourth, we net out the effects of other unobservable time-invariant local confounders (i.e., municipal policies and regional health systems) at the finest and more disaggregated geographical level through municipality fixed effects. Last, we use instrumental variable approaches to mitigate endogeneity concerns.

Our empirical findings show a positive and significant relationship between particulate matter and excess deaths in Italian municipalities and provinces during the first pandemic wave. Our results are robust and confirmed when using instrumental variables and when controlling for heterogeneous epidemics dynamics and spatial correlation. In terms of economic significance, we find that if conclusions from our IV estimates pointing at causality hold, a 1 μg/m3 increase in PM2.5 causes a 10–20% surge in excess deaths in Italian municipalities, depending on the model used. This is equivalent to an overall 3–6% increase in mortality rate.

Data

Our first data source is the Italian National Statistical Institute (ISTAT), which provides information on daily deaths in each municipality. We use the difference between daily deaths during the pandemic and the 2015–2019 five-year average of the corresponding days as the dependent variable (see Table A1 in the Appendix). This measure overcomes two well-known problems that arise when using official COVID-19 registered deaths. First, it is not always possible to ascertain whether victims died because of COVID-19 or with COVID-19, with Italy's independent local health systems, different municipalities, and regions that interpret this distinction differently. The resulting heterogeneity in death registration, therefore, creates an implicit measurement error. Second, at the peak of the pandemic in Italy—March and April 2020—hospitals were overcrowded and several COVID-19 victims could not access hospital care and died without a proper diagnosis.

Table A1

Variable legend.

Variable	Description	Source
PM10	11-day (from t−10 to t) moving average of particulate matter with diameter < 10 μm (μg/m3)	Copernicus Atmospheric Monitoring Service (CAMS) -
PM2.5	11-day (from t−10 to t) moving average of particulate matter with diameter < 2.5 μm (μg/m3)	Copernicus Atmospheric Monitoring Service (CAMS) -
Excess Deaths	Daily difference of total deaths in 2020 and the 2015–19 average total deaths at municipality level	Italian National Statistical Institute
Mobility	Number of people in transit in subway, bus, train stations, sea port, taxi stand, highway rest stop and car rental agencies in the given Italian province (Change compared to the baseline of the median value, for the corresponding day of the week, during the previous 5-week period).	Google: Community Mobility Report
Population	Number of residents in 2011 at municipality level per 1000 inhabitants.	Italian National Statistical Institute
Employees	Number of employees operating in all economic sectors at municipality level per 1000 inhabitants.	Italian National Statistical Institute
Employees in Essential Sectors	Number of employees operating in essential economic sectors (as defined by the Decree of the Italian President of the Municipality of Ministers, released on March 22nd and revised on March 25th), at municipality level per 1000 inhabitants.	Italian National Statistical Institute
Density	Population per municipality area per 1000 inhabitants.	Italian National Statistical Institute
Over 65	Share of people aged 65 or above per 1000 inhabitants.	Italian National Statistical Institute
Income	Total municipality gross income (billion euros)	Italian National Statistical Institute
Rain	11-day (from t−10 to t) moving average of total precipitation in mm at municipality level	Copernicus Climate Change Service (CCCS)
Temperature	11-day (from t−10 to t) moving average of air temperature measure at the height of 2 m above ground, at municipality level.	Copernicus Climate Change Service (CCCS)
Region	Italian regions.
Days since lockdown	Days since the start of the national lockdown considering the different starting days based on subsequent government decisions (see footnote 7).

The second and third data sources are the Copernicus Atmospheric Monitoring Service (CAMS) and the Copernicus Climate Change Service (C3S), which provide data on air quality data and weather conditions.2 The analysis conducted in this paper exploits the following C3S datasets: (i) the C3S ERA5-Land hourly data from 1981 to present, which contains measures on various land variables over several decades at a global scale3 ; (ii) the CAMS European Air Quality Forecast, which provides information on air quality in Europe.4

We use the following variables:

[C3S – ERA5-Land] 2 m temperature [K]: Air temperature at 2 m (height) obtained by interpolating the lowest model level and the Earth's surface; we converted temperature values to Celsius degrees by subtracting 273.15.

[C3S – ERA5-Land] Total precipitation [m]: Accumulated precipitation in millimeters, including rain and snow, that falls to the Earth's surface considering possible steps during a single day; we converted precipitation values to millimeters by dividing by 1000.

[CAMS – European Air Quality] Surface (average individual's height level) particulate matter d < 2.5 μm [μg/m3]: Fine solid or liquid particles in the atmosphere emitted by natural and anthropogenic sources with a diameter less than 2.5 μm (PM2.5);

[CAMS – European Air Quality] Surface (average individual's height level) particulate matter d < 10 μm [μg/m3]: Fine solid or liquid particles in the atmosphere emitted by natural and anthropogenic sources with a diameter less than 10 μm (PM10).

The grid cover consists of points where the information is recorded spanning 0.1 degrees in latitude and longitude, i.e., the grid was made by the four points of the vertices. We use the Python-3 high-level programming language to download and process Copernicus data to extract the final dataset at the municipal level. Moreover, hour-specific filters were applied to get daily mean values averaging data available at 4-step hours during the day (8:00, 12:00, 20:00, 00:00) for the variables taken into account.

We identify the municipal polygon's centroid for each municipal area and define each municipal through its centroid. Since both weather and pollution variables are available according to a regular 10-km square grid, each centroid is associated with a variable's value based on minimum distance. By using this procedure, we obtain an average mean distance of less than 4 km. In many cases, the distance between the centroid and the closest grid is less than 1 km, which represents a sharp characterization of the real observations in the proximity of the spatial coordinates. In the section that follows, we explain how we use inverse distance weights to account for these differences.

Fig. 1 A–C provide clear descriptive evidence of the geographical distribution of PM2.5 and PM10 concentrations in relation to the geographical distribution of the high death rates. These figures show that the northern macroregion of Pianura Padana was the most affected area due to its higher economic activity and peculiar geographical conformation, i.e. a large plain surrounded by high mountains, where air tends to stagnate more than in other areas of the country. In general, the PM2.5 map shows that a large part of the country has dark brown areas above the World Health Organization's threshold (10 μg/m3).

Fig. 1

A–C. Excess deaths, PM2.5 and PM10 in Italian municipalities.

Note: Excess deaths is the daily difference of total deaths in 2020 and the 2015–19 average total deaths at municipality level (Source: Istat); PM2.5 is the 11-day (from t−10 to t) moving average of particulate matter with diameter < 2.5 μm (μg/m3); PM10 is the 11-day (from t−10 to t) moving average of particulate matter with diameter < 2.5 μm (μg/m3).

Econometric model

To test the impact of particulate matter controlling for potential concurring factors, we estimate the following equation:where our dependent variable (Excess Deaths ) is the difference between total deaths in 2020 in municipality m on day t and the 2015–19 total average deaths in the corresponding municipality and day of the year. The main independent variable of interest is Pollution(MA), calculated as a moving average from day t - 10 to day t of PM10 or PM2.5, measured in municipality m on day t.5 We introduce linear, quadratic, and cubic time trends (t, t , t ) starting with the disease outbreak, which is conventionally fixed as February 24, 2020 (the beginning of our sample period) among control variables. These trends capture part of the deterministic evolution of the pandemic consistently with standard epidemiological modeling approaches (further robustness checks for heterogeneous pandemic dynamics are presented and discussed in section 6).6 Among other controls, Days_Since_Lockdown counts the days since the national lockdown, taking into account the three government decisions that progressively introduced mobility restrictions in Italian municipalities.7 Population is the number of residents in municipality m from the last Italian census (2011) (per 1000 inhabitants); Density is the population density in municipality m (per 1000 inhabitants); Over65 is the proportion of people aged 65 or above and living in municipality m (per 1000 inhabitants); Income is the total before-tax income in municipality m (in billion euros); Employees and Essential-Employees are the number of employees operating in all sectors and in essential sectors only (per 1000 inhabitants) at the municipal level. The essential sectors are those on a list of activities that the Italian government allowed to operate during the lockdown.8 These last two variables capture lockdown-induced local differences in job commuting flows due to the different incidences of essential and non-essential sectors in each municipality. Temperature(MA) is the 11-day (from t− 10 to t) moving average of daily air temperature in each municipality. Last, we control for time-varying human interactions with a variable (Mobility) measuring transit in the subways, bus and train stations, seaorts, taxi stands, highway rest stops, and car rental agencies in Italian provinces. The variable is calculated in first differences, that is, as a change in the number of people in the above-mentioned transit areas compared to the baseline of the median value, for the corresponding day of the week, during the previous 5-week period. We also add region dummies (DRegion) to control for time-invariant features of the regions, such as urbanization rate or health system characteristics. In fact, health policies in Italy are run at the regional level, thereby making health capital endowments highly heterogeneous across regions.9 Standard errors are clustered at the municipal level. A detailed description of variables and their sources is given in Table A1 in the Appendix.

Table 1 presents descriptive statistics of the variables used in our econometric specifications. As expected, the moving averages used in the estimates smoothen extreme values of pollution and atmospheric indicators, with maxima of moving averages of particulate concentration reaching 45.71 and 61.72 μg/m3, respectively. Nonetheless, the mean value of the PM2.5 moving average (14.02 μg/m3) during the sample period is above the average yearly threshold suggested by the World Health Organization (10 μg/m3).10 Our sample period covered the end of winter and spring, and therefore we did not observe extreme hot temperature events (the single-day maximum is 28.86, while the moving average maximum 24.47).

Table 1

Descriptive statistics.

Variable	Obs	Mean	St. Dev.	Min	Max
Excess deaths	713244	0.011	0.252	−14.708	14.925
PM2.5⁎	713244	14.022	7.408	1.862	45.714
PM10⁎	713244	18.679	8.948	2.522	61.715
Rain⁎	713244	2.390	2.449	0.0002	18.331
Days since lockdown	713244	35.307	26.684	0	98
Population	713244	7.715	41.816	0.034	2617.175
Density	713244	0.307	0.659	0.001	12.924
Over 65	713244	1.79	10.70.	0.007	638.523
Income	713244	0.108649	0.783832	0.542792	49314.36
Employees	713244	2.248	17.862	0.001	1023.890
Employees in essential sectors	711480	0.985	9.607	0.001	547.307
Mobility	713244	−3.092	15.927	−79	65
Temperature(MA)⁎	713244	11.376	4.862	−11.15	24.47014

Contains modified Copernicus Climate Change Service Information [2017–2020], DOI: 10.24381/cds.e2161bac’. Contains modified Copernicus Atmospheric Monitoring Service Information [2017–2020].

Results

Before running our estimates, we perform panel stationarity tests and find that all our series are stationary. More specifically, we perform the Lein-Lin-Chu (2020) test for unit roots in panel datasets and find that the null of non-stationarity is rejected in all cases (p < 0.001 for all series)11 .

In Table 2 , we present the results of our main econometric specification. Columns 1 and 2 display (unweighted) pooled OLS estimates of the effects of PM2.5 and PM10, respectively. In columns 3 and 4, observations are weighted for the inverse of the distance from the centroid to give more importance to municipality centroids that lie closer to the geographical point of our meteorological observation of PM.

Table 2

Pooled OLS estimates.

Variables	(1)	(2)	(3)	(4)
PM2.5	0.000966***		0.00113***
	(6.86e-05)		(9.95e-05)
PM10		0.000590***		0.000660***
		(3.74e-05)		(5.59e-05)
T (linear day trend)	0.00282***	0.00287***	0.00286***	0.00294***
	(0.000129)	(0.000129)	(0.000185)	(0.000186)
T2 (quadratic day trend)	−7.88e-05***	−8.58e-05***	−8.10e-05***	−9.05e-05***
	(5.35e-06)	(5.14e-06)	(1.10e-05)	(1.12e-05)
T3 (Cubic day trend)	4.80e-07***	5.11e-07***	4.87e-07***	5.30e-07***
	(3.11e-08)	(3.01e-08)	(5.84e-08)	(5.91e-08)
Days since lockdown	0.000686**	0.00100***	0.000830	0.00123*
	(0.000281)	(0.000272)	(0.000645)	(0.000652)
Population	−8.03e-05	−7.14e-05	−0.000135*	−0.000122*
	(5.58e-05)	(4.99e-05)	(7.27e-05)	(6.37e-05)
Density	−0.00274***	−0.00236***	−0.00300***	−0.00249***
	(0.000457)	(0.000429)	(0.000560)	(0.000527)
Over 65	0.000978***	0.000875***	0.00134***	0.00118***
	(0.000357)	(0.000313)	(0.000484)	(0.000418)
Income	−0.0136***	−0.0123***	−0.0180***	−0.0159***
	(0.00329)	(0.00298)	(0.00510)	(0.00456)
Employees	−0.000453**	−0.000382**	−0.000460	−0.000323
	(0.000206)	(0.000192)	(0.000311)	(0.000295)
Employees in Essential Sectors	0.00130***	0.00112***	0.00152**	0.00120**
	(0.000395)	(0.000367)	(0.000614)	(0.000578)
Temperature	−0.000333**	−0.000154	−0.000263	−1.02e-05
	(0.000162)	(0.000163)	(0.000337)	(0.000342)
Mobility	7.99e-05***	6.80e-05***	7.16e-05	5.84e-05
	(2.46e-05)	(2.45e-05)	(6.21e-05)	(6.20e-05)
Region dummies	Yes	Yes	Yes	Yes
Constant	−0.0218***	−0.0184***	−0.0223***	−0.0181***
	(0.00228)	(0.00214)	(0.00363)	(0.00349)
Observations	685,451	685,451	685,451	685,451
Log Likelihood	36,835	36,827	6369	6339

Note: Columns (1) and (2) do not weight observations, while columns (3) and (4) use as weight the inverse distance of municipality centroids from the meteorological point of observation. Standard errors clustered at municipality level in parentheses. Contains modified Copernicus Climate Change Service Information [2017–2020], DOI: 10.24381/cds.e2161bac. Contains modified Copernicus Atmospheric Monitoring Service Information [2017–2020]; *** p < 0.01, ** p < 0.05, * p < 0.1.

Our empirical findings show that the high mortality in 2020 is significantly and positively related to both air pollution measures. In terms of magnitude, the effect of PM2.5 is larger than that of PM10, with results from weighted and unweighted estimates for the same pollution variable being quite similar. The estimated pollution effect in column 3 implies that 1 μg/m3 of additional PM2.5 concentration creates an approximately 10% increase in the average value of the dependent variable, that is 0.113 extra deaths per day per 100,000 inhabitants, which corresponds to a 3.32% increase in mortality rate. The total effect over the 94 days of the pandemic considered in our sample is 1.07 extra excess deaths per 100,000 inhabitants. This implies that the effect over the entire Italian population is about 647.96 extra deaths per μg/m3. Based on our coefficient magnitude, we estimate that a difference of about 19 μg/m3between the municipalities with the highest and lowest PM2.5 average concentration in the sample would generate a difference of 1231.13 more deaths in the overall sample period.

The linear, quadratic, and cubic trends are strongly significant among the control variables and with the expected sign, displaying non-linear pandemic dynamics during the first phase. The share of employees in essential sectors is positive and significant and likely to capture the positive effect of high death rates on economic activity in industries that could not stop their operations during the lockdown. Time-varying mobility is, as expected, positive and significant as an increase in people in transit stations has a positive and significant effect on excess deaths. The negative sign of the density variable can be explained by the fact that, with population among the regressors, the variable captures the positive effect of municipality surface on excess deaths, likely to be explained by how far inhabitants are from institutions and less accessible health services. The positive and significant effect on excess deaths of the share of the elder population at the municipal level is also expected.

We implement an instrumental variable approach from omitted variables and reverse causality to mitigate a possible estimation bias deriving from measurement error in the dependent variable (Table 3 ). We instrument the PM moving averages in Eq. (1) with the four-day lagged corresponding 11-day moving average of daily rainfalls controlling in our estimates for the mobility variable. We can confidently argue that the chosen instrument is relevant since rainfalls have a strongly significant and negative effect on PM concentrations in the first-stage estimation.12 The exclusion restriction is also likely to be satisfied in our case since—apart from its direct effects on pollution—it is implausible that four-day lagged rainfall moving averages significantly affect the difference in deaths between 2020 and the previous years on a given day. Rainfall may discourage mobility and reduce contagion or increase car vs. public-transport mobility (again reducing contagion), thereby potentially invalidating the exclusion restriction. However, these potential threats to the exclusion restriction can be excluded since we condition for mobility in all estimates. Furthermore, most of the mobility decisions made during the lockdown that cover most of our sample period were forced, with little impact on atmospheric conditions. The pairwise correlation between rain and mobility during lockdown is 0.03, which supports our hypothesis. This is a positive (yet low in magnitude) and non-statistically significant correlation, which goes against the prediction of a potential negative association between the two variables. Furthermore, the instrument is not statistically significant if we introduce it in our baseline non-instrumented specification (Eq. (1)), further supporting its validity hypothesis. Our main findings remain unchanged and the coefficient magnitude is remarkably close to the non-instrumented estimates.

Table 3

Pooled IV estimates.

Variables	(1)	(2)
PM2.5	0.00117***
	(0.000297)
PM10		0.000522***
		(0.000133)
T (linear day trend)	0.00303***	0.00329***
	(0.000227)	(0.000189)
T2 (quadratic day trend)	−6.52e-05***	−8.14e-05***
	(1.01e-05)	(6.49e-06)
T3 (Cubic day trend)	3.98e-07***	4.78e-07***
	(5.41e-08)	(3.68e-08)
Days since lockdown	−0.000155	0.000398
	(0.000386)	(0.000291)
Population	−8.73e-05***	−7.46e-05***
	(1.95e-05)	(1.89e-05)
Density	−0.00291***	−0.00227***
	(0.000392)	(0.000281)
Over 65	0.00104***	0.000878***
	(0.000121)	(0.000102)
Income	−0.0143***	−0.0121***
	(0.00149)	(0.00113)
Employees	−0.000497***	−0.000369***
	(0.000121)	(0.000105)
Employees in Essential Sectors	0.00140***	0.00110***
	(0.000246)	(0.000203)
Temperature	−0.000512*	−6.94e-05
	(0.000304)	(0.000213)
Mobility	1.01e-05	2.57e-05
	(2.79e-05)	(2.66e-05)
Region dummies	Yes	Yes
Constant	−0.0289***	−0.0244***
	(0.00317)	(0.00268)
Observations	670,865	670,865
Log Likelihood	35,659	35,648

Note: Contains modified Copernicus Climate Change Service Information [2017–2020], DOI: 10.24381/cds.e2161bac. Contains modified Copernicus Atmospheric Monitoring Service Information [2017–2020]; *** p < 0.01, ** p < 0.05, * p < 0.1.

In Table 4 , we re-estimate our benchmark specification through OLS panel fixed effects. This allows for capturing unobservable time invariant idiosyncratic factors at the finest geographical unit, e.g. the quality of local majors or local health governance at the municipality level. The significance of the PM2.5 and PM10 variables is also confirmed in this model. In Table 5 , we instrument the PM moving average as the fixed-effect estimates with the instrument used in Table 3. Our results are again confirmed.

Table 4

OLS panel fixed-effects estimates.

Variables	(1)	(2)	(3)	(4)
PM2.5	0.00123***		0.00145***
	(7.60e-05)		(0.000144)
PM10		0.000644***		0.000713***
		(3.80e-05)		(5.87e-05)
T (linear day trend)	0.00295***	0.00298***	0.00302***	0.00307***
	(0.000131)	(0.000131)	(0.000195)	(0.000196)
T2 (quadratic day trend)	−5.59e-05***	−7.10e-05***	−5.33e-05***	−7.35e-05***
	(5.52e-06)	(5.28e-06)	(8.85e-06)	(9.41e-06)
T3 (Cubic day trend)	3.54e-07***	4.28e-07***	3.33e-07***	4.34e-07***
	(3.17e-08)	(3.07e-08)	(5.01e-08)	(5.13e-08)
Days Since Lockdown	−0.000670**	0.000134	−0.000794	0.000251
	(0.000303)	(0.000285)	(0.000505)	(0.000535)
Temperature	−0.000409*	−0.000404*	−0.000425	−0.000375
	(0.000219)	(0.000221)	(0.000419)	(0.000427)
Mobility	0.000107***	8.00e-05***	9.97e-05*	6.81e-05
	(2.49e-05)	(2.46e-05)	(5.70e-05)	(5.88e-05)
Constant	−0.0271***	−0.0207***	−0.0317***	−0.0234***
	(0.00247)	(0.00227)	(0.00592)	(0.00503)
Number of municipalities	7260	7260	7260	7260
Observations	685,451	685,451	685,451	685,451
Log Likelihood	41,862	41,837	11,200	11,146

Note: Columns (1) and (2) do not weight observations, while columns (3) and (4) use as weight the inverse distance of municipality centroids from the meteorological point of observation. Standard errors clustered at municipality level in parentheses. Contains modified Copernicus Climate Change Service Information [2017–2020], DOI: 10.24381/cds.e2161bac. Contains modified Copernicus Atmospheric Monitoring Service Information [2017–2020]; *** p < 0.01, ** p < 0.05, * p < 0.1.

Table 5

IV panel fixed-effect estimates.

Variables	(1)	(2)
PM2.5	0.00158***
	(0.000270)
PM10		0.000736***
		(0.000126)
T (linear day trend)	0.00380***	0.00385***
	(0.000194)	(0.000192)
T2 (quadratic day trend)	−3.02e-05***	−5.36e-05***
	(1.06e-05)	(7.15e-06)
T3 (Cubic day trend)	2.10e-07***	3.26e-07***
	(5.67e-08)	(4.01e-08)
Days since lockdown	−0.00282***	−0.00161***
	(0.000563)	(0.000409)
Temperature	−0.000800***	−0.000583**
	(0.000249)	(0.000233)
Mobility	4.84e-05*	4.80e-05*
	(2.63e-05)	(2.63e-05)
Constant	−0.0318	−0.0231
	(3.473e+10)	(3.473e+10)
Wald χ2	1618.25	1618.07
Observations	670,865	670,865
Number of municipalities	7260	7260

Note: Contains modified Copernicus Climate Change Service Information [2017–2020], DOI: 10.24381/cds.e2161bac. Contains modified Copernicus Atmospheric Monitoring Service Information [2017–2020]; robust standard errors in parentheses. *** p < 0.01, ** p < 0.05, * p < 0.1.

Prolonged exposure vs. short term effects: decomposing the total PM effect

The 11-day moving average used so far mainly captures the time-varying effect of PM on excess deaths. However, it is reasonable to assume that this measure is also influenced by a long-term component capturing long-term, pre-COVID exposure to PM. This component is regarded as the crucial factor affecting the negative consequences of COVID-19 infection, according to the first strand of the literature described in the introduction.

To disentangle the effects deriving from these two—long-term structural and short-term time–varying—components, we propose the following decomposition. First, since the time-varying component can be correlated with historical levels of PM concentration, we regress the PM 11-day moving average, that is, Pollution(MA), on (time-invariant) average PM concentration in the two years before the pandemic, i.e., PM(2018–2019). More specifically, we estimate the following model:

Then, we compute the time-varying residuals , which can be interpreted as a “cleaner” measure of the time-varying effect, i.e., the variation Pollution(MA) that is not explained by the variation in the long-term PM component.

We, therefore, run our benchmark model as in Eq. (1) by replacing Pollution(MA) with its time-varying residual component , and the two-year (time-invariant) average of PM concentration. The estimating model reads as:

The results from the OLS pooled estimates of Eq. (3) show that the coefficients of both components (ß1 and ß2) are positive and statistically significant (columns 1 and 2 of Tables 6 −7 for PM2.5 and PM10, respectively). Our interpretation is that long-term exposure and time-varying effect significantly predict excess mortality.

Table 6

Decomposition of the long-term and short-term effects.

Variables	(1)	(2)	(3)	(4)
PM2.5 (short term component)	0.000850***	0.000872***	0.00110***	0.00119***
	(6.79e-05)	(9.71e-05)	(7.65e-05)	(0.000114)
PM2.5 (ex-ante component)	0.000945***	0.00127***
	(0.000177)	(0.000264)
T (linear day trend)	0.00312***	0.00327***	0.00347***	0.00368***
	(0.000157)	(0.000226)	(0.000162)	(0.000237)
T2 (quadratic day trend)	-7.67e-05***	−8.34e-05***	−5.21e-05***	−5.56e-05***
	(4.95e-06)	(1.11e-05)	(5.30e-06)	(8.95e-06)
T3 (cubic day trend)	4.60e-07***	4.91e-07***	3.27e-07***	3.42e-07***
	(2.82e-08)	(5.71e-08)	(2.98e-08)	(4.85e-08)
Population	−8.02e-05	−0.000135*
	(5.33e-05)	(7.12e-05)
Employees	−0.000434**	−0.000430
	(0.000200)	(0.000307)
Density	−0.00252***	−0.00285***
	(0.000451)	(0.000550)
Employees in Essential sectors	0.00124***	0.00144**
	(0.000383)	(0.000605)
Income	−0.0129***	−0.0174***
	(0.00319)	(0.00507)
Over 65	0.000944***	0.00131***
	(0.000341)	(0.000475)
Days since lockdown	0.000318	0.000598	−0.00135***	−0.00131**
	(0.000290)	(0.000709)	(0.000324)	(0.000571)
Temperature	−0.000252	−0.000304	−0.000515**	−0.000492
	(0.000182)	(0.000368)	(0.000219)	(0.000407)
Mobility	3.96e-05	2.86e-05	6.27e-05**	5.33e-05
	(2.49e-05)	(6.56e-05)	(2.50e-05)	(6.28e-05)
Constant	−0.0251***	−0.0275***	−0.0172***	−0.0132***
	(0.00291)	(0.00469)	(0.00251)	(0.00473)
Municipality fixed effects	No	No	Yes	Yes
Observations	685,385	685,385	685,385	685,385
Log Likelihood	37,948	7547	42,993	12,422
Number of Municipalities			7260	7260

Columns (1) and (2) pooled estimates, columns (3) and (4) fixed effect estimates. Standard errors clustered at municipality level in parentheses. Contains modified Copernicus Climate Change Service Information [2017–2020], DOI: 10.24381/cds.e2161bac. Contains modified Copernicus Atmospheric Monitoring Service Information [2017–2020]; *** p < 0.01, ** p < 0.05, * p < 0.1.

We also re-estimate Eq. (3) through an OLS fixed-effects model. The results are in columns 3 and 4 of Table 6, Table 7 for PM2.5 and PM10, respectively. Given the nature of this regression model, the effect of pre-COVID time-invariant exposure to PM is now absorbed by municipality fixed effects. The rationale of this last estimate is to check whether the time-varying PM component is statistically significant when local unobserved time-invariant characteristics are accounted for. Our findings confirm that this is the case.

Table 7

Decomposition of the long-term and short term effects.

	(1)	(2)	(3)	(4)
PM10 (short term component)	0.000494***	0.000489***	0.000563***	0.000577***
	(3.60e-05)	(5.24e-05)	(3.79e-05)	(5.28e-05)
PM10 (ex-ante component)	0.000666***	0.000955***
	(0.000149)	(0.000226)
T (linear day trend)	0.00313***	0.00329***	0.00342***	0.00363***
	(0.000158)	(0.000230)	(0.000161)	(0.000236)
T2 (quadratic day trend)	−8.24e-05***	−8.95e-05***	−6.60e-05***	−7.22e-05***
	(4.80e-06)	(1.07e-05)	(5.00e-06)	(9.19e-06)
T3 (cubic day trend)	4.84e-07***	5.18e-07***	3.94e-07***	4.23e-07***
	(2.76e-08)	(5.53e-08)	(2.85e-08)	(4.89e-08)
Population	−7.63e-05	−0.000132**
	(5.01e-05)	(6.68e-05)
Employees	−0.000394**	−0.000361
	(0.000193)	(0.000303)
Density	−0.00233***	−0.00266***
	(0.000438)	(0.000550)
Employees in Essential sectors	0.00114***	0.00129**
	(0.000368)	(0.000597)
Income	−0.0122***	−0.0165***
	(0.00303)	(0.00486)
Over 65	0.000892***	0.00125***
	(0.000317)	(0.000443)
Days since lockdown	0.000630**	0.000925	−0.000532*	−0.000364
	(0.000280)	(0.000692)	(0.000304)	(0.000582)
Temperature	−0.000225	−0.000293	−0.000423*	−0.000354
	(0.000187)	(0.000361)	(0.000220)	(0.000414)
Mobility	3.53e-05	2.51e-05	5.26e-05**	4.41e-05
	(2.49e-05)	(6.58e-05)	(2.50e-05)	(6.30e-05)
Constant	−0.0220***	−0.0247***	−0.0155***	−0.0132***
	(0.00285)	(0.00450)	(0.00250)	(0.00473)
Municipality fixed effects	No	No	Yes	Yes
Observations	685,385	685,385	685,385	685,385
Log Likelihood	37,941	7532	42,969	12,384
Number of Municipalities			7260	7260

Columns (1) and (2) pooled estimates, columns (3) and (4) fixed effect estimates. Standard errors clustered at municipality level in parentheses. Contains modified Copernicus Climate Change Service Information [2017–2020], DOI: 10.24381/cds.e2161bac. Contains modified Copernicus Atmospheric Monitoring Service Information [2017–2020]; *** p < 0.01, ** p < 0.05, * p < 0.1.

The first robustness check we perform features the use of an alternative instrument calculated as the residual from the following regression:

The residual η is, by construction, exogenous when used as an instrument in our benchmark estimate in Eq. (1). The advantage of this instrument is that through Eq. (4), we control for the complex pattern of relationships through which rain and mobility can affect the relationship between pollution and excess deaths.

The new IV findings confirm that this instrument is also relevant since first-stage regression coefficients are significant. Moreover, the falsification exercise of introducing the instrument in non-instrumented estimates confirms that the former has no significant direct impact on the dependent variable. In terms of magnitude, we note, however, that the coefficient size of the instrumented variable is much higher in the new IV estimates than in the non-IV ones.

We also test whether the short-term effect estimated in our decomposition exercise presented in Table 6, Table 7 remains significant when instrumented under our two different IV approaches. We find this to be the case (Panel 8.5, columns 1 and 2).

There are two additional potential concerns in our estimates: (i) heterogeneity of the pandemic dynamics at the municipal level; and (ii) spatial dependence of the pandemic.

With regard to the first concern, we take two approaches. First, we estimate the Pesaran and Smith (1995) mean group estimator model where slope coefficients are separately calculated for each municipality and averaged across all municipalities. Our main variables of interest remain strongly significant. However, this approach corrects more for heterogeneity of PM impact than of the virus spread net of the PM effect. We, therefore, estimate this model with a mean group estimator specification allowing for province-specific trends. Again, our main results are unchanged (Table 8 , panel 8.5, column 4). Second, we test whether our findings are confirmed when data are aggregated at the province level as the problem of heterogeneous infection dynamics is particularly severe at the municipal level, but less so at the province level. Our main findings are confirmed in non-instrumented and instrumented specifications with province-level data (Table 9 ). Finally, we check for the contemporaneous presence of PM between and within effects to test whether particulate matter has an impact through both effects at the municipal level. This is another way to address the heterogeneity of pandemic dynamics problem since PM between-effects cannot be affected by such a problem. To this purpose, we estimate hybrid models that split the effect of particulate matter into within- and between-municipality effects (Schunck, 2013; Schunck and Perales, 2017) using a Mundlak (1978) random-effects approach. The estimated findings show that both between and a within municipality variation in PM2.5 and PM10 significantly matter in explaining variation in excess deaths. The within-effect, however, has more power since it accounts for three-fourth of the overall effect in the decomposition estimated in the hybrid model (Table 8, panel 8.5, column 3). Note that this decomposition allows us to disentangle contemporary between and within effects; this is a different approach from that proposed in section 4 Eq. (3), where the between effect is long term, lagged, and aims to capture previous long-term exposure to particulate matter.

Table 8

Robustness checks: day effects and alternative specifications.

Panel 8.1	(1)	(2)	(3)	(4)	(5)	(6)
	Pooled OLS	Pooled IV	Pooled IV (second instrument*)	OLS panel fixed effects	IV panel fixed effects	IV panel fixed effects (second instrument*)
Day effects
PM2.5	0.00141***	0.00174***	0.00284***	0.00178***	0.00203***	0.00692***
	(8.20e-05)	(0.000281)	(0.000478)	(9.39e-05)	(0.000238)	(0.000911)
PM10	0.000907***	0.000645***	0.00182***	0.00101***	0.000824***	0.00271***
	(4.60e-05)	(0.000104)	(0.000290)	(4.84e-05)	(9.72e-05)	(0.000367)
Panel 8.2
Day effects (with inverse distance weights)
PM2.5	0.00158***			0.00201***
	(0.000123)			(0.000185)
PM10	0.00101***			0.00111***
	(7.40e-05)			(9.20e-05)

Second instrument: instrument built as explained in Eq. (4) section 5. Robust standard errors in parentheses. *** p < 0.01, ** p < 0.05, * p < 0.1.

Table 9

Provincial specifications and robustness checks: day effects and other specifications.

Panel 9.1	(1)	(2)	(4)	(5)
	Pooled OLS	Pooled IV	OLS panel fixed effects	IV panel fixed effects
PM2.5	0.00113***	0.00123***	0.000783***	0.00131***
	(0.000189)	(6.73e-05)	(0.000108)	(6.39e-05)
PM10	0.000766***	0.000936***	0.000516***	0.000991***
	(0.000127)	(5.30e-05)	(6.85e-05)	(4.94e-05)

For the second concern, that is, spatial correlation, we run a spatial Durbin model for our panel with the province level data following the approach proposed by Belotti et al. (2017).

Furthermore, to account for other possible endogeneity issues, we build a spatial panel IV model. First, we run the fixed-effects quasi-maximum likelihood estimator on the endogenous regressor against both the instruments and the exogenous covariates of the main model. Then, after getting the control function, ie., the prediction of the overall error component from this regression, we run the full spatial model again, controlling for this component. This allows us to further mitigate the remaining endogeneity of the PM variables (Table 9, panel 9.5, columns 1–3).

To check whether our findings are robust to a more flexible control for the aggregate pandemic dynamics that do not assume any particular functional form, we repeat our estimates by introducing day fixed effects (Table 8, panels 8.1–8.4). Our main results remain significant and the coefficient magnitude do not vary significantly.

In an additional robustness check, we calculate COVID-19 non-synchronous regional trends by assuming independent regional pandemic dynamics. To this purpose, we set the regional contagions at n = 100 and use this conventional number as the starting point of the pandemic trends in each region. This approach allows us to account for unobserved time-varying region-level characteristics. Our main findings do not change after attributing a specific regional trend to each municipal (Table 10 , panels 3–5,8–10, 13 and 16).13

Table 10

Robustness checks.

Panel 10.1 Pooled estimates excluding super-spreader events	(1)	(2)	(3)	(4)
PM10		0.000561***		0.000619***
		(3.89e-05)		(5.54e-05)
PM2.5	0.000972***		0.00113***
	(7.24e-05)		(0.000101)
Observations	650,491	650,491	650,491	650,491
Log Likelihood	43,826	43,812	10,945	10,911
Standard errors clustered at municipality level in parentheses
*** p < 0.01, ** p < 0.05, * p < 0.1

The table reports synthetic statistics for PM coefficients in three different robustness checks. In the first panel we remove provinces of Milan and Bergamo to account for the super-spreader event of the Champions League match Atalanta-Valencia. In the second panel we exclude from the sample observations where moving average rainfalls are above 95th centile. In the third panel we replace the national (linear, quadratic, cubic) trend variables measuring contagion dynamics with nonsynchronous regional trend variables starting from the day of the 100th contagion in the given region. For tables 10.1 to 10.3 and 10.6–10.8, columns (1) and (2) do not weight observations, while columns (3) and (4) use as weight the inverse distance of municipality centroids from the meteorological point of observation. For tables 10.11 to 10.16, columns (1) and (2) presents pooled estimates, while columns (3) and (4) fixed effect estimates.

We further refine our main instrument by ruling out episodes of extreme rainfalls from the sample. More specifically, we eliminate observations where the instrument (rainfall moving average) is above the 95th centile and can be suspected to directly affect excess deaths (Table 10, panels 2,4-5,7,9–10, 12 and 15).

To test whether our findings are robust for “super-spreader” events during the pandemic, we consider the UEFA Champions League match between Atalanta and Valencia that took place February 19, 2020, when around 40,000 Atalanta supporters gathered in the San Siro stadium in Milan for the match.14 To do so, we repeat our estimates by removing data for the Bergamo and Milano provinces. Again, the results remain unchanged in terms of magnitude (the pooled estimate coefficient changes only at the fifth decimal digit) and statistical significance (Table 10, panels 1,4-5,6,9–10, 11 and 14).

Discussion

To compare the magnitude of our results with those of the existing literature, we calculate what our coefficients imply in terms of the impact of 1 μg/m3 of PM on mortality. For the magnitude of the PM effects, the estimated PM coefficients vary between different estimates that look at different sources of variability. However, presenting all of them at least as a robustness check is important to evaluate the robustness and extension of the significance of our findings. For example, the fixed-effect coefficient compared with its pooled estimated counterpart captures only the within-effect controlling for unobserved time-invariant municipality effects. The IV effect depends in turn on the quality of the instrument and corrects for endogeneity problems.

Based on all our different estimates, we conclude that the overall non-instrumented PM2.5 effect can be reasonably estimated in a range between 0.001 and 0.002. The highest coefficient is that of the fixed effect estimates augmented for day fixed effects. The same numbers for the PM10 effect are between 0.0006 and 0.001, also when considering estimates of provincial data in Table 8 and robustness checks in Table 9. Given the average daily mortality rate in Italy in the last four years, the effect implies that one additional μg/m3 of PM2.5 is associated with an increase in mortality rate by 2.9 to 5.59%. This effect is in the range of findings made in other studies, slightly above that estimated in the Netherland (Cole et al., 2020) and below that obtained in the US (Wu et al., 2020) and Northern Italy (Cocker et al. 2020) (see introduction).

Note that the severe lockdown measures adopted at the beginning of March 2020 significantly contributed to air quality. The lockdown, therefore, reduced the short term effect of PM on high mortality rate during the pandemic. To understand to what extent this occurred, we calculated the difference between the average daily PM concentration during the pandemic's first wave (February 2020 to May 2020) and during the corresponding days in the previous two years (2018–2019 average). If we limit our analysis to the Northern regions,15 the difference is above 1 μg. Hence, if we can interpret estimates in section 4 as causal, we may conclude that lockdown measures saved between 1 and 2 extra deaths per 100,000 inhabitants.

Extracting the PM concentration differential using fixed-effects estimates do not change the significance and magnitude, thereby confirming our previous analysis.

This study has a number of caveats and limitations. First, our instruments are relevant, but we cannot test their validity and have to prove that on logical grounds. However, the instruments' lack of significance when included as explanatory variables in the main specification, the control for time-varying mobility, the robustness in the sensitivity analysis when excluding extreme rainfalls supports our exclusion restriction.

Second, data at the municipal level are rarely available and, when available, come from the last census in 2011. Consequently, while our analysis controls for many variables like the share of population aged above 65 and the number of employees, we cannot control for other possible factors influencing COVID-19 contagion. For instance, we cannot control for the number of doctors in a given municipality. However, this characteristic is likely captured by municipality fixed effects at the finest level of geographical disaggregation. In our robustness checks described in section 5, we account for municipal heterogeneous pandemic dynamics looking at between effects, using province trends in mean group estimators, and aggregating data at the province level. Note that the two above-mentioned and all other unobserved components were invariant during our 3-month sample period but updated in time to the 2011 census. Third, similarly to other papers (see Cocker et al. 2020), our dependent variable measures total deaths and does not discriminate between COVID-19 deaths and deaths caused by other diseases. This is because of the heterogeneity of COVID-19-related deaths registration, both over time and across regions, that we have explained when motivating the choice of our dependent variable.

A final consideration relates to the interpretation of our findings on the decomposition between the two-year average and the time-varying component of PM. The first fixed component captures the ex-ante long-term exposure effect, while the second the effect of changes in PM during the pandemic. We do not explicitly and exclusively identify this last component in the “short term effect.” Therefore, it can be questioned whether the time-varying effect derives from the PM capacity to increase survival outside the human body (short term effect) or it may further weaken the capacity of lungs and alveoli to resist respiratory and pulmonary diseases on top of the long-term exposure. While further research could clarify this point, this paper is the first, to the best of our knowledge, to show that historical pre-COVID and contemporary time-varying effects matter.

Conclusions

Two research hypotheses in the literature on the impact of particulate matter on COVID-19 contagion and deaths argue that prolonged pre-COVID exposure to, and contemporary levels of particulate matter can play a positive and significant role.

To test these two hypotheses, we evaluate the impact of particulate matter concentration in Italian municipalities on daily deaths between the first COVID-19 outbreak in Italy and the previous five years. The specific contribution of this study to the literature hinges on the use of the geographically finest controls for concurring factors through municipality fixed-effects, instrumental variable estimates to tackle endogeneity issues within a model taking spatial correlation into account, and the decomposition of the two effects, that is, long-term pre-COVID exposure and time-varying effect.

Our findings show that the impact of both components is positive and significant. Our estimates control for standard time-trend components accounting for the non-linear deterministic evolution of the pandemic, the effects of lockdown measures, and several other controls, such as time-varying mobility. Additional results from province-level data accounting for spatial correlation and instrumental variable estimates addressing endogeneity problems further underline the robustness of our findings. More specifically, taking an average PM2.5 effect estimated across all our different models, we find that particulate matter concentration predicts more than 1231 more deaths if we consider the difference between the municipalities with the highest and lowest average PM2.5 concentration during the first pandemic wave.

Indeed, our instrumental-variable estimates are inevitably subject to discussion and limitations. However, if they can be interpreted as causal, our empirical results have relevant policy implications. They highlight an additional and important reason to contrast particulate matter beyond those already known. For example, beyond the impact of atmospheric phenomena, the sources of PM depend on around 90% of human choices such as domestic heating systems, mobility, agriculture, and industrial production processes. Therefore, urgent steps should be taken to accelerate the transition to frontier technology, reducing each source's contribution to PM.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.