Inference for Indirect Inference Estimators

Why Inference Requires Both Auxiliary Estimates

Indirect inference estimates a structural parameter by matching an auxiliary model’s estimate from observed data to the corresponding estimate from data simulated under candidate structural parameters. The auxiliary model is estimated rather than known, so its sampling variation contributes to uncertainty in the structural estimate. The simulated auxiliary estimate can also be random, and it can be statistically dependent on the data-side estimate under some simulation designs.

Inference therefore starts with the difference between the two auxiliary estimates, not with either estimate in isolation. The key result below gives the first-order effect of that difference on the structural parameter estimator. It also makes an important distinction: a covariance formula does not, by itself, establish that the limiting distribution is normal.

Definition: Let $\theta_0$ be the true $p$-dimensional structural parameter, and let $\hat b_d$ be a $q$-dimensional auxiliary estimate from the observed data. Let $\hat b_s(\theta)$ be the auxiliary estimate from data simulated at candidate value $\theta$. The binding function is $b(\theta)$, the probability limit of $\hat b_s(\theta)$. An indirect inference estimator minimizes $$ [\hat b_d-\hat b_s(\theta)]'W_n[\hat b_d-\hat b_s(\theta)] $$ over $\theta$, where $W_n$ is a weighting matrix. The matched population value satisfies $b(\theta_0)=b_0$, the probability limit of $\hat b_d$.

Usually $q\geq p$: the auxiliary model can provide more matching conditions than there are structural parameters. Define the $q\times p$ derivative matrix $D=\partial b(\theta_0)/\partial\theta'$. Local identification requires $D$ to have full column rank. We take $W_n$ to converge to a symmetric positive-definite matrix $W$. These conditions make $A=D'WD$ nonsingular.

First-Order Distribution and Covariance

The result uses a joint limit for the two auxiliary estimates. Assume that, as the sample size $n$ grows,

\begin{pmatrix} \sqrt n(\hat b_d-b_0)\\ \sqrt n(\hat b_s(\theta_0)-b_0) \end{pmatrix} \xrightarrow{d} \begin{pmatrix} Z_d\\ Z_s \end{pmatrix}.

Assume the limit has finite second moments and that both components have mean zero. Write $V_d=\operatorname{Var}(Z_d)$, $V_s=\operatorname{Var}(Z_s)$, and $C=\operatorname{Cov}(Z_d,Z_s)$. The covariance of the auxiliary-estimate difference is

\Omega =\operatorname{Var}(Z_d-Z_s) =V_d+V_s-C-C'.

The cross-covariance terms matter: omitting them is correct only when the two limiting errors are uncorrelated. For the expansion below, also assume that the simulated auxiliary estimate has a locally stable first-order behavior around $\theta_0$, so its local change is given by the derivative $D$, and that the minimizer is consistent and satisfies the usual local minimization conditions.

Theorem: Under these conditions, with $A=D'WD$, $$ \sqrt n(\hat\theta-\theta_0) = A^{-1}D'W(Z_d-Z_s)+o_p(1). $$ Its limiting covariance matrix is $$ \operatorname{Avar}\!\left(\sqrt n(\hat\theta-\theta_0)\right) =A^{-1}D'W\Omega WD A^{-1}. $$ The limiting distribution is the distribution of $A^{-1}D'W(Z_d-Z_s)$. It is centered because $Z_d-Z_s$ has mean zero, but it is not necessarily normal. If the joint limit $(Z_d,Z_s)$ is mean-zero Gaussian, then the estimator limit is Gaussian with the displayed covariance.

Here $\operatorname{Avar}$ denotes the covariance of the limiting root-$n$ distribution. Accordingly, the approximate covariance matrix of $\hat\theta$ itself is this matrix divided by $n$. Normal Wald intervals require the additional normal-limit condition (or an appropriate normal approximation); the covariance formula alone is not enough.

Deriving the Expansion

Let $r_n(\theta)=\hat b_d-\hat b_s(\theta)$. In a neighborhood of $\theta_0$, the local approximation gives $\hat b_s(\theta)\approx \hat b_s(\theta_0)+D(\theta-\theta_0)$. Thus, writing $h=\sqrt n(\theta-\theta_0)$, the scaled discrepancy is approximately $\sqrt n r_n(\theta_0)-Dh$. The limiting random part at $\theta_0$ is $Z_d-Z_s$.

The limiting quadratic objective for $h$ is therefore $[Z_d-Z_s-Dh]'W[Z_d-Z_s-Dh]$. Differentiating with respect to $h$ and setting the derivative to zero gives the first-order condition below. Full column rank of $D$ and positive definiteness of $W$ imply that $A$ is invertible, so this condition has a unique solution.

-2D'W(Z_d-Z_s-Dh)=0, \qquad h=(D'WD)^{-1}D'W(Z_d-Z_s).

The covariance formula follows by applying the covariance rule for a linear transformation to this expression for $h$. In particular, the middle covariance is $\operatorname{Var}(Z_d-Z_s)=\Omega$, which produces $A^{-1}D'W\Omega WD A^{-1}$. This derivation does not turn an arbitrary limiting distribution into a normal one; linear transformations preserve normality only when the relevant input combination is normal.

Worked Example: A Scalar Wald Interval

Suppose a hypothetical scalar model has $D=2$ and $W=1$. Assume the centered joint limit of the data and simulated auxiliary errors is Gaussian, with $\operatorname{Var}(Z_d)=0.36$, $\operatorname{Var}(Z_s)=0.16$, and $\operatorname{Cov}(Z_d,Z_s)=0.06$. For a sample of $n=100$, suppose the estimate is $\hat\theta=1.20$. Calculate an approximate 95% Wald interval.

First include the dependence between the two auxiliary errors: $$ \Omega=0.36+0.16-0.06-0.06=0.40. $$ Here $A=D'WD=2(1)(2)=4$. The variance of the root-$n$ estimator is $$ A^{-1}D'W\Omega WD A^{-1} =\frac{1}{4}(2)(1)(0.40)(1)(2)\frac{1}{4} =0.10. $$

The estimated standard error for $\hat\theta$ is therefore $$ SE(\hat\theta)\approx\sqrt{\frac{0.10}{100}}=0.03162. $$ Using the standard normal 0.975 quantile, 1.96, gives $$ 1.20\pm1.96(0.03162) =1.20\pm0.06198 =[1.138,\ 1.262]. $$

Under the stated Gaussian-limit assumption, this is an approximate 95% confidence interval for the structural parameter, in the same units as $\theta$. Its repeated-sampling coverage is approximately 95% under the model and regularity conditions; it does not mean there is a 95% probability that the fixed parameter lies in this realized interval.

Estimating the Covariance in Practice

In an application, $D$, $W$, and the covariance components are unknown. They are replaced by consistent estimates. Let $\hat D$ estimate the binding-function derivative, $\hat W$ estimate the limiting weight, and let $\hat V_d,\hat V_s,\hat C$ estimate the relevant covariance matrices. Set

\hat\Omega=\hat V_d+\hat V_s-\hat C-\hat C', \qquad \hat A=\hat D'\hat W\hat D.

The plug-in estimate of the root-$n$ covariance is $\hat A^{-1}\hat D'\hat W\hat\Omega\hat W\hat D\hat A^{-1}$; divide it by $n$ to estimate the covariance of $\hat\theta$. The derivative can be estimated analytically or by numerical differentiation of the binding function. The covariance estimates must correspond to the actual estimation and simulation design. In particular, treating the data-side auxiliary estimate as fixed discards its estimation variance and generally understates uncertainty.

Worked Example: Two Auxiliary Statistics

Consider a hypothetical model with two auxiliary statistics and one structural parameter. Let $D=(1,2)'$ and $W=I_2$, where $I_2$ is the two-dimensional identity matrix. Suppose the estimated limiting covariance components are $$ V_d= \begin{pmatrix}0.8&0.1\\0.1&0.4\end{pmatrix}, \quad V_s= \begin{pmatrix}0.3&0.05\\0.05&0.2\end{pmatrix}, \quad C= \begin{pmatrix}0.1&0\\0.02&0.05\end{pmatrix}. $$ Assume the joint limit is Gaussian, $n=400$, and the estimate is $\hat\theta=0.70$.

The difference covariance is $$ \Omega=V_d+V_s-C-C' =\begin{pmatrix}0.9&0.13\\0.13&0.5\end{pmatrix}. $$ Since $A=D'D=1^2+2^2=5$, the root-$n$ variance is $$ A^{-1}D'\Omega D A^{-1} =\frac{(1,2) \begin{pmatrix}0.9&0.13\\0.13&0.5\end{pmatrix} \binom{1}{2}}{25} =\frac{0.9+0.52+2.00}{25} =0.1368. $$

Thus $SE(\hat\theta)\approx\sqrt{0.1368/400}=0.01849$, and the approximate 95% Wald interval is $0.70\pm1.96(0.01849)=[0.6638,\ 0.7362]$. The off-diagonal terms contribute to the calculation: both covariance among the auxiliary statistics and cross-covariance between their data and simulation estimates are part of the uncertainty calculation.

Bootstrap Confidence Intervals

A bootstrap can approximate the sampling distribution without relying directly on a plug-in covariance formula. For a parametric bootstrap, generate a bootstrap data set under the fitted structural model at $\hat\theta$. In each replication, estimate the auxiliary model on that bootstrap data set, repeat the simulation and estimation steps used by the indirect inference procedure, and re-estimate the structural parameter. Denote the resulting estimates by $\hat\theta^{*(1)},\ldots,\hat\theta^{*(B)}$.

The auxiliary model must be re-estimated in each replication. Holding the observed-data auxiliary estimate fixed would fail to reproduce its sampling variation, and would not generally approximate the sampling distribution of the original estimator. The bootstrap should also reproduce the estimator’s weighting and simulation protocol closely enough for its conditional distribution to represent the intended procedure.

A percentile interval uses empirical quantiles of the bootstrap estimates. A basic interval reflects the bootstrap errors around $\hat\theta$ back across the estimate. A bootstrap standard-error interval instead uses the standard deviation of the bootstrap estimates, usually with an appropriate critical value. These methods have different finite-sample behavior; their validity depends on regularity and on the bootstrap consistently approximating the estimator’s sampling distribution.

Worked Example: Percentile and Basic Bootstrap Intervals

Suppose a hypothetical parametric bootstrap with $B=2{,}000$ replications gives $\hat\theta=0.70$, a 2.5th percentile of $0.662$, and a 97.5th percentile of $0.741$. These are stated as empirical quantiles from the full set of bootstrap estimates.

The 95% percentile interval is directly $$ [0.662,\ 0.741]. $$ For the basic interval, reflect the bootstrap quantiles around the original estimate: $$ [2(0.70)-0.741,\ 2(0.70)-0.662] =[0.659,\ 0.738]. $$

Both intervals are on the structural parameter scale. They differ because they use the bootstrap distribution in different ways; neither is guaranteed to be better in every setting. Their interpretation remains a repeated-sampling coverage statement under the conditions supporting the bootstrap, not a probability statement about the fixed parameter after observing the interval.

When the Normal Approximation Is Not Justified

A covariance calculation summarizes second moments; it does not identify the shape of a distribution. For instance, suppose in a scalar setting $D=W=1$ and $Z_d-Z_s$ has a mean-zero, variance-one nonnormal limiting distribution. The expansion then gives $\sqrt n(\hat\theta-\theta_0)$ the same nonnormal limit, with variance one. Reporting that limit as $N(0,1)$ would be unjustified. A nonzero mean in the joint limit would also shift the estimator’s limiting distribution, so a centered confidence interval would require further analysis.

For conventional analytical Wald intervals, a joint central limit result for the auxiliary estimates is a common sufficient condition for the Gaussian limit. If only a general joint weak limit is available, the linear-transform result still describes the estimator’s limit, but normal critical values need not give the claimed coverage. Bootstrap inference can be useful when its approximation is valid, but it is not an automatic remedy for nonregularity or a misspecified bootstrap design.

Common Mistakes and Practical Checks

Calling every limit normal. A covariance formula gives second moments, not distributional shape. Establish a Gaussian joint limit, or report the linear-transform limiting distribution without labeling it normal.
Ignoring estimated auxiliary parameters. The data-side auxiliary estimate is random. Its variability belongs in $V_d$, and dependence with simulation error belongs in $C$.
Dropping the cross-covariance without justification. The difference variance is $V_d+V_s-C-C'$, not automatically $V_d+V_s$. The latter follows only if the limiting errors are uncorrelated.
Bootstrapping while fixing the auxiliary estimate. Re-estimate the auxiliary model in each replication to reproduce its estimation uncertainty.
Mixing up root-$n$ covariance and estimator covariance. The covariance formula is for $\sqrt n(\hat\theta-\theta_0)$; divide by $n$ before calculating the standard error of $\hat\theta$.

Key takeaway: Indirect inference uncertainty is the transformed uncertainty in the difference between the data and simulated auxiliary estimates. Estimate that difference’s covariance with its dependence structure intact, and use normal Wald inference only when a normal limiting distribution is justified.

Check Your Understanding

Work through these questions using the definitions and assumptions above.

Why does estimating the auxiliary model on observed data contribute to the variance of an indirect inference estimator?
For scalar auxiliary estimates with variances $0.5$ and $0.2$, and covariance $0.1$, calculate the variance of their difference.
What additional condition justifies describing the linear-transform limit as normal?
In a parametric bootstrap for indirect inference, which quantities or models need to be re-estimated within each replication, and why?
Explain why the covariance of $\sqrt n(\hat\theta-\theta_0)$ must be divided by $n$ to obtain an approximate covariance matrix for $\hat\theta$, and how to obtain a scalar standard error from the estimated variance.