{'fraction_non_alphanumeric': 0.044704276451676635, 'fraction_numerical': 0.016003017403954954, 'mean_word_length': 4.3198929868144464, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Estimation in a Generalization of Bivariate Probit Models with Dummy Endogenous Regressors

The purpose of this paper is to provide guidelines for empirical researchers who use a class of bivariate threshold crossing models with dummy endogenous variables. A common practice employed by the researchers is the specification of the joint distribution of the unobservables as a bivariate normal distribution, which results in a bivariate probit model. To address the problem of misspecification in this practice, we propose an easy-to-implement semiparametric estimation framework with parametric copula and nonparametric marginal distributions. We establish asymptotic theory, including root-n normality, for the sieve maximum likelihood estimators that can be used to conduct inference on the individual structural parameters and the average treatment effect (ATE). In order to show the practical relevance of the proposed framework, we conduct a sensitivity analysis via extensive Monte Carlo simulation exercises. The results suggest that the estimates of the parameters, especially the ATE, are sensitive to parametric specification, while semiparametric estimation exhibits robustness to underlying data generating processes. We then provide an empirical illustration where we estimate the effect of health insurance on doctor visits. In this paper, we also show that the absence of excluded instruments may result in identification failure, in contrast to what some practitioners believe.


Introduction
The purpose of this paper is to provide guidelines for empirical researchers who use a class of bivariate threshold crossing models with dummy endogenous variables. This class of models is typically written as follows. With the binary outcome Y and the observed binary endogenous treatment D, we consider where X denotes a vector of exogenous regressors that determine both Y and D, and Z denotes a vector of exogenous regressors that directly affect D, but not Y (i.e., instruments for D). Since Y does not appear in the equation for D, this model forms a triangular model, as a special case of a simultaneous equations model, with the binary endogenous variables. In this paper, we investigate the consequences of the common practices employed by empirical researchers who use this class of models. As an important part of this investigation, we conduct a sensitivity analysis on the specification of the joint distribution of the unobservables (ε, ν). This is the component of the model that practitioners have the least knowledge about, and thus typically impose a parametric assumption. To address the problem of misspecification, we propose a semiparametric estimation framework with parametric copula and nonparametric marginal distributions. The semiparametric specification is an attempt to ensure robustness while achieving point identification and efficient estimation.
The parametric class of models (1.1) includes the bivariate probit model, in which the joint distribution of (ε, ν) is assumed to be a bivariate normal distribution. This model has been widely used in empirical research, including the works of Evans and Schwab (1995), Neal (1997), Goldman et al. (2001), Altonji et al. (2005), Bhattacharya et al. (2006), Rhine et al. (2006) and Marra and Radice (2011) to name a just few. The distributional assumption in this model, however, is made out of convenience or convention, and is hardly justified by underlying economic theory and thus susceptible to misspecification. With binary endogenous regressors, the objects of interest in model (1.1) are the mean treatment parameters, in addition to the individual structural parameters. Because the outcome variable is also binary, the mean treatment parameters such as the average treatment effect (ATE) are expressed as the differential between the marginal distributions of ε. Therefore, the problem of misspecification when estimating these treatment parameters can be even more severe than that when estimating individual parameters.
To one extreme, a nonparametric joint distribution of (ε, ν) can be used in a bivariate threshold crossing model, as in Shaikh and Vytlacil (2011). Their results, however, suggest that the ATE is only partially identified in this fully flexible setting. Instead of sacrificing point identification, we impose a parametric assumption on the dependence structure between the unobservables using copula functions that are known up to a scalar parameter. At the same time, in order to ensure robustness, we allow the marginal distribution of ε (and ν), which is involved in the calculation of the ATE, to be unspecified. Our class of models encompasses both parametric and semiparametric models with parametric copula and either parametric or nonparametric marginal distributions.
This broad range of models allows us to conduct a sensitivity analysis on the specification of the joint distribution of (ε, ν).
The identification of the individual parameters and the ATE in this class of models is established in Han and Vytlacil (2017, hereafter, HV17). They show that when the copula function for (ε, ν) satisfies a certain stochastic ordering, identification is achieved in both parametric and semiparametric models under an exclusion restriction and mild support conditions. Building on these results, we consider estimation and inference in the same setting. For the semiparametric class of models (1.1) with parametric copula and nonparametric marginal distributions, the likelihood contains infinite-dimensional parameters (i.e., the unknown marginal distributions).
To estimate this model, we consider the sieve maximum likelihood (ML) estimation method for the finite-and infinite-dimensional parameters of the model, as well as their functionals. The estimation of the parametric model, on the other hand, is within the standard ML framework.
The contributions of this paper can be summarized as follows. Through these contributions, this paper is intended to provide a guideline to empirical researchers. First, we establish the asymptotic theory for the sieve ML estimators in a class of semiparametric copula-based models. This result can be used to conduct inference on the functionals of the finite-and infinitedimensional parameters, such as inference on the individual structural parameters and the ATE.
We show that the sieve ML estimators are consistent and that their smooth functionals are root-n asymptotically normal.
Second, in order to show the practical relevance of the theoretical results for empirical researchers, we conduct a sensitivity analysis via extensive Monte Carlo simulation exercises. We find that the parametric ML estimates, especially those for the ATE, can be highly sensitive to the misspecification of the marginal distributions of the unobservables. On the other hand, the sieve ML estimates perform well in terms of the mean squared error (MSE) as they are robust to the underlying data generating process. Moreover, their performance is comparable to that of the parametric estimates under a correct specification. We also show that copula misspecification does not have a substantial effect in estimation, as long as the true copula is within the stochastic ordering class of the identification. As copula misspecification is a problem common to both parametric and semiparametric models considered in this paper, our sensitivity analysis suggests that a semiparametric consideration may be more preferable in estimation and inference.
Third, we formally show that identification may fail without the exclusion restriction, in contrast to the findings of Wilde (2000). The bivariate probit model is sometimes used in applied work without instruments (e.g., White and Wolaver (2003) and Rhine et al. (2006)). We show, however, that this restriction is not only sufficient but also necessary for identification in paramet-ric and semiparametric models when there is a single binary exogenous variable common to both equations. We also show that under joint normality of the unobservables, the parameters are, at best, weakly identified when there are common (and possibly continuous) exogenous variables. 1 We also note that another source of identification failure is the absence of restrictions on the dependence structure of the unobservables, as mentioned above.
The sieve estimation method is a useful nonparametric estimation framework that allows for a flexible specification, while guaranteeing the tractability of the estimation problem; see Chen (2007) for a survey of sieve estimation in semi-nonparametric models. The estimation method is also easy to implement in practice. The sieve ML estimation has been used in various contexts: Chen et al. (2006, hereafter, CFT06) consider the sieve estimation of semiparametric multivariate distributions that are modeled using parametric copulas; Bierens (2008) applies the estimation method to the mixed proportional hazard model; and Hu and Schennach (2008) and Chen et al. (2009) use the method to estimate nonparametric models with non-classical measurement errors.
The asymptotic theory developed in this paper is based on the results established in the sieve extremum estimation literature (e.g., CFT06; Chen (2007); Bierens (2014)). A semiparametric version of bivariate threshold crossing models is also considered in Marra and Radice (2011) and Ieva et al. (2014). In contrast to our setting, however, they introduce flexibility for the index function of the threshold, and not for the distribution of the unobservables.
The remainder of this paper is organized as follows. The next section reviews the identification results of HV17, and then discusses the lack of identification in the absence of exclusion restrictions and in the absence of restrictions on the dependence structure of the unobservables. Section 3 introduces the sieve ML estimation framework for the semiparametric class of models defined in (1.1), and Section 4 establishes the large sample theory for the sieve ML estimators. The sensitivity analysis is conducted in Section 5 by investigating the finite sample performance of the parametric ML and sieve ML estimates under various specifications. Section 6 concludes the paper.
2 Identification and Failure of Identification
1 HV17 only show the sufficiency of this restriction for identification. Mourifié and Méango (2014) show the necessity of the restriction, but their argument does not exploit all information available in the model; see Section 2.2 of the present paper for further details. Assumption 1 X and Z satisfy that (X, Z) ⊥ (ε, ν), where "⊥" denotes statistical independence.
Assumption 2 (X , Z ) does not lie in a proper linear subspace of R k+l a.s. 2 Assumption 3 There exists a copula function C : (0, 1) 2 → (0, 1) such that the joint distribution where F ε and F ν are the marginal distributions of ε and ν, respectively, that are strictly increasing and absolutely continuous with respect to Lebesgue measure. 3 Assumption 4 As scale and location normalizations, α 1 = β 1 = 1 and α 0 = β 0 = 0.
A model with alternative scale and location normalizations, V ar(ε) = V ar(ν) = 1 and E[ε] = E[ν] = 0, can be viewed as a reparametrized version of the model with the normalizations given in Assumption 4; see, for example, the reparametrization (2.1) below. For x ∈ supp(X) and z ∈ supp(Z), write a one-to-one map (by Assumption 3) as (2.1) Take (x, z) and (x,z), for some x ∈ supp(X|Z = z) ∩ supp(X|Z =z), where supp(X|Z) is the conditional support of X, given Z. Then, by Assumption 1, model (1.1) implies that the fitted probabilities are written as p 11,xz = C(r 1,x , s xz ), p 11,xz = C(r 1,x , s xz ), p 10,xz = r 0,x − C(r 0,x , s xz ), as the basis for the identification and estimation of the model. Depending upon whether one is willing to impose an additional assumption on the dependence structure of the unobservables (ε, ν) via C(·, ·), the underlying parameters of the model are either point identified or partially identified.
We first consider point identification. The results for point identification can be found in HV17, which we adapt here given Assumption 4. The additional dependence structure can be characterized in terms of the stochastic ordering of the copula parametrized with a scalar parameter.
Such copulas are referred to as stochastically increasing (SI) or stochastically decreasing (SD), respectively. Then,C is strictly more SI (or less SD) than This ordering is equivalent to having a ranking in terms of the first order stochastic dominance.
Let (U 1 , U 2 ) ∼ C and (Ũ 1 ,Ũ 2 ) ∼C. WhenC is strictly more SI (less SD) than C is, then Pr[Ũ 2 > u 2 |Ũ 1 = u 1 ] increases even more than Pr[U 2 > u 2 |U 1 = u 1 ] does as u 1 increases. 5 Assumption 5 The copula in Assumption 3 satisfies C(·, ·) = C(·, ·; ρ) with a scalar dependence parameter ρ ∈ Ω, is twice differentiable in u 1 , u 2 and ρ, and satisfies The meaning of the last part of this assumption is that the copula is ordered in ρ in the sense of the stochastic ordering defined above. This requirement defines the class of copulas that we allow for identification. Many well-known copulas satisfy (2.3): the normal copula, Plackett copula, Frank copula, Clayton copula, among many others; see HV17 for the full list of copulas and their expressions. Under these assumptions, we first discuss the identification in a fully parametric model.
The proof of this theorem is a minor modification of the proof of Theorem 5.1 in HV17.
Although the parametric structure on the copula is necessary for the point identification of the parameters, HV17 show that the parametric assumption for F ε and F ν are not necessary. In addition, if we make a large support assumption, we can also identify the nonparametric marginal distributions F ε and F ν .
Assumption 7 (i) The distributions of X j (for 1 ≤ j ≤ k) and Z j (for 1 ≤ j ≤ l) are absolutely continuous with respect to Lebesgue measure; (ii) There exists at least one element X j in X such that its support conditional on (X 1 , ..., X j−1 , X j+1 , ..., X k ) is R and α j = 0 and β j = 0, where, without loss of generality, we let j = 1. and F ν (·) are identified up to additive constants.
An interesting function of the underlying parameters that are point identified under the parametric and semiparametric distributional assumptions is the conditional ATE: (2.4)

The Failures of Identification
In this section, we discuss two sources of identification failure in the class of models (1.1): the absence of exclusion restrictions and the absence of restrictions on the dependence structure of the unobservables (ε, ν).

No Exclusion Restrictions
There are empirical works where (1.1) is used without excluded instruments; see, e.g., White and Wolaver (2003) and Rhine et al. (2006). Identification in these papers relies on the results of Wilde (2000), who provides an identification argument by counting the number of equations and unknowns in the system. Here, we show that this argument is insufficient for identification. We show that without the excluded instruments (i.e., when γ = 0), the structural parameters are not identified, even with a full parametric specification of the joint distribution (Assumptions 5 and 6). The existence of common exogenous covariates X in both equations is not very helpful for identification in a sense that becomes clear below.
Before considering the lack of identification in a general case with possibly continuous X 1 in X = (1, X 1 ), we start the analysis with binary X 1 . Mourifié and Méango (2014) show the lack of identification when there is no excluded instrument in a bivariate probit model with binary X 1 . They, however, only provide a numerical counter-example. Moreover, their analysis does not consider the full set of observed fitted probabilities, and hence possibly neglects information that could have contributed to the identification. Here, we provide an analytical counter-example in a more general parametric class of model (1.1) that nests the bivariate probit model. We show that there exist two distinct values of (δ 1 , ρ, µ ε , σ ε ) that generate the same observed fitted probabilities, even if the full set of probabilities are used. Note that the reduced-form parameters (µ ν , σ ν ) are always identified from the equation for D, and α = β = (0, 1) as a normalization using scalar X 1 .
In showing this lack-of-identification result, we find a counter-example where the copula density induced by C(u 1 , u 2 ) is symmetric around u 2 = u 1 and u 2 = 1 − u 1 , and the density induced by F ε is symmetric. Note that the bivariate normal distribution, namely, the normal copula with normal marginals, satisfies these symmetry properties. That is, in the bivariate probit model with a common binary exogenous covariate and no excluded instruments, the structural parameters are not identified.
One might argue that the lack of identification in Theorem 2.4 is due to the limited variation of X. Although this is a plausible conjecture, this does not seem to be the case in the model considered here. 6 We now consider a general case with possibly continuous X 1 , and discuss what can be said about the existence of two distinct sets of (β, δ 1 , ρ, µ ε , σ ε ) that generate the same observed data. To this end, define Then, Similar to the proof strategy for the binary X 1 case, we want to show that, given (α, µ ν , σ ν ), there are two distinct sets of parameter values (β, δ 1 , ρ, µ ε , σ ε ) and (β * , δ * 1 , ρ * , µ * ε , σ * ε ) that generate the same observed fitted probabilities p yd,x for all (y, d) ∈ {0, 1} 2 and x ∈ supp(X) under some choices of C(u 1 , u 2 ) and F ε .
Let t(x) ≡ Fε(x β) ∈ (0, 1) for all x and for some β. Also, choose δ 1 = 0 and some ρ ∈ Ω. For ρ * > ρ, we want to show that there exists (β * , δ * 1 ) such that, for t * (x) ≡ Fε(x β * ), The question is whether we find (β, δ 1 , ρ) and (β * , δ * 1 , ρ * ) such that (2.5)-(2.7) hold simultaneously. First, note that, since ρ * > ρ, we have t * > t and hence β * = β by the assumption that there is no linear subspace in the space of X. Now, choose C(·, ·; ρ) to be a normal copula and choose ρ = 0 and ρ * = 1. Then, using arguments similar to those of the binary case (found in the Appendix), (2.8) and s † (x) = q(x)t(x). Then, (2.7) can be rewritten as The complication here is to ensure that this equation is satisfied for all x. Note that (2.8) and (2.9) are consistent with the definition of a distribution function of a continuous r.v.: Fε(+∞) = 1, Fε(−∞) = 0, and Fε(ε) is strictly increasing. We can then numerically show that a distribution function that is close to a normal distribution satisfies the conditions with a particular choice of (β * , δ * 1 ); see Figure 1. Although no formal derivation of the counterexample is given, this result suggests the following: (i) In the bivariate probit model with continuous common exogenous covariates and no excluded instruments, the parameters will be, at best, weakly identified; (ii) This also implies that, in the semiparametric model considered in Theorem 2.3, the structural parameters and the marginal distributions are not identified without an exclusion restriction, even if X 1 has large support.

No Restrictions on Dependence Structures
When the restriction imposed on C(·, ·) (i.e., Assumption 5) is completely relaxed, the underlying parameters of model (1.1) may fail to be identified, regardless of whether the exclusion restriction holds. That is, a structure describing how the unobservables (ε, ν) are dependent on each other is necessary for identification. This is closely related to the results in the literature that the treatment parameters (which are lower dimensional functions of the individual parameters) in triangular models similar to (1.1) are only partially identified without distributional assumptions; see Bhattacharya et al. (2008), Chiburis (2010), Shaikh and Vytlacil (2011), and Mourifié (2015).
Suppose Assumptions 1-4 hold. Then the model becomes a semiparametric threshold crossing model in that the joint distribution is completely unspecified. Then, as a special case of Shaikh and Vytlacil (2011), one can easily derive bounds for the ATE F ε (x β + δ 1 ) − F ε (x β). The sharpness of these bounds is shown in their paper under a rectangular support assumption for (X, Z), which is, in turn, relaxed in Mourifié (2015). In addition, using Assumption 6, one can also derive bounds for the individual parameters x β and δ 1 , as shown in Chiburis (2010). When there are no excluded instruments in the model, Chiburis (2010) shows that the bounds on the ATE do not improve on the bounds of Manski (1990), whose argument applies to the individual parameters.

Sieve and Parametric ML Estimations
Let ψ ≡ (α , β , δ 1 , γ, ρ) denote the vector of the structural individual parameters. Let f and f ν be the density functions associated with the distribution functions F and F ν , respectively, of the unobservables. Then, (ψ , f , f ν ) is the set of parameters in the semiparametric version of the model. The model becomes fully parametric, once the infinite-dimensional parameters f and f ν are fully characterized by some finite-dimensional parameters, i.e., f (·; η ) and f ν (·; η ν ) for η ∈ R dη and η ν ∈ R dη ν . This yields (ψ , η , η ν ) to be the set of parameters in the parametric version of the model. For either case, the parameter of the model is denoted as θ for convenience.
That is, θ ≡ (ψ , f , f ν ) in the semiparametric model and θ ≡ (ψ , η , η ν ) in the parametric model. For the rest of this paper, we explicitly express θ 0 to be the true parameter value for θ.
This applies to all the other parameter expressions.
LetΨ be the parameter space for ψ. For the parametric model, the spaces for the finitedimensional parameters η and η ν are denoted as H ⊆ R dη and H ν ⊆ R dη ν , respectively. Then, the parameter spaceΘ for θ ≡ (ψ , η , η ν ) becomes a Cartesian product ofΨ, H , and H ν , i.e., Θ ≡Ψ × H × H ν ⊆ R d ψ +dη +dη ν , in the parametric model. 7 For the semiparametric model, we consider the following function spaces as the spaces for f and f ν : where j ∈ { , ν} and F is a space of functions, which we specify later. Then, the parameter spacẽ Θ of θ ≡ (ψ , f , f ν ) can be written asΘ ≡Ψ × F × F ν in the semiparametric model. Note that the function spaces F and F ν contain functions that are nonnegative.
We adopt the ML method to estimate the parameters in the model.
For example, if one imposes Assumption 6, then η = (µ , σ ) and ην = (µν , σν ) . i = 1, 2, ..., n} be the random sample. For both parametric and semiparametric models with corresponding θ, we define the conditional density function of (Y i , D i ) conditional on (X i , Z i ) as where p yd (x, z; θ) abbreviates the right hand side expression that equates p yd,xz in (2.2). Then, the log of density l(θ, w) ≡ log f (y, d|x, z; θ) becomes Consequently, the log-likelihood function can be written Now, the ML estimatorθ n of θ 0 ≡ (ψ 0 , η 0 , η ν0 ) in the parametric model is defined as For the semiparametric model, let F εn and F νn be appropriate sieve spaces for F ε and F ν , respectively, and let f n (·; a n ) and f νn (·; a νn ) be the sieve approximations of f and f ν on their sieve spaces F n and F νn , respectively. Then, we define the sieve ML estimatorθ n of θ 0 ≡ (ψ 0 , f 0 , f ν0 ) in the semiparametric model as follows: whereΘ n ≡Ψ × F n × F νn is the sieve space for θ.
With the parameter spaces F and F ν in (3.1), we are interested in a class of "smooth" univariate square root density functions. Specifically, we assume that √ f and √ f ν belong to the class of p-smooth functions and we restrict our attention to linear sieve spaces for F and F ν . 8 In this case, the choice of sieve spaces for F and F ν depends on the supports of and ν. If the supports are bounded, then one can use the polynomial sieve, trigonometric sieve, or cosine sieve.
When the supports are unbounded, then we can use the Hermite polynomial sieve or the spline wavelet sieve.
Note that the sieve ML estimatorθ n can be equivalently obtained from the following uncon-strained optimization problem: max ψ,an n i=1 l(ψ, a n , W i ) − λ n Pen(a n ) + τ ε 1 −ˆf ε (t, a n )dt + τ ν 1 −ˆf ν (t, a νn )dt , where a n = (a n , a νn ) , l(W i , ψ, a n ) is the log likelihood with sieve approximations f ε (·, a n ) and f ν (·, a νn ); Pen(a n ) is a penalization term that imposes, for example, the properties of the class of functions, F; and the remaining penalization terms impose the properties of a density function.
In this paper, we implicitly assume that the copula function is correctly specified. As mentioned earlier, using a parametric copula may lead to model misspecification. It is well known that when the model is misspecified, the ML estimator converges to a pseudo-true value which minimizes the Kullback-Leibler (KL) divergence (e.g., White (1982)). This result applies to a semiparametric model (e.g., Chen and Fan (2006a) and Chen and Fan (2006b)) as in our semiparametric case. We, however, do not investigate the asymptotic properties of the sieve estimators under copula misspecification, as it is beyond the scope of this paper. Instead, later in simulation, we investigate how the copula misspecification affects the performance of estimators. 9

Asymptotic Theory for Sieve ML Estimators
In this section, we provide the asymptotic theory for the sieve ML estimatorθ n of θ ≡ (ψ , f , f ν ) in the semiparametric model. This theory will be useful for practitioners to conduct inference.
The asymptotic theory for the ML estimatorθ n of θ ≡ (ψ , η , η ν ) in the parametric model is relatively standard and can be found in, e.g., Newey and McFadden (1994). The theory establishes that the parametric ML estimator is consistent, asymptotically normal, and efficient under some regularity conditions. To investigate the asymptotic properties of the sieve ML estimator, we slightly modify our model as follows.
Let G (·) and G ν (·) be strictly increasing functions mapping from R to [0, 1]. We further assume that G and G ν are differentiable and that their derivatives are bounded away from zero on R. Then, we consider the following specification: where H 0 (·) and H ν0 (·) are unknown distribution functions on [0, 1]. It is possible that G (·) and G ν (·) are different, but this is not crucial because the main difficulty with the estimation is in estimating the unknown functions H 0 and H ν0 . Therefore, we assume that G (·) = G ν (·) ≡ G(·) to avoid complexity with the specification and write the derivative of G as g (i.e., g(x ≡ dG(x) dx ). The transformation in equation (4.1) can be found in the literature (e.g., Bierens (2014)), and we do not have any loss of generality. For the known distribution function G, we can choose the standard normal distribution function or the logistic distribution function. Since we assume that the distribution functions of and ν admit density functions, we require that H 0 and H ν0 (·) be differentiable, and write their derivatives as h 0 (·) and h ν0 (·), respectively. For each j ∈ { , ν}, let H j ≡ {h j = g 2 : g ∈ F } for some function space F. With this modification, we redefine the We first verify that there exist H 0 and H ν0 that satisfy (4.1) for given F 0 , F ν0 , and G.
Since G is assumed to be strictly increasing, there exists an inverse function G −1 . Letting , it is straightforward to show that H 0 and H ν0 are mappings from [0, 1] to [0,1] and that satisfy the relations in (4.1). Note too that this transformation does not change the identification results. That is, 1]. Assuming that g is bounded away from zero on R and bounded above, the unknown density function h j0 can be written as 1]. In addition, we can see that h 0 and h ν0 are identified if and only if the unknown marginal density functions f 0 and f ν0 are identified. In what follows, we assume that G is strictly increasing from R to [0, 1] and has a derivative, which is bounded away from zero and bounded above on R.
We redefine the parameter space to facilitate developing the asymptotic theory. The identification requires that the space of the finite-dimensional parameterΨ be open and convex (see Theorems 2.2 and 2.3), and thusΨ cannot be compact. We introduce an "optimization space" that contains the true parameter ψ 0 and consider it as the parameter space of ψ. Formally, we restrict the parameter space for estimation in the following way.

Assumption 8 There exists a compact and convex subset
With the optimization space, we define the parameter space as Θ ≡ Ψ × H × H ν , and the corresponding sieve space is denoted by Θ n ≡ Ψ × H n × H νn . Then, the sieve ML estimator in equation (3.4) is also redefined as follows: θ n ≡ arg max θ∈Θn Q n (θ). (4.2)

Consistency of the Sieve ML Estimators
We begin by showing the consistency of the sieve ML estimator. Since the parameter involves both finite-and infinite-dimensional objects, we establish the consistency of the sieve ML estimators with respect to a pseudo distance function d c on Θ × Θ. 10 All of the norms and the definitions of function spaces in this paper are provided in the Appendix.
We present the following assumptions, under which the sieve ML estimator in equation (4.2) is consistent with respect to the pseudo-metric d c (·, ·).
Assumption 9 There exists a measurable function p(X, Z) such that for all θ ∈ Θ and for all being a Hölder ball with radius R; (iii) the density functions h 0 and h ν0 are bounded away from zero . The derivatives C j (·, ·; ·) and C ρ (·, ·; ·) are uniformly bounded for all j = 1, 2.
Assumption 9 guarantees that the log-likelihood function l(θ, W i ) is well defined for all θ ∈ Θ and that Q 0 (θ 0 ) > −∞. Assumption 10 restricts the data generating process (DGP), and assumes the existence of moments of the data. Assumption 11 defines the parameter space and implies that the infinite-dimensional parameters are in some smooth class called a Hölder class. Note that conditions (i) and (ii) in Assumption 11 together imply that h 0 and h ν0 belong to Λ p R ([0, 1]), whereR ≡ 2 m+1 R 2 < ∞. 11 Thus, we may assume that h 0 and h ν0 belong to a Hölder ball with smoothness p under Assumption 11. 12 The condition that H and H ν are the same can be relaxed, but it is imposed for simplicity. The first part of Assumption 12 restricts our choice of sieve spaces for H and H ν to linear sieve spaces with order k n . This can be relaxed so that the choice of k n is different for h and h ν . The latter part of Assumption 12 requires that the sieve space be chosen appropriately so that the unknown parameters can be well-approximated.
Because the unknown infinite-dimensional parameters belong to a Hölder ball and are defined on bounded supports, we can choose the polynomial sieve, trigonometric sieve, cosine sieve, or spline sieve. 13 For example, if we choose the polynomial sieve or the spline sieve, then one can show Lorentz (1966)). Assumption 13 imposes the boundedness of the derivatives of the copula function.
The following theorem demonstrates that under the above assumptions, the sieve estimator θ n is consistent with respect to the pseudo metric, d c .

Convergence Rates
In this section, we derive the convergence rate of the sieve ML estimator. The convergence rate provides information on how fast the estimator converges to the true parameter value. Heuristically, the faster the convergence rate, the larger the effective sample size is for estimation. The next theorem demonstrates the convergence rate of the sieve ML estimator with respect to the Theorem 4.2 Suppose that Assumptions 1-5 and 7-13 hold. If Assumption 18 in the Appendix additionally holds, then we have The convergence rate given in (4.3) is standard in the literature and slower than the parametric rate (n −1/2 ). In addition, it depends on two components. The first is related to the convergence 11 See the Appendix for details. 12 These conditions implicitly define the strong norm (Hölder norm). 13 Refer to Chen (2007) or CFT06 for details on the choice of sieve spaces.
rate of the "variance term," and this rate increases with the complexity of the sieve space, k n . In contrast, the second component, which reflects the convergence rate of the deterministic approximation error ||θ 0 − π k θ 0 || 2 , decreases as k n increases. The choice of k n ∝ n 1 2p+1 yields the optimal convergence rate, and we can see that with this choice, the convergence rate of the sieve estimatorθ n increases with the degree of smoothness, p.

Asymptotic Normality of Smooth Functionals
We now establish the asymptotic normality of smooth functionals. The parameters in our model contains both finite-and infinite-dimensional parameters, and many objects of interest are written as functionals of both types of the parameters. Therefore, we focus on the asymptotic normality of functionals rather than the parameters themselves. The results of this section can be used to calculate the standard error of the estimate of a functional of interest, or to conduct inference (i.e., testing hypotheses and constructing confidence intervals) based on normal approximation.
Before presenting the results on the asymptotic normality of smooth functionals, we strengthen the smoothness condition in Assumption 5. Let C ij (u 1 , u 2 ; ρ) denote the second-order partial derivative of a copula function C(u 1 , u 2 ; ρ) with respect to i and j, for i, j ∈ {u 1 , u 2 , ρ}.

Assumption 14
The copula function C(u 1 , u 2 ; ρ) is twice continuously differentiable with respect to u 1 , u 2 , and ρ, and its first-and second-order partial derivatives are well defined in a neighborhood of θ 0 .
provided the right hand side limit is well defined. The following assumption characterizes the smoothness of the functional T .
Assumption 15 The following conditions hold: (i) there exist constants w > 1 + 1 2p and a small 0 > 0 such that for any v ∈ V with ||v|| ≤ 0 , Assumption 15 defines a smooth functional T and guarantees the existence of v * ∈V such that Here, we call v * the Riesz representer for the functional T .
The next assumption requires that the Riesz representer be well approximated over the sieve space and that it converges at a rate with respect to the Fisher norm.
The following proposition states that the plug-in sieve ML estimator T (θ n ) of T (θ 0 ) is √ nasymptotically normally distributed under certain conditions. The technical conditions (Assumptions 19 and 20) can be found in the Appendix.

Example 1: Asymptotic Normality for the Finite-Dimensional Parameter ψ 0
The finite-dimensional parameter ψ 0 is a special case of the smooth functionals. Here, we demonstrate the asymptotic normality of the sieve estimator of the finite-dimensional parameter ψ 0 . For any arbitrary λ ∈ R d ψ with |λ| ∈ (0, ∞), let T : Θ → R be a functional of the form T (θ) = λ ψ.
Theorem 4.3 Suppose that Assumptions 1-5, 7-14, and 16-20 hold. Then, we have The covariance matrix in (4.7) needs to be estimated. To do so, CFT06 adopt the covariance estimation method proposed by Ai and Chen (2003). Since an infinite-dimensional optimization is involved in calculating S ψ 0 , we provide a sieve estimator of I * (θ 0 ) −1 . The sieve spaces for b and b ν can be the same as those for h and h ν , respectively. As in Ai and Chen (2003), we first estimate b * j 's in equation (4.6) by solving the following minimization problem: for all k = 1, 2, ..., d ψ , to obtain a consistent estimator of I * (ψ 0 ). We now summarize this result as follows: Theorem 4.4 Suppose that assumptions in Theorem 4.3 hold. Then, The proof of the theorem can be found in Theorem 5.1 in Ai and Chen (2003).

Example 2: Asymptotic Normality for the Conditional ATE
We now consider the conditional ATE, Proposition 4.1, we provide the asymptotic normality of the sieve plug-in estimator of the conditional ATE: Theorem 4.5 Let x ∈ supp(X) be given. Suppose that the conditions in Proposition 4.1 hold ||v|| , and the form of ∂AT E(θ 0 ;x)

∂θ
[v] is given by Furthermore, the asymptotic variance in (4.8) can be estimated as follows:

Monte Carlo Simulation and Sensitivity Analysis
It is important to show the practical relevance of the theoretical results developed in this paper.
In this section, we conduct a sensitivity analysis via Monte Carlo simulation exercises to provide guidance for empirical researchers. To this end, we investigate the finite sample performance of the sieve ML estimator of θ 0 and compare it with the performance of the parametric ML estimator under various DGPs and model specifications. Specifically, we are interested in the estimators of ψ 0 and the ATE. Here, ψ 0 is the finite-dimensional parameter and the ATE involves ψ 0 and the marginal distribution function of , which is an infinite-dimensional parameter in the case of semiparametric models.
The finite-dimensional parameter ψ 0 and the marginal effect of the treatment at some value of the covariate (i.e., ATE) are objects of interest in this class of models. For the ATE, we focus on the value evaluated at the means of covariates X. As a performance measure of the estimators, we consider the root mean squared errors (RMSEs) in our simulation.

Simulation Design
We consider various DGPs for both copulas and marginal distributions. The identification of the model requires that the copula function satisfy the stochastic increasing property. HV17 provide several examples of copulas, including the Gaussian, Frank, Clayton, and Gumbel copulas. We use these copulas to generate the sample. We also compare the performance of the parametric and semiparametric estimators when the marginal distributions are misspecified. To do so, we consider two marginal distributions for and ν: the standard normal distribution and a mixture of normal distributions. To estimate the parametric models, we specify normal distributions with unknown mean and variance parameters for the marginal distributions, owing to their popularity.
We refer to the parameters characterizing the marginal distributions as nuisance parameters. Note In the simulation, the data are generated as follows: There are several ways to normalize the location and scale of the model. We impose a restriction that X has no constant for the location normalization, and set the coefficient on X 1 to one for the scale normalization. As mentioned in the previous section, this normalization allows us to easily compare the performance of the parametric and semiparametric estimators, and can be used in more generalized models. To apply this normalization to our simulation design, α and β are assigned the value -1 in our simulation design.
The dependence structure between and ν is characterized by one-dimensional parameters in all copulas considered here, but the interpretation of the dependence parameter differs across the copulas. To resolve the difficulty in comparing the degree of dependence between and ν, we report the Spearman's ρ corresponding to the estimated dependence parameter in each copula specification. We also estimate the models with several values of Spearman's ρ to examine whether the performance of the estimators varies with the degree of dependence.
We consider two marginal distributions for the DGPs: (i) the standard normal distribution; and (ii) a mixture of normal distributions. We set the parameters so that E[ ] = E[ν] = 0 and V ar( ) = V ar(ν) = 1 in both cases. Specifically, and ν are generated from 0.6N (−1, σ 2 ) + 0.4N (1.5, σ 2 ) for some σ > 0 when the true marginal distribution is a mixture of normal distributions.
Although we assume that the copula is correctly specified, economic theory does not provide a justification for the choice of copula. In this simulation study, we examine the effect of copula misspecification on the performance of the estimators. 14 To examine the effect of copula misspecification, we only consider the family of copulas that satisfy the stochastic ordering property (Assumption 5) in order to ensure the identification of the model in our simulation.

Estimation of Parametric Models
The parametric models can be estimated by the standard ML method. Since it is common to use bivariate probit models for parametric estimation, we specify normal distributions for the marginal distributions. If we additionally choose the Gaussian copula, the model becomes a bivariate probit model. For the location and scale normalization of the parametric models, we adopt the same normalization as that of the semiparametric models. Thus, we can easily compare the performance of the estimators between the parametric and semiparametric models.

Estimation of Semiparametric Models
Recall that we assume that h j ∈ Λ p ([0, 1]). Therefore, we approximate h j to k=0 is the set of approximating functions for h j (·), and k nj is the number of approximating functions. Note that the approximation in (5.1) guarantees that´1 0 h j (x)dx = 1 by construction. We take the space of the polynomials as the sieve space for h and h ν . The 14 Misspecification problems in copula-based models have been documented using Monte Carlo simulations in the statistic literature (e.g., Kim et al. (2007a,b); Lawless and Yilmaz (2011)). In particular, Lawless and Yilmaz (2011) compare the performance of the parametric and semiparametric ML estimators in a copula-based model and show that the semiparametric two-step method outperforms the parametric estimation method when the copula function is misspecified. orders of the polynomials (k n and k nν ) are set to be proportional to n 1/7 . To incorporate the specification given in (4.1), we choose the standard normal distribution function for function G.

Simulation Results
We begin by examining the simulation results under a correct specification (i.e., the true marginal distributions are standard normal and the copula is correctly specified). Table 1 shows the simulation results for n = 500. We find that the ML estimators of ψ and the ATE perform well in the parametric models. The biases of the estimators are negligible and the variances are small. In addition, the performance of the sieve ML estimators of ψ and the ATE in the semiparametric models is as good as that in the parametric models. 15 For example, the RMSE of the estimator of δ 1 in the semiparametric model with the Gaussian copula is 0.4181, and that in the corresponding parametric model is 0.3982. We also find that the estimators of the ATE in both the parametric and semiparametric models perform well. Table 7 presents the simulation results with 1000 observations. The results show that the RMSEs decrease as the sample size increases and thus the estimators perform better with a larger sample size in both models. Now, we consider those cases where the marginal distributions are misspecified in the parametric models. Table 2 considers the case where the true marginal distributions are a mixture of normal distributions, but a researcher specifies them as normal distributions. Table 8 is obtained from simulation under the same situation, but with 1000 observations. From these tables, we can find that the RMSEs of the parametric ML estimators are larger than those of the sieve ML estimators. This implies that the parametric ML estimators could be misleading and that the sieve ML estimators are robust to this misspecification. Moreover, the parametric estimators of the ATE are substantially distorted under this misspecification. Note that the poor performance of the parametric estimators is attributed not only to bias, but also to variance. Specifically, the bias of the parametric estimator of the ATE with the Gaussian copula is 0.1377, which is about eight times larger than that of the corresponding sieve estimator, when the sample size is 500. These biases of the parametric estimators of the ATE are substantial in that they do not disappear, even when we increase the sample size. Comparing Tables 2 and 8, we find that the decreases in the RMSEs of the parametric estimators of the ATE with a larger sample size are mostly due to smaller variances. Therefore, the simulation results demonstrate that when the marginal distributions are misspecified, the sieve estimators outperform the parametric estimators in terms of the RMSE. In particular, the parametric estimates of the ATE can be significantly misleading when the marginal distributions are misspecified.
Finally, we examine the simulation results when both the copula and the marginal distributions are misspecified, reported in Tables 3-6 and 9-12. If both the copula and the marginal distributions are misspecified, the performance of the parametric ML estimators are comparable to, or slightly worse than that under marginal misspecification. Consider, for example, the case where the true copula function is the Frank copula and the sample size is 500. The estimators of ψ under both the copula and marginal misspecification (Table 4) have slightly larger RMSEs than the corresponding estimators under the marginal misspecification (Table 2). On the other hand, the performance of the estimators of the ATE varies across copula specifications. In particular, when the true DGP is based on the Gumbel copula, the copula and marginal misspecification has a significant effect on the performance of the parametric estimators of the ATE. The RMSEs of the estimators of the ATE under the copula and marginal misspecification (Table 6) are larger than those under the marginal misspecification (Table 2). Specifically, the RMSE of the parametric estimator of the ATE under the marginal misspecification is 0.1637 (Table 2) (Table 6). On the other hand, there is no clear evidence that the performance of the sieve ML estimators under both the copula and marginal misspecification is worse than that under misspecification of the marginal distributions. For example, when the true copula belongs to the Frank family but the copula is specified as the Gaussian or Gumbel copula, we can see that the RMSEs of the sieve ML estimators of the finite-dimensional parameters other than γ and the ATE under the copula and marginal misspecification (Table 4) are lower than those under the marginal misspecification (Table 2). In contrast, we can see from the same tables that the Clayton copula specification draws the opposite conclusion when the true copula is the Frank.
In general, no matter whether the copula is misspecified, we find that the sieve ML estimators outperform the parametric estimators in terms of the RMSE when the marginal distributions are misspecified.
Overall, the simulation results suggest that researchers are recommended to use the semiparametric models and the sieve ML estimation proposed in this paper when they are concerned about model misspecification. The following is the summary of the main findings from our simulation study: (i) When the model is correctly specified, the performance of the sieve ML estimators is comparable to that of the parametric ML estimators.
(ii) When the marginal distributions are misspecified, the sieve ML estimation is recommended in order to improve the performance.
(iii) The performance of the parametric ML estimators under both the copula and marginal misspecification can be worse than that under the marginal misspecification, and thus the semiparametric models are preferred to the parametric models.
(iv) Especially for the ATE, whenever the marginal distributions are misspecified, the parametric ML estimates can be significantly misleading.

Conclusions
In this paper, we propose semiparametric estimation and inference methods for generalized bivariate probit models. Specifically, we develop the asymptotic theory for the sieve ML estimators of semiparametric copula-based triangular systems with binary endogenous variables. We show that the sieve ML estimators are consistent and that their smooth functionals are √ n-asymptotically normal under some regularity conditions. This semiparametric estimation approach allows for flexibility in the models and thus provides robustness in estimation and inference.
We conduct a sensitivity analysis to examine how sensitive the estimation results are to model specifications. The results show that, overall, the semiparametric sieve ML estimators perform well in terms of both bias and variance. When the marginal distributions are misspecified, the sieve ML estimators substantially outperform the parametric ML estimators and the latter exhibit substantial bias. In particular, we find that the parametric estimates of the parameters involving the misspecified marginal distributions, such as the ATE, are highly misleading. When the model is correctly specified, we find that the performance of the sieve ML estimators is comparable to that of the parametric ones. When the copula is also misspecified, the distortion of the parametric estimates under the marginal misspecification can become even more severe, whereas the semiparametric estimates do not seem to be affected by this misspecification as long as the copula of the true DGP is within the stochastic ordering class. A related and interesting question is how the results would change if the data are not generated from this class of copulas.
Suppose that the following observational equivalence holds: One can easily show that these equations yield that t * 0 = 5/9, t † 0 = 1/9, t * 1 = 8/9, and t † 1 = 4/9. Consider the equation (7.10), which can be rewritten as t † , which is, in fact, the value of δ * 1 . In sum, the values of parameters that give the observationally equivalent fitted probabilities are This argument can be made slightly more general, and thus the counterexample more realistic, by relaxing Fε ∼ U nif (0, 1) and ρ * = 1. We show that a similar argument goes through with Fε being a general distribution function with a symmetric density function, and −1 ≤ ρ * ≤ 1 as long as the copula density is symmetric around u 2 = u 1 (i.e., the 45 degree line) and u 2 = 1 − u 1 .
The above counter-example to identification involves a parameter on the boundary of the parameter space (ρ * = 1), while the identification results in the paper assume that the parameter space is open and thus that ρ ∈ (−1, 1). We now show that the key idea of the argument remains the same with −1 < ρ * < 1. Suppose that the copula density is symmetric around u 2 = u 1 and u 2 = 1 − u 1 . The normal copula satisfies this condition for any ρ ∈ (−1, 1). Because of this condition, the symmetry of s † 0 and t * 1 (and of s † 1 and t * 0 ) around u 1 = 1/2 does not break at a different value of ρ * , even though the values of s † 0 , t * 1 , s † 1 , and t * 0 themselves change. Therefore, (7.13) continues to hold with ρ * = 1.

Hölder Norm and Hölder Class
Let C m (X ) be the space of m-times continuously differentiable real-valued functions on X . Let ζ ∈ (0, 1] and, given a d-tuple ω, let [ω] = ω 1 + ... + ω d . Denote the differential operator by D and let D ω = Letting p = m + ζ, the Hölder norm of h ∈ C m (X ) is defined as follows: where ζ is the Hölder exponent.

Sup-norm and Pseudo-metric d c
For any h ∈ H (or H ν ), define the sup-norm on H (or H ν ) as follows: 1] |h(t)|.

Fisher inner product and Fisher norm
Recall that V is the linear span of Θ − {θ 0 }. Define the Fisher inner product on the space V as for given v,ṽ ∈ V. Then, the Fisher norm for v ∈ V is defined as

Relationship between the Fisher norm and L 2 -norm
Note that for any θ 1 , θ 2 ∈ Θ, we have The following proposition is a modification of Theorem 3.1 in Chen (2007) and establishes the consistency of sieve M-estimator. 18 Proposition 7.1 Letθ n be the sieve extremum estimator defined in (4.2). Suppose that the following conditions hold : (ii) Θ is compact under d c (·, ·), and Q 0 (θ) is upper semicontinuous on Θ under d c (·, ·); (iii) The sieve spaces, Θ n , are compact under d c (·, ·) ; (iv) Θ k ⊆ Θ k+1 ⊆ Θ for all k ≥ 1, and there exists a sequence π k θ 0 ∈ Θ k such that We show that the conditions in Theorem 4.1 imply those in this proposition to prove consistency of the sieve estimator. We first need to verify that (i) the true parameter θ 0 is the unique maximizer of Q 0 (·) over Θ and that (ii) the sample log-likelihood function Q n (·) uniformly converges to Q 0 (·) over the sieve space in probability to establish the consistency of the sieve ML estimator. The following lemma shows that if the model with unknown marginal distributions are identified and some additional conditions are satisfied, then the true parameter θ 0 is the unique maximizer of Q 0 (·) over Θ.
Proof. By Theorem 2.3, the model parameter is identified. Under Assumption 9, we can see that for any θ ∈ Θ, |Q 0 (θ)| ≤ E|l(θ, W i )| ≤ y,d∈{0,1} E| log(p yd,XZ (θ))| < ∞, and thus the function Since the model is identified, it implies that for θ = θ 0 , there exists a set E ⊂ supp(X, Z) such that´E dP XZ > 0 and for some y, d ∈ {0, 1}, where P XZ is the distribution function of (X, Z). Thus, we have where the strict inequality holds by the fact that p yd,xz (θ) = p yd,xz (θ 0 ) on E and Jensen's inequality. Hence, θ 0 is the unique maximizer of Q 0 (·).
For any ω > 0, let N (ω, Θ n , d c ) be the covering numbers without bracketing of Θ n with respect to the pseudo-metric d c . We now establish the uniform convergence of Q n (·) to Q 0 over the sieve space.
By Lemma 2.5 in van de Geer (2000), we have log N ω 4 , H n , || · || ∞ ≤ k n log 1 + 32R ω under Assumption 12-(i); and hence if k n /n → 0. Since the condition k n /n = o(1) is imposed by Assumption 12-(i), the last condition in Condition 3.5M is also satisfied. In all, we have the uniform convergence of Q n to Q 0 over Θ k .
To finish proving Theorem 4.1, we verify the conditions in Proposition 7.1. By Lemmas 7.2 and 7.3, the conditions (i) and (v) in Proposition 7.1 are satisfied. Using (7.23) and Jensen's inequality, we can see that, for any θ,θ ∈ Θ, for some B < ∞. Thus, Q 0 (·) is continuous with respect to d c . Note that since the parameter space of the finite-dimensional parameter ψ, Ψ, is assumed to be compact in Assumption 8, the original parameter space Θ is compact under the d c , by Theorems 1 and 2 in Freyberger and Masten (2015), and thus the conditions (ii) and (iii) are satisfied with the specified parameter space and the norm. Since the condition (iv) is directly imposed, we have d(θ n , θ 0 ) = o p (1) by Proposition 7.1.

Proof of Theorem 4.2
To establish the convergence rate with respect to the norm || · || 2 , we consider the following assumption: Then, there exist B 1 , B 2 > 0 such that Assumption 18 implies that the L 2 -norm || · || 2 and the square-root of the KL divergence are equivalent.
We derive the convergence rate of the sieve M-estimator with respect to the norm || · || 2 by checking the conditions in Theorem 3.2 in Chen (2007). Since {W i } n i=1 is assumed to be i.i.d by Assumption 10, Condition 3.6 in Chen (2007) is satisfied. For Condition 3.7 in Chen (2007), we note that for a small 1 > 0 and for any θ ∈ Θ n such that ||θ − θ 0 || 2 ≤ 1 , we have By the same logic in (7.23), we have By Lemma 2 in Chen and Shen (1998), we have with 1 ≤ 1 and some constant B 1 , and thus Condition 3.7 in Chen (2007) is satisfied.

Proof of Proposition 4.1
We first provide some technical assumptions for the asymptotic normality. Let µ n (g) = 1 n n i=1 {g(W i )− E[g(W i )]} be the empirical process indexed by g. Let the convergence rate of the sieve estimator be δ n (i.e., ||θ n − θ 0 || = O p (δ n )).

Proof of Theorem 4.3
Take any λ ∈ R d Since λ was arbitrary, we obtain the result by Cramer-Wold device.