{'fraction_non_alphanumeric': 0.03340276557833714, 'fraction_numerical': 0.015659899114111475, 'mean_word_length': 4.848818435344195, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian time‐varying autoregressive models of COVID‐19 epidemics

Abstract The COVID‐19 pandemic has highlighted the importance of reliable statistical models which, based on the available data, can provide accurate forecasts and impact analysis of alternative policy measures. Here we propose Bayesian time‐dependent Poisson autoregressive models that include time‐varying coefficients to estimate the effect of policy covariates on disease counts. The model is applied to the observed series of new positive cases in Italy and in the United States. The results suggest that our proposed models are capable of capturing nonlinear growth of disease counts. We also find that policy measures and, in particular, closure policies and the distribution of vaccines, lead to a significant reduction in disease counts in both countries.


INTRODUCTION
Time series of observed counts arise in a variety of contexts in which a number of events are recorded over time, including studies of incidences of a particular disease (Davis, 2003;Zeger, 1988). The reported occurrences are often expressed with daily, weekly, or monthly frequency. Disease modeling has significant practical importance for governments since it provides information that can improve decisionmaking. In addition, studying the correlation between counts and policy variables allows to better gauge the impact of control measures. This paper provides a Bayesian analysis of a time series of counts with the goal of providing reliable short-term predictions of disease counts and to assess their correlation with explanatory variables. The time series of interest is the incidence of COVID-19 infectious disease in Italy and Unites States, two paradigmatic countries, and the explanatory variables are the varying Non-pharmaceutical interventions (NPIs) implemented by most countries to stop the virus propagation.
Since COVID-19 first emerged in early 2020, it has infected nearly every country and claimed the lives of several thousands of people worldwide. So far, many countries have been locked-down and stringent social distancing measures have been implemented to contrast the high incidence of the infection over time and its high economic consequences. Previous studies, mostly focused on the first wave of the pandemic (before July 2020), have shown the effectiveness of NPIs in pandemic mitigation (Cowling et al., 2020;Hao et al., 2020). They have also shown that the impact of NPIs may differ among countries with various country characteristics, like health capacity, population density, air temperature, etc. (Esteve et al., 2020;Marcon et al., 2020). Existing studies also indicate that vaccination is the most promising strategy for overcoming the epidemic. However, the unequal distribution and allocation of vaccines among nations and demographic groups may obstruct the path to herd immunity. It is crucial to understand how vaccination combined with NPIs reduces COVID-19 transmission (Giordano et al., 2021).
The Poisson AutoRegressive (PAR) model is one of the most commonly used models for analyzing disease time series data, as it can well capture time dependence among observations. It can be seen as a special case of the generalized linear model (GLM) methodology (Fahrmeir et al., 2013;Kedem & Fokianos, 2002;Nelder & Wedderburn, 1972) where the link function is a linear function of the covariates, with coefficients that do not change over time. The advantage of this specification is the ease of model implementation and interpretation. However, when the impact of the epidemic changes over time, especially under the effect of policy measures, this type of model may be restrictive, and time-varying coefficients may be needed.
A flexible way to introduce time-varying parameters in a GLM model is through spline basis functions, as in the seminal work of Cleveland and Loader (1996). Spline basis functions have also been proposed to specify time-varying coefficients in AutoRegressive (AR) time series models (Amorim et al., 2008;Cai et al., 1999;Huang & Shen, 2004). For this type of models, Zeger (1988) has suggested a generalized estimating equation methodology, which still is the most commonly used estimation method. Within the Bayesian framework, Chan and Ledolter (1995) proposed a Monte Carlo expectation maximization (MCEM) technique. More recently, Hastie et al. (2001) compared different estimation procedures. Moving to the more complex autoregressive moving average (ARMA) models, the only work that considers time-varying parameters is that of Andrade et al. (2015).
None of these papers focuses on the time-varying character of the coefficients. We contribute to fill the gap with a model that extends Fokianos et al. (2008) within a time-varying coefficient framework. In doing so, we extend the timedependent contagion models recently proposed by Agosto et al. (2021), Geir et al. (2022), and Roy and Karmakar (2021) with the inclusion of time-varying policy covariates. And we also extend the causal contagion models recently proposed by Bo et al. (2021) and Islam et al. (2020) with the inclusion of time-varying covariate effects, which also allows to compare alternative time lags of policy measures. Note also that our proposed model allows to account for overdispersion of counts, as in the Markov-switching Poisson GARCH model proposed by Roberts et al. (1998).
More specifically, We consider two classes of time series generalized linear models. The first is a class of AR models, in which the conditional mean of the Poisson process is linked to the past observations and to the potential covariate effects. The second one adds the dependence of the past values on the conditional mean. Since the mean of a Poisson model is the same as the variance, the latter can be seen as an Integer-Valued Generalized Autoregressive Conditional Heteroskedasticity (INGARCH) model. Originally introduced by Ferland et al. (2006), these models have been thoroughly analyzed in Zhu (2012aZhu ( , 2012bZhu ( , 2012c and Roberts et al. (1998).
We thus compare alternative models along three dimensions: AR versus INGARCH structures; time-constant versus time-varying coefficient, and models with or without covariates while testing for different time lags.
The rest of the paper is organized as follows. Section 2 describes the proposed models, Section 3 describes the Bayesian specification and inference. Section 4 describes the data. Section 5 presents the main empirical findings. Section 6 summarizes the highlights of the paper and suggests future research directions.

MODEL SPECIFICATION
The time dependence of COVID-19 transmission represents a critical factor to understand the epidemic dynamics. The reference Susceptible-infected-recovery (SIR) model assumes an exponential growth model, which assumes that the logarithm of the counts follows a linear trend, summarized by 0 , a function of the ratio between the new cases originating on consecutive days. However, the evolution of the epidemics is characterized by time heterogeneity, which leads to a nonlinear growth. This derives, for example, from time-varying testing and measurement procedures, among others. From the previous consideration, it appears that the exponential growth assumption needs to be relaxed in favor of a time-dependent modeling framework. Recently, Agosto et al. (2021) proposed a PAR model for epidemic counts that can take into account both short-term and long-term dependence on past contagion shocks. The model is specified by a Poisson conditional distribution for the COVID-19 counts, along with a logarithmic link function which allows to model the effect of covariates without constraining them to assume positive values, as in a linear link model.
While the model has been quite useful to monitor contagion in the first wave of the pandemic, it has become less performing during the following times since, even if there is a nonlinear growth assumption, it is based on coefficients that do not change over time and, therefore, cannot easily adapt to regime switches.
In this paper, we aim to extend the fixed coefficients PAR model of Agosto et al. (2021) into a time-varying coefficient framework, more flexible in model specification. Our modeling strategy is motivated to adapt to the rapid changes of the COVID spread. Hence, we propose a nonstationary time-varying model, which is suitable for our analysis. The constraints of the parameter space in Equations (7) and (18) are motivated by the stationary conditions for the time-constant case. This is fairly common in the time-varying AR and INGARCH model literature (Dahlhaus & Rao, 2006;Davis & Mikosch, 2009;Fryzlewicz et al., 2008;Ferreira et al., 2017;Rohan & Ramanathan, 2013).
In detail, we consider two alternative specifications: one based on INGARCH models for counts, equivalent to the PAR models and other one based on the more parsimonious AR models.

AR models
The time series of counts reported at time (day) is assumed to follow a Poisson distribution, conditional on the information up to − 1: where  −1 is the -field that contains all relevant information up to time . Let = ( ,1 , … , , ) be a time-varyingdimensional covariate vector and let = log( ). A Poisson autoregressive model can then be specified, following Agosto et al. (2021), by the following link function: where ∈ ℝ, ∈ ℝ with = 1, … , and = ( 1 , … , ) ∈ ℝ are unknown parameters to be estimated from the data. From an interpretational viewpoint, is the intercept, describes the dependence of the expected count of cases on past observations and describes a dependence of the expected counts on explanatory covariates at time point . The above model, as discussed in Agosto et al. (2021), allows to take into account time dependency of COVID-19 counts, thereby reducing overdispersion. However, it remains the issue of a strong time heterogeneity of the policy measures taken by different countries, which can lead to an incorrect measurement of the effects of the covariates. To solve this issue, we propose to extend the AR model into a time-varying coefficient version, in which the link function for becomes: where ( ∕ ) and ( ∕ ) represent, for = 1, … , , time-varying splines. We recall that a spline is a continuous function that coincides with a polynomial on every subinterval in which it is defined. The coefficients of the polynomial differ from interval to interval, but the order of the polynomial is the same. We consider, for parsimony, cubic B-splines. The introduction of splines allows the effect of covariates to vary along different time intervals, trying to better capture the effect of time-varying policy measures.

INGARCH models
The PAR model of Agosto et al. (2021) is an INGARCH model which allows to specify the logarithmic link not only as a function of the past observations but also as a function of the past conditional mean of the process.
More formally, the INGARCH model can be specified as follows: where the additional parameter ∈ ℝ, = 1, … , represents the dependence of the logarithmic link on the conditional mean of the Poisson process.
As done for AR models, we propose to consider a time-varying extension of INGARCH models, described by the following: The constraints introduced in Equation (7) on ( )s are motivated by the time constant models where these conditions lead to stationarity and ergodicity. In the time-varying AR literature, they are not uncommon (Dahlhaus & Rao, 2006;Fryzlewicz et al., 2008).
The prior distributions for the time-varying Bayesian AutoRegressive (tvBAR) models can be specified, for each B-spline basis expansion as follows: , ∼ (0, 100), 1 ≤ ≤ , 1 ≤ ≤ 2 , , ∼ (0, 100), 1 ≤ ≤ , 1 ≤ ≤ 3, where are B-spline basis functions. The subscripts i and k refer to the time-constant dimension of the parameters (i.e., = 1, … , and = 1, … , ), while s refers to the dimension of the B-splines functions. For the unconstrained parameters and , (which respectively refer to the intercept and the effect of explanatory covariates) we assume independent Gaussian priors in Equations (11) and (12) with variance equal to 100 to make the prior weakly informative. With reference to the constrained autoregressive parameter, we may specify a uniform prior in order to satisfy the stationarity conditions of time-constant models.
However, from a computational viewpoint, this implies that the autoregressive parameters ( ) are to be sampled in a bounded space between -1 and 1, reducing the efficiency of the Markov chain Monte Carlo (MCMC) posterior approximation.
where 1 ≤ ≤ . We therefore propose a parameter transformation, mapping the bounded parameter vector Θ to an unbounded parameter vector to improve sampling efficiency. The proposed mapping function is the following: where ∼ (−1, 1). In detail, the final autoregressive parameter ( , ) to be sampled in Langevin Monte Carlo (LMC) algorithm is obtained as follows: where 1 ≤ ≤ and 1 ≤ ≤ 2 . In this way, we put normal priors on the transformed B-spline coefficients to be sampled within the LMC algorithm and then transformed back the parameters for the likelihood function ( −1 ( , ), where −1 is the inverse of function ). This method results in better mixing of the posterior samples. The log-likelihood function of the model, combined with the proposed prior distributions, gives rise to the following log posterior distribution ( 1 ) for the tvBAR model: For each B-spline basis expansion, prior distributions can be specified as follows: where are B-spline basis functions. The subscripts i, j, and k refer to the time-constant dimension of the parameters (i.e., = 1, … , , = 1, … , , and = 1, … , ), while s refers to the dimension of the B-splines functions.
Priors for and can be recovered from Equations (11) and (13). Moreover, we need to additionally estimate the parameter 0 to start the recursion of the conditional mean (Equation 27).
Similar to what is done in Section 3.1, we propose a transformation of the time-dependent parameters to ensure the condition about the absolute value of their sum (Equations 24 and 25 ), taking the logarithmic function (15) to map into an unbounded parameter space for sampling. The starting values for time-dependent parameters ( , ), together with their transformed version ( , , , ) can be summarized as follows: , , ∼ (0, 100), where 1 ≤ ≤ , 1 ≤ ≤ , and 1 ≤ ≤ with = ( 2 , 3 ).
Specifically, , , , , , , , and 0 are the parameters to be sampled in MCMC iterations. After sampling, the parameters , and , are transformed back ( −1 ( , ), −1 ( , )) for the likelihood function. The posterior distribution  2 for the Bayesian INGARCH (BINGARCH) model can then be obtained as follows for the time-varying case: To calculate the posterior distribution from the previous expressions, we employ the gradient-based. The Hamiltonian Monte Carlo (HMC) algorithm is similar to Roy and Karmakar (2021). Tuning of HMC parameters is necessary to obtain a good acceptance rate. HMC has two parameters, the number of leapfrog steps and step length. Since gradient computation is expensive, we consider setting the number of leapfrog steps to one. In this way, HMC becomes LMC (Neal, 2011). According to Roberts et al. (1998), the optimal acceptance rate for LMC is around 0.57. We thus maintain an acceptance rate between 0.45 and 0.7 by adaptively tuning the step-length parameter. Based on our results, this produces excellent mixing of the posterior samples (see Figure 3 for reference).
We want to remark that, due to the increasing complexity of the parameter space for the INGARCH model, , , , , and 0 are sampled together.

DATA
The proposed models have been applied to two types of data sources: COVID-19 impact data (taken from the European Centre for Disease Prevention and Control (ECDC): www.ecdc.europa.eu/en/covid-19/data) and policy intervention data (taken from Oxford Covid-19 Government Response Tracker (OxCGRT): www.bsg.ox.ac.uk/research/research-projects/ covid-19-government-response-tracker). From ECDC, we collected daily frequency COVID-19 count data from two representative countries: Italy and the United States. We consider, as target response variable, the change in total positive cases, defined as the difference between the total counts of positives for the current day and those the day before. The starting point of the time series was chosen as the date when the country recorded more than 100 recorded cases.
We also obtained data on the taken policy measures from OxCGRT (Hale et al., 2020), a project from Oxford's Blavatnik School of Government, in which each government response to the pandemic is measured according to an ordinal scale that ranges from no measures to the most severe form of measure. Our primary interventions of interest are the measures defined as C-containment and closure policies as well as those defined as H-health system policies.
We remark that the epidemic and the policy measures taken to tackle the epidemics have been quite heterogeneous across countries, with NPI measures taken at the local authority level, varying across regions and across time. This heterogeneity suggests that these covariates have an intrinsic measurement error.
Tables 1 and 2 show a preliminary analysis of the OxCGRT variables of interest. Table 1 indicates sparse frequency data, mostly in correspondence with health system policies, which have less frequent changes over time, compared to containment and closure policies. Table 2 indicates high collinearity between policy measures, which can be explained by the fact that NPIs have not been implemented independently of each other but, rather, in a multivariate policy mix approach.
We have addressed these issues by grouping NPIs into homogeneous categories, in line with Cowling et al. (2020) and Davies et al. (2020) who concluded that combined interventions have been the most effective strategy. Bo et al. (2021) also showed that social distancing in combination with other NPIs has led to a significant decrease in COVID-19 transmission.
More specifically, following , the original policy variables have been transformed from the ordinal to the binary scale, taking a value of "0" in case of no intervention or recommended intervention and a value of "1" in case of required intervention. We also grouped policies into five variables, following Islam et al. (2020): closure policies (X1_closures), restrictions on mass gatherings and events (X2_events), lockdown-type policies (X3_lockdown), testing, tracing, and masks (X4_masks), and vaccination policy (X5_vaccines). The precise mapping between the OxCGRT variables and ours is shown in Table 3.

EMPIRICAL FINDINGS
We first implemented the frequentist AR and INGARCH models using the tsglm function from the R package tscount (Liboschik et al., 2017). As shown in Table 4, the order of the INGARCH structure has been selected according to the TA B L E 3 Transformation of OxCGRT NPIs. Variables have been transformed from ordinal to binary scale taking the value of 0 in case of "no intervention" or "recommended intervention" and 1 in case of "required intervention." Policy variables have then been grouped based on the thematic area of reference absolute error of the median predictions with respect to the true values (median absolute percentage error (MAPE))). Based on the mean incubation period of the disease, equal to 1, 7, and 14 have been tested. Results show that equal to 1 reports the lowest MAPE in both Italy and United States. Then, the time-dependent model structure has been extended adding exogenous covariates for policy effect. Additionally, we have implemented the same model structures within the Bayesian framework for which we have developed specific R functions. Since we extend the frequentist framework models by introducing B-splines, we discuss about the selection of the number of knots for the B-splines (Table 7) before presenting parameter estimates for both time-constant and time-varying coefficients.
Finally, we compared and discussed the results obtained in terms of predictive performance, using the last week of data as a test set. Point forecasts have been evaluated in terms of MAPE, whereas the whole distribution has been evaluated using scoring rules (Gneiting & Raftery, 2007), implemented in the R package scoringutils (Bosse et al., 2020).
The source code to reproduce the results is available as Supporting Information on the journal's web page.

Frequentist analysis
The results from the frequentist Poisson AR/INGARCH models, grouped by country, are shown in Table 5. The AR estimates in Table 5 show a strong positive 1 for both countries, which indicates a positive dependence of counts over time. The INGARCH estimates in the same table show that, for Italy, 1 is much greater than 1 . This indicates that short-term dependence on past observed disease counts prevails over long-term dependence: a result that could be explained by the containment and closure measures introduced very promptly in response to the COVID-19 spread. Differently, in the United States, the estimates for 1 and 1 have a similar magnitude, revealing the persistence of long memory effects.
The impacts of policy measures are summarized in Table 6. Table 6 shows that, for Italy, closure policies (X1_closures) and vaccines (X5_vaccines) have significant negative coefficients, meaning that these policy measures have significantly reduced infection counts. In contrast, restrictions on gatherings and events (X2_events), lockdown measures (X3_lockdown), and testing, tracing, and masks (X4_masks) reveal a significant positive impact on counts. This could be explained by the fact that these policy measures are typically imposed when the number of cases is high, leading to a "reverse effect" that leads to stronger policies when counts are higher.
The results for the United States also indicate negative significant effects for closure policies (X1_closures) and vaccines (X5_vaccines), but also for restrictions on gatherings and events (X2_events) and for testing, tracing, and masks policies (X4_masks). Only lockdown-type policies (X3_lockdown) have a positive coefficient, which can be attributed to the "reverse effect" found for Italy. These results are in line with Hsiang et al. (2020) who have found effectiveness of most anticontagion policies in slowing down the exponential growth of the pandemic in the United States.

Bayesian analysis
To implement the proposed Bayesian time-varying models, a preliminary selection of the number of spline knots is necessary. A large number of knots can cause overfitting, whereas a small number can result in underfitting. We select  Table 7. From Table 7, the best number of knots, which minimizes the MAPE for the tvBAR model is 4 for Italy and 6 for the United States; for the tvBINGARCH model, it is 4 for both countries. Figures 1 and 2 show the fitted cubic B-splines of the daily new observed cases of COVID-19 for Italy and for the United States, applying tvBAR models with, respectively, four and six knots. Figure 2 shows that the peaks and the evolution of COVID-19 data are well captured by a spline with six knots for the United States, which identifies three main peaks. Figure 1 shows that a simpler four knots spline is a good approximation for Italy, which identifies two main peaks, smoothing the two observed peaks in the second wave of the pandemic.
To derive the posterior distributions, the LMC algorithm has been run for 10,000 iterations, with the first 5000 discarded as burn-in. All computations have been validated with convergence diagnostic assessments. For illustration, we report here the trace plots of the average likelihood function deviations along the MCMC chain for the coefficient functions ( ), where ( Tables 8 and 9 for both countries. From Tables 8 and 9, note that, for the United States, the frequentist findings in Section 5.1 are confirmed by the Bayesian estimates: both short-term and long-term dependence of counts are present. The same result is found for Italy, which is different from what occurred with the frequentist estimates, and consistent with the similarity between the two countries. We now examine the posterior distributions of the time-varying coefficients. Figures 4 and 5 depict the time evolution of the estimated and 1 from the Bayesian AutoRegressive (BAR) models for Italy and the United States. Figure 4 suggests that, for Italy, the estimated autoregressive coefficient of order 1 has a high impact on disease counts at the beginning of the period of analysis: the first wave of the pandemic. The impact starts to decrease when the number of cases decrease and increases in correspondence to the second wave of the pandemic. An almost opposite behavior is shown by the mean trend, which depends on the stock of cases, rather than on their flows. In April 2021, decreasing autoregressive coefficient estimates emerge, together with a low mean trend, resulting in a decreasing magnitude of the pandemic.  Figure 5 suggests that, for the United States, no evident turning points in the trend of estimated functions appear; the autoregressive coefficient is almost always close to 1 with a high mean trend. In the last period of analysis, the autoregressive component decreases and the mean trend becomes the prevailing factor. Figures 6 and 7 report the time evolution of the , 1 , and 1 estimated parameters from the BINGARCH models. For Italy, Figure 6 confirms the results from the BAR model, also showing an increasing contribution of the conditional mean in the initial period of analysis, reaching the peak in the summer, where the contagion was decreasing and quite stable.
For the United States, Figure 7 reports a higher contribution of 1 over 1 during the first wave of the pandemic and the opposite behavior in the second wave. A growing mean trend is still present in April 2021.
The further step of the analysis involves adding NPIs covariates within the proposed modeling framework, comparing the results from the BAR and BINGARCH models. Tables 10 and 11 show the results, when time constant parameters are considered.
The frequentist results in Table 6 are confirmed by the Bayesian posterior estimates in Tables 10 and 11. Table 10 suggests that for Italy the closures of schools and workplaces (X1_closures) and vaccination policy (X5_vaccines) are both effective in reducing COVID-19 counts; whereas other variables have a positive sign. Table 11 confirms that, in the United States, most policies (with the exception of lockdown measures) reduce COVID-19 counts.
When time-varying parameters are considered, the median point estimates in Tables 10 and 11 are replaced by a range of estimates, corresponding to different time points. Figures 8 and 9 show the results for Italy and the United  States, comparing the median estimates of the Bayesian AutoRegressive eXogenous (BARX) models with the estimates of the time-varying Bayesian AutoRegressive eXogenous (tvBARX) models. In addition, we consider different lag length windows of 1, 7, and 14 days to identify a possible time lag between policy interventions and their epidemiological effect. From Figure 8, we can gauge that a time lag of about 1-2 weeks seems to increase the impact of covariates representing containment and closure policies (X1_closure), confirming that a large number of days are needed between the activation of this type of policy and the realization of its epidemiological effect. Figure 9 suggests a more immediate time impact of policies for the United States, as they lead to more negative estimates without time lag or with a delay of at most 1 day.
Both Figures 8 and 9 show a large variation in policy effects within the considered period. In particular, the effects of lockdowns (X3_lockdowns) and testing, tracing, and masks (X4_masks) vary much more than any of the other policy TA B L E 1 1 Posterior median and 5%, 25%, 75%, and 95% posterior quantiles of NPI covariates using time-constant BARX/BINGARCHX models for the United States  measures. This can be attributed to the presence of reverse causality effects, as policies have been implemented in response to the rise in pandemic counts. Overall, it can be concluded that Bayesian estimates are in line with frequentist ones since the model structure is the same in the time constant case, but, additionally, incorporate greater uncertainty by means of the posterior distribution and the time-varying coefficient estimates.

Model comparison
All models have been assessed in terms of the absolute error of the median forecasts, calculated as the absolute value of the difference between the true values and median posterior forecasts obtained from MCMC iterations. A high MAPE means that the median forecasts tend to be far away from the true values. The accuracy of model prediction intervals has also been evaluated in terms of overprediction, underprediction, and bias, calculated from the scoring rules applied to the predictive Monte Carlo samples. They are defined with respect to the upper bound (overprediction) or lower bound (underprediction) of the prediction interval. The bias is a measure between -1 and 1 that expresses the tendency to underpredict (-1) or overpredict (1) the true values. All posterior predictions The results show that most of the proposed models fit the data quite well. For Italy, the tvBARX model seems to be the best, with a small MAPE (equal to 12%) on the training set and a bias almost equal to 0. The model also performs quite well on test data, with a 15% MAPE and a bias equal to 0.057. Similar conclusions are obtained for the United States, although the bias on the test set is higher (0.436), consistent with the higher number of cases with respect to Italy.
To further demonstrate the accuracy of our proposed models, Figures 10 and 11 plot the posterior median of the fitted and predicted values, together with the number of observed cases of COVID-19 for Italy and the United States. The interval range (50-90%) has been added to both figures.
From Figures 10 and 11 both models appear to perform quite well. We finally remark that our proposed TvBINGARCHX model was employed within the ECDC European Forecast Hub initiative and compared with many other models that attempted to predict counts of positive cases from 1 to 4 weeks ahead. The latter model was chosen since allows to account for both short-term and long-term dependence together with

DISCUSSION
Motivated by the need to deliver reliable estimates of the COVID-19 contagion dynamics, in this paper we have provided novel Bayesian Poisson autoregressive models, with time-varying coefficients. While time-varying coefficients allow to better take into account of nonlinear growth patterns, the Bayesian framework allows to incorporate the uncertainty of model outputs. The addition of policy intervention variables allows to better understand the timing, sign, and magnitude of the policy measures undertaken by the governments to face the pandemic.
The application of the model to two paradigmatic cases as Italy and the United States reveals that some combinations of policy interventions together with the first phase of the vaccination campaign are associated on average with a reduction in the incidence of COVID-19. In detail, findings suggest beneficial effects of the closure of schools and workplaces for both countries. We also found that earlier implementation of restrictions results in a greater reduction in the incidence of COVID-19 for Italy. The same conclusion seems not to be shared by the United States.
Accuracy evaluation of our model prediction intervals shows good in-sample and out-of-sample performance. With respect to the latter, since our tvBINGARCHX model has also been employed within the ECDC European Forecast Hub initiative, low error with respect to 1-4 weeks ahead forecasts confirm the previous results.

6.1
Literature comparison Agosto et al. (2021) extended the log-linear Poisson autoregression (Geir et al., 2022) to epidemiologic contagion. The frequentist framework of our paper (Section 5.1) refers to this model, being the starting point of our analysis and the model against which comparisons have been made with the Bayesian framework. Unlike Agosto et al. (2021), policy covariates have been introduced. The Bayesian framework instead introduces the time-varying models of Roy and Karmakar (2021) and evaluates their performances against the extended structure with time-varying policy covariates. Time constant coefficient models have also been added to have a direct comparison with Agosto et al. (2021) in the Bayesian framework. So, the model structure has been recovered and extended from these two contributions in terms of policy covariates.
On the other hand, in terms of policy evaluation, Bo et al. (2021) and Islam et al. (2020) have been recalled as a starting point for policy definition and evaluation. The effectiveness of NPIs was evaluated by Bo et al. (2021) in terms of timevarying reproduction number (Rt) while Islam et al. (2020) evaluated policies in terms of counts, more similar to our approach. We confirm results of Bo et al. (2021) and Islam et al. (2020) on the effectiveness of distancing policies but we extend their work by introducing and quantifying the effectiveness of vaccines. Similar Islam et al. (2020), we introduce lag effects and a combination of policies, resulting in a better evaluation of policy effects.

Strengths and limitations of this study
The strength of our study relies on both methodological extension and epidemic monitoring. The structure of our proposed methods shows both extensions of previous works in terms of time-dependence structure (i.e., time-varying coefficients) and exogenous components of the model (i.e., policy covariates). Our model formulation is more flexible since it allows to consider both time dependence structure in terms of autocorrelation with previous counts but also with exogenous covariates. Alternative time lags allow to capture the timing of the effect of public interventions on the trajectory of the pandemic. Furthermore, the time-varying coefficient framework allows both to model rapid changes of COVID-19 spread and address overdispersion issues. At the same time, our study is also subject to several limitations. First, data quality could be an issue for forecasting models, especially in terms of the measurement of policy covariates. We relied on OxCGRT variables to represent government interventions, but it is really challenging to collect information on the exact date, nature, and extent of policies undertaken by the different countries. Measurement with error could result in higher model bias. The discussion could be extended in the same way with regard to pandemic counts since day-to-day reports are usually made up of errors. In the end, real-time forecasts of cases are really challenging, given data issues in COVID-19 databases and intrinsic model error together with changing policies, population behavior, and testing procedures.

Implications for policy-makers
Quantifying the impact of NPIs is a critical step in ensuring a strong public health response to COVID-19 in the first phase of the vaccination campaign. Risk factors for cluster formation are expected to be similar across countries, including closed and poorly ventilated locations together with crowded places where people have limited possibility for physical distancing. High-risk environments include venues, events, and activities with certain environmental conditions. As a result, the findings suggest beneficial effects of the closure of schools and workplaces for both countries, the latter representing a possible option for reducing transmission but with a large-scale economic and social impact.
Mobility plays an important role in human-to-human transmission. However, given the difference in terms of population size, limiting people's movement and gathering by requiring public transport closures in a combination of stay-at-home orders, and internal movement restrictions allows to decrease municipal (and regional) mobility fluxes and, consequently, the infection rate in Italy. On the other hand, requiring canceling gatherings and events at the stricter level (between 11 and 110 people or less) is effective in large countries such as the United States.
We also found that earlier implementation of restrictions results in a greater reduction in the incidence of COVID-19 in Italy. The impact on the trajectory of infection changes for the United States, being more immediate. However, more factors need to be accounted for to obtain more accurate conclusions for both countries.
While NPIs reduce transmission rates, vaccination primarily reduces symptomatic disease, hospitalization, and death from infection, as well as infection rates to a lesser amount. Our findings show that vaccination is effective on epidemic evolution but in combination with some NPIs, advocating for the need to keep NPIs in place during the first phase of the vaccination campaign.
However, our results regarding several health system interventions need to be cautiously interpreted. The noneffectiveness of some NPIs can be attributed to several factors. For example, if the outcome of the study is related to the number of cases, as case detection improves through more effective testing and tracing, the number of cases reported will increase without representing a real rise of cases. Finally, for both countries, we have found a reverse effect of some policies, consistent with the fact that policy measures are triggered by the increase in COVID-19 counts.

Further studies
Further extensions of the work in the paper should consider the inclusion of finer policy measurements, as new databases become available; and, possibly, the inclusion of spatial dependence in the model, especially for large countries such as the United States, as longer time series become available. For all countries, including Italy and the United States, it would also be important to include policy and contagion data at the subnational and regional level, if available, to better capture spatial variation.

A C K N O W L E D G M E N T S
This research has received funding from the European Union's Horizon 2020 research and innovation program "PERISCOPE: Pan European Response to the ImpactS of COvid-19 and future Pandemics and Epidemics" [101016233, H2020-SC1-PHE_CORONAVIRUS-2020-2-RTD]. The project is the result of a close collaboration among the three authors; in addition, Paolo Giudici's main role has been in ideating and supervising the work; coding and model specification, and data analysis is mainly due to Barbara Tarantino together with Arkaprava Roy with regard to coding and revision of the paper. The Authors thank the European COVID-19 Forecast Hub of the European Center for Disease Prevention and Control, to which the forecasts based on this paper are sent every week, for useful comments, feedback, and discussions. Open Access Funding provided by Universita degli Studi di Pavia within the CRUI-CARE Agreement.

C O N F L I C T O F I N T E R E S T
The authors have declared no conflict of interest.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data from the case study are available as Supporting Information on the journal's web page.

O P E N R E S E A R C H B A D G E S
This article has earned an Open Data badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available in the Supporting Information section. This article has earned an open data badge "Reproducible Research" for making publicly available the code necessary to reproduce the reported results. The results reported in this article could fully be reproduced.