Mustapha Er-Riani
email: m.erriani@fstt.ac.ma

K Chetehouna

On the critical behaviour of exothermic explosions in Class A geometries

Keywords: Thermal criticality, chemical kinetics, Homotopy perturbation method, Domb Sykes method

The aim of this work is to apply the homotopy perturbation method for solving the steady state equations of the exothermic decomposition of a combustible material obeying Arrhenius, Bimolecular and Sensitised laws of reaction rates. These equations are formulated on some Class A geometries (an infinite cylinder, an infinite slab and a sphere). We also investigate the effect of Frank-Kamenetskii parameter on bifurcation and thermal criticality by means of the Domb-Sykes graphical method.

The safety in transport and storage of combustible materials is a key issue in pyrotechnic applications. These materials are often subjected to self-ignition. This internal heating occurs when an explosive substance is brought to a sufficient temperature so that the process of decomposition begins to produce significant exothermic effects. This involves a thermal runaway phenomenon accompanied with an increase of the temperature producing a rapid thermal decomposition. The understanding of the factors that control this phenomenon is of fundamental importance in many industrial processes.

This phenomenon was first introduced in the 1930s by Semonov, Zeldovith and Frank-Kamentskii and their pioneering contributions were summarized in [START_REF] Frank | Diffusion and Heat Transfer in Chemical Kinetics[END_REF]. Furthermore, Frank-Kamentskii developed the steady-state theory of thermal explosion. In this theory, the Frank-Kamenetskii approximation allows us to determine the critical values, which constitute limit values not to be exceeded to avoid the phenomenon of self-ignition. Some studies deal with a chain of several reactions in this phenomenon of self-ignition. Several works in the literature have applied this theory to different combustible materials geometries. Boddington et al. [START_REF] Boddington | Thermal explosions and the disappearance of criticality at small activation energies: exact results for the slab[END_REF][START_REF] Boddington | Thermal Explosions, Criticality and the Disappearance of Criticality in Systems with Distributed Temperatures. I. Arbitrary Biot Number and General Reaction-Rate Laws[END_REF][START_REF] Boddington | Theory of Thermal Explosions with Simultaneous Parallel Reactions. I. Foundations and the One-Dimensional Case[END_REF] have considered the special case of two-step parallel exothermic reactions for the infinite slab where solutions by quadrature were possible. Graham-Eagle and Wake [START_REF] Graham-Eagle | Theory of Thermal Explosions with Simultaneous Parallel Reactions. II. The Two-and Three-Dimensional Cases and the Variational Method[END_REF] Recently Ajadi and Gold'shtein [START_REF] Ajadi | Critical behaviour in a three-step reaction kinetics model[END_REF] employed a three-step reaction kinetics model (initiation, propagation and termination steps). The calculation of the criticality was made by means of a variational method for (infinite slab, infinite cylinder, and sphere) under Arrhenuis laws of reaction rates by using an effective activation energy approximation. Balakrishnan et al. [START_REF] Balakrishnan | Critical values for some nonclass A geometries in thermal ignition theory[END_REF] calculated the critical values for some non Class A geometries (infinite square rod and cube).

Their critical values were found using the finite difference method.

In applied mathematics or engineering problems, numerical methods commonly used such as finite difference, finite element or characteristics method, need large size of computational works due to discretization and usually the effect of round-off error causes loss of accuracy in the results. In addition to this drawback, these methods with limited precision include slow runtimes, numerical instabilities and difficulties in handling redundant constraints.

Analytical traditional methods commonly used for solving these problems can be very useful, especially for the calibration of numerical calculations. Among these approaches, the classical perturbation method is based on the existence of small parameters but the overwhelming majority of linear and nonlinear problems have no small parameters at all.

To overcome this shortcoming, the homotopy perturbation method (HPM) was first introduced by He [START_REF] He | A coupling method of a homotopy technique and a perturbation technique for non-linear problems[END_REF][START_REF] He | Homotopy perturbation method: A new nonlinear analytical technique[END_REF][START_REF] He | Homotopy perturbation method for solving boundary value problems[END_REF][START_REF] He | An elementary introduction to the homotopy perturbation method[END_REF]. The main idea of this method is to introduce an embedding parameter p [0,1] to construct a homotopy and give solutions of the deformed problem as power series in p . When p = 0 the system of equations is reduced to a simplified form which admits exact solutions and at p = 1, the system takes its original form and gives the desired solutions. This method has been extensively employed by numerous authors to solve a large variety of linear and nonlinear problems. We can cite the works [START_REF] Ariel | The homotopy perturbation method and analytical solution of the problem of flow past a rotating disk[END_REF][START_REF] Hosseini | Solution of temperature distribution in a radiating fin using homotopy perturbation method[END_REF][START_REF] Rajabi | Application of homotopy perturbation method in nonlinear heat conduction and convection equations[END_REF][START_REF] Siddiqui | Homotopy perturbation method for heat transfer flow of a third grade fluid between parallel plates[END_REF], the list is not exhaustive.

In this paper, we examine the steady-state solutions for the strongly exothermic decomposition of a combustible material of a symmetric class A geometries, uniformly heating, under Arrhenius, Biomolecular and Sensitised kinetics, neglecting the consumption of the material.

The contribution of the present work is twofold. First, we calculate the temperature field using the HPM in a symbolic computational language. The second aim is to study the thermal criticality conditions of the problem.

Critical values for different geometries are found by using the Domb-Sykes technique [START_REF] Domb | On the susceptibility of a ferromagnetic above the Curie point[END_REF] as a useful tool to extract singularities and to perform analytic continuation. This method has shown to be very accurate and simpler to implement, when compared with variational method [START_REF] Graham-Eagle | The theory of thermal explosions with simultaneous parallel reactions.III. Disappearance of critical behaviour with one exothermic and one endothermic reaction[END_REF][START_REF] Ajadi | Critical behaviour in a three-step reaction kinetics model[END_REF] or Hermite-Padé approximants method [START_REF] Makinde | Strongly exothermic explosions in a cylindrical pipe: a case study of series summation technique[END_REF][START_REF] Makinde | Exothermic explosions in a slab: A case study of series summation technique[END_REF][START_REF] Makinde | Exothermic explosions in symmetric geometries-An exploitation of perturbation technique[END_REF].

The structure of this paper is as follows. Section 2 presents the boundary value problem governing the ignition of a viscous combustible material for symmetric class A geometries. In section 3, we will apply the homotopy perturbation method to this nonlinear boundary value problem. The fourth section is assigned to a brief description of the Domb and Sykes method and its application to calculate the bifurcation points. Section 5 is devoted to the results obtained by the proposed method and comparison with other works will be performed.

Conclusions will appear in Section 6.

Mathematical model

We consider the steady-state solutions for the strongly exothermic decomposition of a viscous combustible material. Neglecting the reactant consumption, the equation for the temperature T(r) may be written in terms of physical variables together with the boundary conditions, as:

0 0 2 2           RT E m e h KT k A QC dr dT r j dr T d  (1) 0 ) ( , 0 ) 0 ( T a T dr dT   (2)
with T the absolute temperature, T 0 the wall temperature, k the thermal conductivity of the material, Q the heat of reaction, A the rate constant, E the activation energy, R the universal gas constant, C 0 the initial concentration of the reactant species, h the Planck's number, K the Boltzmann's constant, thevibration frequency, a the geometry half width and r the radial distance in the normal direction. m is the numerical exponent, such that 2

1 , 0 , 2   m
represent numerical exponent for Sensitised Arrhenius, and Bimolecular kinetics (Boddington et al. [START_REF] Boddington | Thermal Explosions, Criticality and the Disappearance of Criticality in Systems with Distributed Temperatures. I. Arbitrary Biot Number and General Reaction-Rate Laws[END_REF], Makinde [START_REF] Makinde | Exothermic explosions in a slab: A case study of series summation technique[END_REF], Bowes [START_REF] Bowes | Self-Heating, evaluating and controlling the hazard[END_REF], Bebernes and Eberly [START_REF] Bebernes | Mathematical Problems from Combustion Theory[END_REF], Okoya [START_REF] Okoya | On the behaviour of solutions to a system of ordinary differential equations modelling branched chain reaction[END_REF]) and j = 0, 1, 2 is the geometry factor representing respectively Slab, Cylindrical and Spherical geometries.

The following dimensionless variables and parameters are introduced in Eq. (1):

0 2 0 0 2 0 2 0 0 , , , ) ( RT E m m m m e Rk h T K C QEAa E RT a r r RT T T E           
to obtain the dimensionless governing equation together with the corresponding boundary conditions. We get the following equations:

0 1 exp ) 1 ( 2 2                  m r d d r j r d d (3) 0 ) 1 ( , 0 ) 0 (     r d d ( 4
)
where ,  represent the Frank Kamenetskii and activation energy parameters respectively.

Hereafter we will suppress the bar symbol for clarity.

Homotopy perturbation method solution

In this section, we will apply the homotopy perturbation method (HPM) to nonlinear ordinary differential equation [START_REF] Boddington | Thermal Explosions, Criticality and the Disappearance of Criticality in Systems with Distributed Temperatures. I. Arbitrary Biot Number and General Reaction-Rate Laws[END_REF]. According to this method, we can construct a homotopy as follows:

  0 1 exp ) 1 ( ) 1 ( 2 2                               m dr d r j dr d p p L (5)
where

  dr d r j dr d      2 2 L ( 6
)
Applying the perturbation technique [START_REF] Nayfeh | Problems in Perturbation[END_REF], we can assume that the solution of Eq. ( 5) can be expressed as a series in p

        ... , 2 2 1 0     r f p r pf r f p r  (7)
where p  [0,1] is an embedding parameter. When p = 0, we can obtain the initial guesses; when p = 1, Eq. ( 5) turns out to be the original one.

Setting p = 1, we obtain an approximate solution of (3):

        ... , lim 2 1 0 1       r f r f r f p r p   (8)
One has to substitute relation [START_REF] Gelfand | Calculus of Variations[END_REF] into the governing equation ( 5), collect the powers of p and obtain a sequence of differential equations and boundary conditions. The solution for the temperature field for Sensitised, Arrhenius and Bimolecular reaction rates for different geometries are given as

           3 2 2 2 2 2 1 5 1 24 1 1 2 1 0 , , p O m r r p r p j p r              (9)            3 2 2 2 2 2 1 3 1 64 1 1 4 1 1 , , p O m r r p r p j p r              (10)            3 2 2 2 2 2 1 7 3 1 360 1 1 6 1 2 , , p O m r r p r p j p r              (11)
According to Eq. ( 8) we get the solution for the temperature field for the three reaction rates.

Bifurcation study by Domb Sykes method

The modelling of physical phenomena often results in nonlinear problems for some unknown function, say f() depending on a parameter . Usually, the problems cannot be solved exactly. The solutions of these nonlinear systems are dominated by their singularities which must be real and positive in order to have a physical sense.

We suppose that up to the point  0 , the solution is analytic in the interval ]- 0 ,  0 [, then one can solve the problem by expanding the solution in a power series

  n n n c f       0 (12)
If an infinite number of c n is known, the radius of convergence, which is the distance from the origin to the nearest singularity, limiting the range of validity of the series [START_REF] He | Homotopy perturbation method for solving boundary value problems[END_REF], can be calculated using D'Alembert's ratio.

However, for most nonlinear problems it is rare to find an unlimited number of terms of the power series [START_REF] He | Homotopy perturbation method for solving boundary value problems[END_REF]. The nearest singularity cannot therefore be determined precisely.

Using symbolic calculus codes, it is now possible to calculate a sufficient number of terms in the series, to study precisely the solution, and there exist a variety of methods devised for extracting the required information of the singularities from a finite number of series coefficients. The most frequently used methods are the ratio-like methods, such as the Domb-Sykes method [START_REF] Gutmann | Asymptotic analysis of power series expansions[END_REF][START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers[END_REF], Neville-Aitken extrapolation [START_REF] Gutmann | Asymptotic analysis of power series expansions[END_REF], and semi numerical approximant methods, such as Padé approximants [START_REF] Baker | Padé Approximants[END_REF] or Hermite-Padé approximants [START_REF] Makinde | Strongly exothermic explosions in a cylindrical pipe: a case study of series summation technique[END_REF][START_REF] Makinde | Exothermic explosions in a slab: A case study of series summation technique[END_REF][START_REF] Makinde | Exothermic explosions in symmetric geometries-An exploitation of perturbation technique[END_REF].

Herein we are concerned with the bifurcation analysis by analytic continuation as well as with the dominant behavior of the solution by using partial sum [START_REF] Gelfand | Calculus of Variations[END_REF]. This may be done calculating the nearest singularity to the origin. A useful tool for extracting this singularity is the Domb-Sykes method. This technique is easier to implement than others and when it is improved, it gives a precise value of the singularity. This method, successfully applied in various problems [START_REF] Aziz | Perturbation Methods in Heat Transfer[END_REF][START_REF] Er-Riani | Shapes of liquid drops obtained using symbolic computation[END_REF][START_REF] Schwartz | A semi-analytic solution for nonlinear standing waves in deep water[END_REF][START_REF] Ohkitania | Numerical study of singularity formation in a class of Euler and Navier-Stokes flows[END_REF], is presented here.

According to Fuchs (Bender and Orszag [START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers[END_REF]), there are two possible forms:

                           , , , 1 , 0 for ln , , , 1 , 0 for 0 0 0 0 n n c f n n n               (13)
The radius of convergence 0  of the series (12) may be found by d'Alembert ratio:

n n n c c 1 0 lim      (14)
Of course only a finite number of coefficients c n are known, so that it is difficult to obtain precisely the limit. Domb and Sykes have suggested that the inverse ratio

1  n n c c
has the following expansion: ), which is appropriate for this kind of sequences and can be defined as follows.

) / 1 ( 1 1 1 1 0 1 n h n o n c c a n n n                      (
If n a can be written in the following form:

     2 2 1 0 n S n S S a n , then the Richardson extrapolation is:             m k m k n m k mn k m k k n a s 0 )! ( ! 1 ( 16
)
This new sequence has a quicker convergence to the limit S 0 .

Results and discussion

Using a computer algebra system, we obtained the first 30 terms of the solutions series ( 9)- [START_REF] He | Homotopy perturbation method: A new nonlinear analytical technique[END_REF].

In order to verify numerically whether the proposed methodology leads to high accuracy, we evaluate the numerical solutions using a collocation method proposed by Shampine et al. [START_REF] Shampine | Solving Boundary Value Problems for Ordinary Differential Equations in MATLAB with bvp4c[END_REF].

In Figures1-3 below, we have plotted the numerical and the HPM solutions in all class A geometries and for all numerical exponent m. These figures show a perfect agreement of both solutions. The obtained solutions, in comparison with the numerical solutions admit a remarkable accuracy. A clear conclusion can be drawn from the numerical results that the HPM provides highly accurate numerical solutions for nonlinear differential equations.

For a further information, Figure 4 illustrates a comparison between the exact solutions known for the special cases (m = 0,  = 0, j = 0, j = 1) when a closed form is available [START_REF] Ajadi | Critical behaviour in a three-step reaction kinetics model[END_REF] and the series solutions obtained by using the homotopy perturbation method. The obtained results are found to be in good agreement with the exact solutions. We can notice that the deviation between HPM and exact solution do not exceed 0.2 % for both slab and cylindrical geometries. In the following we will focus on the calculation of radius of convergence  cr of the series solutions ( 9)- [START_REF] He | Homotopy perturbation method: A new nonlinear analytical technique[END_REF]. Up to the point  c , the solution  is analytic, and has a real singularity at the value  cr. It would be interesting to know the exact nature of the singularity, and if this singularity corresponds to a bifurcation point. Let the maximum temperature

      , , 0 max 
be a characteristic quantity which qualifies the solution. We rearrange max



in a series of  to write:

     30 0 max n n n c    (17)
where different coefficient c n for  = 0 are given in table 1. Seeing Figure 5, one can say that the ratio , m = 0, j = 1.

.

, m = 0, j = 1.
It is possible to calculate the limited series defining max  = max  () by a division procedure.

This function is plotted in Figure 7.

The graph of max  , as shown in Figure 6(a), is not compatible with the expected singularity (13). The curve max  versus should be in the form of Figure 6(b). To analyse the paradox, the inverse function ( max  ) is considered. Its series is given by inverting the series max  .

One can see in Figure 7 that this function has a horizontal tangent for the values in the case j = We summarize the results of all Class A geometries in The critical values ( cr and  cr ) were computed using the improved Domb-Sykes applied to the temperature field obtained by HPM in all Class A geometries.

In order to verify the accuracy of this approximate method, we compare the critical values obtained using this method with the exact methods [START_REF] Frank | Diffusion and Heat Transfer in Chemical Kinetics[END_REF][START_REF] Boddington | Thermal explosions and the disappearance of criticality at small activation energies: exact results for the slab[END_REF], variational methods [START_REF] Graham-Eagle | Theory of Thermal Explosions with Simultaneous Parallel Reactions. II. The Two-and Three-Dimensional Cases and the Variational Method[END_REF][START_REF] Ajadi | Critical behaviour in a three-step reaction kinetics model[END_REF] and 20

Hermite-Padé approach [START_REF] Makinde | Strongly exothermic explosions in a cylindrical pipe: a case study of series summation technique[END_REF][START_REF] Makinde | Exothermic explosions in a slab: A case study of series summation technique[END_REF][START_REF] Makinde | Exothermic explosions in symmetric geometries-An exploitation of perturbation technique[END_REF] for the special case of  = 0.The method used in the present work gives results approximately equal to those discussed in literature, as shown in Table 4. Tables 5-7 depict the variations of δ cr and θ cr with  In these tables, we observe that the magnitude of thermal criticality conditions for a viscous combustible material with high activation energy ( = 0) is lower than the one for moderate values of activation energy ( = 0.01,  = 0.1). This implies that for moderate values of activation energy, the criticality varies from one reaction to another, as shown in Tables 567 7. Variation of δ cr and θ cr with respect to  (Sphere).

Exact method Variational method

Conclusion

We studied the problem of exothermic explosion of a viscous combustible in Class A geometries under Arrhenius, Bimolecular and Sensitised laws of reaction rates with the homotopy perturbation method. The results show that this method provides excellent approximation of the solution of this nonlinear system with high accuracy.

A bifurcation study is performed with the Domb-Sykes graphical method to calculate the critical value of this problem. The results show that these critical values increase with the activation energy.

2 , 3 , 4 ]

 234 considered a system of simultaneous exothermic reactions. They extended the investigations of Boddington et al. [to the other two geometries of the infinite circular cylinder and the sphere. A variational method was used to evaluate the critical values of Frank-Kamentskii parameterand the maximum temperature. The same authors[START_REF] Graham-Eagle | The theory of thermal explosions with simultaneous parallel reactions.III. Disappearance of critical behaviour with one exothermic and one endothermic reaction[END_REF] extended the treatment of simultaneous reactions to the case where one reaction is exothermic and the other endothermic, leading to the phenomenon of the disappearance of criticality or transition, which can also happen with a single reaction in the case of very low activation energy. In this case, the Frank-Kamenetskii approximation is no longer valid. They used a variational method [7] to determine numerically the values of the parameter which characterizes the transition in the parameter space.

15 )

 15 Domb and Sykes have pointed out that it is more reliable to plot 1  n n c c versus 1/n, i.e. to bring n→∞ to the origin, rather than plotting 1/n is known as Domb-Sykes plot. In this plot, the intersection of the straight line (of this line gives the exponent . Unfortunately, n a is often a slowly converging series. A good way to improve the convergence is to use the Richardson extrapolation (Bender and Orszag [27]

Figure 1 .

 1 Figure 1. Comparison of slab numerical and HPM results of the boundary value problem (3-4) for ε = 0.1,  = 0.5. Solid lines: Numerical solutions; Symbols: HPM solutions.

Figure 2 .

 2 Figure 2. Comparison of cylinder numerical and HPM results of the boundary value

Figure 3 .

 3 Figure 3. Comparison of sphere numerical and HPM results of the boundary value

Figure 4 .

 4 Figure 4. Comparison of HPM results and exact solutions of the boundary value

Figure 5

 5 Figure 5. Sykes-Domb plot for 0 , max

Figure 6 .

 6 Figure 6.  max versus ; (a) calculated, (b) expected form.

Figure 7 .

 7 Figure 7. Plot of the function  = ( max ).

Table 1 . Coefficients c n of 0 ,

 1 (i.e. turning point) as shown in a sketch of bifurcation diagram in Figure6.

		n a		1  n n c c	in the simple geometry of infinite circular cylinder converges;
	however, it is extremely difficult to conclude to a precise value of radius of convergence. It is
	necessary to use some process of convergence acceleration. The Richardson extrapolation is
	appropriate for this kind of convergence. The results of Table 2 show that the radius of
	convergence is				2	. Two solution branches are therefore identified with a bifurcation
				cr
	point at	 .Two solution branches (type I and II) are identified with a bifurcation point at
		cr
	n	c n				n	c n	n	c n	n	c n	n	c n
	1	0.25				7	2.337e-4	13	1.455e-6	19	1.290e-8	25	1.338e-10
	2	4.687 e-2 8	9.588e-5	14	6.515e-7	20	5.978e-9	26	6.310e-11
	3	1.302 e-2 9	4.024e-5	15	2.939e-7	21	2.778e-9	27	2.982e-11
	4	4.272e-3		10	1.720e-5	16	1.334e-7	22	1.296e-9	28	1.412e-11
	5	1.538e-3		11	7.465e-6	17	6.096e-8	23	6.064e-10 29	6.700e-12
	6	5.874e-4		12	3.279e-6	18	2.798e-8	24	2.845e-10 30	3.184e-12

cr  max   

Table 2 .

 2 Values of the Richardson extrapolation s nm for

	 cr	,  	0

table 3 ,

 3

	4 and 5 while		cr 		 0 cr  ,	,			is

calculated by improving the series (17) by a suitable Padé approximant

[START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers[END_REF]

.

Table 3 .

 3 Computation showing criticality for Sensitised, Arrhenius and Bimolecular reactions for  = 0 and for all Class A geometries.

Table 4 .

 4 Comparison of critical values obtained by different methods.

Table 5 .

 5 . Explosion in Bimolecular reactions seems to occur faster than in Arrhenius and Sensitised reactions. Variation of δ cr and θ cr with respect to (Slab).

	m		 cr 	 cr 
	-2	0.01	0.90624	1.24566
		0.1	1.31389	2.22239
	0	0.01	0.8878	1.19841
		0.1	0.98819	1.52434
	0.5	0.01	0.88331	1.19541
		0.1	0.93221	1.42024

Table 6 .

 6 Variation of δ cr and θ cr with respect to  (Cylinder).

	m		 cr 	 cr 
	-2	0.01	3.45641	1.67123
		0.1	4.32871	2.1177
	0	0.01	3.42112	1.62756
		0.1	4.22245	2.21856
	0.5	0.01	3.15167	1.63331
		0.1	4.04567	2.23811
	Table

Acknowledgements

This work is supported by the French Ministry of Research and Education through invited professor contract. The authors are highly grateful to "ENSI de Bourges" for providing excellent research environment and facilities and to Professor J.H. for his valuable comments.


Zeynep Akata
email: zeynep.akata@xrce.xerox.com

Christian Thurau
email: christian.thurau@iais.fraunhofer.de

Non-negative Matrix Factorization in Multimodality Data for Segmentation and Label Prediction

With the increasing availability of annotated multimedia data on the Internet, techniques are in demand that allow for a principled joint processing of different types of data. Multiview learning and multiview clustering attempt to identify latent components in different features spaces in a simultaneous manner. The resulting basis vectors or centroids faithfully represent the different views on the data but are implicitly coupled and they were jointly estimated. This opens new avenues to problems such as label prediction, image retrieval, or semantic grouping. In this paper, we present a new model for multiview clustering that extends traditional non-negative matrix factorization to the joint factorization of different data matrices. Accordingly, the technique provides a new approach to the joint treatment of image parts and attributes. First experiments in image segmentation and multiview clustering of image features and image labels show promising results and indicate that the proposed method offers a common framework for image analysis on different levels of abstraction.

Motivation and Background

The rise of the social web and the user generated content movement have turned the Internet into a virtually limitless repository of annotated and rated multimedia data. For example, as of this writing, there are more than 4.5 billion images available on flickr most of which are tagged, rated, categorized, and appraised by the community. This development offers tremendous possibilities for research on image understanding but also calls for methods that allow for an integrated processing of different types of data.

Our goal is a principled joint treatment of image features and image tags. We present a new technique for multiview clustering that simultaneously determines latent dimensions or centroid vectors in different feature spaces. In contrast to ad hoc methods such as, say, concatenating different types of features into a single descriptor, multiview clustering is faithful to the different characteristics of different descriptors. Since latent components or centroids are jointly estimated, multiview techniques allow for advanced inference. Since for every centroid in one feature space there is a corresponding centroid in another space, transitions between different views are straightforward. This offers auspicious new approaches to segmentation, automatic image tagging, or tag-based image retrieval.

Although they have a long and venerable tradition, there is a renewed interest in multiview learning and multiview clustering. The canonical example of a method that simultaneously uncovers latent components in different spaces is Hotelling's canonical correlation analysis (CCA) [12,2] for which kernelized and probabilistic extension have been proposed as of late [START_REF] Fukumizu | Statistical Consistency of Kernel Canonical Correlation Analysis[END_REF]11,3]. Other recent developments consider extensions of spectral clustering to multiple graphs that encode different types of similarities [START_REF] Zhou | Spectral Clustering and Transductive Learning with Multiple Views[END_REF][START_REF] Tang | Clustering with Multiple Graphs[END_REF].

Our new approach to multiview clustering extends non-negative matrix factorization (NMF) [START_REF] Lee | Learning the Parts of Objects by Non-negative Matrix Factorization[END_REF][START_REF] Lee | Algorithms for Non-negative Matrix Factorization[END_REF] to the joint factorization of several data matrices. It is motivated by the following considerations: i) Similar to principal component analysis (PCA) [START_REF] Chaudhuri | Multi-View Clustering via Canonical Correlation Analysis[END_REF] or singular value decomposition (SVD) [START_REF] Golub | Matrix Computations[END_REF] CCA does not necessarily do justice to purely non-negative data such as color histograms or term frequency vectors. Non-negative matrix factorization, however, typically yields results that can be seen as part-based representations and accommodate human perception.

ii) Methods based on spectral clustering of similarity matrices scale quadratically with the number of data and are therefore prohibitive in modern, largescale data and image analysis problems.

iii) For NMF, on the other hand, there exist efficient algorithms that factorize matrices of billions of entries [START_REF] Thurau | Convex Non-Negative Matrix Factorization in the Wild[END_REF] which may apply to the multiview setting.

In the next section, we clarify the relation between matrix factorization and clustering. Then, in section 3, we briefly review NMF according to [START_REF] Lee | Learning the Parts of Objects by Non-negative Matrix Factorization[END_REF][START_REF] Lee | Algorithms for Non-negative Matrix Factorization[END_REF] and extend this approach toward the joint factorization of different data matrices. In section 4, we present experiments on using multiview NMF in image segmentation, label prediction, and image retrieval. A conclusion will end this contribution.

Matrix Rank Reduction and Clustering

In this section, we briefly review how matrix rank reduction applies to the problem of clustering or vector quantization.

Consider a data matrix X = [x 1 . . . x n ] ∈ R m×n of rank r ≤ min(m, n) whose column vectors x i correspond to feature vectors obtained from some measurement process. Using the singular value decomposition (SVD) [START_REF] Golub | Matrix Computations[END_REF] any matrix X ∈ R m×n can be written as

X = U ΣV T = r j=1 σ j u j v T j (1)
where

U = [u 1 . . . u m ] ∈ R m×m and V = [v 1 . . . v n ] ∈ R n×n
are orthogonal matrices and Σ = diag(σ 1 , . . . , σ r ). The SVD is a popular tool in data analysis because it is known that the optimal solution to the rank reduction problem

min rank( X)=k<r X -X 2 (2)
is given by

X = k j=1 σ j u j v T j = Ũ Σ Ṽ T . (3)
Substituting W = Ũ ∈ R m×k and H = Σ Ṽ T ∈ R k×n , we recognize that X ≈ W H is approximated as a product of a matrix of basis vectors and a matrix of coefficients. This allows for dimensionality reduction, since

x i ≈ k j=1 w j h ji (4)
so that to every data vector

x i ∈ R m there is a coef- ficient vector h i ∈ R k where k < m.
Depending on which constraints are imposed on W and H, one obtains different dimensionality reduction schemes when solving the general matrix factorization problem

min W ,H X -W H 2 . (5)
For instance, principal component analysis (PCA) [START_REF] Chaudhuri | Multi-View Clustering via Canonical Correlation Analysis[END_REF] is recovered from

min W ,H X -W H 2 s.t. W T W = I. (6)
Casting matrix factorization in a yet more general form reveals a connection to vector quantization and clustering. For example, running the k-means algorithm is tantamount to solving

min G,H X -XGH 2 s.t. g T j 1 = 1 (7)
g j 0

h i = [0 . . . 010 . . . 0] T .
Due to the convexity constraints on the columns of G, the resulting basis vectors in W = XG are convex combinations of certain data points in X and since the coefficient vectors in H are unitary vectors, every data point x i in X will be represented by exactly one centroid w j in W .

NMF for Multiview Clustering

In this section, we first summarize non-negative matrix factorization (NMF) and then introduce our generalization of NMF toward multiview clustering.

Factorization of Data via NMF

Orthogonal basis vectors such as determined by PCA or SVD are not always the best choice for dimensionality reduction or clustering [START_REF] Lee | Learning the Parts of Objects by Non-negative Matrix Factorization[END_REF][START_REF] Lee | Algorithms for Non-negative Matrix Factorization[END_REF]25,[START_REF] Donoho | When Does Nonnegative Matrix Factorization Give a Correct Decomposition into Parts?[END_REF]15,[START_REF] Klingenberg | Nonnegative Matrix Factorization: Ill-posedness and a Geometric Algorithm[END_REF]. In particular data that consist exclusively of non-negative measurements cannot be guaranteed to retain non-negativity after projection onto lowerdimensional subspaces that are spanned by its dominant eigenvectors. As an alternative that is true to the non-negative nature of certain data Lee and Seung popularized the idea of non-negative matrix factorization [START_REF] Lee | Learning the Parts of Objects by Non-negative Matrix Factorization[END_REF][START_REF] Lee | Algorithms for Non-negative Matrix Factorization[END_REF]. In computer vision where image data typically consists of non-negative values, NMF was observed to yield superior results in segmentation, feature extraction, motion-, or pose estimation [START_REF] Zass | A Unifying Approach to Hard and Probabilistic Clustering[END_REF]10,[START_REF] Cheriyadat | Non-Negative Matrix Factorization of Partial Track Data for Motion Segmentation[END_REF][START_REF] Thurau | Pose Primitive Based Human Action Recognition in Videos or Still Images[END_REF].

Viewed as a constrained least squares optimization problem, NMF amounts to solving

min W ,H X -W H 2 s.t. W , H 0. (8)
Although ( 8) is convex in either W or H, the simultaneous estimation of basis vectors and coefficients in (8) does not admit a closed form solution and is known to suffer from many local minima. A unique optimum provably exists [25], however, algorithms that are guaranteed to find it are not known to date (see the discussions in [25, [START_REF] Donoho | When Does Nonnegative Matrix Factorization Give a Correct Decomposition into Parts?[END_REF]15,[START_REF] Klingenberg | Nonnegative Matrix Factorization: Ill-posedness and a Geometric Algorithm[END_REF]).

In the work presented here, we consider multiplicative fixed point iterations to find a solution to NMF because their extension to multiview clustering is immediate. In the following, A B ∈ R m×n denotes the Hadamard product of two matrices A, B ∈ R m×n where A B ij = a ij • b ij . The Hadamard division is defined accordingly but for better readability we write A B = A/B.

Concerned with the problem in (8), Lee and Seung [START_REF] Lee | Learning the Parts of Objects by Non-negative Matrix Factorization[END_REF][START_REF] Lee | Algorithms for Non-negative Matrix Factorization[END_REF] randomly initialize the matrices W and H. They derive the following update rules

W ← W XH T W HH T and H ← H W T X W T W H (9)
and prove their convergence using an expectation maximization argument. Next, we will extend this approach to multiview data.

Simultaneous Factorization of Multiview Data via NMF

Our main motivation behind the work presented in this paper is to cluster entities for which there are different types of data available. For instance, images retrieved from flickr can be characterized by means of different abstract image features but at the same time there are user generated tags or labels available that describe their content or formation. We hypothesize that simultaneous clustering of such different views on the data will yield more meaningful clusters and may provide a tool to fill in missing information. In particular, multiview clustering of image features and image tags may provide a way to predict a set of tags given an image or to retrieve relevant images from a database given a set of query tags.

Assuming a set of n different images, it can be characterized by an m × n image-feature matrix X as well as by an l × n term-by-image matrix Y . Our basic idea is to uncover suitable bases W and V for the image-and text features, respectively, which are implicitly coupled via a common coefficient matrix H. In other words, we aim at finding two low rank approximations

X ≈ W H and Y ≈ V H (10)
where W ∈ R m×k , V ∈ R l×k , and H ∈ R k×n . Our solution is to formalize this idea as a convex combination of two constrained least squares problems

min W ,V ,H (1 -λ) X -W H 2 + λ Y -V H 2 s.t. W , V , H 0 (11)
where λ ∈ [0, 1] is user specified constant that allows for expressing preferences for either of the two feature types. Just as with the original NMF problem in (8), the extended problem in (11) does not admit a closed form solution. We therefore adapt the Lee and Seung type fixed point iteration to our case. For the matrices of basis vectors W and V , the update rules immediately carry through and read:

W = W XH T W HH T and V = V Y H T V HH T . ( 12
)
Since the coefficient matrix H now couples two bases, its update is slightly more involved. The simplified version of the fixed point iteration for the coefficients is:

H = H (1 -λ)W T X + λV T Y (1 -λ)W T W + λV T V H . (13)

Discussion

Our choice of a convex combination of the individual optimization problems in (11) is not an arbitrary decision. There is a known close relation between non-negative matrix factorization and probabilistic latent semantic analysis [8,[START_REF] Ding | NMF and PLSI: Equivalence and a Hybrid Algorithm[END_REF]. Assuming an appropriate normalization, NMF can be understood as learning the parameters of a joint probability distribution which is expressed as a product of marginal distributions. By choosing a convex combination of two NMF problems, this analogy may be lifted to the level of learning a distribution of distributions. This is akin to Latent Dirichlet Allocation [18,[START_REF] Blei | Latent Dirichlet Allocation[END_REF] but we will leave possible implications to future work.

We note that by setting λ = 0 or λ = 1 our model and its updates reduce to the original form of NMF. Moreover, the model is not confined to the case of two different types of views. Its extension to convex combinations of p different views is straightforward:

min W i ,H p i=1 λ i X i -W i H 2 s.t. W i , H, λ 0 (14)
λ T 1 = 1
Finally, as with with all alternating least squares schemes, convergence of the extended update algorithm for multiview NMF is guaranteed. We omit the formal proof but sketch the argument: Given H, none of the updates in (12) will increase either term in (11); given W and V , the update in (13) cannot increase the expression in (11).

Experiments

In the following subsections we present first experimental results obtained from using multiview NMF for image segmentation, label prediction, and image retrieval. Note that, so far, these are preliminary experiments intended to validate the approach. We are currently working on extended experimental evaluations to compare the proposed approach to other methods in the literature.

Image Segmentation via Joint Non-negative Matrix Factorization

In a first series of experiments, we apply simultaneous NMF to the problem of image segmentation. We consider color images of natural scenes downloaded from flickr. We convert the RGB pixel values into the LUV color-space because of its alleged perceptual uniformity which ensures that equally distant colors in the color space would be also equidistant perceptually.

In order to segment an image into homogeneous regions, we sample 1000 pixels from each image and build two feature matrices, one containing 1000 three dimensional column vectors of color information and one containing 1000 two dimensional column vectors containing pixel coordinates. This way, we separate color from location and run simultaneous NMF to obtain centroid vectors W and V in the respective spaces that are coupled via the common coefficients H.

We conduct several experiments where we vary the number of centroids k = {4, 10, 20} and the weighting parameter λ = {0.1, 0.5, 0.9}. When λ is larger, more weight is given to the color descriptor of the pixels and when it is smaller more weight is given to the location of the pixels. After random initialization to positive values sampled from a Gaussian distribution, we run the update rules for the matrices W , V and H until convergence but at most 100 times.

Given the results of the training phase, the test phase in these experiments consist in assigning every pixel x of an image to one of the k resulting cluster centroids. Given W and V , we solve min(1-λ) x-W h 2 +λ x-V h 2 for the coefficients h and determine the cluster index c according to c = argmax

j h j . ( 15
)
Figure 1 shows examples of images we considered in our segmentation experiments. The accuracy of the segmentation appears to improve with an increasing value of the weighting parameter λ. This corresponds to intuition because assigning more weight to color information should yield image segments grouped together based on color rather than on spatial proximity. However the result of segmentation seems best for λ = 0.5 where location and color values of the pixels contribute equally to the resulting matrix factors. This resembles the behavior of a bi-Figure 1. A sample image and its segmentation results obtained from computing cluster centroids using multiview NMF applied to pixel location-and color information. For a smaller λ, more weight will be assigned pixel location information, for a larger λ, more weight will be assigned to pixel color information. With larger weights on location information, small regions of rather homogenous color disappear in the segmentation process. For larger weights on color information, we observe a tendency towards over-segmentation and noisy segment boundaries. For the case where color and location information contribute equally, small regions are preserved and segment boundaries are smoother. lateral filter [START_REF] Tomasi | Bilateral Filtering for Gray and Color Images[END_REF] which also incorporates color-and location information and is known to yield smooth segment boundaries.

Label Prediction and Image Retrieval via Joint Factorization of Image-and Text-Features

This series of experiments aims at exploring whether or not multiview NMF is capable of filling in missing information. We considered a training set of natural images retrieved from the "most interesting" category at flickr. This set of training images contains 10 different (clouds, moonlight, beach, ship, bridge, mountain, forest, city, church, castle) of motives and we considered 300 images per class.

In these experiments, the feature vectors are calculated using local self similarity (SSIM) [START_REF] Shechtman | Matching local selfsimilarities across images and videos[END_REF] feature extraction scheme. The feature vectors are then clustered into a visual vocabulary of k = 750 visual words. For each image in the dataset, a histogram of this vocabulary is created. The individual histograms of all the images in the dataset are then collected in an image-feature matrix F ∈ R k×n .

Textual descriptors for the tag list of the images are created by using the well known Bag of Features [19] approach. Firstly, the most frequent tags in the dataset are collected and the textual vocabulary or the dictionary is generated by filtering the irrelevant tags such as foreign names, flickr group names, and abbreviations. Secondly, all the tag lists corresponding to the respective images are compared with the dictionary and according to the presence (1) or absence (0) of the dictionary words in the tag list of an image, a binary text feature vector is formed. Finally, the feature vectors are stored in a matrix X ∈ R m×n with n being the number of images in the dataset and m = 1000 being the size of the textual dictionary. the matrices W , V and H were initialized to random positive values sampled from a Gaussian and we ran the multiview NMF update algorithm until convergence but at most 100 times, to obtain coupled factorizations (k = 10, λ = 0.5) of the image-and text-feature matrices X and Y , respectively. In the test phase of these experiments, we considered two different settings.

Label Prediction

Given an image that was not part of the training set, we compute its image-feature vector x and solve min x -W h 2 s.t.h 0 for h. Given h, we plug it into y = V h to obtain a corresponding vector y in the text-feature space.

Given y, we search for that column vector y i of the training data matrix Y for which yy i is minimal. We use y to predict a ranked list of tags. To this end, we determine and rank those words in the lexicon that correspond to the 20 basis vectors t i in the original text-by-image space for which the projection y T t i is maximal. The 10 highest ranked tags are selected to be the tag list of the test image. In Figure 1, the retrieved tags for some of the images are shown.

Image Retrieval

In this setting, we queried random words such as bridge, sea, sky individually or as a group to retrieve the best corresponding images. The text feature vector y of the random words are created the same way as training tag lists of the images. We then solve min y -V h 2 s.t.h 0 for h. Given h, we plug it into x = W h to obtain a corresponding vector x in the image-feature space.

Given x, we search for that column vector x i of the training data matrix X for which xx i is minimal. The four most similiar images are shown in Table 2

Conclusion and Future Work

The work presented in this paper aims at the analysis of images for which there is additional information available. We introduced a new model for multiview clustering that extends the idea of non-negative matrix factorization (NMF) towards the joint analysis of different types of features. We cast multiview NMF as a convex combination of individual optimization problems and adopt the well known multiplicative fixed point algorithm for NMF to this case. The approach avoids ad hoc combinations of different types of features and thus stays true to the nature of different descriptors. The individual optimization problems in our multiview NMF formulation are coupled via a common coefficient matrix. Due to this coupling, the resulting basis vectors or cluster centroids allow for inferring one type of descriptor (e.g. image labels) from another type of descriptor (e.g. image features).

In preliminary experiments we validated the applicability of the proposed approach in image segmentation, tag prediction, and tag-based image retrieval. Our first results suggest that multiview clustering can provide a framework for image analysis that applies to different levels of abstraction. Image parts could be identified by combining pixel-color and -location information in the principal manner that is provided by the multiview approach. Information as diverse as color histograms and text-by-image vectors were coupled using our framework and we found it to be capable to predict missing information from what data was available.

Currently, we are conducting more extensive experiments to provide a more quantitative analysis as well as to compare the proposed approach to other multiview methods such as (kernelized) canonical component analysis. In contrast to related methods from the literature, we expect that highly efficient implementations of multiview NMF will be possible. To this end, we are currently adopting techniques such as convex-hull NMF to our model. We will also further explore how multiview NMF relates to LDA and whether it offers an alternative approach to hierarchical latent topic models. Finally, we envision further applications of the proposed method, for instance in the area of hyperspectral imaging.

  Figure 2. The 3 most relevant images retrieved by querying with the word or the group of words below them. The retrieved images tend to be more specific with the increasing number of words used in the queries.

Table 1 .

 1 that correspond to the words below.

	high	blue	water
	travel	water	sky
	cruise	trees	clouds
	holiday	bridge	waves
	morning	reflections	rocks
	cityscape	grass	sea
	daybreak	yellow	ocean
	tower	woods	seascape
	sea	railing	raining
	land	waterscape	waterscape
	water	sky	nightscape
	beach	clouds	blue
	sunrise	blue	stars
	outdoors	holiday	afterdark
	nature	red	sky
	reflection	castle	night
	landscape	bluesky	landscape
	raining	raining	moonlight
	walking	disneyland	yellow
	yellow	middleages	city

Results of automatic image annotation. The taglist corresponds to the first ranked 10 tags retrieved by querying an unknown image.


Environmental Tax and the Distribution of Income among Heterogeneous Workers

Keywords:

This paper analyzes the environmental tax policy issues when labor is heterogeneous. The objective is to assess whether an environmental tax policy could be Pareto improving, when the revenue of the pollution tax is recycled by a change in the labor tax properties. We show that, depending on the heterogeneity characteristics of labor and on the initial structure of the tax system, a policy mix could be designed in order to leave each class of workers unharmed. It consists of an increase in progressivity together with a decrease in the ‡at rate component of the wage tax.

I. Introduction

One of the most publicly debated issues about environmental regulation concerns its e¤ects on the poorest agents and the fairness of such a policy. This question is being raised again because of the relative failure observed until now of the emissions trading system chosen as an instrument for the Greenhouse Gases (GHG) regulation, which explains the emergence of some new arguments in favor of a carbon tax. What would be the inequality consequences of the European Carbon Tax Project or of the Climate-Energy Contribution planned by the "Grenelle de l'environnement" in France? How could these consequences be smoothed or even eliminated? More generally, how is it possible to implement an environmental tax without inducing further inequalities? The economic literature does not contribute much to this debate and our paper aims to …ll the existing gap.

The still open question "Can an environmental tax reform be designed without negatively a¤ecting the economic welfare?" has given rise to a huge literature on the double dividend. [START_REF] Terkla | The E¢ ciency Value of E-uent Tax Revenues[END_REF], [START_REF] Parry | Pollution Taxes and Revenue Recycling[END_REF], or [START_REF] Poterba | Global Warming Policy: A Public Finance Perspective[END_REF] …rst had the intuition that the recycling of the revenue of an environmental tax could reduce or even eliminate the gross cost of its implementation. As governments use the revenues from pollution taxes to decrease other distortionary taxes, environmental taxes may lead to a double dividend, according to Goulder's de…nition, by improving the environmental quality and achieving a less distortionary tax system [START_REF] Goulder | Environmental Taxation and the "Double Dividend": A Reader's Guide[END_REF]). [START_REF] Baumol | The Theory of Environmental Policy[END_REF], [START_REF] Pearce | The Role of Carbon Taxes in Adjusting To Global Warming[END_REF] and [START_REF] Oates | Pollution Charges As a Source of Public Revenues[END_REF] suggested that these e¢ ciency gains could be a powerful argument in favor of environmental taxation. After some prior refutations of the double dividend hypothesis [START_REF] Bovenberg | Environmental Levies and Distortionary Taxation[END_REF] and their followers), [START_REF] Goulder | Environmental Taxation and the "Double Dividend": A Reader's Guide[END_REF] and [START_REF] Ligthart | The Macroeconomic E¤ects of Environmental Taxes : A Closer Look at the Feasibility of Win-Win Outcomes[END_REF] showed that the existence of the double dividend essentially depends on the possibility of transferring the global tax burden from the wage earners to some …xed production factors or to other consumers, thus emphasizing the role of heterogeneity. Following this line, Chiroleu-Assouline and Fodha [2005] and [2006] studied the existence conditions of a long term double dividend, taking into account the distinction between wage earners and retired consumers, by means of OLG models. Surprisingly enough, although this large body of literature has deeply analyzed the double dividend issue, it has neglected the distribution issue of the welfare gain, although it is usually obtained at the expense of some groups of agents. It is well known that the costs and bene…ts of environmental policies are unequally distributed among agents. The rich and the poor seem to assign di¤erent degrees of priority to environmental protection [START_REF] Baumol | The Theory of Environmental Policy[END_REF]), since wealthier individuals would accept paying for a higher level of environmental quality (considered as a normal good). Beyond this fact, distributive elements also matter when we consider how the costs of a policy of environmental protection are likely to be distributed among individuals with di¤ering incomes. Most existing studies are empirical works dealing mainly with the distribution of the bene…ts of the environmental policies among income classes [START_REF] Christainsen | Distributional and Macroeconomic Aspects of Environmental Policy[END_REF], [START_REF] Elliott | Toward a Theory of Statutory Evolution: The Federalization of Environmental Law[END_REF], [START_REF] Harrison | Climate Change, Economic Instruments and Income Distribution[END_REF], [START_REF] Peskin | Environmental Policy and the Distribution of Bene…ts and Costs[END_REF]), and neglecting the cost side of the policies. One can infer from some empirical studies emphasizing the regressivity of the indirect taxes, that any environmental policy is likely also to be regressive. In particular, in the French case, a tax on energy or transport consumption harms the lowest wage households three times more than the highest wage households [START_REF] Ruiz | Le caractère régressif des taxes indirectes : les enseignements d'un modèle de microsimulation[END_REF]). Moreover, the usual recycling of the environmental tax revenues through a decrease in the labor tax rate could also be regressive [START_REF] Metcalf | A Distributional Analysis of Green Tax Reforms[END_REF]). This point is of interest in a world in which inequality has assumed high priority among social issues.

Policies and programs that are perceived as unfair will stand little prospect of passage even if they enhance the prospects for e¢ ciency and sustainability.

Is it then possible to design an environmental policy in order to ensure a non-decreasing welfare for each class of workers? Compared to the standard double dividend literature which pursues two objectives in a second-rank framework -improve the environment by increasing an environmental tax (…rst dividend) and the economic welfare by decreasing another distorsive tax (second dividend) -we add the third objective of Pareto improvement. It corresponds to an acceptability/unanimity criterion for the policy when agents are heterogenous. Three instruments are consequently needed: we argue that the distributive properties of the tax policy could be one of the instruments of fair internalization of the environmental externalities1 . The aim of this paper is therefore to design a balanced environmental tax reform able to correct these regressive properties of taxes and to leave all classes of worker better o¤.

We assume capital-intensive sectors to be polluting sectors. To keep things simple, as we deal only with a one sector aggregate output, we assume that capital input is the main source of pollution (even if it is only a proxy). The key factor is that capital accumulation favors the production of pollution-intensive goods. For industrial pollution, this is consistent with evidence [START_REF] Brown | A North American Free Trade Agreement: Analytical Issues and a Computational Assessment[END_REF], Gale and Mendez [1998], Antweiler et al. [2001]). Indeed, Brown et al. [1992] show that a 10% change in the Mexican capital stock induces a 9% increase in electric utility output. The conclusion they draw is simply that electricity generation is highly capital intensive, and this creates a strong ling between capital accumulation and sulfur dioxide pollution. In the same way, the empirical work of [START_REF] Gale | The Empirical Relationship between Trade, Growth, and the Environment[END_REF] investigate the role that capital abundance may play in predicting cross-country di¤erences in pollution levels. It suggests a strong link beween factor endowments and pollution level. Finally, Antweiler et al. [2001] …nd that a 1% increase in a nation's capital-to-labor ratio, holding scale, income, and other determinants constant, leads to perhaps a one percentage point increase in pollution.

Our model shares the mean features of Chiroleu-Assouline and Fodha [2005] and [2006]. As in Chao and Peck [2000] or [START_REF] Williams | Environmental Tax Interactions when Pollution A¤ects Health or Productivity[END_REF] or [2003], we assume that the degradation of environmental quality has a negative impact on the total productivity of factors. This assumption is justi…ed by the results of an increasing number of empirical studies measuring the health e¤ects of pollution (OECD [2008]) and the impact of the health of workers on labor productivity [START_REF] Bloom | The E¤ects of Health on Economic Growth: A Production Function Approach[END_REF], in a sample consisting of both developing and industrial countries, found that good health, proxied by life expectancy, has a sizable, positive e¤ect on economic growth). Since [START_REF] Ostro | The E¤ects of Air Pollution on Work Loss and Morbidity[END_REF], many papers have emphasized the loss of productivity caused by the health e¤ects of pollution, e.g. [START_REF] Pervin | Societal Costs of Air Pollution-Related Health Hazards: A Review of Methods and Results[END_REF] for air-pollution, and also [START_REF] Bosello | Economy-wide Estimates of the Implications of Climate Change: Human Health[END_REF] or [START_REF] Hübler | Costs of Climate Change -The E¤ects of Rising Temperatures on Health and Productivity[END_REF] for the health e¤ects of climate change. For example, according to [START_REF] Bosello | Economy-wide Estimates of the Implications of Climate Change: Human Health[END_REF], strong heat stress causes a productivity loss of 3% and extreme heat a loss of 12%. Another source of productivity loss originates in the impact of pollution on the quality of natural resources [START_REF] Gollop | From Total Factor to Total Resource Productivity: An Application to Agriculture[END_REF] for the agricultural sector) 2 .

Samakovlis et al. [2005], or

We assume that the production technology is a function of capital and heterogeneous labor. Heterogeneous workers live two periods (young and old) and earn wages correspond- 2 Our model adresses the speci…c issue of pollutions due to industrial, highly capital-intensive sectors, that emit …ne particles, NOx or SOx. These pollutants are harmmful for the health of all agents, and especially of the workers employed by these sectors. As an anonymous referee pointed out, this stylized model can not describe the GHG case, since the dioxyde emissions have no direct e¤ect on the health and productivity of workers and they are due to the use of fossil energies which can be susbsituted by carbon free inputs, such as capital or knowledge.

ing to their skill and consequently to their productivity. The labor tax is a very general one that could be either a progressive or proportional tax. Our demographic assumptions allow us to take into account several income classes; indeed, we consider (i) the heterogeneity characteristics of the labor market (high wages -skilled workers, middle wages, low wages -non skilled workers...), (ii) the heterogeneity of the individual income sources (wages for workers, savings for retirees). The environmental policy consists of increasing the environmental tax on savings, in a second-rank framework. We then characterize the necessary conditions for the obtaining of a double dividend, i.e. an improvement of the environmental quality and an improvement of the welfare when the revenue of the pollution tax is recycled by a change in the labor tax rates. Previous studies show that the obtaining of a double dividend requires economic conditions such that the double dividend hypothesis seems unrealistic. Conversely, we show that the conditions for the obtaining of a double dividend lie in the distributive properties of the labor taxes. Even when the double dividend is not possible, the cost of the pollution regulation can be minimized by a new designing of the progressivity of the labor tax instead of an homogenous cut in the labor tax rates. The results are dependent on the initial tax system.

The paper is organized as follows. Section 2 presents the model and Section 3 presents the welfare analysis of a tax policy. Section 4 gives the speci…cation of the balanced tax reform. In Section 5 we present the environmental e¤ects of the tax reform and Section 6 examines the welfare e¤ects of such a reform. Section 7 presents some empirical illustrations and the last section concludes.

II. The Model

We consider an overlapping generations economy with polluting technology of production i.e. with polluting capital. Heterogeneous workers live two periods (young and old) and earn wages corresponding to their skill and consequently to their productivity 3 . We assume that N t individuals are born in period t and that the population remains constant, so we can normalize N t to unity. Each household is characterized by its labor class (or skill) i. There are I classes of agents in the economy (i and I 2 N), and the size of each class i is equal to q i (with

P I i=1 q i = 1)
. Each agent supplies one unit of labor 4 when she is young and earns a wage w i t ; she divides her labor income between consumption and savings s i t . The sub-optimal initial tax system is composed by an income tax and an environmental tax intended to …nance the public expenditures G: The income (or labor) tax is a very general one that could be either a progressive or proportional tax as the labor tax rate is i 2 [0; 1[. In the second period the household consumes her savings and the interest she earns. The welfare of an individual born at t is measured with the intertemporal separable utility function 5 :

U c iy t ; c io t+1 = (1 ) ln c i;y t + ln c i;o t+1
with c y t denoting the …rst-period consumption of the agent born at t, c o t+1 her secondperiod consumption.

2 ]0; 1[ represents the individual propensity to save. The two instantaneous components of the utility function obviously exhibit the usual properties:

they are increasing in their argument, strictly concave and satisfy the Inada conditions.

The real interest rate is r t+1 . As the capital is polluting, the environmental policy consists in a tax on savings e 2 [0; 1[. The household's budget constraints writes:

1 i t w i t = c iy t + (1 + e t ) s i t c io t+1 = (1 + r t+1 ) s i t (1)
The household's problem is to choose her consumption path to maximize her lifetime utility subject to the intertemporal budget constraint.

This yields the optimal consumption and savings path of the representative household, within the Diamond framework [START_REF] Diamond | National Debt in a Neoclassical Growth Model[END_REF]) with a homothetic utility function:

8 > > > > > > < > > > > > > : c iy t = (1 ) 1 i t w i t c io t+1 = (1 + r t+1 ) (1 + e t ) 1 i t w i t s i t = 1 i t (1 + e t ) w i t ( 2
)
4 Our long term view allows us to assume full employment. Moreover, we focus on the e¢ ciency double dividend (according to [START_REF] Goulder | Environmental Taxation and the "Double Dividend": A Reader's Guide[END_REF]) and not on the employment double dividend. 5 We do not introduce any direct e¤ect of pollution on the household's welfare, but only the indirect one through the consequences on productivity of the degradation of environmental quality (see below). Indeed the direct e¤ect would have no impact on the welfare distribution among heterogenous agents while the indirect one a¤ects the wage gaps, because agents do not di¤er here by their preferences.

The production sector consists of many …rms, each of them being characterized by the same Cobb-Douglas production function F ( ) with constant returns to scale. They use di¤erent kinds of labor (high wages -skilled workers, middle wages, low wages -non skilled workers) and the total productivity of factors A (P t ) is negatively a¤ected by pollution P t because pollution degrades the health of workers or the quality of natural resources (A (P ) = P e with e > 0):

Y t = A (P t ) F (K t ; L i;t ) = A (P t ) K t I Y i=1 L i i;t
where > 0 and i 0 stand for the shares of the input factors in production, + i = 1.

The maximization problem of the representative …rm is (taking the output price as numeraire):

M ax K;fL i g t = A (P t ) K t I Y i=1 L i i;t I X i=1 w i t L i;t (1 + r t ) K t
with t the current net revenue, w i t the real wage rate of class i of workers. The depreciation rate of capital is equal to unity.

Since markets are competitive, capital and labor earn their marginal products: 8 > > < > > :

Y t K t = 1 + r t i Y t L i;t = w i t (3)
This yields, at the equilibrium of the labor markets:

Y t = A (P t ) K t I Y i=1 q i i (4)
The ratio of wages is (using 3):

w i t w 1 t = i 1 q 1 q i
Assumption 1: The di¤ erent labor classes are ranked by growing skills, i.e. by growing wages:

w i t > w i 1 t , i q i > i 1 q i 1
We assume that government spending (G) is entirely …nanced by current taxes. The government's budget constraint states that, at each period, its purchases must be equal to its tax revenues generated by the pollution tax and the labor tax:

I X i=1 q i i t w i t + e t I X i=1 q i s i t = G t (5)
We de…ne a progressivity index of the labor tax, a i = ab(i), such as:

i = 1 + a i = 1 + ab (i)
where 1 is the ‡at component of the tax rate and a is a positive parameter. Assume b(1) = 0 and b(i + 1) > b(i) 8i > 1: We consider the general case for the characteristics of the tax progressivity 6 . For example, the design of progressivity …ts well the characteristics of the French tax system when b(i

+ 2) b(i + 1) < b(i + 1) b(i):
This yields

I X i=1 q i i w i = I X i=1 i i Y = Y 1 (1 ) + a I X i=1 b(i) i !
Let us de…ne B I as the sum

I P i=1 b (i) i .
B I > 0 is constant for any given I and progressivity characteristics fb(i)g 1 i I .

The pollution ‡ow is due to the capital stock, and we assume that the welfare is a¤ected by the stock of pollution only through its e¤ect on the total factor productivity.

The dynamics of pollution are described by the following equation:

P t = (1 h) P t 1 + K t 1 ( 6
)
where h is the constant rate of natural absorption of pollution (0 < h < 1) and > 0 stands for the emission rate of pollutants.

The equilibrium condition of the capital market, meaning that the capital stock in period t+1 is the amount saved by young individuals in period t, is obtained by substituting the zero-pro…t condition, the government's budget constraint (eq. 5) and the household's budget constraints (eq. 1) into the equilibrium of the output good market. It writes:

K t = I X i=1 q i s i t 1
By substituting eq. 2:

K t = 1 + e t 1 q 1 1 1 1 t 1 (1 ) aB I w 1 t 1 (7)
and using eq. 3, we …nally obtain:

K t = 1 + e t 1 1 1 t 1 (1 ) aB I Y t 1 (8)
The steady-state A steady state equilibrium is an allocation where capital and pollution are stationary, i.e. K and P are such that (using eq. 6 and by substitution of eq. 4 in eq. 8):

K = " 1 + e XA (P ) I Y i=1 q i i # 1 1 (9) P = h K (10)
with X = 1 1 (1 ) aB I , constant for a given tax system.

It is easy to verify that X is always strictly positive, which ensures the existence of a positive steady-state equilibrium. X > 0 reduces to b ({) < 1 1 a where { is the particular class of worker 7 such that b ({) =

I 1 b(i) i I 1 i
. Hence b ({) is the average of the coe¢ cients of progressivity, weighted by the shares of each class of worker in the output. X > 0 is then equivalent to { = 1 + ab ({) < 1, meaning that the revenue of this particular class is not totally con…scated. This is always veri…ed since we assumed that i < 1 8i.

The following proposition characterizes the steady state of this model. 7 More precisely, b ({)

I 1 b(i) i I 1 i and b ({ 1) < I 1 b(i) i I 1 i
, with { I.

Proposition 1 For given initial conditions fK 0 ; K 1 ; P 1 g and for a given triplet of tax policy instruments, e ; 1 ; a , the equilibrium conditions can be summed up by eq.

(2) ; (3) ; (5) ; (6) ; (8). Thus, there exists a locally stable steady-state de…ned by (K ; P ) such that:

K = 1 + e X 1 1 ( e)
"

h e I Y i=1 q i i # 1 1 ( e) (11)
P = h 1 + e X 1 1 ( e)
"

h e I Y i=1 q i i # 1 1 ( e)
Proof. See above for the de…nition of the capital market equilibrium. Moreover, (i) the utility function is homothetic, (ii) the utility and the production functions verify the Inada conditions, (iii) Using Proposition 1 and eq. 9, one obtains the output value at the steady-state equilibrium:

F 00 K < 0, (iv) lim K!0 (F 0 K :K=F (K)) 2 [0; 1[. Conditions (i) (iv)
Y = 1 + e X e 1 ( e)
"

h e I Y i=1 q i i # 1 1 ( e) (12)
Assumption 2: e > 0:

It means that the direct contribution of capital to output ( ) is greater than the indirect one through pollution and productivity (e)8 .

In the rest of the paper, we will only consider the consequences of variations of the tax rates on the steady-state equilibrium de…ned in Proposition 1. The existence of the new steady-state is ensured if and only if X remains positive, i.e. { + d { < 1. Notice that X may remain positive even if some high wages classes bear a tax rate equal or greater than unity. A su¢ cient condition to ensure the positivity of X after the tax reform is that dX > 0, i.e. that the variations of the labor tax components verify da d 1 < 1 B I , 8d e > 0.

III. Welfare Analysis

Like Chiroleu-Assouline and Fodha [2006], we examine here the welfare e¤ects of the tax change for a generation during its life-cycle, once the …nal steady-state equilibrium is reached. In this section, the welfare issue is thus a long term one.

One can measure the welfare e¤ects of small tax changes by the marginal excess burden.

It corresponds to the additional income that needs to be provided to the representative household to keep her utility at its initial level: this is the compensatory income variation, denoted dR. It stands for the excess welfare loss of the consumers over and above the tax revenues collected by the government and can be interpreted as the hidden costs of …nancing public spending: a positive value for the marginal excess burden indicates a loss in welfare after the tax reform.

Let us determine the compensatory income variation which, after the tax reform, would leave the level of life-cycle utility unchanged (dU = 0): @U @c iy dc iy + @U @c io dc io = 0

() (1 + r) (1 + e ) dc iy + dc io = 0
The intertemporal budget constraint of household i writes:

1 i w i = c iy + (1 + e ) (1 + r) c io
Remember that

c io = (1 ) (1 + r) (1 + e ) c iy
We use the …rst-order conditions of the representative household's program and the de…nition of the compensatory income variation dR:

1 i dw i w i d 1 + b (i) da + dR i = dĉ iy ĉio (1 + r) 2 dr + ĉio (1 + r) d e + (1 + e ) 1 + r dĉ io = (1 + e ) ĉio (1 + r) 2 dr + ĉio (1 + r) d e
this leads to:

dR i = 1 i dw i + w i d 1 + b (i) da (1 + e ) ĉio (1 + r) 2 dr + ĉio (1 + r) d e (13)
Unlike Bovenberg and de Moiij [1994] and the greater part of the literature on this subject, it is here impossible to distinguish an environmental component and a nonenvironmental one, because pollution and production a¤ect each other. In this paper, we are thus constrained to depart from the usual de…nition of the double dividend [START_REF] Goulder | Environmental Taxation and the "Double Dividend": A Reader's Guide[END_REF]) because of this non-separability: a double dividend will be characterized by the simultaneous decrease of pollution (which stands for the usual …rst dividend) and increase of economic welfare (which depends here also on the pollution level).

The increase of welfare is veri…ed for any variations of consumptions such that dc y =c y +

(1 )dc o =c o > 0. This property plays an important role in our results. The …rst dividend requires a decrease of the capital stock inducing a decrease of the output, but this does not prevent the obtaining of the economic dividend. This originates in the opposite impacts of dK on w and r: dw=dK > 0 and dr=dK < 0.

Proposition 2 The increase of the environmental tax is regressive i.e. it more heavily penalises the lowest wage-earners.

Proof. Let us compute the steady-state value of the compensatory income variation for agent i:

dR i i q i Y = Z 1 i (1 + e ) d e (14)
+ (Z ) (1 ) X + 1 (1 i ) 1 i d 1 + (Z ) B I X + b (i) (1 i ) 1 i da with Z = e
1 ( e) > 0. Notice that from the optimality conditions of the …rm (eq. 3) and the equilibrium condition of the labor market, we have i q i Y = w i . Hence, when d 1 = da = 0, one obtains, for any environmental tax increase d e > 0 :

dR i w i = Z 1 i (1 + e ) d e > 0
then, as Z > 0, dR i > 0 for every i. But dR i =w i is smaller when i is higher (i.e. i is also high) and when w i is higher. When the environmental policy harms the consumers, the relative burden is greater for the less skilled9 .

The consequences of the increase of the environmental tax e are clearly harmful for all the agents as it implies a decrease of the second-period consumption without any positive e¤ect on the …rst-period one. Concerning the labor tax, its impact is ambiguous, for the ‡at component as well as for the progressivity index. This is because the decrease of the income of the agents could be compensated by an increase of the interest rate raising the consumption of the old agents. This latter e¤ect increases with .

IV. The Speci…cation of the Balanced Tax Reform

We assume an exogenous increase of the environmental tax rate, imposed by the government in order to control pollution. The amount of government spending is assumed ex post invariant, and the tax policy has to maintain the amount of the tax revenues constant 10 . This increase d e of the environmental tax rate can be accompanied by a variation of the labor tax rates d i by two potential means: a homogenous variation of all labor tax rates through the ‡at rate component (d 1 ) or a variation in the progressivity of the labor tax (da i through a variation of a). At the steady-state equilibrium, the government's budget constraint is:

1 (1 ) + aB I Y + e K = G
The link between the variations of the pollution tax and of the characteristics of the labor tax is obtained through the di¤erentiation of this constraint, which is quite direct (using eq. 12).

Any balanced tax reform is then characterized by the following relationship between d e , d 1 and da (with dG = dq i = d i = 0):

(1 )

d 1 + B I da = ( + ) (1 + ) X d e 1 + e (15)
where

= 1 (1 ) + aB I + K Y 1 ( e)
e Z X 0 (with X > 0), constant.

1 0 As noticed previously, this tax reform does not constitute the optimal one. Seeking for optimality would require to take into account the intergenerational externality on capital accumulation (dynamic e¢ ciency) and the environmental externality. These two externalities would be corrected by the use of our two taxes. Concerning the distribution of welfare among the agents, we have to consider intergenerational and intragenerational characteristics. Intergenerational optimality considerations turn to the very general debate about the choice of the discount rate. But in our very simple model, the heterogeneity of workers is not an externality and would not deserve any correction. Intragenerational welfare would be the sum of the welfare of each class weighted by its size.

We will study two polar cases for balancing this environmental tax reform: uniform variation of all labor tax rates, with invariant progressivity (da = 0):

d 1 = ( + ) (1 + ) X (1
)

d e 1 + e = d e (16)
variation of the progressivity, with invariant ‡at rate component of the tax rate

(d 1 = 0): da = ( + ) (1 + ) X B I d e 1 + e = d e (17)
Proposition 3 The sign of the balanced tax reform multipliers and is a priori undetermined and depends on the initial tax rates and on the values of the various elasticities.

Proof. (i) The numerator measures the e¤ect of the change in pollution tax rate on its revenue. There are both a value e¤ ect (the tax revenue increases with the tax rate, for unchanged stock of capital), summed up by , and a tax base e¤ ect (the stock of capital decreases, and so does the output; as the tax rate increases, so the tax base erodes), conditioned by , which work in opposite ways. As a result, this term might be positive

(if > ) or negative.
(ii) The denominator measures the e¤ect of the change in labor tax rate on its revenue.

There are also both a value e¤ ect proportional to (the tax revenue increases with the tax rate, for unchanged wages) and a tax base e¤ ect (the wages decrease as the tax rate increases, thus the tax base erodes) which work in opposite ways. As a result, this term too might be positive or negative (if > 1).

As the signs of the numerator and of the denominator are undetermined, the sign of the necessary change in the labor tax components is also undetermined.

Assumption 3:

We only consider the La¤ er-e¢ ciency case, where the sign of the balanced tax reform multipliers and is positive.

Under eq. 15, Assumption 3 implies

+ 1 + > 0 , > 1 or < .

V. The Environmental E¤ects of the Tax Reform

As steady-state pollution depends only on steady-state capital, we …nd straightforward that the …rst dividend (i:e: a decrease of pollution) is reached if dP = h dK < 0: This condition rewrites (see Proposition 1):

dK K = 1 1 ( e) d e 1 + e + (1 ) d 1 + B I da X < 0
As 1 ( e) > 0, this condition reduces to

d e 1 + e + (1 ) d 1 + B I da X > 0:
We then have the two cases corresponding to the two alternative policies:

invariant progressivity (da = 0) :

d e 1 + e > (1
) d e X variation of the progressivity d 1 = 0 :

d e 1 + e > B I d e X :
We obtain a …rst dividend if the multipliers are not too high:

8 > > < > > : < X (1 + e ) (1 ) < X (1 + e ) B I () + 1 + < 1 (18)
This means that the decrease of the labor tax (either of the ‡at rate component or of the progressivity index) should be not too great in order not to jeopardize the environmental bene…t of the reform.

Assumption 4:

We only consider the case where the tax policy decreases pollution.

Proposition 4 Under Assumptions (1-4), the tax policy is e¢ cient and reduces pollution if and only if the initial environmental tax rate is not too high: e < 1=Z.

Proof. We have to examine the two cases which verify the tax e¢ ciency assumption (Ass.

3).

If (i) : 1 < ; then + 1 + < 1 , > 1 : which con ‡icts with the basic assumption on : The environmental bene…t can never be achieved.

If (ii) :

< ; then + 1 + < 1 , < 1 : always veri…ed.
Finally, such that to reach the …rst dividend under the tax e¢ ciency assumption, one needs to impose < ; which is a necessary and su¢ cient condition. This condition is met if and only if: Proof. Let us compare the compensatory income variations associated with each policy.

1 (1 ) + aB I + 1+ e X 1 ( e) e Z X < , 1 (1 )
First case (dR i ( 1 ) ): uniform variation of all labor tax rates (with invariant progressivity): da = 0 and

d 1 = d e dR i ( 1 ) i q i Y = Z 1 i (1 + e ) d e (Z ) X (1 ) + 1 (1 i ) 1 i d e
Second case (dR i (a) ): variation of the progressivity (with invariant ‡at rate component):

d 1 = 0 , da i = b (i) da and da = d e dR i (a) i q i Y = Z 1 i (1 + e ) d e (Z ) X B I + b (i) (1 i ) 1 i d e i %

Agents prefer da < 0

Agents prefer1 < 0 The signs of these compensatory income variations depend on the values of the characteristics of the economy ( , e, b(i), i ) and of the initial tax rates ( e , i , a).

0 i b(i) i
We now compare the relative e¤ects on welfare of the two tax policies. Consequently, we have to determine the sign of i = dR i

( 1 ) dR i (a) for the workers of class i.

i i q i Y = (Z ) X (1 ) 1 i + 1 + (Z ) X B I 1 i + b (i) T 0
Using eq. 16 and 17 ( = 1 B I ), we …nd that this sign is equal to the sign of:

b (i) 1 B I 1 = b (i) i b (i) i 1
Therefore the particular "average" workers' class { de…ned above is indi¤erent between both policies ( { = 0). Each class of workers of a higher skill would prefer a decrease in the progressivity index, da < 0, and, conversely, each less skilled class of workers would prefer a decrease in the ‡at rate component of the wage tax, d 1 < 0 (Figure 1 corresponding to the French case of a progressive labor tax). 11

The result above suggests that, in the case where some workers'classes would be suffering from a deterioration of their welfare after the above tax reforms, an appropriate policy mix could be designed in order to leave each workers'class unharmed by the environmental tax reform : it would consist in an increase of the progressivity index together with a decrease of the ‡at rate component of the wage tax rate.

Assume that the balanced tax reform is de…ned by d 1 < 0 and da = d 1 with > 0 hence da > 0. Such a compensation for the increase in the environmental tax rate will imply a greater decrease of 1 than above because the degree of progressivity is raised.

Proposition 6 For economies where a decrease in the minimal wage rate does not su¢ ce to re-establish the welfare of some class of worker, one can choose in order to ensure that all classes will be better o¤ with the environmental tax reform.

Proof. Let us specify the link implied by such a policy between the increase in the environmental tax rate and the decrease in the ‡at component of the wage tax.

d e = 1 d 1 1 da = 1 d 1 = 1 1 B I d 1
To ensure the existence of the new steady-state equilibrium, we have shown above (Ass.

2) that we need da d 1 = < 1 B I . This condition also ensures the decrease of the ‡at rate component of the labor tax.

The compensatory income variation of the balanced mix policy is:

dR i ( ) i q i Y d e = 1 B I " 1 dR i 1 i q i Y d e + " (Z ) B I X 1 i + b (i) B I 1 1 Z 1 i (1 + e ) ##
with Z > 0; X > 0; constants.

-dR i ( ) S 0; 8i = 1; :::I:

-dR i ( ) is not monotonous in b (i), nor in , 8i.
-If dR i ( 1 ) < 0; 8i = 1; :::I, all classes will be better o¤ even with = 0:

If 9i; dR i ( 1 ) > 0; as the function dR i ( ) is bound for each , there is one { which maximizes it (see Fig. 2 One can choose > 0 in order to nullify dR { ( ) : it ensures that all classes will be better o¤ (their compensatory income variations are all negative or null).

= Z 1 1 ab ( {) (1 + e ) h 1 1 ab ( {) (Z ) (1 ) X + 1 i " B I 1 1 Z 1 1 ab ( {) (1 + e ) (Z ) B I X (1 1 ab ( {)) b ( {)
# This policy mix aims to reduce the pre-existing distortions of the tax system, in a second-best world, but not to reach an optimal outcome. Therefore comes from our Pareto-improving criterion but not from any optimality criterion.

It may be useful now to sum up the main mechanisms at work. The primary consequence of the environmental policy is the decrease of the pollution allowed by the decrease of the capital stock which also implies a reduction of the output. The cost of this environmental bene…t consists in a fall in the gross wages of all agents. Nevertheless, for some groups of agents, this decline of the gross wages may be slowed by a decrease of their labor tax rate (for any group i such that b(i) < 1= ). The harmful e¤ect of the fall in net wages on the total life-cycle consumption can be somehow counterbalanced by the drop in the capital stock, implying a rise of the interest rate which bene…ts the second-period consumption.

even the …rst class would bene…t from the environmental tax policy. The increase in the welfare of the upper classes will be reduced, but is still high, comparatively to the lower classes.

For some intermediate levels of the labor tax rate ( = 0:5 and 1), it is impossible to de…ne a policy compatible with the Pareto-improvement criterion. This con…guration could be explained by the initial characteristics of the tax system. For these cases, the initial values are e¢ cient, close to their optimal level. Hence, some agents would su¤er a big decrease in their revenues that could not be compensated by the small gains of the agents who bene…t from the tax policy. When the progressivity part of the labor tax rate is too high ( = 3), the tax system is obviously characterized by the ine¢ ciency properties. Finally, for some reasonable values of the progressivity index = 1:5 (resp. 2), close to the French economic data, the alternative policies are unambiguously harmful for all the classes of workers. But, by …xing to 0:26 (resp. 0:23), the environmental policy would bene…t all the agents. The environmental tax policy is then Pareto-improving and acceptable.

VIII. Conclusion

In this paper, we have shown that a budget-neutral environmental tax reform may result in a double dividend (de…ned as a decrease in pollution and an increase in the global economic welfare), even when the economy is characterized by heterogenous agents (old and young) and many worker classes (heterogenous labor). We have also emphasized that the conditions for the obtaining of a double dividend depend on the distributive properties of the labor taxes. Hence, we have shown that (i) an increase of the environmental tax deteriorates the welfare of all and is regressive, (ii) the low paid workers prefer an environmental tax reform balanced by a decrease in the ‡at rate component but the high paid workers prefer a decrease of the progressivity.

Even when the double dividend is not possible, the cost of the pollution regulation can be minimized by a new designing of the progressivity of the labor tax instead of a homogenous cut in the labor tax rates. We conclude that the distributive properties of the tax policy could be one of the instruments of fair internalization of the intergenerational externalities.

To a certain extent, our paper highlights the gap between economic e¢ ciency and vertical equity and illustrates the problem of the aggregation of positive and negative as: 2 + (P (h(2 h) (1 h) 1) K(1 h+ ))+ P (1 h(2 h))+ K( (1 h)+e) P (1 h)+ K i

We know that P =K = =h: So, we have the following properties P (0) = 0; P (1) =

(1 +e)

h (1 h)+ > 0 and P ( 1) = h ((2 h)(1+ )+he)

h (1 h)+ > 0:

For su¢ cient conditions to local stability, we need that the modulus of each eigenvalue to be lower than one. We can decompose the characteristic polynomial into two parts:

P ( ) = P 1 ( ) :P 2 ( ), with P 1 ( ) = 2 . Basically, the characteristic polynomial has two eigenvalues equal to 1 = 2 = 0: Then we have to study the second part of P ( ) i.e. P 2 ( ) = 2 + (P (h(2 h) (1 h) 1) K(1 h+ ))+ P (1 h(2 h))+ K( (1 h)+e) P (1 h)+ K :

We already know that P (1) > 0 and P ( 1) > 0; which induce that P 2 (1) > 0 and P 2 ( 1) > 0: Regarding P 2 (0) = P (1 h(2 h))+ K( (1 h)+e) P (1 h)+ K = h (1 h(2 h))+ ( (1 h)+e)

h (1 h)+

=

(1 h) + he > 0: We now need P 2 (0) < 1 which turns to (1 h) + he < 1; that is always true under Assumption 2. Hence, the steady state equilibrium (K ; P ) is locally stable.

  ensure the existence of the steady-state equilibrium (see de la Croix and Michel [2002] and Konishi and Perera-Tallo [1997]). Concerning local stability, we show that the steady state associated with the dynamic system, (6; 8), is a sink (see Appendix).

Figure 1 :

 1 Figure 1: An Illustration of Workers Class Threshold in the French Case

Figure 2 :

 2 Figure 2: An Illustration of the Policy-Mix Criterion

  t+1 = P t andz t+1 = K t : The characteristic polynomial of is P ( ) = 2

  h

  + aB I + The environmental tax reform preferred by the workers of the lowest classes consists of a decrease in the ‡at rate component of the tax rate 1 (ii) but the one preferred by the workers of the highest classes consists of a decrease in the progressivity of the wage tax rate.

	which implies
	e 1+ e X 1 ( e) < X Z
	, (	e) < 1 1+ e
	,	e < 1 ( e) e
	QED.
	VI. The Welfare E¤ects of the Tax Reform
	As each policy does not a¤ect all classes equally, one can wonder which one will be
	preferred by each class of worker.
	Proposition 5 (i)
		e 1+ e X 1 ( e) < X Z

In line with the double dividend literature, we are not seeking for the optimal tax system, which would require to use several instruments in order to internalize the externalities. Nevertheless, according to our Pareto-improving criterion, our results give conditions such that the economy gets closer to its Pareto optimal equilibrium.

In our very stylized model, the heterogeneity of workers derives, by assumption, from their skills. In our case, because the utility function is logarithmic, this model is equivalent to a model with endogenous work supply and heterogenous desutilities of work. We thank an anonymous referee to have drawn our attention to this equivalence.

Unlike the usual speci…cation of the tax progressivity in the literature, we assume here that i depends on the class i of the agent and do not depend directly on the wage w i : As assumption 1 implies a ranking between i and w i , our de…nition of the tax progressivity is nearly, but not exactly equivalent to the usual one, since it relies here on an ordinal concept rather than on a cardinal one, as usual.

The reverse ( < e) would simplify the analysis but this case would be of negligible economic interest because any rise in e would result in a decrease of K allowing an increase of Y .

This result does not rely on the usual channels causing the regressivity of an environmental tax. Indeed, as preferences are identical among agents, their savings rate and their budget share for polluting goods are the same. Our regressivity property results from the progressivity characteristics of the labor tax.

1 In order to simplify the graphics (especially Figure3), we have chosen to plot …gures as if all functions were continuously de…ned over a great number of workers classes.

The heterogeneity characteristics of the economy and the progressivity of the labor tax allow the government to obtain more revenues from the environmental policy. The decrease of the ‡at rate component can be greater without compromising the environmental bene…t.

Moreover, it enables the redistribution of welfare between the agents in order to ful…ll our Pareto-improvement criterion.

Observe that even if the productivity is insensitive to pollution (e = 0) and/or if the pollution would a¤ect the household's welfare in an additive way, our results would still be robust and contribute to the standard literature.

VII. The French Case: A Simple Empirical Illustration

In order to illustrate the economic and welfare consequences of the environmental policy-mix, we choose some realistic values for the parameters of our theoretical model. Then, under these parameters'con…guration, we compute the variations of compensatory income with respect to the class of agents i: Actually, we want to test the sensitivity of with respect to the initial level of the labor tax rates. In more detail, we compute the values of the dR i ( 1 ) ; dR i (a) and dR i ( ) for di¤erent realistic values of i :

The initial tax system is supposed to be de…ned by the following values: e = 0:01; cases: = (0:01; 0:1; 0:5; 1; 1:5; 2; 3), see Fig. 3. In particular, the French case corresponds to = 1:5:

The share of capital in total output is very common for developed countries ( = 0:35) and we assume no particular preferences for the future ( = 0:5).

The environmental parameters are not easy to de…ne. We …xed these values to some reasonable levels for this literature: e = 0:05; = 10; h = 0:3. Nevertheless, our results are robust to any changes in the values of these parameters. q i 0:1 0:1 0:1 0:1 0:1 0:1 0:1 0:1 0:1 0:1 i 0:01 0:02 0:03 0:04 0:05 0:06 0:07 0:1 0:12 0:15

Table I. presents the labor market characteristics. We assume an equi-distribution of the workers among the classes, such that q i = 0:1; 8i and I = 10; hence each i represents a decile of the corresponding class. We calibrated the shares in total output for workers of class i ( i ) such that the ratios of wages …t the French data. This ratio is equal to 15 for the highest class. Prog. Case: 1 2 3 4 5 6 7 0:01 0:1 0:5 1 1:5 2 3 0:1 0:11 ? ? 0:26 0:23 Tax ine¢ ciency

We have numerically computed the value of required by the policy mix according to the di¤erent cases of progressivity (see Table II.). Concerning the low progressivity cases = 0:01 (resp. 0:1); the consequences on welfare of the alternative policies, measured by dR i

( 1 ) and dR i (a) , show that the …rst class would be damaged by the policies, the second class is una¤ected, and all the other classes would be better o¤. These policies are not Pareto improving. But, if the government can combine variations of the two parts of the labor tax rate and design a policy-mix system, he has to …x to 0:1 (resp. 0:11). Therefore, compensatory variations: the usual method of aggregation gives a higher weight to the wealthiest classes and introduces a bias when assessing the desirability or the acceptability of any environmental tax reform.

Aknowledgements

This paper bene…ts from the …nancial support of the French National Research Agency (Grant ANR-09-BLAN-0350-01). We would like to thank two anonymous referees for precious suggestions. We would like to thank also participants to the 2008 AFSE Thematic To prove the local stability part of the proposition, we linearize the dynamic system (6; 8) around the steady state. The linear approximation can be written in matrix form