Carmelo J León

Juan M Hernández

Miguel León-Santana

THE EFFECTS OF WATER TEMPERATURE IN AQUACULTURE MANAGEMENT

Keywords: Expected returns, Farm aquaculture, Water temperature, Canary Islands, Seabream Expected returns, Farm aquaculture, Water temperature, Canary Islands, Seabream JEL Classification: C61, L11, Q12

published or not. The documents may come

Introduction

The cultivation of marine species is expected to be an increasingly fast growing industry in the next decades as techniques for off-shore production develop. New species are being introduced and there is need for appropriate models that are useful for evaluating management decisions. In the last decades, several models have been developed for the optimal management of the most relevant control variables in the farm.

The early model by [START_REF] Bjorndal | Optimal Harvesting of farmed Fish[END_REF] analyzed how the optimal harvesting time schedule is influenced by the price of fish and the most important costs involved in production, i.e. the feed and harvesting costs. [START_REF] Arnason | Optimal Feeding Schedules and Harvesting Time in Aquaculture[END_REF] incorporated the feeding schedule in the optimization problem, and considered various functional forms for the fish growth function. [START_REF] Heaps | The Optimal Feeding of a Farmed Fish[END_REF][START_REF] Heaps | Density Dependent Growth and the Culling of Farmed Fish[END_REF] included the possibility of previous culling as dependent on density, while [START_REF] Mistiaen | Optimal Feeding and Harvest Time for Fish Weight-Dependent Prices[END_REF] assumed a step-wise influence of fish's size on price. The general conclusion of these results is that an increase in the interest rate would result in shorter optimal harvesting time and weight. These models are based on fish growth functions which abstract from the influence of water temperature on fish growth.

In this paper we propose a model for aquaculture management that incorporates a fish growth function influenced by water temperature. This model is utilized to derive the optimal decisions on the feeding rate and the harvesting time, which have relevant economic implications in terms of expected returns for the farm.

In general, previous economic models in aquaculture have made basic assumptions about water temperature. The only work which has explicitly considered water temperature in the growth function is [START_REF] Cacho | Optimal control of fish growth[END_REF]. This model 3 introduces water temperature in an optimization framework for catfish farming, but the analysis is limited to the feeding rate and the costs of feeding.

Water temperature is a major determinant of fish growth and can not be controlled in off-shore farming. Rather, it is determined by environmental and geographical conditions. In addition, the variability of water temperature adds a risk component to farming decisions, which makes theoretical predictions of standard models to depart from empirical findings [START_REF] Pascoe | Optimal harvesting strategies: practice versus theory[END_REF].

Our results indicate that the input costs have a negative impact on the optimal harvesting size if a particular relation between the economic and biological factors is given. However, this is not the case for the output transportation costs. Since output transportation costs are incurred at the end of the culture cycle they have an ambiguous impact on the optimal harvesting size. That is, a high output transportation cost could be overcome by higher revenues to be obtained with larger harvesting sizes.

The next sections present a model for a fish farm, with the fish weight as the state variable and the ration size as control, and with harvesting time and harvesting weight optimally determined. An empirical case study is presented for two scenarios of seabream culture in off-shore cages. The first represents a farm in Mediterranean waters, close to population demand centers in Europe, and the second is placed in the Canary Islands, an archipelago in the Atlantic Ocean located 2.000 kms. from the European coast but with favorable water temperature for seabream cultivation. The simulation results indicate that product differentiation by the optimal size could be an efficient strategy to overcome the economic disadvantages of the more distant farms.

The Model

Let us assume the producer has the objective to maximize profits, defined by the difference between revenues and costs. In aquaculture management, the production process is based on the growth process of the organisms that determines the amount of biomass generated in a period of time. Let us consider a growth model based on fish physiology proposed by [START_REF] Brett | Environmental Factors and Growth[END_REF] and applied by [START_REF] Muller-Feuga | Modélisation de la croissance des poissons en élevage[END_REF], which has been shown to perform appropriately to represent the growth process of a large number of fish species. The model equations and their specifications are presented in Table 1.

The growth function is ( ) ( ( ), ( ), ( )), w t G w t r t t θ = & (1)
where w is the fish's weight, r is the ration size given to the fish, θ is water temperature and t stands for time. The ration size r(t) is defined in normalized terms, i.e., in relation to the maximum ration which would lead to satiation. r is equal to 0 when there is no feeding and approaches 1 at satiation level.

As in [START_REF] Cacho | Optimal control of fish growth[END_REF], water temperature is modeled by a cyclical function depending on the time of the year, , 365

) ( 2 sin ) ( ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + + = θ π θ θ t t A t ( 2
)
where θ is the average value of temperature during a year, A is the maximum level of temperature, and

θ θ t t t - = 0
, where 0 t is the single day of the year when fingerlings are stocked and θ t is the day of average temperature ( The influence of these factors on fish growth can assume a multiplicative form, i.e.

( ) ( ) ( )

1 2 3 ( ( ), ( ), ( )) ( ) ( ) ( ) , G w t r t t G w t G r t G t θ θ = Γ (3)
where functions G i (⋅), i=1,2,3, represent the effect of each variable on fish weight and Γ is a correction parameter.

The firm grows a given number of species N 0 =N(t=0) with a given initial weight w 0 until they reach a commercial size at t=T. However, not all fingerlings reach maturity because of some mortality. Thus, N(t)=N 0 e -Mt , where M is the mortality rate.

On the other hand, the price of the product p is assumed to be dependent on the size of the fish harvested, i.e. p(w(t)). Without lack of generality, p(⋅) is assumed increasing with weight ( 0 w p > ) because larger sizes have larger market prices 1 .

The total cost function includes the costs of feed, labor, transportation, fingerlings and other costs. Because the efficiency of ration varies across different sizes of fish and is affected by water temperature [START_REF] Brett | Environmental Factors and Growth[END_REF], the costs of feed (C f (⋅)) depend on the conversion ratio, i.e. the quantity of food which is necessary for the fish weight to grow by 1 gr. Therefore we assume the following expression:

( , , , ) ( ) ( ( ), ( ), ( )) ,

f f C N w r c N t f w t r t t w θ θ = & (4)
where c f is the per gram price of feed and f(w,r,θ) is the conversion rate for given levels of fish weight w(t), ration size r(t) and temperature θ(t). The amount of feed in a day is

given by multiplying the growth rate by the number of surviving individuals in the cages (N(t)). The fingerling costs c s are incurred at the beginning of the culture, while the output transportation costs C τ ( ) are seen at the end of the cycle and are given by Let us consider that the farm is interested in optimally managing both the ration size and the final weight of the fish which is going to be sold in the market, i.e. the size of the product. The problem involves the maximization of discounted profits, and can be formulated as follows:

( ) ( ), ) ( ) ( ), ( T w T N c T w T N C τ τ = (
F o r P e e r R e v i e w ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ), , ( ) 0 1 2 3 0 ( ) , , , . . ( ) ( ) ( ) , (0) , , ( ) , 2 (
) ( ) sin , 365 0 ( ) 1, ( ) 0, T hT hT ht Mt S f o r

t T w T e p w T w T c e w T c e c f w r w t e K c w dt s t w G w t G r t G t w w

T free w T free

t t t A r t w t Max τ θ θ θ π θ θ - - - -- - + + = = + ⎛ ⎞ = + ⎜ ⎟ ⎝ ⎠ ≤ ≤ ≥ ∫ & & ( 8
)
where i is the interest rate, h=i+M is the sum of the interest and mortality rates, and K=c L L/N 0 is the daily cost of labor per individual. Thus, the state variable of this non autonomous problem is represented by the fish weight w(t), which is controlled by the ration size r(t). Both harvesting time and size are freely determined. The solution will lead us to the optimal ration size r * (t) for each time period (i.e. per day) as well as to the optimal market size w * (T) for the product that the firm should deliver to the market.

This size is also grown in an optimal time span T * .

The current value Hamiltonian is:

1 2 3 1 2 3 1 2 ( ) ( ) ( ) ( , , ) ( ) ( ) ( ) ( ) (1 ), Mt f o H c G w G r G f w r e K c w t G w G r G r r θ θ ψ θ µ µ = -Γ - - + Γ + + - (9)
where ψ(t) is the costate variable and µ 1 (t), µ 2 (t) are non-negative multipliers with properties µ 1 (t)r(t)=0 and µ 2 (t)(1-r(t))=0 respectively. The specification of functions G i , i=1,2,3 is included in Table 1. In particular, G 2 (r)=Z(r)/Z(r c ), where Z(r) is the normalized growth function and r c is a constant ration of reference or culture ration. The conversion rate for different ration sizes f(w,r,θ) is defined by means of the theoretical ( )

1 3 1 1 2 1 ( ) ( ) ( , ) ( ) ( ) (1 ), ( ) Mt f o c H G w G c r f w t Z r e K c w r r Z r θ θ ψ µ µ = Γ - + - - + + - (10)
By applying the Pontryagin's maximum principle we have ( )

1 3 1 1 2 1 ( ) ( ) ( , ) ( ) ( ) 0. ( ) r f c H G w G c f w t Z r Z r θ θ ψ µ µ ′ = Γ -⋅ + + -= (11)
The following results are obtained from the conditions on interior and corner solutions,

1 1 1 1 0 ( ) 1 ( , ) ( ) 0, ( ) 0 ( , ) (0) 0, ( ) 1 ( , ) (1) 0 ( , ) 0. f f f f if r t c f w Z r if r t c f w Z if r t c f w Z c f w θ ψ θ ψ θ ψ θ ′ < < ⇒-⋅ + = ′ = ⇒ -⋅ + ≤ ′ = ⇒ -⋅ + ≥ ⇔ -⋅ ≥ (12)
The last condition is not verified since both c f and 1 f (w,θ) are positive functions. Thus, it is never optimal to feed at satiation 2 .

The first condition shows that the optimal ration size is determined such that the marginal benefits of weight (ψZ´(r)) have to be equal to the relative cost of the 100% ration size (c f ⋅ 1 f (w,θ)). If the marginal benefits are less than the relative price for any ration size then the optimal ration is null (r=0). The latter case is difficult to occur, since it would involve high feeding costs or high conversion rates for which commercial culture would not be viable. Applying the first order optimal conditions we obtain the following differential equations which govern the state and co-state variables:

( )

3 1 1 1 1 1 1 2 3 ( ) ( , ) ( ) ( , ) ( ) ( ) , ( ) ( ) ( ) ( ) ( ). f o w c f c h G f w G w f w G w G w r c Z r c Y r w G w G r G ψ ψ ψ θ θ θ θ ⎛ ⎞ Γ ′ ′ = + + - ⋅ + ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ = Γ & & (13)
On the other hand, if we impose the transversality conditions for w(T) and T (see Appendix 1 for details), we have

( ) ( ) ( ) ( ) , ( ) ( ) ( ) ( ) ( ) ( ( ), , ) ( ) ( ). w MT w f o T p wT p c p w T w T p c w T i M p c w T c f w T r w T e K c w T τ τ τ ψ θ = + - + - = + - + + + (14)
From the first equation in ( 14), it is found that in the last period T the marginal benefits In this condition were fulfilled, fish weight is inversely related to most of the input costs considered, as it is shown in Table 2. Therefore, the optimal weight decreases as the The main implication is that adverse cost conditions -i.e. higher wages or feed priceslead to earlier harvesting, with the result of lower optimal market sizes. The opposite would be the case for those conditions which result in an increase in prospective revenues for the farm. For the interest and mortality rate, the impact on harvesting weight is also negative for those values which satisfy condition (15).

Condition ( 15) is more likely to be satisfied as the growth rate of fish price is positive, i.e. p(w) is convex, or the labor cost per unity of biomass is very small. In this case, the influence of output transportation costs on the optimal size is positive if .

w w i M < + &
For growth rates of fish at the harvesting time larger than the interest rate plus the mortality index, a rise in transportation costs would reduce the optimal weight. Consequently, transportation costs show an ambiguous effect on the optimal management of fish weight. That is, the combination of the specific shape of the fish price schedule and the growth rate makes the influence of these costs uncertain.

With respect to water temperature, Appendix 1 proves that its impact on the harvesting weight at time T is always positive for values of the parameters in condition (15). This means that higher temperature would induce larger market sizes because of the advantage of a higher growth rate of fish weight. The implication is that favorable environmental conditions could compensate for adverse cost conditions in production and in the transportation of output.

Application

In this section we illustrate the model above using an empirical example in the Canary Islands, where water temperature has played an important role in the choice of the optimal harvesting weight. The particular expressions in Table 1 have been calibrated for seabream culture utilizing statistical techniques with empirical data (Appendix 2).

More details of the calibration can be found in [START_REF] Hernández | A growth model for gilthead seabream (Sparus aurata)[END_REF]. For the price schedule, we assume a logistic-type function, which reflects the large degree of continuity and decreasing rate of the empirical prices published for the European market 3 , from €3.15 per kgm. for sizes around 200 g. to €6.00 per kgm. for 700 g. (see Figure 1). The empirical conversion rate was statistically fitted based on results from fish physiology and growth data.

We consider two alternative scenarios which vary in terms of the costs of feed and the average water temperature. Table 3 presents the assumptions made with respect to the parameters of the model for both scenarios. Scenario A represents average values for the Mediterranean waters while scenario B represents average values for the Canary Islands waters, off the cost of Northwestern Africa, in the Atlantic Sea. The latter incurs in somewhat higher costs of feed because the input has to be imported from the mainland. These input transportation costs can be incorporated in parameter c f. It has also higher average water temperature and incurs in larger output transportation costs for access to the main consumer markets. However, local farms in the Canary Islands receive full compensatory subsidies for these costs (Gasca-Leyva, León and Hernández 2003), so we can assume identical parameters for the output transportation net costs for both Table 4 presents the optimal results for the harvesting time (T), final fish weight (w T ) and the conversion rate (f) for both scenarios respectively. The stocking is assumed in May 4 . The solutions were obtained with the Runge-Kutta numerical method with order five and one day step time. The algorithms were validated by finding very low local errors. The results for the conversion rate respond to the accumulated rate over all the production cycle 5 . The optimal weight and harvesting time are larger for scenario B than for scenario A.

These results did not vary significantly with the capacity of the farm. The conversion rate is also larger for the Canary Islands waters. Thus, producers in the Canary Islands should optimally choose a product differentiation strategy based on larger sizes which would compensate for its adverse cost conditions and would take advantage from higher growth rates which are possible because the higher water temperature. The Mediterranean scenario specializes in lower market sizes which are harvested at the beginning of winter, when the conversion rate reaches its highest levels, and feeding is more inefficient.

The optimal trajectories for the ration size under the assumptions of each scenario for farm capacity of 200 tm. are presented in Figure 2. The optimal ration size is always higher for the Canary Islands because of the advantageous water temperature. While optimal ration size in the Mediterranean waters varies between a maximum of 84% and a minimum of 64%, the range in scenario B is 92% when fish is stocked to 70% at the end of the culture. The pattern is declining as the fish grows and enters the colder approaches the summer peaks. This can be observed in scenario B after one year of culture and the second summer of culture has begun. Thus the optimal feeding schedule is critically influenced by the assumptions of the water temperature.

The differentiation strategy according to the optimal fish weight is not affected by changes in the relevant parameters of the model. Figures 3 and4 show respectively the impacts of the output transportation costs and the interest rate on the optimal fish weight. The theoretical model predicts that large labor cost per unity of biomass and a large price response to weight could lead to an ambiguous effect of transportation costs on market weight. For the numerical application, the output transportation costs have a negative impact on the harvest size, but it is very small and does not produce significant changes in the optimal choice across the relevant range. However, when the output transportation costs are very high the optimal weight tends to zero because the negative returns. These results are valid for both scenarios and match observed behavior, with producers in the Canary Islands selecting a larger size than in the Mediterranean waters.

The impact of the interest rate on the optimal size is also negative but very small. The lack of sensitivity found in the calibrated results is due to the particular assumptions made with respect to the set of parameters. Nevertheless, the interest rate does not have an impact on the differentiation strategy based on the choice of a larger fish size under favorable conditions of water temperature. The theoretical results show that most common input costs have a negative effect on the optimal harvesting weight if a particular condition over these costs and the other biological and economic factors is fulfilled. The particular condition that could lead to a positive impact of these parameters on harvesting weight would be a large rate of growth of the fish price with respect to weight and a large labor cost per unity of biomass. Nevertheless, the effect of output transportation costs can be still ambiguous.

Conclusions

Water temperature has a positive effect on the optimal weight, indicating that those farms with favorable environmental conditions could improve their performance by producing larger sizes. That is, the results suggest that product differentiation on the basis of fish's weight could be an optimal strategy for those distant scenarios that enjoy advantageous conditions of water temperature. Past research has shown that the location of firms could be a major determinant of technical efficiency (see for instance, Gumbau-Albert and Maudos (2002) for the Spanish industry). Thus, further research is needed on the impact of advantageous environmental conditions on the location of firms' investment.

The application of the model to two alternative scenarios in the Mediterranean and the Canary Islands indicate that water temperature is the main variable influencing the choice of the optimal size. The optimal size for each scenario was not significantly

Notes

1. As will be shown below this assumption is appropriate for the case of some aquaculture products, such as gilthead seabream. The assumption of a stepwise function for the price schedule, as in [START_REF] Mistiaen | Optimal Feeding and Harvest Time for Fish Weight-Dependent Prices[END_REF] can be considered a degenerated case which would limit the range of quantitative results obtained with our model. On the other hand, the stepwise assumption becomes somewhat intractable with more than two steps, and complicates unnecessarily the analysis. Nevertheless, the general assumption of a nonlinear continuous function reflects observed prices for some aquaculture species. From a market perspective, [START_REF] Bjorndal | A Dominant Firm Model of Price Determination in the US Fresh Salmon Market: 1985-88[END_REF] find out that the oligopoly pricing models perform well in estimating the price of farmed salmon in the US market, suggesting that price discrimination and product differentiation are viable strategies among firms.

2. However, feeding at satiation is a common practice (Azevedo et al. 1998, Glasser and[START_REF] Glasser | High stocking densities reduce Oreochromis niloticus yield: model building to aid the optimisation of production[END_REF].

3. http://www.globefish.org/marketreports/Groundfish/seabass-seabream.htm 4. The assumption about the stocking date affected the optimal harvesting weight in our model, but did not change the conclusions obtained.

( ) ( )
( ) ( ) ( ) ( ) ( ) ( ( ), , ) ( ) ( ).
2 ( ) 2 1 ( ) , 2
- + ∂ = < ∂ - -Feed cost (c f ): ( ) ( ) ( , , ) 0. f w w m f wr w w c L R γ θ θ + ∂ = < ∂ - & -Average water temperature ( ) θ : ( ) ( ) 1 ( ) ( ) ( ) 0. f f w w m cY r f c f w L R θ θ γ θ γ θ θ + + ∂ = > ∂ - -Labor cost per individual (K): 0. MT w w w e K L R ∂ = < ∂ - -Other costs (c o ): 0. o w w w w c L R ∂ = < ∂ - Appendix 2 (Please, insert

Growth model equations

Parameters description: m, growth function parameter, specific for each species, 0<m<1; θ M , maximum lethal temperature; α and β, thermal function parameters, specific for each species; r c , culture ration, recommended levels of feed suppliers; r m , maintenance ration; r o , optimal ration, where a minimum conversion rate is reached.

(*) While r> r m , and r m ,<0.5, generally accepted.

Description Expression

Fish growth function

1 2 3 ( ( ), ( ), ( )) ( ( )) ( ( )) ( ( )) G w t r t t G w t G r t G t θ θ = Γ Weight function 1 ( ) m G w w = Ration function 2 ( ) ( ) ( ) c Z r G r Z r = Thermal function ( ) ( ) ( ) 3 ( ) M M G D e e α θ θ β θ θ θ - - = - Normalized conversion rate ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - - - - - + = ) )( 2 3 )( 1 ( ) 1 ( ) ( 1 ) ( 2 2 2 m o m o m o o m r r r r r r r r r r r r Y Normalized growth rate (*) 0 , 0 , ) ( ) ( < ′ ′ > ′ = Z Z r Y r r Z Conversion rate 1 ( , , ) (

  function (C L ) depends on the number of workers L and the daily cost of labor, which is assumed to be constant c L . Thus, .

  c τ is the output transportation cost per unit of weight and s τ is the subsidy to transportation cost per unit of weight.The function of other costs (C o (⋅)) includes sanitary expenditure, maintenance, etc. which depend on the cumulated biomass, i.e. represents the daily costs.

  Y(r) and function f 1 (w,θ), which represents the conversion rate for a 100% ration size. Introducing these definitions and after some simplifications the Hamiltonian is transformed into

  of allowing the individual to grow by one gram (ψ(T)) are directly proportional to the increase in revenue to be obtained with this weight gain, i.e. (p+ p & w). The output transportation net costs reduce the potential benefits which would be obtained by an extra growth of the individual. The second equation in (14) gives the harvesting weight w(T) as dependent on some economic and biological factors in the model. The influence of these factors on the optimal harvesting weight is determined by the following condition (see Appendix 1),

  or other costs becomes larger. Similarly, if the wage rate increases then the weight for the individual is lower. On the other hand, if the farm's size -as measured by the number of fingerlings N 0 -rises, then labor cost per individual would be lower, which implies a larger optimal weight.

  rest of data have been obtained from previous market analysis and observations of real culture in Canary Islands.

  , in both scenarios the ration size should rise as the cycle

  focused the management of the key variables of the ration size and the harvesting weight in off-shore aquaculture, taking into account the impact of water temperature as a critical environmental condition. The model allows us to evaluate some of the potential economic trade-offs farms can face between adverse economic conditions and advantageous environmental factors.

  in the output transportation costs and the interest rate. The differentiation according to the market size was neither affected by changes in other parameters of the model.

  side and the right hand side in equation (A.5) respectively. For the sake of simplicity, we omit the argument of variable w. The influence of any factor

Figure 2 .Figure 3 .Figure 4 .

 234 Figure 1. Fish price function:

  Table A2 here)

	Page 21 of 30		Submitted Manuscript
		F o
		r
	3.5	P e
	3		e
	2.5		r
	0 2	200	400	e R 600 800	1000	1200
				v i
				e w
	Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK 21

Table 1

 1

Table 2

 2 Effect of changes in unitary costs over harvesting weight if condition (15) is fulfilled.

	Submitted Manuscript
	Sign (+) indicates positive influence; sign (-) indicates negative influence; 0 indicates
	no influence.
	Unitary cost	Effect over harvesting weight
	c f (feed price)		-
	c L (wage rate)		-
	c S (fingerling cost)		0
	c o (other costs)		-
	L (number of employees)		-
	N 0 (number of fingerlings)		+
	i (interest rate)		-
	M (mortality index)		-
	c τ (output transportation net costs)	+/-
	θ (Average water temperature)	+

Page 27 of 30 Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Page 30 of 30 Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript

& & &

Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UKSubmitted Manuscript

& & & &

5. The algorithm used in the numerical solution was implemented with MATLAB® and consisted on determining the trajectory r(t) for the state variables in (13) to satisfy the transversality conditions ( 14) with a tolerance of one thousandth. This is available from authors upon request.

. We deduce from here that

The second transversality condition is obtained by solving expression ( )

where H 1 is the Hamiltonian introduced in (9) and removing the last two terms. The calculations are


Alessandra Ferrari
email: a.ferrari@reading.ac.uk

Econometric and Statistical Methods: Special Topics < C -Mathematical and Quantitative Methods, I11 -Analysis of Health Care Markets < I1 -Health < I -Health, Education, and Welfare

Keywords:

The internal market and hospitals efficiency: a stochastic distance function approach.

Introduction

The UK was one of the first European countries to reform its health services 1 . As regards hospital services, a clear distinction was introduced between purchasers (District Health Authorities and GP fund-holders) and providers (hospital trusts).

By becoming trusts hospitals were given more autonomy in the management of their resources and had to sell their services to the purchasers via contracts on the so-called internal market [START_REF] Bartlett | The performance of trusts[END_REF]. Competition for contracts was thought to give hospitals incentives to efficient behaviour. Even though the "quasi-market" was formally eliminated in 1997, its main features remained in place. Other countries are now following in the same direction, and the debate about the efficiency and effectiveness of this kind of reform is a major policy issue (see for example [START_REF] Gerdtham | Health system effects in cost efficiency in the OECD countries[END_REF].

Not surprisingly the reform generated scientific interest, a lot of which of a theoretical nature, on the efficiency properties of different contracts or the credibility of competitive conditions on the market. The empirical work has been instead more restricted, especially when it comes to the analysis of efficiency changes [START_REF] Glennerster | Competition and Quality in Health Care: the UK Experience[END_REF]. Among the main contributions, [START_REF] Propper | Market Structure and Prices: The Responses of Hospitals in the UK NHS to Competition[END_REF] and [START_REF] Propper | Competition in the NHS Internal Market: an Overview of its effects on Hospital Prices and Costs[END_REF] studied the effects of competition on hospital prices and costs, and found some evidence of price reduction eventually mirrored in the costs, but mainly a great variability in prices that bore no respect to the average cost pricing rule advocated by the Government. The effects of competition on the quality of hospital services were the focus of [START_REF] Propper | Does Competition between Hospitals Improve the Quality of Care? Hospital Death Rated and the NHS Internal Market[END_REF] who found evidence of a decrease in quality associated to the working of competition. Studies of hospitals efficiency on more recent, post-reform data sets can be found for example in Jacobs (2001) and [START_REF] Street | Relative performance evaluation of the English acute hospital sector[END_REF], who 3 focused on the (lack of) significance of the differences in inefficiency between hospitals. [START_REF] Jacobs | Hospital Efficiency Targets[END_REF] showed the non-effectiveness of hospitals efficiency targets; these have become increasingly popular with policy makers, and their feasibility, as well as that of general efficiency analysis of public services, have been thoroughly analysed by Smith and Street (2003).

Not many contributions exist instead about the effects on hospitals efficiency during the years of the reform, when competition was actively promoted as the key feature of the system. Among these studies, Soderlund et al. (1997) estimated a classical linear regression model on a sample of NHS hospitals in England for the years 1992 to1994 which revealed a general productivity improvement however not associated to trust status. Due to its easier availability, others have used the acute Scottish hospitals data set used in this paper 2 . Scott and Parkin (1995) estimated a translog cost function for the years 1992/93 which highlighted the prevalence of constant returns to scale and economies of scope between different kinds of outputs. [START_REF] Maniadakis | The impact of the internal market on hospital efficiency, productivity and service quality[END_REF][START_REF] Maniadakis | Assessing productivity changes in UK hospitals reflecting technology and input prices[END_REF] used Data Envelopment Analysis (DEA) to calculate Malmquist indexes of total factor productivity (TFP) for the period 1991/92-1995/96, finding a general improvement in TFP mainly attributable to shifts of the frontier, and a worsening of the quality level 3 . The aim of this paper is to analyse the changes in technical efficiency and performance of hospitals during the years of the reform by means of a stochastic distance function. As opposed to DEA, the stochastic parametric approach allows the statistical testing of hypotheses, making the results more reliable 4 . It also allows analysing the characteristics of the production process that a nonparametric method by definition does not identify. Furthermore the chosen model [START_REF] Coelli | Efficiency measurement, multiple output technologies and distance functions: with application to European railways[END_REF] for technical efficiency, which is a frontier model as opposed to the classical linear regression one; as will be seen in Section 2, the frontier model specifically separates the noise in the data from the estimation of inefficiency, which is the aim of the exercise. The choice of a distance function form is due to the multiple output nature of the production process that rules out the direct estimation of a production function 5 .

Unfortunately, data for the whole of the NHS are not available; the analysis is performed on a sample of 52 acute hospitals in Scotland observed between 1991/92 and 1996/97, thus covering the whole duration of the reforms.

The rest of the paper is structured as follows. The estimation of stochastic frontiers and the distance function model are discussed in Section 2. Section 3 describes the data set. The model selection process and the results are in Section 4, and general conclusions in Section 5. A final appendix details the additional statistical analysis performed to check on the robustness of the model.

The stochastic frontier and the distance function.

Following Debreu and Farrell seminal papers, the efficiency of a firm can be defined and measured as the distance of its actual performance from a frontier 6 .

If this frontier is the production function, i.e. the maximum attainable output from a given set of inputs, the distance will measure technical inefficiency.

Following Shephard, this can be measured as a radial expansion D o for the output vector(s) in order to reach the frontier, where The efficiency measure u i must come from a positively skewed distribution; for instance if this is a half normal (which will be used later), then

0 < D o 1. If D o =1
u i = |U i | U i N(0, u 2 )
Due to the presence of a composite error term OLS gives consistent but inefficient slope parameters' estimates, and the use of ML is to be preferred if the distribution of u i is known or an assumption can be made about it.

Following Battese and Corra (1997) the influence of the inefficiency component can be measured by a parameter = 2 u / 2 , where 2 = 2 u + 2 . The significance of can be tested with an LR test which, if the null hypothesis H 0 : = 0 is true, follows a mixed 2 distribution. If the null hypothesis is true and inefficiency is not significant, the model is equivalent to a standard "average" production function, and its log-likelihood is the same as that of the linear model estimated by OLS.

As (1) can be estimated only for the single output case, an alternative model has been proposed by Coelli andPerelman (1996, 2000) to deal with the multiple

= + + + + + + = = = = = = = = = (2)
Using the fact that D o is linearly homogeneous in outputs, adding the stochastic component v i and setting lnD oi =-u i ,, (2) can be rearranged as

9 i i mi ki km M m K k li ki kl K l K k ki k K k ni mi mn M n M m mi m M m Mi u v y x x x
x y y y y

+ + + + + + + = = = = = = = = = * ln ln ln ln 2 1 ln * ln * ln 2 1 * ln ln 1 1 1 1 1 1 1 1 1 1 1 1 0 (3)
where y M is the M-th output and y* m =y m /y M . Equation (3) can be now estimated as a usual production frontier, by regressing (the log of) one output on (the logs of) the inputs and (the logs of) the outputs ratios. Note that the coefficients of a production frontier correspond to the negative of the coefficients of a distance function.

The data

The data are a sample of 52 acute hospitals in Scotland in the years 1991/92 to 1996/97 (from now on referred to as 1992 and 1997 respectively), making a concerns [START_REF] Scott | Investigating hospital efficiency in the new NHS: the role of the translog cost function[END_REF], revealed that the information on some of the cross sectional units was not reliable. Furthermore, as suggested also in [START_REF] Smith | Measuring the Efficiency of Public Services: the Limits of Analysis[END_REF], additional statistical analysis of the residuals has been performed to check for the reliability of the results and the robustness of the model.

Coming to the definition of the inputs and outputs variables, the following choices have been made 11 . Output is usually measured as the total number of cases treated. Since cases are very heterogeneous they are usually divided into various specialty (or casemix) categories, which qualify the hospital as a multiple-output unit. The use of index numbers can overcome the problem of the trade-off between heterogeneity and degrees of freedom as long as one can define weights that correctly represent the differences between cases. On the assumption that more difficult illnesses are more input demanding than the less serious ones, measures of their reference average cost are produced by the NHS and customarily used by the literature.

This paper follows this approach and the weights reflect the average costliness of a case in each category, calculated for the whole of Scotland in order not to bias them with some measure of inefficiency of each hospital. Furthermore, to avoid average cost changes over time due to possible inefficiency changes, they have been normalised each year to sum to 1. In detail, if q jit is the total number of cases treated in category j by hospital i at time t, and c jit is the average cost per case in category j at hospital i at time t; then is the average cost (across hospitals) per case in category j at time t The weight for each category is calculated as

= = J j jt jt jt c c w 1 .
. and so finally the index for each hospital in each year is

= = J j jit jt it q w y 1
The many output categories have then been summarised in two indexes: one for the inpatients and one for the outpatients, day patients and day cases. The main difference between the two categories is that in the former patients spend several days in the hospital and in the latter no more than one day, sometimes without even using a bed (and the staff associated to it). As a substitution between the two kinds of services could have taken place, it was preferred to keep them separated, so the two final output indexes are y 1 (index of inpatients) and y 2 (index of outpatients, day patients and day cases).

Finally, 5 variables identify the inputs:

x 1 = total capital charges (£000)

x 2 = medical staff FTE (full time equivalent);

x 3 = nursing staff FTE;

x 4 = other staff FTE;

x 5 = total number of beds.

The capital measure 12 was deflated using the "Hospital and Community Health Services pay and price inflation values". The "other staff" input includes 1.

[Table 1 here] 4.

The model

The main problem encountered in the estimation of the translog output distance function was the modelling of time effects. A series of LR tests was performed to choose among alternative time specifications; this initially led to the estimation of model with intercept dummy variables for each year in the panel, like the In the above, (the log of) the index of outpatients, day patients and day cases (y 2 ) is regressed over (the log of) the five inputs (x) and the outputs ratio (y* = y1/y2), plus five dummy variables to allow for a different intercept each year, and a dummy variable for teaching hospitals. A dummy variable for trust status could not be introduced because of the implicitly assumed correlation with the inefficiency component (the issue is discussed more in detail later). Following [START_REF] Battese | Frontier production functions, technical efficiency and panel data: with application to paddy farmers in India[END_REF], the inefficiency component was modelled as a function of time as

following 13 : it it teach t t t kit k k lit kit kl l k kit k k it u v D D y x x
u it = u i exp[-(t-T)]
(5) so that a value of > 0 (<0) implies increasing (decreasing) efficiency over time.

A value of = 0 implies no time effect, and the hypothesis can be tested by 4) however produced unsatisfactory results; in particular, even though the model with the intercept dummies had been statistically selected against other time specifications, not all the dummies were significant, their overall pattern was counterintuitive and contradicting existing evidence (see for example [START_REF] Maniadakis | The impact of the internal market on hospital efficiency, productivity and service quality[END_REF], 2000[START_REF] Ferrari | Market oriented reforms of health Services: a non-parametric analysis[END_REF] and finally one of the inputs elasticities was negative. All this raised the suspicion that the effect of time had not been adequately captured, namely that a pooled model might not be appropriate and that all the parameters of the distance function, not just the intercepts, might have changed over time.

As the use of Chow tests for parameters stability was ruled out for lack of degrees of freedom an alternative approach was used instead. This consists of estimating several times the distance function with a time interaction dummy instead of the intercept dummies. In particular, a time dummy d is introduced, which takes a value of 1 for a particular year(s), and 0 else, and this is multiplied to all the variables in the translog distance function, as: The inputs and outputs of (6) are the same as those of (4), D teach is a dummy variable for teaching hospitals and d is the time interaction dummy. When d=0 the parameters of the function are the s, s and s; when d=1 they are the respective ( s+ s), ( s+ s) and ( s+ s). The inefficiency component u it varies The interaction dummy is first set equal to 1 for 1992 (and 0 else), then for 1992 and 1993 (and 0 else) and so on. In this way 5 different distance functions are estimated, each with a different time effect which is captured by the parameters of the interaction dummy. The likelihood results of the five estimations of ( 6) are reported in Table 2.

it it teach kit k k lit kit kl l k kit k k kit k k lit kit kl l k kit k k it u v D d y x d x x
[Table 2 here] In each case the significance of the time interaction parameters is tested for by means of LR tests against a restricted model with no time effects and as expected the null hypothesis is always rejected. On the grounds of the Akaike 14 information criterion the model as specified in ( 6) is also to be preferred to the pooled one in (4).

Since the models are not nested in one another but they all have the same number of parameters, selection by minimisation of any standard information criterion is equivalent to selection of the model with the greatest maximised log likelihood. This happens when separating 1992 and 1993 from all other years, as the model has about 27 points of difference in the log-likelihood from its closest alternative. This therefore points to the fact that the parameters of the distance (and production) function might have changed after 1993.

The model in which 1992 and 1993 are separated from the following years is therefore analysed. The main results are shown in Table 3, which reports in order: the ML and OLS log-likelihood values, the estimated values of and , the inputs elasticities and then the elasticity of y 2 with respect to the outputs ratio 15 . All elasticities are calculated at the sample mean and their significance is tested by means of an LR test that in all cases leads to reject the null hypothesis (i.e. all are significant) 16 .

[Table 3 here]

The results of the estimation can be summarised as follows. The parameter is significant (LR test is 337.62) meaning that so is inefficiency. Very interestingly however is not: the result of the LR test (2.46) leads not to reject the null hypothesis that = 0, so that no significant difference in (in)efficiency appears to have taken place over time.

The most notable difference between the two time periods is the change in the elasticity of substitution between y 2 and y 1 : the absolute value increases by 60%, implying that the opportunity cost of treating someone as an inpatient over time has increased considerably. A pattern has emerged towards treating patients more and more as day patients/cases or directly as outpatients. This is confirmed when looking at Table 1, that shows a very big rate of increase in the value of y 2 as opposed to a relatively small increase in that of y 1 .

As regards the inputs, two variables improve their elasticity after 1994 (namely the nursing staff and the beds), and the other three lower it. Whether this is due to a change in the productivity of an input or just to a change in its levels can be revealed by testing for the significance of the relevant dummy variable parameters. This reveals that the increased elasticity of the beds input is a consequence of the reduction in the levels of the variable, well known also to the 13 general public via the news. This is also consistent with the aforementioned trend towards day-based treatment, as beds would in that case be used more intensively.

The reduced elasticities of capital and of other staff are taken to be a direct consequence of the reforms. The increase in capital levels could be due to the investment in information technology that hospitals made in order to deal with the new contracting issues [START_REF] Fattore | Cost containment and health care reforms in the British NHS[END_REF]. As this activity is not directly linked to the treatment of patients (the output variable) this might explain the reduced elasticity of the input. However, increased capital levels are also the consequence of accountancy changes related to the change to trust status, which made the hospitals owners of their assets, so concluding that there was overcapitalisation would be misleading. The data did not offer any other measure of capital but the one used, which has its limitations although its entire removal from the equation was rejected from the data; keeping this in mind, since no further detail is available, the result has to be taken with caution.

Similarly, the "other staff" variable increases in level and its elasticity decreases from 0.33 to 0.12. One reasonable explanation is that administrative staff increased to deal with the new contracting issues. Another possibility is that the lower increase in nursing staff might have led to the transfer of some of their duties to cheaper but less qualified (and therefore less efficient) staff. A pattern towards the use of cheaper labour inputs in Scottish hospitals was revealed by [START_REF] Gray | Factor input in NHS hospitals[END_REF], and this might have been reinforced by the financial concerns of the reforms. Finally, the results also show a positive difference in the intercepts indicating a possible improvement in average productivity, that is a shift upwards of the frontier. Considering that the shape of the frontier has changed, the higher intercept could indicate that the improvement is mainly in the production of the dependent variable, i.e. again y 2 .

The fact that 1994 is the first year in which some hospitals start to change to trust status naturally leads as to think that is the reason behind the structural break.

The relevance of trust status in explaining changes in technology and inefficiency is therefore analysed more in detail. One approach to do so is the estimation of a "one-step model" (Wang and Schmidt, 2002), which however failed to converge to a maximum. More sophisticated models are currently discussed by the literature and will be worth exploring for future research. For this paper the interaction dummy specification remains the preferred one. Using the estimates of ( 6), the elasticities of trusts and non trusts are calculated and compared to one another, as reported in Table 4 17 . This shows that the pattern revealed by Table 3 after 1994 seems to be more marked for the trusts sample than for the other hospitals, which confirms the hypothesis that the change in technology is related to the change in status. However, no significant link between trust status and efficiency can be detected: a t-test on the equality of the mean inefficiency score, computed as E[u i i ] 18 of trusts and non trusts is performed and the null hypothesis cannot be rejected (the p-value is 0.34).

[Table 4 here]

As mentioned in Section 3, to check on the reliability of the results and the robustness of the model specification, two additional pieces of analysis have . We find that there is very little overlap and thus that inefficiency is not only significant overall but also significantly different between hospitals.

Concluding remarks.

This paper estimated a stochastic distance function to analyse the changes in (technical) efficiency and performance of acute hospitals in Scotland during the years of the reform. A structural break was detected in 1994, the year of the first trust wave, after which hospitals change not only the way in which they provide their services, but also the mix of services they provide. The opportunity cost of inpatients increases as hospitals tend to treat patients more and more on a day basis. This view is supported by the fact that both the number and the costliness of outpatients, day patients and day cases increase quite significantly, indicating a possible "swap" between the two categories of output considered. This could be the result of an attempt of hospitals at reducing their costs by reducing the length of stay, especially if the contracts constrained them to provide minimum levels of treatment (as it was the case especially with the widely used 'block contracts'). The involvement in the new contracting activity, and the financial concerns that it brings with itself, appears also to translate into reduced inputs productivity. This at least seems to be true for the capital and other staff variables This increase in the day-basis treatment is also behind the shift upwards of the frontier. However, although technical inefficiency remains significant, it does not show any significant improvement.

There are clearly limitations to this work, among which the non-availability of a specific quality measure and the relatively small sample size. With these in mind, the general conclusion that can be drawn from this work is that the reform did not unambiguously produce the expected efficiency improvements. The move towards day-based care might be going hand in hand with the reduction in the quality of treatment found elsewhere in the literature, and this makes the overall scenario not particularly optimistic.

Analysis of the residuals of equation (4).

In this appendix we provide a graphical comparison of the theoretical distribution of the composed error term, with an estimate of the density of the actual residuals. The theoretical distribution, resulting from a convolution, may be found in [START_REF] Greene | Frontier production functions[END_REF]. For the density estimator we use a kernel method (see for example [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF], with the standard plug-in bandwidth. For the relevance of this kind of analysis see for example [START_REF] Smith | Measuring the Efficiency of Public Services: the Limits of Analysis[END_REF] and [START_REF] Street | Relative performance evaluation of the English acute hospital sector[END_REF] for a similar application to hospitals efficiency.

Figure A1 shows in order: the theoretical density of the composed error term i = v i -u i (with v i ~N(0, v 2 ), u i = |U i | and U i N(0, u 2 )), the estimated density of the residuals from the estimation of equation ( 4), and the distribution of the 312 residuals.

The residuals seem to be in broad agreement with their theoretical distribution and there are no serious outliers, indicating that the results are robust.

[Figure A1 here] 3 Measured as the survival rate 30 days after discharge. 4 See Jacobs (2001) for a discussion of the relative advantages of the two methodologies for the hospital sector, and Parkin and Hollingsworth for an application to the Scottish data set. 5 This problem arises when estimating a production function as opposed to a cost function. See for example [START_REF] Conrad | A multiple-output, multiple-input model of the hospital industry of North Carolina[END_REF][START_REF] Conrad | A multiple-output, multiple-input model of the hospital industry of North Carolina[END_REF]Wagstaff (1989). 6 For an overview of the many techniques of frontier estimation, including more recent developments, see for example [START_REF] Sena | The Frontier Approach to the Measurement of Productivity and Technical Efficiency[END_REF]. 7 Many comprehensive reviews of the stochastic frontier model exist in the literature; see for example [START_REF] Greene | Frontier production functions[END_REF] and [START_REF] Kumbhakar | Stochastic Frontier Analysis[END_REF]. 8 A more parsimonious specification like the Cobb-Douglas should not be used because, apart from other restrictions, it is not concave in the output dimensions. 9 The detailed explanation is omitted for reasons of space and can be found in Coelli andPerelman (1996, 2000). 10 The data are from the Scottish Health Service Costs statistics. 11 See in particular [START_REF] Mcguire | A production function analysis of acute hospitals[END_REF]; [START_REF] Tatchell | Measuring hospital output: a review of the case mix and service mix approaches[END_REF][START_REF] Mcguire | Methodological considerations of hospital production and cost functions: relationships to efficiency[END_REF]; Butler (1995).

12 This was the only available capital measure in the data, and it comprises: a) depreciation on fixed assets; b) interest paid on money borrowed to finance any of the projects in a); c) a 6% return on capital (trusts only). Even though not ideal its entire removal from the estimations was rejected by the data. 13 Since this first model was eventually discarded, for reasons of space the details are omitted, but they are available in a previous longer version of this work [START_REF] Ferrari | NHS reforms and efficiency of hospital services: a stochastic distance function approach[END_REF]). 14 The Akaike information criterion is used to compare models that are non nested in one another.

It is specified as AIC = -2 + 2n where is the value of the maximised log-likelihood and n is the number of parameters. The preferred model is the one with the lowest AIC value. 15 The details of the calculation of the elasticity of substitution are in Ferrari 2004.. 16 The presence of the squared and interaction terms makes the translog prone to multicollinearity. As a consequence it is usually advisable to test for joint parameters' significance rather than relying on their individual ones. 17 In particular these are calculated at the average sample values in 1994 and 1995, which are the years where a reasonable mix of trusts and non trusts exists. 1992 and 1997 in fact have 0 in one category, and 1993 and 1996 have 7 or less in one category. 18 See [START_REF] Greene | Frontier production functions[END_REF] for a discussion.

  idea is to directly express D o as a function of the K inputs and M outputs of each of the N firms. Using a log-linear translog function specification 8 this is:

  of 312 observations 10 . The data set was actually reduced from the originally available one after a detailed analysis, motivated also by others'

  , administrative, clerical and all other staff. The descriptive statistics of the data are reported in Table

  LR test. Equation (

  in (5) and is modelled as a half normal distribution, since the LR test on the 0-mean of a truncated normal distribution leads not to reject the null hypothesis in all cases.

  e y* ), the total input elasticity (or elasticity of scale) and finally the elasticity of substitution between y 2 and y 1

  are reported in detail in the appendix. First we compare the theoretical density function of it with that of the estimated residuals; this shows that the model fits the data well and rules out the presence of outliers. Secondly for each hospital we computed the 2.5%, 50% and 97.

  increase over time), and for the medical staff. Nurses and beds instead improve their productivity, and are also associated with the lowest increases in levels (with the latter strongly negative in fact).

  (1989) "Estimating efficiency in the hospital sector: a comparison of three statistical cost frontier models", Applied Economics, 21, pp.659-672.Wang, H and Schmidt, P (2002) "One-step and two-step estimation of the effects of exogenous variables on technical efficiency levels", Journal of ProductivityAnalysis, 18,[129][130][131][132][133][134][135][136][137][138][139][140][141][142][143][144]

Figure

  FigureA2for shows, for the 52 hospitals, the quantiles corresponding to

Table 1 : Descriptive statistics of the inputs and outputs variables (standard deviation into brackets) Average 92-97 Yearly rate of growth

 1

	Inpatients index	141746 (125124)	0.3
	Outpatients et al. index	30342 (28349)	9.0
	Capital (£000)	1513 (1396)	1.6
	Medical staff (WTE)	88 (87)	4.1
	Nursing staff (WTE)	457 (360)	0.7
	Other staff (WTE)	302 (256)	4.0
	Beds	357 (272)	-2.7

Table 2 : Log-likelihood of the translog distance function with time interaction dummy. 92 92 -93 92 -94 92 -95 92 -96

 2

	177.99 218.54 190.93 191.68 170.19

Table 3 : results from the estimation of equation (4).

 3

	Variable		Estimate
			218.54
	OLS		49.73
			0.955*
			0.02
		1992-1993	1994-1997
	e cap	0.028*	0.005*
	e med	1.049*	0.672*
	e nurs	0.156*	0.254*
	e oth	0.334*	0.121*
	e bed	0.043*	0.289*
	e y*	-0.493*	-0.613*
	e tot	1.611	1.341
	e y1	-0.973	-1.582

* = significant at 5% (or less) ** = significant 1 at 10%.

Table 4 : partial elasticities of trusts and non trusts hospitals.

 4

		Trusts	Non trusts
	e cap	-0.001	0.022
	e med	0.667	0.687
	e nurs	0.252	0.259
	e oth	0.121	0.128
	e bed	0.308	0.207
	e y1	-1.594	-1.506

comparison of the theoretical density function of i with that of the estimated residuals. Figure A2: quantiles corresponding to

  White Paper Working for Patients (December 1989), NHS and Community Care Act (June 1990, effective from April 1991). 2 ISD Scotland, Scottish Health Service Costs, NHS in Scotland.

	Page 25 of 26						Submitted Manuscript
	1
		1.5	Theor. density of
		1.0
		0.5
				-1.4	-1.2	-1.0	-0.8	-0.6		-0.4	-0.2	0.0	0.2	0.4
		1.5	Estimated density of residuals.
		1.0
		0.5
			0	-1.4 Residuals. F -1.2 o	-1.0	-0.8		-0.6	-0.4	-0.2	0.0	0.2	0.4
		-1		r
			0	25		50	75 P 100 e	125		150		175	200	225	250	275	300
								f r e	( u	i	i	1	,...	iT	)	, and
	to 1.00 f 1.25	(	oi D	) 97.5 % quantile ,... 1 iT i 2.5 % quantile	.	50 % quantile (median) R e v i
	0.25 0.50 0.75												e w
			0	5		10	15	20		25			30	35	40	45	50
	100			2.5 % quantile 97.5 % quantile		50 % quantile (median)
	75
	50
			0	5		10	15	20		25			30	35	40	45	50
													25

Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK

Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UKSubmitted Manuscript

Acknowledgements.

Very useful comments and help were received from Wiji Arulampalam, Simon Burke, Mike Pitt, and from David Spiegelhalter as discussant of the paper at the RSS 2003 one-day meeting on "Performance Monitoring and Surveillance".

Thanks go also for the advice of an anonymous referee. Usual disclaimers apply.

Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK

Submitted Manuscript


Jorge E Arana

Carmelo León

Modeling Unobserved Heterogeneity in Contingent Valuation of Health Risks

Keywords: Contingent Valuation, Flexible Distributions, Heterogeneity Contingent Valuation, Flexible distributions, Heterogeneity, Health preferences, Normal mixture, Willingness to pay

  L'archive ouverte pluridisciplinaire

Introduction

The economic valuation of health risks is relevant for adopting efficient health policies, and is increasingly approached with non-market valuation methods such as contingent valuation (CV). This is a survey approach involving a market transaction carried out on a questionnaire posed to the relevant population. Past studies by [START_REF] Viscusi | Pricing Environmental Health Risks: Survey Assessments of Risks-Risks and Risks-Dollar Trade-Offs for Chronic Bronchitis[END_REF] and [START_REF] Johannesson | The Value of Private Safety Versus the Value of Public Safety[END_REF] among others have shown that people are sensitive to the level of risk of illness.

However, the values people place on policies for reducing health risks could be subject to unobserved heterogeneity. This can be explained because of a number of related reasons: i) subjects can have different degrees of risk aversion, ii) they could perceive and/or be sensitive to health risks in different ways, and iii) respondents could present different preferences for alternative health statuses. These sources of variability in the data can rarely be observed by the researcher.

Thus, sample data on the values of health risks elicited following the CV method can be conditioned by unobserved heterogeneity, breaking down the assumption of the same distribution across the sample. This assumption is maintained in parametric and semiparametric models for the dichotomous choice elicitation method. This is a method which requests binary responses from the subject to a given price for the proposed policy involving the reduction in health risks. Most applications in health valuation have fitted a parametric distribution such as the probit model to these binary responses.

In this paper we consider a mixture distribution approach which allows us to represent heterogeneity across the sample in contingent valuation surveys utilizing the dichotomous 1 F o r P e e r R e v i e w choice format. The proposed approach has both theoretical and practical advantages. From the statistical point of view, it provides flexibility and with enough components virtually any empirical distribution can be approximated arbitrarily well (Zubov (1985)). In particular, multiples modes are possible and fatter tails than the normal can also be accommodated by mixing in normal components with large variance. From the economic point of view, unlike most of the other methods, an estimator of the link function may be strictly monotonic. The violation of this assumption implies that consumers are not necessarily utility maximizing (Coppejans, 1999).

There are other approaches which have been proposed for flexible modeling of this type of data. For instance, Horrowitz (1996) proposes a consistent and asymptotically normal nonparametric estimator for single index models with unknown transformation of the dependent variables. However, this estimator it is not necessarily monotonic and these results do not extend to the binary choice model. Another possible way to incorporate flexibility in the model is by using a exponential generalized beta type 2 [START_REF] Nagler | Scobit: an alternative estimator to logit and probit[END_REF]. Of particular interest, is the work of [START_REF] Mcdonald | An application and comparison of some flexible parametric and semi-parametric qualitative response models[END_REF]McDonald ( , 2000)), who extended this family of beta distributions to permit heterogeneity in the scale parameter. McDonalds´ results suggest that a random coefficient probit performs well compared with semiparametric alternatives.

In the context of modeling discrete choice contingent valuation (DCCV) data (see [START_REF] Hanemann | The Statistical Analysis of Discrete-Response CV Data[END_REF] for a detailed review). This is the case of the flexible semi-nonparametric model by [START_REF] Creel | Semi-Nonparametric Distribution-Free Dichotomous Choice Contingent Valuation[END_REF], which is similar to [START_REF] Chen | Semi-nonparametric Estimation of Binary Response Models with an Application to Natural Resource Valuation[END_REF]. These approaches allow the researcher to deal with the misspecification of rigid distribution functions imposed on irregular empirical data, but do not address the problem of unobserved heterogeneity across the sample observations. In addition, by using 2

F o r P e e r R e v i e

w

Monte Carlo simulations, Araña and León (forthcoming) show that the mixture of normals approach reduces bias and improves performance with respect to the Creel and Loomis semi- nonparametric model, particularly when the sample is characterized by heterogeneous preferences. Another nice feature about the mixture of normals is that it can be viewed as a natural extension of the paradigm case for purely parametric setting -the probit model.

The mixture model is applied to data obtained from the valuation of the reduction in the probability of an episode of respiratory illness. These data are likely to show heterogeneity because of the different degrees of risk aversion and different risk perceptions that people can have with respect to potential health outcomes. The model for sample heterogeneity based on a mixture of normal distributions is shown to improve performance with respect to a rigid probit model, which is also rejected according to Horowitz and Hardle's (1994) test.

The outline of the paper is as follows. Section 2 presents the modeling approach to unobserved heterogeneity and the estimation algorithm using Bayesian Markov chain Monte Carlo. Section 3 discusses an extension of the Horowitz and Hardle´s (1994) test for probit specifictation. Section 4 presents the data obtained from a CV survey on the reduction of the risk of an episode of respiratory illness, while section 5 discusses the results. Finally, Section 6 summarizes the main conclusions.

Model of sample heterogeneity

The dichotomous choice format requests from the subject to answer a yes/no willingness to pay question. The subject is offered a bid price (B j ) for a reduction in the level of health risk 3 F o r P e e r R e v i e w from z 0 to z 1 . Let e i (.) be the expenditure function for individual i, that is, the inverse of the indirect utility function with respect to income. If the subject seeks the largest level of utility then the answers would be yes if and no otherwise, where

V ) , ( ) , ( * 0 * 1 V z e B V z e i j i

≤ +

* is some fixed level of utility, and we assume other arguments of the expenditure function, such as prices, as constants.

The expenditure difference could be seen as the individual's willingness to pay (WTP) for the reduction in health risk, that is,

WTP i = . ( 1
) ) , ( ) , ( * 0 * 1 V z e V z e i i -
Therefore, the observed answer {y i } to bid price {B i }takes the value one if WTP i is higher or lower than the bid price, and zero otherwise. We can assume that the latent variable WTP is a function of two components, a deterministic µ and a random component ξ. In general, we can write , where µ is the mean of WTP and

) , ( WTP WTP i i i i ξ µ = i
ξ is a random error term, which is assumed normally distributed, with zero mean and σ standard error. Assuming independent answers and fixed covariates, the probability of a positive answer to bid price

B i is [ ] [ ] ) , B ( WTP B ) , ( WTP ob Pr ) 1 ( ob Pr i 1 - i i i i i i i F y µ ξ µ ξ ξ = > = =
(2) where is the linear predictor associated with a k x 1 regression parameter vector β and a covariate vector x β µ The error distribution can be specified as some parametric form, and the model estimated by maximum likelihood [START_REF] Hanemann | The Statistical Analysis of Discrete-Response CV Data[END_REF]). However, if there are different groups of individuals across the sample this assumption cannot be maintained, since each group could have a specific distribution, and unimodality does not hold. Thus, the assumption of a common distribution can be overcome by considering that the stochastic terms adopts the following form:

' i i x = i ,
) ( ... ) ( ) ( 2 2 2 2 1 1 1 1 k k k ik i i i e e e σ η α σ η α σ η α ξ + + + + + + = (3)
where α j and σ j (∀j=1,2,…,k) represent the mean and the standard deviation of each of the normal forms in the mixture 1 , and and

k k k k R R + ∈ = ∈ = ) ,..., ( ; ) ,..., ( 1 ' 1 ' σ σ σ α α α X i η ∼ N(0,
= = = ij j e π .
Without further restrictions, the model is clearly unidentified in the sense that more than one set of values of the parameters imply the same p(y|X). In particular, to prevent interchanging the components of the mixture some labelling restrictions are needed. Here we impose i) α 1 < α 2 < ... < α m . The additional restrictions require are the followings 2 : ii) rank (X)=k and a′X′ ≠ (1,…,1) for any k x 1 vector a; (ii) p j > 0 ∀ j; (iii) the support of is a set of positive Lebesgue measure; (v) σ β ' i x j = 1 for some j.

1 The mean and the variance of the error terms are ; respectively, and where [START_REF] Geweke | Mixture of Normals Probit Models[END_REF] for specific details about the specification. The inclusion of a mixture of normals is also a way of introducing flexibility in the distribution assumption, since it deviates from the assumption of a single parametric distribution for the entire sample. Thus, the mixture model approaches a semi-parametric model, since it eliminates the normality assumption. In addition, [START_REF] Heckman | A Method for Minimizing the Impact of Distributional Assumptions in Econometric Models for Duration Data[END_REF] and [START_REF] Zubov | Standardization of the analytic representation of probability distributions and of stochastic matrices[END_REF] show that any distribution can be approximated arbitrarily well by increasing the number of mixtures k in the model.

[ ] ∑ = = k j j j i E 1 α π ξ [ ] ∑ = = k j j j i Var 1 2 2 σ π ξ ∑ = = k j j 1 1 π 2 See
The likelihood function across the sample under the mixture assumption is:

( ) ( ) [ ] ∏ ∑ = = + - - Φ = N i j i j i k j j i x B y X p y p 1 ' 1 , , , , σ β α π σ α β ( ) [ ] ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ - - Φ - - ∑ = + j i j i k j j i x B y σ β α π ' 1 1 ) 1 ( (4)
Maximum likelihood estimators (MLE) of the parameters can be obtained by accordingly optimizing this function using some non-linear method. However, sampling distributions of MLE cannot be derived in closed form so that estimator properties rely almost exclusively on asymptotic or large sample situations. In particular, any model requires extremely large (as much as 1000 observations per parameter) samples to insure the adequacy of asymptotic approximates [START_REF] Mcculloch | An Exact Likelihood Analysis of the Multinomial Probit Model[END_REF]). In addition, asymptotic properties of ML estimators do not need to be maintained with small and finite samples 3 . [START_REF] Anderson | Logistic Discrimination and Bias Correction in Maximun Likelihood Estimation[END_REF] and [START_REF] Griffiths | Small Sample Properties of Probit Models Estimators[END_REF] found out relevant biases with numerical simulations of probit and logit models with small samples, while [START_REF] Copas | Binary Regresión Models for Contaminated Data[END_REF] 3 See [START_REF] Amemiya | Advances Econometrics[END_REF] or [START_REF] Huber | Robust Statistics[END_REF] for a more detailed discusión of the asymptotic proporties of ML estimators.

6

Page utilize Taylor series expansions to define the bias obtained for a logit model with small samples.

In these cases it could be utlilized at Bayesian approach (see for example Araña and León (2002)). Bayesian methods, as developed by [START_REF] Chib | Bayes Inference in the Tobit Censored Regression Model[END_REF] and [START_REF] Albert | Bayesian Analysis of Binary and Polichotomous Response Data[END_REF], adhere to the likelihood principle and are conducted using formal rules of probability theory. This means that under mild conditions Bayes estimators are consistent, asymptotically efficient and admissible in small samples. As a practical matter, Bayesian inference is free from the use of asymptotic approximations and delivers exact, finite sample inference [START_REF] Rossi | Bayesian Statistics and Marketing[END_REF]). The prior distributions are defined as follows:

β ∼ ( ) 1 , - β β H N α ∼ ( ) 1 , - α α H N j j s σ 2 ∼ ( ) j v 2 χ ( ) [ ] ( ) ( ) ∏ = - - Γ = k j j j v j v j j v s s v p j j j 1 2 2 / ) 2 ( 2 / 2 2 / 5 . 0 exp 1 2 / 2 1 ) 1 ( σ σ σ ) (π p ∼ Beta(r)
where ,

( ) k r r r ,..., 1 = ) ,..., ( 1 k α α α =
and α H is an k x k matrix. The last two matrixes are positive defined. The posterior distribution is derived by using a Gibbs sampling algorithm, similar to [START_REF] Geweke | Mixture of Normals Probit Models[END_REF]. This involves sampling from the conditional posterior distributions when they are known in a feasible form. Thus, even though WTP = (WTP 1 , WTP 2 , ..., WTP n )' is not observed, it is possible its simulation from available information. Taking the starting value for θ i.e. θ (0) = (β (0) , α (0) , σ (0) , π (0) ), the Gibbs sampling algorithm obtains iterated samples from each of the posterior conditional distributions. The algorithm is carried on t times leading to the simulated vector (WTP (t) , ,

β (t) , α (t) , σ (t) , π (t) ) obtained from the joint distribution (WTP, β, α, σ, π)|Y.
β h (t) , α h (t) , σ h (t) , π h (t)
)] H h=1 . The moments of interest are obtained from these simulated values.

A formal test of the probit model

Before applying the flexible modeling approach of a mixture distribution it is convenient to test whether a rigid model such as probit is appropriate for the binary choice data obtained in a CV experiment. Horowitz and Hardle (1994) propose a specification test (H-H) for the distributional assumption of the error term of the probit model. In this section we adapt this test to the parameterization of a latent variable WTP model.

The test is based on the difference between the parametric fit of the model ( )

[ ][ ] ) ( ) ( ~ ) ( ) ( ' ' 1 ' ' β β β β i i i n i i i i n X X F X Y X u h T Φ - Φ - = ∑ = (5)
which is asymptotically normal with zero mean and variance (

1 ) 1 ( 16 15 ) ( 2 2 ≤ - = v I v v K ) (6)
where I is the indicator function. This estimator is asymptotically unbiased. Thus there is no trade-off between asymptotic bias and variance as is commonly the case with standard techniques of determining the bandwidth. Proença (1993) used Monte Carlo simulation to show that the test can be influenced by the choice of h with finite samples. These results 5 Appendix 1 presents the development of the statistics in more detail. 6 See for instance, [START_REF] Martins | Parametric and Semiparametric Estimation of Sample Selection Models: An Empirical Application to the Female Labour Force in Portugal[END_REF] for a demonstration. suggest that under-smoothing could leads to the rejection of H 0 . Thus, it is convenient to utilize large values of the bandwidth in order to improve the potency of the test.

Application

The data for the application of the mixture distribution model come from a contingent valuation study conducted in the island of Gran Canaria (Canary Islands) to measure the economic benefits from the reduction in the probability of an episode of respiratory illness caused by air pollution. The field work was carried out in 1999. The questionnaire was improved after two pre-tests studies and two focus groups where the key design aspects Heterogeneity in the sample could emerge because some people might not be interested in reducing the possibilities of a respiratory illness. This behavior can be explained because they could not feel as potentially affected by this health problem. That is, they could perceive that their risk of getting the health problem is low because of their current good health conditions or because they have already taken precautionary measures. In addition, the proposed solution might involved unpopular measures such as the restrictions in the use of private transportation that some people could prefer against the reduction in the probability of illness in the general population.

As a provision rule, the subject was asked to assume that the proposed measures were going to be undertaken only if the majority of the population agreed on them. The subject was also recalled about her budget constraint and was asked about her willingness to pay for the reduction in the probability of having an episode of respiratory illness. The suggested risk reduction levels were presented with the aid of a risk ladder and were defined as [START_REF] Copas | Binary Regresión Models for Contaminated Data[END_REF]40,60 and 70 percent reductions. The risk card was tested in focus groups and was found effective as a risk communication device. randomly ordered with the aim of avoiding order bias effect. The price vector was designed with the distribution of the responses to an open ended question in the pretest study and following [START_REF] Cooper | Optimal Bid Selection for Dichotomous Choice Contingent Valuation Surveys[END_REF] design method for a predetermined number of bids. The bid price vector was 6, 12, 18 and 24 Euros for the lowest risk level, which was proportionately raised for the higher levels of risk reduction.

Results

The H-H test allows us to determine whether there is need to utilize a more flexible modeling approach to the empirical data. This could be convenient in order to avoid misspecification errors as coming from the utilization of rigid parametric structures, such as probit or logit. The use of these models could result in biased estimates of the structural parameters defining the mean and median willingness to pay for the reduction of the risks of respiratory illness.

The H-H test statistics was computed for different values of h, in order to compare the sensitivity of the results to the choice of this parameter. The variance was estimated by (A.5) and the weighting function u(v) was assumed to take the unitary value inside the confidence interval and zero in other case (Proenςa, 1993). The limits of the intervals were taken as the 10 th and 90 th percentiles, although the results were not sensitive to changes in these limits. there is need to consider a more flexible approach to model the empirical data. The left tail of the estimated distribution is always fatter for any value of h than the non-parametric approach. Thus, the probit model particularly fails to represent those subjects on the tails of the distribution. A more flexible approach might be able to model more precisely these extreme responses to the WTP question, establishing different classes of individual positions with respect to the value of the health prevention policy.

Since we have no prior information on model parameters for the Bayesian models, we assume very non-informative diffuse priors with a large variance. The starting values for Gibbs sampling were taken from ML estimation, although the results were quite robust to changes in these parameters. The posterior results were generated by running the Markov chains for a burn-in period of 10000 draws and the retaining every 10 th draw of the next 150000 draws. Convergence checks as in [START_REF] Raftery | How many iterations in the Gibbs sampler? in Bayesian Statistics[END_REF] and [START_REF] Geweke | Evaluating the Accuracy of Sampling-Based Approaches to the Calculation of Posterior Moments" en Bayesian Statistics[END_REF] did not indicate any problems, and this was corroborated by visual inspection of trace plots and rerunning the samplers many times from different starting values.

The estimated WTP model for the reduction of the risk of respiratory illness includes some covariates which allowed us to improve the predictive power and the goodness of fit of the model. Table 2 presents the description of the covariates which were significant in explaining WTP, obtained after probing with various potential specifications. The results of the flexible Bayesian model for heterogeneous samples are presented in Table 3.

For comparison purposes we also present the results of the Bayesian probit model (BPM) estimation following [START_REF] Albert | Bayesian Analysis of Binary and Polichotomous Response Data[END_REF] approach applied by Araña and León (2002) in DCCV, which disregards modeling heterogeneity. There are three main reasons for this 13 F o r P e e r R e v i e w decision. Firstly, simulation studies show that under non informative priors the Bayesian probit model (BPM) produces asymptotically equivalent results that the random coefficients probit model (RCPM) [START_REF] Geweke | Mixture of Normals Probit Models[END_REF]. The random coefficients models have been shown to be an interesting alternative to semiparametric and nonparametric approaches to incorporate heterogeneity [START_REF] Mcdonald | An application and comparison of some flexible parametric and semi-parametric qualitative response models[END_REF]. Secondly, the asymptotic properties of the RCPM estimated by maximum likelihood are not guaranteed with small and finite samples, while the BPM was shown to improve the results with respect to maximum likelihood estimation, especially for small samples. Finally, since the BPM is the most simply version of the mixture of normals model, the comparison of its fit with the more complex mixture models can be used as an additional test of the relevance of the unobserved heterogeneity in a specific application [START_REF] Bianchi | Testing for Convergence: Evidence from Non-Parametric Multimodality Tests[END_REF].

WTP is significantly influenced by income, as expected by theory. The effect of income is quadratic, indicating that WTP raises with income but eventually at a decreasing rate. On the other hand, WTP is also explained by the number of days of the potential episodes of respiratory illness and by the level of risk reduction. Both variables have positive effects on WTP. That is, those subjects posed with larger number of days of the episodes and larger risk reduction levels responded with higher latent WTP values. In addition, WTP becomes smaller for those subjects that have some type of chronic disease. This result could seem surprising, but can be explained because the marginal benefit of a health risk reduction measure could be smaller for the already ill. The results with the BMNP model improve substantially over the basic model. It can be seen that the t-ratio statistics become larger when the latent classes are modeled. In order to compare both models we calculate the Bayes' factor. This is defined by the ratio of the predictive probabilities of the alternative models, and can be interpreted as the times that one model is more likely than an alternative model.

The likelihood of alternative models can be evaluated with the predictive distribution or marginal likelihood. In general, for a model k with data density distribution π k (y|θ k ) and prior density π k (θ k ), θ k ∈Θ k , the predictive probability is π k (y)= ∫π(y|θ k ) π(θ k ) dθ k , which is the denominator of Bayes' theorem. If we want to compare two alternative models, for instance a probit model with k=1 versus a four mixture model with k=4, then the Bayes factor is defined as the ratio of their respective predicted probabilities, i.e. isπ (k4) /π (k1) . In the case of equal prior probabilities for each model, the Bayes' factor becomes B k4,k1 = π k4 (y|θ k14 )/ π k1 (y|θ k1 ). For any particular value θ * of the parameters, the predictive likelihood can be obtained from Bayes' theorem as Ln (π(y)) = Ln (π(θ * )) + Ln (π(y|θ * )) -Ln (π(θ * |y)). Hence, the predictive likelihood can be evaluated by combining the values of the prior, the likelihood and the posterior density at any value of θ * .

The results of the computed marginal likelihoods or predictive distributions for both the probit model and the mixture model are presented in Table 3. The larger probability for the mixture model means that the Bayes' factor favors the mixture model against the more restricted probit model. Thus, the mixture model adds to the explanation of the latent WTP data, although it increases the parameterization of the model. The implication is that modeling unobserved heterogeneous preferences results in better predictions and more accurate estimation of health risk benefits. That is, health preferences for the risk of a respiratory illness are subject to sample heterogeneity. The assumption of a common distribution structure is not maintained with these data and would lead to biased results of the welfare estimates of the health policy measures.

Conclusions

Health preferences could be heterogeneous because of the different perceptions people in society might have about their health status, their environment and their probability of becoming ill. The level of risk aversion to particular health states could also vary across individuals, leading to heterogeneity in preferences. The result is that the values elicited with heath valuation methods such as contingent valuation can not be assumed to be The model has been applied to data on the valuation of the reduction of the probability of a respiratory illness potentially caused by air pollution or natural factors. The mixture model selects four different distributions embedded within the data, which might respond to different preferences for health risks. The group with the largest proportion is given by those individuals who are not willing to pay any amount or might be on the negative side when asked for an amount of money for the reduction of the risk of respiratory illness.

These people might have developed negative preferences based on their perceptions of the potential risks associated with environmental pollution. Thus, the mixture distribution approach is particularly appropriate to model zero and negative values in situations where there might be potential conflicts of interest or different opinions about the implementations of the health prevention measures. Subjects who were willingness to pay positive amounts also showed different degrees of intensity in their preferences, with the largest proportion represented by the group with the higher positive values. The mixture model improves prediction against the model which imposes a unique distribution on the data. The latter structure is also rejected by utilizing the H-H test for probit misspecification, calling for a more flexible approach to model the data. Thus, the mixture approach reduces bias in the estimation of the benefits from reducing health while providing a more accurate representation of the different preferences within society. Further research should focus on the development of practical methods to identify the characteristics of these groups of individuals. In addition, the mixture model 18 Horowitz and Hardle (1994) propose a procedure to test for rigid structures for the error distribution in endogenous dichotomous choice models. In general, the problem is to test the specification of the following model:

[ ] ( ) β , x f x X Y E = = (A.1)
where Y ∈ R, x ∈ R k and f is a function in R which is known and associated to the parameter vector β ∈ Θ ⊂ R k . Let us consider a random sample{(x 1 ,y 1 ), … , (x n ,y n )} where (x i ,y i ), i = 1,…n, is jointly distributed ( )

β y x g , defined in R k × R. The general regression function ( ) β , x f
can be considered a simple index model, such as probit,

where [ ] ( )

β x x X Y P x X Y E ′ Φ = = = = = ) 1 ( . When equation (A.1) is correct, we have that ( ) [ 0 , = - x x f Y E β ]
) . Thus, the functions of the explanatory variables should be correlated with the errors

( β , x f Y - . If we assume that function ( ) β , x f can be written as a function F of ) , ( β x v
, where F and v are known functions defined in R, then equation (A.1) takes the following form:

[ ] ( ) [ β , x v F x X Y E = = ] (A.2) If we have the simple index model β β x x v ′ = ) , ( and function
( ) [ ] ( ) β β x F x v F ′ = ,
, where F is the link function, then the test involves the comparison of the following hypotheses:

[ ] [ ] ) ( : ) ( : 1 0 v H v X Y E H v F v X Y E H = = ′ = = ′ β β
where H is an unknown function. When the link function is probit, the test statistic is as follows

: [ ][ ] ) ( ) ( ~ ) ( ) ( ' ' 1 ' ' β β β β i i i n i i i i n X X F X Y X u h T Φ - Φ - = ∑ = (A.3) 24 F o r P e e r R e v i e w Appendix 2

. Gibbs Sampler

The conditional posterior densities over which Gibbs sampling is applied are the following: and , represents chi-square and beta pdf's respectively. Then, the vector of means are In the next question we want to know how much people value a reduction in the possibility that the population of Las Palmas have of becoming ill because of a respiratory problem. The following card shows the symptoms of an episode of respiratory illness, which could be associated with some diseases such as asthma of bronquitis.

( ) ( )[ ] ( )[ ] ⎪ ⎩ ⎪ ⎨ ⎧ = ∞ + ′ = ∞ + ′ = - - 0 , 1 , 1 ' 1 ' i i L i i i i i L i i i i si y , B - I e x WTP si y , B I e x WTP Y,θ WTP f i i σ α β φ σ α β φ (A2.1) ( ) ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = -1 , , , , , β β β φ π σ α β π H WTP Y (A2.2) ( ) ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = -1 , , , , , α α α φ π σ β α π H WTP Y s.t. α 1 < α 2 < ... < α k 7 (A2.3) ( ) ( ) j j j j j s f WTP Y s f υ σ π α β σ χ 2 2 2 , , , , = s.t. σ 1 < σ 2 < ... < σ k 6 (A2.4) ( ) ( ) k k Beta T r T r T r f WTP Y f + + + = ..., , , , , , 2 2 1 1 π σ α β π (A2.
∑ = = N i ij j e T 1 2 χ f Beta f ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = ∑ = - N i L i i i WTP x H H 1 1 σ β β β β ; and ( ) ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - + = ∑ = - β α α α α ' 1 1 i i N i i x WTP e H H . Thus, L i = (i:e ij =1); ( ) ( ) [ ] j i j i j i x WTP j L P σ β α π 2 exp 2 ' - - - = = and the covariance matrix are ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = ∑ = N i i i L x x h H H i 1 ' β β and ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = ∑ = N i i i e e H H 1 ' α α . Then, ∑ - - + = i i j i ij j j x WTP e s s
These types of episodes are possible in the city of Las Palmas because of traffic pollution, and tend to increase in some times of the year when there are adverse climatic conditions, such as when the North Atlantic Alisios winds stop blowing through the Canary Islands. A possible solution which could reduce the incidence of the episodes of respiratory illness would be to implement traffic restrictions in the days that traffic emissions could cause more problems. The reduction in the possibility of having an episode of respiratory illness depends on the type and effectiveness of the measures which are undertaken. Let us considered that this possibility could be reduced in 20%, 50% or 70% as shown in the following card.

Assume that measures involving traffic restrictions in some areas are undertaken only if they benefit to the majority of the population. Now consider how much you value in money terms the reduction in the possibility of becoming ill of a respiratory illness.

Also recall all the personal damages that these symptoms would imply to you, including the costs on medicines and working time lost, as well as your limited income and your necessary expenses. Would you be willing to pay x euros for the reduction in the possibility of become ill in x %?

  including z, and -1 i WTP ξ is the inverse of the willingness to pay function with respect to ξ i . The linear predictor is linked to the probability of a positive response by a known cumulative distribution function {F ξ }or link function.

5

  Page

7

  Page Y = (y 1 , y 2 , ..., y n )' and θ = (β, α, σ, π), the posterior distribution following data augmentation π(θ|Y,WTP) and the conditional density of the latent variable f(WTP i |Y,θ) are known in a manageable form 4 .

  value of h determines the bandwidth of the confidence interval in the kernel regression. u(.) is a weighted function which gives smaller relevance to extreme observations, and (.) ~i F is a kernel regression on the data with order r ≥ 2. Parameter β is estimated under the null assumption using ML for the probit model. (.) ~i F should have the properties of being asymptotically unbiased and uncorrelated with Y . The first property is obtained by using Bierrens (A.6) correction. The second follows by eliminating observation i from the computation of the kernel regression or "leave-one-out" estimation. The asymptotic variance is replaced by the consistent estimator. We utilize a quartic kernel, which satisfies the requirements for Horowitzi 2 T σ 6 proposition 1 to hold:

9

  Page 11 of 33Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK

Page

  were discussed. The final sample was obtained from 368 in-person interviews addressed to the general population of the island. The sample was taken randomly in the households subject to quotas for representative age and sex groups. Protest responses as indicated by individuals who opposed the valuation process because the government should pay for the risk prevention program amounted to 10.5 percent of the sample. Those subjects who were not willing to pay any positive quantity represented another 11.2 percent.The valuation scenario presented the subject with a hypothetical reduction in the possibility of becoming ill with an episode of respiratory illness. The wording of the scenario is shown in Appendix 3. The symptoms of this type of health problem were described by a card which specified also the number of days of their duration. The subjects were informed that the episodes could be associated with some diseases such as asthma and bronquitis. The possible causes of the episodes were attributed to traffic pollution in some areas of the city of Las Palmas. The possibilities of having an episode could also increase due to whether 10 in the days of the year when Alisios winds stop blowing through the Canary Islands. The possibility of developing the episode of respiratory illness could be reduced by implementing policy measures involving traffic restrictions in the days of higher risks.

Page

  The elicitation format was single bounded dichotomous choice involving a single binary yes/no question to a given price for a concrete risk reduction level. Each subject received up to four binary choice questions for the alternative risk reduction levels which were 11

Page

  The value of k for the number of normal mixture distributions was chosen by the best Akaike and Bayesian information criteria statistics. The result was k=4, i.e. the best model was obtained with a mixture of four distributions reflecting four latent segments in the 14 preferences for the health risks prevention measures. The four latent classes differ particularly in mean WTP. There are three classes with postive means while there is one class with a large negative mean WTP. The class with the largest mean value (latent segment 3) is willing to pay 41.70 Euros for the risk prevention measures. The other two classes with positive values (latent segments 1 and 2) are quite close in their mean WTP, although there is less uncertainty in segment 2 due to its shorter confidence interval.Negative WTP is reflected in latent segment 4 with a mean WTP of -31.78 Euros. These are people who might object the adoption of preventive measures for the risk of respiratory illness, probably because they oppose the traffic resctrictions or because they think they are not affected by the contaminated air. Overall, mean WTP with the BMNP model is practically equal to the one obtained with the basic model that does not take into account unobserved sample heterogeneity. However, it is clear that the policy implications are quite different. Whereas in the basic model every segment of the population would be treated the same, in the flexible model the population could be grouped into different classes, each charged with its corresponding mean WTP.The flexible model also gives us the sample proportions for each latent segment. It can be seen in Table3that the largest proportion is obtained for the segment with negative WTP (38% of the sample) while the smallest proportion is found for those subjects who are willing to pay the smallest positive amount (12%). The subjects who are willing to pay the largest amount to prevent respiratory illnesses represent 29% of the sample.

Page

  attributed to a unique representative individual. Rather, there might be different patterns of behavior across the sample which are reflected in different health values. Rigid parametric approaches to model health valuation empirical data impose strong assumptions, since all individuals are treated as emerging from the same data generation process. In this paper we have considered a mixture of normal distributions to model unobserved heterogeneity with the dichotomous choice contingent valuation method. The model is flexible because any empirical distribution can be represented by increasing the number of mixture distributions. Flexible approaches can accommodate different patterns of behavior resulting from the dispersion of individuals' preferences. In order to simplify computations, the mixture model is estimated utilizing Bayesian methods involving data augmentation and Gibbs sampling. The results of the mixture model are compared with the probit model 17 ' factor, showing that it leads to a larger predictive probability under the assumption of a four mixture distribution.

  5) where φ(.) I[a,b] is the normal pdf truncated in the interval [a,b],

Table 1

 1 presents the results of the test statistics for different values of h and their p-values.The probit model is not rejected for low values of h. However, for values of h above 0.55 the test leads to the rejection of this rigid specification at the 5 percent level, suggesting that 12

	Page 14 of 33

Table 1 .

 1 Horowitz and Härdle (H-H) test results

	Yes	No	Do not know/ Refuse answer

Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK

Table 2 .

 2 Description of explanatory variables in WTP models.

	RP	= 1 for a 20% risk reduction, =2 for a 40% risk reduction, =3 for a 60%
		risk reduction and =3.5 for a 70% risk reduction
	HS	Self-reported health status index for the individual (1-4), where 1 is very
		good and 4 is very bad.
	AGE
	INCOME	Annual personal income of the subject in Euros.
	EXERCISE	= 1 if individual usually does exercise at least twice a week, 0 otherwise
	CHRONIC	= 1 if individual has a chronic disease, 0 otherwise.
		28

Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK

Table 3 .

 3 Descriptive Statistics

	Variable	Mean	Standard Deviation	Min	Max
	RP	2.3315	0.8213	1	4
	CHRONIC	0.6277	0.4840	0	1
	DAYS	5.6311	2.8543	4	7
	INCOME	1116.07	436.14	0.4808	2.4040
	HS	2.9103	0.7158	1	4
	AGE	36.5815	10.6616	14	72
	EXERCISE	0.56521	0.4964	0	1

Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK

Table 4 .

 4 WTP valuation functions (standard error in parentheses)

	Variables	Bayesian normal (k=1)	BMNP(k=4)
	DAYS	31.90 (12.85)	23.87 (4.87)
	RP	20.0472 (9.53)	18.67 (7.10)
	CHRONIC	20.35 (12.08)	18.93 (10.53)
	INCOME	0.0355 (0.0243)	0.0203 (0.01)
	HS	-5.77 (8.81)	-6.06 (6.99)
	AGE	-2.80 (1.37)	-2.11 (0.55)
	EXERCISE	11.69 (11.58)	10.45 (9.83)
	σ	96.10 (46.06)	80.71 (11.99)

Table 5 .

 5 Model Selection based on Marginal Likelihood values (the preferred model is highlighted in bold)

		Marginal Likelihood
	K=1	-2366.83
	K=2	-1898.14
	K=3	-1836.56
	K=4	-1822.68
	K=5	-1832.37

Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK

Submitted Manuscript

This restriction guarantees the identificability of the model[START_REF] Geweke | Mixture of Normals Probit Models[END_REF]

Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript Editorial Office, Dept of Economics, Warwick University, Coventry CV4 7AL, UK Submitted Manuscript

Under the null hypothesis and some regularity conditions, it can be shown that T n is

The weighted function u(v) can be set to one if it is inside a bounded interval and zero in other case.

In order to obtain Bierrens's (1987) property of consistency, ) ( ~v F can be defined as a linear combination of Kernel regressors with different bandwidth (h and s respectively). This correction relies on the following result:

where is the density function for v estimated by

25