Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
64,171 characters · 10 sections · 40 citation commands
Estimating Dyadic Treatment Effects with Unknown Confounders
Dyadic data are ubiquitous in our society. International trade, travels, population flows, military alliances, partnerships between firms, research collaboration, and many others can be represented as dyadic data, where each dyad represents a pair of countries, firms, or individuals, depending on the context. Dyadic data analysis is particularly prevalent in the literature of international trade, where regression-based analysis, the so-called gravity model, serves as a primary analytical approach in these fields since the pioneering work by tinbergen1962shaping (see also, e.g., anderson1979theoretical, anderson2011gravity, head2014gravity and references therein). For reviews of recent econometric literature on dyadic data analysis in general, see, for example, graham2020dyadic,graham2020network.
Despite the popularity of dyadic data, there are only a few causal inference methods tailored specifically for dyadic data analysis, with some exceptions such as baier2009estimating, arpino2017implementing, and nagengast2023staggered. This may be due to the non-standard and complex endogeneity structure often encountered in typical applications of dyadic data. For example, suppose we are interested in the impacts of free trade agreements (FTA) on trade flows between countries. The treatment variable, FTA, should be considered endogenous because both the decision to enter into FTA and the trade outcome should be influenced by each country's economic factors and the economic and political relationship between the countries involved. Thus, if one tries to resolve the endogeneity issue by using the instrumental variables (IV) method, for instance, then he/she needs to prepare at least three different types of IVs: those accounting for confounding factors at the “origin” country, those at the “destination”, and pair-specific factors.\footnote{ In the specific context of estimating the impacts of FTA on trade flows, this intrinsic difficulty of finding appropriate IVs is discussed in depth in baier2007free. } In addition, not limited to international trade, dyadic data typically emerge in observational studies where it is difficult to find exogenous variations that help the researcher alleviate the endogeneity issue.\footnote{ Although there are a certain number of experimental studies utilizing the network structure in the data, they usually consider assigning a treatment to each individual, but not to each pair of individuals. } Even in cases where the treatments are exogenous, dyadic data naturally entail complex dependence of variables arising from pairwise interactions, which pose additional challenges in developing suitable inferential methods.
In this paper, we aim to develop an easy-to-implement causal inference method specialized for dyadic data, where the treatment variable may potentially be endogenous, correlating with unobserved confounders in an arbitrary manner. We suppose a situation where researchers have access to non-experimental data, comprising solely dyad-level treatments and outcomes. We do not rely on any natural-experimental variations commonly exploited in the literature to identify causal effects, such as IVs, unexpected policy changes, policy-induced discontinuities, etc.
Our causal inference method is built on the following two requirements: (i) the dyadic treatments follow an exchangeable distribution, and (ii) only unit-level factors (rather than dyad-level factors) can cause endogeneity in the treatment choice. The first condition is a common requirement in the literature of network statistics concerning the estimation of graphons, where a graphon is a symmetric nonparametric function that generates graphs (i.e., networks) of any sizes (a more formal definition will be given later). As an important special case, if the treatments are independent and identically distributed (IID) across all dyads, they are automatically exchangeable. Since requiring the exchangeability can be restrictive in some empirical situations, researchers should verify it on a case-by-case basis. Nonetheless, it still encompasses a wide range of network formation models considered in the literature, several examples of which are given in the next section. Fully leveraging the exchangeability of treatments enables us to nonparametrically estimate the propensity score at each dyad without using any additional covariate information. The exchangeability plays also an important role in establishing our inference method.
For the second assumption, which limits the sources of endogeneity to unit-level attributes, it should be noted that variables representing (dis)similarities or social/economic distances between two units are not excluded as sources of endogeneity, as long as they can be expressed as functions of unit-level variables. Even variables such as unobserved strategic partnerships between firms do not violate our requirement only if they are consequences from the interaction of unit-level factors. This requirement should be somewhat similar to an assumption commonly made in panel data analysis with endogeneity, which stipulates that unobserved time-invariant individual heterogeneity is the source of endogeneity, while time-varying individual errors are considered pure idiosyncratic noise.
With these assumptions, we study statistical inference methods for several treatment parameters. Among them, the most basic component is the dyadic average treatment effect (DATE), defined as the difference in the conditional means of potential outcomes given a propensity score. We demonstrate that DATE can be estimated by regressing the outcome on the propensity score separately for the treatment and control groups. Given this result, we propose a nonparametric kernel smoothing method for estimating DATE, termed neighborhood kernel smoothing (NKS). The NKS method extends the neighborhood smoothing approach introduced by zhang2017estimating, which was originally developed for estimating probabilities of network edges (i.e., propensity scores in our context). Once DATE is estimated for all dyads, we can then estimate the individual average treatment effect (IATE) and the global average treatment effect (GATE). IATE represents the unit-level average treatment effect obtained by averaging DATEs over the “destinations”, while GATE is simply the average of DATEs over all dyads. For statistical testing, we develop a permutation test for the sharp null hypothesis of no treatement effects, which involves permuting the rows and columns of the treatment network matrix based on the exchangeability assumption (which ensures that the permuted treatment matrix is equally likely as the actual one). Regarding the theoretical properties of our method, under a set of regularity conditions, we derive the rate of convergence of the NKS estimator. Additionally, for the cases when the endogeneity is caused by selection-on-returns, we demonstrate that the permutation test exhibits a desirable size control property.
As an empirical application of our method, we investigate the impact of FTAs between two countries on export and import volumes. Our dataset consists of 37 selected countries and regions in Asia, Oceania, and North America. The outcome variables of interest are the share of export amount to the partner country among total exports to all countries and the share of import similarly defined, both for the year 2021. For the treatment variable, we consider all FTAs related to our sample countries that were active as of 2021. Our empirical results suggest that IATE of FTAs is positive for almost all countries (except Singapore) for both exports and imports, and that enactment of an FTA increases bilateral trade flows by approximately 2% on average. The results of the permutation test reveal $p$-values sufficiently small for both exports and imports, demonstrating that the treatment effects exist significantly.
\paragraph{Literature.}
The key idea in this study is to model the allocations of the dyadic treatments across units as a network generated from a graphon. A graphon is a nonparametric function of node-specific latent variables that determines the probability of link connection. Estimation of graphon has been one of the central themes in network statistics. Prior studies proposed various approaches for estimating graphons, including modularity maximization (bickel2009nonparametric), nonparametric block-model approximation (e.g., olhede2014network, gao2015rate, klopp2017oracle), sorting (e.g., chan2014consistent, yang2014nonparametric), neighborhood smoothing (e.g., zhang2017estimating, su2020network), among many others. In particular, we extend the neighborhood smoothing method of zhang2017estimating to the estimation of propensity score in our setting in combination with nonparametric kernel smoothing.
The use of the neighborhood smoothing method, or its variants, has been gaining some attention in the econometrics literature. auerbach2022identification examined a partially linear regression model with network-induced unobserved individual heterogeneity. To identify individuals with similar networking types, he proposed a neighborhood smoothing method conceptually similar to zhang2017estimating. Another intriguing application is by zeleneev2020identification, who considered estimating a dyadic regression model with nonparametric unobserved heterogeneity. To accommodate dyad-specific fixed effects, he extended the neighborhood smoothing method to a regression framework. As far as we are aware, our study would be the first to apply the neighborhood smoothing method in the context of causal inference with dyadic data.
As mentioned earlier, only a few studies have explored causal inference methods specifically designed for dyadic data. While numerous studies, particularly in the field of international trade, have attempted to address endogeneity issues in identifying the effects of a treatment, many of them rely on some IV-based approaches or panel data analysis with a fixed-effects approach (e.g., baier2007free, kohl2014we, jochmans2022instrumental, among others), which were not necessarily articulated formally within the framework of causal inference (i.e., the potential outcomes framework). As one notable exception, arpino2017implementing considered extending the propensity score matching method to a dyadic data framework to estimate the causal effect of the General Agreement on Tariffs and Trade (GATT) on bilateral trade. They formulated the probability of being a member of GATT as a function not only of economic and demographic covariates of countries but also of the network structure of the world trade network. Another example is nagengast2023staggered, who also examined the causal impacts of trade agreements on trade flows. Considering that trade agreements were enacted at different time points across different pairs of countries, and that the impacts of the agreements are likely heterogeneous across countries, the authors developed a novel causal inference method by combining the gravity model with the staggered difference-in-differences approach. Compared to these studies, our approach offers significant advantages in terms of analytical simplicity and minimal data requirement.
Our paper also contributes to the literature on permutation inference. As previously mentioned, our identification and estimation approach relies on the exchangeable distribution of treatments. Another advantage of exchangeability is that it naturally leads to a permutation inference method based on the permutation of treatments (e.g., pesarin2010permutation,lehmann2022testing,zhang2023randomization). To the best of our knowledge, no prior studies have considered permutation tests in causal inference with dyadic data. In the context of dyadic regression analysis, erikson2017dyadic explored a permutation inference approach similar to our test. They illustrated the usefulness of permutation inference in a real data analysis for international relations but did not investigate its theoretical properties.
\paragraph{Paper organization.}
In Section (ref), we formally describe our setup and treatment parameters to be estimated. In Section (ref), we first present the characterization and estimation of our target parameters. Then, under certain conditions, the rate of convergence of our estimator is established. After that, we present our permutation inference method for testing the sharp null hypothesis and prove its size control property. We report the results of numerical simulations and empirical data analysis on international trade in Sections (ref) and (ref), respectively. Section (ref) concludes the paper. The proofs of technical results and supplementary numerical results are provided in the appendix.
\paragraph{Notation.}
Throughout the paper, $C$ (possibly with subscripts) denotes a global constant that does not depend on other variables, whose value may differ in different contexts. We write a matrix whose $(i,j)$-th element is $A_{ij}$ as $A = (A_{ij})$ and its $i$-th row as $A_{i \cdot}$. For random variables $X$ and $Y$, we write $X \overset{d}{=} Y$ if $X$ and $Y$ follow the same probability distribution, $X \lesssim Y$ if $X = O(Y)$ almost surely (a.s.), and $X \lesssim_{\mathbb{P}} Y$ if $X = O_{\mathbb{P}}(Y)$.
Consider a population of interest composed of individuals, households, firms, or countries, depending on the context, that are interacting to each other pairwisely. We can observe $n$ units $[n] = \{ 1, 2, \dots, n \}$ drawn from the population, and we treat each pair of units (i.e., dyad) as one observation. We are interested in estimating the causal effect of a dyadic treatment variable on a dyadic outcome variable. Let $A_{ij} \in \{ 0, 1 \}$ and $Y_{ij} \in \mathbb{R}$ respectively denote an observed binary treatment and an observed outcome variable for each dyad $(i, j)$. Throughout this paper, we focus on the case of undirected treatment such that $A_{ij} = A_{ji}$ for all $(i, j)$. As opposed to the treatment, the outcome variable can be specific to each unit in each dyad; that is, $Y_{ij} \neq Y_{ji}$ in general. Let $A = (A_{ij})$ and $Y = (Y_{ij})$ denote the $n \times n$ matrices of treatments and outcomes, respectively, where the diagonals are normalized to zero: $A_{ii} = 0$ and $Y_{ii} = 0$ for all $i \in [n]$.
The above setup encompasses many empirical situations of interest. For instance, in our empirical analysis, $A_{ij}$ represents the indicator for the presence of FTA between countries $i$ and $j$, and $Y_{ij}$ corresponds to the amount of exports from $i$ to $j$ or imports from $j$ to $i$. In other examples, $A_{ij}$ may denote whether firms $i$ and $j$ participate in a researchers exchange program, and $Y_{ij}$ represents the R&D collaboration output; $A_{ij}$ may be the indicator of political alliance between countries $i$ and $j$ with $Y_{ij}$ representing population migration; and so forth.
The dyadic data $\{ (Y_{ij}, A_{ij}) \}_{i, j \in [n]}$ generally exhibit complex dependence structure due to interactions within and across the dyads, which poses highly non-trivial statistical challenges, if no restrictions are imposed. Following the literature on graphon estimation, we address this problem by imposing the exchangeability on the probability distribution of the treatment matrix $A$. A graphon is a symmetric measurable function $f: [0,1]^2 \to [0,1]$, where the value of graphon $f(U_i, U_j)$ is considered to be the edge probability between $i$ and $j$. We say that $A$ is exchangeable if $(A_{ij}) \overset{d}{=} (A_{g(i)g(j)})$ for any permutation $g$ on $[n]$. The exchangeability plays an essential role in our analysis in combination with the celebrated Aldous--Hoover representation theorem (aldous1981representations; hoover1979relations): a symmetric random matrix $A = (A_{ij})$ is exchangeable if and only if there exists a symmetric function $f:[0, 1]^2 \to [0, 1]$ such that
where $\{U_i\}$ and $\{U_{ij}\}$ are mutually independent such that $U_i, U_{ij} \overset{\mathrm{IID}}{\sim} \mathrm{Uniform}[0, 1]$ with $U_{ij} = U_{ji}$. Obviously, the above claim is equivalent to $A_{ij} \mid U_i, U_j \sim \mathrm{Bernoulli}(f(U_i, U_j))$ with a symmetric function $f$. That is, the Aldous--Hoover theorem says that if the treatment matrix $A$ is exchangeable, then its associated graphon takes this particular form (and the converse is also true). We can interpret $U_i$ and $U_{ij}$ as the unit- and dyad-level latent composite attributes, respectively.
Formally, we assume the following:
As just mentioned, Assumption (ref) is equivalent to the exchangeability of $A$. Note that an obvious sufficient condition for exchangeability is the IID-ness of $\{A_{ij}\}$, which is commonly assumed in the network formation literature. Below, we provide three examples of treatment models that can be analyzed in our framework.
One notable limitation of Assumption (ref) is that the resulting network will automatically be a dense network. It is easy to see that the expected degree for each unit is a constant times $n-1$, which diverges as $n$ increases. In order to accommodate the sparsity, klopp2017oracle proposed the so-called sparse graphon estimation method. Extending our method in this direction would be important, but is left for future study.
Now, let $P_{ij} \coloneqq \mathbb{E}[ A_{ij} \mid U_i, U_j ] = f(U_i, U_j)$ denote the conditional probability of the treatment take-up given the pair of unit-level latent attributes, which we call the propensity score. Note that $P_{ij}$ is a function of unobserved $(U_i, U_j)$, differently from the standard propensity score defined with respect to observed covariates. Also notice, for all $(i, j)$, that $P_{ij} = P_{ji}$ by Assumption (ref) and that $P_{ij} = \mathbb{E}[ A_{ij} \mid P_{ij} ]$ by the law of iterated expectations. Since $A_{ii} = 0$, we normalize $P_{ii} = 0$ for all $i \in [n]$.
Next, we let $Y_{ij}(a)$ denote the potential outcome when $A_{ij} = a$. By construction, $Y_{ij} = A_{ij} Y_{ij}(1) + (1 - A_{ij}) Y_{ij}(0)$. The treatment effect of $A_{ij}$ at dyad $(i, j)$ is defined as $Y_{ij}(1) - Y_{ij}(0)$, which can be heterogeneous across dyads. This potential outcomes framework implicitly subsumes the stable unit treatment value assumption (rubin1974estimating) that rules out treatment interference or spillovers across dyads. Similarly as above, we set $Y_{ii}(0) = Y_{ii}(1) = 0$ for all $i \in [n]$.
As the main causal parameter of interest, we define the dyadic average treatment effect (DATE):
DATE helps us understand the nature of heterogeneity in treatment effects across dyads. For example, if $\tau_{ij}$ tends to take a positively large value for dyads with a higher $P_{ij}$, this suggests that the treatment would be more effective for dyads that are more likely to take the treatment (i.e., selection-on-returns). Other causal parameters of interest would include the individual average treatment effect (IATE) for unit $i$ and the global average treatment effect (GATE), which are respectively defined as
Since both IATE and GATE can be estimated once the DATEs are estimated, in the following, we mainly focus on DATE.
In this section, we first discuss the estimation of the DATE parameter. Then, the rate of convergence for the proposed estimator with respect to a certain norm is derived. After that, we introduce our permutation test for the sharp null hypothesis of no treatment effects and prove its size control property.
We first introduce the following two assumptions.
Assumption (ref) requires that only the unit-level latent factors $(U_i, U_j)$ can be the source of endogeneity. That is, the potential outcomes and treatment can be correlated only through $(U_i, U_j)$. To be more specific, we consider the following potential outcome equation:
where $y_a$ is an unknown function, and $\xi_{ij}(a)$ is a dyad-level composite attribute affecting the potential outcome independently of $(U_i, U_j)$. Then, Assumption (ref) is satisfied if the dyad-level composite factor $U_{ij}$ in the treatment choice equation is independent of those $(\xi_{ij}(0), \xi_{ij}(1))$ in the potential outcomes.
Assumption (ref) is violated when there is a dyad-level factor that affects both the potential outcomes and treatment. For instance, in the context of international trade, the connectedness of two countries would be such a factor, as it may be inherently difficult to explain the locations of two countries solely from a combination of independent unit-level attributes. However, since countries' locations are readily available data, this particular type of endogeneity can be addressed by restricting the analysis to pairs with similar locational patterns. In our empirical study, we mitigate the impacts of regional confounders by focusing on countries in the Asia-Pacific region. If there are unobservable dyad-level confounders, we might need additional data such as IV to address the endogeneity.
Assumption (ref) requires that each dyad $(i, j)$ has a propensity score strictly bounded away from zero and one. Even with this assumption, it is possible in practice for some of the estimated propensity scores to be very close to zero or one. In our numerical analysis, we simply censor such propensity scores at some pre-specified minimum or maximum value (e.g., 0.01 and 0.99).
Under these assumptions, we obtain the following simple result, which serves as the basis of our estimator.
An implication of Proposition (ref) is that we can estimate the DATE parameter by conducting a nonparametric regression of the outcome on the propensity score separately for the treatment and control groups. However, since $\{P_{ij}\}$ are unknown, this is not feasible. Alternatively, if we could identify a sub-sample satisfying $\{i' \in [n]: A_{i'j} = a, P_{i'j} \approx P_{ij}\}$, then computing the average of $Y_{i'j}$'s over this subset would give us a valid estimate of $m_{ij}(a)$. We pursue this task by adopting the approach proposed by zhang2017estimating. Let $d(i, i')$ denote the squared $\ell_2$ distance between graphon slices:
Intuitively, if $d(i,i')$ is close to zero, we would have $P_{ij} \approx P_{i'j}$ for all $j \neq i,i'$. Thus, if $d(i,i')$ were computable, we can estimate $m_{ij}(a)$ by computing the average outcome over the subsample $\{i' \in [n]: A_{i'j} = a, d(i,i') \approx 0\}$. Of course, directly estimating the value of $d(i,i')$ is not feasible. However, for the purpose of selecting a neighborhood of $i$, it is not necessary to know the precise value of $d(i,i')$; rather, it suffices to consider a tractable upper bound. zhang2017estimating have demonstrated that the following pseudometric $\widetilde d(i,i')$ serves as a viable approximate upper bound, in the sense that $\widetilde d(i,i') \approx 0$ implies $d(i,i') \approx 0$ with high probability:
Then, using this pseudometric in place of $d(i,i')$, nonparametric regression estimators for $m_{ij}(0)$ and $m_{ij}(1)$ are given by
respectively, where
and $b > 0$ is a bandwidth parameter such that $b \to 0$ as $n \to \infty$. In particular, $\widecheck P_{ij}$ is equivalent to the neighborhood smoothing estimator studied in zhang2017estimating.
One can see that the estimator of $m_{ij}(a)$ presented above corresponds to a uniform kernel regression estimator such that the neighborhood observations are weighted equally. As a slight generalization, we can consider the following general kernel estimation procedure, which we call the neighborhood kernel smoothing (NKS) estimator.
For a kernel weighting function $\mathcal{K}: \mathbb{R}_+ \to \mathbb{R}_+$, let
Notice that this estimator is not symmetric; that is, $\widetilde P_{ij}$ is generally different from $\widetilde P_{ji}$. Since we have assumed $P_{ij} = P_{ji}$ in Assumption (ref), we take the average of the two estimates for $(i, j)$ and $(j, i)$:
Then, with this propensity score estimator, the NKS estimators for $m_{ij}(0)$ and $m_{ij}(1)$ are given by
respectively. Finally, our proposed DATE estimator is defined as
Moreover, the estimators for IATE and GATE are obtained by
respectively.
In this subsection, we derive the convergence rate of the NKS estimator. We impose the following assumptions on the kernel weighting function $\mathcal{K}$ and the bandwidth $b$.
Assumption (ref) requires the kernel function $\mathcal{K}$ to be supported on the unit interval and bounded away from zero and above by a constant on the support. The uniform kernel function is a typical example that satisfies this assumption. Furthermore, commonly used compact support kernel functions can satisfy Assumption (ref) if we make slight modifications on their functional forms. For instance, in line with Assumption (ref), the triangular and Epanechnikov kernel functions can be respectively modified as
with some constant $\underline{C}_{\mathcal{K}} > 0$. The motivation behind this modification stems from the discrete nature of $\widetilde d(i, i')$. That is, due to the limited variation in the value of $\widetilde d(i,i')$, particularly for small $n$, it is often the case that a certain number of observations who are exactly $b$ distant away from $i$ exist. Consequently, when using standard kernel weight functions whose value degenerates to zero at the boundary, the resulting number of observations involved in the estimation can unintentionally be small, leading to a large variance. We numerically check this issue in Appendix (ref) -- see Table (ref) for more detail.
Assumption (ref) is a technical condition that is needed to derive the rate of convergence of our NKS estimator. The same assumption is introduced in zhang2017estimating. In our numerical simulations, we confirm that the estimator with $C_0 = 1$ performs reasonably well.
To derive the convergence rate of the NKS estimator, we assume that graphon $f$ is a piecewise-Lipschitz function. The following definition is due to Definition 2 of zhang2017estimating.
Lastly, we impose the following set of assumptions on the potential outcomes.
In Assumption (ref)(ii), we assume that the potential outcomes are bounded, which is in fact stronger than necessary, but significantly facilitates theoretical investigations. Assumption (ref)(iii) is a high-level condition. This can be satisfied, for instance when $\{Y_{ij}(a)\}$ are identically distributed such that $\mathbb{E}[Y_{ij} \mid P_{ij}] = g(P_{ij})$ and $g$ is $L_Y$-Lipschitz on $[0,1]$.
Now, we are ready to present our main theorem, which gives the convergence rate of the KNS estimator for the DATE parameter with respect to the following matrix norm: for $n \times n$ matrices $A$ and $B$ with zero diagonals, we define the normalized $(2, \infty)$ matrix norm as
where $\| \cdot \|_2$ denotes the Euclidean norm.
From Theorem (ref) and Jensen's inequality, we can obtain a bound for the convergence rate of the IATE estimator as follows:
Similarly, we can easily observe that $|\widehat \tau^\text{GATE} - \tau^\text{GATE}|^2 \lesssim_{\mathbb{P}} \sqrt{ (\log n) / n }$ holds.
Note that these bounds including the one given in Theorem (ref) are potentially very crude. It is straightforward to see that the same bound as in Theorem (ref) applies to $[n(n-1)]^{-1}\left\| \widehat \tau - \tau \right\|_F^2$, where $\left\|\cdot\right\|_F$ denotes the Frobenius norm. In the context of graphon estimation for a certain smooth graphon class, it is known that the minimax rate with respect to the norm $[n(n-1)]^{-1}\left\| \: \cdot \: \right\|_F^2$ is $(\log n)/n$ (gao2015rate,zhang2017estimating,gao2021minimax), which is square times faster than the derived bound. One example of a graphon estimator that achieves the optimal rate is the combinatorial least squares method proposed by gao2015rate. However, this method is often impractical due to its extreme computational complexity. It remains an open question whether a more computationally efficient algorithm exists that achieves the optimal rate.\footnote{ It is worth noting that the sorting algorithm proposed by chan2014consistent does not involve combinatorial optimization and can achieve the minimax rate, but assuming additionally that the graphon is monotonic. It is certainly possible to use these alternative graphon estimation methods to estimate the propensity scores and incorporate them into the DATE estimation, which should be an important topic of future research. In the current manuscript, considering its algorithmic simplicity, we advocate using the neighborhood smoothing method. }
As shown in the previous subsection, the rate of convergence of the NKS estimator is non-standard. Also, recalling that the construction of the pseudometric $\widetilde d(i,i')$ involves all information of $A$, the neighborhood $\mathcal{N}_i$ of each $i$ is endogenously determined. For these reasons, it is very challenging to develop statistical inference methods by deriving a tractable limiting distribution. Alternatively, in this subsection, we propose an easy-to-implement permutation inference method making full use of the exchangeability assumption.
We consider testing the following sharp null hypothesis:
Under $\mathbb{H}_0$, the observed outcome $Y_{ij}$ satisfies $Y_{ij} = Y_{ij}(0) = Y_{ij}(1)$, and hence the potential outcomes schedule $W \coloneqq \{ Y(0), Y(1) \}$ is imputable from the observed outcomes $Y$, where $Y(0) = (Y_{ij}(0))$, $Y(1) = (Y_{ij}(1))$, and $Y = (Y_{ij})$. Consequently, $\mathbb{H}_0$ implies that $\tau_{ij} = 0$ for all $(i, j)$, so that the rejection of $\mathbb{H}_0$ indicates the significant presence of $\tau_{ij}$ for some $(i, j)$.
For testing $\mathbb{H}_0$, let $T(A, W)$ denote any real-valued test statistic, which can be computed under $\mathbb{H}_0$. For example, we can consider the following test statistic:
Denote the set of all permutations of $[n]$ as $\bm{G}$, such that $|\bm{G}| = n!$. For a given permutation $g \in \bm{G}$, let $gA = (A_{g(i)g(j)})$ be the permuted treatment matrix. Then, the $p$-value for testing $\mathbb{H}_0$ is defined by
With a given nominal level $\alpha \in [0, 1]$, we reject $\mathbb{H}_0$ at the $100 \cdot \alpha\%$ significance level if $P(A, W) \le \alpha$.
The next theorem demonstrates the size control property of our test under the additional assumption that the pre-treatment outcomes $Y(0)$ and the treatments $A$ are independent.
Note that the additional independence condition in Theorem (ref) does not preclude the endogeneity arising from the dependence between $A$ and the treatment effect $Y_{ij}(1) - Y_{ij}(0)$. For example, the assumption allows selection-on-returns, where units with higher $Y_{ij}(1) - Y_{ij}(0)$ are more likely to take the treatment. In the presence of a more general form of endogeneity, there is no guarantee that (ref) holds, and our permutation test may not be valid. In such a case, an alternative testing procedure should be developed, but this task is left for future research.
In this section, we conduct a set of numerical simulations to evaluate the performance of our method. To save space, only the main results are reported here, and supplementary simulation results are provided in Appendix (ref).
\paragraph{Setup.}
We consider two data generating processes (DGPs) for propensity scores. The first is a stochastic block model with a covariate. Letting $X_{1i} \overset{\mathrm{IID}}{\sim} \mathrm{Uniform}[0, 1]$ be a unit-level covariate, the units are classified into three groups according to the following group assignment probability:
where $Z_i \in \{ 1, 2, 3 \}$ indicates the group to which agent $i$ is assigned. We then generate the propensity score at each dyad as follows:
For the second DGP to generate the propensity score, we borrow the same network formation model as Design A.1 in graham2017econometric:
where $F$ denotes the standard logistic cumulative distribution function, $X_{2i} \in \{ -1, 1 \}$ is a unit-level covariate such that $\mathbb{P}(X_{2i} = -1) = \mathbb{P}(X_{2i} = 1) = 0.5$, and $\beta_i$ represents an individual-level degree heterogeneity defined as $\beta_i = -0.5 + V_i$ with a centered Beta random variable $V_i \overset{\mathrm{IID}}{\sim} ( \mathrm{Beta}(1, 1) - 0.5 )$.
For both DGPs, the potential outcomes are generated by
where $\xi_{ij}, \zeta_{ij} \overset{\mathrm{IID}}{\sim} \mathrm{Normal}(0, 1)$ independent of $\{(X_{1i}, X_{2i}, V_i)\}$. We consider the following three patterns for the values of $\gamma$'s: (i) $\gamma_0 = \gamma_1 = \gamma_2 = 0$, (ii) $\gamma_0 = \gamma_1 = 1$, $\gamma_2 = 0$, and (iii) $\gamma_0 = \gamma_1 = \gamma_2 = 1$. Under these three setups, the DATE parameter becomes (i) $\tau_{ij} = 0$, (ii) $\tau_{ij} = P_{ij}$, and (iii) $\tau_{ij} = P_{ij} + P_{ij}^2$, respectively. Also notice that $Y_{ij}(0) = Y_{ij}(1)$ in case (i). Hence, we refer to these as the null setup and the linear and quadratic DATE setups, respectively.
For each setup, we consider two sample sizes: $n \in \{ 40, 80 \}$. For the choice of kernel weighting function, in view of Assumption (ref), we consider three kernel functions with $\underline{C}_{\mathcal{K}} = 1$: (i) uniform kernel, (ii) triangular kernel, and (iii) Epanechnikov kernel. In line with Assumption (ref), bandwidth $b$ is set to the $h = \sqrt{ (\log n) / n }$-th quantile of $\{ \widetilde d(i, j) \}_{j: j \neq i}$.
For each simulation scenario, we perform 1,000 Monte Carlo repetitions to assess the performance of the NKS estimators for estimating DATE, IATE, GATE, and propensity score (i.e., (ref), (ref), and (ref)) in terms of the average bias (AB) and average mean squared error (AMSE):
where the superscript $(r)$ indicates that it is obtained from the $r$-th replicated dataset. The AB and AMSE for the IATE, GATE, and propensity score estimates are defined analogously. For evaluating the permutation test, we report the rejection frequency (RF) among the 1,000 repetitions with nominal levels of $\alpha \in \{ 0.1, 0.05, 0.01 \}$. Considering the computational difficulty of calculating the $p$-value in (ref) due to the too large size of $|\bm{G}|$, we approximate it by randomly drawing 1,000 permutations from all $n!$ possible permutations. The test statistic used is as given in (ref).
\paragraph{Results.}
Tables (ref) and (ref) summarize the simulation results. For all simulation setups, AB and AMSE for the DATE and propensity score estimates are satisfactorily small, especially when $n$ is relatively large. As expected from the fact that IATE and GATE are obtained by aggregating DATEs, AMSEs for the IATE and GATE estimation are significantly smaller than those for the DATE estimation.
When comparing the kernel functions, the three kernels demonstrate quite similar performances. Particularly, the choice of kernel has almost no recognizable impact on the estimation of propensity scores. However, for the estimation of DATEs, overall, the triangular and Epanechnikov kernels tend to produce slightly smaller bias than the uniform kernel. In contrast, in terms of AMSE, the estimator with the uniform kernel outperforms the other two.
For the performance of our permutation test, the RFs in the null setup are sufficiently close to their nominal levels, which corroborates the result in Theorem (ref). In contrast, in the linear and quadratic DATE setups, the RFs are sufficiently high and tend to 100% as $n$ increases, suggesting the consistency of our permutation test.
As an empirical illustration of our method, we investigate the effects of free trade agreements (FTAs) on bilateral trade flows. This has been a classical and one of the most important empirical questions in international economics (e.g., tinbergen1962shaping). In the literature, there seems to be an agreement that FTAs indeed help increase exports and imports between member countries; however, it is still somewhat unclear to what extent they are actually beneficial; see, for example, the discussion in nagengast2023staggered. In addition, as mentioned above, since most papers on this topic use a regression-based gravity model approach, which is not necessarily based on a causal inference framework, our study would contribute to the literature in this respect.
In this empirical analysis, we focus on the trade flows across 37 selected countries and regions mainly in Asia, Oceania, and North America.\footnote{ The list of countries included in our dataset: Indonesia, Cambodia, Singapore, Thailand, Philippines, Brunei, Vietnam, Malaysia, Myanmar, Laos, Mongolia, India, Pakistan, Sri Lanka, Korea, China, Japan, Hong Kong, Australia, New Zealand, United Arab Emirates, Bahrain, Oman, Qatar, Kuwait, Saudi Arabia, Jordan, Lebanon, Israel, Canada, United States, Russia, Turkey, Kazakhstan, Tajikistan, Kyrgyz, and Uzbekistan. } The outcomes of interest are bilateral export and import amounts. Considering huge variation and unbalance of these values across countries, we transform them into percents of total exports and imports with all countries. The data of exports and imports for the year 2021 are obtained from the World Integrated Trade Solution (WITS), World Bank (\url{https://wits.worldbank.org}).
For the treatment variable, we consider all FTAs related to our sample countries that were enacted as of 2021. The data for the FTAs are obtained from the website of Japan External Trade Organization (JETRO) (\url{https://www.jetro.go.jp/theme/wto-fta/ftalist.html}). We do not differentiate between the contents of FTAs. If a pair of countries jointly participates in at least one FTA, we consider that dyad as treated.
Figure (ref) shows the FTA network of our data. Among the total of 666 dyads, the number of treated dyads is 197. In our dataset, Singapore has the largest number of FTA partners (degree 27), and the country with the smallest is Mongolia (degree 1).
Based on this dataset, we estimate the treatment effects of FTAs using the same estimation procedure as in the simulation experiment in Section (ref). For the choice of kernel function, we use the uniform kernel. The estimated IATEs on exports and imports for all countries are depicted in Figure (ref). We can observe that for both exports and imports, the impacts of FTAs are positive for almost all countries. In line with this, the permutation test of no treatment effects yields $p$-values of 0.052 and 0.018 for exports and imports, respectively, which would indicate the significant impacts of FTAs.
It appears that there is a positive correlation between the impacts on exports and those on imports, suggesting that FTAs are mutually beneficial for both “origin” and “destination” countries. For both exports and imports, the presence of an FTA increases the trade volume between countries by up to approximately 4 percent points. Averaged over the countries, the estimated GATE on exports is 1.938 and that on imports is 2.161, suggesting that the impact of establishing an FTA, on average, increases trade inflow and outflow by about 2 percent of the total. Meanwhile, we could not observe the impact of FTA only in Singapore. As mentioned above, Singapore has FTAs with almost three-quarters of the countries in the dataset, which might have made it more difficult to identify the impact of FTAs compared with other countries.
To investigate the nature of treatment heterogeneity, we draw scatter plots of IATE against individual average propensity score (IAPS), defined in a similar manner to IATE, in Figure (ref). Interestingly, the scatter plots reveal a certain negative correlation between the treatment effect and propensity score. This observation might indicate the possibility of decreasing marginal returns to the number of FTAs.
In this paper, we developed a statistical inference method for assessing the effects of a dyadic treatment on a dyadic outcome. Our main focus was on estimating the dyadic average treatment effect (DATE), which represents the difference in mean potential outcomes conditional on the propensity score with respect to unit-level latent attributes. By assuming the exchangeability for the distribution of dyadic treatments, we proposed the neighborhood kernel smoothing (NKS) estimator for estimating the DATE parameter and derived its convergence rate. For testing the null hypothesis of no treatment effects, we introduced a permutation inference method based on the exchangeability of treatments and the exogeneity of pre-treatment outcomes. To demonstrate the empirical usefulness of our method, we conducted an analysis on the impacts of FTAs on bilateral trade flows. Our findings suggested that enacting an FTA increases both exports and imports between two countries by approximately two percent points, on average.
The validity of our proposed approach relies heavily on the exchangeability assumption, which introduces important restrictions on the data. For instance, it implies a dense network structure for the treatment network and imposes limitations on the dependence structure within and across dyads. To address these issues, we can consider extending our method by incorporating sparse graphon estimators, such as the one in klopp2017oracle. Additionally, we could explore incorporating observable dyad-level covariates into our framework to accommodate some dependencies more effectively. These would be promising directions for future research.
This work was supported by JSPS KAKENHI Grant Numbers 20K01597 and 24K04817.