Pick-freeze Sobol ratios are consistent and have a joint delta-method CLT

For each row, Cauchy–Schwarz [F9] and $f\in L^4$ make $f(X)f(X^{[i]})$ and $f(U)f(V)$ square-integrable; $(a-b{)}^4\le 8(a^4+b^4)$ makes $T_i$ and $D$ square-integrable as well. Thus every coordinate of $Z$ has finite second moment. By [F10]–[F11], $\mathbb{E}[f(U)f(V)]={\mu }^2$ and $\mathbb{E}D=\frac{1}{2}(2\mathbb{E}[f(X{)}^2]-2{\mu }^2)=V$ by [F8]. The covariance identity [F2] gives $\mathbb{E}[f(X)f(X^{[i]})]-{\mu }^2=V_{\{i\}}$ and $\mathbb{E}{T}_i={\tau }_i$, hence $\mathbb{E}{N}_i=V_{\{i\}}$ and $\mathbb{E}Z=\eta$. Independent rows make these row vectors iid.

Apply [F3] to each of the finitely many coordinates of $Z$. Then $\bar{Z}\to \eta$ almost surely. Since $V>0$, almost surely its last coordinate is eventually positive, so the continuous ratio formula gives ${\tilde{S}}_n=g(\bar{Z})\to g(\eta )$ almost surely and $\mathbb{P}(\bar{D}=0)\to 0$. Every fixed extension on the zero-status event differs only on an event whose probability tends to zero; for any bounded continuous test function the expectation difference is at most twice its supremum norm times that probability. Hence it has the same distributional limit, and convergence in probability is also unchanged.

The finite second moment from step 1.1 and [A1] satisfy [F6]. Hence $\sqrt{n}(\bar{Z}-\eta )\Rightarrow N(0,\mathrm{\Sigma })$. The map $g$ in the Statement is Borel, being the continuous ratio map on the open half-space $s>0$ and zero on its complement. Since $V>0$, it is smooth near $\eta$ with derivative $J$. Applying [F7] gives the joint delta-method limit, and step 2.1 transfers it to every stated fixed extension.

For every two coordinates $Z_a,Z_b$, Cauchy–Schwarz gives $\mathbb{E}\mathrm{\mid }{Z}_a{Z}_b\mathrm{\mid }<\mathrm{\infty }$. Each product is the same measurable function of every iid row; [F4] gives independence and the common row law gives identical distributions. Thus [F5] gives convergence in probability of the sample means of all finitely many products $Z_a{Z}_b$. Step 2.1 gives convergence of the row means. Expanding the sample covariance around those means by [F11] proves $\widehat{\mathrm{\Sigma }}\to \mathrm{\Sigma }$ in probability. Continuity of the derivative of $g$ near $\eta$ then gives $\widehat{J}\widehat{\mathrm{\Sigma }}{\widehat{J}}^{\mathsf{T}}\to J\mathrm{\Sigma }{J}^{\mathsf{T}}$ in probability. Since the algorithm returns $M_n=\widehat{J}\widehat{\mathrm{\Sigma }}{\widehat{J}}^{\mathsf{T}}\mathrm{/}n$, this proves the stated scaled covariance convergence. An arbitrary fixed finite extension on the zero-status event preserves convergence in probability because that event has probability tending to zero by step 2.1. The argument does not claim an unbiased ratio.

The loops over $n$ rows and $d$ coordinates are finite, so the procedure terminates after the $n$ rows. In each row [F1] evaluates $f(X)$ once and reuses it in every hybrid expression, evaluates the $d$ first-order and $d$ total-effect hybrids, and evaluates the four additional values $f(U),f(V),f(W),f(W^{\mathrm{\prime }})$. This requires $1+d+d+4=2d+5$ evaluations per row, hence at most $n(2d+5)$ overall. Forming the covariance of the $(2d+1)$-dimensional row vectors costs $O(n{d}^2)$ scalar operations, and the plug-in matrix products cost $O(d^3)$.