Zahra Ismaïl

Olivier Jamet

Keywords: Geoid, Remove-Compute-Restore, Stokes' integration, Integration radius, Gravity anomaly grid resolution, EGM2008, GRAVSOFT

Geoid determinations by the Remove-Compute-Restore (R-C-R) technique involve the application of Stokes' integral on reduced gravity anomalies. Numerical Stokes' integration produces an error depending on the choice of the integration radius, grid resolution and Stokes' kernel function.

In this work, we aim to evaluate the accuracy of Stokes' integral through a study on synthetic gravitational signals derived from EGM2008 on three different landscape areas with respect to the size of the integration domain and the resolution of the anomaly grid. The influence of the integration radius was studied earlier by several authors. Using real data, they found that the choice of relatively small radii (less than 1 • ) enables to reach an optimal accuracy. We observe a general behaviour coherent with these earlier studies. On the other hand, we notice that increasing the integration radius up to 2 • or 2.5 • might bring significantly better results. We also note that, unlike the smallest radius corresponding to a local minimum of the error curve, the optimal radius in the range 0 • to 6 • depends on the terrain characteristics. We find also that the high frequencies, from degree 600, improve continuously with the integration radius in both semi-mountainous and mountain areas.

Finally, we note that the relative error of the computed geoid heights depends weakly on the anomaly spherical harmonic degree in the range from degree 200 to 2000. It remains greater than 10% for any integration radii up to 6 • . This result tends to prove that a one centimetre accuracy cannot be reached in semi-mountainous and mountainous regions with the unmodified Stokes' kernel.

Accuracy of unmodified Stokes' integration in the R-C-R procedure for geoid computation

Introduction

Different approaches for regional geoid and quasi-geoid determination have been proposed, such as least-squares collocation [START_REF] Krarup | A Contribution to the Mathematical Foundation of Physical Geodesy[END_REF][START_REF] Yildiz | Comparison of Remove-Compute-Restore and Least Squares Modification of Stokes' Formula Techniques to Quasi-geoid Determination over the Auvergne Test area[END_REF], wavelet modeling [START_REF] Panet | Wavelet Modelling of the Gravity Field by Domain Decomposition Methods: An Example Over Japan[END_REF] and the remove-compute-restore (R-C-R) techniques.

In an operational context, one widely used method is the R-C-R procedure along with the Residual Terrain Model (RTM) reduction [START_REF] Duquenne | A New Solution for the Quasigeoid in France QGF98[END_REF][START_REF] Kilicoglu | Regional Gravimetric Quasi-Geoid Model and Transformation Surface to National Height System for Turkey (THG-09)[END_REF][START_REF] Srinivas | Gravimetric Geoid of a Part of South India and its Comparison with Global Geopotential Models and GPS-Levelling Data[END_REF][START_REF] Yildiz | Comparison of Remove-Compute-Restore and Least Squares Modification of Stokes' Formula Techniques to Quasi-geoid Determination over the Auvergne Test area[END_REF][START_REF] Valty | Quasi-Geoid of New Caledonia: Computation, Results and Analysis[END_REF]. This method was introduced by Forsberg [START_REF] Forsberg | Modelling the Fine Structure of the Geoid: Methods, Data Requirements and Some Results[END_REF] in order to combine rigorously the information of three frequency ranges of the gravity field -the short, intermediate and long wavelength signals -through three steps, including a Stokes' integration, to obtain a geoid model.

The precision of the Stokes' integration phase has been studied for long, in particular to investigate the effects of the size of the integration domain, of the computation method and of the improvement brought by modifying the integration kernel itself. Kearsley (1986) assessed the so-called ring integration method (discrete approximation of the integral) on free-air anomalies by comparing integration results to geoid heights measured as the difference between GPS ellipsoidal heights and levelling on a test network in Ohio, USA [START_REF] Kearsley | The Determination of Precise Geoid Height Differences Using Ring Integration[END_REF][START_REF] Kearsley | Data Requirements for Determining Precise Relative Geoid Heights from Gravimetry[END_REF]. He found that increasing the cap size degraded the evaluation of N, and concluded that the integration radius to compute the contribution of gravimetric data should be in the order of 0.5 • for the degree and order of the reference GGM used.

A more detailed study by the same author investigated, within a R-C-R procedure using the OSU81 geopotential model, the relation between the maximum degree n max of the geopotential model and the optimal integration cap size over two test areas in Canada [START_REF] Kearsley | Tests on the Recovery of Precise Geoid Height Differences from Gravimetry[END_REF]. A strong correlation between 180/n max and the optimal cap radius was found, with an optimal cap radius of about 1 • for n max = 180. He also observed that the choice of the optimum radius was not critical around this value.

In 1998, Higgins et al. compared ring integration to Fourier-based computation of the Stokes' integral in a R-C-R procedure on a data set in Queensland, Australia. They used the EGM96 geopotential model for the removal of low frequencies [START_REF] Higgins | The Effects of Varying Cap Sizes on Geoid Computations: Experiences with FFT and Ring Integration[END_REF]. They tested for cap sizes up to 1 • and found that the optimum was reached for values ranging between 0.2 • and 0.5 • .

However, some other works, including recent ones, choose larger values of this parameter, ranging from 1 • to 3 • depending on the terrain characteristics [START_REF] Ellmann | The Geoid for the Baltic Countries Determined by the Least Squares Modification of Stokes Formula[END_REF][START_REF] Hwang | Terrain Correction Computation Using Gaussian Quadrature[END_REF][START_REF] Valty | Auvergne Dataset: Testing Several Geoid Computation Methods[END_REF].

In this paper, we address the question of the Stokes' integral accuracy from a complementary point of view. Recent satellite missions dedicated to the observation of the Earth's gravity field considerably improved the accuracy and resolution of geopotential models, in particular since the publication of the EGM2008 model complete to degree and order 2159 [START_REF] Nikolaos | The Development and Evaluation of the Earth Gravitational Model 2008 (EGM2008)[END_REF]. This new data gives us the ability to evaluate the computational techniques on synthetic gravimetric data with a reasonably realistic spectral content. The use of synthetic data provides indeed a reference to control the accuracy of the computation process independently of the gravity data and interpolation quality. A similar approach was adopted for instance by Featherstone [START_REF] Featherstone | Tests of Two Forms of Stokes's Integral Using a Synthetic Gravity Field Based on Spherical Harmonics[END_REF], using the EGM96 geopotential model, in order to compare the effect of several Stokes' kernel modifications. However, to our knowledge, synthetic experiments were not performed to revisit the question of the Stokes' integration parameters.

The experiments presented here thus intend to study the dependency of the precision of the Stokes' integration on the gravimetric anomalies grid resolution and on the integration cap size, in the ideal situation when the gravimetric data and the reference data are free of any error. They also intend to assess the dependency of this precision to the spectral content of the gravity signal. Gridded gravimetric data as well as reference data are derived from EGM2008 taken as the reference gravitational field. The Stokes' integration study is limited to the unmodified kernel and to a common discretized space-domain integration.

Background

The geoid is defined as the equipotential surface of the Earth's gravity field corresponding to the mean sea surface. High-resolution geoid models can be used in numerous applications such as geodesy, geophysics, oceanography, etc. In the last decades, improving the quality of regional geoid models toward a 1 cm accuracy has become one of the major goals of various research groups and mapping agencies [START_REF] Esan | Spectral Analysis of Gravity Field Data and Errors in View of Sub-Decimeter Geoid Determination in Canada[END_REF][START_REF] Featherstone | The AUSGeoid09 model of the Australian Height Datum[END_REF][START_REF] Krynski | Centimetre Quasigeoid Modelling in Poland using Heterogeneous Data[END_REF].

The R-C-R method with RTM reduction [START_REF] Forsberg | Modelling the Fine Structure of the Geoid: Methods, Data Requirements and Some Results[END_REF] combines rigorously the information of three frequency ranges of the gravity field -the short, intermediate and long wavelength signals -through three steps to obtain a geoid model.

1. First, one computes the residual gravity anomalies ( g res ) -on the geoid in case of geoid computation and on the topography in case of quasi-geoid computation -defined as :

g res = g -g M -g RT ( 1
)
where g is the free-air gravity anomaly, g M is the long wavelength part computed from a global gravity model, g RT is the short wavelength part computed as the attraction of the residual topography. g res anomalies are interpolated into a regular grid.

2. The next step is the integration, where the N res geoid or ζ res quasi-geoid residual heights are derived from the residual gravity anomalies through the computation of Stokes' integral. The original Stokes' formula for gravimetric geoid model is given by Stokes [START_REF] Stokes | On the Variation of Gravity at the Surface of the Earth[END_REF]:

N res = R 4πγ σ S(Ψ)∆g res dσ ( 2
)
where R is the mean Earth radius,γ is the normal gravity on the reference ellipsoid, dσ is the integration surface element of the unit sphere and S(Ψ) is Stokes' kernel as function of the spherical distance Ψ between the computation and data point [START_REF] Heiskanen | Physical geodesy[END_REF].

S(Ψ) = 1 sin( Ψ 2 ) -6sin Ψ 2 + 1 -5cosΨ -3cos(Ψ)ln sin Ψ 2 + sin 2 Ψ 2 (3)
3. The final step is the "restore" phase, where the long and short wavelengths, which have been removed earlier from the gravity signal, are restored to compute the geoid or quasi-geoid heights.

N = N M + N RT + N res (4)
In step 2, Stokes' formula assumes an integration over the whole sphere. For practical considerations, the domain of computation is usually limited to a reduced area around the computation point, and computed on discretized data. Thus, Equation 2 is approximated and evaluated as a limited summation, for Ψ ≤ Ψ 0 , of products of Stokes' kernel and gravity anomalies over a grid, where Ψ 0 is called the integration radius or the integration cap size according to the authors.

Several authors have proposed to modify the Stokes' kernel in order to reduce the truncation error of the remote zone, either in a deterministic way [START_REF] Ellmann | The Geoid for the Baltic Countries Determined by the Least Squares Modification of Stokes Formula[END_REF] or though a stochastic approach [START_REF] Sjoberg | A General Model for Modifying Stokes's Formula and its Least-Squares Solution[END_REF], but the unmodified kernel remains in wide use.

In the following, we study the effect of the grid resolution and of the cap size Ψ 0 on the precision of the integration. In section 3, we present the methodology and data used. In section 4 we discuss the experimental results.

3 Methodology and data

General principles

In order to evaluate the precision of the sole integration step, we use the global gravity field model EGM2008 [START_REF] Nikolaos | The Development and Evaluation of the Earth Gravitational Model 2008 (EGM2008)[END_REF] to generate a set of synthetic data composed of gridded gravity anomalies and reference geoid heights on selected test areas. The low frequencies of EGM2008, up to degree n max , play the part of the geopotential model used in a R-C-R method. The higher frequencies, from degree n max to 2000, play the part of the residual gravity signal.

We consider the case of a quasi-geoid computation. Free-air anomalies are computed on the topographical surface given by a digital terrain model. Since the geopotential model we used has a limited resolution of about10 km, and since we can generate gridded anomalies at any resolution, we did not include, in our test, a removal of the high frequencies of the gravimetric signal. An error analysis is performed on a set of points chosen within the test areas. In order to observe the effect of the characteristics of the terrain, test areas are of limited size (1

• × 1 • ).

Test data synthesis and evaluation

Synthesis

More precisely, let us note W M and g M respectively as the gravitational potential and acceleration computed from the degrees 0 to n max of EGM2008 (the geopotential model in our experiment), W and g the total potential and acceleration of the considered gravitationnal field given by EGM2008 up to degree and order 2000, and γ the normal acceleration.

For a given

1 • × 1 • test area, let {P i (ϕ i , λ i , H i )} i=1.
.p be a set of p points, where (ϕ i , λ i ) are the geographical coordinates of the point P i and H i its normal altitude. Let Q i be the corresponding point on the telluroid.

The quasi-geoid height ζ i at point P i is decomposed into

ζ i = ζ M,i + ζ res,i (5)
where ζ M,i is the contribution of the geopotential W M and ζ res,i is the residual quasi-geoid height. At point P i , the reference quasi-geoid residual height will be given by the Brun's formula [START_REF] Heiskanen | Physical geodesy[END_REF]

ζ ref res,i = W (Q i ) -W M (Q i ) γ (Q i ) (6)
The estimated quasi-geoid height ζ comp res,i is obtained by Stokes' integration of gridded residual anomalies ∆g res , computed at each grid point, on the topographical surface given by the DTM, by the spherical approximation of the fundamental equation of geodesy [START_REF] Heiskanen | Physical geodesy[END_REF]

∆g res = - ∂ (W -W M ) ∂r -2 W -W M R ( 7
)
where R is the average radius of the Earth.

The geographical coordinates of the points P i are randomly chosen within the test area, while their normal altitudes are obtained by a bilinear interpolation in a digital terrain model (DTM) of the area.

The anomaly grids were computed up to 6.5 • from the center of each test area. These data sets allow us to test integration radii up to 6 • . Although results may be sensitive to larger values of this parameter, 6 • is already far beyond the possible values in practice, either because of the computational effort, or because of the lack of gravity data far from the area of interest.

Evaluation

The evaluation of the Stokes' integration results is performed though the comparison between the estimated quasi-geoid heights ζ comp res,i and their reference values ζ ref res,i . Knowing that Stokes' integration over a limited spherical cap is likely to be biased -since the convolution kernel has no longer a zero mean -, we use the same criterion as Higgins [START_REF] Higgins | The Effects of Varying Cap Sizes on Geoid Computations: Experiences with FFT and Ring Integration[END_REF] by computing the root mean square error of the differences of quasi-geoid heights along the baselines between every pair of points of the area. For each test area, the precision of the Stokes' integration is thus evaluated as

= 2 p (p -1) p-1 i=1 p j=i+1 ∆H comp i,j -∆H ref i,j 2 (8)
where, for computed and reference values

∆H i,j = ζ res,i -ζ res,j (9)
Due to the sizes of the zones tested, the baselines taken into account in the computation of the error range from a few kilometers to about 100 km. We did not separate them into several classes of distance as done in [START_REF] Higgins | The Effects of Varying Cap Sizes on Geoid Computations: Experiences with FFT and Ring Integration[END_REF], considering that, in the latter study, baselines lesser than 50 km and lesser than 100 km showed a very similar behaviour.

We found that this criterion led to the same conclusions as an evaluation by the standard deviation σ of the difference ζ comp res,i -ζ ref res,i , and will also use the standard deviation in the following.

Data and software

As mentioned before, synthetic data are generated using the EGM2008 model [START_REF] Nikolaos | The Development and Evaluation of the Earth Gravitational Model 2008 (EGM2008)[END_REF]. This model can be used for different tasks concerning geoid computation due to its relatively high frequencies and its good fit to the gravity field over different areas [START_REF] Eshagh | A Strategy Towards an EGM08-based Fennoscandian Geoid Model[END_REF][START_REF] Nikolaos | The Development and Evaluation of the Earth Gravitational Model 2008 (EGM2008)[END_REF]. We used the digital model ETOPO5 [START_REF]Digital Relief of the Surface of the Earth[END_REF] as a model of the topographical surface considered conventionally as giving the normal heights over the processed areas. We considered that the resolution of 5 arc minutes of the DTM was sufficient with respect to the resolution of the gravity data used. Stokes' integrals were computed with the stokes function of the GRAVSOFT package [START_REF] Forsberg | Geodetic Gravity Field Modelling Programs[END_REF]. Synthetic data are produced by software developed by the authors.

Test areas

Three different areas were selected in France: a plain, a semi-mountainous and a mountainous area (Figure 1), in order to investigate the effect of the landscape on the choice of the parameters, especially in the case of the mountainous regions. The accuracy of one centimetre is known to be hard to reach in mountainous areas because of the insufficient gravity data coverage and strong topography signal at short wavelengths [START_REF] Flury | Short-Wavelength Spectral Properties of the Gravity Field from a Range of Regional Data Sets[END_REF][START_REF] Yildiz | Comparison of Remove-Compute-Restore and Least Squares Modification of Stokes' Formula Techniques to Quasi-geoid Determination over the Auvergne Test area[END_REF]. Topography and gravity anomalies are illustrated in Figure 2 and Figure 3. Statistical information on the test areas are displayed in Table 1 andTable 2. 4 Results

Grid resolution

Since we use the EGM2008 up to degree 2000, a resolution of about 0.1 • should be sufficient to represent the content of the gravimetric signal to perform the integration. Therefore, it is expected that the grid resolution will not have much influence over the precision if it is smaller than 0.1 • . This conclusion is confirmed by the experiments. We tested resolutions from 0.005 • to 0.5 • at a fixed integration radius of 2 • and for residual anomalies in several frequency ranges. The value of 2 • was chosen in accordance with the results of other experiments (see section 4.2). We decomposed the residual anomalies into frequency ranges in order to test sparse grids fulfilling the requirements of the signal theory.

Figure 4 shows the evolution of the standard deviation of Stoke's integral errors as a function of the grid resolution, averaged on the 3 test areas, and for spherical harmonic degrees from 200 to 2000. Each curve corresponds to a given frequency range and is drawn in the resolution interval allowed by the sampling theory.

For all the frequency ranges, one can observe that choosing a resolution fulfilling strictly the sampling theory (the right extremity of each curve in Figure 4 is not sufficient to ensure an optimal precision. The precision loss at raw grid resolution is likely to come from the numerical accuracy of the integration algorithm. As expected, after a small precision improvement when going toward finer resolutions, the error reaches its minimal value.

All degrees being considered, the degradation of the resolution from 0.005 • to 0.075 • does not significantly alter the precision. Experiments with other integration radii confirmed these results.

Integration radius

Findings

Figure 5 shows the error defined by Equation 8 as a function of the integration radius Ψ 0 for n max = 199, (a) for Ψ 0 varying from 0 to 6 • , and (b) with a focus on the range 0 to 3 • . Our choice to focus the study on the integration of the 200-2000 spherical harmonic degree range was driven by the RCR-RTM context: due to the limits of spatial gravity models and of the residual terrain effect removal, we suppose that the residual anomalies will keep a significant energy in these frequencies. Ψ 0 was sampled at resolution 0.1 • , with an oversampling at resolution 0.01 • in the neighbourhood of the first local minimum of the curves.

When averaged on the 3 test areas (Figure 5, continuous curve) the evolution of the error with respect to Ψ 0 shows a general behaviour coherent with previous studies. Kearsley, for instance, already noticed the steep decreasing of the error for small values of Ψ 0 , followed by large oscillations [START_REF] Kearsley | Tests on the Recovery of Precise Geoid Height Differences from Gravimetry[END_REF]. The first minimum of the error occurs for Ψ 0 = 0.24 • . The absolute minimum of the error is reached for Ψ 0 = 1.9 • , but with a value of the error smaller than the first minimum by only 7%.

However, this behaviour shows a significant variability according to the area. In the flatter region (Figure 5, Zone 1, dotted curve with the symbol ×), the difference between the minimal error ( = 0.022 at Ψ 0 = 2 • ) and the first minimum ( = 0.025 at Ψ 0 = 0.25 • ) is weak. In the semi-mountainous area (Figure 5, Zone 2, dashed line with the symbol +), the curve shows a much poorer precision at Ψ 0 = 0.25 • and a real gain to process the integral with an integration radius of 2 • (error nearly divided by 3). In mountain (Figure 5, Zone 3, finely dotted line with the symbol ), the first minimum is the best choice if one cannot use integration radii larger than 1.5 • . However, a much better minimum of is reached for Ψ 0 = 2.6 • , with an error decrease of more than 40%. Another fact that must be underlined is that, for all the 3 areas, the error of the integration appears to be very sensitive to the choice of Ψ 0 around the value yielding the best results. Such a behaviour was observed with geopotential models of low degree (up to degree 90) [START_REF] Kearsley | Tests on the Recovery of Precise Geoid Height Differences from Gravimetry[END_REF], but not at all with models of higher degree.

These findings have some significance regarding operational applications. The precision of the 1 st order of the French levelling network is estimated to be 2 mm/ √ km, meaning that one will expect the rms value of the height differences to be in the order of 2 cm for baselines of 100 km. This value is barely achieved by the best minima on the 3 zones.

Interpretation

Several facts may explain the differences observed with previously published results. The main ones probably concern the data used. Our experiments rely on error-free sets of data. Previous works used real data, with the consequence that both gravimetric data and reference data may be biased. We do not believe that the errors of the reference data play a great part: Kearsley uses for instance an evaluation criterion that eliminates the effects of biases and drifts in the levelling, and more importantly, systematic errors in the GPS heighting [START_REF] Kearsley | Tests on the Recovery of Precise Geoid Height Differences from Gravimetry[END_REF] . On the contrary, the quality of the gridded anomalies is likely to play a major part. Previous studies rely on gravity measurements with a spacing in the order of 10 km, and an interpolation of this data over a grid -after the removal of the residual terrain contribution for some works.

However, other factors may have to be taken into consideration. First, in our synthesis of the gravimetric anomalies, we used the spherical approximation of the fundamental equation of geodesy (Equation 7), for which the Stokes' integration is the theoretical solution. With real measurements, anomalies are computed as the difference between the module of the measured acceleration on the topography and the module of the normal acceleration on the telluroid. We verified that the synthetic anomalies from Equation 7were not biased with respect to anomalies computed from the acceleration modules (rms 0.06 mGal on a 7 • × 7 • area around our second test area). We thus do not think that this choice influences our results.

Secondly, the results in Figure 5 correspond to n max = 199, while other works used different values (180 in [START_REF] Kearsley | Tests on the Recovery of Precise Geoid Height Differences from Gravimetry[END_REF], and 360 in [START_REF] Higgins | The Effects of Varying Cap Sizes on Geoid Computations: Experiences with FFT and Ring Integration[END_REF]). The dependency of the optimal radius to the degree of the geopotential model has been well documented as mentioned in section 1. However, this dependency was supposed to be more or less independent of the landscape, and to be characterised by smaller optimal integration radii as the degree of the geopotential model increase. Figure 6 shows the evolution of the error as function of the integration radius, for the semi-mountainous area and for a few values of n max . The dependency of the location of the first minimum of the curve conforms to the one observed in previous works. One also observes that, as n max increases, the oscillations of the curves become smaller, and nearly fade out for n max = 360 (Figure 5, continuous curve with the symbol ). However, even for this latter value, the integration over a larger domain, up to Ψ 0 2 • , brings some precision improvement.

Finally, the accuracy of the geopotential model itself may be a cause. In our tests, the geopotential model is chosen as the low frequencies of the model used for synthesising the gravity data, which corresponds to an unrealistic situation of perfect removal of the low frequencies. In order to evaluate the possible induced bias, we compared the use of the geopotential models OSU81 [START_REF] Rapp | The Earth's Gravity Field to Degree and Order 180 Using Seasat Altimeter Data, Terrestrial Ggravity Data and other Data[END_REF] and EGM96 [START_REF] Lemoine | The Development of the Joint NASA GSFC and NIMA Geopotential Model EGM96[END_REF] at degree and order 180 in place of the low degrees of EGM2008 for the values of W M in the computation of the residual anomalies (Equation 7). Figure 7 shows this comparison for the three areas. In the plains, and with the OSU81 model (Figure 7, left, dotted curve with the symbol ×), we find a behaviour similar to the one observed by Kearsley [START_REF] Kearsley | Tests on the Recovery of Precise Geoid Height Differences from Gravimetry[END_REF]. The first minimum is reached for Ψ 0 0.3. Then the error oscillates and grows with the integration radius. However, with EMG96 and EGM2008 (dashed curve with the symbol + and continuous curve with the symbol ), while one observes the same oscillations and location of the minima, the error decreases with the radius. This example is interesting since Kearsley's data correspond also to a rather flat area. It tends to show that the improvement of the geopotential models has changed the conditions of use of Stokes' integral in the plain. On the other test areas (Figure 7, center and right), even when using OSU81, there is a significant gain in choosing a larger integration radius, especially in the semi-mountainous area (center). The quality of the geopotential model seems less of a determining factor in those cases.

Dependency on the spectral content

While the previous experiments bring some information about the behaviour of Stokes' integration, the very error values observed may not be significant. Indeed, the spectral content of the signals we processed corresponds to the full gravity signal in the corresponding spectral bands. In real contexts, this spectral content will vary according to the quality of the geopotential model in the locality of the test areas, and according to the choice made for the removal of the residual terrain effect.

In order to investigate whether general conclusion could be drawn from this study as far as the absolute error value is concerned, we computed the integration on separated spectral bands -defined as intervals of spherical harmonics degrees. Figure 8 shows the dependency of the standard deviation of ζ comp res -ζ ref res as function of Ψ 0 to these spectral bands. In the plains (Zone 1), the error is clearly dominated by lowest degrees (200 -359), with a minimum for Ψ 0 2 • . In the semi-mountainous area (Zone 2), increasing the integration radius from the first minimum does not improve the quality for the lowest degrees (200 -359). However, it does improve the results for intermediate degrees (200 -599; 200 -799). In the mountains, the behaviour in much more erratic, but, for Ψ 0 ≥ 1 • , the error is also dominated by the degrees up to 599.

These later spectral ranges are also obviously the ones that model the major part of the energy of the signal, and would be largely modified by a terrain effect removal. The Figure 9 presents the relative precision ε rel of the Stoke's integration values, defined as the ratio between the standard deviation of the error and the standard deviation of the reference signal computed on the test points:

ε rel = σ ζ comp res,i -ζ ref res,i σ ζ ref res,i (10)
where the standard deviation σ (v i ) of a quantity v i at the p test points P i is defined by

σ (v i ) = 1 p p i=1 v 2 i - 1 p p i=1 v i 2
One observes that, except for the plain test area, the quality of the estimate of the high frequencies (from degree 600) of the geoid improves continuously with the integration radius. The exception of the plain is not surprising since the energy of the corresponding signal in this area is very low.

Secondly, we note that, from degree 600, the tested spherical harmonic frequency ranges show a similar behaviour, with a relative error keeping between 10% and 20%-25% for all integration radii greater than 1 • , with the exception of the degree range 600-799 of the semi-mountainous area. Moreover, disregarding this exception, 10% appears to give a lower bound of the relative error for all the frequency ranges. This result tends to show that independently of the spectral content of the gravity anomalies or of the geoid, for integration radii lesser than 6 • , the precision of the Stoke's integration cannot be better than 10% of the standard deviation of the computed residual geoid.

This gives a basis for a rough estimate of the best achievable accuracy in real cases. For instance, Duquenne [START_REF] Duquenne | A Data Set to Test Geoid Computation Methods[END_REF] provided, in the Auvergne data set, a geoid model covering our second test area, and based on the R-C-R technique using a global gravity field up to 360 • . The standard deviation of residual geoid height for this geoid model is about 11 cm. Based on our conclusion, the best accuracy we can expect from a standard Stokes' integration is cm, larger values being in fact expected due the other error sources.

Conclusion and Recommendations

In order to evaluate the accuracy of integration step in R-C-R procedure, we generated a set of synthetic data based on EGM2008 to degree and order 2000. We studied the influence of the integration radius and of the grid resolution using a standard Stokes' kernel, through an analysis of the accuracy of the computed height differences and of the behaviour the results in different frequency ranges.

As expected, the grid resolution does not affect the results providing that it is kept small enough, as suggested by the sampling theory. Considering the overall errors, the choice of very large integration radii, up to 6 • , seems to lack any benefit. Beyond 2 • to 2.5 • , the precision does not improve significantly. However, the relative precision of the high frequencies (degree > 600) seems to improve gradually with the integration radius for the whole studied interval, suggesting that the local variations of the geoid be potentially better represented with large integration radii. This question would deserve a deeper study dedicated to the evaluation of the accuracy of the height differences along small baselines.

These results partly contradict previous studies led on the basis of real data. First, the choice of relatively small radii (lesser than 1 • ) does not allow to reach the optimal accuracy on synthetic data. A local minimum of the error is always observed for radii around 0.5 • as in previous works. Its very location depends mainly on the degree of the geopotential model as noticed by other authors. Nevertheless, this first local minimum might yield far from optimal results. The choice of the best radius and the gain it brings with respect to the first minimum seems to be highly dependent on the characteristics of the terrain. In addition, the quality of the results appears to be very sensitive to the chosen radius around the optimal location. These findings do not seem to depend significantly on the quality of the geopotential model used to remove the low frequencies of the gravity signal.

These differences with previous works are very likely to come from the lower quality of the gravity anomalies available in real cases. However, this tends to indicate that, with the improvement of the gravity surveying, one has to choose the Stokes' integration radius with caution, possibly by checking the process results on control data.

Lastly, our results tend to show that, independently of the landscape and from the spectral content of the residual gravity anomalies, the unmodified Stokes' kernel cannot reach a precision better than 10% of the geoid signal to be retrieved. The lower bound of the accuracy is probably even higher since our study does not take into account other error sources, such as the interpolation errors, the spherical approximation errors, the aliasing error, etc. This work, however, poses new questions. The experiments presented here apply the Stokes' integral to anomalies computed on the topography. Our synthetic data do not in that sense fulfil the strict requirements of the Stokes' theory but should be processed in the frame of Molodensky's theory [START_REF] Moritz | Advanced Physical Geodesy. Sammlung Wichmann: Neue Folge[END_REF], or at least take into account the effect of the flattening of the ellipsoid.

Moreover, this study was limited to the evaluation of the unmodified Stokes' kernel and on a restricted number of test areas, regardless of the accuracy of the interpolation of the observed gravity anomalies. A similar methodology could be used to evaluate other integration kernels and interpolation methods.

Graphics

Figure 1 : 13 Figure 2 :

 1132 Figure 1: Location of test zones in France

Figure 3 :

 3 Figure 3: Residual gravity anomaly of the test zones in France for wavelength bands between degree 600 and 800. Radius of Stokes' integration= 2 • . Gravity anomaly grid resolution = 0.050 • . The black dots represent the test points in each area.

Figure 4 :

 4 Figure 4: Accuracy of the geoid heights in meters computed using the standard Stokes' kernel as a function of the gravity anomaly grid resolution averaged on the three test areas. Integration radius=2 • . Colors correspond to different frequency ranges(see legend)

Figure 5 :

 5 Figure 5: Error in m (Equation 8) as function of the integration radius in angular degrees for a geopotential model at degree 199. Left (a): evolution up to 6 • ; right (b): zoom of (a) on small integration radii. Curves correspond to different test areas (see legend)

Figure 6 :

 6 Figure 6: Error in m (Equation 8) as function of the integration radius in angular degrees for the semi-mountainous test area (Zone 2). Curves correspond to different values of maximum degree n max of the EGM2008 geopotential model (see legend)

Figure 7 :

 7 Figure 7: Error in m (Equation8) as function of the integration radius in angular degrees for the 3 test areas. Left: plain; center: semi-mountainous; right: mountain. Curves correspond to the use of different geopotential models for the computation of W M as defined in section 3.2 (see legend)

Figure 8 :

 8 Figure 8: Accuracy of the geoid heights in meters computed using the standard Stokes' kernel as a function of the integration radius in three test zones. Gravity anomaly grid resolution =0.05 • . Colors correspond to different frequency ranges(see legend)

Figure 9 :

 9 Figure 9: Relative accuracy of the geoid heights computed using the standard Stokes' kernel as a function of the integration radius in three test zones. Gravity anomaly grid resolution =0.05 • .Colors correspond to different frequency ranges(see legend)

Table 1 :

 1 • ) 47.50-48.49 45.51-46.49 44.00-44.98 Statistics of heights in test areas derived from the digital model ETOPO5 [1]

	Long( • )	1.00-1.99	2.50-3.49	3.00-3.99
	Hmin(m)	84.00	247.00	142.50
	Hmax(m)	242.00	1379.70	1394.40
	Hmean(m) 131.96	523.70	958.54
	Rms(m)	134.24	576.22	989.43
	Std(m)	24.47	240.43	245.39

Table 2 :

 2 • . The black dots represent the test points in each area. Statistics of Residual gravity anomaly in test areas derived from the EGM2008 model[START_REF] Eshagh | A Strategy Towards an EGM08-based Fennoscandian Geoid Model[END_REF]

	g res	Zone(1) Zone(2) Zone(3)
	Min(mgal)	-22.8	-23.33	-27.92
	Max (mgal) 26.86	26.86	26.86
	Mean(mgal) -0.002	-0.003	-0.074
	Std (mgal)	4.63	6.12	7.34

Acknowledgements:

This work was funded by a doctoral scholarship from the University of Tichrine, Latakia, Syria. Figure 1 was drawn with the Generic mapping Tool [ref].


Stéphane Goutte

Thomas P Éran

Thomas Porcher

Social Inequalities and Vulnerability of population facing the COVID-19: the case of Seine-Saint-Denis in Ile-de-France

Keywords: Pandemic, COVID-19, Social distancing, Health system, Territorial vulnerabilities, Poverty, Housing JEL classification:I14, I18, J14, H12, R11

The vast majority of research focuses on the individual factors leading to coronavirus mortality. Numerous studies have shown that the age of the population is the dominant factor explaining mortality. Other more recent work has added gender, comorbidity, ethnicity and obesity. Based on the most populous and dense region of France -Île-de-France, grouping 8 heterogeneous departments in terms of wealth -our study seeks to identify whether economic and financial or structural factors related to housing can explain a faster circulation of the virus during social distancing like lockdown, and therefore lead to excess mortality. We show that agglomerations with higher precariousness indicators (unemployment benefit income, poverty rate, social minima in income, little or no graduate in the workforce) and less suitable housing (potentially unworthy housing, household size, overcrowded housing) are more at risk, including if their population is younger. Our study therefore provides political leaders with a number of indications allowing them to take effective measures in the event of a second wave of COVID-19 or forthcoming coronavirus pandemics.

Introduction

On the 11th of March 2020, the Director General of the World Health Organization (WHO) declared the COVID-19 epidemic as an ongoing pandemic, while more than 110 states were heavily affected worldwide, totalizing approximately 120,000 cases of coronavirus (WHO, 2020a). Very quickly, researchers from around the World devoted their work to the study of this new virus, by mainly using three different approaches. First, a race against the clock was launched by epidemiologists to find a vaccine [START_REF] Shoenfeld | Corona (COVID-19) time musings: Our involvement in COVID-19 pathogenesis, diagnosis, treatment and vaccine planning[END_REF][START_REF] Cohen | Vaccine designers take first shots at COVID-19[END_REF][START_REF] Le | The COVID-19 vaccine development landscape[END_REF] and reach in the earliest possible delay a satisfactory level of collective immunity [START_REF] Altmann | What policy makers need to know about COVID-19 protective immunity[END_REF]. Second, the medical profession devoted itself to studying the effects of the virus on the health of individuals. Finally, the major part of researchers has tried to identify the most effective ways to staunch this global scourge. Relative to this last axis in which this current study fully fits, we aim to understand what are the factors of transmission of the virus [START_REF] Li | Risk factors for severity and mortality in adult COVID-19 inpatients in Wuhan[END_REF] and worsening of the situation health (Di [START_REF] Lorenzo | Coronavirus Disease (COVID-19) in Italy: Analysis of Risk Factors and Proposed Remedial Measures[END_REF], but also to assess the extent of the economic consequences that the health crisis inevitably causes [START_REF] Mckee | If the world fails to protect the economy, COVID-19 will damage health not just now but also in the future[END_REF][START_REF] Yue | Risk Prediction and Assessment: Duration, Infections, and Death Toll of the COVID-19 and Its Impact on China's Economy[END_REF].

Fundamental elements to understand the COVID-19 phenomenon are immediately revealed by the first studies. It is thus demonstrated that the source of the pandemic is to be located in the Chinese city of Wuhan inside Hubei Province and that in the category of elderly individuals (NHS England, 2020) the highest mortality rate is to be deplored. Extremely fast, these results are refined to allow a precise identification of other comorbidity factors [START_REF] Bacon | OpenSAFELY: factors associated with COVID-19-related hospital death in the linked electronic health records of 17 million adult NHS patients[END_REF]). Thus, for example, it seems very likely that patients suffering from other pathologies such as diabetes [START_REF] Klonoff | COVID-19 in patients with diabetes: risk factors that increase morbidity[END_REF] or asthma [START_REF] Abrams | Managing Asthma during COVID-19: An Example for Other Chronic Conditions in Children and Adolescents[END_REF][START_REF] Who | Clinical management of severe acute respiratory infection when COVID-19 is suspected[END_REF] are more affected than healthiest patients, but also that the rhesus of the blood group and the ethnic origin of the patients [START_REF] Mihm | COVID-19: Possible Impact of the Genetic Background in IFNL Genes on Disease Outcomes[END_REF]NHS England, 2020;[START_REF] Hooper | COVID-19 and Racial/Ethnic Disparities[END_REF] could constitute a medical field fostering the mortality of the virus. In other words, a standard "robot portrait" of the patient most at risk in front of COVID-19 is already drafted.

On the other hand, geographical studies show that if not a single continent is sheltered1 [START_REF] Hopman | Managing COVID-19 in Low-and Middle-Income Countries[END_REF][START_REF] Gilbert | Preparedness and vulnerability of African countries against importations of COVID-19: a modelling study[END_REF], we observe within the West the highest morbidity rate, according to the following top Global ranking: United States of America (¿ 79,500 deaths), United Kingdom (¿ 31,900 deaths), Italy (¿ 30,500 deaths), Spain (¿ 26 600), France (¿ 26,300) 2 . Therefore, in order to better understand the vectors of the virus transmission as well as the morbidity factors, it seems interesting to conduct comparative studies at the three continental, regional and State levels. However, the first observations establish that biases exist and that it therefore could be reasonable to limit comparative analyzes to territorial units with the same human and climatic characteristics [START_REF] Desjardins | Rapid surveillance of COVID-19 in the United States using a prospective space-time scan statistic: Detecting and evaluating emerging clusters[END_REF][START_REF] Liu | Impact of meteorological factors on the COVID-19 transmission: A multi-city study in China[END_REF]. For example, the observation is made that in Sub-Saharan Africa, the contamination and transmission/propagation rates are extremely low compared to the countries in the North and West of the Globe [START_REF] Martinez-Alvarez | COVID-19 pandemic in west Africa[END_REF][START_REF] Nuwagira | Is Sub-Saharan Africa prepared for COVID-19?[END_REF]. The positive effects of various factors are discussed in this regard, such as the protective role of previous injections of Malaria vaccine on populations exposed to COVID-19 [START_REF] Sargin | Is Coronavirus Disease 2019 (COVID-19) seen less in countries more exposed to Malaria?[END_REF]. Moreover, the global death Report show that the counting of national deaths is not uniform. Some countries report exclusively deaths in hospital (like France at the early stage of the pandemic) while others merge deaths in hospital, domestic and nursing homes (like Germany). Accordingly, an international study seems to be unrealistic at the moment.

Compared to the other research papers mentioned above, our study proposes an approach to the structural factors of contamination/contagion/mortality, rather than a focus on human factors. In other words, it is a question of asking whether, in addition to genetics and geography, new elements can be put forward to explain the excess mortality in certain populations. For this, we have limited our research to Île-de-France (see Figure (2.1)), a French region built of eight departments3 , which has the unique characteristic of not constituting a cluster of contamination due to an identifiable and outstanding event.

Indeed, certain territorial units have formed clusters of contamination largely responsible for the degree of contamination at the National level. For example, in the French city of Mulhouse belonging to the Grand Est region, a major gathering of faithful evangelists is considered to be responsible for a nonnegligible part of the whole excess mortality linked to COVID-19 among the State. In addition, the Île-de-France region is highly populated with 12 174 880 million inhabitants (19% of the whole French population) and socially heterogeneous (ethnicity, professional qualification of workers, graduate of higher education, quality of the health system...) but its boarders stand inside a small geographic area with no climatic ecosystems effects. Under these conditions, Îlede-France constitutes a relevant field of study for the various structural factors other than individual ones like age or comorbidities promoting the contamination/contagion/mortality tryptic.

Data and Approach

Starting with a large panel of 66 variables4 which were representative of the economic, financial or structural factors related to housing in Ile de France and its population. This data collects, analyses and disseminates information on the French economy and society. We finally select using a Principal Component Analysis (PCA) a closer panel of 30 variables which are all more significant in term of segmentation of our Departments (and in particular the Seine-Saint-Denis).

Principal Components Analysis

Thus in order to characterize the best set of discriminant variables we proceed with a principal component analysis. This approach allows us to best capture the explanatory and segmenting power of the variables we have. If we want to go more deeply in the explanation of this segmentation, we can take a look on the weights of each variables given by each two axis. These results are provided in Table 1. We can see that a positive value on the first axis (i.e. horizontal) characterizes:

• a high share of graduates of higher education in pop. out of school 15 years or more at a level of 96.50%

• Average household size at a level of 95.96%

• a high value of the aging index at a level of 93.60%

• a high average hourly net salary at a level of 85.45%

This means that a department with a high coordinate in this Factor 1 exhibits all this points; and higher is its coordinate higher is these factors. The projection of IDF Departments on the these two main PCA axis are stated in Table 2.

The Departments of Paris and Haut de Seine which take a high value in this axis projection are so fundamentally and intrinsically characterized and determined by a population with a high level of education, with a higher salary than the other departments and an older population than the other one. This last factor is of course the main reason why the mortality are important in both departments.

Conversely, the Seine-Saint-Denis department which takes the most negative value in this projection is so characterized by a younger population with a less level of education and so an average value of salary at the end. But, as we show previously, its excess of surmortality dues to COVID-19 is the higher. We so have to look the interpretation of the second axis (i.e. vertical) since the Seine-Saint-Denis is again isolated than the other departments in the upper region (i.e. positive values).

We can see that a positive value in this factor characterizes:

• a high number of main residences overcrowded at a level of 96.00%.

• a high share of private park accommodation potentially unworthy (PPPI) at a level of 91.44%

• a high number of people living in an apartment in a household of at least 4 people at a level of 88.83%

• a high poverty rate at a level of 88.20%

• a high value of share of unemployment benefits in the rev. avail. at a level of 72.74%

This means that the Seine-Saint-Denis is highlighted by very difficult economic and health conditions: overcrowded housing, potentially unworthy; a low-income population, mostly from income linked to unemployment benefits. Hence, conditions socio-economics which bring surmortallity in period of pandemic.

Results

The list of these variables can be see in Tables 3 to 7 in Appendix. To compare the values of these set of variables we decided to evaluate the variation in percentage of each value for each department with respect to the average of the region. It means that a value of 10% represents that this department has a value 10 % higher than we average of all departments in Ile de France. Our study provides interesting results. We note in Figure 3 that the link between the population over age of 75 and excess mortality is not absolute.

Potentially unworthy housing

Household size

Overcrowded housing Two departments with the highest population deltas over age of 75, Paris (+14.38%) and Yvelines (+14.38%), are among the departments with the lowest excess mortality (respectively +73,90% and +66.60%). Conversely, while Seine-Saint-Denis department displays the lowest delta on the population over 75 (-38.51%) it shows the highest excess mortality (+128.10%). Theoretically, the standard observation would have been the opposite side: The high mortality rate observed among people over 75s in France, representing 78.3% of deaths with an average age of 81.2 (Santé Publique France, 2020), should have led to a negative ranking on such departments. The Val-d'Oise is also a department with a negative delta regarding the population over 75 years old (-15.54%) but with the fourth excess mortality in Île-de-France (+88.6%). Seine-et-Marne department has also a smaller population of over 75 (-7.41%) associated to an excess mortality rate of +71.70%. Our study allows us to identify a broader number of factors, first, and to point out the specificity of each department with a significant excess mortality despite their more advantageous demography compared to others, second. By taking economic, social and financial variables that can reveal the insecurity of department populations, such as unemployment benefits, poverty rate, minimum social benefits or level of education, and by adding variables specific to the structure of housing, we offer a chance to implement tailor-made structural policies. For instance, in regard to unemployment benefit income, we observe that Seine-Saint-Denis and Val-d'Oise are the only departments to have a positive delta (respectively +39.13% and +7.51%), with a very clear demarcation for Seine-Saint-Denis (see Figure 3). Among the cluster, all the other departments have negative delta (see unemployment benefit income in Figure 3). We come to the same observation according to the poverty rate, where four departments have positive deltas with a clear demarcation of the Seine-Saint-Denis (+84.07%), then we rank Val-d'Oise (+8.12%) and Val-de -Marne (+7.48%) (see Poverty rate in Figure 3). We find similar result at the observation of social minima where three departments (Seine-Saint-Denis, Vald'Oise and Val-de-Marne) have all positive deltas with a clear demarcation for Seine-Saint-Denis (+118.20%), Val-d'Oise and Val-de-Marne tied (+5.79%) (see Social minima in income in Figure 3). Finally, concerning the part of individuals without diploma into the workforce, Seine-Saint-Denis still occupies the first place with a delta of +46.12% compared to the average of the cohort, followed by Val-d'Oise (+11.60%), Seine-et-Marne (5.84%) and Val-de-Marne (+4.69%) (see graph little or no graduate in the work force in Figure 3). Based on the analysis of economic and financial variables, the first conclusion that can be drawn is that there are therefore several common points between Seine-Saint-Denis and Val-d'Oise, two departments having a smaller population of 75+ but suffering a significant excess mortality, despite social distancing measures implemented by the French Government. Indeed, since the promulgation of the Law 2020-290 of March 23, 2020 said 'Emergency to face the epidemic of COVID-19', prolonged by the Law 2020-546 of May 11, 20205 , the French Government is authorized to rule into legislative matters by decree when it comes to fight against the epidemic in France. In addition, regarding inequalities related to the structure of housing, specifically the unworthily housing, the two departments with positive deltas are Seine-Saint-Denis (+104.77%) and Paris (+77.47%) (see Potentially unworthly housing in Figure 3). For the average size of households, five departments have a positive delta: Seine-Saint-Denis (+9.47%), Seine-et-Marne (+5.26%), Yvelines (+1.05%), Essonne (+5.26%) and Val-d'Oise (+9.47%) (see Household size in Figure 3). Finally, four departments have positive deltas at the variable "overcrowded main residences": Seine-Saint-Denis (+68.16%), Paris (+29.79%), Hauts-de-Seine (+6.12%) and Val-de-Marne (+ 11.83%), with a delta far above for Seine-Saint-Denis (see Overcrowded housing in Figure 3).

Conclusion and opening to future work

Seine-Saint-Denis differs from other departments according to a number of important variables. These variables, on the main field of financial economic poverty, on the one hand, and structural factors linked to housing, on the other hand, may explain the excess mortality during social distancing and lockdown put in place by French Government. Six of these seven variables are also significant in Val-d'Oise, another department which, like Seine-Saint-Denis, has a significant excess mortality with a lower proportion of people over the age of 75.

Our study therefore provides political leaders with a number of inputs allowing them to better implement effective measures in the event of a second wave of COVID-19 or new pandemics due to viruses inside the coronavirus family. Agglomerations with higher precariousness indicators (unemployment benefit income, poverty rate, social minima in income, little or no graduate in the workforce) and less suitable housing (unworthy housing, household size, overcrowded housing) are more at risk, including when their population is younger. It is therefore a requirement to set up new health policies facilitating an accurate monitoring of the inhabitants and their environment in these departments or agglomerations, with the main objective of breaking human-to-human transmission chains more quickly and efficiently. It would be interesting in an additional work to see if, in other countries and other continents, the results obtained at

Figure 1 :

 1 Figure 1: Map of Île-de-France

Figure 2 :

 2 Figure 2: Projection of IDF Departments on the two main PCA axis

Figure 3 :

 3 Figure 3: Cumulative Variability obtained for each factor axis

Figure 4 :

 4 Figure 4: Link between Age and Excess mortality

Figure 6 :

 6 Figure 6: Inequalities linked to housing

Table 1 :

 1 Explication Weights of each variables on the two main axis factor

	F1	F2

Table 3 :

 3 Values of the delta percentage of our panel of data for each departments -Part I

	Departments	Excess mortality	0 to 19y	20 to 39y	40 to 59y	60 to 74y	75y and more
	Paris	73,90% -37,41% 18,28% -4,97%	7,17%	14,38%
	Seine-et-Marne	71,70%	5,55%	-8,82%	2,29%	3,93%	-7,41%
	Yvelines	66,60%	2,34%	-15,93% 3,63%	6,78%	14,38%
	Essonne	88,20%	4,76%	-6,47%	0,98%	-0,54%	2,45%
	Hauts-de-Seine	127,80% -6,54%	3,56%	1,31%	-3,09%	8,48%
	Seine-Saint-Denis 128,10% 10,44%	4,51%	-2,28% -15,50% -38,51%
	Val-de-Marne	96,50%	-1,70%	1,44%	-0,13% -0,93%	2,71%
	Val-d'Oise	88,60%	9,18%	-4,52% -1,36% -1,22% -15,54%
	6 Appendix

Table 4 :

 4 Values of the delta percentage of our panel of data for each departments

	-Part II	2019-Q4 quarterly unemployment rate	Population density	Average household size	Median of standard of living	Share of taxable households	Share of unemployment benefits in disposable income
	Paris	-14,04% 270,31% -20,00% 15,85%	8,98%	-5,14%
	Seine-et-Marne	-7,11%	-95,79%	5,26%	-2,43%	-1,78%	-5,14%
	Yvelines	-14,04% -88,83%	1,05%	11,60%	10,40% -14,62%
	Essonne	-9,88%	-87,29%	5,26%	-0,06%	2,81%	-8,30%
	Hauts-de-Seine	-12,65% 62,69%	-7,37%	14,83%	10,72% -11,46%
	Seine-Saint-Denis 44,19%	21,21%	9,47%	-26,55% -24,40% 39,13%
	Val-de-Marne	-0,17%	0,22%	-3,16%	-5,11%	-1,46%	-1,98%
	Val-d'Oise	13,69% -82,52%	9,47%	-8,13%	-5,26%	7,51%

Table 5 :

 5 Values of the delta percentage of our panel of data for each departments -Part IIIShare of social minima in disposable incomePoverty rate Share of those with little or no education in the out-of-school population aged 15 and over Share of higher education graduates in the out-of-school population of 15 years or more Share of apartments in total housing Share of houses in total housing

	Paris	-27,27%	1,69%	-30,97% 52,54%	44,56% -97,14%
	Seine-et-Marne	-7,44%	-25,34%	5,85%	-25,87% -39,52% 84,59%
	Yvelines	-47,11% -37,57% -15,24% 10,18% -16,22% 35,35%
	Essonne	-14,05% -16,98% -0,29% -10,75% -22,49% 49,01%
	Hauts-de-Seine	-33,88% -21,48% -21,76% 31,61%	29,62% -63,15%
	Seine-Saint-Denis 118,18% 84,07%	46,12% -33,69% 10,51% -22,48%
	Val-de-Marne	5,79%	7,48%	4,70%	-2,43%	13,50% -29,47%
	Val-d'Oise	5,79%	8,13%	11,60% -21,59% -19,96% 43,29%

It also seems that not a single country has been sheltered and that the few localities where no deaths have been recorded have chosen not to report the cases. Refer to figures from Johns Hopkins University which are widely accepted by the Global scientific Community. Available at: https://coronavirus.jhu.edu/map.html (accessed 10 May

2020).2 Figures updated toMay 11, 2020.

In France, administrative levels in order of importance (ascending order) are municipalities/agglomerations (35 357 units), departments (101 units) and regions (18 units).

taken from the statistical database of The National Institute of Statistics and Economic Studies (INSEE)

Refer to https://www.legifrance.gouv.fr/.


Donia Aloui

Stephane Goutte

Khaled Guesmi

Rafla Hchaichi

COVID 19's impact on crude oil and natural gas S&P GS Indexes

Keywords: COVID-19, oil price war, health crisis, drop oil price, S&P GS commodities Indexes, TVP-SVAR model, crude oil, natural gas. JEL Classification:

 to asses the impact of COVID-19 shocks on the energy futures markets, particularly on crude oil and natural gas S&P GS Indexes. The findings confirm that energy commodities S&P GS Indexes respond to COVID-19 shock that varying over time due to fundamentals factors as well as behavioral and psychological factors.

Introduction

For a long time, oil has been considered the engine of the world's economy. Crude oil so-called "black gold" is considered one of the world's most precious commodities that affect the global economy as well as world trade. Indeed, the sudden price drops or unexpected can send the global world into uncertainty and fears. Crude oil prices are subject to booms and busts and it depends on 1 This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=3587740 the economic conjuncture, business cycle, and crisis. Previously, supply shortages triggering by political events have caused change dramatically in oil prices, such as the Iranian revolution, Iran-Iraq war, Arab oil embargo, and Gulf wars. Besides, the Asian financial crisis and the global economic crisis of 2008 have been generating also price fluctuations [START_REF] Baumeister | Forty years of oil price fluctuations: Why the price of oil may still surprise us[END_REF]).

In 2014, after reaching monthly peaks of $112, crude oil prices fell sharply following a production that exceeded demand [START_REF] Baumeister | Understanding the decline in the price of oil since june 2014[END_REF]). Currently, we are facing the biggest oil price drop since the Gulf War in 1991 [START_REF] Opec | Organization for the petroleum exporting countries[END_REF]). This is illustrated by the dramatically falling since Mars 12, 2020 in reaction to the coronavirus pandemic that led an unexpectedly sharp drop in demand for oil [START_REF] Opec | Organization for the petroleum exporting countries[END_REF]). The current drop in oil demand is essentially caused by the quarantine of countries particularly restriction and containment measure that causes a drop in consumption. Besides, an oil price war went off following a spat about oil production between Russia and the Saudiled Organization of Petroleum Exporting Countries (OPEC). OPEC wanted to slash production to match a drop off in demand and the global economic turndown, while Russia was not cooperative (OPEC (2020b)).

It seems that Russia, by refusing to reduce production quotas, undermine the extraction of shale oil in the United States. According to Rosneft, which is a Russian state-owned company specializing in the extraction, processing, and distribution of petroleum, if Russia had accepted a drop in oil production to maintain the price, American shale oil would have found takers on the world energy market, which would not have been favorable for the flow of Russian oil. However, Saudi Arabia and Russia failed to reach an agreement on production cuts that exacerbate the crisis (OPEC (2020b)). As a result, the price of U.S. crude oil was fluctuating just around $30 per barrel. As we can see, the crude oil prices can vary greatly from $150 per barrel in 2014 to $30 in 2020 [START_REF] Alkhathlan | Analysis of saudi arabia's behavior within opec and the world oil market[END_REF]).

It turns out that crude oil prices react to many variables, including, the current health crisis caused by coronavirus outspread, economic news, overall supplies, consumer demand, and speculators and hedgers behavior. Besides, OPEC as an international oil-producing cartel plays a considerable role in determining global oil supplies that influence oil prices [START_REF] Kaufmann | Determinants of opec production: Implications for opec behavior[END_REF]; [START_REF] Kaufmann | Does opec matter? an econometric analysis of oil prices[END_REF]). The oil supply is also determined by the new technology that allows extracting oil from the rock so-called shale oil his allowed the USA to become the world's largest producer of oil in 2018 and a major source of global oil supplies [START_REF] Frondel | The us fracking boom: Impact on oil prices[END_REF]).

Presumably, oil demand fell sharply as a result of the global COVID-19 pandemic (OPEC (2020b)). In December 2019, a virus designed as a severe acute respiratory syndrome (SARS) has spread in Wuhan, a Chinese city. The first contaminated case was a 55-year-old patient reported on November 17, 2019. The cases infected has increased to twenty-seven, on December 15 and reached to sixty in five days. It has been identified as a viral zoonotic disease later named the coronavirus (COVID-19) that was spread rapidly in China and was declared, on 30 January 2020, as being a public health emergency of international concern that causing 811 deaths within 3373 confirmed cases, on 8 February 2019 (WHO (2020)).

Too soon, the World Health Organization (WHO (2020)) announced that the disease can be transmitted between humans by droplets ejected, direct contacts or via contaminated objects and an early detection process must be put in place to avoid the spread of this virus with exponential growth. Unfortunately, the coronavirus has, speedily, emerged as an alarming contagious disease increasing dramatically the number of death that requires a rigorous containment measure around the world. On 11 March 2020, the epidemic was qualified as a pandemic by the WHO (2020). The COVID-19 pandemic has spread rapidly to many other countries and territories, causing border closures, a sudden slowdown in the global economy and a stock market crash on March 12, 2020, due to the failure of negotiations of OPEC to reach a production quotas agreement.

Although, there is a vast oil market literature [START_REF] Frondel | The us fracking boom: Impact on oil prices[END_REF]; [START_REF] Chuliá | Volatility spillovers in energy markets[END_REF]; [START_REF] Charlier | From residential energy demand to fuel poverty: income-induced non-linearities in the reactions of households to energy price fluctuations[END_REF]; [START_REF] Cleveland | Oil supply and oil politics: Deja vu all over again[END_REF]; [START_REF] Baumeister | Understanding the decline in the price of oil since june 2014[END_REF]; [START_REF] Kilian | How the tight oil boom has changed oil and gasoline markets[END_REF]), the health crisis only recently has become a topic of interest, and the economic impact of the COVID-19 on oil price is not widely analyzed yet. To fill this gap in literature this study performs the TVP-SVAR model to asses the COVID-19 chocks on the energy market. The reason we choose the energy market, particularly crude oil and natural gas, is due to the following considerations. Commodity markets are closely related to economic conjuncture, in contrast to financial assets. They play a considerable role in related risk-sharing. As a result, energy commodities, particularly crude oil and natural gas have become an asset class, which is widely used for diversification, hedging or speculation prospects [START_REF] De Roon | Hedging pressure effects in futures markets[END_REF]; [START_REF] Gorton | Facts and fantasies about commodity futures[END_REF]; [START_REF] Jorion | Firm value and hedging: Evidence from us oil and gas producers[END_REF]). So, investors, speculators, hedgers, and policymakers consider commodity markets as an alternative investment instrument for hedging against the equity market risk [START_REF] Gatfaoui | Diversifying portfolios of us stocks with crude oil and natural gas: A regime-dependent optimization with several risk measures[END_REF]).

The remainder of this paper is organized as follows. Section 2 provides a related literature review. Section 3 introduces the methodology, data, and the model investigating the effect of coronavirus on crude oil and natural gas price. In Section 5, we discuss the estimation results and we draws the main conclusions.

2 Related literature 2.1 The coronavirus's impact on world economy

The outbreak of coronavirus has disrupted china's economy ensuring uncertainties that slowed world trade. China's demand for raw materials has declined which penalized the commodity-exporting countries such as Saudi Arabia, Russia, and Australia. The shutdown of many Chinese companies has resulted in a shortage of spare parts in other countries' industries, such as the Hyundai company that has shut down factories in South Korea due to a lack of supplies. Pharmaceutical and manufactured supply shortages get worse due to the disruption of Chinese factories. In some affected countries such as Italy and Hong Kong, the pandemic has generated panic purchases and food shortages OPEC (2020b).

Besides, the impact on tourism and airlines has worsened so Thai tourism will be affected by the number of Chinese tourists that could drop by 20% in 2020. According to the OECD (2020), global growth could drop at 2.4% this year, compared to 2.9% expected and the central banks could drop key rates to revive the economy by adopting expansionary monetary policy as did by several central banks in several countries. Several studies have highlighted the impact of economic activity on the oil price. [START_REF] Baumeister | Understanding the decline in the price of oil since june 2014[END_REF] shows that a slowing on the global economy has reduced the Brent price per barrel by $49 between June and December 2014. [START_REF] Kilian | How the tight oil boom has changed oil and gasoline markets[END_REF] explained this decline by a weakening of the demand for crude oil. [START_REF] Yardimcioglu | Opec ülkelerinde hollanda hastaligi: Petrol fiyatlan ve ekonomik büyüme iliskisinin ekonometrik bir analizi/dutch disease in the opec countries: An econometric analysis of the relationship between oil prices and economic growth[END_REF] shows that there is a strong cointegration relationship between oil prices and economic growth in the 10 OPEC countries over the period from 1970 to 2011 and he confirms the impact of Dutch Disease on price oil fluctuations. Several empirical studies have been focused on understanding the macroeconomic consequences of oil price shocks on the US economy. [START_REF] Ghalayini | The interaction between oil price and economic growth[END_REF] investigates whether economic world growth can be explained by changes in the oil price. For the oil-importing countries, oil price increase and economic growth are negatively correlated while the relationship is positively correlated for oil-exporting countries. This study covered the G-7 group, OPEC countries in addition to Russia, China and India.

Recently many studies investigated COVID-19's impact on global economic activity. [START_REF] Yilmazkuday | Coronavirus disease 2019 and the global economy[END_REF],by using daily data on the global coronavirus disease 2019 (COVID-19) deaths, shows that the effect of COVID-19 deaths is insignificant, and explained the plunge of oil price on March 2020 by the OPEC disagreement. In the same line of work, [START_REF] Maijama'a | Corona virus outbreak and the global energy demand: A case of people's republic of china[END_REF] investigates the impact of coronavirus outbreak on the global energy demand by using daily data on china's population, currency Chinese exchange rate and international crude oil prices, from 23rd January to 8th February 2020. The result revealed that the total population has a positive and significant impact on total coronavirus infected cases while the crude oil price is negative and significantly related to the coronavirus infected cases.

Determinants of oil price fluctuation

The fundamental driver of oil's price is supply and demand in the market. So, oil supply is controlled by a cartel of oil-producing countries called OPEC [START_REF] Alkhathlan | Analysis of saudi arabia's behavior within opec and the world oil market[END_REF]) and Oil demand is driven by production and transport activity such as electricity, gasoline, and kerosene. A large literature on the economic determinants of oil price fluctuations has emerged. Initially, oil price fluctuations were explained by the disruptions of global oil production associated with political events such as wars and revolutions in OPEC countries [START_REF] Hamilton | What is an oil shock?[END_REF]). Subsequent research has shown that oil price fluctuations since 1973 are explained by shifts in the crude oil demand [START_REF] Kilian | The economic effects of energy price shocks[END_REF]; Kilian (2017); [START_REF] Gong | Time-varying effects of oil supply and demand shocks on china's macro-economy[END_REF]; [START_REF] Kilian | Oil price shocks: Causes and consequences[END_REF]; [START_REF] Baumeister | Forty years of oil price fluctuations: Why the price of oil may still surprise us[END_REF]; [START_REF] Kilian | Did unexpectedly strong economic growth cause the oil price shock of 2003-2008[END_REF]; [START_REF] Kilian | Quantifying the speculative component in the real price of oil: The role of global oil inventories[END_REF]).

According to [START_REF] Hamilton | What is an oil shock?[END_REF], the oil price shock of 1973-74 has been explained by the lack of supply that led to the rise of its price. [START_REF] Hamilton | What is an oil shock?[END_REF] attributes the Arab oil production cuts to the Arab oil embargo against selected Western countries, which lasted from October 1973 to March 1974. This crisis was followed by a second oil crisis in 1979-80, when the West Texas Intermediate price rose from less than $15 in September 1978 to almost $40 in April 1980. [START_REF] Hamilton | What is an oil shock?[END_REF] explained this oil price surge by the reduction in Iranian oil production following the Iranian Revolution.The price rise was, also, explained by an increased inventory demand in anticipation of oil future shortages and high future oil demand from a booming global economy. From 1980until 1990[START_REF] Hamilton | What is an oil shock?[END_REF] attributes the increase in the West Texas Intermediate price of oil from $36 in September to $38 in January 1981 to an oil supply disruption following the outbreak of the Iran-Iraq war causing the destruction of Iranian oil facilities. However, the rise of the U.S. interest rates and the global recession that emerged lowered the oil demand and hence its prices. Hence, this economic conjuncture caused a sell-off of the oil inventories accumulated in 1979. The most surge of the oil price since 1979 has occurred between mid-2003 and mid 2008, the West Texas Intermediate price climbed from $28 to $134. The price surge was caused by an increase in crude oil demand. [START_REF] Kilian | The economic effects of energy price shocks[END_REF], [START_REF] Hamilton | International energy agency. An unprecedented global health and economic crisis[END_REF], and [START_REF] Kilian | Did unexpectedly strong economic growth cause the oil price shock of 2003-2008[END_REF] have shown that this dramatic increase in demand was explained by an unexpected expansion of the global economy driven by a strong rise in oil demand from Asiatic emerging countries. This surge in oil prices in the physical market can be explained as the result of speculator behavior in the future market. Recently the financial crisis of 2008 illustrates the powerful effects of a sharp drop in commodities demand worldwide that was sharply curtailed in the second half of 2008 as crude oil that plummeted, causing a fall in the price of oil from $134 in June 2008 to $39 in February 2009.

Currently, as a result of the COVID-19 pandemic, global industrial production is reduced causing lower oil consumption that decreases the barrel price. As announced, on February 15, 2020, by the International Energy Agency (IEA (2020)), global demand could decrease by 425,000 barrels per day following a dramatic drop in Chinese consumption. Global oil demand drop compared to 2019: China, which consumes about 14% of world oil production and air carriers, which consume 8% of world demand, have greatly reduced their activities. As a consequence, OPEC decides to regulate the supply of oil to set the price on the world market while other oil-producing countries, such as Russia, have braked OPEC's ability to control supply.

The OPEC decision consists of reducing their production by 1.5 million barrels per day during the second quarter of the year by negotiating with Russia as well as other producer countries. On March 6, 2020, Russia refused to collaborate causing the fell by 10% of the price of oil. Hence, Saudi Arabia cuts the price of its oil and decides to increase its production, to force Russia to return to negotiating with them (OPEC (2020a)). Previous studies have focused on OPEC behavior. [START_REF] Kaufmann | A model of the world oil market for project link integrating economics, geology and politics[END_REF], [START_REF] Frondel | The us fracking boom: Impact on oil prices[END_REF], [START_REF] Kaufmann | Does opec matter? an econometric analysis of oil prices[END_REF], and [START_REF] Dees | Modelling the world oil market: Assessment of a quarterly econometric model[END_REF] have investigated OPEC's influence on oil prices over the medium-and long-run.They found that the key decision variables are the OPEC production quota and the OPEC ability to add production capacity. [START_REF] Dees | Modelling the world oil market: Assessment of a quarterly econometric model[END_REF] emphasizes that if the OPEC cartel fails to cooperate, it results in a sharp drop in the crude oil price, as the competition among OPEC countries leads to supply quantities that surpass demand which will be stored.

Oil Prices and future market

As a component of derivative markets, futures markets are hedging markets allowing economic agents to protect themselves against price risk. In recent years, traders have significantly increased their investment in the commodities futures market that has led to a significant development of speculation in these markets. Several researchers have studied the importance of commodity speculation and have analyzed the bidirectional relationship between spot prices and futures prices, [START_REF] Buchanan | Which way the natural gas price: an attempt to predict the direction of natural gas spot price movements using trader positions[END_REF]; [START_REF] Sanders | Hedgers, funds, and small speculators in the energy futures markets: an analysis of the cftc's commitments of traders reports[END_REF]; [START_REF] Bu | Price dynamics and speculators in crude oil futures market[END_REF]). On the one hand, the assessment of future contracts is done by using the cost-of carry model, which requires that the futures price should depend on spot price and the cost of carrying or storing goods [START_REF] Nicolau | Are Futures Prices Influenced by Spot Prices or Vice-Versa?: An Analysis of Crude Oil, Natural Gas and Gold Markets[END_REF]).On the other hand, transaction and speculation in the futures market could affect oil spot prices [START_REF] Kaufmann | The role of market fundamentals and speculation in recent price changes for crude oil[END_REF]; [START_REF] Wen | Exploring the dynamic effects of financial factors on oil prices based on a tvp-var model[END_REF]).They assume that oil price fluctuations are not only correlated to the real fundamentals of the economy but also follow speculative logic. According to behaviour theory, it turns out that investors' sentiment can influence oil prices. Thus, when they anticipate a decrease in oil demand in the future, investors made short selling oil futures contracts that leads to a dramatic drop in prices, this is meaning that prices can depend on market psychology [START_REF] Kaufmann | Oil prices, speculation, and fundamentals: Interpreting causal relations among spot and futures prices[END_REF]). The behavioural finance theory emphasize the importance of psychology in the evolution of commodity market prices [START_REF] Deaton | On the behaviour of commodity prices[END_REF]; [START_REF] Deaton | Competitive storage and commodity price dynamics[END_REF]; [START_REF] Chambers | A theory of commodity price fluctuations[END_REF]; [START_REF] Shiller | From efficient markets theory to behavioral finance[END_REF]).

Other researches suggest that the futures market movements involve relevant information on commodity spot price evolution. They are based on the fact that fluctuations of the oil prices are difficult to predict and extremely uncertain that's why economic agents use the crude oil price of the futures markets as a measure to predict future spot prices. In other words, oil futures prices are used to determine oil price expectationsBaumeister and [START_REF] Baumeister | Understanding the decline in the price of oil since june 2014[END_REF]. This relationship between spot prices and futures prices is based on the wisdom of the participants in the financial market, which are considered as the best-informed agent about the financial and economic environment of oil price determinants.

Therefore, it is interesting to study the dynamics of futures prices during this health crisis and to assess the investor's behavior in the face of this unpredictable crisis. This paper discusses whether the health crisis of COVID-19 could affect crude oil and natural gas in futures markets.

Methodology

Due to the rapid and unpredictable development of this health phenomenon, there have been structural changes in the commodity market, which will cause changes in investor behavior on oil prices at different times. For that reason, it is interesting to examine the time-varying effects of COVID-19 on commodity indices. Therefore, we opt for a Structural VAR model with time-varying coefficients and stochastic volatility (TVP-SVAR model) developed by [START_REF] Primiceri | Time varying structural vector autoregressions and monetary policy[END_REF]. This model was recently adopted by [START_REF] Gong | Time-varying effects of oil supply and demand shocks on china's macro-economy[END_REF] and [START_REF] Wen | Exploring the dynamic effects of financial factors on oil prices based on a tvp-var model[END_REF] to examine macroeconomic and financial determinants of oil price fluctuations. This model allows us to measure variables' responses to structural shocks and capture their changes over time.

The model

In order to asses the impact of Covid-19 chocks on the energy market, we use a Structrual VAR model with time varying coefficients and stochastic volatility (TVP-SVAR model) as presented by [START_REF] Primiceri | Time varying structural vector autoregressions and monetary policy[END_REF]. The model can be expressed as follows:

Y t = c t + b 1,t Y t-1 + +b p,t Y t-p + u t (1)
Where Y t is a vector composed of two dependent variables, namely: the SP Goldman Sachs indices of Crude Brent and Natural Gas and a shock variable related to covid-19. c t is n × 1 vector of the time variable constant , b p,t is n × n matrix of coefficients varying over time and u t is n × 1 vector of structural innovations which are normally distributed.

The covariance variance matrix of the residues u t ,which is varying over time, is expressed as follows:

V ARu t = Ω t = A -1 t H t (A -1 t )
(2)

H t = Σ t Σ t (3) Σ t =       σ 1,t 0 .... 0 0 σ 2,t . . . . . . . . . . . . . . . 0 0 ... 0 σ n,t       (4) A t =       1 0 .... 0 a 21,t 1 . . . . . . . . . . . . . . . 0 a n1,t ... a nn-1,t 1       (5)
The TVP-SVAR model can be summarized as follows:

y t = X t C t + A -1 t Σ t ε t (6)
Where :

X t = I ⊗ [1, y t-1 , ..., y t-l ] (7)
C t = vec([c t , b 1,t , ..., b p,t ]) (8) and V AR(ε t ) = I n (9)
This model supposes the variation of the matrix A t . This econometric specification allows a structural shock to generate effects varying overtime on the variables studied and thus to capture the change in the relationships between the variables over time. In order to obtain the impulse response functions and stochastic volatility, we estimate our model according to the method of [START_REF] Primiceri | Time varying structural vector autoregressions and monetary policy[END_REF].

Estimation

In order to obtain the impulse response functions and stochastic volatility, we estimate our model according to the method of [START_REF] Primiceri | Time varying structural vector autoregressions and monetary policy[END_REF]. His method consists of estimating the covariance variance matrix Ω t which is varying over time in three blocks. The first block consists of the vector of the coefficients C t = [c t , b 1,t , ..., b p,t ] the second block consists of the non-zero and non-unitary elements of the matrix A t namelyk t where k t = (a 21,t , . . . , a n1,t , . . . a n1,t ) and the third block h t includes the standard deviations σ 1,t ; σ 2,t ; . . . σ n,t which are located on the diagonal of the matrix Σ t with h t = ln(σ t ) is the stochastic volatility of innovations. This process can replace the ARCH specification. [START_REF] Primiceri | Time varying structural vector autoregressions and monetary policy[END_REF] In addition, the coefficients of the model follow an autoregressive random walk as follows:

C t = C t-1 + ω t (10)
Similarly, the non-zero and non-unitary elements of the matrix A t follow an autoregressive random process:

k t = k t-1 + ϕ t (11)
While the block representing volatility follows a geometric random walk:

t = lnσ t = lnσ t-1 + t ( 12
)
h t = h t-1 + t (13)
We assume that the residual vectors follow the normal distribution as follows:

    t ω t ϕ t ε t     ∼ N (0, V ) (14)
V t =     I n 0 0 0 0 Q 0 0 0 0 S 0 0 0 0 W     (15)
It is the identity matrix of dimension n. Q, S, W are matrices that represent the hyper-parameters. According to [START_REF] Primiceri | Time varying structural vector autoregressions and monetary policy[END_REF] the most suitable method for estimating these parameters is the Bayesian simulation method. Following the work of [START_REF] Primiceri | Time varying structural vector autoregressions and monetary policy[END_REF] and [START_REF] Korobilis | Assessing the transmission of monetary policy shocks using dynamic factor models[END_REF], we use the Gibbs sampler to estimate the parameters in question. It is a Bayesian estimation method based on the Markov chain Monte Carlo simulation method (MCMC). This technique allows the determination of the posterior distribution of the parameters and hyper-parameters.

Priors

The estimation by the TVP-SVAR model requires the fixing of several hypotheses beforehand. We follow the choices adopted by [START_REF] Primiceri | Time varying structural vector autoregressions and monetary policy[END_REF] to establish these hypotheses. We assume, first of all, that the coefficients, the parameters as well as the hyper-parameters are independent. The first 30 observations represent the training sample from which we can extract prior distributions. Based on the work [START_REF] Primiceri | Time varying structural vector autoregressions and monetary policy[END_REF] the coefficient matrix is normally distributed, the average is obtained by the OLS ordinary least squares method from the pre-sample.

C 0 ∼ N ( C OLS , 4.V ( C OLS ) (16)
The previous distribution of the standard deviation as well as the lower elements of the covariance variance matrix A t also follow the normal distribution with an average estimated by the OLS ordinary least squares method:

ln σ0 ∼ N ( σ OLS , I n ) (17) A 0 ∼ N ( A OLS , 4.V ( A OLS ) (18)
Q, W and S are assumed to be distributed according to the inverse Wishart distribution expressed as follows:

Q ∼ IW(k 2 Q .30.V ( C OLS , 30) (19)
The value 30 corresponds to the size of the pre-sample.

W ∼ IW(k 2 w .4.I n , 4) ( 20)

S 1 ∼ IW(k 2 S .2.V ( A 1OLS , 2) (21) S 2 ∼ IW(k 2 S .3.V ( A 2OLS , 3) (22)
S 1 and S 2 are two blocks that constitute S and A 1 , A 2 two blocks of A. k 2 Q = 0, 01, k 2 W = 0, 01 and k2 S = 0, 1. This model provides flexible estimate of the mechanism for transmitting structural shocks.

Data and Results

Data

The aim of this paper is to analyze the impact of COVID 19 shocks on the energy market. We use the number of deaths to measure the evolution of COVID 19 and S& P Goldman Saches indices. According to the Commodity Future Trading Commission (CFTC), the value of the Goldman Sachs indices increases with the futures prices of the selected contracts and vice versa. We use daily data from January 2, 2020, to April 9, 2020. The objective is to study the impact of a COVID 19 shock on the oil and natural gas futures indices, distinguishing 3 periods: (i)the first period, from February 22 to March 10, is before WHO (2020) announces a global pandemic, (ii)the second period from March,10 to March 20: during this period the WHO (2020) announces a global pandemic and (iii) the third period from 20 March to April, 9: this period is characterized by a high spread of the Virus in the USA. All variables are studied in the first difference.

Results

Stochastic volatility

Figure 1 shows the stochastic volatility of the three studied variables: deaths from COVID-19, S&P GSCI Crude Oil and S&P GSCI. We can observe the high volatility of the number of deaths between March, 12 and April, 9. This period is characterized by a very rapid worldwide spread of the virus accompanied by a significant increase in the number of deaths especially in Italy, France, Spain, and the USA. The volatility curve of the Goldman Sachs crude oil index shows high volatility during February to attend a peak during 6 March and 10 March. This period precedes the announcement of the World Health Organisation (WHO) that qualifies the virus as a global pandemic on 11 Mars, 2020. Investors changed their investment strategy during the first period of the coronavirus spread when the oil spot price registered the first drops leading to high volatility in the future commodity market. We can notice a peak in volatility during the same period for the natural gas future index. However, the stochastic volatility of natural gas remains moderate compared to the other two series. Our results are consistent with the theory of expectations in the futures market which suggests that speculators considered as economic agents who are able to predict future price trends [START_REF] Kaldor | The role of commodity prices in economic recovery[END_REF]).

Impulse responses

The charts show that the Goldman Sachs indices for the two commodities react during the first two days after the shock. During the first period namely before WHO announces a global pandemic, we observe an insignificant decrease in the oil index and a significant decrease of only 5 units in a natural gas index. In the Second period which is featured by a stock market crash due to COVID-19, we notice a slight increase in the two commodities indices. This result could be explained as follows: during periods of the financial crisis when stock market prices , speculators move their investments from the equity market to the commodity market. Indeed, investors tend to move towards the commodity futures market since it is considered as a risky market with higher expected profitability compared to the equity market. In the third period, the figure shows a significant decrease of two indices, especially, the oil index which registered a decline of about 150 units. This result is mostly linked to the spread of the virus in the USA. During this health crisis, the USA suffers from a high decrease in dollars' exchange rates and an economic recession. Commodity futures indices are not only linked to the imbalance between supply and demand of oil but also the exchange rate of the dollar, the state of the global economic system and equity market.

Conclusion

Using a TVP-SVAR with a stochastic volatility model, we investigate the impact of COVID-19 on the commodity futures market. First, we show the importance of investors expectations in future commodity price dynamics. Speculators in the futures markets anticipate price movements which caused high stochastic volatility during the period that precedes the huge increase in the number of deaths. Then, we show that commodity responses to COVID-19 shock varied from one period to another. The determinants likely to explain the dynamics of the SP GSCI Crude oil index and SP GSCI natural gas index face the COVID-19 include fundamental factors as well as financial factors. In fact, the fall in oil futures prices is partly due to a structural imbalance between supply and demand, particularly a strong decreasing of energy demand following the COVID 19 in large emerging countries such as China and an oil shock caused by disagreements between the Organization of the Petroleum Exporting Countries (OPEC) and Russia, triggering the drop in oil spot prices. In addition,the spread of the virus in the USA impacted negatively the commodity futures market. On the other side, based on behavioral finance theorist, the commodity indices fluctuations are linked to the investors' expectations. In fact, during the crisis , speculators operating on a purely financial logic, turned to the commodity markets to increase their returns which lead to a slight increase in commodity futures index during stock market crash. The short-term outlook for the commodity futures market will depend on the efficiency of government policy to contain the COVID-19 outbreak, and to weaken the impact of the global health crisis has on economic activity.

This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn.com/abstract=3587740

This preprint research paper has not been peer reviewed. Electronic copy available at: