Parametric Bootstrap Calibration With Estimated Margins

Why estimated margins matter

A copula separates marginal distributions from dependence. This is useful for testing independence, but it does not make marginal estimation disappear. If the marginal distributions are unknown, the test statistic can depend on how their parameters are estimated and on how the data are transformed using those estimates. A calibration that treats estimated margins as though they were known can therefore mimic the wrong procedure.

The parametric bootstrap addresses this by simulating data under a fitted independence model and repeating the analysis from the beginning on each simulated sample. In particular, each replicate must include the same marginal refitting and copula fitting used to calculate the observed statistic. The resulting p-value is a plug-in calibration for a composite null, not generally an exact finite-sample test.

Model and test statistic

Definition: Let $(X_i,Y_i)$, $i=1,\ldots,n$, be independent and identically distributed pairs with continuous parametric margins $F_1(\cdot;\theta_1)$ and $F_2(\cdot;\theta_2)$, densities $f_1$ and $f_2$, and a copula density $c_\eta$. The joint density is $$ f(x,y;\theta_1,\theta_2,\eta) =c_\eta\!\left(F_1(x;\theta_1),F_2(y;\theta_2)\right) f_1(x;\theta_1)f_2(y;\theta_2). $$ Independence is the null hypothesis $H_0:\eta=\eta_0$, where $c_{\eta_0}(u,v)=1$ for $0<u,v<1$. The marginal parameters are nuisance parameters: they are unknown but unrestricted by the independence hypothesis.

For a likelihood-ratio test, write the full log-likelihood as

\ell(\theta_1,\theta_2,\eta) =\sum_{i=1}^n \left[ \log c_\eta\!\left(F_1(X_i;\theta_1),F_2(Y_i;\theta_2)\right) +\log f_1(X_i;\theta_1)+\log f_2(Y_i;\theta_2) \right].

The likelihood-ratio statistic is

T_{\mathrm{obs}} =2\left\{ \sup_{\theta_1,\theta_2,\eta}\ell(\theta_1,\theta_2,\eta) -\sup_{\theta_1,\theta_2}\ell(\theta_1,\theta_2,\eta_0) \right\}.

The suprema are over the parameter spaces specified by the model. Under $H_0$, the copula density is one, so the log-likelihood reduces to the sum of the two marginal log-likelihoods. Thus the null estimates $\hat\theta_{1,0}$ and $\hat\theta_{2,0}$ are obtained by fitting the margins separately, provided their parameter spaces and likelihoods factor in this way. Under the alternative, the copula term generally depends on both marginal parameters too, so the full alternative likelihood must be maximized over $\theta_1,\theta_2,\eta$.

The bootstrap procedure

The bootstrap null distribution is generated from the fitted null model, not from the fitted alternative. Each simulated pair has independent components, with each component drawn from its corresponding fitted margin. Then the entire statistic calculation is repeated on that pair of bootstrap samples.

Fit the null model.
Estimate each marginal parameter under $H_0$, obtaining $\hat\theta_{1,0}$ and $\hat\theta_{2,0}$. Calculate the observed statistic using the prescribed null and alternative fits.

Generate a replicate under fitted independence.
For $i=1,\ldots,n$, independently generate $X_i^*\sim F_1(\cdot;\hat\theta_{1,0})$ and $Y_i^*\sim F_2(\cdot;\hat\theta_{2,0})$, with independence also across pairs.

Refit the replicate.
Fit both margins under the replicate null. Then fit the alternative copula model and its marginal parameters using the same method, constraints, and optimization rules as for the observed data.

Recalculate the statistic.
Compute $T_b^*$ from the replicate's null and alternative fits. Repeat the generation, refitting, and calculation for $b=1,\ldots,B$.

Calculate the Monte Carlo p-value.
For an upper-tail statistic, count $R=\sum_{b=1}^B \mathbf{1}\{T_b^*\geq T_{\mathrm{obs}}\}$ and use $\hat p=(R+1)/(B+1)$.

The refitting step is essential. For example, the transformed values $F_j(X_i^*;\hat\theta_{j,0})$ used to define a copula observation are based on fitted parameters; the corresponding transformed values in a replicate must use that replicate's fitted margins if that is how the observed statistic was constructed. Holding the original estimates fixed instead changes the procedure whose null distribution is being approximated.

Why the Monte Carlo correction is used

First consider the simpler case where the null parameters are known and the observed sample and $B$ simulated samples are generated independently from exactly the same null distribution. If the statistic has a continuous distribution, the rank of the observed statistic among the $B+1$ statistics is uniform. Therefore, the upper-tail rank divided by $B+1$ is a valid Monte Carlo p-value. With ties, counting simulated statistics greater than or equal to the observed one gives a conservative version. The $+1$ convention prevents a reported p-value of zero just because none of a finite number of replicates exceeded the observation.

Result: If the null parameters are known, the observed sample and bootstrap replicates are exchangeable under $H_0$, and ties are counted using $\geq$, then $(R+1)/(B+1)$ is super-uniform: its probability of being at most $\alpha$ is no greater than $\alpha$.

To see this, let $K=B+1$ and rank the observed statistic among all $K$ statistics, counting ties conservatively in the upper tail. Exchangeability makes each observation's position equally likely in the no-tie case; the upper-tail rank divided by $K$ is then uniform on $\{1/K,\ldots,1\}$. With ties, the conservative tail count can only move the p-value upward. Thus the probability of a p-value at most $\alpha$ is bounded by $\alpha$. $\square$

When the null parameters are estimated, the observed data are not exchangeable with replicates generated from the fitted model: the latter use $\hat\theta_{1,0},\hat\theta_{2,0}$ as if they were the true values. The same Monte Carlo formula is still used, but it no longer gives general finite-sample exactness. Its justification is typically asymptotic and depends on the model and statistic.

Worked Example: Gaussian Margins and a Likelihood-Ratio Statistic

Suppose the model has normal margins and a Gaussian copula, equivalently a bivariate normal distribution. For four observed pairs, take $x=(-1,-1,1,1)$ and $y=(-2,0,0,2)$. Both sample means are zero. The maximum-likelihood marginal variances under independence are

\hat\sigma_x^2=\frac{1+1+1+1}{4}=1,\qquad \hat\sigma_y^2=\frac{4+0+0+4}{4}=2.

The maximum-likelihood covariance estimates use divisor $n$, so

\hat\sigma_{xy}=\frac{2+0+0+2}{4}=1,\qquad \hat\rho=\frac{\hat\sigma_{xy}}{\hat\sigma_x\hat\sigma_y} =\frac{1}{\sqrt{2}}.

The null covariance determinant is $1\cdot2=2$; the unrestricted fitted covariance determinant is $1\cdot2(1-\hat\rho^2)=1$. Consequently,

T_{\mathrm{obs}} =4\log\!\left(\frac{2}{1}\right) =-4\log(1-\hat\rho^2) =4\log 2 \approx 2.7726.

For each bootstrap replicate, generate four independent $X_i^*\sim N(0,1)$ values and four independent $Y_i^*\sim N(0,2)$ values. Re-estimate both means and variances, calculate the replicate correlation using the same maximum-likelihood covariance convention, and evaluate $-4\log(1-\hat\rho_*^2)$. If $B=999$ and 63 replicate statistics are at least 2.7726, the bootstrap p-value is $(63+1)/(999+1)=0.064$. The original estimates are used to generate samples; replicate-specific estimates are used to calculate each replicate statistic.

Worked Example: Refitting Exponential Margins in a Replicate

Suppose the null model has independent exponential margins, with rate parameters $\lambda_1,\lambda_2$, and the observed samples are $x=(1,2,1,4)$ and $y=(2,1,3,2)$. The exponential maximum-likelihood rate estimate is $n/\sum_i x_i$, so both null estimates equal $4/8=0.5$. Generate each bootstrap sample independently from $\operatorname{Exp}(0.5)$ in each margin.

Consider one possible replicate $x^*=(0.5,1,2,1.5)$, $y^*=(1,2,0.5,1.5)$. Both sums equal 5, so the replicate null estimates are

\hat\lambda_{1,0}^*=\frac{4}{5}=0.8,\qquad \hat\lambda_{2,0}^*=\frac{4}{5}=0.8.

For an exponential rate $\lambda$, the distribution function is $F(t;\lambda)=1-e^{-\lambda t}$. The replicate's transformed values therefore include

u^*=\left(1-e^{-0.4},1-e^{-0.8},1-e^{-1.6},1-e^{-1.2}\right), \quad v^*=\left(1-e^{-0.8},1-e^{-1.6},1-e^{-0.4},1-e^{-1.2}\right).

These are calculated with the replicate estimates $0.8$, not the original estimates $0.5$. To obtain $T^*$, fit the alternative copula and its marginal parameters to this replicate by the same full-likelihood procedure used for the observed data, fit its null margins, and take twice the alternative-minus-null maximized log-likelihood. The transformed values alone are not the bootstrap statistic; they are inputs to that refitting procedure.

What plug-in calibration can guarantee

A standard sufficient route to asymptotic validity is that the conditional bootstrap distribution consistently estimates the null limiting distribution of the statistic. Let $G$ be a continuous limiting cumulative distribution function, and let $G_n^*$ be the conditional bootstrap cumulative distribution function given the observed data. Suppose, under the null,

T_{\mathrm{obs}}\xrightarrow{d}G,\qquad \sup_t|G_n^*(t)-G(t)|\xrightarrow{p}0.

Here $T_{\mathrm{obs}}\xrightarrow{d}G$ means convergence in distribution to a random variable with cumulative distribution function $G$. The bootstrap upper-tail probability $1-G_n^*(T_{\mathrm{obs}})$ is then asymptotically uniform. Indeed, its difference from $1-G(T_{\mathrm{obs}})$ is bounded by $\sup_t|G_n^*(t)-G(t)|$, which converges in probability to zero. By continuity of $G$, the probability integral transform gives $G(T_{\mathrm{obs}})\xrightarrow{d}\operatorname{Uniform}(0,1)$; hence the upper-tail probability also converges in distribution to a uniform random variable. $\square$

This result states a sufficient condition, not a guarantee that every fitted parametric bootstrap satisfies it. The model must be correctly specified enough for the null simulation to represent the data-generating process, the estimates and optimizations must behave appropriately, and the statistic's bootstrap law must track its null law. Constraints, boundary parameters, weak identification, or a misspecified margin model can invalidate a routine plug-in argument. In practice, the p-value also has Monte Carlo error because only $B$ replicates are used.

Worked Example: Reporting a Finite-Simulation p-Value

Suppose a fitted null model has been used to generate $B=499$ replicates, and the observed likelihood-ratio statistic is $T_{\mathrm{obs}}=3.10$. Assume 18 replicate statistics are at least 3.10. Then $R=18$, and

\hat p=\frac{R+1}{B+1} =\frac{19}{500}=0.038.

The estimated tail proportion without the correction would be $18/499\approx0.0361$, but the prescribed Monte Carlo p-value is 0.038. The conclusion should be reported as a plug-in parametric-bootstrap result based on 499 replicates, rather than as an exact p-value for the composite null. More replicates can reduce simulation variability, but they do not fix a poorly specified null model or omitted refitting.

Common pitfall and practical checks

A common mistake is to generate from fitted independent margins and then evaluate every replicate using the original data's fitted margins. This treats estimated quantities as fixed in the replicate statistic, even though their estimation contributed to the observed statistic. The correction is to rerun the complete analysis for each replicate: null margin fit, alternative fit, and statistic calculation. Also ensure that the same optimizer, parameter restrictions, treatment of failed fits, and definition of the statistic are applied to observed and simulated data.

Before interpreting the result, check that the simulated margins are generated from the null fit, not the alternative; that each pair is generated independently under the null; and that the replicate margins are refitted. State $B$, the exceedance count, and the corrected p-value. These details make clear both what was calibrated and the limits of that calibration.

Check Your Understanding

Use these questions to test whether the bootstrap procedure is clear.

Under a copula independence null with unknown marginal parameters, which parameters are estimated from the observed data to generate bootstrap samples?
Why must the marginal parameters be refitted in each bootstrap replicate when the observed statistic uses fitted margins?
If $B=999$ and 24 replicate statistics are at least as large as the observed statistic, what is the corrected Monte Carlo p-value?
Why does the $+1$ p-value formula not make a plug-in parametric bootstrap exact for unknown null parameters?
State one condition sufficient for the asymptotic bootstrap tail probability to be uniform, and explain why it may fail in a misspecified model.