Dario Slejko
email: dslejko@ogs.trieste.it

Alessandro Caporali
email: alessandro.caporali@unipd.it

Mark Stirling
email: m.stirling@gns.cri.nz

Salvatore Barba
email: barba@ingr.it

Occurrence probability of moderate to large earthquakes in Italy based on new geophysical methods

Keywords: Seismic hazard, Strong earthquakes, Probability, Italy

We develop new approaches to calculating 30-year probabilities for occurrence of moderate-to-large earthquakes in Italy. Geodetic techniques and finite-element modelling, aimed to reproduce a large amount of neotectonic data using thin-shell finite element, are used to separately calculate the expected seismicity rates inside seismogenic areas (polygons containing mapped faults and/or suspected or modelled faults). Thirty-year earthquake probabilities obtained from the two approaches show similarities in most of Italy: the largest probabilities are found in the southern Apennines, where they reach values between 10% and 20% for earthquakes of M W ≥ 6.0, and lower than 10% for events with an M W ≥ 6.5.

Introduction

Plate tectonics processes produce deformation within the crust of the Earth, and the measurement of such surface deformation is an important boundary condition constraint on the interaction of crustal blocks through time. The distribution in space and time of geodetically derived strain rate will correlate, to some degree, with the geometry and activity of the underlying seismogenic sources, as well as with aseismic processes. In this context, many studies around the world have used geodetic strain rates to estimate earthquake recurrence and probability, providing a valuable supplement or alternative to parameters derived from geologic and seismic catalogue data.

During the period spring 2005 to summer 2007 the Department of the Italian Civil Protection funded several seismological and volcanological projects. One of the seismological projects was entitled "Assessing the seismogenic potential and the probability of strong earthquakes in Italy" (designated S2) and its main goals were: (1) to identify of the seismic sources capable of generating destructive earthquakes (i.e. events with a magnitude larger than 5.5) and (2) to assess the occurrence probability of these events for the sources. The S2 project was organized in four tasks, and Task 4 was dedicated to the actual computation of the occurrence probabilities. One of the goals was to assess the occurrence probability of strong earthquakes using seismological information, and to calibrate the results with geodetic data. Our working hypothesis is that large earthquakes occur along major faults according to the characteristic earthquake model [START_REF] Schwartz | Fault behavior and characteristic earthquakes: examples from the Wasatch and San Andreas fault zones[END_REF]. According to this hypothesis, faults show a tendency to generate earthquakes of similar characteristics (magnitude, slip, rupture length, etc.). While considered overly simplistic by the Italian research community, the characteristic earthquake model is useful as a basis for the broad data-based and methodological comparisons made in this study.

Our efforts are strongly motivated by the need to find ways to augment the incomplete coverage of fault mapping in Italy to date. While the fault data set is considered to be 95% complete in the southern Apennines, it is at best only 70% complete in the rest of the country, and 50% complete or less in the offshore (Valensise, personal communication).

Considerable international literature is available regarding the evaluation of the occurrence probability of strong earthquakes on well-defined faults (e.g. [START_REF] Wgcep | Working Group on California Earthquake Probabilities[END_REF] and the assessment of time-independent (e.g. [START_REF] Frankel | USGS National Seismic Hazard Maps[END_REF] and time-dependent seismic hazard [California: 2007[START_REF] Wgcep | The uniform California earthquake rupture forecast, version 2 (UCERF 2)[END_REF], Italy: [START_REF] Peruzza | Earthquake probabilities and probabilistic shaking in Italy in 50 years since 2003: trials and ideas for the 3rd generation of Italian seismic hazard maps[END_REF]]. However, two main problems thwart these efforts in Italy:

(1) the knowledge of the seismogenic faults in Italy is incomplete in terms of the number of faults and their geometric and seismic characteristics and (2) constraints from geodetic data are problematic because the number of permanent global positioning satellite (GPS) stations in Italy is small, the fact that they have only been in operation a short time interval (about 5 years), and the campaign measurements provide velocities with a large associated uncertainty. Not surprisingly, areal seismogenic sources were defined for areas suspected to contain fault sources as yet unmapped in the Database of the Italian Seismogenic Sources [see [START_REF] Basili | The database of individual seismogenic sources (DISS), version 3: sum-marizing 20 years of research on Italy's earthquake geology[END_REF] for details]. In the S2 project, these seismogenic areas were used to estimate earthquake recurrence parameters from geodetically measured strain rates. Specifically, the geodetic strain rate would be converted into seismic moment rate (M 0 R) and then used to give an upper limit to the seismic potential of the seismogenic sources. Critical to this approach would be the assessment of M 0 R from GPS data, the association of this M 0 R to a geographical area, and the exact definition of a seismogenic source within the geographical area.

Our paper summarizes work undertaken thus far in the framework of the S2 project for the assessment of the occurrence probability in 30 years of moderate-to-large earthquakes within seismogenic areas defined in the Database of the Italian Seismogenic Sources. We develop estimates of M 0 R from: (1) observations from permanent GPS stations (geodetic constraint) and ( 2) from a 3D geophysical model that incorporates stateof-the-art knowledge on faults, and rheology, and is calibrated from GPS observations [geophysical constraint: [START_REF] Barba | Numerical modelling of strain rates in Italy[END_REF]].

Basic ideas

Two types of seismogenic sources are defined in the Database of the Italian Seismogenic Sources (Fig. 1). These are mapped faults and seismogenic areas. The mapped faults are generally well constrained from geological and geophysical data, in that a complete geometric and seismic parameterization (length, dip, slip rate, slip-per-event, etc.) are available, along with an evaluation of the uncertainties associated with the source parameters.

Seismogenic areas do not contain mapped faults, but are assumed to produce earthquakes of magnitude 5.5 or greater based on other geological, geomorphological, and geophysical data. Features such as linear valleys along strike from the mapped faults are assumed to be fault controlled, but the lack of field mapping prevents the definition of fault sources at the present time. For the seismogenic areas, the polygon defining the overall source is given in the Database of the Italian Seismogenic Sources, along with the associated parameters (depth, strike, dip, rake, and slip rate, expected maximum magnitude) and uncertainty bounds.

The research documented here is aimed at defining the seismic potential of the seismogenic areas from the present strain rate in Italy and in consideration of regional seismicity patterns. As not all the existing faults in the seismogenic areas are known, a statistical procedure was designed [START_REF] Stirling | Seismotectonic modelling in northeastern Italy[END_REF]) to fill the empty space of the seismogenic areas with modelled faults of rupture lengths similar to those of the known fault sources.

The working hypotheses of the present study are as follows:

1. the regional geodetic strain is proportional to the seismic potential of the region (i.e. strain is released by earthquakes and aseismic creep); 2. the earthquakes occur on a pre-defined set of faults (a combination of mapped faults and modelled faults: the total number of faults in each seismogenic area is given by the sum of the mapped and modelled faults); 3. faults produce earthquakes according to the characteristic earthquake model (i.e. a tendency to produce a narrow range of earthquakes at or near the maximum size possible from physical constraints such as fault length) and the total regional M 0 R is released as the sum of characteristic earthquakes; 4. the general magnitude-frequency distribution of a region is described by the Gutenberg-Richter behavior, i.e. at the regional scale, the frequencies of the characteristic earthquakes form a Gutenberg-Richter distribution whose b value is in agreement with the regional b value, obtained by the past seismicity.

According to the above hypotheses, we assume that in the long-term (10 4 -10 6 years) the majority of regional seismic release (proportional to the regional geodetic strain) will occur on the mapped or modelled faults, each acting according to the characteristic earthquake model and all together representing a Gutenberg-Richter behavior.

The data

The network of some 160 permanent GPS stations has been in operation for a decade (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) and as part of the European Permanent Network of EUREF and the CERGOP 2 Project of the European Union. Additional local densification stations provide a valuable contribution to the estimate of the average surface strain rate. The strain rate budget for central Europe, determined from GPS observations, is of the order of 20-40 nanostrain per year in a circular area of 150-200-km radius of the eastern Alps (see Fig. 2 where a 100-km ray is shown), corresponding to a velocity range of a few millimeters/year over distances of some hundreds of kilometers [START_REF] Caporali | Geokinematics of Central Europe: new insights from the CERGOP-2/Environment Project, J. of Geodynamics[END_REF]. We compute velocity gradients by least squares co-location, which is a minimum variance algorithm capable of rigorously taking into account the stochastic properties of the input velocities [START_REF] Caporali | Average strain rate in the Italian crust inferred from a permanent GPS network. Part 2: strain rate vs. seismicity and structural geology[END_REF]. For this purpose a covariance function is needed in order to represent the fall off of the correlation coefficient with the lag distance (average distance between the stations). Once the covariance function has been assigned then the velocity and the associated uncertainties can be computed at any point. For deformation analyses it is crucial to know how velocity changes spatially. The horizontal velocity gradient can be split into a symmetric and an anti-symmetric part. The symmetric part represents strain rate, whereas the anti-symmetric part represents a rigid rotation and is, hence, ignorable for deformation studies as the rigid rotation has no associated deformation. The symmetric part can eventually be diagonalized, yielding eigenvectors or principal directions of strain rate. The uncertainty in the components of the strain rate tensor can be quantified from the formal uncertainties of velocities at the actual stations. A final question relates to the method used to compute the strain rates. There exist two schools of thought. One school computes the strain rate on a regular grid, and propagates the uncertainty to account for the loss of accuracy as one moves away from the data points. The other school is more conservative, in the sense that the strain rates are computed only at those points where the estimates are sufficiently well constrained (i.e. where a significant number of stations are close to the site of interest). Hence the strain rate map is patchy, but well constrained where the calculations are made. We adopt this latter approach in our analysis. Specifically, we compute the strain rates at the location of those stations which are surrounded by four or more stations with known velocity within a search radius comparable to the decorrelation distance, that is, the distance at which the average correlation of horizontal velocity pairs drops of 50%.

The horizontal velocities at the GPS stations (see Table 1) are inverted into maximum geodetic strain rate using a search radius variable from 100 to 300 km, according to the number of stations Apennines, and a 300-km radius in the rest of Italy. This approach is based on that applied by [START_REF] Ward | Methods for evaluating earthquake potential and likelihood in and around California[END_REF] in California. The maximum geodetic strain rate was translated into geodetic M 0 R at the same station locations by the application of the [START_REF] Kostrov | Seismic moment and energy of earthquakes, and seismic flow of rock[END_REF] formula, considering the volume to which that strain rate is related. This passage from strain rate to M 0 R is not trivial. In fact, the mean strain rate ε is equal to the sum of the moment (M 0 ) tensors of all earthquakes occurring per unit time t in a unit volume V [V = AH, where A is the area and H is the seismogenic thickness; [START_REF] Kostrov | Seismic moment and energy of earthquakes, and seismic flow of rock[END_REF]], and is given by:

. ε = M 0 2μV t (1)
in which μ is the shear modulus. The crustal volume is obviously a crucial parameter. In our case, the average surface geodetic strain rate is computed considering a circle of 100-to 300-km radius and, consequently, it is assumed to represent the volumetric strain for the seismogenic zone beneath that circle. Equation 1 should define the surface strain rate which can be accommodated in earthquake production in the volume V. A similar definition of the strain volume was already applied by [START_REF] Ward | A multidisciplinary approach to seismic hazard in southern California[END_REF][START_REF] Ward | On the consistency of earthquake moment rates, geological fault data, and space geodetic strain: the United States[END_REF] and [START_REF] Savage | Surface strain accumulation and seismic moment tensor[END_REF], although these authors considered the area to be limited strictly by the GPS stations, while our circles are slightly larger. We adopt an average value of 10 km for the seismogenic thickness because the majority of earthquakes in Italy do not exceed this depth. A shear modulus of 3.5 × 10 10 N/m 2 is also adopted. The resulting estimates of regional M 0 R from GPS geodesy show a considerable range across the different stations within a given domain (see Fig. 2); for example, in the eastern Alps the M 0 R of different stations varies between 0.03 × 10 18 and 2.12 × 10 18 N•m/year. The M 0 R of the central station in each domain was chosen as the M 0 R of the whole domain (see Table 2).

Considering the approximations and uncertainties introduced in the method above, we use an additional data set in our analyses. This second data set comes from finite-element modelling conducted to address the misfits between model predictions and a large amount of neotectonic data for Italy [START_REF] Barba | Numerical modelling of strain rates in Italy[END_REF][START_REF] Barba | Experimental evidence for mantle drag in the Mediterranean[END_REF][START_REF] Barba | Active stress field modelling in the southern Apennines (Italy) using new borehole data[END_REF]. In these studies, model predictions were compared to three independent data sets: geodetic horizontal velocities from temporary and permanent GPS stations [START_REF] Serpelloni | Combination of permanent and non permanent GPS networks for the evaluation of the strain rate field in the Central Mediterranean area[END_REF][START_REF] Serpelloni | Kinematics of the Western Africa-Eurasia plate boundary from focal mechanisms and GPS data[END_REF][START_REF] Caporali | Geophysical characterization of the main seismogenic structures[END_REF]; stress regime data, based on relative stress magnitudes, and the directions of maximum horizontal compressive stress [START_REF] Montone | An improved stress map for Italy and surrounding regions (central Mediterranean)[END_REF]). The model incorporates faults and realistic rheology in a two-layer grid (crust and lithospheric mantle) with laterally varying seismogenic thickness, heat flow, and topography. The horizontal components of the momentum equation [START_REF] Kong | Shells: a thin-plate program for modeling neotectonics of regional or global lithosphere with faults[END_REF] were solved to predict long-term horizontal velocities, anelastic strain rates, verti- The index 1 refers to the GPS observations, while the index 2 refers to the results of the geophysical modelling. M 0 R1 is computed for the reference station while M 0 R2 is given by the sum of the M 0 Rs calculated by modelling for the seismogenic areas belonging to each domain. SAN1 represents the number of seismogenic areas inside the search circle and can be larger than the actual number of seismogenic areas inside the domain (the same seismogenic area can belong to more than one domain if it is located in the overlapping areas of search circles). SAN2 is the number of seismogenic areas inside the domain for which M 0 R was possible to compute by the geophysical modelling (no overlapping areas as the seismogenic areas are associated to the pertinent domain only)

cally integrated stresses, and fault slip rates. Seismogenic areas were simplistically represented as a single fault trace in these prior studies. In our analysis, we compute the strain and slip rates on the model grid and transform them into M 0 R for the different seismogenic areas considering the volume represented by the length, width, and thickness of crust represented by the seismogenic area. A combination of fault slip rate data (where available) and strain rate collectively accounts for known as well as unknown faults (see Table 3).

Figure 3 shows the M 0 Rs obtained for the seismogenic areas. We note that the transformation of strain rate into M 0 R by Eq. 1 is simplistic in the sense that it ignores the possibility that some of the geophysical strain rate may be released aseismically. We consider this issue specifically later in the paper. Table 2 gives details about the M 0 Rs in the domains. The index 1 of the table refers to the GPS observations, while the index 2 refers to the results of the geophysical modelling. In the case of the geodetic constraint, the domain M 0 R(M 0 R1 in Table 2) corresponds to that calculated for the central GPS station and represents the sum of the M 0 R of each seismogenic area in the domain, plus the M 0 R of the distributed seismicity (earthquakes of M W less than 5.5 and, consequently, outside the seismogenic areas), plus that released as aseismic creep. In the case of the geophysical constraint, the domain M 0 R(M 0 R2 in Table 2) is given by the sum of the M 0 Rs of the seismogenic areas calculated by the geophysical modelling. In the case of the geodetic constraint, the number of seismogenic areas inside the search circle (SAN1 in Table 2) can be larger than the actual number of seismogenic areas inside the domain because the same seismogenic area can belong to more than one domain if it is located in the overlapping areas of search circles. As it was not possible to compute the M 0 R by the geophysical modelling for all the seismogenic areas, the number of seismogenic areas inside a domain in the case of the geophysical constraint (SAN2 in Table 2) can be less than the actual number of seismogenic areas inside that domain defined in the Database of the Italian Seismogenic Sources. This explains the differences between the numbers in Table 2 and what shown by Figs. 2 and3.

Figure 2 displays the M 0 Rs calculated for the GPS stations from the geodetic observations (circles): the squares highlight the reference stations, whose M 0 Rs are associated to the domain (M 0 R1 in Table 2) and used as input data in the following computations. The circles in Fig. 3 quantify the M 0 R computed by the geophysical modelling for each seismogenic area: the sum of the M 0 Rs inside each domain gives the value reported in Table 2 (M 0 R2). A direct comparison between the M 0 R estimates obtained with the two different methods is not possible. In fact, according to the geodetic constraints, the value reported in Table 2 overestimates the M 0 R of the domain because some seismogenic areas are counted more than once as they appear in several domains, depending on the overlapping of the search radii. According to the geophysical constraints, the value reported in Table 2 underestimates the M 0 R of the domain because it is not possible to compute the strain rate (and, consequently, the M 0 R) for all the seismogenic areas by the geophysical modelling. The two estimates differ by a factor from 5 to 12, and those calculated with geodetic constraints are higher than those computed by the geophysical modelling (see Table 2). This discrepancy is motivated by the fact that the M 0 R in each domain from GPS observations is given by the sum of the M 0 R released as characteristic earthquakes plus that released as distributed seismicity and as aseismic creep, while the M 0 R from geophysical modelling refers only to the contribution of the characteristic earthquakes.

Comparison of moment rate derived from geodetic and seismicity data

To compare the observed geodetic moment rate (M 0 R1) with the seismic one, domain catalogues have been extracted from the Italian earthquake catalogue (Gruppo di lavoro CPTI 2004). These catalogues collect all the events in a circle centered on the reference GPS stations (Faloria, Pavia, Medicina, and Tito) and with the same radius as that used for the computation of M 0 R1: 100 km for Faloria, 200 km for Pavia and Medicina, and 300 km for Tito. The observed seismic moment rate (M 0 R oss ) in each of the four domains has been computed considering all the earthquakes which have occurred since the beginning of the eighteenth century, because this period can be considered complete for M W 5.5 and over. The [START_REF] Hanks | A moment magnitude scale[END_REF] relation was used to compute seismic moment from M W and the observed seismic moment rate of events with M W 5.5 and over (M 0 R oss5.5 ) was computed as well (Fig. 4).

As the above M 0 R oss does not represent the total seismic moment rate, the Gutenberg-Richter relation was calculated for the four domains. First, the seismicity rates for each magnitude class were computed by the [START_REF] Albarello | Seismic hazard estimates using ill-defined macroseismic data at site[END_REF] approach, where the whole catalogue is divided into time intervals for each of which the probability that it is complete is evaluated and the related rate is weighted accordingly. The Gutenberg-Richter parameters (a and b values) were calculated by the application of the maximum likelihood approach according to the [START_REF] Weichert | Estimation of the earthquake recurrence parameters for unequal observation periods for different magnitudes[END_REF] procedure. From the Gutenberg-Richter distribution, the individual annual rates were derived and, from them and the M W s (from the maximum observed magnitude to 0.1), the calculated seismic moment rate (M 0 R cal ) was obtained. As this M 0 R cal was calculated considering all the M W s, it should represent the actual seismic moment rate released by each domain. In addition, the computed seismic moment rate for events with M W 5.5 and over (M 0 R cal5.5 ) was also computed. It can be seen from Table 4, where all the results are reported, that the M 0 R cal is larger than the M 0 R oss and the difference varies from one domain to the other from less than two times in the eastern Alps to more than three times in the western Alps. It is notable the situation of the southern Apennines, where M 0 R oss is larger than M 0 R cal because of the large number of moderate-to-large earthquakes which occurred there in the last three centuries. Increasing the time period over which the calculation is done the value ofM 0 R oss decreases notably (e.g.: it is 12.0 × 10 17 when calculated considering the last five centuries). Moreover, it is quite interesting to observe which is the contribution in terms of M 0 R cal given by the strong earthquakes (M W 5.5 and over) in each domain. Also in this case the ratio spans over a large interval: strong earthquakes contribute largely in the eastern Alps and in the southern Apennines, where their presence is frequent, while their contribution is limited in the other two domains.

The final comparison refers to the ratio between M 0 R cal and M 0 R1 and quantifies the amount of strain which is supposed to be released seismically. This ratio is quite constant around 20-30% with the exception of the western Alps, where it is about 60%. [START_REF] Bressan | Seismic deformation pattern in the Friuli-Venezia Giulia region (north-eastern Italy) and western Slovenia[END_REF] have found that a significant part of the deformation occur aseismically in the eastern Alps.

In summary, the comparison among the three estimates of M 0 R shows that: (1) M 0 R cal is much larger than M 0 R oss , with the exception of the southern Apennines; (2) M 0 R1 is larger (about four times) than M 0 R cal , with the exception of the western Alps. These values will be introduced in the following computations. More precisely, the part of the M 0 R which is supposed to be released aseismically (obtained from M 0 R cal /M 0 R1 in Table 4) will be subtracted from M 0 R1 in the analysis that refers to the domains. The ratio between the total moment rate from the finite element model and the seismic moment rate for M w ≥ 5.5 is 1.06. So, this correction is applied, subtracting 6% from M 0 R2 in the analysis referring to the seismogenic areas. The M 0 R that refers to the distributed seismicity (events with an 2), M 0 R oss is the value observed in the last three centuries (see Fig. 4) according to the Italian earthquake catalogue (Gruppo di lavoro CPTI 2004), M 0 R oss5.5 is the same as M 0 R oss but referring to earthquakes with an M W 5.5 and over, M 0 R cal is the value calculated from the Gutenberg-Richter relation for all M W s, M 0 R cal5.5 is the same as M 0 R cal but calculated for earthquakes with an M W 5.5 and over M W < 5.5), which is not considered in M 0 R2 (Barba, personal communication), can be calculated from the ratio between the M 0 R released by large earthquakes and that released by all earthquakes (M 0 R cal5.5 /M 0 R cal in Table 4) and it is added to M 0 R2 before the computations. Conversely, M 0 R1 does not need any addiction, as it is comprehensive of the earthquakes of all magnitude. As the geophysical model is satisfactorily constrained only for the Apennines, we will restrict the following elaborations only to the seismogenic areas along the Apennines.

b σb M 0 R1 M 0 R oss M 0 R oss5.5 M 0 R cal M 0 R cal5.5 M 0 R cal5.5 M 0 R cal value (×10 17 ) (×10 17 ) (×10 17 ) ( ×10 17 ) (×10 17 ) / M 0 R cal /M 0 R1

Definition of faults inside the seismogenic areas

Our knowledge of the seismogenic faults in Italy is far from complete. Hence an ability to account for unmapped faults in seismic hazard models is highly advantageous. Seismogenic areas in Italy are often empty (no mapped faults within) or partially filled with mapped faults. A statistical procedure was developed by [START_REF] Stirling | Seismotectonic modelling in northeastern Italy[END_REF] to fill the empty space of the seismogenic areas with modelled faults with parameters equivalent to those of the mapped faults. The procedure for defining the modelled faults is of Monte Carlo type, and therefore takes account of all the uncertainties in parameters from the mapped faults. In the end, a distribution of modelled fault source parameters is defined, enabling probabilities of specific rupture parameters to be obtained from the distribution.

Entering into detail, each of the four domains contains a certain number of seismogenic areas. Each seismogenic area contains a certain number of mapped faults; the empty space is the difference between the length of the seismogenic area and the sum of the lengths of the mapped faults. This empty space is filled by modelled faults, whose dimensions are simplistically assumed to mimic those [mean and standard deviation (σ )] of the mapped faults present in the domain. The average value of the length of the mapped faults in each domain is, then, considered representative also of the length of the modelled faults in that domain, and the length of these modelled faults is allowed to vary in the range given by the average length of the mapped faults ± 1 σ . The procedure for filling each empty space of the seismogenic areas is of the Monte Carlo type, in our case with 1,000 repetitions. The algorithm samples the rupture length repeatedly and randomly between the minimum (mean -1 σ ) and maximum (mean + 1 σ ) values, and each sample is used to calculate a set of realistic earthquake rupture lengths that fill the empty space of the seismogenic area. The rupture lengths are assumed to be positioned endon-end, meaning that no overlapping ruptures are assumed. A random number generator producing uniformly distributed random numbers between 0 and 1 [subroutine Random1 in Press et al. (1992)] was the basis for sampling the range, using the extreme endpoints ± 1-σ values for sampling the average rupture.

In the end, we obtain a set of modelled faults, which mimic the distribution of the mapped faults and fill the empty space of each seismogenic area of each domain.

The number of faults of a specific rupture length in each seismogenic area is given, then, by the sum of the numbers of mapped (N1 i ) and modelled (N2 i ) faults, and the total number of faults of a specific rupture length in each domain (N * i ) is obtained summing up the number of faults in each seismogenic area:

N * i = N SA j=1 N1 i j + N2 i j (2)
where N SA is the number of seismogenic areas in the considered domain.

An example of the approach followed for the definition of the modelled faults and the associated parameters is shown for the eastern Alps domain in Table 5 and Fig. 5. This domain contains eight seismogenic areas (purple boxes in Fig. 5a), each of which can or cannot contain mapped faults (red boxes in Fig. 5a); the number of mapped faults in each seismogenic area is reported in Table 5 according to its characteristic magnitude derived by its rupture length [START_REF] Wells | New empirical relationship among magnitude, rupture length, rupture width, rupture area, and surface displacement[END_REF].

The mean length of the mapped faults is 11.96 km with a σ of 4.65 km. After the Monte Carlo simulations, the resulting number of modelled faults is reported in Table 5. Figure 5b shows the combined result for the modelled and mapped fault sources. The distribution of the modelled faults shows a Gaussian pattern with central tendency around magnitude 6.1. The number of modelled faults is larger than that of the mapped faults around magnitude 6.1 while it is minimal at the edges of the distribution. Some peculiar features are given by the number of mapped faults: no faults referring to magnitude 6.2 and 6.3 have been identified yet in this domain while the number of those of exceeding magnitude 6.4 and larger is notable.

Geodetic constraints on seismicity rates

Our approach is to convert M 0 R for each domain into earthquake magnitude and frequency. The earthquakes are assumed to be described by the characteristic earthquake model on a given fault, and these characteristic earthquakes collectively describe a Gutenberg-Richter distribution at a regional scale (i.e. within each domain, or Italy as a whole). We also calculate distributed earthquake recurrence parameters for each domain from a combination of catalogue seismicity and geodetic observations. Gutenberg-Richter b values of each domain are computed from the seismicity observed in the domain itself, while the a value of each domain can be obtained from the geodetic observations, i.e. from the M 0 R in the domain. In fact,

log N = a -b log M 0 ( 3
)
where N is the number of earthquakes with seismic moment larger than, or equal to, M 0 , where

a = a + 9.1 1.5 b (4) and b = 1 1.5 b (5)
considering the relation between magnitude, M W , and M 0 in N•m [START_REF] Hanks | A moment magnitude scale[END_REF]:

log M 0 = 9.1 + 1.5M W (6)
Knowing the value of the M 0 R in the domain under study, in our case obtained from the geodetic observations, and fixing the maximum value for M 0 , derived from the maximum magnitude for that domain, we can write:

M 0 R = i N i M 0 i (7)
where N i is the unknown annual non-cumulative number of earthquakes with M 0i , and the index i represents all the classes of M 0 in the domain.

From Eq. 3 we have

N i = 10 a -b log M 0 i (8)
where N i is the cumulative number of earthquakes with M 0i and above.

From the cumulative number we can compute the non-cumulative number N i :

N i = N i -N i+1 = 10 a -b log M 0 i -10 a -b log M 0 i+1 = 10 a 1 10 b log M 0 i - 1 10 b log M 0 i+1 (9)
Substituting the value of N i given by Eq. 9 into Eq. 7, we obtain an equation where the only unknown term is the a value. In such a way we can compute the a value in the domain under study from strain instead of from seismicity and, consequently, we can obtain the a value from Eq. 4, and the number of earthquakes, N i , for all classes of seismic M 0i from Eq. 9.

In agreement with the identification of faults (mapped and modelled) inside the seismogenic areas introduced in the previous section and considering the existing equivalence between M 0 and fault rupture length, obtained by the application of the [START_REF] Wells | New empirical relationship among magnitude, rupture length, rupture width, rupture area, and surface displacement[END_REF] and [START_REF] Hanks | A moment magnitude scale[END_REF] scaling laws, in each domain, the number of earthquakes, N i , of a specific class of M 0i is given by the product of the number of faults of a given rupture length [N * i in Eq. 2] and the related rate of occurrence, ν i

N i = N * i • ν i (10)
In such a way, we can compute the rate of occurrence for each class of M 0 . This rate will determine a regional Gutenberg-Richter distribution (i.e.: in the domain), while the Gutenberg-Richter distribution will not be reflected inside the seismogenic areas.

In summary, our procedure defines a suite of mapped and modelled fault sources inside each seismogenic area, each with a defined fault source length and, consequently, with a characteristic M 0 . From the M 0 R obtained from geodetic observations, we constrain the occurrence rate for each class of M 0 . The earthquake recurrence for fault sources within seismogenic areas is defined according to the characteristic earthquake model, and the distributed seismicity parameters for the surrounding domains are modelled according to the Gutenberg-Richter relationship.

Geophysical constraints on seismicity rates

We also use strain rates produced through numerical modelling by [START_REF] Barba | Numerical modelling of strain rates in Italy[END_REF]. Such strain rates are those corresponding to the least misfit between geodetic [START_REF] Caporali | Geophysical characterization of the main seismogenic structures[END_REF][START_REF] Serpelloni | Kinematics of the Western Africa-Eurasia plate boundary from focal mechanisms and GPS data[END_REF]) and stress data [START_REF] Montone | An improved stress map for Italy and surrounding regions (central Mediterranean)[END_REF]) and the numerical predictions produced by trialand-error, using thin-shell finite element code SHELLS [START_REF] Bird | Thin-plate and thin-shell finite-element programs for forward dynamic modeling of plate deformation and faulting[END_REF], in a large number of iterations which mostly differ by boundary conditions, fault friction coefficient [as for example in [START_REF] Barba | Experimental evidence for mantle drag in the Mediterranean[END_REF][START_REF] Barba | Active stress field modelling in the southern Apennines (Italy) using new borehole data[END_REF]].

The technique gives a value of the M 0 R for each seismogenic area and, consequently, the computation refers to each seismogenic area instead of to each domain. In the geodetic constraint approach the characteristic earthquakes of the mapped and modelled faults form a Gutenberg-Richter distribution in the domains but not necessarily inside the seismogenic areas. Differently, in the geophysical constraint approach, the characteristic earthquakes of the mapped and modelled fault sources defined inside each seismogenic area collectively produce a Gutenberg-Richter distribution. The observation of a Gutenberg-Richter distribution for large regions, countries, and the globe is well documented (e.g. WGCEP 1995; [START_REF] Stirling | Probabilistic seismic hazard analysis of New Zealand[END_REF].

Application

As a working hypothesis it is assumed that only earthquakes with M W 5.5 and above occur in the seismogenic areas, and while the M 0 R calculated from geodetic observations represents the regional total value (i.e.: earthquakes in the seismogenic areas plus distributed seismicity plus aseismic creep), the M 0 R obtained from the geophysical modelling refers to the seismicity of the seismogenic areas only. In the geodetic constraint case, the percentage of M 0 R accommodated as an aseismic creep was taken from Table 4 (column 10), while the distributed seismicity is represented by the Gutenberg-Richter relation for events of M W less than 5.5. [Eq. 7 applies to the whole range of M 0 (i.e., M W in the range 0-M max ) although only events related to M W 5.5 and larger are treated in the computation of probabilities]. In the geophysical modelling case, the percentage of moment rate related to the distributed seismicity is derived from Table 4 (column 9) and added to the value calculated from the modelling (that refers only to events with an M W 5.5 and larger). In both cases a regional b value is assumed for the earthquake rate calculations (see Table 4). The Poisson distribution is then used to estimate the probability of exceeding magnitude m in t years, F Mt (m):

F M t (m) = P [M t > m] = 1 -e -νt ( 11
)
where ν is the annual rate of magnitude m and over. This approach is applied to four domains covering Italy and for which the regional M 0 R can be computed. The number of events for each M 0 class was determined according to Eq. 10 and it was then scaled in the seismogenic areas according to Eq. 2. The cumulative rate for each magnitude class is obtained in such a way, and the associated probabilities were computed for exceeding different magnitude classes inside each seismogenic area.

The approach based on the geophysical constraints was applied directly considering the seismogenic areas.

Results

We show our 30-year probability for M W ≥ 6.0 and ≥ 6.5 in Table 6 and Figs. 6 and7, obtained by way of the geodetic constraints and derived from the geophysical modelling. It is interesting to note that there is no direct proportionality between the forecasts for M W 6.0 (Fig. 6) and 6.5 (Fig. 7): this is due to the fact that the seismogenic areas are supposed to be composed of faults which behave according to the characteristic earthquake model. Consequently, the probability distribution of each seismogenic area is strongly peaked at specific magnitude values, e.g.: the characteristic magnitudes of the faults existing there. A clear example is given by the probabilities of the two seismogenic areas along the northern coast of Sicily according to the geodetic constraint: the eastern one shows a higher probability for M W ≥ 6.0 (Fig. 6a), while the western one does for M W ≥ 6.5 (Fig. 7a). This aspect simply shows that events with an M W ≥ 6.5 are more likely to occur along the western sector, while events with an M W between 6 and 6.5 are largely more probable along the eastern one. The 30-year probability for M W ≥ 6.0 based on our geodetic modelling approach (Fig. 6a) is low everywhere and varies in one seismogenic area to another mainly from 1% to 5%. Most of the most 'hazardous' seismogenic areas (in terms of earthquake probabilities) are located in the southern Apennines (black large dots in Fig. 6a, corresponding to a probability greater than 5%). The equivalent probabilities in the case of the results obtained from the geophysical modelling approach (Fig. 6b) are very similar to those with geodetic constraints (Fig. 6a) showing almost the same seismogenic areas. Some disagreements can be noted along the Adriatic coast, where the forecasts from the geophysical modelling are higher than those from the geodetic model in the central sector and lower in the southern one. In the southern sector, the geophysical model has been mainly constrained through [START_REF] Serpelloni | Kinematics of the Western Africa-Eurasia plate boundary from focal mechanisms and GPS data[END_REF] geodetic data, showing a moment rate that is lower than the moment rate relative to the geodetic data used to derive the earthquake rates within this work. We consider the probabilities derived by the current geodetic data more reliable than those derived by the geophysical model in the southeastern side of the peninsula as the current data are greater in number and more accurate. The largest probabilities for M W ≥ 6.0 are found in the southern Apennines, where they reach values between 10% and 20% (Fig. 6b).

The 30-year probability for M W ≥ 6.5 (Fig. 7) is obviously very low again and varies mainly between 0.4% and 3%. The agreement between the estimates according the two approaches is again fairly good, although the estimates under geophysical constraints are slightly lower. Again a disagreement can be noted in a small sector of the southern Apennines. The largest probabilities are found again in the southern Apennines, where they reach values lower than 10%.

Figure 8 shows the differences between the estimates obtained by considering the two different constraints and quantifies how much the results with the geodetic constraints differ from those from geophysical modelling. Considering M W ≥ Fig. 8 Comparison of the probability estimates (P1 indicates the exceedance probability in 30 years computed with geodetic constraints, P2 indicates the exceedance probability in 30 years computed with geophysical constraints): a P2 vs. P1 for M W 6.0, b P2 vs. P1 for M W 6.5, c P1 for M W 6.5 vs. P1 for M W 6.0, d P2 for M W 6.5 vs. P2 for M W 6.0 6.0 (Fig. 8a), it can be seen that the great majority of points remain within the lines indicating the ratios 1:2 and 2:1, that is where the results with one approach are less than double those with the other approach. There are a few among the remaining points where the estimates with geodetic constraints are by far larger (more than four times) than those taken from geophysical modelling. The situation is worse when events with an M W ≥ 6.5 are considered (Fig. 8b). In this case, almost all forecasts with geodetic constraints are larger than those from geophysical modelling and even values 10 times and larger are obtained in a few cases.

The only region where the estimates with the two approaches, referring to M W ≥ 6.0 and M W ≥ 6.5 as well, differ greatly is the promontory in the southern Adriatic sea (Gargano promontory). There, the forecasts from the geophysical modelling are much lower than those with geodetic constraints: this difference is notably larger than that in the rest of Italy.

Some explanations for the differences shown by the two approaches can be suggested. Both methods suffer some limitations: the geodetic one because it was possible to compute the strain rate over wide regions (the four domains, see Appendix A) and some peculiar differences are, consequently, lost. The geophysical modelling is not yet perfectly tuned and a satisfactory agreement with all the boundary conditions is not reached yet (see Appendix B).

The lack of proportionality between the estimates referring to a different magnitude threshold has been already justified, and is more evident in the case of a geodetic constraint (Fig. 8c) than when data from the geophysical modelling have been used (Fig. 8d). The latter, in fact, display quite a nice alignment along the line 10:1, indicating that the probability of M W ≥ 6.0 is about 10 times larger than that of M W ≥ 6.5.

Conclusions

We have developed two new approaches to estimating earthquake probabilities for moderateto-large earthquakes in Italy, using geodetic and geophysical modelling methods. The two ap-proaches show a good agreement in the resulting probabilities, with noticeable differences only in limited areas of Italy. Thirty-year probabilities both based on geodetic and geophysical constraints are less than 5% for M W ≥ 6.0 with the exception of the southern Apennines, where they reach values between 10% and 20% in a very few seismogenic areas. In the same areas, 30year probabilities for M W ≥ 6.5 remain lower than 10%. Future work will be focused on improving the methodologies developed, in an effort to constrain better the derived earthquake probabilities.

SINEX format and with removable constraints are combined into one network solution by an independent combination center. The constraints imposed by the individual analysis centers are removed and new constraints are imposed, so that the final adjustment is properly aligned with the IGS reference frame. The comparison of the individual sub-network solutions with the final network solutions yields a quantitative estimate of the mutual consistency of the processing strategies of the LACs. Typical discrepancies between solutions are at the sub millimeter level in translation, fraction of milliarcsec in rotation and few parts per billion in scale. Furthermore, local analysis centers process regional networks, e.g. at a national level, of permanent GPS stations using the same standards as for the EPN sub-networks. This results in additional SINEX files, which may be combined with the EPN SINEX files at corresponding epochs for network densification. Several software packages are available for such combination work: Bernese's ADDNEQ or ADDNEQ2, CATREF, GIPSY, GLOBK are well-known examples. The SINEX format has the advantage that is a software independent format. Hence SINEX files, if generated by a standardized processing, can be considered a form of metadata of higher level than the RINEX files containing raw phase and pseudorange data.

The rest of this appendix analyses velocities and derived products (velocity field, strain rate) of permanent GPS stations resulting from a combination of the weekly SINEX files concerning the EPN, an Italian network processed by the University of Padua (UPA) and an Austrian network processed by the Astronomical Observatory in Graz (GP_). Both Padova and Graz are EPN Local Analysis Centers.

The SINEX used in the combination analysis are summarized as follows: Beginning GPS week 995, the three normal equations are combined with the program ADDNEQ of the Bernese Software v.4.2 and the appropriate constraints are imposed on those stations with position and velocities listed in the ITRF2000 solution. Because our combination scheme fully considers the variance covariance matrix of the individual network solutions, also the non-ITRF2000 stations have coordinates consistently defined with that system. To ensure that the EPN solution is the backbone, the weight of its SINEX files is larger than for the UPA and GP_ solutions A total of 372 permanent GPS stations are present in the combined network, although only for a fraction of them a reliable estimate of the velocity can be made. The ITRF2000 [START_REF] Altamimi | ITRF2000: a new release of the International Terrestrial Reference Frame for earth science applications[END_REF] constraints in position and velocity of the datum defining stations are available at http:// itrf.ensg.ign.fr/ITRF_solutions/2000/sol.php. The datum defining stations are listed below, and are chosen on the basis of their continuous tracking for several years. Figure 9 shows the velocities of stations with sufficiently reliable time series (2.5 years minimum), which were used in the subsequent strain rate analysis.

A2 Statistical properties of the estimated velocities

The horizontal velocities estimated in the ITRF2000 frame exhibit a dominant NE trend of the order of 2 cm/year. Most of this signal can be accounted for with a rigid rotation about an Eulerian pole and can be filtered out. The residual velocities (Fig. 9) are spatially correlated: the likelihood that two stations have similar velocities decreases with increasing distance between the two stations. This likelihood function is shown in Fig. 10 and forms the basis for computing a velocity field and strain rates out of the observed velocities. According to the analytical model shown in Fig. 10, the characteristic distance d 0 is such that the likelihood of velocities of sites at such distance is reduced to 50% that at zero lag. More details are given in [START_REF] Caporali | Average strain rate in the Italian crust inferred from a permanent GPS network. Part 2: strain rate vs. seismicity and structural geology[END_REF].

Once the velocities and their uncertainties are given at each GPS station, and the correlation function has been specified, the velocities and their uncertainties can be interpolated at other points P by a minimum variance (or 'optimal') algorithm known as least squares collocation:

v n v e P = s C d P,s s C d s,s + W ss -1 • v n v e s s, s = station indices (12)
with

W ss = 1 σ 2 s s 1 σ 2 s δ ss (13)
The mean velocity is removed prior to interpolation and added back to the interpolated value.

Likewise, the velocity derivatives [Eq. 14] yield a velocity gradient, or strain rate tensor:

v n,n v n,e v e,n v e,e P = s ⎡ ⎢ ⎢ ⎣ ∂C ∂n ∂C ∂e ∂C ∂n ∂C ∂e ⎤ ⎥ ⎥ ⎦ P,s × s C d s,s + W ss -1 • v n v e s s, s = station indices (14)
The fact that the individual velocities have different uncertainties is embodied into the algorithm by augmenting the correlation matrix C between the sites with a weight matrix W built with the variances of the individual velocities. This set up ensures that the resulting velocity field will be smooth. The velocity field, when computed at the location of a GPS station, will be closer to the observed velocities, the smaller the associated uncertainty. In general the weight matrix W acts as a smoother or low pass filter. This ensures that local velocity anomalies, e.g. induced by monument instabilities, will not affect the velocity field or the strain rate.

The components of the strain rate matrix [Eq. 14] are expressed in a geographical frame. The eigenvectors and their azimuth can be computed by standard matrix diagonalization. We consider as extensional the eigenvector with positive eigenvalue.

The estimate of the strain rate involves a differentiation and must be done with great care. While the interpolation of the velocity reduces the noise by smoothing (integration), differentiation tends to amplify the noise, especially over short distances. If two stations separated by a distance d have each an uncorrelated uncertainty s in the velocity, then a quick estimate of the uncertainty in the resulting strain rate is 2 1/2 s/d. If s = 0.5 mm/year and d = 100 km, the expected uncertainty is 7 nstrain/year (1 nstrain = 10 -9 ). Shorter distances will proportionally increase the uncertainty to levels above the sought for signal, which is of the order of several tens up to 100 of nstrain/year, typically.

It is crucial to identify the optimal locations, from the point of view of minimal variance, where the strain rate can be reliably computed. A first criterion for optimal choice is that there exist uniformly distributed stations in the neighborhood of the computation point. A second criterion is that the nearest stations should be given higher weight than more distant stations, if their velocities have the same variance. This last criterion is automatically satisfied by the collocation algorithm. To ensure that the strain rate is computed at points such that it is sufficiently well constrained by data from neighboring stations we have adopted the following rule: the strain rate is computed at the geographical location of those permanent stations of known velocity, such that at least one additional station with known velocity exists in each quadrant within a radius equal to the correlation distance d 0.

The adopted rule is exemplified in Fig. 11. We require that within a circle of radius d 0 five velocities are defined: the velocity at the center of the circle and at four additional sites. The least squares collocation represents in a continuous form a finite difference between the center velocity and the peripheral velocities. In this manner we ensure that the resulting strain rate reliably describes the rate of deformation in the circled area. Unfortunately the strain rate computed by this rather conservative approach not always relates to areas of seismological interest. Mapping the strain rate to other locations is mathematically feasible but not recommended, due to the fast increase of the variance. The numerical results for northern, central, and southern Italy obtained the velocities of permanent GPS stations resulting from the combined analysis described above were used as geodetic constraints in the computation of the earthquake occurrence probabilities (Table 2). For southern Italy we have a number of stations insufficient to meet the adopted criterion for strain rate computation. Therefore we have used a set of velocities coming from campaign and permanent stations, available from the literature and processed by analysis centers outside the EUREF network. For this reason we have thought it appropriate not mixing the sets of velocities, although the methodologies of computation are, in a broad sense, consistent with each other. strain rates and stresses. Deformation occurs by frictional sliding or nonlinear dislocation creep. Given a strain rate, the deviatoric stress is evaluated separately for each of three flow laws: frictional faulting, dislocation creep (power law), and Newtonian creep (linear). At each point, the flow law that provides the lowest maximum shear stress is selected. The rheological parameters are different for the crust and mantle lithosphere. The rheological parameters impose the lithospheric rigidity and the coupling between the crust and the lithospheric mantle. The chosen parameters represent an average lithosphere. Frictional faulting stress is evaluated under the assumption of hydrostatic pore pressure. Faults are distinguished from continuum elements because of double "slipping" nodes, and a lower frictional coefficient with respect to the continuum medium.

The method incorporates some 3-D characteristics since volume integrals of density and strength are performed numerically in a lithosphere model with laterally varying crust and mantle-lithosphere layer thicknesses, heat flow, and elevation. Crustal and lithospheric mantle structure, and steady state thermal regime have been derived from literature data.

Fig. 1

 1 Fig. 1 Seismogenic areas (red) and mapped faults (yellow) in the Database of the Italian Seismogenic Sources [modified from Basili et al. (2008)]

Fig. 2

 2 Fig. 2 M 0 R (in N•m/year) computed from GPS observations. Blue small dots for M 0 R < 1.0 × 10 18 ; red medium dots 1.0 × 10 18 ≤ M 0 R < 2.0 × 10 18 ; purple medium circles 2.0 × 10 18 ≤ M 0 R < 3.0 × 10 18 ; black large circles M 0 R ≥ 3.0 × 10 18 . The central GPS station of each domain is marked by a square with size and color according to its M 0 R. The numbered large circles identify the four domains: 1 eastern Alps, 2 western Alps, 3 central Apennines, 4 southern Apennines

Fig. 3

 3 Fig. 3 M 0 R (in N•m/year) computed by geophysical modelling for the seismogenic areas. Blue small dots for M 0 R < 5.0 × 10 15 ; red medium dots 5.0 × 10 15 ≤ M 0 R < 1.0 × 10 16 ; purple medium circles 1.0 × 10 16 ≤ M 0 R < 3.0 × 10 16 ; black large circles M 0 R ≥ 3.0 × 10 16 . The numbered large circles identify the four domains: 1 eastern Alps, 2 western Alps, 3 central Apennines, 4 southern Apennines

Fig. 4

 4 Fig. 4 M 0 release since 1700 in the four domains: a eastern Alps, b western Alps, c northern Apennines, d southern Apennines. Solid line for all earthquakes of the Italian

Fig. 5

 5 Fig. 5 Mapped and modelled faults in the eastern Alps domain: a seismogenic areas (purple areas marked by ITSA) and mapped faults (red areas marked by ITGG),

Fig. 6

 6 Fig. 6 Thirty-year probability for M W ≥ 6.0 according to a Poisson model: a computed by geodetic constraints, b computed by geophysical constraints. Blue small dots for

Fig

  Fig. 9 Horizontal velocities of GPS permanent stations in the Alpine Mediterranean area, after removing a common rigid rotation which approximates the rigid rotation of the Eurasian plate in the ITRF2000 frame. Error ellipses are 2σ. Velocities have been computed only for stations with at least 2.5 years of continuous tracking. Not shown are the velocities of EPN stations falling outside the plotting box

Fig. 11

 11 Fig. 11Example of strain rate calculation. On the left figure five sites are present: one is at the center of a circle of radius equal to the typical correlation length of the ensemble of velocities, and four additional velocities are uniformly distributed about that point within the circle. On the right we have the corresponding strain rate eigenvectors. The structure of the velocity distribution is such that one expects an extension in the NE-SW direction. This deformation regime is then a mean property of the entire circular area

Table 1

 1 Strain rate (ε ) and moment rate (M 0 R1), both with related standard deviations, calculated for the GPS stations

	STA	Lon. E	Lat. N	Radius	ε	σ ε	M 0 R1	σ M 0 R1
				(km)	(ns/year)	(ns/year)	(N•m/year • 10 18 )	( N •m/year • 10 18 )
	ACOM	13.51	46.55	100	1.73	7.50	0.033	0.141
	AFAL	12.17	46.53	100	29.99	6.39	0.565	0.120
	AQUI	13.35	42.37	300	31.25	0.25	5.299	0.042
	ARDE	10.20	46.78	200	32.88	9.10	2.478	0.686
	ASIA	11.53	45.87	200	6.88	1.64	0.518	0.124
	BRAS	11.11	44.12	200	17.05	0.33	1.285	0.025
	BRIX	10.23	45.56	200	35.90	2.16	2.705	0.163
	BZRG	11.34	46.50	100	112.73	4.32	2.124	0.081
	CAME	13.12	43.11	300	36.68	3.06	6.219	0.519
	ELBA	10.21	42.75	300	8.97	2.00	1.521	0.339
	GENO	8.92	44.42	200	22.57	0.18	1.701	0.014
	GRAS	6.92	43.75	300	8.28	0.07	1.404	0.012
	LEC1	9.41	45.86	200	49.63	6.22	3.740	0.469
	MDEA	13.44	45.92	200	29.31	0.87	2.209	0.066
	MEDI	11.65	44.52	200	26.30	0.84	1.982	0.063
	MERA	11.16	46.67	100	66.98	1.05	1.262	0.020
	MPRA	12.99	46.24	200	27.65	3.27	2.084	0.246
	NOVA	8.61	45.45	200	37.21	3.34	2.804	0.252
	PADO	11.90	45.41	200	31.58	4.56	2.380	0.344
	PAVI	9.14	45.20	200	8.24	3.28	0.621	0.247
	PRAT	11.10	43.89	200	10.69	1.13	0.806	0.085
	ROVE	11.04	45.89	200	18.09	0.46	1.363	0.035
	ROVI	11.78	45.09	200	44.33	7.14	3.341	0.538
	SARG	9.51	46.98	200	101.92	2.79	7.681	0.210
	TITO	15.72	40.60	300	45.24	1.91	7.671	0.324
	TORI	7.66	45.06	300	34.56	2.14	5.860	0.363
	TREN	11.12	46.07	200	18.69	0.13	1.408	0.010
	TRIE	13.76	45.71	200	17.44	2.95	1.314	0.222
	UNPG	12.36	43.12	300	31.95	0.51	5.417	0.086
	VLCH	13.85	46.61	100	8.12	9.18	0.153	0.173
	ZIMM	7.47	46.88	300	33.26	6.63	5.640	1.124
	ZOUF	12.97	46.56	200	3.38	1.60	0.255	0.121
	available in the circle. More precisely, five seismic
	domains are considered: eastern Alps, western
	Alps, northern Apennines, central Apennines,
	and southern Apennines (Figs. 2 and 3). A 100-km
	search radius was applied in the eastern Alps, a
	200-km radius in the western Alps and northern

Table 2 M

 2 0 Rs of the domains

	Domain	Reference	Search	SAN1	M 0 R1	SAN2	M 0 R2
		station	ray (km)		(N•m/year)		(N•m/year)
	E Alps	AFAL-Faloria	100	10	0.56 × 10 18	5	1 . 0 1 × 10 17
	W Alps	PAVI-Pavia	200	20	0.62 × 10 18	4	1 . 5 3 × 10 17
	N Apennines	MEDI-Medicina	200	47	1.98 × 10 18	18	3.29 × 10 17
	S Apennines	TITO-Tito	300	34	7.67 × 10 18	26	4.77 × 10 17

Table 3

 3 Strain rate (ε )

	and moment rate (M 0 R2), modelled for the seismogenic areas	Source ITSA002 ITSA003	Region Cent. South. Alps Ripabottoni	Lon. E 9.859 15.029	Lat. N 45.460 41.685	ε (ns/year) 0.90 1.60	M 0 R2 (N•m/year • 10 15 ) 6.554 1.544
		ITSA004	Ascoli Satriano	15.677	41.309	1.40	2.450
		ITSA005	Picerno-Massafra	16.328	40.634	21.40	7.068
		ITSA006	Sciacca-Gela	13.374	37.407	4.00	5.314
		ITSA008	Conero onshore	13.668	43.506	7.80	6.391
		ITSA010	Copparo-Comacchio	12.035	44.743	5.90	7.420
		ITSA012	Portomaggiore	12.001	44.629	5.40	7.486
		ITSA013	Aremogna	14.037	41.822	30.60	16.483
		ITSA014	South. Tyrrhenian	14.072	38.460	6.90	53.857
		ITSA015	Crati Valley	16.285	39.186	45.20	55.998
		ITSA016	Aspromonte	15.516	38.133	118.90	16.699
		ITSA017	Scicli-Catania	14.908	37.161	57.90	79.656
		ITSA019	Crotone-Rossano	17.022	39.287	14.60	11.859
		ITSA020	Southern Marche	13.523	43.229	13.50	46.704
		ITSA021	Marsala-Belice	12.921	37.764	5.30	2.142
		ITSA024	Castelpetroso	14.487	41.468	51.30	60.454
		ITSA025	In. C. Apennines	13.166	42.543	29.10	49.059
		ITSA027	Out. C. Apennines	12.401	43.630	12.50	30.959
		ITSA028	Colfiorito	12.889	43.024	81.60	15.391
		ITSA029	Gela-Catania	14.577	37.324	63.60	28.045
		ITSA031	Conero offshore	13.687	43.578	7.50	8.296
		ITSA032	Pesaro-Senigallia	13.251	43.650	4.00	18.134
		ITSA033	Mt. Pollino South	16.215	39.809	22.40	8.153
		ITSA034	Irpinia	15.464	40.663	37.40	54.261
		ITSA035	Ragusa-Palagonia	14.766	37.151	13.90	18.460
		ITSA037	Mugello	11.237	44.004	23.60	4.128
		ITSA038	Mercure Basin	15.929	40.015	39.70	25.133
		ITSA040	Castelluccio	13.350	42.561	28.40	57.716
		ITSA041	Selci Lama	12.226	43.502	37.60	8.261
		ITSA042	Patti-Eolie	14.974	38.298	30.40	9.523
		ITSA043	Pesaro-Senigallia	13.019	43.955	7.00	10.580
		ITSA051	Mirandola	11.354	44.801	3.10	4.685
		ITSA053	Southern Calabria	16.220	38.616	79.70	37.045
		ITSA054	Porto San Giorgio	13.924	42.889	13.70	8.434
		ITSA055	Bagnara	15.928	38.233	15.70	15.197
		ITSA056	Gubbio Basin	12.611	43.210	59.40	23.835
		ITSA057	Pago Veiano	15.163	41.238	7.30	2.808
		ITSA058	Mattinata	15.830	41.710	1.30	0.969
		ITSA059	Tremiti	14.757	42.180	2.20	5.313
		ITSA060	Montello	12.202	45.864	15.90	4.273
		ITSA061	Cansiglio	12.558	46.117	16.20	8.883
		ITSA062	Maniago-Sequals	12.850	46.202	14.60	25.187
		ITSA063	Andretta-Filano	15.415	40.874	49.10	26.406
		ITSA064	Tramonti-Kobarid	13.378	46.294	11.70	33.937
		ITSA066	Gemona-Tarcento	13.265	46.210	28.20	19.751
		ITSA068	Catanzaro Trough	16.365	38.861	30.22	24.122
		ITSA075	Pietracamela	13.704	42.487	9.50	5.804

Table 4

 4 Annual M 0 Rs (in N•m/year) and b values (with related standard deviation) for the four domains Domain

Acknowledgements

The present study has been devel-

oped

J Seismol (2010) 14:27-51 5.5 0.000 1 0.000 0 0.000 0 0.000 0 0.000 0 0.000 0 0.000 0 0.000 0 0.000 1 5.6 0.000 0 0.000 0 0.000 0 0.000 0 0.000 0 0.000 0 0.000 0 0.000 0 0.000 0 5.7 0.000 0 0.000 0 0.000 0 0.000 0 0.000 0 0.000 0 0.000 1 0.000 0 0.000 1 5.8 0.000 0 0.000 0 0.000 0 0.000 0 0.000 1 0.069 0 0.000 0 0.181 0 0.250 1 5.9 0.000 0 0.000 0 0.000 0 0.000 1 0.306 0 0.267 0 0.000 0 0.263 0 0.836 1 6.0 0.256 0 0.015 0 0.127 0 0.000 0 0.826 0 0.250 0 0.000 1 0.270 0 1.744 1 6.1 0.606 0 0.182 1 0.238 1 0.000 0 0.977 0 0.332 0 0.000 0 0.299 0 2.634 2 6.2 0.138 0 0.302 0 0.283 0 0.000 0 1.047 0 0.083 0 0.000 0 0.122 0 1.975 0 6.3 0.000 0 0.319 0 0.326 0 0.000 0 0.355 0 0.000 0 0.000 0 0.000 0 1.000 0 6.4 0.000 0 0.260 1 0.054 0 0.000 0 0.000 0 0.000 1 0.000 0 0.000 0 0.314 2 6.5 0.000 0 0.062 1 0.000 0 0.000 1 0.000 0 0.000 0 0.000 1 0.000 0 0.062 3 6.6 0.000 2 0.000 0 0.000 0 0.000 0 0.000 0 0.000 0 0.000 0 0.000 0 0.000 2

The number of seismogenic areas refers to Fig. 5a.

Mod. indicates the modelled faults (i.e.: resulting number from the Monte Carlo simulations divided by 1000), M. the mapped ones

Appendix A: Geodesy

A1 Validation of GPS velocity data

The analysis of GPS data normally rests on the IGS standards: these prescribe consistent orbits and Earth Rotation Parameters, and recommend models for Phase Center Variations of the antennas, the elevation cutoff angle and a set of datum defining coordinates and velocities. The final product of the adjustment must be available in the SINEX format and the constraints adopted in the adjustment must be explicitly given, so that further analyses can be made with possibly different constraints. An example of this procedure is given by the weekly maintenance of the European Reference Frame done by the European Permanent Network (http://www.epncb.oma.be): 16 local analysis centers (LACs) process partially overlapping sub-networks of several tens of permanent GPS stations in Europe. The 16 weekly solutions in

Appendix B: Finite-element modelling

The finite-element modelling aimed to reproduce a large amount of neotectonic data for Italy by means of trial-and-error approach, using thin-shell finite element code SHELLS [START_REF] Bird | Thin-plate and thin-shell finite-element programs for forward dynamic modeling of plate deformation and faulting[END_REF], in a large number of iterations which mostly differ by boundary conditions, fault friction coefficient [as for example in [START_REF] Barba | Numerical modelling of strain rates in Italy[END_REF] and [START_REF] Barba | Experimental evidence for mantle drag in the Mediterranean[END_REF][START_REF] Barba | Active stress field modelling in the southern Apennines (Italy) using new borehole data[END_REF]]. In all these studies, model predictions were compared to three independent data sets: geodetic horizontal velocities from tem-porary and permanent GPS stations [START_REF] Serpelloni | Kinematics of the Western Africa-Eurasia plate boundary from focal mechanisms and GPS data[END_REF][START_REF] Caporali | Geophysical characterization of the main seismogenic structures[END_REF]; stress regime data, based on relative stress magnitudes, and the directions of maximum horizontal compressive stress [START_REF] Montone | An improved stress map for Italy and surrounding regions (central Mediterranean)[END_REF]. Residuals between model predictions and the calibration data sets allowed to characterize of the degree of realism of the numerical results. The model results that exhibited the lowest misfits were averaged, thereby accounting for uncertainties in boundary conditions and model parameters. Misfits were standardized to account for the different calibration datasets.

The model incorporates faults and realistic rheology in a two-layer grid (crust and lithospheric mantle) with laterally varying seismogenic thickness, heat flow, and topography. The horizontal components of the momentum equation [START_REF] Kong | Shells: a thin-plate program for modeling neotectonics of regional or global lithosphere with faults[END_REF] were solved to predict long-term horizontal velocities, anelastic strain rates, vertically integrated stresses, and fault slip rates.

In the thin-shell neotectonic modelling program SHELLS [START_REF] Bird | Thin-plate and thin-shell finite-element programs for forward dynamic modeling of plate deformation and faulting[END_REF] and references therein] uses isostasy and vertical integration of lithospheric strength to reduce the threedimensional problem of deforming lithosphere to two dimensions, where the horizontal velocity components do not depend on the depth. The components of the momentum equation are vertically integrated using 1-km steps at each of seven Gauss integration points in each finite element (continuum or fault) and the horizontal components of velocity predicted. The vertical component of the momentum equation comes from the assumption of isostatic equilibrium [START_REF] Bird | New finite element techniques for modeling deformation histories of continents with stratified temperature-dependent rheology[END_REF]). The model is composed of two layers with variable thicknesses, plus the topography. Temperature, strength, and the shear stress tensor are depth dependent. By integrating the shear stress along z the problem is reduced to the plane strain. By computing the corresponding strain solution using only the in-plane terms makes the problem bidimensional. SHELLS solves the momentum equation in two dimensions, and the normal vertical stress, assumed as lithostatic, is then added to the result.

The rheology has the same analytical form at all points of the model. The code neglects all elastic strain accumulation and release and solves for velocities, fault slip rates, and anelastic


Greqam

Clément Bosquet
email: clement.bosquet@univmed.fr

Hervé Boulhol

Joao Santos

Thierry Magnac

Pierre-Philippe Combes

SCALE-DEPENDENCE OF THE NEGATIVE BINOMIAL PSEUDO-MAXIMUM LIKELIHOOD ESTIMATOR

Keywords: C13, C21, F10 pseudo-maximum likelihood methods, negative binomial estimator, Poisson regression, gamma PML

, various studies have used the Poisson Pseudo-Maximum Likelihood to estimate gravity specications of trade ows and noncount data models more generally. Some papers also report results based on the Negative Binomial estimator, which is more general and encompasses the Poisson assumption as a special case. This note shows that the Negative Binomial estimator is inappropriate when applied to a continuous dependent variable which unit choice is arbitrary, because estimates articially depend on that choice.

Introduction

Pseudo-Maximum Likelihood (PML) methods were introduced and then derived for Poisson models by Gourieroux, Monfort and Trognon (1984a,b). Following these seminal works, the Poisson PML (PPML) estimator, which assumes proportionality between the conditional variance and the conditional expectancy of the dependent variable, has often been used for count data models. However, beyond count data, [START_REF] Gourieroux | Pseudo Maximum Likelihood Methods: Theory[END_REF] note that "the pseudomaximum likelihood method with Poisson family may be applied even if the dependent variable is any real number". Santos [START_REF] Silva | The Log of Gravity[END_REF] highlight the advantages of this estimator for gravity equations of bilateral trade ows specied in levels, relative to the common practice of estimating these equations in log-levels by Ordinary Least Squares. Indeed, these authors show that the log-linear specication leads to biased estimates following Jensen's inequality, due to heteroskedasticity in trade levels. 1 Moreover, they provide some evidence that the PPML estimator is more ecient than the nonlinear least squares estimator of the trade specication in level.

As a result, a number of empirical studies of trade ows apply the PPML estimator. As an extension, some researchers consider other PML estimators based on non-Poisson distributions such as gamma according to which the variance is proportional to the square of the conditional mean. The Negative Binomial (NB) PML estimator has also been increasingly used recently in trade as well as mergers and acquisitions studies, including [START_REF] Head | How Remote is the Oshoring Threat?[END_REF], [START_REF] Burger | On the Specication of the Gravity Model of Trade: Zeros, Excess Zeros and Estimation[END_REF], [START_REF] Briant | Product Complexity, Quality of Institutions and the Pro-Trade Eect of Immigrants[END_REF], [START_REF] Westerlund | Estimating the Gravity Model Without Gravity Using Panel Data[END_REF] and [START_REF] Garita | Countries of a Feather Flock Together Mergers and Acquisitions in the Global Economy[END_REF]. The NB distribution assumes that the conditional variance is a linear combination, to be estimated, of the conditional mean and of its square. The NB PML estimator is appealing because it encompasses both PPML and gamma PML as special cases.

This note shows that the NB PML estimator is inappropriate when applied to continuous dependent variables, such as trade or M&A ows, for which the choice of the unit measure is arbitrary. For example, in the case of trade equations, the NB PML estimated parameters depend articially on whether trade ows are measured in thousands of dollars, in billions of dollars or in millions of euros. More precisely, when ows are measured in small units (e.g. thousands of dollars), the NB PML converges towards the gamma PML estimator. In contrast, when ows are measured in large units (e.g. trillions of dollars), the NB PML converges towards the Poisson PML estimator. This scale dependence has been unnoticed so far.

As an interesting case, [START_REF] Garita | Countries of a Feather Flock Together Mergers and Acquisitions in the Global Economy[END_REF] use the NB PML estimator with either the value or the number of mergers and acquisitions as the dependent variable. According to this note, the estimator based on the value will articially depend on the choice of unit, while in principle the estimator based on the number is immune to this problem. However, even for count data, the NB estimator is sensitive to whether the dependent variable is measure in the actual number, in hundreds, in thousands, etc.

Section 2 provides the proof of the scale-dependence of the NB estimator, and section 3 illustrates this proposition with an application based on the trade gravity equation.

Proof

The specication is y i = exp(X i β + u i ) where u i is the residual. The rst-order conditions for PPML, NB PML and gamma PML are, respectively (Gourieroux, Monfort and Trognon ;1984a,b):

PPML: i y i -exp(X i β) X i = 0 (1)
NB PML:

i 1 + α exp(X i β) -1 y i -exp(X i β) X i = 0 (2) gamma PML: i exp(-X i β) y i -exp(X i β) X i = 0 (3)
Whereas the underlying assumption of the PPML and gamma PML is that the conditional variance is proportional to the conditional expectancy and to its square, respectively, the NB PML assumes that V ar(y|X) = E(y|X) + α E 2 (y|X), where α is a constant, generally considered to be positive. Eq. (1-3) conrm that when α → 0, NB PML → PPML, while when α → ∞, NB PML → gamma PML.

This note focuses on the impact of using ỹ = λ y as the dependent variable instead of y where λ is a scalar that can be either very small or very large depending on the unit choice.

The rst-order conditions indicate that both the Poisson and gamma estimator are independent of scale, as only the constant, denoted β 0 , is aected by the linear transformation according to β0 (λ) = β 0 + Log λ, such that exp( β0 (λ)) = λ exp(β 0 ). This implies that exp(X i β(λ)) = λ exp(X i β) ∀i, and the FOC (1) and ( 3) are unaected by scale 2 :

     β0 (λ) = β 0 + Log λ βk (λ) = β k ∀k = 0 ⇒ (4) PPML: i ỹi -exp(X i β) X i = 0 ⇔ λ i y i -exp(X i β) X i = 0 (5) gamma PML: i exp(-X i β) ỹi -exp(X i β) X i = 0 ⇔ i exp(-X i β) y i -exp(X i β) X i = 0 (6)
In contrast, the rst-order condition for NB PML (eq. 2) is sensitive to λ. When ỹ is the dependent variable, that condition is:

i 1 + α exp(X i β) -1 ỹi -exp(X i β) X i = 0 (7)
For β to be unaected (except the intercept) by the transformation, the FOC must be:

i 1 + α λ exp(X i β) -1 λ y i -exp(X i β) X i = 0 (8)
The comparison of ( 2) and ( 8) implies that the condition under which the NB PML estimator

is independent of λ is α(λ) = α/λ, with α ≡ α(λ = 1). 3
However, this condition is violated in general. Let us rst consider the two-step estimator implemented in various econometric softwares. The rst step consists in computing a consistent 2 β ≡ β(λ = 1).

3 Another way to see this is as follows. The NB PML assumption is

V ar(ỹ(λ)|X) = E[ỹ(λ)|X] + α(λ) E 2 [ỹ(λ)|X].
Under the condition that β is independent from λ (except β0 ), this becomes:

λ 2 V ar[y|X] = λ E[y|X] + λ 2 α(λ) E 2 [y|X] ⇔ V ar[y|X] = 1/λ E[y|X] + λ α(λ) E 2 [y|X]
. Independence of the estimator with respect to λ implies that λ α(λ) = α.

estimator, e.g. PPML which is used in most softwares (Stata, SAS), and ŷi denotes the rst-step estimated observations. In a second step, α is estimated by OLS from the following regression:

(y i -ŷi ) 2 -ŷi = α ŷi 2 + i ( 9
)
where i is a residual, which yields: What happens to this estimator when the linear transformation ỹ = λ y is used as the dependent variable? As shown before, the rst-step PPML estimator is unaected by scale, hence ŷi (λ) = λ ŷi . It follows that:

α = i [(y i -ŷi ) 2 -ŷi ] ŷi
λ α(λ) = λ i [( ỹi -ŷi ) 2 -ŷi ] ŷ2 i i ŷ4 i = λ i [λ 2 (y i -ŷi ) 2 -λ ŷi ] λ 2 ŷi 2 λ 4 i ŷi 4 (11) = i [λ (y i -ŷi ) 2 -ŷi ] ŷi 2 i ŷi 4 = α + (λ -1) i (y i -ŷi ) 2 ŷi 2 i ŷi 4
which proves that:

• λ α(λ) = α as soon as λ = 1 ;

• when λ → +∞, λ α(λ) → +∞ and NB PML → gamma PML ;

• when λ → 0, λ α(λ) → α -i (y i -ŷi ) ŷi 2 i ŷi 4 = -i ŷi 3 i ŷi 4 .
When the software constrains the estimated value to be positive (e.g. Stata), λ α(λ) → 0, and NB PML → PPML.

Section 3 provides an empirical example illustrating these results.

This problem is not an artefact of using a two-step estimator. Even the theoretical NB PML estimator exhibits scale dependence. To see this, calculate instead of eq. ( 9) the FOC on α. The full log-likelihood is: 4

LogL = i Log Γ(α -1 + y) Γ(α -1 )Γ(y + 1) -α -1 Log[α] + y i X i β (12) -α -1 + y i Log[1 + α exp(X i β)] -Log[α]
where Γ is the standard Gamma function. 5 While dierentiating this expression with respect to β leads to (2), dierentiating with respect to α entails: 6

α -2 i Log[1 + α exp(X i β)] + y i -exp(X i β) α -1 + exp(X i β) - +∞ k=0 y i (k + α -1 )(k + α -1 + y i ) = 0 (13)
With the linear transformation, the FOC becomes:

α-2 i Log[1 + α exp(X i β)] + ỹi -exp(X i β) α-1 + exp(X i β) - +∞ k=0 ỹi (k + α-1 )(k + α-1 + ỹi ) = 0 (14)
As seen above, the non-scale-dependence of β (except the intercept) is equivalent to exp(X i β(λ)) = λ exp(X i β) ∀i and α(λ) = α/λ. By absurd reasoning, the FOC is then written:

λ 2 α -2 i Log[1 + α exp(X i β)] + λ y i -λ exp(X i β) λ α -1 + λ exp(X i β) - +∞ k=0 λ y i (k + λ α -1 )(k + λ α -1 + λ y i ) = 0 (15) ⇔ i Log[1 + α exp(X i β)] + y i -exp(X i β) α -1 + exp(X i β) = i +∞ k=0 y i ( k λ + α -1 )( k λ + α -1 + y i ) (16)
which is absurd given that α is dened according to (13) by:

⇔ i Log[1 + α exp(X i β)] + y i -exp(X i β) α -1 + exp(X i β) = i +∞ k=0 y i (k + α -1 )(k + α -1 + y i ) (17)
That is, the FOC on α is scale-dependent as the right hand side term

i +∞ k=0 y i ( k λ +α -1 )( k λ +α -1 +y i )
in eq. ( 16) depends on λ. The empirical example developed in the following section illustrates 4 Calculations details are provided in Appendix A.

5 Γ(x) = +∞ 0 t x-1 e -t dt.

6 Calculations details are provided in Appendix B.

that the direct estimation of α(λ) based on eq. ( 14) violates the condition that λ α(λ) is independent of scale.

3 Application to the trade gravity equation

Data

Trade ow data are taken from the IMF Direction Of Trade Statistics (DOTS) database. For the year 2000, there are 21,543 non-zero ows between 196 trading partners. 7 The geographical variables (distance between countries, common border, common language and colonial linkage dummies) are provided by the CEPII database8 , and the FTA data are based on [START_REF] Fontagné | A Re-evaluation of the Impact of Regional Agreements on Trade Patterns[END_REF] who improve those used by [START_REF] Baier | Do Free Trade Agreements Actually Increase Members International Trade?[END_REF].

Specication

The bilateral trade equation is estimated following [START_REF] Anderson | Gravity with Gravitas: A Solution to the Border Puzzle[END_REF] according to:

x ij = exp(β 0 + β 1 Log d ij + β 2 B ij + β 3 L ij + β 4 C ij + β 5 F T A ij + F X i + F M j )u ij (18)
where x ij is the value of export from country i to country j, F X i and F M j are exporting and importing countries xed eects, respectively. B ij , L ij and C ij are the traditionnal control covariates: common border, common ocial language and colonial linkage dummies, respectively.

The u ij are the multiplicative error terms of the nonlinear estimates. Following [START_REF] Anderson | Gravity with Gravitas: A Solution to the Border Puzzle[END_REF], importer and exporter xed eects are used to control for multilateral resistance terms as well as for the income levels of both importers and exporters.

Results

Table 1 shows the estimates of equation ( 18 when the unit is small such as thousands of US$), the NB PML are very close to the gamma PML estimates. This illustrates that the NB PML estimator is inappropriate as the estimates depend arbitrarily on the unit choice of the dependent variable.

Table 2 compares the dispersion parameter α(λ) estimated by various NB PML estimators, and shows that for all of them the condition under which these estimators are not sensitive to scale (i.e. λ α(λ) does not depend on λ) is violated. The NB PML estimators that are compared are: that computed by Stata, that by SAS, the two-step estimator using PPML in the rst step according to eq. (9-11), the two-step estimator using the geometric estimator in the rst step according to Head et al. (2009), and the one-step estimator computing α such that likelihood is maximized (eq. 14, using Newton algorithm and PPML for β). Both Stata and SAS compute an iterated estimator, reestimating eq. (9-11) at each step and starting with PPML, with α being the nal iterated value, and yield similar results. For all these NB PML estimators, λ α depend on λ, converging towards zero when λ becomes small and innity when λ becomes large. That is, all NB PML estimators converge towards PPML and GPML, respectively.

Conclusion

Although it is being increasingly used, the NB PML is not appropriate when the unit choice of the dependent variable is arbitrary.

Rewriting:

B i = y i X i β -α -1 + y i Log[1 + α exp(X i β)] (27)
leads to:

LogL = i A i + B i -α -1 Log[α] + α -1 + y i Log[α] = i A i + B i + y i Log[α] (28)
Consistent with i B i as the objective function with respect to β in Gourieroux, Monfort and Trognon (1984a).

B

First-order condition with respect to α

For derivating the Gamma function, the digamma function, denoted ψ is used:9

ψ(x) = D Log[Γ(x)] = Γ (x) Γ(x) = -γ - 1 x + x +∞ k=1 1 k(k + x) ( 29
)
where D is the dierential operator and γ the Euler-Mascheroni' constant.

Dropping terms without α from (28) leads to the objective function with respect to α:

i Log Γ(α -1 + y i ) Γ(α -1 ) -α -1 + y i Log[1 + α exp(X i β)] + y i Log[α] (30)
Rewriting:

C i = Log Γ(α -1 + y i ) Γ(α -1 ) (31)
and:

D i = -α -1 + y i Log[1 + α exp(X i β)] + y i Log[α] (32)
the FOC would be i C i + D i = 0. Dierentiating C i leads to:

C i = -α -2 -γ- 1 α -1 + y i +(α -1 +y i ) +∞ k=1 1 k(k + α -1 + y i ) --γ- 1 α -1 +α -1 +∞ k=1 1 k(k + α -1 ) (33)
C i = -α -2 1 α -1 - 1 α -1 + y i + (α -1 + y i ) +∞ k=1 1 k(k + α -1 + y i ) -α -1 +∞ k=1 1 k(k + α -1 ) (34)
C i = -α -2 y i α -1 (α -1 + y i ) +y i +∞ k=1 1 k(k + α -1 + y i ) +α -1 +∞ k=1 1 k(k + α -1 + y i ) - 1 k(k + α -1 ) (35)
Rewriting E i , the last term within the parenthesis:

E i = α -1 +∞ k=1 1 k(k + α -1 + y i ) - 1 k(k + α -1 ) = α -1 +∞ k=1 -y i k(k + α -1 + y i )(k + α -1 ) (36) = -y i +∞ k=1 α -1 k(k + α -1 + y i )(k + α -1 ) (37)
leads to:

C i = -α -2 y i α -1 (α -1 + y i ) +y i +∞ k=1 1 k(k + α -1 + y i ) -y i +∞ k=1 α -1 k(k + α -1 + y i )(k + α -1 ) (38)
⇔ C i = -α -2 y i α -1 (α -1 + y i ) + y i +∞ k=1 1 (k + α -1 )(k + α -1 + y i ) (39)
⇔ C i = -α -2 +∞ k=0 y i (k + α -1 )(k + α -1 + y i ) (40)
Dierentiating D i leads to:

D i = α -2 Log[1 + α exp(X i β)] -α -1 + y i exp(X i β) 1 + α exp(X i β) + y i α (41) ⇔ D i = α -2 Log[1 + α exp(X i β)] + y i -exp(X i β) α 1 + α exp(X i β) (42)
⇔ D i = α -2 Log[1 + α exp(X i β)] + y i -exp(X i β) α -1 + exp(X i β) (43)
which yields the FOC with respect to α:

α -2 i Log[1 + α exp(X i β)] + y i -exp(X i β) α -1 + exp(X i β) - +∞ k=0 y i (k + α -1 )(k + α -1 + y i ) = 0 (44)

  the Quasi-generalized PML estimator for NB proposed by Gourieroux et al. (1984b), renamed "two-step NB" by Wooldridge (1999), and for which Head et al. (2009) provide a Stata code.

  ) from the PPML and gamma PML estimators, in the rst and last columns, respectively. The columns in between report the NB PML estimates based on dierent unit values for trade ows. Is is computed either by the Stata nbreg (or glm with nbinomial family) estimator, the SAS proc genmod procedure or the two-step Head et al. (2009)'s code, which yield identical estimated values. When ows are measured with a very large unit (trillions of US$), which means that ow values are very small, NB PML and PPML estimates are visually identical. At the other extreme, when ows take very large values (i.e.

Table 2 :

 2 .1, ** p<0.05, *** p<0.01 ; PML = Pseudo-Maximum Likelihood, PPML = Poisson PML, NB = Negative Binomial, GPML = gamma PML ; USD = United States Dollars, Tr.=Trillions, B.=Billions, M.=Millions, Th.=Thousands. Fixed eects are importer and exporter country xed eects. Estimation of the dispersion parameter α = λ Stata = nbreg procedure, SAS = genmod procedure, HMR = Head et al. (2009) ; USD = United States Dollars, Tr.=Trillions, B.=Billions, M.=Millions, Th.=Thousands.

	Tables
	Table 1: Scale-dependence of the Negative-Binomial Estimator
			equation ; 2000
		PPML		Negative Binomial PML		GPML
		PPML	Tr. USD B. USD	M. USD Th. USD	GPML
		b/se	b/se		b/se	b/se	b/se	b/se
	Distance	-0.606*** -0.606*** -0.758*** -1.259*** -1.232*** -1.231***
		(0.032)	(0.032)		(0.035)	(0.023)	(0.024)	(0.024)
	Contiguity dummy	0.644*** 0.644*** 0.500*** 0.880*** 0.959*** 0.959***
		(0.099)	(0.099)		(0.094)	(0.086)	(0.090)	(0.090)
	Common-language dummy	0.154*	0.154*		0.167**	0.513*** 0.541*** 0.541***
		(0.082)	(0.082)		(0.072)	(0.047)	(0.051)	(0.051)
	Colonial-tie dummy	0.175*	0.175*	0.536*** 1.296*** 1.290*** 1.289***
		(0.103)	(0.103)		(0.125)	(0.089)	(0.100)	(0.100)
	Free-trade agreement dummy 0.479*** 0.479*** 0.322*** 0.173*** 0.249*** 0.250***
		(0.088)	(0.088)		(0.073)	(0.060)	(0.061)	(0.061)
	Fixed eects	Yes	Yes		Yes	Yes	Yes	Yes
	Observations	21543	21543		21543	21543	21543	21543
	Notes: * p<0α
	unit		Tr. USD B. USD M. USD Th. USD
	Stata			0	0.9e-4	1.75	2,510
	SAS			0	0.9e-4	1.75	18,186
	two-step (rst step=PPML, eq. 9-11)	0	0.2e-4	0.02	22
	two-step (rst step=geometric, HMR)	0	3.6e-4	0.88	900
	one-step (eq. 14)		1.4e-8	0.6e-4	2.50	3,734
	Notes:

Because E(Log x) = Log E(x) the expected value of the logarithm of trade ows depends on higher moments, including the variance. Since the variance of the residuals is likely to depend on explanatory variables, estimators using the log specication are biased.

Focusing on non-zero ows is sucient for illustration purposes. Including zero ows or focusing on other years unsurprisingly leads to the same conclusion as the proof in section 2 is general.

http://www.cepii.fr/anglaisgraph/bdd/distances.htm, Centre d'Etudes Prospectives et d'Informations Internationales.

[START_REF] Andrews | Special Functions[END_REF]

Appendix

A Log-likelihood of the Negative Binomial estimator

Negative Binomial density:

Log-likelihood:

with:

Replacing µ i by exp(X i β) leads to:

Rewriting:

yields: