Monte Carlo permutation Shapley effects have a joint CLT

Fix $i$ and $S\subseteq [d]\setminus \{i\}$. Exactly $\mathrm{\mid }S\mathrm{\mid }!(d-\mathrm{\mid }S\mathrm{\mid }-1)!$ permutations place precisely $S$ before $i$: order $S$ first, place $i$, and order the remaining coordinates last. Dividing by $d!$ gives the coefficient of $v(S\cup \{i\})-v(S)$ in the expectation of ${\mathrm{\Delta }}_i$. This is exactly the Shapley weight in [F1], so $\mathbb{E}{\mathrm{\Delta }}_i=V{\varphi }_i$ and the estimator is unbiased.

The range and monotonicity of $v$ in [F3] give $0\le {\mathrm{\Delta }}_i\le V$ for every permutation. Thus the increment vector is bounded and has finite covariance. Each increment coordinate and each pairwise product is the same measurable function of each sampled permutation; [F5] gives independence, and the common input law gives identical distributions. Apply [F4] to their coordinate and product sequences; the sample mean and covariance converge almost surely to their population values.

The bounded increment vector satisfies the finite second-moment hypothesis of [F6]. The multivariate CLT gives $\sqrt{n}(\bar{\mathrm{\Delta }}-\mathbb{E}\mathrm{\Delta })\Rightarrow N(0,\mathrm{Cov}(\mathrm{\Delta }))$. Scaling by the fixed positive $V$ and using step 1.1 gives the stated joint limit for $\widehat{\varphi }$. Step 1.2 gives consistency of the covariance and standard-error estimates.

For each of the $n$ rows, the procedure traverses $d$ entries and requests two coalition values per entry. The loops are finite, so it terminates after $n$ rows and makes at most $2nd$ oracle calls. Computing the $d$ sample means and the $d\times d$ sample covariance from these rows takes $O(n{d}^2)$ arithmetic operations; this excludes permutation generation and oracle costs.