The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
54,051 characters
Time-Varying Poisson Autoregression
\begin{frontmatter}
\title{\LARGE Time-Varying Poisson Autoregression}
\author[bo]{Giovanni Angelini}
\author[bo,ex]{Giuseppe Cavaliere}
\author[vu]{Enzo D'Innocenzo}
\author[bo]{Luca De Angelis}
\address[bo]{Department of Economics, University of Bologna}
\address[vu]{School of Business and Economics, Vrije Universiteit Amsterdam}
\address[ex]{Essex Business School, University of Essex}
\begin{abstract}
In this paper we propose a new time-varying econometric model, called Time-Varying Poisson AutoRegressive with eXogenous covariates (TV-PARX), suited to model and forecast time series of counts. {We show that the score-driven framework is particularly suitable to recover the evolution of time-varying parameters and provides the required flexibility to model and forecast time series of counts characterized by convoluted nonlinear dynamics and structural breaks.} We study the asymptotic properties of the TV-PARX model and prove that, under mild conditions, maximum likelihood estimation (MLE) yields strongly consistent and asymptotically normal parameter estimates. Finite-sample performance and forecasting accuracy are evaluated through Monte Carlo simulations.
The empirical usefulness of the time-varying specification of the proposed TV-PARX model is shown by analyzing the number of new daily COVID-19 infections in Italy and the number of corporate defaults in the US.
\end{abstract}
\end{frontmatter}
\section{Introduction}
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 \citep{fokianos2011log}, corporate defaults \citep{Agosto2016}, environmental economics \citep{Wang2014},
epidemiology \citep{davis2003}, COVID-19 infections and deaths \citep{khismatullina2020nonparametric,li2021will}, sports \citep{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 \cite{davis2016} for an overview of recent developments on the econometric models for discrete values.
In particular, \cite{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. \cite{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. \cite{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.
\cite{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 \cite{creal2013} and \cite{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.
\cite{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.
\cite{gorgi2020beta} uses the beta–negative binomial distribution for the narcotics trafficking reports in Sydney and the euro–pound sterling exchange rate, \cite{koopman2019} use the bivariate Poisson distribution for the number of goals in football matches and the Skellam distribution for the score difference, and \cite{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 \cite{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 \cite{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 \cite{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{sec:model} describes the new
TV-PARX model and discusses the properties of TV-PARX. Section \ref{sec:ml_est} presents the MLE for the TV-PARX. Section \ref{sec:monte_carlo} shows the finite sample properties and Section \ref{sec:empirical} reports the empirical illustrations on COVID-19 data and on corporate defaults. Section \ref{sec:conclusion} concludes the paper. Details on notation and all the proofs are provided in
\ref{proofs}.
\section{Time-Varying Poisson Autoregression}
\label{sec:model}
Consider a time series of counts $\{ y_t \}_{t \in \mathbb{Z}}$ with a Poisson conditional distribution,
\begin{align}
\label{cond_distribution}
y_t | \mathcal{F}_{t-1} \sim \mathcal{P}\textit{ois}(\lambda_t) \hspace*{1cm}t \in \mathbb{Z},
\end{align}
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 \cite{fokianos2009poisson} and \cite{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.\ \cite{creal2013} and \cite{Harvey2013} for a detailed review $-$ thus allowing for time variation in the parameters of the dynamic equation for $\lambda_t$.
\cite{Koopman2008} show that a score-driven model for Poisson time series of counts encompasses most of the observation-driven models considered by \cite{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 \eqref{cond_distribution}, which implies that the conditional $\log$-density of $y_t$ is given up to a constant by
\begin{align}
\label{cond_dens}
\log p(y_t | \mathcal{F}_{t-1})
=
{y_t} \log \lambda_t -\lambda_t.
\end{align}
We propose the following nonlinear Poisson autoregressive model with exogenous covariates and deterministic components defined as
\begin{eqnarray}
\log \lambda_{t+1} &=& \omega + \beta \log \lambda_t + \alpha_{t+1} (y_t - \lambda_t)\lambda_t^{-1} + \gamma_{t+1}^{\prime}x_t + \psi^\prime d_t \label{eq:tv_lambda_parx}
\end{eqnarray}
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 \eqref{eq:tv_lambda_parx} 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 \eqref{eq:tv_lambda_parx} is that $\alpha_t$ and $\gamma_t$ are time-varying.
Similarly to \cite{Blasques2019}, we characterize $\alpha_t$ by the following updating equation
\begin{eqnarray}
\alpha_{t+1}
&=
\delta_\alpha + \phi_\alpha \alpha_t + \kappa_\alpha (y_t - \lambda_t)\lambda_t^{-1} (y_{t-1} - \lambda_{t-1})\lambda_{t-1}^{-1},\label{eq:tv_alpha_parx}
\end{eqnarray}
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 \cite{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$, \cite{Harvey2013} shows that the law motion for $\log \lambda_t$ is given by the first-order autoregressive process in \eqref{eq:tv_lambda_parx} 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
\begin{eqnarray}
\gamma_{t+1} &=& \delta_{\gamma} +\phi_{\gamma}\gamma_{t} +\kappa_{\gamma} (y_t - \lambda_t)\lambda_t^{-1} x_t \label{eq:tv_gamma_parx}
\end{eqnarray}
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 \eqref{eq:tv_lambda_parx} 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 \cite{Agosto2016} (see Section \ref{sec:defaults} for an empirical illustration).
\medskip
\noindent {\bf Remark 1.}
Several models discussed in the literature are special cases of the model proposed here.
For instance, the model proposed by \cite{Harvey2013} is a special case of the model in \eqref{eq:tv_lambda_parx} 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 \cite{davis2003} for the case of time-invariant coefficients of the exogenous covariates. \hfill $\square$
\medskip
Despite the time-varying specification of the model in \eqref{eq:tv_lambda_parx}-\eqref{eq:tv_gamma_parx} is in the same vein of \cite{Blasques2019}, i.e.\ the coefficient $\alpha$ is allowed to change at each time period by means of the score-driven approach, with respect to \cite{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 \cite{fokianos2009poisson} and static PARX model by \cite{Agosto2016} by allowing for time-varying parameters, including the coefficients for the exogenous covariates.
\medskip
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 \eqref{cond_dens} can be expressed as
\begin{align*}
y_t = N_t(\lambda_t),
\end{align*}
where $\lambda_t$ follows the dynamics given in \eqref{eq:tv_lambda_parx} and \eqref{eq:tv_gamma_parx}, see e.g.\ \cite{fokianos2009poisson}. We begin by introducing the first result.
\begin{prop}
\label{prop_ergodic_lnLambda}
Consider the model in \eqref{eq:tv_lambda_parx} and \eqref{eq:tv_gamma_parx} under the condition $|\phi_\alpha|<1$. Then, $\{ \alpha_t \}_{t\in\mathbb{Z}}$ admits a strictly stationary and ergodic solution, with $\bar{\alpha} := \mathbb{E}[\alpha_t]= \delta_{\alpha} / (1 - \phi_\alpha) < \infty$. Moreover, under the additional conditions that $0<\beta<1$ and $ \beta| \beta + \bar{\alpha}| <1$, then, $\{ \log \lambda_t \}_{t\in\mathbb{Z}}$ is geometrically ergodic and, moreover, $\{ \lambda_t \}_{t\in\mathbb{Z}}$ is geometrically ergodic with $\mathbb{E}[\lambda_t]<\infty$.
\end{prop}
\noindent {\bf Remark 2.}
In line with the $\log$-linear Poisson autoregressive model of \cite{fokianos2009poisson}, by considering the full model given by equations \eqref{eq:tv_lambda_parx} and \eqref{eq:tv_gamma_parx} 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{prop_ergodic_lnLambda}. 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{prop_ergodic_lnLambda} are remarkably easy to check. Moreover, a direct consequence of Proposition \ref{prop_ergodic_lnLambda} 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 \cite{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.
\begin{prop}
\label{prop_sec_mom}
Under the assumptions of proposition \ref{prop_ergodic_lnLambda} $\{ \log \lambda_t \}_{t\in\mathbb{Z}}$ and $\{ y_t \}$ have finite moments of any order, i.e.\ for any $k>0$, $\mathbb{E}[|\log \lambda_t|^k]< \infty$ and $\mathbb{E}[y_t^k]< \infty$.
\end{prop}
It is worth to note that Proposition \ref{prop_sec_mom} 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.
\section{Maximum likelihood estimation}
\label{sec:ml_est}
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.
\begin{eqnarray}
\log \lambda_{t+1} &=& \omega + \beta \log \lambda_t + \alpha_{t+1} (y_t - \lambda_t)\lambda_t^{-1}
\label{ln_tv_lambda}
\end{eqnarray}
which corresponds to \eqref{eq:tv_lambda_parx} 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 \cite{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 \eqref{cond_distribution}, \eqref{ln_tv_lambda} and \eqref{eq:tv_alpha_parx}, evaluated at $\boldsymbol{\theta}_0$. The $\log$-likelihood function is given by
\begin{align}
\label{approx_log_lik_tot}
\hat{L}_T(\boldsymbol{\theta})
=
\frac{1}{T}
\sum_{t=1}^T \hat{l}_t(\boldsymbol{\theta}),
\end{align}
where
\begin{align}
\label{approx_log_lik}
\hat{l}_t(\boldsymbol{\theta})
=
y_t \log\hat{\lambda}_t(\boldsymbol{\theta}) - \hat{\lambda}_t(\boldsymbol{\theta})
\end{align}
and
\begin{eqnarray}
\label{ln_lambda_theta}
\log\hat{\lambda}_{t+1}(\boldsymbol{\theta})
&=&
\omega + \beta \log\hat{\lambda}(\boldsymbol{\theta})
+ \hat{\alpha}_{t+1}(\boldsymbol{\theta})
(y_t - \hat{\lambda}_t(\boldsymbol{\theta}) )\hat{\lambda}^{-1}_t(\boldsymbol{\theta})\\
\label{alpha_theta}
\hat{\alpha}_{t+1}(\boldsymbol{\theta})
&=&
\delta_\alpha + \phi_\alpha \hat{\alpha}_t(\boldsymbol{\theta})
+ \kappa_\alpha (y_t - \hat{\lambda}_t(\boldsymbol{\theta}) )\hat{\lambda}_t^{-1}(\boldsymbol{\theta}) (y_{t-1} - \hat{\lambda}_{t-1}(\boldsymbol{\theta}) )\hat{\lambda}_{t-1}^{-1}(\boldsymbol{\theta}).
\end{eqnarray}
Note that the time-varying parameters in \eqref{ln_lambda_theta} and \eqref{alpha_theta} 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
\begin{align}
\label{MLE}
\hat{\boldsymbol{\theta}}_T
=
\operatorname*{argmax}_{\theta \in \boldsymbol{\Theta}}
\hat{L}_T(\boldsymbol{\theta}).
\end{align}
The conditions stated in Proposition \ref{prop_ergodic_lnLambda} 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 \eqref{ln_lambda_theta} and \eqref{alpha_theta}, and therefore the approximate $\log$-likelihood function in \eqref{approx_log_lik_tot}, 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 \cite{Straumann_Mikosch2006}. Stationarity and ergodicity of $\{\log \lambda_t \}$ (as established in Proposition \ref{prop_ergodic_lnLambda}), and Proposition \ref{prop_invertibility} below, are the key ingredients for establish asymptotic properties of the MLE.
\begin{prop}
\label{prop_invertibility}
Consider the model in \eqref{ln_tv_lambda} and \eqref{alpha_theta} under the conditions of Proposition \ref{prop_ergodic_lnLambda}. Let $\boldsymbol{\Theta}$ be compact
with $\kappa_\alpha>0$ and assume that
\begin{align}
\mathbb{E}\Big[\log \sup_{\boldsymbol{\theta} \in \boldsymbol{\Theta}}
\Big|
\beta \exp\Big\{
\omega + \bar{\alpha}_{t+1}(\boldsymbol{\theta})(y_t e^{-\ell} - 1) - \ell(1-\beta)
\Big\} \Big| \Big] < 0, \label{eq:contraction}
\end{align}
where
\begin{align*}
\bar{\alpha}_{t+1}(\boldsymbol{\theta})
&:=
\delta_\alpha + \phi_\alpha \bar{\alpha}_{t}(\boldsymbol{\theta})
+ \kappa_\alpha(y_t e^{-\ell} - 1)(y_{t-1} e^{-\ell} - 1),\\
\ell
&:=
\frac{\omega - \frac{\delta_\alpha + \kappa_\alpha}{1 - \phi_\alpha}}{1 - \beta}
\end{align*}
Then,
\begin{align}
\label{contr_cond}
\sup_{\boldsymbol{\theta} \in \boldsymbol{\Theta}}|
\hat{\alpha}_t(\boldsymbol{\theta}) - {\alpha}_t(\boldsymbol{\theta})
| \xrightarrow[]{e.a.s.} 0
\hspace*{1cm}\sup_{\boldsymbol{\theta} \in \boldsymbol{\Theta}}|
\hat{\lambda}_t(\boldsymbol{\theta}) - {\lambda}_t(\boldsymbol{\theta})
| \xrightarrow[]{e.a.s.} 0
\,\,\,\,\,\text{as} \,\,\,\,\, t \rightarrow \infty,
\end{align}
where $\{ {\alpha}_t(\boldsymbol{\theta}) \}_{t\in\mathbb{Z}}$ and $\{ {\lambda}_t(\boldsymbol{\theta}) \}_{t\in\mathbb{Z}}$ are stationary and ergodic.
Moreover, if
\begin{align}
\label{contr_cond_mom}
\sup_{(\boldsymbol{\theta}\times y) \in
(\boldsymbol{\Theta} \times \mathbb{N}_0)}
\Big|
\beta \exp\Big\{
\omega + \bar{\alpha}_{t+1}(\boldsymbol{\theta})(y e^{-\ell} - 1) - \ell(1-\beta)
\Big\} \Big| <1,
\end{align}
then $\exists m >0$ such that $\mathbb{E}[\sup_{\boldsymbol{\theta} \in \boldsymbol{\Theta}} |\lambda_t(\boldsymbol{\theta})|^m]<\infty$.
\end{prop}
Condition \eqref{eq:contraction} is a sufficient condition which ensures the invertibility of the TV-PARX model by an application of Proposition 5.2.12 in \cite{Straumann2005}, or equivalently Theorem 2.8 in \cite{Straumann_Mikosch2006}. It is different from the conditions for stationarity and ergodicity given in Theorem \ref{prop_ergodic_lnLambda}, and moreover, it is only concerned with the filtering recursion relative to $\lambda_t(\boldsymbol{\theta})$ because the stationarity condition imposed in Proposition \ref{prop_ergodic_lnLambda} for $\alpha_t(\boldsymbol{\theta})$ is also sufficient for filter invertibility. The complex form of the contraction condition in \eqref{eq:contraction} is due to the exponential function. A similar problem with related discussion can be found in \cite{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 \eqref{eq:contraction} as in \eqref{contr_cond_mom}, Proposition \ref{prop_invertibility} 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 \cite{white_1994}.
First, for consistency we consider the almost sure uniform convergence of the approximate $\log$-likelihood function in \eqref{approx_log_lik_tot}. In particular, Lemma \ref{lemma_as_approxLogLik} below ensures that the average $\log$-likelihood $\hat{L}_T(\boldsymbol{\theta})$ converges uniformly to a function $L(\boldsymbol{\theta})$.
\begin{lemma}
\label{lemma_as_approxLogLik}
Let the conditions of Propositions \ref{prop_ergodic_lnLambda}-\ref{prop_invertibility} hold true. Then,
\begin{align*}
\sup_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}
| \hat{L}_T(\boldsymbol{\theta}) - {L}(\boldsymbol{\theta}) |
\xrightarrow[]{a.s.} 0
\,\,\,\,\,\text{as} \,\,\,\,\, T \rightarrow \infty,
\end{align*}
where ${L}(\boldsymbol{\theta}) := \mathbb{E}[l_1(\boldsymbol{\theta})] = \mathbb{E}[y_1\log\lambda_1(\boldsymbol{\theta}) - \lambda_1(\boldsymbol{\theta})]$.
\end{lemma}
Second, we verify the identifiability of the model parameterization, which is the argument of the next Lemma.
\begin{lemma}
\label{lemma_identifiability}
The true parameter vector $\boldsymbol{\theta}_0$ is the unique maximizer of ${L}(\boldsymbol{\theta})$ in $\boldsymbol{\Theta}$; that is, for any $\boldsymbol{\theta} \neq \boldsymbol{\theta}_0$, then ${L}(\boldsymbol{\theta}) < {L}(\boldsymbol{\theta}_0)$.
\end{lemma}
As a consequence, given the results obtained in Lemma \ref{lemma_as_approxLogLik} and Lemma \ref{lemma_identifiability}, an application of Theorem 2.11 of \cite{white_1994} implies the strong consistency of the MLE.
\begin{theorem}
\label{thm_consistency}
Let the conditions of Propositions \ref{prop_ergodic_lnLambda}-\ref{prop_invertibility} hold true. Then,
$\hat{\boldsymbol{\theta}}_T
\xrightarrow[]{a.s.}
\boldsymbol{\theta}_0$ as
$T \rightarrow \infty.$
\end{theorem}
Next, we discuss the asymptotic distribution of the MLE.
Consider the normalized score evaluated at $\boldsymbol{\theta}_0$
\begin{align}
\label{score_fun}
\boldsymbol{\bar{\eta}}_{T}
=
\sqrt{T}
\frac{\partial L_T(\boldsymbol{\theta}_0)}{\partial \boldsymbol{\theta}}
=
\frac{1}{\sqrt{T}} \sum_{t=1}^{T} \boldsymbol{\eta}_t
\ \ , \ \
\boldsymbol{\eta}_t = e_t \lambda_t \log\lambda^{\boldsymbol{\theta}_0}_t
,
\end{align}
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{lemma_deriv_proc} in the Appendix.
Asymptotic normality then follows by a standard central limit theorem for martingale difference sequences, see
e.g.\ Corollary 3.1 of \cite{Hall1980}, by showing that
\begin{align}
\label{condition_hall1}
\sum_{t=1}^T
\mathbb{E}[\boldsymbol{\eta}_t \boldsymbol{\eta}_t^\prime | \mathcal{F}_{t-1}]
\xrightarrow[]{P} \boldsymbol{V},
\end{align}
where
\begin{align}
\label{V}
\boldsymbol{V}
=
\lim_{T\rightarrow\infty}
\frac{1}{T} \sum_{t=1}^T \lambda_t (\log\lambda^{\boldsymbol{\theta}_0}_t
\log\lambda^{\boldsymbol{\theta}_0^{\prime}}_t),
\end{align}
and, the Lindberg condition holds, i.e., $\forall \epsilon > 0$,
\begin{align}
\label{condition_hall2}
\mathbb{E}[\boldsymbol{\eta}_t \boldsymbol{\eta}_t^\prime
\mathbbm{1}(\| \boldsymbol{\eta}_t \| > \epsilon)
| \mathcal{F}_{t-1}]
\xrightarrow[]{P} \boldsymbol{0} .
\end{align}
These conditions are verified in the following proposition.
\begin{prop}
\label{prop_norm_score}
Let the conditions of Propositions \ref{prop_ergodic_lnLambda}-\ref{prop_invertibility} hold true. Moreover, assume that $\mathbb{E}[|A_t|^2]<1$ and $0<|\phi_\alpha|<1$, then \eqref{score_fun}, \eqref{condition_hall1} and \eqref{V} hold and the score function in \eqref{score_fun}
satisfies
\begin{align*}
\boldsymbol{\bar{\eta}}_{T}
\overset{d}\to
\mathcal{N}(\boldsymbol{0}, \boldsymbol{V}).
\end{align*}
\end{prop}
Denote the information matrix at the true value by $\boldsymbol{J}$, i.e.
\begin{align}
\label{J_matrix}
\boldsymbol{J}
=
\mathbb{E}\bigg[
\frac{1}{\lambda_t}
\frac{\partial \lambda_t(\boldsymbol{\theta}_0)}{\partial \boldsymbol{\theta}}
\frac{\partial \lambda_t(\boldsymbol{\theta}_0)}{\partial \boldsymbol{\theta}^\prime}
\bigg].
\end{align}
Then, we finally state the following theorem.
\begin{theorem}
\label{thm_asy_norm}
Let the conditions of Propositions \ref{prop_ergodic_lnLambda}-\ref{prop_norm_score} hold true. Moreover, assume that the true parameter vector $\boldsymbol{\theta}_0$ is in the interior of $\boldsymbol{\Theta}$. Then, the MLE $\hat{\boldsymbol{\theta}}_T$ is asymptotically normal,
\begin{align*}
\sqrt{T}(\hat{\boldsymbol{\theta}}_T - {\boldsymbol{\theta}}_0)
\Rightarrow
\mathcal{N}(\boldsymbol{0}, \boldsymbol{J}^{-1}).
\end{align*}
\end{theorem}
\noindent {\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 \cite{white1982regularity} and \cite{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
\begin{align*}
\sqrt{T}(\hat{\boldsymbol{\theta}}_T - {\boldsymbol{\theta}}_0)
\Rightarrow
\mathcal{N}(\boldsymbol{0}, \boldsymbol{\Sigma}),
\end{align*}
where $\boldsymbol{\Sigma} = \boldsymbol{J}^{-1}\boldsymbol{I}\boldsymbol{J}^{-1}$, with $\boldsymbol{J}$ as in equation \eqref{J_matrix} 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 \cite{Ahmad2016} or \cite{Aknouche2021}.
\section{Finite sample performance}
\label{sec:monte_carlo}
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
\begin{equation}
y_t = N_t(\lambda_t^0), \label{eq:MC1}
\end{equation}
\noindent where $N_t(\cdot)$ denotes a sequence of independent Poisson processes with unity intensity. Similarly to \cite{Blasques2019}, for the time-varying $\lambda_t^0$, we consider a deterministic step-function
\begin{align}
\lambda_t^0
=
\begin{cases}
2 & \hspace*{1cm}\text{if} \hspace*{1cm}\sin(\gamma^{-1} 10^{-2} (\pi t - 1) \geq 0 \\
2+\delta & \hspace*{1cm}\text{if} \hspace*{1cm}\sin(\gamma^{-1} 10^{-2} (\pi t - 1) < 0
\end{cases}
\label{step}
\end{align}
\noindent 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{fig:lambda_0} 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 \eqref{ln_tv_lambda} and the static PAR, i.e.\ where $\alpha_t=\alpha$ in \eqref{ln_tv_lambda}.
Table \ref{tab:monte_carlo_RMSE} reports the Root Mean Squared Errors (RMSE) for all the cases considered. The results in Table \ref{tab:monte_carlo_RMSE} 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{fig:mc_est} 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{fig:mc_est}, 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 \eqref{step}.
\section{Empirical illustrations}
\label{sec:empirical}
In this section, we show the usefulness of the time-varying specification of the TV-PARX model by considering two empirical illustrations.
Section \ref{sec:covid} deals with the daily counts of infections of SARS-COV-2 virus in Italy until the end of May 2021.
Section \ref{sec:defaults} analyzes the monthly corporate default counts in the US as previously analyzed by \cite{Agosto2016}.
\subsection{COVID-19}
\label{sec:covid}
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 \cite{batista2020estimation} for a logistic growth model and \cite{giordano2020modelling} and \cite{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{fig:emp_fitted_actual}), 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,
\cite{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 \cite{khismatullina2020nonparametric} and \cite{dong2020time} using non-parametric approaches and by \cite{palmer2021count} using Bayesian hierarchical models.
Therefore, we fit the TV-PARX model in \eqref{eq:tv_lambda_parx}-\eqref{eq:tv_gamma_parx} 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{fig:emp_fitted_actual}).
The counts dynamics in Figure \ref{fig:emp_fitted_actual} shows a seasonal component mainly due to reporting.
Hence, we consider the following TV-PARX model
\begin{eqnarray}
\log \lambda_{t+1} &=& \omega + \beta \log \lambda_t + \alpha_{t+1} (y_t - \lambda_t)\lambda_t^{-1} + \psi^\prime d_t
\nonumber
\\
\alpha_{t+1} &=& \delta_\alpha + \phi_\alpha \alpha_t + \kappa_\alpha (y_t - \lambda_t)\lambda_t^{-1} (y_{t-1} - \lambda_{t-1})\lambda_{t-1}^{-1}
\label{covid}
\end{eqnarray}
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 \eqref{covid}.
Table \ref{tab:estimates} shows the estimates and Table \ref{tab:empirical} 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{fig:emp_fitted_actual} 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.
\subsection{Corporate defaults}
\label{sec:defaults}
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.
\cite{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 \cite{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 \cite{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{fig:emp_fitted_actual_defaults} shows the default count time series. As discussed in \cite{Agosto2016}, Figure \ref{fig:emp_fitted_actual_defaults} 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{fig:emp_fitted_actual_defaults}) 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 \eqref{eq:tv_lambda_parx}-\eqref{eq:tv_gamma_parx}, 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 \cite{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 \cite{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
\begin{eqnarray}
\log \lambda_{t+1} &=& \omega + \beta \log \lambda_t + \alpha_{t+1} (y_t - \lambda_t)\lambda_t^{-1} + \gamma_{1,t+1}RV_t + \gamma_{2,t+1}LI_t
\nonumber
\\
\alpha_{t+1} &=& \delta_\alpha + \phi_\alpha \alpha_t + \kappa_\alpha (y_t - \lambda_t)\lambda_t^{-1} (y_{t-1} - \lambda_{t-1})\lambda_{t-1}^{-1} \nonumber \\
\gamma_{1,t+1} &=& \delta_{\gamma_1} +\phi_{\gamma_1}\gamma_{1,t} +\kappa_{\gamma_1} (y_t - \lambda_t)\lambda_t^{-1} RV_t \nonumber\\
\gamma_{2,t+1} &=& \delta_{\gamma_2} +\phi_{\gamma_2}\gamma_{2,t} +\kappa_{\gamma_2} (y_t - \lambda_t)\lambda_t^{-1} LI_t. \label{eq:emp_default}
\end{eqnarray}
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 \eqref{eq:emp_default}. The maximum likelihood estimates are reported in Table \ref{tab:estimates}.
According to the information criteria and the RMSE reported in Table \ref{tab:empirical}, 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 \cite{Agosto2016}, namely (non-logs) PARX(2,1).}
As can be noted from Figure \ref{fig:emp_fitted_actual_defaults}, 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{fig:emp_alphagamma_defaults} 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{fig:emp_alphagamma_defaults} 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 \eqref{eq:tv_lambda_parx}-\eqref{eq:tv_gamma_parx} provides a more flexible alternative to the time-invariant counterpart.
\section{Conclusion}
\label{sec:conclusion}
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.
\newpage
\bibliographystyle{apalike}
\bibliography{references}
\newpage
\begin{table}[h]
\centering
\begin{tabular}{c c c c c c c}
\hline \hline
& \multicolumn{2}{c}{$\gamma=1$} & \multicolumn{2}{c}{$\gamma=1.5$} & \multicolumn{2}{c}{$\gamma=2$} \\
& PARX & TV-PARX & PARX & TV-PARX & PARX & TV-PARX \\
\hline
\multicolumn{7}{c}{$T=250$} \\
$\delta=0$ & \textbf{0.0077} & 0.0086 & \textbf{0.0077} & 0.0086 & \textbf{0.0077} & 0.0086 \\
$\delta=2$ & 0.4052 & \textbf{0.3686} & 0.3317 & \textbf{0.2630} & 0.3358 & \textbf{0.2776} \\
$\delta=4$ & 0.6331 & \textbf{0.5928} & 0.5034 & \textbf{0.3749} & 0.5002 & \textbf{0.3842} \\
\hline
\multicolumn{7}{c}{$T=500$} \\
$\delta=0$ & \textbf{0.0033} & 0.0049 & \textbf{0.0033} & 0.0049 & \textbf{0.0033} & 0.0049 \\
$\delta=2$ & 0.4097 & \textbf{0.3914} & 0.3729 & \textbf{0.3477} & 0.3371 & \textbf{0.3115} \\
$\delta=4$ & 0.6324 & \textbf{0.6092} & 0.5726 & \textbf{0.5359} & 0.5198 & \textbf{0.4810} \\
\hline
\multicolumn{7}{c}{$T=1000$} \\
$\delta=0$ & \textbf{0.0026} & 0.0031 & \textbf{0.0026} & 0.0031 & \textbf{0.0026} & 0.0031 \\
$\delta=2$ & 0.4245 & \textbf{0.4121} & 0.3784 & \textbf{0.3644} & 0.3373 & \textbf{0.3212} \\
$\delta=4$ & 0.6529 & \textbf{0.6363} & 0.5863 & \textbf{0.5626} & 0.5224 & \textbf{0.4931} \\
\hline \hline
\end{tabular}
\caption{Root Mean Squared Error (RMSE) where the error is between the true $\lambda_t^0$ and the filtered parameter $\lambda_t$ from TV-PARX and PARX models for different true values of $\delta$ and $\gamma$.
\label{tab:monte_carlo_RMSE}}
\end{table}
\newpage
\begin{table}[h]
\centering
\begin{tabular}{c c c c c}
\hline \hline
& \multicolumn{2}{c}{\textbf{SARS-COV-2}} & \multicolumn{2}{c}{\textbf{Defaults}} \\
& TV-PARX & PARX & TV-PARX & PARX \\
\hline
$\omega $& $\underset{(0.0002)}{0.0045}$ & $\underset{(0.0002)}{0.0039}$ & $\underset{(0.0323)}{-0.0227}$ & $\underset{(0.0056)}{0.0219}$\\
$\beta $& $\underset{(0.0001}{0.9990}$ & $\underset{(0.0001)}{0.9990}$ & $\underset{(0.0386)}{0.9423}$ & $\underset{(0.0038)}{0.9441}$ \\
$\psi_1 $& $\underset{(0.0197)}{-0.0514}$ & $\underset{(0.0002)}{-0.0416}$ & - & - \\
$\psi_2 $& $\underset{(0.0071)}{-0.1310}$ & $\underset{(0.0001)}{-0.1246}$ & - & - \\
$\psi_3 $& $\underset{(0.0004)}{-0.2767}$ & $\underset{(0.0001)}{-0.2724}$ & - & - \\
$\psi_4 $& $\underset{(0.0106)}{0.1739}$ & $\underset{(0.0005)}{0.1794}$ & - & - \\
$\psi_5 $& $\underset{(0.0065)}{0.1578}$ & $\underset{(0.0002)}{0.1602}$ & - & - \\
$\psi_6 $& $\underset{(0.0100)}{0.1043}$ & $\underset{(0.0006)}{0.1083}$ & - & - \\
$\delta_{\alpha} $& $\underset{(0.0017)}{-0.0013}$ & $\underset{(0.0021)}{0.7033}$ & $\underset{(0.0061)}{-0.0012}$ & $\underset{(0.0021)}{0.2369}$ \\
$\phi_{\alpha} $& $\underset{(0.0004)}{0.9984}$ & - & $\underset{(0.0227)}{0.9946}$ & - \\
$\kappa_{\alpha} $& $\underset{(0.0250)}{0.5000}$ & - & $\underset{(0.0798)}{0.0157}$ & -\\
$\delta_{\gamma_1} $& - & - & $\underset{(2.6856)}{1.1954}$ & $\underset{(0.1994)}{11.8267}$ \\
$\phi_{\gamma_1} $& - & - & $\underset{(0.1418)}{0.9498}$ & - \\
$\kappa_{\gamma_1} $& - & - & $\underset{(273.6449)}{917.6204}$ & - \\
$\delta_{\gamma_2} $& - & - & $\underset{(0.0062)}{-0.0004}$ & $\underset{(0.0018)}{0.0010}$ \\
$\phi_{\gamma_2} $& - & - & $\underset{(0.3673)}{0.9766}$ & - \\
$\kappa_{\gamma_2} $& - & - & $\underset{(0.0021)}{0.0070}$ & - \\
\hline \hline
\end{tabular}
\caption{Estimated parameters, standard errors in brackets.
\label{tab:estimates}}
\end{table}
\newpage
\begin{table}[h]
\centering
\begin{tabular}{c c c}
\hline \hline
& \textbf{COVID-19} & \textbf{Defaults} \\
& TV-PARX $-$ PARX & TV-PARX $-$ PARX \\
\hline
log-likelihood & 422.53 & 15.11 \\
AIC & -841.05 & -18.21 \\
HQC & -837.79 & -8.94 \\
BIC & -832.78 & 5.10 \\
RMSE & -14.29 & -0.01 \\
\hline \hline
\end{tabular}
\caption{Difference between the information criteria and the Root Mean Squared Error (RMSE) of the estimated TV-PAR and PAR models for the empirical analysis of the daily number of new COVID-19 cases in Italy and the number of monthly corporate defaults in the US.
\label{tab:empirical}}
\end{table}
\newpage
\begin{figure}[h]
\centering
\begin{tabular}{c c c}
\multicolumn{3}{c}{$T=250$} \\
\includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelA_11.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelA_12.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelA_13.pdf} \\
\includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelA_21.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelA_22.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelA_23.pdf} \\
\includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelA_31.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelA_32.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelA_33.pdf} \\
\multicolumn{3}{c}{$T=500$} \\
\includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelB_11.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelB_12.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelB_13.pdf} \\
\includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelB_21.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelB_22.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelB_23.pdf} \\
\includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelB_31.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelB_32.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelB_33.pdf} \\
\multicolumn{3}{c}{$T=1000$} \\
\includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelC_11.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelC_12.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelC_13.pdf} \\
\includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelC_21.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelC_22.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelC_23.pdf} \\
\includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelC_31.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelC_32.pdf} & \includegraphics[width=3cm]{Figure/MonteCarlo/Figure1_PanelC_33.pdf} \\
\end{tabular}
\caption{Dynamic of $\lambda_t^0$ for $\delta \in \{0,2,4 \}$ and $\gamma \in \{1, 1.5, 2\}$. \label{fig:lambda_0}}
\end{figure}
\newpage
\begin{figure}[h]
\centering
\begin{tabular}{c}
\multicolumn{1}{c}{$T=250$} \\
\includegraphics[width=10cm]{Figure/MonteCarlo/Figure2_PanelA.pdf} \\
\multicolumn{1}{c}{$T=500$} \\
\includegraphics[width=10cm]{Figure/MonteCarlo/Figure2_PanelB.pdf} \\
\multicolumn{1}{c}{$T=1000$} \\
\includegraphics[width=10cm]{Figure/MonteCarlo/Figure2_PanelC.pdf} \\
\end{tabular}
\caption{Estimation of $\lambda_t$ for $\delta = 4 $ and $\gamma = 1$ with the 95 \% confidence bands obtained with the TV-PARX (in green) and the static PARX (in red). \label{fig:mc_est}}
\end{figure}
\newpage
\begin{figure}[h]
\includegraphics[width=16cm]{Figure/Empirical/emp_fitted_actual.pdf}
\caption{Fitted vs actual daily number of cases of Covid-19 in Italy. \label{fig:emp_fitted_actual}}
\end{figure}
\newpage
\begin{figure}[h]
\includegraphics[width=16cm]{Figure/Empirical/emp_alpha_95.pdf}
\caption{Estimated time-varying parameter $\alpha_t$ in \eqref{covid} for the daily number of cases of Covid-19 in Italy. Shaded areas are the 95\% confidence intervals. \label{fig:emp_alpha}}
\end{figure}
\newpage
\begin{figure}[h]
\includegraphics[width=16cm]{Figure/Empirical/emp_fitted_actual_defaults_2cov.pdf}
\caption{Fitted vs actual monthly number of corporate defaults in US. \label{fig:emp_fitted_actual_defaults}}
\end{figure}
\newpage
\begin{figure}[h]
\begin{center}
\includegraphics[width=10cm]{Figure/Empirical/emp_alpha_defaults_2cov.pdf} \\
\includegraphics[width=10cm]{Figure/Empirical/emp_gamma1_defaults_2cov.pdf} \\
\includegraphics[width=10cm]{Figure/Empirical/emp_gamma2_defaults_2cov.pdf} \\
\caption{Estimated time-varying parameter $\alpha_t$, $\gamma_{1,t}$ and $\gamma_{2,t}$ in \eqref{eq:emp_default} for the monthly number of corporate defaults in US. Shaded areas are the 95\% confidence intervals. \label{fig:emp_alphagamma_defaults}}
\end{center}
\end{figure}
\clearpage