Doubly constrained gravity models for interregional trade estimation

Correspondence Mattia Cai, European Commission – Joint Research Centre, Edificio Expo, Calle Inca Garcilaso, 3, 41092 Sevilla, Spain. Email: mattia.cai@ec.europa.eu Abstract This paper discusses a family of methods grounded in the doubly constrained gravity model (DCGM) for the estimation of product-specific origin–destination matrices of interregional trade. We argue that several estimation procedures documented in the literature, although outwardly unrelated, can be conceptualized as applications of the DCGM framework. We show that DCGM estimation requires less restrictive assumptions and fewer data than commonly thought. We demonstrate that, in ideal conditions, the unknown trade flows can be recovered exactly through a standard application of the RAS algorithm. Finally, we examine how a DCGM estimator might fare in applications using Monte Carlo techniques and give an example of its use with real-world data.

survey-based statistics on interregional trade. This paper is concerned with how to bridge that information gap in applied work. Specifically, we discuss a family of methodological approaches to the estimation of product-specific origin-destination (OD) matrices of interregional trade.
The question of how to estimate bilateral trade flows at the sub-territorial level has attracted relatively little attention from economists since Schwarm, Okuyama, and Jackson (2006, p. 84) complained about the "dearth of examples in the literature."Understandably, regional analyses that concentrate primarily on a single territorial unit generally show little interest in constructing a full set of trade flows explicitly linking all parts of the country. Instead, they focus on estimating only those economic aggregates that are strictly necessary for their study. Thus, for example, input-output studies often rely on non-survey regionalization of model coefficients estimated from country-level data rather than attempt estimating interregional trade (Hermannsson & McIntyre, 2014;Loizou, Chatzitheodoridis, Polymeros, Michailidis, & Mattas, 2014;Morrissey, 2016;PwC, 2014;Yu, Hubacek, Feng, & Guan, 2010). In recent years, non-survey methods have been at the centre of a large and rapidly growing body of research (Bonfiglio, 2009;Bonfiglio & Chelli, 2008;Flegg & Tohmo, 2013Flegg & Webber, 2000;Flegg, Webber, & Elliott, 1995;Kowalewksi, 2015;Kronenberg, 2009;Lamonica & Chelli, 2018;Többen & Kronenberg, 2015).
Although comparatively limited in size, the literature on interregional trade estimation has already experimented with a number of approaches. Frequently, the starting point is provided by data on freight transport. These are typically subjected to numerous and often complex adjustments in order to translate the observed flows from physical to monetary units (Többen, 2017), account for multiple modes of transportation (Llano, Esteban, Pérez, & Pulido, 2010), or address problems of statistical reliability (Schwarm et al., 2006), missing data and mismatching commodity classifications (Park, Gordon, Moore, & Richardson, 2009).
The estimation of interregional trade has also been approached using the gravity model as a framework. At its simplest, the gravity model posits that the intensity of trade between two regions is a function of their respective economic masses and the distance between them. The main difficulty with these methods is that it is not at all clear how to parametrize the model (Simini, González, Maritan, & Barabási, 2012), as that seems to require precisely the type of interregional trade data whose unavailability one is trying to overcome. In applications, model parameters are usually determined through a combination of calibration with freight transport data (Boero, Edwards, & Rivera, 2018;Lindall, Olson, & Alward, 2006), econometric estimation from international trade statistics (Fingleton, Garretsen, & Martin, 2015;Riddington, Gibson, & Anderson, 2006), and guesswork (Johansen, Egging, & Ivanova, 2018;Sargento, Ramos, & Hewings, 2012).
As these approaches have fairly heavy data and labour requirements, some authors have attempted to develop simpler alternatives by extending the logic of non-survey input-output regionalization to a multiregional context (Gallego & Lenzen, 2009;Haddad, Samaniego, Porsse, Ochoa, & de Souza, 2011;Jahn, 2017). In this case, the estimates generated by non-survey techniques are often used in conjunction with a balancing procedure (e.g., the RAS technique) to incorporate additional information available to the analyst.
For all its diversity, the existing literature on interregional trade estimation is skewed towards heavily applied contributions. Overcoming the limitations of the data often requires ad-hoc assumptions and estimation procedures are seldom linked to a broader theoretical framework. As a result, it is not always clear how to generalize beyond the specific data availability conditions under which they were developed, or how to assess the relative strengths and weaknesses of alternative approaches. It adds to the difficulty of drawing lessons from past experiences that-again, due to a lack of suitable data-very few attempts have been made to evaluate the accuracy of the trade estimates (Distefano, Tuninetti, Laio, & Ridolfi, 2019;Fournier Gabela, 2020).
In this paper, the problem of estimating interregional trade is examined through the lens of one specific flavour of the gravity model, the doubly constrained gravity model (DCGM). The DCGM is popular among transportation planners as a tool to estimate (or forecast) OD matrices of passenger and, less commonly, freight flows (Barbosa et al., 2018;Roy & Thill, 2004). Despite the similarity in objectives, that body of work has inspired surprisingly few applications-most notably, Lindall et al. (2006)-and no broader methodological discussion in the context of interregional trade estimation. Conversely, the freight models of transportation planners place great emphasis on issues (e.g., modal split, transshipment and warehousing) that in the context of interregional trade estimation represent nuisances rather than objects of primary interest (Davydenko & Tavasszy, 2013).
Here, we use the term DCGM with reference to a framework with two defining features (Section 2). First, the unobserved bilateral flows to be estimated behave according to some version of the gravity equation. Second, the total flows into and out of each region (e.g., overall supply and use) are known or can at least be estimated. It is worth noting that the DCGMs of the transportation literature have one additional characteristic: they use a specific form of the gravity equation derived via maximum entropy arguments and accordingly calibrate the key parameter of the model to cost data (Wilson, 1967(Wilson, , 1970. In our case, by contrast, the gravity equation is simply assumed to hold. It is stated, however, in a form general enough to accommodate maximum entropy formulations, specifications derived from economic theory, as well as a variety of empirical implementations encountered in the applied trade literature. This flexibility expands the range of empirical techniques that can be used to determine the model's parameters. For example, a gravity specification based on microeconomic principles (e.g., Anderson & Van Wincoop, 2004;Chaney, 2018) can help establish a link between the unknown parameters of the interregional model and those of a behavioral equation that can be estimated econometrically from international trade data. Furthermore, in cases where the analyst has no other option but to make numerical assumptions about parameter values, it ensures that those assumptions can at least be informed by a large body of empirical literature.
This paper contributes to the methodological discussion on interregional trade estimation in four ways. First, it calls attention to the fact that estimation procedures based on the DCGM actually impose weaker restrictions on the data-generating process than it is commonly assumed. Secondly, it shows that DCGM estimation of interregional trade need not be as data and computationally intensive as the few existing applications. In fact, under ideal conditions, a very simple procedure-which we refer to as the "plain vanilla" DCGM estimator-is enough to recover an unknown OD matrix of trade with perfect accuracy. Thirdly, we make a case that several experiences described in the literature, although seemingly unrelated, can be thought of as applications of the DCGM approach. Finally, we examine how the DCGM estimator might behave in real-world applications through a series of Monte Carlo experiments (Section 3). Focusing on the plain vanilla DCGM, we assess how sensitive the bilateral trade estimates are to various types of model misspecification and data quality problems (Section 4). This exercise provides useful indications regarding not only the likely magnitude of the estimation error in applications, but also any systematic biases in the trade estimates. The paper also includes an example of how the plain vanilla DCGM estimator can be operationalized in a real world application (Section 5).

| METHODOLOGICAL FRAMEWORK
Suppose that an economy of interest consists of n > 1 regions and an unspecified number of traded commodities.
Our analysis will only look at one of those commodities (a 'widget'), but the same line of reasoning applies to any other traded product. For i, j = 1,…, n, let x ij denote the overall amount of widgets produced in region i and used in region j during a certain time period of interest. Taken together, this collection of bilateral flows defines an n × n OD 'trade matrix' , X. The analyst's problem is that X, which contains information necessary for applied work, cannot be observed and must therefore be estimated.

| A general doubly-constrained gravity model
As a basic premise, the family of estimation procedures considered here posits that interregional trade behaves according to a gravity model. Specifically, it is assumed that, for all i and j, bilateral widget flows can be written in the form: The terms a i (Á), b j (Á) and f ij (Á)-which are sometimes referred to as origin, destination and separation factors, respectively-represent strictly positive real functions of a vector of explanatory variables, z ij . Note that, under these assumptions, x ij is also strictly positive. As a statement of the gravity model, Equation 1 is more general than the specifications commonly used in applications. For example, consider setting: and Here, the Greek letters denote parameters; the variables s i , u j and d ij measure, respectively, widget supply in the origin region, widget use in the destination region and distance between the two regions. In this case, Equation 1 becomes: In empirical work, gravity model specifications like Equation 5  show it generally provides a good fit for bilateral trade data (e.g., Sargento et al., 2012). Trade models similar to Equation 5 can also be derived from economic theory (Anderson & Van Wincoop, 2004). In our framework, what justification is provided for using an equation like Equation 5 is not important. In fact, it is only necessary that the gravity model hold in the much less restrictive version of Equation 1. Also, we emphasize that Equation 1 is assumed to apply not only to inter-but also to intra-regional flows (i.e., for i = j as well).
With a view to estimating X, a useful starting point is represented by the following special case. With bilateral trade defined according to Equation 1, suppose that the value of the separation factor f ij can be observed for all i, j pairs. Suppose further that the marginal totals of the trade matrix are also available. Thus, region i's total widget production m i* = P j x ij and region j's total widget demand m *j = P i x ij are known for all i and j. Under these assumptions, standard arguments imply that the origin and destination factors-the a i 's and the b j 's of equation 1-can be solved for (up to a constant factor) (Bacharach, 1970;Idel, 2016;Miller & Blair, 2009;Sen & Smith, 1995). Effectively, this means that the true trade matrix X can be recovered exactly. To do so, a simple iterative procedure known as the RAS algorithm can be used. In a nutshell, RAS rescales bi-proportionally a matrix of initial values (the "prior" matrix) so that the resulting matrix (the "scaled" matrix) matches a set of pre-specified row and column totals (the "target" totals). In our case, X emerges as the scaled matrix when RAS is applied to the prior matrix F = [f ij ] with the target totals provided by the m's. A more formal discussion of this point can be found in the Appendix. This idealized setup in which the true trade matrix can be recovered exactly will be referred to as the DCGM.

| Doubly-constrained gravity model estimation
The information necessary to recover X itself is hardly ever available in the real world. In most applied settings, neither the separation factors nor the marginal totals can be observed. Even so, estimates of those quantities should not be difficult to construct. Throughout the paper, we will adopt the following notational convention: we will use a superimposed tilde to distinguish an estimate from the corresponding true value. Thus, suppose that a complete set of estimatesf ij ,m iÃ andm Ãj can be assembled for i,j = 1,…,n. As long as those are reasonably accurate, we should be able to approximate X using the DCGM logic with the estimates in place of the unknown true values. In this spirit, an estimate of x ij can be constructed as: withã i andb j calculated by RAS so that the adding up restrictions P jxij =m iÃ and P ixij =m Ãj are verified for all i and j. Below, this approach to interregional trade estimation is referred to as the DCGM estimator. Under very mild conditions (e.g., if P imiÃ = P jmÃj andf ij ,m iÃ ,m Ãj > 0 ), the DCGM estimateX =x ij Â Ã exists and is unique (Bacharach, 1970). Obviously, iff ij ,m iÃ andm Ãj match their true counterparts for all i and j, we are back in the ideal conditions andX = X.
Consider how the DCGM estimator of Equation 6 may work out in practice. Generally speaking, constructing estimates of the marginal totals is not especially problematic. For the degree of product detail typically used in applications, the literature offers several examples of how to estimate this type of aggregates by combining a country-level input-output (or supply-and-use) table and industry-level data on regional employment (e.g., Eding, Oosterhaven, de Vet, & Nijmeijer, 1999;Jackson, 1998;Kronenberg, 2009;Lahr, 2001;Madsen & Jensen-Butler, 1999, 2005Többen & Kronenberg, 2015;Yamada, 2015). Let us now turn to the separation factors. A natural way to approach the estimation of the fs is by specifying a parametric model. The simplest option is perhaps to postulate that Equation 4 holds, so that trade decays with distance according to a power law. Under this assumption, given that data on the interregional distance d ij are widely available, the fs are effectively known up to the parameter θ (the "distance elasticity"). Leaving aside for the time being the issue of how to obtain an estimate of θ in applied work, suppose that such estimate is indeed available. Then, an estimate of f ij is given by: In passing, it should be noted that this specification of the separation factors requires that distances be strictly positive. This rules out several simple metrics like flat-earth distances between region centroids (or capital cities), as those imply d ij = 0 whenever i = j.
The simple distance decay process of Equation 7 could be enriched with additional explanatory variables.
Suppose, for example, there is reason to believe that there is something systematically different about intraregional trade. This can be accommodated in our framework by assuming that function, that is, I(i = j) takes the value one for intraregional flows and zero in all other cases. Naturally, flexibility comes at the price of increased data requirements, as more elaborate models inevitably have more parameters (e.g., not only the distance elasticity θ, but also the home bias parameter λ) for which estimates have to be obtained.
Irrespective of what assumptions and data sources are employed to construct thefs, oncef ij ,m iÃ andm Ãj are available for all i and j, estimating the trade matrix amounts to a standard application of the RAS algorithm: to findX, one only needs to scale the matrixF =f ij h i to the estimated marginal totals. It is worth emphasizing that DCGM estimation does not require any information nor make any assumptions about the origin and destination factors.
Identifying the functional form, explanatory variables and parameters of the as and the bs in Equation 1 is not necessary.
The remainder of this paper will focus predominantly on the case in which the separation factors are specified according to the power law of Equation 7. To distinguish it within the broader class of DCGM estimators, this approach will be referred to as the "plain vanilla" DCGM estimator. The reason of our emphasis on the plain vanilla estimator is that Equation 7 represents the standard way by which distance decay is modeled in applied trade economics. Thus, its single unknown parameter sits in the middle of a considerable body of theoretical and empirical literature. Yamada (2015). In these analyses, initial estimates of bilateral trade are first predicted using an empirical gravity equation (e.g., one of an econometric nature) and then improved upon using the RAS algorithm.

| Relationship with previous applied studies
In terms of our framework, all those experiences can be described as follows. First of all, it is assumed that Equations 1 to 4 hold, so that bilateral trade flows are as in Equation 5. Then, using estimates to replace the unobserved s i , u j , γ, δ and θ, one computes:g for i,j = 1,…,n. The details of hows i andũ i are specified (e.g., measures of total supply and demand, per capita GDP, etc.) vary across studies. So does the way the values ofγ,δ andθ are determined. Typically, parameter values are selected on the basis of theoretical considerations (Lindall et al., 2006), econometric analyses (Riddington et al., 2006;Yamada, 2015) or educated guesses (Fournier Gabela, 2020; Johansen et al., 2018;Sargento et al., 2012). In some instances,θ is calibrated to freight transport data using an iterative procedure (Lindall et al., 2006;Riddington et al., 2006). Occasionally Equation 8 is augmented with a scaling factor (Fournier Gabela, 2020; Sargento et al., 2012). Sometimes the model includes additional explanatory variables (Johansen et al., 2018;Sargento et al., 2012) or is assumed not to apply to intra-regional flows (Fournier Gabela, 2020; Sargento et al., 2012;Yamada, 2015). We can dispense with such idiosyncrasies here, as they do not affect our main argument. In all cases, the estimation process goes on to use the resulting matrixG =g ij Â Ã as the prior for a RAS procedure. Given a set of marginal totals estimates for X, the RAS algorithm yields estimates of the α's and the β's.
Finally, trade flows are calculated to be: How does this trade matrix estimateX 1 =x 1 ij h i relate to the plain vanilla estimate? LettingX o denote the plain vanilla estimate, it is easy to show that-provided that the same marginal total and distance elasticity estimates are used in the estimation process-the matricesX o andX 1 are equal. From Equations 7 and 8, it is indeed apparent that g ij =sγ iũδ jfij . In other words,G can be obtained by proportionally scaling each row and each column ofF by a strictly positive factor (and vice versa). As long as the target totals do not change, such transformations of the prior matrix do not affect the scaled matrix obtained by RAS (lemma A1 in the Appendix). Hence,X o =X 1 .
As a corollary, this argument implies that the values selected forγ,δ,s i andũ j in Equation 8 do not affect the value ofX 1 . Even if the corresponding true values were known, the accuracy of the trade estimates would not improve. In the same spirit, augmenting the model of Equation 5 with an additional explanatory variable does not affect the trade estimates as long as that variable enters the equation multiplicatively and its value varies only by origin or by destination (e.g., degree of specialization in Sargento et al., 2012). This is true whether or not the variable in question is a meaningful determinant of interregional trade. Similar considerations apply to variables that vary only across product categories (e.g., the tax and margin rates in Johansen et al., 2018). Just likeX o , the estimated trade matrixX 1 is compatible with an entire family of gravity model specifications of which Equation 5 is just a member.
In fact, the relevance of this analysis extends beyond the estimation methods that are explicitly based on the gravity model. Specifically, there are interesting implications for those estimation procedures that, although building on a theoretical framework other than the gravity model, still make use of RAS scaling. For example, one of the techniques Distefano, Tuninetti, et al. (2019) use to reconstruct an unknown bilateral trade network is to apply the RAS algorithm to a prior matrix with generic entry d − 1 ij . This procedure can be conceptualized as the plain vanilla DCGM estimator withθ set equal to 1. In turn, this provides an intuitive interpretation for the RAS scaling factors: they represent estimates of a gravity model's origin and destination factors. More generally, it is common practice for analysts to 'balance' their interregional trade estimates-however obtained-to ensure consistency with the broader accounting framework in which they are embedded. Often, this is done using some version of the RAS method (Gallego & Lenzen, 2009;Haddad et al., 2011;Zhao & Squibb, 2019). When looked at through the lens of the DCGM, such estimation approaches comprise a full, though implicit, specification of the separation factors. In addition, our analysis of Equation 8 implies that seemingly important variables used in the construction of the prebalancing estimates may end up having no effect whatsoever on the post-balancing estimates. In particular, this is the case of any variable whose inclusion only augments the model with a multiplicative term that scales proportionally all flows from a certain origin (likesγ i ) or to a certain destination (likeũδ j ).

| SIMULATION ANALYSIS
Section 2 has shown that, under certain assumptions an unknown OD matrix of interregional trade can be recovered with perfect accuracy using the plain vanilla DCGM estimator. Namely, those assumptions are that: (i) the bilateral flows represented in the trade matrix follow the deterministic gravity model of Equation 1 with the separation factors given by the distance decay function in Equation 4; (ii) the elasticity θ of trade to distance is known; and (iii) the row and column totals of the matrix, the m's, are observed without uncertainty. In real-world applications, it is highly unlikely that any of these assumptions will hold exactly. On the one hand, any empirically feasible model specification can only approximate the actual trade flows. On the other, both the distance elasticity and the marginal totals of the trade matrix will be estimated with some degree of error. How can these departures from the assumptions be expected to reflect on the accuracy of the estimated trade matrix?
This type of sensitivity analysis is complicated by the fact that DCGM estimates are obtained using the RAS method. There is unfortunately no analytical expression either for the scaled matrix or for its derivatives with respect to the prior matrix and the target totals. Thus, to investigate how plain vanilla DCGM estimation might perform in real-world applications, we use a Monte Carlo simulation approach. In each simulation run, we construct a hypothetical trade matrix with known properties and imagine having to recover it from partial information. We examine a variety of scenarios that differ in the type and quality of the information that is presumed available. All analyses are carried out in the R environment for statistical computing (R Core team, 2016).

| Random trade matrix generation
To generate the true OD matrix X, we use a generalized gravity model in the spirit of Equation 1, with trade decaying with distance according to Equation 4. In this case, however, the model is augmented with a random multiplicative error term, ϵ ij . Thus, the flow of widgets from region i to region j is simulated as: for i,j = 1,…,n. Unless explicitly stated otherwise, all of the analyses below assume that n = 25, a realistic choice given the observed administrative structure of many European economies. The as and bs are generated as random draws from a uniform distribution. Note that the arguments of Section 2 imply that modeling the economic determinants of those origin and destination factors explicitly is not necessary. The ϵ's are sampled independently from a lognormal distribution with location parameter zero and shape parameter σ ϵ . The latter varies across simulations. This distributional assumption is broadly consistent with the pattern of residuals generally observed in econometric analyses.
When it comes to generating the ds, we seek to mimic the main features of real-world data using the following procedure. Using a k-d tree approach, a hypothetical square country is randomly partitioned in the desired number of regions. In each partition, a handful of points are then placed at random. These can be thought of as the region's main centres of economic activity. On the basis of this simulated geography, d ij is computed as the mean Euclidean distance between the points in region i and those in region j. The elasticity of trade to distance, θ, is set equal to 0.9, which Disdier and Head (2008) report as the mean estimate in the existing international trade literature. Finally, k is a normalizing constant whose value is chosen so that P i,j x ij = 1. In our experiments, it will be necessary to simulate imperfect knowledge of the trade matrix's row and column totals. Again, the true marginal totals of X are given by m i* = P j x ij and m *j = P i x ij for all i and j. We assume that the analyst cannot observe these quantities, but only the corresponding estimatesm iÃ = m iÃ η iÃ andm Ãj = m jÃ η Ãj , where the ηs represent stochastic multiplicative error terms. In any given region i, the pair (η i* , η *i ) is drawn independently from a bivariate normal distribution with mean one, standard deviation σ η and correlation ρ η . Here, ρ η is presumed to be non-negative, reflecting the fact that-given how estimations are likely to be carried out in in practice-row and column totals that refer to the same region will generally be biased in the same direction. 1 Thems are balanced to ensure that the resulting DCGM estimation problem is well-behaved (i.e., that P imÃi = P jmÃj ).

| Simulation experiments
To assess how the plain vanilla DCGM estimator performs when the assumptions that guarantee its accuracy are relaxed, we conduct four types of analyses. Initially, the empirical significance of each of the assumptions is investigated in isolation. Thus, a first set of simulations examines the implications of misspecifying the distance elasticity, while retaining the assumptions the trade matrix follows a deterministic gravity model and that its marginal totals are known (experiment 1). A second set of simulations considers the consequences of inaccurately estimated row and column totals in the context of a deterministic gravity equation with known distance elasticity (experiment 2). A third set of simulations allows bilateral trade flows to randomly depart from the deterministic gravity equation, but still posits that the distance elasticity and the marginal totals are accurate (experiment 3). Eventually, we try to determine 1 In the existing literature, estimates of unobserved regional-level variables are typically obtained by apportioning a known country total in proportion to some other economic variable that can indeed be observed at the regional level. For example, national widget output might be disaggregated among regions in proportion to employment in the widget industry. Such approaches will tend to overstate economic activity in some regions and to understate it in others. In any given region, it is likely that both the supply and the use of widgets will be distorted in similar ways.
the likely magnitude and pattern of the estimation errors under more realistic scenarios in which several assumptions are violated at the same time (experiment 4). All simulation experiments consist of one thousand replications for each combination of parameter values.
In each simulation run, a trade matrix estimateX has to be compared with its true counterpart X. The difference between corresponding elements,x ij −x ij , will be referred to as a "residual." Following Sargento et al. (2012), we condense the residuals into a single scalar measure of estimation error using the standard total percentage error, 2 Although a useful gauge of overall discrepancy, the STPE does not convey any information as to what trade flow estimates are biased more severely and in what direction. In the analyses below, it is therefore complemented by several other dissimilarity measures. For example, we will often be interested in assessing whether there is a general tendency for intraregional flows to be over-or under-estimated. For that purpose, we will compute the diagonal total percent error, DTPE = P ixii −

| Experiment 1: inaccurate distance elasticity
In applied work, an accurate estimate of the distance elasticity θ may not be easy to obtain. How is the trade estimates' validity affected ifθ differs from the true θ of the gravity model? To explore this issue, we examine what happens when we try to recover a true matrix generated with θ = 0.9 under various alternative choices ofθ. Specifically, the analysis considers values ofθ at 0.1 intervals in the [0.3, 1.5] range. This range contains approximately 90% of the distance elasticity estimates reviewed by Disdier and Head (2008). Throughout, it is assumed that both σ ϵ and σ η are zero. Table 1 describes the distribution over simulation runs of several measures of dissimilarity between the true and the estimated trade matrix. The symbols M and SD respectively refer to the mean and the standard deviation of the simulated distributions. As expected, whenθ is set equal to the actual distance elasticity of the model, the trade matrix is recovered exactly. Indeed, it can be easily verified that in this case the STPE (column (1)) is always zero. As long as the assumed distance elasticity remains reasonably close to the truth, the estimated matrix exhibits tolerable levels of error. Asθ moves further away from θ in either direction, both the mean STPE and its standard deviation increase.
Whetherθ is over-or under-stated, the effects on the STPE are analogous. The corresponding estimates, however, deviate from the true trade matrix in radically different ways. This is apparent in Figure 1, where the same simulated trade matrix is estimated four times for different values ofθ. Each panel plots element-by-element per- × 100 against bilateral distances. In relative terms, a large value ofθ implies a steep distance decay process. Thus, whenθ > θ , the estimated matrix overstates trade between nearby locations and understates trade across larger distances. Forθ < θ, the reverse happens.
A natural way of measuring the intensity of biases of the type displayed in Figure 1 is through the slope of the ordinary least squares regression of PE ij on logd ij . In Table 1, column (2) summarizes its simulated distribution under varying distance elasticity assumptions. Asθ increases, the relationship between the relative estimation error associated with a certain trade flow and the distance between its origin and its destination goes from steeply positive to 2 By construction, the entries ofX add up to the same grand total as those of X. Each widget thatX assigns to a wrong origin-destination pair contributes +2 to the numerator of the STPE. In this sense, the STPE can be interpreted as twice the share of total trade that has been allocated to an incorrect link. steeply negative. In column (3), the R-squared of the regression suggests that the relationship is relatively tight at all levels ofθ. From these results, it is clear that the single-simulation exercise of Figure 1 is just an instance of a much more general pattern.
The implications for intra-regional trade estimation are worth noting. By nature, intra-regional flows take place over comparatively short distances. In the simple setup of this simulations, this means that wheneverθ is specified incorrectly, the resulting intra-regional trade estimates are all biased in the same direction: upwards when the F I G U R E 1 Percent estimation error versus distance for a sample 50-region simulation run under varying distance elasticity assumptions distance elasticity is set too high and downwards when it is set too low. This can indeed be seen in column (4), which reports what share of the trade matrix's diagonal entries is overestimated for varying choices ofθ.
Finally, we examine how the results of Table 1 would be affected if a number of regions other than n = 25 were assumed. This sensitivity analysis considers two alternative values of n, namely 10 and 50. For reasons of space, the full results of this exercise are only displayed in the supplemental materials (Table S1). As expected, in qualitative terms they are entirely analogous to those of Table 1. From a quantitative perspective-everything else being the same-the bias associated with an inaccurate choice ofθ becomes more severe as the number of regions grows. This makes intuitive sense: an increase in n makes the problem of recovering the trade matrix more difficult, as the number of unknown matrix entries rises more rapidly than the number of constraints provided by the row and column totals. Loosely speaking, when the number of regions becomes larger, the estimator leans more heavily on the distance elasticity parameter.
As the number of regions in the system grows, how fast does the accuracy cost of a bad of choice ofθ increase?
Suppose, for example, thatθ is set equal to 1.2. If n = 10, this leads on average to a STPE of 15% (with SD 2.1). By contrast, when n = 50 the average STPE becomes 19 (SD 1.5). It is tempting to interpret these figures as indicating that, even in the face of a considerable increase in the number of regions, the accuracy of the estimates deteriorates only modestly. Once again, however, focusing solely on the STPE misses the fact that a bad choice ofθ biases introduces a very systematic source of bias. Consider, for instance, the estimation error associated with total intraregional trade as summarized by the DTPE. Withθ = 1:2, when the system consists of 10 regions, intraregional trade is overestimated by 20% (SD 2.2) on average. When the number of regions increases to 50, the average DTPE becomes as large as 46% (SD 2.4). Abstracting from this specific example, Figure 2 represents the simulated distribution of the DTPE for various combinations ofθ and n. It is indeed apparent that intra-regional trade estimates are biased much more severely by distance elasticity misspecification if the number of regions is large than if it is small.

| Experiment 2: measurement error in the marginal totals
In reality, the row and column totals of the target matrix are generally unknown. Instead, they have to be estimated indirectly from the available data. Thus, one will typically have to work with target totals that contain at F I G U R E 2 Bias in estimated intraregional trade for various choices ofθ and n least some amount of measurement error. How will this reflect on the accuracy of the estimated trade matrix? In our simulation framework, the degree of reliability of the available marginal totals is controlled by the parameter σ η . We assess the repercussions of measurement error in row and column totals for four alternative values of σ η , namely 0.025, 0.05, 0.10 and 0.15. Given that the ηs are normally distributed, setting σ η equal to, say, 0.05 implies that approximately 95% of the time the estimated marginal totals will be off by up to 10%. Thus, our choices for σ η denote degrees of measurement error that range from modest to severe. In this analysis, it is assumed that the distance elasticity is known (i.e.,θ = 0:9) and that the underlying gravity model is deterministic (σ ϵ = 0). Initially, ρ η is set equal to zero.
The main results are presented in Table 2. Naturally, as increasing amounts of measurement error are introduced in the marginal totals, the estimated matrix gradually departs from its true counterpart. Thus, the simulated STPE distribution (column (1)) is centred at 3% (SD 0.3) or σ η = 0.025 and at 17% (SD 2.2) for σ η = 0.15. At the elementby-element level, inaccurate row and column totals give rise to a distinctive pattern of residuals. Depending on its sign, measurement error in a certain marginal total tends to push or pull all entries of the corresponding row or column in the same direction. Evidence of this can be found in column (2), which describes the simulated distribution of the 'row-wise rate of error concordance' . By row-wise rate of error concordance we refer to the percent share of the trade matrix's entries for which the estimated flowx ij is biased in the same direction as the estimated totalm iÃ of the row it belongs to. Formally, it is given by F I G U R E 3 Percent error heatmap for a sample sevenregion simulation run respectively denoting the indicator and the sign function. On average, the rate of concordance is about three quarters. The tendency for measurement error in marginal totals to distort entire rows and columns of the matrix in a similar way is also visible in Figure 3, which uses a heatmap to represent the percent errors associated with the entries of a trade matrix estimate for a sample 7-region simulation with σ η = 0.10.
Figure 3 also helps to visualize what can happen when the measurement error in the total of a row total is positively correlated with the measurement error in the total of the corresponding column. Consider, for example, region 1. In this case, the row and the column totals are overestimated by 12 and 7%, respectively. As a result, the element that lies at the intersection of row 1 with column (1) is overstated by a non-negligible 16%. Something similar occurs in region 6, where intraregional trade is underestimated by 16%. Figure 3 was generated under the assumption that ρ η = 0.9. As ρ η approaches unity, it becomes increasingly likely that a large positive (or negative) error in any given row total will be matched by another large error of the same sign in the corresponding column total. This, in turn, tends to distort the elements that lie along the diagonal of the calibrated matrix.
The general pattern can be appreciated from Figure 4. For various values of σ η , the graph displays how the distribution of the standard percentage error along the diagonal, STPED = P ixii − x ii j j = P i x ii À Á × 100 , tends to move upwards as ρ η becomes larger. For less than extreme levels of measurement error, however, the effect seems far from dramatic.
Finally, sensitivity analysis with respect to n (table S2 in the supplemental materials) suggests that the results of this section are fairly unresponsive to the number of regions in the system.

| Experiment 3: idiosyncratic errors
Until now, σ ϵ has remained fixed at zero in all simulation analyses. Effectively, this forces trade to adhere to a deterministic gravity model. What happens if interregional trade follows the gravity equation only approximately? Assuming that both the distance elasticity and the marginal totals are specified correctly (i.e.,θ = 0:9 and σ η = 0), we run three sets of simulations using different values of σ ϵ . In each case, we select σ ϵ such that fitting a log-linearized version of Equation 10 to the simulated data by ordinary least squares (i.e., regressing logx ij on loga i , logb j and logd ij ) yields a certain pre-specified R-squared (in expectation). The target R-squared values that define our three scenarios are 0.9, 0.8 and 0.7. This range of values is compatible with the typical findings of econometric analyses. In each case, the average STPE is 34% (SD 3.8), 51% (SD 5.8) and 65% (SD 7.5), respectively. As expected, with the error term ϵ ij independent across origin-destination pairs, there is no noticeable pattern in the residuals. The idiosyncratic error acts as a potentially loud but random source of noise.
F I G U R E 4 Correlated measurement error in marginal totals and bias in intraregional trade estimates

| Experiment 4: realistic scenarios
In applications, one is typically faced with a situation in which-all at the same time-the gravity equation describes interregional trade only loosely, and both the distance elasticity and the marginal totals are uncertain. To assess what degree of accuracy should be expected of the plain vanilla DCGM estimator in such cases, we examine three scenarios characterized by empirical problems of increasing severity. The first scenario (labelled LOW for convenience) posits moderate idiosyncratic error variance (σ ϵ such that the expected R-square of the log-linear model is 0.9) and known marginal totals (σ η = 0). In the second scenario (MEDIUM), σ ϵ is increased to make the expected R-squared 0.8 and a certain amount of measurement error is introduced in the marginal totals (σ η = 0.05 and ρ η = 0.6). Finally, the most extreme scenario (HIGH) assumes σ η = 0.10, ρ η = 0.9 and an expected R-squared of 0.7. In each case, several values ofθ are contemplated.
The distribution of v error measures over simulation runs is summarized in Table 3. Predictably, the overall accuracy of the estimates deteriorates substantially as empirical difficulties become more pervasive. For any given choice of the distance elasticity parameter, the STPE (column (1)) increases substantially moving from the LOW to the HIGH scenario. Misspecification of the distance elasticity tends to distort the relative sizes of long-versus short-distance trade flows in ways similar to what has been observed in the deterministic setting of experiment 1. This is typically be reflected in the DTPE (column (2)): in general, intraregional trade is underestimated whenθ is set too low and overestimated when it is set too high. For a given choice ofθ, the average DTPE over replications is indeed not far from its counterpart in Table 1. The standard deviations, however, indicate that there is now quite a significant amount of variation across simulation runs, especially in the more challenging scenarios. When it comes to measurement error in the marginal totals, the row-wise rate of error concordance in column (3) imply its repercussions on the accuracy of our estimates are no longer as clear-cut as in the earlier deterministic simulations. It is still true that, when a certain margin is measured inaccurately, each element in the corresponding row ofX gets dragged in the T A B L E 3 Simulated distribution of some dissimilarity metrics in various realistic scenarios same direction. However, whether that element turns out to over or underestimate the corresponding entry of X also depends on the magnitude of the randomly drawn error term associated with it.
In this case as well, the analysis was repeated for n = 10 and n = 50 (supplemental materials, Table S3). Generally speaking, the main results of this section hold up. Also, changes in n appear interact with the choice ofθ in a way consistent with the findings of experiment 1. Finally, the range of STPE values observed in this simulation exercise seems compatible with the findings of a series of conceptually similar analyses carried out by Sargento et al. (2012) on real-world data.
With a view to applying the plain vanilla DCGM to real world problems, these results suggest that-unless the functional form of the gravity model is particularly ill-suited to represent the unknown trade flows we are trying to work out-the approach is indeed capable of producing fairly accurate estimates. It is crucial, however, that we are able to initialize the estimation algorithm with a good approximation of the true distance elasticity, the more so if the number of regions in the system is large.

| DCGM ESTIMATION WITH REAL DATA: AN EXAMPLE
To demonstrate how DCGM estimation may play out in a real world setting, this section presents an example based on data from the World Input-Output Database (WIOD) (Timmer, Dietzenbacher, Los, Stehrer, & De Vries, 2015). Produced by harmonizing macroeconomic information from a variety of sources, WIOD is an internally coherent multi-country accounting framework detailing the economic interdependences among 56 branches of economic activity ("products," for short) and more than 40 geographical entities ("countries"). Given the general lack of data at the sub-national level, it is not uncommon for methodological research of the type conducted here to experiment with international trade datasets like WIOD (e.g., Lamonica & Chelli, 2018;Sargento et al., 2012).
For our purposes, WIOD presents two advantages over raw trade statistics (e.g., UN Comtrade): first, it contains no trade asymmetries, as all discrepancies between mirror flows were removed in the harmonization process (Distefano, Chiarotti, Laio, & Ridolfi, 2019;Gehlhar & Pick, 2002); second, it includes information on how much each country trades with itself.
Our starting point is the WIOD world input-output table for 2014. We aggregate it into a series of productspecific OD matrices of international trade. Only one of those matrices will be considered here, namely, the one relating to agricultural products (identified in WIOD by the code A01). The same estimation procedure, however, can be applied separately and independently to each of the remaining products. Also, the scope of the analysis will be limited to EU member states (including the United Kingdom), as they are all part of the same customs union. Thus, the trade matrix at the centre of this example is 28 × 28. We will treat it as if it were the unknown true X matrix that needs to be estimated.
Once again, we focus on the plain vanilla flavour of the DCGM. Implementing the estimator requires three types of inputs. First, we need information about the row and column totals of X. Here, we will use the true marginal totals, as in the simplified setting of this example they are known. In more realistic circumstances, the marginal totals would have to be estimated beforehand. Second, we necessitate data on bilateral distances between EU member states.
For that, we turn to the GeoDist database (Mayer & Zignago, 2011). Specifically, as our measure of distance d ij we use the GeoDist variable distwces, which is a population-weighted aggregation of bilateral distances between major cities. As required, distwces is strictly positive (including for i = j). Finally, the third necessary piece of information to implement our estimator is a distance elasticity estimateθ. Initially we setθ = 1 like Distefano, Tuninetti, et al. (2019), Fournier Gabela (2020), Johansen et al. (2018) and Sargento et al. (2012). Accordingly, the prior matrixF is computed asf ij = d − 1 ij . Using the RAS algorithm, we scaleF to the target marginal totals and obtain our plain vanilla estimate of X as a result. Let this estimated trade matrix be denotedX jθ = 1 . Because in this illustrative example the true matrix X is known, we can actually assess how closely our estimates approximate the true values from WIOD. The STPE associated withX jθ = 1 turns out to be a disappointing 87%. What is wrong with the estimates is apparent from the left-hand panel of Figure 5, which plots corresponding entries of X andX jθ = 1 against each other. Clearly,X jθ = 1 has a marked tendency to overestimate the small flows and underestimate the large ones. In addition, the DTPE indicates that total intra-national trade is underestimated by about a half.
In light of the earlier simulation analysis, this pattern of error strongly suggests that we have setθ too low. When that happens, our estimator has indeed been observed to artificially inflate the flows between distant locations (which, all things equal, the gravity model predicts to be relatively small) at the expense of those between shorter distances (which the model predicts to be relatively large).
The hypothesis that settingθ = 1 understates how steeply trade declines with distance in our true trade matrix can actually be tested. To this end, we econometrically fit the empirical counterpart of Equation 5 to the entries of X. (More information on this analysis can be found in the supplemental materials.) Estimation by log-linear ordinary least squares yields a distance elasticity parameter of 1.85 (with a robust standard error of 0.067), which is indeed well in excess of one. Incidentally, the R-squared of the regression is 0.78, which reassures us that even a very simple specification like Equation 5 captures the main features of the true trade matrix data reasonably well. 3 How would the plain vanilla DCGM estimator perform if, instead of assuming thatθ = 1, we were able to guess the distance elasticity of agricultural trade more closely? Suppose, for example, that we setθ = 2. The resulting trade matrix estimate,X jθ = 2 , has an STPE of 22%. This represents a remarkable improvement onX jθ = 1 . As an aside, it also compares favourably with the results of a similar exercise in Sargento et al. (2012). 4 As can be seen in the right-hand panel of Figure 5, comparatively small flows can still be substantially misrepresented inX jθ = 2 , but there is no systematic pattern to the estimation error. On the other hand, the largest flows are estimated with remarkable accuracy.
Many of those larger flows relate to within-country trade and correspond to the main diagonal of the X matrix.
F I G U R E 5 Estimated versus actual bilateral trade in agricultural products in the EU (2014, million US dollars) Finally, we consider what happens if we initialize the plain vanilla DCGM using a distance elasticity estimate obtained econometrically from the true trade matrix itself. Hence, on the basis of our earlier regression analysis, we setθ = 1:85. The STPE and DTPE of the ensuing estimate are 27 and −11%, respectively. In other words, the accuracy of the estimates does not improve relative toX jθ = 2 .
By and large, these findings are consistent with our earlier simulation analysis. At least within the context of this example-and as long as one is able to make an informed guess about the relevant distance elasticity-the plain vanilla DCGM estimator recovers the unknown trade matrix with a degree of accuracy that would be deemed sufficient for most applications.

| CONCLUDING REMARKS
As a rule, national statistical programmes do not collect data on the amount of goods and services traded among the various regions of a country. The resulting dearth of interregional trade data often constitutes a major obstacle to the construction of policy analysis models with sub-national detail. This paper has examined a family of methods for estimating interregional trade on the basis of a DCGM gravity model. Our analysis has highlighted several appealing characteristics of the approach. First, the DCGM estimator makes only mild assumptions about the causal mechanism behind interregional trade. The gravity equation that lies at its core-one of the most robust empirical regularities in economics-can be justified using a number of alternative economic (Baldwin & Taglioni, 2006) or statistical (Wilson, 1967) arguments. There has even been speculation that "just about any plausible model of trade would yield something like the gravity equation" (Deardorff, 1998, p.12). That the interregional trade estimation process remains as agnostic as possible about the nature of the underlying datagenerating process is indeed desirable. The more assumptions we introduce in the estimation, the more likely it becomes that some of those assumptions will conflict with those of the policy analysis models in which we intend to use the ensuing estimates. In our framework, the gravity equation is assumed to hold in a very general (i.e., unrestrictive) form and there is no need for the analyst to specify-let alone measure-all the economic drivers of interregional trade.
Another attractive feature of DCGM methods is represented by their fairly modest data and computational needs. In practice, estimating the OD matrix of interregional trade for a certain product only amounts to a standard application of the RAS algorithm. Typically, the constraints on the row and column totals can be constructed from readily available information using standard methods. The prior matrix is also very easy to obtain, at least in an uncomplicated version of the model. The simplest implementation of the estimator, the plain vanilla DCGM, assumes that interregional trade decays with distance according to a power law. In this case, only data on pairwise interregional distances and an estimate of the distance elasticity of trade are necessary to compute the starting values for RAS. Such data availability and processing requirements are not significantly more demanding than those of common input-output regionalization methods. In spite of its simplicity, the plain vanilla DCGM estimator is demonstrably equivalent to several seemingly more elaborate approaches described in the literature.
Finally, a simulation study has been conducted to assess how well the plain vanilla DCGM might perform in realworld applications. Using Monte Carlo techniques, we have examined the potential implications of various types of model misspecification and data quality problems. Overall, the method seems to produce reasonably accurate estimates under a range of realistic circumstances, even though-especially when the number of regions in the system is large-getting the crucial distance elasticity parameterθ badly wrong can cause severe amounts of bias. A demonstrative example based on actual data has given analogous findings.
Admittedly, our analysis does not discuss how to identify the appropriate value ofθ in applications. Assessing the relative merits of alternative approaches to the estimation ofθ , and more generally finding satisfying ways of specifying the distance decay process for applied DCGMs, will require additional work. Even so, DCGM estimators appear to be useful tools for cost-effectively estimating interregional trade at the typically coarse level of territorial and commodity aggregation that is typical of input-output and computable general equilibrium analysis.

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section at the end of this article.

APPENDIX A
This appendix states and proves an invariance property of the RAS algorithm that lies at the root of various assertions made in the methodological discussion of Section 2. Below, we use angled brackets hÁi to denote a square matrix that has the relevant vector along the main diagonal and zeros everywhere else.
Lemma A1. Consider a generic nonnegative m × n matrix, A o . Let v and w be two strictly positive vectors of appropriate length and define a second m × n matrix as A 1 = hviA o hwi. Suppose that, using the RAS algorithm, we scale the matrix A o to match a target m-vector of row sums, m R , and a target n-vector of column totals, m C .
Assume that this scaling problem is feasible and let the resulting scaled matrix be denoted B. Then, the same scaled matrix B is obtained if A 1 is used as the prior for RAS instead of A o .
One way of deriving lemma A1 is a follows. First, consider the problem of scaling A o to row sums m R and column sums m C by RAS. A solution takes the form of a pair of suitably sized vectors r o and s o such that the scaled matrix B = hr o iA o hs o i satisfies both Bi = m R and B 0 i = m C . The notation i refers to a vector of ones with the appropriate length.
If they exist, such r o and s o are unique up to a multiplicative constant (Bacharach, 1970). Now suppose that, instead of A o , we use A 1 as the prior for RAS. The target margins m R and m C , on the other hand, are left unchanged. It is easy to verify that a pair of scaling vectors that solves the modified problem is given by r 1 = hvi −1 r o and s 1 = hwi −1 s o . Indeed, hr 1 iA 1 hs 1 i = hr o ihvi −1 hviA o hwihwi −1 hs o i = B, which is known to match the prescribed marginal totals. Just like before, the solution given by r 1 and s 1 is unique up to a constant factor. Thus, for target margins m R and m C , the scaled matrix obtained through the RAS algorithm is B irrespective of whether A o or A 1 is used as the prior. Lemma A1, for example, justifies the claim in subsection 2.1 that X can be computed via the RAS algorithm using F as the prior matrix and the m's as the target totals. To see this, express Equation 1 in matrix form and rearrange to obtain F = hai −1 Xhbi −1 , with a = [a i ] and b = [b j ]. Trivially, scaling X to its own marginal totals (i.e., the m's) returns X itself as the scaled matrix. Then, scaling F to those same marginal totals must also yield X.
The equivalence betweenF andG as priors for RAS in subsection 2.3 is also a direct implication of lemma A1.