Severine Rosat
email: severine.rosat@oma.be

Toshio Fukushima

Tadahiro Sato

Yoshiaki Tamura

Application of a Non-Linear Damped Harmonic Analysis method to the normal modes of the Earth

Keywords: non-linear harmonic analysis, normal modes, superconducting gravimeters, Slichter triplet

The search for weak oscillations in noisy records requires the development of analysis tools more sophisticated than the Fourier spectrum analysis. We present in this paper a new non-linear damped harmonic analysis (NLDHA) method that enables us to detect any damped harmonic signal hidden in noise and to simultaneously estimate the amplitude, frequency, phase and damping factor of the harmonic signal. We also compare this newly developed NLDHA method with a conventional autoregressive method called ESPRIT. The application of the NLDHA method to the seismic normal modes 0S2 and 0S3 excited after the 2004 Sumatra-Andaman earthquake highlights its ability to precisely recover the frequencies and quality factors of the singlets that constitute these two seismic modes. The application of the NLDHA in the sub-seismic frequency band reveals the presence of some unmodelled waves but no detection of the translational oscillation of the inner core has been claimed.

Introduction

The density structure in the deep Earth is poorly constrained, especially inside the core. Knowledge of the density structure inside the core can be achieved by the observation of the surface gravity effect of the core modes and of the translation of the solid inner core in the fluid outer core, the so-called Slichter triplet [START_REF] Slichter | The fundamental free mode of the Earth's inner core[END_REF]. The period of the Slichter triplet is directly dependent on the density jump at the inner core boundary (ICB) as the main restoring force is the Archimedean force. The PREM [START_REF] Dziewonski | Preliminary Reference Earth Model[END_REF] reference model gives a density jump at the ICB of 600 kg/m 3 while recent observations by [START_REF] Masters | On the resolution of the density within the Earth[END_REF], and [START_REF] Koper | Observations of PkiKP/PcP amplitude ratios and implications for Earth structure at the boundaries of the liquid core[END_REF], reached diverging conclusions, respectively that the density jump at ICB should be greater or smaller than the PREM predicted value. Moreover a convincing detection of this harmonic degree one normal mode could also constrain the viscosity at the ICB (e.g. [START_REF] Smylie | The Inner Core Translational Triplet and the Density Near Earth's Center[END_REF][START_REF] Smylie | Viscous and rotational splitting of the translational oscillations of Earth's solid inner core[END_REF]2000;[START_REF] Rieutord | Slichter modes of the Earth revisited[END_REF]. However, many attempts using spectral methods and stacking processes have not yet led to the reliable detection of the surface gravity signature of the free oscillation of the inner core [START_REF] Jensen | Noise limitations in the core-mode band of superconducting gravimeter data[END_REF][START_REF] Hinderer | A search for the Slichter triplet in superconducting gravimeter data[END_REF][START_REF] Rosat | The search for the Slichter mode: comparison of noise levels of superconducting gravimeters and investigation of a stacking method[END_REF][START_REF] Rosat | Optimal Seismic Source Mechanisms to Excite the Slichter Mode[END_REF]Guo et al. 2006a and2006b) while the unique observation by [START_REF] Smylie | The Inner Core Translational Triplet and the Density Near Earth's Center[END_REF] and [START_REF] Smylie | The product spectra of gravity and barometric pressure in Europe[END_REF] has been confirmed by [START_REF] Courtier | Global superconducting gravimeter observations and the search for the translational modes of the inner core[END_REF]. However, [START_REF] Rieutord | Slichter modes of the Earth revisited[END_REF] and [START_REF] Rogister | Splitting of seismic free oscillations and of the Slichter triplet using the normal mode theory of a rotating, ellipsoidal Earth[END_REF] have demonstrated that this unique claim of detection is irrelevant to the theoretical predictions.

As the problem is to detect weak signals buried in larger amplitude geophysical phenomena and with unknown parameters (frequency, damping factor, and excitation source), we have to use more sophisticated tools than the usual spectral methods. As the period of the Slichter triplet is poorly constrained, the frequency must not be fixed a-priori. The process of excitation of the free oscillation of the inner core is also poorly known. Therefore the damping factor must also be an estimated parameter and not fixed a-priori. That highlights our motivation to develop a non-linear damped harmonic analysis (NLDHA), which is based on the non-linear harmonic analysis (NLHA) developed by [START_REF] Harada | A method of non-linear harmonic analysis and its application to dynamical astronomy[END_REF]. The NLHA method was successfully applied to the estimation of the time ephemeris of the Earth [START_REF] Harada | Harmonic decomposition of time ephemeris TE405[END_REF] as well as to the determination of planetary precession [START_REF] Harada | A new determination of planetary precession[END_REF].

Here we propose a generalization of the method to damping oscillations and a first application to the seismic normal modes of the Earth that induced time-varying gravity signals recorded on the Earth's surface by Superconducting Gravimeters (SGs) of the Global Geodynamics Project (GGP; [START_REF] Crossley | Network of superconducting gravimeters benefits a number of disciplines[END_REF] after the 2004 Sumatra-Andaman earthquake. SGs have proved to be well-suited for the study of the long-period seismic and sub-seismic modes (e.g. [START_REF] Freybourger | Comparative study of superconducting gravimeters and broadband seismometers STS-1/Z in subseismic frequency bands[END_REF][START_REF] Van Camp | Measuring seismic normal modes with the GWR C021 superconducting gravimeter[END_REF][START_REF] Widmer-Schnidrig | What can Superconducting Gravimeters contribute to normal mode seismology?[END_REF][START_REF] Rosat | Performance of superconducting gravimeters from tides to long-period seismology[END_REF]. The precise estimation of the frequencies and damping of the seismic modes of the Earth below 1 mHz, is fundamental as they are very sensitive to the density structure inside the Earth. In particular, their frequency splitting is directly linked to the 1D-density structure [START_REF] Widmer-Schnidrig | What can Superconducting Gravimeters contribute to normal mode seismology?[END_REF]. Moreover their splitting possesses high sensitivity to the 3D-density structure in the Earth's mantle and core.

Therefore these modes are a mean to improve the density profile inside the Earth without any tradeoff with elastic parameters. Finally, an application of the NLDHA to the SG gravity residuals at Strasbourg (France) and Canberra (Australia) in the sub-seismic frequency range is performed in the search for the surface gravity effect of the Slichter modes.

Description of NLDHA method

We extended Harada's NLHA [START_REF] Harada | A method of non-linear harmonic analysis and its application to dynamical astronomy[END_REF][START_REF] Harada | A method of non-linear harmonic analysis and its application to dynamical astronomy[END_REF][START_REF] Harada | Harmonic decomposition of time ephemeris TE405[END_REF]2004) so as to be applicable to time series containing unknown numbers of damping oscillations with unknown frequencies and damping factors. The extended method tries to minimize Φ, the weighted sum of the square of residuals of the given data after subtracting a polynomial trend and multiple damping oscillations, defined by:

  2 1 0 1 ) cos sin )( exp( 2 1 ) , , , , (                   M k i i k k i k k i k L j j i j N i i f t c t s t t p w c s p      ,
where N is the number of data, fi is the i-th datum obtained at the time ti, and wi is its weight; L is the order of the trend polynomial, pj is its j-th coefficient; M is the number of damping oscillations, sk and ck are the coefficients of the sine and cosine parts of the k-th oscillation component, the angular frequency and the damping factor of which are ωk and λk, respectively. Note that the introduction of the data weights enables NLDHA to deal with records containing some gaps by attributing a weight of zero to the missing samples.

The outline of the method is roughly the same as that of Harada's method [START_REF] Harada | A method of non-linear harmonic analysis and its application to dynamical astronomy[END_REF].

Once the non-linear parameters, ωk and λk, are all fixed, the other linear parameters, pj, sk , and ck, are optimally determined by the method of weighted linear least squares. In this sense, the optimal linear parameters, pj, sk , and ck, can be regarded as functions of the non-linear parameters. Then we regard the objective function Φ as a function of the non-linear parameters only:

            , ), , ( ), , ( ), , ( ) , ( c s p    .
Thus our problem reduces to a non-linear least squares problem, which is effectively solved by a quasi-Newton method such as the well-known BFGS (Broyden-Fletcher-Goldfarb-Shannon) method [START_REF] Press | Numerical Recipes in FORTRAN: the Art of Scientific Computing[END_REF]. Therefore the issue is reduced to the questions of how to determine the optimal number of damping oscillations, M, and how to find suitable initial guess of the non-linear parameters.

We solve the first issue by the step-by-step approach, namely increasing M one by one starting from 0 as was done in Harada's NLHA [START_REF] Harada | A method of non-linear harmonic analysis and its application to dynamical astronomy[END_REF]. This iterative approach is also similar to the iterative Fourier transform method used by [START_REF] Guo | A spectral search for the inner core wobble in Earth's polar motion[END_REF] in the search for the inner core wobble in Earth's polar motion data. As for the second issue, we pick up a candidate pair of frequency and damping factor which maximizes the decrement of the objective function in the two-dimensional search of parameters in the range pre-assigned. Note that we simplified the procedure of the search by conducting not the full two-dimensional but the alternative search in the sense of frequency-first. The initial frequency is obtained from the power spectrum of the signal for the zero-damping factor. When the frequency has been finely estimated, the damping factor is then searched. The initial damping value is obtained from the maximum of the damping power spectrum estimated at the frequency determined in the previous fine search of frequency. The power spectrum in the damping factor domain is uni-modal, namely containing only one maximum (cf.

Fig. 1). For each direction, in the frequency or in the damping factor, we followed a three-stage approach (1) a coarse but comprehensive search, (2) a refined but still comprehensive search, both based on the quasi-Newton method implemented with the BFGS algorithm, and (3) a parabola fitting near the optimum. The final stage to find the global minimum of Φ is conducted by the Brent's method using the golden rule section [START_REF] Brent | Algorithms for Minimization without Derivatives[END_REF]. The full algorithm is sketched in Fig. 2.

For the optimization subroutines used in NLDHA, the reader may refer to the Numerical Recipes [START_REF] Press | Numerical Recipes in FORTRAN: the Art of Scientific Computing[END_REF] Chapters 10.1 and 10.7, respectively for the golden section search of a local minimum and for the BFGS optimization of the quasi-Newton method.

The NLDHA method is different from the well-known auto-regressive Prony's method, discovered by Prony in 1795 (Prony, 1795) during the long-time work to compute numerical tables of sine and cosine, and its modifications (Osborne 1975[START_REF] Smyth | Coupled and Separable Iterations in Nonlinear Estimation[END_REF][START_REF] Osborne | A modified Prony algorithm for exponential function fitting[END_REF] in the sense that (1) the time of sampling can be arbitrary, (2) the secular trend can be determined simultaneously, (3) the effect of data weight is considered, and (4) the optimal number of damping oscillations is also estimated. Note that Prony's method is applicable only to evenly-sampled data.

Prony's method and its modern variation were re-discovered independently many times, for example, the "Sompi" method [START_REF] Hori | A new method of spectral analysis and its application to the Earth's free oscillations: the "Sompi" method[END_REF][START_REF] Kumazawa | A theory of spectral analysis based on the characteristic property of a linear dynamic system[END_REF][START_REF] Imanishi | Observation of seismic core modes from a superconducting gravimeter record[END_REF],

which is now popular in geophysics. Recently the ESPRIT method [START_REF] Roy | ESPRIT-Estimation of signal parameters via rotational imvariamce techniques[END_REF][START_REF] Roy | ESPRIT-Estimation of Signal Parameters via Rotational Invariance Techniques[END_REF], which is a variation of Prony's method and originally established in the field of acoustic signal processing, seems to be better than the modified Prony's method in analyzing transient phenomena in electric power systems [START_REF] Dafis | Proceedings of the 14 th PSCC held at Sevilla on 24-28 th[END_REF]. The ESPRIT method requires fixing apriori the model order, i.e. the number of harmonic signals to be analyzed. However, based on the estimation error [START_REF] Badeau | Selecting the modeling order for the ESPRIT high resolution method: an alternative approach[END_REF]) for instance, it is possible to evaluate a-posteriori the optimal model order. The characteristics of the NLDHA and the auto-regressive ESPRIT methods are compared in Table 1.

A first application of the NLDHA method to the seismic normal modes excited after the 2004 Sumatra-Andaman earthquake is proposed in the next section. The seismic normal modes are well-known decaying waves; therefore they are useful to test our NLDHA method and to compare it with more conventional analyses.

3. Application to the seismic modes 0 S 2 and 0 S 3 after the 2004 Sumatra-Andaman earthquake

SGs have proved to be well suited for the study of the seismic modes below 1 mHz and the seismic mode analysis after the Mw = 9.3 Sumatra-Andaman earthquake on 2004 December 26 has confirmed their high quality [START_REF] Rosat | High resolution analysis of the gravest seismic normal modes after the 2004 Mw=9 Sumatra earthquake using superconducting gravimeter data[END_REF]. Contrary to the Slichter triplet, the frequencies and Q-values of the seismic modes are well-known from many previous studies based on seismometer data (e.g. [START_REF] He | Normal-mode constraints on the structure of the mantle and core[END_REF] or SG data (e.g. [START_REF] Roult | New determinations of Q quality factors and eigenfrequencies for the whole set of singlets of the Earth's normal modes 0S0,0S2,0S3 and 2S1 using SG data from the GGP network[END_REF]. Therefore the analysis of the seismic modes presents a great opportunity to test the NLDHA method.

In the following, we use 22-day SG data recorded at the Strasbourg (ST, France) station after the 2004 Sumatra event. The gravity records have been corrected for the local tides and atmospheric pressure effect through a nominal barometric admittance of -3 nm/s²/hPa. The atmospheric pressure reduction is necessary to reduce the noise level at these frequencies (e.g. [START_REF] Zürn | On noise reduction in vertical seismic records below 2 mHz using local barometric pressure[END_REF].

In order to reduce the computation time of NLDHA, we limit the analysis to the frequency band 0.28-0.34 mHz by pre-processing the data with a band-pass filter. This frequency range A comparison of the frequency estimates after the 2004 December 26th Andaman-Sumatra earthquake using a standard Lorentzian fitting of an amplitude spectrum [START_REF] Rosat | High resolution analysis of the gravest seismic normal modes after the 2004 Mw=9 Sumatra earthquake using superconducting gravimeter data[END_REF] with the results of the NLDHA and ESPRIT methods shows a perfect similitude within the error bars (cf.

Table 2). The errors computed by [START_REF] Rosat | High resolution analysis of the gravest seismic normal modes after the 2004 Mw=9 Sumatra earthquake using superconducting gravimeter data[END_REF] are based on the signal to noise ratio of each singlet as proposed by [START_REF] Dahlen | The effect of data windows on the estimation of free oscillations parameters[END_REF], and thus do not reflect the method uncertainty. So we have also estimated the uncertainties of the Lorentzian fitting. The estimated quality factors (Table 3) are also in good agreement between the amplitude decay measurement and NLDHA methods. It is worth to note that the NLDHA method gives lower uncertainties on the estimated parameters of 0S2 seismic mode. The predicted frequencies and Q values of 0S2 for a PREM model are also given in Table 3.

We have also applied the NLDHA method to the same SG data but band-pass filtered between the corner frequencies 0.2 and 0.5 mHz. The wider band enables to estimate precisely the 5 singlets of 0S2 and also to recover the 6 singlets of 0S3 (Fig. 4 (a)). Indeed there should be 7 singlets for this harmonic degree 3 seismic mode, but because of the geographical position of ST site with respect to the earthquake source, the central singlet corresponding to m = 0 at the PREM frequency 0.4686 mHz, cannot be observed at ST. In that case, the RMS residuals do not indicate any drop up to the final step (Fig. 4 (b)). The process stops at the step 17 by itself as it cannot find any more harmonic signal to analyze. Although the NLDHA succeeds to find the harmonic signal corresponding to m = 2 and to estimate the frequency and damping factor, it fails computing the uncertainties of these non-linear parameters. The reason is that the frequency resolution is too poor to distinguish the singlet m = 2 from the singlet m = 1, so the Cholesky decomposition used to compute the inverse matrix of the uncertainties fails (the matrix becomes singular). We have resolved that problem by using a Singular Value Decomposition instead of the Cholesky decomposition when the latter fails. The results of NLDHA obtained for 0S3 are compared with the PREM predictions by [START_REF] Millot-Langet | Normal modes and long period seismograms in a 3D anelastic elliptical rotating Earth[END_REF] in Table 2 for the frequencies and in Table 3 for the quality factors.

We have demonstrated with this example that NLDHA can analyze a very small oscillation hardly distinguishable from a larger harmonic signal and give a good estimate of its parameters. Of course to increase the quality of the parameter estimates we could improve the frequency resolution of the 0S3 seismic mode by taking longer data series, when they are available, in a similar way than when doing a spectral analysis.

To conclude this part, the NLDHA method performs well in the analysis of the seismic normal modes of the Earth. In the following section, we will show its application to the search for the Slichter triplet in the sub-seismic frequency band.

Application to the search for the Slichter modes in the sub-seismic frequency band

The mechanism of excitation of the Slichter modes is not well known. It could be excited by a strong earthquake (e.g. [START_REF] Rosat | Optimal Seismic Source Mechanisms to Excite the Slichter Mode[END_REF], randomly excited by some turbulent flow in the fluid outer core and/or by the variable motions of the atmosphere and oceans. We suppose that the Slichter modes are continuously excited or are slowly damping, so that they are statistically well represented in a year of SG gravity records, in order to be detectable by the NLDHA. We consider the SG gravity records at the Canberra (Australia) and Strasbourg (France) SG sites. These SG sites have been chosen because of their low noise level in the frequency range considered [START_REF] Rosat | Performance of superconducting gravimeters from tides to long-period seismology[END_REF] and because of their nearly opposite geographical location. Before applying the NLDHA method,

we have corrected the gravity data for the local tides and the atmospheric pressure effect. The obtained residuals have then been band-pass filtered between the corner frequencies 0.02 mHz and 0.08 mHz. That frequency range covers all the possible periods for the Slichter modes predicted or observed by the previous works (e.g. [START_REF] Smylie | The Inner Core Translational Triplet and the Density Near Earth's Center[END_REF][START_REF] Courtier | Global superconducting gravimeter observations and the search for the translational modes of the inner core[END_REF][START_REF] Rieutord | Slichter modes of the Earth revisited[END_REF][START_REF] Rogister | Splitting of seismic free oscillations and of the Slichter triplet using the normal mode theory of a rotating, ellipsoidal Earth[END_REF].

The NLDHA analysis detects the tidal waves as well as the diurnal harmonics of the atmospheric mode S1 remaining in the SG data. The reconstructed signal and the final residuals are plotted in Fig. 5 The amplitude spectra of the residuals at every five steps from 1 to 50 are also plotted in Fig. 5 (c).

Notice a spectral peak in the Canberra data found at the frequency 0.0736 mHz with a negative Qfactor equal to -37051 and an amplitude of 0.008 nm/s², and a spectral peak in the Strasbourg SG data found at the frequency 0.0774 mHz with Q = 20810 and an amplitude of 0.003 nm/s² with a damping-rate of about 9 years. Note that, as discussed by [START_REF] Zürn | Investigation of the "Nearly diurnal free wobble" -resonance in individual tidal records[END_REF], it is possible to obtain a negative Q-value when the system has a very large Q and the signal amplitude is small.

Another spectral peak with amplitude of 0.01 nm/s² emerges clearly above the noise level in the Canberra residuals at the frequency 0.0655 mHz with an estimated Q-factor of 119860 giving a damping rate of about 58 years. This could be also an actual geophysical signal. However we have not enough evidence for the moment to conclude.

The frequencies detected by [START_REF] Courtier | Global superconducting gravimeter observations and the search for the translational modes of the inner core[END_REF] and attributed to the Slichter triplet were 0.069186, 0.073767 and 0.077544 mHz, with quality factors between 100 and 400 [START_REF] Smylie | The Inner Core Translational Triplet and the Density Near Earth's Center[END_REF]. [START_REF] Crossley | On the excitation, detection and damping of core modes[END_REF] have shown that the damping due to the seismic anelasticity of the inner core and mantle has a Q-value of the order of 5000 with a corresponding damping time of 400 days. The damping due to the outer core viscosity has been formulated by [START_REF] Smylie | Viscous and rotational splitting of the translational oscillations of Earth's solid inner core[END_REF] and also by [START_REF] Rieutord | Slichter modes of the Earth revisited[END_REF]. The estimates of the dynamic viscosity range from 1.6 10 -2 Pa.s using laboratory experiments [START_REF] Rutter | Towards evaluating the viscosity of the Earth's outer core: An experimental high pressure study of liquid Fe-S (8.5 wt.% S)[END_REF], giving a Q-value of the order of 10 7 , to 1.2 10 11 Pa.s [START_REF] Smylie | The inner core as a dynamic viscometer[END_REF] using the claimed Slichter modes of [START_REF] Courtier | Global superconducting gravimeter observations and the search for the translational modes of the inner core[END_REF] giving a Qvalue less than 10. [START_REF] Mathews | Visco-electromagnetic coupling in precession-nutation theory[END_REF] have proposed an upper limit of 1.7 10 5 Pa.s using nutation data corresponding to a Q-value of 5000. The magnetic damping of the inner core has been studied by Buffet and Goertz (1995) who have shown that the Q-value should be between 5.8 10 5 and 2200 for a magnetic field ranging from 0.0005 to 0.001 T.

Since the expected triplet for the Slichter modes has not been commonly detected at the two SG sites Strasbourg and Canberra, we cannot say that the Slichter triplet has been clearly detected from the present analysis. However, the obtained signals close to the proposed frequencies of the Slichter triplet show a very large Q-value (between 20000 and 120000) suggesting a very low viscosity of the fluid core -if these signals can be attributed to the Slichter modes and if we consider that the main damping is due to the viscosity of the fluid outer core -and a very small amplitude of these waves as suggested from the theoretical studies.

Besides the S1 wave and its over-tones observed in the residual time series data indicate a necessity to improve the accuracy of the atmospheric pressure correction. To strengthen the reliability of the detection of the core modes and the Slichter modes, further study is needed;

particularly using an extension the NLDHA to combine several datasets should be considered (e.g.

expanding the present NLDHA algorithm so that it can treat many SG data sets simultaneously).

Conclusion

The advantages of the NLDHA method reside in the complete analysis of the signal without fixing a-priori the number of waves, a highly precise estimation of the parameters values and the application to unevenly spaced datasets.

We have demonstrated that the NLDHA is well suited to the precise analysis of the Earth's normal modes. The NLDHA will be also very useful to detect and analyze the time-varying gravity record in the search for the weak oscillations of the core and inner core. In particular, a multistation extension of the NLDHA would be required to stack several SG records and search for the global signature of the translational motion of the inner core. This research will be the aim of a future work. Frequency (mHz) m = -2 m = -1 m = 0 m = 1 m = 2 Lorentzian fit (Rosat et al.2005) 0.29997 0.30462 0.30929 0.31389 0.31845 SNR errors (Rosat et al.2005) 6.6 10 -5 1.1 10 -4 6.1 10 -5 9.6 10 -5 3.3 10 -5 Uncertainties 1.6 10 -5 1.9 10 -5 1.5 10 -5 2.1 10

Figure caption

  contains the spheroidal seismic mode 0S2 only. The second subscript corresponds to the spherical harmonic degree l=2 and the first figure is the radial number n=0 in the notation nSl for the spheroidal modes. Because of the Earth's rotation and ellipticity, each mode is split into 2l+1 singlets having distinct frequencies. Thus 0S2 is split into 5 frequencies. The amplitude spectrum of 0S2 is represented in the upper plot of Fig.3 (a). The reconstructed signals and the residuals obtained with NLDHA and with ESPRIT are plotted in Fig.3 (a). We have also plotted the RMS residuals after each step of the NLDHA in Fig.3 (b). The sudden and clear drop after the 5 th step, i.e. after the analysis of five decaying waves, indicates that all the meaningful harmonic signals have been detected and that only noise is remaining.

  (a). The diurnal harmonics of S1 and its overtones (from S2 to S7) are clearly observed. Moreover, the main frequencies of the non-linear tides (MN4, M4, SN4, MS4 and MK4 for the quart diurnal band; 2MN6, M6, MSN6, 2MS6, 2MK6, 2SM6 and MSK6 for the tidal band around 6h) are indicated. The RMS amplitudes of the residuals at each step are plotted in Fig.5 (b).

Fig. 1

 1 Fig. 1 Example of a damping power spectrum. Fig. 2 NLDHA algorithm.

Fig. 3

 3 Fig. 3 NLDHA and ESPRIT method applied to the time-varying gravity record of the SG at Strasbourg after the Sumatra event. (a) From top to bottom, amplitude spectrum of the SG data, amplitude spectrum of the reconstructed signal and amplitude spectrum of the residuals (original signal minus the reconstructed damped oscillations) obtained from the NLDHA (solid line) and ESPRIT (dotted line) analysis results. The data have been band-pass filtered between the corner frequencies 0.28 and 0.34 mHz before applying the NLDHA and ESPRIT. (b) Residual amplitude (RMS) after each step of the NLDHA. Note that the step number corresponds to the number of selected waves.Fig.4 Results of the NLDHA method applied to the time-varying gravity record of the SG at Strasbourg after the Sumatra event. (a) From top to bottom, the amplitude spectra of the SG data, of the reconstructed signal and of the residuals (original signal minus the reconstructed damped oscillations) obtained from the NLDHA analysis results are plotted. The data have been band-pass filtered between the corner frequencies 0.2 and 0.5 mHz before applying the NLDHA. (b) Residual amplitude (RMS) after each step of the NLDHA. Note that the step number corresponds to the number of selected waves.Fig. 5 Results of the NLDHA method applied to the one-year time-varying gravity record of the Canberra (Australia) and Strasbourg (France) SG sites, respectively on the left-hand and right-hand sides. (a) From top to bottom, the amplitude spectra of the SG signal, of the reconstructed signal and of the residuals (original signal minus the reconstructed damped oscillations) obtained from the NLDHA analysis results are plotted. (b) Residual amplitude (RMS) after each step of the NLDHA. (c) Amplitude spectrum of the residuals after every five steps from 1 to 50. The diurnal harmonics of the thermal atmospheric wave S1 from S2 to S7 are indicated in vertical dashed and dotted lines and the non-linear tides around 4 and 6h are indicated in vertical dotted lines.

Fig. 1

 1 Fig.1

  Fig. 2

Table 1

 1 Characteristics of the NLDHA and ESPRIT methods.

	NLDHA	ESPRIT

Table 2

 2 Estimated frequencies in mHz of 0S2 singlets from the 22-day Strasbourg SG record after the 2004 Sumatra earthquake using a Lorentzian fitting in amplitude spectrum, NLDHA and ESPRIT methods. Note that in[START_REF] Rosat | High resolution analysis of the gravest seismic normal modes after the 2004 Mw=9 Sumatra earthquake using superconducting gravimeter data[END_REF] the errors attributed to the frequency estimates of 0S2 are based on the signal-to-noise ratio (SNR) of the signal and are not formal uncertainties.

Table 3

 3 Estimated quality factors of 0S2 singlets from the 22-day Strasbourg SG record after the 2004 Sumatra earthquake using an amplitude decay measurement method, NLDHA and ESPRIT methods.

	Q	m = -2	m = -1	m = 0	m = 1	m = 2
	Amplitude decay measurement	440	493	545	484	555
	Uncertainties	6	11	14	11	5
	NLDHA	449.3	481.5	506.7	457.7	518.7
	Uncertainties	0.1	0.1	0.4	0.3	0.4
	ESPRIT	453.2	433.2	528.1	481.4	539.1
	Uncertainties	5	2.7	4.3	17.8	12.5
	PREM (Millot-Langet et al. 2003)	494.6	501.8	509.3	517.0	525.0

Table 4

 4 Estimated frequencies and quality factors of 0S3 singlets from the 22-day Strasbourg SG record after the 2004 Sumatra earthquake using NLDHA.

		0 S 3	m = -3	m = -2	m = -1	m = 0	m = 1	m = 2	m = 3
	NLDHA	frequency (mHz) Q	0.4615728 ± 4.9 10 -6 380.3 ± 0.04	0.4642270 ± 4.9 10 -6 412.7 ± 0.1	0.4664168 ± 4.6 10 -7 447.5 ± 0.5	--	0.4707937 ± 3.8 10 -6 446.1 ± 0.1	0.4727150 ± 6.9 10 -6 477.2 ± 0.2	0.4747854 ± 1.1 10 -6 328.2 ± 0.3
	Rosat et al. (2005)	frequency (mHz)	0.46167	0.46424	0.46639	-	0.47084	0.47266	0.47474
		frequency (mHz)	0.4618	0.4641	0.4664	0.4686	0.4707	0.4728	0.4748
	PREM	Q (Millot-Langet et al. 2003)	411.7	413.6	415.5	417.4	419.3	421.3	423.3

Acknowledgments

S. Rosat was supported by the Japan Society for the Promotion of Science during this study. This research is partly funded by the Japan Society for the Promotion of Science (JSPS) Grants-in-Aid for scientific research (B) 16340134.


Ion Lapteacru Larefi
email: ion.lapteacru@u-bordeaux.fr

Murphy-Topel adjustment of the variance-covariance matrix of a two-step panel data model: Evidence from competition-fragility nexus in banking

Keywords: JEL classification: Banking, Murphy-Topel adjustment, competition, risk C33, G21

We develop the Murphy-Topel adjustment of the variance-covariance matrix for two-step panel data models. We apply it on the competition-fragility nexus in banking with different samples for two equations. Indeed, this issue is often observed in this field of research. A competition measure of banks is constructed for each country (first equation), whereas a risk measure is regressed on the entire sample of countries (second equation). Any statistical adjustment will only provide approximate results for the second equation, because of possible correlations between the results of both models. The Murphy-Topel method eventually seems to be more appropriate.

Introduction

Many papers in economic empirical literature treat a subject with a two-step econometric model. The authors regress a first equation whose results are then applied in the computation of an explanatory variable that is used in the regression of a second equation. And many of them make these regressions independently of each other, considering therefore that this explanatory variable has no estimation error. Nevertheless, the errors of parameters from the first regression provide to this variable a random character and the values of parameters allow the determination of its expected value, which cannot be considered as an observed data unlike other explanatory variables. Consequently, since it depends on the regression results of the first model, the precision (the standard errors) of parameters of the second model must be adjusted.

Given this embedding of regressions and results, few authors merely bootstrap the standard errors of the second equation (for instance, [START_REF] Schaeck | Are competitive banking systems more stable?[END_REF][START_REF] Buch | Do Banks Benefit from Internationalization? Revisiting the Market Power-Risk Nexus[END_REF]). However, this statistical procedure does not take into account the results of the first equation, neither its parameters nor its variance-covariance matrix, and does not consider a possible correlation between these results and those of the second equation, thereby providing approximate results. Murphy andTopel (1985, 2002) propose an adjustment of the variancecovariance matrix for the second model and therefore of its standard errors. The aim of the present paper is to develop this procedure for panel data models with different samples for the first and second equations and to apply it to the estimation of the competition-fragility relationship in banking.

This subject is widely studied in empirical banking literature and all studies are made on panel data. However, very few of them bootstrapped the standard errors and none made a precise adjustment with the Murphy-Topel approach. The importance of the subject and the difficulty to estimate correctly the standard errors of the second equation, namely different samples for regression equations, lead us to develop the Murphy-Topel adjustment of the variance-covariance matrix of two-step panel data models, which seems more appropriate.

Two-step econometric model and Murphy-Topel variance-covariance matrix

Two-step econometric model consists of two equations, where the parameters of the first equation will appear in the second equation and the precision of their estimation, that is, the variance-covariance matrix, interferes therefore in the precision of the second equation's parameters.

Let Model 1: 𝑦𝑦 1 = 𝑓𝑓 1 (𝐱𝐱 1 , 𝛂𝛂) be the first equation, whose econometric results, that is, coefficients and variance-covariance matrix, are involved in the regression of the following second equation Model 2:

𝑦𝑦 2 = 𝑓𝑓 2 (𝐱𝐱 1 , 𝐱𝐱 2 , 𝛂𝛂 �, 𝛃𝛃),
where 𝐱𝐱 𝑖𝑖1 and 𝐱𝐱 𝑖𝑖2 are vectors of explanatory variables of the first and second models, respectively, and 𝛂𝛂 and 𝛃𝛃 are vectors of parameters of model 1 and model 2, respectively. Vector 𝛃𝛃 includes also the coefficient 𝛄𝛄 that connects to the second equation the explanatory variable computed with the parameters of the first equation.

Unlike many empirical studies that use two-step econometric models and consider the variables with coefficient 𝛄𝛄 as they contain observed data, one might consider them as estimators and adjust therefore the variance-covariance matrix of the second equation. Murphy andTopel (1985, 2002) proposes a full information maximum likelihood estimation with two-step procedure that consists to estimate, firstly, the parameters 𝛂𝛂 of model 1 by maximizing its likelihood function

ln𝐿𝐿 1 (𝛂𝛂) = � ln𝑓𝑓 1 (𝑦𝑦 1𝑖𝑖 |𝐱𝐱 1𝑖𝑖 , 𝛂𝛂) 𝑖𝑖 (1)
and, secondly, the parameters 𝛃𝛃 of model 2 by maximizing the respective likelihood function

ln𝐿𝐿 2 (𝛂𝛂 �, 𝛃𝛃) = � ln𝑓𝑓 2 (𝑦𝑦 2𝑖𝑖 |𝐱𝐱 1𝑖𝑖 , 𝐱𝐱 2𝑖𝑖 , 𝛂𝛂 �, 𝛃𝛃) 𝑖𝑖 , ( 2
)
where the embedded consistent estimator 𝛂𝛂 � is treated as given. Greene (2012, p. 576) explains the reasons why one might proceed with two-step estimations. Among them, the easiness to formulate and to maximize the two separate log-likelihoods leads to choose this approach, instead of the maximisation of a complicated joint distribution function. Thus, Murphy andTopel (1985, 2002) show that if the regression results can be obtained with MLE, then the variance-covariance matrix of the second equation must be adjusted in the following way

𝑉𝑉 2 * = 𝑉𝑉 2 + 𝑉𝑉 2 [𝐶𝐶𝑉𝑉 1 𝐶𝐶 ′ -𝑅𝑅𝑉𝑉 1 𝐶𝐶 ′ -𝐶𝐶𝑉𝑉 1 𝑅𝑅′]𝑉𝑉 2 ,
(3) where

𝑉𝑉 1 = �� �- 𝜕𝜕 2 ln𝑓𝑓 1,𝑖𝑖 𝜕𝜕𝛂𝛂 �𝜕𝜕𝛂𝛂 �′ � 𝑖𝑖 � -1 = �� � 𝜕𝜕ln𝑓𝑓 1,𝑖𝑖 𝜕𝜕𝛂𝛂 � � � 𝜕𝜕ln𝑓𝑓 1,𝑖𝑖 𝜕𝜕𝛂𝛂 �′ � 𝑖𝑖 � -1 (4)
is the variance-covariance matrix of model 1,

𝑉𝑉 2 = �� �- 𝜕𝜕 2 ln𝑓𝑓 2,𝑖𝑖 𝜕𝜕𝛃𝛃 � 𝜕𝜕𝛃𝛃 � ′ � 𝑖𝑖 � -1 = �� � 𝜕𝜕ln𝑓𝑓 2,𝑖𝑖 𝜕𝜕𝛃𝛃 � � � 𝜕𝜕ln𝑓𝑓 2,𝑖𝑖 𝜕𝜕𝛃𝛃 � ′ � 𝑖𝑖 � -1 (5)
is the variance-covariance matrix of model 2 when no adjustment is made, that is, when the constructed explanatory variable is considered with given observed data instead of computed econometric results, and

𝐶𝐶 = � �- 𝜕𝜕 2 ln𝑓𝑓 2,𝑖𝑖 𝜕𝜕𝛃𝛃 � 𝜕𝜕𝛂𝛂 �′ � 𝑖𝑖 = � � 𝜕𝜕ln𝑓𝑓 2,𝑖𝑖 𝜕𝜕𝛃𝛃 � � � 𝜕𝜕ln𝑓𝑓 2,𝑖𝑖 𝜕𝜕𝛂𝛂 �′ � 𝑖𝑖 , (6)
𝑅𝑅 = � � 𝜕𝜕ln𝑓𝑓 2,𝑖𝑖 𝜕𝜕𝛃𝛃 � � � 𝜕𝜕ln𝑓𝑓 1,𝑖𝑖 𝜕𝜕𝛂𝛂 �′ � 𝑖𝑖 (7)
are the matrices of cross products of first derivative vectors. Since the variance of the first derivative vector of the log-likelihood is the negative of the expected second derivative matrix, the two formulas for matrices 𝑉𝑉 1 , 𝑉𝑉 2 and 𝐶𝐶 provide the same result. As the covariance of the two first derivative vectors, it is very likely that the matrix 𝑅𝑅 has all elements very close to zero.

Application: competition-fragility nexus in banking

There is some difficulty to develop the Murphy-Topel adjustment for panel data models with different samples, which it is often the case of studies on the competition-fragility relationship.

Two-step panel data model

At the first stage, one might estimate the bank competition measure. Let us take the following examples for the first equation, which are widely applied in banking empirical literature. 1 st stage: Lerner index

The Lerner index is the mark-up of price over the bank's marginal cost:

𝐿𝐿 𝑖𝑖𝑖𝑖 = 𝑝𝑝 𝑖𝑖𝑖𝑖 - 𝜕𝜕𝐶𝐶 𝑖𝑖𝑖𝑖 𝜕𝜕𝑦𝑦 𝑖𝑖𝑖𝑖 𝑝𝑝 𝑖𝑖𝑖𝑖 , ( 8a
)
where 𝑝𝑝 𝑖𝑖𝑖𝑖 = 𝑇𝑇𝑅𝑅 𝑖𝑖𝑖𝑖 �𝑦𝑦 1,𝑖𝑖𝑖𝑖 + 𝑦𝑦 2,𝑖𝑖𝑖𝑖 + 𝑦𝑦 3,𝑖𝑖𝑖𝑖 � ⁄ is the average price of a bank i's output at time t, 𝐶𝐶 𝑖𝑖𝑖𝑖 is the bank's total cost and 𝑦𝑦 𝑖𝑖𝑖𝑖 is its output. Higher is the Lerner index, higher is the market power and less competitive banking institutions are.

The econometric task is thus to estimate the marginal costs of banks, 𝜕𝜕𝐶𝐶 𝑖𝑖𝑖𝑖 𝜕𝜕𝑦𝑦 𝑖𝑖𝑖𝑖 ⁄ . As in the most literature, we take three inputs, prices of which are 𝑤𝑤 𝑗𝑗 , and three outputs 𝑦𝑦 𝑚𝑚 :

ln�𝐶𝐶 𝑖𝑖𝑖𝑖 /𝑤𝑤 3,𝑖𝑖𝑖𝑖 � = 𝑐𝑐 0 + 𝛼𝛼 𝑖𝑖 + 𝜇𝜇 𝑖𝑖 + � 𝛽𝛽 𝑗𝑗 ln�𝑤𝑤 𝑗𝑗,𝑖𝑖𝑖𝑖 /𝑤𝑤 3,𝑖𝑖𝑖𝑖 � 2 𝑗𝑗=1 + � 𝛾𝛾 𝑚𝑚 ln𝑦𝑦 𝑚𝑚,𝑖𝑖𝑖𝑖 3 𝑚𝑚=1 + 1 2 � � 𝜑𝜑 𝑘𝑘𝑘𝑘 ln�𝑤𝑤 𝑘𝑘,𝑖𝑖𝑖𝑖 /𝑤𝑤 3,𝑖𝑖𝑖𝑖 �ln�𝑤𝑤 𝑘𝑘,𝑖𝑖𝑖𝑖 /𝑤𝑤 3,𝑖𝑖𝑖𝑖 � 2 𝑘𝑘=1 2 𝑘𝑘=1 + � � 𝜒𝜒 𝑚𝑚𝑘𝑘 ln�𝑤𝑤 𝑘𝑘,𝑖𝑖𝑖𝑖 /𝑤𝑤 3,𝑖𝑖𝑖𝑖 �ln𝑦𝑦 𝑚𝑚,𝑖𝑖𝑖𝑖 2 𝑘𝑘=1 3 𝑚𝑚=1 + 1 2 � � 𝜌𝜌 𝑚𝑚𝑚𝑚 ln𝑦𝑦 𝑚𝑚,𝑖𝑖𝑖𝑖 3 𝑚𝑚=1 ln𝑦𝑦 𝑚𝑚,𝑖𝑖𝑖𝑖 3 𝑚𝑚=1 + � 𝜔𝜔 1𝑚𝑚 ln𝑦𝑦 𝑚𝑚,𝑖𝑖𝑖𝑖 ln𝑔𝑔𝑔𝑔𝑝𝑝 𝑖𝑖 3 𝑚𝑚=1 + � 𝜔𝜔 2𝑚𝑚 ln𝑦𝑦 𝑚𝑚,𝑖𝑖𝑖𝑖 𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑝𝑝 𝑖𝑖 3 𝑚𝑚=1 + 𝜀𝜀 𝑖𝑖𝑖𝑖 , (8𝑏𝑏)
To control the effects of the business cycle, the real GDP, 𝑔𝑔𝑔𝑔𝑝𝑝 𝑖𝑖 , and/or the real GDP growth, 𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑝𝑝 𝑖𝑖 , can be considered in levels and/or in cross-product terms. The estimated coefficients of this equation are applied to compute the marginal cost:

𝑚𝑚𝑐𝑐 𝑖𝑖𝑖𝑖 = 𝜕𝜕𝐶𝐶 𝑖𝑖𝑖𝑖 𝜕𝜕�𝑦𝑦 1,𝑖𝑖𝑖𝑖 + 𝑦𝑦 2,𝑖𝑖𝑖𝑖 + 𝑦𝑦 3,𝑖𝑖𝑖𝑖 � = � � 𝛾𝛾 𝑚𝑚 3 𝑚𝑚=1 + � � 𝜒𝜒 𝑚𝑚𝑘𝑘 ln�𝑤𝑤 𝑘𝑘,𝑖𝑖𝑖𝑖 /𝑤𝑤 3,𝑖𝑖𝑖𝑖 � 2 𝑘𝑘=1 + � � 𝜌𝜌 𝑚𝑚𝑚𝑚 ln𝑦𝑦 𝑚𝑚 3 𝑚𝑚=1 3 𝑚𝑚=1 3 𝑚𝑚=1 + � 𝜔𝜔 1𝑚𝑚 ln𝑔𝑔𝑔𝑔𝑝𝑝 𝑖𝑖 3 𝑚𝑚=1 + � 𝜔𝜔 2𝑚𝑚 𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑝𝑝 𝑖𝑖 3 𝑚𝑚=1 � × 𝐶𝐶 𝑖𝑖𝑖𝑖 𝑦𝑦 1,𝑖𝑖𝑖𝑖 + 𝑦𝑦 2,𝑖𝑖𝑖𝑖 + 𝑦𝑦 3,𝑖𝑖𝑖𝑖 . (8𝑐𝑐)
1 st stage: Boone indicator The Boone indicator, proposed by [START_REF] Boone | A new way to measure competition[END_REF], is the elasticity of market share to marginal cost. Based on the idea that more efficient banks, that is, with lower marginal costs, gain higher market shares, this effect will be stronger in more competitive banking market. One needs to run the following econometric model:

ln𝑚𝑚𝑚𝑚 𝑖𝑖𝑖𝑖 = 𝜋𝜋 0 + 𝜋𝜋 1𝑖𝑖 + 𝜈𝜈 𝑖𝑖 + 𝜃𝜃 1 ln𝑚𝑚𝑐𝑐 𝑖𝑖𝑖𝑖 + 𝜃𝜃 2 ln𝑚𝑚𝑐𝑐 𝑖𝑖𝑖𝑖 × ln𝑚𝑚𝑐𝑐 𝑖𝑖𝑖𝑖 + 𝜃𝜃 3 ln𝑚𝑚𝑐𝑐 𝑖𝑖𝑖𝑖 × ln𝑔𝑔𝑔𝑔𝑝𝑝 𝑖𝑖 +𝜃𝜃 4 ln𝑚𝑚𝑐𝑐 𝑖𝑖𝑖𝑖 × 𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑝𝑝 𝑖𝑖 + 𝜀𝜀 𝑖𝑖𝑖𝑖 , (9
) where 𝑚𝑚𝑚𝑚 𝑖𝑖𝑖𝑖 stands for earning assets share for bank i in year t. As the elasticity of market share to marginal cost, the Boone indicator is calculated as following:

𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵 𝑖𝑖𝑖𝑖 = 𝜃𝜃 1 + 2𝜃𝜃 2 ln𝑚𝑚𝑐𝑐 𝑖𝑖𝑖𝑖 + 𝜃𝜃 3 ln𝑔𝑔𝑔𝑔𝑝𝑝 𝑖𝑖 + 𝜃𝜃 4 𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑝𝑝 𝑖𝑖 .
(10) Since 𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵 might take negative values, higher values (or lower absolute values) imply less competitive behaviour. 1 st stage: H-statistic

The H-statistic of [START_REF] Panzar | Testing for "monopoly" equilibrium[END_REF] is the aggregate elasticity of banks' revenue to input prices. It can therefore be computed regressing the following translog revenue function of banks:

ln𝑇𝑇𝑅𝑅 𝑖𝑖𝑖𝑖 = a 0 + ϑ i + 𝜐𝜐 𝑖𝑖 + � 𝜆𝜆 𝑗𝑗 ln𝑤𝑤 𝑗𝑗,𝑖𝑖𝑖𝑖 3 𝑗𝑗=1 + 𝜓𝜓 1 ln𝑇𝑇𝑇𝑇 𝑖𝑖𝑖𝑖 + 𝜓𝜓 2 ln𝐿𝐿𝐵𝐵𝐿𝐿𝐵𝐵 𝑖𝑖𝑖𝑖 + 1 2 � � 𝜅𝜅 𝑘𝑘𝑘𝑘 ln𝑤𝑤 𝑘𝑘,𝑖𝑖𝑖𝑖 ln𝑤𝑤 𝑘𝑘,𝑖𝑖𝑖𝑖 3 𝑘𝑘=1 3 𝑘𝑘=1 + � 𝜙𝜙 1𝑘𝑘 ln𝑤𝑤 𝑘𝑘,𝑖𝑖𝑖𝑖 ln𝑇𝑇𝑇𝑇 𝑖𝑖𝑖𝑖 3 𝑘𝑘=1 + � 𝜙𝜙 2𝑘𝑘 ln𝑤𝑤 𝑘𝑘,𝑖𝑖𝑖𝑖 𝐿𝐿𝐵𝐵𝐿𝐿𝐵𝐵 𝑖𝑖𝑖𝑖 3 𝑘𝑘=1 + � 𝜙𝜙 3𝑘𝑘 ln𝑤𝑤 𝑘𝑘,𝑖𝑖𝑖𝑖 ln𝑔𝑔𝑔𝑔𝑝𝑝 𝑖𝑖 3 𝑘𝑘=1 + � 𝜙𝜙 4𝑘𝑘 ln𝑤𝑤 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑝𝑝 𝑖𝑖 3 𝑘𝑘=1 + 𝜀𝜀 𝑖𝑖𝑖𝑖 , ( 11
)
where 𝑇𝑇𝑅𝑅 𝑖𝑖𝑖𝑖 stands for total revenue of bank i in year t, 𝑇𝑇𝑇𝑇 𝑖𝑖𝑖𝑖 for total assets and 𝐿𝐿𝐵𝐵𝐿𝐿𝐵𝐵 𝑖𝑖𝑖𝑖 for ratio of loans on total assets. As for the Lerner index, we consider three inputs for bank intermediation function. Usually in banking empirical literature, the use of 𝑇𝑇𝑇𝑇 and 𝐿𝐿𝐵𝐵𝐿𝐿𝐵𝐵 variables is to take into account the differences in size and business model, respectively, akin to De Bandt and Davis (2000), [START_REF] Bikker | Competition, concentration and their relationship: an empirical analysis of the banking industry[END_REF], [START_REF] Weill | Bank competition in the EU: How has it evolved?[END_REF]. And as for previous two measures, the economic state variables are considered to control the business cycles. Thus, the H-statistic is equal to

𝐻𝐻 𝑖𝑖𝑖𝑖 = � 𝜕𝜕𝑇𝑇𝑅𝑅 𝑖𝑖𝑖𝑖 𝜕𝜕𝑤𝑤 𝑘𝑘,𝑖𝑖𝑖𝑖 × 𝑤𝑤 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑖𝑖𝑖𝑖 3 𝑘𝑘=1 = � 𝜆𝜆 𝑗𝑗 3 𝑘𝑘=1 + � � 𝜅𝜅 𝑘𝑘𝑘𝑘 ln𝑤𝑤 𝑘𝑘,𝑖𝑖𝑖𝑖 3 𝑘𝑘=1 3 𝑘𝑘=1 + � 𝜙𝜙 1𝑘𝑘 ln𝑇𝑇𝑇𝑇 𝑖𝑖𝑖𝑖 3 𝑘𝑘=1 + � 𝜙𝜙 2𝑘𝑘 𝐿𝐿𝐵𝐵𝐿𝐿𝐵𝐵 𝑖𝑖𝑖𝑖 3 𝑘𝑘=1 + � 𝜙𝜙 3𝑘𝑘 ln𝑔𝑔𝑔𝑔𝑝𝑝 𝑖𝑖 3 𝑘𝑘=1 + � 𝜙𝜙 4𝑘𝑘 𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑝𝑝 𝑖𝑖 3 𝑘𝑘=1 . ( 12
)
H-statistic takes values between zero, monopoly or oligopoly situation, and one, perfectly competitive behaviour.

All regressions can be made individually on each country 𝑘𝑘 of the sample. However, since the risk measure cannot be computed for many of banks, namely for non-listed banks if a market-based measure is employed, the second equation is often regressed on entire sample of 𝑁𝑁 𝑐𝑐 countries.

At the second stage, in order to estimate the relationship between bank competition and risk-taking, two different panel data models are considered in the banking literature: 2 nd stage: linear relationship where 𝑅𝑅𝑅𝑅𝑚𝑚𝑘𝑘 𝑖𝑖,𝑗𝑗,𝑖𝑖 is a risk measure, n is the number of bank-specific variables, m-n is the number of country-specific variables, 𝐵𝐵𝐶𝐶 𝑖𝑖,𝑗𝑗,𝑖𝑖 and 𝑋𝑋 𝑘𝑘,𝑖𝑖,𝑗𝑗,𝑖𝑖 are, respectively, bank competition measure and bank-specific factor that controls for the k-th characteristic of bank i at year t in country j, and 𝑌𝑌 𝑘𝑘,𝑗𝑗,𝑖𝑖 is the country-specific factor that controls for the l-th characteristic of country j.

𝑅𝑅𝑅𝑅𝑚𝑚𝑘𝑘 𝑖𝑖,𝑗𝑗,𝑖𝑖 = 𝛽𝛽 0 + 𝛾𝛾 1 𝐵𝐵𝐶𝐶 𝑖𝑖,𝑗𝑗,𝑖𝑖 + � 𝛽𝛽 𝑘𝑘 𝑛𝑛 𝑘𝑘=1 𝑋𝑋 𝑘𝑘,𝑖𝑖,𝑗𝑗,𝑖𝑖 + � 𝛽𝛽 𝑘𝑘 𝑚𝑚 𝑘𝑘=𝑛𝑛+1 𝑌𝑌 𝑘𝑘,𝑗𝑗,𝑖𝑖 + 𝜇𝜇 𝑖𝑖 + 𝜋𝜋 𝑖𝑖 + 𝜀𝜀 𝑖𝑖,𝑗𝑗,𝑖𝑖 (
In order to determine the adjusted variance-covariance matrix 𝑉𝑉 2 * of our risk equation, the variance-covariance matrices of the first and second equation 𝑉𝑉 1 and 𝑉𝑉 2 , respectively, are extracted from the respective regressions and we must compute only the matrices 𝐶𝐶 and 𝑅𝑅.

The matrix 𝑉𝑉 1 is of dimension 𝑙𝑙 1 𝑁𝑁 𝑐𝑐 × 𝑙𝑙 1 𝑁𝑁 𝑐𝑐 and composed on the diagonal by variancecovariance matrices 𝑉𝑉 1,𝑘𝑘 for each country 𝑘𝑘 (𝑘𝑘 = 1, … , 𝑁𝑁 𝑐𝑐 ) and all other elements are zeros:

𝑉𝑉 1 = � 𝑉𝑉 1,1 0 ⋮ 0 0 𝑉𝑉 1,2 ⋮ 0 … … ⋱ … 0 0 ⋮ 𝑉𝑉 1,𝑁𝑁 𝑐𝑐 �.
The estimations of competition measures for all countries are therefore taken into consideration in the adjustment of the variance-covariance matrix 𝑉𝑉 2 , which is of dimension 𝑙𝑙 2 × 𝑙𝑙 2 . 𝑙𝑙 2 and 𝑙𝑙 1 are the size of vectors 𝐱𝐱 2 and 𝐱𝐱 1 , respectively. These vectors of explanatory variables include also the intercept and the time fixed effects.

The difficulty of computation of matrices 𝐶𝐶 and 𝑅𝑅 is the difference between samples of two equations. They are determined as cross "products", for matrix 𝑅𝑅, or cross derivations, for matrix 𝐶𝐶, between the single sample of the second equation and the bank sample of each country. The computation of the elements 𝑔𝑔[𝑅𝑅, 𝑗𝑗] and 𝑐𝑐[𝑅𝑅, 𝑗𝑗] of matrices 𝑅𝑅 and 𝐶𝐶, respectively, are explained in next sub-section.

Murphy-Topel adjustment of the variance-covariance matrix of the risk equation

For the first panel data model (eq. 8b, for the Lerner index, eq. 9 for the Boone indicator, and eq. 11 for the H-statistic), the likelihood is formulated for each country 𝑘𝑘 of the sample.

ln𝑓𝑓 1,𝑘𝑘,𝑖𝑖 = - 1 2 � 1 𝜎𝜎𝜎𝜎 1,𝑘𝑘 2 � � 𝜀𝜀 1,𝑘𝑘,𝑖𝑖𝑖𝑖 2 𝑇𝑇 1,𝑘𝑘,𝑖𝑖 𝑖𝑖=1 - 𝜎𝜎𝐵𝐵 1,𝑘𝑘 2 𝜎𝜎𝐵𝐵 1,𝑘𝑘 2 + 𝑇𝑇 1,𝑘𝑘,𝑖𝑖 𝜎𝜎𝜎𝜎 1,𝑘𝑘 2 � � 𝜀𝜀 1,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 1,𝑘𝑘,𝑖𝑖 𝑖𝑖=1 � 2 � + ln �1 + 𝑇𝑇 1,𝑘𝑘,𝑖𝑖 𝜎𝜎𝜎𝜎 1,𝑘𝑘 2 𝜎𝜎𝐵𝐵 1,𝑘𝑘 2 �� - 1 2 �𝑇𝑇 1,𝑘𝑘,𝑖𝑖 ln�2π𝜎𝜎𝐵𝐵 1,𝑘𝑘 2 �� ( 14
)
𝑘𝑘 ranges from 1 to 𝑁𝑁 𝑐𝑐 . 𝜎𝜎𝐵𝐵 1,𝑘𝑘 2 and 𝜎𝜎𝜎𝜎 1,𝑘𝑘 2 are the variances of residuals and individual effects, respectively, for each panel data regression 𝑘𝑘. 𝜀𝜀 1,𝑘𝑘,𝑖𝑖𝑖𝑖 are residuals of this regression and 𝑇𝑇 1,𝑘𝑘,𝑖𝑖 number of observations for the bank 𝑅𝑅. The log-likelihood of the second model takes therefore the following form

ln𝑓𝑓 2,𝑖𝑖 = - 1 2 � 1 𝜎𝜎𝜎𝜎 2 2 �� 𝜀𝜀 2,𝑖𝑖𝑖𝑖 2 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 - 𝜎𝜎𝐵𝐵 2 2 𝜎𝜎𝐵𝐵 2 2 + 𝑇𝑇 2,𝑖𝑖 𝜎𝜎𝜎𝜎 2 2 �� 𝜀𝜀 2,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 � 2 � + ln �1 + 𝑇𝑇 2,𝑖𝑖 𝜎𝜎𝜎𝜎 2 2 𝜎𝜎𝐵𝐵 2 2 �� - 1 2 �𝑇𝑇 2,𝑖𝑖 ln(2π𝜎𝜎𝐵𝐵 2 2 )�, ( 15
)
where 𝜎𝜎𝐵𝐵 2 2 and 𝜎𝜎𝜎𝜎 2 2 are the variances of residuals and of individual effects, respectively. 𝜀𝜀 2,𝑖𝑖𝑖𝑖 are residuals of this regression and 𝑇𝑇 2,𝑖𝑖 number of observations for the bank 𝑅𝑅.

To compute the elements of matrices 𝑅𝑅 and 𝐶𝐶, we determine the following derivatives: For the first equation: Consequently, the elements of 𝑅𝑅 matrix for each country 𝑘𝑘 are

𝜕𝜕ln𝑓𝑓 1,𝑘𝑘,𝑖𝑖 𝜕𝜕𝛂𝛂 � 𝑘𝑘 = 1 𝜎𝜎𝐵𝐵 1,𝑘𝑘 2 � � 𝜀𝜀
𝑔𝑔0[𝑘𝑘, 𝑙𝑙 2 , 𝑙𝑙 1 ] = � � 1 𝜎𝜎𝐵𝐵 2 2 �� 𝜀𝜀 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝐱𝐱[𝑙𝑙 2 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 - 𝜎𝜎𝐵𝐵 2 2 𝜎𝜎𝐵𝐵 2 2 + 𝑇𝑇 2,𝑖𝑖 𝜎𝜎𝜎𝜎 2 2 �� 𝜀𝜀 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 � �� 𝐱𝐱[𝑙𝑙 2 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 �� 𝑁𝑁 2 𝑖𝑖=1 × 1 𝜎𝜎𝐵𝐵 1,𝑘𝑘 2 � � 𝜀𝜀 1,𝑘𝑘,𝑖𝑖𝑖𝑖 𝐱𝐱[𝑙𝑙 1 ] 1,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 1,𝑘𝑘,𝑖𝑖 𝑖𝑖=1 - 𝜎𝜎𝐵𝐵 1,𝑘𝑘 2 𝜎𝜎𝐵𝐵 1,𝑘𝑘 2 + 𝑇𝑇 1,𝑘𝑘,𝑖𝑖 𝜎𝜎𝜎𝜎 1,𝑘𝑘 2 � � 𝜀𝜀 1,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 1,𝑘𝑘,𝑖𝑖 𝑖𝑖=1 � � � 𝐱𝐱[𝑙𝑙 1 ] 1,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 1,𝑘𝑘,𝑖𝑖 𝑖𝑖=1 ��� , ( 18
)
where 𝑁𝑁 2 is the number of banks of the second model. Finally, the matrix 𝑅𝑅 has elements 𝑔𝑔[𝑅𝑅, 𝑗𝑗] = 𝑔𝑔0[𝑘𝑘, 𝑙𝑙 2 , 𝑙𝑙 1 ], where 𝑅𝑅 = 𝑙𝑙 2 and 𝑗𝑗 = 𝑘𝑘 × 𝑙𝑙 1 .

Besides the matrix 𝑅𝑅, the results of the first equation are also considered in the second one with the matrix 𝐶𝐶. Since we have two versions of the risk equation, this second derivative considers only 𝐵𝐵𝐶𝐶 variable for the linear relationship equation (eq. 13a) and both 𝐵𝐵𝐶𝐶 and 𝐵𝐵𝐶𝐶 2 variables for non-linear relationship equation (eq. 13b) as follows: For a liner model: -The derivative of the 𝐵𝐵𝐶𝐶 variable with respect to parameters of the first model

𝑐𝑐0[𝑘𝑘, 1, 𝑙𝑙 1 ] = -� � 1 𝜎𝜎𝐵𝐵 2 2 ���𝜀𝜀 2,𝑘𝑘,𝑖𝑖𝑖𝑖 -𝛾𝛾 1 𝐵𝐵𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 �𝐱𝐱′[𝑙𝑙 1 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 𝑁𝑁 2 𝑖𝑖=1 - 𝜎𝜎𝐵𝐵 2 2 𝜎𝜎𝐵𝐵 2 2 + 𝑇𝑇 2,𝑖𝑖 𝜎𝜎𝜎𝜎 2 2 �� 𝜀𝜀 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 -𝛾𝛾 1 � 𝐵𝐵𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 � �� 𝐱𝐱′[𝑙𝑙 1 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 ���, (19𝐿𝐿)
-The derivative of other variables of the second equation with respect to parameters of the first model

𝑐𝑐0[𝑘𝑘, 𝑙𝑙 2 , 𝑙𝑙 1 ] = -� � 1 𝜎𝜎𝐵𝐵 2 2 �� -𝛾𝛾 1 𝐱𝐱′[𝑙𝑙 1 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝐱𝐱[𝑙𝑙 2 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 𝑁𝑁 2 𝑖𝑖=1 - 𝜎𝜎𝐵𝐵 2 2 𝜎𝜎𝐵𝐵 2 2 + 𝑇𝑇 2,𝑖𝑖 𝜎𝜎𝜎𝜎 2 2 �� -𝛾𝛾 1 𝐱𝐱′[𝑙𝑙 1 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 � �� 𝐱𝐱[𝑙𝑙 2 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 ���, (19𝑏𝑏)
For a non-linear model:

-The derivative of the 𝐵𝐵𝐶𝐶 variable with respect to parameters of the first model

𝑐𝑐0[𝑘𝑘, 1, 𝑙𝑙 1 ] = -� � 1 𝜎𝜎𝐵𝐵 2 2 ���𝜀𝜀 2,𝑘𝑘,𝑖𝑖𝑖𝑖 + �-𝛾𝛾 1 -2𝛾𝛾 2 𝐵𝐵𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 �𝐵𝐵𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 �𝐱𝐱′[𝑙𝑙 1 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 𝑁𝑁 2 𝑖𝑖=1 - 𝜎𝜎𝐵𝐵 2 2 𝜎𝜎𝐵𝐵 2 2 + 𝑇𝑇 2,𝑖𝑖 𝜎𝜎𝜎𝜎 2 2 ����-𝛾𝛾 1 -2𝛾𝛾 2 𝐵𝐵𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 �𝐱𝐱′[𝑙𝑙 1 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 � �� 𝐵𝐵𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 � + �� 𝜀𝜀 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 � �� 𝐱𝐱′[𝑙𝑙 1 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 ����, (20𝐿𝐿)
-The derivative of the 𝐵𝐵𝐶𝐶 2 variable with respect to parameters of the first model where 𝑅𝑅 = 𝑙𝑙 2 and 𝑗𝑗 = 𝑘𝑘 × 𝑙𝑙 1 , in both cases.

𝑐𝑐0[𝑘𝑘, 2, 𝑙𝑙 1 ] = -� � 1 𝜎𝜎𝐵𝐵 2 2 ���2𝐵𝐵𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝜀𝜀 2,𝑘𝑘,𝑖𝑖𝑖𝑖 + �-𝛾𝛾 1 -2𝛾𝛾 2 𝐵𝐵𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 �𝐵𝐵𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 2 �𝐱𝐱′[𝑙𝑙 1 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 𝑁𝑁 2 𝑖𝑖=1 - 𝜎𝜎𝐵𝐵 2 2 𝜎𝜎𝐵𝐵 2 2 + 𝑇𝑇 2,𝑖𝑖 𝜎𝜎𝜎𝜎 2 2 ����-𝛾𝛾 1 -2𝛾𝛾 2 𝐵𝐵𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 �𝐱𝐱′[𝑙𝑙 1 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 � �� 𝐵𝐵𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 2 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 � + �� 𝜀𝜀 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 � �� 2𝐵𝐵𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝐱𝐱′[𝑙𝑙 1 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇 2,𝑖𝑖 𝑖𝑖=1 ����, (
The difficulty to construct the matrices 𝑅𝑅 and 𝐶𝐶 is the difference in sample size of two equations: the first equation is regressed for each country 𝑘𝑘 and the second equation for entire sample of countries. To overcome it, we transformed the explanatory variables, 𝐱𝐱[𝑙𝑙 2 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 , and the error term, 𝜀𝜀 2,𝑘𝑘,𝑖𝑖𝑖𝑖 , of the second equation as that each country of the sample has its own variables. Each 𝑘𝑘 variable (or error term) takes therefore its own values if they correspond to country 𝑘𝑘 and zero otherwise.

The vector 𝐱𝐱′[𝑙𝑙 1 ] 2 is the derivative of bank competition measure 𝐵𝐵𝐶𝐶, estimated in the first step, with respect to parameters 𝑙𝑙 1 of the first equation. The index 2 signifies that this vector is obtained from the derivative of the log-likelihood of the second equation. The derivation is firstly made with respect to parameters of the explanatory variables and then with respect to parameters of fixed time effects and intercept. 𝐱𝐱′[𝑙𝑙 1 ] 2 vector is therefore individual for each of our three non-structural competition measures: -For the Lerner index, based on equations (8a)-(8c), and the definition of the average price of bank output The size of this vector is the same as that of the explanatory variables of equation (8b), that is, 27+ 𝑡𝑡𝐵𝐵, and includes 26 explanatory variables, 𝑡𝑡𝐵𝐵 time fixed effects and 1 intercept.

𝐱𝐱′[𝑙𝑙 1 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 = �0,
-For the Boone indicator, based on equations ( 9) and ( 10)

𝐱𝐱′[𝑙𝑙 1 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 = �1,2ln𝑚𝑚𝑐𝑐 𝑘𝑘,𝑖𝑖𝑖𝑖 , ln𝑔𝑔𝑔𝑔𝑝𝑝 𝑘𝑘,𝑖𝑖𝑖𝑖 , 𝑔𝑔𝑔𝑔𝑝𝑝𝑔𝑔𝑔𝑔 𝑘𝑘,𝑖𝑖𝑖𝑖 , 0, … ,0 ��� 𝑖𝑖𝜎𝜎 , 0�
The size of this vector is the same as that of the explanatory variables of equation ( 9), that is, 5+ 𝑡𝑡𝐵𝐵, and includes 4 explanatory variables, 𝑡𝑡𝐵𝐵 time fixed effects and 1 intercept.

-For the H-statistic, based on equations ( 11) and ( 12) The size of this vector is the same as that of the explanatory variables of equation ( 11), that is, 24+ 𝑡𝑡𝐵𝐵, and includes 23 explanatory variables, 𝑡𝑡𝐵𝐵 time fixed effects and 1 intercept.

𝐱𝐱′[𝑙𝑙 1 ] 2,𝑘𝑘,𝑖𝑖𝑖𝑖 = �1,1,1,0,0,

Conclusion

In this paper, we developed the Murphy-Topel adjustment of the variance-covariance matrix of a two-step panel data model and for different samples of two equations, and we applied it on the example of competition-stability relationship in banking. This subject is widely studied in banking empirical literature and this econometric issue is therefore present. Because of possible correlations between the results of both models, the Murphy-Topel method eventually seems to be more appropriate than any statistical (numerical adjustment) as the bootstrapping procedure.

RECENT LAREFI WORKING PAPERS

  20𝑏𝑏) and the derivative of other variables of the second equation with respect to parameters of the first model is the same equation (19b). The matrix 𝐶𝐶 has the elements 𝑐𝑐[𝑅𝑅, 𝑗𝑗] = 𝑐𝑐0[𝑘𝑘, 𝑙𝑙 2 , 𝑙𝑙 1 ],

For

  earlier LAREFI Working Papers, please go to http://lare-efi.org.u-bordeaux4.fr All Discussion Papers can be downloaded free of charge CR16-EFI07 Faure, P-H., "Le bitcoin peut-il être assimilé à une monnaie ? Un examen à partir des différentes grilles de lecture de la science économique" CR16-EFI06 Lapteacru, I., "Convergence of bank competition in Central and Eastern European countries : Does ownership matter ?" CR16-EFI05 Lapteacru, I., "On the consistency of the Z-score to measure the bank risk" CR16-EFI04 Lapteacru, I., "Income and funding structures, banking regulation and bank risk-taking : The role of ownership in Central and Eastern European banks" CR16-EFI03 Humblot, T., "Distance et accès au crédit sous Bâle III" CR16-EFI02 Berthoumieu, J. and V. Lamani, "Vertical Differentiation, Uncertainty, Product R&D and Policy Instruments in a North-South Duopoly" CR16-EFI01 Raffestin, L., "Foreign exchange investment rules and endogenous currency crashes" CR15-EFI05 Berthoumieu,J., "Technology Diffusion via Patent Collaborations : The Case of European Integration" CR15-EFI04 Bouet, A. and A-G. Vaubourg, "Financial Constraints and International Trade with Endogenous Mode of Competition" CR15-EFI03 Lahet, D. and A-G. Vaubourg, "Banks' shareholding in multilateral trading facilities : A two-sided market perspective" CR15-EFI02 Berthoumieu, J., "Policy Instruments, Patents and International Technology Diffusion in a North-South Duopoly" CR15-EFI01 Dupuy, L., "International trade and structural change : a dynamic model of weak sustainability"

  𝑘𝑘,𝑖𝑖𝑖𝑖 + ln𝑦𝑦 1,𝑘𝑘,𝑖𝑖𝑖𝑖 �, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄ × �ln𝑦𝑦 1,𝑘𝑘,𝑖𝑖𝑖𝑖 + ln𝑦𝑦 2,𝑘𝑘,𝑖𝑖𝑖𝑖 �, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄ × �ln𝑦𝑦 1,𝑘𝑘,𝑖𝑖𝑖𝑖 + ln𝑦𝑦 3,𝑘𝑘,𝑖𝑖𝑖𝑖 �, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄ × �ln𝑦𝑦 2,𝑘𝑘,𝑖𝑖𝑖𝑖 + ln𝑦𝑦 2,𝑘𝑘,𝑖𝑖𝑖𝑖 �, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄ × �ln𝑦𝑦 2,𝑘𝑘,𝑖𝑖𝑖𝑖 + ln𝑦𝑦 3,𝑘𝑘,𝑖𝑖𝑖𝑖 �, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄ × �ln𝑦𝑦 3,𝑘𝑘,𝑖𝑖𝑖𝑖 + ln𝑦𝑦 3,𝑘𝑘,𝑖𝑖𝑖𝑖 �, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄ × ln𝑔𝑔𝑔𝑔𝑝𝑝 𝑘𝑘,𝑖𝑖𝑖𝑖 , -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄ × ln𝑔𝑔𝑔𝑔𝑝𝑝 𝑘𝑘,𝑖𝑖𝑖𝑖 , -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄ × ln𝑔𝑔𝑔𝑔𝑝𝑝 𝑘𝑘,𝑖𝑖𝑖𝑖 , -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄ × 𝑔𝑔𝑔𝑔𝑝𝑝𝑔𝑔𝑔𝑔 𝑘𝑘,𝑖𝑖𝑖𝑖 , -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄ × 𝑔𝑔𝑔𝑔𝑝𝑝𝑔𝑔𝑔𝑔 𝑘𝑘,𝑖𝑖𝑖𝑖 , -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄ × 𝑔𝑔𝑔𝑔𝑝𝑝𝑔𝑔𝑔𝑔 𝑘𝑘,𝑖𝑖𝑖𝑖 , 0, … ,0 ���

	0, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄	, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄	, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄	, 0,0, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄
	× ln�𝑤𝑤 1,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑤𝑤 3,𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄	�, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄	× ln�𝑤𝑤 1,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑤𝑤 3,𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄	�, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄
	× ln�𝑤𝑤 1,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑤𝑤 3,𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄	�, 0, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄	× ln�𝑤𝑤 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑤𝑤 3,𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄	�, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄
	× ln�𝑤𝑤 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑤𝑤 3,𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄	�, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄	× ln�𝑤𝑤 2,𝑘𝑘,𝑖𝑖𝑖𝑖 𝑤𝑤 3,𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄	�, -𝐶𝐶 𝑘𝑘,𝑖𝑖𝑖𝑖 𝑇𝑇𝑅𝑅 𝑘𝑘,𝑖𝑖𝑖𝑖 ⁄
	× �ln𝑦𝑦 1,

𝑖𝑖𝜎𝜎

, 0�

  ln𝑤𝑤 1,𝑘𝑘,𝑖𝑖𝑖𝑖 + ln𝑤𝑤 1,𝑘𝑘,𝑖𝑖𝑖𝑖 , ln𝑤𝑤 1,𝑘𝑘,𝑖𝑖𝑖𝑖 + ln𝑤𝑤 2,𝑘𝑘,𝑖𝑖𝑖𝑖 , ln𝑤𝑤 1,𝑘𝑘,𝑖𝑖𝑖𝑖 + ln𝑤𝑤 3,𝑘𝑘,𝑖𝑖𝑖𝑖 , ln𝑤𝑤 2,𝑘𝑘,𝑖𝑖𝑖𝑖 + ln𝑤𝑤 2,𝑘𝑘,𝑖𝑖𝑖𝑖 , ln𝑤𝑤 2,𝑘𝑘,𝑖𝑖𝑖𝑖 + ln𝑤𝑤 3,𝑘𝑘,𝑖𝑖𝑖𝑖 , ln𝑤𝑤 3,𝑘𝑘,𝑖𝑖𝑖𝑖 + ln𝑤𝑤 3,𝑘𝑘,𝑖𝑖𝑖𝑖 , ln𝑇𝑇𝑇𝑇 𝑘𝑘,𝑖𝑖𝑖𝑖 , ln𝑇𝑇𝑇𝑇 𝑘𝑘,𝑖𝑖𝑖𝑖 , ln𝑇𝑇𝑇𝑇 𝑘𝑘,𝑖𝑖𝑖𝑖 , 𝐿𝐿𝐵𝐵𝐿𝐿𝐵𝐵 𝑘𝑘,𝑖𝑖𝑖𝑖 , 𝐿𝐿𝐵𝐵𝐿𝐿𝐵𝐵 𝑘𝑘,𝑖𝑖𝑖𝑖 , 𝐿𝐿𝐵𝐵𝐿𝐿𝐵𝐵 𝑘𝑘,𝑖𝑖𝑖𝑖 , ln𝑔𝑔𝑔𝑔𝑝𝑝 𝑘𝑘,𝑖𝑖𝑖𝑖 , ln𝑔𝑔𝑔𝑔𝑝𝑝 𝑘𝑘,𝑖𝑖𝑖𝑖 , ln𝑔𝑔𝑔𝑔𝑝𝑝 𝑘𝑘,𝑖𝑖𝑖𝑖 , 𝑔𝑔𝑔𝑔𝑝𝑝𝑔𝑔𝑔𝑔 𝑘𝑘,𝑖𝑖𝑖𝑖 , 𝑔𝑔𝑔𝑔𝑝𝑝𝑔𝑔𝑔𝑔 𝑘𝑘,𝑖𝑖𝑖𝑖 , 𝑔𝑔𝑔𝑔𝑝𝑝𝑔𝑔𝑔𝑔 𝑘𝑘,𝑖𝑖𝑖𝑖 , 0, … ,0 ���

𝑖𝑖𝜎𝜎

, 0�


Mickaël Royer

Alain Pirovano

Garcia Fabien

Fabien Garcia

Performance of Context-Aware Publish/Subscribe systems for AANET

Keywords: Vehicular Ad Hoc Network (VANET), Aeronautical Ad Hoc Network (AANET), Publish / Subscribe communication systems

Introduction

Mobile Ad-Hoc Networks (MANET) are self-conguring networks of mobile nodes connected by wireless links without a xed infrastructure. Vehicular Ad-Hoc Networks (VANET) are a subclass of MANET in which nodes are vehicle (cars, trains...). This sub-category shows a lot of specic characteristics, the most important being the high mobility of nodes and the predictability of movements. As detailed in [START_REF] Royer | Survey on Context-Aware Publish/Subscribe Systems for VANET[END_REF], AANET (Aeronautical Ad hoc NETworks, which consist in aircraft acting as nodes of a MANET) can be seen as a subcategory of VANET.

In a context of increasing air trac and new needs in term of air ground communication systems, AANET represents a relevant solution as explained in [START_REF] Besse | Wireless Ad Hoc Networks Access For Aeronautical Communications[END_REF]. Content-based information dissemination enjoys wide applicability in AANET where resources (like bandwidth) are limited. Notication of dangerous weather events to aircraft which will navigate close to the concerned geographical area or important messages from airline operational center to aircraft belonging to the company are good examples. Indeed, this approach allows a space decoupling between source and recipients since the sender of a message does not know which nodes will receive it. The aim of the network is then to nd the recipients, based on the message's content. To accomplish content-based routing in a network, the Publish / Subscribe communication paradigm is the most widespread solution. In this model, event distribution from publisher (event producer) to subscriber (event consumer) is based on the subscriber interest and the event's content. Many Publish / Subscribe systems exist for wired network ( [START_REF] Antony | SCRIBE: The Design of a Large-Scale Event Notication Infrastructure[END_REF], [START_REF] Pietzuch | Hermes: a distributed event-based middleware architecture[END_REF]) but are not suitable for MANET with low bandwidth and frequent network partition. Moreover, the applications described above need two extra features compared to classical Publish / Subscribe systems. The rst one is the persistence of events in the network to reduce the problem of network partitions (similar principle as Delay Tolerant Networks, [START_REF] Zhang | Routing in intermittently connected mobile ad hoc networks and delay tolerant networks: overview and challenges[END_REF]). The second is the introduction of node context to reduce the scope of events/subscriptions and use more eciently the network resources. As explained in [START_REF] Royer | Survey on Context-Aware Publish/Subscribe Systems for VANET[END_REF], this context can be assimilated to a location context (an area or just a position). For example, with a classical Publish / Subscribe solution, aircraft can subscribe to "dangerous weather situation". However, this subscription will be too wide and will lead to the reception of useless events since aircraft are interested especially by "dangerous weather situation close to my position or route".

In the past years, some solutions have been proposed by the research community to provide Publish / Subscribe systems adapted to VANET (especially for cars), introducing the persistence of events and the use of node context in event delivery decision process. The aim of this article is to evaluate some solutions using dierent approaches against a typical scenario (notication of weather events) in a realistic AANET based on real trac of commercial aircraft.

The next section presents a brief state of the art on Context-Aware Publish/Subscribe systems adapted for VANET. Section 3 compares some existing solutions based on a typical use case for AANET. The scenario and all implementation choices made to simulate aircraft mobility and Publish / Subscribe systems are explained in this section. Finally, section 4 concludes this article and introduces further work on this topic.

2

Context-Aware Publish/Subscribe solutions in VANET As explained in [START_REF] Royer | Survey on Context-Aware Publish/Subscribe Systems for VANET[END_REF], existing Context-Aware Publish/Subscribe systems that are adapted for VANET can be classied in three categories: geographic routing, proximity routing and overlay network based solutions. Geographic routing based solutions like [START_REF] Leontiadis | Persistent content-based information dissemination in hybrid vehicular networks[END_REF] or [START_REF] Frey | Context-aware publish subscribe in mobile ad hoc networks[END_REF] are based on the assumption that each event / subscription is linked to a geographical area. The publish / subscribe system uses a geographic routing protocol to dispatch events from the publisher to the subscribers. Actually, an event is sent to the nodes located in the geographical area dened as the context of the event. In the same way, a subscription is sent to the geographical area dened in the context. Then, as described in gure 1a, when a subscription matches with an event, nodes located at the intersection can act as middleman between publisher and subscribers.

Proximity routing based solutions like [START_REF] Meier | STEAM: event-based middleware for wireless ad hoc networks[END_REF] or [START_REF] Eugster | Design and Implementation of the Pervaho Middleware for Mobile Context-Aware Applications[END_REF] provide many likeness with the previous solution. The main dierence relates to the area where events and subscriptions are sent. This kind of solutions disseminate events/subscriptions only within a radius close to the sender (see gure 1b). Overlay network based solutions introduce a new component in P/S architecture, the broker. Brokers are responsible for dispatching events from publishers to subscribers. They are interconnected and thus form an overlay network. In most of these solutions, like [START_REF] Cugola | REDS: a recongurable dispatching system[END_REF], brokers form an acyclic graph. Subscriptions are thus sent to all brokers of the graph (or just a part of it if context is used to lter recipient brokers). Based on these subscriptions, they can populate the table used to forward events to interested subscribers. After a brief description of existing context-aware P/S solutions, this section studies the performance of these solutions. Section A describes the scenario used to compare the solutions. This application relies on an AANET. Section B explains all the choices made during the implementation of the scenario and the P/S solution in a simulation model. The last section explains the results obtained for each solutions.

Scenario

The scenario used to compare publish/subscribe solutions is a common problematic for civil aviation: how to inform aircraft of new weather situation ? Indeed, some phenomena like storm or wind shear are very dangerous and thus aircraft must be aware of these phenomena. Nowadays, these information are broadcasted by ground entities (like Air Trac controller or airlines operational center) using radio, VDL (VHF Data Link) or SATCOM means.

A new approach to solve this problem can be found with the use of a context aware P/S system over an AANET composed of commercial aircraft. In this scenario, aircraft aware of a dangerous phenomena, using its weather radar for example, sends an event over the network to announce the new weather situation. This event contains at least the type of event and the estimated duration of the event. The context associated to this event will be the geographical area cover by the weather phenomena. Meanwhile, aircraft know their short term trajectory and then can subscribe to weather events which happen near their route. The context associated to this subscription will be a geographical area that covers a part of the trajectory (see Figure 3). Fig. 3: principle of subscription to weather events 3.2 Simulation / Implementation choices All the simulations have been performed with Omnet++ ( [START_REF] Varga | An overview of the OMNeT++ simulation environment[END_REF]). This paragraph deals with all the implementation choices that has been made to simulate the scenario and the Publish / Subscribe systems with Omnet++.

mobility The mobility model used for the simulations is based on real trac traces from the French civil aviation authority and Eurocontrol based on either radar data or position report from the aircraft itself. Concretely, for a given day of the year, all position reports of commercial ights around the world (in spherical coordinate with an interval of 30min) can be obtained in a le following a specic format.

This le is then parsed by a tool developed at the ENAC, acparse. It oers many options to lter aircraft, interpolate trajectories and export the result to common format. In our case, acparse is used to get an XML le which contains two parts for each aircraft. The rst one (surround by the tag mobility) contains position reports with an interval of 1 minute (thanks to the interpolation that acparse oers) and is used by a new homemade mobility module for Omnet++. The second one (surrounded by the tag plan) can be assimilated to a ight plan. In the simulation model, this ight plan is used by nodes in order to subscribe to weather event close to their short term trajectory.

Due to logistical limitations, it was dicult to simulate all the world trac in reasonable time. As AANET is especially interesting for area where no ground infrastructures are available, as over oceans, we choose to reduce the scope to the transatlantic trac. That is why our simulations are taking into account aircraft whose latitude belongs to the interval [23,5 • -70 • ] and the longitude belongs to the interval [-90

• -10 • ].
Physical/Mac Studies presented in this article deal with Publish / Subscribe service and not with the AANET feasibility which has been demonstrated in other studies like [START_REF] Besse | Wireless Ad Hoc Networks Access For Aeronautical Communications[END_REF]. That's why it is not mandatory to have realistic physical and MAC layer. We choose to use IdealWireless module provided by the inet library [12], which simulates a simple physical and MAC layer. When a node sends a message, all nodes whose distance from the sender is below the parameter transmissionRange receive the message. In the same way, bit rate is a congurable constant. All messages are sent at this bit rate. However the values used for the simulations are deduced from a previous studies evaluating the feasibility of AANET built with commercial aircraft. In particular, [START_REF] Besse | Interference estimation in an aeronautical ad hoc network[END_REF] evaluates the use of CDMA as the layer 2 to interconnect aircraft and obtains an available bit rate of several Mbits/s for a transmission range of several hundred of kilometers.

Events / Subscriptions The generation of events meets the following rules. We x the number of simultaneous events at the beginning of the simulation. We choose to keep this number low due to the size of the simulated AANET (limited to the transatlantic trac). Each event occurs close to a randomly chosen aircraft. These events are then published by the aircraft close to it through the Publish / Subscribe system. When an event ends, a new event is created close to another aircraft and so on until the end of the simulation. The geographical area linked to an event is a cylinder which starts at the ground and nishes at an altitude of 20km (the cruising altitude for a commercial aircraft is around 10km). The cylinder radius is a parameter of the simulation. The generation of subscriptions by aircraft is based on ight plan described in the section 3.2. This ight plan can be broken up to segments. Each segment t 10 minutes of a ight and a subscription is generated for each segment. The time period associated to this subscription is the time interval between two position reports that form the segment. The geographical area linked to the subscription is the segment itself. Consequently, a subscription matches with an event if two conditions are veried. Firstly, the time period of the event ts the subscription's one. Secondly, the segment associated to the subscription intersects or is inside the cylinder which represents the geographical area of the event.

P/S Systems Three Publish / Subscribe systems have been simulated. The rst one is the most simple. Events are disseminated to all nodes of the network. Subscriptions never leave the source node and are recorded in a table. When an event is received by a node, the event is checked against recorded subscriptions. If one matches, the event is forwarded to the application. This solution serves as a reference to compare performance of P/S solutions and will be called "broadcast event solution" in the rest of this article. Intuitively, this solution gives the best result in term of delivery since all events are received by every node but is not resource ecient since even nodes that are not interested by an event receive it anyway.

The second one is a geographic routing based solution. LBM (Location Based Multicast) [START_REF] Ko | Geocasting in mobile ad hoc networks: Locationbased multicast algotithms[END_REF] is the geographical routing protocol used to dispatch events and subscriptions. LBM allows to send a packet to all nodes within a target geographic area (called a multicast region). In this way, the sender of the message denes a forwarding area according to his position and the multicast region. Only nodes within the forwarding area retransmit the message. Thus, the size of the forwarding area has a direct eect on the delivery ratio of the protocol (better with a large region) and the resulting load on the network (lower with a small region). The home-made implementation of LBM oers two options to dene this forwarding region. The rst one is explained in [START_REF] Ko | Geocasting in mobile ad hoc networks: Locationbased multicast algotithms[END_REF] and is based on the distance. A node A forwards the message received from node B for the region R if dist(A, R) < dist(B, R) + δ where dist(A, R) is the distance from A to R and δ is a parameter of the protocol. The second option is a conic approach. In this case, A forwards the message if A ∈ cone(B, R, α)) where cone(B, R, α) is a cone whose the vertex is B, the axis is dened by the points B and the center of the region R, α is the opening angle. This family of Publish / Subscribe solutions also requires a unicast protocol in order to send an event from the node used as an intermediary to subscribers interested by the event. In our case, we choose AODV [START_REF] Perkins | Ad hoc On-Demand Distance Vector (AODV) Routing[END_REF] (Ad hoc On-Demand Distance Vector). This choice was motivated by the nature of the network (hundreds of highly mobile nodes) and the quality of service expected by the application. Our main objective is to be as ecient as possible in order to save bandwidth, furthermore, in this scenario, the delay to disseminate event to all subscribers is not a performance criteria. Thus a reactive routing protocol seems adequate and AODV is the reference in this category.

The third one is a proximity routing based solution. For each event/subscription, a proximity perimeter is dened (parameter of the simulation). The message is sent to all nodes inside this area. As with previous solution, when a node receives an event from which it knows subscribers, AODV is used to send this event to recipients.

Results

For all the results introduced in this section, the table 1 gathers the value of main parameters used for the simulations. The transmission range and bit rate chosen for the simulations are coherent with previous studies on AANET like [START_REF] Vey | Aeronautical Ad Hoc Network for Civil Aviation[END_REF].

For each simulation, two results will be observed: Delivery ratio, for each event Ratio = Number of subscribers that received the event Total number of subscribers Total cost of network, which is the total number of bytes sent in the network (including control message) to dispatch events to subscribers. It is calculated as the sum of sent bytes (at the level of MAC layer) for all nodes of the network. The rst paragraph details some general results, valid for all the simulations and useful to understand other results. The three following paragraphs show and explain results for each simulated pub/sub solution.

General results Table 2 lists the number of events and subscriptions generated during the simulation. We observe there is far more subscriptions message than event. This is the consequence of the chosen scenario. Indeed, only one event is active at a given time, then few events have been generated during the simulation. On the other hand, aircraft generate subscription every 10 minutes. As there are hundreds of aircraft ying for several hours, the result is more than 20000 generated subscriptions.

Number of events 150 Number of subscriptions 23978

Table 2: Number of events and subscriptions Figure 4 shows the number of ying aircraft in the targeted area during the simulated day. This number has a direct impact on the network connectivity and thus on the delivery ratio possible by the pub/sub solution. We see number of aircraft during a part of the night is low compared to rest of the day. Consequently, it is dicult to obtain a delivery ratio close to 100% since it seems obvious that during this period, network partition exists. Previous studies on AANET like [START_REF] Besse | Wireless Ad Hoc Networks Access For Aeronautical Communications[END_REF] found a reachable packet delivery ratio around 90%. It appears to be the maximum possible value of delivery ratio reachable with a publish / subscribe solution. 3 shows the delivery ratio and the network load obtained with this solution.

Delivery ratio 64,02% Network load 3, 70 • 10 7

Table 3: Standard broadcast event solution

The observed result is not as good as estimated. However, it can be easily explained. The event has been disseminated in the network just after its creation. As the life of the event is between 30 and 60 min, it is possible that an aircraft appears after the generation of the event but it is interested by it anyway. With the implemented solution, it can never received the event.

Figure 5 highlights this phenomena. We introduced a new parameter in the simulation: retransmiTime. If its value is strictly higher than 0, it allows the retransmission of the event every retransmiTime seconds. The graphs show the obtained delivery ratio and network load depending on the retransmit time. We see delivery ratio can be increased up to 25 percent if the retransmit time is short (less than 600 seconds). The counter part is the exponential increase of the network load. As expected, this solution can give decent results at the cost of important network load which is not acceptable for AANET where bandwidth resources are very limited. Several reasons can explain this result. Firstly, events and subscriptions are only sent to nodes belonging to the area dened with the geographical context. In our case, it corresponds to the cylinder which represents the weather region for events and a cylinder that surrounds the aircraft trajectory for subscriptions. These regions are small thus few or no aircraft are present to act as middleman between publishers and subscribers. To highlight this problem, the model has been modied to extend geographic context area thanks to a new parameter, proximitySize which denes the radius of cylinder used as context region, the center being the center of the weather situation for events and the center of the segment for subscriptions. Figure 6 shows obtained results according to proxim-itySize's value. We see delivery ratio can be increased if the proximity size is extended but stays below 50% and thus this solution is far less ecient than the broadcast. The high network load can be explained by the characteristic of the application using the Publish / Subscribe service. As seen in section 3.3, aircraft send far more subscriptions than events and this solution forwards subscriptions compared to the previous one. Thus, even if subscriptions' size is lower than events' one, the important number of subscriptions increases substantially the load of the network. 5 gives rst results obtained with the proximity routing based solution (the perimeter has been set to 250km). These results are consistent with those obtained in the previous section. The delivery ratio is lower than with the geographic solution and the network load is roughly equivalent. In the same way than the geographic solution, we verify the inuence of the delivery perimeter on the performance of the solution. Figure 7 shows results according to size of proximity area. As expected, the delivery ratio and the network load increase with the size of the proximity region. The maximum delivery ratio obtained is close to the one obtained with the broadcast event solution. However, the network load is higher. The reason is the same as earlier with the dissemination of subscriptions in the network. In this article, after a brief survey about context aware Publish / Subscribe solutions adapted for VANET, we compare some solutions on AANET based on real commercial trac, especially with a typical application for this type of network. We show that interesting approaches for VANET like geographical routing protocol are not suitable in this situation and lead to worse results than a basic approach based on the broadcast of the event in the network. Proximity based solution gives better results in terms of delivery ratio but with higher load of the network, which is not acceptable for AANET with low bandwidth resources. Further work will concentrate on two axes. Firstly, this article shows that the development of a new publish / subscribe system is necessary to be suited to AANET and specic applications used on it. A promising approach is based on the denition of cluster in the network with dierent forwarding strategy inside and outside cluster groups. For example, we can imagine a solution where subscriptions are broadcasted inside a cluster group and events are sent to all cluster heads. It allows to take into account applications which generate far more subscriptions than events (as the one described in this article). The second axes is the research and the simulation of new application for which content based routing will be interesting. For example, in civil aviation domain, some research team studies ways to build real-time wind map and use of a Publish / Subscribe service over an AANET is a very promising solution to support this application.

Fig. 1 :

 1 Fig. 1: Principle of geographic and proximity routing based P/S

Fig. 2 :

 2 Fig. 2: principle of overlay network based P/S

Fig. 4 :

 4 Fig. 4: of simultaneous aircraft during the day

Fig. 5 :

 5 Fig. 5: Broadcast event solution depending on retransmit time

Fig. 6 :

 6 Fig. 6: Geographic routing based solution

Fig. 7 :

 7 Fig. 7: Proximity routing based solution

Table 1 :

 1 Common parameter values for all simulations

	Duration of the simulation	24h
	Day	2012/12/25
	Transmission range	300 km
	Physical layer bit rate	1 Mb/s
	Duration of an event	uniform(1800s, 3600s)
	Duration of subscriptions	30 min
	Size of zone covered by an event	250km
	Size of a subscription message	32 Bytes
	Size of a event message	4096 Bytes
	Number of simultaneous events	5

Table 4 :

 4 Table 4 gives rst results obtained with this solution. They are worse than expected. The delivery ratio appears lower than with the broadcast solution and most importantly the network load is higher. In conclusion, the broadcast event solution is more ecient than this one. Standard geographic routing based solution

	Delivery ratio 30,66%
	Network load 7, 16 • 10 8

Table 5 :

 5 Standard proximity routing based solution


Fabien Mickaël Royer

Alain Pirovano

AN ENHANCED 1-HOP CLUSTERING ALGORITHM FOR PUBLISH / SUBSCRIBE SYSTEMS IN AANETS

Keywords: Vehicular Ad Hoc Network (VANET), Aeronautical Ad Hoc Network (AANET), Clustering Algorithm, Publish / Subscribe communication systems

Content-based Publish / Subscribe communication paradigm offers a new approach to disseminate messages in the network, where the message content determines the recipients. Many applications used on AANETs, which are a subclass of VANETs, could be more efficient using this paradigm. Many Publish / Subscribe systems suitable for VANETs have been developed, however they are not efficient for some AANET applications. A promising approach is to build a Publish / Subscribe system over a cluster structure to reduce the control overhead and to offer a good scalability. However, the efficiency of this approach strongly depends on the performance of the clustering algorithm. The aim of this article is to propose a new clustering method, named CAPS, which will be the basis for a future content-based Publish / Subscribe system for AANETs. To validate our approach, a simulation model has been developed. Our algorithm has been compared to some other solutions in a modeled AANET context based on real air traffic traces. We show that CAPS gives better results than other solutions in terms of stability while maintaining at a low level the number of cluster groups.

Introduction

The content-Based Publish / Subscribe paradigm offers many advantages over the classical communication paradigm. From the point of view of applications, such a paradigm allows sending messages without using destination addresses. Instead, based on the interests of nodes defined through subscriptions to specific content, the network deduces the recipients of a message from its content and forwards this message as efficiently as possible to all relevant nodes. Publish / Subscribe systems have been initially defined for fixed networks but may also offer advantages for MANETs (Mobile Ad-Hoc Network -selfconfiguring networks of mobile nodes connected by wireless links) and VANETs (Vehicular Ad-Hoc Network, MANETs composed of vehicles). Innovative Applications for VANETs like accident or traffic jam information work naturally over a Publish / Subscribe system. In the same way, AANETs (Aeronautical Ad hoc NETworks which consist in MANETs with aircraft as nodes) can be seen as a subcategory of VANETs [START_REF] Royer | Survey on Context-Aware Publish/Subscribe Systems for VANET[END_REF]. In the context of increasing air traffic and new needs in term of air ground communication systems such as the introduction of the 4D trajectory data link service in a near future [START_REF] Jackson | The 4D trajectory data link (4DTRAD) service -Closing the loop for air traffic control[END_REF], AANETs represent a relevant solution to supplement existing systems like satellite links. This solution does not need a heavy ground infrastructure and can offer an additional bit rate of several hundred of kbits/s for air/ground or air/air communications. Content based Publish / Subscribe application enjoys wide applicability in AANET where resources are limited compared to wired network. Notification of dangerous weather events (e.g.: storms, dust plumes due to volcanic activity) to aircraft whose trajectory comes close to the concerned geographical area is a good example of application. Some Publish / Subscribe systems for VANETs, generally adapted to be aware of the geographical position of the node for the dissemination of events and subscriptions, have already been proposed. [START_REF] Royer | Survey on Context-Aware Publish/Subscribe Systems for VANET[END_REF] establishes a survey of these solutions, categorizing them in three classes: geographic based routing solutions, proximity based routing solutions and overlay based solutions. The first two categories give good results in term of delivery ratio and network load for some specific applications but are too simplistic when events or subscriptions are applicable to wide geographical area. The third category of solutions, based on brokers to disseminate events and subscriptions in the network is the most promising to suit all kind of applications. However, in the context of large and highly dynamics network like AANETs, managing the connectivity and the synchronization between brokers imply that Figure 1. Principle of a Publish / Subscribe system over clusters a large part of the bandwidth has to be dedicated to the overlay management. A solution to mitigate this problem and to offer a better scalability is to build a Publish / Subscribe system over a hierarchical structure such as a clustered structure. In a such a structure, the mobile nodes are organized in a number of virtual groups called clusters. Each cluster is composed of a cluster-head, serving as a local coordinator and in our case a broker for the Publish / Subscribe system, and regular nodes. Cluster-heads record all events sent by nodes belonging to their cluster and the subscriptions are sent to all cluster-heads using a multicast routing protocol as shown in Figure 1. When a cluster-head receives a subscription matching a recorded event, it sends it to the source of the subscription. The efficiency of this kind of Publish / Subscribe system strongly depends on the performance of the clustering method and particularly on the stability of the clusters created. The aim of this article is to propose a new 1-Hop clustering algorithm dubbed CAPS (Cluster Algorithm for Publish / Subscribe) adapted to AANETs in terms of efficiency and stability. The next section presents the related work on this domain. Section details the CAPS algorithm. The validation of CAPS compared to other existing solutions is presented in Section . Finally, Section concludes this article and introduces further work on this topic.

Related work

Clustering algorithms have been proposed for MANET since a number of years. [START_REF] Yu | A survey of clustering schemes for mobile ad hoc networks[END_REF] proposes a survey and a classification for these solutions. A well known algorithm is LCC (Least Cluster Change) [START_REF] Chiang | Routing in clustered multihop, mobile wireless networks with fading channel[END_REF] which is an enhancement of Lowest ID Clustering (LIC) [START_REF] Ephremides | A design concept for reliable mobile radio networks with frequency hopping signaling[END_REF] or Highest Connectivity Clustering (HCC) [START_REF] Ephemrides | A Design Concept for Reliable Mobile Radio Networks with Frequency Hopping Signaling[END_REF]. To reduce the control overheads, LCC dissociates two phases: cluster formation and cluster maintenance. The latter being the key point since the modification of the cluster structure happens only when two cluster-heads move into communication range or when a mobile node looses its clusterhead. Another solution is 3hBAC (3-hop Between Adjacent Cluster-heads) [START_REF] Yu | 3hBAC (3-hop between adjacent clusterheads): a novel non-overlapping clustering algorithm for mobile ad hoc networks[END_REF] which forms a 1-hop non-overlapping cluster structure. To do so, 3hBAC introduces a new role, named cluster-guest which is used to avoid ripple effects in case of re-clustering.

Some other propositions use the mobility as first criteria to build the cluster structure. A good example is MOBIC [START_REF] Basu | A mobility based metric for clustering in mobile ad hoc networks[END_REF] which supposes that nodes with low speed relative to their neighbours are the best candidate to become cluster-heads and thus proposes an aggregate local mobility metric for the cluster head election process. On-Demand WCA (Weighted Clustering Algorithm) [START_REF] Chatterjee | An on-demand weighted clustering algorithm (WCA) for ad hoc networks[END_REF] combines 4 metrics for this election: degree-difference, sum of the distance with all neighbours, average moving speed and clusterhead serving time. Another approach to minimize the control overhead consists in optimizing re-clustering conditions. α-SSCA (Stability Structure Clustering Algorithm) [START_REF] Guizani | Hierarchical cluster-based link state routing protocol for large self-organizing networks[END_REF] proposes to add a condition to merge clusters when two cluster-heads move into communication range of each other. One cluster-head gives up its role only if the stability factor (detailed in Section ) is below α whose value is between 0 and 1.

VANETs being a subclass of MANETs with specific characteristics like constraints trajectory, specific clustering algorithms have been proposed for VANETs. They share the same aim which is to reduce the control overhead and use efficiently the distinctive features of VANETs. [START_REF] Bali | Clustering in vehicular ad hoc networks: taxonomy, challenges and solutions[END_REF] proposes a survey on such techniques, many solutions it details are specific to VANETs composed by car (like MAC based clustering or Driver behavior clustering). However some solutions seem interesting for AANETs like a new stability-based clustering algorithm (SBCA) [START_REF] Ahizoune | A new stability based clustering algorithm (SBCA) for VANETs[END_REF] which aims to reduce communication overhead and to increase the lifetime of the cluster. The election of a cluster-head is based on a merit value that combines the velocity difference between nodes and the number of neighbouring nodes. The specificity of this solution is the election of a secondary clusterhead which serves as a backup in case of problem with the primary cluster-head.

The proposed clustering algorithm

The main idea is to enhance solutions introduced in the previous section using efficiently the additional information available in AANETs nodes, especially its current position and speed. Therefore, as most of the 1-Hop clustering algorithm, this new solution is based on 3 steps: neighbours discovery, thanks to periodic hello messages, a procedure to elect the clusterhead at the beginning and a maintenance procedure.

The key points of this new solution are following.

• Hello messages contain more information, especially the current position and the current speed of the node. • Elections of cluster-head are based on a combined-metric like WCA or SBCA with two metrics: one based on the number of neighbours, the other one based on relative speed. • The maintenance procedure follows the α-SSCA strategy, introducing two conditions before merging clusters. The first one is based on relative speed of cluster-heads and the second one is the same as α-SSCA. To go deeper in the cluster head election process , the first metric uses the number of neighbour to calculate a coefficient D coe f f which is calculated thanks to the formula 1.

D coe f f =              0 if N < 1 N-1 N L thresh -1 if 1 ≤ N ≤ N L thresh 1 if N L thresh < N ≤ N H thresh 1 + N H thresh -N N L thresh ) if N H thresh < N (1)
N is the number of neighbour. N L thresh and N H thresh are parameters of the algorithm. Figure 2 gives a graphical representation of this equation. This coefficient allows the selection of nodes whose the number of neighbour is between N L thresh and N H thresh to act as cluster-head. Hence, a good cluster-head is a node with a minimum number of neighbour to reduce the number of clusters, yet with limited number of neighbours since, in our case, clusterheads have to record all the events of this group and all the active subscriptions of the network. Too many neighbours imply the management of a lots of events and subscriptions and thus some resources or synchronization issue when two clusters have to be merged. The formula 1 translates this consideration into a merit value for the election process. A study has been conducted to optimize the parameters of the algorithm which conclude to

N L thresh = 20 N H thresh = 25
for all the simulations detailed in the next section.

The second metric considers the node velocity. If, for a node n, v n is the vector representing the velocity of the node and N(n) represents the set of n neighbours, then this metric is calculated thanks to the following equation.

S coe f f = ∑ i∈N(n) v n . v i v n v i
It selects preferably nodes which have a low relative velocity with respect to their neighbours. The weight used to elect cluster-heads is a combination of these two metrics.

Score = β D coe f f + (1 -β )S coe f f
β being chosen between 0 and 1. The other important characteristic of this new algorithm is the condition used to determine whether two cluster-heads moving into communication range of each other need to merge their clusters. In such a case, one cluster-head has to give up its role if two conditions are fulfilled. The first one verifies that the two cluster-heads move in the same general direction. The speed vector is the best information to test this condition:

head1. head2 > 0
The second condition is the same as α-SSCA algorithm. If u and v are two neighbour cluster-heads, u gives its role up if

|{X i ∈ Cl(u)/X i ∈ N(v)}| |Cl(u)| >= α
Where N(v) is the set of neighbours of v and Cl(u) represents the set of cluster members of u .

Validation

In order to validate our new clustering algorithm, a dedicated simulation model has been developed on OMNET++ ( [START_REF] Varga | An overview of the OMNeT++ simulation environment[END_REF]). We also implemented three others 1-Hop clustering algorithms to assess the performances of our new solution: LCC-HCC, α-SSCA (with α equals to 0.5) and MOBIC.

Scenario

These 4 clustering algorithms have been compared in an modeled AANET context based on real air traffic traces taken from Eurocontrol. As AANETs are especially interesting in areas where no ground infrastructures are available, as over oceans, we choose to focus on the north transatlantic traffic. Therefore, our simulations are taking into account aircraft whose latitude are in the [23,5 • N -70 • N] interval and the longitude belongs to the [90 • W -10 • E] interval. Furthermore, in order to evaluate the solutions in different situations, the traffic of several days has been used as input for the simulations with similar results. Since most of the north transatlantic traffic is composed by regular flights, any two days close on the calendar are very similar in terms of number of aircraft and trajectories. The difference lies mainly in the period of the year. For example, on average, there is much more traffic during the summer. That's why, in this article, we choose to present the results for two days. The first one (2012/12/25) is a day where the commercial traffic is very low as opposed to the second one (2012/10/05). Table I summarizes the configuration used for the simulations. The layer 2 protocol has been simulated using the IdealWireless module form INET1 with fixed transmission range and bit rate. However the values used for the simulations are deduced from a previous study evaluating the feasibility of AANETs with commercial aircraft over the north transatlantic corridor. In particular, [START_REF] Besse | Interference estimation in an aeronautical ad hoc network[END_REF] evaluates the use of CDMA in the layer 2 to interconnect aircraft. Depending on the configuration such as the used bandwidth or modulation, the obtained link capacity between aircraft varies from hundreds of Kbits/s to two Mbits/s for a transmission range of several hundred of kilometers. The choice of one Mbits/s is a good trade off for our simulations.

To evaluate the efficiency of the clustering algorithms, three parameters are measures during the simulations:

1) The number of cluster-head at a given time.

2) The total number of cluster-head changes, a lower value indicates a better stability of the clustering algorithm.

3) The lifetime of cluster-heads which is also an indicator of the cluster stability.

Results

Figures 3 and4 show all the results obtained from the simulations. To facilitate the analysis, the cluster head duration is plotted as a cumulative distribution probability. In the same way, the first graph represents the number of aircraft during the considered day to easily make a correlation between this number and the number of cluster head. As expected, most of the traffic in the North Atlantic Airspace is organized in two major flows (peaks on the graph): a westbound flow departing Europe in the morning and an eastbound flow departing North America in the evening (UTC time). The results show that the clustering algorithms taking into account the mobility of nodes offer better stability since the cluster-head duration is higher. The counterpart is a higher number of cluster head. In detail, CAPS gives the best performance compared to others solutions in terms of stability. Table II shows that the mean lifetime of clusters for CAPS is 10% higher than for MOBIC while the number of clusterheads stays close to the MOBIC one. CAPS seems to be a solid basis for a Publish /Subscribe system adapted to AANETs since it will offer a very stable overlay.

Conclusion

In this article, after emphasizing the interest of a cluster structure as an overlay of an efficient VANET adapted Publish / Subscribe systems, a brief Our further work will concentrate on the OM-NET++ implementation of this Publish / Subscribe system and its evaluation in comparison to other solutions like proximity or geographic routing based Publish / Subscribe.

Figure 2 .

 2 Figure 2. Calculation of D coe f f (example with N L thresh = 10 and N H thresh = 15)

Figure 3 .Figure 4 .

 34 Figure 3. Results for the December 25 2012

Table I .

 I Simulations parameters

	Duration of the simulation	24h
	Day	2012/12/25 and 2012/10/05
	Transmission range	300 km
	Physical layer bitrate	1 Mb/s
	Hello interval	1 second

Table II .

 II Detailed results for the 2012/12/25 algorithms adapted to MANETs is proposed, focused on low maintenance algorithm. While interesting, these solutions don't use all the information available to the nodes in an AANET. Particularly, to exploit the velocity of aircraft in addition to the number of neighbours, we created a new clustering algorithm, and then compared it to other solutions through simulations on OMNET++ with a realistic AANET mobility model. Simulation results show that this new algorithm provides better performance than other simulated solutions in terms of stability. It minimizes the number of cluster-head changes and increases the lifetime of clusters. This new algorithm represents a solid basis for a hierarchical Publish / Subscribe systems adapted to AANET.

		CAPS MOBIC
	Number of cluster-heads (mean)	80.7	78.9
	Cluster-head duration (mean)	4775s	4152s
	Number of cluster-head changes	1981	2298
	survey on clustering

http://inet.omnetpp.org/