Vincenzo Verardi
email: vverardi@fundp.ac.be

Nicolas Debarsy
email: ndebarsy@fundp.ac.be

n-consistent semiparametric regression estimator in Stata

Keywords: semipar, Semiparametric estimation JEL Classification: C14, C21

 double residual semiparametric regression estimator and Hardle and Mammen's (1993) specification test implementation in Stata. Some simple simulations illustrate how this newly coded estimator outperforms the already available semiparametric plreg command.

Introduction

The objective of this paper is to present the implementation in Stata of [START_REF] Robinson | Root-N-consistent semiparametric regression[END_REF] double residual semiparametric regression estimator. Also, to check if the nonparametric part of the relation may be approximated by a polynomial functional form, we introduce [START_REF] Hardle | Comparing nonparametric versus parametric regression fits[END_REF] specification test as an option in the programmed estimator. We also briefly describe this test.

The structure of the paper is the following: in Section 2, [START_REF] Robinson | Root-N-consistent semiparametric regression[END_REF] semiparametric regression estimator and [START_REF] Hardle | Comparing nonparametric versus parametric regression fits[END_REF] specification test are described. In Section 3, the implemented Stata command (semipar) is presented. Some simple simulations assessing the performance of the estimator and of the test are performed in Section 4. In Section 5, we illustrate the use of the semipar command with an empirical application. Section 6 concludes.

2 Estimation method 2.1 [START_REF] Robinson | Root-N-consistent semiparametric regression[END_REF]

semiparametric regression estimator

Consider a general model of the type

y i = θ 0 + x i θ + f (z i ) + ε i , i = 1, ..., N (1)
where y i is the value taken by the dependent variable for individual i, x i is the row vector of characteristics for individual i, θ 0 is a constant term and ε i is the disturbance assumed to have zero mean and constant variance σ 2 ε . Variable z is an explanatory variable that enters the equation nonlinearly according to a non-binding function f . This model can be estimated using [START_REF] Robinson | Root-N-consistent semiparametric regression[END_REF] double residual methodology that starts by applying a conditional expectation to both sides of (1). This leads to

E (y i |z i ) = θ 0 + E (x i |z i ) θ + f (z i ) i = 1, ..., N (2)
By subtracting (2) from (1), we have

y i -E (y i |z i ) = (x i -E (x i |z i )) θ + ε i i = 1, ..., N (3)
If the conditional expectations are known, parameter vector θ can easily be estimated by fitting (3) by ordinary least squares. If they are unknown, they have to be estimated by calling on some consistent estimators

y i = m y (z i ) + ε 1i and x ki = m x k (z i ) + ε 2ki ,
where k = 1, ..., K is the index of the explanatory variables entering the model parametrically. [START_REF] Robinson | Root-N-consistent semiparametric regression[END_REF] double residual estimator is hence the OLS estimation of model

y i -my (z i ) = (x i -mx (z i )) θ + ε i i = 1, ..., N (4)
where x imx (z i ) is the row-vector of the differences between each explanatory variable x ki and the fitted conditional expectation of x ki given z i .

2

The estimated coefficients vector is therefore

θ = i (x i -mx (z i )) ′ (x i -mx (z i )) -1 i (x i -mx (z i )) ′ (y i -my (z i ))
(5) with variance (if errors are i.i.d)

V ar θ = σ 2 ε i (x i -mx (z i )) ′ (x i -mx (z i )) -1 (6)
where σ 2 ε is the variance of the error term. If errors are non i.i.d., standard sandwich and cluster variance formulas can be used.

Having estimated parameter vector θ, it is now possible to fit the nonlinear relation between z i and y i by simply estimating equation ( 7) presented below nonparametrically.

y i -x i θ = θ 0 + f (z i ) + ε i , i = 1, ..., N (7)

Hardle and Mammen's (1993) test

It is sometimes suggested that nonparametric functions may be approximated by some parametric polynomial alternative. To test for the appropriateness of such an approximation, [START_REF] Hardle | Comparing nonparametric versus parametric regression fits[END_REF] propose a statistic which compares the nonparametric and parametric regression fits using squared deviations between them. The test-statistic is:

T n = N √ h N i=1 f (z i ) -f (z i , θ) 2 π(•) (8)
where f (z i ) is the nonparametric function estimated in (7), f (z i , θ) is an estimated parametric function and h is the bandwidth used. π(•) is a weight function. To obtain critical values for the test, [START_REF] Hardle | Comparing nonparametric versus parametric regression fits[END_REF] suggest to call on simulated values obtained by wild bootstrap. Obviously, an absence of rejection of the null (i.e. "accepting" the parametric model) means that the polynomial adjustment is at least of the degree that has been tested.

We implemented this estimator in Stata under the command semipar which is described below.

The semipar command

The semipar command fits Robinson's double residual estimator in the case of a unique variable entering the model nonparametrically. The default kernel regression used for all stages is a gaussian kernel weighted local polynomial fit. This kernel is of order 2. The optimal bandwidth used minimizes the conditional weighted mean integrated squared error.

The general syntax for the command is: The first option, nonpar, is compulsory and necessary to declare which variable enters the model nonparametrically. All other choices are optional. The first one, generate, reproduces the "nonparametrically" fitted dependent variable. The user chooses the name of this new variable by defining it in parentheses. Similarly, the partial option is needed to generate a new variable that contains the parametric residuals (i.e. the left-hand side of equation ( 7)). The option nograph will be used if the user does not want the graph of the nonparametric fit of the variable set in nonpar to appear. The ci option allows to visualize the confidence interval around the nonparametric fit, while the level option sets the level of confidence for inference (by default set to 95%). The options title, ytitle and xtitle are used to indicate respectively the title and the labels of the axes of the graph illustrating the nonparametric relation between the dependent variable and the variable defined in the nonpar option. The robust and cluster options call for standard errors of the estimated parameters that are respectively resistant to heteroskedasticity and clustered errors. The test option implements [START_REF] Hardle | Comparing nonparametric versus parametric regression fits[END_REF] statistic to test whether the nonparametric fit could be approximated by a polynomial fit, the order of which must be set by the user. For the sake of clarity, we rescaled the statistic in such a way that it can be compared with the quantile of a Normal distribution. Note however that the test is not normally distributed. The nsim option defines the number of bootstrap replicates used to get inference. Its default value is set to 20. The weight option allows to introduce a weighting vector for the test (i.e. π(•) in equation ( 8)). By default, this weighting vector is set to ι N /N with ι N a unit vector of dimension N . Finally, the degree option allows the user to specify the degree of the local polynomial fit used to nonparametrically estimate the regressions. By default, it is set to 1.

semipar varlist if in ,
To assess the performance of the programmed estimator, in the next section we present some simple simulations in which we compare this estimator with the already available plreg command. The latter implements [START_REF] Yatchew | Nonparametric regression techniques in economics[END_REF] difference estimator where the nonparametric part in (1) is partialled out by differencing rather than by removing the conditional expectations. Since the highest efficiency of [START_REF] Yatchew | Nonparametric regression techniques in economics[END_REF] estimator is attained by a differencing of order 10, we will use this differencing order as a benchmark.

Simulations

The simulation setup is the following. To begin, we generate (for a sample of 300 observations) two explanatory variables x 2 and x 3 from two independent N (0, 1). An additional random variable x 1 is generated from a discrete Uniform distribution on [-10, 10]. This sample design remains unchanged for all simulations. Then, for each replication, we generate an error term e from a standard normal and create variable y according to DGP y = x 1 + x 2 1 + x 2 + 3 + e. We run the semipar and plreg estimators for each replicate and Table 1 reports both the bias and MSE of the coefficients associated with x 2 and x 3 . We carry out 1000 simulations. The variable that enters the equation nonparametrically is generated from a discrete Uniform distribution on purpose to illustrate the fragility of plreg with respect to this kind of data. [START_REF] Robinson | Root-N-consistent semiparametric regression[END_REF] estimator, that is based on partialling out the nonparametric part removing conditional expectations rather than by differencing, behaves much better.

In this setup, [START_REF] Robinson | Root-N-consistent semiparametric regression[END_REF] estimator leads to smaller biases than Yatchew's (1998) differencing estimator. From equation ( 7), this also implies that the nonparametric fit is better estimated by semipar than by plreg.

To illustrate the fitting performance of the proposed estimation procedure, we generate four samples according to the following DGPs:

a) y = x 2 + x 3 + e b) y = x 1 + x 2 + x 3 + e c) y = x 1 + x 2 1 + x 2 + x 3 + e d) y = x 1 -x 2 1 -x 3 1 + x 2 + x 3 + e.
In Figure ( 1), we present the scatter plots, the non-parametric fit (thick plain line) and the true DGP (red dashed line) related to the four DGPs described above. As expected, the results are unambiguous.

In the absence of any relation between x 1 and y (panel a), no clear pattern emerges and the non-parametric curve lies close to the horizontal line (the true DGP). In the three other cases (panels b, c and d), the nonparametric estimation of the relation matches the true functional form quite well.

As mentioned in the previous section, the T n statistic assesses the adequacy of a polynomial adjustment compared to a nonparametric fit. Table 2 presents the performance of the test for the DGPs described above. The rows indicate the order of the generated polynomial while the columns specify the order of the polynomial that has been tested. Thus, the diagonal (and the upper triangle) elements are the simulated sizes of the test while elements below the main diagonal are some measure of power. To construct this table we replicated the DGPs 1000 times. Each time a new error term is randomly drawn and a new dependent variable is generated (the design space remains unchanged). Inference for the test is based on 100 bootstrap replications. We observe that the test has good rejection rates when the order of the polynomial adjustment tested is lower than the generated one.

Besides, the size of the test (whose theoretical value is set at 5%), is very close to its nominal value.

Example

To illustrate the usefulness of this semiparametric model in empirical applications, we call on a dataset used by [START_REF] Wooldridge | Econometric analysis of cross section and panel data[END_REF] that studies the effects of an incinerator location on housing prices. The data are for houses that were sold during the year 1981 in North Andover, MA; 1981 was the year construction began on a local garbage incinerator. The dependent variable is the log of price of houses (lprice) and the variable of interest is the distance from the house to the incinerator measured in feet and expressed in logs (ldist). To control for confounding effects, the author suggest to include the log of interstate distance (linst), the log of the square footage of the house (larea), the log of the lot size in square feet (lland), the number of rooms (rooms), the number of bathrooms (baths), and the age of the house (age) as additional covariate. However, he also asserts that the effect of the log of the interstate distance is not linear and proposes to consider it squared. In this application we carry out this exercise again but do not impose any functional form to the log of interstate distance and estimate the model semiparametrically. We then check if the square approximation is appropriate. More precisely, we run the following command lines:

. The results of the parametric part (see Stata output above) show that the distance from the incinerator does not seem to be significant (the t-stat associated with the coefficient is smaller than the critical value of 1.96).

As far as the effect of the log of the interstate distance is concerned, Figure (2) shows that it is clearly nonlinear. Indeed, when the interstate distance increases, the effect of house prices first increases and then decreases. When we check if the quadratic approximation proposed by [START_REF] Wooldridge | Econometric analysis of cross section and panel data[END_REF] is appropriate, it turns out that this assumption is clearly rejected by [START_REF] Hardle | Comparing nonparametric versus parametric regression fits[END_REF] test (see below). However, when we compare it with a polynomial adjustment of degree 3, the null is no longer rejected which means that instead of a semiparametric model, a pure parametric model with a polynomial fit of degree 3 of linst could be used.

The two Stata outputs below summarize results of the Hardle and Mammen (1993) test when the polynomial adjustment tested is of order 2 or 3 respectively. These outputs do not present the results concerning the parametric part since they are the same as in the output presented above. H0: Parametric and non-parametric fits are not different -------------------------------------------------------Test statistic T: .9705025 Critical value (95%): 1.959964 P-value: .3

Conclusion

In econometrics, semiparametric regression estimators have become standard tools for applied researchers. In this paper, we present [START_REF] Robinson | Root-N-consistent semiparametric regression[END_REF] double residual semiparametric regression estimator and [START_REF] Hardle | Comparing nonparametric versus parametric regression fits[END_REF] specification test. We then present the Stata codes we created to implement them in practice. Some simple simulations and an empirical application to illustrate the usefulness of the procedure are also

Figure 1 :

 1 Figure 1: Non-parametric fit of the four DGPs

Figure 2 :

 2 Figure 2: Nonlinear link between the price and interstate distance (in logs)

.

  use http://fmwww.bc.edu/ec-p/data/wooldridge/HPRICE3 . semipar lprice ldist larea lland rooms baths age, nonpar(linst) nograph test(2) > nsim(100) Simulation the distribution of the test statistic bootstrap replicates ............................................... 50 .................................................. 100 H0: Parametric and non-parametric fits are not different -------------------------------------------------------Test statistic T: 2.762337 Critical value (95%): 1.959964 P-value: 0 . use http://fmwww.bc.edu/ec-p/data/wooldridge/HPRICE3 . semipar lprice ldist larea lland rooms baths age, nonpar(linst) nograph test(............................................... 50 .................................................. 100

Table 1 :

 1 Comparison between semipar and plreg Bias x 2 Bias x 3 MSE x 2 MSE x 3

	plreg	-0.4695 -0.1039 0.2208	0.0111
	semipar -0.0435 -0.0183 0.0022	0.0007

Table 2 :

 2 Performance of the comparison Test T n

				Order tested
			0	1	2	3
		0 0.053 0.06 0.055 0.039
	True	1	1	0.064 0.055 0.021
	Order 2	1	1	0.06 0.062
		3	1	1	1	0.066
	Figures correspond to rejection rates of the
	test.

  use http://fmwww.bc.edu/ec-p/data/wooldridge/HPRICE3

	. semipar lprice ldist larea lland rooms baths age, nonpar(linst) xtitle(linst) ci
						Number of obs =	321
						R-squared	= 0.4437
						Adj R-squared = 0.4331
						Root MSE	= 0.2646
	lprice	Coef.	Std. Err.	t	P>|t|	[95% Conf. Interval]
	ldist	.1083941	.0640184	1.69	0.091	-.0175636	.2343519
	larea	.4887243	.0668208	7.31	0.000	.3572527	.6201959
	lland	.0866459	.036037	2.40	0.017	.0157423	.1575495
	rooms	.0436451	.0221781	1.97	0.050	9.12e-06		.087281
	baths	.0806555	.0335251	2.41	0.017	.014694		.146617
	age	-.003481	.0005436	-6.40	0.000	-.0045506	-.0024114

shown.

Acknowledgments

We would like to thank all our colleagues at CRED, ECARES, and CERPE as well as Emmanuel Flachaire, Marjorie Gassner, Whitney Newey, Christopher Parmeter and Sébastien Van Bellegem. Vincenzo Verardi is an associated researcher of the National Science Foundation of Belgium and gratefully acknowledges their financial support. Nicolas Debarsy is a research fellow of the National Science Foundation of Belgium and gratefully acknowledges their financial support.


Luc Gillibert

Dominique Jeulin

Stochastic Multiscale Segmentation Constrained by Image Content

Keywords: Granular media, granulometry, multiscale 3D segmentation, stochastic watershed, constrained segmentation

We introduce a noise-tolerant segmentation algorithm efficient on 3D multiscale granular materials. The approach uses a graphbased version of the stochastic watershed and relies on the morphological granulometry of the image to achieve a content-driven unsupervised segmentation. We present results on both a virtual material and a real X-ray microtomographic image of solid propellant.

Introduction

The stochastic watershed segmentation was first introduced by Angulo and Jeulin in [START_REF] Angulo | Stochastic watershed segmentation[END_REF]. The approach is based on using a large number of realizations of random markers to build a probability density function (pdf) of contours, starting from a standard watershed algorithm producing oversegmentations.

The stochastic watershed was proved to be efficient for unsupervised segmentation [START_REF] Noyel | Random germs and stochastic watershed for unsupervised multispectral image segmentation[END_REF] [START_REF] Faessel | Segmentation of 3D microtomographic images of granular materials with the stochastic watershed[END_REF]. The two parameters used for its construction are k, the number of random markers used in each realization, and R, the number of realizations. From the law of large numbers, the pdf converges when increasing R. The parameter k needs to be proportional to the number of desired regions in the segmented image. Therefore, in the case of granular materials, k needs to be proportional to the number of grains contained in the image, that can be automatically estimated. In [START_REF] Faessel | Segmentation of 3D microtomographic images of granular materials with the stochastic watershed[END_REF], the authors use the covariance for this estimation, deduced from the average radius of the grains, and a Boolean model assumption [START_REF] Jeulin | Modèles morphologiques de structures aléatoires et de changement d'échelle[END_REF].

For multiscale images, the covariance is not so efficient. In this paper, we introduce a new approach for stochastic segmentation which relies on the full granulometry of the image. Using morphological openings, this granulometry can be automatically computed from the image and is used as a constraint during iterations of segmentation steps.

Stochastic Segmentation

Stochastic Watershed

The first method introduced for computing the stochastic watershed is based on a large number of realizations of random markers to estimate a pdf of con-tours, or of surface boundaries in 3D. The random markers are generated with an uniform distribution corresponding to a constant intensity. In the case of granular materials, a constant background marker, extracted by an automatic thresholding, is added to each set of random markers. For each set of markers, a constrained watershed is computed. Then, the Parzen window is used to estimate the pdf of contours.

For a good estimation of the stochastic watershed, 100 to 200 realizations are required [START_REF] Angulo | Stochastic watershed segmentation[END_REF]. However, using λ-flat zones, a stochastic watershed segmentation can be achieved with only 50 realizations [START_REF] Faessel | Segmentation of 3D microtomographic images of granular materials with the stochastic watershed[END_REF]. This number is low, but the computation of 50 watersheds is very time-consuming, especially on large 3D data sets.

From the pdf, it is possible to obtain the segmentation. The first approach uses this pdf as a gradient for a new watershed [START_REF] Angulo | Stochastic watershed segmentation[END_REF]. A more efficient approach uses λ-flat zones to overcome the fact that the estimated pdf is not constant over each branch of contour [START_REF] Faessel | Segmentation of 3D microtomographic images of granular materials with the stochastic watershed[END_REF]. Illustration of the pdf of contours is on Fig. 1 (A). The resulting segmentation is illustrated on Fig. 2 (A).

Graph-Based Stochastic Watershed

Computing a large number of watersheds from simulations provides good results but is a slow process, mainly in 3D. A more efficient solution for computing stochastic watersheds is to use a graph-based approach. Probability of boundaries is directly computed with a good approximation without the use of any realisation [START_REF] Jeulin | Remarques sur la segmentation probabiliste[END_REF][START_REF] Meyer | Tailored hierarchies for segmentation[END_REF]. A first watershed is computed from the local minima of the gradient, as in standard segmentation, but restricted to the complementary set of the background extracted by an automatic thresholding. For this purpose, the background is used as a marker. In the present case, a very strong oversegmentation is obtained as a result of the presence of noise. From this watershed, an adjacency graph is constructed. Vertices of the graph are associated to each basin of the watershed, connecting edges between adjacent regions. A vertex is associated to the background too. Values are given to the vertices corresponding to the volumes of the regions (or to the areas, in 2D). Each edge is labelled with the minimum of the gradient function on the boundary between the corresponding regions. From this valued graph, a minimum spanning tree is extracted [START_REF] Meyer | Tailored hierarchies for segmentation[END_REF]. Then the regions in the minimum spanning tree are merged, starting with the edge of lowest value. After each merging, the probability p of the boundary between the corresponding regions is estimated using the following equation [START_REF] Jeulin | Remarques sur la segmentation probabiliste[END_REF]:

p = 1 -(1 -V 1 /V ) k -(1 -V 2 /V ) k + (1 -(V 1 + V 2 )/V )) k (1)
This is the probability of obtaining at least one random marker in each of the two regions, knowing the volumes V 1 and V 2 of the two regions, the total volume V of the image and the number of markers k. It is seen from equation 1 that the probability p increases with the volume of the grain, that may cause a bias towards largest grains if there is a wide distribution of sizes in the image.

After the merging of all the nodes in the original minimum spanning tree, the probability of all the edges of the tree is known. The result is projected from the tree on the graph and from the graph on the image.

This approach provides uniform probability on each part of boundary between two regions, as illustrated on Fig. 1 (B). Therefore, the λ-flat zones are useless and a simple threshold can be used for the segmentation Fig. 2 (B).

Multiscale Image Segmentation

The multiscale image segmentation process is based on a simple idea: estimate the full granulometry of the image, using morphological openings, then use multiple stochastic watersheds with different numbers of markers for each size, and finally combine them to get a segmentation that is pertinent for each size of grains. For achieving this goal, a hierarchy on the boundaries of the stochastic watershed is required. Many hierarchical segmentation algorithms have been studied, as the waterfalls [START_REF] Beucher | Watershed, hierarchical segmentation and waterfall algorithm[END_REF] or the P algorithm [START_REF] Beucher | P algorithm, a dramatic enhancement of the waterfall transformation[END_REF]. Here we introduce an approach similar to the one used for the computation of the waterfall: the merging of the watersheds basins using a minimum spanning tree [START_REF] Marcotegui | Fast implementation of waterfall based on graphs[END_REF].

Hierarchy on Boundaries

Each boundary is an edge in the adjacency graph. For a non-multiscale graphbased stochastic watershed, a fixed threshold is efficient. So, a first approach for constructing a hierarchy on the edges of the adjacency graph is to use a threshold. All the edges with a probability less than a given value t are removed from the graph. With t ∈ [0, 1], a complete hierarchy is obtained [START_REF] Jeulin | Remarques sur la segmentation probabiliste[END_REF].

This approach is fast and easy to compute, but, for a large value of t the remaining edges are not exploitable as boundaries for a segmentation. The removal of edges implies the removal of boundaries and therefore the merging of regions in the images. This process changes the probability of the boundaries of the remaining regions if we iterate the segmentation. This phenomenon is ignored with a fixed threshold. As seen on Fig. 2 (B), some over-segmented grains will disappear by a simple trheshold.

Working with a graph, it is possible to update the probability of the boundaries of the remaining regions at each edge deletion with the following iterative algorithm [START_REF] Jeulin | Remarques sur la segmentation probabiliste[END_REF]:

-While the MST has at least 2 vertices.

• In the MST, the edge e with the lowest probability is chosen.

• The edge e is deleted and the adjacent vertices are merged.

• The probability of the edges incident to the merged vertex is updated, using Eq. 1 and the new volume of the vertex.

This leads to better hierarchy on boundaries. As seen on Fig. 3, with the threshold approach some over-segmented grains disappear easily, while this problem is solved with the merging algorithm.

Granulometry-Driven Multiscale Approach

In order to get a pertinent segmentation starting from a highly over-segmented image, we can introduce constraints in the merging process. In the present case, the first step of the approach is to estimate the full granulometry of the image, using morphological openings. For fast computations, we use a rhombicuboctahedron as structuring element and we work on the binarized (thresholded) image. Gaussian noise, border effects and the shape of the structuring element leads to a few errors on the granulometry, as seen on Fig. 4, but they will not induce errors in the segmentation process. From the granulometry, a small number of classes are chosen. A good approach for this choice is maximising the interclass variance. For our sample material R 1 , we uses 3 classes: radius [13,18], radius [START_REF] Faessel | Segmentation of 3D microtomographic images of granular materials with the stochastic watershed[END_REF]13] and radius [START_REF] Angulo | Stochastic watershed segmentation[END_REF][START_REF] Marcotegui | Fast implementation of waterfall based on graphs[END_REF]. The total volume of the grains in each class will be denoted v(x). The number of grains in each class is deduced from v(x) and is denoted n(x). It is used to generate the corresponding number of markers in the calculation of the probability.

(A) (B) (C) (D) (E) (F) (G) (H)
Then, the standard watershed is computed from the local minima of the gradient. From this watershed, the adjacency graph is constructed and a minimum spanning tree is extracted.

The first class is chosen. The stochastic watershed is computed with a number of markers equals to n(1). In our 3D example (material R 1 ), 417 markers are needed. Based on this stochastic watershed, a first hierarchy on boundaries is computed with the merging algorithm. For each step i of the hierarchy, the granulometry of the corresponding segmentation is computed (v i (1)).

In the full hierarchy, there is a size step i * which minimises the difference |v i (1) -v(1)|. In our example, i * = 40 (Fig. 11 and Fig. 5). This step is used for the segmentation of the grains in the first class. All the segmented grains are removed from the image and added to the background mask. The minimum spanning tree is updated and the next class is chosen.

The same process is applied to all the classes. When no more class are left, all the segmentations are combined together. This provides the final result illustrated on Fig. 2 (C) and Fig. 12. The main lines of the multiscale image segmentation process are summarized on Fig. 6.

Validation

Random Model

For the validation, we use a simple random model. We generate 7560 random discrete spheres. The distribution of the radius of the spheres is fixed and is given on Fig. 4. The centres of the spheres are uniformly distributed in a cubic volume of size 300 × 300 × 300 with a rejection sampling algorithm (acceptance-rejection method). Two spheres cannot overlap more than a given threshold (fixed to 2 voxels). Spheres have non-uniform grey values, and show a visible boundary generating information in the gradient image. The volume fraction of the grains on the sample is 0.432. A strong Gaussian noise is finally added. A slice of the simulated material and of the corresponding binary mask is illustrated on Fig. 7.

Microtomographic Images

Our second sample is an X-ray microtomographic image provided by the CEA Gramat. The image was generated at with a Skyscan 1172 high-resolution micro-CT. We are operating on a 954 × 243 × 243 subimage. It is a solid propellant sample with fragmented grains due to a mechanical impact. A voxel is 1.80µm. The original diameter of the grains is 400 µm, but there are many small fragments. A slice is shown on Fig. 8. Results of the segmentation are illustrated on Fig. 10 and Fig. 9, where it is clear that the multiscale approach overcomes the standard stochastic watershed algorithm.

Conclusion

On both the simulated image and the X-ray microtomographic image, our multiscale approach provides a good segmentation (it is essential to use a "volume" weighted granulometry, all our attempts with a "number" granulometry being unsuccessful, due to its inherent sensitivity to noisy data). On the simulation, there is no over-segmentation for large grains, and the small grains are present as illustrated on Fig. 2 (C). On the microtomographic image, some of the large grains are over-segmented, but most of the small fragments are present and the results are better than the original stochastic watershed approach as illustrated on Fig. 10. The granulometry of the segmented image fits to the morphological granulometry of the input images, since it was used as a constraint in the merging process (Fig. 12 and Fig. 9). This new segmentation technique, combining an iterative probability-based merging of boundaries and a size distribution constraint, is very robust with respect to the noise contained in the image, without the necessity to apply a filter that would destroy the smallest grains. In addition, the segmentation is fully non-parametric, since at every step the required parameters are automatically deduced from the image.

(A) (B)
Due to the graph approach, the process is much faster than the original stochastic watershed, as summarized on Tab. 1.

Finally, it is possible to generalize the process to multi-criteria segmentation, using any other measurement tool as a constraint during the multiscale segmentation.

Fig. 1 .

 1 Fig. 1. (A) Pdf of contours on a simulated 3D material R1, estimated with 50 realizations of a Poisson point process (slice). (B) Pdf of contours computed with a graph-based approach (slice).

Fig. 2 .

 2 Fig. 2. (A) Stochastic watershed segmentation (slice). The λ-flat zones are used. A few grains are over-segmented. Small grains are missing (B) Graph-based stochastic watershed segmentation with a threshold of 0.5 on the probability (slice). A few oversegmented grains, as well as small grains are missing. (C) Final granulometry-driven multiscale segmentation result for the random material R1 (slice of a 3D image).

Fig. 3 .

 3 Fig. 3. (A) Probability density function of contours on a 3D material R1, computed with a graph-based approach (detail of a slice). (B), (C), (D): hierarchy on boundaries using increasing thresholds. (A), (E), (F), (G), (H): hierarchy on boundaries using the merging algorithm, step 0 is (A).

Fig. 4 .Fig. 5 .

 45 Fig.4. Input granulometry of the random material R1 and granulometry measured on the thresholded image from morphological openings. On x axis, the radius, and on y axis the volume fraction of grains with this radius (in percent).

Fig. 6 .

 6 Fig.6. The main steps of the multiscale image segmentation process.

Fig. 7 .

 7 Fig. 7. (A) Random material R1 with a strong Gaussian noise and a wide grain size distribution (slice). (B) The binary mask for material R1 (slice), the threshold is calculated via the maximization of the interclass variance.

Fig. 8 .

 8 Fig.8. Image E1: a 3D X-ray microtomographic image of a fragmented granular material and its binarization by maximisation of the interclass variance (slice).

Fig. 9 .Fig. 10 .

 910 Fig.9. Morphological granulometry of the X-ray microtomographic image of a fragmented granular material and granulometry of the final segmented image. On x axis, the diameter (µm), and on y axis the volume fraction of grains with this diameter.

Fig. 11 .

 11 Fig. 11. Morphological granulometry of the random material R1 and granulometry obtained at different steps of the hierarchy. The best result is step 40.

Fig. 12 .

 12 Fig. 12. Morphological granulometry of the random material R1 and granulometry of the final segmented image.

Table 1 .

 1 Computational cost of the stochastic watershed, the graph-based stochastic watershed and our multiscale image segmentation approach. Times are give on a 3.00 GHz Pentium 4.

	Algorithm	Image	Time
	Stochastic Watershed	Random material R1 1h 36m 36s
		Propellant (microtomographic) 3h 9m 21s
	Graph-Based Stochastic Watershed	Random material R1	5m 11s
		Propellant (microtomographic)	14m 9s
	Multiscale Image Segmentation	Random material R1	9m 22s
		Propellant (microtomographic)	25m 7s

Acknowledgments: This work was supported by a grant fom DGA (contract 2009 34 0006). The authors are grateful to Alain Fanget (CEA Gramat) for his advice during this study.


Marko Kryvobokov
email: mark.kryvobokov@let.ish-lyon.cnrs.fr

Is it worth identifying service employment (sub)centres when modelling apartment prices?

Keywords: service employment centres, centrality, accessibility, apartment price, hedonic modelling service employment centres, centrality, accessibility, apartment price, hedonic modelling

The use of the attributes of the central business district and several subcentres instead of the characteristics of all the land parcels or zones can be seen as a higher level of analysis in real estate valuation. However, old technological limitations on considering smaller territorial units are being successfully overcome. The question is whether or not we still need generalisation, i.e. to identify urban centres when modelling real estate prices, or whether it is preferable to operate at a lower spatial level. The application of the traditional approach of identifying centres is compared with an "objective" centrality index and a "subjective" accessibility index calculated for each zone. The purpose is to find out, which of the three concepts best fits a regression model of apartment prices and provides the best prediction. Both global and geographically weighted ordinary least squares regressions are used as well as spatial lag and spatial error models. We conclude that if a model is spatially weighted or the spatial effects are controlled, it is not that important which of the concepts is applied. Nevertheless, in most cases the highest predictive capacity is obtained with duocentric models.

Introduction

The tradition to consider urban centres in different aspects of urban study is as old as urban modelling itself. Begun by von Thünen, it was strongly theoretically developed by Alonso and many others. With growth of secondary urban centres, Wingo, Wendt, Harris and Ullman and others shifted the focus from monocentric to polycentric models (e.g. [START_REF] Merlin | Méthodes quantitatives et espace urbain[END_REF][START_REF] Harvey | Urban Land Economics[END_REF]. Secondary centres (subcentres) are supposed to be sufficiently large to significantly influence the urban structure, including such crucial components as travel patterns and real estate values.

The use of the attributes of the central business district (CBD) and several subcentres instead of the characteristics of all land parcels/neighbourhoods/districts/zones can be seen as an attempt to achieve a higher level of analysis. In the real estate terms [START_REF] Grissom | Real estate valuation: guide to investment strategies[END_REF] it corresponds to the transition from the second level of location (the relationship of a site to its surroundings) to the third level (the overall urban structure and the interrelationships of a community's land use pattern). In terms of urban geography (e.g. [START_REF] Sanders | Agent Models in Urban Geography[END_REF] it corresponds to a transition from a meso-geographical level to a macro-geographical level of a city. Thus, the identification of urban centres should be not an excessive simplification, but a reasonable generalisation of the description of reality.

With the rapid development in GIS and transportation analysis software, the old technological limitations on considering smaller territorial units are being successfully overcome, and research efforts can be directed to more detailed analyses, where, in principle, each land parcel in a city can be taken into account. Thus, we ask whether we still need generalisation, i.e. identifying urban centres, or whether we can operate with integral accessibility measures at meso-or even micro-geographical levels when modelling real estate prices. In the current study we address this question at the geographical level of Traffic Analysis Zones (TAZs). Our study was motivated by the existence of several important centres in Lyon that are difficult to identify formally and include in a hedonic price model. These difficulties stimulated us to apply alternative methods that avoid explicitly considering centres. The purpose of the paper is to find out, which of the three concepts: travel time to urban centres, an "objective" centrality index or a "subjective" accessibility index, best fits the hedonic model and provides the best prediction. In the hedonic modelling, we apply four approaches: global ordinary least squares (OLS), geographically weighted OLS, spatial lag, and spatial error.

Similarly to [START_REF] Sivitanidou | Do Office-Commercial Firms Value Access to Service Employment Centers? A Hedonic Value Analysis within Polycentric Los Angeles[END_REF], instead of using employment centres in general, we focus on service employment centres. We wish to avoid a site without commercial services, but with a large industrial enterprise, being identified as a subcentre. Thus, we follow the suggestion of [START_REF] Mcdonald | Urban Transportation and Land Use[END_REF] to consider employment density by industry sector.

The paper is organised as follows. Section 2 analyses the relevant studies from the urban economics, transport planning, and real estate valuation literature. The service employment centres in the Lyon Urban Area are identified in the Section 3 by applying residual analysis. Section 4, uses the other two approaches, namely the centrality index based on travel time as a result of transportation modelling, and the accessibility index based on a travel survey. Section 5 creates hedonic regression models of residential real estate prices, exploiting the identified service employment centres and the indices.

Conclusions are then drawn.

Literature review

Theoretical models of formation of non-monocentric patterns developed in the urban economics literature during the last three decades include the studies by [START_REF] Fujita | Multiple Equilibria and Structural Transition of Nonmonocentric Urban Configurations[END_REF], [START_REF] Fujita | A Monopolistic Competition Model of Spatial Agglomeration: Differentiated Products Approach[END_REF], [START_REF] Anas | General Equilibrium Models of Polycentric Urban Land Use with Endogenous Congestion and Job Agglomeration[END_REF] and [START_REF] Fujita | On the Endogenous Formation of Secondary Employment Centres in a City[END_REF]. In this paper the focus is on identifying existing subcentres and their influence on real estate prices, rather than on the formation of new subcentres.

Subjectivity in the identification of urban centres has been recognised and criticised [START_REF] Mcdonald | The Identification of Urban Employment Subcenters[END_REF][START_REF] Mcmillen | Evolving subcenters: employment and population densities in Chicago, 1970-2020[END_REF]. Several formal identification procedures have been developed. Thus, McDonald (1987) proposes identifying employment subcentres as secondary peaks in the gross employment density (employment divided by total land area) and the employment-population ratio. A peak means that all the adjacent zones outside a subcentre have a smaller density or ratio. [START_REF] Giuliano | Subcenters in the Los Angeles region[END_REF] and [START_REF] Small | Population and Employment Densities: Structure and Change[END_REF] apply a similar definition of a centre, which is a continuous set of zones, selected with cut-offs for density and total employment. In a more formal approach, employment density was found to be a function of the distance from the CBD as well as from subcentres. Applying a monocentric analysis of employment density, [START_REF] Mcdonald | Suburban Employment Centres: The Case of Chicago[END_REF] define subcentres as locations with significantly positive residuals. However, McMillen and McDonald (1994) noted that such an approach may overlook fairly large subcentres populated by firms with large internal scale economies because the employment effects may be highly localised. [START_REF] Mcmillen | Nonparametric Employment Subcenter Identification[END_REF] identified the potential subcentres as sites, which have statistically significant residuals of locally weighted regression of employment density on the distance from the CBD (the first stage). He then checked whether they provide significant explanatory power in a semiparametric employment density regression estimation (the second stage). [START_REF] Mcmillen | Nonparametric Employment Subcenter Identification[END_REF] applied the proposed procedure for six American metropolitan areas. [START_REF] Mcmillen | The number of subcenters in large urban areas[END_REF] applied this procedure to 62 large metropolitan areas in the USA. [START_REF] Craig | Using Quantile Smoothing Splines to Identify Employment Subcenters in a Multicentric Urban Area[END_REF] used a nonparametric employment density function, namely quantile smoothing splines.

The CBD is usually the primary focus in hedonic price models. Although the a priori CBD identification can be seen as a weak point, there are a few papers in the real estate domain where the CBD is not simply taken as the area usually referred to as the CBD. Söderberg and Janssen (1999) re-estimated their regression model changing the precise location of the CBD of Stockholm by a step of 50 metres; as a result, the best model of apartment prices has been obtained for a location one kilometre east of the place commonly viewed as the city centre. [START_REF] Sivitanidou | Do Office-Commercial Firms Value Access to Service Employment Centers? A Hedonic Value Analysis within Polycentric Los Angeles[END_REF]) used the McDonald (1987) definition of urban centres. At the same time, a formal approach to secondary centres in the real estate literature is rare. Among the relatively few examples are: [START_REF] Mcdonald | Employment subcenters and land values in a polycentric urban area: The case of Chicago[END_REF], [START_REF] Mcmillen | One Hundred Fifty Years of Land Values in Chicago: A Nonparametric Approach[END_REF] and [START_REF] Sivitanidou | Do Office-Commercial Firms Value Access to Service Employment Centers? A Hedonic Value Analysis within Polycentric Los Angeles[END_REF].

As [START_REF] Mcmillen | Evolving subcenters: employment and population densities in Chicago, 1970-2020[END_REF] note, it is important to operate with an optimal number of subcentres. On the one hand, listing too many centre sites produces inefficient estimates and can influence other estimated coefficients when distance to a subcentre is highly correlated with other explanatory variables. On the other hand, incorrectly omitting subcentres causes other estimates to be biased. [START_REF] Ross | Inconsistency in Welfare Inferences from Distance Variables in Hedonic Regressions[END_REF] have highlighted the common inability to fit more than two distance variables arguing that two points in space triangulate the optimal position by fundamental geometry.

In the discussion above, the identification of urban centres and their accessibility measures as distance or travel time were addressed. However, the concept of accessibility is not limited to urban centres. Despite being the focus of research in transport planning for a long time (see [START_REF] Hansen | How Accessibility Shapes Land Use[END_REF][START_REF] Morris | Accessibility indicators for transport planning[END_REF], accessibility remains a rather illusive concept (Miller, 2008). As [START_REF] Morris | Accessibility indicators for transport planning[END_REF] note, there is a critical distinction between the derivation of "objective" indicators of accessibility, and perceived measures. The former refers to location of opportunities and potential access to them [START_REF] Morris | Accessibility indicators for transport planning[END_REF][START_REF] Krizek | Perspectives on Accessibility and Travel[END_REF] and is related to the concept of centrality [START_REF] Samaniego | Cities as organisms: Allometric scaling of urban road networks[END_REF]. The latter concerns realisation of this potential in terms of actual travel. In this respect Des Rosiers and [START_REF] Des Rosiers | Mass Appraisal, Hedonic Price Modelling and Urban Externalities: Understanding Property Value Shaping Process[END_REF] define accessibility as the ease with which persons, living at a given location, can move to reach activities and services which they consider to be the most important. This, mostly behavioural and subjective, concept of accessibility is quite distinct from centrality which relies on structural features and relates to the proximity to urban amenities. These concepts of accessibility and centrality are exploited in the present paper. [START_REF] Anas | General Equilibrium Models of Polycentric Urban Land Use with Endogenous Congestion and Job Agglomeration[END_REF] have analysed the urban economics studies which did not pre-specify any centres but used rent gradients and land use density peaks around the most accessible place(s) in urban space. [START_REF] Thériault | Modelling accessibility to urban services using fuzzy logic. A comparative analysis of two methods[END_REF] and Des Rosiers and [START_REF] Des Rosiers | Mass Appraisal, Hedonic Price Modelling and Urban Externalities: Understanding Property Value Shaping Process[END_REF] have provided examples of hedonic modelling of real estate prices without explicitly considering the CBD and secondary centres, but with integral measures of centrality and accessibility. They found that the perceptual index of accessibility, based on interview and fuzzy logic criteria, far outweighs the centrality index in the hedonic model of housing prices in Quebec City.

Identification of service employment centres

The Lyon Urban Area1 (Figure 1) is the second largest by population in France. The The historic centre of the city of Lyon, which was founded in the first century BC, is located three kilometres to the North of the confluence of the Rhône and the Saône. Further development, in the medieval period and later, mainly occurred close to the historic centre, in the Peninsula (the area between the two rivers close to their confluence). Nowadays the Peninsula, which contains the city hall, other administrative buildings and a large shopping district, is traditionally considered to be the city centre of Lyon. Its main transportation junction, Bellecour-Sala, is located in the middle of the Peninsula. Part-Dieu, the largest shopping centre in the Lyon Urban Area, was built in the 1970s and is located to the east of the Peninsula.

There are 812 TAZs in the Lyon Urban Area (see Figure 1). A zone corresponds to a French statistical unit IRIS (les îlots regroupés pour l'information statistique). The average zonal population slightly exceeds two thousand inhabitants. In the current paper, a zone is used as the spatial unit of data collection and analysis.

A peculiarity of French statistics is that the data about commercial employment itself are not available. Instead, INSEE supplies the number of employees in the tertiary economic sector in 2006. This sector contains "services" in a very broad sense including commerce, education, medicine, transport and other spheres. The highest number of service employees (above 14,000) is in Part-Dieu, the second highest number (above 12,500) is in Bellecour-Sala.

We use travel times between the centroids of the zones. The origin-destination (O-D) matrix of travel time for 2007 for this study was obtained from the MOSART4 transportation model for the Lyon Urban Area. Though public transport is well developed in the area, the travel time by car provides better regression results. We use the travel times estimated for the travel by car in the morning peak. When the same zone is both the origin and the destination, the travel time is equal to the minimum among all the other cases, thus null values are avoided.

As [START_REF] Mcdonald | Suburban Employment Centres: The Case of Chicago[END_REF], we run a simple regression model of service employment density on travel time to Bellecour-Sala, which is considered to be the CBD5 , in order to find positive significant residuals. [START_REF] Mcmillen | Suburban Subcenters and Employment Density in Metropolitan Chicago[END_REF] argue that employment density functions are biased if only non-zero densities are included. In our case, however, the zones are sufficiently large and only five of them, located in different districts, have zero densities. Of the 812 zones, fifteen have positive standardised residuals higher than 3.3, i.e. very significant (see Table 1). In further analysis we will use all these pre-identified service employment centres in order to avoid the rather subjective step of selecting a cut-off point.

The centres are described in Table 1 andshown in Figure 2. Twelve of them are situated in the 1 st , 2 nd , 3 rd , 6 th and 7 th arrondissements of Lyon, and the other three (Stalingrad, Charles Hernu and Gratte Ciel est) are in Villeurbanne. The service employment density and service employment to population ratio are the best measures with which to identify urban centres [START_REF] Mcdonald | The Identification of Urban Employment Subcenters[END_REF][START_REF] Mcdonald | Urban Transportation and Land Use[END_REF][START_REF] Sivitanidou | Do Office-Commercial Firms Value Access to Service Employment Centers? A Hedonic Value Analysis within Polycentric Los Angeles[END_REF]. In the current study, the former measure provides better results and so is used.

Centrality and accessibility measures

A service centrality index for zone i is calculated with a simple gravity-like model:

∑ = = N j ij j i tt A CI 1 ,
where j A -the attraction of zone j; ij tt -the travel time from zone i to zone j6 ;

N -the number of zones.

The attraction of a zone is its service employment density. The normalised service centrality indices calculated for the fifteen pre-identified centres are shown in Table 1, where i CI for each zone is divided by the maximum value and multiplied by 100, as in [START_REF] Thériault | Modelling accessibility to urban services using fuzzy logic. A comparative analysis of two methods[END_REF]. Figure 3 shows the clusters7 of the normalised centrality indices grouped into five classes. As in [START_REF] Thériault | Modelling accessibility to urban services using fuzzy logic. A comparative analysis of two methods[END_REF], we estimate a suitability index applying fuzzy membership and using the 50 th percentile and 90 th percentile satisfaction thresholds from the survey for travel times from the O-D matrix. Any travel time less than the observed 50 th percentile is seen to be acceptable. A travel time larger than the 90 th percentile is unsatisfactory. Intermediate cases yield satisfaction levels obtained by linear interpolation. A suitability index ij S for travelling from zone i to zone j is calculated by the following formula from the aforementioned source:

1 = ij S ∀ 50 C tt ij ≤ , ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - - = 50 90 50 1 C C C tt S ij ij ∀ 90 50 C tt C ij < < , 0 = ij S ∀ 90 C tt ij ≥ ,
where ij tt -the travel time from zone i to zone j; Table 2 includes the values of percentiles and the number of cases in the O-D matrix, where the suitability index has a value of unity, between zero and unity, or zero.

∑ = = N j j ij i A S AI 1 ,
where ij S -the suitability index for travelling from zone i to zone j; j A -the attraction of zone j;

N -the number of zones.

The attraction of a zone is its service employment density. In contrast to Thériault et al. (2005), we do not multiply the suitability index by the population of the zone. This is because in our study we are not analysing how many people can reach a particular zone, but rather we are analysing how attractive a zone is, taking into account the service employment of those zones, which can be reached from this particular zone.

The normalised service accessibility indices calculated for the fifteen preidentified centres are shown in Table 1, where form belts around the central part of the Lyon Urban Area with peak values in the city core (see Figure 3). The spatial configuration of the accessibility indices is more complex; however, their structure also resembles belts (see Figure 4). For remote locations, where distances from the identified centres are much longer than distances between the centres themselves, the centres are "merged" into a city core, like a whole city becomes a point on a smaller-scale map. Thus, to understand the individual influences of the identified centres, it is better to focus on real estate prices in the central part of the Lyon Urban Area.

Hedonic model of apartment prices

Data and model specification

In a hedonic price model, the dependent variable is price and the independent variables are real estate attributes and location attributes. The estimated parameters in the OLS can be interpreted as the willingness to pay for different attributes [START_REF] Rosen | Hedonic prices and implicit markets: product differentiation in pure competition[END_REF].

Hedonic regression analysis is widely used in investigations of real estate around the world. Examples of its application to apartment prices include: [START_REF] Asabere | Thoroughfares and Apartment Values[END_REF], [START_REF] So | Estimating the influence of transport on house prices: evidence from Hong Kong[END_REF], [START_REF] Watkins | Are new entrants to the residential property market informationally disadvantaged[END_REF], [START_REF] Brañas-Garza | The North-South divide and house price islands: the case of Córdoba (Spain)[END_REF], and [START_REF] Björklund | Correlation between Negotiated Rents and Neighbourhood Quality: A Case Study of Two Cities in Sweden[END_REF].

The data on sale prices and apartment attributes were provided by Perval, which collects information about real estate transactions in France. Data on approximately 10,000 apartment sales selected randomly from all sales in the central part of the Lyon Urban Area in the period 1997-2008 were obtained. After deleting observations with incomplete data, 4,362 apartments remained. The apartments are mainly located in Lyon and Villeurbanne and also in the surrounding urbanised area (Figure 5). The location variables in Table 3 include the percentage of middle-income households, the percentage of high-income households in zones, and dummies for proximity to water and location in one of four ad hoc districts. The percentages of middle-and high-income households in zones were obtained from the INSEE data. The middle-income group includes households in the middle 60% of the income range and the high-income group is composed of the 20% households with the highest income. A dummy for location within a 100 metre buffer created for rivers and lakes is a proxy for a water view, though we admit that in densely built areas water is not necessarily visible from each apartment. The four ad hoc districts, created as proxies for submarkets, are quite large, but relatively homogenous territories, divided by water frontiers and the boundaries of the urbanised area. District 1 is the Peninsula and the urbanised area to the north of it, between the Rhône and the Saône. District 2 is an urbanised area on the left bank of the Rhône. District 3 is an urbanised area on the right bank of the Saône.

District 4 is the less urbanised territory, which occupies most of the area in Figure 1.

The following attributes are used as default values in the hedonic model: Year97, Bath1, Park0, FloorGr, Constr1981_1991, CondGood, ViewGood, Cellar0 and District1.

Global and GWR OLS models

The relative importance of variables in regressions is quite often discussed in the literature. As the choice of one or another concept of relative importance often affects conclusions [START_REF] Kruskal | Concepts of Relative Importance in Recent Scientific Literature[END_REF], it is important to select a meaningful measure.

We can use a contribution to adjusted R 2 and an unstandardised regression coefficient.

The latter is appropriate to compare variables, which have the same unit of measurement, e.g. travel times to different centres measured in minutes. The comparison of travel times with indices is more complicated. One could consider a standardised regression coefficient for this purpose. However, this beta coefficient has been much criticised in the statistics literature (e.g. [START_REF] Darlington | Regression and Linear Models[END_REF][START_REF] Bring | How to Standardize Regression Coefficients[END_REF]. As [START_REF] King | How Not to Lie With Statistics: Avoiding Common Mistakes in Quantitative Political Science[END_REF] noted, this measure is a mixture of the estimated effect and the standard deviation, which should be analysed separately.

Perhaps a better approach to the evaluation of different models is to define the in-sample estimates and then use the result for an ex-sample prediction. For this purpose, as e.g. in [START_REF] Bourassa | Do housing markets really matter[END_REF], we randomly select 80% of observations as the in-sample and the other 20% are the ex-sample. Using the latter, we calculate the percentages of predictions that deviate by less than 10% and 20% from the actual sale prices.

To investigate the existence of multicollinearity, we estimate the maximum of variance inflationary factors (VIF). The principle that a VIF in excess of 10 indicates multicollinearity is usually used in the literature (e.g. [START_REF] Seiler | The Impact of World Class Great Lakes Water Views on Residential Property Values[END_REF][START_REF] Thériault | Modelling accessibility to urban services using fuzzy logic. A comparative analysis of two methods[END_REF]. We measure global spatial autocorrelation in the error term with Moran's I [START_REF] Anselin | Local Indicators of Spatial Association -LISA[END_REF][START_REF] Dubin | Spatial Autocorrelation: A Primer[END_REF]; see also the empirical examples in De [START_REF] De Graaff | A general misspecification test for spatial regression model: dependence, heterogeneity, and nonlinearity[END_REF][START_REF] Munroe | Exploring the determinants of spatial pattern in residential land markets: amenities and disamenities in Charlotte, NC, USA[END_REF], which is a weighted correlation coefficient ranging between -1 and 1. In this paper, it is calculated with the row-standardised weight matrix of inverse squared distances.

The influence of the pre-identified centres was examined in the following way.

First, the OLS model with only apartment variables was estimated. This model with the adjusted R-square of 0.781 does not contain the dummies for construction periods, which are correlated with location. Its significant attributes are presented in the Appendix. We then added the location attributes and travel time to each of the fifteen pre-identified centres, one at a time. Thus we obtained fifteen global models with adjusted R-squared varying in the range 0.819-0.849. We do not report standard errors of regression and F-values, because, as shown in Söderberg and Janssen (1999), as the number of variables is the same for each run, the case which has the minimum standard error of regression, has maximum R-squared, maximum F-value and maximum adjusted R-squared. Sorting the adjusted R-squared from high to low, we added all fifteen variables to the equation and then excluded them one by one from the bottom in order to obtain a model with an acceptable VIF. Due to high multicollinearity, an acceptable VIF can only be achieved with a small number of centres in a model.

The two best global OLS models: with travel times to the CBD (Centre 3) and Les Belges (Centre 10), and with travel times to Les Belges (Centre 10) and Jussieu (Centre 6) are named Duocentric 1 and Duocentric 2, respectively, in Table 4 8 . It should be noted that in the former model, the travel time to Les Belges outweighs that to the CBD. However, in the latter model, the travel time to Jussieu outweighs that to Les

Belges. The estimates of the former model that are significant at the 5% level are presented in the Appendix. Other location variables behave as follows. Water does not significantly influence apartment prices. %MidIncome is more significant than %HighIncome, perhaps because of the relatively small percentage of high-income households. District1, used as the default ad hoc district, is the most attractive one, whereas District2 is the least attractive urbanised ad hoc district, and less urbanised District4 has the highest negative coefficient. The extracted results for the alternative models: the model with longitude X and latitude Y9 , the monocentric model, and those, which use a centrality index or an accessibility index instead of travel times, are also presented in Table 4. For all the models, a Jarque-Bera test and a Breusch-Pagan test indicate no rejections of the assumptions of normality and heteroskedasticity. The goodness-of-fit and prediction of the global models with two centres is better than those of the alternative models. The model with coordinates predicts prices better than other alternative specifications, but is worse than the duocentic models. Among the global models, Moran's I for residuals is the lowest for the duocentric models and the highest for the monocentric model.

To detect local peculiarities, we apply a geographically weighted regression (GWR) [START_REF] Brunsdon | Geographically Weighted Regression: A Method for Exploring Spatial Nonstationarity[END_REF] as an alternative to global regression modelling. The GWR model used in the study is OLS and the error term is Gaussian. The fixed kernel type is used; the kernel bandwith is determined by the Akaike Information Criterion. In the current application of the GWR, a separate equation is solved for each observation, and the overall results are reported in this paper. The extracted GWR results are shown in Table 4. For them, the goodness-of-fit of the models with two centres are also higher than those of alternative models, though the differences are less than in the global OLS. Moran's I for the indices is 1% better than that for the monocentric model, but 1% worse than that for the duocentric models.

In all cases, the GWR models outperform the global OLS models with respect to goodness-of-fit. Naturally, the spatial autocorrelation is also better controlled by the GWR. The predictive capacity measured using average estimates in each zone also demonstrates the superiority of the duocentric models, although the GWR monocentric model predicts not only better than the models with indices, but also better than Duocentric 1 within 20% of the sale price.

Spatial lag and spatial error models

Moran's I indicates that quite high spatial autocorrelation still exists10 . Spatial econometrics can be applied to take account of spatial dependence and spatial heterogeneity. Discussions of spatial models can be found in [START_REF] Anselin | Spatial Econometrics: Methods and Models[END_REF], [START_REF] Dubin | Spatial Autocorrelation: A Primer[END_REF] and [START_REF] Lesage | Introduction to Spatial Econometrics[END_REF]. We control for spatial effects with two models: a spatial lag model, where a dependent variable is not only a function of independent variables, but also of the dependent variables in nearby areas; and a spatial error model, where the error term is a function of the errors in neighbouring areas. In both cases, the weight matrix includes row-standardised binary weights, which are equal to unity for neighbouring observations and zero otherwise. The band is defined in such a way that there is at least one neighbour for all observations. Such a weighting scheme avoids the correlation of spatial structure with travel time variables (this problem, with respect to a distance variable, is discussed in [START_REF] Wilhelmsson | Spatial Models in Real Estate Economics[END_REF].

Parameter estimates in global OLS, GWR and a spatial error model have a straightforward interpretation as partial derivatives of the dependent variable with respect to the explanatory variable in question. In a model containing the spatial lag of the dependent variable, the interpretation of parameters becomes richer and more complicated [START_REF] Lesage | Introduction to Spatial Econometrics[END_REF]. [START_REF] Kim | Measure the benefits of air quality improvement: A spatial hedonic approach[END_REF] have shown that if a unit change was induced at every location and a weight matrix is row-standardised, the marginal effect can be calculated as a spatial lag estimate multiplied by (1 -ρ) -1 . LeSage and Pace ( 2009) have called this simple form a summary measure of total impacts.

The estimates of the duocentric spatial models with travel times to the CBD and Les Belges are presented in the Appendix. For most non-location variables the estimates of the spatial models are slightly lower than those of the OLS. The coefficients and significance of the travel time to Les Belges (Centre 10) are higher than those to the CBD (Centre 3) in the spatial error model and they (as well as their marginal effects) are practically equal in the spatial lag model.

The extractions of all the spatial models are presented in Table 4. A Breusch-Pagan test indicates that the assumption of heteroskedasticity is not rejected. Lower rho, i.e. lower spatial dependency, is observed for the duocentric models and for the models with indices. Lambda shows a similar tendency, though to a lesser degree; it is lowest for the duocentric models. For the spatial models, Table 4 contains pseudo R 2 . Both spatial models explain the variation in price better than the global OLS, but worse than the GWR OLS. Interestingly, in most cases the spatial autocorrelation of the residuals was either not decreased or decreased by only 1% in comparison with the global OLS.

Only the monocentric model demonstrates a larger decrease in Moran's I. This may be due to the weighting scheme which was applied.

Comparing the spatial lag results with other models, we can see that while the marginal effect of an accessibility index (0.005) is the same as the estimate in other models, the marginal effects of travel time variables are quite large, for example in Duocentric Model 1 it is equal to -0.203 for both of the centres that were examined.

The ex-sample predictions are calculated with the average in-sample dependent variable for the spatial lag model and with an average in-sample prediction error in each zone for the spatial error model. In most cases the best predictions are obtained with the duocentric models. Only in the spatial error model within the 10% interval of sale price the prediction is the best with an accessibility index, while it is second best with the monocentric model.

In general, when we control for spatial effects, the predictive capacity usually increases, especially with the spatial error model, though there are exceptions, the most visible being the model with a centrality index. The spatial models are no more helpful than the OLS when determining the leaderships or ranking among the three selected centres.

Conclusions

The first conclusion in this study concerns the importance of particular centres when modelling apartment prices in the Lyon Urban Area. The best results were obtained with travel times to three centres: Bellecour-Sala traditionally viewed as the CBD, Les

Belges, and Jussieu. However, due to high correlation between them, it was impossible to include all three in a single model (this can be done in a future study with a factor analysis), and two alternative duocentric models were created instead. We acknowledge, as in Söderberg and Janssen (1999), that we did not determine the alternative CBDs in the traditional CBD sense. The two subcentres found to be important are not the leaders according to their service employment attributes (see Table 1), but rather are highly desirable residential locations11 .

The second conclusion is related to the three concepts: identified urban centres, an "objective" centrality index, and a "subjective" accessibility index. In most cases both duocentric models have the highest predictive capacity. There are no large differences in the predictions for the GWR and the spatial lag models with one or another index, although an accessibility index provides much better prediction with the spatial error model. Thus, the answer to the question asked in the title is positive.

However, if a model is spatially weighted or the spatial effects are controlled, it is usually less important which of the concepts is applied. This is especially noticeable with the GWR methodology, where the differences between the results are rather marginal.

Appendix

  data on population from INSEE 2 refer to 2005 3 . The central part of the area with a population of 613 thousand people consists of the cities of Lyon and Villeurbanne. These cities, which have a common planning structure and transportation network are shown in Figure 1 (in white in the centre) divided by small zones and almost surrounded by other urbanised areas shown in dark grey. The population of these areas is 580 thousand inhabitants. The rest of the territory with 711 thousand inhabitants is less urbanised and is represented in light grey.

Figure 1 .

 1 Figure 1. The Lyon Urban Area and surrounding areas

Figure 2 .

 2 Figure 2. The pre-identified centres

Figure 3 .

 3 Figure 3. Clusters of centrality indices

50C

  -the 50 th percentile of the observed travel time; 90 C -the 90 th percentile of the observed travel time.

  i AI for each zone is divided by the maximum value and multiplied by 100. Figure 4 shows the clusters of the normalised accessibility indices grouped into five classes.

Figure 4 .

 4 Figure 4. Clusters of accessibility indices

Figure 5 .

 5 Figure 5. Location of apartments

Table 1 .

 1 Standardised residuals, description and indices of the pre-identified centres

						Tertiary
	Centre number	Centre name	Standardised residual	Tertiary employment	Population	100 sq. m employment density per	Centrality index	Accessibility index
		Stalingrad	14.47	4,113	1	1,103 98.74	89.49
		Louis	10.73	1,867	56	841 100.00	95.96
		Pradel
		Bellecour-	9.47	12,522	2,629	762 90.32	92.70
		Sala
		Victor	6.45	6,077	3,536	531 85.39	97.59
		Bach
		Molière	5.51	4,692	2,577	463 85.62	97.87
		Jussieu	5.35	4,726	1,832	450 86.12	98.69
		Saxe-	5.13	2,930	2,207	430 92.24	100.00
		Bossuet
		Mutualité-	4.76	6,232	3,001	415 84.94	98.05
		Liberté
		Charles	4.46	5,488	4,581	375 69.23	87.26
		Hernu
		Les	4.25	5,458	2,226	364 70.30	93.39
		Belges
		Villette	4.24	7,434	2,836	359 66.24	92.88
		Gare
		Gratte	3.82	3,525	4,020	320 55.84	78.36
		Ciel est
		Terreaux-	3.82	7,782	3,727	340 82.37	95.16
		Bat
		d'Argent
		Part-Dieu	3.48	14,205	2,869	311 71.49	98.46
		Marechal	3.35	3,978	3,086	300 88.30	98.62
		Lyautey

Table 2

 2

		. Description of suitability index
	50 C , minutes	90 C , minutes	S	ij	=	1	Number of cases where 1 0 < < ij S S	ij	=	0
	5.69	18.49	11,530	157,079	490,735
	A service accessibility index for zone i is calculated as follows:

Table 3

 3

also gives the travel times to each of the pre-identified centres and includes the definition and descriptive statistics of the normalised centrality index and accessibility index for the zones, where the analysed apartments are located. The logarithmic transformations of Area, %MidIncome, %HighIncome and travel times are used. The dependent variable is the logarithm of Price.

Table 3 .

 3 Definition of variables and descriptive statistics

	Variable	Description	Mean	Minimum	Maximum	Std. deviation
	Price	Transaction price, Euros	123,635.86	12,196	1,120,000	76,838
	Year97-	Dummies for year of	0.02-0.13	0	1	0.15-0.34
	Year08	transaction
	Area	Apartment area, square	69.05	8	301	28.14
		metres
	Bath1-Bath3	Dummies for number of	<0.01-0.93	0	1	0.05-0.26
		bathrooms
	ParkUn	Dummy for cases with no	0.26	0	1	0.44
		data about parking places
	Park0-Park3	Dummies for number of	<0.01-0.50	0	1	0.06-0.50
		parking places
	FloorGr	Dummy for ground floor	0.13	0	1	0.33
	Floor1	Dummy for storey 1	0.19	0	1	0.39
	Floor2_4	Dummy for storey 2 to 4	0.49	0	1	0.50
	Floor5_8	Dummy for storey 5 to 8	0.18	0	1	0.38
	Floor9+	Dummy for storey 9 or	0.02	0	1	0.14
		more
	Constr<1850-	Dummies for period of	0.03-0.34	0	1	0.17-0.48
	Constr1992<	construction
	CondGood	Dummy for good state	0.81	0	1	0.39
	CondMed	Dummy for state when	0.16	0	1	0.37
		some maintenance is
		needed
	CondBad	Dummy for state when	0.03	0	1	0.17
		renovation is needed
	ViewNo	Dummy for cases with no	0.60	0	1	0.49
		data about view
	ViewGood	Dummy for view	0.38	0	1	0.48
		increasing value
	ViewBad	Dummy for view	0.02	0	1	0.13
		decreasing value
	Cellar0-	Dummies for number of	0.02-0.33	0	1	0.13-0.47
	Cellar2	cellars
	Garden	Dummy for existence of	0.05	0	1	0.22
		garden
	Terrace	Dummy for existence of	0.09	0	1	0.29
		terrace
	%MidIncome	Percentage of middle-	57.96	42.70	66.20	3.30
		income households
	%HighIncome	Percentage of high-	12.55	4.34	28.77	2.91
		income households
	Water	Dummy for location	0.03	0	1	0.18
		within a 100 m buffer of
		water
	District1-	Dummies for location in	0.01-0.16	0	1	0.12-0.37
	District4	districts
	Travel time to	Travel time to centre 1 -	9.57-11.67	0.45	24.43-31.28	4.85-5.67
	Centre 1 -	travel time to centre 15
	Travel time to
	Centre 15
	Centrality	Centrality index	42.19	13.45	100.00	17.21
	Index
	Accessibility	Accessibility index	63.23	0.34	100.00	27.39
	Index

Table 4 .

 4 The extracted results of hedonic models Estimates

	Model	Travel time to	Centre 3	Travel time to	Centre 10	Travel time to	Centre 6	Centrality	index	Accessibility	index	Adjusted R 2 /	Pseudo R 2	Maximum VIF	Moran's I	Rho/lambda	Within 10% of	sale price	Within 20% of	sale price

. Duocentric 1: global OLS without location attributes, global OLS, and spatial models

		Global OLS, no location attributes		Global OLS		Spatial lag	Spatial error
	Variable	Coeffi-cient	t-value	Coeffi-cient	t-value	VIF	Coeffi-cient	Assym-ptotic t-value	Coeffi-cient	Assym-ptotic t-value
	Constant	6.982	106.34	4.364	11.34	-	0.522	0.66	4.584	11.92
	Travel time to Centre 3	-	-	-0.148	-12.63	2.129	-0.137	-11.71	-0.123	-8.56
	Travel time to Centre 10	-	-	-0.174	-19.26	1.678	-0.137	-12.06	-0.159	-13.15
	Year99	0.140	3.90	0.099	3.32	4.303	0.095	3.20	0.092	3.13
	Year00	0.197	5.43	0.158	5.24	4.154	0.156	5.23	0.155	5.25
	Year01	0.268	7.59	0.247	8.42	4.901	0.241	8.29	0.243	8.41
	Year02	0.319	9.04	0.307	10.48	5.036	0.299	10.33	0.302	10.48
	Year03	0.480	13.91	0.458	16.01	5.979	0.452	15.93	0.456	16.21
	Year04	0.656	18.66	0.638	21.88	5.117	0.630	21.80	0.637	22.19
	Year05	0.803	22.87	0.791	27.17	5.127	0.783	27.14	0.788	27.48
	Year06	0.932	26.11	0.908	30.59	4.496	0.901	30.66	0.906	31.03
	Year07	1.027	28.57	1.003	33.56	4.332	0.997	33.70	1.002	34.11
	Year08	0.968	24.21	0.969	29.18	2.632	0.963	29.32	0.966	29.60
	Area	0.902	66.23	0.922	80.84	1.442	0.927	81.91	0.927	82.24
	Bath2	0.139	6.65	0.077	4.40	1.247	0.075	4.36	0.075	4.37
	ParkUn	0.228	14.78	0.073	5.08	2.427	0.070	4.95	0.057	3.99
	Park1	0.283	19.93	0.150	11.54	2.674	0.145	11.28	0.138	10.69
	Park2	0.366	16.70	0.197	10.13	1.793	0.192	9.98	0.182	9.50
	Park3	0.154	2.03	-0.012	-0.19	1.090	0.191	3.639	-0.026	-0.42
	Floor1	0.101	5.10	0.068	4.12	2.620	0.068	4.19	0.065	4.00
	Floor2_4	0.149	8.31	0.094	6.26	3.576	0.094	6.35	0.092	6.21
	Floor5_8	0.160	7.97	0.115	6.88	2.652	0.113	6.80	0.112	6.82
	Floor9+	0.118	3.04	0.112	3.47	1.238	0.109	3.42	0.114	3.59
	CondMed	-0.165	-12.01	-0.108	-9.22	1.170	-0.109	-9.41	-0.106	-9.20
	CondBad	-0.244	-8.20	-0.213	-8.55	1.106	-0.216	-8.74	-0.216	-8.81
	ViewNo	-0.030	-2.56	-0.037	-4.23	1.130	-0.036	-4.18	-0.037	-4.30
	ViewBad	-0.127	-3.34	-0.111	-3.49	1.089	-0.109	-3.47	-0.108	-3.47
	Cellar1	-0.042	-3.76	0.027	2.69	1.391	0.023	2.37	0.018	1.81
	Cellar2	0.056	1.41	0.070	2.10	1.147	0.066	2.02	0.062	1.90
	Garden	0.099	3.71	0.057	2.58	1.534	0.061	2.75	0.058	2.65
	Terrace	0.111	5.77	0.054	3.33	1.370	0.054	2.37	0.058	3.64
	Constr1850_1913	-	-	-0.108	-4.91	1.948	-0.099	-4.49	-0.102	-4.66
	Constr1914_1947	-	-	-0.078	-3.68	1.892	-0.071	-3.36	-0.071	-3.38
	Constr1948_1969	-	-	-0.135	-8.30	2.912	-0.129	-7.94	-0.130	-8.07
	Constr1970_1980	-	-	-0.107	-6.50	2.295	-0.100	-6.12	-0.100	-6.16
	Constr1992<	-	-	0.142	9.55	3.163	0.143	9.73	0.142	9.75
	%MidIncome	-	-	0.790	7.57	2.339	0.785	7.61	0.718	6.95
	%HighIncome	-	-	0.114	4.12	2.736	0.103	3.78	0.086	3.14
	District2	-	-	-0.090	-7.43	2.152	-0.064	-5.00	-0.076	-4.44
	District3	-	-	-0.030	-1.96	2.210	-0.025	-1.70	-0.040	-2.19
	District4	-	-	-0.107	-2.89	1.197	-0.054	-1.43	-0.031	-0.61
	Rho	-	-	-	-	-	0.324	5.58	-	-
	Lambda	-	-	-	-	-	-	-	0.881	16.42

In this paper, three areas outside of the traditional boundaries of the Lyon Urban Area are also included. A considerable number of commuters reside in these more distant areas with local centres in Villefranchesur-Saône, Bourgoin-Jallieu, and Vienne.

Institut national de la statistique et des études économiques (National Institute for Statistics and Economic Studies).

Due to the methods used by INSEE, only the estimated population is available after 1999.

Modélisation et Simulation de l'Accessibilité aux Réseaux et aux Territoires (Modelling and Simulation of Accessibility to Networks and Territories).

Using Part-Dieu as the CBD leads to a much worse regression performance.

Travel time squared or the square root of travel time as the denominator does not lead to significant changes in the result.

Hereafter, clusters are created with the method based on fuzzy equality proposed in[START_REF] Kryvobokov | Formalisation of urban land value zoning with application of GIS-model[END_REF]. The advantage of the method is that the number of clusters is not determined a priori. Indices for clusters presented in the figures are calculated as a weighted average in a cluster with areas of zones used as weights.

In Table4, t-values for global OLS and asymptotic t-values for spatial models are in parentheses. VIF always has its maximum for Year03. There are median values for GWR coefficients.

This model suggested by[START_REF] Ross | Inconsistency in Welfare Inferences from Distance Variables in Hedonic Regressions[END_REF] as an alternative to models with distance variables is tested in our OLS group. It was impossible to include in this model X-squared and Y-squared due to enormous multicollinearity.

We should admit, however, that Moran's I is highly dependent on a weighting scheme. With the rowstandardised binary weights the lowest Moran's I is 0.04 and the highest is 0.10. Nevertheless, we focus on Moran's I calculated with the row-standardised weight matrix of inverse squared distances.

For example, high quality schools are located in Les Belges and it is adjacent to Park de la Tête d'Or, which is regarded as the best urban park in France.

Acknowledgements

The study is part of project PLAINSUDD sponsored though French ANR (number ANR-08-VD-00). The author is grateful to Perval and Pierre-Yves Peguy for providing data on apartment prices and attributes and to Nicolas Ovtracht and Valérie Thiebaut for calculating the O-D matrices of travel times and the coordinates of apartments. The author also thanks Alain Bonnafous, Fabrice Marchal, Dominique Bouf, Gunther Maier and the anonymous referees for valuable suggestions and John Mc Breen for proofreading the paper. The GWR software was distributed by its developers, Chris Brunsdon, Stewart Fotheringham and Martin Charlton. The publically released GeoDa TM software by Luc Anselin was used for spatial models. The availability of both software products is gratefully acknowledged.