Control variates, antithetics, and stratification

For $Z_i=f(X_i)-{\beta }^{\mathsf{T}}(h(X_i)-\theta )$, the expectation convention [F6] and linearity [F3] give $E{Z}_i=Pf-{\beta }^{\mathsf{T}}(Ph-\theta )=\mu$. Thus $E{\widehat{\mu }}_{\beta }=\mu$.

The $Z_i$ are the same measurable function of the IID variables $X_i$, so [F2] and [F4] make them IID.

By [F7]–[F8], the variance as a function of $\beta$ is ${\mathrm{Var}}_P(f)-2{\beta }^{\mathsf{T}}c+{\beta }^{\mathsf{T}}\mathrm{\Sigma }\beta$.

Expanding the centered square of $n^{-1}{\sum }_i{Z}_i$, the diagonal terms contribute $n^{-2}n\mathrm{Var}(Z_1)$. For $i\mathrm{\ne }j$, $Z_i-E{Z}_i$ and $Z_j-E{Z}_j$ are independent integrable variables, so [F5] factors their product expectation into zero. Hence $\mathrm{Var}({\widehat{\mu }}_{\beta })=n^{-1}\mathrm{Var}(Z_1)=n^{-1}{\mathrm{Var}}_P(f-{\beta }^{\mathsf{T}}h)$.

If $\mathrm{\Sigma }$ is positive definite, substituting $c=\mathrm{\Sigma }{\beta }^{\ast }$ and completing the square gives ${\mathrm{Var}}_P(f)-c^{\mathsf{T}}{\mathrm{\Sigma }}^{-1}c+(\beta -{\beta }^{\ast }{)}^{\mathsf{T}}\mathrm{\Sigma }(\beta -{\beta }^{\ast })$. The last term is nonnegative and vanishes only at $\beta ={\beta }^{\ast }$; for $q=1$ the displayed scalar formula follows.

In dimension one, let $\mu$ be the law of $U_1$. By [F7], the expectations defining the covariance are the corresponding $\mu$-integrals. Expanding those integrals using [F8] and applying [F3] gives $2\mathrm{Cov}(G(U_1),H(U_1))=\int \text{\hspace{0.17em}⁣}\int (G(x)-G(y))(H(x)-H(y))\mu (dx)\mu (dy)$. The integrand is nonpositive for every $x,y$ because the two functions have opposite monotonicity, so the covariance is nonpositive. The two differences are square-integrable on $({\mathbb{R}}^2,\mu \times \mu )$, so [F4] makes their product integrable and justifies the double integral.

Assume the claim holds for $d-1$ and first suppose $G,H$ are bounded. By [F2] the joint law is $\mu \times \nu$, where $\mu$ is the first $d-1$ coordinate laws and $\nu$ is the law of $U_d$. The change-of-variables formula [F7] identifies expectations under the joint law with integrals against $\mu \times \nu$. Define $G_0(x)=\int G(x,z)\nu (dz)$ and $H_0(x)=\int H(x,z)\nu (dz)$. Fubini [F3] makes these measurable and bounded; they retain their respective coordinate monotonicities. For each fixed $x$, step 1.1 applied to the last coordinate gives ${\mathrm{Cov}}_{\nu }(G(x,U_d),H(x,U_d))\le 0$, and the induction hypothesis gives ${\mathrm{Cov}}_{\mu }(G_0,H_0)\le 0$.

Expanding both sides by the covariance formula [F5], using expectation linearity [F8], the joint law identity [F7], and Fubini [F3] yields $\mathrm{Cov}(G(U),H(U))=\int {\mathrm{Cov}}_{\nu }(G(x,U_d),H(x,U_d))\mu (dx)+{\mathrm{Cov}}_{\mu }(G_0,H_0)\le 0$. Thus the bounded case holds in dimension $d$, completing the induction.

For general square-integrable $G,H$, clamp each function to $[-m,m]$. The clamps remain coordinatewise monotone and bounded, with $\mathrm{\mid }{G}_m\mathrm{\mid }\le \mathrm{\mid }G\mathrm{\mid }$, $\mathrm{\mid }{H}_m\mathrm{\mid }\le \mathrm{\mid }H\mathrm{\mid }$, $\mathrm{\mid }{G}_m-G{\mathrm{\mid }}^2\le G^2$, and $\mathrm{\mid }{H}_m-H{\mathrm{\mid }}^2\le H^2$. Applying [F6] to the squared differences gives $L^2$ convergence, while $\mathrm{\parallel }{G}_m{\mathrm{\parallel }}_2\le \mathrm{\parallel }G{\mathrm{\parallel }}_2$ and $\mathrm{\parallel }{H}_m{\mathrm{\parallel }}_2\le \mathrm{\parallel }H{\mathrm{\parallel }}_2$. Cauchy–Schwarz [F4] shows that $E\mathrm{\mid }{G}_m-G\mathrm{\mid }$ and $E\mathrm{\mid }{H}_m-H\mathrm{\mid }$ tend to zero, and that $E\mathrm{\mid }{G}_m{H}_m-GH\mathrm{\mid }\le \mathrm{\parallel }{G}_m-G{\mathrm{\parallel }}_2\mathrm{\parallel }{H}_m{\mathrm{\parallel }}_2+\mathrm{\parallel }G{\mathrm{\parallel }}_2\mathrm{\parallel }{H}_m-H{\mathrm{\parallel }}_2$ tends to zero. Thus expectation linearity [F8] and the covariance identity [F5] imply $\mathrm{Cov}(G_m(U),H_m(U))\to \mathrm{Cov}(G(U),H(U))$. Each bounded covariance is nonpositive by step 3.1, so the limit is nonpositive.

For a uniform coordinate, $\mathbb{P}(U_{1j}=s)=0$ for every $s\in [0,1]$, since $\mathbb{P}(U_{1j}\le s)=s$ and the left limit of this cdf at $s$ is also $s$. Thus for $0\le t\le 1$, $\mathbb{P}(1-U_{1j}\le t)=\mathbb{P}(U_{1j}\ge 1-t)=1-\mathbb{P}(U_{1j}<1-t)=t$. By [F3], the reflected coordinates remain independent; [F4] then shows $1-U_1$ has the same product-uniform law as $U_1$. Therefore both $T(U_1)$ and $T(1-U_1)$ have law $P$. Applying [F10] to $f$ and $f^2$ identifies their means and variances as $\mu$ and ${\sigma }^2$, so $Y$ and $\tilde{Y}$ have the same law, mean $\mu$, and variance ${\sigma }^2$.

The pair value $Z_i=[f(T(U_i))+f(T(1-U_i))]\mathrm{/}2$ is a measurable function of $U_i$. Thus the $Z_i$ are IID by [F2]–[F3], each has mean $\mu$ by [F9], and their average is unbiased.

By [F5]–[F6], $\mathrm{Var}(Z_i)=[\mathrm{Var}(Y)+\mathrm{Var}(\tilde{Y})+2\mathrm{Cov}(Y,\tilde{Y})]\mathrm{/}4=[{\sigma }^2+\mathrm{Cov}(Y,\tilde{Y})]\mathrm{/}2$. For distinct $i,j$, the centered $Z_i,Z_j$ are independent, so [F7] gives zero covariance; expanding the centered square of their average therefore gives the stated variance divided by $n$.

The mean of $2n$ independent $P$-draws has variance ${\sigma }^2\mathrm{/}(2n)$ by the same centered-square expansion and [F7]. Comparing this with step 2.1 gives the iff conditions for weak and strict variance reduction.

For each coordinate choose a monotonicity direction for $g$: define $S_j(v)=v_j$ in a nondecreasing coordinate and $S_j(v)=1-v_j$ in a nonincreasing coordinate; a constant coordinate may be assigned either direction. Then $S$ is a coordinatewise reflection map with $S(S(v))=v$ and $S(1-v)=1-S(v)$. By [F3] and [F4], $V=S(U_1)$ again has independent uniform coordinates, and $\varphi =g\circ S$ is coordinatewise nondecreasing on $[0,1{]}^d$. Define the coordinatewise clamp $C:{\mathbb{R}}^d\to [0,1{]}^d$ by $C_j(x)=\min \{1,\max \{0,x_j\}\}$. It is continuous and nondecreasing in every coordinate, and $C(1-v)=1-C(v)$. Hence $\bar{\varphi }:=\varphi \circ C$ is a Borel coordinatewise nondecreasing function on ${\mathbb{R}}^d$, while $v\mapsto \bar{\varphi }(1-v)$ is Borel and coordinatewise nonincreasing. On the cube these extensions equal $\varphi (v)$ and $\varphi (1-v)$, respectively; the identities for $S$ give $g(U_1)=\varphi (V)$ and $g(1-U_1)=\varphi (1-V)$. By step 1.1, these functions are square-integrable, so [F8] applies to $\bar{\varphi }(V)$ and $\bar{\varphi }(1-V)$ and gives $\mathrm{Cov}(Y,\tilde{Y})\le 0$. In the stated ordered-subset case, $u\le v$ coordinatewise implies $T(u)\le T(v)$ and then $f(T(u))\le f(T(v))$, which makes $g$ nondecreasing and is the special case $S(v)=v$.

For each $h$, [F1] gives ${\int }_{\mathsf{X}}{f}^{2}d{P}_h=p_h^{-1}{\int }_{A_h}{f}^{2}dP<\mathrm{\infty }$. Since $\mathrm{\mid }f\mathrm{\mid }\le 1+f^2$, $f$ is integrable under $P$ and every $P_h$. Thus [F3] gives ${\mu }_h=p_h^{-1}{\int }_{A_h}fdP$, and [F4] plus the finite partition identity yield $E{\widehat{\mu }}_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{t}}={\sum }_h{p}_h{\mu }_h={\sum }_h{\int }_{A_h}fdP=\int fdP=\mu$.

For fixed $h$, expand the centered square of ${\bar{f}}_h$. The diagonal terms give ${\sigma }_h^2\mathrm{/}{n}_h$; for distinct draws the centered product expectation is zero by independence and [F5]. Thus $\mathrm{Var}({\bar{f}}_h)={\sigma }_h^2\mathrm{/}{n}_h$.

Splitting $\int (f-\mu {)}^{2}dP$ over the partition and writing $f-\mu =(f-{\mu }_h)+({\mu }_h-\mu )$ on $A_h$, the cross term integrates to zero since ${\int }_{A_h}(f-{\mu }_h)dP=p_h{E}_{P_h}(f-{\mu }_h)=0$. Therefore ${\mathrm{Var}}_P(f)={\sum }_h{p}_h{\sigma }_h^2+{\sum }_h{p}_h({\mu }_h-\mu {)}^2$.

By [F1] and [F7], the sigma-algebras generated by all individual draws are independent. Grouping them by stratum and applying [F8] shows that the sigma-algebras generated by the disjoint stratum blocks are independent. Each ${\bar{f}}_h$ is measurable with respect to its block sigma-algebra because it is a finite sum of measurable values of $f$, using [F9]; hence the stratum averages are independent. Applying [F5] to their centered products and expanding the variance of ${\sum }_h{p}_h{\bar{f}}_h$ gives ${\sum }_h{p}_h^2{\sigma }_h^2\mathrm{/}{n}_h$.

Under $n_h=n{p}_h$, step 3.1 becomes $n^{-1}{\sum }_h{p}_h{\sigma }_h^2$. The ordinary mean of $n$ IID $P$-draws has variance ${\mathrm{Var}}_P(f)\mathrm{/}n$, by the same centered-square expansion as step 2.1. Step 2.2 shows the difference is $n^{-1}{\sum }_h{p}_h({\mu }_h-\mu {)}^2\ge 0$, proving the comparison.

The entries $h_j(X_i)$, $f(X_i)$, $h_j(X_i)h_k(X_i)$, and $h_j(X_i)f(X_i)$ are measurable by [F11]. Since each $X_i$ has law $P$, applying the same measurable function to each gives a common law; [F10] preserves their mutual independence, so each resulting sequence is IID by [F2]. The fourth-moment hypothesis and [F8] make every product entry integrable, and [F9] identifies its expectation with the corresponding integral under $P$. Applying [F4] to these IID sequences gives ${\bar{h}}_n\to \theta$, ${\bar{f}}_n\to Pf$, $n^{-1}{\sum }_{i}h(X_i)h(X_i{)}^{\mathsf{T}}\to E_P[h{h}^{\mathsf{T}}]$, and $n^{-1}{\sum }_{i}h(X_i)f(X_i)\to E_P[hf]$ almost surely. The centered formulas therefore give ${\widehat{\mathrm{\Sigma }}}_n\to \mathrm{\Sigma }$ and ${\widehat{c}}_n\to c$ almost surely; positive definiteness of $\mathrm{\Sigma }$ and continuity of determinant imply ${\widehat{\mathrm{\Sigma }}}_n$ is eventually invertible almost surely. The fallback definition makes ${\widehat{\beta }}_n$ measurable: on the open nonsingular set the inverse entries are ratios of polynomial functions, and on its closed complement the value is constant.

Put $W_i=h(X_i)-\theta$ and $V_i=f(X_i)-Pf$. The finite fourth-moment assumptions imply finite fourth moments for every component of $W_i$ and for $V_i$; for example, $\mathrm{\mid }x-a{\mathrm{\mid }}^4\le 8(\mathrm{\mid }x{\mathrm{\mid }}^4+\mathrm{\mid }a{\mathrm{\mid }}^4)$. Then [F8] gives finite second moments for every product $W_{ij}{W}_{ik}$ and $W_{ij}{V}_i$. Each centered coordinate and product is a fixed measurable function of $X_i$, so [F2] and [F10] make each sequence IID. By [F9] and [F7], $E{W}_i=0$, $E{V}_i=0$, and the product means are $\mathrm{\Sigma }$ and $c$. Independence and [F6] give variance equal to the one-observation variance divided by $n$ for each sample average. Applying [F5] with $a=M\mathrm{/}\sqrt{n}$ and taking a finite union over the fixed $q$ coordinates and matrix entries proves ${\bar{h}}_n-\theta =O_P(n^{-1\mathrm{/}2})$, ${\bar{f}}_n-Pf=O_P(n^{-1\mathrm{/}2})$, $n^{-1}{\sum }_i{W}_i{W}_i^{\mathsf{T}}-\mathrm{\Sigma }=O_P(n^{-1\mathrm{/}2})$, and $n^{-1}{\sum }_i{W}_i{V}_i-c=O_P(n^{-1\mathrm{/}2})$. The identities ${\widehat{\mathrm{\Sigma }}}_n=\frac{1}{n}{\sum }_i{W}_i{W}_i^{\mathsf{T}}-({\bar{h}}_n-\theta )({\bar{h}}_n-\theta {)}^{\mathsf{T}},{\widehat{c}}_n=\frac{1}{n}{\sum }_i{W}_i{V}_i-({\bar{h}}_n-\theta )({\bar{f}}_n-Pf)$ show that the centering corrections are $O_P(n^{-1})$. Thus ${\widehat{\mathrm{\Sigma }}}_n-\mathrm{\Sigma }$ and ${\widehat{c}}_n-c$ are jointly $O_P(n^{-1\mathrm{/}2})$.

On the event that ${\widehat{\mathrm{\Sigma }}}_n$ is invertible, ${\widehat{\beta }}_n-{\beta }^{\ast }={\widehat{\mathrm{\Sigma }}}_n^{-1}[({\widehat{c}}_n-c)-({\widehat{\mathrm{\Sigma }}}_n-\mathrm{\Sigma }){\beta }^{\ast }]$, since $c=\mathrm{\Sigma }{\beta }^{\ast }$. Step 1.1 and the adjugate-over-determinant formula show that the inverse converges to ${\mathrm{\Sigma }}^{-1}$ almost surely on the eventually nonsingular set; extending it by zero on the singular event gives an $O_P(1)$ matrix sequence. Given $\epsilon >0$, choose bounds for this inverse norm and for the bracket in step 1.2 after scaling by $\sqrt{n}$, with total failure probability less than $\epsilon \mathrm{/}2$. On their intersection the product is bounded after scaling by $\sqrt{n}$, so ${\widehat{\beta }}_n-{\beta }^{\ast }=O_P(n^{-1\mathrm{/}2})$ there. The singular-event probability tends to zero by step 1.1; the fallback value $0$ therefore does not change this order in probability.

By the plug-in definition in the Statement and [F1] for the fixed coefficient ${\beta }^{\ast }$, subtracting the two estimators cancels ${\bar{f}}_n$ and leaves exactly $-({\widehat{\beta }}_n-{\beta }^{\ast }{)}^{\mathsf{T}}({\bar{h}}_n-\theta )$. Given $\epsilon >0$, choose finite bounds for each factor after scaling by $\sqrt{n}$ so that the probability either bound fails is less than $\epsilon \mathrm{/}2$; on their intersection, Cauchy--Schwarz in ${\mathbb{R}}^q$ bounds $n$ times the absolute inner product by the product of those bounds. Thus the inner product is $O_P(n^{-1})$ by definition.

The production sample means are ${\mathcal{X}}_n$-measurable and satisfy $E\bar{f}=E_{P}f$ and $E(\bar{h}-\theta )=0$ by [F2], [F5], and [F6]. Each ${\widehat{\beta }}_j$ is $\mathcal{P}$-measurable; for any Borel sets, its preimage event belongs to $\mathcal{P}$ while the corresponding production-mean event belongs to ${\mathcal{X}}_n$, so [F7] makes it independent of each production mean. Since $E\mathrm{\parallel }\widehat{\beta }\mathrm{\parallel }<\mathrm{\infty }$, [F4] applied to the absolute product makes every ${\widehat{\beta }}_j({\bar{h}}_j-{\theta }_j)$ integrable; the full estimator is therefore integrable.

Let $A\in \mathcal{P}$. The random variable $\bar{f}$ is ${\mathcal{X}}_n$-measurable, so its Borel preimage events are independent of $A$ by [F7]; factorization [F4] gives $E[1_A\bar{f}]=P(A)E_{P}f$. Also $1_A{\widehat{\beta }}_j$ is $\mathcal{P}$-measurable and independent of ${\bar{h}}_j-{\theta }_j$, which is ${\mathcal{X}}_n$-measurable, so [F4] gives $E[1_A{\widehat{\beta }}_j({\bar{h}}_j-{\theta }_j)]=E[1_A{\widehat{\beta }}_j]E({\bar{h}}_j-{\theta }_j)=0$. By [F1] and [F5], $E[1_A{\widehat{\mu }}_{\widehat{\beta }}]=P(A)E_{P}f$.

The constant $E_{P}f$ is $\mathcal{P}$-measurable and its integral on every $A\in \mathcal{P}$ equals the integral of the estimator on $A$ by step 2.1. The defining characterization [F3] therefore gives $E[{\widehat{\mu }}_{\widehat{\beta }}\mid \mathcal{P}]=E_{P}f$ almost surely.