I W Martin
email: i.martin@physics.gla.ac.uk

R Bassiri

R Nawrodt

M M Fejer

A Gretarsson

E Gustafson

G Harry

J Hough

I Maclaren

S Penn

S Reid

R Route

S Rowan

C Schwarz

P Seidel

J Scott

A L Woodcraft

Effect of heat treatment on mechanical dissipation in Ta 2 O 5 coatings

Keywords: numbers: 05.40.Jc, 61.43.Er, 68.60.Bs, 04.80.Nn, 95.55.Ym, 42.60.Da

Thermal noise arising from mechanical dissipation in dielectric reflective coatings is expected to critically limit the sensitivity of precision measurement systems such as high-resolution optical spectroscopy, optical frequency standards and future generations of interferometric gravitational wave detectors. We present measurements of the effect of post-deposition heat treatment on the temperature dependence of the mechanical dissipation in ion-beam sputtered tantalum pentoxide between 11 K and 300 K. We find the temperature dependence of the dissipation is strongly dependent on the temperature at which the heat treatment was carried out, and we have identified three dissipation peaks occurring at different heat treatment temperatures. At temperatures below 200 K, the magnitude of the loss was found to increase with higher heat treatment temperatures, indicating that heat treatment is a significant factor in determining the level of coating thermal noise.

Introduction

Several long baseline gravitational wave detectors have been constructed and used in searches for gravitational radiation from astrophysical sources [1,2,3,4]. These detectors use interferometric sensing to search for displacements, induced by gravitational waves, of test masses which are suspended as pendulums at the ends of perpendicular arms up to 4 km in length. The test masses are coated to form mirrors which are highly reflective at 1064 nm. The thermal displacement noise associated with the test mass mirrors is a significant limit to the sensitivity of these detectors [5,6,7]. The magnitude of the thermal noise is directly related to the mechanical dissipation of the materials used to construct the mirrors. In particular, the mechanical dissipation of the amorphous dielectric mirror coatings has been identified as an important noise source which is likely to limit the sensitivity of the generation of detectors currently being constructed, such as Advanced LIGO [START_REF] Harry | for the[END_REF] and Advanced Virgo [START_REF] Flaminio | Advanced Virgo White Paper[END_REF], in their most sensitive frequency band [START_REF] Levin | [END_REF]11,12,13]. Coating thermal noise is also expected to significantly limit the performance of highly frequency-stabilized lasers for use in high resolution optical spectroscopy [14], optical frequency standards [15,16] and fundamental quantum measurements [17].

The mirror coatings used in current gravitational wave detectors are formed from alternating layers of silica (SiO 2 ) and tantalum pentoxide (Ta 2 O 5 ) deposited by ionbeam sputtering. Research has shown that the mechanical dissipation of these multilayer coatings is dominated by the dissipation of the Ta 2 O 5 layers [18,19,20] and that the total dissipation of a multi-layer coating can be reduced by up to ∼ 40 % by doping the Ta 2 O 5 with titania (TiO 2 ) [21]. However, neither the process responsible for mechanical dissipation in Ta 2 O 5 nor the mechanism by which TiO 2 doping reduces the dissipation is as yet understood.

In a recent paper we reported data showing a low temperature dissipation peak in a Ta 2 O 5 coating doped with TiO 2 and heat treated at 600 • C [22]. There is evidence that this peak arises from a thermally activated dissipation mechanism, possibly related to the reorientation of Ta-O bonds within a double well potential, similar to the mechanism believed to occur in fused silica [START_REF] Braginsky | Systems with small dissipation[END_REF][START_REF] Bömmel | [END_REF]25]. Further work has suggested that the effect of the TiO 2 doping is to increase the activation energy of the dissipation process [26], thereby reducing the mechanical loss of the material. A barrier height distribution function analysis revealed that the effect of TiO 2 doping was to shift the distribution of potential barriers to a higher energy and to broaden the peak in the distribution. Heat treatment is another possible method of altering the distribution of potential barrier heights. Ion-beam sputtered coatings are often heat treated after deposition to reduce the stress in the film and to reduce the optical absorption [START_REF] Netterfield | Proc. SPIE[END_REF]. Several studies have shown that post-deposition heat treatment of ion-beam sputtered amorphous oxide layers can increase the thickness of the film, with a corresponding reduction in density [START_REF] Netterfield | Proc. SPIE[END_REF][START_REF] Tilsch | Proc. SPIE[END_REF][START_REF] Brown | [END_REF]30]. Furthermore it is known heat treatment can significantly reduce the mechanical loss of bulk fused silica [31]. This is believed to be a result of a reduction of internal stresses in the silica, possibly altering the distribution of the potential barrier heights [START_REF] Lunin | Physical and Chemical Bases for the Development of Hemispherical Resonators for Solid-State Gyroscopes[END_REF]. Thus investigating the effects of heat treatment on the mechanical loss of Ta 2 O 5 is of particular interest for understanding, and perhaps reducing, the mechanical dissipation. In addition to providing insights into the mechanisms of energy dissipation in Ta 2 O 5 , measurements of the temperature dependence of the coating dissipation are critical for the design of future gravitational wave detectors such as ET [START_REF] Punturo | Einstein Gravitational Wave Telescope, proposal to the European Commission, Framework Programme 7[END_REF] and LCGT [START_REF] Kuroda | Status of LCGT[END_REF] which may use cryogenic cooling to reduce thermal noise.

Sample preparation and experimental technique

The Ta 2 O 5 coatings were deposited onto rectangular silicon cantilever substrates. These substrates consist of a 0.5-mm-thick clamping block and a thinner flexing portion of length 34-mm and thickness 50 µm, as shown in Figure 1. The sample geometry was designed to reduce frictional energy losses into the clamp used to support the cantilever during mechanical loss measurements [35,36]. A thermal oxide layer of approximately 30 nm thickness was grown on the cantilevers prior to coating to ensure proper adhesion of the Ta 2 O 5 . A Ta 2 O 5 film (0.49 ± 0.1) µm in thickness was deposited on one face of the flexing part of each cantilever by ion-beam sputtering using argon as the sputtering ion. Four fused silica disks, 25.4 mm in diameter, were coated in the same coating run to serve as witness samples for microscopic analysis of the coating structure. A coated cantilever, a nominally identical uncoated cantilever and a coated silica disk were heat treated in air at each of the following temperatures: 300, 400, 600 and 800 • C. The samples were heated to the desired temperature at 2 • C per minute and held at this temperature for 24 hours after which they were allowed to cool naturally. The coating and heat treatment process was carried out by CSIRO †. Ion-beam sputtered Ta have been observed to crystallize when heated to approximately 600 to 700 • C [37], to the detriment of their optical properties. The highest heat treatment temperature of 800 • C was chosen to allow study of the mechanical loss on either side of this transition region.

The mechanical dissipation factors of several resonant modes of each coated cantilever and the associated uncoated control cantilever were measured between 11 and 300 K in a cryostat using a ring-down technique, in which a resonant mode of the sample of angular frequency ω 0 was excited to an amplitude A 0 and the vibration allowed to decay freely. The time dependence of the amplitude ring-down A(t) is given by:

A(t) = A 0 e -φ(ω 0 )ω 0 t/2 , ( 1
)
where φ(ω 0 ) is the mechanical dissipation factor. The cantilever was mounted in a stainless steel clamp within the vacuum chamber of the cryostat and the bending modes excited by an electrostatic actuator located a few mm below the cantilever. The frequency of the n th bending mode can be calculated as follows [START_REF] Mclachlan | Theory of Vibrations[END_REF]:

ω n = α 2 n t 4π √ 3L 2 Y ρ 1/2 (2)
where ρ is the density, Y the Young's modulus, t the thickness and L the length of the cantilever. α n is a numerical factor which takes the values 1.875, 4.694, 7.855 and 10.996 for n = 1 to 4 respectively. For n ≥ 5, α n can be approximated as α n (2n -1)π/2. The vibration amplitude was measured by means of a laser beam reflected from the cantilever surface and directed onto a photodiode sensor outside the cryostat. The sensor consisted of a rectangular photodiode masked so that a thin triangular strip of the surface was exposed. Vibrations of the cantilever resulted in motion of the reflected laser beam along the triangular sensor, producing a signal proportional to the vibration amplitude. This sensor arrangement was found to be more practical than a quadrant photodiode due to the large amplitude of motion of the laser spot. Several measurement cycles, in each of which the sample temperature was increased incrementally from 11 to 300 K, were carried out. Repeated ring-down measurements at each temperature generally showed a variation in loss factor of less than 5 % for each mode. The variation in the loss factor between the various temperature cycles was typically less than 10%. The temperature of the cantilever was measured using a silicon-diode sensor (Lakeshore DT-670 A) mounted within the clamp immediately below the fixed end of the cantilever.

Results

Figure 2 shows, as an example, the loss factors measured for the coated and uncoated cantilevers heat treated at 300 • C. The increase in loss due to the presence of the Ta 2 O 5 film can be clearly observed throughout the temperature range studied. Of particular note is the peak in the loss of the coated cantilever at approximately 35 K -this will be discussed in more detail in the following section. The dissipation of the uncoated cantilever also exhibits a peak at between 25 to 30 K. A similar peak was observed for all of the uncoated samples studied. Peaks have been observed at similar temperatures in silica films grown on silicon by wet thermal oxidation [START_REF] White | [END_REF], and it seems likely therefore that these peaks are due to the thermal oxide layer grown on all of the cantilevers. The mechanical loss of the Ta 2 O 5 film was calculated from the difference in the loss factors of the coated and uncoated cantilevers using the following equation [40]:

φ(ω 0 ) coating = Y s t s 3Y c t c (φ(ω 0 ) coated -φ(ω 0 ) substrate ), (3)
where ω 0 is the angular frequency of the mode, φ(ω 0 ) coated is the loss factor of the coated cantilever, φ(ω 0 ) substrate is the loss factor of the cantilever prior to coating, t s and Y s are the thickness and Young's modulus of the substrate respectively and t c and Y c are respectively the thickness and Young's modulus of the coating. The Young's moduli of silicon and Ta 2 O 5 were taken to be 166 GPa [START_REF] Touloukian | Thermo-physical Properties of Matter[END_REF] and (140 ± 15) GPa [START_REF] Martin | Materials Research Society Symposium Proceedings[END_REF] respectively. Assuming that the temperature dependence of the Young's modulus of Ta 2 O 5 is typical of other amorphous oxides [START_REF] Marx | [END_REF], then its effect on the calculation over the temperature range studied here is negligible. Figure 3 shows the calculated mechanical dissipation of the Ta 2 O 5 film heat treated at 300 • C, obtained from the data in Figure 2 using Equation 3. The same procedure was followed to find the loss of the coatings which had been heat treated at 400, 600 and 800 • C. Data was obtained at several mode frequencies for each sample, and Figures 4 to 6 show comparisons of the loss of the series of coatings at 1000 Hz, 1900 Hz, 3000 Hz and 350 Hz respectively. The results clearly show that the temperature of the post-deposition heat treatment had a significant effect on the temperature dependence of the mechanical loss factor of the Ta 2 O 5 coatings. There are several notable features in the data. Firstly, a broad peak in the dissipation of the Ta 2 O 5 film heat treated at 300 • C was observed at approximately 35 K. The data for the coated cantilever which was heat treated at 400 • C is incomplete as the cantilever cracked in the clamp during a measurement cycle. The breakage was indicated by a sudden change in the resonant frequencies of the cantilever and an increase in the measured mechanical dissipation factors. However, below 100 K the loss of this coating follows a similar trend to the coating heat treated at 300 • C, with a broad dissipation peak at 35 K.

Secondly, a sharper additional peak in the dissipation of the Ta 2 O 5 film heat treated at 600 • C was observed at approximately 20 K. This is in qualitative agreement with previous measurements of Ta 2 O 5 films supplied by LMA † and heat treated at the same temperature [22,26]. This peak at 20 K is significantly narrower and higher than the peak at 35 K observed in the coatings heat treated at 300 and 400 • C. At temperatures above 40 K, the dissipation of the coating heat treated at 600 • C shows a similar trend in shape to that for the 300 and 400 • C coatings, although the absolute level of dissipation is somewhat higher below 250 K. This suggests that the 35 K peak may also be present in the film heat-treated at 600 • C, partially underlying the tail of the 20 K peak.

Thirdly, the Ta 2 O 5 coating treated at 800 • C exhibits a large and very broad dissipation peak at 80 to 90 K. A small plateau in the low temperature edge of this peak is apparent at each measurement frequency, suggesting that the 20 K dissipation peak is also present. At temperatures between 200 and 250 K, this coating has a lower dissipation than those heat treated at lower temperatures, with the minimum dissipation occurring at approximately 225 K. The increase in the dissipation between 225 K and 300 K may be indicative of an additional dissipation peak occurring above room temperature. However, further investigation was beyond the scope of the current study, as the sensitivity with which the loss factor of the coating can be measured is significantly reduced at room temperature and above due to the relatively high dissipation factor of the silicon cantilever substrates at these temperatures. As noted above, ion-beam sputtered Ta 2 O 5 has been observed to crystallize at temperatures between 600 and 700 • C. The structure of the Ta 2 O 5 films was investigated by electron diffraction measurements of the the witness samples which were coated and heat treated at the same time as the cantilevers. The results, shown in Figure 8, indicate that the coatings heat treated at 600 • C and below are amorphous, while the coating heat treated at 800 • C has a crystalline structure. The difference in structure is also apparent from comparison of transmission electron microscope (TEM) images of the † Laboratoire des Matériaux Avancés, LMA, CNRS-IN2P3, France, http://lma.in2p3.fr. coatings heat treated at 600 and 800 • C, as shown in Figure 9. It is likely therefore that the large dissipation peak centered on 80 to 90 K is associated with the poly-crystalline structure of the Ta 2 O 5 induced by the heat treatment at 800 • C. One possible dissipation mechanism in poly-crystalline materials is phonon scattering at the grain boundaries [START_REF] Nowick | Anelastic Relaxation in Crystalline Solids[END_REF].

Analysis

For both the dissipation peak near 20 K and the peak near 35 K, the temperature at which the peak occurred was found to vary with frequency. This behaviour, which can be seen clearly in Figure 10 for the coating heat treated at 600 • C, is a typical characteristic of a thermally activated dissipation process. These processes can be characterized by a rate constant, τ 0 , and an activation energy, E a , which are related by the Arrhenius equation [START_REF] Nowick | Anelastic Relaxation in Crystalline Solids[END_REF]:

τ = τ 0 e E a /k B T ( 4
)
where τ is the relaxation time associated with the dissipative system returning to equilibrium after a perturbation. Analysis of these dissipation processes shows that the temperature of the dissipation peak, T peak , at angular frequency ω 0 is related to the activation energy and rate constant as follows [START_REF] Nowick | Anelastic Relaxation in Crystalline Solids[END_REF]:

ω 0 τ 0 e E a /k B T peak = 1. (5)
Thus a plot of log ω 0 against 1/T peak should yield a straight line, from which the activation energy and rate constant for the dissipation process can be obtained. Figure 11 shows this analysis for the peak near 20 K in the coating heat treated at 600 • C. The activation energy and rate constant, calculated from the linear fit, were found to be (35.6 ± 2.5) meV and (9.9 ± 0.5)×10 -13 s respectively. The activation energy is somewhat higher than a previous measurement of (28.6 ± 1.6) meV for Ta 2 O 5 deposited by LMA. The difference in activation energy could arise from differences in the detailed structure of the amorphous coating induced by variations in the ion-beam sputtering deposition process between coating vendors. The Arrhenius plot for the peak near 35 K observed for the coating heat treated at 300 • C is shown in Figure 12. A straight line could be fitted to the data for four of the modes studied. However, the point on the Arrhenius plot corresponding to the mode at 3200 Hz clearly does not follow this linear trend. It is possible that a torsional mode close in frequency to the 5 th bending mode was measured unintentionally, which may explain the difference in the loss peak for this mode. If this point is neglected the fit for the remaining modes gives an activation energy of (66 ± 10) meV and a rate constant of (9.4 ± 0.9)×10 -14 s. In many amorphous solids, the dissipation at temperatures above 10 K is thought to result from thermally activated reorientations of bonds over potential barriers of height V [START_REF] Gilroy | [END_REF]. If the dissipation is modelled as arising from transitions in an asymmetric double-well potential in some local configuration coordinate, it can be shown that the measured mechanical loss is related to the distribution of energy barrier heights g(V ) by the following expression [46]:

φ = πγ 2 f 0 C ii k B T g(V ), (6)
where

V = k B T ln( 1 ωτ 0 ). ( 7
)
and C ii is the appropriate elastic constant for the mode-shape under consideration, γ is the elastic coupling constant which represents the coupling between the defect (e.g. a bond re-orienting within the double well potential) and the applied strain and f 0 is a constant representing the asymmetry between the depth of the two potential wells in each double-well system. The relaxation time τ associated with a barrier height V is given by Equation 4. The assumptions made in this model are discussed in detail by Topp and Cahill [46]. Thus the measured mechanical loss can be used to extract the distribution of barrier heights (i.e. activation energies) for the loss mechanism associated with a particular dissipation peak. The extracted barrier height distributions for the Ta 2 O 5 heat treated at 300 • C and 600 • C are shown in Figure 13. At all barrier energies, the magnitude of the function g(V )f 0 is larger for the coating heat treated at 600 • C. An interpretation of this is that there are more barriers contributing to the dissipation in this coating, resulting in a higher mechanical dissipation. The feature in the distribution function at a barrier height energy of 35 meV corresponds to the average activation energy calculated from the Arrhenius plot. This feature appears to be superimposed on a background distribution similar to that found for the coating heat treated at 300 • C.

The shape of the barrier height distribution for the 35 K loss peak in the coating heat treated at 300 • C is broadly similar to the distribution found in bulk fused silica [46], with g(V )f 0 decreasing smoothly with increasing barrier energy. We postulate that the dissipation mechanisms in Ta 2 O 5 may be similar to that in fused silica and other amorphous solids, and be related to the reorientation of oxygen atoms, or groups of tantalum and oxygen atoms, in a double-well potential. Several models of the dissipation mechanism in fused silica have been proposed, including motion of Si-O bonds between stable angles [47], elongation of Si-O bonds [25], rotation of SiO 4 tetrahedral structures [48,49] or of groups of these structures [50]. If similar processes are responsible for the loss in Ta 2 O 5 , it is possible that the peaks observed at 20 K and 35 K correspond to two different types of motion, one of which is only activated by changes in the local ordering induced by heat treatment at temperatures close the the crystallization temperature. Further detailed studies of the short range structure in Ta 2 O 5 coatings using electron diffraction data and computer modelling are ongoing [51].

Conclusions

The temperature dependence of the mechanical dissipation in Ta 2 O 5 films has been found to be strongly dependent on the temperature at which post-deposition heat treatment is carried out. We have identified three dissipation peaks arising from dissipation mechanisms which occur after heat treatment at different temperatures. The previously observed dissipation peak at 20 K has been shown to arise from heat treatment close to the temperature at which the amorphous Ta 2 O 5 crystallizes. The dissipation mechanism responsible for this peak may therefore be related to pre-crystallization changes in the local ordering of the Ta 2 O 5 film. A dissipation peak at approximately 35 K was observed in Ta 2 O 5 films heat-treated at 300 and 400 • C. There is evidence that this peak is also present in Ta 2 O 5 heat treated at 600 • C, underlying the larger dissipation peak at 20 K. We postulate that the two low temperature dissipation peaks may arise from two different types of thermally activated transitions. The Ta 2 O 5 film heat-treated at 800 • C, which was found to have crystallized, had a very broad dissipation peak centered on 90 -100 K. Below 200 K the dissipation of the Ta 2 O 5 coatings was found to increase at higher heat treatment temperatures, suggesting that minimizing the temperature at which post-deposition heat treatment is carried out may allow a significant decrease in the coating thermal noise in future cryogenic gravitational wave detectors, thus extending their astronomical reach. As heat treatment is also known to affect the optical absorption of dielectric coating, studies of the optical properties of Ta 2 O 5 heat-treated at various temperatures would be of interest.

The dissipation of Ta 2 O 5 coatings at room temperature appears to be related to the broad tail of the dissipation peak at 35 K. However there is some evidence, particularly for the coating heat-treated at 800 • C, of an increase in dissipation close to room temperature. Thus further study of the dissipation in Ta 2 O 5 at temperatures above 300 K would be of significant interest to determine if additional dissipation peaks occurring at temperatures greater than 300 K may contribute to the room temperature dissipation.

Further work is underway to investigate the relationship between the short range structure of Ta 2 O 5 films and the mechanical loss. In addition, it would be of significant interest to study the emergence of the various dissipation peaks in Ta 2 O 5 by studying coatings heat treated at more temperature intervals between 400 and 800 • C.

Figure 1 .

 1 Figure 1. A schematic diagram of one of the silicon cantilever samples used for the coating loss measurements. The cantilever length is parallel to the 110 crystal axis.

Figure 2 .Figure 3 .

 23 Figure 2. Measured mechanical loss for a cantilever coated with Ta 2 O 5 and an uncoated control cantilever heat treated at 300 • C

Figure 4 .Figure 5 .Figure 6 .Figure 7 .

 4567 Figure 4. The calculated mechanical loss factors of the various Ta 2 O 5 coatings at approximately 1000 Hz

Figure 8 .

 8 Figure 8. (a)Electron diffraction pattern from the Ta 2 O 5 coating heat treated at 600 • C. The diffuse rings indicate an amorphous structure. (b) Electron diffraction pattern from the Ta 2 O 5 coating heat treated at 800 • C, indicating crystalline structure.

Figure 9 .

 9 Figure 9. (a) Bright field TEM image of a typical amorphous Ta 2 O 5 film (darker region), heat treated at 600 • C.(b) Dark field TEM image of a crystallized Ta 2 O 5 film (lighter region) which was heat treated at 800 • C.

Figure 10 .

 10 Figure 10. Mechanical dissipation of the coating heat treated at 600 • C at several frequencies, showing the frequency dependence of the temperature of the peak in the low temperature dissipation.

Figure 11 .Figure 12 .

 1112 Figure 11. Arrhenius plot for loss peak at 20 K observed in the Ta 2 O 5 coating heat treated at 600 • C

Figure 13 .

 13 Figure 13. Comparison of the calculated potential barrier height distribution function g(V )f 0 as a function of barrier height for Ta 2 O 5 heat treated at 300 • C and 600 • C, using the data measured at ∼ 1000 Hz.

† Industrial Physics Division, Commonwealth Scientifc and Industrial Research Organisation, West Lindfield, NSW, Australia

Acknowledgments

We are grateful for the financial support provided by STFC, the University of Glasgow, EPSRC, the German Science Foundation under contract SFB Transregio 7 and the National Science Foundation under grants PHY-07 57896 and PHY-05 02641 (Stanford) and NSF-0653590 (HWS). I. Martin and S. Reid are supported by an STFC Postdoctoral Fellowship and a Royal Society of Edinburgh Research Fellowship respectively. R. Nawrodt acknowledges the support of the FP7 EU Einstein Telescope design study project under contract number 211743. LIGO was constructed by the California Institute of Technology and Massachusetts Institute of Technology with funding from the National Science Foundation and operates under cooperative agreement PHY-0107417. We would like to thank our colleagues in the LSC and VIRGO collaborations and within SUPA for their interest in this work. This paper has LIGO Document Number LIGO-P1000074.


Christophe Chesneau
email: chesneau@math.unicaen.fr

Maher Kachour
email: maher.kachour@univ-nantes.fr

A parametric study for the first-order signed integer-valued autoregressive process

Keywords: Integer-valued time series, INAR models, SINAR models, Rademacher(p) -N class, Skellam distribution. 2000 Mathematics Subject Classification 62M10, 62M20

In recent years, many attempts have been made to find accurate models for integer-valued times series. The SINAR (for Signed INteger-valued AutoRegressive) process is one of the most interesting. Indeed, the SINAR model allows negative values both for the series and its autocorrelation function. In this paper, we focus on the simplest SINAR(1) model under some parametric assumptions. Explicitly, we obtain the form the probability mass function of the innovation when the marginal distribution of the process is known. Moreover, we give an implicit form of the stationary distribution for a known innovation. Simulation experiments as well as analysis of real data sets are carried out to attest the models performance.

Introduction

Integer-valued times series occur in many practical situations and several application fields. The observed data can be on N, e.g. a counting process, or on Z, e.g. the meteorological data. For both cases, the analysis of such a times series may present some difficulties, especially stochastic models must reflect the integer peculiarity of the observed series. For example, a standard AR(p) process does not guarantee this peculiarity. In the past decades, many attempts have been made to propose models adapted to integer-valued times series.

Prior art on integer-valued times series on N. The class of INAR (for INteger-valued AutoRegressif) models is based on the thinning operator, denoted by •, defined as follows. Let Z be non-negative integer-valued random variable and a ∈ [0, 1]. According to [START_REF] Steutel | Discrete analogues of self-decomposability and stability[END_REF],

a • Z = Z k=1 ξ k ,
where the so-called counting sequences (ξ k ) k∈Z are independent and identically distributed (i.i.d.) such that ξ 1 ∼ Bernoulli(a) i.e. P(ξ 1 = 1) = a = 1 -P(ξ 1 = 0), and are independent of Z. The sequence (ξ k ) k∈Z is called a counting sequences.

Based on this thinning operator, Al-zaid and Al-Osh (1987) defined an INAR(1) model mimics a real AR(1) process. In [START_REF] Al-Zaid | An integer-valued pth-order autoregressive structure (INAR(p)) process[END_REF], the INAR(p) model has a correlation structure similar to the one of an ARMA(p, p-1) process. [START_REF] Du | The integer-valued autoregressive (INAR(p)) model[END_REF] introduced another INAR(p) process with a correlation structure identical to the one of a real AR(p) process. [START_REF] Gauthier | Convergence forte des estimateurs des paramètres d'un processus GENAR(p)[END_REF] considered a more general version of the model, the GINAR(p) model, based on the generalized Steutel & van Harn operator, where the i.i.d. sequence (ξ k ) k∈Z (the counting sequences) follows a discrete distribution with a finite mean a ∈ [0, 1] and a finite variance µ. Hence, INAR(p) becomes a special case of GINAR(p) model. For a recent review on discrete variate time series models, we refer for example to [START_REF] Kedem | Regression models for time series analysis[END_REF], [START_REF] Mckenzie | Discrete variate time series[END_REF], Jung and Tremayne (2006) and [START_REF] Weiß | Thinning operations for modeling time series of counts a survey[END_REF].

However, all the previous models have the same serious limitations. Thus, because of the thinning operators, all the coefficients of the models must be positive. Therefore, the modeling of series with possible negative autocorrelations is excluded. Moreover, these models are defined on N, and so they cannot fit a time series with negative observations.

Prior art on integer-valued times series on Z. Integer-valued times series with negative values are very frequent in practice. But also there are other applications when values on Z is the domain of interest. For example, when analyzing a non-stationary integer-valued time series with non-negative values, i.e. observations with a time trend or seasonality. In such case we need to use the differencing operator in order to achieve stationarity. Therefore, we obtain a time series on Z. Another example, let (X t ) be an observation of a counting process where the values are very high, and we want to reduce the variations while keeping the integer peculiarity of the observed series. In some way, we are going to analyze the centered series (Y t ) instead of (X t ). Note that, the empirical mean m of (X t ) is real rather than integer-valued in all but very rare cases. A mapping into the discrete support of the series is readily obtained by applying Gaussian brackets (integer part of), denoted by [•], or by using the rounding operator, denoted by • . Hence, one can define Y t = X t -[m] or Y t = X t -m . The observations of (Y t ) have a weak variation and they oscillate around 0 and consequently we have a time series on Z.

In spite of the importance of such a modeling, it seems that few models have been proposed. A first model, we owe it to [START_REF] Kim | A non-stationary integer-valued autoregressive model[END_REF], it is based on the signed binomial thinning operator and denoted by INARS. Then, [START_REF] Zhang | Inference for INAR(p) processes with signed generalized power series thinning operator[END_REF] propose a generalized version of the previous model, where the counting sequences have a generalized power series as common distribution. Recently, [START_REF] Kachour | A p-order signed integer-valued autoregressive (SINAR(1)) model[END_REF] introduced a more general class, denoted by SINAR, based on a modified version of the generalized thinning operator, called also the signed thinning operator which is different than that introduced by [START_REF] Kim | A non-stationary integer-valued autoregressive model[END_REF]. Recall that all these processes have a autocorrelation structure identical to that of a real AR. Furthermore, consistency and the asymptotic properties of estimators proposed to estimate the parameters of these models have been well studied in the papers cited above. Note that in a different context from the thinning operation, [START_REF] Kachour | First-order rounded integer-valued autoregressive (RI-NAR(1)) process[END_REF] and [START_REF] Kachour | On the rounded integer-valued autoregressive process[END_REF] proposed the RINAR model, based on the rounding operator. They showed that the stationarity and ergodicity of their model is ensured under similar stationary condition as in a AR process. On the other hand, the rounding operation lead to technical difficulties: the autocorrelation structure of the RINAR(p) process is not identified, in spite of the similarity of their model with a real AR process and no asymptotic distribution is available for parameter estimates.

Contributions and relation to prior work. In the past several years, the study of the marginal distribution for INAR(1) (or GINAR(1)) type processes has become very popular. Indeed, by choosing appropriate distributions for the counting sequences and the innovation it is possible to construct a stationary first-order autoregressive integer-valued processes (INAR(1) or GI-NAR(1) processes) with discrete marginal distributions such as the Poisson, negative binomial, geometric, and many others (see [START_REF] Al-Zaid | First-order integer-valued autoregressive (INAR(1)) process[END_REF], [START_REF] Mckenzie | Some simple models for discrete variate time series[END_REF], [START_REF] Mckenzie | Some ARMA models for dependent sequences of Poisson counts[END_REF] and Ristic et al. (2009a,b)). For example, consider a INAR(1) stationary process, with Bernoulli(α) as common distribution of the counting sequences, and suppose that the stationary distribution of the process is a Poisson with λ/(1 -α) as parameter. Then, we can deduce that the noise ε has also a Poisson distribution with λ as parameter. This model is called Po-INAR(1) process.

To the best of our knowledge, the question of marginal distribution for models of time series on Z has often been avoided. Recently, [START_REF] Freeland | True integer value time series[END_REF] proposed the TINAR(1) process which can be interpreted as the difference between two independent Po-INAR(1) processes. Explicitly, let

Z t = X t -Y t , where X t = α • X t-1 + t , Y t = α • Y t-1 + η t , with t ∼ Poisson(λ 1 ) and η t ∼ Poisson(λ 2 ). Hence, the marginal distribution of Z t is a Skellam (λ 1 /(1 -α), λ 2 /(1 -α))
, which is the difference between two independent Poisson with respectively λ 1 /(1-α) and λ 2 /(1-α) as parameters. In particular, when λ 1 = λ 2 = λ, the marginal distribution of Z t is symmetric around zero on Z. Furthermore, note that with large values for µ 1 + µ 2 , a Skellam(µ 1 , µ 2 ) distribution can be sufficiently approximated by a gaussian distribution.

Indeed, the literature on classes of discrete distributions defined on Z is very limited. The Skellam distribution, originally introduced in [START_REF] Skellam | The frequency distribution of the difference between two poisson variates belonging to different populations[END_REF], is one of the most interesting member. Skellam distribution has found some applications recently (see Karlis andNtzoufras (2006, 2009) and Al-Zaid and Omair ( 2010)). Some other distributions defined as the difference of two discrete variables are given in [START_REF] Ong | A class of discrete distributions arising from difference of two random variables[END_REF]. One class of existing models are provided by the discretization of some continuous distribution such that the Gaussian distribution (see, e.g. [START_REF] Kemp | Characterizations of a discrete normal distribution[END_REF], [START_REF] Szablowski | Discrete normal distribution and its with jacobi theta functions[END_REF], [START_REF] Roy | The discrete normal distribution[END_REF] and [START_REF] Inusah | A skew Laplace distribution on integers[END_REF]).

In the present article we introduce a new class of distribution on Z, denoted by Rademacher(p) -N. This latter can be interpreted as the extension on Z of any discrete distribution defined on N, and defined as follows.

Definition 1 (Rademacher(p) -N class) Let p ∈ (0, 1), R and W be two independent random variables such that R ∼ Rademacher(p) i.e. P(R = 1) = p = 1 -P(R = -1), and support(W ) ⊆ N. We say that a random variable X belongs to the Rademacher(p) -N class, if and only if

X d = RW.
Note that, from this class and by specifying the distribution of W , we can extract many subclasses. Two of them are presented below.

Definition 2 (R-Po(p, λ) and SPo(λ) classes) Let p ∈ (0, 1) and λ > 0. We define the Rademacher-Poisson class, and we denote it by R-Po(p, λ), all the random variables X belonging to the Rademacher(p) -N class where W ∼ Poisson(λ). The R-Po(1/2, λ) is called the symmetric-Poisson class and denoted by SPo(λ).

The main purpose of this paper is to introduce and develop a new family of SINAR(1) process, by considering a parametric assumption on the common distribution of the counting sequence of the model. The new model will be denoted by P-SINAR(1). Thus, using the characteristic function, we give a full form of innovation distribution when the marginal distribution of the new process belongs to the Rademacher(p) -N class. Then, we extend the previous result by giving an implicit form of the innovation distribution when the marginal distribution of P-SINAR(1) is known (not necessarily a member of the Rademacher(p)-N class, e.g. Skellam). Furthermore, we obtain an implicit form of the marginal distribution of P-SINAR(1) process when the innovation distribution is known (e.g. symmetric Poisson, Skellam). Note that these results are very useful for the for prediction purposes and the goodness of fit of the P-SINAR(1) model. Namely, identifying the noise distribution allows to calculate the transition probability of the process, which is simply the one step ahead predictive distribution. On the other hand, knowing the marginal distribution of P-SINAR(1) allows to verify if given real data can be adjusted by the model.

Paper organization. The rest of this paper is organized as follows. In section 2, the first-order Parametric Signed INteger-valued AutoRegressive (P-SINAR(1)) model is presented together with some of its properties. Our main results are set in Section 3. A one-step ahead forecasting approach is described in Section 4. Simulation experiments as well as analysis of real data are discussed in Section 5 and 6. The proofs of all forthcoming results are postponed to Section 7.

The P-SINAR(1) process

The so-called signed thinning operator, originally proposed by [START_REF] Latour | An integer-valued bilinear type model[END_REF], seems to be the natural extension of the classical Steutel & Van Harn operator to Z-valued random variables and is defined as follows.

Definition 3 (Signed thinning operator ) Let (Y i ) i∈Z be a sequence of i.i.d. integer-valued random variables with F as common distribution, independent of an integer-valued random variable X. The signed thinning operator, denoted by F •, is defined by

F • X =            sign(X) |X| i=1 Y i , if X = 0, 0 otherwise, (1)
where for an integer x, sign(x) = 1 if x > 0 and -1 if x < 0. The sequence (Y i ) i∈Z is referred to as a counting sequence.

Based on the above operator, [START_REF] Kachour | A p-order signed integer-valued autoregressive (SINAR(1)) model[END_REF] introduced the SINAR (for Signed INteger-valued AutoRegressive) process, which can be considered as a natural extension on Z of the well-known INAR (or GINAR) process. In this paper, we focus on the SINAR(1) model defined as follows.

Definition 4 (SINAR(1) process) A sequence (X t ) t∈Z is said to be a SINAR(1) process if it admits the following representation

X t = F 1 • X t-1 + ε t , t ∈ Z, (2)
where (ε t ) t∈Z is a sequence of i.i.d. integer-valued random variables, with mean µ ε and variance σ 2 ε . It is understood that the counting sequence Y

(1) t associated to the operator F 1 • is independent of ε t . Moreover, ε t and the sigma-algebra

G t-1 = σ (X t-1 , • • • ), t ∈ Z, are supposed to be mutually independent.
Next, we summarize some theoretical results concerning stationarity ergodicity and autocorrelation structure of the SINAR(1) process (for details, we refer to [START_REF] Kachour | A p-order signed integer-valued autoregressive (SINAR(1)) model[END_REF]). First, we recall the following assumptions.

(A1) F 1 charges zero, i.e. F 1 {0} = P(Y 1 = 0) > 0.
(A1') There exist an even integer a ∈ Z and an odd integer b ∈ Z such that

P (ε 1 = a) > 0, P (ε 1 = b) > 0. (A2) We have |E(Y 1 )| < 1.
Suppose that (A1) (or (A1')) and (A2) hold. Then the equation ( 2) has a unique stationary solution X satisfying E(X 2 ) < ∞. Finally, note that the process X is also ergodic.

Remark 1 A stationary SINAR(1) process has the same autocorrelation structure than a classical AR(1) process

Y t = α 1 Y t-1 + η t , t ∈ Z,
where α 1 denotes the mean of the distribution F 1 (which is not necessarily positive) and η t denotes a white noise. In other words, we have cor

(X t , X t+k ) = α |k| 1 .
Now, we introduce a new family of SINAR(1) process, by specifying F 1 the common distribution of the counting sequence of the process. The new model will be denoted by P-SINAR(1) process (for Parametric Signed INteger-valued AutoRegressive).

Definition 5 (P-SINAR(1) process) A sequence (X t ) t∈Z is said to be a P-SINAR(1) process if it admits the same representation of a SINAR(1) process, i.e.

X t = F 1 • X t-1 + ε t , t ∈ Z, (3)
with F 1 the common distribution of the i.i.d. counting sequence (Y i ) i∈Z defined as follows:

P(Y 1 = -1) = (1 -α) 2 , P(Y 1 = 0) = 2α(1 -α), P(Y 1 = 1) = α 2 , ( 4
)
where α ∈ (0, 1).

Remark 2 We have Y 1 d = Z 1 -1, where Z 1 ∼ Binomial(2, α).
Moreover, for all k ∈ N * , we have

P k i=1 Y i = l = P(Z k = k + l), l ∈ {-k . . . , k}, (5)
where

Z k ∼ Binomial(2k, α).
In the next theorem, we derive some properties of the P-SINAR(1) process.

Theorem 1 Consider a P-SINAR(1) model, as defined in (3). Then, 1. The process has a unique stationary solution, denoted by X, and it is also ergodic.

2. (X t ) forms an homogeneous Markov chain with state space E = Z and its transition probability function is given by

P(X t = a | X t-1 = 0) = P(ε = a), a ∈ Z, (6)
and, for all a ∈ Z and b ∈ Z * ,

P(X t = a | X t-1 = b) = |b| l=-|b| P(ε = a -l)P Z |b| = |b| + sign(b)l , ( 7
)
where

Z |b| ∼ Binomial(2|b|, α).
Moreover, under second order stationarity conditions, we obtain

E(X t ) = µ ε 1 -α 1 , V(X t ) = β 1 E(|X t |) + σ 2 ε 1 -α 2 1 , (8)
where

α 1 = E(Y 1 ) = 2α -1, β 1 = V(Y 1 ) = 2α(1 -α) = P(Y 1 = 0).

Results

Main results

First, some recalls and notations are necessary. For any random variable U , let φ U be the characteristic function of U and G U its probability generating function i.e.

φ U (t) = E(e itU ), t ∈ R, G U (z) = E(z U ), |z| < 1,
Remark 3 When U has a distribution on Z, only the characteristic function is useful to characterize the distribution: it exists for all probability distributions and it determines uniquely the law of a random variable in the sense that φ U = φ V is equivalent to "U and V have the same distribution".

In Theorem 2 below, we give the form of the characteristic function of the noise ε d = ε 1 of a P-SINAR(1) process with a stationary distribution belonging to the Rademacher(p) -N class.

Theorem 2 Consider a P-SINAR(1) process as defined in (3). Suppose that its marginal stationary distribution, denoted by X, belongs to the Rademacher(p)-N class, i.e. X d = RW , where R and W are two independent random variables such that R ∼ Rademacher(p), p ∈ (0, 1), and support(W ) ⊆ N. Then, the characteristic function of the noise ε is

φ ε (t) = pφ W (t) + (1 -p)φ W (-t) pG W (φ Y (t)) + (1 -p)G W (φ Y (-t)) , t ∈ R, where φ Y (s) = 2α(1 -α) + α 2 e is + (1 -α) 2 e -is , s ∈ R.
Moreover, we have

P(ε = k) = 1 2π π -π φ ε (t)e -ikt dt, k ∈ Z. ( 9
)
Remark 4 Unfortunately, even if W is known (e.g. W ∼ Poisson(λ)), we can not confirm that ε has a known distribution of Z. Indeed, this is not a real surprise, since the known classes of distributions on Z are very limited and often such classes are rather built.

In the next theorem, we give an implicit form of the distribution of the noise ε of the P-SINAR(1) process (depending on the marginal stationary distribution of (3)).

Theorem 3 Consider a P-SINAR(1) process as defined in (3). Suppose that X its marginal distribution is known (not necessarily a member of the Rademacher(p) -N class). Then, we have

P(ε = k) = 1 1 -P(X = 0) (Q 1 (k) + Q 2 (k)), k ∈ Z, (10)
where

Q 1 (k) = ∞ m=1 P(X = m) k+m j=k-m m l=-m P(ε = j-l)P(Z m = m+j-k)P(Z m = m+l) (11) and Q 2 (k) = ∞ m=1 P(X = -m) k+m j=k-m m l=-m P(ε = j-l)P(Z m = m-j+k)P(Z m = m-l), (12) with Z m ∼ Binomial(2m, α).
Remark 5 Theorem 3 seems more general than Theorem 2, since we do not specifying the distribution of X. On the other hand, the implicit form of the distribution of the noise ε, given by the equations ( 10), ( 11) and ( 12), is difficult to exploit even numerically. However, supposing that the support of X is finite (which means that the stationary process (3) is an homogeneous Markov chain with a finite state space E = support(X)), the calculation becomes feasible. For example, let support(X) = {-M, . . . , M }, with M ∈ N * . From the definition of signed thinning operator (1) and under Assumption [H], we have support(ε) = {-2M, . . . , 2M }.

Therefore, ( 10), ( 11) and ( 12) lead to a general system of 2M linear equations with 2M unknowns. Moreover, since 2M k=-2M P(ε = k) = 1, this system has a unique solution (even if M is very large) and then the distribution of ε can be calculated.

In the sequel, we give an implicit form of the marginal stationary distribution X of the P-SINAR(1) process (depending on the distribution of the noise ε).

Theorem 4 Consider a P-SINAR(1) process as defined in (3). Suppose that the distribution of the noise ε is known. Then, the marginal stationary distribution, denoted by X, of (3) is given as follows. For all k ∈ Z, we have

P(X = k) = P(F 1 • X = 0)P(ε = k) + j∈Z * P(F 1 • X = j)P(ε = k -j), (13)
where, for all l ∈ Z,

P(F 1 • X = l) = ∞ k=1 (P(Z k = k + l)P(X = k) + P(Z k = k -l)P(X = -k)) 1 {l∈{-k,...,k}} + P(X = 0)1 {l=0} , with Z k ∼ Binomial(2k, α).
Remark 6 Note that the implicit form of the distribution of X given by the Theorem 4 is also difficult to handle. However, additional assumptions may reduce the calculation, e.g. a noise with a finite support, without making it really feasible.

Estimation

In this section, we describe estimation for the parameters of the P-SINAR(1) process under some parametric additional assumptions. Here, we recall that, from (4), we have

E(Y 1 ) = α 1 = 2α -1, V(Y 1 ) = β 1 = 2α(1 -α). ( 14
)
Assumption on X the marginal distribution of the process. Suppose that X follows a symmetric-Poisson distribution SPo(λ), i.e. X d = RW where R ∼ Rademacher(1/2) and W ∼ Poisson(λ). Since E(X) = 0 and from (8), we have that µ ε the mean of the noise ε equals zero. Let X 0 , X 1 , • • • , X n be observations from the P-SINAR(1) process with the previous assumption. Thus, the conditional least squares (CLS) estimator of α 1 (mean of F 1 ) is defined by

αn = n t=1 X t-1 X t n t=1 X 2 t-1 . ( 15
)
Moreover, from ( 14), one can deduce that the CLS estimator of α (the actual parameter of the distribution F 1 ), is given by αn = 1 + αn 2 .

Theorem 5 Consider a P-SINAR(1) process as defined in (3). Assume that X ∼ SPo(λ). Then, αn is strongly consistent, i.e.

lim n→∞ αn = α 1 , a.s. and √ n (α n -α 1 ) d -→ N (0, Ψ ), ( 16
)
where

Ψ = 2α(1 -α) 2λ 3 + 4λ 2 + 4λ + 1 λ(1 + λ) 2 .
On the other hand, since P(X = 0) = e -λ , the parameter λ of the marginal stationary distribution can be estimated by

λn = log n n i=1 1 {Xi=0} . ( 17
)
Note that, by using the ergodicity of the process (3), one can prove that λn is strongly consistent, and has a gaussian asymptotic behavior, i.e √ n λn -λ where μn denotes the CLS estimator of µ ε (the noise mean), which has also a gaussian asymptotic behavior and defined by μn = Xn,0 -αn Xn,1 , with Xn,j = (1/n) n t=1 X t-j , j ∈ {0, 1}.

Remark 8 Consider now that X follows a symmetric Skellam distribution, i.e. X = Z 1 -Z 2 , where Z 1 and Z 2 are i.i.d. with Z 1 ∼ Poisson(λ). The conditional least squares estimator αn still has gaussian asymptotic behavior. In this case, the asymptotic variance is more complicated to determine analytically, since it involves calculating E(|X|) and E(|X| 3 ), which can be deduced as in [START_REF] Katti | Moments of the Absolute Difference and Absolute Deviation of Discrete Distributions[END_REF] but there are not in a helpful form. On the other hand, since the marginal variance is 2λ, the parameter λ can be estimated by the sample variance as follows:

λn = 1 2n n t=1 X 2 t . ( 18
)
Assumption on the noise distribution. Suppose that ε in (3) follows a symmetric-Poisson distribution SPo(λ). Therefore, we have µ ε = 0 and σ 2 ε = λ(1 + λ). Hence, the conditional least squares will yield only an estimate for α 1 and not for λ. Let X 0 , X 1 , • • • , X n be observations from the P-SINAR(1) process with the previous assumption. Here, we will use the moment estimation which is simple and straightforward based on (8) and the fact that λ > 0. Let ρn (1) denotes the first order sample autocorrelation function, s 2 the empirical variance, and

| x| = 1 n n t=1 |X t |.
Then, it follows that

α1 = ρn (1), α = 1 + α1 2 , (19) and λ
= √ 1 + 4κ n -1 2 , ( 20
)
where

κ n = 2α(1 -α)(2s 2 -| x|). ( 21
)
Remark 9 Suppose now that ε follows a symmetric Skellam distribution SKellam(λ). Thus, we have µ ε = 0 and σ 2 ε = 2λ. Therefore, the moment estimates of the parameters are defined as follows: α1 = ρn (1),

and λ = κ n 2 , ( 22
)
where κ n is given by ( 21).

Remark 10 Since the distribution of ε is known, others estimation methods can be considered. Namely, the conditional maximum likelihood (CML) estimators of the parameters are obtained by maximizing the function

L θ (x) = n t=2 P(X t = x t | X t-1 = x t-1 ), where θ = (α 1 , λ), x = (x 1 , • • • , x n ) ∈ Z n
, and the conditional distribution X t | X t-1 is given by ( 6) and ( 7) of Theorem 1. Note that, the CML estimation is more complicated (but tractable) and their results are slightly better than those obtained by the moments estimation method described above.

Forecasting

In this section, we focus on one step ahead predictions for the P-SINAR(1) process, as defined in (3). The one-step least squares ahead forecast XT +1 of

X T +1 equals XT +1 = E(X T +1 | X T ) = α 1 X T + µ ε ,
where α 1 = 2α -1 denotes the mean of F 1 and µ ε that of the noise. In general XT +1 is real-valued. In order to get integer values we can use the median from the transition probability function of the process, defined in ( 6) and ( 7), together with its 95% confidence limits, i.e.

N j=M P(X T +1 = j | X T = i) ≥ 0.95
such as N -M to be as small as possible. Note that, the conditional distribution X T +1 | X T depends on the noise distribution.

Remark 11 Another mapping into the discrete support of the series can be obtained by rounding to the nearest integer. Thus, we find

XT +1 = α 1 X T + µ ε , (23)
where • denotes the rounding operator. Note that this method is not very efficient. For instance, suppose that we consider two P-SINAR(1) models with different centered noises (e.g. ε t ∼ SPo(λ) and η t ∼ Skellam(λ)) to fit a given real data. According to (23), in spite of their difference, both models have the same prediction performances.

Numerical illustration

In the sequel, we consider a P-SINAR(1) process, as defined in (3), with a fixed α = 0.8. It follows that, P(Y 1 = -1) = 0.04, P(Y 1 = 0) = 0.32, P(Y 1 = 1) = 0.64.

First numerical study: on the distribution of ε. We suppose that support(X) = {-2, -1, 0, 1, 2}, and its distribution is given by P(X = 0) = 0.7, P(X = 1) = P(X = -1) = 0.14,

P(X = 2) = P(X = -2) = 0.01. ( 24
)
According to Remark 5, we have support(ε) = {-4, . . . , 4} and its distribution is identified by solving the following system of equations: for all k ∈ {-4, . . . , 4},

P(ε = k) = 2 m=1 P(X = m) k+m j=k-m m l=-m P(ε = j -l)H(j, l, m, k) ,
where

H(j, l, m, k) = P(Z m = m + j -k)P(Z m = m + l) +P(Z m = m -j + k)P(Z m = m -l) =       2m m + l             2m m + j -k       G(j, l, m, k),
with G(j, l, m, k) = 0.8 2m+(j-k+l) 0.2 2m-(j-k+l) + 0.8 2m-(j-k+l) 0.2 2m+(j-k+l) .

In Figure 1, we depict the probability mass function of ε distribution for this case.

Remark 12

The histogram of the real data, studied in Section 6 below, is consistent with the specific distribution of X defined in (24).

Second numerical study: on the distribution of X. Now, we suppose that the distribution of ε is known. Hence, by Theorem 4, the marginal stationary distribution X of the process (3) has an implicit form given by ( 13). Next, we consider two types of distribution for ε:

-First, we suppose that ε follows a symmetric Skellam distribution distribution with parameter λ, i.e. X = Z 1 -Z 2 , where Z 1 and Z 2 are i.i.d. with Z 1 ∼ Poisson(λ). Recall that with large values for λ, a symmetric Skellam(λ) distribution can be sufficiently approximated by a centered gaussian distribution.

In Figure 2, we depict the probability mass function of the marginal stationary distribution of the P-SINAR(1) for λ ∈ {1, 2, 3, 7, 10, 30}. Hence, we remark that with large values for λ, X can be also sufficiently approximated by a centered gaussian distribution. -Now, we suppose that ε follows a symmetric-Poisson distribution. Recall that, Symmetric-Poisson is a subclass of Rademacher(p) -N class (see Definition 2). Moreover, for large values of λ, a R-Po(p, λ) distribution can be sufficiently approximated by two gaussian distributions, one on N and the other on Z -= Z \ N (details are given in Annex).

In Figure 3, we depict the probability mass function of the marginal stationary distribution of the SINAR(1) for λ ∈ {1, 2, 3, 7, 10, 30}. Thus, we find that X has a gaussian approximation with large values of λ.

Application to real data

The models introduced in this paper are applied to 597 counts of the price change for the Australian firm Broken Hill Proprietary (BHP) Limited, measured in cents, at each transaction time during the day of October 1, 2001. Note that, the observations vary from -2 to 2, the sample mean equals -0.00335 and the variance 0.41274.

In Figure 4, we have the time series plot, histogram, sample autocorrelation function and sample partial autocorrelation function. Thus, the correlogram shows that the first sample autocorrelation, is more significant than the others and there exists a cutt-off after lag 1 in the partial autocorrelations. That justify an AR(1) process fitting. Recall that a SINAR(1) process is defined on Z, has a correlation structure similar to a real AR(1), and allows negative autocorrelation. Thus, for the previous reasons, one can consider a SINAR(1) to analyse the data. On the other hand, according to the nature of the observations, assuming that the distribution F 1 is defined by (4) does not seem abusive. Therefore, in the sequel we will use a P-SINAR(1) to fit the data.

Furthermore, since the observations are centered, it is natural to consider a mean zero noise. For our first attempt to fit the data, we propose a P-SINAR(1) model such that the noise follows a symmetric Poisson distribution, i.e. ε ∼ SPo(λ 1 ). First, we estimated the model parameters using ( 19), ( 20) and ( 21). It follows that α1 = -0.2623 α = 1 + α1 2 = 0.3688 , and λ1 = 0.1859.

Recall that, according to Theorem 4, the marginal stationary distribution of the process has an implicit form given by ( 13), which is difficult to handle. Thus, in order to assess the goodness of fit we followed a computational approach. Explicitly, using the above estimated parameter values, we simulated 1000 series of length 597 and we kept the expected relative frequencies for each series. The reported frequencies are the mean over the 1000 series.

The Table 1 below contains the observed and expected relative frequencies based on the proposed model. One can see that the model fits very well the data. In Figure 5, we depict the one ahead predictive distribution X t+1 | X t = i for values of i = 1, -1, 2, -2. Now, we will calculate through simulation the standard errors of the moment estimation method considered to estimate the parameters of the above proposed model. Again, by using the estimated parameters, i.e. α = 0.3688 and λ = 0.1859, we simulate a serie of length 597, and we find the moment estimates. After K = 1000 independent replications, we calculate the average and standard deviation of {α i , i = 1, • • • , K} and { λi , i = 1, • • • , K} the sequences of the obtained estimates. This resulted in α = 0.3677 and λ1 = 0.1867 with standard errors of 0.0296 and 0.0257, respectively.

Note that, one can propose a P-SINAR(1) model with a symmetric Skellam distribution for the noise, i.e. ε ∼ Skellam(λ 2 ). Then, we estimated the model parameters using ( 19), ( 21) and ( 22). Hence, we obtain α1 = -0.2623, and λ2 = 0.1103. Again, in order to access the goodness of fit, we proceed in the same manner as the previous case.

Table 2 contains the observed and expected relative frequencies based on the new proposed model. One can see that the new model fits also very well the data. We depict the one ahead predictive distribution X t+1 | X t = i for values of i = 1, -1, 2, -2, in Figure 6.

Remark 13

The estimated values of λ 1 and λ 2 are quite close and weak (both < 0.2). Furthermore, remark that Skellam and Symmetric Poisson distributions have similar behavior with weak parameter values. This explain that the one ahead predictive distributions for both above models (displayed in Figure 5 and 6) are also very close. Finally, we note that the one step ahead prediction (for both models) will often leads to zero.

Remark 14 The histogram of the data (see Figure 4) is consistent with a symmetric Skellam distribution. Therefore, one can also propose a P-SINAR(1) process, where its marginal distribution X ∼ Skellam(λ 3 ), to fit the data. The parameters of the model, i.e. α 1 (or α) and λ 3 , can be estimated by ( 15) and ( 18). According to Theorem 3, the noise distribution of the process has an implicit form given by (10) ( 11) and ( 12), which is difficult to exploit analytically. However, using the data, one can simulated the noise distribution. Stationarity of P-SINAR(1). Recall that, we have

F 1 {0} = P(Y 1 = 0) = 2α(1 -α) > 0 and |E(F 1 )| = |2α -1| < 1.
Thus, assumptions (A1) and (A2) of Theorem 1 of [START_REF] Kachour | A p-order signed integer-valued autoregressive (SINAR(1)) model[END_REF] are verified (they have been recalled in Section 2), and then the P-SINAR(1) process (3) has a unique stationary solution, denoted by X, and it is also ergodic. Markovian properties of the process. Since the noise (ε t ) t∈Z in (3) is an i.i.d. sequence of integer-valued random variables, then the P-SINAR(1) process, denoted by (X t ) t∈Z , is Markovian. Thus, from the definition of the signed thinning operator, for all a ∈ Z, we have (with

ε d = ε t ) Γ (a, 0) = P(X t = a | X t-1 = 0) = P(ε = a). (25)
Furthermore, for any a ∈ Z and b ∈ Z * , we obtain

Γ (a, b) = P(X t = a | X t-1 = b) = P   sign(b) |b| k=1 Y k + ε t = a   = l∈Z P(ε = a -l)P   |b| k=1 Y k = sign(b)l   .
Recall that, from (5), for any k ∈ N * , we have

P k i=1 Y i = l = P(Z k = k + l), l ∈ {-k, . . . , k},
where Z k ∼ Binomial(2k, α). Therefore, we find

Γ (a, b) = |b| l=-|b| P(ε = a -l)P Z |b| = |b| + sign(b)l . ( 26
)
Theorem 1 is proved.

Proof of Theorem 2

Recall that for any random variable U , we denote by φ U its characteristic function and G U its probability generating function.

By the stochastic representation of the model (3), we have (with

X d = X t ) φ ε (t) = φ X (t) φ F1•X (t) . ( 27
)
Calculation of φ X (t). We have

φ X (t) = P(X = 0) + E e itX 1 {X>0} + E e itX 1 {X<0} .
Observe that

E e itX 1 {X>0} = ∞ k=1 e itk P(X = k) = p ∞ k=1 e itk P(W = k) = p ∞ k=0 e itk P(W = k) -P(W = 0) = pφ W (t) -pP(X = 0)
and

E e itX 1 {X<0} = -1 k=-∞ e itk P(X = k) = (1 -p) -1 k=-∞ e itk P(W = -k) = (1 -p) ∞ k=0 e -itk P(W = k) -P(W = 0) = (1 -p)φ W (-t) -(1 -p)P(X = 0). Therefore φ X (t) = pφ W (t) + (1 -p)φ W (-t). ( 28
)
Calculation of φ F1•X (t). By the definition of the signed thinning operator F 1 •, we have

φ F1•X (t) = P(X = 0) + E e it X i=1 Yi 1 {X>0} + E e -it -X i=1 Yi 1 {X<0} .
It follows from the i.i.d. nature of (Y i ) i∈Z , the independence of (Y i ) i∈Z and X, and

X d = RW that E e it X i=1 Yi 1 {X>0} = E 1 {X>0} E e it X i=1 Yi | X = ∞ k=1 (φ Y (t)) k P(X = k) = p ∞ k=1 (φ Y (t)) k P(W = k) = p ∞ k=0 (φ Y (t)) k P(W = k) -P(W = 0) = pG W (φ Y (t)) -pP(X = 0).
Similarly, we have

E e -it -X i=1 Yi 1 {X<0} = E 1 {X<0} E e -it -X i=1 Yi | X = -1 k=-∞ (φ Y (-t)) -k P(X = k) = (1 -p) -1 k=-∞ (φ Y (-t)) -k P(W = -k) = (1 -p) ∞ k=0 (φ Y (-t)) k P(W = k) -P(W = 0) = (1 -p)G W (φ Y (-t)) -(1 -p)P(X = 0). Hence φ F •X (t) = pG W (φ Y (t)) + (1 -p)G W (φ Y (-t)). (29)
Combining ( 27), ( 28) and ( 29), we obtain

φ ε (t) = pφ W (t) + (1 -p)φ W (-t) pG W (φ Y (t)) + (1 -p)G W (φ Y (-t)) , t ∈ R.
Finally, from the definition of the characteristic function, P(ε = k) represents the kth Fourier coefficient, and then (9) follows. Theorem 2 is proved.

Proof of Theorem 3

We have, for any k ∈ Z,

P(ε = k) = m∈Z j∈Z P (ε t = k, X t = j, X t-1 = m) = P (k) + Q(k),
where

P (k) = j∈Z P (ε t = k, X t = j, X t-1 = 0)
and

Q(k) = m∈Z * j∈Z P (ε t = k, X t = j, X t-1 = m) . Let Γ (a, b) = P(X t = a | X t-1 = b), (a, b) ∈ Z 2
, be the transition probability function of the process (3).

Observe that, thanks to (25) and the stochastic representation of the model

(3), (with

X d = X t ) P (k) = P(X t = k, X t-1 = 0) = Γ (k, 0)P(X = 0) = P(ε = k)P(X = 0). Therefore P(ε = k) = 1 1 -P(X = 0) Q(k).

Now remark that

Q(k) = m∈Z * j∈Z P   X t -sign(X t-1 ) |Xt-1| i=1 Y i = k | X t = j, X t-1 = m   × P(X t = j, X t-1 = m) = m∈Z * j∈Z P   sign(m) |m| i=1 Y i = j -k   Γ (j, m)P(X = m) = Q 1 (k) + Q 2 (k),
where, by ( 5) and ( 26), for Z m ∼ Binomial(2m, α),

Q 1 (k) = ∞ m=1 P(X = m) j∈Z P m i=1 Y i = j -k Γ (j, m) = ∞ m=1 P(X = m) k+m j=k-m P(Z m = m + j -k)Γ (j, m) = ∞ m=1 P(X = m) k+m j=k-m m l=-m P(ε = j -l)P(Z m = m + j -k)P(Z m = m + l) and Q 2 (k) = -1 m=-∞ P(X = m) j∈Z P -m i=1 Y i = -(j -k) Γ (j, m) = ∞ m=1 P(X = -m) j∈Z P m i=1 Y i = -(j -k) Γ (j, -m) = ∞ m=1 P(X = -m) k+m j=k-m P(Z m = m -j + k)Γ (j, -m) = m∈N * P(X = -m) k+m j=k-m m l=-m P(ε = j -l)P(Z m = m -j + k)P(Z m = m -l).
Theorem 3 is proved.

Proof of Theorem 4

Implicit form for the stationary distribution. By the stochastic representation of the model ( 3) and the independence between ε t and the sigma-algebra

G t-1 = σ (X t-1 , • • • ), for any k ∈ Z, we have (with X d = X t ) P(X = k) = j∈Z P(F 1 • X = j)P(ε = k -j) = P(F 1 • X = 0)P(ε 1 = k) + j∈Z * P(F 1 • X = j)P(ε = k -j).
Calculation of F 1 • X distribution. For all l ∈ Z, we have

P(F 1 • X = l) = k∈Z P(X = k, F 1 • X = l) = k∈Z P(F 1 • X = l | X = k)P(X = k) = k∈Z P   sign(k) |k| i=1 Y i = l   P(X = k) = A(l) + B(l) + C(l), where
A(l) = P(X = 0)1 {l=0} , by (5), for Z k ∼ Binomial(2k, α), B(l) = ∞ k=1 P k i=1 Y i = l P(X = k) = ∞ k=1 P(Z k = k + l)P(X = k)1 {l∈{-k,...,k}} and
C(l) = -1 k=-∞ P -k i=1 Y i = -l P(X = k) = ∞ k=1 P k i=1 Y i = -l P(X = -k) = ∞ k=1 P(Z k = k -l)P(X = -k)1 {l∈{-k,...,k}} .
Combining the values of A(l), B(l) and C(l), we end the proof of Theorem 4.

Proof of Theorem 5

Note that we consider a consider a P-SINAR(1) process, as defined in (3). Thus, the existence of a unique stationary solution, denoted by X, and the ergodicity of the process are ensured. Therefore, consistency and asymptotic normality of conditional least squares estimator follow from the result of [START_REF] Klimko | On conditional least squares estimation for stochastic process[END_REF], as it is explained in [START_REF] Du | The integer-valued autoregressive (INAR(p)) model[END_REF] for the case of INAR processes. Now, it remains to calculate the finite variance Ψ of the normal distribution in (16), under the assumption X ∼ SPo(λ), i.e. X d = RW , where R ∼ Rademacher(1/2) and W ∼ Poisson(λ). Hence, we note that

√ n (α n -α 1 ) = n -1 2 n t=1 X t-1 (X t -α 1 X t-1 ) n -1 n t=1 X 2 t-1
and further that the denominator converges almost surely to E(X 2 ) = λ(1+λ). So, by the law of large numbers and the continuous mapping theorem, this has the same asymptotic distribution as

1 λ(1 + λ) n -1 2 n t=1 X t-1 (X t -α 1 X t-1 ).
Since this sum is a martingale and martingale differences are stationary ergodic with finite moments, it converges to a normal distribution with mean zero and variance

Ψ = 1 λ 2 (1 + λ) 2 E(X 2 t-1 (X t -α 1 X t-1 ) 2 ) = 1 λ 2 (1 + λ) 2 E(X 2 t-1 V(X t | X t-1 )) = 1 λ 2 (1 + λ) 2 β 1 E(X 2 t-1 |X t-1 | + σ 2 ε EX 2 t-1 ) = 1 λ 2 (1 + λ) 2 β 1 E(X 2 |X|) + σ 2 ε E(X 2 ) ,
where α 1 = 2α -1 and β 1 = 2α(1 -α) denotes respectively the mean and variance of F 1 , and σ 2 ε the variance of ε. Now, from (8), we have

σ 2 ε = (1 -α 2 1 )V(X) -β 1 E(|X|) = β 1 (2V(X) -E(|X|)).
Recall that W ∼ Poisson(λ), and

V(X) = E(X 2 ) = E(W 2 ) = λ(1 + λ). Thus E(|X|)) = E(W ) = λ,
and then

σ 2 ε = 2αλ(1 -α)(2λ + 1). Note also that E(X 2 |X|) = E(W 3 ) = λ(λ 2 + 3λ + 1).
Therefore,

Ψ = 1 λ 2 (1 + λ) 2 β 1 E(W 3 ) + σ 2 ε E(W 2 ) = 2α(1 -α) 2λ 3 + 4λ 2 + 4λ + 1 λ(1 + λ) 2 .
Theorem 5 is proved.

Annex: on the behavior of a R-Po(p, λ) distribution with large λ

In the sequel, by sake of simplicity, we suppose that λ = n ∈ N * . Let p ∈ (0, 1) and X n = RW n , where R and W n are two independent random variables such that R ∼ Rademacher(p) and W n ∼ Poisson(n). Hence, for any x ∈ R, we have

P(X n ≤ x) = (1 -p)P(-W n ≤ x) + pP(W n ≤ x) = (1 -p)(1 -P(W n ≤ -x)1 {x≤0} ) + pP(W n ≤ x)1 {x≥0} = (1 -p) 1 -P W n √ n - √ n ≤ - x √ n - √ n 1 {x≤0} + pP W n √ n - √ n ≤ x √ n - √ n 1 {x≥0} .
By the central-limit theorem, for n large enough, one can approximate W n / √ n-√ n by a random variable U ∼ N (0, 1). Now, let f U be the density function of U and F U its cumulative distribution function. Therefore,

P(X n ≤ x) ≈ (1 -p) 1 -F U - x √ n - √ n 1 {x≤0} + pF U x √ n - √ n 1 {x≥0} = G n (x),
where G n (x) denotes a cumulative distribution function of density

g n (x) = G n (x) = (1 -p) 1 √ n f U - x √ n - √ n 1 {x≤0} + p 1 √ n f U x √ n - √ n 1 {x≥0} , x ∈ R.
Thus, one can deduce that with large values for λ, a R-Po(p, λ) distribution can be sufficiently approximated by two gaussian distributions: one on N and the other on Z -= Z \ N.

In Figure 7, we depict the probability mass function of a R-Po(1/2, λ) distribution (called also Symmetric-Poisson distribution and denoted by SPo(λ)) for λ ∈ {4, 6, 8, 10}. For each case, the blue curve represents the function g n ≡ g λ given above.

  d-→ N (0, Φ) where Φ > 0.Remark 7 Let X be a random variable following a R-Po(p, λ) distribution, i.e. X d = RW where R ∼ Rademacher(p) and W ∼ Poisson(λ). Then, the parameter λ can be again estimated by(17). Moreover, from (8), the parameter p can be estimated

Fig. 1

 1 Figures

Fig. 2 Fig. 3 Fig. 4 Fig. 5

 2345 Fig.2The probability stationary distribution X, when ε ∼ Skellam(λ) with different values of λ.

Fig. 6

 6 Fig.6Predictive distributions for Xt | X t-1 = i for values of i = 1, -1, 2, -2, with ε ∼ Skellam(λ).

Fig. 7

 7 Fig.7The probability distribution of a SPo(λ) with different values of λ.

Table 1

 1 Observed versus expected probabilities, when ε ∼ SPo

	Values	Observed Expected
	≤ -2	0.0134	0.0217
	-1	0.1507	0.1557
	0	0.67839	0.6700
	1	0.1407	0.1373
	≥ 2	0.0167	0.01507

Table 2

 2 Observed versus expected probabilities, when ε ∼ Skellam

	7 Proofs
	Proof of Theorem 1