\documentclass[12pt]{article} \usepackage[hang,flushmargin]{footmisc} \setlength{\pdfpagewidth}{8.5in} \setlength{\pdfpageheight}{11in} \usepackage[utf8]{inputenc} \usepackage{csquotes} \usepackage{amsmath} \usepackage{setspace} \usepackage{geometry} \usepackage{booktabs} \usepackage{amssymb} \usepackage{threeparttable} \usepackage{graphicx} \usepackage{caption} \usepackage{subcaption} \usepackage{bbm} \geometry{margin=1in} \usepackage[backend=biber,style=apa,sorting=nyt, doi=false, url=false]{biblatex} \usepackage[labelfont=bf]{caption} \usepackage{titling} \usepackage[flushmargin]{footmisc} \addbibresource{refs.bib} %\AtBeginRefsection{\GenRefcontextData{sorting=ynt}} %\AtEveryCite{\localrefcontext[sorting=ynt]} \usepackage[hidelinks]{hyperref} \usepackage{tikz} \usepackage{pgfplots} \makeatletter \renewcommand{\@makefntext}[1]{% \setlength{\parindent}{0pt}% no paragraph indent in footnotes \noindent\@makefnmark\ #1% } \makeatother \begin{document} \thispagestyle{empty} { % set font to helvetica (arial) to make it 508-compliant \fontfamily{phv}\selectfont \begin{center} {\LARGE \textbf{Structural gravity and nonparametric trade costs}} \\ \vspace{0.25in} %{\LARGE \textbf{Second Line if Needed}} \\ \vspace{1.75in} {\Large Kevin D. Duncan} \\\vspace{.25in} {\Large Ward M. Reesman} \vspace{1.5in} {\large ECONOMICS WORKING PAPER SERIES}\\ Working Paper 2026--05--A \\ \vspace{0.5in} U.S. INTERNATIONAL TRADE COMMISSION \\ 500 E Street SW \\ Washington, DC 20436 \\ \vspace{0.5in} May 2026 \\ \end{center} \vfill %\noindent Optional: Acknowledgements, thanks, etc.\\\vspace{1em} \noindent The views expressed solely represent the opinions and professional research of the authors. The content of the working paper is not meant to represent the views of the U.S. International Trade Commission, any of its individual Commissioners, or the United States government. \newpage \thispagestyle{empty} % remove headers, footers, and page numbers from cover page \begin{flushleft} Structural gravity and nonparametric trade costs \\ Kevin D. Duncan and Ward M. Reesman \\ Economics Working Paper 2026--05--A\\ May 2026 \\~\\ \end{flushleft} \vfill \begin{abstract} \noindent This paper introduces a flexible approach to estimating nonlinear trade cost effects within the structural gravity framework using generalized additive models (GAMs). Unlike standard log-linear gravity models with constant elasticities, our method allows continuous policy and geographic trade cost variables to enter the trade cost vector nonparametrically while preserving theoretical consistency. Using bilateral trade data, we document significant nonlinearities: trade is relatively inelastic at low and high levels of trade costs and more elastic in central ranges of trade cost distributions. We show how these patterns are masked under conventional specifications, and discuss how estimates of state-dependent elasticities can inform general equilibrium counterfactuals for richer policy analysis. \end{abstract} \vfill \begin{flushleft} Kevin D. Duncan \\ kevindduncan@gmail.com \\ \vspace{2em} Ward M. Reesman \\ Research Division, Office of Economics\\ U.S. International Trade Commission\\ ward.reesman@usitc.gov \\ \vspace{3em} Suggested Citation: \\ \vspace{2em} Duncan, K. and Reesman, W. (2026). "Structural gravity and nonparametric trade costs," U.S. International Trade Commission. Economics Working Paper 2026-05-A. \\ \end{flushleft} } % end of helvetica (arial) font \clearpage \newpage \doublespacing \section{Introduction} Understanding how bilateral trade costs shape international trade flows is central to quantitative study of trade policy. The structural gravity framework delivers a theoretically grounded mapping from trade costs to bilateral trade flows while consistently accounting for general equilibrium forces through multilateral resistance terms \parencite{anderson_gravity_2003,head_chapter_2014}. Empirically, most applications rely on parametric and log-linear functional forms for bilateral trade costs, typically modeled as iceberg and with constant elasticity of substitution (CES) preferences \parencite{eaton_technology_2002, anderson_trade_2004,head_chapter_2014}. This approach is tractable, widely used, and remarkably successful at explaining variation in trade flows, especially when paired with pseudo-Poisson maximum likelihood (PPML) estimation and high-dimensional fixed effects \parencite{silva_log_2006,fally_structural_2015,yotov_advanced_2016}. However, these functional form assumptions rule out nonlinearities in the interaction between parameterized trade costs and trade flows, and they imply constant elasticities across the support of the distribution of observed trade costs. There are several reasons why trade elasticities may not be constant. On the demand side, CES preferences impose constant elasticity by assumption, but more flexible demand systems such as translog \parencite{feenstra2003homothetic, novy_international_2013} or \textcite{melitz2008market} formulations imply substitution patterns that vary with relative prices or market size. On the supply side, the marginal response of trade flows to trade cost shocks may differ across the level of trade costs due to extensive margin adjustment, heterogeneous firm-level productivity, or the presence of resource or capacity constraints or fixed entry costs. These mechanisms suggest that the relationship between changes in trade costs and bilateral trade may naturally vary across the support of observed frictions. This paper proposes a practical method to flexibly capture trade cost nonlinearities while maintaining the core structure and interpretability of gravity. We formulate the estimation problem as a generalized additive model (GAM), allowing some bilateral trade costs to enter nonparametrically. The GAM specification nests the standard parametric model and can be estimated with a Poisson distribution with typical high-dimensional exporter-time and importer-time fixed effects to capture multilateral resistance. By modeling trade cost covariates with smooth functions, we recover local (gradient) elasticities of trade that can vary across the support of the variable. This yields state-dependent partial trade effects and trade elasticities, providing a natural lens for testing and quantifying nonlinearities in how trade policy and geographic frictions shape trade flows. Methodologically, we leverage established tools from the GAM literature for estimation, penalized smoothing, and inference \parencite{hastie1986generalized, hastie1990generalized, wood2011reml, wood2015generalized}, adapted to the setting of structural gravity. We demonstrate the approach with an empirical application that combines bilateral trade flows with measures of tariffs and geographic trade costs. The empirical exercise is intentionally simple: we estimate gravity under (i) a conventional log-linear PPML specification and (ii) a GAM-based specification with nonparametric functions over bilateral distance and tariffs. We then compare implied elasticities and partial effects across specifications. The results document economically and statistically meaningful nonlinearities at empirically relevant ranges of tariffs and distance -- trade is relatively inelastic at both low and high levels of bilateral tariffs and distance, and relatively more elastic in the central region of the trade cost distributions. These patterns are masked or compressed by the log-linear specification but emerge clearly under the GAM estimator. Importantly, the GAM design leaves the fixed effects structure and identification intact, isolating partial effects while properly controlling for multilateral resistance. Our contribution is threefold. First, we provide a method to estimate gravity with nonparametric trade costs without sacrificing most core features (e.g. Poisson estimation, multilateral resistance, domestic trade, and interpretability of coefficients and partial effects).\footnote{Two core features we cannot include are zero trade flows and bilateral fixed effects.} Second, we show how to recover non-constant trade elasticity gradients with respect to individual continuous cost covariates. Estimated gradients provide a quantitative summary of nonlinearities that can be directly read from the smooth functions and their derivatives. Third, we discuss how the estimated gradient elasticities can be embedded in general equilibrium frameworks to analyze policy changes when the marginal effect of a cost shock depends on its initial level, where standard constant-elasticity models may misstate aggregate consequences \parencite{arkolakis2012new, egger_heterogeneous_2024}. Other approaches can generate nonlinear trade elasticities by imposing specific structural features (such as translog demand systems, firm heterogeneity with fixed costs, or supply constraints), and these models deliver valuable insights but often require stronger functional form assumptions or additional data \parencite{feenstra2003homothetic, novy_international_2013, melitz2008market}. Our contribution is complementary: we recover nonlinearities directly from observed variation in trade costs without the addition of structural restrictions beyond those already embedded in the gravity framework.\footnote{Additionally, our estimator is unique in its ability to identify heterogeneity over the support of the trade cost distribution. Approaches like translog gravity \parencite{novy_international_2013} only model heterogeneity bilaterally.} \section{Empirical method} For reference, we start with the canonical structural gravity model, which assumes a constant elasticity of substitution. It is given by % \begin{align}\label{eq:strucGrav} X_{ij} &= \left( \frac{\tau_{ij}}{P_{j} \Pi_{i}} \right)^{1-\sigma} Y_{i} E_{j} \\ \Pi_{i}^{1-\sigma} &= \sum_{j} \left( \frac{\tau_{ij}}{P_{j}} \right)^{1-\sigma} E_{j} \\ P_{j}^{1-\sigma} &= \sum_{i} \left( \frac{\tau_{ij}}{\Pi_{i}} \right)^{1-\sigma} Y_{i} \end{align} % where bilateral trade flows between exporter $j$ and importer $j$ ($X_{ij}$) are a function of exporter income $Y_{i}$, importer expenditure $E_{j}$, and bilateral trade frictions. In the model, trade frictions take three distinct forms. First, bilateral trade frictions ($\tau_{ij}$) reflect factors that alter the ability of any two countries to trade -- factors such as distance, contiguity, bilateral tariff and non-tariff barriers, and bilateral agreements. Second and third, indexes of aggregate trade costs for the exporter ($\Pi_{i}$) and the importer ($P_{j}$) capture outward and inward multilateral resistance, which have been shown to be important to theoretically consistent structural gravity models \parencite{anderson_gravity_2003}. All trade frictions are augmented by the elasticity of substitution ($\sigma$), which is typically assumed to be constant. It is important to note that while the presence of $\tau_{ij}$ as a generalized iceberg trade cost is guided by theory, the structure and content of it are not. In empirical applications, the researcher must choose a parameterization in order to incorporate a set of appropriate covariates to estimate trade elasticities, subject to the constraint $\tau_{ij} \geq 1$. The most commonly used structure for this is multiplicative \parencite{anderson_trade_2004}. Formally, for a $1 \times K_1$ vector of continuous or discrete variables $Z_{ij} = \{z_{ij}^{k}\}_{k=1}^{K_1}$ and a $1 \times K_2$ vector of binary variables $W_{ij} = \{w_{ij}^{k}\}_{k=1}^{K_2}$, this is % \begin{equation}\label{eq:mtc} \tau_{ij}^{1-\sigma} = \left( \; \prod_{k=1}^{K_1}\left(z_{ij}^{k} \right)^{\beta_{k}} \; \prod_{k=1}^{K_2}\exp(\beta_k w_{ij}^{k} ) \; \right). \end{equation} % This structure permits the following log transformation, % \begin{equation}\label{eq:log_mtc} (1-\sigma) \ln (\tau_{ij}) = \left(\sum_{k=1}^{K_1}\beta_{k} \ln z_{ij}^{k}+\sum_{k=1}^{K_2} \beta_k w_{ij}^{k} \right) \end{equation} % which lends itself well to substitution into the common empirical specification % \begin{equation*}\label{eq:empGrav} \ln X_{ij} = (1-\sigma) \ln (\tau_{ij}) + \ln (Y_{i}) + \ln (E_{j}) + (\sigma-1) \ln (\Pi_{i}) + (\sigma-1) \ln (P_{j}) + \varepsilon_{ij} \end{equation*} % where $(1-\sigma)\ln(\tau_{ij})$ is specified as outlined in (\ref{eq:log_mtc}), and we note income, expenditure, and multilateral resistance terms as defined in (\ref{eq:strucGrav}) can be addressed through the use of importer and exporter fixed effects:\footnote{Multilateral resistance terms can also be computed or approximated through a variety of methods \parencite{anderson_gravity_2003, baier_bonus_2009, fally_structural_2015, freeman2025unlocking}.} % \begin{equation}\label{eq:empGravFE} \ln X_{ij} = (1-\sigma) \ln (\tau_{ij}) + \phi_{i} + \phi_{j} + \varepsilon_{ij}. \end{equation} % This can then be estimated either via OLS, or further exponentiated and estimated via Poisson pseudo-maximum likelihood (PPML) as widely recommended \parencite{yotov_advanced_2016, silva_log_2006, fally_structural_2015}. In the latter case, the estimating equation is % \begin{equation*} X_{ij} = \exp \Big( (1-\sigma) \ln (\tau_{ij}) + \phi_{i} + \phi_{j} \Big) \xi_{ij} \end{equation*} % where we again note $(1-\sigma)\ln(\tau_{ij})$ is defined as in (\ref{eq:log_mtc}), which yields a tractable log-linear design for estimation. Estimates on trade cost covariates have a structural representation: $\hat{\beta}_k = (1-\hat{\sigma}) \rho_k$ where $\rho_k$ is the elasticity of \textit{trade costs} with respect to covariate $z^k$. We note for tariffs (i.e., $z_{ij} = (1+t_{ij})$ where $t_{ij}$ is the ad-valorem tariff rate), the estimated coefficient $\hat{\beta}$ exactly identifies the trade elasticity $\sigma$.\footnote{Provided the trade data used in estimation are expressed at CIF prices. This property emerges as tariffs are an exact price shifter. See \textcite{yotov_advanced_2016}, \textcite{HEID201670}.} It is clear this parameterization and estimation strategy imposes a direct log-linear link between trade policy, embodied in $\tau_{ij}$, and bilateral trade. However, there are reasons to suspect there may be nonlinear effects of trade policy or natural trade costs on bilateral trade.\footnote{\textcite{egger_heterogeneous_2024} find broad nonlinear effects of tariff and nontariff barriers on bilateral trade costs. Other examples include \textcite{herman2024too} and \textcite{larch_estimating_2024} on free trade agreement provisions, and \textcite{eaton_technology_2002, henderson_is_2008, hillberry_trade_2008} on natural (geographic) trade costs.} Such a hypothesis cannot be tested using the canonical design. Rather than assume a direct log-linear relationship between trade costs and trade flows as in (\ref{eq:log_mtc}), we consider specifying a more flexible form, % \begin{equation} \tau_{ij}^{1-\sigma} = \prod_{k=1}^{K_1} f_k(z_{ij}^{k}) \prod_{k=1}^{K_2} \exp( \beta_{k} w_{ij}^{k}) \end{equation} that generalizes the relationship between continuous trade cost covariates $z_{ij}^{k}$ and total trade costs $\tau_{ij}$ with unknown smooth functions $f_k(\cdot)$, indexed by $k$.\footnote{One could also consider a parameterization that includes $K_1$ continuous or discrete covariates modeled with smooth functions, and $K_2$ continuous or discrete covariates and $K_3$ binary or factor covariates modeled linearly as in (\ref{eq:log_mtc}). We also note binary or factor variables can be added as interactions, allowing a smooth function $f_k(\cdot)$ to vary by group, which may be insightful based on the research question.} Applying the log linearization as before, bilateral trade costs take the shape % \begin{equation} (1-\sigma) \ln (\tau_{ij}) = \sum_{k=1}^{K_1} \ln f_{k}(z_{ij}^{k}) + \sum_{k=1}^{K_2} \beta_k w_{ij}^{k}. \end{equation} % Then, we can estimate (\ref{eq:empGravFE}) as a generalized additive model (GAM), % \begin{equation}\label{eq:gam} X_{ij} = \exp \left( \sum_{k=1}^{K_1} g_{k}(z_{ij}^k) +\sum_{k=1}^{K_2} \beta_k w_{ij}^{k} + \phi_i + \phi_j \right)\xi_{ij}, \end{equation} % where $g_k(\cdot) \equiv \ln f_k(\cdot)$ generalizes the relationship between trade cost covariates $z_{ij}^{k}$ and bilateral trade flows.\footnote{It is assumed that the response variable $X_{ij}$ is distributed via some exponential family and is related to smooth functions of predictor variables through a link function. For a Poisson or Gamma model, this is a natural logarithm function. For a Gaussian model, this is the identity function.} Under this design, the marginal effect of a change in trade cost measure $z_{ij}^{k}$ on aggregate bilateral trade costs is characterized by $f_k'$, and the marginal effect of a change in trade cost measure $z_{ij}^{k}$ on bilateral trade flows is characterized by $g_k'$, which are structurally linked.\footnote{In our empirical application, we report estimates of $g_k(\cdot)$ and thus characterize the full trade flow effects.} Importantly, these marginal effects can take any shape and may not be constant for all values of $z_{ij}^{k}$. The smooth functions $g_k(\cdot)$ are unknown and must be estimated. In principle, one could estimate each $g_k$ using a full spline, fitting a nonparametric curve to the data. However, this is prone to overfitting and quickly becomes computationally intensive.\footnote{Particularly so for large samples, as the computational cost of a full spline approach is $O(n^{3})$.} Instead, we approximate each function $g_k$ using a basis expansion: a weighted sum of simpler known functions (splines or polynomials) defined over subsets of the covariate space. Formally, we approximate the functions as % \begin{equation} g_k(z) = \sum_{m=1}^{M_k} \theta_{km} B_{km}(z) \end{equation} % where the $B_{km}(\cdot)$ are known as basis functions, and $\theta_{km}$ are coefficients to be estimated in the model fitting process for each subset $m$ of the support of $z$. The basis dimensions $M_k$ are selected to be sufficiently large to avoid bias from model oversimplification, but sufficiently small to retain computational efficiency.\footnote{For $p = \sum M_{k}$, the computational cost of estimation is $O(np^{2})$. For $p < n$ this is significantly more efficient than a full spline approach.} To avoid overfitting, we choose to model the basis functions $B_{km}(z)$ as penalized regression splines which impose a penalty on smoothness of the function. We estimate the model via restricted maximum likelihood (REML), which maximizes a penalized log-likelihood function given by % \begin{equation*} l(\boldsymbol{\theta}) - \sum_k \lambda_k \int \left(g''_{k} (z)\right)^{2} dz \end{equation*} % where $l(\boldsymbol{\theta})$ is the log-likelihood function and the $\lambda_k$ are penalty terms \parencite{wood2011reml}. Practically, the objective function is % \begin{equation} l(\boldsymbol{\theta}) - \sum_{k=1}^{K} \lambda_k \boldsymbol{\theta}_{k}^{T} \mathbf{S}_k \boldsymbol{\theta}_k \end{equation} % where $\boldsymbol{\theta}_k$ is the vector of basis coefficients for the smooth function $g_k$, $\mathbf{S}_k$ is the penalty matrix derived from basis functions (i.e., approximating the integral of $(g''_{k} )^2$), and the penalty terms $\lambda_k$ are estimated jointly with all other parameters in the model.\footnote{REML estimation treats the smooth terms as Gaussian random effects, and the penalty therefore becomes a variance component in a mixed model, which is what permits this matrix representation. Other common GAM estimation approaches do not jointly estimate $\lambda_k$, rather selecting these through cross validation.} \section{Demonstration} In this section, we demonstrate the approach outlined in the previous section with a simple application featuring both policy-driven trade costs (tariffs) and geographic barriers (distance). We briefly describe the data used before outlining an illustrative specification and discussing some estimation issues. We then present results of the estimation, and conclude with a short discussion of estimation implications. \subsection{Data} To estimate the model, we leverage a panel of bilateral trade and trade cost data. Bilateral trade data is sourced from the International Trade and Production Database for Estimation (ITPD-E), which document trade flows across 170 industries for over 200 countries \parencite{borchert2021international}. We use an aggregated country-level sample derived from these data covering 100 importing countries and over 220 exporting countries, representing over 90\% of global trade flows. For estimation, we first consider a limited sample over the period 2010 to 2019 for computational reasons, but in our final results expand this sample to cover the period 2000 to 2019. For both panels, we use data for consecutive years on the recommendation of \textcite{egger2022gravity}. The longer panel is helpful in introducing a wider distribution of bilateral tariff rates to aid in identification, but this does comes at a computational cost. An important feature of the ITPD data is the inclusion of domestic trade data derived from administrative data on domestic production, which are necessary for proper identification in structural gravity models. However, there is a concern related to the presence of missing domestic trade observations in the ITPD-E dataset (due to data availability issues in underlying primary sources) leading to significantly understated domestic market sizes when industries are aggregated to the country level.\footnote{This is not an econometric issue when using ITPD-E data at the disaggregated industry level. We choose to aggregate for computational feasibility, given the nonparametric estimator, and as such we must deal with the missing data issue.} This in turn leads to biased estimates. To address this concern, we use domestic trade data from the International Trade and Production Database for Simulation (ITPD-S) database, which uses a variety of methods to fill in missing domestic trade values \parencite{borchert2024itpds}. While an imperfect solution to the missing data problem, the inclusion of proxied production data delivers country-level domestic trade flows that are closer to accurate.\footnote{Estimation was also carried out using aggregated ITPD-E data without adjusting for missing domestic trade values and thus overstating international trade relative to domestic trade. Results were not significantly different but estimates do suffer from the measurement error of understated aggregate domestic trade.} Natural trade cost variables, such as distance and contiguity, as well as additional policy variables such as joint membership in preferential trade agreements (PTAs) and the World Trade Organization (WTO), are sourced from version 2.1 of the Dynamic Gravity Database \parencite{gurevich2018dgd}. We also add an indicator variable to denote international versus domestic trade flows, as recommended in gravity applications with domestic trade data. Bilateral tariff rates are sourced from a beta version of Feodora Teti's Global Tariff Database (v\_beta1-2024-12) from \textcite{teti2024missing}. These data are based on tariff data collected and organized by Trade Analysis Information System (TRAINS) and distributed by World Integrated Trade Solutions (WITS), but correct for erroneous interpolation and selection bias issues.\footnote{For more information on these data, see \textcite{teti2024missing}.} Some summary statistics are presented in Table \ref{tab:sumstats} for our sample of aggregated bilateral trade flows, tariff rates, and standard gravity variables, covering the shorter sample period in panel (a) and the longer sample period in panel (b). \subsection{Methodology} Our illustrative specification is of the class outlined in (\ref{eq:gam}), with the addition of a panel dimension to introduce additional variation in bilateral tariff rates. We estimate \begin{equation} X_{ijt} = \exp \Big( g_1 (t_{ijt}) + g_2 (d_{ij}) + W_{ijt}\beta + \phi_{it} + \phi_{jt} \Big)\xi_{ijt} \end{equation} where $t$ denotes years and we consider $K=2$ continuous trade cost covariates to be modeled flexibly: bilateral tariffs $t_{ijt}$ (defined as 1 plus the ad-valorem tariff rate), and log bilateral distance $d_{ij}$.\footnote{We log transform distance $d_{ij}$ in this application to scale outliers, as the distribution of bilateral distance in the sample is very skewed. This is not structurally required. The results we present are very similar to those for unscaled distance in the ranges of the distribution that are well-identified.} We use penalized thin plate regression splines to estimate $g_k(\cdot)$, and set a maximum number of knots at six.\footnote{The number of knots determines the effective degrees of freedom of the smooth fits.} As suggested by theory, we include (linear) exporter-year and importer-year fixed effects to address multilateral resistance and control for exporter income and importer expenditure.\footnote{Pair fixed effects are omitted due to computational constraints and the tendency of (too many) high-dimensional fixed effects in GAMs to inflate variance in smooth terms, straining the penalized likelihood and rendering basis functions unidentifiable.} $W_{ijt}$ contains other trade cost variables that are modeled linearly; namely, these are binary indicators such as contiguity and common language, as well as bilateral membership in PTAs, the WTO, and the EU. We estimate the model assuming a quasi-Poisson distribution for $X_{ijt}$, which is consistent with the recommended use of PPML in standard structural gravity estimation due to favorable treatment of heteroskedasticity, overdispersion in the conditional mean-variance relationship, and zero trade flows. To address additional concerns related to heteroskedasticity and possible sensitivity of the estimated smooth terms in the tails of covariate distributions, our standard errors are bootstrapped. The bootstrap resamples at the country-pair level to allow for panel interdependence. While zero trade flows can be appropriately modeled with the use of the Poisson distribution, mass points in the data pose a challenge for estimation of the smooth terms $g_k(\cdot)$. Zero trade flows form a significant mass point is our bilateral trade data and are thus problematic for identification. Therefore, we do not include zero trade in the sample. This nullifies one of the advantages of the Poisson model, and is a drawback of the approach.\footnote{This is also a drawback of the nonparametric method of \textcite{egger_heterogeneous_2024}.} Methodological designs to capture flexible nonlinear trade cost effects in a theoretically consistent manner while still accounting for zero trade flows are worth further exploration in future research. \subsection{Results} % benchmark To provide a benchmark for the parametric coefficients of the GAM estimation, we first present standard PPML estimates in Table \ref{tab:ppml}. Columns (1) and (2) present estimates using the shorter sample over 2010-2019, with column (1) reporting results for the closest specification to our GAM design, and column (2) including bilateral fixed effects as another reference. Columns (3) and (4) report the same but for the fuller sample over 2000-2019. As discussed before, the point estimate on $\ln(1+t_{ijt})$ is an estimate of $\sigma$, which is largely within ranges identified in previous trade literature. % parametric coefs The parametric coefficients of the GAM estimation are presented in Table \ref{tab:gam}. In column (1), we report estimates using the 2010-2019 sample, and in column (2) we report estimates using the 2000-2019 sample. Point estimates are largely the expected sign and of similar magnitude to the benchmark PPML results. Note, however, the standard errors cannot be clustered or otherwise adjusted for heteroskedasticity, so reported standard errors are somewhat deflated and inference must be conducted with caution.\footnote{This is another drawback of the naive GAM approach, and provides additional motivation for the bootstrap.} For the smooth terms, estimates are not easily reported in table form. Rather, Figure \ref{fig:gam_raw} plots the raw estimated smooth functions over tariffs and distance, respectively, with standard errors for the 2010-2019 sample. Figure \ref{fig:gam_raw_full} presents the same for the 2000-2019 sample. While inference should be cautious in the absence of a bootstrap for the standard errors, the point estimates display significant nonlinearities. Panel (a) of both figures displays the estimated function over bilateral tariffs. Given the use of CIF import values, the gradient can be thought of as an estimate of the trade elasticity $\hat{\sigma}$ over the support of the distribution of tariff rates. The trade elasticity is relatively flat at low tariff rates but turns negative at around 7.5\% ad valorem tariff rates for both samples. Along the right tail of the distribution, the trade elasticity flattens. In the 2010-2019 sample, it eventually slopes upward but this is likely poorly identified due to sparsity in this region of the distribution. Panel (a) of Figure \ref{fig:gam_raw_full} shows with broader support in the tariff distribution over the 2000-2019 period, the trade elasticity indeed flattens around 15\% ad valorem tariff rates, before turning more negative in the far right tail where again, the distribution becomes much sparser.\footnote{These identification issues become clearer when we turn to the bootstrapped results below.} Panel (b) of both figures displays the estimated function over bilateral distance. In both cases, at low values of distance, the marginal effect on trade is flat, and even weakly positive. The elasticity becomes sharply negative around the median of the distribution, before flattening in the right tail, with the slope halving in the far right tail. Given the log transformation, the right tail has been compressed and the estimates over this region of the distribution are well-identified. Illustrative results using raw distance (in kilometers) replicate these diminishing marginal effects at high values, though with much poorer identification due to sparsity and skew in the distribution of unscaled bilateral distances. % bootstrapped results Figure \ref{fig:gam_boot} presents results of the estimation with the bootstrap.\footnote{Bootstrap samples are drawn ensuring that all importer-exporter combinations exist in each bootstrap sample. This corresponds to pair-wise block-bootstrapping.} We use the full sample period over 2000-2019 to saturate the distribution of bilateral tariff rates to aid in identification.\footnote{Due to the computational burden of estimating the model on the full sample, we limit the number of bootstrap samples to 100.} Point estimates are largely preserved, though at low tariff values the trade elasticity is much flatter, reflecting the importance of the bootstrap in addressing sensitivity in the point estimates of $g_k(\cdot)$. Additionally, steepening elasticities along the right tail are poorly identified due to the sparsity of high-tariff trading pairs in the sample, illustrated by ballooning confidence intervals. With respect to bilateral distance, statistically significant nonlinearities are still identified despite some identification issues in the left tail of the distribution. Partial trade effects are relatively flat until the median of the distribution, when it turns negative before flattening in the right tail. Taken together, our results suggest trade is relatively inelastic when trading partners are closer and tariffs are low, but becomes much more elastic as both natural and policy-driven trade costs increase, with some diminishing marginal effects -- in the case of distance, partial effects on trade flatten in the right tail of the distribution, and for tariffs, the trade elasticity flattens again around 15\% ad valorem bilateral rates. Tariff effects may grow larger as tariffs surpass 20\% ad valorem rates, but this is poorly identified given the lack of high tariff observations even in the broader sample period. % ad hoc squared terms for comparison To provide a point of comparison for these nonlinear effects, we estimate log-linearized structural gravity with the inclusion of ad hoc squared terms for tariffs and distance. Structurally, these do not make sense given the parameterization and log linearization of continuous covariates in bilateral iceberg trade costs.\footnote{That is, these are modeled as $(1-\sigma)\ln(\tau_{ij}) = \beta_1 \ln(z_{ij}) + \beta_2 \ln(z_{ij})^{2}$, with the square applied ad hoc after log linearization. This design cannot be traced back to a sensible model of multiplicative iceberg trade costs. If you start instead with $\tau_{ij}^{1-\sigma} = z_{ij}^{\beta_1} * (z_{ij}^{2})^{\beta_2}$, then log linearization yields $(1-\sigma)\ln(\tau_{ij}) = \beta_1 \ln(z_{ij}) + 2*\beta_2 \ln(z_{ij}) = (\beta_1 + 2*\beta_2) \ln(z_{ij})$ which does not capture a quadratic effect. Another option is $\tau_{ij}^{1-\sigma} = Z_{ij}^{\beta_1}$ where $Z_{ij} = z_{ij} + \beta_2 z_{ij}^{2}$, yielding $(1-\sigma)\ln(\tau_{ij}) = \beta_1 \ln(z_{ij} + \beta_2 z_{ij}^{2})$ where the squared term is not separately identifiable.} Because these quadratic terms are ad hoc and lack structural grounding in multiplicative trade cost settings, we discuss them here only as a rough comparison for our nonparametric results. Estimates for these specifications are reported in Table \ref{tab:ppml_quad}. For each sample period, we first square tariffs alone in columns (1) and (2), then distance alone in columns (3) and (4), and then jointly square tariffs and distance in columns (5) and (6). In each case, the first column reports results using the 2010-2019 sample, and the second column reports results using the 2000-2019 sample. The results suggest an ``inverse U shape" for both tariffs and distance, where as bilateral trade costs rise, the trade elasticities become larger. While these are consistent with some regions of the estimated smooth plots, the quadratic structure does not permit the identification of further nonlinearities, e.g. the diminishing marginal effects in the right tail of the distance distribution, and the flattening of the trade elasticity around 15\% ad valorem tariffs. The use of quadratic terms verifies some of the nonlinear effects we identify with the GAM, but masks other components -- particularly in the distributional tails, even when well-identified. \subsection{Discussion} \subsubsection*{Practical implications for policy modeling} Our findings reveal that trade elasticities vary substantially across the distribution of trade costs, challenging the common assumption of constant elasticities in structural gravity models. This has important implications for policy analysis. For example, tariff reductions at very low initial levels may yield smaller trade responses than predicted by log-linear models, while cuts from moderate levels can have much larger effects. Conversely, at high tariff levels, marginal effects diminish again, suggesting that deep liberalization may not always deliver proportional trade increases. These nonlinear patterns matter for welfare analysis and counterfactual simulations. Standard models compress marginal effects into a single elasticity, which can bias predictions. This is most influential in a setting with large shocks or heterogeneous bilateral relationships. Incorporating state-dependent elasticities into general equilibrium frameworks would allow policy simulations to reflect the true shape of trade responses, improving the accuracy of welfare estimates and trade forecasts. While we do not present a full simulation here, the estimated gradients provide a natural input for such exercises and highlight the potential for richer policy insights.\footnote{For an example of how nonparametrically estimated elasticities impact simulation outcomes, see \textcite{egger_heterogeneous_2024}.} % \subsubsection*{Methodological extensions} Several methodological issues merit further exploration. First, our approach excludes zero trade flows, which limits its applicability, especially if zero trade flows are economically meaningful. Developing semi-parametric or zero-inflated GAM estimators could address this limitation while preserving theoretical consistency. Second, while we adopt a Poisson specification for comparability with PPML, alternative distributional assumptions may offer computational advantages and different robustness properties. A systematic comparison of these alternatives would clarify trade-offs between flexibility and efficiency. Third, the nonlinearities uncovered here suggest that standard PPML estimates may suffer from functional form misspecification. Quantifying the bias introduced by constant-elasticity assumptions would help assess the practical importance of adopting nonparametric methods. Finally, while ad hoc quadratic terms are commonly used to approximate nonlinearities, as we offer as a point of rough comparison, they lack structural grounding and cannot capture complex patterns. Future work could explore spline-based parametric approximations or finite exponential curve-fitting as interpretable, computationally efficient alternatives to fully nonparametric designs. \section{Conclusion} In this paper, we propose a new method to flexibly estimate nonlinearities in trade costs within structural gravity using GAMs. The approach maintains the core structure and interpretability of gravity, and can be estimated with a Poisson distribution and with high-dimensional fixed effects to address multilateral resistance, while allowing for continuous bilateral trade costs to enter nonparametrically. We estimate nonparametric functions over both policy-driven and geographic trade cost variables, and show how the resulting estimates exhibit strong nonlinearities over the support of trade cost distributions. The estimated gradients can be easily extracted for use in quantitative trade models to study the effects of trade policy shocks in general equilibrium while accounting for a rich structure of heterogeneity. However, while this approach is practical and easy to implement, there remain estimation issues (e.g., concerning zero trade flows and bilateral fixed effects) and fruitful methodological extensions worth exploration that we leave to future research. \printbibliography \pagebreak \begin{table}[!htbp] \centering \renewcommand*{\arraystretch}{1.1}\caption{Summary Statistics}\resizebox{\textwidth}{!}{ \input{tables/sumstats} \label{tab:sumstats} }\end{table} \begin{table}[!htbp] \centering \caption{PPML estimation} \input{tables/ppml_baseline} \label{tab:ppml} \begin{tablenotes} \footnotesize \centering \item Clustered (pair) standard errors in parentheses. \item * $p<0.01$, ** $p<0.05$, * $p<0.1$ \end{tablenotes} \end{table} \begin{table}[!htbp] \centering \caption{Parametric GAM coefficients} \input{tables/gam_parametric_coef} \label{tab:gam} \begin{tablenotes} \footnotesize \centering \item Standard errors in parentheses. \item * $p<0.01$, ** $p<0.05$, * $p<0.1$ \end{tablenotes} \end{table} \begin{table}[!htbp] \centering \caption{PPML estimation with quadratic trade costs} \input{tables/ppml_quadratic} \label{tab:ppml_quad} \begin{tablenotes} \footnotesize \centering \item Clustered (pair) standard errors in parentheses. \item * $p<0.01$, ** $p<0.05$, * $p<0.1$ \end{tablenotes} \end{table} \pagebreak \begin{figure} \centering \caption{Poisson GAM, 2010-2019 sample} \begin{subfigure}[t]{0.8\linewidth} %\centering \caption{Estimate of $g(1+t_{ijt})$} \begin{tikzpicture} \begin{axis}[ width=\linewidth, axis lines=none, ticks=none, enlargelimits=false, clip=false ] % Insert external graphic \addplot graphics[xmin=0, xmax=1, ymin=0, ymax=1] {figures/smooth_tariff.png}; % Overlay axis labels \node at (axis description cs:0.5,-0) {$1+t_{ijt}$}; \node[rotate=90] at (axis description cs:-0.01,0.5) {$g(1+t_{ijt})$}; \end{axis} \end{tikzpicture} \label{fig:gam_raw_tariff} \end{subfigure} \begin{subfigure}[t]{0.8\linewidth} %\centering \caption{Estimate of $g(d_{ij})$} \begin{tikzpicture} \begin{axis}[ width=\linewidth, axis lines=none, ticks=none, enlargelimits=false, clip=false ] \addplot graphics[xmin=0, xmax=1, ymin=0, ymax=1] {figures/smooth_dist.png}; \node at (axis description cs:0.5,-0) {$d_{ij}$}; \node[rotate=90] at (axis description cs:-0.01,0.5) {$g(d_{ij})$}; \end{axis} \end{tikzpicture} \label{fig:gam_raw_dist} \end{subfigure} \label{fig:gam_raw} \end{figure} \begin{figure} \centering \caption{Poisson GAM, 2000-2019 sample} \begin{subfigure}[t]{0.8\linewidth} %\centering \caption{Estimate of $g(1+t_{ijt})$} \begin{tikzpicture} \begin{axis}[ width=\linewidth, axis lines=none, ticks=none, enlargelimits=false, clip=false ] % Insert external graphic \addplot graphics[xmin=0, xmax=1, ymin=0, ymax=1] {figures/smooth_tariff_full.png}; % Overlay axis labels \node at (axis description cs:0.5,-0) {$1+t_{ijt}$}; \node[rotate=90] at (axis description cs:-0.01,0.5) {$g(1+t_{ijt})$}; \end{axis} \end{tikzpicture} \end{subfigure} \begin{subfigure}[t]{0.8\linewidth} %\centering \caption{Estimate of $g(d_{ij})$} \begin{tikzpicture} \begin{axis}[ width=\linewidth, axis lines=none, ticks=none, enlargelimits=false, clip=false ] \addplot graphics[xmin=0, xmax=1, ymin=0, ymax=1] {figures/smooth_dist_full.png}; \node at (axis description cs:0.5,-0) {$d_{ij}$}; \node[rotate=90] at (axis description cs:-0.01,0.5) {$g(d_{ij})$}; \end{axis} \end{tikzpicture} \end{subfigure} \label{fig:gam_raw_full} \end{figure} \begin{figure} \centering \caption{Bootstrapped Poisson GAM, 2000-2019 sample} \begin{subfigure}[t]{0.8\linewidth} \caption{Bootstrapped estimate of $g(1+t_{ijt})$} \begin{tikzpicture} \begin{axis}[ width=\linewidth, axis lines=none, ticks=none, enlargelimits=false, clip=false ] % Insert external graphic \addplot graphics[xmin=0, xmax=1, ymin=0, ymax=1] {figures/smooth_tariff_boot.png}; % Overlay axis labels \node at (axis description cs:0.5,-0) {$1+t_{ijt}$}; \node[rotate=90] at (axis description cs:-0.01,0.5) {$g(1+t_{ijt})$}; \end{axis} \end{tikzpicture} \end{subfigure} \begin{subfigure}[t]{0.8\linewidth} \caption{Bootstrapped estimate of $g(d_{ij})$} \begin{tikzpicture} \begin{axis}[ width=\linewidth, axis lines=none, ticks=none, enlargelimits=false, clip=false ] \addplot graphics[xmin=0, xmax=1, ymin=0, ymax=1] {figures/smooth_dist_boot.png}; \node at (axis description cs:0.5,-0) {$d_{ij}$}; \node[rotate=90] at (axis description cs:-0.01,0.5) {$g(d_{ij})$}; \end{axis} \end{tikzpicture} \end{subfigure} \label{fig:gam_boot} \end{figure} \end{document}