Policy Research Working Paper 8985 Generating Gridded Agricultural Gross Domestic Product for Brazil A Comparison of Methodologies Timothy S. Thomas Liangzhi You Ulrike Wood-Sichra Yating Ru Brian Blankespoor Erwin Kalvelagen Development Economics Development Data Group August 2019 Policy Research Working Paper 8985 Abstract This paper examines two new methods to generate gridded regression approach for all the key crops and contributors to agricultural Gross Domestic Product (GDP) and com- agricultural GDP. However, the issue of degrees of freedom pares the results with a traditional method. In the case of is an important limiting factor, as the approach requires Brazil, these two new methods of spatial disaggregation sufficient subnational data. The cross-entropy method with and cross-entropy outperform the prediction of agricultural readily available spatially distributed crop, livestock, forest, GDP from the traditional method that distributes agricul- and fish allocation far outperforms the traditional method, tural GDP using rural population. The paper finds that the at least in the case of Brazil, and can operate with national- best prediction method is spatial disaggregation using a and/or subnational-level data. This paper is a product of the Development Data Group, Development Economics. It is part of a larger effort by the World Bank to provide open access to its research and make a contribution to development policy discussions around the world. Policy Research Working Papers are also posted on the Web at http://www.worldbank.org/prwp. The authors may be contacted at tim.thomas@cgiar.org (corresponding author) and bblankespoor@worldbank.org. The Policy Research Working Paper Series disseminates the findings of work in progress to encourage the exchange of ideas about development issues. An objective of the series is to get the findings out quickly, even if the presentations are less than fully polished. The papers carry the names of the authors and should be cited accordingly. The findings, interpretations, and conclusions expressed in this paper are entirely those of the authors. They do not necessarily represent the views of the International Bank for Reconstruction and Development/World Bank and its affiliated organizations, or those of the Executive Directors of the World Bank or the governments they represent. Produced by the Research Support Team Generating Gridded Agricultural Gross Domestic  Product for Brazil: A Comparison of Methodologies1  Timothy S. Thomas2  Liangzhi You3  Ulrike Wood‐Sichra3  Yating Ru3  Brian Blankespoor4  Erwin Kalvelagen3  Keywords: agriculture, regional economic activity, gridded GDP, cross‐entropy, spatial disaggregation  JEL: Q10, R12  1  The authors thank Juha Siikämaki (IUCN), Tim Robinson (FAO), Jean‐François Arvis (World Bank), Uwe Deichmann  (World Bank) and Siobhan Murray (World Bank) for earlier discussions. We also thank the participants of the  Workshop on Data Integration: Realising the Potential of Statistical and Geospatial Data that took place in  Belgrade, Serbia and was organised in the context of the partnership agreement between UNECE and UN‐GGIM:  Europe, with the support of Eurostat and EFTA. We appreciate Juha Siikämaki sharing the non‐timber data with us.  We appreciate the support of our institutions and financing from DFID under the World Bank’s Strategic Research  Program ‘Big Data’. The findings, interpretations, and conclusions expressed in this paper are entirely those of the  authors. They do not necessarily represent the views of the International Bank for Reconstruction and  Development/World Bank and its affiliated organizations, or those of the Executive Directors of the World Bank or  the governments they represent.   2  Corresponding author. International Food Policy Research Institute (IFPRI), Email: tim.thomas@cgiar.org  3  International Food Policy Research Institute (IFPRI)  4  Development Data Group, World Bank, Email: bblankespoor@worldbank.org  Introduction  Doll, Muller, and Elvidge (2000) were the first to use satellite data to spatially distribute GDP on a  geographic grid. Their methodology relied on night‐time lights, which generally gave an indication of  where cities were, while the intensity of the light gave an indication of both population density and  intensity of production, with the more intense lights being associated with higher GDP per unit area. The  methodology they used was especially helpful in countries for which finely distributed subnational data  on economic activities were not available. It allowed researchers and policy makers to have a better idea  as to which cities and locales were of high importance to the national economy, and which were under‐ developed.  However, this methodology performed poorly for countries with relatively large agricultural GDPs, since  the night‐time lights were more of an indication of manufacturing and urban population rather than  farming and rural population. In 2006, Nordhaus (2006) proposed using population to distribute GDP. If  GDP per capita were relatively evenly distributed within administrative areas in which the values were  aggregated, this could provide a reasonably accurate estimate of the geographic distribution of GDP.  Several authors since then have used this idea and similar methodology to distribute agriculture GDP  using rural population5 or land cover (e.g. Ghosh, Powell, and Elvidge, et al., 2010; Gunasekera et al.  2015; Rafa, Moyer, and Wang et al., 2017).   In this paper, we propose two new methods for distributing agricultural GDP across a geographic grid.  The first is a spatial allocation methodology which relies on several gridded data sets showing location of  agricultural production. It starts with several components of agricultural GDP such as crop, livestock,  fishery and timber production to provide a prior estimate of gridded agricultural GDP. The estimate is  then reconciled with national and sub‐national agricultural GDP values using an entropy‐based data  fusion method. The second is a “spatial disaggregation” method that relies on imposing a data‐ generating process at the gridcell level with the variable of interest (in our present case, agricultural  GDP, though the technique can be applied to many types of questions) – which is not measured at the  gridcell level ‐‐ being defined by a combination of known gridded explanatory variables and a set of  parameters to be estimated. The equations are aggregated up to the level at which the variable of  interest is known, and a multiple regression is used to solve for the parameters, which in turn can be  used to make predictions at the gridcell level.  In this paper, we use Brazilian agricultural GDP data to test these two new techniques against the  traditional technique of using rural population to distribute the agricultural GDP. We have data available  at the municipio level, with 5,564 units. However, we instead use the data at a more aggregated  microregion level, with 558 units. Then we make predictions at the gridcell level and aggregate up to the  municipio, at which level we compute two measures of predictive accuracy.  Review of Previous Work  Previous analyses use two primary geographically‐referenced data types to locate GDP. The one that has  been used the longest is nighttime light emissions captured by satellite imagery (Doll, Muller, and                                                               5  It should be pointed out that gridded population data sets are technically not satellite data sets but GIS raster  (gridded) data sets, but they may include satellite data sets as inputs to their construction.  2    Elvidge 2000; Sutton and Costanza 2002; Ebener, Murray, Tandon et al. 2005,  Doll, Muller, and Morley  2006; Henderson, Storeygard, and Weil 2012; Sutton, Elvidge, and Ghosh 2007; Chen and Nordhaus,  2011; Ghosh, Anderson, and Powell et al. 2009). Night lights are seen as indicators of where people –  especially urban people – are concentrated, and more importantly, where economic activity (particularly  non‐agricultural) is concentrated.   The second is population. Nightlights seemed to capture industrial and commercial activity reasonably  well, but it performed poorly in capturing significant portions of total GDP in countries that have  relatively large agricultural sectors.   Both types of methods use national or subnational GDP measures, and distribute the measures to each  pixel based on the intensity of the indicator: night lights or population density. In the latter case, the  assumption is normally that GDP per capita is uniform across people in the national or subnational unit  for which data are available.  Often the distribution happens using separate measures for rural and non‐rural. If the two are  distinguished, often nightlights are used with non‐rural GDP distribution while population is used to  distribute to rural areas the agricultural GDP. Because this is a mechanistic approach, in whatever way  the population data set has an error, there will be a corresponding error in the gridded agricultural GDP.  Furthermore, in whatever way the true distribution deviates from the assumption of perfectly even  distribution of agricultural GDP per person this will also produce errors.   Nordhaus (2006) uses subnational GDP data to distribute the GDP based on population per pixel. He  then does post‐distribution analysis at a pixel level with climate and geographic variables. He does not  distinguish between urban and rural populations, even though data for most countries show a gap  between rural and urban income. Kummu, Taka, and Guillaume (2018) use a similar methodology, and  Yetman, Gaffin, and Xing (2004) use it to generate a future gridded GDP based on a future gridded  population. The method to grid GDP by UNEP and the World Bank distinguishes between rural and  urban population assuming the latter to have a higher GDP per capita and allocates GDP using an  estimated ratio of urban to rural GDP per capita by country (UNISDR 2011).6  Ghosh, Powell, Elvidge, et al. (2010) use LandScan population grids (e.g. Dobson et al. 2000) to distribute  national level agricultural GDP into the rural area, then use night lights to distribute the remaining GDP  to urban areas. Rafa, Moyer, and Wang et al. (2017) also use a hybrid method for distributing rural and  urban GDP in Uganda. Likewise, Gunasekera et al. (2015) use separate models for the non‐agricultural  and agricultural sectors with refinements of the non‐agricultural model in Blanchard et al. (mimeo).   Turning to the first of two new methodologies for generating gridded agricultural GDP, we note that  previous work of spatial allocation models includes GDP models (Ghosh et al. 2010; Nordhaus 2006;  UNISDR 2011; and Gunasekera et al. 2015), population models (Tobler et al. 1997) and agricultural  models such as MapSPAM (You et al. 2006) and Global Agro‐Ecological Zones (GAEZ) (FAO/IIASA 2012).  GAEZ has a total crop production data set for the year 2000 that combines actual yield and production  and FAO crop statistics (international) prices for 23 major commodities.                                                                6   Available at Global Risk Data Platform: http://preview.grid.unep.ch.  3    Methodology  Distribution Based on Rural Population  The methodology for distributing rural population used in this analysis is relatively straightforward. First,  as with all the methodologies compared here, the agricultural GDP is aggregated from the municipio  level up to the microregion level.  Then, the rural population grid was produced using the Gridded  Population of the World (GPW) version 4 (CIESIN 2017) for 2010 adjusted to the United Nation’s World  Population Prospects (UN WPP), and masking out urban area using the Global Human Settlement grid  for 2015 (Pesaresi and Freire 2016).7    Once these steps are done, it is straightforward to aggregate population to the 5 arc minute pixel used  by all three methodologies, then for each microregion, distribute the agricultural GDP per rural person  to each pixel based on the rural population in each pixel.  Spatial Allocation Using Cross‐Entropy  The methodological framework is a spatial allocation model that estimates the subnational agricultural  GDP to 5 arc minute pixels (or any other spatial resolution desired). The probabilistic model includes the  control (constraint) of national or subnational agricultural GDP available from various sources of  national ministries by country or reports. The five components of agricultural GDP (crops, livestock,  fisheries, timber, and non‐timber forest products) are used to provide a prior estimate of agricultural  GDP at the pixel level.   We define our spatial AgGDP allocation problem in a cross‐entropy framework following You et al.  (2014). The first step is to transform all real‐value parameters into a corresponding probability form.  Let  Si be the percentage of the total agricultural GDP allocated to pixel i within a certain country (e.g.  Country X) except the crop value for that pixel. AgGDPi is the agricultural GDP allocated to pixel i in  Country X while TotalAgGDP is the total agricultural GDP for that country.  Therefore:    Let PreAlloci be the prior allocation of AgGDP share from our best guess. One possibility for the first  approximation is from pixel level components of AgGDP where we have gridded data sets:  ℎ   As we have the major components of AgGDP, and the sum of these component should be close to the  official value. At least, we make sure that the official AgGDP value (TotalAgGDP) should be no less than  the sum of all five components of agricultural GDP. Therefore, we first sum them up    ∈                                                              7  We select GHS Settlement grid as an example of a geospatial definition of urban population to exclude from the  model.  The definition of urban is an active topic of research (e.g. Dijkstra and Poelmann 2014; Roberts et al. 2017;  Angel et al. 2018; Bosker et al. 2018).  4    Then rescale the prior Ag GDP to be consistent with the official AgGDP value:    Finally, we could calculate our prior for Si in a probability form by normalizing PriorAgGDP:    ∑∈ We could calculate livestock, forestry, fishing and hunting values by their productions multiplied by the  corresponding prices, either national average price or better local farmgate prices. Then, we could  formulate our cross‐entropy model in the following mathematical optimization framework:  log log   Subject to  1  ∀  ∈ 0 1 ∀   Where  i: i=1,2,3,… pixel identifier within the allocation unit (e.g. Brazil)  k: k=1,2,3, … identifiers for sub‐national geopolitical units (e.g. state) where AgGDP values  (SubAgGDPk) are available  The objective function of the spatial allocation model is the cross‐entropy of AgGDP shares and their  prior. The first equation of the constraints is the adding‐up constraint for AgGDP which simply specifies  that the sum of all allocated AgGDP values is equal to the total AgGDP of the country. The second  equation sets the sum of all allocated AgGDP within those subnational units with available data to be  equal to the corresponding subnational AgGDP values. The last equation is simply the natural constraint  of the percentage of AgGDP, which is the probability in the cross‐entropy model.  The modeling framework is flexible to add more constraints if the data exist or we have some  reasonable assumptions on how AgGDP should be spatially disaggregated (e.g. population density or  market access may play a role in determining the spatial distribution or spatial structure of AgGDP).   After the model run, we need to add back crop value to get the pixel level AgGDP. The final AgGDP by  pixel is calculated by:    5    Spatial Disaggregation  The idea of spatial disaggregation is relatively simple. We assume there is a data‐generating process at  the fine spatial level (pixel), and that we only have data on that value at a more aggregated level. Even  though we do not have a value (i.e., a dependent variable for the process) at the disaggregated level, we  have data for variables that are part of the data‐generating process at that level which interact with a  set of parameters with undetermined values. The sum of the estimates of the statistic at the  disaggregated level must equal the aggregated value. Therefore, we estimate the parameters used at  the disaggregated level that give the best fit at the aggregated level, based on some criterion such as  minimizing the sum of squared errors or maximizing a likelihood function.  We can make the idea more explicit. Let i index the aggregated units, and yi be the statistic of interest. In  our case, yi would be the agricultural GDP in municipio i.  Let j be a spatial unit inside i. In our case, we  will use 5 arc minute pixels as the unit, but there is no reason that the unit would have to be a regular  shape or of uniform size. Variables of interest are placed in the vector Xij. The disaggregated statistic of  interest would be given by     (1)  where e, the residual, is a measure of influence of unmeasured variables that we assume has a mean of  0.  We can sum over all j for each i, to get    which can be written as  ∗ ∗   (2)  where  ≡ , ≡ , ≡   and where Xi* is simply Xi with the first column of ones (the constant term) removed, and * is simply   with the first parameter associated with the constant removed.  While the first equality in equation (2) is written like an ordinary least squares (OLS) regression, there  are a few key differences to be noted. First, assuming the disaggregated regression has a constant term  in it, Xi does not. Instead of a constant for each aggregated unit i, there is now a value Ni, since that  many disaggregated units were added together. Second, while adding together a number of stochastic  elements with mean 0 still gives mean 0, the variance of each summed ei can become quite complicated.   Even in the easiest case – where each eij is independent of all of the other residuals – the variance is still  Ni2, which means that there is a special form of heteroskedasticity that needs to be corrected for,8                                                               8  In Stata, this is done simply by using the “aweights” option together with the number of disaggregated units  within each aggregated unit, which here we are calling Ni.  6    unless Ni = Nj for all i and j, in which case there is no need to make any correction, because it is  homoscedastic.  Controlling for spatial relationships  Because of the nature of data based on location, there is often some kind of spatial correlation that e  affects the statistical analysis. Sometimes there are omitted variables that might be spatially correlated,  or the variables that are included might have some type of spatial error themselves due to the locations  not being perfectly known. And it could very well be that the dependent variable is influenced directly  by the neighboring dependent variables.  It is possible – though far from easy ‐‐ to control for two kinds of spatial relationships at the level of the  pixel. The first is a spatial lag in which the dependent variable is influenced by the neighboring  dependent variable, and the second is a spatial error in which there is a spatial correlation in the  residual.  Writing (1) in matrix form, and adding in these two spatial effects, we get    (3) where u is iid normal with mean 0 and variance 2. Both  and  are constrained to be between ‐1 and  1, and unless one is able to argue the reasonableness of a negative value for these parameters, it is  better to constrain them to be between 0 and 1. W is the weights matrix, a generally non‐symmetric  square matrix with the rows summing to 1 and the non‐zero elements measuring the degree of  influence one unit has on another. We could use a different W for each type of spatial relationship, but  it is much more common to use the same one for both. It is typical for there to only be a small number  of non‐zero elements in each row. If, for example, we used a weights matrix with the four nearest  neighbors, each row would only have four non‐zero numbers. A common functional form of W is inverse  distanced weighted, which we use. W is typically programmed as a sparse matrix, because as a regular  (full) matrix, the memory requirements can be huge.  For the analysis here with Brazil data, we have just over 100,000 disaggregated units. With an even  100,000, a full matrix would have 10 billion elements. If it takes 4 bytes per element, that would mean  we would need 40 gigabytes of RAM to hold W. With a sparse matrix, it would require much less.  Depending on how the data are coded, it might take only around 12 megabytes of memory.  If we did not need to aggregate the data to the level for which we have agricultural GDP computed –  which for Brazil is at the level of the municipio, our task would be easier, because we could focus on an  easier way to estimate version of equation (3), given by    However, because we need to sum a function of the explanatory variables and the variance matrix, we  unfortunately need to compute       (4) 7    Here, T is the transpose operator. The inverse of a very efficient sparse matrix might not be sparse  anymore, and computing the inverse not only can take a large amount of RAM, but also can take a long  time to compute.9   One solution to this dilemma is to take advantage of the following identity   ⋯ where the right‐hand side is an infinite series. The matrix W2 is the weights matrix for the neighbors of  neighbors, W3 is for neighbors of the neighbors of neighbors, etc. We can see from the identity that  ultimately each pixel is related to many if not all the others. With  being less than 1 and with each row  of W summing to 1 and therefore each element of W being less than or equal to 1, each successive term  on the right‐hand side is smaller in magnitude than the one before it. We thus can take an  approximation for the inverse by cutting off this series at a pre‐determined number of terms. For this  analysis, we chose to have the last term be W6.  To compute the variance matrix for the aggregated units from the variance matrix at the pixel level, we  simply sum all of the variance and covariance terms in each submatrix associated with each aggregated  unit i. This is based on the well‐known formula  2 ,   and another, less‐known formula  , , , , ,   For  , we do not need to compute the inverse. We simply need to know the product of the  inverse and X, which can be found relatively quickly using Gaussian elimination. Therefore, we do not  need to use the approximation for the inverse, but rather can compute   exactly, once we  specify the value for .  As we did in the simpler case without any spatial relationships to control for, our next step is to  aggregate. We can aggregate the variance matrix in (4) by summing over rows and columns for each  aggregate unit. For the explanatory variables, instead of aggregating X, we aggregate  .  In the aggregate, we have,  ̃, where ̃ is distributed normally with mean of 0 and variance  given by  .   We can then write the log likelihood function (Hauke and Kossowski, 2011) as  2 | |   where N is the number of aggregated units.   There are a number of ways to compute a determinant, but one of the most reliable ways to solve it is  through computing the eigenvalues,10 since a determinant is equal to the product of its eigenvalues.                                                               9  We tried to estimate this on a computer with 64 GB of RAM, taking the inverse of both L and U from the LU  decomposition in a 64‐bit version of R. The computation of each inverse took around 4 hours.  10  We tried using the “det” command in R, and it kept giving the value of infinity as the answer. Therefore, we  computed the eigenvalues in R with the “eigen” command, and then took the log and summed the results.  8    Because we are taking the log of the determinant, the log is equal to the sum of the log of each of the  eigenvalues. If we let the eigenvalues be given by I, the log likelihood function can then be written as  2 ∑   (5) With the aggregated variance matrix for Brazil having 5,434 rows and columns, the calculation can take  several minutes. One of the ways to speed this up is to compute a grid of eigenvalues, say for equal to  0 to 0.9 with a step of 0.1, and the same for . Then we could use bilinear interpolation on this grid of  eigenvalues to approximate the log of the Jacobian for any pair of  and . In our program written in R,  we create this as an option, but generally computed the exact Jacobian instead of the approximate  Jacobian.  Once the parameters are estimated, it is then possible to compute the predicted value at the aggregated  level of the dependent variable, and subtracting that from the actual value of the dependent variable  leaves the residual. We can use this residual to make the prediction at the disaggregated level perfectly  sum to the actual value of the dependent variable at the aggregated level. The residual value can be  treated as a “fixed effect” for the aggregated unit: a characteristic of the unit that was not measured by  variables of the regression but is by definition a characteristic of the aggregated unit.  The fixed effect at the aggregated level needs to be distributed to the disaggregated units. The easiest  way to do that is simply divide the fixed effect by the number of disaggregated units in that aggregated  unit, and this now becomes the fixed effect at the disaggregated level. However, this was to some  extent an ad hoc method, and other methods of rescaling could be developed.  To solve the likelihood function relatively quickly, we would maximize the concentrated log likelihood  function, concentrated for  and . Noting that the Jacobian and the aggregated variance matrix V  (along with V‐1) are functions of  and , but not for  or 2. We can solve for  and 2 as functions of J  and V‐1 and therefore ultimately of  and . They are    ⁄  Once the optimal values of  and  are found, in principle the values of all the parameter estimates can  be inserted into the full likelihood to find the proper standard errors of not only of  and , but also   and 2. In practice, this requires inverting a large numerically computed hessian, which may fail. An  alternative, used in this paper, is to loop through estimates of the log likelihood function computed for a  regression which excludes only one variable. This allows for use of the likelihood ratio test against the  full regression, which gives a probability for the parameter being non‐zero. This probability can be used  to create a pseudo‐z‐statistic, which can be used with the parameter estimate to give a pseudo‐standard  error. In practice, the standard errors computed this way were similar to the standard errors found  through the concentrated likelihood function that assumed that the spatial parameters were known  rather than estimated.  9    Data  The agricultural GDP data at the municipio level come from official data (IBGE 2016a), as does the  shapefile used to spatially locate the data (IBGE 2010b). We use the agricultural GDP data for 2010,  which are shown in Figure 1.   Figure 1. Brazilian agricultural GDP at the municipio‐level in 2010 (thousands of reais per  square kilometer)     We also have available several gridded databases that either show production of agricultural  commodities or show proxies of commodities that we can treat like production measures. the crop  maps from Spatial Production Allocation Model (SPAM) is one that shows crop harvested area and yields  (You et al. 2014). SPAM has production of multiple crops (in metric tons), gridded with pixels that are 5  arc‐minutes on edge. We also have similar data sets that cover livestock, fish, and forestry as measures  that summarize all the agricultural production in each pixel (though the fish and one of the forestry  measures are best seen as proxies of production).11  Another issue that is important to keep in mind is that the SPAM data set and the Gridded Livestock  data set present outputs of statistical models rather than real data. Ultimately, any statistical deficiency                                                               11  More specifically, we use measures derived from the following data sets: livestock data in Gridded Livestock of  the World (Gilbert et al. 2018), forestry measures of non‐wood from Siikämaki 2015 and wood from MODIS Land  Cover products (Friedl et al. 2010) and fishery measures from FISHSTAT and European Space Agency Climate  Change Initiative (2016).  10    in the models producing the data will contaminate our analysis which is based on the data. This might  lead to large errors from the regression if, for example, a crop in SPAM is in large quantity in a given  pixel, but in reality it is not in the pixel at all, or only in small quantity.  Results and Discussion  Spatial disaggregation  The regression results from using crop production from SPAM ‐‐ along with livestock, forest, and  fisheries data – are found in Table 1.   For crops, the units are in metric tons, and the dependent variable is in thousands of reais. So, if we  were to try to interpret the parameters as prices, the value for the parameter estimate for wheat is  0.456, or R$456 per metric ton (MT), which at the end of 2010 would be equal to USD274 per MT. This is  within the price range it could have sold for. Maize price from the regression, on the other hand, is low  at R$101 per metric ton. Given the explanatory variables are all modeled values, it is not entirely  surprising that the parameters do not reflect prices as well as we might hope. They do, however, allow  us to use these parameters to project agricultural GDP to each pixel, which was the goal.  Table 1 shows that most of the agricultural outputs that are important for the Brazilian agricultural  sector have statistically significant parameters. In addition to key cereals, cassava, beans, soybeans, and  a number of cash crops can be found, along with chickens, cattle, and pigs, and non‐timber forest  products.  Table 1. Regression using pixel‐level production data  Std  Variable Param  error  t‐stat  Prob  Wheat  0.456  0.083  5.473  0.000  Rice  0.061  0.028  2.167  0.030  Maize  0.101  0.025  3.982  0.000  Barley  ‐0.296  0.612  0.483  0.629  Sorghum  0.321  0.311  1.035  0.301  Other cereals  0.572  0.552  1.035  0.301  Potatoes  0.167  0.074  2.263  0.024  Sweet potatoes  2.786  0.526  5.298  0.000  Yams  ‐1.037  1.619  0.640  0.522  Cassava  0.081  0.029  2.761  0.006  Beans  0.504  0.212  2.381  0.017  Other pulses  11.812  9.957  1.186  0.236  Soybeans  0.043  0.019  2.204  0.028  Groundnuts  1.128  0.664  1.699  0.089  Coconuts  0.446  0.077  5.808  0.000  Oil palm  0.079  0.098  0.802  0.422  Sunflower  ‐1.708  3.836  0.445  0.656  11    Std  Variable Param  error  t‐stat  Prob  Rapeseed  71.528  197  0.363  0.717  Sesame  ‐726.78  1,899  0.383  0.702  Other oils  ‐0.214  1.348  0.159  0.874  Sugar cane  0.009  < 0.001  > 10  0.000  Cotton  0.384  0.133  2.882  0.004  Other fibers  0.120  0.391  0.308  0.758  Coffee, arabica  0.367  1.111  0.330  0.741  Coffee, robusta  1.547  1.703  0.908  0.364  Cocoa  1.072  0.768  1.395  0.163  Tea  32.869  20.193  1.628  0.104  Tobacco  1.369  0.233  5.883  0.000  Bananas  0.379  0.052  7.244  0.000  Fruit, tropical  0.050  0.011  4.681  0.000  Fruit, temperate  0.146  0.047  3.089  0.002  Vegetables  0.172  0.048  3.623  0.000  All others  0.990  0.539  1.838  0.066  Chickens  0.003  0.000  7.299  0.000  Cattle  0.045  0.008  5.591  0.000  Goats  ‐0.007  0.064  0.103  0.918  Pigs  0.069  0.015  4.537  0.000  Sheep  0.087  0.050  1.735  0.083  Inland water  0.107  0.092  1.159  0.246  Forest loss  0.211  0.117  1.801  0.072  Non‐timber products  0.000  0.000  3.181  0.001  Constant  ‐7.538  27.877  0.270  0.787    Log likelihood  ‐7199.47    Lam  0.6223  0.169  3.677  0.000  Rho  ‐0.5343  0.419  1.274  0.203  Source: Authors.  Notes: We were unable to compute the inverse of the Hessian matrix. So we iterated through the parameters,  assigning the value of 0 to one at a time, and recomputing the value of the maximized likelihood function. We then  applied the likelihood ratio test to get a probability. We treated the probability as if it were from a normal  distribution to get the t‐statistic (z‐statistic). Then we used the t‐statistic with the parameter estimate to derive a  pseudo‐standard error. The pseudo‐standard error was very close to the uncorrected standard errors that were  computed as part of the concentrated likelihood function.  Comparing to cross‐entropy estimate  We compare the predictive accuracy of the three methods by making predictions at each pixel and then  aggregating up to the municipio level. We have the true agricultural GDP measures at the municipio  12    level, even though we only used in the regressions the values at the microregion level. This allows us to  test how well each method works, and it gives us a way to compare the methods.   We see from Table 2 that the spatial disaggregation method performed best in terms of both mean  absolute deviation and root means square error. The cross‐entropy method performed second best, and  not that much less successively, with a mean absolute deviation that is only 15 percent higher than that  for the spatial disaggregation method, and a root mean square error that is only 10 percent higher. Both  methods performed better than that of the traditional method of using rural population to map  agricultural GDP. In fact, the traditional method had a mean absolute deviation that was 4 times higher  than that of the spatial disaggregation method. The root mean square error was only 50 percent higher,  meaning that it was much less likely to generate large outliers.  Figure 2 shows the complete distribution of actual and predicted values for municipios for the three  methods. The spatial disaggregation method had several municipios with predicted values of 0 but  actual values greater than 0. In fact, there are several large outliers using this method in which the  predicted value was much higher than the actual. This may suggest that the somewhat ad hoc method  for adjusting the fixed effects was less than optimal and an alternative method should be found.  Generally, though, the spatial disaggregation method cluster near the diagonal line is narrower than that  of the cross‐entropy spatial allocation method, which in turn is narrower than the one for the rural  population density method.  Table 2. Deviations of pixel agricultural GDP summed to municipio from actual value    Mean absolute  Root mean      deviation    square error  Rural population density  28,744  25,397  Cross‐entropy spatial allocation  8,249  18,347  Spatial disaggregation from regression on agricultural  7,214  16,673  production  Source: Authors.  Figure 3 shows the pixel‐level results produced by the three methods. The difference between the white  and the orange should mostly be overlooked. It shows different assumptions that fixed certain  municipios and their underlying pixels to be valued at 0 (the white areas). Ignoring a somewhat artificial  contrast between the white and orange, the cross‐entropy method and the spatial disaggregation  method appear to be very similar in most locations. The results using the population density are less  smooth, looking more speckled, based simply on the nature of the underlying rural population data set.  Conclusions  When there are enough data to allow high enough degrees of freedom to use a regression approach for  all the key crops and contributors to agricultural GDP, the spatial disaggregation approach is superior. In  principle, we could have used many more parameters – though the hedonic price interpretation would  13    potentially no longer be valid. We could, for example, include population density as an explanatory  variable, in addition to all the other variables. We could have used market distance, elevation, distance  to roads, soil type, climate variables, etc.   Figure 2. Predictions and actual values for the three methodologies  Rural population density (correl = 0.8105)  Cross‐entropy spatial allocation (correl = 0.9066)      Spatial disaggregation (correl = 0.9211)    Source: Authors.  Note: For the sake of graphing using a logarithmic scale, values of 0 for agricultural GDP were converted to values  of 1. The correlation of actual and predicted (“correl”) is given in the title of each graph.  Furthermore, the spatial disaggregation approach can be used for problems unrelated to GDP. It could  be used for any kind of spatial allocation problem for which gridded data can help explain at a pixel level  the value associated with any aggregated total or average.  14    Figure 3. Comparison of projections of agricultural GDP          Notes: Top left, spatial disaggregation; top right, cross entropy; bottom left, rural population density.  Source: Authors.  However, the issue of degrees of freedom is an important limiting factor for that method as it applies to  GDP. There are not many countries that have enough subnational units with computed agricultural GDP.  Many developing countries only have national‐level agricultural GDP available. That is why the cross‐ entropy spatial allocation method is of such importance. It can operate with only national‐level  agricultural GDP – just as the traditional method of using rural population density can operate with only  national‐level agricultural GDP. But the cross‐entropy approach on readily available spatially distributed  15    crop, livestock, forest, and fish allocation far outperforms the method of using rural population to  distribute the agricultural GDP, at least in the case of Brazil.  References  S. Angel, P. Lamson‐Hall, B. Guerra, Y. Liu, N. Galarza, and A. M. Blei (2018). Our Not‐So‐Urban World,  Working Paper No. 42 The Marron Institute of Urban Management, New York University.  Arvis, J. F., & Shepherd, B. (2013). The Poisson quasi‐maximum likelihood estimator: a solution to the  ‘adding up’ problem in gravity models. Applied Economics Letters, 20(6), 515‐519.  Arvis, J.F; Blankespoor, B. Downscaling Economic Variables: An Entropy Estimation. Mimeo.  Blanchard, P.B.; Blankespoor, B. Rivera‐Fuentes, J.; Gunasekera, R.; Ishizawa O. and Jimenez L.F. Spatial  Disaggregation of Gross Domestic Product for applications in Disaster Risk Management. Mimeo.  Bosker, M.; Park, J.; Roberts, M. (2018). Definition Matters : Metropolitan Areas and Agglomeration  Economies in a Large Developing Country (English). Policy Research working paper; no. WPS 8641.  Washington, D.C. : World Bank Group.  Center for International Earth Science Information Network ‐ CIESIN ‐ Columbia University. (2017).  Gridded Population of the World, Version 4 (GPWv4): Population Count Adjusted to Match 2015  Revision of UN WPP Country Totals, Revision 10. Palisades, NY: NASA Socioeconomic Data and  Applications Center (SEDAC). https://doi.org/10.7927/H4JQ0XZW.  Chen, X., & Nordhaus, W. D. (2011). Using luminosity data as a proxy for economic statistics.  Proceedings of the National Academy of Sciences, 108(21), 8589‐8594.  Dobson, J. E., Bright, E. A., Coleman, P. R., Durfee, R. C., & Worley, B. A. (2000). LandScan: a global  population database for estimating populations at risk. Photogrammetric engineering and remote  sensing, 66(7), 849‐857.  Dijkstra, L., & Poelman, H. (2014). A harmonised definition of cities and rural areas: the new degree of  urbanisation. Working Paper, 1, 2014.  Doll, C. H., Muller, J. P., & Elvidge, C. D. (2000). Night‐time imagery as a tool for global mapping of  socioeconomic parameters and greenhouse gas emissions. AMBIO: a Journal of the Human  Environment, 29(3), 157‐163.  Doll, C. N., Muller, J. P., & Morley, J. G. (2006). Mapping regional economic activity from night‐time light  satellite imagery. Ecological Economics, 57(1), 75‐92.  Ebener, S., Murray, C., Tandon, A., & Elvidge, C. C. (2005). From wealth to health: modelling the  distribution of income per capita at the sub‐national level using night‐time light imagery. International  Journal of Health Geographics, 4(1), 5.  Elvidge, C. D., Erwin, E. H., Baugh, K. E., Ziskin, D., Tuttle, B. T., Ghosh, T., & Sutton, P. C. (2009, May).  Overview of DMSP nightime lights and future possibilities. In 2009 Joint Urban Remote Sensing  Event (pp. 1‐5). IEEE proceedings of the 7th international urban remote sensing conference, Shanghai,  China.  16    Friedl, M. A., Sulla‐Menashe, D., Tan, B., Schneider, A., Ramankutty, N., Sibley, A., & Huang, X. (2010).  MODIS Collection 5 global land cover: Algorithm refinements and characterization of new datasets.  Remote sensing of Environment, 114(1), 168‐182.  IIASA/FAO. (2012). Global Agro‐ecological Zones (GAEZ v3.0). IIASA, Laxenburg, Austria and FAO,  Rome, Italy.  Gilbert, M., Nicolas, G., Cinardi, G., Van Boeckel, T. P., Vanwambeke, S. O., Wint, G. W., & Robinson, T. P.  (2018). Global distribution data for cattle, buffaloes, horses, sheep, goats, pigs, chickens and ducks in  2010. Scientific data, 5, 180227.  Ghosh, T., Anderson, S., Powell, R., Sutton, P., & Elvidge, C. (2009). Estimation of Mexico’s informal  economy and remittances using nighttime imagery. Remote Sensing, 1(3), 418‐444.  Ghosh, T., L Powell, R., D Elvidge, C., E Baugh, K., C Sutton, P., & Anderson, S. (2010). Shedding light on  the global distribution of economic activity. The Open Geography Journal, 3(1).  Gunasekera, R., Ishizawa, O., Aubrecht, C., Blankespoor, B., Murray, S., Pomonis, A., & Daniell, J. (2015).  Developing an adaptive global exposure model to support the generation of country disaster risk  profiles. Earth‐Science Reviews, 150, 594‐608.  Hauke, J., & Kossowski, T. (2011). Comparison of values of Pearson's and Spearman's correlation  coefficients on the same sets of data. Quaestiones Geographicae, 30(2), 87‐93.  Henderson, J. V., Storeygard, A., & Weil, D. N. (2012). Measuring economic growth from outer  space. American economic review, 102(2), 994‐1028.   Instituto Brasileiro de Geografia e Estatística (IBGE), Coordenação de Contas Nacionais. (2016a). Produto  interno bruto dos municípios : 2010‐2014. Rio de Janeiro : IBGE, 2016.  https://sidra.ibge.gov.br/tabela/5938#notas‐tabela.  Instituto Brasileiro de Geografia e Estatística (IBGE). (2016b). Malha Municipal 2010. Rio de Janeiro:  IBGE. https://mapas.ibge.gov.br/bases‐e‐referenciais/bases‐cartograficas/malhas‐digitais.html  Kummu, M., Taka, M., & Guillaume, J. H. (2018). Gridded global datasets for gross domestic product and  Human Development Index over 1990–2015. Scientific data, 5, 180004.2018/02/06/online.  https://www.nature.com/articles/sdata20184.  Nordhaus, W. D. (2006). Geography and macroeconomics: New data and new findings. Proceedings of  the National Academy of Sciences, 103(10), 3510‐3517.  Pesaresi, Martino; Freire, Sergio (2016): GHS Settlement grid following the REGIO model 2014 in  application to GHSL Landsat and CIESIN GPW v4‐multitemporal (1975‐1990‐2000‐2015). European  Commission, Joint Research Centre (JRC) [Dataset] PID: http://data.europa.eu/89h/jrc‐ghsl‐ ghs_smod_pop_globe_r2016a  Rafa, Mickey, Jonathan D. Moyer, Xuantong Wang, and Paul Sutton. (2017). Estimating District GDP in  Uganda, report for Frederick S. Pardee Center for International Futures, Josef Korbel School of  International Studies, University of Denver, November.  17    Roberts, M.; Blankespoor, B.; Deuskar, C.; Stewart, B. (2017). Urbanization and Development : Is Latin  America and the Caribbean Different from the Rest of the World?. Policy Research Working Paper;No.  8019. World Bank  Siikamäki, F. J. Santiago‐Ávila, and P. Vail. (2015). GLOBAL ASSESSMENT OF NONWOOD FOREST  ECOSYSTEM SERVICES, Resources For the Future, Washington, D.C.  Sutton, P. C., Elvidge, C. D., & Ghosh, T. (2007). Estimation of gross domestic product at sub‐national  scales using nighttime satellite imagery. International Journal of Ecological Economics &  Statistics, 8(S07), 5‐21.  Sutton, P. C., & Costanza, R. (2002). Global estimates of market and non‐market values derived from  nighttime satellite imagery, land cover, and ecosystem service valuation. Ecological Economics, 41(3),  509‐527.  Tobler, W., Deichmann, U., Gottsegen, J., & Maloy, K. (1997). World population in a grid of spherical  quadrilaterals. International Journal of Population Geography, 3(3), 203‐225.  UNISDR (2011) Global Assessment Report on Disaster Risk Reduction. Geneva, Switzerland: United  Nations International Strategy for Disaster Reduction.  Yetman, G., S.R. Gaffin, and X. Xing. (2004). Global 15 x 15 Minute Grids of the Downscaled GDP Based  on the SRES B2 Scenario, 1990 and 2025. Palisades, NY: NASA Socioeconomic Data and Applications  Center (SEDAC). http://dx.doi.org/10.7927/H4NC5Z4X.  You, L., Wood, S., Wood‐Sichra, U., & Wu, W. (2014). Generating global crop distribution maps: From  census to grid. Agricultural Systems, 127, 53‐60.  You, L., & Wood, S. (2006). An entropy approach to spatial disaggregation of agricultural  production. Agricultural Systems, 90(1‐3), 329‐347.7  18