Miklos Molnar

Exact Algorithm to Solve the Minimum Cost Multi-Constrained Multicast Routing Problem

Keywords: Multicast Routing, Quality of Service, Multi-Constrained Steiner Problem, Hierarchy, Partial Minimum Spanning Hierarchy, Branch and Bound

The optimal solution of the multi-constrained QoS multicast routing problem is a tree-like hierarchical structure in the topology graph. This multicast route contains a feasible path from the source node to each of the destinations with respect to a set of QoS constraints while minimizing a cost function. Often, it is a tree. In other cases, the hierarchies can return several times to nodes and links of the topology graph. Similarly to Steiner problem, finding such a structure is an NP-hard problem. The usual tree and topology enumeration algorithms applied for the Steiner problem cannot be used to solve the addressed problem. In this paper, we propose an exact algorithm based on the Branch and Bound principle and improved by the Lookahead technique. We show relevant properties of the optimum hierarchy permitting efficient pruning of the search space. To our knowledge, our paper is the first to propose an exact algorithm for this non-trivial multi-constrained optimal multicast route computation. Simulations illustrate the efficiency of the proposed pruning operations. The analysis of the execution time shows that in simple topologies and with tight QoS constraints the exact algorithm requires relatively little execution time. With loose constraints the computation time cannot be tolerated even for off-line route computation. In these cases, the solution is close to a Steiner tree and heuristics can be applied. These results can serve as basis for the design of efficient, polynomial-time routing algorithms.

Introduction

Quality of Service (QoS) multicasting needs the computation of a multi-constrained multicast route connecting a source node to a set of destinations. The objective of network operators (and transitively of users) is to minimize the network resource consumption. To find an appropriate solution within the feasible solutions, a cost function can be used. The minimum cost solution can minimize for instance the hop count or an additive cost on the links. As it is discussed in [START_REF] Masip-Bruin | Research Challenges in QoS Routing[END_REF], guaranteeing QoS and optimizing resource utilization may be two conflicting interests. We are interested in finding the best cost solution with respect to a set of QoS constraints.

Without QoS constraints, the minimum cost multicast route computation corresponds to the well known NP-hard Steiner problem. In our study, we consider the constrained multicast routing with multiple QoS constraints. As it is shown in [START_REF] Mieghem | Concepts of Exact QoS Routing Algorithms[END_REF], the multi-constrained routing problems are NP-hard even if there is only one destination.

Often, in works talking about QoS aware multicast routing, authors suppose that the optimal solution is a partial spanning tree. Section 2.2 gives an overview of the most important related works. In some cases, the optimal multicast route is neither a tree, nor a set of trees, nor a set of optimal QoS paths.

In [START_REF] Kuipers | MAMCRA: A Constrained-Based Multicast Routing Algorithm[END_REF], it has been mentioned that a feasible/optimal solution of the problem may be different from a tree. Indeed, the authors suggest that it corresponds to a sub-graph containing feasible paths and eventually minimizing a length/cost function. Unfortunately, the sub-graph concept is not sufficient to define the optimal solution. Paper [START_REF] Molnar | The Cost Optimal Solution of the Multi-Constrained Multicast Routing Problem[END_REF] shows that a generalization of the tree concept called hierarchy defines accurately the optimum. Furthermore, hierarchies can also describe feasible solutions approaching the optimal one. To our knowledge, there is no algorithmic solution proposed in the literature to compute the optimal route, which is a hierarchy. In this work, we investigate on finding the optimal solution (the partial minimum spanning hierarchy) of the corresponding multi-constrained partial spanning problem. The problem is NP-hard. Since hierarchies are different from trees and do not obligatory correspond to any sub-graph, the known tree computation algorithms cannot be applied to find the optimal hierarchy. Indeed, hierarchies can contain the same graph node/edge multiple times, and the algorithms based on tree and node enumeration are not suitable for hierarchy computation.

Therefore, we propose a particular Branch and Bound algorithm permitting returns to the nodes several times.

Our paper is organized as follows. Section 2 presents the multi-constrained partial spanning problem for multicast QoS routing and provides an overview of the previously proposed approaches to solve the routing problem. Section 3 briefly presents the hierarchy concept generalizing spanning trees as well as the relevant properties of the multi-constrained minimum cost hierarchies. These properties as well as the Lookahead concept (cf. [START_REF] Mieghem | Concepts of Exact QoS Routing Algorithms[END_REF]) are useful to design the exact algorithm described in Section 4.

The computation results can be found in Section 5.

The Multi-Constrained QoS Multicast Routing Problem

Often, multimedia applications require multi-constrained multicast routes. For each destination of a given multicast group, the route should meet a set of QoS requirements such as limited delay, jitter, bandwidth, loss rate etc. Indeed, generally the source based QoS multicast routing aims to compute routes by satisfying the given QoS constraints at the destinations and minimizing the network resource consumption or cost.

Problem Formulation

The network topology is represented by an undirected graph ( )

, G V E =
, with the set of nodes V (routers) and the set of edges E (links). The source corresponds to a node s V ∈ and the destination node set is given by The multicast QoS routing aims at finding a multicast route ( )

{ } { } \ , 1
, M W F =
, where W is a multi-set of node occurrences (the visited nodes) and F is a multi-set of edge occurrences (the used edge occurrences) in the route 1 . If W and F are simple sets, the route is a sub-graph. In our constrained routing problem, M is not always a sub-graph and due to the constraints, the route can visit nodes and edges multiple times [3] [4].

Consequently, a node or an edge may be present multiple times in the multi-sets W and F respectively. We propose to refer to occurrences of nodes and edges/arcs in spite of referring to nodes and edges themselves. This multicast route ( )

, M W F = must
contain at least one feasible path ( )

, j p s d
from the source node s to each destination , 1, ,

j d j r =
. The network resource usage can be expressed by an adequate cost function. The cost may be dependent or independent from the m QoS metrics. Without loss of generality, in the following we suppose an arbitrary additive cost. Therefore, we consider the cost of the multicast communication as the total cost of the forwarded data by using the edges of the multicast structure ( )

, M W F = : ( ) ( ) e F c M c e ∈ = ∑ (2)
Since the structure M may be different from a sub-graph (the sets W and F are multi-sets with possible repetitions of graph elements), if an edge e E ∈ of G is present twice in F (because it is used twice to forward the multicast messages), its cost must be added twice to ( )

c M .
The following formulation of the optimal QoS multicast routing has been proposed in [START_REF] Molnar | The Cost Optimal Solution of the Multi-Constrained Multicast Routing Problem[END_REF]. To facilitate the analysis, we use undirected graph models despite the fact the route is directed. We consider the links can be used in both directions and the weights and cost values are the same in both directions.

Multi-Constrained Minimum Cost Multicast problem (MCMCM)

satisfies the given constraint vector L :

( ) ( ) ( ) * , , 1, , and is minimum
. ≤ = d j p s d j r c M w L (3)
As it has been demonstrated, the solution is a directed minimum cost partial spanning hierarchy [START_REF] Molnar | The Cost Optimal Solution of the Multi-Constrained Multicast Routing Problem[END_REF]. Section 3 recalls the definitions. Finding the optimum is NP-hard.

Related Works

For solving the QoS multicast routing problem, several propositions aim to compute spanning trees. The reason for basing multicast structures on multicast trees is that the tree structure avoids redundancies (thus minimizes the cost for instance).

Without constraint, the basic problem is known as the Steiner problem. Two large overviews on this well known NP-hard problem are presented in [START_REF] Winter | Steiner Problem in Networks: A Survey[END_REF] and [START_REF] Hwang | Steiner Tree Problems[END_REF]. The most known exact algorithms are based on the enumeration of the Steiner nodes 2 and Steiner trees. The solution of the Steiner problem is always a sub-graph (a tree).

If only one constraint (usually on the end-to-end delay) is given, the addressed problem is the Constrained Steiner Tree (CST) problem that is also NP-hard [START_REF] Aggarwal | Minimal Spanning Tree Subject to a Side Constraint[END_REF]. To solve it, the authors propose a Branch and Bound algorithm by using the Lagrangian Relaxation and heuristics to get lower and upper bounds in the Branch and Bound tree. A compact Mixed-Integer Program formulation with efficient reduction techniques can be found in [START_REF] Leggieri | The Steiner Tree Problem with Delays: A Compact Formulation and Reduction Procedures[END_REF]. In [START_REF] Gouveia | Modeling Hop-Constrained and Diameter-Constrained Minimum Spanning Tree Problems as Steiner Tree Problems over Layered Graphs[END_REF] the diameter-constrained minimum spanning tree variant (DM-STP) is defined. In this problem a natural number D is given, and the goal is to find a spanning tree of the graph with minimum total cost such that the unique path from any node i to any node j has no more than D edges. A Branch and Cut algorithm is proposed using heuristics to cutting plane generation. Several heuristics have also been proposed to solve approximately the problem. In [START_REF] Kompella | Multicast Routing for Multimedia Communication[END_REF], a heuristic that constructs a low cost spanning tree, with respect to a bounded delay on each multicast destination is described. In [START_REF] Kumar | The MASC/BGMP Architecture for Inter-Domain Multicast Routing[END_REF], an algorithm that approaches the minimum cost spanning tree solution with respect to the delay constraint has been presented. The authors in [START_REF] Zhu | A Source-Based Algorithm for Delay-Constrained Minimum-Cost Multicasting[END_REF] proposed a heuristic algorithm called The Bounded Shortest Multicast Algorithm (BSMA). The BSMA algorithm always finds a delay constrained tree, if such a tree exists, since it begins with the least-delay spanning tree.

When more than one QoS constraint are considered, the problem becomes more complex, and different tree based solutions are proposed [START_REF] Zhu | A Source-Based Algorithm for Delay-Constrained Minimum-Cost Multicasting[END_REF] [13] [14] [START_REF] Korkmaz | Multi-Constrained Optimal Path Selection[END_REF]. A rendezvous point based tree computation has been presented recently in [START_REF] Stachowiak | Rendezvous Point Based Approach to the Multi-Constrained Multicast Routing Problem[END_REF]. Three different formulations can be found in [START_REF] Kuipers | MAMCRA: A Constrained-Based Multicast Routing Algorithm[END_REF]. The Multiple Constrained Multicast (MCM) problem aims to compute a sub-graph { } ( )

, , M s D H =
that contains a feasible path from the source to each destination. Since the sub-graph can be arbitrarily large and the paths are not explicitly identified in this sub-graph, this formulation is not interesting. The Multiple Parameter Steiner Tree (MPST) problem searches for a partial spanning tree minimizing an arbitrary length function. Since the solution of this problem is a tree, it can be configured and used for multicasting. The drawback of this formulation is that in some cases this kind of solution does not exist. The authors say that the MPST, although optimal in terms of resource utilization, does not always satisfy the constraints. A third formulation is given by the combination of the two cited ones above; the Multiple Constrained Minimum Weight Multicast (MCMWC) problem searching for a sub-graph with minimum cost. Meta-heuristics have also been applied for solving QoS multicast routing. The papers in [START_REF] Sun | Multiple Constraints QoS Multicast Routing Optimization Algorithm in MANET Based on GA[END_REF] and [START_REF] Lu | Genetic Algorithm for Energy-Efficient QoS Multicast Routing[END_REF] propose genetic algorithm based solutions even for networks with uncertain parameters and for energy minimization respectively. The thesis work in [START_REF] Xu | Metaheuristic Approaches for QoS Multicast Routing Problems[END_REF] presents several metaheuristics (e.g., swarm optimization). An original computation of multicast trees based on a Chemical Reaction Optimization Algorithm for QoS Multicast Routing can be found in [START_REF] Sahoo | A Tree Based Chemical Reaction Optimization Algorithm for QoS Multicast Routing[END_REF]. Let us notice that these algorithms compute trees and not the optimal solution presented in our paper.

In the literature, few proposed algorithms allow solutions that are different from spanning trees. The Multicast Adaptive Multiple Constraints Routing Algorithm (MA-MCRA) [START_REF] Kuipers | MAMCRA: A Constrained-Based Multicast Routing Algorithm[END_REF] is one of the most relevant algorithms that looks for multicast routes that are different from spanning trees. MAMCRA solves the multi-constrained mul-ticast routing problem by computing a special routing structure in two steps:

1) In the first step, the algorithm computes a set of optimal (multi-constrained) paths 2) In the second step, it tries to eliminate the useless redundancies that are produced in the first step.

The result may be a hierarchy (cf. its definition in the next section).

Instead of multiple constraints, different QoS objectives can be established on the QoS values: minimizing the delay, the jitter, the packet losses, etc. A survey on multiobjective minimum spanning tree problem is available in [START_REF] Ruzika | A Survey on Multiple Objective Minimum Spanning Tree Problems[END_REF]. Our approach is based on multiple constrained formulation.

Hierarchies to Solve the Optimal QoS Multicast Routing

Usually, spanning trees are considered as the minimum cost solutions for partial spanning in graphs. However, they have some limitations. Trees can not always satisfy the end-to-end constraints and in some cases there is no tree solution for the problem.

As it has been demonstrated in [START_REF] Molnar | The Cost Optimal Solution of the Multi-Constrained Multicast Routing Problem[END_REF] the solution is a connected, graph related structure called hierarchy.

Hierarchies

Trees are connected (sub-) graphs without cycles. To solve our problem, a structure permitting the graph nodes to be visited several times has been proposed. A hierarchy is a graph-related structure defined as follows [START_REF] Molnár | Hierarchies to Solve Constrained Connected Spanning Problems[END_REF]. For multicast routing, the route is rooted at the source and we consider rooted hierarchies (the tree T is rooted). A non-empty rooted hierarchy can be considered as a connected route containing occurrences of the graph elements, where each node occurrence has at most one parent node.

Trivially, a hierarchy is not a sub-graph, but (similarly to walks) a graph-related structure. In the example of Figure 1, the nodes c and d are present twice in the hierarchy.

Since the different occurrences of the same element can play different roles, the distinction and the identification of the occurrences is substantial.

Even though a hierarchy is not a sub-graph, it generates a sub-graph in the graph G.

This sub-graph is the image of the hierarchy and may contain cycles. To facilitate the use of the hierarchies, we propose to retain some tree related terms.

A sub-hierarchy of a (rooted) hierarchy H is a hierarchy only containing elements of H.

A branching node occurrence is a node occurrence having at least a degree three in the hierarchy.

A leaf is a node occurrence with a degree one in the hierarchy. In Figure 1, there are three leaves (one occurrence of c, e and f), the second occurrence of d is a branching node occurrence and the reader can easily detect the sub-hierarchy rooted at b.

Hierarchies can be directed or undirected. A directed hierarchy can be related to an undirected graph G. In this case, an arc relies two node occurrences in the hierarchy iff there is an edge between the corresponding nodes in the graph G. If the graph G is not directed, we suppose that each edge can be replaced by a pair of arcs.

Moreover, we suppose that the weights of the edge correspond to the weights of the arcs in both directions. The directed rooted partial spanning hierarchy concept allows an accurate and exact definition of the solution of the multi-constrained partial minimum spanning problem. These rooted hierarchies are directed from the source to the destinations. As it has been proved in [START_REF] Molnar | The Cost Optimal Solution of the Multi-Constrained Multicast Routing Problem[END_REF], the optimal multicast route M, with respect to multiple constraints on positive additive metrics, is always a directed partial spanning hierarchy.

Figure 1. Example of a rooted hierarchy.

In the following, we refer to this solution as the Multi-Constrained Minimum Partial Directed Spanning Hierarchy, abbreviated by MC-MPDSH. Using the hierarchy concept, the MCMCM problem can easily be re-formulated as follows.

The MCMCM problem consists in finding the multi-constrained minimum partial directed spanning hierarchy containing at most one directed path ( ) . The cost of all edges is equal to 1 in this example. The cost minimal hierarchy respecting the QoS constraints is drawn with bold arrows. Another feasible tree also exists (drawn with dotted arrows) but this tree is more expensive.

,
Trees are special hierarchies obtained by injective homomorphism. Thus, a tree is a hierarchy. Some properties of trees are true for hierarchies but not all of them, whereas all properties characterizing hierarchies are true for trees. In [START_REF] Molnar | The Cost Optimal Solution of the Multi-Constrained Multicast Routing Problem[END_REF], relevant properties of the MC-MPDSH have been enumerated. These properties enable to design efficient exact and heuristic algorithms. Thus, we investigate on the properties of the MC-MPDSH.

Properties of the MC-MPDSH

We summarize properties that permit the design of our Branch and Bound algorithm by detecting infeasible "solutions" ASAP. The proofs of these properties are generally trivial and omitted in this paper.

1) The leaves in the optimal spanning hierarchy * M are destinations. Consequently, in an MC-MPDSH there are at most D leaves.

2) The directed path from the source s to any arbitrary node occurrence v is a feasible path.

Figure 2. The cost optimal solution of a multi-constrained multicast routing problem.

3) The edges of the topology graph G can be used multiple times and in both direction. 4) In a directed path from the source s to an arbitrary destination j d in the MC-MPDSH, any node v V ∈ has at most one occurrence (a node can not be repeated in a path, paths are elementary). This property is also true regarding the source node itself:

no directed path in an MC-MPDSH can return to the source node. Consequently, the source is present only once in an MC-MPDSH. ( )

s i s i M M v W d v l + ∈ ≤ ∑ (4)
Remember that the metrics in the graph are positive and additive. The following property enables to establish an important relation between the QoS constraints and gives a sub-optimality property. (  )

, v p s v = - L L w 4 .
10) Several paths from the source to different destinations can use an edge. We say that two paths share an edge in a hierarchy if the same occurrence of the edge belongs to both paths. Let 1 P and 2 P be two paths from the source to two distinct destinations both con- taining an edge a. The edge a is shared by 1 P and 2 P in the MC-MPDSH, iff the pre- fix of 1 P and 2 P from the source to the edge a is the same:

( ) ( )

1 2 P P pref a pref a = .
Figure 3 illustrates a graph, that contains only one feasible path 1 P and another 2 P from the source s to the destinations 1 d and 2 d respectively. In the minimum cost spanning hierarchy, the two paths can share the edge { } , s b (the prefixes are the same, i.e. this prefix is empty), but they cannot share the edge { } This property corresponds to the well known Bellman's principle of optimality [START_REF] Bellman | Dynamic Programming[END_REF] and permits to search the optimal solution with dynamic programming. Unfortunately, the computation of all the smallest and possible sub-hierarchies is expensive. Remember that the problem is NP-hard.

MC-MPDSH Computation by Branch and Bound Algorithm

The multi-constrained multicast routing problem is NP-hard. Unfortunately, the algorithms proposed to solve the Constrained Steiner Problem such as the Steiner Tree Enumeration Algorithms and the Topology Enumeration Algorithms cannot be applied here, because of the allowed multiple occurrences of the graph elements in the solution.

In the following, we present a Branch and Bound algorithm that exactly solves the MCMCM problem. The efficiency of this kind of algorithms depends strongly on the applied pruning operations. In our algorithm, two ways are proposed to reduce the search space and accelerate the computation:

• pruning based on properties of the possible optimum as presented in the previous section.

• pruning based on the Lookahead concept.

We will demonstrate that the proposed pruning operations accelerate significantly the computation.

Main Algorithm

Despite the fact that there are repeated nodes and edges in a hierarchy, the Branch and Bound algorithm can easily be adapted for the spanning hierarchy computation. The proposed framework (cf. Algorithm 1) is a frontier search type Branch and Bound algorithm. It allows the enumeration of the partial spanning hierarchies, which are candidates to lead to the solution. Each node in the virtual search tree corresponds to a hierarchy rooted at the source node. Generally, these hierarchies are sub-hierarchies, potential feasible germs of the solution spanning the destinations. The frontier is the set of leaves in the search tree. It is initialized with a hierarchy that contains only the source node s. At each step of the algorithm a simple greedy strategy is applied: the hierarchy H with the best cost is selected from the frontier. The first hierarchy in the course of the enumeration, which covers all destinations and satisfies the end-to-end constraints corresponds to the optimal solution and is returned. If the destinations are not yet covered by the selected hierarchy, possible continuations of this selected hierarchy H are computed to create larger hierarchies. These successors are then added to the frontier. The particularity of the problem is that repeated elements may occur in the continuations. If there is no hierarchy in the frontier satisfying the end-to-end constraints, the problem has no solution. The first steps of the search tree exploration is illustrated in Figure 4.

Let K be the frontier (the set of leaves in the search tree). The meta-code of the main algorithm for computing the MC-MPDSH is given by Algorithm 1.

The brute enumeration of all possible hierarchies is very expensive. However, the properties of the optimal solution and the Lookahead technique permit to realize pruning operations. We analyze the efficiency of these accelerations in the following section. The crucial question is how the valid successors of a selected hierarchy H are computed to continue the enumeration of hierarchies in the search tree. Duplication of spanning hierarchies should be avoided and all possible candidates should be enumerated.

Computation of the Successors

To enumerate the candidate hierarchies, we propose the following construction. Like trees, hierarchies can always be decomposed into layers. In our case, we consider that the source corresponds to layer 0. When new edges are added to a selected hierarchy in

K ≠ ∅ do select H from K with minimal cost ( ) c H if D H ⊂ then return H {it is the optimal solution} S ← successor_hierarchies_of( H ) {cf. Algorithm 2} for all s H S ∈ do {( )} s K K H ← ∪ {add s H to K } end for \ ( ) K K H ← {delete H from K }
end while STOP without solution the algorithm, these edges are added always at the highest layer of the hierarchy. This layer contains only leaves. If the selected minimum cost hierarchy have n layers, the new successor hierarchies will have 1 n + layers. Following the example in Figure 4, the successors of the root hierarchy (containing only the source) are hierarchies having two layers (and never three). The successors of the hierarchy selected in the example are hierarchies having three layers and the new edges are added only to the leaves of layer 1.

Let H be the selected hierarchy and let H L be the set of the leaves in the last layer of H. The successors of H in the search tree are the hierarchies enlarged by adjacent edges from H L . We call this kind of edges fringe edges, expression borrowed from [24].

Since a node can not be repeated in any path (Property 4), edges returning to the parent node of a leaf are not in the fringe edges. The fringe edges of a hierarchy are illustrated by Figure 5.

The computation of the set of successors of a given hierarchy H is described by Algorithm 2. Only the hierarchies, which correspond to the end-to-end constraints and to the properties of the optimal solution and which are not excluded by the Lookahead (cf. Subsection 4.3) can be sub-hierarchies of any solution. To simplify, we call them valid candidate hierarchies. The algorithm proceeds in two steps.

In the first step, it selects the set of fringe edges (indicated by A) of the hierarchy H. (Non empty combinations of these edges can be added to H to create successors in the second step.) Trivially, if a fringe edge f leads to a new leaf that violates the QoS constraints, f is not added to the set A (the QoS constraints must be respected in the new leaves of the successor hierarchies, cf. Properties 2 and 9). Moreover, a node can not be repeated twice in a path from the source in the optimal solution (Property 4). Consequently, a fringe edge from a given leaf a is not added to the candidate set A if its other extremity already belongs to the set of parents of a.

t ← opposite_node_of ( , ) n f ( ) ( ) ( ) w t w n w f ← + {accumulated weight at t } if ( ( ) ) & ( ) n w t L t V ≤ ∉ then { } A A f ← ∪ end if end for end for
Step 2: Combination of the fringe edges computed in Step 1 for all combination ( )

C A of the edges in A do ( ) H H C A ′ ← ∪ {add the edges of ( ) C A to H } ( ) a a DR D i ← ∪ {set of reachable destinations} if (number_of_leaves ( ) | | H D ′ ≤ and max-degree_of_nodes ( ) | | /2 H D ′  ≤  and DR D ⊆ ) then { } S S H′ ← ∪ {add H ′ as new successor to S} end if end for
In the second step, the algorithm combines the selected fringe edges (enumerated in the first step) for construction of valid successors. Since the optimal hierarchy can have at most D leaves, only successors corresponding to this limitation are enumerated (Property 1). Following Property 7, the branching node occurrences of a given node cannot exceed 2

D      
. Property 8 permits to verify the sum of the degrees of the occurrences of the nodes present in a new configuration. Moreover, if a new combination does not permit to reach all the destinations (see also in the next sub section), the combination is deleted and not added to the successor set.

In the worst case, the algorithm constructs at most ( )

min 2 1, 2 1 D k - - valid candi- date hierarchies,
where k is the number of fringe that are added to A in step one.

Many other not valid hierarchies (which do not match any solution) can be detected using the Lookahead technique. In the following we summarize the adaptation of the Lookahead technique in our case.

Application of the Lookahead

The Lookahead is originally used in Artificial Intelligence [25] [26] to accelerate the computation by applying upper bounds and reducing the search space. The concept was successfully applied in the computation of optimal paths under multiple constraints [START_REF] Mieghem | Concepts of Exact QoS Routing Algorithms[END_REF]. In our case, in a node reached on a path from the source, one can consider the remaining information of any path toward the destinations and decide if this node can be usefull and considered for the computation or not. More precisely, the Lookahead concept can be applied in the computation of the optimal hierarchy as follows.

At 2) Only one m-dimensional vector

( ) v LV
is considered in each node v with the smallest values from this node to the closest destination concerning each metric. For this vector, ( ) ( )

min j j i d D i LV v LV v ∈ =
. This vector is constant and cannot be reconsidered.

Trivially, the first method consumes more memory but it is more precise and efficient for pruning. In this study, we implemented the simpler second solution: only one Lookahead vector per node was used.

It is important to state that the computation of the Lookahead vectors is held off-line, before the construction of the candidate hierarchies. Then these vectors can be used at the computation of valid successors. For each fringe link selection (first step of Algorithm 2), the following additional test can be performed.

• Let f be a fringe edge from node a in the hierarchy to a newly examined node b.

Trivially, if

( ) ( ) ( ) , d j p s d b + ≥ w LV L
, there is no further destination which can be reached from b respecting the constraints. The fringe arc can not be used to create successors (there is no feasible path to any destination).

In Step 2, the algorithm combines these edges and generates new candidates. The usefulness of each combination can be tested by using the Lookahead information, if the Lookahead vector ( )

j v LV is available for each destination j d D ∈ . • Let { } , 1, , k B b k l = =
be the set of leaves of the candidate hierarchy. (B contains the extremities of the newly added fringe edges and also leaves, which are not concerned by the new fringe edges.) Let k D be the set of "reachable" destinations from k b . For a "reachable" destination, (

,

d j k j p s b d + ≤ w LV L 5 .
For each combination, the set of destination which can be achieved by the combination can be computed. If the set does not cover all the destinations, this new hierarchy combination is not added to the candidate set.

Exactness of the Proposed Branch and Bound Algorithm

It is important to state that Algorithm 2 computes all valid hierarchies without skipping any of these hierarchies, even if it considers the proposed pruning operations.

Lemma 1. Algorithm 2 does not create duplication.

Proof. Remember, the successors are created following the layers of the hierarchies.

Starting from a hierarchy 1 H selected by Algorithm 1, non-empty combinations of fringe edges are added to create successors. A given combination is added only once to 1 H . Two identical successors can not be produced from the hierarchy 1 H . Proof. Let us suppose that there is a hierarchy H ′ with less cost than the optimum and this hierarchy is skipped. By deleting the fringe edges from H ′ , the hierarchy H ′′ is obtained, which is also valid and of less cost. Since H ′ is not enumerated, H ′′ can not be enumerated because there is no filter from H ′′ to H ′ . In this way, recursively we can arrive in the route hierarchy, which can not be enumerated by the algorithm. It is absurd.

Theorem 1. The proposed Branch and Bound algorithm returns the MC-MPDSH if it exists.

Proof. The algorithm enumerates all the valid candidate hierarchies before arriving at the solution (Lemma 2). The valid candidates are selected by Algorithm 1 in an increasing order of costs. If a hierarchy H is selected, then there is no solution H ′′ with less cost than the cost of H (otherwise H ′′ is selected to examine). The first selected hierarchy spanning all the destinations is the MC-MPDSH.

Performance Evaluation of the Branch and Bound Algorithm

To analyze its performance, the algorithm have been executed in two benchmark networks. The first is the DARPA CORONET CONUS topology with 75 nodes and 99 links (illustrated by the first topology in Figure 7). The second is the NTT (Nippon Telephone Telegraph of Japan) network topology. This network contains 55 nodes and 144 links (the second topology in the figure).

In the presented simulations, three or four metrics were associated with each link and the link values were randomly generated from the same integer set (for example from { } 1, 2,3, 4,5 ) for each metric. We also considered the cost as a random value (corresponding to the first metric). Different constraint vectors are generated from tight to loose ones. For a given destination d, we can consider the tightness of constraints as follows. For each metric i, a shortest path i p with length ( ) i i l p can be computed. This length is a lower bound for the feasible solutions from the source to this destination regarding this metric. A vector of constraints can be considered as tight, when the constraints are close to the upper bounds but there is at least a feasible solution. Contrarily, when the constraints are far from the upper bounds, feasible solutions are easier to find, we consider the constraint vector as loose.

The execution time of the algorithm is quasi proportional to the number of computed hierarchies in the search tree. In this study, we use this number as a fundamental performance metric. A second metric can be the number of iterations before obtaining the solution.

Relevance of the Pruning and Lookahead

Due to the possible repetitions of nodes and edges in the spanning hierarchies, their enumeration by the Branch and Bound exact algorithm is more expensive than the enumeration of the more restricted spanning trees. However, the presented properties of the optimal hierarchy allow to implement efficient pruning. In addition to the pruning operations, we expect that the application of the Lookahead concept will also reduce considerably the search space.

To measure the impact of both pruning and the Lookahead technique, three variants of the algorithm have been executed: 1) BBB: the basic Branch and Bound algorithm without pruning and Lookahead 2) BBwP: the algorithm improved by pruning based on the properties of the solution 3) BBwPLA: the algorithm with both the pruning and Lookahead operations.

At first, we executed the three algorithms in CORONET network for a multicast group of five arbitrarily chosen members. With randomly generated link values and moderately tight constraints 6 ten executions were launched, and only the first 20000 iterations of the algorithms were registered. Significant differences between the performances were observed. Figure 8 illustrates the progression of the algorithms showing the average number of computed hierarchies in the ten executions and the corresponding table indicates the average number of the computed hierarchies at the end of the 20,000 iterations.

The results show the efficiency of the propositions. In the example, the proposed pruning operations and the Lookahead technique divide the number of enumerated hierarchies in the search tree by more than ten.

The Impact of the Tightness

The value of the constraint vector influences strongly the computation time. The tightness of the constraints is relative to metrics and destinations. With a given constraint vector and for some destinations, one can find several QoS routes satisfying the constraints whereas it is not possible for some other destinations. Since the solution (the spanning hierarchy) should cover all the destinations, we consider that the tightness always concerns the most critical destination. Intuitively, when the constraints are loose (for all destinations), the optimal solution is a route of minimum cost but practically without constraint. It converges to the minimum Steiner tree. We will see that our Branch and Bound algorithm (which permits largely node repetitions) is not efficient to compute usual spanning trees (which are without node repetitions).

The enumeration of the possible hierarchies with potential node repetitions becomes essential, when the constraints are tight. In these cases both the Lookahead mechanism and the pruning operations are more efficient. The following experiences illustrate this effect. Let us notice that the computation of the exact solution is expensive. We developed only a few examples but the results are qualitatively similar in these examples.

In the DARPA CORONET CONUS topology, we tested the algorithm (the proposed BBPwLA version) for an arbitrarily configured multicast session with different constraint vectors. The (randomly selected) source was in Austin and five destinations were (also randomly) selected at nodes Abilene, Buffalo, Los Angeles, Billings and El Paso. In this computation, three metrics were associated with each link and the values were generated randomly from the integer set { } 1, 2,3, 4 . The cost corresponded to the first metric.

We computed the optimal solution for several QoS requirements. Since the link values were generated from the same range, we used the same end-to-end constraint for the three metrics. Table 1 resumes the results on the number of iterations and com- The progression of the computation is illustrated by Figure 9. This experience shows that the computation becomes longer with loose constraints (curves are stopped quickly in tight environment).

Why not compute the optimal solution with tighter constraints, if the computation is faster? Unfortunately, the response is negative. Tight constraints do not always implicate the same quality and cost that can be obtained with loose constraints. Intuitively, loose constraints permit to compute a solution near the cost minimal Steiner tree but at the expense of the QoS. With tighter constraints (satisfying tighter quality requirement) the optimal solution can have higher cost. This phenomenon is illustrated by the following example. In the CORONET CONUS topology for the same multicast group but using differently generated random link values, the following results were obtained. In this case, there is no solution for the constraint vector [ ] T 16,16,16 and the vector [ ] T 17,17,17 can be considered as tight constraint vector. We computed the exact solution with four constraint vectors as it is shown in Table 2. Figure 10 shows the progression of the algorithm in the different cases. For tight constraints, the optimal solution is obtained rapidly (computing 14,334 hierarchies) but he cost is higher (49) than the cost (37) obtained in the case of loose constraints. The solutions are different.

To obtain the optimal solution for the constraint vector 3).

Figure 9. The progression of the number of computed hierarchies in the first example. With looser constraints the cost diminishes as it is expected. The particularity of this case is that with the constraints [ ] T 14,14,14,14 the number of iterations (the number of the selected hierarchies from the front of the search tree) is less (17,527) than with the constraints [ ] T 13,13,13,13 (24,039), but the number of computed successor hierarchies is higher. The loose constraints implicate less pruning operation and in this way, the number of computed successors increases rapidly. But, as it illustrated in this example, the computation can be terminated using less iterations. The progression is shown by Figure 11. Our summary concerning the results is as follows.

The investigated simulations corroborate our intuitions. In the case of loose constraints, the exact algorithm is very expensive, even if the pruning operations limit the search space. In this case, the solution is close to the Steiner tree. The enumeration of the candidate hierarchies can not be tolerated in real route computations. Steiner heuristics can offer solutions. Contrarily, in the case of tight constraints, due to the pruning, the number of computed candidate hierarchies is lower and the proposed algorithm or a somewhat modified variant of this algorithm can be used.

Conclusions and Perspectives

The exact solution of the multi-constrained QoS multicast routing problem corresponds to a directed rooted hierarchy. However, finding feasible and/or minimum cost Table 3. The number of computed hierarchies for different constraints in the third example. hierarchies are NP-hard problems. Exact algorithms to compute the optimum have not yet been proposed. The algorithms computing heuristic solutions manipulate trees and sets of paths rooted at the source. We argue that it is essential to identify the optimal solution of the addressed problem. It is clear that this solution does not obviously belong neither to the set of minimum cost paths nor to the set of shortest paths computed by using the earlier proposed and known non-linear length. Since hierarchies may contain multiple occurrences of a node or an edge/arc, most of the enumeration algorithms are not appropriate to compute the optimal solution. The main result of our investigations in this paper is an exact Branch and Bound algorithm to compute the optimal hierarchy. We proved that the proposed algorithm finds the optimum if it exists. To accelerate it, the design of the algorithm is based on some important properties of the hierarchy-type solutions. Moreover, the Lookahead concept is applied successfully. The simulation's results highlight the gain on the computation time, when the algorithm considers the proposed pruning operations. Due to the complexity of the problem, the exact algorithm cannot be applied for routing in large networks. Your results illustrate that in the case of loose QoS constraints, the solution can be a tree and in this manner modified Steiner heuristics can be applied. In the case of tight constraints, the number of computed candidate hierarchies is limited and we can hop the design of efficient quasi-exact algorithms. These first results on the optimum are good start points for the future design of routing algorithms.

≤

  end-toend QoS requirements, expressed as constraints from the source to the destinations, are fixed by an m-dimensional constraint vector is the Pareto dominance.

  This problem deals with finding the structure ( ) , M W F = with a minimum cost ( ) c M , containing at least one path ( ) , j p s d from the source node s to each destination j d D ∈ that 1

  de- fines a hierarchy in G.

3 A

 3 homomorphism preserves the adjacencies of nodes; 1 w W ∈ and 2 ∈ w W can be adjacent if the corresponding nodes 1 ∈ v V and 2 ∈ v V are also adjacent.

Figure 1

 1 Figure1illustrates a hierarchy. In our case, the graph G is the topology graph of the network and the tree

Figure 2

 2 Figure2presents the optimal directed hierarchy solution for two destinations (nodes h and i) when two weights are associated with edges and the QoS requirements correspond to

M

  be an MC-MPDSH rooted at s. It satisfies the constraints L in the lea- ves. Let v M be a sub-hierarchy of * M rooted at v s ≠ght vector of the path ( ) , p s v in * M . The sub-hierarchy v M is an MC-MPDSH from v to the destination occurrences that it contains, with respect to the constraints ( )

,

  f g . (The edge { } , f g is present twice in the hierarchy, as it is illustrated by the associated T in the figure).This study of the properties of the optimal solution is very useful to design exact algorithms and heuristics. Our proposal for an exact Branch and Bound algorithm follows.

4

 4

Figure 3 .

 3 Figure 3. Shared and unshared common edges in the optimal hierarchy.

Figure 4 .∈

 4 Figure 4. An example for the evolution of the search tree in a simple topology.

Figure 5 .

 5 Figure 5. The fringe edges of a hierarchy.

1 )

 1 first, we compute the m shortest path trees from each destination to all other nodes by considering the m metrics. m D * shortest path trees (computed in polynomial time) are obtained. To store the different computed values at each node, two propositions can be considered: All values from this node to the destinations are saved in D Lookahead vectors (initially m D * values at each node). In the following, we denote by ( ) j v LV the Lookahead vector computed for the node v regarding the destination j d is the shortest path from v to j d regarding the metric i. When a destination is reached, the vectors are updated (one can remove the corresponding Lookahead vectors).

Figure 6

 6 Figure 6 illustrates the fist Lookahead test. The source is the node s and two destinations are concerned (nodes k and m). Two metrics are used and the link values are indicated. The computed Lookahead vectors are also indicated with bold characters near the nodes. During the first iteration, there are three candidates to create eventual successors of the source ( ) , s g , ( ) , s h and

Figure 7 .

 7 Figure 7. Topology CORONET CONUS [27] and NTT [28].

Figure 8 .

 8 Figure 8. The number of computed hierarchies in the first 20,000 iterations in an example.

  927 hierarchies was needed. A third computation case study was performed in the NTT network topology. The link values were generated randomly between 1 and 5 and for four QoS metrics. A multicast group of four destinations (nodes 8, 17, 25 and 28) and a source (node 14) were also randomly selected. For the constraint vector [ ] T 11,11,11,11 there is no solution in the generated case, and the tight constraints permitting a feasible solution correspond to [ ] T 12,12,12,12 . With different feasible constraint vectors, the following solutions were obtained (cf. Table

Figure 10 .

 10 Figure 10. The evolution of the number of computed hierarchies in the second example.

Figure 11 .

 11 Figure 11. he evolution of the number of computed hierarchies in the NTT example.

  ) contains the fringe edges. By deleting these fringe edges, the previous hierarchies can be obtained. These previous hierarchies are the same by deleting the same fringe edges both from 1 H ′ and from

	Let us suppose that there are two identical hierarchies 1 H ′ and	2 H ′ produced from
	two different hierarchies	1 H and	2 H . The last layer of	1 H ′ (and also in
		2 H ′ . The contradiction is trivial.
	Lemma 2. Algorithm 2 enumerates all valid candidate hierarchies without skipping
	any of them.
			.

6. Lookahead gain when [ ] T 10,10 = L 2 H ′

Table 1 .

 1 The number of computed hierarchies for different constraints in the first example.

	constraints	cost of the opt	# of iterations	# of computed hierarchies
	[	] T 21, 21, 21	no solution	776	1279
	[	22, 22, 22	] T	38	6213	14,334
	[	] T 23, 23, 23	38	6213	14,334
	[	24, 24, 24	] T	38	22,409	78,206
	[	25, 25, 25	] T	38	22,409	78,206
	[	26, 26, 26	] T	38	72,615	267,228
	puted hierarchies to obtain the optimal solution. In the observed case, there is no
	solution with the constraints [	]
	For instance, the solution with minimum cost for the constraint vector [	26, 26, 26	] T
	(with cost 38) was obtained after 72,615 iterations and computing 267,228 hierarchies.
	The solution with the same quality (with cost 38) was also obtained with the tighter
	constraints [	]

T 21, 21, 21 . In this way, we can consider that the vector [ ] T 22, 22, 22 corresponds to a tight constraint vector and the higher values are looser. T 22, 22, 22 after 6452 steps in the Branch and Bound algorithm computing only 14,334 hierarchies.

Table 2 .

 2 The number of computed hierarchies for different constraints in the second example.

	constraints	cost of the opt	# of iterations	# of computed hierarchies
	[ 17,17,17	] T	49	6213	14,334
	[ 18,18,18	] T	42	17,152	20,907
	[ 19,19,19	] T	37	47,574	75,064

A Steiner node is a node that has a degree greater than two in the spanning tree.

It is not sure that a "reachable" destination can be reached using the proposed test. What is sure is that it cannot be reached if the condition is not satisfied.

The selected constraints were neither tight nor loose.

Acknowledgements

I thank Dr. Samer Lahoud, associate professor at the University of Rennes 1, who has worked with me on an earlier version of the proposition. I thank also Dr. Fernando Kuipers, associate professor at Delft University of Technology, who proposed the application of the Lookahead concept. Finally, the algorithm has been implemented by Dimitri Morere, master student at the University of Montpellier.


Giulio Bottazzi
email: giulio.bottazzi@sssup.it

Ugo M Gragnolati
email: ugo.gragnolati@univ-paris1.fr

Fabio Vanni
email: fabio.vanni@sssup.it

Non-linear externalities in firm localization

Keywords: C12, C16, C51, R30 Firm localization, Externalities, Non-linearities

published or not. The documents may come

Introduction

The uneven spatial distribution of some scarcely mobile factors of production may also drive the firms that use such factors to mimic their spatial distribution. For example, extractive industries may be clustered around mines, gushers, or gas fields, while being virtually absent where these natural resources are unavailable. Many other industries, however, are also strongly concentrated in space despite not being bound to immobile factors of production. In these cases, clustering is generally explained through the self-reinforcing dynamics stemming from different types of interdependencies in the localization choices of firms. In particular, these interdependencies can be categorized according to the role that market mechanisms play in their unfolding (see [START_REF] Scitovsky | Two concepts of external economies[END_REF]. When the choice of a firm enters the production problem of other firms by affecting the prices thereby involved, then the interdependence among firms is classified as a pecuniary externality. Otherwise, when the choice of a firm affects directly and exclusively the factor inputs of other firms, then the interdependence is characterized as technological externality. There is a strictly practical reason for why this taxonomy is of interest especially when the interdependencies at stake are spatially bounded. Namely, if such interdependencies stand at the root of the self-reinforcing dynamics that allow the formation of spatial clusters, then policy-makers may well be interested in unleashing localized externalities so as to promote regional development. Yet, the policies that may be thought to govern pecuniary externalities differ from those that could target technological externalities (see [START_REF] Martin | Paul Krugman's geographical economics and its implications for regional development theory: a critical assessment[END_REF]Sunley, 1996, Ottaviano, 2003).

Pecuniary externalities can lead to self-reinforcing dynamics and thus to spatial agglomeration through a process of cumulative causation based on the accumulation of local demand (see [START_REF]Economic Theory and Underdeveloped Regions[END_REF]. As an example, take two perfectly symmetrical regions. If one firm moves from one region to the other, local wages decrease in the region of departure and increase in the region of destination, thus leading workers to move in the same direction as the firm. In turn, workers are also customers, so that consumer demand rises in the destination region. In this sense, the localization choice of a single firm ends up affecting all other local firms through market demand, thus constituting a pecuniary externality. This mechanism may keep on attracting additional firms as long as the gain from an expansion in local demand outstrips the loss associated with fiercer local competition. Hence, under suitable transportation costs, imperfectly competitive markets, and economies of scale, the initial relocation of one firm calls for other firms to move in the same direction and gives rise to spatial concentration (see Krugman, 1991a). In this framework, the geography of production can be characterized in terms of core and periphery. This is especially the case since the attractive pull of demand acts across sectors, thus suggesting that the resulting agglomeration should look like a diversified city rather than as a specialized cluster.

Yet, specialized clusters do exist in areas that are neither particularly populated nor especially well-connected to markets. Only a limited part of these agglomerations can be explained through the dependence of one sector on some immobile factor inputs, as the contribution of broadly-meant natural advantages to spatial concentration is empirically modest across sectors (see [START_REF] Ellison | The geographic concentration of industry: Does natural advantage explain agglomeration?[END_REF]. These specialized clusters could then stem from pecuniary externalities that take place with the co-localization of firms belonging to vertically-integrated industries, as discussed by [START_REF] Venables | Equilibrium locations of vertically linked industries[END_REF]. Also in this case, however, the empirical evidence suggests that the co-agglomeration between industries with strong upstream-downstream ties is limited (see [START_REF] Ellison | Geographic concentration in US manufacturing industries: a dartboard approach[END_REF]. Therefore, something else is likely to be at play within sectors.

Sector-specific technological externalities represent an alternative source of self-reinforcing dynamics that may contribute to explain the riddle of the spatial distribution of firms. A typical example of how these untraded interdependencies may unfold locally is represented by the case of knowledge spillovers. That is, some of the private productive knowledge incorporated in firms spills into a common knowledge pool whose availability is bounded in space, due to the fact that knowledge is at least partly tacit and thus not perfectly transferable. In these circumstances, firms have an incentive to localize where a larger share of other firms are already settled, so as to benefit from the vaster knowledge pool that is available only locally (see Marshall, 1890, bk IV, ch.X). What makes this mechanism sector-specific is that firms are interested in acquiring knowledge that is related to their own production process, so that they will co-localize with similar-rather than with generic -firms. In this sense, localized technological externalities stemming from knowledge spillovers can make sense of specialized clusters.

On the previous special issue of this journal, [START_REF] Bottazzi | Cities and clusters: economy-wide and sector specific effects in corporate location[END_REF] provide a methodology to disentangle and measure the determinants of the spatial distribution of firms. Their work focuses in particular on singling out and comparing the strength of localized externalities acting within sectors against the pull of intrinsic regional characteristics. The idea is that if a particular sector is very concentrated where other sectors are not, this may stem either from particular local advantages or from sector-specific externalities. Once the former are controlled for, the latter can then be properly measured. In particular, the intrinsic attractiveness of locations is defined by a number of controls possibly acting across sectors, such as local demand, industrial variety and infrastructural advantages; yet, the controls also include a sector-specific dummy indicating the presence of an industrial district. In this way, the measurement of within-sector externalities is meant to be cleared from those pecuniary effects that may act across sectors -such as consumer demand -and from those that may be at least partly sector-specific -such as specialized labor market pooling. Hence, the type of effect that remains to be captured on the side of localized externality lends itself to be interpreted as a technological externality. [START_REF] Bottazzi | Cities and clusters: economy-wide and sector specific effects in corporate location[END_REF] base their analysis on Italian plant data for a variety of manufacturing and service sectors as observed in year 2001 at the scale of commuting zones. Their main result is that the pull of externalities and population size are by far the most important determinants of firm localization, and these two drivers are comparable with each other in terms of their magnitudes. In fact, externalities are found to have an even stronger effect than population size across most sectors. If one interprets sector-specific externalities as being mostly technological, such a result contrasts to some extent with the limited weight that technological externalities have received, for instance, in the context of the earlier coreperiphery models (see [START_REF] Fujita | The spatial economy: Cities, regions, and international trade[END_REF]. For example, Krugman (1991b, p.54 and pp-61-62) states:

[W]hile I am sure that true technological spillovers play an important role in the localization of some industries, one should not assume that this is the typical reason-even in the high technology industries themselves [p.54].

[...] An accident led to the establishment of the industry in a particular location, and thereafter cumulative processes took over. [...] What the historical record shows us are two things. First, such cumulative processes are pervasive [...]. And second, Marshall's first two reasons for localizations, labor pooling and the supply of specialized inputs, play a large role even when pure technological externalities seem unlikely [pp.61-62, italics added].

By contrast, the results discussed in [START_REF] Bottazzi | Cities and clusters: economy-wide and sector specific effects in corporate location[END_REF] suggest that technological externalities may hardly be considered as "unlikely". In fact, they should be regarded as the rule in the economy, rather than as an exotic exception concerning only high-tech sectors (see [START_REF] Bottazzi | Gobalizing knowledge: How technological openness affects output, spatial inequality, and welfare levels[END_REF] for a model investigating the role of technological externalities in the context of New Economic Geography core-periphery framework).

There is, however, at least one reason to suspect that such a conclusion would overestimate the strength of technological externalities in determining firm localization. Namely, [START_REF] Bottazzi | Cities and clusters: economy-wide and sector specific effects in corporate location[END_REF] base their estimates on a linear specification of externalities. In their framework, the individual advantage for a firm to choose a location increases proportionally to the number of other firms in the sector that are already placed in the same location. Externalities, however, may well be non-linear. In its "weaker" declination, non-linearity would still preserve monotonicity. In the case of localized knowledge spillovers, for instance, the advantages from having an additional neighbor may increase up to a critical threshold and then keep stable. In this case, externalities are non-decreasing and thus monotone, but they are nonlinear. In a "stronger" declination of non-linearity, instead, monotonicity is lost. For example, the advantages of having an additional neighbor may increase up to a critical threshold and decrease beyond it, for instance because firms may incur into growing search costs to select and exploit the relevant pieces of productive knowledge. In this case, firm localization is subject to congestion and externalities are non-monotone (see [START_REF] Fagiolo | Endogenous neighborhood formation in a local coordination model with negative network externalities[END_REF], for an interpretation based on endogenous network formation).

The present work extends the localization model originally presented by Bottazzi et al. (2007) by allowing for quadratic externalities. In carrying out such an extension, two key characteristics of the original framework are entirely preserved. First, it remains possible to disentangle localized externalities from the effect of other location-specific variables. Second, the localization choices of firms converge in the long-run toward an ergodic invariant distribution. In this sense, the model prescribes the same stochastic equilibrium regardless of the initial distribution of firms across regions. On the other hand, the introduction of non-linear externalities prevents us from deriving the equilibrium distribution in closed form. Hence, we apply numerical simulations to derive the equilibrium distribution of firms across regions. Such a theoretical distribution is then compared with the observed one through χ 2 minimization, so as to estimate its unknown the parameters. More specifically, the main objective of this procedure is to test whether the quadratic externality coefficient is statistically different from zero. If not, localized externalities would result that are sufficiently well-approximated by a linear form; otherwise, they are better approximated by a quadratic shape. Notably, depending on the sign and magnitude of the estimated quadratic coefficient, the model allows one to accommodate both monotone and non-monotone shapes of externalities. If the quadratic term is estimated to be statistically different from zero and sufficiently negative, a region can reach the point in which the addition of one extra firm decreases the probability for that same region to attract others firms.

Model

The localization choices of firms depend on the intrinsic features of regions as well as on the distribution of other firms across regions. A class of stochastic models which is naturally suited to capture both dependencies is that of "generalized" Polya urn schemes. These models have been applied, for instance, to the description of phenomena like technological adoption and diffusion (see the early contributions in [START_REF] Arthur | Path-dependent processes and the emergence of macro-structure[END_REF] and [START_REF] Dosi | Generalized urn schemes and technological dynamics[END_REF]). In those applications, however, the "fixed effects" provided by the intrinsic features of the objects of choice tend to progressively disappear, be it in time or space. They are as such less useful to describe cases of persistent fluctuations. Moreover, they often present non-ergodic dynamics, which imply a substantial impossibility to estimate their parameters from the data. For these reasons, Bottazzi and Secchi (2007) and [START_REF] Bottazzi | Sectoral and geographical specificities in the spatial structure of economic activities[END_REF] modified the original urn scheme framework toward a Markov process known in the physics literature as Ehrenfest-Brillouin (E-B) model. The E-B model can be obtained under rather general conditions from discrete utility theory, assuming a linear externality effect very similar to the one adopted by Brock and Durlauf (2001a,b). The equilibrium distribution entailed by the E-B model is in fact equivalent to a finite time Polya process; but, being ergodic, it can be directly estimated from the data.

Notably, ergodicity allows one to overcome the concerns on identification presented by [START_REF] Blume | Identification of social interactions[END_REF]. The present paper introduces an extension of the E-B model to include quadratic externalities while still allowing one to estimate them in high dimensional problems.

A population of N firms has to choose among L regions, and each generic firm i chooses location l with probability

p l ∼ a l + bn l + cn 2 l , (1)
where n l is the number of firms already present in l, a = (a 1 , . . . , a L ) is an L-dimensional vector assumed with positive components and b ≥ 0 and c are scalars. The probabilistic outcome of the choice is meant to capture the heterogeneity of firms and the consequently idiosyncratic nature of their preferences. One has to think of these probabilities as an average across all different firms composing the economy. Expression (1) has two distinct components: on the one hand, terms such qs a l capture the intrinsic attractiveness of the different locations; and, on the other, the term n l establishes an interdependence between the choice of firm i and those made by other firms, thus capturing externalities. When c = 0, the probability in ( 1) is linear in n l and the invariant distribution of the model is the multivariate Polya distribution (see Bottazzi and Secchi, 2007). In this case, the existence of a stochastic equilibrium can be demonstrated either by proving the attainment of a detailed balance condition or by studying the convergence of sequential choices (see [START_REF] Scalas | Statistical equilibrium in simple exchange games[END_REF]. But when c = 0, (1) is not linear and one cannot analytically exploit the detailed balance condition to obtain the explicit expression for the invariant distribution. Therefore, the second approach is used instead by assuming a sequence of choices and then studying numerically their convergence properties.

The sequence of individual choices is structured so that, at each time step, one firm is selected at random to revise its choice. In general, alternative non-random rules could be adopted to select who is called on to operate a revision, thus affecting the dynamics of the model. In the present case, however, the aim is to keep the structure of selection as agnostic as possible precisely by attributing to all firms an equal probability to be selected for choice revision. Under this premise, the evolution of the system from configuration n at time t to configuration n at t + 1 with n m = n m -1 , n l = n l + 1 and n k = n k for any k = l, m is defined in terms of the transition probability P {n t+1 |n t }. Such probability corresponds to the intersection between the event "Firm revises its previous choice of region m" (event B) and the event "Firm chooses region l" (event A). That is, the conditional probability reads P {n t+1 |n t } = Pr{Firm revises its choice m}•

• Pr{Firm chooses region l | Firm revise its choice m} .

Since the firm that is called on to revise its current choice m is selected at random, it follows that Pr{B} = n m /N . This probability must be then multiplied by the probability p l to select region l, as given in (1), conditional to the fact that the firm is no longer among those opting for m, that is without considering self-interaction

P {n |n} = n m N a l + b(n l -δ l,m ) + c(n l -δ l,m ) 2 L i=1 {a i + b(n i -δ i,m ) + c(n i -δ i,m ) 2 } , (2)
where the Kronecker term δ l,m is 1 if l = m and 0 otherwise. The term δ l,m is what makes the above expression a conditional probability. With the previous assumptions one has the following:

Proposition 2.1. Assume that a l , b and c are such that for any occupancy vector n the probability of choosing any location is positive, that is p l (n) > 0, ∀l. Then, the Markov chain is irreducible.

Proof. For any system configuration n, let m be a region such that n m > 0 and consider the configuration n = n + δ l -δ m obtained by adding one firm among those who chose l while removing one firm from those who chose m. Then, the one-step probability P {n |n} to move from configuration n to configuration n is given by equation 2. Since by hypothesis P {n |n} > 0, any configuration can be reached from any other configuration in a finite number of steps. The statement follows.

Since the Markov chain associated with the non-linear model is irreducible, it is also ergodic. Namely, there exists a unique distribution π(n; a, b, c) which depends on the set of parameters (a, b, c), such that for any realization {n t } of the process and irrespectively of the initial configuration n 0 it is:

lim T →∞ 1 T T t=0 δ nt,n = π(n, a, b, c) , (3)
where the Kronecker term δ nt,n is 1 if n t = n and 0 otherwise. Equation 3 is precisely the method adopted here to compute the distribution π(n; a, b, c).1 Notably, Proposition 2.1 implies a careful investigation of the region of the parameter space in which the positivity condition of the transition probability is fulfilled. This is indeed the only region in which it is guaranteed that numerical time averages converge asymptotically to the invariant distribution of the process. In this respect, it is important to clarify why the parameter b can be safely assumed to be non-negative in the context of the present analysis. In fact, the model will be used here to test the presence of quadratic externalities in firm localization. Such externalities, however, have already been estimated to be strongly positive under a linear specification. In particular, [START_REF] Bottazzi | Cities and clusters: economy-wide and sector specific effects in corporate location[END_REF] have done so using exactly the same data that will be used below in Section 4. Crucially, their analysis compares two alternative nested models which differ from each other only for the presence of the linear externality term. In their procedure of model selection, the model with linear externalities outperform systematically the model without externalities, even after penalizing parametric numerosity. For this reason, it is safe to assume that b ≥ 0 in the specific context of the present work. Conversely, more attention needs to be placed on the behavior of the model for varying values of c.

To this purpose, Figure 1 shows the behavior of p l as a function of n l for varying values of c. Notwithstanding the constraints imposed on c to guarantee p l ≥ 0, a wide range of behaviors can still be observed. When c = 0 as in Figure 1a, the probability p l grows linearly in n l according to equation (1), thus delineating the standard EhrenfestâĂŞBrillouin model. Conversely, as c moves away from 0, as in Figures 1b-1e, the probability p l becomes a non-linear function of n l . At least two further aspects of Figure 1 are worth pointing out. First, having c < 0 is not a sufficient condition for p l (n l ) to be non-monotone. In fact, there exist negative values of c for which p l is increasing in n l , as shown in Figures 1d-1e. Second, the magnitude and sign of the derivative ∂p l /∂n l will depend on n l . These two aspects imply that it is not trivial to determine whether individual choices are affected by a congestion effect. Defining congestion as ∂p l /∂n l < 0, such condition occurs only for sufficiently negative values of c and for sufficiently high values of n l , as shown in Figure 1f. In this sense, if occurring at all, congestion is not a general condition affecting all regions, but rather a condition that would be specific to some regions according to their actual occupancy.

Before moving on to the inferential analysis, it is worth making a remark on the applicability of the present model to empirical data. A time tick in the model corresponds to a localization choice rather than to a unit of real time; hence, the equilibrium distribution π(n; a, b, c) describes the state of the system after the occupancy vector has been revised a sufficiently large number of times. Practically, this occurs both when incumbents move their establishments and in the case of genuine entry and exit. Considering that the gross yearly turnover rate in any sector is typically between 15% and 20% in OECD countries (see Bartelsman et al., 2005, p.378), the real time needed to reach the invariant distribution may well be relatively short.

Inferential analysis

Under the assumption that the model described above represents the true data-generating process, its estimation consists in selecting the set of parameters (a, b, c) that generate the invariant distribution for the occupancy vector n(a, b, c) which is most statistically "compatible" with the observed occupancy vector n o . Specifically, the parameters are estimated by minimizing the chi-square statistic of the implied invariant distribution of occupancy vectors π(a, b, c) with respect to the observed occupancy. The "categorization" required by the chi-square statistics is obtained by binning the locations according to their occupancy to obtain a binned occupancy distribution. More precisely, this is how the procedure works. 2 Let f (n) denote the number of regions chosen by n firms. This can be either the number observed in the data or the expected number implied by the invariant distribution of the process. Thus, f (0) is the number of regions chosen by zero firms, f (1) is the number of regions selected exactly by one firm, and so on. The sum of these quantities is equal to the number of available regions, that is N n=0 f (n) = L. Given a finite partition of the integers, that is a set of increasing integer numbers C = {y 1 , y 2 , . . . , y J }, let h j with j = 0, 1, . . . , J denote the number of regions 2 The use of maximum likelihood might appear as more appropriate. However, it poses enormous numerical difficulties as the simulation length necessary to assign reliable probabilities to all system configurations soon becomes soon unfeasible with the increase of the number of locations and firms. This problem is efficiently overcome by considering a binned statistics and reverting to chi-square minimization. Dubious readers are reminded that chi-square minimization in general possesses an efficient asymptotic behavior and, in any case, it is not more severely affected by small sample biases than maximum likelihood.

with occupancy in the class (or bin) [y j , y j-1 ), that is:

h j = y j+1 -1 n=y j f (n)
with y 0 = 0 and y J+1 = +∞. The set {h 0 , h 1 , . . . , h J } defines the binned occupancy distribution of the regions with respect to their occupancy. Notice that J j=0 h j = L. Keeping fixed the partitions C, one can compute the observed occupancy distribution h o using the observed distribution of firms across regions and also the theoretical occupancy h π (a, b, c) implied by the invariant distribution of the occupancy vector π(a, b, c). Then the parameters can be estimated solving the problem

(â, b, ĉ) = argmin a,b,c χ 2 (h π (a, b, c), h o ) = J j=0 (h π j (a, b, c) -h 0 j ) 2 h π j (a, b, c) . ( 4
)
To improve the efficiency of the method, we define the partition C in such a way that the different bins contain a roughly constant number of regions when the observed occupancy is considered. This implies that the histogram of the observed binned occupancy distribution is roughly flat (see Figure 2b). When c = 0, the estimation procedure can be carried out more easily because the theoretical occupancy distribution can be directly obtained from the analytic expression of the invariant distribution π [START_REF] Bottazzi | Cities and clusters: economy-wide and sector specific effects in corporate location[END_REF]. In the general case, however, the invariant distribution is not known and one has to revert to numerical methods, specifically Monte Carlo techniques, to obtain an estimate of the theoretical occupancy distribution. 3 Since the main interest here is to test whether or not the null hypothesis H 0 : c = 0 can be statistically rejected, the following exposition we will make explicit reference only to the parameter c while assuming that a and b are known, so as to simplify the notation. Hence, the invariant distribution predicted by the model will be labeled as π(c) and h π (c) is the implied occupancy distribution. In fact, the values of a and b are initialized according to the estimates obtained under c = 0 and are successively re-estimated if the null hypothesis H 0 : c = 0 is rejected (see below).

As usual in inferential analysis, the estimation of the unknown parameter c rests on two steps. First, one searches for the particular value ĉ that generates the occupancy distribution which is closest, on average, to the observed occupancy h o according to the chi-square statistics. The search takes place within a predetermined set [c min , c max ], which contains 0. 4 Second, the null hypothesis H 0 : c = 0 undergoes statistical testing against the alternative H 1 : c = 0. This means one has to evaluate whether the null hypothesis c = 0 can be rejected, at a given level of statistical significance, based on the observation n o . To understand better the estimation procedure, it is convenient to go through its illustration in Figure 3.

The first step consists in deriving the point estimate ĉ as illustrated in the upper-left panel of Figure 3 3 The Monte Carlo techniques adopted here are characterized by very poor performances when a maximum likelihood approach is adopted. The problem arises every time the attractiveness of two locations is similar. In this case any redistribution of firms across locations generates a variation in the likelihood function which is, however, spurious. Washing away this spurious effect requires a computation time which is unfeasible for practical applications. For more details, see [START_REF] Bottazzi | Non-linear externalities: A computational estimation method[END_REF]. 4 To guarantee the convergence of the algorithm and improve its speed, the initial search interval should be tuned on the initial values of the parameters a and b. Essentially, one has to avoid that the search algorithm probes the objective function outside its domain of definition. For more details, see [START_REF] Bottazzi | Non-linear externalities: A computational estimation method[END_REF]. Then, each configuration is further evolved for a sufficiently large number of steps T . Hence, K independent trajectories {n(t min ; c k ), n(t min + 1; c k ), . . . , n(T ; c k )} are obtained for each c ∈ {c 1 , . . . , c K }. Along each trajectory the average occupancy distribution is computed:

h j (c k ) = 1 T t min +T t=t min h t,j (c k ) , (5)
where h t,j (c k ) is the occupancy distribution obtained from the occupancy vector n(t min + 1; c k ).

If t min and T are large enough, this provide a good numerical approximation of the theoretical distribution h π (c k ). Finally, the value of ĉ is chosen solving:

ĉ(n o , t min , T ) = argmin c∈(c 1 ,...,c K ) χ 2 (h(c), h o ) . (6)
By successive refinements of the grid one can achieve the desired precision.5

The second step consists in a statistical test of the null hypothesis H 0 : c = 0 as illustrated in the right panels of Figure 3(b-c). Initially, the model is evolved with c = 0 for s ∈ (1, . . . , S) different realizations, thus obtaining the set of replicas {n 1 , . . . , n S }. Then, the same search procedure described above is performed by considering each replica as the initial observed occupancy. In this way a set of independent estimates {c 1 , . . . , cS } is obtained. From these estimates, we build the empirical distribution function F0 (c). This is an approximation of the theoretical distribution of estimated parameters under the null c = 0, so that the two-sided p-value of ĉ is given by F0 (-|ĉ|) + 1 -F0 (|ĉ|).

If focusing on the estimate of a single parameter might have eased the exposition thus far, the present estimation method can indeed rely on a multistage procedure to estimate multiple parameters (for the details on the coordinate descend method, see Bertsekas, 1999, sect. 1.2.1). Basically, such a procedure consists in estimating each unknown parameter conditional on an initial value of the other parameter(s), thus producing an iterative cycle that stops as soon as

observation n o n(t min ; c 1 ) . . . n(t min ; c k ) . . . n(t min ; c K ) n(t min + T ; c 1 ) . . . n(t min + T ; c k ) . . . n(t min + T ; c K ) ĉ(n o ) = argmin c∈(c 1 ,...,c K ) {χ 2 (h(c), h o )} ĉ model n 1 . . .
c(n 1 ) c(n s ) c(n S ) (c = 0) c 1 c k c K χ 2 χ 2 χ 2
Hypothesis testing

0 ĉ c F (a) (b) (c)
Figure 3: Estimation approach (grid method).

convergence is reached. As an example, imagine that the vector of intrinsic advantages a were known, while the two externality parameters (b, c) were unknown. Starting with initial values (b 0 , c 0 ), a first estimate ĉ1 is obtained keeping the value of b fixed. If ĉ1 is significantly different from c 0 , one estimates a new value b1 keeping c = ĉ1 fixed. If b1 is significantly different from b 0 , one re-estimates c with b = b1 and so on until the procedure convergences to a couple of values ( b, ĉ).

Application

This section shows how the model described and the related estimation method detailed above can be applied to detect the presence of a quadratic component in the localized externalities that drive the spatial distribution of firms. More precisely, the aim is to test whether or not the null hypothesis H 0 : c = 0 is to be rejected in the different sectors. This inferential exercise is carried out at the sectoral level using Italian census data for year 2001 (see ISTAT, 2006). To make the present results directly comparable with those of [START_REF] Bottazzi | Cities and clusters: economy-wide and sector specific effects in corporate location[END_REF], the number of plants by commuting zone is used to measure n l and their same 2-3 digit NACE sectoral disaggregation is adopted. If the null hypothesis H 0 : c = 0 is not rejected, then the model presented in Section 2 reduces to the Polya model with linear externalities estimated by [START_REF] Bottazzi | Cities and clusters: economy-wide and sector specific effects in corporate location[END_REF]. In this case, their estimates will hold as being sufficiently accurate and quadratic externalities would turn out to be redundant to explain the observed spatial distribution of plants. Otherwise, if the null hypothesis H 0 : c = 0 is rejected, the results obtained under linearity would need to be revised.

According to the model in Section 2, the observed distribution of plants across commuting zones is interpreted as the equilibrium outcome of discrete localization choices on the side of firms. Consequently, the effect of non-linear externalities in shaping the observed spatial distribution of plants is detected by estimating the unknown parameter c. The analysis operates at the sectoral level, so that each economic sector taken into account is characterized by a varying number N of plants (reported in the second column of Table 1) located across the L = 686 commuting zones that compose Italy. In this sense, two of the parameters of the model in Section 2, namely N and L, are given directly by the data. Other parameters, instead, need to be tuned.

On the one hand, S is a precision parameter. To recall, S defines the number of stochastic realizations that determine the distribution of c. Hence, higher values of T and S will ensure greater precision, but with an increasing computational cost. For the purpose of the present application, the authors set S = 1000, which ensures a level of precision in the estimate of p-values in the order of 10 -2 .

On the other hand, t min and T are tuned in relation to the size of N . To recall, t min defines the minimum number of time steps needed to guarantee the convergence of numerical simulations to the invariant equilibrium distribution π = (a, b, c), whereas T defines the number of time periods over which equilibrium time averages are computed in problem (6). Therefore, t min can be interpreted as indicating how many times, on average, the N firms have to revise their individual choices in order to reach the equilibrium. As shown by [START_REF] Garibaldi | Finitary probabilistic methods in econophysics[END_REF], the rate of approach to equilibrium as a function of N is

r = A/b N (A + N -1) . (7)
Since the time to reach the equilibrium is of the order of 1/r, one can accordingly set t min = 1/r time steps. In parallel, the number of time periods on which to compute time averages is set to T = 2/r, which guarantees the use of the ergodic property of the process. For estimation to be feasible, it is also necessary to know the vector of intrinsic features a = (a 1 , . . . , a L ). Otherwise, according to equation (1), the model would have L + 2 unknown parameters, thus making estimation unfeasible. To escape this problem, a l can be characterized as a function of H L variables that are meant to describe the regions at stake along some relevant dimensions. As a consequence, the number of unknown parameters reduces to H +2 L. To make the results fully comparable with those of [START_REF] Bottazzi | Cities and clusters: economy-wide and sector specific effects in corporate location[END_REF], their same Cobb-Douglas form is assumed:

a(β, x l ) = exp H h=1 β h log(x h,l ) + β 0 = H h=1 x β h h,l e β 0 . (8)
When substituting the specification for a in (8) into the expression of p l in (1), one obtains the probability for commuting zone l to be chosen at given time step, which serves to evolve the system as described in Section 2 and obtain the simulated trajectories that are necessary for the numerical optimization described in Section 3. In the present application, rather than to estimate fully each parameter of the model, our primary objective is to test whether the addition of a quadratic externality term is statistically meaningful as compared with the linear case. Hence, one looks particularly at c. To this purpose, one introduces a useful normalization: by dividing numerator and denominator by b the expression for the probability p l of selecting

Results

Considering that the focus here is on the detection of non-linearity, one looks at the estimates of c. To provide a visual summary of the data on which estimation is based, Figure 4 illustrates the spatial distributions for two different sectors together with the spatial distribution of a key control variable such as population. The results of the estimation procedure as applied on a comprehensive set of manufacturing and service sectors are reported in Table 1. As a premise, it should be noted that the estimates can be derived for high values of both N and L, thus allowing one to apply the present methodology at virtually any spatial scale or sectoral disaggregation.

The basic result we obtained is a negative one. The estimate of c turns out to be statistically different from zero only in a very small minority of the 43 sectors under analysis. More precisely, only in four sectors is the null hypothesis H 0 : c = 0 rejected at a 90% confidence level, and only in two sectors it is rejected at a 95% confidence level. These numbers are consistent with the overall validity of the null hypothesis, as they correspond to the 10% and 5% of the independent estimates, respectively. Hence, at the scale of commuting zones, the geography of production across the various sectors of the Italian economy in year 2001 is not generally characterized by quadratic localization externalities. In fact, a linear specification of the Polya model is already accurate enough to capture most of the effect of externalities on firm localization. Therefore, adding an extra firm to a location l that hosts n l = 10 or to a location m that hosts n m = 100 firms has approximately the same incremental effect on the probability for each location to further attract other firms at the next localization round. Nonetheless, all other things being equal, positive localization externalities make m more attractive than l in absolute terms, precisely because m hosts more firms to begin with.

Moreover, in the rare instances in which it is significantly different from zero, ĉ is positive. This result rules out the presence of congestion effects at the scale of commuting zones, especially if one takes into account that the condition ∂p l /∂n l < 0 would need a sufficiently negative value of c in order to be met (see the related discussion above in Section 2). In fact, albeit rarely, firm localization might be affected by more-than-linear positive feedbacks, thus making the overall effect of localization externalities at the scale of commuting zones positive.

That is why having an extra firm generally increases the attractiveness of a commuting zone.

Relatedly, the multistage procedure that would serve to estimate the other unknown parameters (β, b) is almost never triggered. Since c is normally not statistically different from zero, the multistage estimates of the other unknown parameters turn out to be never statistically different from the initialization values obtained by estimating the linear Polya model. In this sense, the original marginal effects put froward by [START_REF] Bottazzi | Cities and clusters: economy-wide and sector specific effects in corporate location[END_REF] are not falsified. It also follows that their conclusion about the relative strength of localization and urbanization economies result to be robust to a changing functional specification of localization externalities. Similarly, other results in the literature that were based on a linear specification of externalities may also possibly represent a sufficiently accurate estimate of the determinants of firm localization (see [START_REF] Black | Spatial evolution of population and industry in the United States[END_REF][START_REF] Desmet | Employment concentration across US counties[END_REF][START_REF] Devereux | The geographic distribution of production activity in the UK[END_REF][START_REF] Dumais | Geographic concentration as a dynamic process[END_REF][START_REF] Duranton | Testing for localization using micro-geographic data[END_REF][START_REF] Ellison | The geographic concentration of industry: Does natural advantage explain agglomeration?[END_REF][START_REF] Henderson | Marshall's scale economies[END_REF][START_REF] Maurel | A measure of the geographic concentration in French manufacturing industries[END_REF][START_REF] Rosenthal | The determinants of agglomeration[END_REF].

Conclusion

This paper has presented an empirical analysis on the functional shape of localization externalities. So far, most of the econometric studies regarding the determinants of the localization of economic activities have adopted a linear specification of externalities. Such an approach may produce an inaccurate measurement of the actual strength of localization externalities, par-ticularly if congestion effects are present. Hence, before jumping to strong conclusions about the weight of externalities relative to other determinants of firm localization, one may want to investigate the effect of an alternative function specification of externalities. In particular, the present work has adopted a quadratic specification. The main result is of a negative kind. The quadratic externality coefficient is almost never statistically different from zero. Therefore, localization externalities are approximated accurately enough by a linear specification and, as already discussed, by a specification that suggests the presence of positive technological externalities that increase with the number of co-located firms. If one accepts the assumption that pecuniary agglomerating factors act more or less commonly across the entire economy while technological ones are somehow limited in scope inside a specific sector, the analysis confirms a strong presence of the latter in basically all sectors under scrutiny, reinforcing the conclusions put forward by [START_REF] Bottazzi | Cities and clusters: economy-wide and sector specific effects in corporate location[END_REF].

Notably, the entire analysis is structured so as to allow for a direct comparison with [START_REF] Bottazzi | Cities and clusters: economy-wide and sector specific effects in corporate location[END_REF]. On the one hand, this approach allows one to test their results further, which however cannot be generally falsified. On the other hand, the focus here on comparison has led to the adoption of the same data. This meant estimating the quadratic Polya model at the scale of commuting zones. Nonetheless, non-linearities in firm localization could possibly be occurring at other spatial scales. In particular, even in the very few sectors in which the quadratic externality coefficient is not statistically equal to zero, its size is extremely small. This may signal that commuting zones are too coarse at a spatial scale to make non-linearities detectable at a sufficient level of statistical precision. For this reason, it could be interesting to apply the analytical framework presented here on finer spatial scales, which may allows to capture better, for instance, the effect of spatial congestion.

Figure 1 :

 1 Figure 1: The probability p l as a function of n l for varying values of c. Note: The common parameters of these examples are L = 2, N = 100, a 1 = 1, a 2 = 2, and b = 0.5.

  (a). Starting from the observed configuration n o , the model is simulated with K different values of c ∈ [c min , c max ]. For each value of c, the model is run for an initial number of steps t min to eliminate transient effects (see Section 4 below for further discussion on t min ).

Figure 2 :

 2 Figure 2: Binning procedure

Table 1 :

 1 Number of plants, estimate of c, and corresponding p-value.

	Sectors	N	ĉ	p-value
	15-Food products	73680 6.77E-006	0.460
	17-Textiles	31984 -4.63E-006	0.300
	18-Apparel	46377 -7.94E-006	0.360
	19-Leather products	24195 1.91E-004	0.230
	20-Wood processing	50250 -1.73E-005	0.360
	21-Pulp and paper	5175 6.82E-004	0.190
	22-Publishing and printing	29166 1.91E-005	0.500
	23-Coke,petroleum and nuclear fuel	913 1.45E-003	0.110
	24-Organic and inorganic chemicals	7721 -3.99E-004	0.210
	25-Rubber and plastic products	15115 7.87E-005	0.025
	26-Non metallic mineral products	31177 1.10E-004	0.160
	27-Basic metals	3984 5.22E-004	0.410
	28-Fabricated metal products	102295 -7.67E-006	0.420
	29-Industrial machinery	46481 -5.68E-005	0.570
	30-Office machinery	1715 -7.84E-004	0.540
	31-Electrical machinery	20282 -7.12E-005	0.590
	32-Radio TV and TLC devices	9677 9.27E-005	0.330
	33-Precision instruments	26244 -6.38E-005	0.320
	34-Motor vehicles and trailers	2229 5.41E-004	0.690
	35-Other transport equipment	4951 -7.25E-005	0.730
	361-Furniture	35784 5.47E-005	0.680
	362-Jewelry	10906 2.02E-005	0.790
	363-Musical instruments	695 -4.96E-003	0.081
	36R-Residual of sector 36	6728 5.97E-008	0.710
	40-Electricity and gas	4159 -6.97E-004	0.730
	41-Water	1408 -5.66E-003	0.220
	45-Construction	529757 1.50E-002	0.015
	50-Sale and services of motor vehicles	164079 2.94E-006	0.810
	51-Wholesale and commission trade	404278 3.75E-007	0.840
	52-Retail trade	772730 1.80E-005	0.011
	55-Hotels and restaurants	261304 -2.61E-005	0.320
	60-Land transport	135135 -1.56E-005	0.690
	61-Water transport	1319 -3.88E-004	0.840
	62-Air transport	457 6.22E-004	0.910
	63-Auxiliary transport activities	33765 -1.14E-004	0.088
	64-Post and telecommunications	18056 -1.98E-004	0.540
	65-Financial intermediation	30587 -1.04E-004	0.350
	66-Private insurance and pensions	1771 -4.46E-004	0.860
	67-Auxiliary financial activities	84677 4.73E-006	0.920
	70-Real estate activities	149990 -1.57E-005	0.790
	71-Renting of machinery and equipment	13291 4.41E-004	0.110
	72-Computer and related activities	84100 -3.39E-005	0.210
	74-Business services	216883 -2.32E-005	0.120

Note: For each sector, N is the number of plants and ĉ is the estimate of c, which is obtained with the corresponding p-value. The number of regions is fixed to L = 686 commuting zones.

Differently from the analytical solution of the model, the numerical approach does not require that the Markov chain is reversible and does not exploit the detailed balance condition.

The precision is however limited by the finite computation time. As a robustness check, a minimization algorithm based on successive parabolic interpolations was also employed, obtaining consistent results.

commuting zone l becomes

where (β, c) are the H + 1 unknown parameters to be estimated and b = 1 is set without loss of generality.

Given the focus on non-linearity, one can then initialize the stochastic process using the estimates (â 1 / b, . . . , âL / b) obtained under the linear Polya model to check whether the numerical search converges elsewhere relative from where it started. Hence, the first step of the multistage estimation procedure estimates c conditional on the initialization values (â 1 / b, . . . , âL / b) obtained under the linear Polya and b = 1. Only if the null hypothesis H 0 : c = 0 is rejected will the multistage estimation procedure be triggered.

In general, there are multiple intrinsic features of a commuting zone that may be relevant to firm localization, thus entering among the H regressors that shape a(β, x l ). Factors such as the local size of final goods and labor markets, infrastructural endowments, the extent of local productive variety, as well as the presence of especially conducive local institutions may all play a role in the localization choices of firms. Indeed, all these different factors were taken into account in the original analysis of [START_REF] Bottazzi | Cities and clusters: economy-wide and sector specific effects in corporate location[END_REF], so that the present study will use exactly the same controls in order to guarantee comparability. Namely, this paper considers population, spatial extension, the presence of transportation infrastructures, aggregate consumption and average labor productivity of the location, an index of production variety plus two dummies for the presence of industrial districts and metropolitan areas (see [START_REF] Bottazzi | Cities and clusters: economy-wide and sector specific effects in corporate location[END_REF] for further details). Nonetheless, it is worth recalling that population plays by far the most important role, serving as a proxy for the local size of both final goods and labor markets. This fact is discussed at length by Bottazzi and Gragnolati (2015, esp. pp. 9-11, tab. 2), where the marginal effect of population in shaping the attractiveness of commuting zones is shown to be orders of magnitude stronger than any other intrinsic feature. To facilitate visual inspection of this key variable, Figure 4a illustrates the spread of population among Italian commuting zones in 2001.