**3. Results**

In this section, we present and comment on the empirical results relating both to the differences between pre-and-post policy averages and medians and to the policy intervention analysis for the Milan Area B case study. Section 3.1 shows the variations in concentration levels of NOx and NO2 for all the stations installed in Milan and for the other seven cities around it. Section 3.2 presents the model selection results, the values of the selection criteria, and the final model specifications. Section 3.3 reports the empirical estimates of the policy effects obtained through the basic structural model augmented by the policy intervention.

#### *3.1. Average and Median Differences*

Empirical differences of concentrations levels for all the considered stations are presented in Table 1. For both nitrogen oxides and nitrogen dioxide, using the difference of mean and median respectively, it reports the estimates of the difference between the average (the median) concentrations for the year 2019 and the average (the median) concentrations of the same period for the years 2014–2018.


**Table 1.** Differences between the average concentration level of the sub-period 2014–2018 and the treatment period 2019. Differences are expressed in μg/m3.

The estimates highlight large negative differences in oxides concentrations between 2019 and the period 2014–2018, both in the metropolitan area of Milan and in almost all the surrounding towns. Particularly heavy reductions, and similar to those in Milan, were recorded in the cities of Bergamo (East) and Pavia (South). The simultaneous abatement inside and outside Milan confirms the presence of a general decreasing trend in the aggregate levels of pollutants for the Lombardy Po basin as already indicated by the previous figures.

The differences registered in Milan are relevant both in suburban districts, such as the stations Liguria (West) and Marche (North), and in the historical centre at the Senato station. For those monitoring stations, the reductions are larger than 16 μg/m<sup>3</sup> for NOx and 10 μg/m<sup>3</sup> for NO2. In general, the differences between the averages and between the medians are quite similar, but in many stations, the reductions for the medians are stronger than the average differences. This fact is related to the skewed and non-symmetric characteristics of the distribution involved, as shown also in Figure 4. The above considerations on average and median pollution abatement are valid for both pollutants, in fact, the stations where the greatest differences are recorded for nitrogen oxides are the same for nitrogen dioxide.

These preliminary results do not allow for identifying the causes of the reductions and to state if they depend on common causes related to the environment and climatic factors or if they have been generated by the introduction of the new policy in Milan. The next section will attempt to investigate the variations through the modeling of possible environmental and anthropogenic factors able to influence the air quality of the city.

#### *3.2. Model Selection*

#### 3.2.1. Step 1: Detection of the Seasonal Components

Results relative to the first step of model selection are summarized in Figures 5 and 6, which show the evaluation criteria for all the stations. For each station, the plots are organized in paired-panels; the left panel represents the 10-steps-ahead MSFE as a function of the forecast horizon and the number of harmonics modeling the seasonality (scale colour); the right panel shows the AICc–BIC pairs for each model. The optimal number of harmonics to model the seasonality is identified as the one that evaluates the minimum pair of AICc and BIC and returns the lowest MSFE curve. Both the estimates for NOx and NO2 for the city of Milan agree unanimously in the selection of the model in which the seasonality is composed by a single harmonic (*k* = 1); therefore, it can be rewritten as

$$\text{Optimal\\_seasonal}: \quad \gamma\_t = \gamma\_{1,t} \,. \tag{10}$$

where *γ*1,*t* is

$$
\begin{bmatrix}
\gamma\_{1,t} \\
\gamma\_{1,t}^\*
\end{bmatrix} = \begin{bmatrix}
\cos(2\pi/365) & \sin(2\pi/365) \\
\end{bmatrix} \begin{bmatrix}
\gamma\_{1,t} \\
\gamma\_{1,t}^\*
\end{bmatrix} + \begin{bmatrix}
\omega\_{1,t} \\
\omega\_{1,t}^\*
\end{bmatrix} \tag{11}
$$

*ωt* ∼ *WN*(0, *σ*<sup>2</sup> *ω*) and *ω*<sup>∗</sup> *t* ∼ *WN*(0, *σ*<sup>2</sup> *ω*).

**Figure 5.** Model selection-Step 1-NOx. Seasonal component selection for the five nitrogen oxides stations in Milan. Left panel: 10-steps-ahead MSFE in log-scale as function of the number of harmonics. Right panel: AICc and BIC pairs for each model.

**Figure 6.** Model selection-Step 1-NO2. Seasonal component selection for the 5 nitrogen dioxide stations in Milan. Left panel: 10-steps-ahead MSFE in log-scale as function of the number of harmonics. Right panel: AICc and BIC pairs for each model.

#### 3.2.2. Step 2: Detection of the Counter-Factual Component

After selecting the seasonal component, we proceed to the selection of the counter-factual term. Estimates are summarized in Figures 7 and 8, which show the results for NOx and NO2. The plots are graphically organized like those related to step one, with the difference that the MSFEs and the AICc-BIC pairs are functions of one of the seven counter-factual candidates instead of the number of harmonics. The selection criteria follow the same rules used for the previous step.

The search for the optimal counter-factual term requires greater attention and detail than in the previous step as the minimizers are not unique. According to the plots, there is a restricted set of stations that are good candidates for the counter-factual role. The set includes the following cities: Treviglio, Pavia, Saronno, and Cremona. In particular, Pavia's station achieves one of the best forecast and fitting performances for almost all the stations in Area B for both NO2 and NOx. Based on this last consideration, we select as the counter-factual term for future models the air quality monitoring station of Pavia, located South to Milan. Therefore, the final specification of the basic structural model will include as counter-factual term the logarithm of the concentrations in Pavia, *xt* = *log*(*Paviat*).

**Figure 7.** Model selection-Step 2-NOx. Counter-factual term selection for the five nitrogen oxides stations in Milan. Left panel: 10-steps-ahead MSFE in log-scale as function of the candidate. Right panel: AICc and BIC pairs for each model.

**Figure 8.** Model selection-Step 2-NO2. Counter-factual term selection for the five nitrogen dioxide stations in Milan. Left panel: 10-steps-ahead MSFE in log-scale as function of the candidate. Right panel: AICc and BIC pairs for each model.

3.2.3. Step 3: Detection of the Weather and Calendar Factors

The last step of model selection aims to select the optimal subset of local weather and calendar regressors after having selected the optimal unobservable components, common weather, and socio-economic factors, captured by the counter-factual. For each station and pollutant, the best models are reported in Tables 2 and 3.

**Table 2.** Model selection-Step 3-NOx: Best subset of covariates using backward-forward stepwise algorithms for NOx.


*Note*: symbol indicates that the regressor is selected within the best subset of covariates.

**Table 3.** Model selection-Step 3-NO2: Best subset of covariates using backward-forward stepwise algorithms for NO2.


*Note*: symbol indicates that the regressor is selected within the best subset of covariates.

As expected, BIC-based models, being more parsimonious, retain fewer covariates than AIC-based models. Following this fact, we will use the BIC-selected models, but we now discuss some details about the selection process. Concerning the calendar, both criteria include in almost all cases the holidays, week-end, and Sunday holidays effects. AIC suggests adding also the interaction term between Sunday and holidays. Even if the interaction between Saturday and holidays is not always included, the final model will take into account the full set of calendar events and their interactions.

Regarding weather covariates, except for the wind speed, none of the others is included within the final models. Winds blowing from Southwest (QSW) are always selected and those coming from Northwest (QNW) are often included, and hence kept in the final model. Moreover, temperature, rainfall, solar radiation, and humidity are considered only by the AIC. This fact can be explained by the presence of the counter-factual, which captures not only common characteristics in terms of human behaviour and air quality conditions but also homogeneous climatic conditions common to all the areas considered.

#### 3.2.4. Final Model Specification

Based on the results of the three-step model selection procedure, the final specification of the BSM augmented by the policy intervention can be expressed using the following model:

 ,

$$\begin{aligned} Macassurement: \quad y\_t &= \mu\_t + \gamma\_t + \theta \mathbf{x}\_t + \phi\_1 Holidays + \phi\_2 \mathbf{W}eckEnd \\ &+ \phi \mathbf{y}Saturday: Holidays + \phi\_4 Sunday: Holidays \\ &+ \phi \mathbf{y}WindSpeed\_{\text{SW}} + \phi\_6 \mathbf{W}indSpeed\_{\text{NW}} \\ &+ \delta\_1 D\_{1t} + w\_t + \varepsilon\_{t,\*} \end{aligned} \tag{12}$$

where *εt* ∼ *N*(0, *σ*<sup>2</sup> )

$$\text{Sasonal component}: \quad \gamma\_t = \gamma\_{1,t}, \tag{13}$$

$$\text{Level}: \quad \mu\_t = \mu\_{t-1} + \beta\_{t-1} + \eta\_t \qquad \eta\_t \sim \text{NN}(0, \sigma\_\eta^2) \; , \tag{14}$$

$$\text{Slope}: \quad \beta\_t = \beta\_{t-1} + \zeta\_t \qquad \zeta\_t \sim \text{WN}(0, \sigma\_{\zeta}^2) \text{ }, \tag{15}$$

$$Transtory\ policy:\quad w\_t = \lambda w\_{t-1} + \delta\_0 D\_{2,t} \, . \tag{16}$$

where *yt* is the logarithm of pollution concentrations in one of the stations located in Milan and *xt* = *log*(*Paviat*) is the logarithm of pollution concentrations in Pavia.

#### *3.3. Basic Structural Model and Policy Intervention*

In this section, we show the numerical results obtained using state space modeling to estimate both permanent and transitory effects generated by the introduction of Area B controlling for local weather conditions, anthropogenic effects, and common areal trends.

The maximum likelihood estimates of the coefficients and the components' variances for the five air quality monitoring stations installed in Milan are reported in Tables 4 and 5. The results appear to be coherent both for nitrogen oxides and nitrogen dioxide. First, the models identify statistically significant and positive coefficients for the counter-factual term, highlighting its capability to capture socio-economic and climatic factors common to neighbouring areas and coherent with the expected sign. Second, weekends and holidays exert a negative effect on concentration levels probably due to the reduction in the movements and productive activities of the city in those days. Their interactions are almost everywhere not statistically significant but with a positive sign and always less than the sum of the individual effects of the weekend and holidays. This fact underlines how the holiday weekends enjoy more contained effects of emission reductions compared to generic weekends of the year. Third, as to be expected, winds blowing from the West (QSW and QNW) greatly reduce the amount of pollutants all over the city with peaks over 40% to 50%.


**Table 4.** ML estimates of BSM parameters and variances for NOx.

*Note 1*: values in parenthesis are standard errors. *Note 2*: \* *p* < 0.10, \*\* *p* < 0.05, \*\*\* *p* < 0.01.


**Table 5.** ML estimates of BSM parameters and variances for NO2.

*Note 1*: values in parenthesis are standard errors. *Note 2*: \*\* *p* < 0.05, \*\*\* *p* < 0.01.

The short-term impacts adjusted for common anthropic and weather factors are summarized in Table 6, which shows the estimated permanent and transitory effects for each station in Milan expressed in logarithmic scale, hence interpretable as relative variations in concentrations levels. None of these two coefficients identifies an improvement of the considered pollutant concentrations. Moreover, the permanent effect (*δ*1) is always positive and in some cases moderately statistically significant. This means that, compared to the generally decreasing areal trend, Milan air quality went worst. It is worth observing that the most significant results are obtained at the Senato station, which is located in the already existing Area C, hence already subject to some car traffic restrictions.

Such a result could be justified by the presence of multiple causes. As a first justification, we are approaching the initial phase of a progressive policy and the time elapsed since its outset may be too short to assess any significant impacts on pollutant levels. This fact is consistent with the forecasts expected by the municipality of Milan about the reductions in nitrogen oxide levels; in fact, the first significant reductions should be observed starting from 2022 [11]. Furthermore, since we are dealing with limitations to human behavior and social perception of new norms, it is not always clear how agents adapt to changes. The deterioration in the air quality of the centre could be linked to new traffic congestions in that area or to a panic shock of drivers, who need time to understand the functioning of the restrictions and adapt their behavior, exactly as in situations mismanagement of individual and organizational changes [41,42]. Eventually, the recent climate changes and the extreme weather conditions that affected the Po valley, such as temperatures higher than the seasonal average and extreme atmospheric events, could increase the noise present in the data and thus mask the real repercussions of the limitations.


**Table 6.** Estimated permanent and transitory effects in log scale on NOx and NO2 for each station.

> *Note 1*: \* *p* < 0.10, \*\* *p* < 0.05.
