\section{Introduction}

Estimating the marginal treatment effects (MTE) is essential in treatment evaluation with an endogenous treatment.
The MTE can provide rich information on how treatment effects vary across economic agents in terms of their observed and unobserved characteristics.
Furthermore, after estimating the MTE, researchers can identify many treatment parameters of interest, such as the average treatment effects (ATE), local ATE (LATE), and policy-relevant treatment effects (PRTE) as some weighted averages of the MTE (\citealpmain{heckman1999local, heckman2005structural}).
Prior studies clearly demonstrate the usefulness of MTE methods in various empirical fields such as labor economics, health economics, and education economics (e.g., \citealpmain{basu2007use}; \citealpmain{carneiro2011estimating}; \citealpmain{cornelissen2016late}; \citealpmain{felfe2018does}, among others).

An important but often neglected issue in studies of treatment effects is the presence of potential interference between agents.
For example, when evaluating the effect of smoking behavior on health outcomes for couples, one partner's smoking behavior would affect the health outcome of the other.
That is, the well-known stable unit treatment value assumption (SUTVA) would not hold due to the treatment spillover.
In addition, it is natural to expect that one partner's smoking behavior interacts ``strategically'' with the other partner's smoking behavior, in the sense that one partner's action can directly affect the utility of smoking of the other partner, and vice versa.
This example suggests the presence of two different types of interference that need to be addressed: (i) treatment spillover and (ii) strategic interaction in the treatment decisions.

In reality, the co-existence of treatment spillover and strategic interaction should be fairly common.
For example, it is frequently observed that peer behavior such as delinquent activities among close friends has a significant impact on students' academic performance, invalidating the SUTVA when viewing delinquency as a treatment variable.
Additionally, many prior studies empirically show that strategic (social) interaction between friends is a primary cause of delinquency.
For another example, consider two large-scale rival retailers, which are deciding whether to enter a given local market.
Here, we would be interested in measuring the effects of the entry of these large retailers on other local businesses and communities, similar to \citemain{jia2008happens}.
This example is viewed as a special case of endogenous treatment models with the both types of interference.
\bigskip

This study aims to develop identification and estimation procedure for MTE models that allow both treatment spillover and strategic interaction in the treatment decisions.
In particular, we focus on an empirically relevant setup in which interactions occur between a pair of agents (e.g., couples, best friends, twins, or duopoly firms).
We postulate that they make decisions on their treatment status simultaneously in a binary game of complete information.
Within this framework, we formulate a set of sufficient conditions under which it is possible to point-identify the MTE parameters of interest.
To this end, we need to address the following two issues.
The first issue is the non-applicability of standard (i.e., ordered) {\it monotonicity} and {\it unordered monotonicity} (\citealpmain{imbens1994identification}; \citealpmain{heckman1999local, heckman2005structural}; \citealpmain{heckman2018unordered}).
The monotonicity conditions require that shifts in instrumental variables (IVs) determine the direction of changes in the treatment choices uniformly in all agents.
As \citetmain{vytlacil2002independence} and \citetmain{heckman2018unordered} show, these monotonicity conditions are equivalent to assuming a model in which the effects of the IVs on the treatment decisions are additively separable from unobservable error terms.
However, in our case, an agent's treatment choice nonlinearly depends on another agent's treatment choice through strategic interaction.
Consequently, a shift in one agent's IV can interact with both his/her and the partner's unobservables in a non-additive way, implying the failure of the monotonicity conditions.
The second issue is the possibility of multiple equilibria in the treatment decisions.
The presence of multiple equilibria leads to an {\it incomplete} econometric model (e.g., \citealpmain{tamer2003incomplete}; \citealpmain{lewbel2007coherency}; \citealpmain{ciliberto2009market}; \citealpmain{chesher2020structural}, among others) in the sense that the model-consistent treatment assignment is not unique.
The issue of incompleteness is a common problem in the literature on game model estimation; however, it is not yet well understood in the context of treatment evaluation.

Our identification strategy, which solves these two issues simultaneously, is a combination of extending the local IV (LIV) method in \citemain{heckman1999local, heckman2005structural} and imposing a stochastic equilibrium selection rule in the treatment decision game.
The key idea is to use local variations of the pair of IVs that satisfy a relevance condition and an exclusion restriction.
Specifically, we require the existence of player-specific continuous IVs that alter players' treatment status, but do not directly affect their outcomes.
Although this is a natural extension of the LIV method to a multidimensional space, it is still insufficient to point-identify the MTEs due to the presence of multiple equilibria.
We overcome this issue by explicitly introducing an equilibrium selection rule.

We present the following novel identification results.
First, we can point-identify the following MTE parameters: the {\it direct MTE}, in which only one agent's treatment status switches from untreated to treated, whereas the partner's remains unchanged; the
{\it indirect MTE}, in which only the partner's treatment status switches from untreated to treated, and the focal agent's status remains unchanged;
and the {\it total MTE}, in which the treatment status of both players switches from untreated to treated.
In contrast to the conventional MTE framework, our MTEs can reveal treatment effect heterogeneity in terms of the ``pair'' of unobservables.
Second, we show that the presence of multiple equilibria results in over-identification of MTE parameters.
Importantly, this over-identifiability provides an opportunity to improve the efficiency of the MTE estimation and to derive a testable implication for the equilibrium selection assumption.
Third, based on the identification of the MTE, we present the identification of the other treatment parameters, including the LATE and PRTE.
Finally, while treatment evaluation is our main concern, we also provide new supplementary results for the identification of binary game models of complete information.

Our identification is constructive in that we can estimate the MTEs directly by following the identification strategy.
We propose a two-step semiparametric procedure for estimating the MTEs.
In the first step, we estimate the parameters in the treatment decision game using a maximum likelihood (ML) approach.
Using the ML estimates, we estimate the MTEs in the second step by employing semiparametric series (sieve) techniques to estimate the outcome equation.
The proposed estimator is uniformly consistent with the optimal convergence rate and is asymptotically normally distributed.
In addition, the estimator possesses an oracle property in the sense that its limiting distribution is the same as that of the infeasible estimator where the parameters in the treatment decision game are treated as known.

To illustrate our methods empirically, we investigate the impacts of an opposite-gender best friend's delinquency (e.g., smoking and drinking alcohol) on adolescents' academic outcome (GPA: grade point average).
Following the literature (e.g., \citealpmain{card2013peer}), we model the decision to participate in such risky activities as a complete information game.
Our methods reveal new empirical evidences that the direct treatment effect of risky behaviors on GPA is significantly negative for female students and that the total treatment effect is significant for both genders.
In particular, we find that an opposite-gender friend's delinquency more easily influences male students' GPAs than those of female students.
\bigskip

To our best knowledge, few studies address treatment evaluation in the presence of strategic interaction modeled explicitly as games.
One important exception is \citetmain{balat2019multiple}, who investigate the partial identification of the ATE, differently from our focus on the point identification of the MTEs.
The key identification assumptions in \citetmain{balat2019multiple} are nonparametric shape restrictions and variations in possibly discrete IVs, while we rely principally on a stochastic equilibrium selection rule and the existence of player-specific continuous IVs to achieve the point identification.

Another closely related study is \citetmain{lee2018identifying}, who show that it is possible to identify MTE parameters by modeling the treatment selection with a set of threshold-crossing rules.
Their identification result is obtained by computing the changes in the expected outcome with respect to local variations of ``all'' threshold variables, even including those irrelevant to a particular treatment status.
Appendix C of \citetmain{lee2018identifying} illustrates how this approach can be used in the case when the treatments are determined in a binary game.
In contrast, our approach uses only a set of the threshold variables which are needed to well-define a treatment status of interest.
This means that our identification result holds with a fewer number of IVs than theirs, and this distinction would be crucial from a practical perspective (see Remark \ref{remark:leesalanie} for more details).
\bigskip

The rest of the paper is organized as follows.
In Section \ref{sec:model}, we introduce our model and review the incompleteness problem for discrete game models.
In Section \ref{sec:identification}, we establish the identification of the MTE parameters.
We develop our two-step semiparametric estimators for the MTEs and show their asymptotic properties in Section \ref{sec:estimation}.
Section \ref{sec:numeric} presents the numerical illustrations, including a set of Monte Carlo simulations and the empirical analysis of adolescents' academic performance.
Section \ref{sec:conclusion} concludes the paper.
The proofs of all theorems are provided in Appendix \ref{sec:proof}, and other supplementary technical results in Appendices \ref{appendix:supptech} and \ref{appendix:game}.
In Appendix \ref{sec:severalparameter}, we discuss the identification of several treatment parameters besides the MTE.
Appendix \ref{sec:MC} presents detailed information about the Monte Carlo experiments in Section \ref{sec:numeric}.
Finally, we provide the supplementary tables and figures for the empirical analysis in Appendix \ref{sec:supp}.

\section{Model} \label{sec:model}

In this section, we introduce our treatment effect model with strategic interaction.
We denote a player by $j \in \{1, 2\}$ and his/her partner (or opponent) by $-j$.
We aim to evaluate the effects of player $j$'s treatment $D_j \in \{0, 1\}$ and/or partner's treatment $D_{-j} \in \{0, 1\}$ on player $j$'s outcome $Y_j$ and/or the partner's outcome $Y_{-j}$.
The outcomes may or may not be common to both players; that is, we allow both $Y_j = Y_{-j}$ and $Y_j \neq Y_{-j}$.

Let $Y_j^{(d_j, d_{-j})}$ be the potential outcome for player $j$ when his/her own treatment status is $D_j = d_j$ and the partner's is $D_{-j} = d_{-j}$.
Then, the observed outcome can be written as $Y_j = \sum_{d_j = 0}^{1}\sum_{d_{-j} = 0}^{1} I^{(d_j, d_{-j})} Y_j^{(d_j, d_{-j})}$, where $I^{(d_j, d_{-j})} \coloneqq \mathbf{1}\{(D_j, D_{-j}) = (d_j, d_{-j})\}$.
Suppose that player $j$'s potential outcome can be written as
\begin{align}\label{eq:model-Y1}
Y_j^{(d_j, d_{-j})} = \mu_j^{(d_j, d_{-j})} \left( X_j, U_j^{(d_j, d_{-j})} \right),
\end{align}
where $X_j \in \mathbb{R}^{\mathrm{dim}(X)}$ is a vector of the observable covariates, $U_j^{(d_j, d_{-j})} \in \mathbb{R}$ is an unobservable random variable, and $\mu_j^{(d_j, d_{-j})}$ is an unknown structural function.
The covariates $X_1$ and $X_2$ may contain common elements as well as some player-specific elements.
For simplicity, we assume that the dimensions of $X_1$ and $X_2$ are both equal to $\mathrm{dim}(X)$, and the same assumption will be made for the other variables.
The model specification \eqref{eq:model-Y1} is fairly general in that $U_j^{(d_j, d_{-j})}$ and $\mu_j^{(d_j, d_{-j})}$ can depend on the player type $j$ and treatment status $(d_j, d_{-j})$.

\subsection{Strategic interaction in the treatment decisions}

To account for the strategic interaction between players, we extend the latent index model of \citetmain{heckman1999local, heckman2005structural}.
Suppose that player $j$'s treatment is determined by
\begin{align}\label{eq:model-D1}
	D_j = \mathbf{1} \left\{ \pi_j (D_{-j}, W_j) \ge \varepsilon_j \right\},
\end{align}
where $W_j \coloneqq (X_j^\top, Z_j^\top)^\top$ is a vector including $X_j$ and the instruments $Z_j \in \mathbb{R}^{\mathrm{dim}(Z)}$, $\varepsilon_j \in \mathbb{R}$ is a continuously distributed unobserved random variable, and $\pi_j$ is an unknown function.
We do not restrict the dependence structure between $\varepsilon_j$ and $U_j^{(d_j, d_{-j})}$ for $(d_j, d_{-j}) \in \{0, 1\}^2$, which is the source of endogeneity.
Without endogeneity, one may relatively easily identify some causal parameters by standard approaches based on an unconfoundedness assumption, even in the presence of strategic interaction.
We assume that the IVs include at least one player-specific variable so that the distribution of $Z_{j}$ is non-degenerate given $Z_{-j}$.
The instruments $Z_j$ are assumed to be independent of the unobservables $\varepsilon_j$ and $U_j^{(d_j, d_{-j})}$ for $(d_j, d_{-j}) \in \{0, 1\}^2$ conditional on $X = (X_1, X_2)$.
In addition, when using an LIV method to identify treatment effects, we require that not all the elements of $Z_j$ are discrete.

For the subsequent analysis, it is convenient to transform \eqref{eq:model-D1} as follows.
Firstly, since $D_{-j}$ is a binary variable, we can write
\begin{align}\label{eq:pi}
	\pi_j (D_{-j}, W_j) = \pi_j^0(W_j) + D_{-j} \left( \pi_j^1 (W_j) - \pi_j^0(W_j) \right),
\end{align}
where $\pi_j^0(\cdot) \coloneqq \pi_j (0, \cdot)$ and $\pi_j^1 (\cdot) \coloneqq \pi_j (1, \cdot)$.\footnote{
	The expression in \eqref{eq:pi} clearly indicates that the additive separability imposed on the treatment decision model \eqref{eq:model-D1} implies the interaction effect to be a function of only the observed variables.
	In the context of our empirical setting, this requires that the impact of friends' delinquency is independent of own unobserved factors, such as his/her latent attitude towards risky activities.
	With few exceptions (e.g., \citealpmain{kline2015identification}), models where the strategic interaction effects can depend on some unobservables have not been studied in detail in the literature, and how to cope with this issue is an important future task.
 	}
Let
\begin{align}\label{eq:definition-P}
	V_j \coloneqq F_{\varepsilon_j} ( \varepsilon_j | X ),
	\qquad
	P_j^0 \coloneqq F_{\varepsilon_j} \left( \pi_j^0(W_j) \Big| X \right),
	\qquad
	P_j^1 \coloneqq F_{\varepsilon_j} \left( \pi_j^1(W_j) \Big| X \right),
\end{align}
where $F_{\varepsilon_j}(\cdot | X)$ is the conditional cumulative distribution function (CDF) of $\varepsilon_j$ given $X$.
By construction, $V_j $ is distributed as $\text{Uniform}[0,1]$ conditional on $X$.
Combining these and \eqref{eq:pi}, we can rewrite \eqref{eq:model-D1} as
\begin{align}\label{eq:model-D2}
	D_j = \mathbf{1} \left\{ P_j^0 + D_{-j} (P_j^1 - P_j^0) \ge V_j \right\}.
\end{align}
The treatment decision (\ref{eq:model-D2}) is equivalent to $D_j(D_{-j})$ almost surely (a.s.), where $D_j(d_{-j}) = \argmax_{d_j \in \{0, 1\}} u_j(d_j, d_{-j})$, and  $u_j(d_j, d_{-j}) \coloneqq d_j \cdot (P_j^0 + d_{-j}  (P_j^1 - P_j^0) - V_j)$.
That is, we can view $u_j(d_j, d_{-j})$ as the payoff function for player $j$ and $D_j(d_{-j})$ as the best response given $D_{-j} = d_{-j}$.\footnote{
	If players $1$ and $2$ are collaborating partners such that their payoffs are transferable to each other, the players action may be jointly characterized by $(D_1, D_2) = \argmax_{(d_1,d_2) \in \{0, 1\}^2} [u_1(d_1, d_2) + u_2(d_2, d_1)]$.
	This type of cooperative decision-making model is proposed in \cite{lewbel2007coherency}.
	As shown in \cite{lewbel2007coherency}, this model is \textit{complete} (i.e., no multiple model-consistent solutions) and, thus, is analytically simpler than non-cooperative models including \eqref{eq:model-D2}.
	}
Here, we implicitly assume a complete information setup, where $(W_j, W_{-j}, \varepsilon_j, \varepsilon_{-j})$ is common knowledge to both players.
The set of realized treatments $(D_j, D_{-j})$ generated by model \eqref{eq:model-D1} is a fixed point of the best responses; that is, it is a pure strategy Nash equilibrium of the game with the payoff matrix given in Table \ref{tab:payoff}.\footnote{
	Most of studies in the game econometrics literature focus on pure strategy Nash equilibria as the solution concept for $2 \times 2$ games, except when there are no  Nash equilibria in pure strategies, for example, because of asymmetric strategic interaction (cf., e.g., \citealpmain{bjorn1984simultaneous}; \citealpmain{tamer2003incomplete}; \citealpmain{depaula2013econometric}).
	For some motivations of this choice, we refer readers, for example, to Section 6 of \citemain{bjorn1984simultaneous}.
}

\begin{table}[!ht]
	\caption{Payoff matrix}
	\label{tab:payoff}
	\centering
	\begin{tabular}{c|c|c}
		\hline
		& $D_2 = 0$ & $D_2 = 1$ \\ \hline
		$D_1 = 0$ & $(0, 0)$ & $(0, P_2^0 - V_2)$ \\ \hline
		$D_1 = 1$ & $(P_1^0 - V_1, 0)$ & $(P_1^1 - V_1, P_2^1 - V_2)$\\ \hline
	\end{tabular}
\end{table}

For simplicity, we assume that whether the treatment decisions are strategic complements or substitutes is known to us a priori.
Note that this assumption is empirically testable.\footnote{
	For example, \citetmain{aradillas2019nonparametric} has developed nonparametric tests for the presence and the direction of strategic interaction effects in $2 \times 2$ games of complete information.
	For another example, one may use a Vuong-type model selection test for non-nested alternatives (e.g., \citealpmain{hsu2017model}; \citealpmain{schennach2017simple})

	}
The treatment decisions are complements (resp. substitutes) if and only if the strategic interaction effect $P_j^1 - P_j^0$ is positive (resp. negative).
We hereafter restrict our main focus to the case of complementarity.
\begin{assumption}\label{as:complement}
	For both $j = 1$ and $2$, $P_j^1 - P_j^0 > 0$, a.s.
\end{assumption}
This assumption is in line with our empirical illustration.
That is, there is a consensus in the literature that friends' delinquency tend to promote one's own delinquency.
It is not difficult to modify our approach for the case of strategic substitutes.

\subsection{Incompleteness}

A major difficulty in our model is the {\it incompleteness} of the treatment decision model.
Under strategic complementarity, we have the following relationship between $(D_1, D_2)$ and $(V_1, V_2)$:
\begin{equation} \label{eq:equilibrium}
	\begin{array}{lcl}
		(D_1, D_2) = (1, 1) & \Longrightarrow & V_1 \le P_1^1, \; V_2 \le P_2^1, \\
		(D_1, D_2) = (1, 0) & \Longleftrightarrow & V_1 \le P_1^0, \; V_2 > P_2^1, \\
		(D_1, D_2) = (0, 1) & \Longleftrightarrow & V_1 > P_1^1, \; V_2 \le P_2^0, \\
		(D_1, D_2) = (0, 0) & \Longrightarrow & V_1 > P_1^0, \; V_2 > P_2^0.
	\end{array}
\end{equation}
Figure \ref{fig:nash-comp1} visually summarizes these relationships.
As shown in the figure, the space of $V = (V_1, V_2)$ cannot be partitioned into non-overlapping regions associated with the four alternative realizations of $D = (D_1, D_2)$.
Both $D = (1, 1)$ and $D = (0, 0)$ can occur when $P_1^0 < V_1 \le P_1^1$ and $P_2^0 < V_2 \le P_2^1$, and in this case, the value of $D$ is not uniquely determined (i.e., multiple equilibria).
This non-uniqueness of model-consistent decisions is called incompleteness and has been extensively studied in the literature on simultaneous equation models for discrete outcomes.

\begin{figure}[!h]
	\centering
	\fbox{\includegraphics[width = 9cm, bb = 0 0 730 550]{fig_nash_comp1.png}}
	\caption{Nash equilibrium under strategic complementarity.}
	\label{fig:nash-comp1}
\end{figure}

To handle this incompleteness problem, there are essentially three approaches to achieve point identification in the literature.\footnote{
	Partial identification of game econometric models has been considered, for example, in \citemain{ciliberto2009market} and \citemain{aradillas2019inference}.
	See \citetmain{depaula2013econometric} for an excellent survey on this topic.
}
The first approach is to focus only on the outcomes that can occur as unique equilibria (e.g., \citealpmain{bresnahan1990entry}; \citealpmain{berry1992estimation}).
The second approach is to use an identification-at-infinity strategy based on large support regressors and exclusion restrictions (e.g., \citealpmain{tamer2003incomplete}; \citealpmain{bajari2010identification}; \citealpmain{kline2015identification}).
The third approach is to explicitly introduce a stochastic (or possibly deterministic) equilibrium selection mechanism (e.g., \citealpmain{bjorn1984simultaneous}; \citealpmain{kooreman1994estimation}; \citealpmain{soetevent2007discrete}; \citealpmain{bajari2010identification}; \citealpmain{card2013peer}; \citealpmain{khan2018information}).
Especially among them, the last approach allows us to identify the choice behavior in the region of multiplicity and to perform counterfactual exercises.
Since the estimation of treatment effects inherently involves counterfactual inference, it is natural to employ the stochastic equilibrium selection approach in our situation.

\begin{remark}[Monotonicity]
	From Figure \ref{fig:nash-comp1} and \eqref{eq:equilibrium}, we can find why the standard monotonicity and unordered monotonicity conditions are not applicable to our situation.
	Consider the pairs of agents in the region of multiple equilibria and those in the region of $D = (1,0)$ for a given $Z = z$, where $Z = (Z_1, Z_2)$.
	Suppose that a shift of $Z$ from $z$ to $z'$ induces the former set of pairs to $D = (0, 1)$ and the latter to the multiple equilibria.
	In this case, some of the former pairs would switch their treatment statuses from $D = (0, 0)$ to $D = (0, 1)$, while some of the latter pairs would switch from $D = (1, 0)$ to $D = (0, 0)$.
	Thus, the change in the IV values results in both inflows into $D = (0, 0)$ and outflows from it, implying the violation of the monotonicity conditions.

\end{remark}

\section{Identification}\label{sec:identification}

In this section, the identification of the MTE parameters is discussed.
We here refer to the treatment decision game \eqref{eq:model-D2} as the ``first stage'' and the realization of the outcome \eqref{eq:model-Y1} as the ``second stage''.
To focus on the identification of the parameters in the second stage, we treat the parameters in the first stage as known to us.
More specifically, we assume the following:
\begin{assumption}
	\hfill
	\label{as:PV}
	\begin{enumerate}[(i)]
		\item For both $j = 1$ and $2$, $P_j^0$ and $P_j^1$ are identified.
		\item The conditional CDF and density of $V$ given $X = x$, denoted respectively by $H(v_1, v_2 | x)$ and $h(v_1, v_2 | x)$, are identified (or known) functions.
	\end{enumerate}
\end{assumption}
Assumption \ref{as:PV} can be justified by the existing identification results in the literature with minor modifications (e.g., \citealpmain{bresnahan1990entry}; \citealpmain{berry1992estimation}; \citealpmain{tamer2003incomplete}; \citealpmain{kline2015identification}).
In particular, note that this assumption can be satisfied without imposing auxiliary assumptions on equilibrium selection.

For what follows, we introduce additional notations: $Y = (Y_1, Y_2)$, $W = (W_1, W_2)$, $\varepsilon = (\varepsilon_1, \varepsilon_2)$, and $U^{(d_1, d_2)} = ( U_1^{(d_1, d_2)}, U_2^{(d_2, d_1)} )$ for $(d_1, d_2) \in \{0, 1\}^2$.
We here introduce the exclusion restriction and relevance condition for the IVs.
\begin{assumption}
	\hfil
	\label{as:IV}
	\begin{enumerate}[(i)]
		\item For all $(d_1, d_2) \in \{0, 1\}^2$, the instruments $Z$ are excluded from the structural functions in \eqref{eq:model-Y1} and are independent of the unobservables $(\varepsilon, U^{(d_1, d_2)})$ given $X$.
		\item For both $j = 1$ and $2$, the instruments $Z_j$ contain a player-specific continuous variable such that $\pi_j^0(W_j)$ and $\pi_j^1(W_j)$ are non-degenerate and continuously distributed given $X$.
	\end{enumerate}
\end{assumption}

Assumption \ref{as:IV}(i) requires that the instruments $Z$ do not directly affect the outcomes $Y$ and are conditionally independent of the error terms.
By construction, the transformed error $V$ is also conditionally independent of $Z$ given $X$.
Thus, the conditional distributions of $V$ and $U^{(d_1, d_2)}$ given $W$ are identical to those given $X$.
Assumption \ref{as:IV}(ii) requires that the instrument $Z_j$ must contain a continuous non-trivial determinant of $D_j$ specific to player $j$.

We provide below a series of identification results only for player 1 (the results for player 2 are symmetric and thus omitted).
We first discuss the identification of the conditional mean of the potential outcome:
\begin{align}\label{eq:m}
	m_1^{(d_1, d_2)}(x, p_1, p_2)
	\coloneqq E[Y_1^{(d_1, d_2)} | X = x, V_1 = p_1, V_2 = p_2].
\end{align}
We call this function the marginal treatment response (MTR) function, as in \citetmain{mogstad2018using}.
Since the number of treatment patterns is four (i.e., $D \in \{ (0, 0), (1, 0), (0, 1), (1, 1) \}$), there are six distinct MTE parameters for each player.
From the point of view of player 1, the parameters of interest would be the \textit{direct MTE}, \textit{indirect MTE}, and \textit{total MTE}, which are defined as
\begin{align*}
	\text{MTE}_{\text{direct}}^{(d_2)}(x, p_1, p_2) &\coloneqq m_1^{(1, d_2)}(x, p_1, p_2)  - m_1^{(0, d_2)}(x, p_1, p_2)   \;\; \text{ for $d_2 \in \{0,1\}$},\\
	\text{MTE}_{\text{indirect}}^{(d_1)}(x, p_1, p_2) &\coloneqq m_1^{(d_1, 1)}(x, p_1, p_2)  - m_1^{(d_1, 0)}(x, p_1, p_2)    \;\; \text{ for $d_1 \in \{0,1\}$},\\
	\text{MTE}_{\text{total}}(x, p_1, p_2) &\coloneqq m_1^{(1, 1)}(x, p_1, p_2)  - m_1^{(0, 0)}(x, p_1, p_2),
\end{align*}
respectively.
Notably, these MTE parameters are informative about the variation of the treatment effects in terms of the pair of unobservables $V_1$ and $V_2$, in contrast to the conventional MTE framework.
The identification of the MTE parameters is straightforward once the MTR function for each $(d_1,d_2)$ is identified.

\subsection{Identification of the marginal treatment response functions}\label{subsec:condmean}

We first consider the cases of $D = (1, 0)$ and $D = (0, 1)$.
In these cases, $V$ resides in the regions with a unique equilibrium.
In the following, we investigate the identification of the MTR functions:
\begin{align*}
\begin{array}{ll}
	m_1^{(1, 0)}(x, p_1^0, p_2^1) & \text{ for $(p_1^0, p_2^1) \in \text{supp}[P_1^0, P_2^1| X = x, D = (1, 0)]$,}\\
	m_1^{(0, 1)}(x, p_1^1, p_2^0) & \text{ for $(p_1^1, p_2^0) \in \text{supp}[P_1^1, P_2^0| X = x, D = (0, 1)]$,}
\end{array}
\end{align*}
where $ \text{supp}[P | A = a]$ denotes the conditional support of random variable(s) $P$ given the event $A =a$.
The former is the MTR function for player 1 when $D = (1,0)$ is realized and $V$ is located at point A in Figure \ref{fig:nash-comp2}.
At this point, players 1 and 2 are at the ``margin'' of the actions $D = (1, 0)$, $(1, 1)$, and $(0, 0)$ given that $(P_1^0, P_2^1) = (p_1^0, p_2^1)$ holds.
Thus, only a small deviation from $(p_1^0, p_2^1)$ may result in different treatment decisions from $D = (1,0)$.
The latter MTR function can be interpreted similarly at point D.
Further, we define the following functions:
\begin{align*}
\begin{array}{ll}
	\psi_1^{(1, 0)}(x, p_1^0, p_2^1) \coloneqq E[ I^{(1, 0)} Y_1| X = x, P_1^0 = p_1^0, P_2^1 = p_2^1] & \text{ for $(p_1^0, p_2^1) \in \text{supp}[P_1^0, P_2^1| X = x, D = (1, 0)]$},\\
	\psi_1^{(0, 1)}(x, p_1^1, p_2^0) \coloneqq E[ I^{(0, 1)} Y_1| X = x, P_1^1 = p_1^1, P_2^0 = p_2^0] & \text{ for $(p_1^1, p_2^0) \in \text{supp}[P_1^1, P_2^0| X = x, D = (0, 1)]$},
\end{array}
\end{align*}
which are directly identified from the data.
Hereafter, when there is no confusion, we suppress the subscript $1$ from $m_1^{(d_1, d_2)}$ and $\psi_1^{(d_1, d_2)}$ for notational simplicity.

\begin{figure}[!ht]
	\centering
	\fbox{\includegraphics[width = 9cm, bb = 0 0 730 550]{fig_nash_comp2.png}}
	\caption{Points at which the MTR functions are identified.}
	\label{fig:nash-comp2}
\end{figure}

The next theorem presents formal identification results for the MTR functions for $D = (1,0)$ and $(0,1)$.\footnote{
	Similarly, replacing $Y_1$ in $\psi^{(1, 0)}(x, p_1^0, p_2^1)$ with $1$ leads to $h(p_1^0, p_2^1| x) = - \partial^2 \Pr[D = (1, 0) | X = x, P_1^0 = p_1^0, P_2^1 = p_2^1] /(\partial p_1^0 \partial p_2^1)$.
	Similar results hold for the other realizations of $D$; as such, they can be used to check the specification of the distribution of $V$ imposed in Assumption \ref{as:PV}.
}

\begin{theorem}\label{thm:unique}
	Suppose that Assumptions \ref{as:complement}, \ref{as:PV}, and \ref{as:IV} hold.
	Then, if $m^{(1, 0)}(x, v_1, v_2)$, $m^{(0,1)}(x, v_1, v_2)$, and $h(v_1, v_2 | x)$ are continuous in $(v_1, v_2)$, the MTR functions for $D = (1, 0)$ and $D = (0, 1)$ can be identified as
	\begin{align*}
		m^{(1, 0)}(x, p_1^0, p_2^1)
		&= - \frac{1}{h(p_1^0, p_2^1 | x)} \frac{\partial^2 \psi^{(1,0)}(x, p_1^0, p_2^1)}{\partial p_1^0 \partial p_2^1} ,\\
		m^{(0, 1)}(x, p_1^1, p_2^0)
		&= - \frac{1}{h(p_1^1, p_2^0 | x)} \frac{\partial^2  \psi^{(0,1)}(x, p_1^1, p_2^0)}{\partial p_1^1 \partial p_2^0}.
	\end{align*}
\end{theorem}

A sketch of the proof for $m^{(1, 0)}(x, p_1^0, p_2^1)$ would be helpful to illustrate the identification strategy used here.
It can be shown that
\begin{align*}
	\psi^{(1, 0)}(x, p_1^0, p_2^1)
	&= E[Y_1^{(1, 0)}| X = x, V_1 \le p_1^0, V_2 > p_2^1] \cdot \Pr[V_1 \le p_1^0, V_2 > p_2^1 | X = x]\\
	&=  \int_{p_2^1}^1 \int_0^{p_1^0} m^{(1, 0)}(x, v_1, v_2) h(v_1, v_2 | x) \mathrm{d}v_1 \mathrm{d}v_2.
\end{align*}
Partially differentiating both sides with respect to $p_1^0$ and $p_2^1$ yields the desired result.
In this derivation, we need Assumption \ref{as:complement} so that the treatment status $D = (1, 0)$ may be uniquely linked with the region of $V$ (i.e., the upper left region in Figures \ref{fig:nash-comp1} and \ref{fig:nash-comp2}).
Assumption \ref{as:IV}(i) is required so that the first equality may be derived and the conditions by $P_1^0$ and $P_2^1$ be ignored.
Assumption \ref{as:IV}(ii) ensures that the partial differentiation with respect to $p_1^0$ and $p_2^1$ is well defined.
\bigskip

We move on to the identification of the MTR functions for $D = (0, 0)$ and $D = (1, 1)$.
In contrast to the above cases, the treatment statuses $D = (0, 0)$ and $D = (1, 1)$ are not uniquely selected when $V$ is in the multiple equilibria region.
Hence, we introduce an assumption on the equilibrium selection that leads to the point identification of the parameters of interest.
Write the region of $V$ where multiple equilibria occur as $\mathcal{V}_{\text{mul}}(\mathbf{P}) \coloneqq \{ (v_1, v_2) \in [0, 1]^2 : P_1^0 < v_1 \le P_1^1, P_2^0 < v_2 \le P_2^1 \}$, where $\mathbf{P} = (P_1^0, P_1^1, P_2^0, P_2^1)$.
Further, denote the region of $V$ where $D = (0, 0)$ is uniquely selected by $\mathcal{V}_{\text{uni}}^{(0, 0)}(\mathbf{P}) \coloneqq \{ (v_1, v_2) \in [0, 1]^2 : P_1^0 < v_1, P_2^0 < v_2 \} \backslash \mathcal{V}_{\text{mul}}(\mathbf{P})$.
For notational convenience, we often write $\mathcal{V}_{\text{mul}} = \mathcal{V}_{\text{mul}}(\mathbf{p})$ and $\mathcal{V}_{\text{uni}}^{(0, 0)} = \mathcal{V}_{\text{uni}}^{(0, 0)}(\mathbf{p})$ by suppressing their dependence on $\mathbf{P} = \mathbf{p}$.

\begin{assumption}
	\hfil
	\label{as:multiple}
	\begin{enumerate}[(i)]
		\item There exist a constant $\lambda \in [0, 1]$ and an unobserved random variable $\epsilon$ distributed as $\text{Uniform}[0, 1]$ independent of $(W, \varepsilon, U^{(d_1, d_2)})$ for $(d_1,d_2)\in\{0,1\}^2$ such that when $V \in \mathcal{V}_{\text{mul}}(\mathbf{P})$, $D = (0, 0)$ occurs if and only if $\epsilon \le \lambda$; that is, $D = (0, 0) \Longleftrightarrow V \in \mathcal{V}_{\text{uni}}^{(0, 0)}(\mathbf{P}) \lor (V \in \mathcal{V}_{\text{mul}}(\mathbf{P}) \land \epsilon \le \lambda)$.
		\item $\lambda \in [0, 1]$ is an identified (or known) parameter.
	\end{enumerate}
\end{assumption}
This assumption states that $D = (0, 0)$ is observed with probability $\lambda \in [0, 1]$ in the multiple equilibria situation.
Similarly to Assumption \ref{as:PV}, we assume that $\lambda$ is identifiable or known.
For example, some authors assume that under multiple equilibria, one of model-consistent actions is selected uniformly at random (e.g., \citealpmain{bjorn1984simultaneous}; \citealpmain{soetevent2007discrete}).
In particular, \citemain{card2013peer} introduce this assumption in their empirical application on risky behaviors.
If we adopt the same assumption, we can set $\lambda = 0.5$ a priori.
For another example, one may assume that the realized treatment status corresponds to the ``largest'' Nash equilibrium, as in \citemain{xu2015estimation}.
In this case, since $u_j(1,1) > u_j(0,0)$ in the multiple equilibria region for both players, $D = (1,1)$ holds almost surely (i.e., $\lambda = 0$).
Although we introduce Assumption \ref{as:multiple} for tractability, we can instead consider a more general structure on the equilibrium selection.
For example, we can establish point identification results even when  the assumption is relaxed to allow $\lambda$ to depend on $X$.
In more general case, $\lambda$ may depend also on some unobservable factors that are possibly correlated with other variables; nevertheless, it would be difficult to identify such very general selection from data alone (cf. \citealpmain{jun2019counterfactual}).
Even without any assumptions on the equilibrium selection, it would be possible to partially identify some causal parameters, as in \citemain{balat2019multiple}; see also Remark \ref{remark:partial} below.

To state the next theorem, define the following function which can be identified from the data:
\begin{align*}
	\psi^{(d_1, d_2)} (x, \mathbf{p})
	\coloneqq E[ I^{(d_1, d_2)} Y_1 | X = x, \mathbf{P} = \mathbf{p}]
	\quad
	\text{for $\mathbf{p} \in \text{supp}[\mathbf{P}|X = x, D = (d_1, d_2)]$.}
\end{align*}

\begin{theorem}\label{thm:multiple}
	Suppose that Assumptions \ref{as:complement}, \ref{as:PV}, \ref{as:IV}, and \ref{as:multiple} hold.
	Then, if $m^{(0, 0)}(x, v_1, v_2)$, $m^{(1, 1)}(x, v_1, v_2)$, and $h(v_1, v_2 | x)$ are continuous in $(v_1, v_2)$, the MTR functions for $D = (0, 0)$ and $D = (1, 1)$ can be identified as
	\begin{alignat*}{2}
		 m^{(0, 0)}(x, p_1^0, p_2^0) &= \frac{1}{\lambda h(p_1^0, p_2^0 | x)} \frac{\partial^2 \psi^{(0, 0)}(x, \mathbf{p})}{\partial p_1^0 \partial p_2^0}  && \quad \text{for $\lambda > 0$,}\\
		 m^{(0, 0)}(x, p_1^{d_1}, p_2^{d_2}) &= \frac{1}{(1 - \lambda) h(p_1^{d_1}, p_2^{d_2}|x)} \frac{\partial^2 \psi^{(0, 0)}(x, \mathbf{p})}{\partial p_1^{d_1} \partial p_2^{d_2}}  && \quad \text{for $\lambda < 1$ and $d_1 \neq d_2$,} \\
		 m^{(0, 0)}(x, p_1^1, p_2^1) &= - \frac{1}{(1 -\lambda) h(p_1^1, p_2^1|x)} \frac{\partial^2 \psi^{(0, 0)}(x, \mathbf{p})}{\partial p_1^1 \partial p_2^1}  && \quad \text{for $\lambda < 1$},
	\end{alignat*}
	and
	\begin{alignat*}{2}
		 m^{(1, 1)}(x, p_1^0, p_2^0) &= - \frac{1}{\lambda h(p_1^0, p_2^0 | x)} \frac{\partial^2 \psi^{(1, 1)}(x, \mathbf{p})}{\partial p_1^0 \partial p_2^0}  && \quad \text{for $\lambda > 0$,}\\
		 m^{(1, 1)}(x, p_1^{d_1}, p_2^{d_2}) &= \frac{1}{\lambda h(p_1^{d_1}, p_2^{d_2}|x)} \frac{\partial^2  \psi^{(1, 1)}(x, \mathbf{p})}{\partial p_1^{d_1} \partial p_2^{d_2}}  && \quad \text{for $\lambda > 0$  and $d_1 \neq d_2$,} \\
		 m^{(1, 1)}(x, p_1^1, p_2^1) &= \frac{1}{(1 - \lambda) h(p_1^1, p_2^1|x)} \frac{\partial^2 \psi^{(1, 1)}(x, \mathbf{p})}{\partial p_1^1 \partial p_2^1}  && \quad  \text{for $\lambda < 1$}.
\end{alignat*}
\end{theorem}

These identification results can be visually understood by Figure \ref{fig:nash-comp2}.
For example, $m^{(0, 0)}(x, p_1^0, p_2^0)$ is the mean of $Y_1^{(0, 0)}$ at point C, where the players are at the margin of the actions $D = (0, 0)$ and $(1, 1)$ given that $(P_1^0, P_2^0) = (p_1^0, p_2^0)$ holds.
Thus, a small fluctuation around $(p_1^0, p_2^0)$ can change the treatment status of a certain proportion of the players from $D = (0,0)$ to $(1,1)$ provided $\lambda > 0$, and we  can utilize this local exogenous variation for the identification of the MTR functions.
Note that, if $\lambda = 0$, the identification at point C fails because $D = (0, 0)$ is never chosen at this point.
\begin{comment}
To be more precise, observe that
\begin{align*}
	\psi^{(0, 0)}(x, \mathbf{p})
	&=  E[Y_1^{(0, 0)} | X = x, V \in \mathcal{V}_{\text{uni}}^{(0, 0)}] \cdot \Pr[V \in \mathcal{V}_{\text{uni}}^{(0, 0)} | X = x] \\
	& \quad + \lambda \cdot E[Y_1^{(0, 0)} | X = x, V \in \mathcal{V}_{\text{mul}}] \cdot \Pr[V \in \mathcal{V}_{\text{mul}} | X = x]\\
	&= \int_{p_1^0}^1 \int_{p_2^0}^1 m^{(0, 0)}(x, v_1, v_2) h(v_1, v_2 | x) \mathrm{d}v_1 \mathrm{d}v_2  - (1 - \lambda) \int_{p_1^0}^{p_1^1} \int_{p_2^0}^{p_2^1} m^{(0, 0)}(x, v_1, v_2) h(v_1, v_2 | x) \mathrm{d}v_1 \mathrm{d}v_2,
\end{align*}
where the second equality follows by direct calculation.
Thus, partially differentiating both sides with respect to $p_1^0$ and $p_2^0$ leads to the result for $m^{(0, 0)}(x, p_1^0, p_2^0)$.
\end{comment}

The above identification results are obtained by conditioning the values of all four variables $(P_1^0, P_1^1, P_2^0, P_2^1)$.
Recall that Theorem \ref{thm:unique} required only two out of four values as the conditioning variables, since the point at which identification is achieved is well characterized as the upper left or the lower right corner in the space of $V$.
By contrast, we should fix the values of all four $(P_1^0, P_1^1, P_2^0, P_2^1)$ under multiple equilibria because these four points are required for characterizing the multiple equilibria region.
Nonetheless, the parameter to be identified here, namely $m^{(0, 0)}(x, p_1^0, p_2^0)$, relates to only two values, $p_1^0$ and $p_2^0$, and it is thus irrelevant to $p_1^1$ and $p_2^1$.
As shown in the next section, the computational effort required to estimate $m^{(0, 0)}(x, p_1, p_2)$ and $m^{(1, 1)}(x, p_1, p_2)$ is no greater than that for $m^{(1, 0)}(x, p_1, p_2)$ or $m^{(0, 1)}(x, p_1, p_2)$.

\begin{remark}[\citet{lee2018identifying}]\label{remark:leesalanie}
Here, to contrast the treatment decision model and the identification results in this paper with those of \citet{lee2018identifying}, we review Appendix C of their paper.
Suppose that each player $j$ has payoff $\pi_j^0 = Q_j^0(\mathbf{Z}_j^0) - V_j^0$ when $(D_j, D_{-j}) = (1, 0)$ and payoff $\pi_j^1 = Q_j^1(\mathbf{Z}_j^1) - V_j^1$ when $(D_j, D_{-j}) = (1, 1)$, where $Q_j^0(\mathbf{Z}_j^0)$ and $Q_j^1(\mathbf{Z}_j^1)$ are identifiable functions of IVs, and $V_j^0$ and $V_j^1$ are unobserved variables.
Under strategic complementarity $\pi_j^0 < \pi_j^1$ for both $j =1$ and $2$,\footnote{
	\citet{lee2018identifying} originally assume strategic substitutes, but we focus on the case of strategic complements in line with our setting.
}
if only player $1$ takes the treatment, for example, then we know that $\pi_1^0 > 0 \;\; \text{and} \;\; \pi_2^1  < 0$, i.e., $V_1^0 >  Q_1^0(\mathbf{Z}_1^0) \;\; \text{and} \;\; V_2^1 >  Q_2^1(\mathbf{Z}_2^1) $.
When $\pi_j^0 < 0 < \pi_j^1$ for both $j = 1$ and $2$, both $D = (0, 0)$ and $(1, 1)$ are Nash equilibria.
Now, assume that under the multiple equilibria, $D = (0, 0)$ is selected if and only if $V_3 \le Q_3(\mathbf{Z}_3)$.
Then, all possible treatment decisions are uniquely characterized by the set of threshold crossing rules determined by $\mathbf{Q} = (Q_1^0(\mathbf{Z}_1^0), Q_1^1(\mathbf{Z}_1^1), Q_2^0(\mathbf{Z}_2^0), Q_2^1(\mathbf{Z}_2^1), Q_3(\mathbf{Z}_3))$ and $\mathbf{V} = (V_1^0, V_1^1, V_2^0, V_2^1, V_3)$.
Theorem 3.1 in \citet{lee2018identifying} implies that
\begin{align*}
	E[Y_1^{(d_1, d_2)} | X=x, \mathbf{V} = \mathbf{q}]
	= \frac{\partial_{\mathbf{q}} E[I^{(d_1, d_2)} Y_1 | X = x , \mathbf{Q} = \mathbf{q}]}{\partial_{\mathbf{q}} \Pr[D = (d_1, d_2) | X=x, \mathbf{Q} = \mathbf{q}]},
\end{align*}
where $\partial_{\mathbf{q}} \coloneqq \partial^5 / (\partial q_1^0 \partial q_1^1 \partial q_2^0 \partial q_2^1 \partial q_3)$ denotes the differential operator with $\mathbf{q} = (q_1^0, q_1^1, q_2^0, q_2^1, q_3)$.
It should be noted that this result
requires at least five distinct IVs, which would be a strong requirement in practice.
Moreover, even if there are rich IVs, the resulting estimator based on the above equality will suffer from the curse of dimensionality.
\end{remark}

\begin{comment}
\begin{remark}[\citet{lee2018identifying}]\label{remark:leesalanie}
One may notice that our identification results are similar to those of \citet{lee2018identifying}.
To clarify the difference between our study and theirs, let $\tilde V \coloneqq (\tilde V_1, \tilde V_2, \tilde V_3, \tilde V_4, \tilde V_5)$, where $\tilde V_1 = V_1$, $\tilde V_2 = V_2$, $\tilde V_3 = V_1$, $\tilde V_4 = V_2$, and $\tilde V_5 = \epsilon$.
Further, let $Q \coloneqq (Q_1, Q_2, Q_3, Q_4, Q_5)$, where $Q_1 = P_1^0$, $Q_2 = P_2^0$, $Q_3 = P_1^1$, $Q_4 = P_2^1$, and $Q_5 = \lambda$, and $\partial_q \coloneqq \partial^5/ \prod_{j=1}^5 \partial q_j$ be the differential operator with $q = (q_1, q_2, q_3, q_4, q_5)$.
Then, the treatment choice for each pair can be uniquely determined by a set of threshold crossing rules defined by $\tilde V$ and $Q$.
Theorem 3.1 in \citet{lee2018identifying} would imply that
\begin{align*}
E[Y^{(0, 0)} | X=x, \tilde V = q]
= \frac{\partial_q E[I^{(0, 0)} Y | X=x ,Q = q]}{\partial_q \Pr[D = (0, 0) | X=x, Q = q]}.
\end{align*}
The above result is based on the partial derivatives with respect to all thresholds $q$.
Thus, when directly based on this, the resulting estimator would suffer from the curse of dimensionality.
Another notable difference is the over-identifiability of the MTR functions, as mentioned in the next remark.
\end{remark}
\end{comment}

\begin{remark}[Over-identification]\label{remark:over}
	We can achieve over-identification for $m^{(0, 0)}(x, p_1, p_2)$ and $m^{(1, 1)}(x, p_1, p_2)$.\footnote{
		Note that, if the treatment decisions are substitutes, the over-identification results for $m^{(0, 0)}(x, p_1, p_2)$ and $m^{(1, 1)}(x, p_1, p_2)$ do not hold; rather, similar results hold for $m^{(1, 0)}(x, p_1, p_2)$ and $m^{(0, 1)}(x, p_1, p_2)$.
	}
	For example, assume that some $(p_1, p_2) \in \text{supp}[P_1^0, P_2^0 | X = x, D = (0, 0)] \cap \text{supp}[P_1^1, P_2^1 | X = x, D = (0, 0)]$ exists.
	Then, if $0 < \lambda < 1$, we can identify $m^{(0, 0)}(x, p_1, p_2)$ in at least two ways:
	\begin{align*}
		m^{(0, 0)}(x, p_1, p_2)
		&= \frac{1}{\lambda h(p_1, p_2 | x)} \frac{\partial^2 \psi^{(0, 0)}(x, p_1, p_1^1, p_2, p_2^1)}{\partial p_1 \partial p_2} \\
		&= - \frac{1}{(1 -\lambda) h(p_1, p_2|x)} \frac{\partial^2 \psi^{(0, 0)}(x, p_1^0, p_1, p_2^0, p_2)}{\partial p_1 \partial p_2},
	\end{align*}
	where $(p_1^1, p_2^1)$ and $(p_1^0, p_2^0)$ can be any value as long as they are consistent with the value of $(p_1, p_2)$.

	This over-identification result enables us to improve the estimation efficiency (see Remark \ref{remark:overidentification:est}).
	In addition, it gives us a testable implication for the validity of Assumption \ref{as:multiple} since there is no guarantee that the over-identification holds without Assumption \ref{as:multiple}.
\end{remark}

\subsection{Identification of the marginal treatment effects} \label{subsec:MTE}

Given Theorems \ref{thm:unique} and \ref{thm:multiple}, the identification of MTE is rather straightforward.
Here, it is sufficient to discuss only the identification of the direct MTE and total MTE since $\text{MTE}_{\text{indirect}}^{(0)}(x, p_1, p_2) = \text{MTE}_{\text{total}}(x, p_1, p_2) - \text{MTE}_{\text{direct}}^{(1)}(x, p_1, p_2)$ and $\text{MTE}_{\text{indirect}}^{(1)}(x, p_1, p_2) = \text{MTE}_{\text{total}}(x, p_1, p_2) - \text{MTE}_{\text{direct}}^{(0)}(x, p_1, p_2)$.
Moreover, for the direct MTE, we focus only on the case of $d_2 = 0$ for exposition purposes.

From Theorems \ref{thm:unique} and \ref{thm:multiple}, assuming that $\lambda < 1$, we can write the direct MTE as
\begin{align*}
	 \text{MTE}_{\text{direct}}^{(0)}(x, p_1^0, p_2^1)
	&=  -\frac{1}{h(p^0_1, p^1_2 | x)} \frac{\partial^2}{\partial p^0_1 \partial p^1_2} \left(  \psi^{(1,0)}(x, p^0_1, p^1_2) + \frac{\psi^{(0, 0)} (x, \mathbf{p})}{1 - \lambda} \right),
\end{align*}
which measures the average direct treatment effect on player 1's outcome for the individuals who are at the margin of the actions $D = (1,0)$, $(1,1)$, and $(0,0)$.
Recall that the value of $\psi^{(1,0)}(x, p^0_1, p^1_2)$ cannot be identified from the data if $(p^0_1, p^1_2) \not\in \text{supp}[P_1^0, P_2^1 | X = x, D = (1, 0)]$.
Similarly, if $\mathbf{p} \not\in \text{supp}[\mathbf{P} | X = x, D = (0, 0)]$, $\psi^{(0, 0)} (x, \mathbf{p})$ is not obtained.
However, as $p_1^1$ and $p_2^0$ are irrelevant to the value of $m^{(0, 0)}(x, p^0_1, p^1_2)$, they can take any values as long as $ (p^0_1, p^1_2) \in \text{supp}[P_1^0, P_2^1|X = x, D = (0, 0)]$.
Consequently, the identification of this MTE parameter requires that the support of $(P_1^0, P_2^1)$ for players satisfying $D = (1, 0)$ overlaps with that of $(P_1^0, P_2^1)$ for players satisfying $D = (0, 0)$ conditional on $X = x$.
In practice, the observations for which $(P_1^0, P_2^1)$ is not contained in the common support should be dropped when the MTE is computed.

For identification of the total MTE, assuming that $\lambda > 0$, we have
\begin{align*}
	\text{MTE}_{\text{total}}(x, p_1^0, p_2^0)
	&= - \frac{1}{\lambda h(p^0_1, p^0_2 | x)}\frac{\partial^2}{\partial p^0_1 \partial p^0_2}  \left( \psi^{(1, 1)}(x, \mathbf{p}) +  \psi^{(0, 0)}(x, \mathbf{p}) \right).
\end{align*}
In this case, the conditional support of $(P_1^0, P_2^0)$ for $D = (1, 1)$ and that for $D = (0,0)$ must overlap.

\begin{remark}[Partial identification of MTE without equilibrium selection]\label{remark:partial}
	As shown above, the MTR functions for $D=(1,0)$ and $(0,1)$ can be point-identified without using the equilibrium selection assumption.
	When no assumptions are imposed on the equilibrium selection, $\lambda$ can take any value in $[0,1]$ non-uniformly (possibly depending on other random variables) over the multiple equilibria region.
	Thus, as can be seen from Theorem \ref{thm:multiple}, the identified set of the MTR functions for $D = (0, 0)$ and $(1, 1)$ are typically unbounded.
	However, since $\lambda$ is non-negative, we can still identify the sign of these MTR functions.
	Therefore, for example, if $m^{(1,0)}(x, p_1, p_2)$ is point-identified as a positive value and the sign of $m^{(0,0)}(x, p_1, p_2)$ is identified as non-positive, then we can obtain an informative lower bound of $\text{MTE}_{\text{direct}}^{(0)}(x, p_1, p_2)$ by the inequality $\text{MTE}_{\text{direct}}^{(0)}(x, p_1, p_2) \ge m^{(1,0)}(x, p_1, p_2)$.
\end{remark}

\section{Estimation and Asymptotics} \label{sec:estimation}

In this section, we propose a two-step semiparametric procedure for estimating the MTE parameters given the data $\{\{(Y_{ji}, D_{ji}, W_{ji})\}_{j=1}^2\}_{i=1}^n$ are observed.
In the subsequent analysis, for empirical tractability, we consider the following parametric treatment decision model:
\begin{assumption}
	\hfil
	\label{as:gamemodel}
	\begin{enumerate}[(i)]
		\item Each player $j$ decides his/her treatment by $D_{ji} = \mathbf{1}\{ W_{ji}^\top \gamma_0 + D_{-j, i} \cdot \eta(W_{ji}^\top \gamma_1) \ge \varepsilon_{ji} \}$, where $\gamma = (\gamma_0^\top, \gamma_1^\top)^\top \in \mathbb{R}^{2\mathrm{dim}(W)}$ is a vector of unknown parameters and $\eta$ is a known positive function.

		\item The random variables $\varepsilon = (\varepsilon_1, \varepsilon_2)$ are independent of $W$ and continuously distributed

		with known marginal CDFs $F_{\varepsilon_1}$ and $F_{\varepsilon_2}$, respectively.
		Their joint distribution is represented by a copula $H_\rho$ such that
		$\Pr(\varepsilon_1 \le a_1, \varepsilon_2 \le a_2) = H_\rho(F_{\varepsilon_1}(a_1), F_{\varepsilon_2}(a_2))$, where $\rho \in [\underline{c}_\rho, \bar c_\rho]$ is a correlation parameter between $\varepsilon_1$ and $\varepsilon_2$, and $\underline{c}_\rho$ and $\bar c_\rho$ are real numbers whose values depend on the choice of the copula function.
		The copula $H_\rho$ has a density function $h_\rho$.
	\end{enumerate}
\end{assumption}

In Assumption \ref{as:gamemodel}(i), we assume a strictly positive interaction effect $\eta > 0$ to ensure strategic complementarity as in Assumption \ref{as:complement}.
In the literature, it has typically been assumed that the interaction effect is a constant value.
However, such a model specification cannot be adopted here because the constancy of the interaction effect implies that the conditional support of $P_j^1$ given $P_j^0$ degenerates to a singleton.
We also assume that the coefficients $\gamma$ are common to both players for simplicity.
In Assumption \ref{as:gamemodel}(ii), we assume that the error terms $(\varepsilon_1, \varepsilon_2)$ are
continuously distributed and that their joint distribution function is known up to a parameter $\rho$.
A typical example satisfying this assumption is a standard bivariate normal distribution.
In this case, the copula $H_\rho$ corresponds to the Gaussian copula $H_\rho(v_1, v_2) = \Phi_2(\Phi^{-1}(v_1), \Phi^{-1}(v_2); \rho)$, where $\Phi_2(\cdot, \cdot; \rho)$ is the standard bivariate normal CDF with correlation $\rho$, and $\Phi(\cdot)$ denotes the standard normal CDF.
As another example, \citetmain{aradillas2010semiparametric} uses the Farlie--Gumbel--Morgenstern (FGM) copula $H_\rho(v_1, v_2) = v_1v_2[1 + \rho(1 - v_1)(1 - v_2)]$ in a context similar to ours.
Importantly, as recently shown by \citetmain{khan2018information}, when the marginal CDFs of $(\varepsilon_1, \varepsilon_2)$ are unknown, it is generally impossible to estimate the interaction effect at the parametric rate $n^{-1/2}$.
For completeness, in Appendix \ref{appendix:game}, we provide primitive conditions to verify Assumptions \ref{as:PV} and \ref{as:multiple}(ii) under the setup in Assumption \ref{as:gamemodel} by extending the identification result for the generalized bivariate probit model in \citemain{han2017identification}, where our model can be seen as a game-econometric extension of theirs.

To facilitate the analysis, for the potential outcomes as well, we assume a linear parametric model.
The following assumption is stronger than necessary for identification of the MTE, but is a popular setup in the literature (e.g., \citealpmain{carneiro2009estimating}; \citealpmain{felfe2018does}).

\begin{assumption}
	\hfil
	\label{as:linear}
	\begin{enumerate}[(i)]
		\item For each player $j$, the potential outcome $Y_{ji}^{(d_j, d_{-j})}$ is determined by $Y_{ji}^{(d_j, d_{-j})} = X_{ji}^\top \beta_j^{(d_j, d_{-j})} + U_{ji}^{(d_j, d_{-j})}$, where $\beta_j^{(d_j, d_{-j})} \in \mathbb{R}^{\mathrm{dim}(X)}$ is a vector of unknown parameters.
		\item The random variables $U^{(d_1,d_2)}$ are independent of $W$ for $(d_1, d_2) \in \{0, 1\}^2$.
	\end{enumerate}
\end{assumption}

\subsection{Two-step estimation} \label{subsec:procedure}

\paragraph{First step: Estimation of the treatment decision game.}

In accordance with \eqref{eq:definition-P}, we write $V_{ji} = F_{\varepsilon_j}(\varepsilon_{ji})$, $P_{ji}^0(\gamma) = F_{\varepsilon_j}(W_{ji}^\top \gamma_0)$, $P_{ji}^1(\gamma) = F_{\varepsilon_j}(W_{ji}^\top \gamma_0 + \eta(W_{ji}^\top \gamma_1))$, and $\mathbf{P}_i(\gamma) = (P_{1i}^0(\gamma), P_{1i}^1(\gamma), P_{2i}^0(\gamma), P_{2i}^1(\gamma))$.
Further, let $\theta = (\gamma^\top, \rho, \lambda)^\top$, and $\theta^* = (\gamma^{* \top}, \rho^*, \lambda^*)^\top$ be the true value of $\theta$.
For a given $\theta$, the conditional probability that the $i$-th pair of players is in the multiple equilibria region is given by
\begin{align*}
	\mathcal{L}_{\text{mul},i}(\theta) \coloneqq H_\rho(P_{1i}^1(\gamma), P_{2i}^1(\gamma)) - H_\rho(P_{1i}^1(\gamma), P_{2i}^0(\gamma)) - H_\rho(P_{1i}^0(\gamma), P_{2i}^1(\gamma)) + H_\rho(P_{1i}^0(\gamma), P_{2i}^0(\gamma)).
\end{align*}
Then, letting the probability that they choose an action $D_i = (d_1, d_2)$ be $\mathcal{L}_i^{(d_1, d_2)}(\theta)$, we have
\begin{equation}\label{eq:probabilities}
\begin{array}{ll}
	\mathcal{L}_i^{(1, 0)}(\theta) = P_{1i}^0(\gamma) - H_\rho(P_{1i}^0(\gamma), P_{2i}^1(\gamma)), &\quad \mathcal{L}_i^{(1, 1)}(\theta) = H_\rho(P_{1i}^1(\gamma),P_{2i}^1(\gamma)) - \lambda \cdot \mathcal{L}_{\text{mul},i}(\theta),\\
	\mathcal{L}_i^{(0, 1)}(\theta) = P_{2i}^0(\gamma) - H_\rho(P_{1i}^1(\gamma), P_{2i}^0(\gamma)), &\quad \mathcal{L}_i^{(0, 0)}(\theta) = 1 - \sum_{(d_1,d_2) \neq (0,0)}\mathcal{L}_i^{(d_1, d_2)}(\theta).
\end{array}
\end{equation}
The ML estimator $\hat \theta_n$ of $\theta^*$ is obtained by maximizing the log-likelihood function $\sum_{i = 1}^n \ell_i(\theta)$ with respect to $\theta$, where $\ell_i(\theta) \coloneqq \sum_{d_1 = 0}^1 \sum_{d_2 = 0}^1 I_i^{(d_1, d_2)} \log \mathcal{L}_i^{(d_1, d_2)}(\theta)$.
Once the estimator $\hat \theta_n$ is obtained, we can estimate $P_{ji}^0 = P_{ji}^0(\gamma^*)$ and $P_{ji}^1 = P_{ji}^1(\gamma^*)$ by $\hat P_{ji}^0 = P_{ji}^0(\hat \gamma_n)$ and $\hat P_{ji}^1 = P_{ji}^1(\hat \gamma_n)$, respectively.
Similarly, the true joint CDF and density of $ (V_1, V_2)$, which we denote by $H = H_{\rho^*}$ and $h = h_{\rho^*}$, can be respectively estimated as $\hat H = H_{ \hat \rho_n}$ and $\hat h = h_{ \hat \rho_n}$.
Moreover, we can estimate the probabilities in \eqref{eq:probabilities} evaluated at the true $\theta^*$, namely, $\mathcal{L}_i^{(d_1, d_2)} = \mathcal{L}_i^{(d_1, d_2)}(\theta^*)$ and $\mathcal{L}_{\text{mul}, i} = \mathcal{L}_{\text{mul},i}(\theta^*) $, by $\hat{\mathcal{L}}_i^{(d_1, d_2)} = \mathcal{L}_i^{(d_1, d_2)}(\hat \theta_n)$ and $\hat{\mathcal{L}}_{\text{mul}, i}  = \mathcal{L}_{\text{mul},i}(\hat \theta_n)$, respectively.

\paragraph{Second step: Estimation of the MTE.}

As discussed in Section \ref{sec:identification}, a variety of treatment effect parameters can be identified.
We here specifically discuss the estimation of $\text{MTE}_{\text{direct}}^{(0)}(x, p^0_1, p^1_2)$.
The estimation of $\text{MTE}_{\text{total}}(x, p^0_1, p^0_2)$ can be performed similarly (see Appendix \ref{subsec:totalMTE}).

Recall that $\text{MTE}_{\text{direct}}^{(0)}(x, p^0_1, p^1_2)$ is the difference of the MTR functions $m^{(1, 0)}(x, p_1^0, p_2^1)$ and $m^{(0, 0)}(x, p_1^0, p_2^1) $.
We first discuss estimating $m^{(1, 0)}(x, p_1^0, p_2^1)$.
Assumption \ref{as:linear} implies that
\begin{align*}
	m^{(1, 0)}(x, p_1^0, p_2^1)
	= x_1^\top \beta_1^{(1, 0)} + E[U_1^{(1, 0)} | V_1 = p_1^0, V_2 = p_2^1].
\end{align*}
Since $U_1^{(1, 0)}$ may be correlated with $(V_1, V_2)$, $E[U_1^{(1, 0)} | V_1, V_2] \neq 0$ in general.
Observe that
\begin{align*}
	E[ I^{(1, 0)} Y_1 | I^{(1, 0)}, X, P_{1}^0, P_{2}^1 ]
	& = I^{(1, 0)} X_{1}^\top \beta_1^{(1, 0)} + I^{(1, 0)} E[ U_{1}^{(1, 0)}| V_{1} \le P_{1}^0, V_{2} > P_{2}^1].
\end{align*}
Letting $g_1^{(1, 0)}(p_{1}^0, p_{2}^1) \coloneqq \int_{0}^{p_{1}^0} \int_{p_{2}^1}^1 E[U_{1}^{(1, 0)}| V_{1} = v_1, V_{2} = v_2] h(v_1, v_2) \mathrm{d}v_1 \mathrm{d}v_2$, the same argument as in the proof of Theorem \ref{thm:unique} can show that $E[ U_{1}^{(1, 0)}| V_{1} \le p_{1}^0, V_{2} > p_{2}^1] = g_1^{(1, 0)}(p_{1}^0, p_{2}^1) / \mathcal{L}^{(1,0)}(p_1^0, p_2^1; \theta^*)$, where $\mathcal{L}^{(1, 0)}(p_1^0, p_2^1; \theta^*) = p_1^0 - H_{\rho^*}(p_1^0, p_2^1)$ as in \eqref{eq:probabilities} and we will use analogous notations such as $\mathcal{L}^{(0, 0)}(\mathbf{p}; \theta^*)$ below.
Then, we have the following semiparametric partially linear regression model:
\begin{align}\label{eq:plm1}
	I^{(1, 0)} Y_1 = I^{(1, 0)} X_{1}^\top \beta_1^{(1, 0)} + T^{(1, 0)} g_1^{(1, 0)}(P_{1}^0, P_{2}^1) + e_{1}^{(1, 0)},
\end{align}
where $T^{(1, 0)} \coloneqq I^{(1, 0)} / \mathcal{L}^{(1, 0)}$, and $E[e_{1}^{(1, 0)} | I^{(1, 0)}, X, P_{1}^0, P_{2}^1] = 0$ by construction.
Here, note that since $\mathcal{L}^{(1, 0)}$ may take an arbitrarily small value close to zero, the weight term $T^{(1, 0)}$ can be extremely large for some observations.
To deal with this issue, we exclude such observations from the analysis.
Specifically, we introduce a non-negative smoothed indicator function $\tau_\varpi(\mathcal{L}^{(1, 0)})$  such that (i) $\tau_\varpi:[0, 1] \to [0, 1]$ is a non-decreasing function and $\tau_\varpi(a) = 0$ if and only if $a < \varpi$ for a small constant $\varpi > 0$, and (ii) $\tau_\varpi$ is continuously differentiable with a bounded derivative.
Multiplying both sides of \eqref{eq:plm1} by $\tau_\varpi(\mathcal{L}^{(1, 0)})$ gives
\begin{align}\label{eq:plm2}
	\tilde I^{(1, 0)} Y_1 = \tilde I^{(1, 0)} X_{1}^\top \beta_1^{(1, 0)} + \tilde T^{(1, 0)} g_1^{(1, 0)}(P_{1}^0, P_{2}^1) + \tilde e_{1}^{(1, 0)},
\end{align}
where $\tilde I^{(1, 0)} \coloneqq \tau_\varpi(\mathcal{L}^{(1, 0)}) I^{(1, 0)}$, $\tilde T^{(1, 0)} \coloneqq \tilde I^{(1, 0)} / \mathcal{L}^{(1,0)}$, and $\tilde e_1^{(1, 0)} \coloneqq \tau_\varpi(\mathcal{L}^{(1, 0)}) e_{1}^{(1, 0)}$.
Since $\mathcal{L}^{(1, 0)}$ is a function of $P_1^0$ and $P_2^1$, $E[\tilde e_{1}^{(1, 0)} | I^{(1, 0)}, X, P_{1}^0, P_{2}^1] = 0$ holds.

Based on \eqref{eq:plm2}, we consider estimating $\beta_1^{(1, 0)}$ and $g_1^{(1,0)}$ using the series (sieve) method.
Let $b_K(\cdot, \cdot) = (b_{1K}(\cdot, \cdot), \dots, b_{KK}(\cdot, \cdot))^\top$ be a $K \times 1$ vector of bivariate basis functions.
We assume that $g_1^{(1, 0)}$ can be well approximated by a linear combination of the basis functions: $g_1^{(1, 0)}(\cdot, \cdot) \approx b_K(\cdot, \cdot)^\top \alpha_1^{(1,0)}$ for some coefficient vector $\alpha_1^{(1, 0)}$ with sufficiently large $K$.
Then, letting $\hat I_i^{(1, 0)} \coloneqq \tau_\varpi(\hat{\mathcal{L}}_i^{(1, 0)}) I_i^{(1, 0)}$ and $\hat T_i^{(1, 0)} \coloneqq \hat I_i^{(1, 0)} / \hat{\mathcal{L}}_i^{(1,0)}$, we have
\begin{align*}
	\hat I_i^{(1, 0)} Y_{1i} \approx \hat I_i^{(1, 0)}  X_{1i}^\top \beta_1^{(1, 0)} + \hat T_i^{(1, 0)} b_K(\hat P_{1i}^0, \hat P_{2i}^1)^\top \alpha_1^{(1,0)} + \tilde e_{1i}^{(1, 0)}.
\end{align*}
Let $(\hat \beta_{1n}^{(1,0)}, \hat \alpha_{1n}^{(1,0)})$ be the least squares (LS) estimator of $(\beta_1^{(1,0)}, \alpha_1^{(1,0)})$ obtained by regressing $\hat I_i^{(1, 0)} Y_{1i}$ on $(\hat I_i^{(1, 0)}  X_{1i}, \hat T_i^{(1, 0)} b_K(\hat P_{1i}^0, \hat P_{2i}^1))$.
Then, the estimator of $g_1^{(1, 0)}(p_1^0, p_2^1)$ can be obtained by $\hat g_1^{(1, 0)}(p_1^0, p_2^1) \coloneqq b_K(p_1^0, p_2^1)^\top \hat \alpha_{1n}^{(1, 0)}$.
For the estimation of $E[U_1^{(1, 0)} | V_1 = p_1^0, V_2 = p_2^1]$, applying the same argument as in Theorem \ref{thm:unique} yields
\begin{align*}
	E[U_1^{(1, 0)} | V_1 = p_1^0, V_2 = p_2^1] =  - \frac{1}{h(p_1^0, p_2^1)} \frac{\partial^2 g_1^{(1, 0)}(p_1^0, p_2^1)}{\partial p_1^0 \partial p_2^1}.
\end{align*}
Replacing the unknown parameters by their estimators, we can estimate it by
\begin{align}\label{eq:EU10}
\begin{split}
	\hat E_n[U_1^{(1, 0)} | V_1 = p_1^0, V_2 = p_2^1]
	& \coloneqq - \frac{1}{\hat h(p_1^0, p_2^1)} \frac{\partial^2 \hat g_1^{(1, 0)}(p_1^0, p_2^1)}{\partial p_1^0 \partial p_2^1} = - \frac{\ddot{b}_K(p_1^0, p_2^1)^\top \hat \alpha_{1n}^{(1, 0)}}{\hat h(p_1^0, p_2^1)} ,
\end{split}
\end{align}
where $\ddot{b}_K(p_1, p_2) \coloneqq \partial^2  b_K(p_1, p_2) /(\partial p_1 \partial p_2)$.
Finally, we can estimate $m^{(1, 0)}(x, p_1^0, p_2^1)$ by
\begin{align}\label{eq:m10}
	\hat m^{(1, 0)}(x, p_1^0, p_2^1) \coloneqq x_1^\top \hat \beta_{1n}^{(1, 0)} + \hat E_n[U_1^{(1, 0)} | V_1 = p_1^0, V_2 = p_2^1].
\end{align}

Next, we describe the estimation of the MTR function $m^{(0, 0)}(x, p_1^0, p_2^1)$.
As above, we have
\begin{align*}
	m^{(0, 0)}(x, p_1^0, p_2^1)
	& = x_1^\top \beta_1^{(0, 0)} + E[U_1^{(0, 0)} | V_1 = p_1^0, V_2 = p_2^1]
\end{align*}
under Assumption \ref{as:linear}.
Observe that
\begin{align*}
	E[ I^{(0, 0)} Y_{1} | I^{(0, 0)}, X, \mathbf{P}]
	& = I^{(0, 0)} X_{1}^\top \beta_1^{(0, 0)} + I^{(0, 0)} E[U_{1}^{(0, 0)} | I^{(0, 0)} = 1, \mathbf{P}].
\end{align*}
Further, if $E[U_1^{(0, 0)}| V_1 = v_1, V_2 = v_2] h(v_1, v_2)$ is continuous in $(v_1, v_2)$,
it turns out that there exist the bivariate real-valued functions $g_{11}^{(0, 0)}$, $g_{12}^{(0, 0)}$, $g_{13}^{(0, 0)}$, and $g_{14}^{(0, 0)}$ satisfying
\begin{align}\label{eq:additive}
\begin{split}
	 & E[U_{1}^{(0, 0)}| I^{(0, 0)} = 1, \mathbf{P} = \mathbf{p}]
	= \frac{g_{11}^{(0, 0)}(p_{1}^0, p_{2}^0) + g_{12}^{(0, 0)}(p_{1}^1, p_{2}^0) + g_{13}^{(0, 0)}(p_{1}^0, p_{2}^1) + g_{14}^{(0, 0)}(p_{1}^1, p_{2}^1)}{\mathcal{L}^{(0, 0)}(\mathbf{p}; \theta^*)}.
\end{split}
\end{align}
Hence, similarly to \eqref{eq:plm2}, we obtain the following semiparametric partially linear additive regression model:
\begin{align}\label{eq:plm00}
\begin{split}
	\tilde I^{(0, 0)} Y_{1}  = \tilde I^{(0, 0)} X_{1}^\top \beta_1^{(0, 0)}
	& + \tilde T^{(0, 0)} g_{11}^{(0, 0)}(P_{1}^0, P_{2}^0) + \tilde T^{(0, 0)} g_{12}^{(0, 0)}(P_{1}^1, P_{2}^0) \\
	& + \tilde T^{(0, 0)} g_{13}^{(0, 0)}(P_{1}^0, P_{2}^1) + \tilde T^{(0, 0)} g_{14}^{(0, 0)}(P_{1}^1, P_{2}^1) + \tilde e_{1}^{(0, 0)},
\end{split}
\end{align}
where $\tilde I^{(0, 0)} \coloneqq \tau_\varpi(\mathcal{L}^{(0, 0)}) I^{(0, 0)}$, $\tilde T^{(0, 0)} \coloneqq \tilde I^{(0, 0)} / \mathcal{L}^{(0,0)}$, and $E[\tilde e_{1}^{(0, 0)} | I^{(0, 0)}, X, \mathbf{P}] = 0$ holds by construction.
Assuming again that each $g_{1l}^{(0,0)}$, $l = 1, \ldots , 4$, can be approximated by $g_{1l}^{(0,0)}(\cdot, \cdot) \approx b_K(\cdot, \cdot)^\top \alpha_{1l}^{(0,0)}$ and replacing $\tilde I^{(0,0)}$, $\tilde T^{(0,0)}$, and $\mathbf{P} = (P_{1}^0, P_{1}^1, P_{2}^0, P_{2}^1)$ with their estimators $\hat I^{(0,0)} \coloneqq \tau_\varpi(\hat{\mathcal{L}}^{(0, 0)}) I^{(0, 0)}$, $\hat T^{(0, 0)} \coloneqq \hat I^{(0, 0)} / \hat{\mathcal{L}}^{(0,0)}$, and $\hat{\mathbf{P}} = (\hat P_{1}^0, \hat P_{1}^1, \hat P_{2}^0, \hat P_{2}^1)$, respectively, $\beta_1^{(0, 0)}$ and $\alpha_{1l}^{(0,0)}$'s can be estimated by an LS regression.\footnote{
	It is possible to use different orders of basis terms to approximate each component of the functions $g_{1l}^{(0,0)}$'s for $l = 1, \ldots, 4$, but we use the same order $K$ for all, for simplicity.
	Also note that the ``locations'' of the functions $g_{1l}^{(0,0)}$'s are not identified without further restrictions.

	To simplify our presentation, we postulate that an appropriate location normalization is made implicitly.
}
Let $\hat \beta_{1n}^{(0, 0)}$ and $\hat \alpha_{1ln}^{(0,0)}$ be the resulting LS estimators.
Then, each $g_{1l}^{(0,0)}(p_1, p_2)$ can be estimated by $\hat g_{1l}^{(0, 0)}(p_1, p_2) \coloneqq b_K(p_1, p_2)^\top \hat \alpha_{1ln}^{(0, 0)}$.
Moreover, Theorem \ref{thm:multiple} and \eqref{eq:additive} imply that
\begin{align*}
	E[U_1^{(0, 0)} | V_1 = p_1^0, V_2 = p_2^1]
	& =\frac{1}{(1 - \lambda) h(p_1^0, p_2^1)} \frac{\partial^2 g_{13}^{(0, 0)}(p_1^0, p_2^1)}{\partial p_1^0 \partial p_2^1},
\end{align*}
and thus we can estimate this by
\begin{align}\label{eq:EU00}
\begin{split}
	\hat E_n[U_1^{(0, 0)} | V_1 = p_1^0, V_2 = p_2^1]
	& \coloneqq \frac{1}{ (1 - \hat \lambda_n )\hat h(p_1^0, p_2^1)} \frac{\partial^2 \hat g_{13}^{(0, 0)}(p_1^0, p_2^1) }{\partial p_1^0 \partial p_2^1} = \frac{\ddot{b}_K(p_1^0, p_2^1)^\top \hat \alpha_{13n}^{(0, 0)}}{ (1 - \hat \lambda_n )\hat h(p_1^0, p_2^1)}.
\end{split}
\end{align}
Consequently, $m^{(0, 0)}(x, p_1^0, p_2^1)$ can be estimated by
\begin{align}\label{eq:m00}
	\hat m^{(0, 0)}(x, p_1^0, p_2^1) \coloneqq x_1^\top \hat \beta_{1n}^{(0, 0)} + \hat E_n[U_1^{(0, 0)} | V_1 = p_1^0, V_2 = p_2^1].
\end{align}
The above estimator corresponds to the identification result at the upper left corner (point A in Figure \ref{fig:nash-comp2}), and the estimation at the other points is analogous.
Finally, $\text{MTE}_{\text{direct}}^{(0)}(x, p^0_1, p^1_2)$ can be estimated by
\begin{align}\label{eq:estim-direct}
	\widehat{\text{MTE}}_{\text{direct}}^{(0)}(x, p^0_1, p^1_2) \coloneqq \hat m^{(1, 0)}(x, p_1^0, p_2^1) - \hat m^{(0, 0)}(x, p_1^0, p_2^1).
\end{align}

\subsection{Asymptotics}\label{subsec:asymptotics}

In this subsection, we derive the asymptotic properties of the proposed estimators.
We mainly investigate the estimator of the MTR function for $D = (1,0)$.
Analogous arguments apply to the other cases.

For subsequent discussion, we introduce the following notations.
Denote $\mathcal{S}^{(1,0)} \coloneqq \text{supp}[P_1^0, P_2^1 | D = (1,0)]$.
For a positive integer $a$, $\mathbf{I}_a$ denotes the $a \times a$ identity matrix.
For positive integers $a_1$ and $a_2$, $\mathbf{0}_{a_1 \times a_2}$ denotes the $a_1 \times a_2$ zero matrix.
For a matrix $A$, $\|A\| = \sqrt{\mathrm{tr}\{A^\top A\}}$ denotes its Frobenius norm, where $\mathrm{tr}\{ \cdot \}$ is the trace.
When $A$ is a square matrix, we use $\chi_\mathrm{max}(A)$ and $\chi_\mathrm{min}(A)$ to denote its largest and smallest eigenvalues, respectively.
We denote a symmetric generalized inverse of a matrix $A$ by $A^{-}$.

\begin{assumption}\label{as:data}
	The data $\{\{(Y_{ji}, D_{ji}, W_{ji})\}_{j=1}^2\}_{i=1}^n$ are independent and identically distributed across $i$.
\end{assumption}

\begin{assumption}
	\hfil
	\label{as:1stage}
	\begin{enumerate}[(i)]
		\item $\|\hat \theta_n - \theta^*\| = O_P(n^{-1/2})$ and $\max_{1 \le i \le n} |\hat P_{ji}^d - P_{ji}^d| = O_P(n^{-1/2})$ for $j = 1, 2$ and $d = 0, 1$.

		\item $H_ \rho(p_1, p_2)$ is Lipschitz continuous with respect to $ \rho$ in the neighborhood of $\rho^*$.
		\item For all $\rho$ in the neighborhood of $\rho^*$, $h_ \rho(p_1, p_2)$ is continuous in $(p_1, p_2)\in [0,1]^2$, bounded away from zero uniformly in $(p_1, p_2) \in \mathcal{S}^{(1,0)}$, and Lipschitz continuous in $\rho$.
	\end{enumerate}
\end{assumption}

Assumption \ref{as:data} is a standard and relatively weak condition for microeconomic applications.
We do not restrict the dependence between the variables for $j$ and those for $-j$.
Assumption \ref{as:1stage}(i) is the $\sqrt{n}$-consistency of the ML estimator.
This is a high-level condition; however, it is standard for parametric ML estimation under mild regularity conditions.
Assumption \ref{as:1stage}(iii) ensures that $\hat h$ is uniformly bounded away from zero and infinity with probability approaching one (w.p.a.1) in conjunction with (i).

We introduce the following notations: $\mathbf{Y}_1^{(1, 0)} = (\tilde I_1^{(1, 0)}Y_{11}, \dots, \tilde I_n^{(1, 0)}Y_{1n})^\top$, $\hat{\mathbf{Y}}_1^{(1, 0)} = (\hat I_1^{(1, 0)}Y_{11}, \dots, \hat I_n^{(1, 0)}Y_{1n})^\top$, $R_{K,i}^{(1,0)} \coloneqq (\tilde I_i^{(1,0)} X_{1i}^\top, \tilde T_i^{(1,0)} b_K (P_{1i}^0, P_{2i}^1)^\top)^\top$, $\hat R_{K,i}^{(1,0)} \coloneqq (\hat I_i^{(1,0)} X_{1i}^\top, \hat T_i^{(1,0)} b_K (\hat P_{1i}^0, \hat P_{2i}^1)^\top)^\top$, $\mathbf{R}_K^{(1, 0)} = (R_{K,1}^{(1,0)}, \dots, R_{K,n}^{(1,0)})^\top$, and $\hat{\mathbf{R}}_K^{(1,0)} = (\hat R_{K,1}^{(1,0)}, \dots, \hat R_{K,n}^{(1,0)})^\top$.
\begin{comment}
$\mathbf{X}_1^{(1, 0)} \coloneqq (I_1^{(1, 0)} X_{11}, \dots, I_n^{(1, 0)} X_{1n})^\top$
define $\mathbf{T}^{(1,0)} = \text{diag}(T_1^{(1,0)}, \dots, T_n^{(1,0)})$, $\hat{\mathbf{T}}^{(1,0)} = \text{diag}(\hat T_1^{(1,0)}, \dots, \hat T_n^{(1,0)})$, and
\begin{align*}
	& \mathbf{b}_K^{(1, 0)}  \coloneqq \left( b_K \left( P_{11}^0, P_{21}^1 \right), \dots, b_K \left( P_{1n}^0, P_{2n}^1 \right) \right)^\top, \qquad
	\hat{\mathbf{b}}_K^{(1, 0)} \coloneqq \left( b_K \left( \hat P_{11}^0, \hat P_{21}^1  \right), \dots, b_K \left( \hat P_{1n}^0, \hat P_{2n}^1 \right) \right)^\top,\\
	& \mathbf{R}_K^{(1, 0)} \coloneqq \left(\mathbf{X}_1^{(1, 0)}, \mathbf{T}^{(1,0)} \mathbf{b}_K^{(1, 0)} \right), \qquad
	\hat{\mathbf{R}}_K^{(1, 0)} \coloneqq \left(\mathbf{X}_1^{(1, 0)}, \hat{\mathbf{T}}^{(1,0)} \hat{\mathbf{b}}_K^{(1, 0)} \right).
\end{align*}
\end{comment}
Then, our estimator of $\delta_1^{(1, 0)} \coloneqq (\beta_1^{(1, 0)\top}, \alpha_1^{(1,0)\top})^\top$ is written as \eqref{eq:fcoef10}:
\begin{align}
	\hat \delta_{1n}^{(1, 0)} = \left( \hat \beta_{1n}^{(1, 0)\top},  \hat \alpha_{1n}^{(1, 0)\top} \right)^\top
	&\coloneqq \left[ \hat{\mathbf{R}}_K^{(1, 0)\top} \hat{\mathbf{R}}_K^{(1, 0)} \right]^{-} \hat{\mathbf{R}}_K^{(1, 0)\top} \hat{\mathbf{Y}}_1^{(1, 0)}, \label{eq:fcoef10}\\
	\tilde \delta_{1n}^{(1, 0)}   = \left( \tilde \beta_{1n}^{(1, 0)\top},  \tilde \alpha_{1n}^{(1, 0)\top} \right)^\top
	& \coloneqq \left[ \mathbf{R}_K^{(1, 0)\top} \mathbf{R}_K^{(1, 0)} \right]^{-} \mathbf{R}_K^{(1, 0)\top} \mathbf{Y}_1^{(1, 0)}, \label{eq:infcoef10}
\end{align}
while the one in \eqref{eq:infcoef10}  is an ``infeasible'' estimator with the true parameters in the first stage being treated as known.
The infeasible estimators of $g_1^{(1, 0)}(p_1^0, p_2^1)$, $E[U_1^{(1, 0)} | V_1 = p_1^0, V_2 = p_2^1] $, and $m^{(1, 0)}(x, p_1^0, p_2^1)$ can be defined similarly, which we denote by $\tilde g_1^{(1, 0)}(p_1^0, p_2^1)$, $\tilde E_n[U_1^{(1, 0)} | V_1 = p_1^0, V_2 = p_2^1] $, and $\tilde m^{(1, 0)}(x, p_1^0, p_2^1)$, respectively.
Furthermore, define $\Psi^{(1, 0)}_K \coloneqq E\left[ R_K^{(1,0)} R_K^{(1,0)\top} \right]$, $\Psi^{(1, 0)}_{nK} \coloneqq \mathbf{R}_{K}^{(1, 0) \top} \mathbf{R}_K^{(1, 0)}/n$, $\hat \Psi^{(1, 0)}_{nK} \coloneqq \hat{\mathbf{R}}_{K}^{(1, 0) \top} \hat{\mathbf{R}}_K^{(1, 0)}/n$, and $\Sigma^{(1, 0)}_K \coloneqq E\left[\left( \tilde e_1^{(1,0)} \right)^2 R_K^{(1,0)} R_K^{(1,0)\top} \right]$.

\begin{assumption}
	\hfil
	\label{as:eigen}
	\begin{enumerate}[(i)]
		\item There exist positive constants $\underline{c}_\Psi$ and $\bar{c}_\Psi$ such that $0 < \underline{c}_\Psi \leq \chi_\mathrm{min}\left( \Psi^{(1, 0)}_K \right) \leq \chi_\mathrm{max}\left( \Psi^{(1, 0)}_K \right) \leq \bar{c}_\Psi < \infty$ uniformly in $K$.
		\item There exist positive constants $\underline{c}_\Sigma$ and $\bar{c}_\Sigma$ such that $0 < \underline{c}_\Sigma \leq \chi_\mathrm{min}\left( \Sigma^{(1, 0)}_K \right) \leq \chi_\mathrm{max}\left( \Sigma^{(1, 0)}_K \right) \leq \bar{c}_\Sigma < \infty$ uniformly in $K$.
	\end{enumerate}
\end{assumption}

\begin{assumption}\label{as:emom}
	$E[ ( e_1^{(1,0)} )^4 | W, D ]$ is bounded.
\end{assumption}

Assumption \ref{as:eigen}(i) ensures the existence of the inverse matrices $[ \Psi^{(1, 0)}_{nK} ]^{-1}$ and $[ \hat \Psi^{(1, 0)}_{nK} ]^{-1}$ w.p.a.1.
Assumption \ref{as:emom} is introduced to conveniently derive the limiting distribution of our estimator.
Note that under this assumption, $E[ ( \tilde e_1^{(1,0)} )^4 | W, D ]$ is also bounded by the definition of $\tau_\varpi$.

To state the next assumption, we introduce the following notation.
For a sufficiently smooth function $g(p_1, p_2)$ and a vector of non-negative integers $\mathbf{a} = (a_1, a_2)$, let $\partial ^{\mathbf{a}} g(p_1, p_2) \coloneqq \partial^{|\mathbf{a}|}g(p_1, p_2) / (\partial ^{a_1}p_1 \partial ^{a_2}p_2)$, where $|\mathbf{a}| = a_1 + a_2$. If $|\mathbf{a}| = 0$, then $\partial ^{\mathbf{a}} g(p_1, p_2) = g(p_1, p_2)$.

\begin{assumption}\label{as:sieve}
	For some integer $s \ge 2$, the functions $g_1^{(1,0)}$ and $b_K$ are at least $s$-times continuously differentiable, and there exists a sequence of vectors $\alpha_1^{(1,0)} \in \mathbb{R}^K$ such that $\sup_{(p_1, p_2) \in [0, 1]^2} | \partial ^{\mathbf{a}} g_1^{(1,0)}(p_1, p_2) -  \partial ^{\mathbf{a}} b_K(p_1, p_2)^\top\alpha_1^{(1,0)} | = O( K^{(|\mathbf{a}| - s)/2} )$.
\end{assumption}
Assumption \ref{as:sieve} clearly restricts the smoothness of $g_1^{(1,0)}$ and choice of the basis functions $b_K$. For example, Lemma 2 in \citetmain{holland2017penalized} shows that when $g_1^{(1,0)}$ is $s$-times continuously differentiable on $[0,1]^2$, Assumption \ref{as:sieve} is satisfied by tensor-product B-splines of order $r$ (degree $r-1$) for $r -2 \geq s$.
For slightly more refined results, when $g_1^{(1,0)}$ is in a H\"{o}lder space of smoothness $s$, B-splines, wavelets, and Cattaneo and Farrell's local polynomial partitioning series can satisfy Assumption \ref{as:sieve} (for details, see \citealpmain{chen2018optimal}, Corollary 3.1, and \citealpmain{cattaneo2013optimal}, Lemma A.2).

\begin{assumption} \label{as:rate}
	As $n \to \infty$,
	(i) $\zeta_0(K)\sqrt{(\log K)/n} \to 0$, and
	(ii) $\zeta_1(K)/\sqrt{n} \to 0$,
	where $\zeta_d(K) \coloneqq \max_{|\mathbf{a}| \leq d} \sup_{(p_1, p_2) \in [0,1]^2 }\| \partial ^{\mathbf{a}} b_K(p_1, p_2) \|$.
\end{assumption}

Assumption \ref{as:rate}(i) is used to prove the convergence of the matrix $\Psi^{(1, 0)}_{nK}$ to $\Psi^{(1, 0)}_K$, and \ref{as:rate}(ii) is additionally introduced to ensure the convergence of $\hat \Psi^{(1, 0)}_{nK}$ to $\Psi^{(1, 0)}_K$.
The bound of $\zeta_0(K)$ is well known for several basis functions, which is typically $\zeta_0(K) = O(\sqrt{K})$ (e.g., \citealpmain{chen2007large}; \citealpmain{belloni2015some}).
By contrast, there are fewer readily available results for the bound of $\zeta_1(K)$.
For example, Cattaneo and Farrell's local polynomial partitioning series satisfies $\zeta_1(K) = O(K)$ (see \citealpmain{cattaneo2013optimal}, Lemma A.1).
In Appendix \ref{subsec:zeta}, we show that the tensor-product B-spline also satisfies $\zeta_1(K) = O(K)$.

Next, consider a generic random vector $Q \in \text{supp}[Q] \subset \mathbb{R}^{\dim(Q)}$ with a finite dimension $\dim(Q)$.
Denote the set of uniformly bounded continuous functions on $\mathcal{D}$ as $\mathcal{C}(\mathcal{D})$.
We define the linear operator $\mathcal{P}_{nK}$ that maps a given function $g \in \mathcal{C}(\text{supp}[Q])$ to the sieve space defined by $b_K$ as follows:
\begin{align*}
\begin{split}
	\mathcal{P}_{nK}g
	= b_K(\cdot, \cdot)^\top \mathbb{S}_K \left[ \Psi^{(1, 0)}_{nK} \right]^{-1} \frac{1}{n} \sum_{i = 1}^n R_{K,i}^{(1,0)} g(Q_i),
\end{split}
\end{align*}
where $\mathbb{S}_K \coloneqq (\mathbf{0}_{K \times \mathrm{dim}(X)}, \mathbf{I}_K)$.
The operator norm of $\mathcal{P}_{nK}$ restricted on $\mathcal{S}^{(1,0)}$ is defined as
\begin{align*}
	\| \mathcal{P}_{nK} \|_{\infty}
	\coloneqq \sup\left\{ \sup_{(p_1, p_2) \in \mathcal{S}^{(1,0)}} \left| \left( \mathcal{P}_{nK} g \right) (p_1, p_2) \right| : g \in \mathcal{C}(\text{supp}[Q]), \sup_{q \in \text{supp}[Q]} |g(q)| = 1 \right\}.
\end{align*}

\begin{assumption} \label{as:opnorm}
	(i) $\| \mathcal{P}_{nK} \|_{\infty} = O_P(1)$.
	(ii) $\zeta_0(K) \zeta_1(K) / \sqrt{n} = O(1)$.
\end{assumption}

\begin{assumption}\label{as:deriv}
	$\sup_{(p_1, p_2) \in \mathcal{S}^{(1,0)}}\left| \partial^{\mathbf{a}} b_K(p_1, p_2)^\top \alpha \right| = O(K^{|\mathbf{a}|/2}) \sup_{(p_1, p_2) \in \mathcal{S}^{(1,0)}} \left|  b_K(p_1, p_2)^\top \alpha \right|$ for any $\alpha \in \mathbb{R}^K$.
\end{assumption}

Assumption \ref{as:opnorm}(i) is a stability condition for the projection operator $ \mathcal{P}_{nK}$.
\citetmain{huang2003local} shows that a similar condition holds true for spline bases under some mild regularity conditions.
In addition, wavelets can satisfy such a condition, as shown in Theorem 5.2 in \citetmain{chen2015optimal}.
For a direct verification of Assumption \ref{as:opnorm}(i) under a particular choice of basis functions, see Appendix \ref{subsec:opnorm}.
Assumption \ref{as:opnorm}(ii) is used to show that the operator norm of the ``feasible'' version of $\mathcal{P}_{nK}$ is also bounded in probability (see Lemma \ref{lem:opnorm}).
In the proof of Corollary 3.1 of \citetmain{chen2018optimal}, it is shown that Assumption \ref{as:deriv} holds for splines and wavelets.

\begin{assumption}
	\hfil
	\label{as:rate2}
	\begin{enumerate}[(i)]
		\item There exist finite constants $c_b > 0$ and $\omega \geq 0$ such that $\| b_K(p_1, p_2) -  b_K(p'_1, p'_2)\| \leq c_b K^\omega \| (p_1,p_2)  - (p'_1, p'_2) \| $ for all $(p_1, p_2), (p_1', p_2') \in [0,1]^2$.
		\item $\zeta_0^2(K)\sqrt{(\log n) /n} \to 0$.
	\end{enumerate}
\end{assumption}

The next theorem establishes the uniform convergence rate of the MTR estimators.

\begin{theorem}\label{thm:cdmean}
	Suppose that Assumptions \ref{as:gamemodel}--\ref{as:1stage},   \ref{as:eigen}(i), and \ref{as:emom}--\ref{as:rate2} hold.
	Then, we have
	\begin{align*}
	\renewcommand{\arraystretch}{2}
	\begin{array}{cl}
	\text{(i)} & \displaystyle \sup_{(p_1, p_2) \in \mathcal{S}^{(1,0)}}\left| \tilde m^{(1, 0)}(x, p_1, p_2)  -  m^{(1, 0)}(x, p_1, p_2)  \right| =  O_P(\zeta_0(K) K \sqrt{\log n / n} ) +  O_P(K^{(2 - s) /2}),\\
	\text{(ii)} & \displaystyle \sup_{(p_1, p_2) \in \mathcal{S}^{(1,0)}}\left| \hat m^{(1, 0)}(x, p_1, p_2)  -  m^{(1, 0)}(x, p_1, p_2)  \right| = O_P(\zeta_0(K) K \sqrt{\log n / n} ) +  O_P(K^{(2 - s) /2}).
	\end{array}
	\renewcommand{\arraystretch}{1}
	\end{align*}
\end{theorem}

The proof of the theorem is straightforward from Lemmas \ref{lem:parametric} and \ref{lem:unifconv2}, and thus it is omitted.
When $\zeta_0(K) \asymp \sqrt{K}$, by choosing $K \asymp (\log n /n)^{-1/(1 + s)}$, we can obtain
\[
\sup_{(p_1, p_2) \in \mathcal{S}^{(1,0)}}\left| \hat m^{(1, 0)}(x, p_1, p_2)  -  m^{(1, 0)}(x, p_1, p_2)  \right| = O_P\left(\left(\log n / n\right)^{\frac{s - 2}{2s + 2}} \right).
\]
The above result implies that our MTR estimator can converge at the optimal uniform rate of \citetmain{stone1982optimal}.
This fact is consistent with the finding recently shown by \citemain{chen2018optimal}.

The next theorem shows that the asymptotic distribution of the feasible estimator for the MTR function $m^{(1, 0)}(x, p_1, p_2)$ is equivalent to that of the infeasible oracle estimator.

\begin{theorem}\label{thm:normal1}
	Suppose that Assumptions \ref{as:gamemodel}--\ref{as:rate2} hold.
	For a given $(p_1^0, p_2^1) \in \mathcal{S}^{(1,0)}$, if in addition $K / \| \ddot{b}_K(p_1^0, p_2^1) \| \to 0$, $\zeta_0(K) \sqrt{K/n} \to 0$, and $\sqrt{n} K^{(2 - s)/2} = O(1)$ hold, then we have
	\begin{align*}
		\text{(i)} & \qquad \frac{\sqrt{n} \left( \tilde m^{(1, 0)}(x, p_1^0, p_2^1)  -  m^{(1, 0)}(x, p_1^0, p_2^1) \right)}{\sigma_K^{(1,0)}(p_1^0, p_2^1)} \overset{d}{\to} N(0, 1),\\
		\text{(ii)} & \qquad \frac{\sqrt{n} \left( \hat m^{(1, 0)}(x, p_1^0, p_2^1)  -  m^{(1, 0)}(x, p_1^0, p_2^1) \right)}{\sigma_K^{(1,0)}(p_1^0, p_2^1)} \overset{d}{\to} N(0, 1),
	\end{align*}
	where
	\begin{align*}
		\sigma_K^{(1,0)}(p_1^0, p_2^1) \coloneqq \frac{\sqrt{ \ddot{b}_K(p_1^0, p_2^1)^\top \mathbb{S}_K \left[ \Psi_K^{(1,0)} \right]^{-1} \Sigma^{(1, 0)}_K \left[ \Psi_K^{(1,0)} \right]^{-1} \mathbb{S}_K^\top \ddot{b}_K(p_1^0, p_2^1)}}{h(p_1^0, p_2^1)}.
	\end{align*}
\end{theorem}
The standard deviation $\sigma_K^{(1,0)}(p_1^0, p_2^1)$ can be easily estimated by a sample analog, replacing the true values and functions with their estimates.
\begin{comment}
Specifically,
\begin{align*}
\hat \sigma_K^{(1,0)}(p_1^0, p_2^1) \coloneqq \frac{\sqrt{ \ddot{b}_K(p_1^0, p_2^1)^\top \mathbb{S}_K \left[ \hat \Psi_{nK}^{(1,0)}\right]^{-1} \hat \Sigma^{(1, 0)}_{nK}  \left[ \hat \Psi_{nK}^{(1,0)}\right]^{-1} \mathbb{S}_K^\top \ddot{b}_K(p_1^0, p_2^1)}}{\hat h(p_1^0, p_2^1) },
\end{align*}
where $\hat \Sigma^{(1, 0)}_{nK} \coloneqq n^{-1} \sum_{i=1}^{n} ( \hat{\hat{e}}_{1i}^{(1,0)} )^2 \hat R_{K,i}^{(1,0)} \hat R_{K,i}^{(1,0)\top}$ with $\hat{\hat{e}}_{1i}^{(1, 0)} \coloneqq \hat I_i^{(1, 0)} Y_{1i} - \hat R_{K,i}^{(1,0)\top} \hat \delta_{1n}^{(1,0)}$.
\end{comment}
\bigskip

Let us briefly discuss the asymptotic properties of the estimators for $D=(0,0)$.
Define $\hat{\mathbf{Y}}_1^{(0,0)} = (\hat I_1^{(0,0)} Y_{11}, \dots, \hat I_n^{(0,0)} Y_{1n})^\top$ and $\hat{\mathbf{R}}_K^{(0, 0)} = (\hat R_{K,1}^{(0,0)}, \dots, \hat R_{K,n}^{(0,0)})^\top$, where
\begin{align*}
	\hat R_{K,i}^{(0,0)} \coloneqq (\hat I_i^{(0,0)} X_{1i}^\top, \hat T_i^{(0,0)} b_K(\hat P_{1i}^0, \hat P_{2i}^0)^\top, \hat T_i^{(0,0)} b_K(\hat P_{1i}^1, \hat P_{2i}^0)^\top, \hat T_i^{(0,0)} b_K(\hat P_{1i}^0, \hat P_{2i}^1)^\top, \hat T_i^{(0,0)} b_K(\hat P_{1i}^1, \hat P_{2i}^1)^\top)^\top.
\end{align*}
\begin{comment}
Define $\mathbf{Y}_1^{(0, 0)}$, $\mathbf{X}_1^{(0, 0)}$, $\mathbf{T}^{(0, 0)}$, and $\hat{\mathbf{T}}^{(0, 0)}$ analogously as above.
Further, let
\begin{align*}
	& \hat{\mathbf{b}}_{1K}^{(0, 0)} \coloneqq \left( b_K \left( \hat P_{11}^0, \hat P_{21}^0 \right), \dots, b_K \left( \hat P_{1n}^0, \hat P_{2n}^0 \right) \right)^\top, \qquad
	\hat{\mathbf{b}}_{2K}^{(0, 0)} \coloneqq \left( b_K \left( \hat P_{11}^1, \hat P_{21}^0 \right), \dots, b_K \left( \hat P_{1n}^1, \hat P_{2n}^0 \right) \right)^\top,\\
	& \hat{\mathbf{b}}_{3K}^{(0, 0)} \coloneqq \left( b_K \left( \hat P_{11}^0, \hat P_{21}^1 \right), \dots, b_K \left( \hat P_{1n}^0, \hat P_{2n}^1 \right) \right)^\top, \qquad
	\hat{\mathbf{b}}_{4K}^{(0, 0)}  \coloneqq \left( b_K \left( \hat P_{11}^1, \hat P_{21}^1 \right), \dots, b_K \left( \hat P_{1n}^1, \hat P_{2n}^1 \right) \right)^\top,\\
	& \hat{\mathbf{b}}_K^{(0, 0)} \coloneqq \left( \hat{\mathbf{b}}_{1K}^{(0, 0)}, \hat{\mathbf{b}}_{2K}^{(0, 0)}, \hat{\mathbf{b}}_{3K}^{(0, 0)}, \hat{\mathbf{b}}_{4K}^{(0, 0)} \right), \qquad
	\hat{\mathbf{R}}_K^{(0, 0)} \coloneqq \left(\mathbf{X}_1^{(0, 0)}, \hat{\mathbf{T}}^{(0,0)} \hat{\mathbf{b}}_K^{(0, 0)} \right).
\end{align*}
\end{comment}
The estimator of $\delta_1^{(0, 0)} = ( \beta_1^{(0, 0)\top}, \alpha_{11}^{(0,0)\top}, \ldots , \alpha_{14}^{(0,0)\top} )^\top$ can be written as
\begin{align}\label{eq:coef00}
	\hat \delta_{1n}^{(0, 0)} = \left( \hat \beta_{1n}^{(0, 0)\top},  \hat \alpha_{11n}^{(0, 0)\top}, \ldots, \hat \alpha_{14n}^{(0, 0)\top} \right)^\top
	&\coloneqq \left[ \hat{\mathbf{R}}_K^{(0, 0)\top} \hat{\mathbf{R}}_K^{(0, 0)} \right]^{-} \hat{\mathbf{R}}_K^{(0, 0)\top} \hat{\mathbf{Y}}_1^{(0, 0)},
\end{align}
and the estimator of the MTR function for $D=(0,0)$ can be obtained by \eqref{eq:m00}.
Under conditions similar to those in Theorem \ref{thm:normal1}, we can show the following asymptotic normality result:
\begin{align*}
	\frac{\sqrt{n} \left( \hat m^{(0, 0)}(x, p_1^0, p_2^1)  -  m^{(0, 0)}(x, p_1^0, p_2^1) \right)}{\sigma_{K}^{(0,0)}(p_1^0, p_2^1)} \overset{d}{\to} N(0, 1),
\end{align*}
where
\begin{align*}
	\sigma_{K}^{(0,0)}(p_1^0, p_2^1) &\coloneqq \frac{\sqrt{ \ddot{b}_K(p_1^0, p_2^1)^\top \mathbb{S}_{3K} \left[ \Psi_K^{(0,0)}\right]^{-1} \Sigma^{(0, 0)}_K  \left[ \Psi_K^{(0,0)}\right]^{-1} \mathbb{S}_{3K}^\top \ddot{b}_K(p_1^0, p_2^1)}}{(1 - \lambda) h(p_1^0, p_2^1) },
\end{align*}
and $\mathbb{S}_{3K} \coloneqq (\mathbf{0}_{K \times \mathrm{dim}(X)}, \mathbf{0}_{K \times 2K}, \mathbf{I}_K,  \mathbf{0}_{K \times K})$.
The definitions of the matrices $\Psi_K^{(0,0)}$ and $\Sigma_K^{(0,0)}$ should be clear from the context.

Finally, the limiting distribution of $\widehat{\text{MTE}}_{\text{direct}}^{(0)}(x, p^0_1, p^1_2)$ can be characterized as follows:
\begin{align*}
	\frac{\sqrt{n} \left( \widehat{\text{MTE}}_{\text{direct}}^{(0)}(x, p^0_1, p^1_2)  -  \text{MTE}_{\text{direct}}^{(0)}(x, p^0_1, p^1_2) \right)}{\sqrt{\left[\sigma_K^{(1,0)}(p_1^0, p_2^1)\right]^2 + \left[\sigma_{K}^{(0,0)}(p_1^0, p_2^1)\right]^2}} \overset{d}{\to} N(0, 1).
\end{align*}
The limiting distributions of the other MTE estimators can be derived similarly.

\begin{remark}[An over-identified estimator]\label{remark:overidentification:est}
	If the intersection of the supports of $(P_1^0, P_2^0)$, $(P_1^1, P_2^0)$, $(P_1^0, P_2^1)$, and $(P_1^1, P_2^1)$ conditional on $D = (0, 0)$ is non-empty, we may estimate the MTR function $m^{(0, 0)}(x, p_1, p_2)$ in four ways, as shown in Theorem \ref{thm:multiple}: $\hat m_l^{(0, 0)}(x, p_1, p_2) \coloneqq x_1^\top \hat \beta_{1n}^{(0, 0)} + \hat E_{ln}[ U_1^{(0, 0)} | V_1 = p_1, V_2 = p_2 ]$ for $l = 1,2,3,4$, where

	$\hat E_{ln} [ U_1^{(0, 0)} | V_1 = p_1, V_2 = p_2 ] \coloneqq \hat \kappa_{ln} \ddot b_K(p_1, p_2)^\top \hat \alpha_{1ln}^{(0, 0)}$ with $\hat \kappa_{1n} \coloneqq 1 / (\hat \lambda_n \hat h(p_1, p_2))$, $\hat \kappa_{2n} = \hat \kappa_{3n} \coloneqq 1/((1 - \hat \lambda_n) \hat h(p_1, p_2))$, and $\hat \kappa_{4n} \coloneqq -1/((1 - \hat \lambda_n) \hat h(p_1, p_2))$.
	Let $w = (w_1, w_2, w_3, w_4)^\top \in [0,1]^4$ be a vector of fixed weights such that $\sum_{l = 1}^4 w_l = 1$.
	Then, $m^{(0, 0)}(x, p_1, p_2)$ can be estimated by the weighted average of the four estimators: $\hat m_{\text{over-id}}^{(0, 0)}(x, p_1, p_2; w) \coloneqq \sum_{l=1}^{4} w_l \hat m_l^{(0, 0)}(x, p_1, p_2)$.
	In the same manner as in the proof of Theorem \ref{thm:normal1}, we can show
	\begin{align*}
		\frac{\sqrt{n} \left( \hat m_{\text{over-id}}^{(0, 0)}(x, p_1, p_2; w) - m^{(0, 0)}(x, p_1, p_2) \right)}{\sigma_{\text{over-id}, K}^{(0, 0)}(p_1, p_2; w) } \stackrel{d}{\longrightarrow} N(0, 1),
	\end{align*}
	where $[ \sigma_{\text{over-id}, K}^{(0, 0)}(p_1, p_2; w) ]^2 = w^\top \Upsilon w$, and $\Upsilon = (\Upsilon_{lm})$ is the $4 \times 4$ symmetric matrix with its $(l, m)$-th element $\Upsilon_{lm} \coloneqq \kappa_l \kappa_m \ddot b_K(p_1, p_2)^\top \mathbb{S}_{lK} [\Psi_K^{(0,0)}]^{-1} \Sigma_K^{(0, 0)} [\Psi_K^{(0,0)}]^{-1} \mathbb{S}_{mK}^\top \ddot b_K(p_1, p_2)$.
	Since the asymptotic variance has a quadratic form, we can easily find the  ``optimal'' weights that attain the minimum asymptotic variance in the class of the weighted average estimators $\hat m_{\text{over-id}}^{(0, 0)}(x, p_1, p_2; w)$ by quadratic programming.
	Although the matrix $\Upsilon$ is unknown in practice, we can estimate the optimal weights by minimizing $w^\top \hat \Upsilon_n w$ with a consistent estimator $\hat \Upsilon_n$ of $\Upsilon$.

\end{remark}

\section{Numerical Illustrations}\label{sec:numeric}

\subsection{Summary of the Monte Carlo simulation results}\label{subsec:MC}

In order to evaluate the finite sample performance of our MTE estimator, we conduct a set of Monte Carlo simulations based on the over-identified estimator in Remark \ref{remark:overidentification:est}.
To save space,  we here briefly  report the summary of the results, and detailed information about the designs and results of the experiments are provided in Appendix \ref{sec:MC}.

The simulation results provide us with several practical guidelines as follows.
First, when comparing bivariate power series and tensor-product B-splines for the basis functions, we find that the power-series-based estimator achieves smaller RMSE (root mean squared error) than the B-splines-based estimator, especially when the sample size is not large.
For this result, note that the power series is a globally-supported series (unlike the B-splines), and thus the power-series estimator tends to have smaller variance at the cost of lower flexibility than the B-splines estimator.
Theoretically speaking, B-splines estimators have a better approximation property than power-series estimators, and in fact the latter may not attain the optimal convergence rate.
Hence, the above result may be due to the chosen sample sizes in this simulation setup.
Second, in terms of RMSE, the MTE parameters can be more precisely estimated by ridge regression with a small regularization parameter.\footnote{
	To further improve the results, one may consider a more sophisticated penalization scheme, however it is beyond the scope of this analysis.
}
Remarkably, using the ridge regression often reduces the bias for the B-splines estimator.
Since the tensor-product B-splines involves a large number of basis terms, which causes a severe multicollinearity problem for the non-regularized estimators, there is a higher risk of generating extremely outlying estimates compared with the regularized estimators.
Thus, it would be desirable to employ a regularized regression in practical situations with moderate sample size.
Note that adding a sufficiently small regularization factor does not alter our asymptotic results.

\subsection{Empirical illustration}\label{subsec:empirical}

As an empirical illustration, we estimate the impacts of risky behaviors such as smoking and drinking alcohol on adolescents' academic performance.
The empirical analysis is performed on the National Longitudinal Study of Adolescent Health (Add Health) data that provide nationally representative information on 7--12th graders in the U.S.
The in-school survey was conducted during the 1994--1995 school year.
The survey elicits information on variables such as the social and demographic characteristics of the respondents, education levels and occupations of their parents, and friendship connections.

In this analysis, similarly to \citemain{card2013peer}, we focus on the interactions among pairs of closest opposite-gender friends.
In the survey, each respondent was asked to list up to five friends of each gender in order from best to 5th best friend.
We first exclude missing nomination data (caused by non-response or because the nominee was not included in the Add Health data) and identify each respondent's closest opposite-gender friend in the available dataset.
Since a student's closest friend's friends might not include the student him/herself, we consider that the pair is formed only when his/her best friend nominates him/her either in the first place or second place, except for missing nominees.
This procedure results in 7,631 opposite-gender pairs of students.

The treatment variable ($D_j$) is defined as the student's participation in risky behaviors including smoking, drinking, truancy, and fighting, and the outcome variable ($Y_j$) is the natural log of the student's GPA score.
Table \ref{table:def} in Appendix \ref{sec:supp} summarizes the explanatory variables used ($X_j$) and their definitions.
For the choice of IVs ($Z_j$), we suppose that the delinquent activities are often group activities.
Indeed, \citemain{lee2014binary} find that adolescents' smoking behavior is significantly influenced by the characteristics of their friends, including non-best friends.
Then, we construct $Z_j$ based on the average of the selected explanatory variables over the 3rd--5th best friends of both genders; note that, at the same time, we are assuming that the 3rd--5th best friends do not have direct influence on the student's own GPA.
After excluding observations with missing values for $(D,X,Z)$, the size of the resulting sample used in the estimation of the treatment decision model is 6,053.
When estimating the MTR functions, observations with missing values for $Y_j$ are further excluded.

\begin{comment}
Table \ref{table:d_dist} presents the empirical distribution of $(D_1, D_2)$, where male students are labeled as ``player 1'' and female students as ``player 2''.
As expected, the case of $D = (0,0)$ represents the largest share of our sample.
There is an interesting asymmetry between the cases of $D = (1,0)$ and $D = (0,1)$; that is, the number of males who solely participate in delinquent behaviors is significantly larger than the number of such females.

\begingroup
\renewcommand{\arraystretch}{1.2}
\begin{table}[!h]
	\caption{Distribution of $(D_1, D_2)$}
	\begin{center}
		\small{\begin{tabular}{|cc|ccc|}
				\hline
				&  & \multicolumn{2}{c}{\underline{Male}}  & \\
				&                 & $D_1 = 0$ & $D_1 = 1$ & Total\\
				\hline
				\multirow{2}{*}{\underline{Female}}& $D_2 = 0$ & 2,707       & 1,351       & 4,058\\
				& $D_2 = 1$ & 736          & 1,259       & 1,995\\
				& Total          & 3,443       & 2,610       & 6,053\\
				\hline
		\end{tabular}}
		\label{table:d_dist}
	\end{center}
\end{table}
\endgroup
\end{comment}

To save space, the detailed estimation results of the treatment decision model are omitted here and are provided in Table \ref{table:1ststage} in Appendix \ref{sec:supp}.
A noteworthy finding is that about 80\% of the pairs choose $D = (0,0)$ in the multiple equilibria situation.
Figure \ref{fig:Phist} in Appendix \ref{sec:supp} shows the histograms of the estimated $(P_j^0, P_j^1)$ for $j=1$ (male) and $2$ (female).
These figures suggest that the overlapping support condition can hardly hold outside the interval $[0.2, 0.8]$, implying that MTE estimates might not be reliable when $V_j$ is outside this interval.

Now, we present our main results of estimating the MTE parameters.
Since our sample size is not very large, following the suggestion from the Monte Carlo results, we employ the over-identified estimator in Remark \ref{remark:overidentification:est} with a third-order bivariate power series for the basis function and use ridge regression for the parameter estimation with penalty  equal to $n^{-1}$.
Figures \ref{fig:directMTE} and \ref{fig:totalMTE} summarize the estimated direct MTE and total MTE, respectively.
For the value of $X_j$, the MTEs are evaluated at the median over all observations of men and women altogether.
Figure \ref{fig:directMTE} shows that the direct MTE is significantly negative for female students only.
For male students, the direct impact of the delinquent activities is less significant and almost independent of the female partner's $V$ (partially due to the regularization).
The direct MTE for the female students tend to have a weak negative slope with respect to the male partner's $V$.
That is, when both of the pair have stronger hesitation in participating in risky activities, the direct MTE on the female students' GPA becomes larger.
Interestingly, this phenomenon is not persistent in the total MTE.
Figure \ref{fig:totalMTE} illustrates that the total MTE is significantly negative for both male and female students.
In particular, the total MTE for male students is clearly larger in magnitude than their direct MTE.
This result would imply that male students' academic performances are more easily influenced by their opposite-gender peers' delinquent behavior than females.

\begin{figure}[!ht]
	\centering
	\includegraphics[width = 17cm]{fig_directmte.png}
	\caption{Estimated direct MTE.}
	\label{fig:directMTE}
\end{figure}
\begin{figure}[!ht]
	\centering
	\includegraphics[width = 17cm]{fig_totalmte.png}
	\caption{Estimated total MTE.}
	\label{fig:totalMTE}
\end{figure}

\section{Conclusion} \label{sec:conclusion}

This study developed treatment effect models that admit both the treatment spillover and strategic interaction in the treatment decisions within a pair of economic agents.
We first demonstrated that the interaction in the treatment decisions can be modeled as a binary game of complete information with multiple equilibria.
Assuming an empirically reasonable equilibrium selection rule, we showed that the MTE can be over-identified using an extended version of the LIV method.
Based on our identification results, we proposed a two-step semiparametric series estimator for MTE parameters.
We showed that the proposed MTE estimator is uniformly consistent and asymptotic normal, and also that its limiting distribution is the same as that of the infeasible oracle estimator.

The results of this study suggest several extensions that would be promising to investigate, some of which have been mentioned above.
In addition, it would be worth extending our results to the case of strategic interaction among more than two players.
In the literature on game econometrics, only a few studies address point identification of game models of complete information with more than two players.
The main reason for this is that the characterization of the equilibrium is extremely complicated compared with the case of two-player games; in general, the number of multiple equilibria regions is increasing in the number of players (see, e.g., \citealpmain{soetevent2007discrete}).
Accordingly, the direct applicability of our point-identification approach to such cases is unclear.
It also would be interesting to extend our approach to situations where not only the individual treatment decision but also the pair formation itself is endogenously determined jointly in a simultaneous game, as in \citemain{badev2013discrete}.
Finally, it may be beneficial to develop treatment evaluation techniques under treatment decision games of incomplete information.
These topics are left for future research.

\section*{Acknowledgments}

The authors are grateful to Yoichi Arai, Sukjin Han, Hiroaki Kaido, Shin Kanaya, Hiroyuki Kasahara, Toru Kitagawa, Tatsushi Oka, Ryo Okui, Yuya Sasaki, Yuta Toyama, Takuya Ura, Haiqing Xu, and the participants of AMES, ESAM, EEA-ESEM, IAAE, SETA, and seminars at various institutions for their valuable comments and suggestions.
This work was supported by JSPS KAKENHI Grant Numbers 15K17039, 17K13715, and 19H01473.
This research uses data from Add Health, a program project directed by Kathleen Mullan Harris and designed by J. Richard Udry, Peter S. Bearman, and Kathleen Mullan Harris at the University of North Carolina at Chapel Hill, and funded by grant P01-HD31921 from the Eunice Kennedy Shriver National Institute of Child Health and Human Development, with cooperative funding from 23 other federal agencies and foundations. Special acknowledgment is due Ronald R. Rindfuss and Barbara Entwisle for assistance in the original design. Information on how to obtain the Add Health data files is available on the Add Health website (http://www.cpc.unc.edu/addhealth). No direct support was received from grant P01-HD31921 for this analysis.

\clearpage

\begin{center}
	{\Large Supplementary Appendix for ``Treatment Effect Models with Strategic Interaction in Treatment Decisions''}

	\textbf{(Not for publication)}
	\bigskip

	{\large Tadao Hoshino$^*$ and Takahide Yanagi$^\dagger$}
	\bigskip

	$^*$ School of Political Science and Economics, Waseda University.

	$^\dagger$ Graduate School of Economics, Kyoto University.
\end{center}