Smith ScholarWorks Smith ScholarWorks More Than a Safety Net: Ethiopia's Flagship Public Works More Than a Safety Net: Ethiopia's Flagship Public Works Program Increases Tree Cover Program Increases Tree Cover

More than one billion people worldwide receive cash or in-kind transfers from social protection programs. In low-income countries, these transfers are often conditioned on participation in labor-intensive public works to rehabilitate local infrastructure or natural resources. Despite their popularity, the environmental impacts of public works programs remain largely undocumented. We quantify the impact on tree cover of Ethiopia ’ s Productive Safety Net Program (PSNP), one of the world ’ s largest and longest-running public works programs, using satellite-based data of tree cover combined with difference-in-differences and inverse probability treatment weighting methodologies. We find that the PSNP increased tree cover by 3.8% between 2005 and 2019, with larger increases in less densely populated areas and on steep-sloped terrain. As increasing tree cover is considered an important strategy to mitigate global warming, our results suggest a win – win potential for social safety net programs with an environmental component.


Introduction
Reducing poverty while addressing climate change and restoring terrestrial ecosystems are critical challenges that lie at the core of the United Nations Sustainable Development Goals (SDGs) (United Nations, 2015). Despite sub-Saharan Africa's impressive economic growth over the past two decades , this region is projected to host the greatest number of poor and undernourished people in the world by 2030 Yonzan et al., 2020). Moreover, rapid population growth combined with climate change are likely to hasten environmental degradation in the region (Bradshaw and Di Minin, 2019;Olsson et al., 2019). To address these challenges, governments and international organizations are turning to social safety net programs that provide cash or in-kind transfers to the poorest and most vulnerable segments of society (Kuriakose et al., 2013;World Bank, 2018). It is estimated that more than one billion people worldwide receive assistance from such programs (Alderman et al., 2017). Since 2000, the number of safety net programs in sub-Saharan Africa has doubled (Hickey et al., 2018) and today, all 46 sub-Saharan countries implement at least one program .
While safety net programs have generally been found to improve food security and increase asset accumulation Hidrobo et al., 2018), the evidence on their environmental impacts remains mixed. Studies linking safety net programs to environmental outcomes have been largely limited to cash transfer programs conditioned on beneficiary households meeting health or education related objectives. For example, Mexico's Oportunidades program (which provided conditional cash transfers for school attendance, health clinic visits, and nutritional support) increased deforestation, with larger impacts found in poorer and more remote communities (Alix-Garcia et al., 2013). In contrast, Indonesia's Keluarga Harapan program (which provided conditional cash transfers if households accessed specific health and educational services) reduced expected deforestation (Ferraro and Simorangkir, 2020). Brazil's Zero Hunger social protection program, which includes a conditional cash transfer component (conditioned against on child school attendance and family health checks) had mixed impacts on natural vegetation cover, which varied by biome (Dyngeland et al., 2020).
Many safety net programs include public works components, which hold particular promise in delivering on both social and environmental objectives. In these programs, beneficiary households receive cash or inkind transfers conditioned on labor-intensive works that aim to build or restore community assets, such as roads, schools, or degraded natural resources including communal lands and forests Subbarao et al., 2012). Public works programs are popular in South-Asia and sub-Saharan Africa, with the largest programs found in India and Ethiopia covering millions of beneficiaries (World Bank, 2018). Globally, it is estimated that more than $10 billion USD are spent annually on public works programs that provide work to almost 70 million people (McCord and Paul, 2019). Despite their popularity, the extent to which public works programs generate public goods that promote development and environmental sustainability remains poorly understood (Beierl and Grimm, 2019;Gehrke and Hartwig, 2018;Ravallion, 2019).
We examine the effects of Ethiopia's Productive Safety Net Program (PSNP) on tree cover between 2005 and 2019 and estimate how potential carbon sequestration benefits may offset the administrative costs of the program and reduce CO 2 emissions. With eight million beneficiaries (World Bank, 2020), the PSNP is the largest public works program in the world outside of India (World Bank, 2018). Being implemented by the government of Ethiopia, its design and success at achieving social protection gains has made it a model for other social protection programs on the African continent (Monchuk, 2013). The purpose of the PSNP is to relieve poverty and food insecurity through cash or in-kind transfers in exchange for labor on public works designed to build sustainable community assets that increase communities' resilience to shocks (MoARD, 2006;Wiseman et al., 2010). The public works projects are implemented exclusively on publicly owned lands and are identified and designed by the communities themselves with technical support from higher administrative levels (MoARD, 2006;Wiseman et al., 2010). The focus of these work projects has largely been soil and water conservation activities like terracing, embankments, gully check dams, water-infiltration trenches, and especially reforestation (MoARD, 2006;Wiseman et al., 2010).
These environmental activities of the PSNP can potentially help to alleviate the negative impacts of climate change, contribute to climate change mitigation, and restore terrestrial ecosystems. Among these, the PSNP's potential to increase tree cover is of particular interest in this paper. Deforestation and land degradation are major environmental problems in Ethiopia (Lemenih and Kassa, 2014) with the former being a substantial source of carbon emissions worldwide (IPCC, 2019). In the last three decades, Ethiopia is estimated to have lost 33,400 km 2 of forest cover (falling from 204,100 km 2 in 1990 to 170,700 km 2 in 2020) (World Bank, 2021). Globally, the urgency to maintain and increase tree cover has launched several initiatives including the Bonn Challenge (International Union for Conservation of Nature, 2021), the New York Declaration on Forests (Climate and Land Use Alliance, 2021), and the African Forest Landscape Restoration Initiative, AFR100 (African Union Development Agency, 2021), to which Ethiopia is a major contributor. Forests and trees play a key role in the regulation of water, energy, and carbon cycles and have climatic and environmental benefits that support adaptation and mitigation strategies for climate change (Ellison et al., 2017). Trees reduce erosion, stabilize water supply, increase soil fertility, and can exert a cooling effect and promote rainfall-making communities more resilient against adverse impacts of climate change (Ellison et al., 2017;Zuazo and Pleguezuelo, 2009). In addition, the relatively high rate of carbon sequestration of trees makes increasing tree cover an important global warming mitigation strategy, among others (Griscom et al., 2017;Masson-Delmotte et al., 2018;Vincent et al., 2021).
Independent evaluations show that the PSNP has been successful in improving household food security, resilience, and asset levels (Berhane et al., 2014;Knippenberg and Hoddinott, 2017). These studies used a large panel dataset representative of all areas of PSNP implementation. The few impact evaluations focused on the PSNP's environmental outcomes, however, have focused on a much smaller geographic area using households surveys. These studies have investigated participants' investments in sustainable land management practices like soil erosion and soil fertility practices in two districts (Adimassu and Kessler, 2015) and household-level livestock and tree holdings in six sub-districts (Andersson et al., 2011). While household soil management and tree planting strategies can have positive environmental effects, these studies did not have sufficiently broad data to quantify the size of a programwide benefit from these activities. In addition, PSNP's potential for climate change adaptation and mitigation is largely unknown, although it is increasingly being recognized (Subbarao et al., 2012). Relevant empirical studies include Conway and Schipper (2011)'s analysis of mainstreaming climate risk adaptation actions into development initiatives using a case study on drought risk financing mechanisms within the PSNP, and Woolf et al. (2018)'s estimation of PSNP's Global Greenhouse Gas (GHG) reduction based on 24 site surveys on sustainable land, soil, and water practices using an IPCC based modeling approach.
This study advances our understanding of the PSNP's environmental outcomes and its potential for climate change mitigation by providing a robust impact evaluation of the PSNP on tree cover. In the context of existing studies of the PSNP, our study is novel in four ways. First, we use satellite imagery on tree cover and other spatial variables allowing us to cover a larger area consistently throughout the study period. Second, we use an econometric method to assess the impact of the PSNP by applying difference-in-differences and inverse probability treatment weighting methods to construct a credible counterfactual (i.e., what the tree cover would have been had the PSNP area not participated in the program). The control areas were identified via statistical matching and can be thought of as locations that appear similar to those participating in the PSNP based on agro-ecological and socio-economic characteristics, but were not treated with the program. Third, we conduct a series of checks to explore the robustness of our results. These checks include altering the way our outcome variable is defined, exploring the sensitivity of changing the way we control for common shocks and district characteristics, recalculating our confidence intervals using a method that is robust to spatial autocorrelation (Conley, 1999) and checking whether data quality issues are driving our findings. Finally, we estimate how the social benefits from the estimated tree growth could offset the administrative costs of the program.

Ethiopia's Productive Safety Net Program (PSNP)
Ethiopia is the second most populous country in Africa, with a population of over 110 million, projected to increase 2.5% annually (World Bank, 2021). Rainfed agriculture is a major component of the national economy providing livelihood to approximately 80% of the population. Ethiopia's history is characterized by catastrophic droughts that triggered the large-scale famines in the 1970s and 1980s. Meanwhile, the 1990s and early 2000s were characterized by localized food shortages that were typically addressed by ad hoc requests for humanitarian food aid (De Waal, 2017). Despite substantial economic growth coupled with major improvements in various domains of health and development over the last two decades, the country remains vulnerable to droughts and flooding with climate change expected to further intensify these adverse weather events (Alemu and Mengistu, 2019;Conway and Schipper, 2011;Federal Democratic Republic of Ethiopia, 2021;Funk et al., 2008).
Launched in February 2005, the PSNP was designed as a multiyear food security program to provide a more sustainable response mechanism than recurring ad hoc humanitarian appeals (Wiseman et al., 2010). The households benefiting from the PSNP receive food or cash payments in return for labor-intensive public works carried out over a six-month period outside of the main agricultural season while a small share of households with limited labor capacity (e.g., pregnant and lactating women, elderly) receive unconditional transfers. The PSNP is largely externally funded (World Bank, 2018), but the program is led and implemented by the government of Ethiopia.
At the onset of the program, there were 192 districts (woredas, 3rdlevel administrative division in the country) with 4.8 million beneficiaries in the four highland regions (Amhara; Oromia; Southern Nations, Nationalities and Peoples' Region; and Tigray), as well as smaller and predominantly urban regions in the east (Dire Dawa and Harar) (World Bank, 2020). Since the launch of the PSNP, caseloads in the original PSNP districts in the highland regions have increased and the program has expanded to Ethiopia's lowland regions (Afar and Somali). By 2019, the PSNP operated in more than 300 districts, providing support for approximately eight million people (World Bank, 2020). So far, none of the districts selected into the PSNP have exited the program (World Bank, 2020).  Table A1 for data sources.

4
The PSNP combines geographic and community level targeting. Districts were initially selected for the program based on the frequency they had requested and received emergency food assistance prior to the launching of the program in 2005 (MoARD, 2006;World Bank, 2020). Communities themselves then select the most food-insecure households as PSNP beneficiaries. Evaluations based on household data collected from the PSNP localities show that the program is relatively well targeted at the community level (Coll-Black et al., 2011). However, a recent assessment of the geographic targeting suggests that many poor and food insecure districts are not included into the PSNP (World Bank, 2020).
Our analysis measures the impact of the PSNP on tree cover in the participating districts of four highland regions (Amhara; Oromia; Southern Nations, Nationalities and Peoples' Region; and Tigray) (Fig. 1A). The reasons for this geographic restriction were threefold. First, the PSNP has operated the longest time in these highland regions, permitting a longer time window to observe impacts on tree cover. We also note that the highland regions did not have a staggered roll-out of the program. Second, while the program has expanded to other regions since its inception, the focus on the highland regions has remained. In 2019, more than 70% of all PSNP beneficiaries originated from the four highland regions. Third, compared to the two lowland regions, the PSNP has been relatively better-implemented in the highland regions (Lind et al., 2022;Sabates-Wheeler et al., 2013).
Collection 6 is the most accurate MODIS fractional cover product to date and has been improved from previous collections with updated input data (DiMiceli et al., 2021). The data are distributed as a global tiled grid in Sinusoidal projection at a spatial resolution of 250 m. We mosaicked the VCF-TC tiles to cover Ethiopia throughout the study period and integrated them with the data described below in Sinusoidal projection using ArcGIS 10.7 (ESRI, 2019) and Terrset (Clark Labs, 2019).
Previous research using VCF-TC has noted that year-to-year variation in tree cover estimates appears to be higher than expected and may be driven by the quality of the underlying remote sensing data and precipitation, among other factors (Gao et al., 2018;Zomer et al., 2016), which discourages its use for inter-annual comparisons. To mitigate this issue we followed Zomer et al. (2016) and calculated three-year averages of tree cover for our analysis (except for the first period which is based on a two-year average): 2000-2001; 2002-2004; 2005-2007; 2008-2010; 2011-2013; 2014-2016; 2017-2019. Each VCF-TC pixel was coded to indicate its participation (or not) in the PNSP (Fig. 1A) and matched to its corresponding district and region using the district boundaries from Ethiopia's Central Statistical Agency in 2007 (unpublished data) and the PSNP district administrative records. In addition, we obtained the annual PSNP beneficiaries at district level by digitizing PSNP's annual planning documents drafted by the Ministry of Agriculture of Ethiopia. We used the year 2007 as a benchmark for the administrative units since it matches with the latest Ethiopian census year and corresponding administrative boundaries. Increases in the number of PSNP-eligible highland districts from the census year onward were due to administrative divisions of the districts The PSNP increased tree cover, particularly in less densely populated areas and steep-sloped terrains. Tree cover also increased in forests and woody areas (see Table A3 for exact aggregation) and in croplands, based on land cover classifications defined at the onset of the program in 2005. Estimates measure % change in tree cover due to the PSNP calculated using pixel-level observations. All estimates are based on a difference-in-differences method combined with an inverse probability treatment weighting. The unit of observation is a pixel observed periodically. The 95% confidence intervals are computed from standard errors clustered at the district level. A: Impact estimates for all pixels in the study region (N = 45,229,114) and for rural areas defined as population density <300 people/km 2 (N = 42,977,984) and <150 people/km 2 (N = 35,702,079). B: Impact estimates by terrain slope quintiles: 1st quintile (Q1): 0.0 to 2.9 degrees (N = 10,403,050); 2nd quintile (Q2): 3.3 to 5.6 degrees (N = 7,954,821); 3rd quintile (Q3): 5.8 to 10.7 degrees (N = 8,781,612); 4th quintile (Q4): 10.7 to 19.1 degrees (N = 9,055,319); 5th quintile (Q5): 19.1 to 78.4 degrees (N = 9,034,312). C: Impact estimates for different land cover types at the onset of the program in 2005: Forests (N = 1,956,304); Croplands (N = 16,631,643); Grasslands (N = 18,452,994); Savannas (N = 8,178,569). (Wiseman et al., 2010), which we dealt with by merging the child districts back to their parent district as of 2007 along with the number of PSNP beneficiaries. After removing pixels flagged as no data or over large bodies of water, the area of our study region was approximately 61.4 million ha (about 11.4 million of pixels), including 617 districts. Of these pixels, 49.5% were located in participating PSNP districts (247 districts).
In addition to VCF-TC, we used several datasets to generate spatial variables as controls and to explore impact heterogeneity. These included: population density per km 2 from the Gridded Population of the World, 2005 (GPW) (CIESIN: Center for International Earth Science Information Network, Columbia University, 2016) (Fig. 1C), elevation from the Shuttle Radar Topography Mission (SRTM), a global digital elevation model (DEM) of the world (USGS, 1996) from which we also derived the slope (Fig. 1D), and land cover type at the onset of the program in 2005 from MODIS (MCD12Q1) (Friedl and Sulla-Menashe, 2015). We used the International Geosphere-Biosphere Programme (IGBP) land cover classification scheme and aggregated the land over classes into eight categories for mapping purposes (Fig. 1E). The Climate Hazards Group InfraRed Precipitation with Station (CHIRPS) annual rainfall data (Funk et al., 2015) (Fig. 1F) from 1995 to 2019 was used to control for rainfall.
We also used the aboveground live woody biomass density dataset (AGBM) for the year 2000 (Zarin et al., 2016) distributed by the (Global Forest Watch GFW, 2021) to estimate the average AGBM (in megagrams biomass ha-1) corresponding to different tree cover percentages as the first step of the carbon sequestration and the benefit-cost analyses.
Lastly, the quality flag information provided with VCF was used to assess the sensitivity of our results to the uncertainty in tree cover estimates associated with data quality. The quality information of the input MODIS surface reflectance data used to predict the vegetation cover is provided as a separate quality band indicating if a pixel in any of the eight input data periods used to generate the annual product is flagged as poor quality due to clouds, high aerosol levels, cloud shadows, or having a view zenith angle higher than 45 • (Townshend et al., 2017). Estimates of vegetation cover with two or more flags in a year may be erroneous and should be used with caution (Townshend et al., 2017). Table A1 in the Appendix summarizes the characteristics and sources of the datasets used. For more details about VCF and the spatial methods see Section A.1 in the Appendix, DiMiceli et al. (2021); Townshend et al. (2017), and Hansen et al. (2003).

Impact of PSNP on tree cover
We evaluated the impact of the PSNP program on tree cover by applying a difference-in-differences method. Specifically, we estimated the difference in tree cover before and after the PSNP program began, and in participating PSNP districts versus non-PSNP participating districts (Fig. 1A). Implementing our difference-in-differences method using a regression approach, we estimated: where TC iwrt is the mean percent of tree cover in pixel i in district w in region r during the three-year period t. PSNP w is a binary variable that is defined at the district level; it equals one if the pixel belongs to a district that was selected into the program in 2005-2006 and equals zero otherwise. The variable POST t equals one if period t occurs after the 2005 launch of the PSNP (i.e., periods 2005-2007, 2008-2010, 2011-2013, 2014-2016, or 2017-2019) and equals zero if the period occurs before the PSNP launch (i.e., 2000-2001 or 2002-2004). We controlled for mean annual rainfall in pixel i in period t (X it ) as well as all period-and region-specific aggregate shocks through period-by-region fixed effects (α rt ). The term u iwrt represents the error term. The impact of the PSNP on the change in log of tree cover is measured by β; the coefficient on the interaction between PSNP w and POST t . We converted these coefficients to percentages using the following formula: (e β − 1) * 100. Finally, we clustered our standard errors at district level; i.e., at the level in which the treatment variable was defined (Abadie et al., 2017).
The key identifying condition of the difference-in-differences method in our application is that tree cover in the pixels within treatment (PSNP) and control (non-PSNP) districts was on a parallel trend before the program began in 2005. Ethiopian highlands are extremely diverse agro-ecologically, ranging from rugged high altitude plateaus in the north and central to arid and semi-arid terrains in the south (see Fig. 1F). The western highlands enjoy reliable and abundant rainfall, while the conditions in the east-where most PSNP districts are located-are generally drier with more erratic rainfall (Fig. 1F). Unsurprisingly then, the parallel trend hypothesis was rejected when we used all non-PSNP pixels in the highland regions as our control areas ( Table A4 in the Appendix). To address this, we first restricted the analysis to PSNP and non-PSNP pixels that had similar agro-ecological conditions before the program was launched in 2005. To do so, we used a propensity score matching algorithm (Rosenbaum and Rubin, 1983) to identify an area of common support (Caliendo and Kopeinig, 2008); a set of PSNP and non-PSNP pixels with sufficient overlap in predicted probability to be included into the program based on selected agro-ecological and socioeconomic characteristics (Figs. A1, A2, A3 in the Appendix). Fig. A2 in the Appendix shows the spatial distribution of the propensity scores and the area of common support. As matching covariates, we considered variables that were likely to capture this agro-ecological heterogeneity and thus correlate with selection into the program in 2005-2006: mean and standard deviation of annual rainfall in 1995-2004 (and their quadratics), population density in 2005, elevation and the slope of terrain (Table A5 in the Appendix). Finally, since the PSNP implementation and targeting benchmarks vary across administrative regions (Wiseman et al., 2010), we also included binary indicators for each region in our matching model. We defined the area of common support as pixels with the estimated propensity score within the [0.1; 0.9] interval (Crump et al., 2009) (Table A6).
We then used these pixel-level propensity scores (PS) to calculate inverse probability treatment weights (IPTW) (Abadie, 2005;Joffe et al., 2004): 1/PS for the treated (PSNP) pixels and 1/(1 − PS) for the untreated (non-PSNP) pixels. After restricting the pixels in our dataset to common support and applying IPTW on our regression model, the parallel trend assumption was satisfied; we cannot reject the null hypothesis that the tree cover in the PSNP and non-PSNP districts areas were on a similar trend before the PSNP was launched in 2005 (Table A4, Col.7 in the Appendix). The matching covariates were also in balance after restricting pixels to the common support and applying IPTW (Table A7 in the Appendix). Once we restricted the area to the common support, the final data used in the analysis had approximately 6.5 million pixels (53% from PSNP districts), coming from 513 districts (227 of which were PSNP districts). Section A.2 in the Appendix provides more information about our impact evaluation approach. Figs. A4,A5,A6 in Section A.3 in the Appendix show the distributions of key variables used in the analyses, after restricting to the area of common support.

Spatial variability of tree cover change
We explored heterogeneity in impact across socio-economic and environmental characteristics at pixel-level to better understand the PSNP impacts through space, particularly in the context of the program design and objectives (see Section A.4 in the Appendix). First, as the PSNP is a rural public works program, ideally we would have restricted the analysis to rural areas only. However, Ethiopia does not have an official definition for rural areas based on population density . Mindful of this ambiguity, we estimated the impacts for all pixels as well as for rural pixels based on two different population density thresholds. Following the recent recommendation made by international organizations (EU and UN-Habitat, 2020), we defined rural areas as pixels that fall below 300 people per km 2 . As an alternative definition for rural areas, we used 150 people per km 2 population density threshold based on previous work in Ethiopia (Schmidt and Kedir, 2009;. Second, the PSNP was specifically designed with the objective of rehabilitating sloped areas in mind (MoA, 2010). Sloped areas are also less suited for agriculture due to the increased risk of erosion and soil degradation, problems that exacerbate at higher slope inclination (Shaxson, 1999). We therefore hypothesized that the PSNP's impact on tree cover was likely to be larger in sloped terrain. To explore this, we sequentially restricted the analysis to quintiles based on the inclination of the slope. Third, we explored whether the impacts varied by the type of land cover at the onset of the program in 2005. We aggregated land cover types derived from MODIS into Forests and woody areas (7.6% of all pixels), Croplands (25.7%), Grasslands (43.6%), and Savannas (22.8%). Pixels categorized as urban, wetland, water, or barren (0.3% in total) were not considered (Tables A2 and A3 in the Appendix). Finally, we explored whether the impacts were larger in districts that had more PSNP beneficiaries relative to their total population compared to districts that had fewer. To do this analysis, we computed the average number of PSNP beneficiaries in each district over the study period and divided this number by its total population. Using this variable as our measure of beneficiary caseload intensity, we split the pixels originating from the PSNP districts into two groups using the median caseload intensity as the threshold. We then replaced our treatment variable with these two binary variables and reran the regression.

Spillover analysis
We estimated the percent tree cover change in non-PSNP districts adjacent to PSNP districts to assess if the PSNP had an spillover effect into neighboring non-PSNP districts that could have affected our results. First, we identified the pixels from 134 districts that did not benefit from the PSNP but shared an administrative border with a PSNP district through spatial analysis. We then appended the estimated model with an additional treatment variable capturing these adjacent non-PSNP districts (see Section A.5 in the Appendix for more details).

Carbon sequestration and benefit-cost analysis
Changes in tree cover were converted to sequestered carbon using the VCF-TC data and average aboveground live woody biomass (AGBM) from the Global Forest Watch GFW (2021) dataset (Zarin et al., 2016) in metric tons of biomass per ha for the year 2000. First we calculated the estimated average of AGBM corresponding to the percent tree cover at pixel level in 2000. Then we used our regression estimates, along with the baseline levels of tree cover in 2000 to estimate the predicted increase in tree cover (and hence AGBM), due to the PSNP and converted this biomass to carbon.
To estimate the benefit-cost ratio of the negative carbon emissions relative to the PSNP's cost we used the social cost of carbon (SCC) estimate from the Interagency Working Group on Social Cost of Greenhouse Gases (2016), along with the report's median assumption of 3% for the discount rate and 2.2% for the average growth rate of the SCC to calculate the benefits. We then used information on the administrative costs of the PSNP program from Drechsler et al. (2017) and benchmarked the social benefits of negative CO 2 emissions against the implementation costs of the program.

Change in tree cover
The PSNP increased tree cover on average by 3.8% (95% CI: 0.0006; 0.0777) in the participating PSNP districts in our study area, relative to what would be expected in the absence of the program ( Fig. 2A). This change represents a 0.54 percentage point increase in tree cover. To put this estimate in context, the mean tree cover in the non-PSNP pixels in the common support increased by 0.77 percentage points (from 14.27% to 15.04%) between 2002-2004 and 2017-2019. The impact of the PSNP on tree cover is additional to this; thus, in PSNP districts tree cover increased by 1.31 percentage points (this can be thought of as a change from approximately 14.27% to 15.58%). When we disaggregated by population density (Fig. 1C), the impact estimates were larger for less densely populated areas. Specifically, the estimated impact was 4.4% (95% CI: 0.0051; 0.0843) in areas with less than 300 people per km 2 and 6.0% (95% CI: 0.0142; 0.1078) when a more stringent threshold of 150 people per km 2 was used ( Fig. 2A).
The positive impacts on tree cover were also larger on steeper sloped land areas (Fig. 2B). For the middle quintile (average slope ranging between 5.8 and 10.7 degrees), we estimated that the PSNP increased tree cover by 5.6% (95% CI: 0.0047; 0.1089). The estimated impact was largest at the 4th quintile (10.7 to 19.1 degrees) of the slope distribution; 7.5% (95% CI: 0.0352; 0.1161). The estimated impacts by terrain slope were consistently larger in magnitude when we restricted the area to less densely populated areas (Figs. A7 and A8 in the Appendix).
In addition, we documented statistically significant impacts in areas classified as forests and woody areas (see Table A3 in the Appendix for land cover types and their aggregation) and cropland, but not in grasslands or savannas (Fig. 2C). In forests and woody areas, we estimated that the PSNP increased tree cover by 11.4% (95% CI: 0.0169; 0.2195) and in croplands by 3.7% (95% CI: 0.0075; 0.0682). The magnitudes of the corresponding impact estimates were sizably larger when we restricted the area of analysis to less densely populated pixels (Figs. A9 and A10 in the Appendix).
Lastly, we also found that the impacts on tree cover were larger in districts that had a large number of PSNP beneficiaries relative to their total population compared to districts that had relatively fewer PSNP beneficiaries ( Fig. A11 in the Appendix). We did not detect statistically significant spillover effects to districts directly adjacent to the PSNP districts (see Fig. A12 in the Appendix).

Carbon sequestration and benefit-cost analysis
We estimated the carbon sequestered by increased tree cover over the period 2005-2019. Changes in tree cover were converted to sequestered carbon using the VCF-TC data and AGBM from the Global Forest Watch GFW, 2021 dataset (Zarin et al., 2016) in metric tons of biomass per Ha for the year 2000.
We found that the average increase in biomass per VCF-TC pixel due to the increased tree cover was 1.12 metric tons per ha, which is equivalent to an estimated 62.4 million metric tons of negative CO 2 emissions (95% CI: 1.1 to 113.5; note that this and subsequent confidence intervals are non-symmetrical relative to the point estimate, due to the nonlinear relationship between tree cover and AGBM). This is equivalent to 4.16 million metric tons annual negative CO 2 emissions.
We next estimated the benefit-cost ratio of the negative carbon emissions relative to the PSNP's cost we using the SCC estimate from the Interagency Working Group on Social Cost of Greenhouse Gases (2016) and information on the administrative costs of the PSNP program from Drechsler et al. (2017). We benchmarked the social benefits of negative CO 2 emissions against the implementation costs of the program under the four scenarios shown in Table 1. We note that other factors that can affect carbon sequestration estimates, including species type, carbon uptake rate, and planting location (Griscom et al., 2017;Holl and Brancalion, 2020;Kirby and Potvin, 2007;Schulp et al., 2008) were not considered.
The benefit-cost ratio of the carbon stored in the tree cover results indicate that the social benefits of the carbon sequestered by the program could offset as much as 49% of the administrative costs of the program (Table 1), although the magnitude of the carbon storage benefits depends heavily on how long the increase in tree cover is preserved.

Robustness checks
We conducted a series of robustness checks to assess the sensitivity of our results (see Section A.6 in the Appendix). First, accounting for the skewed nature of the tree cover data (see Fig. A4 in the Appendix), we used a natural logarithm of the tree cover as our outcome variable. This meant discarding 0.02% of observations with a zero tree cover value. Therefore, we reran our regression applying an inverse hyperbolic sine transformation (Burbidge et al., 1988) as well as using a raw tree cover variable instead of the logged variable. Our findings are robust to these alternative ways of defining our outcome variable, see Table A8 in the Appendix.
Second, in our main analyses, we used three-year averages of tree cover. To explore the sensitivity to the time period aggregation, we reestimated our model using annual tree cover data. We also checked whether our results held if we collapsed the data only to two time periods: pre-PSNP (2000)(2001)(2002)(2003)(2004) and PSNP (2005-2019) (see Table A9 in the Appendix). Our results were robust to these alternative ways of constructing our dataset.
Third, we used alternative ways to control for time trends. Instead of region-specific period fixed effects, we showed that our results are robust to using less data-intensive approaches, such as simple linear time trend (=1 if first period; =2 if second period; and so on) and uninteracted period fixed effects (see Table A10 in the Appendix).
Fourth, our impact estimates were not influenced by additional district level characteristics (Table A11) or time-invariant district characteristics (Table A12).
Fifth, to explore the possibility that our findings were driven by a particular district (e.g., due to its size or because of unusually large changes in tree cover after the launch of the PSNP), we reran our regressions by omitting each district one at a time. The estimates remained stable when individual districts were omitted from the dataset, indicating that the findings were not driven by a particular district (Fig. A13 in the Appendix).
Sixth, our confidence intervals were calculated using clustered standard errors, which may not be appropriate if the error terms are spatially auto-correlated. Fully adjusting for spatial autocorrelation is not computationally feasible in our setup due to the large size of the dataset (see Table A13 in the Appendix). However, using random subsets of our data suggests that our results are robust to adjusting our standard errors and confidence intervals to control for spatial autocorrelation (Conley, 1999) (see Table A14 in the Appendix).
Finally, we used MODIS VCF's data quality band as described in section 3 (Spatial Data) to assess the robustness of our results to the quality of the input surface reflectance data used to estimate tree cover. To this end, we identified the number of data quality flags of each pixel per year and re-ran the analysis after discarding all the pixels that were flagged twice or more in a year during our study period. Our findings were not driven by data quality issues (see Table A15 in the Appendix).

Discussion
The United Nations SDGs (United Nations, 2015) underscore the urgent need to address multiple dimensions of climatic, social, and ecological challenges in an integrated manner (Downing et al., 2021;Gil et al., 2019;Norton et al., 2020;Seddon et al., 2020).
Food security, poverty, and forests are closely linked and are affected by and contributors to climate change (FAO, 2008;IPCC, 2019). While higher food production is necessary to feed an increasingly populated world, the agricultural sector remains an important source of GHG emissions, deforestation, and negative environmental impacts (Agrawal et al., 2014;Bahar et al., 2020;Gil et al., 2019;Godfray et al., 2011;Knoke et al., 2013;IPCC, 2019). Forests support climate change mitigation through carbon sequestration and can also contribute to food security through the provisioning of ecosystem services and increased yields in agroforestry systems (Amadu et al., 2020;Bahar et al., 2020). Deforestation and land degradation contribute to climate change through GHG emissions and reduced rates of carbon uptake (IPCC, 2019), while poverty exacerbates food insecurity and increases vulnerability to climate change by reducing coping and adaptive capacity (FAO, 2008;Paul et al., 2016).
While Ethiopia has made considerable commitment to reduce its vulnerability to climate change, its climate change adaptive capacity is still limited and there is a need to strengthen it across sectors, interventions, and actors (Federal Democratic Republic of Ethiopia, 2019; Federal Democratic Republic of Ethiopia, 2021). The PSNP is designed as a safety net program for households that are chronically food insecure and poor while supporting community development and environmental restoration practices through its public works program (MoA, 2010;Wiseman et al., 2010). As such, it is an example of a program that integrates climate change actions into development programming.
Our results show that the PSNP increased tree cover by 3.8% on average over 15 years in the districts of the Ethiopian highlands that participated in the program. The estimated tree cover increases are larger in less densely populated areas, steep-sloped terrain, and areas Note: This table compares the social benefits of the negative carbon emissions from the PSNP against the PSNP program implementation costs. The table does not factor in any other PSNP program benefits such as poverty alleviation. Each row corresponds to a different assumption of how many years elapse before the trees are cut and carbon is released into the atmosphere.
8 classified as forests and croplands at the onset of the program. These heterogeneous impacts align with the literature on forest conservation policies showing that conservation practices tend to be more successful in areas in which the opportunity cost of converting land to agricultural use is relatively high (Angelsen, 2010;Börner et al., 2020): distance to markets and terrain slope drive agricultural income and costs, making it less profitable to convert forests to cropland in steep-sloped terrains and in areas farther away from urban centers (Sandel and Svenning, 2013;von Thünen, 1826). Therefore, conservation projects and programs often target areas in which the land conversion pressures are low to begin with (Barton et al., 2013;Joppa and Pfaff, 2009). This is also the case with the PSNP as the program targets rehabilitation efforts on communal lands and steep-sloped terrain that are typically not well suited for crop-agriculture (MoA, 2010;Wiseman et al., 2010). While in Ethiopia the land area allocated to crop agriculture has grown over the past two decades, a recent study suggests that the returns to converting more of the highlands to agriculture are limited because further agricultural expansion is not economically profitable (Schmidt and Thomas, 2018). We estimate that the annual negative CO 2 emissions from the increased tree cover are equivalent to 1.5% of Ethiopia's annual emissions reduction pledged by 2030 in its Nationally Determined Contribution for the Paris Agreement (Federal Democratic Republic of Ethiopia, 2021). Our estimate is larger, but on the same magnitude as Woolf et al. (2018) on the PSNP negative carbon emissions using different methods. 1 Our study focuses on a large-scale safety net program with a public works component that is increasingly viewed as an important part of Ethiopia's response to climate change (Federal Democratic Republic of Ethiopia, 2020; Wiseman et al., 2010). Related work in this area has focused on conditional cash transfer programs without explicit environmental goals and documented mixed environmental impacts (Alix-Garcia et al., 2013;Dyngeland et al., 2020;Ferraro and Simorangkir, 2020). Also relevant is research on India suggesting that a public works program led to negative environmental impacts in the form of increased air pollution (Behrer, 2019).
Our findings are complementary to the growing literature on the benefits of Payment for Ecosystem Services (PES) programs for forest restoration and deforestation reduction (Alix-Garcia et al., 2012;Alix-Garcia et al., 2015;Alix-Garcia et al., 2018;Jack and Jayachandran, 2019;Jayachandran et al., 2017;Salzman et al., 2018;Vincent et al., 2021). While PES programs are typically designed with environmental benefits as the primary goal, growing research has demonstrated that these programs can have important social benefits including poverty reduction (Alix-Garcia et al., 2015) and increasing social capital (Alix-Garcia et al., 2018). Our work complements the PES literature in that we find a program with primarily social goals (food security and poverty reduction) can have important environmental benefits. Combined with the PES literature, our findings further cement the importance of considering both social and environmental benefits when evaluating programs. At the same time, the relatively modest gains that we find in tree cover from the poverty-focused PSNP are congruent with the relatively modest reductions in poverty that have been found from PES programs (Alix-Garcia et al., 2015).
The social protection literature has raised concerns about the high implementation costs of public works programs, especially when benchmarked against alternative social safety net programs, such as universal basic income schemes (Ravallion, 2019). However, typically public works programs have not accounted directly for the benefits generated by the public goods produced by these programs (Beierl and Grimm, 2019;Gehrke and Hartwig, 2018;Ravallion, 2019;Subbarao et al., 2012). Our estimates suggest that for Ethiopia's PSNP, the positive impact of tree cover alone (through carbon storage) could offset as much as 49% of the administrative costs of the program on the long term. Our findings show that public works programs can have sizable environmental benefits and should be embedded in benefit-cost calculations to avoid under investing in beneficial programs.

Considerations in realizing environmental benefits in social protection programs
Potential pathways to increase the environmental and climatic benefits of social protection programs include adding or strengthening existing environmental components and incentives to increase tree cover and to perform sustainable land management practices. In addition, there is an opportunity for these types of programs to build on synergies with tree planting, forest conservation, and sustainable forest management initiatives at national (e.g., the African Forest Landscape Restoration Initiative, the Green Legacy Initiative, and the Climate Resilient Green Economy strategy) and sub-national level, such as participatory forest management (Ameha et al., 2014;Siraj et al., 2018), Clean Development Mechanism projects (Brown et al., 2011), and other re-greening initiatives (Lemenih and Kassa, 2014). However, realizing the environmental and climatic benefits of social protection programs that have an environmental component is not without challenges, as it requires the full integration of the programs within their socioecological context. Specifically, the success of tree planting projects rests on careful planning, evaluation of potential trade-offs, and consideration of several social and environmental factors before their implementation (Chazdon and Brancalion, 2019;Holl and Brancalion, 2020). 2 Among these, biophysical aspects such as selecting adequate species and location have received most of the attention (Boissière et al., 2021) partly due to their effect on carbon stock as well as the effect of trees in the environment and overall climatic impact (Anderson et al., 2011;Kirby and Potvin, 2007;Schulp et al., 2008). In addition, the importance of including many native species to increase biodiversity and ecosystem services provisioning, has also been emphasized (César et al., 2021;Ellison et al., 2017;IPCC, 2019;Seddon et al., 2020). It is also critical to consider the complexity of socio-economic aspects of tree planting, including a long-term commitment to land protection, management, and funding (Holl and Brancalion, 2020), as well as the needs, goals, and participation of local communities (Boissière et al., 2021), and land tenure issues (Agrawal et al., 2014;Boissière et al., 2021;Legesse et al., 2018;Unruh, 2008).
The design of the PSNP public porks component paid careful attention to many of these issues. First, the public works projects are integrated into community planning to increase their relevance and improve long-term sustainability (MoA, 2010;Wiseman et al., 2010). Combined with technical support from environmental experts (Wiseman et al., 2010), this community-led approach aimed to ensure that public works projects were tailored to the socio-ecological context. Second, the PSNP public works take place during the agricultural slack season to minimize potential crowding-out effects of on-farm labor and output (Holden et al., 2006). Third, while the PSNP remains largely externally funded (World Bank, 2018), the program is led and implemented by the government of Ethiopia, ensuring long-term commitment to implementation and results.
Lastly, we note that the primary focus of the PSNP is on improving food security, with a secondary focus on generating community assets. Many of the community assets aim to enhance climate change 1 Note that we benchmark our estimate against Ethiopia's most recent 2021 reduction pledge (Federal Democratic Republic of Ethiopia, 2021), whereas Woolf et al. (2018) benchmark against Ethiopia's earlier 2016 reduction pledge UNFCCC (2016), so the percentages of the reductions met are not directly comparable. adaptation and resilience. Climate change mitigation, in contrast, has not been a core focus of the program and, as a result, the program strategies are not designed to optimize climate change mitigation benefits. Given that Ethiopia is a resource-poor country with limited implementation capacity, the burden of mitigation should not fall on Ethiopia. However, if external funding to the PSNP were increased, that might allow for a greater focus on climate change mitigation (Jirka et al., 2015).

Limitations
Many social protection programs in low and middle income countries use geographic targeting Coady et al., 2004) making it difficult to causally assess their environmental impacts. We addressed this by constructing a credible counterfactual, however, in the absence of a randomized allocation of the program, we cannot be sure that our estimates are entirely free from bias (Alpízar and Ferraro, 2020). Another limitation of our methodological approach is that our estimates represent local average treatment effects (Imbens and Angrist, 1994). More specifically, by restricting the study area to common support (a set of PSNP and non-PSNP pixels with similar agro-ecological and socioeconomic characteristics at the onset of the program), our impact estimates are identified from a sub-set of the area covered by the PSNP. This may raise a question whether our local average treatment estimates can be generalized to represent overall impacts (i.e., average treatment effects of the whole program) (Deaton, 2010;Imbens, 2010). We note that the area of common support pixels used in the final analysis was quite large, covering 345,000 km 2 (or 34.5 million ha) and spanning multiple agro-ecological zones. Focusing on such a large land area eases the concerns of applying local average treatment effects to make broader policy relevant statements. To further alleviate these concerns, we also uncover several important sources of impact heterogeneity, such as population density, terrain slope, forest and cropland land cover, that are likely to be highly relevant to policy makers designing similar programs.
While our data did not allow us to explore the mechanisms through which the PSNP increased tree cover, we hypothesize that the tree cover increases are due to the nature of the public works projects which were designed to rehabilitate degraded lands. However, it is also possible that the public works projects 'crowd-in' investments by inducing households to plant trees on their private lands (Andersson et al., 2011;Holden et al., 2006). If so, this means that the effect of the PSNP on tree cover goes beyond the PSNP's public works area. In our study, we did find that the PSNP resulted in small increases of tree cover in areas categorized as cropland at the onset of the program in 2005, while the largest effects were observed in areas classified as forests and woody areas. Additionally, it is also possible that the cash or in-kind transfers themselves could have limited the pressure on households to cut and sell trees for their immediate cash needs during economic hardship, preventing deforestation.
Finally, our analysis was retrospective in nature covering a 15-year period between the onset of the PSNP and 2019. We should be careful in extrapolating our findings beyond the study period. The COVID-19 pandemic and the large-scale conflict in Ethiopia that both erupted in 2020 have caused major disruptions to the PSNP. Moreover, the recent international events (e.g., the pandemic and the war in Ukraine) have resulted in economic downturns and sizable government deficits in high income countries raising a risk of future funding cuts to the PSNP. However, it is very promising that in early 2021 funding for the PSNP was renewed for another five-year phase (U.S. Embassy Ethiopia, 2021).

Conclusion
We have measured the impact on tree cover of the PSNP using several datasets including satellite-based data of tree cover combined with difference-in-differences and inverse probability treatment weighting methodologies. To the best of our knowledge, this is the first assessment of the environmental impacts of a major public works program using broad geographic data coverage and counter-factual analysis. It is an example of a design-based causal inference strategy to empirically evaluate a large sustainability intervention (see Barrett (2021) for a recent call to expand this type of analysis in the broader sustainability science community). Our work also buttresses Norton et al. (2020)'s review on the potential of employment-based social assistance programs to promote ecosystem stewardship, in that we quantify the carbon sequestration benefits of the PSNP showing that large social assistance programs can attain both social and environmental aims.
Our results show that the PSNP increases tree cover and supports climate change mitigation efforts through carbon sequestration, with larger increases in less densely populated areas and on steep-sloped terrain. The PSNP is one of the largest social protection programs in Africa, and our results show the potential that these types of programs can have to support mitigation strategies for climate change by increasing tree cover and reducing CO 2 emissions.

Funding
We acknowledge funding from the CGIAR Research Program on Policies, Institutions, and Markets (PIM). The authors declare that they have no competing interests. The funder had no role in the conceptualization, study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Table A1 describes the spatial datasets used in this study. The main dataset used in this study is the Vegetation Continuous fields (VCF) annual dataset (MODIS44B) L3, Collection 6 derived from the MODerate Resolution Imaging Spectrometer (MODIS) sensor on board the Aqua and Terra satellites. This global dataset has a 250 m spatial resolution and provides an estimation of three ground cover components in each pixel: percent tree cover (VCF-TC), percent non-tree vegetation, and percent non-vegetated (bare) from 0 to 100 (Townshend et al., 2017;DiMiceli et al., 2021). The ground cover components are estimated through a regression tree algorithm using training data from Landsat Geocover data, 16-day surface reflectance composites including bands 1-7 and brightness temperature from bands 20, 31, and 32, and the MODIS Global 250 m Land/Water map (Townshend et al., 2017). In addition to these variables, the VCF dataset also includes a cloud cover band, a data quality band, and two standard deviation bands (percent tree cover and percent non-vegetated) bands (DiMiceli et al., 2015;Townshend et al., 2017).

A.1. Spatial data and methods
The Global Forest Change (GFC) dataset (Hansen et al., 2013) is also widely used to assess forest change (Jain, 2020). We used VCF-TC for two main reasons. First, although GFC covers our study period at a higher spatial resolution (30 m), we did not consider it appropriate for our analysis due to the inconsistencies resulting from differences in data processing between the periods 2000 to 2012 and 2013 to 2019 (University of Maryland, 2019). Second, while GFC provides percent tree canopy cover for 2000, the remaining years are coded as a binary variable (either forest gain or loss). VCF-TC is better aligned to our research objectives because it allows assesing forest change as a continuous process at the pixel level (DiMiceli et al., 2021;Ryan et al., 2017).
Our outcome variable of interest is the percent tree cover (VCF-TC), defined as the "amount of skylight obstructed by tree canopies equal to or greater than 5 m in height" (Hansen et al., 2003). We note that this differs from crown cover, "the amount of the ground which is encompassed by the tree's crown regardless of whether light penetrates." (Townshend et al., 2017).

Table A1
Spatial datasets: data source, time period used in the analysis, and spatial resolution.   (2018).
VCF data are provided as discrete tiles in sinusoidal projection. We mosaicked the four tiles covering Ethiopia (h22.v08; h22v07; h21v08 and h21v07) for all years in our study period. All data except land cover (already in sinusoidal projection) were projected to the sinusoidal projection. There are six land cover type classifications available in the MODIS product MCD12Q1. This is also a global yearly product, but at a 500 m spatial resolution. We chose the Geosphere-Biosphere Programme (IGBP) 17 land cover type classification scheme (Table A2) because it is more closely aligned with our research focus, and aggregated land cover types as described in Table A3 to explore the heterogeneity of the PSNP impacts on the main land cover types.
Each VCF-TC pixel was matched to a district and region using the 2007 Ethiopia's Central Statistical Agency (CSA) (unpublished data) administrative boundaries. The latter was joined to the annual PSNP caseloads at district level for our study period. We identified the district splits that occurred for each year after 2007 and merged back the child districts to their parent districts along with the PSNP beneficiaries to generate a spatially consistent dataset.
Finally, we used the VCF data quality band to test the sensitivity of the results to the uncertainty in vegetation estimates associated with input data quality. This involved processing the VCF-TC quality flags for each year in our study period and reclassifying them to extract the pixels that had two or more flags per year during our student period. The flagged pixels were excluded from the analysis for the robustness check.

A.2. Impact assessment method
The key challenge of any impact assessment is the construction of the counterfactual; what the outcome would have been had the districts not received the program. In randomized controlled trial (RCT) designs, this is solved by randomly allocating the treatment (here PSNP) across eligible districts. When program allocation is random, districts assigned to the control arm are identical-in expectation-to districts in the treated group before the onset of the program, so these control districts provide a credible counterfactual. Impacts of the program can then be measured as differences in outcomes (or differences in changes in outcomes over time) between the randomly assigned treatment and control districts. When an RCT or another experimental design is not feasible or ethical, an identification strategy must be developed in which the counterfactual is constructed using statistical techniques to create a control group of districts that are as similar as possible to the treated group. Most social safety net programs across low and middle income countries are targeted to poor people or poor areas, not randomly allocated (Coady et al., 2004). This is also the case for the PSNP: the program was geographically targeted to chronically food-insecure areas of the country (Wiseman et al., 2010). In the absence of an experimental design, we combined difference-in-differences (DiD) and statistical matching methods to estimate the impact of the PSNP on tree cover. This approach is credible because due to funding constraints or spatial inertia, many poor and chronically food-insecure districts in the highlands are not part of the PSNP and instead make recurring annual requests for emergency food assistance (Clay et al., 1999;Jayne et al., 2002;NDRMC, 2018;World Bank, 2020).

A.2.1. Difference in differences method
DiD is a widely used quasi-experimental method to estimate treatment effects when a randomized allocation of a policy or program is not feasible or ethical (Angrist and Pischke, 2009). DiD requires data before and after the intervention began and from a group that was subject to the treatment (treated group) and a group that was not (control group). A key identifying assumption of DiD is that the two groups were on a similar trend before the treatment began (Ryan et al., 2019). To test this 'parallel trend hypothesis', data from at least two periods before the intervention began is needed.
With data before and after the PSNP began from PNSP and non-PSNP districts, we can use the DiD approach to estimate the impact of the PSNP program on tree cover. The VCF-TC data are available from 2000 onwards, permitting us to test the parallel trends hypothesis.
To begin, we tested the parallel trend hypothesis, restricting our data to two periods prior to the launch of the PSNP in 2005 (2000-2001 and 2002-2004) and defining the binary 'treatment' variable to equal one if the period was 2000-2001, and to equal zero if the period was 2002-2004. We first estimated the equation provided in the main text using an ordinary least squares (OLS) method. Column 1 in Table A4 shows that the null hypothesis of parallel trend (β = 0) was comfortably rejected (p < 0.001). We then attempted to adjust for non-parallel trends using various fixed effects estimators. Columns 2 and 3 show that the parallel trends hypothesis was also rejected when both district-level (p < 0.001) and pixel-level (p < 0.001) fixed effects were used.

A.2.2. Statistical matching method
We used statistical matching estimators that are frequently used in the environmental conservation literature to estimate program impacts Naidoo et al., 2019;Pfaff et al., 2015;Soares-Filho et al., 2010), and which have shown to perform well in reducing bias when combined with DiD in various contexts (Chabé-Ferret, 2017;McKenzie et al., 2010;Ryan et al., 2019). More specifically, we used a propensity score matching algorithm (Rosenbaum and Rubin, 1983) to match the PSNP and non-PSNP pixels based on their pre-program characteristics: mean and standard deviation of annual rainfall in 1995-2004 (and their squared terms), population density in 2005, elevation and slope of land (Table A5). We also included binary indicators for each region. Table A6 presents the results of the propensity score estimation based on a logit estimation method in which the dependent variable is a binary variable that equals one if the pixel belongs to a PSNP district, and zero otherwise. As is common in this literature (Imbens, 2015), we were less interested in interpreting the magnitude or statistical significance of the coefficients reported in Table A6. Instead, we used the predictions from this model to construct the propensity score. Fig. A1 shows the distribution of the propensity score for both PSNP and non-PSNP pixels. As expected, there were a large number of non-PSNP pixels that received a very low score, indicating that they are very unlikely to be selected into the program based on their agro-ecological characteristics. Similarly, there were many PSNP pixels for which the probability of selection was close to one. Spatially, we see that these 'poor matches' are primarily located in the east and west of the study area (Fig. A2). We defined the area of common support as pixels with the estimated propensity score within the interval [0.1; 0.9] (Crump et al., 2009). This meant discarding 5.8 million pixels. In the final dataset used in the analysis, we have 6.5 million pixels and reasonable overlap in the propensity score distributions across PSNP and non-PSNP pixels (Fig. A3). Spatially, this meant focusing on the areas in the middle of the study area; those areas just inside and outside of the 'PSNP boundary' (Fig. A2) where the agro-ecological conditions are comparable (see Fig. 1 in the main text). Restricting the area to common support and rerunning the models based on OLS and the two fixed effects methods resulted in smaller coefficients in absolute terms, but the parallel trend hypothesis is rejected (p < 0.001) across columns 4 to 6 in Table A4. Note: The unit of observation is a pixel observed in two time periods: 2000-2001 and 2002-2004. The outcome variable is the mean percent of tree cover in each period. The standard errors are reported in parentheses and they are clustered at the district level in columns 1, 2, 4, 5, and 7 and at the pixel level in columns 3 and 6. All models include a binary variable capturing pixels belonging to PSNP districts (except the models based on fixed effects methods), region specific period fixed effects and mean annual rainfall over the period. OLS = Ordinary Least Squares; IPTW = Inverse Probability Treatment Weighting. Statistical significance denoted at *** p<0.01, ** p<0.05, * p<0.10.

Table A5
Descriptive statistics of pre-program matching covariates.
(1) PSNP (2)  Note: Mean values followed by standard errors in parentheses. The value displayed for t-tests are the differences in the means across the two groups. Standard errors (SE) are clustered at district level. Statistical significance of the t-test (last column) denoted at *** p<0.01, ** p<0.05, * p<0.10. 0/1 refers to binary variable.

Table A6
Propensity score regression results.    Finally, we used these pixel level propensity scores (PS) to calculate inverse probability treatment weights (IPTW) (Abadie, 2005;Joffe et al., 2004): 1/PS for the treated (PSNP) pixels and 1/(1 − PS) for the untreated (non-PSNP) pixels. After restricting the pixels in our dataset to common support and applying IPTW on our regression model, we cannot reject the null hypothesis that the tree cover in PSNP and non-PSNP areas were on a similar trend before the PSNP was launched in 2005 (Column 7 in Table A4); p = 0.473. Table A7 further shows that the pre-program matching covariates are in balance after we restrict the area to common support and apply IPTW.

A.3. Additional descriptive statistics
Figs. A4,A5,A6 display the distribution of average tree cover, population density, and terrain slope, respectively, after restricting to the area of the common support.

Table A7
Covariate balance after restricting the area to common support and applying inverse probability treatment weights.
(1) PSNP (2)  Note: Mean values followed by standard errors in parentheses. The value displayed for t-tests are the differences in the means across the two groups. Standard errors (SE) are clustered at district level. Observations are weighted using inverse probability treatment weights. Statistical significance of the t-test (last column) denoted at *** p<0.01, ** p<0.05, * p<0.10. 0/1 refers to binary variable.

A.4. Heterogeneity analyses A.4.1. Heterogeneity by population density
In the main text, we provide estimates of the impact of the PSNP on tree cover by terrain slope (Fig. 2B). We replicate that analysis, but restrict the area to rural areas using two population density thresholds. Fig. A7 shows the estimates when rural areas are defined as areas with less than 300 people per km 2 and Fig. A8 shows the estimates when rural areas are defined as areas with less than 150 people per km 2 . In line with Fig. 2 presented in the main text, we see that the slope-specific impacts are larger when we move to less densely populated areas (especially areas with <150 people/km 2 ).

Fig. A7.
Estimates measure % change in tree cover due to the PSNP. Area restricted to pixels containing less than 300 people/km 2 . The unit of observation is a pixel observed periodically. Impact estimates for terrain slope quintiles: 0-20 percentile (0.0 to 2.9 degrees; N = 9,924,180); 20-40 percentile (3.3 to 5.6 degrees; N = 7,496,279); 40-60 percentile (5.8 to 10.7 degrees; N = 8,217,594); 60-80 percentile (10.7 to 19.1 degrees; N = 8,568,854); 80-100 percentile (19.1 to 78.4 degrees; N = 8,771,077). Estimates are based on a difference in differences method combined with an inverse probability weighting. Confidence intervals are computed from standard errors clustered at the district level. Fig. 2C in the main text provides the estimates by land cover type at the onset of the program in 2005. We replicate that analysis, but restrict the area to rural areas using two population density thresholds. Fig. A9 shows the estimates when rural areas are defined as areas with less than 300 people per km 2 and Fig. A10 the estimates based on the 150 people per km 2 threshold. As before, we see that the slope-specific impacts are considerably larger when we move to less densely populated areas (<150 people/km 2 ). This is particularly so for the forests and woody area category and croplands where we find that the PSNP increased tree cover by 15.0 percent and 6.9 percent, respectively.

Fig. A8.
Estimates measure % change in tree cover due to the PSNP. Area restricted to pixels containing less than 150 people/km 2 . The unit of observation is a pixel observed periodically. Impact estimates for terrain slope quintiles: 0-20 percentile (0.0 to 2.9 degrees; N = 8,366,596); 20-40 percentile (3.3 to 5.6 degrees; N = 6,204,891); 40-60 percentile (5.8 to 10.7 degrees; N = 6,670,860); 60-80 percentile (10.7 to 19.1 degrees; N = 7,010,199); 80-100 percentile (19.1 to 78.4 degrees; N = 7,449,533). Estimates are based on a difference in differences method combined with an inverse probability weighting. Confidence intervals are computed from standard errors clustered at the district level. Fig. A9. Estimates measure % change in tree cover due to the PSNP. Area restricted to pixels containing less than 300 people/km 2 . The unit of observation is a pixel observed periodically. Impact estimates for different land cover types at the onset of the program in 2005: Forests and woody areas (N = 1,847,615); Croplands (N = 15,959,503); Grasslands (N = 17,587,542); Savannas (N = 7,574,035). See Table A3 for exact aggregations. Estimates are based on a difference in differences method combined with an inverse probability weighting. Confidence intervals are computed from standard errors clustered at the district level.  Table A3 for exact aggregations. Estimates are based on a difference in differences method combined with an inverse probability weighting. Confidence intervals are computed from standard errors clustered at the district level.

A.4.3. Heterogeneity by caseload intensity
We also explored whether the impacts were larger in districts that had more PSNP beneficiaries relative to total population compared to districts that had fewer. To do this analysis, we computed the average number of PSNP beneficiaries in each district over the study period and divided this number by the total population of the district. Using this variable as our measure of beneficiary caseload intensity, we split the pixels originating from the PSNP districts into two groups using the median caseload intensity as the threshold. We then replaced our treatment variables with these two binary variables and reran the regression. The results of this analysis are shown in Fig. A11. For pixels with a lower caseload intensity, we find positive point estimates of the impact of the program, but these estimates are relatively small, ranging from 1.6 to 3.0 percent depending on the population threshold used, and are not statistically significant. On the other hand, for the pixels with a higher caseload intensity, the estimated increase in tree cover is larger (ranging from 7.1 to 10.9 percent) and statistically significant in all specifications. Wald tests further confirmed that the differences in impact estimates between low and high intensity areas were statistically different from zero in all three regressions (p < 0.05). These results are reassuring in that it is participation in the PSNP, and not some other omitted factor, which is driving our main results.

A.5. Spillover analysis
To assess spillovers from PSNP districts to neighboring non-PSNP districts, we split our control area into two groups: non-PSNP districts directly adjacent to PSNP district and other non-PSNP districts. In total, there were 134 adjacent non-PSNP districts that shared a border with at least one PSNP district. We re-ran our regression using two binary treatment variables: one capturing PSNP and the other capturing adjacent non-PSNP districts. A positive and significant impact estimate on the variable capturing adjacent districts would indicate that neighboring non-PSNP districts also benefit from the program. The estimates reported in Fig. A12 quantify the change in tree cover relative to non-PSNP districts that do not share a border with a PSNP district. The impact estimate for the PSNP districts is statistically significant in all three columns, while the estimate for the adjacent non-PSNP districts is not. Therefore, we conclude that there is no statistical evidence in favor of spillover to non-PSNP districts.

Fig. A11.
Estimates measure % change in tree cover due to the PSNP. Impacts are estimated separately for low caseload intensity versus high caseload intensity pixels. A pixel is defined as high caseload intensity if its average caseload per capita over our study period is above the median. The unit of observation is a pixel observed periodically. The figure displays separate panels for all pixels (N = 45,229,114), pixels containing less than 300 people/km 2 (N = 42,977,984), and pixels containing less than 150 people/km 2 (N = 35,702,079).

A.6. Robustness checks
We conducted a series of robustness checks to assess the sensitivity of our results.

A.6.1. Alternative ways of defining the outcome variable
Our results are not sensitive to the way we define our outcome variables. First, accounting for the skewed nature of the tree cover data (see Fig. A4), we used a natural logarithm of the tree cover as our outcome variable. This meant discarding 0.02% of observations with a zero tree cover value. Therefore, we re-estimated our regression applying an inverse hyperbolic sine transformation (IHS) as well as using a raw tree cover variable instead of the logged variable. The results reported in Panel B of Table A8 show that the estimates are near identical to those reported in the main text (reproduced in Panel A of Table A8) when we use the IHS transformed outcome variable. However, the estimate for 'all pixels' is only significant at the 10% level. The estimates are statistically significant when we use the non-transformed (or raw) tree cover variable (Panel C of Table A8), although as before the 'all pixel' estimate is significant at the 10% level. The magnitudes are also comparable to those reported in the main text. We find that the Fig. A12. Estimates measure % change in tree cover due to the PSNP. Impacts are estimated separately for PSNP districts and non-PSNP districts adjacent to PSNP districts ('Adjacent', N = 134 districts), against other, non-adjacent, non-PSNP districts. The unit of observation is a pixel observed periodically. The figure displays separate panels for all pixels (N = 45,229,114), pixels containing less than 300 people/people/km 2 (N = 42,977,984), and pixels containing less than 150 people/ km 2 (N = 35,702,079). PSNP increased tree cover by 0.497 percentage points. Considering the mean tree cover percent in non-PSNP districts before PSNP was launched (14.27), this estimate corresponds to a 3.48 percent increase in tree cover.
Second, in our main analyses, we used three-year averages of tree cover. To explore the sensitivity in this regard, we re-estimated our model using annual tree cover data. We also checked whether our results hold if we collapsed the data to two time periods: pre-PSNP (2000)(2001)(2002)(2003)(2004) and PSNP (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019). The results reported in Panel B and C of Table A9 show that the estimates are very similar to those reported in the main text (reproduced in Panel A of Table A9) and and statistically significant at least at the 10% level. Note: The unit of observation is a pixel observed periodically. The outcome variable is mean percent of tree cover in the period. In Panel A, the outcome variable is logged. In Panel B, the outcome variable is inverse hyperbolic sine transformed. In Panel C, non-transformed tree cover is used. The standard errors are reported in parentheses and they are clustered at the district level. All models include a binary variable capturing pixels belonging to PSNP districts, region specific period fixed effects and mean annual rainfall over the period. *** p<0.01, ** p<0.05, * p<0.10.

Table A9
Sensitivity analyses: data structure. (1) Area: All <300 ppl/km 2 <150 ppl/km 2 Panel A: Periodic data (Fig. 2 A.6.2. Controlling for common shocks We used alternative ways to control for common shocks. Instead of region-specific period fixed effects, we explored sensitivity by using less dataintensive approaches, such as a simple linear time trend (=1 if first period; =2 if second period; and so on) and un-interacted period fixed effects. The results reported in Panel B and C of Table A10 show that the magnitudes of the estimates are similar to those reported in the main text (reproduced in Panel A of Table A10) and and statistically significant at least at the 10% level. Note: The unit of observation is a pixel observed periodically or annually. The outcome variable is the (log) mean percent of tree cover in each period or year. The standard errors are reported in parentheses and they are clustered at the district level. All models include a binary variable capturing pixels belonging to PSNP districts, region specific period or year fixed effects and mean annual rainfall over the period or mean annual rainfall. Panel A: Estimates reported the main text; see Panel A of Fig. 2. Panel B: Annual data used instead of periodic data. Panel C: Data collapsed to two periods: pre-PSNP (2000)(2001)(2002)(2003)(2004) and PSNP (2005PSNP ( -2019. *** p<0.01, ** p<0.05, * p<0.10.

Table A10
Sensitivity analyses: time trends.

Note:
The unit of observation is a pixel observed periodically. The outcome variable is mean percent of tree cover in each period. In Panel A, region specific period fixed effects are used. In Panel B, these are replaced by simple time trend (=1 if first period; =2 if second period; and so on). In Panel B, these are replaced by period fixed effects. The standard errors are reported in parentheses and they are clustered at the district level. All models include a binary variable capturing pixels belonging to PSNP districts, region specific period fixed effects and mean annual rainfall over the period. In panel C, a binary variable obtaining a value 1 if the period is after the launch of PSNP (i.e., in 2005PSNP (i.e., in -2019, and zero if before (i.e., in 2000-2004). *** p<0.01, ** p<0.05, * p<0.10.

A.6.3. Observed and unobserved district characteristics
Both our propensity score model and the IPTW regression model contain covariates that are only defined at the pixel level. This may raise a concern that our model is not correctly specified if district level characteristicsbeyond pixel level characteristicsinfluence program selection. Our estimation approach offers two ways to address this: adding district level variables to the propensity score model and introducing district fixed effects to the regression model. The 'doubly robust' feature of the IPTW estimator means that our approach is valid if either the PS model or the regression model is correctly specified (Sant'Anna and Zhao, 2020). To explore this, we sequentially adjusted both models to assess whether our impact estimates are sensitive to the addition of additional district level controls. We first appended our propensity score model with additional variables capturing district level means of population density, slope and elevation. Re-estimating the model specified in the main text based on these revised propensity scores yields similar coefficients to those reported in the main text (Table A11). We then used the alternative way of controlling for time-invariant -and unobserved-district characteristics by introducing district fixed effects to the regression model. We implemented the fixed effects by appending the main model with binary variables for each district. This IPTW fixed effects estimator yields identical impact estimates to those estimated by the main IPTW estimator (Table A12) We also verified that our findings were not driven by a particular district (e.g., due to its size or because of unusually large changes in tree cover after the launch of the PSNP). To do this, we reran our regression by omitting one district at a time from the dataset. The results of this analysis are presented graphically in Fig. A13. In this figure, the blue line represents the coefficient estimate when a given numbered district is dropped, the shaded gray area represents the 95% confidence intervals. We ran this district exclusion exercise across all pixels and over all rural pixels (restricted to <300 ppl/km 2 or <150 ppl/km 2 ). As can be seen from the figure, our point estimates remain relatively stable through this sensitivity test.

Note:
The unit of observation is a pixel observed periodically. The outcome variable is mean percent of tree cover in each period. In Panel A, the equation reported in the main text is estimated. In Panel B, the inverse probability treatment weights are based on propensity scores estimated using additional district level variables. The standard errors are reported in parentheses and are clustered at the district level. *** p<0.01, ** p<0.05, * p<0.10.

Table A12
Sensitivity analyses: adding district fixed effects to the IPTW regression model.

A.6.4. Controlling for spatial autocorrelation
The confidence intervals reported in the main text are calculated using clustered standard errors. This may not be valid if the error terms exhibit significant spatial autocorrelation. The standard approach to address this in the literature is to use Conley standard errors that are robust to both spatial autocorrelation and heteroskedasticity (Conley, 1999). The Conley approach is based on a weighting matrix that places more weight on observations located closer to each other. These weights decay to zero after a user-specified distance cutoff. Unfortunately, with more than 6 million pixels and 40 million observations, calculating Conley-type standard errors is not computationally feasible. To demonstrate this, we used the userwritten Stata command acreg (Colella et al., 2019) that computes Conley standard errors while permitting the use of probability weights. We then selected small random subsets of pixels from our data and estimated the duration it takes for a standard laptop in 2021 (Quad core processor, 1.   A13. Estimates measure % change in tree cover due to the PSNP. The blue line represents the estimated percent change in tree cover when a given numbered district is dropped from the dataset, and the shaded gray area represents the 95% confidence interval for this estimate. We use the formula (exp(b) − 1)* 100 to convert from our regression estimates to percent changes; as a result, confidence intervals are slightly asymmetrical. The figure displays separate panels for all pixels, pixels containing less than 300 people/km 2 , and pixels containing less than 150 people/km 2 .

Table A13
Computer processing time when calculating Conley standard errors, by different random samples of pixels.
(1) Note: This table shows the estimated processing times when we estimated regressions based on Conley standard error adjustments (Conley, 1999) using the userwritten acreg command in Stata with small random samples of all pixels in our dataset. The parent dataset is defined as all pixels in the common support with population density <300 people/km 2 . The first column shows the % of pixels selected, the second is the number of pixel-time period observations, and the third is the number of pixels in the given sample. The remaining columns indicate the time it took for a standard laptop to estimate the regression.
GHz; 32 GB RAM) to run the regressions. Table A13 shows the results. First, we see that the processing time increases exponentially with the sample size. Second, using just 1% of all pixels (N = 428,679; 61,254 pixels) took more than 11 h. We also tested this on a slightly larger subset of 5% of all pixels using a high-end computer with more processing power (Quad core processor, 3.60 GHz; 24 GB RAM). The processing time in this case was 15,465 min, or 10 days and more than 17 h. Considering all this, computing Conley standard errors using the full set of data would take several months. We therefore settled for using these random subsets of pixels to gauge how the standard errors change when we use the Conley adjustment compared to when clustered standard errors are used. Focusing on rural areas defined as population density below 300 people/km 2 , we used four different distance cutoff values (50 km, 100 km, 200 km, and 500 km) to calculate the Conley adjusted standard errors. Table A14 shows the results. As expected, the standard errors decrease as the size of the sample increases. Interestingly, the standard errors seem to stabilize already when the subset covers at least as little as 0.25% of all pixels. Focusing on the results based on 0.25% or more pixels, we see that the Conley adjusted standard errors are somewhat larger than the clustered standard errors when 50 km or 100 km cutoffs are applied but similar in magnitude or smaller when the cutoff values are larger. However, the Conley adjusted standard errors do not render the estimates statistically insignificant even when 50 or 100 km cutoff values are applied.

A.6.5. VCF Data quality
The MODIS surface reflectance data used to estimate tree cover comes with quality indicators indicating if a pixel in any of the 8 input data periods used to generate the annual product is flagged as poor quality due to clouds, high aerosol levels, cloud shadows, or having a view zenith angle higher than 45 • (Townshend et al., 2017). It is considered that estimates of vegetation cover with two or more flags in a year may be erroneous and should be used with caution (Townshend et al., 2017). To assess the sensitivity of our results in this regard, we reran the analysis after discarding all pixels that were flagged twice or more any year during our study period. Our findings are robust to restricting the dataset to pixels that never had such data quality concerns during the 2000-2019 study period (Table A15).

Table A14
Estimated Conley standard errors, by different random subset of (<300 people/km 2 ) pixels. (1) (2) (3) Note: Panel A shows estimates based on <300 ppl/km 2 pixels as described in the main text. The subsequent panels use random subsets of pixels as labeled. The first row in each panel ('Clustered standard errors') shows the impact estimates ('coeff'), standard errors ('std err'), t-value (t) and p value (p). The remaining rows in Panels B to F show the impact estimates, standard errors, t-values and p values when Conley standard errors are computed with different distance cutoffs.

A.7.1. Carbon sequestration
We calculated the carbon sequestered resulting from increased tree cover from the pooled impact of the program across all districts that participated in the PSNP over the period 2005-2019 in our study area. We estimated that, controlling for rainfall, tree cover in the PSNP districts increased by 3.8% (95% CI: 0.0006; 0.0777) from 2000-2004 levels, relative to changes in tree cover in non-PSNP districts during the same period. Looking at the PSNP districts, we calculated from our VCF-TC data that their average tree cover in 2005 was 8.76%, and hence the 3.8% increase in tree cover due to the PSNP corresponds to a predicted final tree cover value of 9.10% (95% CI: 0.0864, 0.0944).
To convert these changes in tree cover to changes in carbon emissions, we used data on average aboveground live woody biomass (AGBM) from the Global Forest Watch Data (Zarin et al., 2016). This data was also used by (Ferraro and Simorangkir, 2020). The AGBM is distributed in tiles and provides data on the metric tons of biomass per ha at approximately 30 m spatial resolution for the year 2000. To convert from tree cover to tons of biomass, we analyzed the AGBM and VCF-TC data from 2000 and noted that for pixels with 8.76% tree cover in 2000, the average AGBM for those pixels is 30.8 metric tons of biomass per ha. In addition, the average AGBM for VCF-TC pixels with 9.10% tree cover is 31.9 metric tons of biomass per ha. Thus, we calculated that the average increase in biomass per VCF-TC pixel due to the PSNP was 1.1 metric tons per ha. Next, we converted these changes in biomass to negative CO 2 emissions. To begin, we multiplied the average AGBM by 0.5, because biomass is composed of approximately 50% carbon (Penman et al., 2003). We next multiplied the result by 3.67 to convert from tons of carbon to tons of CO 2 , based on the relative molecular weights of carbon and CO 2 . Finally, we scaled up by the total area of all the districts eligible for the PSNP (30.4 million ha), to calculate that the program resulted in 62.4 million metric tons of negative CO 2 emissions (95% CI: 1.1 to 113.5; note that this and subsequent confidence intervals are non-symmetrical relative to the point estimate, due to the nonlinear relationship between tree cover and AGBM).
Annualizing our estimates over the 15 years in our study period during the PSNP was in effect, we estimated that the program-induced increases in tree cover lead to annual negative CO 2 emissions (4.16 million metric tons). Our estimate is larger, but on the same magnitude, as the estimate of carbon negative emissions due to the PSNP found by (Woolf et al., 2018) using different methods. We also note that our estimate is equivalent to 1.5% of the reduction pledged by Ethiopia in the Paris Agreement (Federal Democratic Republic of Ethiopia, 2021). 3

A.7.2. Estimating cost
We estimated the annual administrative costs of the PSNP to be approximately 302 million USD (2007 dollars) using data from Drechsler et al. (2017), which, in turn, draw on data from the PSNP Interim Financial reports, annual reports and World Bank analyses. Program costs are given in current dollars. We deflated them to 2007 using the US CPI, because the social cost of carbon (SCC) is given in 2007 dollars, see below. Given our estimate that the program induced 4.10 (95% CI: 0.058, 9.28) million metric tons of negative CO 2 emissions annually, we thus calculated that the per unit cost to reduce a ton of CO 2 was USD 72.68 (95% CI: 39.98, 4,318).

A.7.3. Estimating benefits
To estimate the benefits due to negative carbon emissions, we used estimates of the SCC from the Interagency Working Group (Interagency Working Group on Social Cost of Greenhouse Gases, 2016). Specifically, we used the SCC from the year 2015 (close to the midpoint of our treatment period) of 36 USD (2007 dollars) per metric ton of CO 2 , which is the median estimate of the report, corresponding to a 3% discount rate. Since we did not have data on how long the increased tree cover will persist, we calculated benefits for four scenarios: assuming that trees are cut after fifteen years, 30 years, 50 years, or never.
To calculate the value of reducing a million metric tonnes of CO 2 emissions, we followed Jayachandran et al. (2017) and Ferraro and Simorangkir (2020), and used the following formula:

Table A15
Sensitivity analyses: restricting the data to pixels with no quality flags. (1) (2) (3) Area: All <300 ppl/km 2 <150 ppl/km 2 Panel A: IPTW model (Fig. 2 Note: The unit of observation is a pixel observed periodically. The outcome variable is mean percent of tree cover in each period. In Panel A, the equation reported in the main text is estimated. In Panel B, the same equation is re-estimated after omitting all pixels with quality flags. The standard errors are reported in parentheses and are clustered at the district level. The models used in both panels include a binary variable capturing pixels belonging to PSNP districts, region specific period fixed effects and mean annual rainfall over the period. *** p<0.01, ** p<0.05, * p<0.10. In this equation, SCC is the social cost of carbon and r is the effective discount rate, calculated by combining the time discount rate (δ) and the growth rate at which the SCC rises over this time period (g), by the formula in r = (1 + δ)/(1 + g) − 1. Following the estimates from the Interagency Working Group on Social Cost of Greenhouse Gases (2016), we used the median time discount rate of 3%, and an average SCC growth rate of 2.2%, resulting in an effective discount rate of 0.78%. S is the length of storage, which captures the period between deforestation and carbon emission; we assume S to be equal to zero in all our scenarios. Lastly D is the program-induced delay in carbon emissions, measured in years.
We calculated the benefit-cost ratio for four different scenarios: 1. Scenario 1: D = Infinity; S = 0. Under this scenario, we assumed that the increased tree cover due to the PSNP remains permanently in place (D = infinity) and CO 2 is released as soon as the trees are cut (S = 0). Using the formula above, we found that the benefit per metric tonne of negative CO 2 emissions = SSC = USD $36. 2. Scenario 2: D = 50 years; S = 0. Under this scenario, we assumed that the increased tree cover from the PSNP remains in place for 50 years, at which time the trees are cut down and all carbon is immediately released into the atmosphere. Using the formula above, we found that the benefit per metric tonne of negative CO 2 emissions = USD $11.62. 3. Scenario 3: D = 30 years; S = 0. Under this scenario, we assumed that the increased tree cover from the PSNP remains in place for 30 years, at which time the trees are cut down and all carbon is immediately released into the atmosphere. Using the formula above, we found that the benefit per metric tonne of negative CO 2 emissions = USD $7.51. 4. Scenario 4: D = 15 years; S = 0. Under this scenario, we assumed that the increased tree cover from the PSNP remains in place for 10 years, at which time the trees are cut down and all carbon is immediately released into the atmosphere. Using the formula above, we found that the benefit per metric tonne of negative CO 2 emissions = USD $3.97. Table 1 in the main text reports the benefit-cost ratios for the four different scenarios analyzed.