Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
73,402 characters · 16 sections · 107 citation commands
Generalized Poisson Difference Autoregressive Processes
Abstract. This paper introduces a new stochastic process with values in the set $\mathbb{Z}$ of integers with sign. The increments of process are Poisson differences and the dynamics has an autoregressive structure. We study the properties of the process and exploit the thinning representation to derive stationarity conditions and the stationary distribution of the process. We provide a Bayesian inference method and an efficient posterior approximation procedure based on Monte Carlo. Numerical illustrations on both simulated and real data show the effectiveness of the proposed inference.\\
Keywords: Bayesian inference, Counts time series, Cyber risk, Poisson Processes.
MSC2010 subject classifications. Primary 62G05, 62F15, 62M10, 62M20.\\
In many real-world applications, time series of counts are commonly observed given the discrete nature of the variables of interest. Integer-valued variables appear very frequently in many fields, such as medicine (see CarRoyLam1999), epidemiology (see Zeg1988 and DavDunWan1999), finance (see LieMolPoh2006 and RydShe2003), economics (see Fre1998 and FreCab2004), in social sciences (see PedKar2011), sports (see ShaMoy2016) and oceanography (see CunVasBou2018). In this paper, we build on Poisson models, which is one of the most used model for counts data and propose a new model for integer-valued data with sign based on the generalized Poisson difference (GPD) distribution. An advantage in using this distribution relies on the possibility to account for overdispersed data with more flexibility, with respect to the standard Poisson difference distribution, a.k.a. Skellam distribution. Despite of its flexibility, GPD models have not been investigated and applied to many fields, yet. ShaMoy2014 proposed a GPD distribution obtained as the difference of two underling generalized Poisson (GP) distributions with different intensity parameters. They showed that this distribution is a special case of the GPD by Con1986 and studied its properties. They provided a Bayesian framework for inference on GPD and a zero-inflated version of the distribution to deal with the excess of zeros in the data. ShaMoy2016 showed empirically that GPD can perform better than the Skellam model.
As regards to the construction method, two main classes of models can be identified in the literature: parameter driven and observation driven. In parameter-driven models the parameters are functions of an unobserved stochastic process, and the observations are independent conditionally on the latent variable. In the observation-driven models the parameter dynamics is a function of the past observations. Since this paper focuses on the observation-driven approach, we refer the reader to MacZuc1997 for a review of parameter-driven models.
Thinning operators are a key ingridient for the analysis of observation-driven models. The mostly used thinning operator is the binomial thinning, introduced by SteVHa1979 for the definition of self-decomposable distribution for positive integer-valued random variables. In mathematical biology, the binomial thinning can be interpreted as natural selection or reproduction, and in probability it is widely applied to study integer-valued processes. The binomial thinning has been generalized along different directions. Lat1998 proposed a generalized binomial thinning where individuals can reproduce more than once. KimPar2008 introduced the signed binomial thinning, in order to allows for negative values. Joe1996 and ZheBasDat2007 introduced the random coefficient thinning to account for external factors that may affect the coefficient of the thinning operation, such as unobservable environmental factors or states of the economy. When the coefficient follows a beta distribution one obtain the beta-binomial thinning (McK1985, McK1986 and Joe1996). AlOAly1992, proposed the iterated thinning, which can be used when the process has negative-binomial marginals. AlzAlO1993 introduced the quasi-binomial thinning, that is more suitable for generalized Poisson processes. ZhaWanZhu2010 introduced the signed generalized power series thinning operator, as a generalization of KimPar2008 signed binomial thinning. Thinning operation can be combined linearly to define new operations such as the binomial thinning difference (Fre2010) and the quasi-binomial thinning difference (CunVasBou2018). For a detailed review of the thinning operations and their properties different surveys can be consulted: MacZuc1997, KedFok2005, McK2003, Wei2008, ScoWeiGou2015. In this paper, we apply the quasi-binomial thinning difference.
In the integer-valued autoregressive process literature, thinning operations have been used either to define a process, such as in the literature on integer-valued autoregressive-moving average models (INARMA), or to study the properties of a process, such as in the literature on integer-valued GARCH (INGARCH). INARMA have been firstly introduced by McK1986 and AlOAlz1987 by using the binomial thinning operator. DuLi1991 extended to the higher order $p$ the first-order INAR model of AlOAlz1987. KimPar2008 introduced an integer-valued autoregressive process with signed binomial thinning operator, INARS($p$), able for time series defined on $\mathbb{Z}$. AndKar2014introduced SINARS, that is a special case of INARS model with Skellam innovations. In order to allow for negative integers, Fre2010 proposed a true integer-valued autoregressive model (TINAR(1)), that can be seen as the difference between two independent Poisson INAR(1) process. AlzAlO1993 have studied an integer-valued ARMA process with Generalized Poisson marginals while AlzOma2014 proposed a Poisson difference INAR(1) model. CunVasBou2018 firstly applied the GPD distribution to build a stochastic process. The authors proposed an INAR with GPD marginals and provided the properties of the process, such as mean, variance, kurtosis and conditional properties.
RydShe2000 introduced heteroskedastic integer-valued processes with Poisson marginals. Later on, Hei2003 introduced an autoregressive conditional Poisson model and Fer2006 proposed the INGARCH process. Both models have Poisson margins. Zhu2012 defined a INGARCH process to model overdispersed and underdispersed count data with GP margins and AloAlzOma2018 proposed a Skellam model with GARCH dynamics for the variance of the process. KooLitLuc2014 proposed a Generalized Autoregressive Score (GAS) Skellam model. In this paper, we extend Fer2006 and Zhu2012 by assuming GPD marginals for the INGARCH model, and use the quasi-binomial thinning difference to study the properties of the new process.
Another contribution of the paper regards the inference approach. In the literature, maximum likelihood estimation has been widely investigated for integer-valued processes, whereas a very few papers discuss Bayesian inference procedures. CheLee2016 introduced Bayesian zero-inflated GP-INGARCH, with structural breaks. ZhuLi2009 proposed a Bayesian Poisson INGARCH(1,1) and Chen2016 a Bayesian Autoregressive Conditional Negative Binomial model. In this paper, we develop a Bayesian inference procedure for the proposed GPD-INGARCH process and a Markov chain Monte Carlo (MCMC) procedure for posterior approximation. One of the advantages of the Bayesian approach is that extra-sample information on the parameters value can be included in the estimation process through the prior distributions. Moreover, it can be easily combined with a data augmentation strategy to make the likelihood function more tractable.
We apply our model to a cyber-threat dataset and contribute to cyber-risk literature providing evidence of temporal patterns in the mean and variance of the threats, which can be used to predict threat arrivals. Cyber threats are increasingly considered as a top global risk for the financial and insurance sectors and for the economy as a whole EIOPA2. As pointed out in Hass16, the frequency of cyber events substantially increased in the past few years and cyber-attacks occur on a daily basis. Understanding cyber-threats dynamics and their impact is critical to ensure effective controls and risk mitigation tools. Despite these evidences and the relevance of the topic, the research on the analysis of cyber threats is scarce and scattered in different research areas such as cyber security Agra18, criminology Bren04, economics Anderson610 and sociology. In statistics there are a few works on modelling and forecasting cyber-attacks. Xu2017 introduced a copula model to predict effectiveness of cyber-security. Wer17 used an autoregressive integrated moving average model to forecast the number of daily cyber-attacks. Edwards2015HypeAH apply Bayesian Poisson and negative binomial models to analyse data breaches and find evidence of over-dispersion and absence of time trends in the number of breaches. See Hus2018 for a review on modelling cyber threats.
The paper is organized as follows. In Section 2 we introduce the parametrization used for the GPD and define the GPD-INGARCH process. Section 3 aims at studying the properties of the process. Section 4 presents a Bayesian inference procedure. Section 5 and 6 provide some illustration on simulated and real data, respectively. Section 7 concludes.
A random variable $X$ follows a Generalized Poisson (GP) distribution if and only if its probability mass function (pmf) is
with parameters $\theta >0$ and $0\leq \lambda < 1$ Con1986. We denote this distribution with $GP(\theta,\lambda)$. Let $X\sim GP(\theta_{1}, \lambda)$ and $Y\sim GP(\theta_{2}, \lambda)$ be two independent GP random variables, Con1986 showed that the probability distribution of $Z=(X-Y)$ follows a Generalized Poisson Difference distribution (GDP) with pmf:
where $z$ takes integer values in the interval $(-\infty, +\infty)$ and $0<\lambda <1$ and $\theta_{1},\theta_{2} >0$ are the parameters of the distribution. See Appendix (ref) for a more general definition of the GPD with possibly negative $\lambda$.
In the following Lemma we state the convolution property of the GPD distribution since which will be used in this paper. Appendix (ref) provides an original proof of this result.
We use an equivalent pmf and a re-parametrization of the GPD, which are better suited for the definition of a INGARCH model. A random variable $Z$ follows a GPD if and only if its probability distribution is
We denote this distribution with $GPD(\mu,\sigma^2,\lambda)$.
The mean, variance, skewness and kurtosis of a GDP random variable can be obtained in close form by exploiting the representation of the GDP as difference between independent GP random variables.
See Appendix (ref) for a proof.\\
Figure (ref) shows the sensitivity of the probability distribution with respect to the location parameter $\mu$ (panel a), the scale parameter $\sigma^{2}$ (panel b) and the skewness parameter $\lambda$ (different lines in each plot). For given values of $\lambda$ and $\mu$, when $\sigma^{2}$ decreases the dispersion of the GPD increases (dotted and dashed lines in the right plot). For given values of $\lambda$ and $\sigma^{2}$, the distribution is right-skewed for $\mu=8$, which corresponds to $S(Z)=0.7155$, and left-skewed for $\mu=-4$, which corresponds to $S(Z)=-0.3578$, (dotted and dashed lines in the left plot). See Appendix (ref) for further numerical illustrations.
Differently from the usual GARCH$(p,q)$ process (e.g., see FraZak2019), the INGARCH$(p,q)$ is defined as an integer-valued process $\lbrace Z_{t}\rbrace_{t \in \mathbb{Z}}$, where $Z_{t}$ is a series of counts. Let $\mathcal{F}_{t-1}$ be the $\sigma$-field generated by $\lbrace Z_{t-j}\rbrace_{j\geq 1}$, then the GPD-INGARCH$(p,q)$ is defined as
with
where $\tilde{\mu}_{t-j}=\mu_{t-j}(1-\lambda)$, $\alpha_{0}\in \mathbb{R}$, $\alpha_{i}\geq 0$, $\beta_{j}\geq 0$, $i=1,\ldots,p$, $p\geq 1$, $j=1,\ldots, q$, $q \geq 0$. For $q=0$ the model reduces to a GPD-INARCH$(p)$ and for $\lambda=0$ one obtains a Skellam INGARCH$(p,q)$ which extends to Poisson differences the Poisson INGARCH$(p,q)$ of Fer2006. From the properties of the GPD, the conditional mean $\mu_t=E(Z_{t}\vert \mathcal{F}_{t-1})$ and variance $\sigma^2_t=V(Z_{t}\vert \mathcal{F}_{t-1})$ of the process are:
respectively. In the application, we assume $\sigma^{2}_{t}=\vert \mu_{t} \vert \phi$. Following the parametrization defined in Remark (ref), we need to impose the constrain $\phi >(1-\lambda)^{-2}$, in order to have a well-defined GPD distribution. In Fig. (ref), we provide some simulated examples of the GPD-INGARCH$(1,1)$ process for different values of the parameters $\alpha_{0}$, $\alpha_{1}$ and $\beta_{1}$.
Simulations from a GPD-INGARCH are obtained as differences of GP sequences
where
Each random sequence is generated by the branching method in Fam1997, which performs faster than the inversion method for large values of $\theta_{1t}$ and $\theta_{2t}$. We considered two parameter settings: low persistence, that is $\alpha_{1}+\beta_{1}$ much less than 1 (first column in Fig. (ref)) and high persistence, that is $\alpha_{1}+\beta_{1}$ close to 1 (second column in Fig. (ref)). The first and second line show paths for positive and negative value of the intercept $\alpha_{0}$, respectively. The last line illustrates the effect of $\lambda$ on the trajectories with respect to the baselines in Panels (a) and (b). Comparing (I.a) and (I.b) in Fig. (ref) one can see that increasing $\beta_{1}$ increases serial correlation and the kurtosis level (compare (II.a) and (II.b)).
We provide a necessary condition on the parameters $\alpha_{i}$ and $\beta_{j}$ that will ensure that a second-order stationary process has an INGARCH representation. First define the two following polynomials: $D(B)=1-\beta_{1}B-\ldots -\beta_{q}B^{q}$ and $G(B)=\alpha_{1}B+\ldots+\alpha_{p}B^{p}$, where B is the backshift operator. Assume the roots of $D(z)$ lie outside the unit circle. For non-negative $\beta_{j}$ this is equivalent to assume $D(1)=\sum_{j=1}^{q}\beta_{j}<1$. Then, the operator $D(B)$ has inverse $D^{-1}(B)$ and it is possible to write
where $H(B)=G(B)D^{-1}(B)=\sum_{j=1}^{\infty} \psi_{j}B^{j}$ and $\psi_{j}$ are given by the power expansion of the rational function $G(z)/D(z)$ in the neighbourhood of zero. If we denote $K(B)=D(B)-G(B)$ we can write the necessary condition as in the following proposition.
We study the properties of the process by exploiting a suitable thinning representation following the strategy in Fer2006 and Zhu2012 for Poisson and Generalized Poisson INGARCH, respectively. We use the quasi-binomial thinning as defined in Wei2008 and the thinning difference (CunVasBou2018) operators.
We show that the INGARCH process can be obtained as a limit of successive approximations. Let us define:
and
where $\lbrace U_{1t} \rbrace_{t \in \mathbb{Z}}$ and $\lbrace U_{2t} \rbrace_{t \in \mathbb{Z}}$ are sequences of independent GP random variables and for each $t \in \mathbb{Z}$ and $i \in \mathbb{N}$, $\lbrace V_{1t,i,j} \rbrace_{j \in \mathbb{N}}$ and $\lbrace V_{2t,i,j} \rbrace_{j \in \mathbb{N}}$ represent two sequences of independent integer random variables. Moreover, assume that all the variables $U_{s}$, $V_{t,i,j}$, with $s \in \mathbb{Z}$, $t \in \mathbb{Z}$, $i \in \mathbb{N}$ and $j \in \mathbb{N}$, are mutually independent.
It is possible to show that $X_{t}^{(n)}$ and $Y_{t}^{(n)}$ have a thinning representation. We define a suitable thinning operation, first used by AlzAlO1993 and follow the notation in Wei2008, let $\rho_{\theta,\lambda} \circ$ be the quasi-binomial thinning operator, such that it follows a QB($\rho$,$\theta /\lambda$,$x$).
Both $X_{t}^{(n)}$ and $Y_{t}^{(n)}$ in Eq. (ref) and (ref) admit the representation
and
where $\varphi \circ X$ is the quasi-binomial thinning operation. See Appendix (ref) for a definition.
In the following we introduce the thinning difference operator and show that $Z_{t}^{(n)}=X_{t}^{(n)}-Y_{t}^{(n)}$ has a thinning representation.
See CunVasBou2018 for an application of the thinning operation to GPD-INAR processes and Appendix (ref) for further details. Using the new operator as defined in Eq. (ref), we can represent $Z_{t}^{(n)}$ as follows.
The proposition above shows that $Z_{t}^{(n)}$ is obtained through a cascade of thinning operations along the sequence $\lbrace U_{t} \rbrace_{t \in \mathbb{Z}}$. For example:
Since $Z_{t}^{(n)}$ is a finite weighted sum of independent GPD random variables, the expected value and the variance of $Z_{t}^{(n)}$ are well defined. Moreover, it can be seen that $E[Z_{t}^{(n)}]$ does not depend on $t$ but only on $n$, hence it can be denoted as $\mu_{n}$. Using Proposition (ref) and $\mu_{k}=0$ if $k<0$, it is possible to write $\mu_{n}$ as follows
from which it follows $D(B)\mu_{n}=G(B)\mu_{n} + \alpha_{0}\Leftrightarrow K(B)\mu_{n}=\alpha_{0}$, where $K(B)=D(B)-G(B)$. From the last equation it can be seen that the sequence $\lbrace \mu_{n}\rbrace$ satisfies a finite difference equation with constant coefficients. The characteristic polynomial is $K(z)$ and all its roots lie outside the unit circle if $K(1)>0$. Under the assumption $K(1)>0$, the following holds true.
Given Proposition (ref), if we can show that $\lbrace Z_{t}^{(n)} \rbrace$ is a strictly stationary process, for any given $n$, then also its almost sure limit $\lbrace Z_{t} \rbrace_{t \in \mathbb{Z}}$ will be a strictly stationary process. In order to show stationarity for $\lbrace Z_{t}^{(n)} \rbrace$, we follow a procedure similar to the one in Fer2006. Let us define the probability generating function (pgf) $g_{\mathbf{W}}(\mathbf{t})$ of the random vector $\mathbf{W}=(W_{1},\ldots,W_{k})$
where $p(\mathbf{W})=Pr(\mathbf{W}=(W_{1}, \ldots, W_{k})')$ and $\mathbf{t}=(t_{1}, \ldots, t_{k})' \in \mathbb{C}^{k}$. The probability generating function has the following properties.
Using the probability generating function, in the following we know the stationarity of the process.
To verify that the distributional properties of the sequence are satisfied, we will follow the same arguments in Fer2006 adjusted for our sequence. Given $\mathcal{F}_{t-1}=\sigma (\lbrace Z_{u}\rbrace_{u\leq t-1})$, for $t \in \mathbb{Z}$, let
The sequence $\lbrace \mu_{t} \rbrace$ satisfies
Moreover, recalling that $Z_{t}=X_{t}-Y_{t}$, for a fixed t, we can consider three sequences, $\lbrace r_{1t}^{(n)} \rbrace_{n \in \mathbb{N}}$, $\lbrace r_{2t}^{(n)} \rbrace_{n \in \mathbb{N}}$ and $\lbrace r_{t}^{(n)} \rbrace_{n \in \mathbb{N}}$, defined by
and
As claimed by Fer2006, there is a subsequence $\lbrace n_{k} \rbrace$ such that $r_{t}^{(n_{k})}$ converges almost surely to $Z_{t}$. We know that
and
Since $X_{t}^{(n)} \overset{a.s.}{\longrightarrow} X_{t}$ and $Y_{t}^{(n)} \overset{a.s.}{\longrightarrow} Y_{t}$, we know that the first term in both Eq. (ref) and (ref) goes to zero. Therefore, we can write
and, as before, $(Z_{t} - Z_{t}^{(n)})$ goes to zero since we have proven almost sure convergence.\\ We have now to show that the second term in the last line of Eq. (ref) goes to zero, for this purpose we need to find a sequence
that converges almost surely to zero. For this reason it is more suitable to rewrite the previous sequence as follows
Fer2006 show that
therefore, we can conclude that also
Equation (ref) implies that $W_{t}^{(n)}$ converges to zero in $L^{1}$, therefore there exist a subsequence $W_{t}^{(n_{k})}$ converging almost surely to the same limit. From this it follows directly that the distributional properties of $X_{t}$ are satisfied.
Since $r_{1t}^{(n_{k})} \overset{a.s.}{\longrightarrow} X_{t}$ and $r_{2t}^{(n_{k})} \overset{a.s.}{\longrightarrow} Y_{t}$, it is also true $r_{t}^{(n_{k})} \overset{a.s.}{\longrightarrow} Z_{t}$. Hence,
However,
and from Zhu2012 we know that both $r_{1t}^{(n)}$ and $r_{2t}^{(n)}$ have a Generalized Poisson distribution. Since the difference of two GP distributed random variables is GPD distributed, we can write
and conclude that
The conditional mean and variance of the process $Z_{t}$ are
where $\phi = \frac{1}{1-\lambda}$.
The unconditional mean and variance of the process are
From Th. 1 in Wei2009 we know a set of equations from which the variance and autocorrelation function of the process can be obtained. Suppose $Z_{t}$ follows the INGARCH(p,q) model in Eq. (ref) with $\sum_{i=1}^{p} \alpha_{i} + \sum_{j=1}^{q} \beta_{j} < 0$. From Th. 1 part (iii) in Wei2009, the autocovariances $\gamma_{Z}(k) = Cov[Z_{t},Z_{t-k}]$ and $\gamma_{\mu}(k) = Cov[\mu_{t},\mu_{t-k}]$ satisfy the linear equations
In order to have an explicit expression for the variance of $\mu_{t}$ and $Z_{t}$ and for the autocorrelations, we consider two special cases as in Zhu2012 and Wei2009. For a proof of the results in these examples, see Section (ref).
We propose a Bayesian approach to inference for GPD-INGARCH, which allows the researcher to include extra-sample information through the prior choice and allows us to exploit the stochastic representation of the GPD and the use of latent variables to make more tractable the likelihood function.
We assume the following prior distributions. A Dirichlet prior distribution for $\boldsymbol{\varphi} =(\alpha_{1},\ldots, \alpha_{p},\beta_{1},\ldots,\beta_{q})$, $\boldsymbol{\varphi} \sim Dir_{d+1}(c)$, with density:
where $\varphi_{i} \geq 0$ and $\sum_{i=1}^{d} \varphi_{i} \leq 1$. Panel (a) in Fig. (ref) provides the level sets of the joint density function of $\alpha_{1}$ and $\beta_{1}$ with hyper-parameters $c_{0}=3$, $c_{1}=4$ and $c_{2}=3$. We assume a flat prior for $\alpha_{0}$, i.e. $\pi(\alpha_{0}) \propto \mathbb{I}_{\mathbb{R}}(\alpha_{0})$. For $\lambda$ and $\phi$ we assume a joint prior distribution with uniform marginal prior $\lambda \sim \mathcal{U}_{[0,1]}$ and shifted gamma conditional prior $\phi \sim \mathcal{G}a^{*}(a,b,c)$, with density function:
where $c=(1-\lambda)^{-2}$. Panel (b) provides the level sets of the joint density function of $\phi$ and $\lambda$, with hyper-parameters $a=b=5$. The joint prior distribution of the parameters will be denoted by $\pi(\boldsymbol{\theta})=\pi(\boldsymbol{\varphi})\pi(\alpha_{0})\pi(\lambda)\pi(\phi)$.
Denote the probability distribution of $Z_{t}$ with
with $\underline{s}=\max(0,-z)$. Since the posterior distribution
is not analytically tractable we apply Markov Chain Monte Carlo (MCMC) for posterior approximation in combination with a data-augmentation approach TanWon1987. See RobCas2013 for an introduction to MCMC. As in KarNtz2006, we exploit the stochastic representation in Eq. (ref) and introduce two GPD latent variables $X_{t}$ and $Y_{t}$ with pmfs
Let $Z_{1:T}=(Z_{1},\ldots ,Z_{T})$, $X_{1:T}=(X_{1},\ldots ,X_{T})$ and $Y_{1:T}=(Y_{1},\ldots ,Y_{T})$. The complete-data likelihood becomes
where $\delta (z-c)$ is the Dirac function which takes value 1 if $z=c$ and 0 otherwise. The joint posterior distribution of the parameters $\boldsymbol{\theta}$ and the two collections of latent variables $X_{1:T}$ and $Y_{1:T}$ is
We apply a Gibbs algorithm RobCas2013 with a Metropolis-Hastings (MH) steps. In the sampler, we draw the latent variables and the parameters of the model by iterating the following steps:
where $\boldsymbol{\theta}_{-\eta}$ indicates the collection of parameters excluding the element $\boldsymbol{\eta}$.
The full conditional for the latent variables is
We draw from the full conditional distribution by MH. Differently from KarNtz2006, we use a mixture proposal distribution which allows for a better mixing of the MCMC chain. At the $j$-th iteration, we generate a candidate $X_{t}^{*}$ from $GP(\theta_{1t},\lambda)$ with probability $\nu$ and $(X_{t}^{*}-Z_{t})$ from $GP(\theta_{2t},\lambda)$ with probability $1-\nu$, and accept with probability
where $q(X_{t})=\nu f(X_{t} \vert \theta_{1t}, \lambda) + (1-\nu) f(X_{t}-Z_{t} \vert \theta_{2t}, \lambda)$ and $X_{t}^{(j-1)}$ is the $(j-1)$-th iteration value of the latent variable $X_{t}$. The method extends to the GPD the technique proposed in KarNtz2006 for the Poisson differences.
As regards to the parameter $\boldsymbol{\varphi}$, its full conditional distribution is
We consider a MH with Dirichlet independent proposal distribution
where $\textbf{c}^{*} = ({c}_{0}^{*},{c}_{1}^{*}, {c}_{2}^{*})$ and acceptance probability
The full conditional distribution of $\phi$ is
We consider the change of variable $\zeta = \log(\phi -c)$ with Jacobian $\exp(\zeta)$ and a MH step with a random walk proposal
where $\zeta_{j-1} = \log (\phi_{j-1}-c)$, $\phi_{j-1}$ is the previous iteration value of the parameter and $c=\frac{1}{(1-\lambda)^{2}}$. The acceptance probability is
where $\phi^{*}=c+\exp(\zeta^{*})$.
The full conditional distribution of $\lambda$ is
We consider a MH step with Beta random walk proposal
where $s$ is a precision parameter. The acceptance probability is:
The purpose of our simulation exercises is to study the efficiency of the MCMC algorithm presented in Section (ref). We evaluated the Gew92 convergence diagnostic measure (CD), the inefficiency factor (INEFF)\footnote{The inefficiency factor is defined as
where $\rho(k)$ is the sample autocorrelation at lag $k$ for the parameter of interest and are computed to measure how well the MCMC chain mixes. An INEFF equal to $n$ tells us that we need to draw MCMC samples $n$ times as many as uncorrelated samples.} and the Effective Sample Size (ESS).
We simulated 50 independent data-series of 400 observations each. We run the Gibbs sampler for 1,010,000 iterations on each dataset, discard the first 10,000 draws to remove dependence on initial conditions, and finally apply a thinning procedure with a factor of 250, to reduce the dependence between consecutive draws.
As commonly used in the GARCH and stochastic volatility literature Chib2002, casarin2009online, MCO12, Casetal2019, we test the efficiency of the algorithm in two different settings: low persistence and high persistence. The true values of the parameters are: $\alpha=0.25$, $\beta=0.23$, $\lambda=0.4$ in the low persistence setting and $\alpha=0.53$, $\beta=0.25$, $\lambda=0.6$ in the high persistence setting. Table (ref) shows, for the parameters $\alpha$, $\beta$ and $\lambda$, the INEFF, ESS and ACF averaged over the 50 replications before (BT subscript) and after thinning (AT subscript).
The thinning procedure is effective in reducing the autocorrelation levels and in increasing the ESS, especially in the high persistence setting. The p-values of the CD statistics indicate that the null hypothesis that two sub-samples of the MCMC draws have the same distribution is accepted. The efficiency of the MCMC after thinning generally improved. On average, the inefficiency measures (19.05), the p-values of the CD statistics (0.18) and the acceptance rates (0.35) achieved the values recommended in the literature Robetal1997.
Data in this application are the number of accidents near Schiphol airport in The Netherlands during 2001 (Fig. (ref)). They have been previously considered in Brietal2008 and AndKar2014. The time series of accident counts is non-stationary and should be differentiated KimPar2008. We applied our Bayesian estimation procedure, as described in Section (ref).
In Fig. (ref) are presented the histograms for the Gibbs draws for each parameters. Table (ref) presents the parameter posterior mean and standard error and the 95% credible interval for the unrestricted INGARCH(1,1) model (model $\mathcal{M}_1$). In the data, we found evidence of high persistence in the expected accident arrivals, i.e. $\hat{\alpha}+\hat{\beta}=0.8673$ and heteroskedastic effects, i.e. $\hat{\beta}=0.4753$. Also, there is evidence in favour of overdispersion, $\hat{\lambda}=0.5892$ and overdispersion persistsence $\hat{\phi}=179.7905$. We study the contribution of the heteroskedasticy and persistence by testing some restrictions of the INGARCH(1,1) (models from $\mathcal{M}_2$ to $\mathcal{M}_4$ in Tab. (ref)).
Bayesian inference compares models via the so-called Bayes factor, which is the ratio of normalizing constants of the posterior distributions of two different models (see CamPet2014 for a review). MCMC methods allows for generating samples from the posterior distributions which can be used to estimate the ratio of normalizing constants.
In this paper we use the method proposed by Gey1994. The method consists in deriving the normalizing constants by reverse logistic regression. The idea behind this method is to consider the different estimates as if they were sampled from a mixture of two distributions with probability
to be generated from the $j$-th distribution of the mixture. Gey1994 proposed to estimate the log-Bayes factor $\kappa=\eta_{2}-\eta_{1}$ by maximizing the quasi-likelihood function
where $n$ is the number of MCMC draws for each model and $X_{ij}=\log f(Z_{1:T},X_{1:T}^{(i)},Y_{1:T}^{(i)} \vert \boldsymbol{\theta}^{(i)})$ is the log-likelihood evaluated at the $i$-th MCMC sample for each model of Tab. (ref).
We performed six reverse logistic regressions, in which we compare pairwise our models. The approximated logarithmic Bayes factors $BF(\mathcal{M}_i,\mathcal{M}_j)$ are given in Tab. (ref). It is possible to see that our GPD-INGARCH$(1,1)$, $\mathcal{M}_{1}$, is preferable with respect to the other models. Notice that $\mathcal{M}_{2}$ corresponds to an INGARCH$(1,1)$ where the observations are form a standard Poisson-difference model PD-INGARCH$(1,1)$, $\mathcal{M}_{3}$ corresponds to an autoregressive model, GPD-INARCH$(1,0)$, whereas $\mathcal{M}_{4}$ is a standard Poisson difference augoregressive model, PD-INARCH$(1,0)$.
According to the Financial Stability Board FSB1, a cyber incident is any observable occurrence in an information system that jeopardizes the cyber security of the system, or violates the security policies and procedures or the use policies. Over the past years there have been several discussions on the taxonomy of incidents classification ENISA, in this paper we use the classification provided in the \citetalias{Hackweb} dataset. \citetalias{Hackweb} is a well-known cyber-incident website that collects public reports and provides the number of cyber incidents for different categories of threats: crimes, espionage and warfare.
Figure (ref) shows the total and category-specific number of cyber attacks at a daily frequency from January 2017 to December 2018. Albeit limited in the variety of cyber attacks the dataset covers some relevant cyber events and is one of the few publicly available datasets Agra18. The daily threats frequencies are between 0 and 12 which motivates the use of a discrete distribution. We remove the upward trend by considering the first difference and fit the GPD-INGARCH model proposed in Section (ref).
We applied our estimation procedure, as described in Section (ref). As in the previous application, we fix $\alpha_{0}=1.05$ that is coherent with the conditional mean of the time series. We ran the Gibbs sampler for 110000 iterations, where we discarded the first 10000 iterations as burn-in sample. In Fig. (ref) are presented the histograms for the Gibbs draws for each parameters.
Figure (ref) shows that, as before, it is reasonable to fit a GPD-INGARCH process to the difference of cyber attacks since both the autoregressive parameter $\alpha$ and $\beta$, that represent the heteroskedastic feature of the data, are different from zero. Additionally, the value of $\lambda$ suggest the presence of over-dispersion in the data.
Given the importance of forecasting cyber-attacks, in this section we present the results of one-step-ahead forecasting exercise over a period of 120. We follow an approach based on predictive distributions which quantifies all uncertainty associated with the future number of attacks and is used in a wide range of applications Cab05,Cab11. We account for parameter uncertainty and approximate the predictive distribution by MCMC. At tht $j$-th MCMC iteration we draw $Z_{T+h}^{(j)}$ from the conditional distribution given past observations and the parameter draw $\boldsymbol{\theta}^{(j)}$
where $j=1,\ldots,J$, denotes the MCMC draw, $h=1,\ldots,H$ the forecasting horizon and
We introduce a new family of stochastic processes with values in the set of integers with sign. The increments of the process follow a generalized Poisson difference distribution with time-varying parameters. We assume a GARCH-type dynamics, provide a thinning representation and study the properties of the process. We provide a Bayesian inference procedure and an efficient Monte Carlo Markov Chain sampler for posterior approximation. Inference and forecasting exercises on accidents and cyber-threats data show that the proposed GPD-INGARCH model is well suited for capturing persistence in the conditional moments and in the over-dispersion feature of the data.