2.1. The contingent valuation method and payment for environmental services

CV uses questionnaires to elicit people's WTP for a good or service, creating a hypothetical market in which people can declare their preferences for the good. There have been over a thousand applications of CV studies in diverse areas of economics; the main results are summarized in a number of books covering theoretical and empirical issues (Mitchell and Carson, Reference Mitchell and Carson1989; Bateman et al., Reference Bateman, Willis and Arrow2001; Champ et al., Reference Champ, Boyle and Brown2003). Additionally, certain guidelines on best CV practices can be found in Arrow and Solow (Reference Arrow and Solow1993) and in Whittington and Pagiola (Reference Whittington and Pagiola2012). Among these practices, the suggested elicitation method is the closed-ended format for the WTP question in which respondents encounter a randomly assigned amount of money $A as a price for the environmental service (ES), and they decide whether they are willing to pay this amount.

Despite the pervasive presence of CV studies in several areas of economics, their application has been highly controversial. A discussion of these controversies is beyond the scope of this paper; however, an appropriate starting point for the interested reader is presented in the Journal of Economic Perspective’s CV symposium (Carson, Reference Carson2012; Hausman, Reference Hausman2012).

In the context of a PES, CV is used to estimate the maximum amount of money that the beneficiaries of a project are willing to pay to implement a program that ensures the provision of an ES. Certain valuation studies of ES are those by Zhongmin et al. (Reference Zhongmin, Guodong, Zhiqiang, Zhiyong and Loomis2003), Jin et al. (Reference Jin, Wang and Liu2008) and Loureiro and Ojea (Reference Loureiro and Ojea2008), and examples that deal specifically with PES include studies by Ortega-Pacheco et al. (Reference Ortega-Pacheco, Lupi and Kaplowitz2009), Van Hecken et al. (Reference Van Hecken, Bastiaensen and Vásquez2012) and Moreno-Sanchez et al. (Reference Moreno-Sanchez, Maldonado, Wunder and Borda-Almanza2012).

2.2. Estimating a binary quantile regression

Given that the elicitation question is the closed-ended format in which respondents encounter an amount of money $A as a price for the ES and they decide whether they are willing to pay this amount, we will solely have a positive (yes) or a negative (no) answer. We denote a positive answer of individual i as y _i  = 1 and a negative answer as y _i  = 0. As Haab and McConnell (Reference Haab and McConnell2002) described, a CV model assumes that the satisfaction that a consumer perceives can be represented by a utility function, denoted by u _j , which has a deterministic component $v_{j}(p,I,q_{j})$ and a random component $\varepsilon_{j}$ . That is, $u_{j} = v_{j}(p,I,q_{j})+ \varepsilon_{j}$ , where j = 0 denotes the initial situation, j = 1 denotes the new situation, p is a vector of current prices, I is income and q _j represents the environmental quality (we have dropped the individual's subscript i). A respondent will be willing to pay the amount A only if the utility of paying for this project is higher than the utility of the status quo in which he or she does not pay for the project (Hanemann and Kanninen, Reference Hanemann, Kanninen, Bateman and Willis1999) (u ₁>u ₀). Therefore, the probability of a positive (yes) answer is:

(1) $$Pr (yes) = Pr [\Delta v \gt \eta] = F_{\eta} (\Delta v)$$

with $\Delta v = v_{1}(p,I - A,q_{1}) - v_0 (p,I,q_0), \eta = \varepsilon_0 - \varepsilon_1$ and F _η as the distribution function of η.

We used a Bayesian approach based on an asymmetric Laplace distribution (ALD) to estimate a binary quantile regression using a parametric linear utility function given by $v_{j}=\alpha_{j}+\beta I + \varepsilon_{j}$ (for other functional forms, refer to Hanemann, Reference Hanemann1984). It can be shown that for this utility function Δv = α−βA, β>0, α = (α₁−α₀). Other explanatory variables enter into the model through coefficient α. The linear utility function is the most popular functional form used in the stated preference models due to the simplicity of estimating both the random utility model and WTP (Louviere et al., Reference Louviere, Hensher and Swait2000). The Bayesian approach used in this paper does not have convergence problems; it is simple to estimate, valid for small sample sizes, robust to heteroskedasticity and guarantees a global optimum in comparison to classic approaches, such as the Maximum Score (Manski, Reference Manski1985) and the Smoothed Maximum Score (Kordas, Reference Kordas2006; Florios and Skouras, Reference Florios and Skouras2008). Furthermore, the estimators are consistent, efficient, asymptotically normal and admissible (Cameron and Trivedi, Reference Cameron and Trivedi2005; Rossi and McCulloch, Reference Rossi and McCulloch2010).

As Hanemann (Reference Hanemann1984) noted, one needs an assumption regarding the distribution of the difference in the error terms of two indirect utility functions $(\eta = \varepsilon_{1} - \varepsilon_0)$ . In the conventional RUM approach, it is assumed that η has either a logistic or a normal distribution. For the Bayesian estimation, it is convenient to assume that η has an ALD (Hewson and Yu, Reference Hewson and Yu2008; Benoit and Van den Poel, Reference Benoit and Van den Poel2012) because the likelihood function associated with the ALD is directly related to the minimization problem associated with a quantile regression. Consider that a quantile regression will solve the following problem (Benoit and Van den Poel, Reference Benoit and Van den Poel2012).

(2) $$\tilde{\beta} = \arg \min \sum\limits_{i=1}^N \rho_{\tau} (y_i - x_i^{\prime} \beta),$$

in which $\rho_{\tau} (z) = z \lpar \tau -I(z \lt 0) \rpar $ , I(z<0) is an indicator function and τ∈(0, 1). This function assigns the weight τ to positive values of z and (1−τ) to negative values of z. Furthermore, a random variable y _i with ALD has a density function given by

(3) $$f(y\vert \theta,\sigma, \tau)= {\tau (1 - \tau) \over \sigma} \exp \left(-\rho_{\tau} \left({y - \theta \over \sigma} \right) \right),$$

in which τ is the asymmetry coefficient, $-\infty \lt \theta \lt \infty $ (we could define $\theta =x_{i}^{\prime} \beta $ ) is the location parameter, and σ>0 is the scale parameter. The main argument of this distribution (Yu and Zhang, Reference Yu and Zhang2005) is the same argument as equation (2).

Assuming that η has an ALD is convenient but not crucial. Li et al. (Reference Li, Xi and Lin2010) showed that the Bayesian quantile regression is robust to ‘functional form’ assumptions; in fact, many Monte Carlo simulations assume that the original error term is normally distributed (Li et al., Reference Li, Xi and Lin2010; Benoit and Van den Poel, Reference Benoit and Van den Poel2012) and use the likelihood function of a ALD for the estimation. As Yu and Moyeed (Reference Yu and Moyeed2001) noted, it is not necessary to specify the distribution of the error term.

The Bayesian approach relies on a combination of a prior distribution for the coefficients and the likelihood function, which provides a posterior distribution for the coefficients. In accordance with Benoit and Van den Poel (Reference Benoit and Van den Poel2012), if we assume that the difference in an indirect utility function $\Delta v^{\ast} $ is a linear function of the error term η, which is distributed as a Laplace distribution, that is, $\eta \sim {ALD}({0,1,\tau}),$ then the latent variables $\Delta v^{\ast} $ have an ALD of the form $\Delta v^{\ast} \sim {ALD}({\theta = x_{i}^{T} \beta,\sigma =1,\tau})$ , where $\theta_{\tau} =x_{i}^{T} \beta =\alpha_{\tau} -\beta_{\tau} A_{t} $ is the location parameter and σ is the scale parameter. The τ parameter should be specified at the quantile of interest, with τ∈(0, 1).

Given that in a binary CV people's responses are codified as y _i  = 1 for a positive response, that is $\Delta v_{i}^{\ast} \gt 0$ , and y _i  = 0 in other cases ( $\Delta v_{i}^{\ast} \lt 0$ ), then,

$$P\lpar y_i = 1 \rpar = 1 - F_{y_i^{\ast}} \lpar - x_i^{\prime} \beta \rpar ,$$

where $F_{y_{i}^{\ast}} (.)$ is the cumulative distribution function of the asymmetric Laplace latent variable $\Delta v_{i}^{\ast} $ . Therefore, the joint posterior density for β and $\Delta v_{i}^{\ast} $ , given people's responses, is $y = (y_{1},\ldots,y_{n})$ and the quantile τ is:

$$\eqalign{&\pi ({\beta,\Delta v^\ast \vert y,x,\tau})\propto \pi (\beta)\prod_{i=1}^n \{ {I({\Delta v_i^{\ast} \gt 0})I({y_i =1})+I({\Delta v_i^{\ast} \le 0})\,I({y_i =0})} \}\cr &\quad \times F_{y_i^{\ast}} ({\Delta v_i^{\ast}; x_i^{\prime} \beta, 1, \tau}),}$$

where I(.) is an indicator function, and π (β) is a prior distribution. As in Benoit and Van den Poel (Reference Benoit and Van den Poel2012), we assume that the prior of the π (β) parameters follows a normal distribution; that is, $\pi(\beta)\sim N(\bar{\beta}_{0}, V_{0})$ , where $\bar{\beta}_{0}$ is the vector of means of the k×1 order prior and V ₀ is the matrix of variances and covariance of β of k×k order.

Although we cannot sample from this unknown posterior distribution, we could use a MCMC algorithm. Benoit and Van den Poel (Reference Benoit and Van den Poel2012) suggested a Metropolis–Hasting within Gibbs algorithm that splits the posterior distribution into a conditional distribution of β given $\Delta v^{\ast}$ and a distribution of $\Delta v^{\ast}$ given β. They also show that the conditional distribution of $\Delta v^{\ast}$ given is a three-parameter ALD truncated from the left at zero when y = 1 and a three-parameter ALD truncated from the right at zero if y = 0; that is, $\pi (\Delta v^{\ast} \vert \beta, y, \tau)$ . Furthermore, this distribution is known, and one can sample from this distribution in a similar manner to sampling from a normal distribution.Footnote ¹

To complete the simulation, we need the conditional posterior distribution of β given $\Delta v^{\ast}$ represented by

$$\pi \lpar {\beta,\vert \Delta v^\ast,y,\tau} \rpar \propto \pi \lpar \beta \rpar \mathop \prod \limits_{i=1}^n \,F_{y_i^{\ast}} \lpar \Delta v_i^{\ast}; \,x_i^{\prime} \beta, 1, \tau \rpar .$$

The joint posterior function is generated from a block of conditional equations between $\Delta v^{\ast}$ and β; that is, given the data, the prior distribution and the quantile of interest, the join posterior distribution is obtained by sequentially drawing values from the truncated $\hbox{ALD}(x_{i}^{\prime} \beta,1,\tau)$ that characterizes $\pi \lpar \Delta v^{\ast} \vert \beta,y,\tau \rpar $ and the conditional posterior distribution $\pi \lpar \beta, \vert \Delta v^{\ast}, y, \tau \rpar $ . Benoit and Van den Poel (Reference Benoit and Van den Poel2012) incorporated this Gibbs sampler to the bayesQR statistical package of the R program.

2.3. Estimating the willingness to pay

After estimating a linear functional form, we can obtain the mean (and median) of the WTP distribution (known as compensation variation and denoted generally by C) as $\hbox{C} = \hbox{E(WTP)} = \alpha /\beta$ (Haab and McConnell, Reference Haab and McConnell2002). However, our purpose is to estimate the mean at each quantiles regression; therefore, we add a subscript τ, C _τ, to denote the corresponding quantile and, using the same logic suggested by Hanemann (Reference Hanemann1984), we estimate the mean as:

$$C_{\tau}^{+} = {\alpha_\tau \over \beta_\tau}.$$

Notice that we can estimate the ‘median’ or other percentiles of the distribution of each quantile regression using:

(4) $$\eqalign{F_\eta [\Delta v_{\tau} (C_\tau^{\ast})] &= F_\eta [\alpha_\tau - \beta_{\tau} C_\tau^{\ast}] = p, \quad \hbox{where} p\in (0,1). \cr F_\eta^{-1} [p] &= \alpha_\tau - \beta_{\tau} C_{\tau}^{\ast}, \cr C_\tau^{\ast} &= {\alpha_{\tau} - F_{\eta}^{-1} [p] \over \beta_{\tau}}.}$$

p is simply the probability that a person would be willing to pay the amount A in that portion of distribution of the dependent variable. In this case, F _η is the cumulative density function of ALD, α_τ and β_τ are the estimated parameters, and $C_{\tau}^{\ast}$ is the quantile welfare measure that satisfies the equality in equation (4). Furthermore, $F_{\eta}^{-1}[p]$ is the inverse ALD with $F_{\eta}^{-1} [p=0.5]=0$ , $F_{\eta}^{-1} [p=0.6]=-0.455$ and $F_{\eta}^{-1}[p=0.9] = -5.87$ .

It is feasible to estimate the quantile of the distribution of the mean WTP provided by the conventional RUM logit (or probit) model, but this is analogous to estimating the predicted quantile in equation (4) for an arbitrary τ using the α and β instead of α_τ and β_τ. A main difference between our approach and a traditional RUM approach is that we estimate α_τ and β_τ for each τ, whereas a traditional approach only has one overall estimate for α and β. This approach allows us to have a better description of the distribution and the credible intervals at each quantile.

Whether they differ significantly is an empirical question, but even if they do not differ significantly the Bayesian quantile regression will provide a distribution of the WTP at each quantile, not only a point. Furthermore, the Bayesian quantile approach allows us to capture the heterogeneous effect of the explanatory variables across quantiles. For instance, individuals with a higher probability of accepting the payment may also be less affected by the increase in the BID.

3.1. Survey design and implementation

Applying the CV survey includes determining the target population and the sample, designing the questionnaire, validating the survey, selecting and training interviewers and, finally, applying the final version of the survey. The relevant population for our survey included the 252,136 households in the urban area of the city of Santa Cruz that constitute the potential demand for the PES scheme. The sample was of 501 observations and has an error of 4.4 per cent. Our sample was concentrated in middle and lower income levels due to difficulties interviewing people in the high-level income bracket (see table A1 in the appendix). We used a probabilistic polietapic sampling on two levels. First, we randomly selected the neighborhood and blocks, and then we systematically selected the households to be interviewed. This approach means that we selected one household on each block starting in the northern corner, and if we did not have an answer at that house, we skipped the next four houses and knocked on the door of the fifth house.

Empirical evidence has shown that the selection of the bid vector and the size of the subsamples in each value may affect the estimation of the mean WTP (Cooper, Reference Cooper1993). Therefore, we followed a sequential procedure to define the bid vector. We combined the minimization of the mean WTP variance with an observation of the empirical distribution of the WTP after a portion of the sample was collected. From pilot surveys of 100 observations, we obtained values for the WTP using an open-ended question; this information was used to define the bid vector and the size of the subsamples for each bid using the method suggested by Cooper (Reference Cooper1993).