WIXLIB

known to affect cardiovascular and respiratorymorbidity and mortality (Cakmak et al., 2018;Jeong et al., 2017; McGuinn et al., 2017; Yinet al., 2017). Moreover, among 1,596 ItalianCOVID-19 patients who died in the hospitals,and for whom it was possible to analyze cliniccharts, data showed substantial comorbiditiesincluding ischemic heart disease (27.9%); atrialfibrillation (22.4%); heart failure (15.6%); stroke(10.9%); hypertension (70.6%), and chronicobstructive pulmonary disease (17.9%) (IstitutoSuperiore di Sanità, 2020). Biologically, longterm PM exposure may be responsible for achronic inflammation status that causes thehyper-activation of the immune system and thelife-threatening respiratory disorders caused byCOVID-19 (Shi et al., 2020).Some preliminary evidence is now emergingabout COVID-19 that shows a positiverelationship between air pollution and morbidityand mortality. Beyond qualitatively describingthe European Air Quality Index for Northern Italyto argue the causal role of air pollution and therelatively high COVID-19 mortality observedin that region, Conticini et al. (2020) reviewthe most recent existing toxicological andepidemiological literature. Based on existingevidence from other empirical studies, theyclarify the relationship between air pollution,prolonged inflammation and immune systemhyper-activation and immune suppression,and the link between the latter and acuterespiratory distress syndrome, and respiratorymortality. Their paper is important in that itsuggests a clinical and biologically plausibleexplanation to our analysis, but does notprovide statistical evidence in support of thehypothesis. A separate empirical analysisby Becchetti et al. (2020) finds preliminaryempirical evidence that confirms such apositive effect of air pollution on mortality inItaly based on the analysis of death data atthe province level. Similarly, Wu et al. (2020)show a positive association between long termPM exposure and COVID-19-related death inUS counties. Ogen (2020) recently analyseddata from 66 administrative regions in France,Spain, Italy, and Germany, and found that thehighest COVID-19 deaths in these regions wereassociated with five regions of Northern Italythat also corresponded with the highest levelsof atmospheric nitrogen dioxide (NO2). Cole etal. (2020) estimate the same relationship usingNetherlands municipality data and find PM2.5positively associated with COVID-19 cases,hospitalization, and deaths.In this paper, we follow this emerging stream ofthe empirical literature and test the hypothesisthat a higher average long-term exposure toPM2.5 is positively associated with the currentextraordinarily high death toll in NorthernItaly. We decided to focus on PM2.5 because,given the complexity of air pollution, it is quitecommon in air pollution epidemiology studiesto focus the analysis on a single pollutant (Wuet al., 2020), although multipollutant analysesare certainly warranted. We selected PM2.5 fora variety of important reasons, including policyimplications and evidence in the literature interms of chronic health effects. Regardingits policy implications, we selected PM2.5 asopposed to PM10 because the former is morecorrelated with human activities than the latter,and it correlates with stronger health effectsthan PM10 does. With respect to respiratorymortality effects from the existing air pollutionliterature, the most robust evidence pointsto PM2.5 as opposed to other gaseous airpollutants (Bowe et al., 2019).FEEM REPORTS 4

Mortality data are collected at the municipalitylevel for the period January-April 2020. Giventhat mortality data are not disaggregated bymortality cause, death counts are measured asthe difference from the last five-years mean toreflect the abnormal number of deaths causedby the spreading of the pandemic. Since PM2.5can be associated to generic mortality evenin the absence of the pandemic outbreak(Dominici et al., 2003; Katsouyanni et al.,2001; Samet et al., 2000), we also estimatethe impact of PM2.5 on the excess mortality inthe sample using 2019 data, a time in whichthe coronavirus epidemic had presumably notyet begun. Data on PM2.5 concentration at themunicipality level refer to the years prior 2020to account for long-term population exposure.We assign municipality PM2.5 concentrationby a set of different methods of spatialinterpolation (kriging) of monitoring station datarelated to the years 2015-2019.of 2020; and the demographic compositionof the municipal population among others. Inaddition, we consider spatial heterogeneityin the distribution of the number of deathsrelated to regional and local unobservablefactors. We account for region-specific effectsbecause regions, in Italy, are the administrativeunits in charge of managing the healthsystems and the measures taken to traceand contrast the spreading of the pandemicvaried greatly among even contiguous regions.We also account for local effects common tofunctionally linked clusters of municipalities(the Local Labour Systems – LLS). We deemthis part of the identification strategy crucialbecause the relationship between PM2.5 andCOVID-19 mortality may be confounded byseveral other factors, some of which were notobservable or measurable, but are neverthelessintrinsically related to the geographical locationof the observed units.We estimate a negative binomial modelof excessive deaths on historical PM2.5concentrations and a series of controlvariables that may plausibly affect bothPM2.5 concentration and mortality, includingpopulation density; the spatial concentrationof the industrial manufacturing sites; climaticconditions observed during the first quarterThe remainder of the paper is organizedas follows. The next section introduces theempirical strategy and describes the dataset.The results are presented and discussed insection three, considering the total numberof (excess) deaths. Section four draws theconclusions and highlights the limitation of thestudy and the indications for future research.FEEM REPORTS 5

02Empirical strategy and dataOur analysis is restricted to the study area ofNorthern Italy (Figure 1), which encompassesthe sub-regions of Valle D’Aosta, Piemonte,Liguria, Lombardia, Emilia-Romagna, Veneto,Friuli-Venezia Giulia and Trentino-AltoAdige/Südtirol. Official territorial data onCOVID-19 mortality in Italy are available atthe rather aggregate regional or provinciallevel, corresponding to the levels 2 and 3,respectively of the European nomenclatureunits for territorial statistics (NUTS)4. Inaddition, these official data refer to the deathsof patients tested positive for severe acuterespiratory syndrome coronavirus 2 (SARS-CoV2) only and do not include (potential)patients without COVID-19 diagnosis becausethey were not tested and died at home orelsewhere. Hence, the officially reporteddeaths are likely underestimated. Becausetesting policies vary among regions in Italy,the induced measurement error is also nonrandomly distributed among the provinces.Ciminelli & Garcia-mandicó (2020) compare theofficial COVID-19 fatality rates with historicaldeath data and report that deaths were higherthan official fatalities throughout the period ofCOVID-19 epidemic.Figure 1– Italian regions included in the studyRegion:Emilia-RomagnaFriuli V.G.LiguriaLombardiaPiemonteTrentino A.A.Valle ts/backgroundFEEM REPORTS 6

Working under the assumption that COVID-19deaths are underestimated in Italy, the choiceis made in this paper to use the total deathsfrom the official registries, accordingly, andto scale the analysis at the municipality level,the smallest administrative units, to have amore granular representation of the spatialdimension of the phenomenon. Since weare interested in excess deaths, we take thedifference between the number of deaths inthe period January 1 – April 30, 2020, and theaverage number of deaths in the same periodof the previous five years (ExDeaths) and usethis ExDeaths metric as the dependent variablein our statistical model. Figure 2 displaysthe geographical distribution of the abovedescribed data among the 4041 municipalitiesfor which data is available.Figure 2: Spatial distribution of cumulative excess deaths in sample municipalities, Northern Italy, January 1 – April 30, 2020The variable is assumed to follow a NegativeBinomial distribution, a generalization of thePoisson distribution that avoids the restrictivemean-variance equality of the latter, and ismodelled as follows:(1) ExDeathsi log(μi)NB(μi, θ) α βPMi δ' Xi εiwhere θ is the overdispersion parameterto be estimated and μi is the municipalityspecific expectation conditional on the valueof the covariates. Among the covariates, PMis the concentration of fine particulate matterin municipality i and β is the associatedparameter, which we expect positive andstatistically different from zero; X is a vector ofcontrol variables that adjusts for the potentialconfounding effects and includes the (logof) total population as the offset while ε is anormally-distributed error term.Our main source of PM2.5 data is the EuropeanEnvironmental Agency’s (EEA) air monitoringdatabase, which is provided to EEA by theFEEM REPORTS 7

Institute for Environmental Protection andResearch (ISPRA). ISPRA conducts ground-levelair measurements of PM2.5 air concentrations(µg/m3) collected at 268 monitoring sitesthroughout Italy. Specifically, we use the EEA’sE1a and E2a datasets, which are primaryvalidated assessment data and primaryup-to-date assessment data reported bythe European Member States, respectively.Although the measurements come both inhourly and daily averaging formats, we workwith daily values and use them to obtain yearlyaggregates for the years 2015, 2016, 2017,2018, and 2019. However, because model(1) does not include a time component, wefurther compute a six-year averaging time toobtain a metric of long-term (chronic) PM2.5concentration levels throughout differentspatial units of Northern Italy. The number ofsix years for the reference period is sufficientlylong to account for long-term exposure whilebeing not too long to be affected by the mobilityof people among municipalities, and it is inline with existing literature assessing long-termeffects of PM exposure (Yorifuji et al., 2019).Since the air monitoring stations provide onlypartial spatial coverage for municipality-levelPM2.5 concentration data, we impute missingobservations using a spatial interpolationmodel. Specifically, we fill in the gaps usinga mean stationary Ordinary Kriging (seeBivand et al., 2013 p 209) defined through anexponential covariance function with nugget,partial sill and range parameters estimatedthrough (restricted) maximum likelihoodmethods. Figure 2 displays the resulting PM2.5concentration data5.Figure 3: Spatial distribution of PM2.5 concentration levels in the sample municipalities, simple kriging of monitoring stations,average across years 2015-20195We also replicate the analysis using other trend-stationary models (i.e. universal kriging) and different covariance functions; theseextensions are discussed in the robustness check section.FEEM REPORTS 8

Comparing figures 2 and 3, it is possible tovisually appreciate a spatial coincidencebetween higher levels of excess mortalityand higher levels of PM2.5, in particular in theLombardia region which notably is the regionwith both the highest particulate concentrationand the highest number of excess mortality.The hypothesis that PM2.5 concentrationaffected COVID deaths, that is (β̂1 0), is testedamong several possible specifications. Inmodel (2) we include regional effects (λj). Theseeffects are expected to capture the aspectsrelated to the management of the outbreak,which may have systematically influencedCOVID-19 mortality and that are common to allthe municipalities in the same region. Italy hasa national health system that ensures equalaccess to healthcare to all citizens. The systemis managed by regions at the local level, and,in the specific case of this pandemic, regionswere responsible for defining the testing andcontact-tracing protocols and implementing thenecessary measures to contain the outbreak,among which the measure to protect healthcareworkers. In model (3), we include LLS-specificeffects (ek). LLS are spatial clusters ofcontiguous municipalities related by commutingflows that share a common specialization ina specific sector of manufacturing productionand correspond to the conceptualization ofMarshallian districts (Becattini, 2002). Thenumber of LLS clusters per-region and thetotal number of municipalities belonging toclusters are reported in table 1, along with theminimum, maximum, and average cluster size.Table 1: Number of LLS spatial clusters in each region.RegionN. LLSN. of municipalitiesSmallest LLS (N. ofmunicipalities)Largest LLS (N. ofmunicipalities)Average number ofmunicipalities by LLSEmilia-Romagna423281388Friuli-Venezia 7425Piemonte391181110426Trentino-Alto Adige/Südtirol2729113010Valle d'Aosta57432912Veneto4956315211The use of LLS captures the interlinkageswithin neighbouring municipalities that mayhave favoured the geographical spreadingof coronavirus around specific hotspots.Mortality data are then expected to vary amongmunicipalities in different LLS, but differencesare expected to be non-systematic in this case.In model (4) we include both the regional fixedeffects and the LLS random effects.(ek)(2) ExDeathsi log(μij)NB(μij, θ) α βPMij δ' Xij λj εij(3) ExDeathsi log(μik)uikNB(μik, θ) α βPMik δ' Xik uik εik ek(4) ExDeathsi log(μijk)uijkNB(μijk, θ) α βPMijk δ' Xijk λi uijk εijk ekFEEM REPORTS 9

Control variables to be included in the modelwere chosen to avoid any potential spatialconfounding effect and considering as wellthe emerging literature on the impact of PMon COVID-19 related deaths (Cole et al., 2020;Wu et al., 2020). The population density andper-capita income account for urbanisationlevel. The most densely populated and wealthymunicipalities are among the most polluted dueto the spatial concentration of manufacturingand service activities but are also the placeswhere the contagion could have been easier,with a potential impact on mortality. In additionto the density of population, the shares ofmunicipality area occupied by industrial sitesand the average size of manufacturing firmsare included in the regression because theyare related to pollutant concentration andpossibly to mortality. National measuresto stop the spreading of the viral infection(lockdown) involved the service sector to thelargest extent while many manufacturingactivities, being considered necessary, wereleft open and, in the absence of social distanceand individual protection measures, thegeographical concentration of these activitiesin a municipality with their complex logisticsand transport interconnections, and the size ofplants, may have influenced mortality. Averagetemperature, for which an association withCOVID-19 deaths has also been found (Ma etal., 2020), is also included in the regression6.Moreover, COVID-19 incidence has proven tobe higher among men than women and peopleaged 65 or more. Hence these two variablesare considered in the model, even though theseaspects are not necessarily connected withthe average PM2.5 exposure in a municipality.Underlying socioeconomic conditions can also6play a role in COVID-19 related mortality (Goutteet al., 2020). Brandt et al. (2020) and Mukherji(n.d.) have shown that, in the US, COVID-19is more threatening for ethnic minorities,and we believe that the share of migrants,identified as non-EU citizens, can control forthis aspect influencing the observed excessmortality. On the other hand, Mukherji (2020)and Goutte et al. (2020) also find that placeswith a higher share of the population with alow level of education have higher deaths. Inour paper, given the lack of updated data oneducation at the municipal level, we proxy itwith the percentage of university students onthe total population. The distance from theclosest airport is a proxy for the functionaland relational linkage between a municipalityand a place of highly frequent national andinternational connections and potential sourcesof coronavirus spreading. Finally, we considerthe number of hospital beds as a proxy for thesupply of health services to account for the factthat many people died at home without beingdiagnosed for coronavirus due to the shortageof beds in public structures. The full details ofthe variables in the model, including sourcesand summary statistics, are presented in Table2.We omit average humidity because the variable shows a too strong linear correlation (rho 0.98) with temperature in our sampleobservations, and its inclusion would cause severe imperfect collinearity in the model.FEEM REPORTS 10

Table 2: Description of model variables and summary sample statisticsVariableDescriptionExDeathsNumber of deaths in the period January 1 – April 30 2020 – absolutedifference compared to the average of the past five years, source: ISTATMeanMedianSD9.32237.13PM2.5Fine (2.5 µg/m³) particulate matter concentration obtained by spatialinterpolation of monitoring stations, average across the years 20152019, source: European Environmental Agency19.6720.854.15Pop. DensityPopulation density computed as total population in number of inhabitantson January 1 2020 over the total artificial area in Km2, sources: ISTATand European Environmental Agency – Corine land Cover data34.4414.6958.19PC IncomeAverage per-capita income, source: Ministry of Finance, 201915658.4615564.522497.44% Ind. LandShare of industrial area on total municipality surface measures throughsatellite observation, source: European Environmental Agency – Corineland Cover data2.620.025.17% Small EntShare of enterprises with less than 10 employees, source: Registrostatistico delle Unità Locali (ASIA - UL)94.2694.443.80TemperatureAverage mean skin temperature during the death observation period,source: Copernicus ERA5 0.25 x0.25 grid resolution dataset.3.755.344.01Female/MaleRatio between female and male population, source: ISTAT% Over 65Share of population older than 65, source: ISTAT% non-EUShare of non-EU residents, source ISTAT% Univ. Stud.Definition, source: share of University students over total population,source: Ministry of University and ResearchDist. AirportDistance in meters to the closest Airport, source: our computation basedon European Environmental Agency – Corine land Cover dataPC HospitalBedsNumber per-capita hospital beds in the municipality, source: HealthMinistryPopulationTotal population, source: ISTATHaving accounted for the confounding effectdue to the omission of relevant informationfrom the empirical specification, we excludeany other potential source of endogeneityconsidered in similar papers. In particular,we exclude endogeneity due to measurementerror in the outcome variable and the mainindependent variable. Concerning the outcomevariable, the relationship between deaths andcases with fine PM could be spurious becausemore cases could be registered, and moreindividuals tested in highly polluted areas aspeople there are more likely to show COVID-19symptoms due to the chronic inflammationinduced by PM. The high toll of deaths ofpeople diagnosed with COVID-19 would be anatural consequence of that. In contrast, thenumber of deaths in excess, used in this paper,is not affected by testing problems since itconsiders all the potential COVID-19 .32254938814.47Concerning the PM variable, measurementerrors are likely to occur when using satellitedata or modelled data. We preferred to usePM2.5 levels observed from monitoring stationsto avoid such a measurement error. Somecaution is needed in the spatial interpolationbecause the method chosen to fill the missingdata may underestimate the value in locationsfarther from the monitoring stations. With thisconcern in mind, we test the robustness of ourresults using PM2.5 data obtained from differentinterpolation approaches.FEEM REPORTS 11

03ResultsAs in

The opinions expressed in this paper do not necessarily reflect the position of Fondazione Eni Enrico Mattei Corso Magenta, 63, 20123 Milano (I), web site: www.feem.it, e-mail: working.papers@feem.it