\documentclass[11pt]{article}
\def\LEZCleanCopy{1}
\usepackage[T1]{fontenc}
\usepackage{amsmath,amssymb}
\usepackage{graphicx}
\usepackage{booktabs}
\usepackage{longtable}
\usepackage{array}
\usepackage{geometry}
\usepackage{setspace}
\usepackage{placeins}
\usepackage{xcolor}
\usepackage{lineno}
\usepackage[hidelinks]{hyperref}
\ifdefined\LEZCleanCopy
  \newcommand{\rev}[1]{#1}
  \newenvironment{revision}{}{}
\else
  \newcommand{\rev}[1]{\textcolor{blue}{#1}}
  \newenvironment{revision}{\color{blue}}{}
\fi
\geometry{margin=1in}
\setlength{\emergencystretch}{3em}
\renewcommand{\thesection}{S\arabic{section}}
\renewcommand{\thetable}{S\arabic{table}}
\renewcommand{\thefigure}{S\arabic{figure}}
\begin{document}
\begin{center}
{\Large Supplementary Information: Online Resource 1}\\[0.6em]
{\large Monitoring-based assessment of spatial and seasonal air-quality associations with low-emission zones in South Korea}\\[0.8em]
Environmental Monitoring and Assessment\\[0.4em]
Minju Lee and Mijeong Kim\\
Department of Statistics, Ewha Womans University\\
Corresponding author: Mijeong Kim (m.kim@ewha.ac.kr)
\end{center}
\onehalfspacing
\linenumbers

\section{Monthly panel and event-study diagnostics}
The common sample contains 247 districts after excluding Ongjin-gun and Ulleung-gun. The 2012--2024 and 2013--2024 panels contain 156 and 144 months, or 38,532 and 35,568 district-months, respectively.

Both periods use the 12 controls listed in main-text Table 2. All four period--estimator combinations use not-yet-treated comparisons, a universal month $-1$ reference, and 1,999 multiplier-bootstrap iterations.

For district $i$ and pre-adoption contrast $j$, let $\widehat{\psi}_{ij}$ denote the fitted influence function and let
\[
\widehat{s}_j=\frac{\sqrt{\sum_i\widehat{\psi}_{ij}^2}}{n}.
\]
The maximum-$|t|$ statistic is $\max_j|\widehat{ATT}_j/\widehat{s}_j|$. Its multiplier reference distribution uses $\max_j|n^{-1}\sum_i \xi_i\widehat{\psi}_{ij}/\widehat{s}_j|$, where each $\xi_i$ is independently $-1$ or $1$ with equal probability and the same draw is used across all contrasts of district $i$. This preserves their joint dependence without inverting the pre-period covariance. Independent district multipliers do not accommodate cross-district dependence or policy interference. These diagnostic tests and bands are distinct from the score-HAC inference for the SDM.

\begin{table}[htbp]
\centering
\caption{Maximum-$|t|$ pre-adoption diagnostic p-values}
\label{tab:s_pretests}
\small
\setlength{\tabcolsep}{5pt}   
\input{tables/Table S1} 
\begin{flushleft}\footnotesize
Notes: \texttt{reg}: regression adjustment; \texttt{dr}: doubly robust estimation. The 2012 and 2013 starts have 231 and 195 non-reference pre-period restrictions, respectively, with full numerical covariance rank. All joint Wald p-values are below 0.001. Multiplier p-values use the $+1$ correction, giving a minimum of 0.0005. The 2020 cohort contains only 13 districts.
\end{flushleft}
\end{table}

\begin{figure}[p]
\centering
\includegraphics[width=\linewidth,height=0.80\textheight,keepaspectratio]{figure/Fig S1.jpg}
\caption{Doubly robust counterpart of the main event-study diagnostics, 2012--2024. Points aggregate available cohorts at each event month; the 2020 cohort does not contribute at month 60. Shading indicates simultaneous 95\% bands, with month $-1$ as the universal reference. Intervals containing zero do not establish equivalence.}
\label{fig:s_dr}
\end{figure}
\FloatBarrier

\section{Seasonal spatial quantities and trend sensitivity}
\begin{table}[!htbp]
\centering
\caption{Seasonal SDM components and trend-adjusted totals underlying main-text Figure 3}
\label{tab:s_seasons}
\singlespacing
\footnotesize
\setlength{\tabcolsep}{4pt}   
\renewcommand{\arraystretch}{1.02}  
\input{tables/Table S2} 
\begin{flushleft}\footnotesize
Notes: Each fit uses 247 districts and 39 months, separate coefficients and fixed effects, and a common five-neighbor matrix. Parentheses: score-HAC(12) delta-method standard errors, retaining calendar gaps. Gases: ppb; PM$_{10}$: $\mu$g/m$^3$. Winter retains individual December--February months and their designation indicators. The last column adds three cohort trends and their spatial lags. 
\end{flushleft}
\end{table}

The supplementary May--June O$_3$ model uses 6,422 district-months and yields a total of 2.680 ppb (HAC SE 0.741). The main text gives the atmospheric rationale for this overlapping window.
\FloatBarrier

\section{Geographic, control, period, and timing sensitivity}
\rev{Table~\ref{tab:s_sensitivity} changes one main specification feature at a time. Three- and seven-neighbor matrices preserve the full 247-unit sample. Omitting transport controls removes both daily vehicle distance and cars per capita and their spatial lags. The 2013-start comparison retains the original covariates, isolating sensitivity to the observation window. The land-use omission removes industrial area, commercial area, and green area per capita and their spatial lags. The January 2018 timing check moves only the second cohort's designation six months earlier.}

\begin{table}[htbp]
\centering
\caption{Additional sensitivity of total monthly spatial quantities}
\label{tab:s_sensitivity}
\singlespacing
\small
\setlength{\tabcolsep}{4pt}  
\input{tables/Table S3}
\begin{flushleft}\footnotesize
Notes: Score-HAC standard errors with a 12-calendar-month bandwidth are in parentheses. All models include district and calendar year-month fixed effects. \rev{The 2013-start specification has 35,568 observations and other specifications have 38,532.} Values use the outcome units of the main tables. Cohort-exclusion estimates are in the main text.
\end{flushleft}
\end{table}
\FloatBarrier

\section{Full monthly SDM coefficient estimates}
Table~\ref{tab:s_coefficients} reports full primary-model coefficients with HAC(12) standard errors. All fits use 38,532 observations and absorbed district and year-month fixed effects. Continuous controls are standardized before constructing spatial lags; outcomes and designation remain unscaled. The coefficient on $WY$ is $\rho$; the spatial multiplier converts coefficients to the associations reported in the main text.

\clearpage
\begingroup
\singlespacing
\footnotesize
\setlength{\tabcolsep}{5pt}
\renewcommand{\arraystretch}{0.95}
\setlength{\LTpost}{0pt}   
\input{tables/Table S4}
\endgroup
\FloatBarrier

\section{Residual dependence and uncertainty estimation}
Table~\ref{tab:s_residuals} summarizes residual dependence under the within-transformed SDM criterion. Temporal correlations use matched calendar months within districts; monthly Moran's $I$ uses centered residuals and the fitted five-neighbor matrix.

\begin{table}[htbp]
\centering
\caption{Dependence diagnostics for the primary monthly SDM}
\label{tab:s_residuals}
\small
\setlength{\tabcolsep}{4pt}
\input{tables/Table S5}
\begin{flushleft}\footnotesize
Notes: Lag-1 and lag-12 entries summarize 247 district-level correlations at one- and twelve-calendar-month separations. IQR: interquartile range. Moran's $I$ is summarized by its median across 156 months.
\end{flushleft}
\end{table}

Let $\widetilde{y}_t$ and $\widetilde{X}_t$ be the two-way within-transformed outcome and design, with the design containing own and spatially lagged regressors. To match the implemented \texttt{splm} criterion, define $z_t=W\widetilde{y}_t$ and $e_t=\widetilde{y}_t-\rho z_t-\widetilde{X}_t b$. \rev{For $\omega=\log\sigma^2$ (distinct from the spatial-control vector $\eta$ in the main text), its contribution, up to a constant, is}
\[
\rev{\ell_t=\log|I-\rho W|-\frac{n}{2}\omega-\frac{e_t'e_t}{2\exp(\omega)}.}
\]
Writing $Q=(I-\rho W)^{-1}W$, the joint score is
\[
s_t=\begin{pmatrix}
-\operatorname{tr}(Q)+z_t'e_t/\sigma^2\\
\widetilde{X}_t'e_t/\sigma^2\\
-n/2+e_t'e_t/(2\sigma^2)
\end{pmatrix}.
\]
Let $H=-\sum_t\partial^2\ell_t/\partial\vartheta\partial\vartheta'$ be the joint observed information, where \rev{$\vartheta=(\rho,b',\omega)'$}; the bread is $H^{-1}$. We retain its off-diagonal blocks rather than treating $\rho$ as independent of the regression coefficients.

\begin{samepage}
Scores are centered over observed months. With $\bar s_t$ denoting the centered score, the HAC meat is
\[
\begin{aligned}
M_L&=\sum_t\bar s_t\bar s_t'
+\sum_{h=1}^{L}\left(1-\frac{h}{L+1}\right)
\sum_t\left(\bar s_t\bar s_{t-h}'+\bar s_{t-h}\bar s_t'\right),
\\
\widehat V_L&=H^{-1}M_LH^{-1}.
\end{aligned}
\]
\end{samepage}
The primary bandwidth is 12 calendar months; sensitivity choices are 6 and 24. In seasonal or May--June fits, unobserved months contribute zero scores, preserving actual calendar gaps. Cross-sectional aggregation retains contemporaneous score dependence among districts. The construction assumes weak temporal dependence and sufficient time information, with seasonal inference approximate for 39 observed months. Covariance adjustment leaves coefficients unchanged and does not address mean misspecification or interpolation uncertainty.

For any direct, indirect, or total quantity $q(\rho,b)$, its variance is $\nabla q'\widehat V_L\nabla q$, with a zero derivative for \rev{$\omega$}. Approximate intervals use $q\pm1.96\operatorname{SE}(q)$. Table~\ref{tab:s_uncertainty} compares the total SE from $H^{-1}$, denoted OIM, against the three score-HAC choices. All SDM tables and the seasonal figure use HAC(12), except the explicitly labeled covariance comparison. Non-spatial benchmarks use district-clustered errors, while group-time diagnostics use their separate multiplier procedure.

\begin{table}[htbp]
\centering
\caption{Total-association uncertainty under alternative covariance estimates}
\label{tab:s_uncertainty}
\footnotesize
\setlength{\tabcolsep}{3pt} 
\input{tables/Table S6}
\begin{flushleft}\footnotesize
Notes: All columns use the same point estimates. OIM: joint observed information; HAC(6), HAC(12), and HAC(24): calendar-month Bartlett bandwidths. Gases: ppb; PM$_{10}$: $\mu$g/m$^3$. 
\end{flushleft}
\end{table}
\FloatBarrier

\begin{revision}
\section{Cross-season comparisons and multiplicity}
\label{sec:s_season_tests}
Let $q_s=(\beta_s+\theta_s)/(1-\rho_s)$ be the total spatial association in season $s$, $H_s$ its full joint observed information, and $\bar s_{st}$ its centered score at calendar month $t$. For $g_s=\nabla_{\vartheta_s}q_s$, the only nonzero gradient entries are
\[
\frac{\partial q_s}{\partial\rho_s}=\frac{\beta_s+\theta_s}{(1-\rho_s)^2},\qquad
\frac{\partial q_s}{\partial\beta_s}=
\frac{\partial q_s}{\partial\theta_s}=\frac{1}{1-\rho_s}.
\]
Define the scalar influence contribution $a_{st}=g_s'H_s^{-1}\bar s_{st}$ for observed months of season $s$ and zero otherwise. Centering uses the 39 observed months within that seasonal fit, before padding onto the complete 156-month calendar. For seasons $s$ and $r$, the joint covariance estimate is
\[
\widehat V_{sr}=\sum_{t=1}^{156}\sum_{u=1}^{156}
 k_{12}(t-u)a_{st}a_{ru},\qquad
 k_{12}(h)=\max\left(1-\frac{|h|}{13},0\right).
\]
The covariance estimator preserves calendar separation and cross-season score dependence; its diagonal elements are the season-specific HAC(12) variances.

Let $\widehat q=(\widehat q_{\mathrm{winter}},\widehat q_{\mathrm{spring}},\widehat q_{\mathrm{summer}},\widehat q_{\mathrm{autumn}})'$ and $C=[-\mathbf{1}_3\ I_3]$. The omnibus statistic is
\[
T=(C\widehat q)'(C\widehat V C')^{-1}(C\widehat q),
\]
compared with an approximate $\chi^2_3$ reference distribution. Pairwise seasonal differences use their corresponding contrast vector and the same joint covariance. The calculation is repeated for baseline and cohort-trend models. Its target is equality of model-implied seasonal associations conditional on each specification.

Holm adjustment uses three separate families: ten omnibus tests; 60 pairwise comparisons (six pairs, five pollutants, two specifications); and five pooled baseline total-association tests. These are exploratory comparisons. Displayed intervals are marginal normal approximations, not simultaneous or multiplicity-adjusted intervals; other sensitivity analyses remain outside these families.

\begin{table}[htbp]
\centering
\caption{\rev{Equality tests for the four seasonal total spatial associations}}
\label{tab:s_season_tests}
\ifdefined\LEZCleanCopy\else\color{blue}\fi
\small\setlength{\tabcolsep}{5pt}
\input{tables/Table S7}
\begin{flushleft}\footnotesize
Notes: Each test compares winter, spring, summer, and autumn for one pollutant and specification. Covariance retains cross-season dependence using the calendar-month HAC(12) kernel. Holm p-values adjust across all ten rows. df: degrees of freedom.
\end{flushleft}
\end{table}
\FloatBarrier

For the baseline ozone contrasts discussed in the main text, spring minus winter is 2.016 ppb (SE 0.427; marginal 95\% interval $[1.180,2.852]$; Holm pairwise $p=0.000138$), and summer minus winter is 2.532 ppb (SE 0.557; $[1.440,3.624]$; adjusted $p=0.000325$). In the separate pooled-baseline family, the adjusted NO$_2$ p-value is 0.00162 and the remaining adjusted p-values are at least 0.499.


\end{revision}

\section{Interpolation validation and regression sensitivity}
\subsection{Station prediction performance}
For each eligible station--month--variable observation, leave-one-out validation removes that station and predicts its value using the other eligible stations in the same month. IDW uses power 2 in projected meters; nearest prediction uses the closest remaining eligible station. Table~\ref{tab:s_prediction} reports pooled absolute and squared-error summaries for the five urban ambient outcomes and six meteorological controls. Holdouts are station-month predictions, not independent sampling units.

\begin{table}[htbp]
\centering
\caption{Leave-one-station-out prediction errors for analysis variables}
\label{tab:s_prediction}
\singlespacing\footnotesize
\setlength{\tabcolsep}{3pt}
\input{tables/Table S8}
\begin{flushleft}\footnotesize
Notes: MAE: mean absolute error; RMSE: root mean squared error. Methods use identical holdouts within each variable; counts differ with variable-specific eligibility. Errors retain the variable's units and assess point prediction, not district-area accuracy.
\end{flushleft}
\end{table}

Although IDW has lower station-holdout RMSE for all 16 environmental series, its urban NO$_2$ mean signed error is 1.319 ppb, versus 0.513 for nearest prediction. The low-wind equivalent has an IDW RMSE of 7.595 days.

For precipitation only, the additional spatial test excludes all stations in the held-out station's 100-km square grid block. With 14,555 held-out station-months, IDW and nearest RMSEs are 1.700 and 1.769~mm/day, compared with 1.422 and 1.636 under single-station exclusion. Their blocked MAEs are 0.989 and 0.983~mm/day, respectively. This is not a 100-km exclusion radius: nearby stations across block boundaries can remain. The two designs therefore assess station prediction under different local data availability.
\FloatBarrier

\subsection{Common-support model comparison}
Table~\ref{tab:s_interpolation} compares pooled totals using the IDW and nearest-station fields described in the main-text Materials and methods. Eligible station-months, grid, and aggregation weights are shared; annual covariates, designation, estimation sample, and policy-neighbor matrix remain fixed. Outcomes and weather change jointly, so their contributions to measurement sensitivity are not separated.

\begin{table}[htbp]
\centering
\caption{Pooled total associations under common-support interpolation alternatives}
\label{tab:s_interpolation}
\singlespacing\small
\setlength{\tabcolsep}{3pt}  
\input{tables/Table S9}
\begin{flushleft}\footnotesize
Notes: IDW and nearest-station grid predictions use identical area-aggregation weights. All fits use 247 districts and 156 months. Parentheses give score-HAC(12) delta-method standard errors. Gas units are ppb; PM$_{10}$ units are $\mu$g/m$^3$. Intervals are marginal and exploratory; no additional multiplicity-adjusted claim is made for this sensitivity table.
\end{flushleft}
\end{table}

The nearest-panel baseline SO$_2$ and PM$_{10}$ marginal 95\% intervals are $[-0.714,-0.260]$ ppb and $[-3.776,-0.977]$~$\mu$g/m$^3$, respectively, supplementing the point estimates and standard errors in Table~\ref{tab:s_interpolation}. The main text interprets the pollutant-specific interpolation and trend comparisons.
\FloatBarrier
\section{Policy-date evidence and timing sensitivity}
Table~\ref{tab:s_policy_dates} links the assigned permanent-LEZ rollout months to documentary evidence or explicit assumptions. As defined in the main text, this geographic measure is distinct from enforcement intensity and other vehicle-restriction schemes.

\begin{table}[htbp]
\centering
\caption{Source-to-code mapping for the permanent-LEZ rollout}
\label{tab:s_policy_dates}
\singlespacing\footnotesize
\setlength{\tabcolsep}{4pt}
\begin{tabular}{@{}>{\raggedright\arraybackslash}p{0.19\linewidth}>{\raggedright\arraybackslash}p{0.38\linewidth}>{\raggedright\arraybackslash}p{0.37\linewidth}@{}}
\toprule
Geographic cohort & Documentary support & Monthly coding and qualification\\
\midrule
Seoul (25 units) & Government rollout announcement (4 August 2016): 2017. Weekly Gonggam, No. 387 (9 January 2017), printed p. 36: January 2017 Seoul-wide expansion. & January 2017. Expansion of narrower pre-existing restrictions, not first-ever diesel regulation. District-specific first camera operation and compliance are unobserved.\\[3pt]
Incheon and 17 Gyeonggi cities (35 units) & Government announcement: 2018. MECAR permanent-restriction page: integrated operation from 1 July 2018. & July 2018. January 2018 is an alternative coding scenario, allowing for earlier administrative rollout. Neither assignment verifies every district's enforcement history.\\[3pt]
Further 11 Gyeonggi cities (13 units) & Government announcement: expansion in 2020. MECAR describes coverage extending to 28 Gyeonggi cities by 2020. & January 2020 is assumed. July and December 2020 are mid- and late-year sensitivity scenarios, not documented common enforcement dates.\\
\bottomrule
\end{tabular}
\begin{flushleft}\footnotesize
Sources: \href{https://www.korea.kr/news/policyNewsView.do?newsId=148819495}{Ministry of Environment rollout announcement}; \href{https://www.korea.kr/goNewsRes/attaches/2019.05/23/107_20170109_01IozeNoYDGJM000.pdf}{Weekly Gonggam, No. 387}; \href{https://www.mecar.or.kr/dr/info/oldDieselCarAlwaysDr.do}{MECAR permanent-restriction information}. 
\end{flushleft}
\end{table}
\FloatBarrier

Each alternative changes only one cohort's start month. The 247 districts, 156 months, environmental fields, controls, and five-neighbor matrix remain fixed. Cohort-specific trend indicators are defined by geographic cohort and are unchanged across timing scenarios. The second-phase scenario changes designation for 35 districts over six months (210 district-months). The July and December third-phase scenarios change 13 districts over six or eleven months (78 or 143 district-months). All scenarios are reported for all five pollutants, with and without cohort trends; no scenario is selected using its estimate or p-value.

Table~\ref{tab:s_timing} reports every timing scenario. All baseline intervals for pollutants other than NO$_2$ include zero. The marginal positive trend-adjusted CO interval persists under July 2020 but not January 2018 or December 2020; the main text summarizes the NO$_2$ comparison.

\clearpage
{\singlespacing\footnotesize
\setlength{\tabcolsep}{3pt}
\begin{longtable}{@{}>{\raggedright\arraybackslash}p{0.92in}lrrrr@{}}
\caption{Pooled total associations under alternative assigned rollout months}\label{tab:s_timing}\\
\toprule
Scenario & Outcome & Baseline (SE) & 95\% interval & Cohort trends (SE) & 95\% interval\\
\midrule\endfirsthead
\multicolumn{6}{l}{\textcolor{red}{Table S11} continued}\\
\toprule
Scenario & Outcome & Baseline (SE) & 95\% interval & Cohort trends (SE) & 95\% interval\\
\midrule\endhead
\bottomrule\endfoot
\input{tables/Table S11}
\end{longtable}}
\noindent\footnotesize Notes: Assigned months are January 2017, July 2018, and January 2020; each alternative changes only the named phase. Standard errors and marginal intervals use score-HAC(12) and the delta method, without multiplicity adjustment. Gases: ppb; PM$_{10}$: $\mu$g/m$^3$.\normalsize
\FloatBarrier
\end{document}