The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
79,037 characters
Anthropogenic Forcing, Climate Change, and the Shape of Warming: Statistical Inference for Distributional Cointegration
\def\spacingset#1{\renewcommand{\baselinestretch}
{#1}\small\normalsize} \spacingset{1}
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}
\newtheorem{corollary}{Corollary}
\newtheorem{assumption}{Assumption}
\newtheorem{lemma}{Lemma}
\newtheorem{remark}{Remark}
\newtheorem{example}{Example}
\theoremstyle{definition}
\newtheorem{exmp}{Example}[section]
\AtEndDocument{\refstepcounter{theorem}\label{finalthm}}
\AtEndDocument{\refstepcounter{proposition}\label{finalprop}}
\let\proglang=\textsf \let\code=\texttt
\setlength{\abovedisplayskip}{4pt}
\setlength{\belowdisplayskip}{4pt}
\date{}
\if11
{
\title{\bf Anthropogenic Forcing, Climate Change, and the Shape of Warming: Statistical Inference for Distributional Cointegration}
\author{\normalsize
Kyungsik Nam\\ \vspace{-0.8em} Division of Climate Change, Hankuk University of Foreign Studies \\ \vspace{1.2em}
Won-Ki Seo\thanks{The data and computational code used to reproduce the reported results are available in the public repository at \url{https://anonymous.4open.science/r/FRSTAT_Program-5F12/}.} \\
School of Economics, University of Sydney
}
\maketitle
}
\fi
\if01
{
\bigskip
\bigskip
\bigskip
\begin{center}
{\LARGE\bf Anthropogenic Forcing, Climate Change, and the Shape of Warming:\\ Statistical Inference for Distributional Cointegration}
\end{center}
\medskip
} \fi
\begin{abstract}
Anthropogenic forcing components follow different long-run paths, while persistent temperature change can involve distributional changes beyond the mean. Scalar regressions aggregate these components and retain only mean temperature, obscuring how distinct forcing paths relate to persistent distributional change. We develop new testing, estimation, and inference methods for long-run relations between an integrated predictor vector and a density-valued response. These comprise a residual-based test of between-cointegration (whether predictor trends account for all stochastic trends in the response density), a fully modified least-squares estimator of predictor-specific functional responses, and simulation-based inference for interpretable projections. We apply the methods to densities of observed local temperature anomalies and anthropogenic effective radiative forcing divided into CO$_2$ and non-CO$_2$ portfolios. The test results are consistent with persistent movements in these portfolios statistically accounting for the persistent evolution of the anomaly distribution, with no additional stochastic trend detected in the residual. A joint test rejects the common-response restriction imposed by aggregating the two portfolios. The fitted CO$_2$ response mainly shifts mass toward warmer anomalies and increases central concentration, whereas the non-CO$_2$ response produces a smaller shift but greater dispersion and off-center reshaping. Positive fitted mean responses for both portfolios conceal these contrasts, demonstrating the information lost through scalar aggregation.
\end{abstract}
\noindent
{\it Keywords: density-valued time series; fully modified least-squares estimator; cointegration; functional regression; radiative forcing; temperature-anomaly distribution.}
\vfill
\newpage
\spacingset{1.7}
\section{Introduction}\label{sec:intro}
\noindent Standard climate econometric models relate the global mean temperature anomaly to aggregate radiative forcing, reducing both sides of the long-run relation to scalars. This scalar-on-scalar formulation targets persistent movements in mean temperature but cannot recover distributional shape on the response side or forcing composition on the predictor side. Similar mean movements may coexist with changes in dispersion, asymmetry, and tail mass, while forcing components combined in an anthropogenic forcing index may have distinct distributional response profiles. Our analysis therefore focuses on stochastic trends, by which we mean persistent nonstationary components arising from the accumulation of stationary increments. We ask whether the stochastic trends in an anthropogenic forcing vector account for all persistent nonstationary variation in annual distributions of observed local temperature anomalies and how the component-specific response profiles vary across anomaly states.
\indent The scalar time-series literature provides a natural starting point for this analysis. Research on temperature persistence examines whether global and hemispheric anomaly series contain stochastic trends or are stationary around deterministic trends subject to structural breaks \citep{estrada2017extracting,kim2020inference}. Related work studies long-run relations between aggregate temperature and radiative forcing \citep{dergiades2016long}. Despite differing in their representations of persistence, these approaches retain a scalar temperature response. Specifications with multiple forcing variables relax aggregation on the predictor side, but their conclusions remain confined to persistent movements in the mean anomaly.
\indent Distributional evidence shows why the scalar response is restrictive. \citet{rivas2020trends} document changes across moments and quantiles of temperature distributions. Treating cross-sectional anomaly densities as functional time series, \citet{chang2020evaluating} find that stochastic trends operate through location and dispersion and differ across hemispheric distributions. Persistent temperature change may involve changes in shape as well as location. A scalar response cannot determine whether a long-run forcing relation operates through the center of the distribution or through dispersion, asymmetry, and tail mass.
\indent
Aggregation on the predictor side imposes a restriction. Effective radiative forcing accounts distinguish anthropogenic components with different signs and historical paths \citep{IPCC_AR6_WGI_AnnexIII_2021}. Unit-root evidence for related greenhouse-gas, sulfur-forcing, and total-forcing measures motivates an integrated long-run treatment of radiative forcing \citep{kaufmann2002cointegration,pretis2020econometric}. Components may differ in spatial incidence: using station-level data, \citet{magnus2011global} distinguish broadly distributed warming associated with greenhouse gases from more localized aerosol-related cooling. These differences make a common distributional response an empirical restriction rather than a consequence of forcing accounting. Once components are combined into an index, component responses cannot be identified.
\indent
The response-side and predictor-side extensions are therefore complements
rather than substitutes. Retaining the anomaly distribution preserves
information about shape, while retaining a forcing vector preserves information
about composition. Our empirical specification takes the annual distribution
of observed local temperature anomalies as a density-valued response and
separates anthropogenic forcing into carbon dioxide and the net contribution of
the remaining components. The conventional aggregate-mean regression serves as
the scalar benchmark, allowing the information recovered from distributional
shape and forcing composition to be assessed separately and jointly.
\indent We formulate this relation as distributional cointegration between an integrated predictor vector and a density-valued response, which we call \emph{between-cointegration}. The centered log-ratio (CLR) transformation represents the anomaly densities in a linear Hilbert space using the geometry of density-valued functional data \citep{Egozcue2006,petersen2016}. Between-cointegration holds when a vector-to-function long-run relation between the forcing variables and the transformed density process leaves a stationary functional residual. The restriction therefore asks whether persistent variation in the anomaly distribution is spanned by the forcing trends through a stable long-run operator. By contrast, we use the term \emph{within-cointegration} for cointegration arising from stationary linear combinations within a single process. This is the principal object of the existing Hilbert-space cointegration literature \citep{Chang2016152,BSS2017,NSS}, which does not directly address the cross-process testing, operator estimation, and response-function inference required here.
\indent The procedure follows the order in which the long-run relation is assessed and interpreted. We first construct residual-based tests for a functional response and an integrated vector predictor by extending the null-of-cointegration approach of \citet{shin1994residual}. Because their null distributions are nonpivotal, critical values are obtained by plug-in Monte Carlo. We then estimate the vector-to-function operator by fully modified least squares, extending \citet{phillips1995fully} to correct for long-run endogeneity and serial dependence. The limit theory yields simulation-based marginal inference for forcing-specific responses averaged over local anomaly intervals, following the projection approach of \citet{seong2021functional} and \citet{Nam2025}.
\indent
We apply the approach to anomaly densities from the HadCRUT5 monthly gridded temperature record and historical effective radiative forcing reconstructions split into CO$_2$ and non-CO$_2$ portfolios. A forcing-block diagnostic supports their treatment as distinct persistent coordinates. At the 5\% level, the between-cointegration results are consistent with a long-run relationship in which the distinct forcing trends account for the persistent evolution of the anomaly distribution. We then ask how strongly the distribution responds to each portfolio across anomaly states. The estimated response functions differ substantially in magnitude and shape, and inference on their difference rejects the common-response restriction imposed by aggregation, revealing heterogeneity concealed by scalar aggregation.
Substantively, an increase in the CO$_2$ forcing portfolio is associated mainly with mean and location shifts and greater central concentration, whereas an increase in the non-CO$_2$ portfolio is associated with greater dispersion. For CO$_2$, the cold- and warm-tail probability changes nearly offset; for non-CO$_2$, combined extreme-state probability increases primarily on the warm side. The mean response alone therefore misses these contrasting changes in dispersion and the allocation of probability across extreme anomaly states.
\indent The paper proceeds as follows. Section~\ref{sec:empirical_motivation} describes data, documents the distributional evolution of temperature anomalies, and characterizes the information lost through scalar aggregation. Section~\ref{sec:Metric_model} develops testing, estimation, and inference methods. Section~\ref{sec:distributional_responses} presents specification evidence, forcing-specific responses, fixed-basis common-response comparison, and implications for distributional shape and tail mass. Section~\ref{conclude} concludes. The Supplement contains proofs, technical assumptions, implementation details, and robustness checks.
\section{Data and the cost of scalar aggregation}
\label{sec:empirical_motivation}
\noindent This section specifies the two sides of the long-run relation and formalizes what is assumed away when either side is reduced to a scalar. The benchmark throughout is the scalar climate-econometric specification linking the global mean temperature anomaly to aggregate radiative forcing \citep[e.g.,][]{dergiades2016long,pretis2020econometric}. Section~\ref{subsec:response_data} constructs the density-valued response and forcing predictor vector, while Section~\ref{subsec:aggregation_cost} discusses the restrictions implicit in scalar reductions of the response and predictor.
\subsection{Response and predictor construction}
\label{subsec:response_data}
\noindent\textbf{The density-valued response.} We construct annual densities of temperature anomalies from observed grid-cell-month data using the ensemble mean of the non-infilled HadCRUT5 record for 1850--2024 \citep{morice2021updated}. The anomalies are measured in degrees Celsius relative to the 1961--1990 climatology. For each year $t$, we pool the available anomalies with equal weights and estimate $f_t$ by kernel smoothing. Each estimate is restricted to and renormalized on a common, tail-trimmed support. Section~\ref{sec_app_density_error} of the Supplement provides the data and numerical details.
\indent Following \citet{chang2020evaluating}, we use non-infilled observations because our estimand is the observed cell-month anomaly distribution; infilling would make dispersion and tail behavior reflect the spatial reconstruction model as well as the observations \citep{morice2021updated}. Observed spatial coverage expands over time, however, changing the locations represented in the annual densities. Section~\ref{sec_app_coverage} of the Supplement therefore examines sensitivity to this changing coverage.
\indent To use these densities in a functional linear model, we map them into a Hilbert space by the CLR transformation. For any positive density $f$ on the common support $[a,b]$, define
\begin{equation}\label{eq_clr_response} \operatorname{clr}(f)(s)=\log f(s)-\frac{1}{b-a}\int_a^b\log f(u)\,du,\qquad f(s)=\frac{\exp{\operatorname{clr}(f)(s)}}{\int_a^b\exp{\operatorname{clr}(f)(u)}\,du}. \end{equation}
For each year, the empirical functional response is $Y_t=\operatorname{clr}(f_t)\in\mathcal H_{Y}$. The CLR transformation is an isometry from the Bayes Hilbert space of densities onto the closed linear subspace $\{g\in L^2[a,b]:\int_a^b g(s)\,ds=0\}$ \citep{Egozcue2006}. The inverse in \eqref{eq_clr_response} shows that this representation retains all distributional information.
\begin{remark}[Density-estimation error]\label{rem_density_error}
Let $f_t^\circ$ denote the underlying density of temperature anomalies and $Y_t^\circ=\operatorname{clr}(f_t^\circ)$. We represent density-estimation error multiplicatively as $f_t(s)\propto f_t^\circ(s)\exp\{\eta_t(s)\}$, which is general for strictly positive densities. Equation~\eqref{eq_clr_response} then gives $Y_t=Y_t^\circ+e_t$, where $e_t(s)=\eta_t(s)-(b-a)^{-1}\int_a^b\eta_t(u)\,du$. Thus, the error enters additively on the CLR scale and is absorbed into the regression disturbance under the conditions in Section~\ref{sec_app_density_error} of the Supplement.
\end{remark}
\indent Panels~(a) and~(c)--(e) of Figure~\ref{Fig:Data_Desc} show that the evolution of the anomaly distribution is not a pure translation. The densities shift persistently to the right, but re-centering each year at its own cross-sectional mean does not collapse them onto a common shape: later densities are less concentrated than mid-sample densities and carry more mass on the warm shoulder. The difference panel makes the reallocation explicit, with a long-run loss of mass over moderately negative anomalies and a gain over near-zero, positive, and upper-tail states. Location, dispersion, and tail mass therefore all move over the sample, and a specification that tracks only the first moment cannot represent the last two.
\begin{figure}[t]
\centering
\captionsetup[subfigure]{font=small,skip=2pt}
\begin{subfigure}[b]{0.48\textwidth}
\centering
\includegraphics[width=\linewidth,trim={1.2cm 0.1cm 0.2cm 0.1cm},clip]{Figures/GTemp1.png}
\caption{Annual density surface}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.48\textwidth}
\centering
\includegraphics[width=\linewidth,trim={1.2cm 0.1cm 0.2cm 0.4cm},clip]{Figures/Forcing_Graph.png}
\caption{Grouped anthropogenic ERF}
\end{subfigure}
\vspace{0.4em}
\begin{subfigure}[b]{0.327\textwidth}
\centering
\includegraphics[width=\linewidth,trim={1.3cm 0.1cm 0.3cm 0.1cm},clip]{Figures/GTemp2.png}
\caption{Representative years}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.327\textwidth}
\centering
\includegraphics[width=\linewidth,trim={1.3cm 0.1cm 0.3cm 0.1cm},clip]{Figures/GTemp3.png}
\caption{Mean-aligned densities}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.327\textwidth}
\centering
\includegraphics[width=\linewidth,trim={1.3cm 0.1cm 0.3cm 0.1cm},clip]{Figures/GTemp4.png}
\caption{Late-minus-early difference}
\end{subfigure}
\caption{Temperature-anomaly densities and anthropogenic forcing, 1850--2024. F1 is CO$_2$ forcing; F2 is the sum of the non-CO$_2$ anthropogenic components listed in Table~\ref{Tab:ForcingAggregation}.}
\label{Fig:Data_Desc}
\end{figure}
\medskip
\noindent\textbf{The forcing predictor vector.} The forcing series are annual best-estimate reconstructions of Effective Radiative Forcing (ERF) from \citet{forster2024indicators}, as updated by \citet{Smith2025IndicatorsGCC} (v2025.06.25), over 1850--2024. The series follow the disaggregated accounting of IPCC AR6 \citep{IPCC_AR6_WGI_AnnexIII_2021} and are measured in W/m$^2$ relative to the 1750 baseline. Solar variability and volcanic aerosol forcing are excluded, so the predictor vector isolates persistent anthropogenic forcing rather than externally driven or episodic natural variation.
\begin{table}[t]
\centering
\caption{Anthropogenic forcing portfolios}
\label{Tab:ForcingAggregation}
\setlength{\tabcolsep}{5pt}
\small
\begin{tabularx}{\textwidth}{
>{\RaggedRight\arraybackslash}p{3.3cm}
>{\RaggedRight\arraybackslash}p{4.9cm}
>{\RaggedRight\arraybackslash\footnotesize}X}
\toprule
Group & Forcing components & Statistical role in the empirical design \\
\midrule
\makecell[tl]{F1: CO$_2$\\ forcing}
&
\makecell[tl]{CO$_2$}
&
The dominant anthropogenic greenhouse-gas component in the IPCC effective radiative forcing accounting. This portfolio provides the main low-frequency CO$_2$ forcing benchmark.
\\
\midrule
\makecell[tl]{F2: Non-CO$_2$\\ anthropogenic\\ forcing}
&
\makecell[tl]{CH$_4$, N$_2$O,\\ aerosol--radiation interactions,\\ aerosol--cloud interactions,\\ O$_3$, contrails, land-use\\ change, BC on snow,\\ H$_2$O$_{\text{strat}}$, halogenated species}
&
Anthropogenic forcing outside the CO$_2$ block. This portfolio collects non-CO$_2$ greenhouse gases, reactive-chemistry, aerosol, cloud-adjustment, surface-albedo, snow-albedo, and aviation-related forcing channels.
\\
\bottomrule
\end{tabularx}
\end{table}
With $T=175$ annual time points, estimating eleven component-specific functional responses would be imprecise because the forcing series share substantial low-frequency variation. We therefore group the components into the two portfolios in Table~\ref{Tab:ForcingAggregation}. F1 contains CO$_2$ forcing, whereas F2 sums the remaining ten anthropogenic components. Keeping CO$_2$ separate also provides a forcing-based counterpart to the CO$_2$-only greenhouse-gas specification of \citet{magnus2011global}, while F2 preserves a separate non-CO$_2$ forcing margin. This specification allows F1 and F2 to have distinct response functions; it does not identify component-specific responses within F2.
\indent Panel~(b) of Figure~\ref{Fig:Data_Desc} shows that F1 and F2 contain distinct low-frequency information. F1 remains positive and follows a smooth, nearly monotone upward path, whereas F2 is nonmonotone: negative aerosol forcing dominates much of the twentieth century, while non-CO$_2$ greenhouse gases and other positive components dominate in the recent period. Thus, despite a sample correlation in levels of $0.6109$, their historical paths differ substantially. The no-within-cointegration diagnostic reported in Section~\ref{sec_test_stat} formally examines whether F1 and F2 carry two stochastic trends. Aggregating them into an aggregate anthropogenic forcing index would collapse the two coordinates into their sum, netting their paths where their signs differ and imposing a single response function on the temperature density.
\subsection{What scalar aggregation misses: a functional alternative}
\label{subsec:aggregation_cost}
\noindent The matched scalar benchmark applies the aggregate anthropogenic forcing regression to the mean of the same observed temperature anomaly density:
\begin{equation} \bar y_t=\alpha+\beta^{\mathrm{agg}}F_t^{\mathrm{agg}}+u_t,\label{eq:linear_regression_sec6} \end{equation}
where $\alpha$ is an intercept, $u_t$ is a scalar disturbance, $\bar y_t=\int_a^b s f_t(s)\,ds$, $F_t^{\mathrm{agg}}$ is aggregate anthropogenic radiative forcing, and $\beta^{\mathrm{agg}}$ is the long-run scalar slope, measured in $^\circ$C per W/m$^2$, associated with a 1 W/m$^2$ increase in aggregate anthropogenic forcing.
\indent Our formulation generalizes both sides of \eqref{eq:linear_regression_sec6} by replacing the scalar forcing index with a $d_{\mathbf{x}}$-dimensional vector of forcing portfolios and the mean response with the full anomaly density in its CLR representation. Let $Y_t$ denote the CLR-transformed density and $\mathbf{x}_t=(x_{1,t},\ldots,x_{d_{\mathbf{x}},t})'$ the forcing vector, with $d_{\mathbf{x}}=2$ in our application. The model developed in Section~\ref{sec:Metric_model} takes the form
\begin{equation} Y_t(s)=\beta_0(s)+\sum_{k=1}^{d_{\mathbf{x}}}\beta_k(s)x_{k,t}+U_t(s),\qquad s\in[a,b].\label{eq:scalar_to_function_sec6}
\end{equation}
Here $\beta_0$ is a functional intercept, $\beta_k$ is the CLR response function of forcing portfolio $k$, and $U_t$ collects distributional variation not explained by the reconstructed forcing series.
\indent Aggregating the forcing vector replaces $\sum_k\beta_k(s)x_{k,t}$ with $\beta^{\mathrm{clr}}(s)F_t^{\mathrm{agg}}$, where $F_t^{\mathrm{agg}}=\sum_{k=1}^{d_{\mathbf{x}}}x_{k,t}$. The aggregate specification reproduces the vector model for every forcing path if and only if $\beta_k=\beta^{\mathrm{clr}}$ in $\mathcal H_{Y}$ for $k=1,\ldots,d_{\mathbf{x}}$, with equality understood almost everywhere. Hence, predictor aggregation is lossless only if CO$_2$ and non-CO$_2$ forcing share a common CLR response; otherwise, the distributional response depends on forcing composition. Reducing the response to its mean entails a separate loss, retaining only the first-moment response and leaving distributional reshaping unidentified. Section~\ref{sec:distributional_responses} examines this restriction empirically.
The central statistical question is whether the stochastic trends in the CLR-transformed anomaly distribution are spanned by those in the forcing vector, leaving a stationary functional residual. The next section formalizes this long-run restriction as between-cointegration and develops procedures for testing it, estimating the long-run response operator, and conducting inference on the forcing-specific responses.
\section{Statistical Framework}\label{sec:Metric_model}
\noindent We begin from the possibility that both the transformed anomaly distribution and the forcing vector contain stochastic trends, a feature examined empirically in Section~\ref{sec_test_stat}. In that case, the regression relation in \eqref{eq:scalar_to_function_sec6} has a long-run interpretation only if its residual is stationary. In the present application, this condition means that no persistent movement in the anomaly distribution remains after accounting for the anthropogenic forcing vector. We formalize this restriction as \emph{between-cointegration}. The statistical analysis follows the same logic: we first test whether between-cointegration holds, then estimate the long-run response operator under that restriction, and finally conduct inference on the forcing-specific response functions.
\indent Let $Y_t$ denote the CLR-transformed anomaly density, taking values in the Hilbert space $\mathcal H_{Y}\subseteq L^2[a,b]$, and let $\mathbf{x}_t=(x_{1,t},\ldots,x_{d_{\mathbf{x}},t})'$ denote the forcing vector, taking values in the $d_{\mathbf{x}}$-dimensional Euclidean space $\mathcal H_{\mathbf{x}}=\mathbb{R}^{d_{\mathbf{x}}}$. For the theoretical development, we write the empirical specification in \eqref{eq:scalar_to_function_sec6} in operator form as
\begin{equation}\label{eqmodel1}
Y_t=\beta_0+B(\mathbf{x}_t)+U_t,\quad B(\mathbf{v})=\sum_{j=1}^{d_{\mathbf{x}}}\beta_jv_j,\quad \mathbf{v}=(v_1,\ldots,v_{d_{\mathbf{x}}})' \in\mathcal H_{\mathbf{x}},
\end{equation}
where $B:\mathcal H_{\mathbf{x}}\rightarrow\mathcal H_{Y}$ is a linear response operator, $\beta_0\in\mathcal H_{Y}$ is a functional intercept, $U_t\in\mathcal H_{Y}$ is the unexplained component, and $\beta_j\in\mathcal H_{Y}$ is the CLR response function associated with forcing portfolio $j$. Between-cointegration holds when $U_t$ is stationary. Under this restriction, the stochastic trends in the anomaly distribution are fully spanned by those in the forcing vector; otherwise, a persistent distributional component remains unexplained and $B$ cannot be interpreted as a long-run response operator.
\indent The remainder of this section is organized as follows. Section~\ref{sec_model1} states the assumptions and formally defines between-cointegration; Section~\ref{sec_betcointeg} develops the residual-based test; and Section~\ref{sec_esti} presents estimation and inference under the maintained long-run relation. The paper then applies these procedures to the climate data in Section~\ref{sec:distributional_responses}.
\subsection{Model and assumption}\label{sec_model1}
\noindent Recall that $Y_t$ takes values in $\mathcal H_{Y}\subseteq L^2[a,b]$ and $\mathbf{x}_t$ in $\mathcal H_{\mathbf{x}}=\mathbb{R}^{d_{\mathbf{x}}}$, equipped with the $L^2$ and standard Euclidean inner products, respectively. We write $\langle\cdot,\cdot\rangle$ for either inner product. For $v_1\in\mathcal H_1$ and $v_2\in\mathcal H_2$, with $\mathcal H_1,\mathcal H_2\in\{\mathcal H_{\mathbf{x}},\mathcal H_{Y}\}$, define the rank-one operator $v_1\otimes v_2:\mathcal H_1\rightarrow\mathcal H_2$ by $(v_1\otimes v_2)(\cdot)=\langle v_1,\cdot\rangle v_2$. The identity operator on the relevant space is denoted by $I$. For any bounded linear operator $A:\mathcal H_1\rightarrow\mathcal H_2$, let $A^\ast:\mathcal H_2\rightarrow\mathcal H_1$ denote its adjoint, defined by $\langle Au,v\rangle=\langle u,A^\ast v\rangle$ for all $u\in\mathcal H_1$ and $v\in\mathcal H_2$. Further notation and mathematical preliminaries are collected in Section~\ref{Sec_prelim}.
\subsubsection{Stochastic trends and within-cointegration}\label{sec_basic}
Following work on nonstationary functional time series \citep[e.g.,][]{Chang2016152,BSS2017}, we assume a finite-dimensional stochastic-trend decomposition for the CLR-transformed anomaly-density process; Section~\ref{AP_FTS} of the Supplement gives formal conditions. Specifically, let $P^N$ and $P^S=I-P^N$ denote orthogonal projections on $\mathcal H_{Y}$, with $P^N$ of finite rank. We write
\begin{equation}\label{eqtimedecom}
Y_t=\mu_Y+Y_t^N+Y_t^S,\qquad Y_t^N=P^N(Y_t-\mu_Y),\qquad Y_t^S=P^S(Y_t-\mu_Y).
\end{equation}
Here, $\mu_Y$ is a deterministic functional level, representing the common or initial shape of the functional observations. The component $Y_t^N$ captures the persistent, nonstationary dynamics driven by stochastic trends arising from the accumulation of stationary random elements, while $Y_t^S$ denotes the stationary, mean-reverting component. Throughout this paper, we use the terms nonstationarity and persistence interchangeably. Although $Y_t$ resides in an infinite-dimensional space $\mathcal H_{Y}$, following the literature, we assume that its nonstationary component $Y_t^N$ is finite-dimensional (i.e., $P^N$ is a finite-rank projection). This assumption is not only empirically relevant (see e.g., \citealp{NSS}), but theoretically necessary to ensure feasible statistical inference with finite samples based on eigenanalysis to be discussed. Note that for any $v \in \operatorname{ran} P^S$, the inner product $\langle Y_t, v \rangle = \int Y_t(s)v(s)ds$, which can be viewed as a continuous linear combination, consists of a stationary sequence. In this context, $Y_t$ is said to be \textit{within-cointegrated}, reflecting an internal synchronization where individual nonstationary behaviors cancel each other out to reveal a stable long-run component.
\indent For the predictor $\mathbf{x}_t$, we assume that it is an $I(1)$ process of full nonstationary rank with a possibly nonzero deterministic level $\mu_X\in\mathcal H_{\mathbf{x}}$:
\begin{equation}\label{eqtimedecomx}
\mathbf{x}_t=\mu_X+\mathbf{x}_t^N.
\end{equation}
Here $\mathbf{x}_t^N:=\mathbf{x}_t-\mu_X$ denotes its nonstationary component, generated by the accumulation of stationary increments. Specifically, we assume that $\mathbf{x}_t$ admits no stationary linear combination $\langle\mathbf{x}_t,v\rangle$ for any nonzero $v\in\mathcal H_{\mathbf{x}}$. This implies that $\mathbf{x}_t$ is characterized by $d_{\mathbf{x}}$-dimensional persistent dynamics that are potentially aligned with the nonstationary behavior of $Y_t$. We focus on this purely nonstationary case because it matches our empirical application, where the no-within-cointegration diagnostic in Section~\ref{sec_test_stat} supports treating the forcing portfolios as distinct persistent coordinates.
\subsubsection{Statistical formulation of between-cointegration}\label{sec_between}
\noindent We next define the long-run relation between the density-valued response and the forcing vector. Within-cointegration concerns stationary linear combinations within a single multivariate or functional process. By contrast, between-cointegration concerns whether the stochastic trends in one process are accounted for by those in another. In the present setting, this means that the persistent component of $Y_t$ is spanned by the persistent components of $\mathbf{x}_t$.
Formally, we say that $Y_t$ and $\mathbf{x}_t=(x_{1,t},\ldots,x_{d_{\mathbf{x}},t})'$ are between-cointegrated if there exist nonrandom elements $\gamma_1,\ldots,\gamma_{d_{\mathbf{x}}}\in\mathcal H_{Y}$ such that \begin{equation}\label{eq:bet_cointeg} Y_t^N-\sum_{j=1}^{d_{\mathbf{x}}}\gamma_jx_{j,t} \end{equation} is stationary in $\mathcal H_{Y}$. Let $x_{j,t}^N$ denote the $j$th coordinate of $\mathbf{x}_t^N$. Since $\mu_X$ is deterministic, between-cointegration can equivalently be defined using $x_{j,t}^N$ in place of $x_{j,t}$ in \eqref{eq:bet_cointeg}. Under this condition, the long-run relation takes the form \eqref{eqmodel1}, with $\beta_j=\gamma_j$ for $j=1,\ldots,d_{\mathbf{x}}$ and stationary $U_t$. In the climate application, this means that after accounting for the anthropogenic forcing vector, no unexplained persistent component remains in the CLR-transformed temperature-anomaly distribution. Thus, testing between-cointegration provides a specification check for the long-run relation before the forcing-specific response functions are estimated.
\subsection{Statistical test for between-cointegration}\label{sec_betcointeg}
\noindent We develop a residual-based test of between-cointegration. The central idea is to determine whether the stochastic trends in $Y_t$ are fully accounted for by projecting $Y_t$ on the nonstationary forcing vector $\mathbf{x}_t$. Under between-cointegration, the resulting residuals should contain only stationary variation; otherwise, they should retain an unexplained stochastic trend. Section~\ref{sec_app_test} of the Supplement presents the asymptotic theory and examines finite-sample properties by simulation (Section~\ref{sec_sim_between}); here we focus on practical implementation in three steps.
\begin{description}[style=nextline, leftmargin=1em, font=\bfseries,topsep=-2pt, itemsep=0pt]
\item[Step 1: Computing residuals from the least-squares projection] \hfill
We first compute the least-squares projection of $Y_t$ on $\mathbf{x}_t$, using demeaned variables to allow for a nonzero intercept. Let $\bar Y_T=T^{-1}\sum_{t=1}^T Y_t$ and $\bar{\mathbf{x}}_T=T^{-1}\sum_{t=1}^T\mathbf{x}_t$. Define $\widehat{C}_{\mathbf{x}\mathbf{x}}=T^{-1}\sum_{t=1}^T(\mathbf{x}_t-\bar{\mathbf{x}}_T)\otimes(\mathbf{x}_t-\bar{\mathbf{x}}_T)$ and $\widehat{C}_{Y\mathbf{x}}=T^{-1}\sum_{t=1}^T(\mathbf{x}_t-\bar{\mathbf{x}}_T)\otimes(Y_t-\bar Y_T)$. The least-squares projection map is $\check B=\widehat C_{Y\mathbf{x}}\widehat C_{\mathbf{x}\mathbf{x}}^{-1}$. Let $(\widehat\lambda_{\mathbf{x},j},\widehat v_{\mathbf{x},j})$, $j=1,\ldots,d_{\mathbf{x}}$, denote the eigenvalue--eigenvector pairs of $\widehat C_{\mathbf{x}\mathbf{x}}$. The fitted component can then be computed as
\begin{equation}\label{eqprelimest}
\check B(\mathbf{x}_t-\bar{\mathbf{x}}_T)=\frac{1}{T}\sum_{s=1}^T\sum_{j=1}^{d_{\mathbf{x}}}\widehat\lambda_{\mathbf{x},j}^{-1}\langle \mathbf{x}_t-\bar{\mathbf{x}}_T,\widehat v_{\mathbf{x},j}\rangle\langle \mathbf{x}_s-\bar{\mathbf{x}}_T,\widehat v_{\mathbf{x},j}\rangle(Y_s-\bar Y_T).
\end{equation}
Intuitively, the fitted component represents the part of the transformed anomaly distribution that is accounted for by the forcing vector; the projection residuals $\check U_t=(Y_t-\bar Y_T)-\check B(\mathbf{x}_t-\bar{\mathbf{x}}_T)$ summarize the remaining variation. Under between-cointegration, the disturbance $U_t$ is stationary, so $\check U_t$ should contain no unexplained stochastic trend. The fitted residuals need not themselves be stationary in finite samples because they depend on an estimated projection and sample means; the plug-in calibration accounts for these effects.
\item[Step 2: Constructing residual-persistence diagnostics] \hfill
To measure the persistence remaining in $\check U_t$, define the residual partial-sum operator
\[ \widehat{\mathcal K}=\frac{1}{T}\sum_{t=1}^T\left(\sum_{s=1}^t\check U_s\right)\otimes\left(\sum_{s=1}^t\check U_s\right). \]
This operator accumulates residual variation over time. Under between-cointegration, $T^{-1}\widehat{\mathcal K}$ remains stochastically bounded despite the estimation and demeaning effects in $\check U_t$; under an alternative with a remaining stochastic trend, it diverges. We first consider the unnormalized diagnostic
\begin{equation}\label{eqteststat_K_main} \widehat{\mathcal T}_K=\Lambda_{\max}\left(T^{-1}\widehat{\mathcal K}\right), \end{equation}
where $\Lambda_{\max}(\cdot)$ denotes the largest eigenvalue. For a covariance-normalized diagnostic, we use the residual covariance operator $\widehat{\mathcal V}=T^{-1}\sum_{t=1}^T\check U_t\otimes\check U_t$ both to identify the direction of greatest contemporaneous residual variation and to normalize the residual scale. Let $\widehat v_V$ be the unit eigenvector associated with the largest eigenvalue of $\widehat{\mathcal V}$. We define
\begin{equation}\label{eqteststat_V_main} \widehat{\mathcal T}_V=T^{-1}{\langle\widehat{\mathcal K}\widehat v_V,\widehat v_V\rangle}/{\langle\widehat{\mathcal V}\widehat v_V,\widehat v_V\rangle}. \end{equation}
This statistic measures accumulated residual persistence along the direction of greatest contemporaneous residual variation, normalized by the residual variance in that direction. The statistic $\widehat{\mathcal T}_K$ depends on the scale of $\check U_t$. By contrast, $\widehat{\mathcal T}_V$ is invariant to rescaling of the residual process: replacing $\check U_t$ by $c\check U_t$ for any $c>0$ leaves it unchanged. Although the scale invariance of $\widehat{\mathcal T}_V$ can be convenient, both statistics are calibrated against their respective plug-in Monte Carlo null distributions and are reported as complementary diagnostics.
\item[Step 3: Interpreting the diagnostics and implementing the tests] \hfill
For practical interpretation, both diagnostics are right-tailed: large values indicate residual persistence inconsistent with between-cointegration. Under the null, $\widehat{\mathcal T}_K$ and $\widehat{\mathcal T}_V$ converge to finite limits; under the residual $I(1)$ alternative, they diverge at rates $T^2$ and $T$, respectively. The null limits are nonpivotal because they depend on the unknown long-run covariance structure and the estimation effect induced by projecting $Y_t$ on the integrated forcing vector. We therefore use the plug-in Monte Carlo procedure described in Section~\ref{sec_app_test_implementation} of the Supplement. At significance level $\alpha$, this procedure yields $\widehat q_{K,\alpha}$ and $\widehat q_{V,\alpha}$, the simulated $(1-\alpha)$ critical values for $\widehat{\mathcal T}_K$ and $\widehat{\mathcal T}_V$, respectively, together with the corresponding upper-tail $p$-values. The procedure uses kernel estimates of the relevant long-run and one-sided covariance operators, constructed from $\Delta\mathbf{x}_t$ and $\check U_t$, to simulate the joint Brownian processes entering the null limits. The tests reject when $\widehat{\mathcal T}_K>\widehat q_{K,\alpha}$ and $\widehat{\mathcal T}_V>\widehat q_{V,\alpha}$, respectively. Section~\ref{sec_app_test} of the Supplement establishes the validity of this approximation and the consistency of both tests.
\end{description}
\noindent Although these diagnostics use the partial-sum logic familiar from existing stationarity tests for functional time series, their target differs. They do not test whether $Y_t$ is stationary; rather, they test whether any stochastic trend remains after projecting $Y_t$ on the forcing vector. They thus provide a specification check for the long-run relation. When between-cointegration is not rejected, we estimate and conduct inference on the forcing-specific response functions under this maintained restriction. Rejection would indicate that the forcing vector does not account for all persistent variation in the temperature-anomaly distribution.
\subsection{Estimation under between-cointegration}\label{sec_esti}
\noindent We now turn to estimation under the maintained hypothesis of between-cointegration, continuing to assume that $\mathbf{x}_t$ has full nonstationary rank and hence no within-cointegration. This condition is assessed for the forcing data in Section~\ref{sec_test_stat}. Under between-cointegration, $U_t$ in \eqref{eqmodel1} is stationary and $B$ represents the long-run response operator. To remove the functional intercept $\beta_0$, write the model in demeaned form as
\begin{equation}\label{eqmodel1a}
Y_t-\bar Y_T=B(\mathbf{x}_t-\bar{\mathbf{x}}_T)+(U_t-\bar U_T),
\end{equation}
where $\bar Y_T$, $\bar{\mathbf{x}}_T$, and $\bar U_T$ denote the corresponding temporal sample means.
A natural starting point is the least-squares projection estimator $\check B=\widehat C_{Y\mathbf{x}}\widehat C_{\mathbf{x}\mathbf{x}}^{-1}$. The centered analogue of the result in Remark~\ref{rem_ls_limit} of the Supplement gives $\|\check B-B\|_{\operatorname{op}}=O_p(T^{-1})$. The corresponding asymptotic expansion of $T(\check B-B)$ contains nuisance-dependent bias terms arising from endogeneity between $U_t$ and $\Delta\mathbf{x}_t$ and serial dependence in their joint process. We therefore extend the fully modified least-squares estimator of \citet{phillips1995fully} to this setting, removing these terms and obtaining a limit suitable for feasible simulation-based inference.
\subsubsection{Computation of the proposed estimator}\label{sec_compest}
Formal assumptions, consistency results, and the limiting distribution of the estimator are established in Section~\ref{sec_app_est1} of the Supplement, while Section~\ref{sec_app_det2} details the intercept extension used here. We focus on its practical computation, which proceeds in three steps:
\begin{description}[style=nextline, leftmargin=1em, font=\bfseries,topsep=-2pt, itemsep=0pt]
\item[Step 1: Preliminary estimation and residuals] \hfill
We begin with the preliminary least-squares estimator $\check B$ and compute the projection residuals $\check U_t=(Y_t-\bar Y_T)-\check B(\mathbf{x}_t-\bar{\mathbf{x}}_T)$. These residuals are used in place of the unobserved disturbance $U_t$ when estimating the covariance operators for the fully modified correction.
\item[Step 2: Estimation of long-run covariance operators] \hfill
The fully modified estimator uses the long-run and one-sided long-run covariance operators $\Omega_{\mathbf{x}\mathbf{x}}=\sum_{j=-\infty}^{\infty}\mathbb E[\Delta\mathbf{x}_t\otimes\Delta\mathbf{x}_{t+j}]$, $\Omega_{U\mathbf{x}}=\sum_{j=-\infty}^{\infty}\mathbb E[\Delta\mathbf{x}_t\otimes U_{t+j}]$, $\Omega_{\mathbf{x}\mathbf{x}}^+=\sum_{j=0}^{\infty}\mathbb E[\Delta\mathbf{x}_t\otimes\Delta\mathbf{x}_{t+j}]$, and $\Omega_{U\mathbf{x}}^+=\sum_{j=0}^{\infty}\mathbb E[\Delta\mathbf{x}_t\otimes U_{t+j}]$. These operators determine the nuisance-dependent terms in the limiting distribution of $\check B$ and the corrections used by the fully modified estimator, as developed in Sections~\ref{sec_app_est1} and~\ref{sec_app_det2} of the Supplement. We estimate them using
\begin{align}
\widehat{\Omega}_{\mathbf{x}\mathbf{x}}&=\sum_{|j|\leq h}\mathrm{k}(j/h)\widehat{\Gamma}_{\mathbf{x}\mathbf{x}}^{(j)},&
\widehat{\Omega}_{U\mathbf{x}}&=\sum_{|j|\leq h}\mathrm{k}(j/h)\widehat{\Gamma}_{U\mathbf{x}}^{(j)},\label{eqsample1}\\
\widehat{\Omega}_{\mathbf{x}\mathbf{x}}^+&=\sum_{j=0}^{h}\mathrm{k}(j/h)\widehat{\Gamma}_{\mathbf{x}\mathbf{x}}^{(j)},&
\widehat{\Omega}_{U\mathbf{x}}^+&=\sum_{j=0}^{h}\mathrm{k}(j/h)\widehat{\Gamma}_{U\mathbf{x}}^{(j)}.\label{eqsample2}
\end{align}
Under our rank-one-operator convention, the lag-$j$ sample autocovariance and cross-covariance operators are $\widehat{\Gamma}_{\mathbf{x}\mathbf{x}}^{(j)}=T^{-1}\sum_{t=(1-j)\vee1}^{T\wedge(T-j)}\Delta\mathbf{x}_t\otimes\Delta\mathbf{x}_{t+j}$ and $\widehat{\Gamma}_{U\mathbf{x}}^{(j)}=T^{-1}\sum_{t=(1-j)\vee1}^{T\wedge(T-j)}\Delta\mathbf{x}_t\otimes\check U_{t+j}$, respectively, so $\widehat{\Gamma}_{U\mathbf{x}}^{(j)}$ maps $\mathcal H_{\mathbf{x}}$ to $\mathcal H_{Y}$. The function $\mathrm{k}(\cdot)$ is a kernel and $h$ is its bandwidth. Assumption~\ref{assum_test_kernel} of the Supplement gives the corresponding conditions. Proposition~\ref{prop1} establishes operator-norm consistency for the baseline specification, and Section~\ref{sec_app_det2} provides the demeaned extension used here. We use the Parzen kernel in the empirical analysis.
\item[Step 3: Computation of the proposed estimator] \hfill
Using the covariance estimators from Step~2, define the modified response and bias-correction operator as
\begin{equation}\label{eqest01_main}
Z_{1,t}=(Y_t-\bar Y_T)-\widehat{\Omega}_{U\mathbf{x}}\widehat{\Omega}_{\mathbf{x}\mathbf{x}}^{-1}\Delta\mathbf{x}_t,\qquad \widehat{\Upsilon}=\widehat{\Omega}_{U\mathbf{x}}^+-\widehat{\Omega}_{U\mathbf{x}}\widehat{\Omega}_{\mathbf{x}\mathbf{x}}^{-1}\widehat{\Omega}_{\mathbf{x}\mathbf{x}}^+.
\end{equation}
The modification of $Y_t-\bar Y_T$ removes the estimated long-run covariance between the disturbance component and the predictor innovations, while $\widehat{\Upsilon}$ corrects the remaining bias. Under the maintained full-rank condition, $\widehat{\Omega}_{\mathbf{x}\mathbf{x}}^{-1}$ is computed as an ordinary matrix inverse on the finite-dimensional space $\mathcal H_{\mathbf{x}}$. The proposed fully modified estimator is
\begin{equation}\label{eq_B_hat}
\widehat B=(\widehat C_{Z_1\mathbf{x}}-\widehat{\Upsilon})\widehat C_{\mathbf{x}\mathbf{x}}^{-1},\qquad \widehat C_{Z_1\mathbf{x}}=\frac{1}{T}\sum_{t=1}^T(\mathbf{x}_t-\bar{\mathbf{x}}_T)\otimes Z_{1,t}.
\end{equation}
For explicit computation, let $\widehat C_{\mathbf{x}\mathbf{x}}=\sum_{r=1}^{d_{\mathbf{x}}}\widehat\lambda_{\mathbf{x},r}\widehat v_{\mathbf{x},r}\otimes\widehat v_{\mathbf{x},r}$ be the eigendecomposition used in Section~\ref{sec_betcointeg}. Then, for any $v\in\mathcal H_{\mathbf{x}}$,
\begin{equation}\label{eq_B_hat_explicit}
\widehat B(v)=\sum_{r=1}^{d_{\mathbf{x}}}\widehat\lambda_{\mathbf{x},r}^{-1}\langle v,\widehat v_{\mathbf{x},r}\rangle\left\{\frac{1}{T}\sum_{t=1}^T\langle\mathbf{x}_t-\bar{\mathbf{x}}_T,\widehat v_{\mathbf{x},r}\rangle Z_{1,t}-\widehat{\Upsilon}(\widehat v_{\mathbf{x},r})\right\}.
\end{equation}
This expression parallels the preliminary least-squares projection in \eqref{eqprelimest}, with the modified response and bias correction incorporated explicitly. The forcing-specific response functions are obtained as $\widehat\beta_j=\widehat B(e_j)$, where $e_j\in\mathcal H_{\mathbf{x}}$ is the $j$th standard basis vector, with one in its $j$th coordinate and zeros elsewhere.
\end{description}
\noindent Theorem~\ref{thmapp2} in the Supplement establishes that $\|\widehat B-B\|_{\operatorname{op}}=O_p(T^{-1})$ under the specification with an intercept and derives the corresponding nonstandard limiting distribution. The next subsection uses a feasible simulation of this limit to conduct inference on local averages of the forcing-specific response functions.
\subsubsection{Simulation-based inference for local average responses}\label{sec_inference}
\noindent A direct confidence region for the full response operator $B$ would consist of maps from the forcing space $\mathcal H_{\mathbf{x}}$ into the CLR-response space $\mathcal H_{Y}$, while a coefficient-specific region for $\beta_j=B(e_j)$ lies in the infinite-dimensional space $\mathcal H_{Y}$. Neither is readily interpretable in applied work. Constructing a conventional uniform confidence band for $\beta_j$ would require additional mathematical structure and stronger assumptions. We therefore conduct inference on average CLR responses over selected regions of the temperature-anomaly support, which provide directly interpretable scalar summaries of the forcing-specific response functions.
\indent For a subinterval $[a_k,b_k]$ of the full temperature-anomaly support $[a,b]$, define the rectangular weight
\begin{equation}\label{eq_local_weight} w_k(s)=(b_k-a_k)^{-1}\mathbf{1}\{s\in[a_k,b_k]\}. \end{equation}
The weight itself need not be centered: for any $g\in\mathcal H_{Y}$, $\langle g,w_k\rangle=\langle g,w_k-(b-a)^{-1}\rangle$, so only its centered projection enters the inner product. Hence, $\langle\beta_j,w_k\rangle=(b_k-a_k)^{-1}\int_{a_k}^{b_k}\beta_j(s)\,ds$ is the average CLR response to forcing portfolio $j$ over the anomaly range $[a_k,b_k]$ and has a direct interpretation. Its estimator is $\langle\widehat\beta_j,w_k\rangle$. Let $e_j$ be the $j$th standard basis vector of $\mathcal H_{\mathbf{x}}$, so that $\beta_j=B(e_j)$. Applying Theorem~\ref{thmapp2} in the Supplement with $v=e_j$ and the centered projection of $w_k$, and using the preceding equality, gives
\begin{equation}\label{eq_local_limit} T\langle\widehat\beta_j-\beta_j,w_k\rangle\to_d\left\langle\left(\int_0^1 W^c_{\mathbf{x}}(r)\otimes dW_{U|\mathbf{x}}(r)\right)\left(\int_0^1 W^c_{\mathbf{x}}(r)\otimes W^c_{\mathbf{x}}(r)\,dr\right)^{-1}e_j,w_k\right\rangle. \end{equation}
Here $\to_d$ denotes convergence in distribution. The processes $W_{\mathbf{x}}$ and $W_U$ are the joint Brownian limits of the $T^{-1/2}$-scaled partial sums of $\Delta\mathbf{x}_t$ and $U_t$, taking values in $\mathcal H_{\mathbf{x}}$ and $\mathcal H_{Y}$, respectively. The demeaned path is $W^c_{\mathbf{x}}(r)=W_{\mathbf{x}}(r)-\int_0^1W_{\mathbf{x}}(u)\,du$, whereas $W_{U|\mathbf{x}}(r)=W_U(r)-\Omega_{U\mathbf{x}}\Omega_{\mathbf{x}\mathbf{x}}^{-1}W_{\mathbf{x}}(r)$ is the uncentered error Brownian motion independent of $W_{\mathbf{x}}$.
\indent Although the limiting distribution in \eqref{eq_local_limit} does not have closed-form quantiles, Theorem~\ref{thmapp3} in the Supplement shows that it can be consistently approximated using estimated covariance eigenelements and simulated scalar Brownian motions. Let $\widehat q_{\tau}(e_j,w_k)$ denote the simulated $\tau$ quantile of this approximation. An asymptotic $(1-\alpha)$ equal-tailed confidence interval for $\langle\beta_j,w_k\rangle$ is
\begin{equation}\label{eqlocalci} \operatorname{CI}_j(1-\alpha,w_k)=\left[\langle\widehat\beta_j,w_k\rangle-{\widehat q_{1-\alpha/2}(e_j,w_k)}/T,\ \langle\widehat\beta_j,w_k\rangle-{\widehat q_{\alpha/2}(e_j,w_k)}/{T}\right]. \end{equation}
\indent Repeating this calculation over a collection of subintervals produces a local-average confidence display for the forcing-specific response function \citep{seong2021functional,Nam2025}. It shows where the average CLR response over an anomaly range is individually distinguishable from zero at the specified confidence level. Each interval has asymptotic marginal coverage for its corresponding local average; without a multiplicity adjustment, the collection should not be interpreted as a simultaneous confidence band for the entire function.
\smallskip
\noindent\textbf{Monte Carlo implementation.} The feasible quantiles $\widehat q_{\tau}(e_j,w_k)$ are computed using the following three-step simulation. The asymptotic justification for this procedure is provided in Section~\ref{sec_app_det3} of the Supplement.
\begin{description}[style=nextline, leftmargin=1em, font=\bfseries,topsep=-2pt, itemsep=0pt]
\item[Step 1: Spectral decomposition of covariance operators] \hfill
Using the residuals $\check U_t$ defined above, we estimate their long-run covariance operator by
\begin{equation}\label{eq_lrv_UU} \widehat{\Omega}_{UU}=\sum_{|j|\leq h}\mathrm{k}(j/h)\widehat{\Gamma}_{UU}^{(j)},\qquad \widehat{\Gamma}_{UU}^{(j)}=\frac{1}{T}\sum_{t=(1-j)\vee1}^{T\wedge(T-j)}\check U_t\otimes\check U_{t+j}. \end{equation}
We then construct the conditional long-run covariance operator $\widehat{\Omega}_{U|\mathbf{x}}=\widehat{\Omega}_{UU}-\widehat{\Omega}_{U\mathbf{x}}\widehat{\Omega}_{\mathbf{x}\mathbf{x}}^{-1}\widehat{\Omega}_{\mathbf{x}U}$, where $\widehat{\Omega}_{\mathbf{x}U}=\widehat{\Omega}_{U\mathbf{x}}^*$. Let $(\widehat\lambda_j,\widehat v_j)$ and $(\widehat\mu_j,\widehat w_j)$ denote the eigenpairs of $\widehat{\Omega}_{U|\mathbf{x}}$ and $\widehat{\Omega}_{\mathbf{x}\mathbf{x}}$, respectively. We use the decompositions
\begin{equation}\label{eqdecomstep} \widehat{\Omega}_{U|\mathbf{x}}\approx\sum_{j=1}^{M}\widehat\lambda_j\widehat v_j\otimes\widehat v_j,\qquad \widehat{\Omega}_{\mathbf{x}\mathbf{x}}=\sum_{j=1}^{d_{\mathbf{x}}}\widehat\mu_j\widehat w_j\otimes\widehat w_j. \end{equation}
The first decomposition retains the leading $M$ eigenpairs of an operator acting on the infinite-dimensional space $\mathcal H_{Y}$, whereas the second requires no truncation because $\mathcal H_{\mathbf{x}}$ has dimension $d_{\mathbf{x}}$. These eigenelements are obtained by functional eigenanalysis and matrix eigendecomposition, respectively. In finite samples, $M$ is chosen to capture the dominant variation. Asymptotically, $M=M_T$ increases slowly enough with $T$ to control eigenelement estimation error, as required by Theorem~\ref{thmapp3} of the Supplement.
\item[Step 2: Generation of synthetic Brownian paths] \hfill
For each Monte Carlo replication $i=1,\ldots,R_{\mathrm{MC}}$, generate standard scalar Brownian motions $\{W_{1,j,(i)}\}_{j=1}^{M}$ and $\{W_{2,j,(i)}\}_{j=1}^{d_{\mathbf{x}}}$, independently across the two families and across replications, on a fine grid over $[0,1]$. To reproduce the demeaning of the integrated forcing variables, define $W^c_{2,j,(i)}(r)=W_{2,j,(i)}(r)-\int_0^1W_{2,j,(i)}(u)\,du$. The required synthetic paths are
\begin{equation}\label{eq_synthetic_BM} \widehat W_{U|\mathbf{x},(i)}(r)=\sum_{j=1}^{M}\widehat\lambda_j^{1/2}\widehat v_jW_{1,j,(i)}(r),\qquad \widehat W^c_{\mathbf{x},(i)}(r)=\sum_{j=1}^{d_{\mathbf{x}}}\widehat\mu_j^{1/2}\widehat w_jW^c_{2,j,(i)}(r). \end{equation}
Only the forcing path $\widehat W^c_{\mathbf{x},(i)}$ is centered. No analogous centering is applied to $\widehat W_{U|\mathbf{x},(i)}$.
\item[Step 3: Monte Carlo approximation of the quantiles] \hfill
For a given anomaly range with weight $w_k$, define the scalar Brownian process $\widehat Z_{k,(i)}(r)=\langle\widehat W_{U|\mathbf{x},(i)}(r),w_k\rangle$, whose increments satisfy
\begin{equation}\label{eq_simulated_Z} d\widehat Z_{k,(i)}(r)=\sum_{\ell=1}^{M}\widehat\lambda_\ell^{1/2}\langle\widehat v_\ell,w_k\rangle\,dW_{1,\ell,(i)}(r). \end{equation}
Because $\mathcal H_{\mathbf{x}}=\mathbb R^{d_{\mathbf{x}}}$, the operator $\int_0^1\widehat W^c_{\mathbf{x},(i)}(r)\otimes\widehat W^c_{\mathbf{x},(i)}(r)\,dr$ is represented by an ordinary $d_{\mathbf{x}}\timesd_{\mathbf{x}}$ matrix. The simulated realization of the limiting random variable in \eqref{eq_local_limit} is therefore computed directly as
\begin{equation}\label{eq_simulated_local_limit} e_j'\left\{\int_0^1\widehat W^c_{\mathbf{x},(i)}(r)\widehat W^c_{\mathbf{x},(i)}(r)'\,dr\right\}^{-1}\int_0^1\widehat W^c_{\mathbf{x},(i)}(r)\,d\widehat Z_{k,(i)}(r). \end{equation}
The ordinary integrals are evaluated by Riemann sums on the simulation grid, and the stochastic integral is evaluated using Brownian increments and left-endpoint values of $\widehat W^c_{\mathbf{x},(i)}$. Repeating Steps~2 and~3 for $i=1,\ldots,R_{\mathrm{MC}}$ gives the empirical distribution of \eqref{eq_simulated_local_limit}. Its empirical $\tau$ quantile is $\widehat q_{\tau}(e_j,w_k)$; in particular, $\widehat q_{\alpha/2}(e_j,w_k)$ and $\widehat q_{1-\alpha/2}(e_j,w_k)$ are the two quantiles used in \eqref{eqlocalci}.
\end{description}
\section{Empirical Analysis of Distributional Responses}
\label{sec:distributional_responses}
\noindent This section first examines whether the CO$_2$ and non-CO$_2$ forcing portfolios carry distinct stochastic trends (no within-cointegration) and then whether their trends account for all stochastic trends in the anomaly-density process (between-cointegration). It then estimates the forcing-specific responses, uses their implied densities to describe changes in location, dispersion, and tail mass, examines the common-response restriction imposed by an aggregate anthropogenic forcing index, compares the fitted responses with matched scalar benchmarks, and localizes probability-mass reallocations across the anomaly support. Section~\ref{sec_app_coverage} of the Supplement examines whether the results are sensitive to changes in the spatial coverage of the temperature data over time.
\subsection{Anthropogenic forcing and persistent climate change}
\label{sec_test_stat}
The empirical specification involves two distinct questions, assessed in sequence. We first examine whether the forcing vector $\mathbf{x}_t=(\mathrm{F1}_t,\mathrm{F2}_t)'$ contains two distinct stochastic trends. If F1 and F2 are within-cointegrated, their persistent variation is driven by fewer than two stochastic trends, and they cannot be treated as distinct persistent coordinates. Given full nonstationary rank of $\mathbf{x}_t$, we then test whether its stochastic trends span all stochastic trends in the density process $Y_t$, leaving the disturbance $U_t$ in \eqref{eqmodel1} stationary.
Let $d_N(\mathbf{x}_t)$ denote the number of stochastic trends in $\mathbf{x}_t$, equivalently, the dimension of its nonstationary subspace. Since $\mathbf{x}_t$ is a bivariate vector with $d_{\mathbf{x}}=2$, we apply the variance-ratio diagnostic of \citet{Breitung2002} directly to its demeaned observations. We test the null of no within-cointegration against the alternative of fewer than two stochastic trends:
\begin{equation*} H_0:\ d_N(\mathbf{x}_t)=2,\qquad H_1:\ d_N(\mathbf{x}_t)<2.
\end{equation*}
Under the null, no nontrivial linear combination of F1 and F2 is stationary. The variance-ratio statistic is $63.2050$, with a $p$-value of approximately $0.95$, far above conventional significance levels. This supports treating F1 and F2 as distinct persistent coordinates. The maintained full-nonstationary-rank condition also ensures that $\widehat C_{\mathbf{x}\mathbf{x}}$ and $\widehat\Omega_{\mathbf{x}\mathbf{x}}$ are invertible with probability approaching one, as required to define $\check B$ and $\widehat B$. It does not, however, imply that F1 and F2 have different distributional response functions; the corresponding common-response restriction is examined below.
Having retained the specification in which the CO$_2$ and non-CO$_2$ forcing portfolios carry linearly independent stochastic trends, we next ask whether these trends account for the persistent evolution of the density process. In our framework, this is the between-cointegration condition: after projecting the density process on the forcing vector, the remaining disturbance $U_t$ in \eqref{eqmodel1} must be stationary. We thus construct the residuals $\check U_t=(Y_t-\bar Y_T)-\check B(\mathbf{x}_t-\bar{\mathbf{x}}_T)$ from the centered least-squares projection, where $\check B=\widehat C_{Y\mathbf{x}}\widehat C_{\mathbf{x}\mathbf{x}}^{-1}$, and first evaluate the unnormalized diagnostic $\widehat{\mathcal T}_K$ defined in \eqref{eqteststat_K_main}. We also use the covariance-normalized diagnostic $\widehat{\mathcal T}_V$ in \eqref{eqteststat_V_main} as a complementary check. Critical values are obtained from the plug-in Monte Carlo procedure described in Section~\ref{sec_betcointeg}. The relevant long-run and one-sided covariance operators are estimated from $(\Delta\mathbf{x}_t,\check U_t)$ using the Parzen kernel with the bandwidth chosen as the nearest integer to $T^{1/4}$, and the two null distributions are approximated using 5{,}000 replications on a grid of 499 subintervals. The unnormalized statistic is $\widehat{\mathcal T}_K=0.0110$, with a simulated 5\% critical value of $0.0165$ and a Monte Carlo $p$-value of $0.1398$. The complementary covariance-normalized check gives $\widehat{\mathcal T}_V=0.2651$, with a simulated 5\% critical value of $0.3114$ and a Monte Carlo $p$-value of $0.0748$. Both statistics fall below their respective critical values, so neither test rejects between-cointegration at the 5\% level. The response analysis thus proceeds under the specification in \eqref{eqmodel1}.
\noindent Within the stochastic-trend framework adopted here, this finding also has a substantive climate interpretation. Our empirical formulation represents climate change as persistent evolution in the cross-sectional distribution of observed temperature anomalies. At the 5\% level, the tests do not detect an additional stochastic trend in the residual after accounting for the two anthropogenic forcing portfolios. Thus, conditional on the constructed forcing measures and the maintained long-run specification, the results are consistent with their trends spanning the persistent distributional component of observed climate change. This interpretation concerns only the persistent component; stationary natural variability and transitory shocks may remain in the residual.
\subsection{Forcing-specific responses and what aggregation misses}
\label{subsec:responses_and_margins}
\noindent Having retained F1 and F2 as distinct persistent coordinates and found no evidence against between-cointegration, we estimate the long-run response operator using the fully modified estimator in Section~\ref{sec_esti}. The forcing-specific estimates are $\widehat\beta_k=\widehat B(e_k)$ for $k\in\{1,2\}$. For illustration, Figure~\ref{Fig:RF_Response} reports $\Delta x_k\widehat\beta_k$, where $\Delta x_k$ is set equal to the sample standard deviation of forcing portfolio $k$ and serves only as a convenient reporting scale, together with marginal confidence intervals for the corresponding local-average CLR responses.
Because $\widehat\beta_k$ belongs to the CLR space, it describes a response on the centered log-density scale rather than directly on the density scale. To illustrate the corresponding change in density, we apply the fitted CLR response for the chosen increment $\Delta x_k$ to a common reference density and map the result back through the inverse CLR transformation. Following \citet{Nam2025}, we take the arithmetic average $f_0(s)=T^{-1}\sum_{t=1}^Tf_t(s)$ on $[a,b]$, where $a=-5.2527$ and $b=5.6224$, as the reference. Because the density process is nonstationary, $f_0$ is not interpreted as a stationary population mean or long-run equilibrium, but only as a representative data-based reference. Its choice affects the density-scale illustration and descriptive margins, but not the CLR response estimates or their local-average inference. The implied probability density following a change $\Delta x_k$ is
\begin{equation}\label{eq:clr_pushforward}
f_k(\,\cdot\,;\Delta x_k)=\operatorname{clr}^{-1}\!\left[\operatorname{clr}(f_0)+\Delta x_k\widehat\beta_k\right].
\end{equation}
In the finite-basis implementation used here, the inverse-CLR normalization ensures that $f_k(s;\Delta x_k)>0$ and $\int_a^b f_k(s;\Delta x_k)\,ds=1$, so $f_k(\cdot;\Delta x_k)$ is a proper probability density on $[a,b]$.
\begin{figure}[t]
\centering
\includegraphics[height=0.27\textwidth,trim={0.1cm 0.1cm 0.1cm 0.1cm},clip]{Figures/F1_Response0.png}
\includegraphics[height=0.24\textwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{Figures/Figure2_F1_density_95CI.png}
\includegraphics[height=0.27\textwidth,trim={0.1cm 0.1cm 0.1cm 0.1cm},clip]{Figures/F2_Response0.png}
\includegraphics[height=0.24\textwidth,trim={0.0cm 0.0cm 0.0cm 0.0cm},clip]{Figures/Figure2_F2_density_95CI.png}
\caption{Estimated responses to the specified increases in F1 (top) and F2 (bottom): local-average CLR responses with marginal 95\% confidence intervals (left), and the common reference density and fitted implied density (right). Shaded bands show central 95\% pointwise simulation envelopes for the implied densities, based on 50{,}000 joint draws; see Section~\ref{sec_app_density_ci} of the Supplement.}
\label{Fig:RF_Response}
\end{figure}
Figure~\ref{Fig:RF_Response} reports the local-average CLR responses scaled by the increments $\Delta x_k$ in the left panels and the implied densities alongside the reference density in the right panels.\footnote{We use an 11-grid-point window, spanning approximately $0.86^\circ\mathrm{C}$ at interior points, to smooth grid-level variation while retaining local response features.} The CLR panels use different vertical scales. At this scale, the implied F1 density shifts probability mass from colder to warmer anomalies while increasing central concentration and upper-tail mass relative to the reference. The implied F2 density shows a smaller warmward shift, greater dispersion, and a modest increase in upper-tail mass. For F1, the marginal intervals for the local-average CLR response lie entirely below zero over a broad cold-anomaly range and entirely above zero over a warm-anomaly range. For F2, they lie entirely below zero over a central range and entirely above zero over a broad range on the warm side.
\indent We next report several descriptive summaries of the implied density changes shown in the right panels of Figure~\ref{Fig:RF_Response}. These quantities summarize location, dispersion, and tail features and are not treated as separate inferential targets. The distinction between location and dispersion is empirically relevant and is consistent with \citet{chang2020evaluating}, who document stochastic trends in both features.
Let $Q_k(u)$ denote the quantile function of $f_k(\cdot;\Delta x_k)$, and let $Q_0(u)$ denote that of the reference density.\footnote{In practice, $Q_k(u)$ is obtained from the CDF evaluated on the discretized anomaly grid. We numerically enforce monotonicity, retain unique CDF values, and use linear interpolation for inversion.} The implied change in the mean is $\Delta\mu_k=\int_a^b s\{f_k(s;\Delta x_k)-f_0(s)\}\,ds$, and the central-half displacement is $\delta_k^{\mathrm{loc}}=(u_H-u_L)^{-1}\int_{u_L}^{u_H}\{Q_k(u)-Q_0(u)\}\,du$, where $(u_L,u_H)=(0.25,0.75)$. This quantity provides a location benchmark for typical anomaly states that is less sensitive to tail behavior than the mean. We also report the median shift $\Delta\mathrm{Median}_k=Q_k(0.5)-Q_0(0.5)$, the interquartile-range change $\Delta\mathrm{IQR}_k=\{Q_k(0.75)-Q_k(0.25)\}-\{Q_0(0.75)-Q_0(0.25)\}$, and the mean--location gap $\Delta\mu_k-\delta_k^{\mathrm{loc}}$, which compares the overall mean change with the central-half displacement.
For the tail summaries, fix the reference thresholds $r_q=Q_0(q)$ for $q\in\{0.05,0.95\}$. To distinguish reallocation between the two tails from a change in combined tail probability, define the cold-tail change as $\Delta\mathrm{ColdMass}_k=\int_a^{r_{0.05}}\{f_k(s;\Delta x_k)-f_0(s)\}\,ds$, the warm-tail change as $\Delta\mathrm{WarmMass}_k=\int_{r_{0.95}}^b\{f_k(s;\Delta x_k)-f_0(s)\}\,ds$, and their sum as $\Delta\mathrm{Extreme}_k=\Delta\mathrm{ColdMass}_k+\Delta\mathrm{WarmMass}_k$.
\begin{table}[t]
\centering
\caption{Descriptive summaries of the fitted density changes under one-standard-deviation increases. Temperature margins are in $^\circ$C; tail margins are probability changes.}
\label{Tab:MarginDecomp}
\setlength{\tabcolsep}{3pt}
\resizebox{\textwidth}{!}{
\begin{tabular}{lrrrrrrrr}
\toprule
Forcing
& $\Delta\mu$
& $\delta^{\mathrm{loc}}$
& $\Delta\mu-\delta^{\mathrm{loc}}$
& $\Delta\mathrm{Median}$
& $\Delta\mathrm{IQR}$
& $\Delta\mathrm{ColdMass}$
& $\Delta\mathrm{WarmMass}$
& $\Delta\mathrm{Extreme}$ \\
\midrule
F1 & 0.3125 & 0.2646 & 0.0479 & 0.2562 & $-$0.1120
& $-$0.0205 & 0.0223 & 0.0018 \\
F2 & 0.0789 & 0.0883 & $-$0.0094 & 0.0894 & 0.1500
& 0.0003 & 0.0080 & 0.0083 \\
\bottomrule
\end{tabular}}
\end{table}
\indent Table~\ref{Tab:MarginDecomp} provides numerical summaries of the changes in the implied densities shown in Figure~\ref{Fig:RF_Response}. We use their signs and relative magnitudes descriptively, not as separately tested quantities. At the chosen reporting scale, the larger central-half displacement and negative IQR point estimate for F1 correspond to its stronger warmward shift and greater central concentration, whereas the smaller displacement and positive IQR point estimate for F2 correspond to its weaker location shift and wider central range.
The table shows two further distinctions. First, the implied mean--location gap is positive for F1 but close to zero for F2. Thus, F1's mean change is not fully summarized by its central displacement, whereas the implied mean and central shifts for F2 are closely aligned. Second, the sum of the two reference-tail changes distinguishes changes in extreme-state composition from changes in combined tail probability. For F1, the cold-tail decrease nearly offsets the warm-tail increase, leaving $\Delta\mathrm{Extreme}_1=0.0018$; the implied tail composition therefore shifts from the cold to the warm tail with little change in combined tail probability. For F2, the cold tail is nearly unchanged, so the warm-tail increase carries through to $\Delta\mathrm{Extreme}_2=0.0083$ and is concentrated on the warm side.
\indent The estimated CLR profiles in Figure~\ref{Fig:RF_Response} differ visibly in shape and sign patterns, suggesting that aggregation may conceal substantial response heterogeneity. As a formal check using the projection inference developed above, we examine the common-response restriction $\beta_1=\beta_2$ introduced in Section~\ref{subsec:aggregation_cost} within the fixed response space used for estimation. This comparison differs from the no-within-cointegration diagnostic in Section~\ref{sec_test_stat}: that diagnostic assesses the number of stochastic trends in the forcing block, whereas the present comparison asks whether F1 and F2 have the same distributional response.
Because F1 and F2 are measured in the same $\mathrm{W\,m^{-2}}$ units, equality is assessed using the unscaled coefficient functions. Let $\mathbf c=e_1-e_2$ denote the contrast vector, so that $B\mathbf c=\beta_1-\beta_2$. Section~\ref{sec_app_common_response} of the Supplement uses the Cram\'{e}r--Wold device and Theorems~\ref{thmapp2}--\ref{thmapp3} to establish joint validity of the feasible approximation for any fixed finite collection of response directions. Let $\{\phi_m\}_{m=1}^J$, with $J=20$, denote the orthonormal nonconstant Fourier directions spanning the empirical response space, and let $P_J$ denote the associated projection. We test
\begin{equation*}
H_{0,J}:\ \langle B\mathbf c,\phi_m\rangle=0\ \text{for }m=1,\ldots,J,\qquad
\widehat{\mathcal S}_{\mathrm{eq},J}
=T\left\{\sum_{m=1}^J\langle\widehat B\mathbf c,\phi_m\rangle^2\right\}^{1/2}
=T\|P_J\widehat B\mathbf c\|_{L^2}.
\end{equation*}
Here $\|g\|_{L^2}=(\int_a^b |g(s)|^2\,ds)^{1/2}$ denotes the $L^2[a,b]$ norm. The null distribution is approximated using 50{,}000 joint draws from the feasible approximation to the intercept-model limit, with the same Brownian paths used across all $J$ coordinates within each replication and with $M=7$. The same section reports the numerical implementation, retained-variation shares, and sensitivity to $J$ and $M$. The observed statistic is 557.48, exceeding the simulated 95th percentile of 463.85; the corresponding add-one Monte Carlo $p$-value is 0.0199. Thus, the common-response restriction is rejected at the 5\% level within the fixed response space.
\indent We next use two scalar-response benchmarks to distinguish the information lost by reducing the density response to its mean from that lost by aggregating the forcing coordinates. The aggregate scalar specification regresses the density mean on total anthropogenic forcing, thereby collapsing both the response and predictor. The vector-to-mean specification retains F1 and F2 separately but models only the mean response. All three specifications use the same fully modified estimation approach.
For Table~\ref{Tab:BenchmarkComparison}, we construct a matched reporting scenario by setting the increase in aggregate anthropogenic forcing (F1+F2) equal to its sample standard deviation and allocating it between F1 and F2 in proportion to their sample standard deviations. This scaling serves only as a convenient empirical normalization.\footnote{Let $s_1$ and $s_2$ denote the sample standard deviations of F1 and F2, and let $s_{\mathrm{agg}}$ denote that of F1+F2. We set $\Delta_j=s_{\mathrm{agg}}s_j/(s_1+s_2)$ for $j=1,2$, so that $\Delta_1+\Delta_2=s_{\mathrm{agg}}$. The F1, F2, and joint rows apply $(\Delta_1,0)$, $(0,\Delta_2)$, and $(\Delta_1,\Delta_2)$, respectively.} The component rows isolate the two allocated changes, whereas the joint row combines them to represent the matched aggregate increase. This normalization differs from the coordinate-specific scenarios in Figure~\ref{Fig:RF_Response} and Table~\ref{Tab:MarginDecomp}, which change F1 and F2 separately by their respective sample standard deviations. For the matched joint change under the vector-to-density specification, the implied probability density is $f(\,\cdot\,;\Delta_1,\Delta_2)=\operatorname{clr}^{-1}[\operatorname{clr}(f_0)+\Delta_1\widehat\beta_1+\Delta_2\widehat\beta_2]$. Because the inverse CLR mapping is nonlinear, the density-derived component margins need not sum to the joint margin.
\begin{table}[t]
\centering
\caption{Matched responses to a one-standard-deviation increase in aggregate
anthropogenic forcing. Temperature margins are in $^\circ$C; WarmMass is a probability change.}
\label{Tab:BenchmarkComparison}
\setlength{\tabcolsep}{7pt}
\begin{tabular}{llrrr}
\toprule
Specification & Channel & $\Delta\mu$ & $\Delta\mathrm{IQR}$ & $\Delta\mathrm{WarmMass}$ \\
\midrule
Aggregate scalar & F1+F2 & 0.3430 & \multicolumn{1}{c}{--} & \multicolumn{1}{c}{--} \\
\midrule
Vector-to-mean & Matched F1 component & 0.2785 & \multicolumn{1}{c}{--} & \multicolumn{1}{c}{--} \\
& Matched F2 component & 0.0604 & \multicolumn{1}{c}{--} & \multicolumn{1}{c}{--} \\
& Matched joint change & 0.3389 & \multicolumn{1}{c}{--} & \multicolumn{1}{c}{--} \\
\midrule
Vector-to-density & Matched F1 component & 0.2874 & $-$0.1030 & 0.0202 \\
& Matched F2 component & 0.0716 & 0.1371 & 0.0073 \\
& Matched joint change & 0.3678 & $-$0.0098 & 0.0283 \\
\bottomrule
\end{tabular}
\end{table}
\indent
Table~\ref{Tab:BenchmarkComparison} provides a descriptive comparison of fitted values; no separate inference is conducted for its nonlinear margins. The $\Delta\mathrm{IQR}$ and $\Delta\mathrm{WarmMass}$ entries summarize features of the implied densities. The matched mean responses are $0.3430^\circ$C under the aggregate scalar specification and $0.3389^\circ$C under the vector-to-mean specification, compared with $0.3678^\circ$C under the vector-to-density specification. Thus, broadly similar fitted mean responses coexist with IQR and tail changes observable only under the density specification. The matched F1 and F2 components have opposite-signed IQR changes. When the two changes are applied jointly, the resulting density exhibits little change in central spread, with $\Delta\mathrm{IQR}=-0.0098^\circ\mathrm{C}$, while upper-tail probability increases by 2.83 percentage points. The joint implied density exhibits a meaningful change in tail behavior despite its small IQR change.
\subsection{Distributional displacement and reshaping}
\label{subsec:loss_decomposition_interpretation}
\noindent The global margins in Table~\ref{Tab:MarginDecomp} summarize changes in location, dispersion, and tail mass, but not where offsetting gains and losses of probability mass occur. Similar mean or IQR changes can reflect quite different reallocations across the anomaly support. Pointwise density differences are also less directly interpretable because they measure density height rather than probability and may be sensitive to grid-level variation. We therefore examine probability changes over fixed-width neighborhoods, providing an intermediate view between global summaries and pointwise comparisons.
We return to the coordinate-specific reporting scenarios and write $f_k^{\mathrm{full}}(s)\equiv f_k(s;\Delta x_k)$ for the full implied density in \eqref{eq:clr_pushforward}. To distinguish displacement from reshaping, we construct a location-only benchmark by shifting the reference density by the central-half displacement:
\begin{equation}\label{eq:loc_counterfactual_sec6}
f_k^{\mathrm{loc}}(s)\propto f_0\!\left(s-\delta_k^{\mathrm{loc}}\right),\qquad s\in[a,b].
\end{equation}
After extending $f_0$ by zero outside $[a,b]$, the shifted density is restricted to $[a,b]$ and renormalized. Apart from this boundary adjustment, the benchmark preserves the reference shape. The difference between $f_k^{\mathrm{full}}$ and $f_k^{\mathrm{loc}}$ therefore captures implied changes in dispersion, asymmetry, and mass allocation beyond a pure displacement.
For window width $\ell>0$ and center $r\in[a+\ell/2,b-\ell/2]$, let $I_\ell(r)=(r-\ell/2,r+\ell/2]$ and define
\begin{equation}\label{eq:local_loss_functional_sec6}
L_\ell(r;f)=\int_{I_\ell(r)}f(s)\,ds.
\end{equation}
Unlike density height, $L_\ell(r;f)$ is the probability assigned to a neighborhood of $r$. Define $\Delta L_{k,\ell}^{\mathrm{full}}(r)=L_\ell(r;f_k^{\mathrm{full}})-L_\ell(r;f_0)$, $\Delta L_{k,\ell}^{\mathrm{loc}}(r)=L_\ell(r;f_k^{\mathrm{loc}})-L_\ell(r;f_0)$, and $\Delta L_{k,\ell}^{\mathrm{dist}}(r)=\Delta L_{k,\ell}^{\mathrm{full}}(r)-\Delta L_{k,\ell}^{\mathrm{loc}}(r)$. Thus, the full local response is decomposed exactly into displacement and residual reshaping at each $r$. This decomposition is descriptive and benchmark-dependent rather than a unique structural separation. We set $\ell=0.50^\circ\mathrm{C}$ and evaluate the profiles on an overlapping grid with a $0.05^\circ\mathrm{C}$ step. Positive values indicate gains in local probability mass and negative values indicate losses. Because the windows overlap, the profiles are localized diagnostics rather than an additive probability partition.
\begin{figure}[t]
\centering
\includegraphics[height=0.35\textwidth, width=0.49\textwidth]{Figures/Damage1.png}
\includegraphics[height=0.35\textwidth, width=0.49\textwidth]{Figures/Damage2.png}
\caption{Local probability-mass responses for F1 (left) and F2 (right): full fitted responses, location-only components, and residual reshaping components. Shaded bands are central 95\% pointwise simulation envelopes for the full responses, based on the same 50{,}000 joint draws as the implied-density panels of Figure~\ref{Fig:RF_Response}; see Section~\ref{sec_app_density_ci} of the Supplement.}
\label{Fig:BinLossDecomp}
\end{figure}
\indent
Figure~\ref{Fig:BinLossDecomp} reveals different mixtures of location-only displacement and residual reshaping. For F1, the location-only component determines the broad sign pattern: fitted probability mass falls over cold and near-central states and rises on the warm side. The reshaping component adds mass around $r\approx0$ and removes it near $r\approx1.5$, consistent with the negative IQR change reported in Table~\ref{Tab:MarginDecomp}. For F2, the location-only component is relatively small, whereas the reshaping component removes mass near the center and adds it on both flanks, around $r\approx-1$ and $r\approx1.5$, providing the localized counterpart of its positive IQR change. In distributional terms, the fitted CO$_2$ response is therefore dominated by a broad warmward displacement of the cross-sectional distribution of observed local temperature anomalies. By contrast, the fitted non-CO$_2$ response exhibits more pronounced reshaping, with greater dispersion and a higher relative prevalence of both cool and warm off-center states, consistent with the heterogeneous signs and spatial incidence of the components in this portfolio.
\section{The shape of warming: concluding remarks}\label{conclude}
\noindent By linking separate anthropogenic forcing portfolios to the anomaly density rather than its mean alone, our analysis gives the \emph{shape of warming} a specific empirical meaning. The specification evidence supports retaining the CO$_2$ and non-CO$_2$ portfolios as distinct persistent coordinates and is consistent with their trends jointly spanning the persistent evolution of the observed temperature-anomaly distribution. Within the fixed response space, rejection of the common-response restriction further indicates that the distributional response depends on the composition, not merely the aggregate level, of anthropogenic forcing.
\indent On the density scale, warming is not simply a uniform shift to the right. For the CO$_2$ portfolio, the fitted density moves warmward and becomes more concentrated around its center: colder states lose probability, the central peak sharpens, and warmer states gain probability, while the two tail changes nearly offset. For the non-CO$_2$ portfolio, the density moves less but spreads outward: central states lose mass, both flanks gain it, and additional extreme-state probability appears mainly in the warm tail. The two portfolios therefore trace distinct shapes of warming: displacement with concentration versus weaker displacement with widening and reshaping. Scalar benchmarks yield similar aggregate mean responses but compress these patterns into a single number. These fitted patterns represent persistent reduced-form associations and do not establish causality.
\bibliography{VtD_biblio}