A Local Spatial STIRPAT Model for Outdoor NOx Concentrations in the Community of Madrid, Spain

Air pollution control is one of the main challenges facing modern societies. Consequently, the estimation of population, affluence, and technology impacts on air pollution concentrations (STIRPAT modeling) has become the cornerstone of environmental decision-making. Spatial effects are not usually included in STIRPAT modeling of air pollution. However, space matters: accounting for spatial dependencies significantly improves the accuracy of estimates and forecasts, especially (or only) when dealing with small information units rather than with large ones (countries, large regions, provinces in China, counties and states in the USA, etc.). The latter scale is typical in the literature on air pollution due to the difficulties in finding data on its drivers at a true local scale. Accordingly, this paper has a double objective. The first is the estimation of a spatial panel data STIRPAT model, with the spatial units being both very small and also highly autonomous, developed municipalities. The second is to examine whether an environmental Kuznets curve relationship exists between income per capita and NOx concentrations. A case study has been carried out in the Autonomous Community of Madrid, Spain, at the municipal level.


Introduction
Air pollution is currently the most important environmental risk to human health, and it is perceived as the second biggest environmental concern for Europeans, after climate change (European Commission [1]). Therefore, it is no surprise that there is growing political, media and public interest in air quality issues and increased public support for action. As outlined by the European Environment Agency (EEA [2]), growing public engagement around air pollution challenges, including ongoing citizen science initiatives engaged in supporting air quality monitoring (EEA [3]) and initiatives targeting public awareness and behavioral changes, have led to a rise in support and demand for measures to control and improve air quality.
Although there have been substantial reductions in emissions of air pollutants in recent decades (ETC/ACM [4]), a large percentage of the European population continues to be exposed to pollutant emissions, particularly in urban areas. Air pollution has significant impacts on Europeans' health. In 2016-the most up-to-date data from the European Monitoring and Evaluation Programme (EMEP) model, which involves 41 countries-412,000 premature adult deaths were attributable to long-term exposure to particulate matter with a diameter of less than 2.5 micrometers (PM 2.5 ) concentrations while the estimated impacts of exposure to nitrogen dioxide (NO 2 ) and tropospheric ozone (O 3 ) concentrations were around 71,000 and 15,100 premature deaths per year, respectively (EEA [2]). In spite of the reductions in emissions of air pollutants in recent decades, In this paper, the relationship between outdoor air pollution concentrations, economic activity and socio-demographics at local level is analyzed while controlling for potential spatial effects. Analyzing these relationships is important in the Community of Madrid, which has registered a very significant increase in its economic activity over the last two decades, as a result of which it is currently the most densely populated region in Europe. The social and economic development has caused a marked increase in the population and subsequent agglomeration, a remarkable change of the urban structure of the region, and a greater demand for transportation, consequently driving up the levels of energy consumption and air pollutant emissions. Programs to change citizens' travel behavior are essential for reducing transport pollutant emissions and the negative effects of air pollution on human health. However, as a result of spatial spillovers, policies implemented in one municipality may influence policies in neighboring municipalities, thus implying that the various municipalities should interact to balance pro-environmental and economic policy goals. This issue is of crucial importance in order to convince policy-makers to act in a coordinated way. In fact, the lack of coordination among mayors of different political leanings is one of the problems that prevent further progress in air pollution control.
In this article, we specifically focus on NO x (which includes not only NO 2 but also NO), the most relevant nitrogen oxides for air pollution. The reasons for focusing on NO x instead of on other criteria air pollutants (CO, Pb, O 3 , PM, and SO 2 ) are numerous and varied. First, the inevitable "clean air-economic development" trade-off in modern societies implies internal combustion vehicles, and the main source of NO x concentrations in the air in cities is emissions caused by cars (at least until we achieve mass adoption of electric vehicles), especially those powered by diesel engines. Therefore, since economic development entails transport demand and NO x is an immediate consequence of said transport demand (despite the technological advances that have taken place in terms of reducing emissions from traffic), it can be considered a good indicator of pollution due to motorized traffic in urban areas. Unfortunately, traffic will continue to be a major cause of air pollution problems, at least in this decade, and this will undoubtedly be the case in the Autonomous Community of Madrid. Second, NO x is the pollutant that causes the secondlargest number of deaths after PM 2.5 . The most important source of particulate matter (PM 10 and PM 2.5 ) is also emissions generated by road traffic (direct emissions of primary particles from the exhaust pipe of internal combustion vehicles, as well as the resuspension of materials that accumulate on the paved roads: mechanical abrasion products of vehicles, brakes, wheels, emissions derived from construction works). Why then would we not focus on PM 10 or PM 2.5 instead of NO x ? The answer is given by the third, and especially the fourth reasons detailed below.
The third reason is that NO x not only causes highly damaging direct effects on the environment, but also plays a part in various chemical reactions that take place in the atmosphere, giving rise to both the production of O 3 and PM 2.5 , which are the most damaging to the environment. Therefore, when considering the effects of NO x on health, it is not only the direct effects that must be taken into account, but also its status as an indicator of pollution due to traffic (as stated above) and as a precursor to other highly damaging pollutants.
The fourth and most important reason is that NO x is the air pollutant that most often exceeds the legal limits in urban areas. In the specific case of the region of Madrid, NO x is the only criteria pollutant that currently exceeds the standard limits, albeit only in a few monitored sites. Fortunately, the other criteria pollutants, including PM 10 and PM 2.5 are not a problem in the region of Madrid: they are under control, and have been for many years in the case of CO and SO 2 , the levels of which are currently far below the limit set by the legislation for health protection. It is true that, in the summer months, high values of O 3 are recorded because O 3 is a secondary pollutant formed from a series of primary precursors with high insolation and temperature conditions. Among these primary precursors, the main ones are NO x and volatile organic compounds. Therefore, an O 3 episode is ultimately a NO x episode when insolation and temperature are high. However, it must be borne in mind that the O 3 alert threshold for the population has never been exceeded in the region of Madrid.
There is also a fifth (technical) reason. This article methodologically extends the article by Laureti et al. [8] by letting spatial dependencies enter into play in the modeling of air pollution concentration levels. Since Laureti et al. [8] focused on NO x concentrations in the region of Madrid, for the sake of comparison this article also uses NO x concentrations as the response variable and the same covariates data base as Laureti et al. [8].
Therefore, this study extends previous models that use NO x as the response variable (Laureti et al. [8]) by explicitly considering and testing for the types of spatial dependence within the relationship between transport-related emissions and socio-economic drivers at municipal level, where the relationship between socio-economic factors and the environment becomes more complex. Accordingly, this study proposes a spatial STIRPAT (Stochastic Impacts by Regression on Population, Affluence, and Technology) model at local scale which takes into account the effects of spatial dependencies inherent in air pollution concentrations.
In addition, we test the Environmental Kuznets Curve (EKC) hypothesis, which describes an inverted U-shaped relationship between environmental quality and economic development, by considering data at a disaggregated level of territorial unit (in our case, at a very local scale). According to Atwi et al. [9], it is very important to control for the spatial effects in the equation. The inverted U-shaped hypothesis is reinforced when it is estimated in a spatial setting, which means that geography is not neutral because technology and social capital, key elements in managing emissions, are not evenly distributed over space. Strategic interaction among governments is another factor that apparently stimulates spatial interaction but is not particularly relevant in our case. Heterogeneity and structural breaks are two other issues that must be treated properly. These three effects, and especially the first two, make the panel framework an excellent candidate for dealing with such issues.
In recent decades, the STIRPAT framework has been widely adopted to explore the environment-development nexus at macro level, involving cross-national data, and to a lesser extent at local scale, yielding evidence that different socio-economic factors (i.e., urbanization and population density, age structure, income and technology) influence energy consumption and pollutant concentrations.
Nevertheless, from a methodological point of view, while researchers have examined how several econometric issues might explain the lack of robust results in the STIRPAT framework, less attention has been paid to the spatial nature of the data. This is not the case for other air pollution related areas; for example, in the EEA's estimation of health impacts attributable to exposure to air pollution, the methodology involves the use of maps of interpolated air pollutant concentrations, with information on the spatial distribution of concentrations from the EMEP model. Although the STIRPAT model has usually been employed with national level data, over the last few years there has been a growing interest in the analysis of socio-economic drivers of air pollutant emissions at local level (Elliott and Clement [10]). However, moving down in scale-from national to local-also allows us to incorporate and assess spatial dependencies in the STIRPAT approach, since it is reasonable to assume a strong spatial correlation in air quality of geographically close areas.
To date, there have been various empirical studies addressing the spatial dependence issue when exploring the socio-economic drivers of environmental pollution, but the analyses have been mainly carried out at aggregate level (national or regional). Moreover, most of the studies are aimed at testing the EKC hypothesis of the inverted U-shaped relationship between pollutant emissions and income.
However, according to Mosconi et al. [11] and the references therein, empirical evidence suggests that environmental externalities play different roles at the country, regional and local scales.
Consequently, the sub-regional scale is assumed to be a particularly appropriate spatial level to delineate the different forms of the human-environment system, to explain their development in the past, and to suggest new directions for future policy governing people's relationships with the environment (Kairis et al. [12]). The structure of production systems and the interaction between environmental drivers and local communities call for a more comprehensive analysis of the integrated policy response, distinguishing between actions carried out by central governments, regional institutions, and local authorities (Brenner [13]). In this regard, verifying a spatially explicit environment-economy relationship at the different operational scales contributes to a refined understanding of the specific impact of any policy strategy within a range of different 'macro' (cross-country or crossregion) and 'micro' (local community, cross-district) dimensions (Destek and Sarkodie [14]; Destek et al. [15]; Expósito et al. [16]) and can help in policy-making at the decentralized governance levels (Li et al. [17]; Lin et al. [18]; Gill et al. [19]; Luzzati et al. [20]).
It is of note that the above literature considers "local scale" to include very large spatial units such as counties in the USA, provinces in China, etc. While it is true that they are sub-national in scale, it is hard to consider US counties and Chinese provinces as representing "local scale". We should also point out that the abovementioned finding by Mosconi et al. [11] is in accordance with both the main principles of geostatistics and spatial econometrics and the literature on air pollution control: significant spatial autocorrelation in air pollution concentrations and in the drivers of such concentrations appears at relatively short distances and disappear as the distances among the spatial units that make up the area under study increase. Therefore, according to the abovementioned principles, which are summarized in the well-known Tobler's law, we should expect to find a notable difference in the various spatial (and also temporal) scales considered. In addition, the use of spatial models to estimate STIRPAT specifications and test for the existence of EKC is questionable when dealing with large spatial units (countries, large regions, or supposedly local scales such as large provinces in China or states in the USA, among others) as basic information units. Usually, the reason for using these large spatial units is the difficulties in finding data on the drivers of air pollution concentrations at a more local level.
Summarizing: (i) Spatial effects, and more specifically spatial dependencies, are present in the economic activity-air pollution concentrations relationship. Therefore, testing and estimating the impact of spatial dependencies on the above relationship can be considered a hot research topic in the field of pollution control. (ii) Since the STIRPAT framework has been the most popular approach adopted to explore the environmentdevelopment nexus, it is crucial to incorporate the spatial dependencies existing in series of air pollution concentrations and the drivers of these concentrations in a STIRPAT framework. (iii) There have been some attempts at addressing the spatial dependence issue when exploring the socio-economic drivers of environmental pollution, but at aggregate level (countries, large regions, provinces in China, counties in the USA). Some studies claim to have used a local level scale, but what was called "local scale" in reality referred to very large spatial units. In addition, most of the time the real focus of those papers was on testing the EKC hypothesis. (iv) According to Tobler's Law, spatial effects have a notable impact on the economic activity-air pollution concentration relationship at a local scale, since it is reasonable to assume a strong spatial correlation in air quality of geographically close areas. Therefore, we should not be surprised by the lack of robust results obtained from the STIRPAT modeling of air pollution concentrations with and without considering spatial effects, because large scales are incompatible with significant spatial autocorrelation.
Why is there so little literature on STIRPAT modeling of air pollution concentrations (none at all in the case of the spatial STIRPAT approach) at a local scale? The reason is the absence of information about pollutant concentrations and their drivers at that scale. However, local scale and spatial STIRPAT modeling go together. This is precisely the gap we fill with this paper: the estimation of a spatial STIRPAT model for air pollution concentrations at a true local scale. At larger scales, it is more than probable that spatial autocorrelation, if it exists, is either spurious or is capturing other different effects, especially spatial heterogeneity. Therefore, we want to draw attention to the disappointing results yielded by the literature on spatial STIRPAT modeling using large information units. The overarching idea that we want to contribute to the literature on the topic is that spatial STIRPAT modeling goes hand in hand with a local scale approach. Spatial STIRPAT modeling at a true local scale is a challenging task in terms of collecting the necessary data to estimate the model. However, in this paper we show that modern techniques of geostatistics and the collection of socio-economic data at local level/generation of small-scale socio-economic data can be of great help. Our hope is that, in the future, the approach followed in this article will be adopted by many more researchers, so that comparisons can be made and the robustness of results can be checked.
In light of the above arguments, it can be seen why this article has a double objective: first, the estimation of a spatial (panel data) STIRPAT model, where the spatial units are both very small and also highly autonomous and developed municipalities; and second, to examine whether an EKC relationship exists between gross disposable income and NO x concentrations in the air. For this purpose, a case study has been carried out in the municipalities that make up the Autonomous Community of Madrid, Spain, where outdoor NO x concentrations are still one of the most concerning pollution problems.
In brief, this article contributes to filling the "true local scale" gap in the literature by exploring the issue of spatial autocorrelation inherent in air pollution at short-to-medium distances within the STIRPAT framework and at a local scale.
It is reasonable to assume that there is spatial dependence in both pollution concentrations and also in most of the drivers of transport-related emissions and other economic forces that cross municipality borders, which may create spatial spillovers. Consequently, we use a spatial Durbin panel data STIRPAT model to account for and estimate such spatial dependencies. Moreover, for comparison purposes, we analyze this spatial dependence by specifying a spatial lag panel data STIRPAT model and a spatial error panel data STIRPAT specification.
Our analysis focuses on the municipalities that make up the region of Madrid, over the period 2000-2009. The reasons behind the selection of this period are explained in Section 2. It is worth noting that only the highly populated municipalities of the Community of Madrid are equipped with one or more MS. In order to overcome the lack of disaggregated pollution data in the municipalities that are not equipped with MS, we use spatio-temporal kriging, which accounts for spatio-temporal dependencies and enables us to estimate the level of NO x in the municipalities with no MS.
The remainder of the paper is structured as follows. Section 2 provides a description of the data and methods while Section 3 reports the empirical results. In Section 4 some conclusions are drawn.

Data
Our analysis is based on the 179 municipalities of the Community of Madrid over the period 2000-2009. There are several reasons behind the selection of this period. The first is that it is the same as the period analyzed in Laureti et al. [8], which allows for the comparison of aspatial and spatial STIRPAT models. The second and more important one is that traffic emissions are considered as the main source of air pollution in the Community of Madrid (especially NO x concentrations); therefore, the fleet of vehicles per capita is assumed to be the most relevant explanatory variable in the analysis. However, technical advances in fuel efficiency are continuously being incorporated in private vehicles and, fortunately, ever greater reductions in emissions are being achieved. Given this fact, we need information about the "equivalent" fleet of private vehicles, that is, a series of data where the number of cars for year t is expressed in terms of the number of cars of a base year that generate the same amount of pollution. This difficult task was done in Laureti et al. [8] for the period 2000-2009, and we take advantage of this. (iii) The third reason is that in 2010 the Atmosphere Pollution Monitoring System (APMS) of the city of Madrid (the most important municipality of the Community of Madrid; 25 of the 48 MS operating in the Community of Madrid are located in the city of Madrid) was reorganized. The number of MS was substantially reduced and ecologist associations suspected that the municipality of Madrid removed stations from the sites that were potentially the most polluted (Montero and Fernández-Avilés [35]) to manipulate the data. While we could have taken on the task of updating the "equivalent" fleet of private vehicles per capita (a real challenge), we decided against it because of the smaller number of MS after the APMS reorganization and the doubts about the optimality of their location expressed by ecologists associations, which might substantially affect the accuracy of the response variable data. fortunately, ever greater reductions in emissions are being achieved. Given this need information about the "equivalent" fleet of private vehicles, that is, a serie where the number of cars for year t is expressed in terms of the number of cars o year that generate the same amount of pollution. This difficult task was done in L al. [8] for the period 2000-2009, and we take advantage of this. (iii) The third reaso in 2010 the Atmosphere Pollution Monitoring System (APMS) of the city of Mad most important municipality of the Community of Madrid; 25 of the 48 MS ope the Community of Madrid are located in the city of Madrid) was reorganized. The of MS was substantially reduced and ecologist associations suspected that the mu ity of Madrid removed stations from the sites that were potentially the most (Montero and Fernández-Avilés [35]) to manipulate the data. While we could ha on the task of updating the "equivalent" fleet of private vehicles per capita (a r lenge), we decided against it because of the smaller number of MS after the APM ganization and the doubts about the optimality of their location expressed by ec associations, which might substantially affect the accuracy of the response variab    As can be seen in Figures 1 and 2, there are different types of MS (traffic, industrial, and background) depending on their location and the impact (or absence) of nearby emissions (pollutant-specific in 2011/850/EU). They also show the type of area, which refers to the environment on a scale of several kilometers. More specifically, the region of Madrid includes three urban areas (Corredor del Henares, Urbana Sur and Urbana Noroeste), three rural areas (Cuenca del Tajuña, Cuenca del Alberche and Sierra Norte) and a seventh zone managed by the City of Madrid Council, which has its own air quality monitoring network. In the city of Madrid, which is divided into five zones (inner M-30, northwest, northeast, southwest, southeast) the APMS is made up of 24 MS. Half of them are background urban, and they are representative of the exposure of the urban population in general. Nine are near traffic, and are located such that their level of pollution is mainly influenced by emissions from a nearby street or road (but it is necessary to avoid measuring very small micro-environments in their vicinity). And three are suburban, located on the outskirts of the city, in places where the highest levels of ozone are found. Background urban and traffic MS are aimed at the protection of human health, while the suburban ones monitor pollution to protect human health and vegetation.
In the rest of the municipalities of the region, in the zoning for NOx, eight out of the 24 stations of the network are traffic-related, two industrial and 14 background.
It is of note that 155 municipalities out of the total of 179 are not equipped with an air pollution MS. In other words, there are no data on NOx, the response variable of the STIRPAT approach, for the aforementioned 155 municipalities in the period 2000-2009. Consequently, it is necessary to estimate the value of NOx in these municipalities. This could be seen as an insurmountable problem; however, as pointed out in Laureti et al. [8], the solution lies in the fact that air pollution data exhibit spatio-temporal dependence (Fernández-Avilés et al. [36]; Montero-Lorenzo et al. [37]), not just spatial dependencies (as mentioned in Laureti et al. [8]). Based on the registers of the existing MS, we used (i) spatio-temporal kriging of the mean (STKotM) for 2000 to 2009, to predict the response As can be seen in Figures 1 and 2, there are different types of MS (traffic, industrial, and background) depending on their location and the impact (or absence) of nearby emissions (pollutant-specific in 2011/850/EU). They also show the type of area, which refers to the environment on a scale of several kilometers. More specifically, the region of Madrid includes three urban areas (Corredor del Henares, Urbana Sur and Urbana Noroeste), three rural areas (Cuenca del Tajuña, Cuenca del Alberche and Sierra Norte) and a seventh zone managed by the City of Madrid Council, which has its own air quality monitoring network. In the city of Madrid, which is divided into five zones (inner M-30, northwest, northeast, southwest, southeast) the APMS is made up of 24 MS. Half of them are background urban, and they are representative of the exposure of the urban population in general. Nine are near traffic, and are located such that their level of pollution is mainly influenced by emissions from a nearby street or road (but it is necessary to avoid measuring very small micro-environments in their vicinity). And three are suburban, located on the outskirts of the city, in places where the highest levels of ozone are found. Background urban and traffic MS are aimed at the protection of human health, while the suburban ones monitor pollution to protect human health and vegetation.
In the rest of the municipalities of the region, in the zoning for NO x , eight out of the 24 stations of the network are traffic-related, two industrial and 14 background.
It is of note that 155 municipalities out of the total of 179 are not equipped with an air pollution MS. In other words, there are no data on NO x , the response variable of the STIRPAT approach, for the aforementioned 155 municipalities in the period 2000-2009. Consequently, it is necessary to estimate the value of NO x in these municipalities. This could be seen as an insurmountable problem; however, as pointed out in Laureti et al. [8], the solution lies in the fact that air pollution data exhibit spatio-temporal dependence (Fernández-Avilés et al. [36]; Montero-Lorenzo et al. [37]), not just spatial dependencies (as mentioned in Laureti et al. [8]). Based on the registers of the existing MS, we used (i) spatiotemporal kriging of the mean (STKotM) for 2000 to 2009, to predict the response value in the period under study in the municipality of Madrid, which was the only municipality with more than one MS, and (ii) spatio-temporal ordinary kriging (STOK) to predict this value in the non-MS-equipped municipalities.
Although both interpolation procedures will be detailed in the methodological section of the paper, Table 1 shows the NO x mean value (registered or estimated by kriging) for municipality clusters based on their population size (thus ensuring the reader is provided with the information while keeping the Table manageable). As shown in Table 1, the average level of NO x concentrations in the city of Madrid (72.48 µg/m 3 in 2009), although displaying a marked downward trend overall, indicates that the legal standard for health protection (40 µg/m 3 , annual mean) was clearly exceeded at many of the MS located in the city. The same can be said for the industrial municipalities located in the Southern and Eastern Metropolitan area (most of them in the group with 100,000-250,000 inhabitants). However, in the municipalities with fewer than 100,000 inhabitants (most of them located in wealthy, verdant, nicely urbanized municipalities, close to the mountains, in the Northern and Western Metropolitan area and on the outskirts of Madrid) the annual mean level of NO x was much lower (slightly below 30 µg/m 3 ), and only exceeded the legal standard in some specific municipalities in 2002 and 2005.

Independent Variables
The factors that potentially influence NO x concentrations are introduced as explanatory variables. Data for the 179 municipalities over the period 2000-2009 were taken from the ALMUDENA database of the Statistics Office of the Community of Madrid (See http://www.madrid.org/desvan/Inicio.icm?enlace=almudena accessed on 15 February 2021). Considering the specific situation in the Community of Madrid, we carried out the corresponding decomposition and improvements on the population (P) and technology (T) variables in order to understand which are the most important drivers affecting NO x emissions at a municipality scale, taking into account the geographical distribution of the driving forces in the Community of Madrid.
Bearing in mind the availability of data at municipality level, we specified Population (P) by considering total population, the proportion of the male population between 25 and 55 years of age, and the share of population over 65. Indeed, the environmental pressures of the population can vary among different age groups, which can reflect their degree of independence and potential contribution to economic activity (Liddle [38]; O'Neil et al. [39]). One would expect the economically active part of the population aged 14-64 to have a higher impact on emissions than pensioners over 65 or the age group encompassing children and adolescents under 14 (Cole and Neumayer [40]). However, York [41]) stated that there are reasons to expect that populations with a high proportion of elderly people consume more energy than those composed of younger people. For example, older populations tend to have smaller households and consume more energy intensive products, such as cars. Okada [42] found an inverted U-shaped relationship between CO 2 emissions per capita and the share of people over 65 years old in 25 OECD countries, which are mostly European countries.
The affluence variable (A), which indicates the level of prosperity of each municipality, is described through the gross disposable income (GDI) per capita at 1992 constant prices, electric power consumption (KWh per capita), and value added by the industrial sector (% of GDP). With the aim of testing the EKC hypothesis, the squared term of GDI per capita is also included to capture the possible nonlinear relationship between NO x concentrations and income. The presence of an EKC is indicated by a negative and significant coefficient for the quadratic term together with a positive and significant coefficient for GDI.
In the STIRPAT literature, the technology variable has been very difficult to define and convert into operational measures. We thus decompose it into various separate factors which influence environmental impacts, such as vehicle fleet, public transport infrastructure and services (number of lines of public transportation, number of train stations, number of bus stops and number of taxis).
The private vehicle population has increased dramatically in many municipalities in recent years. Indeed, the number of vehicles in Madrid has increased by 25. Since the regional government has not paid enough attention to the main source of emissions of this dangerous pollutant, the effect of the environmental plans implemented in recent years has been partially offset by the increasing use of road transport. As a result, the number of vehicles may have considerable impact on municipalities' NO x concentrations.
On the other hand, an increase in air pollutant emissions is explained not only by the growth of transport activity but also by the shift of transport to more energy-intensive modes and the increase in air pollutant emissions of modes in use. Therefore, in addition to the number of vehicles, the number of lines of public transportation, the number of train stations and the number of bus stops are considered to explain the environmental pressure caused by transportation (Dargay and Gately [43]; Eom and Schipper [44]). Table 2 lists the descriptive statistics of the explanatory variables analyzed as potential drivers of outdoor concentrations of NO x . Our database also contains a number of time-independent regressors, including membership of the Spanish Network of Cities for the Climate, which coordinates the actions carried out by the various administrations to tackle climate change; the transformation of landfills to green areas, which is aimed at mitigating unhealthy sources of contamination and the proliferation of infestations, in order to minimize the environmental impact on the surrounding area and to create new green zones; and others related to location (longitude, latitude, altitude, extension, distance to the city of Madrid, etc.). However, in order to not make Table 2 unnecessarily long, and because their impacts on the response variable cannot be estimated under fixed effects, they are not listed.

Spatio-Temporal Kriging
Kriging is a univariate procedure (which can be extended to cokriging, the multivariate version) which interpolates the values of the target variable at unobserved locations using the available observations of the same variable. This interpolation procedure, which is a minimum mean-squared-error method of spatial prediction, produces the best linear unbiased predictor. For this purpose, kriging uses the covariance or variogram function, which is the spatial (or spatio-temporal) equivalent of the autocorrelation function in timeseries analysis. The kriging strategy is based on the idea that variables follow a stochastic process over space or space-time. Other interpolation alternatives have been considered in the literature, such as Thiessen polygons, inverse distance method, splines, etc. However, kriging is more appropriate when dealing with environmental variables (Anselin and Le Gallo [45]) because it considers the spatial or spatio-temporal dependencies and provides a measure for the goodness of the prediction.
In our case, the spatio-temporal framework of an air pollution monitoring situation, let {Z(s 1 , t 1 ), . . . , Z(s n , t n )} be the values of NO x observed at a set of n spatio-temporal sites {(s 1 , t 1 ), . . . , (s n , t n )} equipped with an air pollution MS across a specific geographical area (in our case, the Community of Madrid). The philosophy of the spatio-temporal geostatistics paradigm when it comes to predicting the concentrations of NO x in a particular location that is not equipped with an MS at a specific instant of time, Z(s 0 , t 0 ), or across an area of interest (municipality, province, region, country) has the same logic as that used by economic agents. An economic agent would observe the concentrations at similar spatio-temporal points around (s 0 , t 0 ) and, more than likely, the prediction provided would be the mean of the observed concentrations. This procedure has an overwhelming logic; the only criticism we can make is that the weight assigned to the concentrations recorded at the sites similar to (s 0 , t 0 ) should not be the same. This is why the spatio-temporal predictor provided by geostatistics, the so-called spatio-temporal kriging (STK) predictor, is: where the weights, λ i , are obtained so that the resulting predictor is optimal, in the sense of unbiasedness and minimal variance of the prediction error; and this in a framework of non-independence of the variables representing the NO x concentrations at sites {(s 1 , t 1 ), . . . , (s n , t n )}, because the spatio-temporal dependence of air pollutant concentrations is already considered a stylized fact of air pollution data. The resulting STK equations depend on the degree of stationarity attributed to the random function that supposedly generates the observed realization. In the case of second-order or intrinsic stationarity, the usual case, the STK equations are as follows: with prediction (kriging) variance: where γ s i − s j , t i − t j is the semivariogram selected to capture the structure of the spatiotemporal dependencies existing in the phenomenon under study (in the second-order case, the covariogram C s i − s j , t i − t j can also be used) and α is a Lagrange multiplier associated with the non-bias condition (see Sherman [46]; Montero et al. [47], for details). In the stationary case-which is our case and the typical case in air pollution databases-the STK procedure detailed above is known as "ordinary" (spatio-temporal ordinary kriging, STOK). In the case of non-stationarity, the corresponding STK procedure is known as "universal" (spatio-temporal universal kriging, STUK). STUK equations for the non-stationary case can also be seen in Montero et al. [47].
The following points should be noted. First, as in the spatial-only case, STK strongly depends on the choice of the autocovariance or semivariogram associated with the spacetime random field under study. This is not a problem in spatial-only geostatistics. It could be said that there is a "standard" list of semivariograms that guarantee positive definiteness, covering almost all practical situations. However, in the spatio-temporal context defining such a list is a challenge, because the construction of non-separable spatio-temporal covariance or semivariogram structures (the space-time interaction is crucial for prediction) associated with stationary or non-stationary, isotropic or anisotropic random fields, which meet the requirements to be a valid covariance or semivariogram model, is not an easy task (it is not easy to prove the positive definiteness of a spatio-temporal dependence structure). Even the construction and visualization of an empirical space-time semivariogram can often be considered a challenging task (see example 6.2 in Montero et al. [47]). The original separable, isotropic and stationary covariance models (the metric model by Dimitrakopoulos and Luo [48]; the sum model introduced by Rouhani and Hall [49]; the product model by Rodriguez-Iturbe and Mejia [50], and De Cesare et al. [51]; the product-sum model popularized by De Iaco et al. [52,53]; and the integrated product and product-sum models also popularized by De Iaco et al. [54], among other initiatives) imply a number of restrictive conditions that makes them infeasible options for use in real applications. Therefore, more complex and realistic covariance models were needed to improve predictive performance. Thus, over the past decade a great effort has been made to deal with separability, stationarity, isotropy and full symmetry in the spatio-temporal structures capturing the spatio-temporal dependencies existing in the phenomena under study. Examples of such efforts include the models proposed by Cressie and Huang [55], Gneiting [56], Ma [57][58][59][60], Stein [61] and Porcu and Mateu [62] to take into account the space-time interaction; the specifications by Ma [63] and Gneiting et al. [64] to overcome full symmetry; the mixture-based Bernstein zonally anisotropic covariance functions by Porcu et al. [65] and the non-separable models constructed by using quasiarithmetic functionals (Porcu et al. [66], which go beyond isotropy; and the non-stationary models developed by Fuentes and Smith [67], Chen et al. [68], Ma [57,69] and Porcu and Mateu [62], among others. Given the long list of "recent" models cited above, one might think that, as in the spatial-only case, there is a "standard" list of permissible models; however, this is currently far from being the case. The new families of space-time covariance models are not exempt from problems and, in practice, the selection of one of them to represent the space-time dependencies existing in the phenomenon under study continues to be a challenge. Given the complexity of these "new" models, their visualization is recommended, as well as an analysis of how they evolve when the model parameters change; this is core information in the selection of the appropriate model to capture the spatio-temporal dependencies. However, visualization is not easy either, and constitutes another of the pending questions regarding spatio-temporal covariogram analysis (see Montero and Fernández-Avilés [70] for details and advances). By way of example, Huang and Sun [71] state that "it is challenging to visualize and assess separability and full symmetry from spatio-temporal observations".
As highlighted by Wikle [72] in his book review of Montero et al. [47], "to date, there are no general spatio-temporal covariance or semivariogram models that can adequately describe the complex dynamics that are exhibited by the multivariate, often nonlinear, spatio-temporal processes that govern most real-world processes".
Second, although it does not apply in our case, when dealing with spatio-temporal data, the solution of STK equations can involve considerable computational burden which might lead to failure in the prediction procedure. More specifically, the calculation of the inverse of covariance matrices becomes a crucial problem. This is known as the "big n problem" (Banerjee et al. [73]), which often arises when dealing with massive spatiotemporal datasets. Accordingly, recent literature emphasizes the use of approximation methods and new methodologies for dealing with massive spatio-temporal data sets.
Third, in case studies and practical applications, it has been found that STK prediction exhibits a very rapid reversion to the mean in locations near the border of the spatial region and/or close to the end of the time period considered (and, of course, in out-ofsample prediction).

Spatio-Temporal Kriging of the Mean
In geostatistics, the problem of estimating the mean value of a variable of interest in a specific area is called the kriging of the mean (KotM). We aim to estimate the unknown mean (in our case for NO x concentrations in the city of Madrid), in the years under study, µ t , t = 2000, . . . , 2009, using the estimator: The weights λ it are obtained so that the kriging prediction of the mean is unbiased and the prediction error variance is minimal. In order to ensure the unbiasedness of the kriging predictor, the sum of the weights must be one, as The resulting spatio-temporal equations, in terms of a permissible spatio-temporal covariance function are: which in the stationary case only depends on the distance between the spatio-temporal points. The estimation variance coincides with the Lagrange multiplier α.
It is worth noting that in some ways the KotM is the same as the "change of support problem" (COSP). Kriging is very often the solution to overcome this mismatch of spatial support (Gotway and Young [74]), particularly when dealing with socio-economic data (Census data, pollutant concentrations, housing prices, etc.) since it takes into account spatial or spatio-temporal dependencies. When dealing with air pollution data, the usual solution to the abovementioned problem is to krige the environmental variable(s) to obtain their interpolated values in the locations where socio-economic data are available. Two kriging alternatives can be used: (i) the only-spatial KotM for every instant of time considered in the database, which means the estimation of a semivariogram for each instant of time; or (ii) the STKotM, where only one spatio-temporal semivariogram needs to be estimated and used whatever the period for which the mean is estimated. The good performance of STK in the in-sample case ("sandwich predictions"), is the reason why we opt for the second alternative.

A Spatial Formulation of the STIRPAT Model
As stated in Laureti et al. [8], when analyzing environmental pressure, the most widely-used research frameworks addressing different driving forces linking economic features with environmental issues, are the STIRPAT model (York et al. [75]) and the EKC approach (Andreoni and Levinson [76]). We also refer to Laureti et al. [8] and the literature therein, for a review of the literature on the STIRPAT approach to the study of the drivers of environmental impacts.
According to Dietz and Rosa [77,78] the basic version of the STIRPAT model, which evaluates various theoretical frameworks on the systemic drivers of environmental change, is given by: where I represents the environmental impact, P is population, A is consumption per capita (or affluence), T is technology (or impact per unit of consumption), e is the residual or error term and the subscript i denotes cross-sectional units. β 0 , β 1 , β 2 and β 3 are the model coefficients or elasticities.
The STIRPAT methodology has mainly been applied at a large scale, using crosssectional and panel national-level data (see Liddle [79], for a review until 2013; and other papers published in the last five years, such as Wang et al. [80], Al-Mulali et al. [81], Abdallh and Abugamos [82], Koçak and Ulucak [83], Lohwasse et al. [84], and Velez-Henao et al. [85] and the references therein for an updated review).
In the last decade, the STIRPAT approach has also been applied at a "local" scale (county, city and province level). Yin et al. [86], Xu and Lin [87] and Xu et al. [88] in China; Squalli [89]; Roberts [90] in the USA; Scholz [91] in Japan; and Anser [92] in Pakistan are some examples. However, only Laureti et al. [8], in the Community of Madrid, Spain, is an application of STIRPAT at genuinely local scale.
In recent years, the spatial argument has started to play a key role in the framework of environmental changes (see Atwi et al. [9] for a review of the literature for and against controlling for spatial effects in environmental equations). The reasons behind this are the ever-growing popularity of the new developments in spatial econometrics and the awareness that the environment and many of its explanatory factors have a clear spatial nature. However, much of this literature focuses on the EKC hypothesis. Some examples include Rupasingha et al. [21], which estimated Bayesian Tobit SE models to determine the drivers of toxic pollutants in US counties, while Maddison [22] used a panel of SO 2 , NO x , CO and volatile organic compounds emissions in 135 countries to explore different spatial econometric techniques in estimating EKCs. According to Maddison [22], in EKCs, spatial relationships arise most obviously as a consequence of countries' strategic response to transboundary pollution flows. However, the analysis is a country-level analysis and the two spatial models used (the spatial lag model and the spatial error model) lose much of their utility. Something similar can be said of Poon et al. [93] for SO 2 in China, and for Burnett et al. [94] in the USA, at a state scale, for the period 1970-2009, when analyzing the relationship between CO 2 emissions, energy consumption, economic activity, and pollution emissions while controlling for potential spatial effects within the data, in a panel data framework. At any rate, all of these studies conclude that the EKC equation should control for spatial impacts.
Wang et al. [23] used a spatial econometric approach (spatial lag, spatial error and Durbin models) and found evidence supporting an inverted U-shape EKC between economic growth and ecological footprint.
Pattison et al. [31], in an exploratory study at county level in the USA for the year 2002, and using a spatial regression analysis, tested the hypotheses that affluence is positively related to carbon emissions from consumption activities but negatively related to emissions from production activities.
Zhao et al. [24] used spatial panel data models to investigate the factors influencing energy-related CO 2 emission intensity among a panel of 30 provinces in China covering the period 1991-2010. Hao and Liu [26] analyzed the socio-economic factors influencing urban PM 2.5 concentrations in 73 Chinese cities in 2013 by using Spatial Lag and Spatial Error (SL and SE, respectively) models. Kang et al. [25] examined the CO 2 EKC hypothesis using a spatial panel data model covering the period 1997-2012 for 30 Chinese provinces to avoid the coefficient estimation error. Atwi et al. [9] studied the case of CO 2 , by comparing the estimates derived from the cross-sectional and the panel data approaches, with data from 182 countries during the period 1992-2011. They found that the EKC hypothesis is acceptable under both approaches, although the estimated turning points in cross-sections seem unreliable.
Ding et al. [95] examined the existence of EKC for PM 2.5 pollution in the Beijing-Tianjin-Hebei region of China. They utilized satellite observations of PM 2.5 and panel data of 13 cities from 1998 to 2016 to investigate the relationship between economic growth and PM 2.5 pollution. Having tested the existence of spatial effects, a Spatial Durbin (SD) model was chosen for the above purpose although SL and SE were also estimated. Huang et al. [96], Yu et al. [97] and Feng et al. [98], among many others, conducted similar research.
In the light of the above literature, it is clear that research using STIRPAT models and the EKC equation has evolved by including the corresponding terms to control for spatial effects, as their existence has been sufficiently proved in the literature.
As can be noted, most of the literature using the STIRPAT approach accounting for spatial effects relies on US counties and Chinese cities or provinces as spatial support. It is true that this spatial support can be considered "local" compared with countries. However, most counties in the USA (and, of course, cities or provinces in China) cover more area than the most important cities in Europe. As such, spatial dependencies make more sense when dealing with "municipalities" in Europe than when working with counties, cities or provinces as spatial units in the USA or China. By way of example, the average area of a US county is 2911.4 km 2 , the average area of a Chinese province is 291,406 km 2 and the average area of a municipality in Madrid, the region under study in this article, is 45.36 km 2 . In addition, when considering spatial effects in the estimation of the impacts of a number of covariates on the environment, and especially on air pollution, the population density and the distribution of that population across the spatial unit considered is a key factor; the higher the population density and the more uniformly it is distributed across the territory, the more important it is to consider spatial effects in the STIRPAT model. This is the rationale behind conducting the analysis at a genuinely local scale and carrying out our practical study in a territory in Europe.
We want to make it clear that in local STIRPAT modeling, when dealing with air pollution concentrations, spatial dependencies must be taken into account because in reality a large percentage of local spatial units are not MS-equipped. Consequently, in our case study, where spatial units are the municipalities of the Autonomous Community of Madrid (clearly local units), we are forced to specify a spatial STIRPAT specification because the dependent variable is an interpolated variable (in the sites where pollution concentrations cannot be measured).
Spatial relationships can be modeled in a variety of ways depending on the relationship between the dependent variable and the explanatory variables. Following LeSage and Pace [99], we consider three spatial panel data models; namely, the spatial autoregressive panel data (SL-PD) model, the spatial Durbin panel data (SD-PD) model and the spatial error panel data (SE-PD) model. They allow us to examine three different types of interaction effects among municipalities (geographical units), that is, endogenous interaction effects among the response variable (Y), exogenous interaction effects among the independent variables (X), and interaction effects among the error terms (ε it ).
The SL-PD model hypothesizes that the value of the dependent variable observed at a particular location is partially determined by a spatially weighted average of neighboring dependent variables (the first way of specifying interaction between spatial units). That is, NO x emissions in municipality i are affected by the NO x emissions in municipality j.
The SL-PD model is expressed as: where y it denotes the dependent variable, I it , (NO x emissions) for unit i at time t. i n is a (n x1) unit vector for the intercept (we have opted to remove the unitary column from X to avoid problems of exact multicollinearity in the estimation process of models including spatially lagged covariates), ρ denotes a scalar spatial autocorrelation coefficient on the dependent variables. The term N ∑ j=1 w ij y jt denotes the interaction effect of the dependent variable y it with the dependent variables y it in neighboring municipalities, where w ij is the ij-th element of a pre-specified nonnegative spatial weighting matrix W of order N 2 describing the spatial arrangement of the units in the sample. ζ t captures the time-period effects, u i captures all time-constant unobserved municipality-specific factors, such as the concern of the population regarding the quality of the environment, which could be viewed as being roughly constant over the period in question. ε it is the disturbance term, which is i.i.d. across i and t with zero mean and variance σ 2 , both in fixed and random effects models (see for example, Baltagi and Liu [100]).
In the SD-PD strategy, the dependent variable can be predicted as a function of spatially lagged values of the explanatory variables as well: where θ is a (K × 1) vector of spatial autocorrelation coefficients on the lagged explanatory variables. If a spatial process of the disturbance (error) term is specified (the second way of specifying the interaction between spatial units), the resulting specification is an SE-PD model, which refers to a situation in which the unobserved shock to municipality i is affected by unobserved shocks in neighboring municipalities: where φ it reflects the spatially autocorrelated error term and λ is called the spatial autocorrelation coefficient.
It should be noted that: (i) both SL-PD and SE-PD models are nested in the SD-PD model; and (ii) that there are spatial specifications more general than the SD-PD model, such as the spatial Cliff-Ord type autoregressive panel data model, also known as the spatial autoregressive panel data model with Autoregressive Disturbances (SARAR-PD model), which was introduced by Kelejian and Prucha [101]. However, the present study focuses on spillovers and, according to LeSage and Page [100], SARAR-PD specifications concentrate on a more elaborate model for disturbances. In addition, SL-PD and SD-PD strategies, which are nested in SARAR-PD models, provide consistent estimates for the majority of spatially correlated data-generating processes (Minguez et al. [102]).
By referring to a panel data context, and expressing Equation (6) in logarithms so that it becomes additive, the specific a-spatial panel data STIRPAT (PD-STIRPAT) model including the specific response and covariates described in Section 2, can be specified as: ln NOx it = α + β 1 ln(POP it ) + β 2 ln(MALE25 − 55 it ) + β 3 ln(POP65+ it ) + β 4 ln(GDI it ) +β 5 (ln(GDI it )) 2 + β 6 ln(I NDUSTR it ) + β 7 ln(ENERGY it ) + β 8 ln(FLEET it ) +β 9 ln(HEAVY it ) + β 10 ln(PUBLIC it ) + β 11 ln(BUSSTOP it ) + β 12 ln(TAXI it ) +β 13 PERIOD t + u i + ε it (10) where the subscripts i and t denote municipalities and years, respectively. As mentioned previously, ε it is the disturbance term, u i captures all time-constant unobserved municipality-specific factors, such as the concern of the population regarding the quality of the environment, which could be viewed as being roughly constant over the period in question.
As outlined above, an additional term ρ N ∑ j=1 w ij y jt leads to the SL-PD-STIRPAT specification (with ρ being the spatial autoregressive coefficient); and another additional term N ∑ j=1 w ij x ijt θ leads to the SD-PD-STIRPAT model. However, the inclusion of the spatially autocorrelated error term φ it = λ N ∑ j=1 w ij φ it + ε it would specify the SE-PD-STIRPAT strategy.
It is important to note that in SL-and SD-PD-STIRPAT models, the spatial impacts are not given directly by any vector of parameters. However, the expressions of direct and indirect effects, as well as of the averaged direct, indirect, and total impacts can be seen in Table 3. (LeSage and Page [100]). Table 3. Total, direct, and indirect impacts on the response variable resulting from a change in covariates: Spatial-PD-STIRPAT models.

Model
Direct Impact Indirect Impact Total Impact SL tr (I n −ρW) −1 I n β k n Total-Direct With the aim of testing the EKC hypothesis, the squared term of GDI per capita is also included to capture the possible nonlinear relationship between NO x concentrations and income.
Given that (i) the assumption of zero correlation between the random effects and the explanatory variables is particularly restrictive, and that (ii) in spatial research it is at least controversial that the number of units is potentially able to go to infinity and that the units of observation are representative of a larger population, which are three conditions that should be satisfied before the random effects model can be implemented (Elhorst [103]), we opt for fixed effects modeling. Some researchers implement Lagrange multiplier tests to inform the above decision. However, these tests have proven to have high power for a large number of data-generating processes in the alternative hypothesis (anyway, the spatial extension of the Haussman test (Mutl and Pfaffermayr [104]) indicates that the data do not support the random effects assumption (chi-sq = 649.4, df = 21, p-value < 2.2 × 10 −16 ; chisq = 456.73, df = 21, p-value < 2.2 × 10 −16 ; and chisq = 5311.6, df = 23, p-value < 2.2 × 10 −16 for the SL-, SE-and SE-PD-STIRPAT models, respectively).
The decision to estimate fixed effects models implies that some location-related covariates cannot be included in the model because they are time-invariant.
As for which lagged drivers of NO x to include in the SD-PD-STIRPAT model, the corresponding nested likelihood ratio tests were implemented.
The estimation method used was Maximum Likelihood. Two-Stage Least Squares estimation could have been implemented if it were suspected that some of the regressors were endogenous, but that is not the case in our study. Therefore, there was no need to use instrumental variables.
Finally, all computations were conducted using the libraries splm (Millo and Piras [105]), spdep (Bivand and Piras [106]) and maptools (Bivand and Lewin-Koh, [107]) in R (R Development Core Team [108]).  Table 1 in Section 2.1 summarizes the main descriptive statistics of their values. The spatio-temporal covariance function used in the prediction process, both in the STK and STKotM equations, was the non-separable spatio-temporal covariance function proposed in Gneiting [56]:

Spatio-Temporal Kriging Results for NO x Predicted Values
where a and c are two positive space and time scale parameters, respectively, and σ 2 is the a priori variance of the NO x stochastic process; δ and α are the space and time smoothing parameters, respectively, with values in the interval (0, 1]; and β is the space-time interaction parameter, which moves in the interval [0, 1]. The larger the β, the more intense the spatiotemporal dependence degree in the stochastic process under study. h and τ are the spatial distance and the time lag, respectively, and θ represents the set of parameters involved in the covariance function. were endogenous, but that is not the case in our study. Therefore, there was no need to use instrumental variables. Finally, all computations were conducted using the libraries splm (Millo and Piras [105]), spdep (Bivand and Piras [106]) and maptools (Bivand and Lewin-Koh, [107] in R (R Development Core Team [108]).  Table 1 in Section 2.1 summarizes the main descriptive statistics of their values. The spatio-temporal covariance function used in the prediction process, both in the STK and STKotM equations, was the non-separable spatio-temporal covariance function proposed in Gneiting [56]:

Spatio-Temporal Kriging Results for NOx Predicted Values
where a and c are two positive space and time scale parameters, respectively, and 2  is the a priori variance of the NOx stochastic process;  and  are the space and time smoothing parameters, respectively, with values in the interval (0, 1]; and  is the spacetime interaction parameter, which moves in the interval [0, 1]. The larger the  , the more intense the spatio-temporal dependence degree in the stochastic process under study. h and  are the spatial distance and the time lag, respectively, and θ represents the set of parameters involved in the covariance function. The reasons behind this decision are (i) the shape of the NOx empirical spatio-temporal covariogram; (ii) it is a non-separable function; (iii) it is sufficiently general for use with air pollution concentrations series. Table 4 lists the values of Moran's I for the period under study, which confirm the existence of spatial autocorrelation in NOx concentrations also from a spatial econometrics perspective.  The reasons behind this decision are (i) the shape of the NO x empirical spatio-temporal covariogram; (ii) it is a non-separable function; (iii) it is sufficiently general for use with air pollution concentrations series. Table 4 lists the values of Moran's I for the period under study, which confirm the existence of spatial autocorrelation in NO x concentrations also from a spatial econometrics perspective.

Coefficients, Direct Impacts and Spillovers
The first finding obtained from the results commented in Section 3.1 and those listed in Table 5 is that space matters and that the model used to specify the spatial PD-STIRPAT representation of the NO x concentrations in the Autonomous Community of Madrid is far from irrelevant. Our results indicate that, although the spatial correlation coefficient in the SE-PD-STIRPAT model, λ, is not significant, space matters. This is shown by the fact that (i) in the SL-and SD-PD STIRPAT models ρ has a significant value (at the 10% level), although it is very moderate; and (ii) some lagged explanatory factors such as W*PUBLIC and W*BUSSTOP show significant (and notable) spatial autocorrelation coefficients at the 1% and 5% level of significance, respectively. It is no surprise that λ is not significant, because the set of factors we have used as explanatory variables is among the best used in the literature. Consequently, as advanced in the introductory section, an SD-PD-STIRPAT model at local scale, which takes into account the effects of spatial dependencies inherent in air pollution concentrations and some of the factors that cause them, is a perfect candidate for estimating the stochastic impacts on NO x of Population, Affluence, and Technology in the Autonomous Community of Madrid. Table 5. Estimates of the spatial autocorrelation parameters for a-spatial and spatial-PD-STIRPAT models. Further exploring spatial effects, the individual effects, which allow for the inclusion of a certain degree of time-invariant unobserved spatial heterogeneity in the model, are shown in Figure 4. Since they capture some idiosyncratic characteristics of the municipalities (including the time-invariant variables longitude, latitude, altitude, extension, distance to the city center, etc.), whose impacts on the response variable are beyond doubt but cannot be estimated when dealing with fixed effects estimation, they might account for some part of the autocorrelation effects.

SL
The first finding obtained from the results commented in Section 3.1 and those listed in Table 5 is that space matters and that the model used to specify the spatial PD-STIRPAT representation of the NOx concentrations in the Autonomous Community of Madrid is far from irrelevant. Our results indicate that, although the spatial correlation coefficient in the SE-PD-STIRPAT model, λ, is not significant, space matters. This is shown by the fact that (i) in the SL-and SD-PD STIRPAT models ρ has a significant value (at the 10% level), although it is very moderate; and (ii) some lagged explanatory factors such as W*PUBLIC and W*BUSSTOP show significant (and notable) spatial autocorrelation coefficients at the 1% and 5% level of significance, respectively. It is no surprise that λ is not significant, because the set of factors we have used as explanatory variables is among the best used in the literature. Consequently, as advanced in the introductory section, an SD-PD-STIRPAT model at local scale, which takes into account the effects of spatial dependencies inherent in air pollution concentrations and some of the factors that cause them, is a perfect candidate for estimating the stochastic impacts on NOx of Population, Affluence, and Technology in the Autonomous Community of Madrid. Table 5. Estimates of the spatial autocorrelation parameters for a-spatial and spatial-PD-STIRPAT models. 0448] ** *, **, ***: significant at the 10%, 5% and 1% significance level, respectively.

SL
Further exploring spatial effects, the individual effects, which allow for the inclusion of a certain degree of time-invariant unobserved spatial heterogeneity in the model, are shown in Figure 4. Since they capture some idiosyncratic characteristics of the municipalities (including the time-invariant variables longitude, latitude, altitude, extension, distance to the city center, etc.), whose impacts on the response variable are beyond doubt but cannot be estimated when dealing with fixed effects estimation, they might account for some part of the autocorrelation effects. The first finding detailed above (an SD-PD-STIRPAT model at local scale is a perfect candidate to estimate the stochastic impacts of P, A and T on NOx) is consistent with the The first finding detailed above (an SD-PD-STIRPAT model at local scale is a perfect candidate to estimate the stochastic impacts of P, A and T on NO x ) is consistent with the second one: the SD-PD-STIRPAT model is the spatial specification that best captures and fits the spatial dependence within the relationship between transport-related NO x concentrations and socio-economic drivers at municipal level, where the relationship between socio-economic factors and the environment becomes more complex. As can be seen in Table 6, the three competing spatial panel data models provide similar results for fixed effects estimation. According to the Akaike and Bayesian Information criteria (AIC and BIC, respectively), the best model is SD-PD-STIRPAT. In addition, the LR test-statistic (277.6051, df = 2, p-value < 2.2 × 10 −16 ) to test whether the SD-PD-STIRPAT model can be simplified to SL-PD-STIRPAT has been rejected. It is also observed that the R 2 of the spatial-PD-STIRPAT strategies is slightly larger than the R 2 of the a-spatial one. This is an expected result because, as stated in Elhorst ([109], p. 9), "researchers frequently find a weak evidence in favor of the spatial interaction effects when the fixed time effects are considered in the model". Table 6. Goodness-of-fit measures for a-spatial and spatial-PD-STIRPAT models. Finally, the last reason to select the SD-PD-STIRPAT model is that the focus of the article is on both local and global effects.

A-Spatial
Third, as can be seen in Table 7, only one of the three covariates representing "Population", POP65+, has significant stochastic impacts on the level of NO x concentrations in the region of Madrid, at the significance level of 1% in the SL model and 5% in the SD strategy. The negative relationship of this Population covariate with NO x concentrations is an expected result, because in the Autonomous Community of Madrid, like in the other Spanish Autonomous Communities, the population living in municipalities in the countryside (also called "empty Spain") or near to the mountain is the oldest. The "empty Spain" factor also explains why in the a-spatial PD-STIRPAT model both POP and MALE25-65 are highly significant. Table 7. Estimates of total impacts: a-spatial model and spatial-PD-STIRPAT models.

COVARIATES
A-Spatial SL SE SD green areas, and the industrial ones are in south and east, where land is cheaper-GDI has a negative impact whereas the impact of ENERGY is positive. This finding is reinforced by the fact that the impact of INDUST on NO x concentrations is positive and significant at the 10% level of significance. It is of note that only ENERGY is significant when spatial dependencies are not considered in the panel data model; in addition its impact on NO x is less than half the impact estimated when space is included in the PD-STIRPAT strategy. Finally, regarding "Technology", all the covariates included under this element turned out to be significant (with the exception of BUSSTOP in the SD model). The one with the greatest (negative) impact is the total number of lines of public transportation per capita (PUBLIC), regardless of the spatial model used. The second-ranking covariate is HEAVY, with a positive impact, as expected. Third is FLEET, with similar positive impacts to HEAVY. It is worth noting that in the a-spatial PD-STIRPAT model, the most influential covariate is FLEET, whereas HEAVY has the least impact. BUSSTOP and TAXI have similar impacts on the level of NO x concentrations, although, evidently, with opposite sign. Surprisingly, in the SE model, TAXI has a negative sign.
Time effects are mostly significant and notable, and in line with the evolution of the air pollutant under study over the time period analyzed. Interestingly, there are no singular differences between the estimates provided by the a-spatial and the spatial PD-STIRPAT strategy.
Given the moderate values of the rho coefficient in both SL-and SD-PD-STIRPAT models, the total effects are mostly direct effects (Table 8). However, POP65+, GDI, ENERGY and PUBLIC indirect effects or spillovers are significant in both SL-and SD-PD-STIRPAT models. FLEET and HEAVY indirect effects are significant in the SL model, and the two lagged variables, W*PUBLIC and W*BUSSTOP, indirectly influence NO x concentrations in the SD strategy. Consequently, a change in the previous non-lagged and lagged covariates in a particular municipality has an average cumulative effect on the corresponding variables in neighboring municipalities. Shaded cells in the indirect impacts indicate significance at the 10% significance level. Direct effects can be obtained by subtracting the indirect effects from the total ones. The level of significance for direct effects is the same as for total effects (Table 7). Table A1 in the Appendix A lists the coefficients for the a-spatial and spatial PD-STIRPAT models considered in our analysis. Tables A2-A4 show the results obtained when estimating an RE model.

Testing the EKC
As can be seen in Table 6, GDI displays a negative sign and GDI (squared) a positive one, irrespective of the spatial-PD-STIRPAT strategy considered, indicating a U-shaped curve instead of the traditional inverted-U.
The result obtained is not in line with the EKC introduced by Grossman and Krueger [110]; that is, an inverted U-shaped relationship between income and pollutant emissions, which assumes that environmental quality worsens in the early stages of development but improves after a certain level of income is reached (see Tzeremes [111,112] for details). According to the empirical evidence, it is indeed shown that for some high-income economies and for certain pollutants the relationship between per capita income and pollution follows an inverted U shape (Falconi et al. [113]). However, as stated in Falconi et al. [113], it must be taken into account that the shape of EKC depends on the pollutant under analysis, the region and the period being studied (to which we would add the length of the period under study). It has been empirically proven that there is an inverted U-shaped ECK for local pollutants. However, for global pollutants, emissions increase with economic growth. In our case study, our database refers only to the period 2000-2009 and to a specific (small) region of a specific country. Using these data, we find a U-shaped relationship, which is probably the consequence of the environmental measures promoted by the Autonomous Community of Madrid at the turn of this century and the spatial distribution of income across the region. In this respect, it is noteworthy that when spatial dependencies are not considered, the shape of the EKC is the traditional inverted U-shaped. Thus, space and spatial dependencies matter (in this case in the short term, at a local scale) when dealing with the EKC.

Policy Implications
The above results yield important economic implications, which are of interest to policy-makers. The significance of the spatial autocorrelation coefficient ρ is a finding of vital importance for the success of political decisions on environmental issues. As outlined in the introductory section, the existence of spatial spillovers indicates that the level of NO x concentrations in a municipality affects and is affected by the level of NO x in neighboring municipalities, which in turn means that environmental governance of a municipality impacts and is impacted by environmental governance in neighboring municipalities. In addition, the existence of significant and notable spatial effects on NO x derived from some lagged explanatory covariates or drivers (we are referring to the theta coefficients), indicates to environmental policy-makers in the region of Madrid that NO x -related policies implemented in a municipality influence the results of policies in neighboring municipalities. From this, we can infer that the various municipalities should work together to balance pro-environmental and economic policy goals. This is especially important in geographical areas composed of highly autonomous spatial units, as is the case of municipalities in the region of Madrid. Therefore, in light of this result, convincing policy makers to act in a coordinated way is of crucial importance. In fact, the lack of coordination among mayors belonging to political parties of different ideologies plays a role in preventing further progress on air pollution control. It is true that, in the period under study, the magnitude of the rho coefficient in the region of Madrid is not very high; however, environmental policy-makers must pay special attention to it because it may increase (or decrease) in a short period of time.
Elaborating on the case of the spatial effects produced by the lagged drivers of NO x (theta coefficients), the "Total number of lines of public transportation per capita (PUBLIC)" and the "Number of bus stops per capita (BUSSTOP)" are the only covariates producing spatial effects on NO x derived from their spatial lags in the region of Madrid for the period under study. This suggests that municipal actions aimed at increasing both the number of buses and bus stops contribute to clean air not only in the implementing municipality but also in the neighboring municipalities. Consequently, in the presence of this type of effect, municipal actions aimed at improving air quality are not only the concern of the municipality that carries them out, but are also a matter of solidarity with neighboring municipalities. This applies especially in the case of the number of bus stops per capita, because its total impact on NO x is not significant for the non-lagged variable, but is for its lagged version.
The above results stand in contrast to the fact that none of the private sector-related covariates has significant lagged effects on the NO x concentrations of neighboring municipalities, which reinforces the need for cooperation at a municipal level in public transportation activities. Undoubtedly, this municipal cooperation will be of great help in significantly reducing the concentrations levels of the only air pollutant still out of control in some areas of the region of Madrid.
An additional effort must be made to reduce the negative impact of some technologyrelated drivers of NO x in the region of Madrid. We are referring to the "equivalent fleet of private vehicles per capita" and the "number of heavy vehicles per capita". Both have a positive and significant impact on the level of NO x concentrations. As such, it is our opinion that efforts should be made to promote green driving in general and electric vehicles in particular. Electric vehicles are expected to significantly reduce road transport emissions, given an increasingly renewable power generation. While outstanding technological issues are continually being solved, the economic viability of a greener driving is still constrained by the large gap currently existing between conventional and green vehicles. Consequently, municipalities should strive to democratize green driving by promoting green transport technology (especially for electric vehicles), subsidizing their purchase and facilitating charging options (especially for private households), while promoting other clean transportation alternatives (especially in urban environments). This is a good moment for the municipal authorities of the region of Madrid (and also for regional and state authorities) to put the above measures into action, because they could play a key role in tackling the crisis currently afflicting the automotive sector in particular, and the economic crisis in general. It is more than likely that, in the future, when conventional vehicles have been replaced by electrical and other green vehicles, the covariate "vehicles fleet" will no longer impact the level of NOx concentrations.
As for EKC-related recommendations to municipal policy makers of the region of Madrid, again, spatial effects matter and the shape of the NO x concentrations-income relationship in Madrid is an issue of concern for all the municipalities. Accordingly, they should all collaborate in addressing it to achieve the desired shape.
Finally, in the current framework in which both the private and public sectors strongly support compliance with the principles of Environmental, Social, and Corporate Governance (ESG), having information on the direct and indirect effects of the main drivers of NO x pollution on the magnitude of its concentrations greatly facilitates governance oriented towards balancing economic activity and clean air. We entirely agree with Uchiyama [114] and Ozcan et al. [115] that the path towards sustainable development can be suggested as optimal policy for (in our case) municipalities, emphasizing the coexistence rather the trade-off between economic growth and the environment. Therefore, the main target of the policies associated with the path of sustainable growth is the creation of a working framework unifying economic development and the environment.
In brief, the spatial panel data Durbin STIRPAT model we contribute for estimating the direct and indirect impacts on NO x concentrations of its main drivers in the region of Madrid, at a municipal level, constitutes a valuable instrument for policy making oriented towards controlling this dangerous pollutant. Perhaps the main lesson to be learned from the estimates obtained is that reducing NO x concentrations is a collaborative task, because environmental measures taken by a municipality of the region affect the concentration levels of the pollutant in neighboring municipalities.

Conclusions
In recent decades, the STIRPAT framework has been widely adopted to explore the environment-development nexus at macro level, involving cross-national data, and to a lesser extent at local scale, finding evidence on the way that different socio-economic factors including population, income and technology influence air pollution concentrations. However, the majority of researchers have addressed this relevant research question (i) without considering the spatial dependencies existing in air pollution concentrations and also in the drivers of said concentrations, and (ii) based on large spatial units due to the difficulties in finding data on the drivers of air pollution concentrations at a true local level. These two facts may well be core issues when it comes to explaining the lack of robust results to date.
The first main conclusion obtained in this article is that space matters; consequently, we observe a significant change in the estimates of the impact of the drivers of NO x concentrations on the level of such concentrations when the spatial dependencies in both NO x and some of its drivers are considered in the STIRPAT model. In addition, the decomposition of the total impacts into direct and indirect (spillovers) impacts is of vital importance for policy makers. This is especially true in the region of Madrid, where the level of NO x concentrations in a few areas remains a pending air pollution issue. Knowledge about the magnitude of such impacts is core information for the environmental decision-makers to determine the magnitude of the environmental measures to be taken, while also ensuring the objectives are achieved not only in their own municipality, province, etc., but also in the neighboring municipalities, provinces, etc. Indeed, as this study has shown, their decisions on NO x pollution drivers also affect other neighboring areas; and vice versa-NO x concentration levels of an area are affected by decisions taken by the policy-makers of neighboring areas. More specifically, knowing the size of the spillovers allows for joint agreements and policies with neighboring areas (municipalities in the case of the region of Madrid) to export the benefits of environmental measures taken in an area (municipality) to neighboring areas (municipalities) and limit negative externalities. Therefore, the existence of significant spatial dependencies in air pollution concentrations suggests we should put an end to autonomous governance of this environmental issue and move towards collaborative governance.
The second main conclusion is that the spatial STIRPAT models used to estimate the relationship of population, income and technology with air pollution concentrations and the magnitude of their impacts (direct and indirect) must be estimated at a local scale. The reason reflects Tobler's first law of geography ("everything is related to everything else, but near things are more related than distant things"), which, according to the literature on air pollution, applies in the field of atmospheric pollution. In addition, Tobler's postulate has been supported by the vast literature on the topic using spatial statistics, especially geostatistics and spatial econometrics. However, the range of significant spatial dependencies, which can be estimated with a semivariogram, rarely exceeds a few kilometers, suggesting that spatial STIRPAT models should not be estimated for large spatial units (regions, countries, etc.). Obviously, this calls into question the results obtained from the literature estimating STIRPAT models at a large scale.
More specifically, in this article we have estimated (using fixed effects estimation given the nature of the problem under study) the a-spatial Durbin panel data STIRPAT model together with the spatial lag, error and Durbin panel data STIRPAT models. As noted above, the main conclusion drawn in this article is that space matters when estimating the stochastic impacts of Population, Affluence, and Technology on NO x in the municipalities of the Autonomous Community of Madrid. Accordingly, the models used to estimate such relationships must be able to capture the spatial dependencies existing in both the response variable and some of the covariates. In line with theoretical premises, the spatial panel data STIRPAT strategy that best captures and represents the spatial dependencies in the NO x concentrations and the drivers of these concentration levels in the region of Madrid is the Durbin model. We cannot state that we contribute a new spatial STIRPAT econometric model capably of capturing spatial dependencies, because the panel data version of the spatial lag, spatial error or spatial Durbin STIRPAT models have been already used in the literature. However, in this article, we contribute to filling the "true local scale" gap in the literature by exploring the issue of spatial autocorrelation inherent in air pollution at short-to-medium distances within the STIRPAT framework and at a local scale.
Dealing with data at a local scale usually involves a lack of information, especially when using air pollution variables as response variables, because there are fairly few atmospheric pollution MS, even in very big cities. In this article we have solved this problem using spatio-temporal kriging techniques, which provide very high-quality estimates. In addition, in spite of the limitations in obtaining information on the drivers of NO x at a very local scale, this article contributes one of the best databases (in fact a proprietary database) feeding STIRPAT models with air pollution concentrations in the response variable.
Although space matters, the magnitude of the spatial autocorrelation coefficient is moderate, both in the spatial lag and spatial Durbin strategies, which can be attributed to the fact that some time-invariant covariates representative of the unobserved spatial heterogeneity have not been included in the models given their fixed effects nature.
It has been shown that POP65+, all the "Affluence" variables (GDI, ENERGY and INDUST; especially GDI and ENERGY) and all covariates coming under "Technology" (with the exception of BUSSTOP) have significant impacts (with the expected sign) on the level of NO x concentrations under an SD structure. These total impacts are mostly direct. In addition, POP65+, GDI, ENERGY and PUBLIC as well as the two lagged variables, W*PUBLIC and W*BUSSTOP, have significant spillovers. In other words, a change in them in a particular municipality has an average cumulative effect on the corresponding variables in neighboring municipalities. Time effects are mostly significant and notable, and in line with the evolution of this air pollutant over the period under study.
In addition, we have included the squared term of GDI in the model in order to test the EKC. Our finding is that when spatial dependencies are not considered, the shape of the EKC is the traditional inverted-U; however, we observe a non-inverted-U when space is accounted for in the model. Thus, space and spatial dependencies matter (in this case in the short term, at a local scale) when dealing with the EKC. Nevertheless, it must be borne in mind that our case study refers to a recent and short period of time in a specific region of a specific country and this finding should be corroborated by further, similar local-scale case studies. In the meantime, we support the statement by Falconi et al. [113]) that the shape of the EKC depends on the pollutant under analysis, the region and the period being studied, especially for short periods of time and at a local scale.
In brief, from a practical perspective, we contribute a spatial panel data Durbin STIRPAT model which can be of great interest for the environmental authorities and policymakers of the Autonomous Community of Madrid, an important Spanish region with levels of NO x concentrations which frequently exceed the legal limits, and whose citizens are deeply concerned about the issue of air pollution. As far as we know, this is the first spatial STIRPAT model estimated for the region of Madrid. In addition, we do not know of any other papers using spatial STIRPAT modeling at a similar local scale to the one used in this article. Consequently, we cannot compare the results obtained for the region of Madrid with results obtained for other areas at local scale. An alternative option would be to compare our results with those obtained in studies that estimate spatial STIRPAT models on the basis of large spatial units (provinces in China, counties in the USA, countries, etc.). However, according to Tobler's law, this comparison makes no sense from a spatial econometric perspective. We strongly encourage other researchers to apply these strategies at a local scale in other regions with high air pollution concentrations, in order to see whether their findings corroborate or contradict our findings.
The spatial panel data Durbin STIRPAT model we contribute, which can be estimated dynamically, is an instrument of particular importance in the current environmental frame-work shaped by the ESG principles and the United Nations 2030 Agenda, as it greatly facilitates the environmental governance demanded by the citizens of modern societies.
Finally, we would be interested in future theoretical research on semiparametric spatial panel data STIRPAT modeling, especially using penalized splines for smoothing purposes. This approach would not only account for spatial autocorrelation and spatial heterogeneity but also the non-linear (smooth) relationships between certain drivers of the response variable and the response variable. We are also interested in seeing the application of functional data analysis to the detection of air pollution episodes and outliers (see Martinez Torres et al. [116] for a recent application in the city of Dublin, Ireland).

Acknowledgments:
The authors are grateful to Roman Minguez, (University of Castilla-La Mancha, Spain) for his valuable comments and his support and advice in the programming phase. We also acknowledge the four anonymous referees for their valuable comments and suggestions, which significantly contributed to improving the quality of the manuscript.