Christophe Macabiau

Lina Deambrogio

Géraldine Artaud

Paul Thevenon

Kalman Filter Based Robust GNSS Signal Tracking Algorithm in Presence of Ionospheric Scintillations

SME (Toulouse) involved in various Radionavigation projects: GPS/EGNOS/Galileo signal processing and critical applications of satellite Radionavigation.

Valentin BARREAU received his electronical engineering degree from ENAC in Toulouse. To validate his degree, he has been doing an internship on GPS occultations at UNB in Canada. He then was working as an engineer at the ISAE, Toulouse. The topic of research was weak GPS signals tracking

more frequently. They also happen over Polar Regions and are associated with the penetration of charged particles along the magnetic field lines. Under the auroral oval region, scintillations can be local or spread over a great part of the oval. According to past observations scintillations occur mainly at the period of the equinox and the solstice. During the rest of the time there are almost no observations. After sunset, the modification in the ionosphere layers generates such variations. Scintillations cause very brutal and fast fades of the received signal, and once these phenomena occur they can last from 30 minutes up to several hours. However, continuous GNSS Carrier Phase Measurements are important observations needed for data demodulation and are increasingly used in GNSS receivers, both at user side (e.g. for precise positioning), and for ground segments to compute the navigation and integrity data. The carrier phase is traditionally tracked in the GNSS receivers using PLL, potentially aided by FLL. Carrier tracking loops may be optimized, depending on the application and environment, by selecting the appropriated loop order, integration time, and bandwidth (trade-off between accuracy and robustness). Phase Loops are however known to be less robust than code tracking loops, and the GNSS receivers may thus suffer from phase tracking loss, for example when tracking low C/N0 signals (attenuations), or fast varying signals such as signals affected by scintillations. This strongly impacts the positioning service availability, as well as the capability to demodulate the navigation message data, in situations where ionospheric scintillations affect the received signal. One thus has to implement innovative techniques and receiver architectures to provide robust carrier phase tracking, either by improving or optimizing the classical tracking loops (optimized parameters of the PLL, potentially dynamically) or by defining different architectures such as Kalman filter-based tracking loops (either scalar or vector architectures, taking benefits of the different tracking channels, constellations and frequencies which may be differently affected by the disturbance). Other techniques coming, for example, from the telecommunication domain can also be interesting to estimate the propagation channel parameters. It can also be taken advantage of the improved structure of the modern GNSS signals, providing in particular a pilot signal component. The improvement technique investigated in this paper consists in replacing the conventional phase lock loop filter by a Kalman Filter PLL (KFP), as inspired by the technique proposed by Psiaki et al.. Kalman loop filters provide the optimal filter gain when the statistical levels of uncertainty of the state and observation vectors are well known. So the Kalman filter is continuously adapting the filter bandwidth to the noise level. KFP tracking loop indeed show a better resistance to weak GNSS signal tracking compared to classical loop filters. We were therefore interested in analyzing the potential of such KFP variants for tracking the GPS L1 carrier phase in presence of scintillations, adapting in particular the feedback control signal. This was conducted during a project financed by CNES. The aim of this paper is therefore to present the development of a GPS L1 phase tracking technique based on a Kalman Filter improving the tracking robustness in presence of ionospheric scintillations, and to present results of its performance using simulations. The paper starts with a description of the phenomenon of ionospheric scintillation, including the possible models for signals affected by ionospheric scintillations, focusing on the model selected for this study which is GISM (Global Ionospheric Scintillation Model) developed by IEEA. The second section presents a review of the state of

INTRODUCTION I. BACKGROUND AND SIGNAL MODEL

Scintillation is a well-known phenomenon when observing stars in a summer night, imaging amplitude and phase fluctuations of electromagnetic waves. The cause can be a diffraction when the waves cross a medium composed of irregular structures with variable refraction indexes. An incident plane wave with a uniform phase gets out the medium with a phase which is no longer uniform in space. The recombination of the resulting fields after propagation can be constructive or destructive and can therefore increase or decrease the wave intensity level. This is illustrated in figure 1. The cause can also be a refraction effect when a wave crosses a medium with a phase velocity larger or lower. The wave keeps its plane structure but the signal phase will be the integration of the phase shifts encountered. Ionospheric scintillations concern signal fluctuations in VHF and UHF radio band up to C band which propagate through the ionized atmospheric layer between 100 and 1000 km, but most particularly in the F2 layer between 250 and 600 km, when it is affected by heterogeneities. The intensity of the effects depends on the signal frequency compared to the plasma frequency and the level of irregularity. In the L band, GNSS is therefore affected too. Several studies showed some effects mainly located along the magnetic equator and at high latitudes in the auroral oval (area of convergence of the Earth magnetic field lines) and in the Polar Cap. These studies also showed a dependency with the solar cycle, as the maximum observed intensities coincide with the maximum UV radiated by the sun and the maximum electronic density in the ionosphere [START_REF] Bishop | [END_REF].

Figure 1: Illustration of diffraction effect from moving ionospheric irregularities leading to ionospheric scintilliation

The physical phenomena that lead to scintillation effects are not well understood yet and are the subject of researches. It is however established that the mechanisms that intervene in the polar region and in the equatorial region have a different nature.

A. Equatorial scintillations

In the case of the equatorial ionosphere, it is largely admitted that the instability mechanism of Rayleigh-Taylor is the main cause of the ionization irregularities that develop after the sunset. In its most frequent form, this phenomenon corresponds to the situation encountered when two non-miscible liquids are in contact, the densest liquid being above the least dense liquid. The equilibrium between the two liquids is unstable and any perturbation causes the least dense liquid to go up and reach the surface. Then, at the interface between the two liquids, some bubbles appear that draw the least dense liquid upwards. In the ionospheric plasma, such bubbles with weak ionization level appear frequently after sunset in the lower part of the ionosphere (typically at altitudes between 200 and 400 km) then go up at higher altitudes and take aigrettes shapes [START_REF] Duchesne | Les effets de la scintillation ionosphérique sur le GPS[END_REF]. Equatorial scintillations exhibit seasonal variations that superimpose to the modulation of the effects with the position in the solar cycle. Observations in America, Africa and India show the occurrence of peaks at the equinoxes and the quasi total absence of scintillations during summer. In the Pacific, the seasonal influence is inversed with observed maximums from April to August [DOHERTY et al., 2000].

B. High latitude Scintillations

Scintillations at high latitudes are associated with large spatial extension plasmas typically of the order of 100-1000 km in the F region. They are more intense in winter when the solar radiations do not smooth the ionization irregularities. For the high latitudes, we distinguish the region of the polar cornet in which the magnetic field lines are open to the solar wind and the auroral oval, band centered on the geomagnetic pole between 60 and 75 degrees of latitude. The auroral oval is the area of concentration of shining auroras, emissions linked to particles precipitations coming from the equatorial plasma leaflet, region of the queue of the magnetosphere and carried by aligned electric currents. The oval is run by electric currents (electrojets) that move towards the west. The electrons precipitations happen at an altitude of about 100 km, mainly on the night side along the field lines during substorms,sporadic phenomena of a few hours whose origin is in the magnetosphere or during intense solar eruptions ejecting very energetic charged particles (CME : Coronal Mass Ejection). In case of a magnetic storm (Kp>5), the solar wind puts a pressure on the magnetosphere triggering sub-storms but also by extending the limits of the oval towards the south, explaining why boreal auroras were sighted as far in the South as in French Normandy (oct./nov. 2003). Practically, the polar cornet is an open magnetic region and phase and amplitude scintillations can be met. The aurora oval is a magnetic region where lines are closed, phase scintillations are predominant.

C. Scintillations at mid-latitudes

After an intense magnetic storm and in particular if its origin is a solar eruption with energetic particles, eg. during coronal mass ejection, the engendered ionospheric perturbation can migrate towards the south.

D. Ionospheric scintillations and effects on GNSS signals

The ionosphere causes a group delay of the modulation and a phase advance. Ionosphere irregularities that cause GNSS scintillations can affect the GNSS signal through refraction and diffraction causing rapid variations of the group delay and phase shift of the carrier [START_REF] Kintner | [END_REF] It can be modeled as the square root of the sum of three terms [START_REF] Conker | [END_REF]:


is the contribution due to thermal noise and can be derived as:

[ ( ) ] ( ) Where:  is the PLL loop bandwidth  is the integration time   is the phase scintillation variance ̃ ∫ | | ̂
Where ̂ is the Power Spectral Density (PSD) of the phase scintillation and | | is the closed loop transfer function of the PLL that depends on k, the loop order, and , the loop natural frequency.

| |

Typical values are k = 3 and = 1.91 Hz [BENIGUEL et al., 2004].

 And is the contribution due to the receiver oscillator noise. These indices can be computed over time, and their estimates are generally based on averages obtained over a few tens of seconds. In order to reflect the effect of the scintillation only, it is usually required to remove low frequency variations due to the change of the distance between the satellite and the receiver and local oscillator impact. Some specific GPS receivers called Ionospheric Scintillation Monitors (ISM) have been developed to remain locked and provide continuous estimates of and , typically every minute. b. GISM The GISM model (Global Ionospheric Scintillation Model) enables to compute the scintillation effects on a signal after the crossing of the ionosphere between a fixed terrestrial receiver and a GPS satellite. It is a climatological model, statistical properties are associated to the medium to take into account in particular the known characteristics of its spectrum, the effect of the solar cycle and of local time. It uses the MPS technique (Multiple Phase Screen) to solve the wave propagation equations and to determine the expected fluctuations. The medium is divided in successive layers, each acting as a phase screen. The line of sight is first computed considering the different indices of the layers crossed and thus the incident and refracted angles, the electron density is given by the model at each layer by the NeQuick model [START_REF] Radicella | The NeQuick model genesis, uses and evolution[END_REF]. The fluctuations and the generated errors are computed by solving the Haselgrove equations. The signal spectrum is characterized by 3 parameters: the slope, the cut-off frequency and the lowest frequency. The slope is comprised between -2 and -3 with -2 as the default value. A detailed model is given in [BENIGUEL, 2010]. The model inputs are:

 The date of the simulations  The solar flux  The geographical coordinates at the observation point

The main outputs of the model are:

 The time series of the amplitude and phase of the signals at the selected frequency.  The estimated S4 and sigma-phi computed every minute.

c. WBMOD

WBMOD is a model developed by NWRA with the support of the US government to estimate a severity level of the effects of the scintillations for a system and of particular environmental conditions. It is a model of the ionosphere that provides the global distribution, a behavior of the electron content irregularities at the origin of scintillations and the system impacts. A scintillation database (Wideband, HiLat and Polar Bear experiences, USAF monitoring network) feeds an empirical model (EDIM) of the irregularities of electron density as the function of date and time, the localization but also the solar activity parameters (number of sunspots) and geomagnetic (Kp). WBMOD allows to estimate the severity of the effects of the scintillations at a given altitude and plots the corresponding map. It does not provide the classical and parameters. It is not adapted for a short term forecast. Like GISM, WBMOD is also a climatological model but while GISM estimates the classical scintillation parameters and provides time series, WBMOD estimates probabilities that a scintillation level is reached. The calibration of the model was done with wideband observations at 447 MHz, in UHF and at 1239 MHz. A detailed model is available in [START_REF] Kintner | [END_REF]. The model inputs are the level (with a maximum of 1) and a correlation time called (from 0.5 s to 2s for GPS frequency). The model allows for example to generate time series for a GPS SPIRENT signal simulator. Based on these simulations, performances of robustness of receiver PLL were obtained.

e. Models limitations

The models described above are only applicable in the equatorial region. The phase screen model is insufficient because a distinction is needed between the 2 affected regions, the polar cornet and the auroral oval, each with different mechanisms. Until recently, no prediction model existed for prediction of scintillations for high latitudes. A few recent publications address a tentative of a polar and equatorial model called WAM (Wernick, Alfonsi and Materassi).

For GISM, another limitation of the scintillation model is that it does not reflect reality, which is the presence of low electron density irregularities concentrated in space but a medium with layers with identical characteristics at a given time all around the observation point. Thus, there is no representativity of the fraction of the satellite constellation affected by the phenomenon; all satellites are affected in a homogeneous way. For Cornell model, its application is also limited to equatorial regions. The other limitations are a maximum value of S 4 of 1 although much higher values can be observed. Also, and are assumed constant during a simulation run, thus generating a stationary pattern. Some authors also consider that this model is susceptible to present biases in GPS frequencies from the fact that some calibration parameters were established in wide band [START_REF] Cigala | Cigala Project, WP200, State of art review[END_REF] .

II. REVIEW OF SCINTILLATION TRACKING TECHNIQUES

In this section, only the techniques at signal processing level will be considered. The aim of the techniques will thus be to improve the robustness of the carrier phase tracking, in order to increase the capability of the loops to follow accurately the phase variations; the objective will be to decrease the number of phase tracking losses and to limit the impact on the provided measurements (including data demodulation). This kind of techniques can improve the tracking capability of each signal (receiver channel) without external aiding, or can take benefits of the information coming from other channels tracking (mono-or multiconstellation, on single or dual frequency), or can finally be aided by external sensors information on the receiver dynamic (Doppler aiding by inertial sensors for example).

A. Conventional Loops

Conventional PLL (Phase Lock Loop) are based on close loop filters, which uses discriminator outputs to perform signal tracking. Those conventional architectures have limited performance for tracking phase measurements in difficult conditions, however, smart adaptations of these tracking loops is a good start to achieve good performance in presence of fading, dynamics and scintillations. Several parameters, such as filter order and loop bandwidth but also the discriminator type, affect the signal tracking behavior of the filter. The higher the order of the filter is, the better the estimation of the residual dynamics will be, but the higher the computational burden will be. Moreover, the loop bandwidth is an important parameter of the loop filter. On one hand, a low loop bandwidth decreases the thermal noise but also slows down the capacity of the filter to follow brutal change of dynamic. On the other hand, a high loop bandwidth increases the thermal noise but the loop filter is able to follow more precisely the changing of dynamic. The predetection time is another important parameter. It corresponds to the integration time of the signal in the correlators. The longer the signal is integrated, the lower the thermal noise will be. But the dynamic variations have to be followed during the integration process in order not to introduce some dynamic tracking error. A FLL (Frequency Lock Loop) can also be used and will be more robust than the PLL as it does not control the exact value of the phase of the local carrier, however this loop will not allow the navigation message demodulation. The FPLL (FLL assisted PLL) combines advantages of both PLL and FLL. This loop filter uses two discriminator inputs, one frequency error and one phase error inputs. This loop is more robust than a classical PLL because of the frequency aid. If the phase lock is lost, the FPLL can keep tracking the signal thanks to the frequency aid contrary to the PLL and then restart the phase tracking as soon as the phase is locked again. All these types of loop are still scalar loops, working in parallel on each channel with a closed-loop filter and having traditionally a constant loop bandwidth. The fine tuning of the loops parameters may increase the robustness of carrier phase tracking but a trade-off is still to be done between accuracy and robustness. Finally, the closed loops architecture may be improved by reducing the dynamic range of the PLL. This can be done by ultra-tight hybridized receiver architectures, which use some external receiver dynamic information to help the tracking loop filters. By aiding the loop with receiver dynamic information such as speed and acceleration or Doppler, the filter loop bandwidth can be reduced as well as the noise level with an improved robustness to dynamics. However, this solution is limited in case of scintillation because the phase variations (due to scintillation) can be strong and not predicted by external sensors. It will finally be of high interest for high dynamics applications.

B. Adaptive Bandwidth Loops

On top of the classical loops, "Fast Adaptive Bandwidth" (FAB) technique [LEGRAND, 2002] has been developed to minimize the total tracking error in presence of noise and dynamics. This technique is interesting because it allows the receiver to track the signal with a loop bandwidth adapted to the dynamics and the noise level. In case of signal fading, the C/N0 can be very low and at the same time the phase variations very fast (high dynamic or scintillation). That means that the loop bandwidth should increase because of the high dynamic imposed by the scintillation and also decrease because C/N0 is rather low. Because of that, FAB is not really effective during scintillation or high dynamics operation with low C/N0 [START_REF] Kondo | Evaluation of tracking performance in the presence of ionosphere scintillation on GPS signal in Japan[END_REF], which is the case most of the time, although it proposes the optimal tuning of the loop bandwidth in this kind of situations.

As the loop is characterized by its transfer function, the discriminator output signal can be predicted and analyzed to perform dynamics and noise estimations. Thus, these estimations are used to build an optimization function that will be minimized by an iterative gradient method in order to provide the corresponding optimal loop filter setting.

C. Kalman Filter-Based Tracking Loops

Other improvement techniques will consist in replacing the conventional loop filter by a Kalman loop filter. Kalman loop filter provides in theory the optimal filter gain when the statistical level of noise of the inputs is well known. In such cases the Kalman filter is continuously adapting the filter bandwidth to the noise level. KFP tracking loop indeed showed a better resistance to weak GNSS signal tracking compared to classical loop filters [ZHANG et al., 2010]. Two options can be considered:

 A 3-state (phase error, Doppler frequency, and Doppler frequency rate) 2-equations (prediction and update) Kalman filter is typically used with a constant gain [START_REF] Psiaki | Tracking L1 C/A and L2C signals through ionospheric scintillations[END_REF], [ZHANG et al., 2009]. In this case the loop bandwidth remains constant but the Kalman filter will allow to track at lower C/N0 and to reduce the cycle slips and loss of lock compared to conventional PLL.  A 5-equations Kalman filter (including the prediction steps) would allow adaptation of the loop bandwidth. However, this requires to estimate the statistical variance of the input model (state noise and observation noise), which is usually difficult to obtain. In addition, a constant-gain Kalman is more stable than a gainvarying Kalman. Indeed, Kalman filtering can only process parameters subject to gaussian noises, so fixing the gain can avoid filter divergences because of non-gaussian noise effects. [ZHANG et al., 2010] proposes a Kalman filter architecture with Kalman gain adaptation according to the C/N0 level, which limits the drift of the phase error. The adaptive loop bandwidth combines a weighted Kalman gain with a weighted discriminator estimate to handle the difficult scenario when the filter has no loop bandwidth working range because, for example, of scintillation.

D. Vectorized Architectures

The above solutions are all scalar architectures with a filtered feedback. An alternative is the use of vectorized architecture, which merges observations from all the channels to use most of the available useful information while averaging the tracking errors on each channel and to improve robustness of tracking with aid from all the individual channels tracking. They are usually closed loop architectures but with only one global state vector used in turn to compute the feedback on each channel filter. NCO's control signals on each channel are thus computed from the residuals of position, velocity and clock offset of the receiver instead of using only code and phase discriminator outputs on each channel. This kind of architecture allows reducing the impact of jamming and interferences, allows tracking at lower C/N0, and offers a greater immunity to receiver's dynamic. Moreover, in case of signal outage on one of the channel, the receiver keeps following signal dynamic and so allows the tracking function to restart as soon as sufficient signal power comes back. In particular, Vectorized Delay Frequency Lock Loop (VDFLL) uses information of position and speed of the receiver to compute the Doppler frequency and the code phase of each satellite. The pooling of dynamic information of all the tracking channels allows the receiver to keep tracking satellite signals even when perturbations disturb some of the channels. Indeed, it is possible to rebuild the Line-Of-Sight (LOS) distance and the Doppler for each satellite from the position and the speed of the receiver, and so to compute values to send to the NCO of each tracking channel.

These vectorized architectures work well as long as a good receiver's position/velocity is estimated. But as soon as the position estimate is wrong, the position error will be introduced in all the tracking channels and eventually could cause the loss of all the channels. Some integrity position monitoring could be realized to ensure the use of good position and velocity information or the use of good pseudorange information. An alternative, the Vectorized Phase Lock Loop (VPLL) has been described in [START_REF] Zhodzishsky | [END_REF]]. The previous concept of vectorized DLL cannot be directly applied to PLL, due to the inherent ambiguity of the carrier phase and its short wavelength. Indeed a direct transformation between the phase error and the estimated position error is much more difficult to find and the effect of atmospheric errors and dynamic would greatly influence the tracking. This technique is efficient in case of short signal blockage. The PLL will continue the tracking as long as there are enough satellites to compute a precise evolution of the receiver coordinates and clock. The concept has been pushed further in [START_REF] Henkel | [END_REF]. The goal here is to merge data for the same satellite but from different frequencies to determine the dynamic of each parameter and then to remove the satellite-related parameters (ionospheric, tropospheric errors, satellite clock). The receiver-related information (position error, receiver clock) are merged and back transformed to a phase variation which is then integrated separately on each channel. However, for both techniques, the merging of all channel observations also has some drawbacks, since some perturbations on one of the channel may be reported on all the channels. This will possibly lead to receiver instability or loss of lock on all satellites. The use of integrity monitoring can help removing some perturbed signal from the solution. In addition, the major drawback is the high processing load and the complexity of the receiver. The Kalman filter or the Least Squares estimation technique used by the vector tracking architecture must be iterated on a time scale commensurate with the integrate-and-dump period used by the algorithm (~ 50 Hz) [PETOVELLO, 2009]. [START_REF] Henkel | [END_REF] proved that vectorized phase tracking is working, but at least dual frequency signals have to be processed to remove the maximum of perturbation errors (ionosphere, clocks). This technique is thus very constraining.

E. Other techniques

Other techniques can be used, inherit from other domains such as telecommunication. For example, wavelet de-noising techniques can help to reduce the noise on tracking channels, based on wavelet transform algorithms. Wavelet transform tends to concentrate the signal energy into a relatively small number of large coefficients. On this basis, a method called wavelet shrinkage to use threshold in wavelet domain was proposed, and it was shown to be asymptotically near optimal for a wide range of signals corrupted by additive Gaussian noise [ZHANG, 2000]. It thus can be useful in case of fading or scintillation to better differentiate noise from dynamic residuals. This technique commonly used for video image de-noising, has been adapted to GNSS field, notably for the acquisition process [START_REF] Tian | [END_REF] or for precise point positioning to remove residual errors [WANG et al., 2009]. The results shown by this technique are interesting in terms of performance but the processing load and complexity is rather high and real time implementation may be compromised. The last type of architecture, very different from the others, is the open loop architecture. Actually, open loop estimators are commonly used in telecommunication (cell phones) in deep urban situations, for example to maintain uplink signal power constant. They have the advantages to be more robust than a closed loop filter since the output is not fed back on the entry of the filter (the errors cannot propagate epoch-to-epoch). Such estimators are used with estimation techniques like Maximum Likelihood Estimation (MLE). Another technique is called Maximum A-Posteriori Estimation (MAP), which uses a priori information about the parameter to evaluate. In the GNSS world, the closed-loop approach has historically been used in GPS receivers [START_REF] Kondo | Evaluation of tracking performance in the presence of ionosphere scintillation on GPS signal in Japan[END_REF]. In this type of receivers, acquisition and tracking are separated. Results from acquisition are used to initialize the tracking. The main advantage of the closedloop sequential processing is to reduce receiver complexity as only three correlators are required to track the signal. An open-loop batch processing technique is proposed in [ [START_REF] Van Graas | [END_REF], which does not separate acquisition and tracking. A batch of GPS signal is processed once and the signal parameters (code phase, Doppler and carrier phase), are estimated independently from one batch to another. Finally, this technique offers higher receiver sensitivity because higher integration time can be used. Nevertheless, a Doppler and Doppler rate aiding has to be provided to the receiver to allow a longer integration of the signal. It also reduces the frequency search domain. Another advantage is the immediate re-acquisition when a GNSS signal becomes available. Also, this type of open-loop processing can include the joint estimation of the characteristics of the propagation channel, such as the characteristics of the multipath in case of a deep fading scenario, and specialized techniques can be implemented for this. The major drawback is the large complexity of the receiver and the wide frequency search domain in case no Doppler aid is provided.

F. Synthesis

The techniques analyzed in the previous sections allow a better resistance to noise, co-operation of the channels with each other, loop aiding, and loop bandwidth adaptation. Each technique has advantages and drawbacks that must be considered.

III. ANALYZED KF TRACKING TECHNIQUE

The KFP algorithm (Kalman Filter Phase lock loop) is based on the models and assumptions presented in [START_REF] Psiaki | Tracking L1 C/A and L2C signals through ionospheric scintillations[END_REF]. This algorithm enables to implement a carrier tracking loop of order 3. The signal processed by the receiver is the received GPS L1 C/A signal. The duration of one bit of the navigation message is thus 20 ms. Let us denote the time samples at which the correlator outputs are provided. We define the duration of the coherent integrations. We define K the number of correlator outputs available within the duration of one bit of the navigation message. We define and the correlator outputs available at time , i.e. the correlator outputs computed during the interval . We define the value of the NCO control signal during the interval . The phase of the PLL NCO is computed as .

We define [ ] the state vector of the Kalman Filter, with  the difference between the phase of the incoming carrier and the phase of the local carrier generated by the NCO



the difference between the doppler shift of the incoming carrier and the Doppler shift of the local carrier generated by the NCO in rad/s (pulsation)



the difference between the jerk of the incoming carrier and the jerk of the local carrier generated by the NCO This algorithm only includes 2 equations of the Kalman Filter among the 5 traditional equations. Indeed, the Kalman gain K is proposed to be fixed, therefore the two covariance equations, during the prediction and during the update, are useless and are not considered. The algorithm also includes an iterative algorithm allowing deciding the value of the bit of the navigation message affecting the correlator outputs et in order to provide the extended arctangent discriminator outputs. The 2 equations of the Kalman Filter to be run every are:

[ ̂ | ̂ | ̂ | ] [ ] [ ̂ | ̂ | ̂ | ] [ ] [ ̂ | ̂ | ̂ | ] [ ̂ | ̂ | ̂ | ] [ round ]
The authors of [START_REF] Psiaki | Tracking L1 C/A and L2C signals through ionospheric scintillations[END_REF] propose that the Kalman gain is not updated, but considered fixed and

defined such that the matrix [ ] [ ] has the eigenvector [ [ √ ] [ √ ]
].

We obtain the values of the gain K, as a function of the desired loop bandwidth and integration time: The measurement prediction is:


̂ | [ ] [ ̂ | ̂ | ̂ |
] .

The considered measurement vector is computed from the discriminator output. If we want to use an ATAN2 discriminator, it is necessary to estimate the value of the data bit. In this case, the measurement is obtained as

atan ( ̂ | ̂ | ).
This discriminator output uses a prediction of the value of the bit of the navigation message under the form:

̂ | { si ̂ | ̂ | si ̂ | ̂ | ∑ floor , ∑ floor( ) And ̂ | ∑ ( ̂ | ) floor , ̂ | ∑ ( ̂ | ) floor With ̂ | [ ] [ ̂ | ̂ | ̂ | ]
Note that we can envisage another technique of prediction of the navigation message bits which is:

̂ | { si ∑ [ ( ̂ | ) ( ̂ | )] floor si ∑ [ ( ̂ | ) ( ̂ | )]
floor If the measurement uses the output of an atan discriminator, the estimation of the data bit is then useless and the measurement is directly atan .

Note that the authors of [START_REF] Psiaki | Tracking L1 C/A and L2C signals through ionospheric scintillations[END_REF] propose to compute the NCO control signal in the following way:

[ ( ̂ | ) ( ̂ | ) ̂ | ] ̂ | ̂ | Where:  is such as ̂ | ̂ | ̂ | 
is a tuning factor that can be chosen such as , or very close to 1 if we want to limit the velocity of convergence of towards Another simpler feedback can also be used. This feedback does not use the parameter « etha » which enables to smooth the rapid variations of the phase:

̂ | ̂ | ̂ |
In the static case, the smoothing feedback using is used. This feedback depends on « etha » a smoothing parameter which enables to limit the reactivity of the feedback and thus stabilizes the loop filter. The tracking configuration of the improved receiver is presented in table 1. In the dynamic case, the feedback is used. Contrary to the static case, parameter « etha » is not used, because the dynamics of the receiver is too difficult to follow. This feedback allows a faster reaction to the strong variations of the control signal contrary to the static case as the smoothing tuned by « etha » is not used.

Nominal Scintillations

Intense Scintillations

Worst case Scintillations

Discri

For nominal and intense scintillations, the discriminator used depends on the estimated value of C/N 0 and on the navigation data bit. Above 35dBHz, atan2 is used but only if the two methods for data bit estimation give the same result. In the contrary, atan is used.

Nominal Scintillations

Intense Scintillations

Worst case Scintillations

Discri

IV. SIMULATIONS DESCRIPTION

This section describes the simulation environment used to model the scintillation phenomena, the GNSS signal processing modules (signal generation and reception architectures) and the classical GNSS signal tracking algorithm.

The simulator selected to run the simulations is the JGNSS simulator, owned by CNES. It models a GPS L1 C/A processing chain including: bit level GPS L1 C/A signal generation, propagation channel model including scintillation and multipath, GPS L1 C/A tracking. Position calculation module is not run in this study that focuses on the robustness of the tracking loops. Signal duration of 100 seconds was selected for the simulations. This duration is sufficiently long for the GNSS tracking loops to be disturbed by the scintillations and for the simulation time to remain acceptable (about 12h).

The architecture of the standard receiver is classical. It is based on FLL aided PLLs pushing the DLL using correlator outputs. The correlators and the loops are set with parameters having classical values (Bl=1 Hz for the DLL and Bl=10Hz for the FLL and PLL, Early minus Late chip-spacing =0.25). Integration time of 20 ms is used.

The RF bandwidth will be assumed to be different between the static case (ground station) and dynamic case (aircraft in approach), for the simulations with the standard receiver or for the simulations with the enhanced receiver.

Two cases are considered:  Static receiver (ground station)  Dynamic receiver in an aircraft in approach. A straight line descent with a constant velocity of 70 m/s is considered.

As we can see, the improved receiver has improvement features compared to the classical receiver:

 In the static case, the improved receiver uses an atan2 phase discriminator with known data bits and a Kalman filter. The integration time is 10ms for the improved receiver compared to the classical receiver.

 In the dynamic case, the improved receiver uses an atan2 phase discriminator when C/N 0 >35dBHz, and a Kalman filter. When the C/N 0 is low, both the classical and the improved receivers use an atan phase discriminator. The multipath model to be considered is described below.

Ground station

In the static case, because the observation interval signal is lower than the period of repetition of the GPS constellation ground tracks, although the geometry between the ground obstacles and the ground station is fixed, the periodicity of the errors due to multipath is not taken into account.

The model proposed is the combination of two submodels: the first part is a model of the parameters of the multipath, the second part is a model of the impact of this multipath on the pseudorange measurements.

The model of the parameters of the multipath is simplified mainly by drawing randomly these parameters. We propose to consider simultaneously 2 reflected rays at each epoch. The 2 reflected rays considered are a ray reflected by the ground underneath the antenna, and a ray reflected by a random obstacle located around the antenna.

In Furthermore, two values of C/N0 will be considered for the undisturbed direct path: a C/N0 of 45dBHz and a C/N0 of 35dBHz.

A. GISM simulated data GISM was selected for this study since it allows to generated amplitude and phase time series of the signal at antenna output. However, we know that the model is limited to the equatorial region and provides an identical effect for all satellites of the constellationexcepted the geometry effectsbecause it does not consider the presence of bubbles of electronic plasma. CLS bought the GISM software from the company IEEA. This software provides the times series of the phase and amplitude fluctuations of the signals with a 50 Hz update rate as well as the scintillations parameters S 4 and sigmaphi for the different scenarios considered. These samples are interpolated in a linear way to bring them to the simulator sampling frequency of 100MHz. The following figures illustrate the phase and amplitude variations of the signal for the scintillation scenarios considered. We can see that power variations are larger compared to the normal scintillation case, except for the absence of a very deep C/N0 peak like in the nominal case. These data correspond to a solar flux of 84. We did not choose to consider the maximum case in this situation, so the variations considered are average variations for the intense scintillation case. We can see that Doppler frequency variations generated by scintillations lie within -800 and +800Hz, meaning that the phase can evolve by 1600*0.02= 32 cycles in 0.02 seconds which is the integration time considered for the standard receiver. Such a Doppler variation is extreme and results will not be good in such cases.

V. SIMULATIONS RESULTS

We compare here the simulation results obtained with the standard receiver with the results obtained with the improved receiver in the static and dynamic cases.

The This table 4 presents the better carrier tracking performance of the KFP compared to conventional tracking. In the nominal scintillation scenario, at 35dBHz, the standard deviation of the phase tracking error is 1.10cm for the KFP while it is 1.76cm for the FPLL. At 45dBHz, this sigma is 1cm for KFP while it is 1.37cm for FPLL.

In the worst scintillation case, tracking with the FPLL leads to too many reacquisitions over the 100s simulation time to be able to estimate a valid phase tracking error standard deviation.

In this worst case, the use of KFP with a narrow bandwidth (3Hz) and a large etha (0.96) instead of conventional FPLL allowed decreasing the number of reacquisitions compared to the FPLL case, and a phase tracking error standard deviation could be estimated. However, the large value of these standard deviations at 35 and 45 dBHz show that the phase is not properly tracked. This tuning has to be adopted because phase variations can reach 1.6 cycle/ms, which is too fast to be followed through 1 ms integration as the discriminator has an operating region of -0.5…+0.5 cycles.

No significant change of the code tracking error is observed if we compare the nominal and intense simulations cases. At 35dBHz, the standard deviation is around 2.8m, and it is 90cm at 45dBHz, whatever the receiver type.

In the worst simulation case, these values of standard deviation increase a lot (2.33m at 45dBHz, 8.63m at 35dBHz for the improved configuration). Indeed, the carrier Doppler and phase being badly followed, correlator outputs are noisier, and all discriminators are affected. Also, as the DLL is pushed by the carrier tracking, the DLL is affected.

The number of re-acquisitions is highly reduced by using the KFP, except in the nominal case where the singular power loss of 23dB cannot be overcome. The loss of lock rate is also improved.

No cycle slip is observed for the nominal and intense cases, whatever the tracking configuration tested. Numerous cycle slips are observed in the worst simulation case but this is due to the large number of reacquisitions.

The table below presents the performance for the dynamic receiver, with the classical loop filter replaced by an equivalent 3 rd order Kalman filter. In this situation, both carrier tracking techniques use an atan discriminator at low C/N0, which has the largest influence on the results. The estimated standard deviations of the phase tracking errors are lower for the KFP than for the FPLL.

Sci

A slight increase of the standard deviation of the code tracking error for the KFP is observed if we compare nominal and intense scintillations. This is thought to be due to the large number of re-acquisitions.

In the worst scintillation case, the standard deviation of the phase tracking error is again large, showing that the phase is not tracked accurately.

As for the static case, the number of re-acquisitions is again low compared to the conventional tracking configuration, except for the nominal scintillation case. It can also be noted that at 35dBHz for the nominal and intense scintillation, even though the number of reacquisitions is lower when using the KFP, the loss of lock rate did not decrease. This is due to the fact that the receiver takes more time to achieve bit synchronization after a re-acquisition and stays in a DLL + FLL mode.

In the worst scintillation case, the rate of loss of lock at 45dBHz goes from 37% for the classical dynamic receiver to 27% for the improved dynamic receiver.

The receiver dynamics sometimes creates problems for reacquisitions and can lead to false peak detections, leading in turn to new re-acquisitions until a meaningful peak is found.

CONCLUSION AND FUTURE WORKS

In this study, it was observed that standard tracking was quite robust for nominal scintillations, but often lost lock in more severe situations. We analyzed an improved tracking technique using a Kalman Filter PLL (KFP) algorithm and proposed enhancements to this technique. We observed that the KFP improved carrier tracking in all situations. However, in our tests, KFP did not improve code tracking.

In the future stages of the study, we can propose the following options for further improvement. First of all, it could be tested a way to adapt the loop bandwidth as a function of detected scintillation intensity. Acceleration input could also be used for tracking loops. It could also be analyzed a way to replace the classical DLL by more robust DLL to improve tracking during strong perturbations. Finally, a frequency discriminator could be added as an additional Kalman filter observable input to maintain carrier tracking even if the phase is not tracked.

  Scintillation Model) is a model for perturbations of scintillations in amplitude and phase developed by Cornell University. It is only applicable in the equatorial region. It is a statistical model based on Rice distribution. The disturbed signal component dominated by scintillation uses an autocorrelation function.

  the dynamic case, the high resolution aeronautical model developed for the European Space Agency (ESA) in 2002 by Joanneum Research, University of Vigo and the DLR, based on measurements is used Erreur ! Source du renvoi introuvable.. This model includes a direct (Path 0), a refracted component (Path 1), a strong fuselage echo changing very slowly (Path 2), and a ground echo changing very rapidly (Path3). The outputs of this model are the multipath parameters: relative delay, relative amplitude and relative carrier phase shift. The values of the parameters of the 3 refracted and diffracted rays can then be used to evaluate the carrier and code tracking errors. Three scintillation scenarios were selected:  A nominal scintillation scenario: it occurred on 23/10/2010, with a solar flux of 78, which is rather favorable. The is set to 0.3.  An intense scintillation scenario: it occurred on 01/12/2006, with a solar flux of 84. The is set to 0.5.  A worst case scintillation scenario: it occurred on october 2003. The solar flux reached a maximum of 279. The is set to 0.9.

Figure 2 :Figure 3 :

 23 Figure 2: Variations of the received signal Doppler offset in the nominal scintillation case

Figure 5 :Figure 6 :

 56 Figure 5: Variations of the received signal Doppler offset in the intense scintillation case

Figure 7 :Figure 8 :

 78 Figure 7: Variations of the received signal Doppler offset in the worst scintillation case

Table 1 : Tracking configurations of improved receiver adopted for scintillation in the static case

 1

		atan2	atan2	atan2
		+ known nav	+ known nav	+ known nav
		bits	bits	bits
	Bl (Hz)	10	10	3
	Integration time (ms)	10	10	1
	Etha ( )	0.77	0.77	0.96

Table 2 : Tracking loop configurations of the improved receiver in the dynamic case

 2

	atan2 if	atan2 if
	(C/N 0 >35dBHz and ̂ ̂ )	(C/N 0 >35dBHz and ̂ ̂ )	atan
	atan if	atan if

Table 3 : Receiver parameters

 3

	Aircraft

Table 4 : Performance obtained for static receiver with KFP and FPLL.

 4 table below presents the performance obtained for the static receiver in the standard and improved version (KFP), with the loop filter replaced by a Kalman filter equivalent to a 3 rd order filter. Standard and improved receivers are named respectively C and I in the table. The results obtained are better when a KFP is used than when a FPLL is used.

	Sci	C/N0 (dBHz)	Rec (cm)	(m)	# of c.s.	# of re-acq	LLR (%)
			I	1.10 2.80	0	1	0.60
		35
			C	1.76 2.72	0	1	1.12
	N
			I	0.99 0.89	0	1	0.40
		45
			C	1.37 0.85	0	1	0.51
			I	1.31 2.92	0	1	0.71
		35
			C	1.88 2.55	0	7	3.41
	I
			I	1.25 0.99	0	0	0
		45
			C	1.30 0.91	0	6	2.24
			I 673.56 8.63	45	15 44.68
		35
			C	N.A. N.A. N.A. 55 49.38
	W
			I 693.49 2.33	40	14 37.32
		45
			C	N.A. N.A. N.A. 51 39.87

Table 5 : Table of performance obtained with KFP and FPLL

 5

	C/N0 (dB.Hz)	Rec (cm) (m)	# of c.s.	# of re-acq	LLR (%)
		I	0.91 2.99	0	5 (3) 3.05
	35
		C 2.06 2.89	0	6	2.76
	N
		I	0.47 0.89	0	1	0.41
	45
		C 1.55 0.89	0	1	0.51
		I	1.21 4.2	0	9 (4) 8.47
	35
		C 2.92 2.49	0	12	9.1
	I
		I	0.60 0.95	0	0	0
	45
		C 1.97 0.88	0	6	2.34
		I 813.3 8.56 47	4 (4) 45.2
	35	C N.A. N.A. N.A.	45 (in 80s)	53.3
	W
		I 419.5 3.26 46	9 (7) 27.34
	45	C N.A. N.A. N.A.	41 (in 80s)	37.37

Willy VIGNEAU is head of the


Ahmed Atil
email: ahmed.atil@esc-rennes.fr

Amine Lahiani
email: amine.lahiani@univ-orleans.fr

Khuong Duc

Nguyen
email: duc.nguyen@ipag.fr

Asymmetric and nonlinear pass-through of crude oil prices to gasoline and natural gas prices

Keywords: oil price shocks, gasoline, natural gas, NARDL model, asymmetric pass-through. JEL classification: C58, Q40, Q47

In this article, we use the recently developed nonlinear autoregressive distributed lags (NARDL) model to examine the pass-through of crude oil prices into gasoline and natural gas prices. Our approach allows us to simultaneously test the short-and long-run nonlinearities through positive and negative partial sum decompositions of the predetermined explanatory variables. It also offers the possibility to quantify the respective responses of gasoline and natural gas prices to positive and negative oil price shocks from the asymmetric dynamic multipliers. The obtained results indicate that oil prices affect gasoline prices and natural gas prices in an asymmetric and nonlinear manner, but the price transmission mechanism is not the same. Important policy implications can be learned from the empirical findings.

Introduction

Understanding the pass-through of oil prices to gasoline and natural gas prices is of paramount importance not only for energy policy decision-making but also for portfolio risk management and optimal hedging issues. Indeed, prices of these energy sources influence the level of capital expenditure of almost all economic agents and the incentives to promote alternative energy sources. Investors in energy markets may find the tendency of the prices to adjust to a certain relative price relationship or parity as the basis for a trading strategy. For example, [START_REF] Brown | What drives natural gas prices?[END_REF] document that there was a 10-1 relationship between the price of one barrel of oil and the price of one million BTUs of natural gas in the 1990s but this relative price relationship declined to a 6-1 ratio more recently as the oil and natural gas prices rise.

The empirical literature on the oil price pass-through is extensive and has started with studies examining the relationship between oil prices and the macroeconomy. The study of [START_REF] Hamilton | Oil and the macroeconomy since World War II[END_REF] shows that rising oil prices are responsible for nine out of ten of the U.S. recessions since the Second World War, thus suggesting a negative impact of oil price changes on U.S. output growth. Most subsequent studies reach similar conclusions and attempt to explain the oil-economic growth nexus by the oil price pass-through to production cost, inflation expectation, monetary policy, and investor confidence (e.g., [START_REF] Mork | Oil and the macroeconomy when prices go up and down: An extension of Hamilton's results[END_REF][START_REF] Hooker | What happened to the oil price-macroeconomy relationship[END_REF][START_REF] Cologni | Oil prices, inflation and interest rates in a structural cointegrated VAR model for the G-7 countries[END_REF]. There is also evidence to suggest the existence of asymmetric and nonlinear links between oil price and economic growth, even though they become weaker in recent years (e.g., [START_REF] Hamilton | What is an oil shock?[END_REF][START_REF] Kilian | Exogenous oil supply shocks: how big are they and how much do they matter for the U.S. economy?[END_REF][START_REF] Lardic | Oil prices and economic activity: an asymmetric cointegration approach[END_REF]. Another major stream of research has focused on how oil prices are passed on to stock prices (e.g., [START_REF] Jones | Oil and the stock markets[END_REF][START_REF] Sadorsky | Oil price shocks and stock market activity[END_REF][START_REF] Ciner | Energy shocks and financial markets: nonlinear linkages[END_REF][START_REF] Arouri | Oil prices, stock markets and portfolio investment: evidence from sector analysis in Europe over the last decade[END_REF][START_REF] Aloui | Assessing the impacts of oil price fluctuations on stock returns in emerging markets[END_REF]. These studies document both linear and nonlinear significant effects of oil prices on equity prices, but the sign and strength of the effects depend on the state of the market, the nature of the economic sectors, and the oil-dependence profile of the considered country.

A number of studies have recently investigated the pass-through of crude oil prices to the prices of refined petroleum products and other commodities. Most of them make use of a cointegration approach and focus on the oil-gasoline price relationships. It has been documented that downstream prices tend to respond to increases in upstream prices more rapidly than their responses to decreases in upstream prices. This asymmetric price interaction, initially discovered by [START_REF] Bacon | Rockets and feathers: the asymmetric speed of adjustment of UK retail gasoline prices to cost changes[END_REF] as the "rockets and feathers" effects, then more or less confirmed by various studies (e.g., [START_REF] Borenstein | Do gasoline prices respond asymmetrically to crude oil price changes?[END_REF][START_REF] Bachmeier | New evidence on asymmetric gasoline price responses[END_REF][START_REF] Chen | A threshold cointegration analysis of asymmetric price transmission from crude oil to gasoline prices[END_REF][START_REF] Grasso | Asymmetric error correction models for the oil-gasoline price relationship[END_REF][START_REF] Blair | The effects of hurricane Katrina on price pass-through for Gulf Coast gasoline[END_REF][START_REF] Meyler | The pass through of oil prices into consumer liquid fuels prices in an environment of high and volatile oil prices[END_REF][START_REF] Honarvar | Asymmetry in retail gasoline and crude oil price movements in the United States: an application of hidden cointegration technique[END_REF].

For example, [START_REF] Borenstein | Do gasoline prices respond asymmetrically to crude oil price changes?[END_REF] report evidence that gasoline prices adjust faster when oil prices increase than when they decrease.

Several studies have investigated the price transmission mechanisms between oil and natural gas, but most of them focus on the investigation of the linear relations between these prices (e.g., [START_REF] Yucel | Fuel taxes and cointegration of energy prices[END_REF][START_REF] Asche | The UK market for natural gas, oil and electricity: are the prices decoupled?[END_REF][START_REF] Villar | The relationship between crude oil and natural gas prices. Energy Information Administration[END_REF][START_REF] Brown | What drives natural gas prices?[END_REF][START_REF] Hartley | The relationship of natural gas to oil prices[END_REF]. In general, past findings indicate evidence of a stable long-run relationship between the two prices as well as a directional influence from oil prices to natural gas prices. For example, [START_REF] Villar | The relationship between crude oil and natural gas prices. Energy Information Administration[END_REF] show, from their bivariate VEC models, that a permanent shock of 20% in the WTI price leads to a 16% increase in the Henry Hub price. Brown and Yucel (2007) show, from an error-correction model that the movements in crude oil prices have a prominent role in shaping natural gas prices when the weather conditions, inventories and other additional factors are taken into account. Two recent studies, [START_REF] Grégoire | Using copulas to model price dependence in energy markets[END_REF][START_REF] Aloui | Extreme risk management in oil and natural gas markets[END_REF], attempt to model the nonlinear links between crude oil and natural gas prices. While the former finds evidence of extreme comovement of these prices using Student-t copulas, the latter relies on an extreme value theory -copula approach and documents the presence of extreme dependence only during the bullish market phases. The joint upward movements are typically explained by the simultaneous effects of increasing industrial and economic activities.

This article examines how oil prices are asymmetrically passed on to gasoline and natural gas prices over both the short and long run. We extend the extant literature by making use of a nonlinear autoregressive distributed lags (NARDL) model which consists in a dynamic errorcorrection representation allowing for the asymmetric responses of gasoline and natural gas prices to crude oil price changes. [START_REF] Shin | Modelling asymmetric cointegration and dynamic multipliers in an ARDL framework[END_REF] show that reliable long-run inferences can be achieved by bounds tests regardless of the integration orders of the variables in the system and that imposing either long-run symmetry or short-run symmetry or both simplifies the model to the traditional error-correction approach.

The investigation of nonlinear and asymmetric price transmissions between these energy commodities is motivated by both theoretical and empirical insights. Indeed, most economic and financial variables exhibit a nonlinear behavior over time, and that they may interact with each other in a nonlinear manner. This nonlinear behavior of time series can be driven by, as well as the complexity of financial markets induced by the heterogeneity of economic agents, the globalization, and the multiplicity of regulations. All these factors may cause unexpected changes in the behavior of economic and financial time series, which particularly induce structural breaks, asymmetric responses to news, and leverage effects. Under these conditions, the prices of crude oil, gasoline, and natural gas are likely to behave nonlinearly, and their joint dynamics seem to be more complex than a simple and stable relationships. [START_REF] Aloui | Extreme risk management in oil and natural gas markets[END_REF] note several periods in which natural gas prices deviate from their historical relationships with oil prices owing to particular market conditions or the discovery of hydraulic fracturing. [START_REF] Serletis | Testing for common features in North American energy markets[END_REF] examine the strength of shared trends and cycles between crude oil prices and Henry Hub natural gas prices, and provide evidence to reject the hypothesis of common and codependent cycles. As to the potential of nonlinear links between crude oil and gasoline prices, it is also affected by the structure of retail markets, which depends on extraction, refinery and distribution. On the other hand, the asymmetric responses of gasoline prices to crude oil prices, which have become popular with the validation of the "rockets and feathers" hypothesis by the majority of studies in the energy economics literature, are often explained by factors such as oligopolistic behavior, inventory levels, production lags, and market competition structure (e.g., [START_REF] Bacon | Rockets and feathers: the asymmetric speed of adjustment of UK retail gasoline prices to cost changes[END_REF][START_REF] Borenstein | Do gasoline prices respond asymmetrically to crude oil price changes?[END_REF][START_REF] Al Gudhea | Do retail gasoline prices rise more readily than they fall? A threshold cointegration approach[END_REF][START_REF] Grasso | Asymmetric error correction models for the oil-gasoline price relationship[END_REF]. The asymmetric transmission between crude oil and natural gas prices can arise from the fact that the latter are not only substitutes, but also rivals in consumption and production of other energy sources [START_REF] Aloui | Extreme risk management in oil and natural gas markets[END_REF]. Moreover, the crude oil market responds more to global economic, financial and geopolitical factors given its global nature, whereas natural gas prices are more sensitive to regional factors and affected by its own supply and demand disruptive shocks (e.g., extreme weather events, seasonality, and supply storage constraints).

Our empirical analysis, based on the spot prices of WTI crude oil, gasoline, and Henry Hub natural gas, provides evidence of short-run asymmetry in the oil-gasoline price relationships and long-run asymmetry in the oil-natural gas price relationships. Moreover, our results from the best-suited models for each pair of commodity prices indicate that the impact of oil price decreases on gasoline and natural gas prices is significantly larger than that of oil price increases. Our empirical evidence of asymmetric price transmission between crude oil and gasoline prices is obviously similar to that of, among others, [START_REF] Borenstein | Do gasoline prices respond asymmetrically to crude oil price changes?[END_REF] and [START_REF] Chen | A threshold cointegration analysis of asymmetric price transmission from crude oil to gasoline prices[END_REF]. However, we are able to show that the asymmetry is short-run and that it takes about 8 months for the two prices to establish a new equilibrium state, given the dynamic adjustments of gasoline prices to negative and positive changes in oil prices.

The rest of the article is organized as follows. Section 2 introduces the empirical framework and describes the data. Section 3 reports the obtained results. Section 4 concludes the article.

Methodology and data

The NARDL model

It is now common that dynamic relationships between I(1) variables can be reproduced via an error-correction model -ECM [START_REF] Granger | Cointegrating variables and error correcting models[END_REF][START_REF] Engle | Cointegration and error correction: representation, estimation and testing[END_REF][START_REF] Johansen | Statistical analysis of cointegration vectors[END_REF]. In its traditional form, the linear ECM specification without asymmetry in short-and long-run dynamics takes the following form

1 1 1 1 1 1 0 p q t y t wti t i t i i t t i i y y wti y wti µ ρ ρ α β ε - - - - - - = = ∆ = + + + ∆ + ∆ + ∑ ∑ (1)
In our study, t wti refers to the WTI oil prices, while t y represents either the gasoline prices or the natural gas prices. The symbol ∆ denotes the price variations. Even though the model in Eq. ( 1) enables the investigation of the short-and long-run relationships between price variables, it becomes impertinent when these relationships are nonlinear and/or asymmetric, which is the case for crude oil, gasoline and natural gas as we pointed out in the introduction. That is why a nonlinear and asymmetric ECM is of great interest. We thus adopt the cointegrating NARDL model of [START_REF] Shin | Modelling asymmetric cointegration and dynamic multipliers in an ARDL framework[END_REF]

( ) ∑ ∑ = = + + ∆ = ∆ = t j j t j j t wti wti wti 1 1 0 , max and
( ) ∑ ∑ = = - - ∆ = ∆ = t j j t j j t wti wti wti 1 1 0 , min .
When the asymmetries in the short-and long-run dynamics are introduced into the standard ECM, [START_REF] Shin | Modelling asymmetric cointegration and dynamic multipliers in an ARDL framework[END_REF] show that Eq. ( 1) is extended to obtain a more general cointegration model as follows ( )

1 1 1 1 1 1 0 p q t y t wti t wti t i t i i t i i t i t i i y y wti wti y wti wti µ ρ ρ ρ α β β ε - - + + - - + + - - - - - - - - = = ∆ = + + + + ∆ + ∆ + ∆ + ∑ ∑ (2)
The superscripts (+) and (-) in Eq. ( 2) denote the positive and negative partial sums decomposition defined above. p and q represent the lag order for the dependent variable and the exogenous variable in distributed lag part, respectively. The long-run symmetry can be tested by using a Wald test of the null hypothesis The short-run symmetry can equally be tested by using a standard Wald test of the null hypothesis

- + = i i β β for all 0,..., 1 i q =
-. The model in Eq. ( 2) reduces to the traditional ECM if both null hypotheses of short-run and long-run symmetry are not rejected. The non-rejection of either the long-run symmetry or the short-run symmetry will yield the cointegrating NARDL model with short-run asymmetry in Eq. ( 3) and with long-run asymmetry in Eq. ( 4), respectively.

( )

∑ ∑ - = - - - + - + - = - - - + ∆ + ∆ + ∆ + + + = ∆ 1 0 1 1 1 1 q i t i t i i t i p i i t i t wti t y t wti wti y wti y y ε β β α ρ ρ µ (3) ∑ ∑ - = - - = - - - - + - + - + ∆ + ∆ + + + + = ∆ 1 0 1 1 1 1 1 q i t i t i p i i t i t wti t wti t y t wti y wti wti y y ε β α ρ ρ ρ µ (4)
When asymmetry is detected in the ARDL model (either in the short run, in the long run or in both), the asymmetric responses to positive and negative shocks (i.e., positive or negative variations of the WTI oil prices) are respectively captured by the positive and negative dynamic multipliers associated with unit changes in one can observe, following a shock affecting the system, dynamic adjustments from the initial equilibrium to the new equilibrium between the system variables.

Data

Our data consist of time series of monthly spot closing prices for WTI crude oil, gasoline, and Table 1 also reports the results of the Zivot-Andrews unit root tests. The use of this test is justified by the potential of structural change in the price series over the study period which is characterized by several episodes of crises and important terrorist and geopolitical events. The obtained results indicate that the three price series are not stationary at the conventional levels, and are thus appropriate for the application of the NARDL model.

Results

We first estimate the specification models developed in Eqs. ( 1)-( 4) for the oil price passthrough to gasoline prices and natural gas prices. The obtained results are reported in Tables 2-3 respectively. Wald tests are then conducted to examine the hypotheses of short-run and long-run symmetry. The results from these tests applied to the estimates of the NARDL model with short-run and long-run asymmetries indicate that the long-run symmetry cannot be rejected for the oil-gasoline price relationships, but does not occur in the case of oil-natural gas price relationships. The short-run symmetry is rejected by the Wald tests for oil-gasoline case at the 5% level, but does occur for oil-natural gas case. Similar conclusions are drawn from the models incorporating either long-run or short-run symmetry. The symmetric shortrun relationship between oil prices and gasoline prices can be partially explained by the fact that gasoline is derived from crude oil and this feedstock accounts for more than 60% of the cost of gasoline. This asymmetry may also come from the public's discretionary spending adjustments in order to keep the gasoline prices stable. Nevertheless, the symmetric response is expected over the long-run because the supply of gasoline can be adjusted more easily. Differently, the symmetric short-run link between crude oil and natural gas prices can be explained by factors such as capacity conversion, natural gas-to-oil substitution, and different costs of transportation justify. 1Overall, these findings, consistent with the empirical evidence of the majority of past studies, clearly suggest that linear models are misspecified and thus do not provide the accurate fit to the underlying data. [START_REF] Frey | Econometric models of asymmetric price transmission[END_REF] review the existing empirical literature on price asymmetries in commodities from 38 studies in terms of econometric models, type of asym-metries and empirical findings. They find that half of past studies report evidence of asymmetry, five of them document the symmetry and the remaining studies report mixed results. In what follows, we detail the price transmission mechanism for each case. Note finally that all the models we use are correctly specified as the residual series do not exhibit ARCH effects at the 5% level, except for the NARDL model with short-run asymmetry for crude oil and natural gas. The non-normality of residual series is only found for the oil-natural gas case.

Pass-through of oil prices to gasoline prices

Table 2 shows that the estimated coefficients of all the four symmetric and asymmetric passthrough models are highly significant at the 1% and 5% levels, thus indicating that the price of crude oil is an important driver of gasoline prices. According to the AIC and SIC information criteria and Wald symmetry tests, the NARDL(3,2) specification with short-run asymmetry is the most suitable model for the oil-gasoline case. 2 The results of the best-suited model show that the long-run coefficient is positive (0.976) and significant at the 1% level, indicating that oil price increases (decreases) cause the price of gasoline to move up (down).

Comparatively, the restricted symmetric model also yields a significant long-run coefficient of similar magnitude (0.979), but definitively masks the short-run asymmetric effects of oil prices on gasoline prices. We indeed find that contemporaneous increases and decreases in the price of crude oil have a positive effect on gasoline prices, while the one-period lagged changes in oil prices are negatively passed through to gasoline prices. This finding corroborates the results of previous studies because not only gasoline prices respond asymmetrically to oil price movements, but there is also a temporal delay in their reaction to changes in crude oil prices. Over the short-run, the relationship between the two prices can be positive or negative depending on other factors that affect the gasoline market. It is, for example, possible that the gasoline market is facing an increase or decrease in gasoline inventories. [START_REF] Borenstein | Sticky prices, inventories, and market power in wholesale gasoline markets[END_REF] argue that refineries can adjust the production of gasoline or the level of gasoline inventory to change the supply and the price of gasoline when crude oil prices change, but these alternative adjustments may be very costly over the short period of time, which leads to a temporal delay in the response of gasoline prices to crude oil prices.

Taking the asymmetric effects into account, we see that in the short-run, an increase of 1% in oil prices entails a cumulative increase of 0.483% in gasoline price and that a decrease of 1% in oil prices leads to a cumulative decrease of 1.095% in gasoline price. Overall, the empirical findings reveal the price of crude oil as an important driver of gasoline prices, reflecting the fact that crude oil is the primary raw material used to derive gasoline as well as other refined petroleum products such as heating oil and diesel fuel. They also suggest that oil price decreases are passed on to gasoline prices over the short run and more importantly than the oil price increases. Our results are thus inconsistent with those of [START_REF] Borenstein | Do gasoline prices respond asymmetrically to crude oil price changes?[END_REF] who find no evidence of asymmetry in wholesale gasoline prices from an error-correction model with daily spot gasoline and crude-oil price data over the period 1985-1998. The results of [START_REF] Chen | A threshold cointegration analysis of asymmetric price transmission from crude oil to gasoline prices[END_REF] are also not supported by ours as they indicate that the observed asymmetry in retail gasoline price transmission primarily occurs downstream of the transmission process but not upstream.

Pass-through of oil prices to natural gas prices

Although the price ratio between crude oil and natural gas has recently changed as compared to its previous level, they still tend to comove through time. Their substitutes in consumption and also complements as well as rivals in production of other energy sources including electricity and heating oil may however induce complicated price links than common expectations of linear causality effects. The results of the Wald tests for symmetries, reported in Table 3, show that the NARDL(2,3) specification with long-run asymmetry is selected as the bestsuited model for the oil-natural gas price case. In the specification allowing for both long-run effects of oil prices are insignificant. The price of natural gas is only affected by the contemporaneous change and the one-month lagged decline in the price of oil. The results of the best-suited model indicate that natural gas prices react more to oil price decreases than to oil price increases in the long-run. It can clearly be seen that an increase of 1% in the price of crude oil causes natural gas prices to rise by 1.247%, while a decrease of 1% in the price of crude oil is associated with a fall of 1.666% in natural gas prices. These findings are also confirmed by those of the model incorporating both short-and long-run asymmetries. Notes: the SIC information criterion selects p = 3 and q = 2 as the optimal lag length. We employ a general-to-specific approach to select the final specifications by fixing max p = max q = 5. Notes: the SIC information criterion selects p = 2 and q = 3 as the optimal lag length. We employ a general-to-specific approach to select the final specifications by fixing max p = max q = 5.

Concluding remarks

Appropriate energy policies usually require the strategic development of alternative energy sources and/or the adjustments of gasoline and natural gas prices so that industrial activities are the least affected possible by expected and/or unexpected price changes in the oil market.

This issue will have ramifications for the investment decisions of private firms as gasoline derived from crude oil is used as an input to the production process and natural gas can serve as a substitute for crude oil in consumption or production of other energy sources.

In this article, a nonlinear autoregressive distributed lags model has been implemented to empirically investigate the dynamic oil price pass-through. Our empirical framework helps discern the effects of strength of crude oil price effects with respect to their positive and negative changes. Using spot prices of WTI crude oil, gasoline, and Henry Hub natural gas, we find that both gasoline and natural gas prices significantly adjust to changes in the price of oil.

Moreover, asymmetric price reaction is observed in the short-run for the gasoline and in the long-run for the natural gas. We interpret this as evidence that the gasoline price absorbs oil price shocks more quickly that the natural gas price does, because crude oil constitutes the main feedstock for the production of gasoline. The differences in nature between crude oil and natural gas markets (i.e., global versus regional markets) also explain this result. As natural gas prices vary more with regional factors (e.g., the demand characteristics of the market, the region's access to different supply basins, pipelines and storage facilities), the changes in oil prices owing principally to global factors (e.g., economic recession, financial crises, geopolitical tension, and wars) are likely to affect natural gas prices over the long-run.

Another important finding is that negative oil shocks have greater effects than positive oil shocks. The greater asymmetric impact of the decreases in oil prices can be attributed to downward price expectation spirals that affect gasoline and natural gas prices during downward economic conditions. This finding implies that the share of consumer and producer budget devoted to purchase gasoline and natural gas is greater in times of declining oil prices than in time of rising oil prices. It is worth noting however that the impact of oil price increases on gasoline prices should be more restrained because gasoline prices come under scrutiny by the public authorities on daily basis, which react to higher gasoline prices by reducing its discretionary budget. To the extent that oil prices experience substantial variations through time owing to both expected (e.g., market states and production quotas) and unexpected events (e.g., natural disasters, conflict, political tension and war), the results of this research are crucial for policymakers to draw the national energy outlook and for companies to hedge the energy risk. The growing financialization of the commodity markets under consideration, which manifests itself in increased volume of derivatives trading and activity of financial and institutional investors (e.g., commodity index funds, hedge funds and large investment banks), also strongly urge financial actors to give heed to the short-and/or long-run asymmetric price relationships. Typically, speculators may have interest to adapt their trading strategies to the complex oil-gasoline price links over the short-run as the reaction of gasoline prices to oil shocks can be negative or positive, depending on the specific factors affecting the gasoline market as well as on the government discretionary spending adjustments. The asymmetry of energy price adjustments also conveys valuable information for commodity investments because financial herding during oil price decreases and backwardation during oil price increases are likely to become stronger in financialized commodity markets.

  among others, successive episodes of economic and financial crises (e.g., the 1997-1998 Asian financial crisis, the 2001 dot-com bubble burst, the 2007 subprime crisis in the United States, and the 2008-2009 global financial crisis), wars and extreme events (e.g., the 11 September 2001 terrorist attack, the Second Gulf war from March 2003 to May 2003, and the recent Arab Spring since 2010), geopolitical tensions, and sudden changes in business cycle

  the positive and negative long-run coefficients. The short-run adjustment to a positive and a negative shock in the price of oil is captured by the parameters + i β and - i β , respectively.

  Based on the estimated multipliers,

  and negative long-run coefficients respectively. AIC and SIC are information criteria. JB and ARCH are the empirical statistics of the Jarque-Bera test for normality and the Engle (1982) test for conditional heteroscedasticity, respectively. LR W denotes the Wald test for long-run symmetry, while SR W refers to the Wald test for short-run symmetry examining the hypothesis . Standard errors of the estimated coefficients are in parenthesis. The p-values of statistical tests are in brackets. Bold numbers indicate the best model.

  and negative long-run coefficients respectively. AIC and SIC are information criteria. LR W denotes the Wald test for long-run symmetry, while SR W . JB and ARCH are the empirical statistics of the Jarque-Bera test for normality and the Engle (1982) test for conditional heteroscedasticity, respectively. Standard errors are in parenthesis. The p-values of statistical tests are in brackets. Bold numbers indicate the best model.

Figure 2

 2 Figure 2 Cumulative asymmetric adjustments of natural gas prices Model with LR and SR asymmetryModel with LR symmetry and SR asymmetry

Table 1 . Stochastic properties of commodity price series

 1 Henry Hub natural gas. They are obtained from Bloomberg and cover the period from January 1997 to September 2012. Commodity prices are expressed in US dollars. Table1presents the descriptive statistics for price series as well as their stochastic properties. The monthly average prices range from $1.4360 per gallon for gasoline to $51.536 per barrel for crude oil. On a JB is the empirical statistics of the Jarque-Bera test for normality. ZA is the empirical statistics of the Zivot and Andrews (1992) unit root test which is robust to structural breaks. + , ++ , and +++ indicate rejection of the null hypotheses of normality and unit root.

		Crude oil	Gasoline price	Natural gas
	Mean	51.536	1.436	4.794
	Maximum	133.880	3.292	13.420
	Minimum	11.350	0.307	1.720
	Standard deviation	30.152	0.831	2.410
	Skewness	0.580	0.572	1.090
	Kurtosis	2.274	2.117	4.248
	JB	14.767 +	16.455 +	49.706 +
	ZA	-4.667	-3.615	-4.437
	Notes:

monthly basis, the WTI oil and gasoline prices reach their maximum value of $133.88 and $1.436 respectively in July 2008, while the highest peak in natural gas price ($13.420) is observed in October 2005. All the series are positively skewed and display significant excess kurtosis, which means that they have fatter tails and longer left tails than a normal distribution. The Jarque-Bera test indicates the non-normality of all the price series.

Table 2 . Estimation results for oil price pass-through to gasoline prices

 2

		Symmetric ARDL			NARDL with LR asymmetry		NARDL with SR asymmetry	NARDL with LR and SH asymmetry
	gasoline	t	1 -	-0.314 *** (0.061)	gasoline	t	1 -	-0.322 *** (0.062)	gasoline	t	1 -	-0.347 *** (0.060)	gasoline	t	1 -	-0.359 *** (0.061).)
	1 -gasoline t wti ∆ ∆ gasoline t wti ∆ 1 -∆ t wti Constant	1 -2 -t t	0.308 *** (0.059) 0.183 ** (0.074) -0.104 ** (0.041) 1.060 *** (0.049) -0.183 ** (0.091) -0.106 *** (0.021)	+ -1 t wti --1 t wti ∆ gasoline ∆ gasoline t wti ∆ 1 -∆ t wti	1 -2 -t t	0.303 *** (0.060) 0.297 *** (0.060) 0.184 ** (0.074) -0.097 ** (0.041) 1.049 *** (0.050) -0.176 * (0.091)	1 -gasoline t wti ∆ ∆ gasoline + ∆ t wti -∆ t wti + -∆ 1 t wti	t	1 -2 -t	0.338 *** (0.058) 0.198 *** (0.072) -0.086 ** (0.040) 0.739 *** (0.096) 1.325 *** (0.086) -0.256 ** (0.115)	+ -1 t wti --1 t wti ∆ gasoline ∆ gasoline + ∆ t wti -∆ t wti	1 -2 -t t	0.334 *** (0.058) 0.327 *** (0.059) 0.198 *** (0.072) -0.078 * (0.041) 0.721 *** (0.097) 1.313 *** (0.086)
							Constant				-0.119 *** (0.024)	--1 t wti ∆		-0.230 * (0.123)	+ -1 t wti ∆	-0.262 ** (0.115)
														Constant		-0.095 *** (0.020)	--1 t wti ∆	-0.207 * (0.214)
																	Constant	-0.112 ***
																	(0.024)
	L	wti				0.979 *** [0.000]	L	wti	+				0.939 *** [0.000]	L	wti		0.976 *** [0.000]	L	wti	+	0.932 *** [0.000]
							L	wti	-				0.921 *** [0.000]				L	wti	-	0.912 *** [0.000]
	AIC				-534.658	AIC				-533.719	AIC
	SIC				-512.078	SIC					-507.912	SIC
	JB				3.009 [0.222]	JB					4.350 [0.114]	JB
	ARCH				6.357 [0.273]	ARCH				6.717 [0.243]	ARCH
							LR W				1.031 [0.311]

SR W -

545.591 -516.559 2.263 [0.322] 7.291 [0.120] 8.906 [0.003]

	AIC	-545.305
	SIC	-513.047
	JB	2.892 [0.235]
	ARCH	7.714 [0.173]
	W	LR	1.662 [0.199]
	W	SR	9.820 [0.002]

Table 3 . Estimation results for oil price pass-through to natural gas prices

 3

					Symmetric ARDL			NARDL with LR asymmetry			NARDL with SR asymmetry	NARDL with LR and SH asymmetry
	NatGas	t	1 -	-0.053 ** (0.024)	NatGas	t	1 -	-0.086 *** (0.028)	NatGas	t	1 -	-0.054 ** (0.025)	NatGas	t	1 -	-0.091 *** (0.028)
	1 -NatGas t wti ∆ t t wti ∆ 1 -∆ t wti 2 -∆ t wti Constant	1 -	0.011 (0.019) 0.051 (0.074) 0.318 *** (0.115) -0.074 (0.120) 0.215 * (0.115) 0.070 ** (0.032)	+ -1 t wti --1 t wti ∆ NatGas t t wti ∆ 1 -∆ t wti 2 -∆ t wti	1 -	0.107 ** (0.044) 0.142 ** (0.058) 0.048 (0.073) 0.394 *** (0.118) -0.122 -0.120) 0.145 (0.117)	1 -NatGas t wti ∆ t + ∆ t wti -∆ t wti + -∆ 1 t wti --∆ 1 t wti	1 -	0.011 (0.019) 0.046 (0.074) 0.501 ** (0.097) 0.198 (0. 207) 0.297 (0.234) -0.361 * (0.215)	+ -1 t wti --1 t wti ∆ NatGas t + ∆ t wti -∆ t wti + -∆ 1 t wti	1 -	0.120 *** (0.045) 0.161 *** (0.060) 0.037 (0.073) 0.602 ** (0.231) 0.286 (0.206) 0.309 (0.230)
						Constant			0.172 *** (0.053)	+ -t wti ∆	2	0.136 (0.231)	--1 t wti ∆	-0.455 ** (0.216)
													--t wti ∆	2	0.350 * (0.204)	+ -t wti ∆	2	0.171 (0.227)
													Constant	0.045 *** (0.037)	--t wti ∆	2	0.189 (0.210)
																	Constant	0.150 ***
																	(0.054)
	wti L				0.208 [0.510]	L	wti	+				1.247 *** [0.002]	L	wti		0.211 [0.497]	L	wti	+	1.320 *** [0.001]
						L	wti	-				1.666 *** [0.003]					L	wti	-	1.764 *** [0.001]
	AIC				-219.079	AIC				-222.854	AIC		-217.609	AIC	-222.756
	SIC				-196.499	SIC					-197.048	SIC		-185.352	SIC	-187.268
	JB				27.554 [0.000]	JB					38.187 [0.000]	JB		1.788 [0.409]	JB	46.383 [0.000]
	ARCH				2.386 [0.793]	ARCH				2.315 [0.804]	ARCH		12.148 [0.033]	ARCH	3.491 [0.625]
						LR W				6.580 [0.011]	W	SR		1.723 [0.191]	W	LR	8.122 [0.005]
																	W	SR	3.427 [0.066]

We are particularly grateful to an anonymous reviewer for these valuable suggestions.

The best model is the one that has the lowest AIC and SIC values.

Asymmetric adjustment paths

The asymmetric adjustments from an initial long-run equilibrium to a new long-run equilibrium after a negative or positive unitary shock affecting the crude oil market can be learned from the dynamic multipliers. Figures 12show the predicted dynamic multipliers for the adjustment of gasoline and natural gas prices under the four NARDL specifications we consider.

The asymmetry curve depicts the linear combination of the dynamic multipliers associated with positive and negative shocks. The positive and negative change curves provide the information about the asymmetric adjustment to positive and negative shocks at a given forecasting horizon respectively. Lower band and upper band for asymmetry indicate the 95% confidence interval.