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.
54,051 characters · 8 sections · 67 citation commands
Time-Varying Poisson Autoregression
In the last years time series of counts has been experiencing an increasing interest among researchers and practitioners. Applications on count data span over several different areas, including finance fokianos2011log, corporate defaults Agosto2016, environmental economics Wang2014, epidemiology davis2003, COVID-19 infections and deaths khismatullina2020nonparametric,li2021will, sports angelini2017parx,koopman2019. Models for time series of counts are generally based on the assumption that the observations follow a Poisson distribution conditional on past observations and explanatory variables, see davis2016 for an overview of recent developments on the econometric models for discrete values. In particular, davis2003 propose an observation-driven Poisson model for which they derive the conditions for stationarity and ergodicity and the large-sample properties of the MLE. They apply their method to daily counts of asthma presentations at a Sydney hospital. fokianos2009poisson develop and establish the asymptotic theory of a Poisson autoregression (PAR) model where the conditional distribution of the observations $\{y_t\}_{t \in \mathbb{Z}}$ given the intensity process $\{\lambda_t\}_{t \in \mathbb{Z}}$ is Poisson with parameter $\lambda_t$ and they apply this new model to describe the number of transactions per minute of the stock Ericsson B. Agosto2016 generalise the PAR model by including a set of exogenous covariates in the dynamic specification of the intensity process (PARX) to improve the model performance for forecasting the number of aggregate corporate default counts in US. Wang2014 propose a self-excited threshold structure in the PAR model specification by considering a two regime-switching approach according to the magnitude of the lagged observations, in order to model the number of major earthquakes in the world.
In this paper we contribute to this literature by introducing a new class of Poisson AutoRegressive models allowing for time-variation in the intensity and in the parameters which is particularly suitable to capture several (non-standard) features, including structural breaks, regime-switching and non-linearities, overdispersion and negative correlation, thus providing better fitting and forecasting performance than the original PARX model. The extension to the time-varying parameter specification of the PARX model, labeled TV-PARX model in what follows, is achieved using a score-driven approach. Originally proposed by creal2013 and Harvey2013 to estimate the evolution of time-varying parameters, the Generalized Autoregressive Score (GAS) framework is now a popular methodology to develop time-varying econometrics models. Blasques2019 propose a new class of GAS models for continuous time series allowing for time-varying parameters and provide the theoretical foundations for their acceleration method by showing optimality in terms of Kullback-Leibler divergence.
The literature on GAS models for time series of counts has been growing remarkably in the very recent years. gorgi2020beta uses the beta–negative binomial distribution for the narcotics trafficking reports in Sydney and the euro–pound sterling exchange rate, koopman2019 use the bivariate Poisson distribution for the number of goals in football matches and the Skellam distribution for the score difference, and gorgi2018 uses the Poisson distribution as well as the negative binomial distribution for the number of offensive conduct reports in the city of Blacktown, Australia.
The TV-PARX model we propose generalizes the accelerated score-driven model by Blasques2019 in the framework of time series of counts and allowing for the inclusion of exogenous covariates and deterministic components in the model specification. We provide conditions for stationarity and ergodicity, as well as asymptotic results for the MLE. Moreover, we show the relevance of our modeling approach by focusing on a Monte Carlo simulation study and {two} main applications. First, we analyze the dynamics of new COVID-19 confirmed cases in Italy and, as in li2021will, we found that the time-varying model specification is required. Indeed, the model parameters substantially change over time due to the different `waves' of the coronavirus infections and their association with the introduction of government's measures to control the spread of the virus, e.g. lockdown, curfew, and mass vaccination. Second, we study the dynamics of US corporate default counts. Using the data set analyzed by Agosto2016, we show that the time-varying specification of the PARX model allows to better capture and forecast the number of corporate defaults and to disentangle the main drivers of the US recessions.
The paper is organized as follows. Section (ref) describes the new TV-PARX model and discusses the properties of TV-PARX. Section (ref) presents the MLE for the TV-PARX. Section (ref) shows the finite sample properties and Section (ref) reports the empirical illustrations on COVID-19 data and on corporate defaults. Section (ref) concludes the paper. Details on notation and all the proofs are provided in (ref).
Consider a time series of counts $\{ y_t \}_{t \in \mathbb{Z}}$ with a Poisson conditional distribution,
where $\mathcal{F}_{t-1}= \sigma\{ y_{t-1}, y_{t-2}, y_{t-3}, \dots \}$ denotes the filtration generated by the count process $\{y_{t-i}\}_{i>0}$ up to time $t-1$, whereas $\lambda_t := \mathbb{E}[y_t | \mathcal{F}_{t-1}]$ is the conditional mean of the count process, which is allowed to vary over time. Examples are the Poisson autoregressive models by fokianos2009poisson and Agosto2016. While these papers consider a time-invariant parameter approach, we assume that the evolution of the conditional mean is described by a time-varying coefficient model representation using the score-driven framework $-$ see e.g.\ creal2013 and Harvey2013 for a detailed review $-$ thus allowing for time variation in the parameters of the dynamic equation for $\lambda_t$. Koopman2008 show that a score-driven model for Poisson time series of counts encompasses most of the observation-driven models considered by davis2003. Then, our contribution is to generalize the score-driven model for the intensity process and formally derive the stochastic properties of our Poisson score-driven model and the asymptotic theory of the maximum likelihood estimator. In particular, we consider a score-driven model with the observation equation in (ref), which implies that the conditional $\log$-density of $y_t$ is given up to a constant by
We propose the following nonlinear Poisson autoregressive model with exogenous covariates and deterministic components defined as
where $x_t \in \mathbb{R}^m$ denotes $m$ exogenous covariates, whereas in $d_t \in \mathbb{R}^d$ are stacked the deterministic components, e.g.\ seasonal dummies, impulse and step dummies, trend functions, etc. In the following, the model in (ref) is called Time-Varying Poisson AutoRegressive with eXogenous covariates, TV-PARX. Indeed, while $(\omega, \beta) \in\mathbb{R}^2$ and $ \psi \in \mathbb{R}^d$ are unknown static parameters, the key feature of model (ref) is that $\alpha_t$ and $\gamma_t$ are time-varying.
Similarly to Blasques2019, we characterize $\alpha_t$ by the following updating equation
where $\delta_\alpha$, $\phi_\alpha$ and $\kappa_\alpha$ are treated as fixed and unknown parameters. We note that, analogously to the models for continuous time series considered in Blasques2019, the law of motion of $\alpha_t$ is driven by the product of current and past innovations. Moreover, by scaling the score innovations with the inverse of information quantity $\mathcal{I}_{t} = \lambda_t$, Harvey2013 shows that the law motion for $\log \lambda_t$ is given by the first-order autoregressive process in (ref) in the case of no exogenous covariates ($x_t=0$) and no deterministic components ($d_t=0$). With the same rationale of $\alpha_t$, we define the score-driven-based updating of the parameters for the exogenous covariates by
where $\delta_\gamma$, $\phi_\gamma$ and $\kappa_\gamma$ are fixed and unknown parameters. Allowing for time-varying coefficients of the exogenous covariates is particularly useful to capture dynamic effects of the covariates on the conditional mean of the count process. In particular, this feature of the model enables a much better fit of the structural breaks and time-varying causality due to the occurrence of exogenous shocks, which is particularly recurrent in economic and financial phenomena.
Note also that in (ref) we adopt an exponential link function, which allows the conditional mean $\lambda_t$ to be always positive, i.e.\ there is no need to impose further restrictions to have $\lambda_t>0$, for any $t$. Moreover, there is no need to transform the $m$ exogenous covariates in $x_t$ to make them positive as done, e.g., in Agosto2016 (see Section (ref) for an empirical illustration).
{\bf Remark 1.} Several models discussed in the literature are special cases of the model proposed here. For instance, the model proposed by Harvey2013 is a special case of the model in (ref) with no covariates and deterministics, i.e.\ $x_t=0$ and $d_t=0$, and where the parameter $\alpha$ is time-invariant. Also, the resulting model is equivalent to the Poisson model of davis2003 for the case of time-invariant coefficients of the exogenous covariates. $\square$
Despite the time-varying specification of the model in (ref)-(ref) is in the same vein of Blasques2019, i.e.\ the coefficient $\alpha$ is allowed to change at each time period by means of the score-driven approach, with respect to Blasques2019 our proposal is developed to model time series of counts instead of continuous responses and we can include $m$ exogenous covariates in the model specification. Moreover, our proposal extends the Poisson autoregressive model by fokianos2009poisson and static PARX model by Agosto2016 by allowing for time-varying parameters, including the coefficients for the exogenous covariates.
We complete this section by obtaining sufficient conditions for geometric ergodicity of the Poisson autoregressive model described above. As usual in the literature, let $\{N_t(\,\cdot \,) \}_{t\in\mathbb{Z}}$ denotes a sequence of independent Poisson processes with unity intensity, so that the process $\{ y_t \}_{t\in\mathbb{Z}}$ in (ref) can be expressed as
where $\lambda_t$ follows the dynamics given in (ref) and (ref), see e.g.\ fokianos2009poisson. We begin by introducing the first result.
{\bf Remark 2.} In line with the $\log$-linear Poisson autoregressive model of fokianos2009poisson, by considering the full model given by equations (ref) and (ref) instead of the restricted on with no covariates and dummy variables $x_t$ and $\delta_t$, respectively, does not crucially affect the conditions for stationarity and ergodicity stated in Proposition (ref). In fact, since $\delta_t$ are deterministic, as soon as the exogenous covariate processes $\{ x_t \}_{t\in \mathbb{Z}}$ are Markov chains, we only need to retrive a set of separate conditions for their transition mechanism, together with need the additional assumption that $|\phi_\gamma|<1$. \\
It is well-known that when dealing with nonlinear time series models it is usually not easy to establish clear and/or simple stationarity conditions. However, we note that the sufficient conditions given in Proposition (ref) are remarkably easy to check. Moreover, a direct consequence of Proposition (ref) is that if the process $\{ \log \lambda_t \}_{t\in\mathbb{Z}}$ is initialized at its stationary distribution, then the process $\{ y_t \}_{t\in \mathbb{Z}}$ will also be stationary and ergodic with a finite first moment $\mathbb{E}[y_t] = \mathbb{E}[\lambda_t] < \infty$; see davis2003.
We conclude this section with a new proposition which gives sufficient conditions for $\mathbb{E}[y_t^k]< \infty$, with $k$ is a strictly positive number.
It is worth to note that Proposition (ref) we do not need stricter contraction conditions in order to obtain the finiteness of the higher-order moments. This feature of our model is particularly useful since it allow us to derive bounds for proving the asymptotic properties of the MLE, which is the main argument of the next Section, without imposing further restrictions on the data generating process.
The TV-PARX model can be easily estimated by standard ML, since the predictive log-likelihood is available in closed form. In the following, for the sake of simplicity, we consider the case of no exogenous covariates and no deterministic components, i.e.
which corresponds to (ref) where $x_t=0$ and $d_t=0$, and is labeled the TV-PAR model. Note that in the case where $x_t \neq 0$ and $d_t \neq 0$, the asymptotic theory developed for our MLE below, could be adapted straightforwardly using partial likelihood theory, see Wong1986.
Denote the parameter vector $\boldsymbol{\theta} = (\boldsymbol{\xi}^\prime, \boldsymbol{\psi}^\prime)^\prime \in \boldsymbol{\Theta} \subset \mathbb{R}^{5}$, where $\boldsymbol{\xi} = (\omega, \beta)^\prime$, $\boldsymbol{\psi} = (\delta_\alpha, \phi_\alpha, \kappa_\alpha)^\prime$, and $\boldsymbol{\Theta}$ is a compact parameter space. The true value of the combined parameter vector is denoted by $\boldsymbol{\theta}_0$, and we assume that $\{y_t\}_{t=1}^T$ is generated according to the TV-PAR process described by equations (ref), (ref) and (ref), evaluated at $\boldsymbol{\theta}_0$. The $\log$-likelihood function is given by
where
and
Note that the time-varying parameters in (ref) and (ref) are obtained recursively by using some fixed starting values $\hat{\lambda}_{1}(\boldsymbol{\theta})\in \mathbb{R}_+$, $\hat{\alpha}_1(\boldsymbol{\theta})\in \mathbb{R}$ and the observations $\{y_t\}_{t=1}^T$.
The MLE $\hat{\boldsymbol{\theta}}_T$ of $\boldsymbol{\theta}$ is defined as
The conditions stated in Proposition (ref) implies that $\alpha_t$ and $\log \lambda_t$ have stationary representations. Now, for the likelihood analysis and the asymptotic properties of the MLE we need to derive the stochastic limit properties of the filtered parameters $\{ \hat{\lambda}_t(\boldsymbol{\theta}) \}_{t \in \mathbb{N}}$ and $\{ \hat{\alpha}_t(\boldsymbol{\theta}) \}_{t \in \mathbb{N}}$, which can be seen as two stochastic functions. In particular, because of the initializations $\hat{\lambda}_{1}(\boldsymbol{\theta})$ and $\hat{\alpha}_{1}(\boldsymbol{\theta})$, the filtered parameters in (ref) and (ref), and therefore the approximate $\log$-likelihood function in (ref), are in general non-stationary. Therefore, in the next proposition, we derive sufficient conditions, which ensure that the effect that $\hat{\lambda}_{1}(\boldsymbol{\theta})$ and $\hat{\alpha}_{1}(\boldsymbol{\theta})$ have on the approximate $\log$-likelihood function $\{ \hat{L}_T (\boldsymbol{\theta}) \}_{t \in \mathbb{N}}$ vanish almost surely and exponentially fast (e.a.s.) uniformly over $\boldsymbol{\Theta}$. This phenomenon is well-known in the literature of nonlinear time series models, which is usually referred to the notion of invertibility, see Straumann_Mikosch2006. Stationarity and ergodicity of $\{\log \lambda_t \}$ (as established in Proposition (ref)), and Proposition (ref) below, are the key ingredients for establish asymptotic properties of the MLE.
Condition (ref) is a sufficient condition which ensures the invertibility of the TV-PARX model by an application of Proposition 5.2.12 in Straumann2005, or equivalently Theorem 2.8 in Straumann_Mikosch2006. It is different from the conditions for stationarity and ergodicity given in Theorem (ref), and moreover, it is only concerned with the filtering recursion relative to $\lambda_t(\boldsymbol{\theta})$ because the stationarity condition imposed in Proposition (ref) for $\alpha_t(\boldsymbol{\theta})$ is also sufficient for filter invertibility. The complex form of the contraction condition in (ref) is due to the exponential function. A similar problem with related discussion can be found in Wintenberger2013 for the EGARCH(1,1) model. Nevertheless, despite its sufficiency, this condition it is still important as it underline that the invertibility region is not degenerate.
Finally, by slightly reinforcing the contraction condition in (ref) as in (ref), Proposition (ref) also ensures the existence of an arbitrary large number of unconditional moments for the stationary and ergodic intensity process $\{\lambda_t(\boldsymbol{\theta})\}_{t\in\mathbb{Z}}$. This is a crucial property which will be useful for proving the consistency and asymptotic normality of the MLE.
To derive the asymptotic properties of the MLE we follow the classic theory in white_1994. First, for consistency we consider the almost sure uniform convergence of the approximate $\log$-likelihood function in (ref). In particular, Lemma (ref) below ensures that the average $\log$-likelihood $\hat{L}_T(\boldsymbol{\theta})$ converges uniformly to a function $L(\boldsymbol{\theta})$.
Second, we verify the identifiability of the model parameterization, which is the argument of the next Lemma.
As a consequence, given the results obtained in Lemma (ref) and Lemma (ref), an application of Theorem 2.11 of white_1994 implies the strong consistency of the MLE.
Next, we discuss the asymptotic distribution of the MLE. Consider the normalized score evaluated at $\boldsymbol{\theta}_0$
where $e_t := (y_t - \lambda_t(\boldsymbol{\theta}_0))\lambda^{-1}_t(\boldsymbol{\theta}_0)$, $\lambda_t := \lambda_t(\boldsymbol{\theta}_0)$ and $\log\lambda^{\boldsymbol{\theta}_0}_t := \frac{\partial\log\lambda_t(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}}|_{\boldsymbol{\theta}=\boldsymbol{\theta}_0}$. Since $\boldsymbol{\theta} = (\boldsymbol{\xi}^\prime, \boldsymbol{\psi}^\prime)^\prime$, we have that $ \log\lambda^{\boldsymbol{\theta}_0}_t=( (\log\lambda^{\boldsymbol{\xi}_0}_t)^\prime, (\log\lambda^{\boldsymbol{\psi}_0}_t)^\prime )^\prime$, where $\log\lambda^{\boldsymbol{\xi}_0}_t = \left( \frac{\partial\log\lambda_t(\boldsymbol{\theta}_0)}{\partial \omega}^\prime,\frac{\partial\log\lambda_t(\boldsymbol{\theta}_0)}{\partial \beta}^\prime \right) ^\prime$ and $\log\lambda^{\boldsymbol{\psi}}_t = \alpha^{\boldsymbol{\psi}}_t e_t$, with $\alpha^{\boldsymbol{\psi}_0}_t = \left(\frac{\partial\alpha_t(\boldsymbol{\theta}_0)}{\partial \delta_{\alpha}}^\prime, \frac{\partial\alpha_t(\boldsymbol{\theta}_0)}{\partial \phi_{\alpha}}^\prime, \frac{\partial\alpha_t(\boldsymbol{\theta}_0)}{\partial \kappa_{\alpha}}^\prime \right)$, as defined in Lemma (ref) in the Appendix.
Asymptotic normality then follows by a standard central limit theorem for martingale difference sequences, see e.g.\ Corollary 3.1 of Hall1980, by showing that
where
and, the Lindberg condition holds, i.e., $\forall \epsilon > 0$,
These conditions are verified in the following proposition.
Denote the information matrix at the true value by $\boldsymbol{J}$, i.e.
Then, we finally state the following theorem.
{\bf Remark 3.} It is worth noting that the results obtained for the MLE $\hat{\boldsymbol{\theta}}_T$ in this section hold true even if the conditional distribution is misspecified. This is a consequence of the fact that the Poisson distribution belongs to the linear exponential family. In particular, as proved by white1982regularity and Gourieroux1984, under high-level assumptions the MLE is a QMLE and maintain the same properties, the only change will occur on the asymptotic variance-covariance matrix of $\sqrt{T}(\hat{\boldsymbol{\theta}}_T - {\boldsymbol{\theta}}_0)$, which will assume the classic sandwitch form, that is
where $\boldsymbol{\Sigma} = \boldsymbol{J}^{-1}\boldsymbol{I}\boldsymbol{J}^{-1}$, with $\boldsymbol{J}$ as in equation (ref) and $\boldsymbol{I}= \mathbb{E}\Big[ \frac{\mathbb{V}[y_t|\mathcal{F}_{t-1}]}{\lambda^2_t(\boldsymbol{\theta}_0)}\frac{\partial \lambda_t(\boldsymbol{\theta}_0)}{\partial \boldsymbol{\theta}} \frac{\partial \lambda_t(\boldsymbol{\theta}_0)}{\partial \boldsymbol{\theta}^\prime}\Big]$. See also Ahmad2016 or Aknouche2021.
In this section we investigate the empirical performance in finite samples of the TV-PAR model compared to the standard PAR model. In particular, the aim of this Monte Carlo exercise is to investigate how quickly the $\lambda_t$ can adapt to (structural) changes. The time series are generated from the following DGP
where $N_t(\cdot)$ denotes a sequence of independent Poisson processes with unity intensity. Similarly to Blasques2019, for the time-varying $\lambda_t^0$, we consider a deterministic step-function
where $\delta \in \{0,2,4 \}$ and $\gamma \in \{1, 1.5, 2\}$. In this Monte Carlo exercise we generate $m=1000$ time series of sample sizes $T \in \{250, 500, 1000\}$. Figure (ref) reports the time varying intensity $\lambda_t^0$ for all the combinations of $\delta$ and $\gamma$ considered. For each series we estimate the time-varying PAR in (ref) and the static PAR, i.e.\ where $\alpha_t=\alpha$ in (ref). Table (ref) reports the Root Mean Squared Errors (RMSE) for all the cases considered. The results in Table (ref) show that the RMSE of the TV-PAR is smaller for all DGPs considered, except for, unsurprinsigly, the case of no steps. However, in the case of no steps, i.e.\ when $\delta=0$, the difference in RMSE between the time-varying and static models is relatively small in all the cases considered.
Conversely, when the $\delta>0$, the RMSE of the TV-PAR model is smaller than the one of the standard PAR model for the all cases considered. In particular, the higher the $\delta$, the better the performance of TV-PAR model compared to the standard PAR model, regardless of the sample size considered. Therefore, with more convoluted DGP dynamics, the time-varying model specification works better than to its static counterpart.
To illustrate, we report in Figure (ref) the cases of the DGPs with $\delta =4$ and $\gamma = 1$ for $T = 250, 500,$ and 1000. This figure shows the mean $\hat{\lambda}_t$ and the corresponding 95% confidence bands for the TV-PAR (in green) and the static PAR (in red) models. From Figure (ref), it clearly emerges that the time-varying model specification provides a much faster reaction and better fit to the jumps created by the step-function in (ref).
In this section, we show the usefulness of the time-varying specification of the TV-PARX model by considering two empirical illustrations. Section (ref) deals with the daily counts of infections of SARS-COV-2 virus in Italy until the end of May 2021. Section (ref) analyzes the monthly corporate default counts in the US as previously analyzed by Agosto2016.
The COVID-19 pandemic, started in 2019 in the city of Wuhan in China's region of Hubei and spread all around the world in 2020, has caused a tremendous public health emergency and a huge negative shock on the business cycle. Since then, researchers of various fields have started to work on that topic in order to help policy makers in their decisions. Are the number of cases (deaths) predictable? Are the evolution of the number of cases (deaths) similar across countries? These questions have increased the interest among researchers and practitioners for new models suited for time series of counts. Epidemiologists, for instance, use reduced-form models to produce forecasts of the number of infections (or deaths), see for instance batista2020estimation for a logistic growth model and giordano2020modelling and hethcote2000mathematics for reviews on several variants of the Susceptible, Infectious and Removed (SIR) model.
Due to the convoluted time evolution of the COVID-19 pandemic (see the number of new cases depicted in blue in Figure (ref)), the time-varying specification of the model is expected to capture the different waves of the pandemic and its non-standard dynamics are likely to outperforms much better time-invariant alternatives. Among others, li2021will found that the model parameters change both over time and across countries. Country-specific behaviour and local, transient outbreaks of COVID-19 is also found by khismatullina2020nonparametric and dong2020time using non-parametric approaches and by palmer2021count using Bayesian hierarchical models.
Therefore, we fit the TV-PARX model in (ref)-(ref) with $\gamma_t=0$ to the daily counts of confirmed new cases in Italy from the beginning of the COVID-19 pandemic outbreak (March 2020) to the end of May 2021, which involves three or four main waves of cases (see Figure (ref)). The counts dynamics in Figure (ref) shows a seasonal component mainly due to reporting. Hence, we consider the following TV-PARX model
where $d_t$ denotes seasonal (daily) dummies.
For comparison purposes, we also fit a time-invariant PARX model where $\phi_\alpha = \kappa_\alpha=0$ in (ref). Table (ref) shows the estimates and Table (ref) reports the differences in the log-likelihoods, information criteria and Root Mean Squared Errors (RMSE) of the estimated models. In terms of in-sample performance, these results show that the time-varying specification of the TV-PARX model is preferred over its time-invariant counterpart. Moreover, albeit both the time-varying model (TV-PARX) and its time-invariant counterpart (PARX) provide a good fit to the observed data, from Figure (ref) it emerges that the TV-PARX model approximates better the different waves of the COVID-19 pandemic in Italy. In particular, the estimated time-varying parameter $\hat{\alpha}_t$ ranges from around -0.2 to 0.8, and the minimum value virtually corresponds to the end of the first wave (end of June 2020), while the maximum values are observed at specific time points which are closely related to the beginning of each wave, e.g.\ March 2020, September-November 2020 and January 2021, as well as abrupt changes in $\hat{\alpha}_t$ for the last wave begun in March 2021.
In this section, we stress the usefulness of the time-varying specification of the TV-PARX model compared to the standard approach based on the constant parameter version of the PARX model by analyzing the count time series of corporate defaults in US. Agosto2016 contributed to the literature on modeling and forecasting time series of corporate defaults, by exploring the stylized fact that defaults tend to cluster over time and investigating the causes of this phenomenon. In particular, using a (constant parameter) PARX model they discriminate between comovements in corporate solvency caused by common underlying macroeconomic and financial risk factors, i.e.\ “systematic risk”, and feedback effects, i.e.\ “contagion”, where, conditionally on the common risk factors, the current number of defaults affects the probability of other firms' future insolvency. As pointed out in Agosto2016, one limitation of their approach is that the measure of contagion ignores feedback effects from defaults to covariates; that is, the case when $x_t$ is affected by past defaults, which in turn will affect future defaults. This issue is softened by the time-varying model specification of the TV-PARX as the time-dependent parameters may allow to capture the feedback effects from past defaults to the covariates. The data set on default counts analyzed by Agosto2016 consists of the number of bankruptices among Moody's rated industrial firms in the US collected at monthly frequency between the 1982 and 2011, for a total of $T = 360$ observations, from the Moody's Credit Risk Calculator (CRC). The blue line in Figure (ref) shows the default count time series. As discussed in Agosto2016, Figure (ref) discloses that default counts is a process characterised by (i) high persistence, (ii) clusters over time, and (iii) overdispersion. These three stylized facts are satisfactorily captured by the PARX model. However, the peak of defaults observed in 2008 (Figure (ref)) as well as the high increase in bankruptices around the year 2000, possibly interpretable as structural breaks, may be better captured by the time-varying specification of the TV-PARX model, which can locally allow for explosive roots rather than standard PARX model where the stationary condition must hold. Therefore, we consider the TV-PARX model, described in (ref)-(ref), with no deterministic componenents, i.e.\ $d_t=0$. For the sake of comparison, in the TV-PARX model specification we include the same covariates considered in the PARX model in Agosto2016.\footnote{Note that, with the specification in logs of the TV-PARX model, we are able to consider the whole time series of LI, and not only the negative values (in modulus to ensure non-negativity) as in Agosto2016. } In particular, we consider $x_t=(RV_t, LI_t)^\prime$, where $RV_t$ is the realized volatility computed as a proxy of the S&P 500 monthly realized volatility using daily squared returns, and $LI_t$ is the Leading Index released by the Federal Reserve. More formally, we consider the following TV-PARX model
We provide an analysis for the full sample 1982-2011 and compare the results for the time-invariant PARX model specification. For comparison purposes, we also fit a time-invariant PARX model where $\phi_{\alpha}=\kappa_{\alpha}=\phi_{\gamma_1}=\kappa_{\gamma_1}=\phi_{\gamma_2}=\kappa_{\gamma_2}=0$ in (ref). The maximum likelihood estimates are reported in Table (ref). According to the information criteria and the RMSE reported in Table (ref), the TV-PARX model provides a better in-sample fit than the PARX model, except for the BIC.\footnote{It is interesting to note that our TV-PARX model outperforms the best model specification in Agosto2016, namely (non-logs) PARX(2,1).} As can be noted from Figure (ref), both models capture the actual default counts dynamics well, but the goodness of fit of the TV-PARX model is superior during the periods where the number of defaults is larger, especially around the year 2009. Figure (ref) shows the filtered dynamic parameters, $\alpha_t$ and $\gamma_t$. From this figure, it can be noted that the estimated time-varying parameters are able to adjust according to the observed default count time series. Interestingly, the time-varying specification allows us to determine the main driver of a specific crisis. In particular, the dynamics of the estimated parameters $\gamma_{1,t}$ and $\gamma_{2,t}$ in Figure (ref) suggests that financial volatility is the main driver of the spike of the defaults in the US during the 2008-09 crisis, while during the early 90s recession the main driver is mainly macroeconomic. Therefore, the time-varying version of the PARX model proposed in (ref)-(ref) provides a more flexible alternative to the time-invariant counterpart.
In this paper, we have developed a novel class of score-driven models for time series of counts, allowing for time-varying parameters and the inclusion of exogenous variables. We have provided conditions for stationarity and ergodicity, as well as asymptotic results for the MLE. Moreover, we have shown the relevance of our modeling approach by focusing on a Monte Carlo simulation study. We have considered two empirical illustrations where we apply our nonlinear time-varying Poisson autoregressive model: the first applies the TV-PARX model with deterministic components to the number of COVID-19 infections in Italy; the second focuses on count time series of corporate defaults in the US by including (time-varying) exogenous covariates in the TV-PARX model specification. For these relevant illustrations, we find that the proposed time-varying framework is capable of improving the in-sample fit. Further issues are left to future research. In particular, in this paper we have specified an univariate TV-PARX model. A multivariate extension could be extremely useful to capture the relationships between different time series of counts (e.g. COVID-19 cases or number of defaults in different countries) and is currently under investigation by the authors. Moreover, the time-varying specification of $\alpha_t$ could be used to construct a new test for the presence of structural breaks for time series of counts.