\section{Introduction}
\label{section:introduction}

    
Modern data science techniques like deep neural networks and ensemble methods such as gradient boosting are well-known for yielding high prediction accuracy. However, it is barely understood how the output depends on the input variables for these techniques.
Consequently, the lack of transparency within these \textit{black boxes} is limiting the practical use of machine learning (ML), since decision-makers require profound reasoning for critical applications such as medical diagnoses, credit scoring or price setting.
One way of navigating this issue is to estimate the importance of the variables at hand.


A variable importance measure that recently has gained in popularity is the Shapley value. Shapley values enable the statistician to explain predictions on a local level, where each variable is assigned an importance measure in explaining the dependent variable for a particular observation. The idea of Shapley values originates from game theory and was proposed for giving a unique solution to a cooperative game \citep{shapley1953s}. Hereby, the goal is to assign credit to each player, according to their contribution to a total payoff in the cooperative game. An analogy is found in statistics, where the interest is in assigning an importance measure for each variable. For example, \citet{lundberg2017unified} define the Shapley values to explain predictions of the conditional mean. Other works make use of this flexibility by defining their version in terms of predictiveness measures such as the variance of the conditional mean \citep{benard2021shaff, owen2017shapley} or even $R^2$ and mean-squared-error \citep{williamson2020efficient}. 



Previous literature is mainly concerned with the practical computation of Shapley values, particularly focusing on computational efficiency \citep{benard2021shaff, lundberg2020local}.
To the best of our knowledge, barely any paper approached the topic from a nonparametric perspective, in which it is possible to deduce precise asymptotic results for variable importance measures. We are aware that an asymptotic result is not in the immediate interest of a ``high prediction accuracy" data scientist.
However, we seek to fill this gap of asymptotics in the interest of explainable ML by considering $d$-dimensional \textit{Shapley curves} instead of point-wise Shapley values, thereby convergence of the whole curve is established. In contrast to the previous literature, the Shapley curves result in a variable importance measure for any arbitrary evaluation point in the $d$-dimensional input space.


In order to establish asymptotic results we a priori define the Shapley curves on population level. It is uniquely determined by the $d$-dimensional conditional mean function and the joint distribution of the covariates. We propose two conceptually different estimation approaches for these curves. One way to obtain Shapley curves is by the direct estimation of the required components based on the corresponding subset of variables. We call this approach \textit{component-based}. In this work, we propose a different approach, inspired by the marginal integration literature for additive models \citep{linton1995kernel, linton1997analysis}, as means to consistently estimate Shapley curves in a more general setting. 
The marginal integration estimator is constructed by integrating out the variables of no-interest using a pilot estimator for the conditional mean function over the empirical distribution of covariates. Naively applying this approach is not desirable since it leads to inconsistent estimation of Shapley curves under variable dependence. By using estimates of the conditional distribution of the covariates instead of the marginal one, we achieve consistency even in the presence of dependence among the variables. We call this approach \textit{integration-based}. 

Our contribution is the following. First, we show that both estimation techniques, the component-based and the integration-based, achieve the minimax convergence rate estimating the Shapley curves. Next, asymptotic normality of Shapley curves is established for both approaches. However, these two estimators differ in terms of their asymptotic properties. The integration-based approach is accompanied by the advantage of only requiring a single pilot estimator. While we show that the asymptotic variance is the same, the bias of the integration-based estimator for the Shapley curve inflates. The general issue is that the reliance on a $d$-dimensional pilot estimator leads to oversmoothing. We extrapolate this technical device to other (modern) smoothing techniques such as forest-based estimators and deep neural networks. While the component-based approach is practically infeasible in a variable-rich environment, the integration-based approach circumvents this computational remedy by only relying on a single $d$-dimensional model fit. This is particularly interesting for modern data scientists working under practical constraints. The said inflated bias of the integration-based approach provides a practical trade-off to the quick explanation of a complex algorithm.
Further, we derive a one-dimensional convergence rate if the true process is an additive model and variables are independent. Our theoretical arguments are empirically confirmed in an extensive simulation exercise.

In a real data application, we estimate Shapley curves for vehicle prices as a function of vehicle characteristics. The estimated confidence intervals for the Shapley curves allow the practitioner to conduct statistical inference. In particular, it enables to distinguish between significant and non-significant intervals of the variables. The application of the said variable importance measure allows price setters to apply a data-driven price-setting alternative to common business standards such as markup calculation or surveying. 
    
    
    
    
    
    
    
    
A general overview of state-of-the art methods and current challenges in the vibrant field of interpretable machine learning is given in \citet{molnar2020interpretable}. 
In contrary to the local SHAP method \citep{lundberg2017unified}, \citet{covert2020understanding} provide the SAGE approach, which seeks to find a global variable importance measure for the covariates. \citet{lundberg2020local} provide a local tree ensemble specific Shapley method.
Further, there exist a few closely related publications having a statistical perspective with accompanied asymptotic arguments. \citet{williamson2020efficient} define \textit{Shapley effects} as a global predictiveness measure and show that it is asymptotically normal. 
\citet{benard2021shaff} show consistency for Shapley effects, where the variance of the conditional mean is used with a random forest as an estimator. 
\citet{verdinelli2021decorrelated} focus on solving a problem of correlation bias with a variable importance measure based on dropping covariates from a semiparametric regression model.
Going further, the integration-based estimation approach of our work has its root in the marginal integration literature for additive models. Major contributions have been made by \citet{auestad1991functional, tjostheim1994nonparametric}, \citet{linton1995kernel} and \citet{newey1994kernel}. \citet{linton1995kernel} use marginal integration to estimate the partial dependence function of a two-dimensional additive and multiplicative model by integration of an initial pilot estimator of the conditional mean.
\citet{chen1996nonparametric} extend the asymptotics for additive models while establishing results for dimensions larger than $2$.
Further, \citet{nielsen1998optimization} give insights into smooth backfitting and marginal integration of additive models. \citet{mammen1999existence} derive the corresponding asymptotic distribution for the backfitting estimation technique. 
Moreover, Shapley-related ideas and the marginal integration approach have also found application in the context of causality. \citet{ernest2015marginal}, show that it can be used to formulate the conditional post-intervention mean. 


    




The paper is structured as follows. First, we introduce the nonparametric setting and notation of our work in Section \ref{section:population}. The Shapley curves are defined on population level and examples are given for building the intuition of the reader. In Section \ref{section:estimation} we propose two estimation approaches for Shapley curves. Our main theorems regarding asymptotics are given and we elaborate on their heuristics. Next, we validate our asymptotic results in a simulation exercise in Section \ref{section:numerical}. In Section \ref{section:application} we estimate Shapley curves on a vehicle price data set and analyze the accompanied confidence intervals. Finally we conclude in Section \ref{section:conclusion}. All codes of this paper are available on \href{http://www.quantlet.de}{quantlet.de}. \hspace*{\fill} \raisebox{-1pt}{\includegraphics[scale=0.02]{qletlogo.png}}\,










\section{A Smoothing Perspective on Shapley Values}
\label{section:population}
\subsection{Model Setup and Notation}
Consider a vector of covariates $X=(X_1,\ldots,X_d)^{\top}\in\mathcal{X}\subset\mathbb{R}^{d}$ and a scalar response variable $Y\in\mathcal{Y}\subset\mathbb{R}$. Let $F$ denote the cumulative distribution function (cdf) of $X$ with continuous density (pdf), $f$. Let $\{(X_i,Y_i)\}_{i=1}^{n}$ be a sample drawn from a joint distribution function $F_{X,Y}$. Consider the following nonparametric regression setting,
\begin{align}\label{eq:dgp}
    Y_i=m(X_i)+\varepsilon_i,\quad i=1,\ldots,n,
\end{align}
with ${\sf E}(\varepsilon_i|X_i)=0$ and $m\in\mathcal{M}$, where $\mathcal{M}$ is a rich class of functions. Consequently, $m(x)={\sf E}(Y|X=x)$ is the conditional expectation.
Let $N\stackrel{\text{def}}{=}\{1,\ldots,d\}$, and let $\mathcal{S}$ denote the power set of $N$. For a set $s\in\mathcal{S}$, $X_s$ denotes a vector consisting of elements of $X$ with indices in $s$. Correspondingly, $X_{-s}$ denotes a vector consisting of elements of $X$ not in $s$. We write $m_s\in\mathcal{M}_s$ to denote functions which ignore arguments with index not in $s$, $\mathcal{M}_s=\{m\in\mathcal{M}:m(u)=m(v)\text{ for all }u,v\in\mathcal{X}\text{ satisfying }u_s=v_s\}$. Similarly, we write $m_{-s}\in\mathcal{M}_{-s}$ for functions which ignore arguments in $s$. Finally, we write $f_{X_{-s}|X_{s}}(x_{-s}|x_s)$ for the conditional density functions of $X_{-s}$ given $X_s=x_s$.
The indicator function $\mathbf{I}\{j \in s\}$ is equal to $1$ if $j\in s$ and $0$ otherwise. The sign function $\text{sgn}\{j \in s\}$ takes the value $1$ if $j\in s$ and $-1$ otherwise.





\subsection{Population Shapley Curves}
We now define the population-level Shapley curves, which are functions, $\phi_j(\cdot):\mathbb{R}^{d}\to\mathbb{R}$ measuring the local variable importance of a variable $j$ at a given point $x\in\mathbb{R}^d$. Opposed to previous literature we have a variable importance measure defined on the whole support instead of the observations $x_i$ for $i=1,\ldots,n$. The functions are defined on the population level, meaning that they are uniquely determined by $m(x)$ defined in (\ref{eq:dgp}) and the joint distribution function of the covariates, $F$. This is opposed to the commonly-used Shapley value of \citet{lundberg2017unified}, which measures the importance of a variable $X_j$ in a predictive model. As a consequence, looking at variable importance on the population level is agnostic of the specific form a corresponding estimator might take.

Originally proposed in a game theory framework, Shapley values are measuring the difference of the resulting payoff for a coalition of players and the same coalition including an additional player \citep{shapley1953s}. By keeping a player fixed, this difference is calculated across all possible coalitions of players and averaged with a combinatorial weight.
From a statistical literature point of view, each player is represented by a variable and the payoff is therefore a measure of contribution for the corresponding subset of variables. We define this measure to be the conditional mean. A recent application of Shapley values is the estimation of the probability of default in credit scoring \citep{li2022data}.
Finally, let us define $\phi_j(x):\mathbb{R}^d\to\mathbb{R}$, as follows,
\begin{align}\label{eq:shapley}
    \phi_j(x)=\sum_{s\subset N\setminus{j}}\frac{1}{d}\binom{d-1}{|s|}^{-1}\left\{m_{s\bigcup{j}}(x_{s\cup{j}})-m_{s}(x_s)\right\},
\end{align}
for $j\in N$, where the components are defined as
\begin{align}\label{eq:component}
\begin{split}
    m_s(x_s)&=\int m(x)f_{X_{-s}|X_s}\left(x_{-s}|x_s\right)dx_{-s}\\
    &={\sf E}\left(Y|X_s=x_s\right),
\end{split}
\end{align}
for $s\in\mathcal{S}$. Note that $m_N(x)=m(x)$ represents the full nonparametric model and $m_{\emptyset}={\sf E}(Y)$ is the unconditional mean of the response variable. A point $\phi_j(x)$ in the Shapley curve measures the difference in the conditional mean of $Y$ from including the variable $X_j$, averaged through the combinatorial weight over all possible subsets. 







Besides the simplicity of the Shapley decomposition, another major reason for its popularity is the satisfaction of several convenient properties. Most importantly, our definition of Shapley curves (\ref{eq:shapley}), satisfies the crucial \textit{additivity} property $m(x)-{\sf E}(Y)=\sum_{j=1}^{d}\phi_j(x)$.
Namely, the conditional expectation of the full nonparametric model (\ref{eq:dgp}) subtracted by the unconditional mean of the dependent variable, $m(x)-{\sf E}(Y)$, can be recovered exactly as a sum of Shapley curves of variables $X_j$ evaluated at the point $x$.






    
    

    
    
    
    







In contrast to the Shapley-based variable importance measures of \citet{williamson2020efficient}, \citet{covert2020understanding} and \citet{benard2021shaff}, Shapley curves provide a local instead of a global assessment of the importance of a variable. Consequently, the Shapley curve will take different values for different points on the $d$-dimensional support $\mathcal{X}$ of the covariates. It might be the case that a variable is redundant in certain areas of $\mathcal{X}$, but indispensable in other areas. Local measures are thus able to paint a more nuanced picture. 

We want to highlight the important role that the dependency in $X$ plays in (\ref{eq:component}), and consequently in the definition of the Shapley curves. For a set $s\in\mathcal{S}$, $m_s(x_s)$ represents the conditional expected value of the dependent variable, with variables not in $s$ integrated out with respect to the conditional density $f_{X_{-s}|X_s}(x_{-s}|x_s)$. In the special case of independent covariates, the expression for the conditional density simplifies to a product of (unconditional) marginal densities, $f_{X_{-s}|X_s}(x_{-s}|x_s)=\prod_{j\in N\setminus{s}}f_{X_j}(x_j)$. The density does not depend on the location of the conditioning variables $X_s=x_s$, the integral can thus be interpreted as taking the unconditional expectation with respect to the covariates $X_j$, $j\in N\setminus s$. In general however, the component (\ref{eq:component}), and consequently $\phi_j(x)$, will depend crucially on the dependency structure of $X$. In particular, the difference in the conditional expected value for a given set $s$, $m_{s\bigcup{j}}(x_{s\cup{j}})-m_{s}(x_s)$, can be a result of the direct effect of variable $X_j$ on $Y$ via the functional relationship described by $m(\cdot)$, or it might be due to the dependence of $X_j$ with another variable $X_k$, which in turn has a direct effect on $Y$. This can be seen as a bias caused by the endogeneity of $X_j$. It is therefore crucial to understand, that Shapley curves, even if they are defined on the population, are a predictive measure of variable importance and not a causal measure. We will demonstrate the role of dependence in a few examples in the next subsection.





\subsection{Examples}

Let us have a look at Shapley curves in a few interesting scenarios. Consider different settings for the regression function, $m(x)$, as well as for the dependence structure of $X$. In particular, we are interested in the difference between the case of independent and dependent regressors.

\subsubsection*{Example 1: Additive Model} \label{ex1}

As a first scenario, we consider the case of additive models \citep{buja1989linear}, where the covariates are assumed to be independent. The conditional expectation of the response variable is a sum of one-dimensional nonparametric functions,
\begin{align}
\label{add_ex}
    m(x)=C+\sum_{j=1}^{d}g_j(x_j),
\end{align}
where $C$ is a constant and $g_j(x_j)$ is the partial dependence function for variable $j$. For the purpose of identification we assume that ${\sf E}\{g_j(X_j)\}=0$ for all $j$, and for simplicity we assume $C=0$. Additive models include no interaction effects between different regressors. Together with the independence of the covariates, one obtains the Shapley curve,
\begin{align*}
    \phi_j(x)=g_j(x_j),
\end{align*}
for all $j\in N$. The Shapley curve of variable $X_j$ is thus identical to the $j$-th nonparametric function. In turn, it is also one-dimensional, since it only depends on $X_j$.
Let us consider the following example from \citet{khowaja2022uniform}, 
\begin{align*}
    m(x) = \operatorname{max}(1 - \frac{\lvert x_1 \rvert}{0.2}) + \operatorname{max}(1 - \frac{\lvert x_2 \rvert}{0.6}),
\end{align*}
where $C=0$ and ${\sf E}\{g_j(X_j)\}\neq0$. The Shapley curve results as

\begin{align*}
    \phi_j(x) = \operatorname{max}(1 - \frac{\lvert x_j \rvert}{\eta_j}),\quad \eta_1=0.2,\ \eta_2=0.6
\end{align*}
Figure \ref{iphone} illustrates the corresponding population Shapley curves for both variables.



\begin{figure}[H]
\vspace{-0.3cm}
\hspace{-0.8cm}
\begin{minipage}{0.5\textwidth}
		\includegraphics[scale=0.20]{iphone1.png}
\end{minipage}
\hspace{0.1cm}
\begin{minipage}{0.5\textwidth}
		\includegraphics[scale=0.20]{iphone2.png}
\end{minipage}\\
\vspace{-1.5cm}
\caption{Population Shapley Curves}
\label{iphone}
  \centering
  \vspace{0.1cm}
  \begin{minipage}{15cm}%
\footnotesize
Population Shapley curve for $j=1$ (left) and $j=2$ (right). \href{https://github.com/QuantLet/Shapley_Smoothing/tree/main/surface_plots}{\hspace*{\fill} \raisebox{-1pt}{\includegraphics[scale=0.02]{qletlogo.png}}\,{Ex\_add}}
\end{minipage}
\end{figure}





\subsubsection*{Example 2: Linear Model under Dependence}

Now consider the linear regression model with $d=2$ and assume that the covariates are dependent,
\begin{align*}
    m(x)=x_1\beta_1+x_2\beta_2.
\end{align*}
It follows from (\ref{eq:component}) that the Shapley curve $\phi_1(x)$ is,
\begin{align}
  
   \phi_1(x) =x_1\beta_1+\frac{1}{2}\left\{{\sf E}\left(X_2|X_1=x_1\right)\beta_2-{\sf E}\left(X_1|X_2=x_2\right)\beta_1-{\sf E}(Y)\right\}.
\end{align}
By symmetry, $\phi_2(x)$ is defined in a similar way. One can see that the Shapley curve is decomposed into a direct effect and an indirect effects. The latter is a consequence of variable dependence and is captured by including a conditional density in (\ref{eq:component}) instead of a marginal density. As a practical example of this phenomenon, \citet{berman2012simpson} consider the relation between price and demand of a good. While the direct effect is negative as predicted by economic theory, the time variable is correlated positively with both price and demand and therefore introduces a positive indirect effect, which dominates the overall effect. This is known as Simpson's paradox \citep{aldrich1995correlations}.
Even in the extreme case of conditional mean independence, $\beta_1=0$, the indirect effects lead to $\phi_1(\cdot)\neq0$. For pure prediction purposes, this is an advantage since all effects are captured. For the goal of interpretation, this byproduct is a burden since we can not distinguish between direct and indirect effects on the sample level.





\subsubsection*{Example 3: Threshold Regression Model}

While the previous two examples consider regression models without interactions, we now focus on a simple parametric model in which the effect of the variables does not enter linearly. We look at the threshold regression model \citep{dagenais1969threshold}, with regressor variable $X_1$ and threshold variable $X_2$. The conditional expectation of $Y$ given $x$ is defined by
\begin{align*}
    m(x)=\left\{\psi+\theta\mathbf{I}(x_{2}\leq C)\right\}x_{1},
\end{align*}
where $C,\psi,\theta$ are scalar parameters. An empirical example for the threshold regression model is given in \citet{hansen2000sample} for the multiple equilibria growth model of \citet{durlauf1995multiple}. The analysis suggests the existence of a regime change in GDP growth depending on whether the initial endowment of a country exceeds a certain threshold.
Then we have the following expression for the Shapley curves given the two variables,
\begin{align*}
    \phi_1(x)=\left\{\psi+\frac{1}{2}\theta\mathbf{I}(x_{2}\leq C)+\frac{1}{2}\theta F_{X_2}(C)\right\}\left\{x_{1}-{\sf E}(X_{1})\right\},
\end{align*}
and
\begin{align*}
    \phi_2(x)&=\frac{1}{2}\theta\left\{\textbf{I}(x_{2}\leq C)-F_{X_2}(C)\right\}\left\{x_1+{\sf E}(X_1)\right\},
\end{align*}
where $F_{X_2}$ is the marginal cdf of $X_2$. The Shapley curve for variable $X_1$ will take a larger value in absolute values whenever the coordinate $x_1$ is far away from the unconditional expectation, ${\sf E}(X_1)$. Similarly, the Shapley curve for the threshold variable $X_2$ will take a large value if the difference $\left\{\textbf{I}(x_{2}\leq C)-F_{X_2}(C)\right\}$ is large, i.e., in situations in which the effect of the variable is not well captured by its unconditional expectation.



\section{Estimators and Asymptotics} 
\label{section:estimation}

In this section, we discuss two approaches for estimating Shapley curves. By considering \textit{curves} instead of a \textit{value} of importance, we are able to gain insights on the whole support of the covariates, instead of being restricted to the observed values. The goal is to establish conditions for the consistent estimation under general conditions. This is a crucial task because else we are not able to determine what is estimated on sample level, i.e., whether the estimate is meaningful in any way. For instance, \citet{scornet2020trees} and \citet{benard2022mean} demonstrate that several measures for variable importance in random forests only have a meaningful population-level target in the restrictive case of independence of regressors and no interactions. Both approaches discussed in this chapter are plug-in estimators using the Shapley formula in (\ref{eq:shapley}). The first approach is based on separate estimators of all component functions, for all subsets $s$. We call this the \emph{component-based} approach. Alternatively, one can obtain estimators of the component $m_s(x_s)$ using a pilot estimator for the full model and integrate out the variables not contained in the set $s$ with respect to the estimated conditional densities using (\ref{eq:component}). This method we denote as the \emph{integration-based} approach.

\subsection{Component-Based Approach}
\label{comp_chapter}

    
The component-based approach for estimating Shapley curves involves the separate estimation of each component. For each subset $s\in\mathcal{S}$, we use a nonparametric smoothing technique to regress $Y_i$ on $X_{s,i}$, i.e., those regressors contained in $s$. This yields in total $2^d$ regression equations,
\begin{align}
    Y_i=m_{s}(X_{s,i})+\varepsilon_{s,i},\quad i=1,\ldots,n;\quad s\in\mathcal{S},
\end{align}
with ${\sf E}(\varepsilon_{s,i}|X_{s,i})=0$ and $m_{s}$ is defined as in equation (\ref{eq:component}). Given an estimator for the components, we obtain a plug-in estimator for the Shapley curve of variable $j$ using (\ref{eq:shapley}). As a result, the estimated Shapley curve follows as

\begin{align}\label{eq:shapley_est}
    \widehat{\phi}_j(x)=\sum_{s\subset N\setminus{j}}\frac{1}{d}\binom{d-1}{|s|}^{-1}\left\{\widehat{m}_{s\bigcup{j}}(x_{s\cup{j}})-\widehat{m}_{s}(x_s)\right\}.
\end{align}
Since $m_s$ belongs to a nonparametric class of functions, we consider local linear estimators \citep{fan1993local} for the component functions. Let $Y=(Y_1,\ldots,Y_n)^\top$, $Z_s=(z_{s1},\ldots,z_{sn})^\top$, $z_{si}=(1,X_{s_1,i},\ldots,X_{s_d,i})$ and $K_s=\operatorname{diag}(\{h_s^{-d_s}\prod_{j=1}^{d_s}k(h_{s}^{-1}(X_{s_j,i}-x_j))\}_{i=1}^n)$, where $k$ is a one-dimensional kernel function and $h_s$ is the bandwidth. Then we have,
\begin{align*}
    \widehat{\beta}_s(x_s)=\begin{pmatrix}\widehat{\beta}_{s,0}(x_s)\\\widehat{\beta}_{s,1}(x_s)\end{pmatrix}=\left(Z_s^\top K_sZ_s\right)^{-1}Z_s^\top K_sY,
\end{align*}
and the local linear estimator is $\widehat{m}_s(x_s)=\widehat{\beta}_{s,0}(x_s)$.

We now look at the asymptotic properties of the component-based estimator. In particular, we establish uniform consistency and asymptotic normality. This has at least two important implications. First, it shows that the estimation of Shapley curves on the sample level is meaningful, i.e., that the estimator will estimate the correct quantity. And second, this allows us to construct confidence intervals and to conduct hypothesis tests on the estimated curves. For this purpose, we impose the following regularity assumptions.
\begin{assumption}\label{assum:density}\
\begin{enumerate}[(i)]
    \item The support of $X$ is $\mathcal{X}$.
    \item The density $f$ of $X$ is bounded, bounded away from zero and twice continuously differentiable on $\mathcal{X}$.
    \item ${\sf Var}(\varepsilon|x)=\sigma^2(x)<\infty$ for all $x\in\mathcal{X}$.
    \item ${\sf E}(|Y|^{2 + \delta} \lvert X=x)< \infty$ for some $\delta > 0$.
\end{enumerate}
\end{assumption}

\begin{assumption}\label{assum:function}
Assume $m(x)$ belongs to $\mathcal{M}_d$, the space of $d$-dimensional twice continuously differentiable functions.
\end{assumption}

\begin{assumption}\label{assum:kernel}
Assume $k(\cdot)$ is a univariate twice continuously differentiable probability density function symmetric about zero and $\int s^2k(s)ds=\mu_2(k)<\infty$ and $\int k^{2+\delta}(s)ds<\infty$ for some $\delta>0$.
\end{assumption}

The following theorem shows the uniform consistency of the component-based estimator in the mean integrated squared error (MISE) sense.
\begin{proposition}
\label{proposition:conv_comp}
Let $\widehat{\phi}_j(x)$ be the component-based estimator with components estimated via the local linear method with bandwidths $h_s\sim n^{-\frac{1}{4+|s|}}$. Then we have under Assumptions \ref{assum:density}, \ref{assum:function} and \ref{assum:kernel}, as $n$ goes to infinity,
\begin{align*}
    \operatorname{MISE}\left\{\widehat{\phi}_j(x),\phi_j(x)\right\} = \mathcal{O}\left(n^{-\frac{4}{4+d}}\right).
\end{align*}
\end{proposition}

The proof of Proposition \ref{proposition:conv_comp} can be found in \ref{glob_proof}. The main idea is to write the difference between the estimator and the true Shapley curve in terms of a weighted sum,
\begin{align}
\label{eq:weightedsum}
    \widehat{\phi}_j(x)-\phi_j(x) &= \sum_{s \subseteq N} \omega_{j,s} \left\{\widehat{m}_s(x_s)-m_s(x_s)\right\},
\end{align}
where the weights are defined as
\begin{align}
    \omega_{j,s}&=\text{sgn}\{j \in s\}\frac{1}{d}\binom{d-1}{|s| - \mathbf{I}\{j \in s\}  }^{-1}
\end{align}
for all $s\in\mathcal{S}$.
We show that the leading term in the MISE of the component-based estimator depends on the MISE of the local linear estimator for the full model, and since the corresponding weight is defined as $\omega_{j,m}=\frac{1}{d}$,
\begin{align*}
    \operatorname{MISE}\left\{\widehat{\phi}_j(x),\phi_j(x)\right\} = \frac{1}{d^2}\operatorname{MISE}\left\{\widehat{m}(x),m(x)\right\}+{\scriptstyle\mathcal{O}}(n^{-\frac{4}{4+d}}).
\end{align*}
Ultimately, the convergence rate is determined by the slowest rate of the components, i.e. the convergence rate of the full model, which is known to be $\mathcal{O}(n^{-\frac{4}{4+d}})$.
Since $\phi_j(x)$ also belongs to $\mathcal{M}_d$, the class of twice continuously differentiable functions, $\widehat{\phi}_j(x)$ is a minimax-optimal estimator for $\phi_j(x)$ by \citet{stone1982optimal}.

The following theorem establishes the point-wise asymptotic normality of the component-based estimator and the corresponding proof is given in \ref{proof_comp_asy}. Denote $h_m$ as the optimal bandwidth of the local linear estimator for the full model.

\begin{theorem}
\label{theorem:comp_normal}
Let the conditions of Proposition \ref{proposition:conv_comp} hold. Then we have, for a point $x$ in the interior of $\mathcal{X}$, as $n$ goes to infinity,
\begin{align*}
    \sqrt{nh^d_m}\left\{\widehat{\phi}_j(x)-\phi_j(x)\right\}=\sqrt{nh^d_m}\frac{1}{d}\left\{\widehat{m}(x)-m(x)\right\}+{\scriptstyle\mathcal{O}}_p(1)\stackrel{\mathcal{L}}{\rightarrow}N\left(B(x),V(x)\right),
\end{align*}
where the asymptotic bias is given as
\begin{align*}
    B(x)=\frac{1}{d} \frac{\mu_2(k)}{2} \sum^d_{j=1} \frac{\partial^2 m(x)}{\partial^2 x_j}
\end{align*}
and the asymptotic variance is given as
\begin{align*}
    V(x)=\frac{1}{d^2} ||k||^2_2 \frac{\sigma^2(x)}{f(x)},
\end{align*}
where $||k||^2_2=\int k^2(s)ds$ denotes the squared $L_2$ norm of $k$.
\end{theorem}
        
Interestingly, the asymptotic distribution of the component-based estimator is the same for all variables, $j=1,\ldots,d$, as neither bias nor variance are variable-specific. Building upon the asymptotic normality result, we are able to construct asymptotically valid confidence intervals, which allow us to conduct inference on the estimated Shapley curves.
        
\begin{corollary}\label{col:conf}
As a consequence of Theorem \ref{theorem:comp_normal}, we can construct the $100(1-\alpha)\%$ confidence interval for $\phi_j(x)$ as 
\begin{align*}
\widehat{\phi}_j(x) - n^{-\frac{2}{4+d}} \widehat{B}(x)  \pm n^{-\frac{2}{4+d}}  z_{1-\frac{\alpha}{2}}\sqrt{\widehat{V}(x)} ,
\end{align*}
where we use the optimal bandwidth $h_m \sim n^{-\frac{1}{4+d}}$ and $z_{1-\frac{\alpha}{2}}$ is the $(1-\frac{\alpha}{2})$-quantile of the normal distribution. The estimators for the bias and the variance, $\widehat{B}(x)$ and $\widehat{V}(x)$ respectively, can be obtained via plug-in estimates as described in \citet{hardle2004nonparametric}.
\end{corollary}








The component-based estimator for the Shapley curves, $\widehat{\phi}_j(x)$, requires the estimation of $2^d$ conditional mean functions, $m_s(x_s)$, for all $s\in\mathcal{S}$. This task is particularly cumbersome when $d$ is large and the chosen estimator is computationally intensive. In such a situation we are in need of an approximation for Shapley curves. Fortunately, \citet{williamson2020efficient} provide an approximation technique which is based on subsampling the subsets $s\in\mathcal{S}$. This subsampling idea is incorporated into a constrained weighted least-squares problem. The authors derive a closed-form solution, which drastically reduces the computational cost. In principle, we can directly incorporate both of our estimation approaches into their closed-form solution and obtain approximated Shapley curves. 
Further, we know that the approximated curves converge to the population counterpart as the number of subsamples goes to infinity. Although the inclusion of this approximation is straightforward for our asymptotic results, it is of no relevance as long as the subsample size approaches infinity.  


The underlying function class in Proposition \ref{proposition:conv_comp} and Theorem \ref{theorem:comp_normal} is quite large, leading to a slow convergence rate in a setting with high-dimensional covariates. In fact, the rate is decreasing exponentially if the number of variables, $d$, increases. This is known as the curse of dimensionality. However, our consistency result can be easily extended to other function classes. The key is to show the consistency and convergence rates for the estimators of the component functions in the respective function class. For instance, \citet{chen1999improved} consider a different smoothness class and show that shallow neural networks achieve a rate of convergence of ${\scriptstyle\mathcal{O}}_p(n^{-1/4})$, irrespective of $d$. Using their neural network estimator to estimate the component functions will lead to the same convergence rate for the Shapley curves. Recently, \citet{kohler2021rate} and \citet{schmidt2020nonparametric} considered function classes, in which the functions are compositions of lower-dimensional functions. It is shown that deep neural networks are able to learn those functions at the rate of the lower-dimensional element. Consequently, we can use deep neural networks to estimate the component functions which leads to the same convergence rate for the estimation of the Shapley curve.





\subsection{Integration-Based Approach}
\label{integ_chapter}

In contrast to the component-based method, the integration-based approach requires the estimation of only one regression function, namely that of the full model. This pilot estimator is obtained by local linear estimation and is thus identical to the estimator for the full model in the component-based approach. To estimate the regression functions for the subsets $s\in\mathcal{S}$, the variables not contained in $s$ are integrated out using an estimate for the conditional density,
\begin{align}
    \widetilde{m}_s(x_s)=\int_{\mathbb{R}^{d-|s|}} \widehat{m}(x)\widehat{f}_{X_{-s}|X_s}(x_{-s}|x_s)dx_{-s}.
    \label{eq:dependent}
\end{align}
In the case of mean-independent regressors, there is no need to estimate conditional densities. In this case, a simplified estimator can be used, which averages over the observations,
\begin{align}
    \widetilde{m}_s(x_s)=\frac{1}{n}\sum_{i=1}^{n}\widehat{m}(X_s=x_s,X_{-s,i}).
    \label{eq:independent}
\end{align}
This estimator is well-known in the literature as the marginal integration estimator for additive models \citep{linton1995kernel,tjostheim1994nonparametric} and is discussed in a more general setting by \citet{fan1998direct}. However, the assumption of mean-independent regressors is not plausible in general, and the estimator described in (\ref{eq:independent}) will lead to inconsistent estimates for the true component, $m_s(x_s)$. In this setting, the estimation of the conditional densities is a necessity and one needs to use the estimator (\ref{eq:dependent}).

To study the asymptotic behavior of the integration-based estimator for Shapley curves, we assume the simplified setting of known conditional densities. Thus, we consider the following estimator for the component function of subset $s$,
\begin{align}
\label{eq:integ}
    \widetilde{m}_s(x_s)=\int\widehat{m}(x)f_{X_{-s}|X_s}(x_{-s}|x_s)dx_{-s}.
\end{align}
In analogy to the component-based estimated Shapley curve (\ref{eq:shapley_est}), we obtain the integration-based estimated Shapley curve  $\widetilde{\phi}_j(x)$ as a plug-in estimate using $\widetilde{m}_s(x_s)$.


The following theorem shows the uniform convergence rate and the asymptotic distribution of the integration-based estimator.
\begin{theorem}
\label{theorem:conv_integ}
Under Assumptions \ref{assum:density}, \ref{assum:function} and \ref{assum:kernel}, let $\widetilde{\phi}_j(x)$ be the integration-based estimator with known density and a pilot estimator based on local linear estimation with bandwidth $h_m\sim n^{-\frac{1}{4+d}}$. Then we have as $n$ goes to infinity,
\begin{align}
    \operatorname{MISE}\left\{\widetilde{\phi}_j(x),\phi_j(x)\right\} = \mathcal{O}\left(n^{-\frac{4}{4+d}}\right). \label{eq:conv_integ}
\end{align}
Further we have, for a point $x$ in the interior of $\mathcal{X}$, as $n$ goes to infinity the asymptotic distribution of the integration-based estimator,
\begin{align*}
    \sqrt{nh^d_m}\left\{\widetilde{\phi}_j(x)-\phi_j(x)\right\} =  \sqrt{nh^d_m}\sum_{s \subseteq N} w_{j,s}  \left\{  \widetilde{m}_s(x_s) - m_s(x_s) \right\}\stackrel{\mathcal{L}}{\rightarrow}N\left(B_{int}(x),V(x)\right),
\end{align*}
where the asymptotic bias term is 
\begin{align}
B_{int}(x)= \frac{\mu_2(k)}{2} \sum_{s \subseteq N} w_{j,s} \left\{ \sum^d_{j=1} \int_{X_{-s}} \frac{\partial^2 m(x)}{\partial^2 x_j} f_{X_{-s}|X_s}(x_{-s}|x_s)dx_{-s}\right\}
\label{eq:bias}
\end{align}
and the asymptotic variance term is
\begin{align*}
    V(x)=\frac{1}{d^2} ||k||^2_2 \frac{\sigma^2(x)}{f(x)}.
\end{align*}
\end{theorem}

The proof of Theorem \ref{theorem:conv_integ} is given in \ref{proof_integ_asy}. Notice that the rate is identical to the convergence rate of the component-based approach and it is also optimal in the minimax sense. The reason is that, once again, the (pilot) estimator $\widehat{m}(x)$ determines the convergence rate of $\widetilde{\phi}_j(x)$.

The asymptotic variance is identical to that of the component-based approach. The difference in the asymptotic distribution is solely due to the bias. Compared with the component-based approach, the bias term defined in (\ref{eq:bias}) is now a sum over $2^d$ elements, instead of a single term. As a consequence, the bias will be larger.


This fact becomes clear when looking at the bias and the variance of the estimated components. Following \citet{linton1995kernel}, we get the following expression for the bias and the variance of $\widetilde{m}_s(x_s)$,
\begin{align*}
    \operatorname{Bias}\left\{\widetilde{m}_s(x_s)\right\}&=\frac{1}{2}h_m^2\mu_2(k)\int\sum_{j=1}^{d}\frac{\partial^2 m(x_s)}{\partial x_{j}^2}f_{X_{-s}|X_s}(x_{-s}|x_s)dx_{-s}+{\scriptstyle\mathcal{O}}(h_m^2),\\
    \operatorname{Var}\left\{\widetilde{m}_s(x_s)\right\}&=\frac{1}{nh_m^{|s|}}||k||^2_2\int \frac{\sigma^2(x)f^2_{X_{-s}|X_s}(x_{-s}|x_s)}{f(x)}dx_{-s}+{\scriptstyle\mathcal{O}}\left(\frac{1}{nh_m^{|s|}}\right).
\end{align*}
The bandwidth of the pilot estimator is chosen as $h_m\sim n^{-\frac{1}{4+d}}$, which balances the squared bias and variance only for the component associated with the full model. However, all components are based on this pilot estimator and thus rely on the same bandwidth. For all other subsets $s$, we have $|s|<d$ which leads to oversmoothing,
\begin{align*}
    \operatorname{Bias}^2\left\{\widetilde{m}_s(x_s)\right\}&=\mathcal{O}(n^{-\frac{4}{4+d}})\\
    \operatorname{Var}\left\{\widetilde{m}_s(x_s)\right\}&=\mathcal{O}(n^{-\frac{(4+d-|s|)}{4+d}}).
\end{align*}
As a consequence, the bias of the lower-dimensional components does not vanish in the bias of the integration-based estimator, as it is of the same order for all the components, namely $\mathcal{O}(n^{\frac{2}{4+d}})$. While the order of the bias term in the integration-based estimator is still identical to that of the component-based approach, it might be substantially larger in finite samples.
However, the advantage of the integration-based curve is that the pilot estimator only needs to be estimated once. Essentially, this constitutes a trade-off between computational complexity and accuracy in the estimation.

\begin{remark}
While the focus in this Section lies on the local linear estimator and the class of twice-continuously differentiable functions, the bias problem of the integration-based approach most likely also arises in other contexts. This is due to the pilot estimator, which typically involves the selection of hyperparameters governing the bias-variance trade-off. For random forests, this could be the depth of the individual trees, for neural networks, the number of layers, and number of nodes. There is no reason to believe that the optimal choice of hyperparameters for the pilot estimator, which is based on the set of all regressors, is optimal for the estimation of the components associated with the lower-dimensional subsets. These components are less complex than the full model, which the integration-based approach cannot accommodate, while the component-based approach can.


\end{remark}
















The previous theorems are affected by the curse of dimensionality. Under the assumption of additivity of the true process and independence of the explanatory variables, we know that the Shapley curve simplifies to the partial dependence function, as introduced in Example 1. 
\begin{assumption}\label{assum:additive}
Assume the regression function $m(x)$ follows an additive structure, s.t. $m(x)=\sum^d_{j=1} g_j(x_j)$ with ${\sf E}\{g_j(x_j)\}=0$ for $j=1,\ldots,d$ and the explanatory variables are independent.
\end{assumption}
It follows directly by \citet{stone1985additive} that the partial dependence function, and thus the corresponding Shapley curve, can be estimated with a one-dimensional rate, see Corollary \ref{col:additive}. The estimation of the partial dependence function is sufficient to obtain an estimator for $\phi_j(x)$. This can be done consistently via marginal integration \citep{linton1995kernel} or backfitting \citep{mammen1999existence}. Further insights into the interpretation of smooth backfitting and marginal integration for additive models is given in \citet{nielsen1998optimization}.
The bias and the variance of the estimated Shapley curves follow by the established asymptotic results for the estimator of the partial dependence function for variable $j$. 


\begin{corollary}\label{col:additive}
Let the partial dependence function $\widehat{g}_j(x_j)$ be an estimator for $\phi_j(x)$ obtained by marginal integration or backfitting. Under assumptions \ref{assum:density}, \ref{assum:function}, \ref{assum:kernel} and \ref{assum:additive} we have that
\begin{align*}
    \operatorname{MISE}\left\{\widehat{g}_j(x),\phi_j(x)\right\} = \mathcal{O}\left(n^{-\frac{4}{5}}\right).
\end{align*}
\end{corollary}


\begin{remark}
The previous results in this Section rely on the knowledge of the true conditional densities for the integration-based estimation of the components. In practice, these densities have to be estimated, which introduces another estimation error of order $\mathcal{O}_p(n^{-\frac{2}{4+d}})$. This is a further disadvantage of the integration-based method of estimating Shapley curves compared with the component-based approach. 
\end{remark}







\section{Numerical Studies}
\label{section:numerical}
In this section, we conduct simulation studies to validate the previous asymptotic results. In particular, we empirically show the consistent estimation of both estimation approaches, the component-based as well as the integration-based. Furthermore, empirical coverage is shown. 
Let the regressors be Gaussian, $X_{1:3}\sim N(0, \Sigma)$, where $\Sigma$ is the covariance matrix with variance $\sigma^2=4$. The correlation $\rho$ is either $0$ or $0.8$ for different setups. The corresponding density functions are assumed to be known.
Let the true relationship between the dependent and independent variables be defined as
\begin{align}
y &= m(x) + \varepsilon\ \text{, where}\ \varepsilon \sim N(0,1).
\end{align}
The first data-generating process (DGP) follows an additive model and the second one includes interactions:
\begin{align*}
&\textbf{\text{DGP 1: Additive}}\ \ \ \ \ \ m(x) = -\text{sin}(2x_1) + \text{cos}(2x_2) + x_3 \\
&\textbf{\text{DGP 2: Non-linear}}\ \ \ m(x) = -\text{sin}(2x_1) + \text{cos}(3x_2) + 0.5x_3 + 2\text{cos}(x_1)\text{sin}(2x_2). 
\end{align*}
To investigate for robustness, we include a variation for $\varepsilon \sim t(5)$.
The local linear estimator is used to obtain component-based and integration-based Shapley curves for all three variables. The bandwidths $h_1,\ldots,h_d$ are allowed to differ in each direction and are chosen via leave-one-out cross validation in each $j$-direction. The second-order Gaussian kernel is used. 
To evaluate the global performance of both estimators, we calculate the MISE for each simulation setup:
\begin{align}
\label{mise_ref}
     \operatorname{MISE}\left\{\widehat{\phi}_j(x),\phi_j(x)\right\} = {\sf E} \left[  \int_{\mathbb{R}^d} \left\{ \widehat{\phi}_j(x) - \phi_j(x)  \right\}^2 dx \right].
\end{align}
The integral is approximated with the \textit{cubature} R package, which adaptively calculates the multidimensional integral over hypercubes. Our simulation results are based on $6000$ Monte Carlo iterations for increasing sample sizes, $n=300, 500, 1000, 2000$.

The population Shapley curves as well as the estimated curves using the component-based estimation approach are illustrated in Figure \ref{comp_plot} for the non-linear DGP 2. Since $x_1$ and $x_2$ do not have interactions with the third variable on the response, we fix $x_3=0$ to reduce the four-dimensional illustration to three dimensions. Ultimately, $\phi_j(x_1,x_2,x_3=0)$ is plotted together with the estimated counterparts for $j=1,2$. Since we consider Shapley \textit{curves}, it is possible to evaluate and illustrate the function $\phi_j(x_1,x_2,x_3)$ (and the estimated curves) at any arbitrary fixed vector $x$.

\begin{figure}[H]
\hspace{-0.6cm}
\begin{minipage}{0.5\textwidth}
		\includegraphics[scale=0.23]{comp1.png}
\end{minipage}
\hspace{-0.3cm}
\begin{minipage}{0.5\textwidth}
		\includegraphics[scale=0.23]{comp1_SE.png}
\end{minipage}\\
\vspace{-1cm}
\caption{Component-based Shapley Curves}
\label{comp_plot}
  \centering
  \vspace{0.1cm}
  \begin{minipage}{13cm}%
\footnotesize
Left: Estimated Shapley curve (in red) and population Shapley curve (in blue) for $j=1$. Right: Squared difference between estimated and population Shapley curve. Local Linear estimator for DGP 2 with $n=1000$ observations. \href{https://github.com/QuantLet/Shapley_Smoothing/tree/main/surface_plots}{\hspace*{\fill} \raisebox{-1pt}{\includegraphics[scale=0.02]{qletlogo.png}}\,{CB\_plot}}
\end{minipage}
\end{figure}



Recall that $\widetilde{\phi}_j(x)$ and $\widehat{\phi}_j(x)$ are denoted as the estimated integration-based Shapley curve and the estimated component-based Shapley curve for variable $j$.
The simulation results are displayed in Table \ref{tab1} for both DGPs and Gaussian error terms. First, we observe that the MISE is shrinking for both estimators as the sample size increases. This result empirically confirms Proposition \ref{proposition:conv_comp} and the first part of Theorem \ref{theorem:conv_integ}.
Second, the component-based approach (almost) always results in a smaller MISE than the integration-based approach. This aligns with the asymptotic bias and asymptotic variance comparison of both estimators in Chapter \ref{section:estimation}. Since we oversmooth within the integration-based estimator, the bias of $\widetilde{\phi}_j$ accumulates over each component. In contrary, the bias of $\widehat{\phi}_j$ is only determined by a single component, namely the one associated with the full $d$-dimensional model. 
Since the asymptotic variance of $\widetilde{\phi}_j(x)$ and $\widehat{\phi}_j(x)$ is equal, the bias is the crucial part for the better performance of the component-based estimator. Further, Table \ref{tab1} underlines that the bias term is more prominent the more complex the true process is.


\begin{table}[ht]
\centering
\setlength{\tabcolsep}{8pt}
\begin{tabular}{ r r r r r r r r r }
 DGP & $\rho$ & $n$ & $\widehat{\phi}_1$ & $\widetilde{\phi}_1$ & $\widehat{\phi}_2$ & $\widetilde{\phi}_2$ &$\widehat{\phi}_3$ & $\widetilde{\phi}_3$  \\
\hline
\multirow{8}{*}{Additive} & \multirow{4}{*}{$0$} & 300 & \bftab8.84 & 14.79 & \bftab8.53 & 8.72 & 3.05 & \fontseries{b}\selectfont{0.34} \\ 
  & & 500 & \fontseries{b}\selectfont{4.86} & 5.37 & \fontseries{b}\selectfont{5.47} & 6.30 & 1.91 & \fontseries{b}\selectfont{0.22} \\ 
  & & 1000 & \fontseries{b}\selectfont{3.04} & 3.29 & \fontseries{b}\selectfont{3.48} & 3.87 & 1.09 & \fontseries{b}\selectfont{0.11} \\ 
  & & 2000 & \fontseries{b}\selectfont{1.83} & 1.91 & \fontseries{b}\selectfont{2.13} & 2.43 & 0.66 & \fontseries{b}\selectfont{0.07} \\ 
\cline{2-9} 
& \multirow{4}{*}{$0.8$} & 300 & \fontseries{b}\selectfont{7.29} & 7.81 & \fontseries{b}\selectfont{7.99} & 8.87 & 2.69 & \fontseries{b}\selectfont{1.27} \\ 
 & & 500 & \fontseries{b}\selectfont{5.09} & 5.41 & \fontseries{b}\selectfont{5.65} & 6.05 & 1.86 & \fontseries{b}\selectfont{0.83} \\ 
 & & 1000 & \fontseries{b}\selectfont{3.37} & 3.38 & \fontseries{b}\selectfont{3.46} & 3.79 & 1.27 & \fontseries{b}\selectfont{0.52} \\ 
 & & 2000 & 2.24 & \fontseries{b}\selectfont{2.10} & \fontseries{b}\selectfont{2.22} & 2.37 & 0.87 & \fontseries{b}\selectfont{0.34} \\ 
\hline
\multirow{8}{*}{Non-linear} & \multirow{4}{*}{$0$} & 300 &  \fontseries{b}\selectfont{9.31} & 12.76 & \fontseries{b}\selectfont{11.45} & 14.15 & 1.74 & \fontseries{b}\selectfont{0.40} \\ 
 & & 500 & \fontseries{b}\selectfont{5.68} & 7.51 & \fontseries{b}\selectfont{6.74} & 7.86 & 0.87 & \fontseries{b}\selectfont{0.23} \\ 
 & & 1000 & \fontseries{b}\selectfont{3.20} & 4.20 & \fontseries{b}\selectfont{3.92} & 4.56 & 0.51 & \fontseries{b}\selectfont{0.13} \\ 
 & & 2000 & \fontseries{b}\selectfont{1.90} & 2.45 & \fontseries{b}\selectfont{2.22} & 2.63 & 0.30 & \fontseries{b}\selectfont{0.08} \\ 
\cline{2-9} 
& \multirow{4}{*}{$0.8$} & 300 & \fontseries{b}\selectfont{10.85} & 12.04 & \fontseries{b}\selectfont{11.89} & 14.34 & 3.16 & \fontseries{b}\selectfont{1.81} \\ 
 & & 500 & \fontseries{b}\selectfont{7.59} & 7.82 & \fontseries{b}\selectfont{7.98} & 9.79 & 2.20 & \fontseries{b}\selectfont{1.27}  \\ 
 & & 1000 & 4.99 & \fontseries{b}\selectfont{4.70} & \fontseries{b}\selectfont{4.79} & 5.94 & 1.50 & \fontseries{b}\selectfont{0.79} \\ 
 & & 2000 & 3.27 & \fontseries{b}\selectfont{2.86} & \fontseries{b}\selectfont{3.00} & 3.57 & 1.05 & \fontseries{b}\selectfont{0.49} \\ 
\hline
\hline
\end{tabular}
\caption{MISE of Shapley Curves for component-based estimator $\widehat{\phi}_j$ and integration-based estimator $\widetilde{\phi}_j$ for each variable. An additive and non-linear DGP is used for $\rho=0$ and $\rho=0.8$. The error terms follow $\varepsilon \sim \operatorname{N}(0, 1)$. \href{https://github.com/QuantLet/Shapley_Smoothing/tree/main/sim_consistency}{\hspace*{\fill} \raisebox{-1pt}{\includegraphics[scale=0.02]{qletlogo.png}}\,{normal\_sim}}}\label{tab1}
\end{table}




The Shapley curves for the third variable result in better performance when the integration-based estimation is applied. This observation might seem counterintuitive at first. As the pilot estimator $\widehat{m}(x)$ oversmooths in each dimension, the bandwidth of the third variable, $h_3$, is tuned to take the whole support of $X_3$. Computationally, this happens in every iteration of our Monte Carlo simulation. As a consequence, a linear model is fitted in $X_3$ direction, which takes the whole support for estimation. This reduces the variance as more observations are available for the fit. In contrast, the component-based estimation is not able to tune $h_3$ such that it captures the whole support in $\widehat{m}_s(x_s)$, for each subset that contains the third variable.

Note that we do not assume knowledge of the true process during the estimation of Shapley curves. If we conduct the estimation knowing that the first DGP is of additive structure, we can use backfitting or marginal integration. In this case, the optimal one-dimensional bandwidth for dimension $j$ would be used, instead of the full model bandwidth. Following Corollary \ref{col:additive} a faster decrease of the MISE will result since the convergence is one-dimensional.
Further, we conclude that the Shapley curves are robust to a lower signal-to-noise ratio as we obtain reasonable results for $\varepsilon \sim \operatorname{t}(5)$ in Table \ref{tab2}.
Both tables show that including dependence between the variables in form of correlation is introducing indirect effects as described in example 2. This leads to an increase in the MISE for (almost) all observation setups.



In the subsequent paragraphs, we are going to investigate the empirical coverage probability of the confidence intervals for the component-based and integration-based Shapley curves.
Let the true process be a non-linear function of $X_1$ and $X_2$:
\begin{align}
\label{bs_dgp}
m(x) = -\text{sin}(2x_1) + \text{cos}(3x_2) + 2\text{cos}(x_1)\text{sin}(2x_2).
\end{align}

Unfortunately, the established convergence rate is known to result in under-coverage for finite samples, due to its slow speed \citep{song2012bootstrap}. For this reason, we estimate the coverage probability by using a standard bootstrap approach. In the nonparametric literature, it is well known that bootstrap sampling yields improved finite sample coverage rather than directly estimating the confidence intervals relying on the asymptotic normal distribution \citep{hardle1991bootstrap, hardle1993comparing, hardle2015tie}.  

The coverage probability is estimated for $\alpha=0.05,0.1$ and $0.15$ for increasing sample size.
First, we simulate data following (\ref{bs_dgp}) with additional Gaussian noise. Then we sample $n$ observation pairs $\{(Y_i^b,X_i^b)\}^{n}_{i=1}$ with replacement from the initial data using the classical bootstrap procedure. On the original data, we estimate $m_s(x_s)$ for $s\in\mathcal{S}$ by local linear estimation, where we use leave-one-out cross validation to obtain the data-driven bandwidths $h_s$.
Next we estimate  $\widehat{m}^{\ast}_{b,s}(x_s)$ on the bootstrapped data set with $h_s$ as fixed bandwidths. 
We estimate $\widehat{\phi}^b_j(x_s)$ and $\widetilde{\phi}^b_j(x_s)$ at the point $(-0.5,-0.5)$. 
This procedure is repeated in $B=1000$ replications, where we estimate the component-based as well as the integration-based Shapley curve for each replication $b=1,\ldots,B$.
Finally, the upper and lower confidence bound is determined as the $1-\frac{\alpha}{2}$ and $\frac{\alpha}{2}$-quantiles, respectively. The estimated coverage probabilities for $M=1000$ Monte Carlo iterations are reported in Table \ref{comp_n_bs} for component-based estimation and $\varepsilon \sim \operatorname{N}(0, 1)$.

We almost obtain the desired coverage ratios. As expected, the local linear estimation is accompanied by a bias term which leads to a slight under-coverage in our simulation setup.
Table \ref{comp_t_bs} contains the estimated coverage probability for $t$-distributed error terms. As it shows, the estimation is robust to increased noise in the true process.
The appendix contains the results for integration-based Shapley curves in Table \ref{int_n_bs} for normally distributed errors and in Table \ref{int_t_bs} for $t$-distributed errors. The estimated coverage probability under-covers the true Shapley curve substantially. This is not surprising, since the bias term is inflated by over-smoothing in each direction.  

\begin{table}[ht]
\centering
\begin{tabular}{rrrrrrr}
  \hline
  $n$ & \multicolumn{3}{c}{Variable 1} & \multicolumn{3}{c}{Variable 2}\\\cmidrule{2-7}
 & 0.15 & 0.1 & 0.05 & 0.15 & 0.1 & 0.05 \\ 
  \hline
300 & 0.81 & 0.85 & 0.91 & 0.75 & 0.80 & 0.87 \\ 
  500 & 0.79 & 0.86 & 0.91 & 0.74 & 0.80 & 0.89 \\ 
  1000 & 0.79 & 0.84 & 0.92 & 0.78 & 0.85 & 0.90 \\ 
  2000 & 0.82 & 0.87 & 0.93 & 0.82 & 0.88 & 0.93 \\ 
  3000 & 0.86 & 0.91 & 0.96 & 0.81 & 0.86 & 0.93 \\ 
  4000 & 0.85 & 0.89 & 0.94 & 0.82 & 0.88 & 0.93 \\ 
   \hline
   \hline
\end{tabular}
\caption{Estimated coverage probability for component-based Shapley curves with $\varepsilon \sim \operatorname{N}(0, 1)$ and significance levels $\alpha=0.15, 0.1, 0.05$.\href{https://github.com/QuantLet/Shapley_Smoothing/tree/main/coverage}{\hspace*{\fill} \raisebox{-1pt}{\includegraphics[scale=0.02]{qletlogo.png}}\,{ncomp\_cov}}}\label{comp_n_bs}
\end{table}

\begin{table}[ht]
\centering
\begin{tabular}{rrrrrrr}
  \hline
  $n$ & \multicolumn{3}{c}{Variable 1} & \multicolumn{3}{c}{Variable 2}\\\cmidrule{2-7}
 & 0.15 & 0.1 & 0.05 & 0.15 & 0.1 & 0.05 \\ 
  \hline
300 & 0.82 & 0.88 & 0.93 & 0.73 & 0.79 & 0.87 \\ 
  500 & 0.80 & 0.84 & 0.91 & 0.74 & 0.81 & 0.89 \\ 
  1000 & 0.81 & 0.87 & 0.94 & 0.76 & 0.83 & 0.89 \\ 
  
  2000 & 0.81 & 0.87 & 0.93 & 0.80 & 0.87 & 0.92 \\ 
  3000 & 0.83 & 0.89 & 0.94 & 0.80 & 0.86 & 0.93 \\ 
  4000 & 0.83 & 0.88 & 0.93 & 0.85 & 0.88 & 0.94 \\ 
   \hline
   \hline
\end{tabular}
\caption{Estimated coverage probability for component-based Shapley curves with $\varepsilon \sim \operatorname{t}(5)$ and significance levels $\alpha=0.15, 0.1, 0.05$. \href{https://github.com/QuantLet/Shapley_Smoothing/tree/main/coverage}{\hspace*{\fill} \raisebox{-1pt}{\includegraphics[scale=0.02]{qletlogo.png}}\,{tcomp\_cov}}}\label{comp_t_bs}
\end{table}


Figures \ref{slice_plot1} and \ref{slice_plot2} illustrate the estimated Shapley curves for fixed $x_2=-0.5$ and the corresponding bootstrap confidence intervals. We use $B=10000$ replications for $n_b=n=2000$ bootstrap samples and the estimation procedure is conducted as before. We observe that component-based estimation results in good coverage. In contrary, the integration-based estimation under-covers for the same reason previously stated. 








\begin{figure}[H]
\vspace{-0.3cm}
\hspace{-0.8cm}
\begin{minipage}{0.5\textwidth}
		\includegraphics[scale=0.4]{bs_one_comp.png}
\end{minipage}
\hspace{0.3cm}
\begin{minipage}{0.5\textwidth}
		\includegraphics[scale=0.4]{bs_one_int.png}
\end{minipage}\\
\vspace{-1cm}
\caption{Bootstrap Confidence Intervals}
\label{slice_plot1}
  \centering
  \vspace{0.4cm}
  \begin{minipage}{15cm}%
\footnotesize
Population Shapley Curve (black) and based estimated Shapley Curve (blue) with $95\%$ confidence intervals (red) for the first variable. $x_1$ varies on (-2,2) and $x_2=-0.5$. Left: Component-based estimation. Right: Integration-based estimation. \href{https://github.com/QuantLet/Shapley_Smoothing/blob/main/coverage/CI_slice_plots/CI_slice_plots.R}{\hspace*{\fill} \raisebox{-1pt}{\includegraphics[scale=0.02]{qletlogo.png}}\,{CI\_slice}}
\end{minipage}
\end{figure}


\begin{figure}[H]
\vspace{-0.3cm}
\hspace{-0.8cm}
\begin{minipage}{0.5\textwidth}
		\includegraphics[scale=0.4]{bs_two_comp.png}
\end{minipage}
\hspace{0.3cm}
\begin{minipage}{0.5\textwidth}
		\includegraphics[scale=0.4]{bs_two_int.png}
\end{minipage}\\
\vspace{-1cm}
\caption{Bootstrap Confidence Intervals}
\label{slice_plot2}
  \centering
  \vspace{0.4cm}
  \begin{minipage}{15cm}%
\footnotesize
Population Shapley Curve (black) and integration-based estimated Shapley Curve (blue) with $95\%$ confidence intervals (red) for the first variable. $x_1$ varies on (-2,2) and $x_2=-0.5$. Left: Component-based estimation. Right: Integration-based estimation. \href{https://github.com/QuantLet/Shapley_Smoothing/blob/main/coverage/CI_slice_plots/CI_slice_plots.R}{\hspace*{\fill} \raisebox{-1pt}{\includegraphics[scale=0.02]{qletlogo.png}}\,{CI\_slice}}
\end{minipage}
\end{figure}

\newpage
\section{Empirical Application}
\label{section:application}
This empirical application illustrates our previous results in a real-world situation of vehicle price setting. Vehicle manufacturers typically rely on price-setting approaches that involve teardown data, surveying or indirect cost multiplier adjustments for calculating a markup. Similar to \citet{moawad2021explainable} we use a data-driven approach as a potential alternative. We use a sample of an extensive vehicle price data set for the U.S. provided by the latter authors and collected by the Argonne National Laboratory. Our data set includes the $20$ most important variables identified by \citet{moawad2021explainable} out of a larger pool of characteristics for the prediction of $38435$ vehicle prices ranging from year 2001 to year 2020.
However, our motivation differs in the sense that we are interested in the visualization of Shapley curves over the whole support of the variables at hand, instead of being restricted to have a variable importance measure only for the observations in the sample.
Our goal is not to obtain the most precise prediction of vehicle prices, but rather to provide an interpretation approach.


For this matter, we make use of the three most important variables of the data set, which are horsepower, vehicle weight (in pounds), and vehicle length (in inches). In principle, the dimension could be larger but does not contribute to an interpretable empirical application. Three univariate density plots of the data are given in Figure \ref{density}. The data is divided into three time intervals, ranging from 2001-2007 (12230), 2008-2013 (11955), and 2014-2020 (14250). The pooled summary statistics can be found in Table \ref{sum_stats}. In the following, the vehicle price is reported in $1000$ USD.
    
\begin{table}[H]
\centering
\begin{tabular}{rrrrrrr}
  \hline
 & Mean & Min & $q_{0.25}$ & $q_{0.5}$  & $q_{0.75}$ & Max \\ 
  \hline
Horsepower & 258 & 70 & 181 & 260 & 310 & 808 \\ 
  Weight [lbs] & 4233 & 1808 & 3322 & 3970 & 5042 & 8039 \\ 
  Length [in] & 199 & 106 & 180 & 192 & 219 & 290 \\ 
  Price [1000 USD] & 36.35 & 8.9 & 24.1 & 31.23 & 41.19 & 191.5 \\
   \hline
   \hline
\end{tabular}
\caption{Summary statistics.}\label{sum_stats}
\href{https://github.com/QuantLet/Shapley_Smoothing/tree/main/application_plots}{\hspace*{\fill} \raisebox{-1pt}{\includegraphics[scale=0.02]{qletlogo.png}}\,{Desc\_stats}}
\end{table}
    
    
   
 
    It is of empirical interest to consider Shapley curves as proposed in this work for several reasons. First, we are able to decompose the estimated conditional expectation locally at any point vector of interest. This type of investigation can contribute to an empirical understanding of prices for U.S. vehicle companies. Second, we obtain asymptotically valid confidence intervals around the estimated Shapley curves. This enables the price setter to differentiate between significant and non-significant support-sections of the covariates.
   
   In this application, we use component-based Shapley curves as proposed in the previous chapters. The bandwidth choice is motivated by a dual objective, on the one hand to get a good fit to the data, and on the other hand to get interpretable and smooth curves.. 
    

    
   
   In order to illustrate the Shapley curves in dependence of two variables, we are required to fix the third variable. The surface plot for the first variable, Horsepower (hp), is illustrated in Figure \ref{3d_emp} in dependence of horsepower and vehicle weight for three different time periods. Several observations can be made from these plots. First, we can illustrate the interactive contribution on the price prediction between the variables. Further, very light cars with high horsepower such as sports cars, are resulting in the highest increase in the price prediction. The lighter the car as well as the higher the horsepower, the higher the contribution to the price. Across the time domain, we infer that horsepower in general is contributing less to an increased price as we move to the second and third time interval. This effect is especially prominent as we compare the first and second time interval. The reason is that as technology develops over time, it is less costly for car manufacturers to produce vehicles with higher horsepower.
    

\begin{figure}[H]
  \begin{minipage}{0.5\linewidth}
    \includegraphics[scale=0.18]{first_x1_3d.png} 
    \vspace{-0.5cm}
  \end{minipage
  \begin{minipage}{0.5\linewidth}
    \hspace{-1.4cm}
    \includegraphics[scale=0.18]{second_x1_3d.png} 
  \end{minipage}
  \begin{minipage}{0.5\linewidth}
    \includegraphics[scale=0.18]{third_x1_3d.png} 
  \end{minipage
  \begin{minipage}{0.5\linewidth}
     \hspace{-1.4cm}
    \includegraphics[scale=0.18]{all_x1_3d.png} 
  \end{minipage} 
\caption{Shapley Curves for Horsepower}
\label{3d_emp}
  \centering
  \vspace{0.1cm}
  \begin{minipage}{13cm}%
\footnotesize 
Estimated Shapley curves for horsepower in dependence of vehicle weight (in lbs) and horsepower (in hp) for a vehicle length of $190$ inches. The time periods are 2001 - 2007 (top left), 2008 - 2013 (top right), 2014 - 2020 (bottom left) and pooled (bottom right). 
\href{https://github.com/QuantLet/Shapley_Smoothing/tree/main/application_plots}{\hspace*{\fill} \raisebox{-1pt}{\includegraphics[scale=0.02]{qletlogo.png}}\,{SC}}
\end{minipage}
\end{figure}    
    
  
    We now want to illustrate the results for the estimated confidence intervals for a variable $j$. Therefore we fix the remaining variables at a value that lies in a dense area of the data or at the median. The results of this exercise are conducted in Figure \ref{hp_slice}, \ref{lbs_slice} and \ref{in_slice}, for the variables horsepower, vehicle weight and vehicle length, respectively. For each \textit{'slice plot'} we illustrate $\widehat{\phi}_j(x)$ against $X_j$. 
    The Shapley curve for horsepower (Figure \ref{hp_slice}) is resulting in similar intuition as described previously. 
    The second Shapley curve (Figure \ref{lbs_slice}) shows that a lighter weight contributes to a price decrease, accompanied by stagnation at around 4000 lbs and an increase shortly after. Further, it shows that the Shapley curve is stretched apart on the tails across the time intervals. 
    The estimated confidence intervals have a larger spread as we move to the tails of the variables. This is intuitive as we have more variation in price for the observations in this area. For example, a very heavy car with more than 5500 lbs can either be a Nissan NV Cargo Van or a Cadillac Escalade SUV. On the lower tails, we observe insignificance for the first and the third period. 
    The third variable, vehicle length (in inches) is illustrated in Figure \ref{in_slice}. We observe that vehicle length (mostly) does not significantly contribute to the price prediction after a length of approximately 230 inches for the first and second time period. However, we observe a downward trend in the prices for longer cars the more recent the time interval is.

Finally, we are interested in the cumulative price contribution of the variables at hand. Therefore we sequentially accumulate the Shapley curves to the unconditional mean prediction, such that the additional contribution of a variable to the prediction is visible. For instance, the green area in Figure \ref{seq_plot} indicates an increase in the price prediction, after adding up a certain variable. We call this the \textit{cumulative Shapley curves}, which are a function of a single variable of interest, with fixed remaining variables.
For example, the first row of Figure \ref{seq_plot} includes the cumulative Shapley curves as a function of horsepower for a car length of 190 inches and 3500 pounds. As we see, adding the curves for vehicle weight and length does barely contribute to the prediction of price. In other words, their importance as a function of horsepower is negligible. 
The second row of Figure \ref{seq_plot} shows that as we include the second and third variable as a function of weight, we see that the price contribution mainly changes as we move closer to the tails. 
The first column of said figure indicates how much more the addition of horsepower explains price differences in comparison to the average prediction. If we plot it as a function of horsepower, it contributes a lot, which means the area is relatively large. A similar argument holds for the Shapley curve for horsepower, as a function of weight. However, for the third variable the corresponding area is rather small. In general, we observe that even adding the Shapley curve for weight and length does not contribute a lot to the price prediction as a function of length.

A qualitative comparison can be made to the additive model. First, note that the first row suggests that the price contribution of weight and length does not depend on horsepower. Assume that the unknown data-generating process indeed follows an additive model as defined in Assumption \ref{assum:additive}. In that case Figure \ref{seq_plot} would have relatively large areas in the diagonal plots and an almost non-existent area in the off-diagonals. We conclude that this graphical analysis of sequential Shapley curves enables a heuristic distinction to the additive model.









\begin{figure}[ht] 
  \begin{minipage}{0.3\linewidth}
    \includegraphics[scale=0.23]{seq_one_x1.png} 
  \end{minipage}
\hfill
  \begin{minipage}{0.3\linewidth}
    \includegraphics[scale=0.23]{seq_two_x1.png} 
  \end{minipage}
  \hfill
  \begin{minipage}{0.3\linewidth}
    \includegraphics[scale=0.23]{seq_three_x1.png} 
  \end{minipage}
  \begin{minipage}{0.3\linewidth}
    \includegraphics[scale=0.23]{seq_one_x2.png} 
  \end{minipage} 
  \hfill
  \begin{minipage}{0.3\linewidth}
    \includegraphics[scale=0.23]{seq_two_x2.png} 
  \end{minipage}
  \hfill
  \begin{minipage}{0.3\linewidth}
    \includegraphics[scale=0.23]{seq_three_x2.png} 
  \end{minipage}
    \begin{minipage}{0.3\linewidth}
    \includegraphics[scale=0.23]{seq_one_x3.png} 
  \end{minipage}
  \hfill
    \begin{minipage}{0.3\linewidth}
    \includegraphics[scale=0.23]{seq_two_x3.png} 
  \end{minipage}
  \hfill
    \begin{minipage}{0.3\linewidth}
    \includegraphics[scale=0.23]{seq_three_x3.png} 
  \end{minipage}
  \caption{Cumulative Shapley Curves}
\label{seq_plot}
  \centering
  \vspace{0.1cm}
  \begin{minipage}{13cm}%
\footnotesize
Green area indicates an increase and red area indicates decrease in the price prediction of updated cumulative Shapley curve (in blue).
Left panel shows $\bar{Y}$ in black and $\widehat{\phi}_1(x)+\bar{Y}$ in blue; Mid panel shows $\widehat{\phi}_1(x)+\bar{Y}$ in black and $\widehat{\phi}_1(x)+\widehat{\phi}_2(x) +\bar{Y}$ in blue; Right panel shows $\widehat{\phi}_1(x)+\widehat{\phi}_2(x) +\bar{Y}$ in black and $\widehat{\phi}_1(x) + \widehat{\phi}_2(x) + \widehat{\phi}_3(x) +\bar{Y}=\widehat{m}(x)$ in blue. The third plot in each row contains the conditional mean prediction in blue.
First row for vehicles with weight of 3500 lbs and length of 190 inches.
Second row for vehicles with horsepower of 250 and length of 190 inches.
Third row for vehicles with horsepower of 250 lbs and weight of 3500 lbs.
 \href{https://github.com/QuantLet/Shapley_Smoothing/tree/main/seq_plots}{\hspace*{\fill} \raisebox{-1pt}{\includegraphics[scale=0.02]{qletlogo.png}}\,{Seq\_plots}}
\end{minipage}
\end{figure}











\begin{figure}[ht] 
  \begin{minipage}{0.5\linewidth}
    \centering
    \includegraphics[scale=0.35]{first_x1.png} 
    \vspace{-0.2cm}
  \end{minipage
  \hspace{0.2cm}
  \begin{minipage}{0.5\linewidth}
    \centering
    \includegraphics[scale=0.35]{second_x1.png} 
    \vspace{-0.2cm}
  \end{minipage}
  \begin{minipage}{0.5\linewidth}
    \centering
    \includegraphics[scale=0.35]{third_x1.png} 
  \end{minipage
    \hspace{0.2cm}
  \begin{minipage}{0.5\linewidth}
    \centering
    \includegraphics[scale=0.35]{full_x1.png} 
  \end{minipage} 
  \caption{Shapley curves for horsepower}
\label{hp_slice}
  \centering
  \vspace{0.1cm}
  \begin{minipage}{13cm}%
\footnotesize
Estimated Shapley curves for horsepower in dependence of horsepower (in hp) for a vehicle length of $190$ inches and vehicle weight of $3500$ pounds. The time periods are 2001 - 2007 (top left), 2008 - 2013 (top right), 2014 - 2020 (bottom left) and pooled (bottom right). \href{https://github.com/QuantLet/Shapley_Smoothing/tree/main/application}{\hspace*{\fill} \raisebox{-1pt}{\includegraphics[scale=0.02]{qletlogo.png}}\,{HP\_slice}}
\end{minipage}
\end{figure}


\begin{figure}[ht] 
  \begin{minipage}{0.5\linewidth}
    \centering
    \includegraphics[scale=0.35]{first_x2.png} 
    \vspace{-0.2cm}
  \end{minipage
  \hspace{0.2cm}
  \begin{minipage}{0.5\linewidth}
    \centering
    \includegraphics[scale=0.35]{second_x2.png} 
    \vspace{-0.2cm}
  \end{minipage}
  \begin{minipage}{0.5\linewidth}
    \centering
    \includegraphics[scale=0.35]{third_x2.png} 
  \end{minipage
    \hspace{0.2cm}
  \begin{minipage}{0.5\linewidth}
    \centering
    \includegraphics[scale=0.35]{full_x2.png} 
  \end{minipage} 
\caption{Shapley curves for vehicle weight}
\label{lbs_slice}
  \centering
  \vspace{0.1cm}
  \begin{minipage}{13cm}%
\footnotesize
Estimated Shapley curves for vehicle weight in dependence of weight (in lbs) for a vehicle length of $190$ inches and $190$ horsepower. The time periods are 2001 - 2007 (top left), 2008 - 2013 (top right), 2014 - 2020 (bottom left) and pooled (bottom right). \href{https://github.com/QuantLet/Shapley_Smoothing/tree/main/application}{\hspace*{\fill} \raisebox{-1pt}{\includegraphics[scale=0.02]{qletlogo.png}}\,{Weight\_slide}}
\end{minipage}
\end{figure}


\begin{figure}[ht] 
  \begin{minipage}{0.5\linewidth}
    \centering
    \includegraphics[scale=0.35]{first_x3.png} 
    \vspace{-0.2cm}
  \end{minipage
  \hspace{0.2cm}
  \begin{minipage}{0.5\linewidth}
    \centering
    \includegraphics[scale=0.35]{second_x3.png} 
    \vspace{-0.2cm}
  \end{minipage}
  \begin{minipage}{0.5\linewidth}
    \centering
    \includegraphics[scale=0.35]{third_x3.png} 
  \end{minipage
    \hspace{0.2cm}
  \begin{minipage}{0.5\linewidth}
    \centering
    \includegraphics[scale=0.35]{full_x3.png} 
  \end{minipage} 
  \caption{Shapley curves for vehicle length}
\label{in_slice}
  \centering
  \vspace{0.1cm}
  \begin{minipage}{13cm}%
\footnotesize
Estimated Shapley curves for vehicle length in dependence of length (in) for a vehicle weight of $3500$ pounds and horsepower of $190$. The time periods are 2001 - 2007 (top left), 2008 - 2013 (top right), 2014 - 2020 (bottom left) and pooled (bottom right). \href{https://github.com/QuantLet/Shapley_Smoothing/tree/main/application}{\hspace*{\fill} \raisebox{-1pt}{\includegraphics[scale=0.02]{qletlogo.png}}\,{Length\_slice}}
\end{minipage}
\end{figure}



\begin{figure}[ht] 
  \begin{minipage}{0.5\linewidth}
    \centering
    \includegraphics[scale=0.35]{density_x1.png} 
    \vspace{-0.2cm}
  \end{minipage
  \hspace{0.2cm}
  \begin{minipage}{0.5\linewidth}
    \centering
    \includegraphics[scale=0.35]{density_x2.png} 
    \vspace{-0.2cm}
  \end{minipage}
  \begin{minipage}{0.5\linewidth}
    \centering
    \includegraphics[scale=0.35]{density_x3.png} 
  \end{minipage
  \caption{Density for three variables}
\label{density}
  \centering
  \vspace{0.1cm}
  \begin{minipage}{13cm}%
\footnotesize
Kernel density estimator for Horsepower, Vehicle Weight and Vehicle Length for the pooled data (2001-2021).   \href{https://github.com/QuantLet/Shapley_Smoothing/tree/main/application_plots}{\hspace*{\fill} \raisebox{-1pt}{\includegraphics[scale=0.02]{qletlogo.png}}\,{density\_desc}}
\end{minipage}
\end{figure}






\clearpage
\section{Conclusion}
\label{section:conclusion}

This paper analyzes Shapley curves as a local measure of variable importance in a nonparametric framework. We give a rigorous definition of Shapley curves on the population level. As for estimation, we discuss two estimation strategies, the component-based approach, and the integration-based approach. We show that both estimators globally converge with the nonparametric rate of \citet{stone1982optimal} to the population counterpart. Asymptotic normality is obtained for both estimators. We show that the asymptotic variance is identical for both approaches. However, the integration-based approach is accompanied by an inflated asymptotic bias. In our simulation exercise, this difference is visible in the results for the MISE and in the empirical coverage probabilities of the confidence intervals. The advantage of the integration-based estimator is that only one pilot estimator is required to estimate the components. Under the assumption of an additive model a one-dimensional convergence rate results for Shapley curves \citep{stone1985additive}.
For an extensive vehicle price data set we show that Shapley curves are a useful tool for the practitioner to gain insight into pricing and the importance of vehicle characteristics. The estimated confidence intervals enable to distinguish between significant and non-significant sections of the variables. Building on our results, empirical researchers as well as practitioners can conduct statistical inference for Shapley curves.











\acks{This research is supported by the Deutsche Forschungsgemeinschaft via IRTG 1792 “High Dimensional Nonstationary Time Series”, Humboldt-Universität zu Berlin; the European Union’s Horizon 2020 research and innovation program “FIN-TECH: A Financial supervision and Technology compliance training programme”under the grant agreement no. 825215 (Topic: ICT-35-2018, Type of action: CSA); the European Cooperation in Science \& Technology COST Action grant CA19130 - Fintech and Artificial Intelligence in Finance - Towards a transparent financial industry; the Yushan Scholar Program of Taiwan; and the Czech Science Foundation’s grant no. 19-28231X / CAS: XDA 23020303.}
































    
    
    
    


    

    







\clearpage
