Inference with Estimated Choice Probabilities

Why First-Stage Uncertainty Matters

Two-step conditional choice probability (CCP) estimation replaces some of the calculations in a dynamic choice model with estimates of choice probabilities from observed choices. The second stage then uses those estimated probabilities to recover structural parameters. This can make estimation computationally attractive, but it creates an inferential requirement: the probabilities supplied to the second stage are estimates, not known inputs.

If second-stage uncertainty is calculated as though the estimated CCPs were fixed, the calculation generally omits a source of sampling variation. The omission can affect standard errors, confidence intervals, and tests of structural parameters. Moreover, the first- and second-stage estimates usually come from the same observations, so their errors may be correlated. Correct inference must account for both the uncertainty passed from stage one and any covariance between stages.

Definition: Let \(\eta\) denote a vector of CCPs or other first-stage quantities, and let \(\theta\) denote the second-stage structural parameter. The two-step estimator is \((\hat{\eta},\hat{\theta})\), where \(\hat{\eta}\) is estimated from the data and \(\hat{\theta}\) solves second-stage conditions that depend on \(\hat{\eta}\). Valid second-stage inference accounts for the sampling variation in both estimates.

The general results below are stated for independent sampling units, such as independent individuals whose choice histories form the observations. Choices within one individual's history need not be independent. If sampling is clustered, the variance calculation must treat the independent cluster, rather than an individual period, as the sampling unit.

Two-Stage Estimating Equations

Write the first-stage estimator as a solution to sample moment conditions and the second stage as a solution to another set of moment conditions:

\frac{1}{N}\sum_{i=1}^{N} h_i(\hat{\eta})=0, \qquad \frac{1}{N}\sum_{i=1}^{N} m_i(\hat{\theta},\hat{\eta})=0.

Here \(h_i\) may encode the choice frequencies used to estimate CCPs. The vector \(m_i\) is a second-stage score or moment condition. At the true values, the population expectations satisfy \(E[h_i(\eta_0)]=0\) and \(E[m_i(\theta_0,\eta_0)]=0\). In applications with weights or conditional moments, the equations should reflect the estimator actually used.

For example, for a binary choice \(A_i\) observed in state \(S_i=s\), the sample CCP \(\hat{p}\) solves

\sum_{i=1}^{N} \mathbf{1}\{S_i=s\}( \mathbf{1}\{A_i=a\}-\hat{p})=0.

This equation makes clear that the first-stage probability is estimated from the observations in the relevant state. When the second stage uses \(\hat{p}\), variation in the state-specific choice frequency can therefore change \(\hat{\theta}\).

Worked Example: A CCP as an Input to a Second-Stage Parameter

Suppose a hypothetical model has one binary choice probability \(p\), and the second-stage parameter is \(\theta=\log\{p/(1-p)\}\). In 200 observed visits to a particular state, 80 choices are of action \(a\). Treat the 200 visits as the sampling units for this illustration.

The first-stage estimate is \(\hat{p}=80/200=0.4000\). The second-stage estimate is \[ \hat{\theta}=\log\left(\frac{0.4000}{0.6000}\right)=-0.4055. \] The second stage is a transformation of the estimated CCP, so its uncertainty cannot be zero merely because the transformation is calculated exactly.

For independent binary choices in this state, the estimated variance of \(\hat{p}\) is \[ \widehat{\operatorname{Var}}(\hat{p}) =\frac{\hat{p}(1-\hat{p})}{200} =\frac{0.4000(0.6000)}{200}=0.001200, \] and its standard error is \(\sqrt{0.001200}=0.03464\). The derivative of the logit transformation is \(1/\{p(1-p)\}\), which at \(\hat{p}=0.4000\) equals \(4.1667\). The delta-method standard error for \(\hat{\theta}\) is therefore \[ SE(\hat{\theta})\approx 4.1667(0.03464)=0.1443. \] A 95% normal-approximation interval is \[ -0.4055\ \mathbin{\pm}\ 1.96(0.1443) =[-0.6883,-0.1227]. \] It describes uncertainty about this hypothetical log-odds parameter under the stated sampling and model assumptions; it is not a probability statement about the fixed parameter. This simple transformation illustrates the first-stage contribution before considering more general second-stage moments.

Influence Function and Sandwich Variance

The key calculation is a first-order expansion of the two estimating equations around the true parameter values. Define the population derivative matrices

H=E\!\left[\frac{\partial h_i(\eta_0)}{\partial \eta'}\right], \qquad J=E\!\left[\frac{\partial m_i(\theta_0,\eta_0)}{\partial \theta'}\right], \qquad D=E\!\left[\frac{\partial m_i(\theta_0,\eta_0)}{\partial \eta'}\right].

Assume these derivatives exist, \(H\) and \(J\) are nonsingular, the estimators are consistent, and suitable finite-moment and regularity conditions justify the expansions. The first-stage equation gives

\sqrt{N}(\hat{\eta}-\eta_0) = -H^{-1}\frac{1}{\sqrt{N}}\sum_{i=1}^{N}h_i(\eta_0)+o_p(1).

Expanding the second-stage equation gives

0= \frac{1}{\sqrt{N}}\sum_{i=1}^{N}m_i(\theta_0,\eta_0) +J\sqrt{N}(\hat{\theta}-\theta_0) +D\sqrt{N}(\hat{\eta}-\eta_0)+o_p(1).

Substitute the first-stage expansion and solve for the second-stage term. This yields

\sqrt{N}(\hat{\theta}-\theta_0) = -J^{-1}\frac{1}{\sqrt{N}}\sum_{i=1}^{N} \left[m_i(\theta_0,\eta_0)-D H^{-1}h_i(\eta_0)\right]+o_p(1).

The expression in square brackets is the adjusted second-stage influence contribution, apart from the premultiplying matrix \(-J^{-1}\). It combines the direct second-stage moment with the effect of first-stage estimation error. In particular, the contribution does not generally have the variance of \(m_i\) alone: the two terms can covary because they are computed from the same sampling unit.

Theorem: Under the regularity conditions above and independent sampling units, the asymptotic covariance matrix of \(\sqrt{N}(\hat{\theta}-\theta_0)\) is \[ J^{-1}\Omega J^{-1'}, \qquad \Omega=\operatorname{Var}\!\left(m_i-DH^{-1}h_i\right). \] A consistent estimated covariance matrix for \(\hat{\theta}\) is \(N^{-1}\hat{J}^{-1}\hat{\Omega}\hat{J}^{-1'}\), using estimated derivatives and the empirical covariance of the adjusted contributions. For clustered observations, form and take the covariance of cluster-level sums of the contributions instead.

The same result can be obtained from a stacked system. The first-stage equation has no \(\theta\) term, while the second-stage equation depends on both parameters. The derivative matrix therefore has a block-triangular form:

A= \begin{pmatrix} H&0\\ D&J \end{pmatrix}.

Stack the moment vector as \((h_i',m_i')'\) and estimate its joint covariance, including the off-diagonal covariance between \(h_i\) and \(m_i\). The usual sandwich calculation for this joint system produces the same second-stage variance as the adjusted influence function. This is often the most direct implementation when the first- and second-stage equations can be written explicitly.

Worked Variance Calculation

Worked Example: Propagating CCP Error and Cross-Stage Covariance

Consider a hypothetical second-stage estimator defined by \(\hat{\theta}=\bar{Y}-2\hat{p}\), where \(\hat{p}=\bar{Z}\) and \(Z\) is a binary indicator of a particular choice. Thus the second-stage moment is \(m_i(\theta,p)=\theta-Y_i+2p\). A hypothetical sample of \(N=1000\) independent decision-makers has \(\bar{Y}=3.200\), \(\hat{p}=0.4000\), \(\widehat{\operatorname{Var}}(Y)=1.000\), \(\widehat{\operatorname{Var}}(Z)=0.2400\), and \(\widehat{\operatorname{Cov}}(Y,Z)=0.1200\). These empirical moments are coherent: the binary indicator mean is 0.4000, and the stated covariance is within the feasible range.

The estimate is \[ \hat{\theta}=3.200-2(0.4000)=2.400. \] Since \(\hat{\theta}=\bar{Y}-2\bar{Z}\), its estimated variance is \[ \widehat{\operatorname{Var}}(\hat{\theta}) =\frac{1}{1000}\left[ \widehat{\operatorname{Var}}(Y) +4\widehat{\operatorname{Var}}(Z) -4\widehat{\operatorname{Cov}}(Y,Z)\right]. \] Substituting the reported moments gives \[ \widehat{\operatorname{Var}}(\hat{\theta}) =\frac{1.000+4(0.2400)-4(0.1200)}{1000} =\frac{1.480}{1000}=0.001480. \] Therefore \(SE(\hat{\theta})=\sqrt{0.001480}=0.03847\), and a 95% normal-approximation interval is \[ 2.400\ \mathbin{\pm}\ 1.96(0.03847) =[2.325,2.475]. \] If the analyst incorrectly held \(p\) fixed, the variance calculation would retain only \(1.000/1000\), giving a standard error of \(0.03162\). That is not the correct standard error here. Notice also that simply adding the variance terms while dropping the covariance would give \(1.960/1000\), not \(1.480/1000\). The negative covariance term in the variance of \(\bar{Y}-2\bar{Z}\) matters; first-stage uncertainty does not always increase the variance relative to a particular fixed-input calculation.

Full-Sample Bootstrap

A full-sample bootstrap provides an alternative way to carry first-stage uncertainty into inference. In each bootstrap replication, resample the independent sampling units, re-estimate the CCPs, and rerun the entire second stage using those re-estimated probabilities. The resulting distribution of second-stage estimates reflects the joint sampling variation of the two stages, including their dependence.

The resampling unit must match the data structure. If observations are individuals with multiple choice periods, resample individuals and retain each selected individual's full history. Resampling individual choice occasions as though they were independent would generally destroy within-person dependence and could give invalid standard errors. If the estimator uses estimated state transitions or other first-stage quantities, re-estimate those in each replication as well.

Choose the sampling unit.
Identify the independent unit implied by the data and sampling design, such as a complete individual history or a cluster.

Resample units.
Draw a bootstrap sample with replacement, preserving each selected unit's relevant observations and history.

Repeat both estimation stages.
Re-estimate the CCPs from the bootstrap sample, then re-estimate the structural parameters using those probabilities.

Summarize the estimates.
Use the variation across the replicated second-stage estimates to estimate uncertainty; use enough replications for stable results and check for failed or boundary estimates.

Worked Example: Bootstrap Standard Error From Replicated Estimates

Suppose a hypothetical full-sample bootstrap is run by resampling complete individual histories and re-estimating both stages each time. Five bootstrap replications, shown here only to make the arithmetic transparent, produce second-stage estimates \(1.82,\ 1.91,\ 2.03,\ 2.08,\ 2.16\). In an actual analysis, substantially more replications would ordinarily be used.

Their mean is \[ \bar{\theta}^{*}=\frac{1.82+1.91+2.03+2.08+2.16}{5}=2.000. \] The bootstrap variance estimate is the sample variance of the replicated estimates: \[ \widehat{\operatorname{Var}}_{\mathrm{boot}}(\hat{\theta}) =\frac{(1.82-2.000)^2+(1.91-2.000)^2+(2.03-2.000)^2+(2.08-2.000)^2+(2.16-2.000)^2}{5-1}. \] The squared deviations sum to \(0.0324+0.0081+0.0009+0.0064+0.0256=0.0734\), so the estimated variance is \(0.0734/4=0.01835\) and the bootstrap standard error is \(\sqrt{0.01835}=0.1355\). This calculation illustrates the standard-error summary; five replications are far too few for reliable final inference. Crucially, each listed estimate must come from repeating both stages, not from holding the original CCP estimates fixed.

Assumptions, Interpretation, and Common Pitfalls

The stacked-equation and influence-function results rely on identification of the relevant CCPs and structural parameters, adequate support for the states used in estimation, and smooth, locally well-behaved estimating equations. In particular, probabilities used in log or logit transformations must be in the interior of their allowable range, and the derivative matrices must be nonsingular. If a state is rarely visited, its estimated CCP can be imprecise; if it is never observed, its CCP is not identified by that sample frequency. A small-sample normal approximation can also be poor when probabilities are near zero or one or when the second-stage objective has weak curvature.

A common mistake is to compute the second-stage Hessian or score-based variance while treating the first-stage CCPs as known. The reasoning is that the second-stage optimization is performed exactly once the estimated probabilities are supplied. The correction is that exact optimization does not remove sampling error from those inputs: \(\hat{\eta}\) changes across samples, and so does \(\hat{\theta}\). Use a stacked sandwich, an equivalent adjusted influence function, or a full-sample bootstrap that re-estimates both stages.

A second pitfall is to assume the two stages are independent just because they are described separately. When both use the same individuals, first-stage and second-stage moment contributions can be correlated. A joint sandwich keeps their covariance; a full-sample bootstrap preserves it through joint resampling. If separate independent samples truly are used for the two stages, the covariance structure differs, but first-stage uncertainty still needs to be propagated.

Key takeaway: A second-stage estimate based on estimated CCPs inherits their sampling error. Valid uncertainty quantification must carry that error through the second-stage mapping and account for dependence between the stages.

Check Your Understanding

Use these questions to check whether the inferential logic is clear.

In a two-stage CCP estimator, why is it generally incorrect to treat the estimated first-stage probabilities as fixed when calculating a second-stage standard error?
In the adjusted influence contribution \(m_i-DH^{-1}h_i\), what role does each of the two terms play?
Why must a stacked sandwich estimate the covariance between the first- and second-stage moments when the same sampling units are used in both stages?
If each individual contributes a multi-period choice history, what should typically be resampled in a full-sample bootstrap, and why?
What can go wrong when a CCP used in a log transformation is estimated at zero?