2.1. VM0007: A REDD+ methodology for emission reduction evaluation

VM0007 (Verra, 2020) is one of the major methodologies that defines the calculation of emission reductions in a REDD+ project. A key concept of VM0007 is the Reference Region for projecting the rate of deforestation (RRD), which is equivalent to a control unit in causal inference. The RRD is chosen to match the project area (PA) as closely as possible in terms of the following variables: deforestation drivers, landscape factors (e.g. forest types, soil types, elevation, etc.), and socioeconomic variables (access to infrastructures, policies, and regulations, etc.). For the 10–12 years preceding the implementation of a project, deforestation rates are aggregated over the RRD and projected as a baseline for the crediting period. The projection method can be a simple historical average or a predefined linear/nonlinear model, where the former is often used. After the projection, spatial modeling is applied for certain types of forest to adjust the baseline.

Several studies have reported that baselines set under VM0007 were overestimated because these baselines failed to consider a counterfactual scenario correctly or to eliminate the effect of external factors, such as policy changes (West et al., Reference West, Börner, Sills and Kontoleon2020, Reference West, Wunder, Sills, Börner, Rifai, Neidermeier, Frey and Kontoleon2023). As an example, Figure 1 describes the spatial distribution and the deforestation trends of the PA and RRD of the Valparaiso project (Eaton et al., Reference Eaton, Dickson, McFarland, Freitas and Lopes2014), where the project started in 2011. The deforestation rates are calculated using MapBiomas Brazil (Souza et al., Reference Souza, Shimbo, Rosa, Parente, Alencar, Rudorff, Hasenack, Matsumoto, Ferreira, Souza-Filho, de Oliveira, Rocha, Fonseca, Marques, Diniz, Costa, Monteiro, Rosa, Vélez-Martin, Weber, Lenti, Paternost, Pareyn, Siqueira, Viera, Neto, Saraiva, Sales, Salgado, Vasconcelos, Galano, Mesquita and Azevedo2020), which will be used in the case study in Section 4. The baseline reported in the Project Design Document (PDD) is also included in the chart. We can see that the baseline (the red dashed line in Figure 1b) might have been overestimated considering the actual deforestation level in the PA. In addition, we can also see that the baseline might have failed to capture the change of trends where the deforestation rates in the RRD kept decreasing from 2002 to 2012 but changed to an upward trend again. The main cause of these changes may be the implementation and change of Brazil’s National Climate Change Plan (Ferrante and Fearnside, Reference Ferrante and Fearnside2019; West et al., Reference West, Börner and Fearnside2019, Reference West, Börner, Sills and Kontoleon2020). The current methodology cannot effectively reflect such impacts to the baseline.

Figure 1. The Valparaiso project.

2.2. Related work

The SCM (Abadie et al., Reference Abadie, Diamond and Hainmueller2010) is one of the most popular methods for causal inference with time-series data. This method is designed for a case with one treatment unit and multiple control units, which is suited for the evaluation of a REDD+ project. Given that the RRD consists of multiple subunits (hereinafter “control units”), SCM finds an optimal weight to match both preintervention deforestation trends and covariates of the synthetic control (i.e., the weighted average of control units) to those of the PA. Note that a baseline estimated by SCM can reflect effects of external factors occurred after intervention while it cannot forecast the baseline, because it is calculated based on observations after intervention. CausalImpact (Brodersen et al., Reference Brodersen, Gallusser, Koehler, Remy and Scott2015) is another popular method based on Bayesian state-space models. It has several ideas in common with SCM, but there are several differences. In contrast to SCM, CausalImpact can forecast a baseline with a little modification of the original model since it is based on state-space modeling, while it does not include a covariate balancing mechanism between a treatment and synthetic controls.

Sparse estimation needs to be considered to deal with small sample problems when applying time-series causal inference models to REDD+ program evaluations. In these cases, the sample size is the number of time points of deforestation data before a project starts. Since deforestation data are often updated annually and available from the 1980s at the earliest, the sample size is usually around 30–40. On the other hand, a considerable number of control units can be necessary to get a good synthetic control, which may lead to a larger number of parameters than the sample size. Several studies consider the use of penalization methods in SCMs in the frequentist perspective (Doudchenko and Imbens, Reference Doudchenko and Imbens2016; Abadie, Reference Abadie2021; Ben-Michael et al., Reference Ben-Michael, Feller and Rothstein2021). In CausalImpact, the spike-and-slab prior (George and McCulloch, Reference George and McCulloch1993) is employed to select influential samples on the baseline. Although the spike-and-slab prior is a natural solution to variable/sample selection problems in Bayesian modeling, it leads to substantial computational costs. The Horseshoe prior is a continuous analogue of the spike-and-slab prior and is computationally more tractable (Carvalho et al., Reference Carvalho, Polson and Scott2010). If some prior knowledge is available about the degree of sparsity, it can be incorporated by specifying a prior on the global shrinkage parameter (Piironen and Vehtari, Reference Piironen and Vehtari2017).

In the context of general carbon crediting schemes including REDD+, the issue of junk carbon credit or overcrediting has been discussed for a long time (e.g., Haya, Reference Haya2009; Aldy and Stavins, Reference Aldy and Stavins2012; Badgley et al., Reference Badgley, Freeman, Hamman, Haya, Trugman, Anderegg and Cullenward2022), and one of the sources has been identified to be baseline setting (Bento et al., Reference Bento, Kanbur and Leard2016; Haya et al., Reference Haya, Cullenward, Strong, Grubert, Heilmayr, Sivas and Wara2020). SCM and CausalImpact have been applied to several studies evaluating REDD+ projects (Roopsind et al., Reference Roopsind, Sohngen and Brandt2019; Correa et al., Reference Correa, Cisneros, Börner, Pfaff, Costa and Rajão2020; Ellis et al., Reference Ellis, Sierra-Huelsz, Ceballos, Binnqüist and Cerdán2020; West et al., Reference West, Börner, Sills and Kontoleon2020, Reference West, Wunder, Sills, Börner, Rifai, Neidermeier, Frey and Kontoleon2023) and other forest conservation activities (Sills et al., Reference Sills, Herrera, Kirkpatrick, Brandão, Dickson, Hall, Pattanayak, Shoch, Vedoveto, Young and Pfaff2015; Rana and Sills, Reference Rana and Sills2018; Simmons et al., Reference Simmons, Marcos-Martinez, Law, Bryan and Wilson2018) in considering a counterfactual scenario. Many of the studies using causal inference methods have reported some degree of overcrediting. To the best of our knowledge, there is no study in the literature concerning ex-ante prediction and ex-post evaluation at once.

3.1. Bayesian state-space SCM

We propose a Bayesian state-space SCM, leveraging the time-series modeling in CausalImpact and the covariate balancing in SCM. Consider a REDD+ project that has one PA and whose RRD is divided as multiple control units. We model the preintervention deforestation rates of the PA and the control units by the following state-space model with a local linear trend:

(1) $$ {y}_t={\tilde{z}}_t^{\prime}\hskip0.5em \beta +{\varepsilon}_t,\hskip1em {\varepsilon}_t\sim N\left(0,\hskip0.5em ,{\sigma}_y^2\right), $$

(2) $$ {z}_t={\tilde{z}}_t+{e}_t,\hskip1em {e}_t\sim N\left(0,\hskip0.5em ,{\sigma}_z^2I\right), $$

(3) $$ {\tilde{z}}_{t+1}={\tilde{z}}_t+{v}_t+{\eta}_t,\hskip1em {\eta}_t\sim N\left(0,\hskip0.5em ,{\sigma}_{\tilde{z}}I\right), $$

(4) $$ {v}_{t+1}={v}_t+{\xi}_t,\hskip1em {\xi}_t\sim N\left(0,\hskip0.5em ,{\sigma}_v^2I\right), $$

$$ t=1,\dots, {T}_0, $$

where $ {y}_t $ is an observed deforestation rate for the PA at time $ t $ , $ {z}_t={\left({z}_{1,t},\dots, {z}_{J,t}\right)}^{\prime } $ is a vector of observed deforestation rates of the control units at time $ t $ , $ J $ is the number of the control units, $ \beta $ is a weight to be applied to the control units defined on the simplex,

(5) $$ \beta \in {S}^{J-1}=\left\{\beta |\hskip0.5em \sum \limits_{j=1}^J{\beta}_j=1,\hskip0.5em ,{\beta}_j\ge 0,\hskip0.5em ,j=1,\dots, J\right\}, $$

$ {\tilde{z}}_t={\left({\tilde{z}}_{1,t},\dots, {\tilde{z}}_{J,t}\right)}^{\prime } $ is a latent state vector, and $ {v}_t={\left({v}_{1,t},\dots, {v}_{J,t}\right)}^{\prime } $ is a vector of slope components of the local linear trend. Time $ t=1 $ denotes the year when observation started and $ t={T}_0 $ the previous year before the intervention. Equations (1)–(2) link the observed data, $ {y}_t $ and $ {z}_t $ , to the latent state, $ {\tilde{z}}_t $ : it assumes that for the control units the deforestation rates are observed as the addition of the latent state and noise, while for the PA, the deforestation rate is written as the addition of the weighted sum of the latent states of the control units and noise (Brodersen et al., Reference Brodersen, Gallusser, Koehler, Remy and Scott2015). The latter relates this model to SCM except for covariate balancing. Equations (3)–(4) define the temporal evolution of the latent state by the simple local linear trend model, which enables us to forecast $ {z}_t $ , and thus the baseline.

We consider a shrinkage prior on $ \beta $ to select important control units for more interpretable and stable results. Sparse estimation is needed because the sample size $ {T}_0 $ tends to be small due to data availability, while we need to take a satisfactory number of the control units $ J $ into account so that there exists some weight that well approximates the PA. Here we introduce the regularized Horseshoe prior proposed by Piironen and Vehtari (Reference Piironen and Vehtari2017): for $ j=1,\dots, J $ ,

(6) $$ {\displaystyle \begin{array}{c}p\left({\beta}_j|{\lambda}_j,\tau, c\right)\propto \hskip0.5em \frac{1}{2{\tilde{\lambda}}_j\tau}\exp \left(-\frac{\beta_j^2}{2{\tilde{\lambda}}_j^2{\tau}^2}\right),\\ {}{\lambda}_j\sim {C}^{+}\left(0,1\right),\hskip1em {\tilde{\lambda}}_j^2=\frac{c^2{\lambda}_j^2}{c^2+{\tau}^2{\lambda}_j^2},\hskip1em {c}^2\sim \mathrm{Inv}-\mathrm{Gamma}\left(1/2,1/2\right),\\ {}\tau \mid {\sigma}_y\sim {C}^{+}\left(0,{\tau}_0\right),\hskip1em {\tau}_0=\frac{m_0}{J-{m}_0}\frac{\sigma_y}{\sqrt{T_0}},\end{array}} $$

where $ {C}^{+}\left(0,1\right) $ is the standard half-Cauchy distribution, $ \tau $ is a global shrinkage parameter, $ {\lambda}_j $ is a local shrinkage parameter, and $ {m}_0 $ is the effective number of nonzero coefficients in $ \beta $ , which is given as a hyperparameter from prior knowledge. The global shrinkage parameter shrinks all coefficients toward zero, while the local shrinkage parameter controls the shrinkage level of each coefficient. Although the global shrinkage parameter has large impacts on the results, the systematic way of specifying a prior for $ \tau $ based on the information about the sparsity has not been well discussed. Piironen and Vehtari (Reference Piironen and Vehtari2017) provide the systematic way of specifying a prior for $ \tau $ based on prior knowledge of the sparsity. We can control the number of control units that show significantly large positive weights through the effective sample size $ {m}_0 $ .

For covariate balancing, we rely on the method of general Bayesian updating (Bissiri et al., Reference Bissiri, Holme and Walker2016). With this, we reflect the distance of the covariates between the PA and synthetic controls as a covariate-dependent prior on $ \beta $ . Together with (6), we set the prior on $ \beta $ as

(7) $$ p\left(\beta |{\left\{{x}_j\right\}}_{j=1}^{J+1},\hskip0.5em ,\lambda, \tau, c\right)\propto \exp \left(- wL\left(\beta; {\left\{{x}_j\right\}}_{j=1}^{J+1}\right)\right)\prod \limits_{j=1}^Jp\left({\beta}_j|{\lambda}_j,\tau, c\right), $$

where $ {x}_j $ is a $ K\times 1 $ vector of the covariates with $ j=J+1 $ for the PA, $ L $ is a loss function that measures a distance between the PA and synthetic controls, and $ w $ is a tuning parameter. Here we choose $ L $ to be a SCM-like loss function for covariate balancing:

(8) $$ L\left(\beta; {\left\{{x}_j\right\}}_{j=1}^{J+1}\right)=1/2\cdot {\left({x}_{J+1}-{X}_0^{\prime}\beta \right)}^{\prime }V\left({x}_{J+1}-{X}_0^{\prime}\beta \right), $$

where $ {X}_0={\left({x}_1,\dots, {x}_J\right)}^{\prime } $ and $ V $ is the inverse of the covariance matrix of $ {X}_0 $ .

Combining equations (1)–(8), we obtain the full posterior as

(9) $$ {\displaystyle \begin{array}{l}\hskip1em p\left(\beta, {\left\{{u}_t\right\}}_{t=1}^{T_0},\Sigma, \lambda, \tau, c|{\left\{{z}_t\right\}}_{t=1}^{T_0},{\left\{{y}_t\right\}}_{t=1}^{T_0},{\left\{{x}_j\right\}}_{j=1}^{J+1},w\right)\\ {}\propto \prod \limits_{t=1}^{T_0}f\left({y}_t,{z}_t,{u}_t|{u}_{t-1},\beta, \Sigma \right)\cdot p\left(\beta |{\left\{{x}_j\right\}}_{j=1}^{J+1},\lambda, \tau, c\right)p\left(\lambda, \tau, c\right)p\left({u}_0\right)p\left(\Sigma \right),\end{array}} $$

where $ f $ is the density function of the model (1)–(4), $ {u}_t={\left(\tilde{z}{'}_t,{v}_t^{\prime}\right)}^{\prime } $ , and $ \Sigma ={\left({\sigma}_y,{\sigma}_z,{\sigma}_{\tilde{z}},{\sigma}_v\right)}^{\prime } $ . After the observation at $ t={T}_1\left(\ge {T}_0\right) $ , we can obtain the posterior predictive distribution of the baseline up to the target period $ t={T}_2\left(\ge {T}_1\right) $ as

(10) $$ {\displaystyle \begin{array}{l}p\left({\left\{{y}_t^{\mathrm{bsl}}\right\}}_{t={T}_0+1}^{T_2}|{\left\{{z}_t\right\}}_{t=1}^{T_1},{\left\{{y}_t\right\}}_{t=1}^{T_1},{\left\{{x}_j\right\}}_{j=1}^{J+1},w\right)=\int \prod \limits_{t={T}_0+1}^{T_2}f\left({y}_t^{\mathrm{bsl}},{z}_t,\hskip0.5em |{u}_t,{u}_{t-1},\hskip0.5em ,\beta, \Sigma \right)\\ {}\cdot \hskip0.5em p\left(\beta, {\left\{{u}_t\right\}}_{t=1}^{T_0},\Sigma, \lambda, \tau, c|{\left\{{z}_t\right\}}_{t=1}^{T_0},{\left\{{y}_t\right\}}_{t=1}^{T_0},{\left\{{x}_j\right\}}_{j=1}^{J+1},w\right)\cdot d\beta \cdot d\Sigma \cdot d\lambda \cdot d\tau \cdot d c\cdot \prod \limits_{t=1}^{T_2}{du}_t\cdot \prod \limits_{t={T}_1+1}^{T_2}{dz}_t,\end{array}} $$

where $ {\left\{{y}_t^{\mathrm{bsl}}\right\}}_{t={T}_0+1}^{T_2} $ is the estimated baseline of the PA from $ t={T}_0+1 $ to $ t={T}_2 $ . As the project proceeds, the baseline can be updated in a unified manner for ex-ante baseline ( $ {T}_1<t\le {T}_2 $ ) and for ex-post baseline ( $ {T}_0<t\le {T}_1 $ ). Note that the estimation of $ \beta $ is based on the data up to $ t={T}_0 $ , because $ \beta $ represents the relation between the PA and the control units without intervention.

3.2. Estimation algorithm

In this section, we discuss an efficient Markov chain Monte Carlo (MCMC) algorithm for sampling the parameters from (9). Although we fit our model using Stan (Stan Development Team, 2016) for the simplicity of implementation in the empirical example below (Section 4), it would be worth discussing the structures of the full conditional distributions and the relevant literature. By the construction of the model, the full conditional distribution of each parameter reduces to a standard distribution, which enables us to obtain a Gibbs sampler easily, except for the weight $ \beta $ and the parameters regarding the regularized Horseshoe prior. Therefore, we focus on the sampling from the full conditional distribution of these parameters.

The model (1) is a Gaussian sparse linear regression model, and if the prior is the basic (without-regularization) Horseshoe prior the full conditional distribution of $ \beta $ becomes also a Gaussian distribution and the Gibbs sampler can be derived by using data augmentation techniques (Carvalho et al., Reference Carvalho, Polson and Scott2010; Bhadra et al., Reference Bhadra, Datta, Polson and Willard2019). However, our formulation has two problems: (1) since we consider the regularized Horseshoe prior, it is hard to find a closed form expression for the full conditional distribution of the shrinkage parameters, $ {\lambda}_j\hskip0.1em \left(j=1,\dots, J\right) $ and $ \tau $ and (2) since the weight $ \beta $ is defined on the simplex space, the full conditional is not a simple Gaussian.

For the first problem, we need to rely on a general sampling algorithm that can be applied densities without a closed form, e.g., Hamiltonian Monte Carlo (HMC) and Langevin Monte Carlo methods (Neal, Reference Neal2011). For the second problem, we can refer to the literature of the sampling from truncated Gaussian distributions. It is known that sampling from a Gaussian distribution defined on the simplex $ {S}^{J-1} $ reduces to sampling from certain truncated Gaussian distribution on $ {\mathrm{\mathbb{R}}}^{J-1} $ (Altmann et al., Reference Altmann, McLaughlin and Dobigeon2014). In our formulation, the covariate balancing term in (7) can be incorporated in this framework because the choice of the quadratic loss function in (8) makes the full conditional distribution have also a Gaussian kernel:

(11) $$ p\left(\beta |\cdot \right)\hskip0.5em \propto \hskip0.5em \exp \left(-\frac{1}{2{\sigma}_y^2}{\left(\beta -{A}^{-1}b\right)}^{\prime }A\left(\beta -{A}^{-1}b\right)\right),\hskip1em \beta \in {S}^{J-1}, $$

where

$$ {\displaystyle \begin{array}{c}A={\tilde{Z}}_0{\tilde{Z}}_0^{\prime }+w{\sigma}_y^2{X}_0V{X}_0^{\prime }+D,\hskip1em D=\operatorname{diag}\left({\tilde{\lambda}}_1^{-2},\dots, {\tilde{\lambda}}_p^{-2}\right)/{\tau}^2,\hskip1em \\ {}b={\tilde{Z}}_0^{\prime }y+w{\sigma}_y^2{x}_{J+1}^{\prime }V{X}_0^{\prime },\hskip1em {\tilde{Z}}_0={\left({\tilde{z}}_1,\dots, {\tilde{z}}_{T_0}\right)}^{\prime },\hskip1em y={\left({y}_1,\dots, {y}_{T_0}\right)}^{\prime }.\end{array}} $$

Equation (11) indicates that the full conditional distribution of $ \beta $ is the Gaussian $ N\left({A}^{-1}b,{\sigma}_y^2{A}^{-1}\right) $ constrained by the linear constraint $ \beta \in {S}^{J-1} $ , which reduces to a truncated normal distribution on $ {\mathrm{\mathbb{R}}}^{J-1} $ .

To sample from a truncated Gaussian distribution, we may use the exact Hamiltonian Monte Carlo (EHMC) method (Pakman and Paninski, Reference Pakman and Paninski2014). A traditional approach to this problem is the combination of the Gibbs sampling from a Gaussian and the rejection of samples that are outside of the domain, but it is known to be inefficient, particularly in a case of high dimension. On the other hand, Pakman and Paninski (Reference Pakman and Paninski2014) propose an efficient HMC algorithm for sampling from a multivariate Gaussian that is defined on a constrained space by linear and quadratic inequalities (or products thereof). In this special case, their algorithm achieves a hundred percent acceptance rate, while HMC methods in general need to reject a proposed sample with some small probability.

In Stan, we implement the following variable transformation so that $ \beta $ is contained in the simplex $ {S}^{J-1} $ :

$$ {\tilde{\beta}}_j={\alpha}_j\cdot {\tilde{\lambda}}_j\cdot \tau, \hskip1em {\beta}_j={\tilde{\beta}}_j/\sum \limits_{j^{\prime }=1}^J{\tilde{\beta}}_{j^{\prime }},\hskip1em j=1,\dots, J, $$

where $ {\alpha}_j $ is defined as a non-negative scalar. In this formulation, $ \alpha $ is sampled based on the full posterior (9) and then transformed into $ \beta $ . We use a truncated normal distribution as the prior for $ \alpha $ .