Inference for Network Count Time Series with the R Package PNAR

We introduce a new R package useful for inference about network count time series. Such data are frequently encountered in statistics and they are usually treated as multivariate time series. Their statistical analysis is based on linear or log linear models. Nonlinear models, which have been applied successfully in several research areas, have been neglected from such applications mainly because of their computational complexity. We provide R users the flexibility to fit and study nonlinear network count time series models which include either a drift in the intercept or a regime switching mechanism. We develop several computational tools including estimation of various count Network Autoregressive models and fast computational algorithms for testing linearity in standard cases and when non-identifiable parameters hamper the analysis. Finally, we introduce a copula Poisson algorithm for simulating multivariate network count time series. We illustrate the methodology by modeling weekly number of influenza cases in Germany.


Introduction
Many data examples are frequently observed as multivariate counting processes recorded over a time-span and with known relations among observations (e.g epidemiological data with geographical distances between different areas).An important objective then is to study the effect of a known network to the observed data.This motivates a great amount of interest to network time series models; see Zhu et al. (2017) who developed continuous Network Autoregressive models (abbreviated as NAR).For these models, the observed variable Y , for the node i at time t, is denoted by Y i,t , and it is assumed to depend on the past value of the variable for the node itself, say Y i,t−1 , and of the past values of the average of all its neighboring variables.The unknown parameters of the model are estimated by Least Squares estimation (LS).This work was advanced by Armillotta and Fokianos (2022a) who developed linear and log-linear Poisson Network Autoregression model (PNAR) for multivariate count distributed data.The joint dependence among different variables is specified by a copula construction (Fokianos et al., 2020, Sec. 2).In addition, Armillotta and Fokianos (2022a) have further established parametric estimation under the framework of quasi maximum likelihood inference (see Wedderburn (1974), Gourieroux et al. (1984)) and associated asymptotic theory when the network dimension increases.For the linear model case and in the context of epidemiology, related applied work has been developed by Held et al. (2005) and it was extended by Paul et al. (2008); Paul and Held (2011); Held and Paul (2012); Meyer and Held (2014) and Bracher and Held (2020).
The previous contributions impose linearity (or log-linearity) of the model, which can be a restrictive assumption for real world applications.For example, existence of different underlying states (e.g.exponentially expanding pandemic/ dying out pandemic) implies that different regime switching data generating processes should be applied and this fits the framework we consider.Recently, Armillotta and Fokianos (2022b) specified a general nonlinear Network Autoregressive model for both continuous and discrete-valued processes, establishing also related theory.In addition, the authors study testing procedures for examining linearity of NAR model against specific nonlinear alternatives by means of a quasi score test statistic.This methodology was developed with and without the presence of identifiable parameters under the null hypothesis.
Even though there exists sufficient statistical methodology for PNAR models, there has been a lack of up-to-date software for implementing their analyses.The aim of this work is to fill this gap by introducing the new package PNAR (Tsagris et al., 2023) and to demonstrate its usefulness for count network data analysis.Related R packages do not provide tools for estimating nonlinear models and applying associated testing procedures.The package GNAR (Leeming et al., 2023;Knight et al., 2020) studies Generalized NAR models (GNAR); this is a linear NAR model which takes into account the effect of several layers of connections between the nodes of the network.This package deals with continuous-valued time series and does not contain tools for testing linearity.PNAR complements GNAR as it provides additional methodology for testing and inference about nonlinear discrete-valued network models.
Package surveillance (Höhle et al., 2022;Meyer et al., 2017) fits only linear models for spatialtemporal disease counts with Poisson or Negative Binomial distribution and with an autoregressive network effect.The package does accommodate various structural break-point tests but it does not contains functions for testing linearity and for log-linear model fitting.Moreover, standard errors of estimated parameters are computed by considering the quasi-likelihood as the true likelihood of the model (Paul et al., 2008, Sec. 2.3).However, if the count time series are cross-sectional dependent, as it is usually the case, the likelihood function is misspecified and the obtained standard errors are not consistently estimated.
The PNAR package provides several advancements to the state of the art software: i) efficient estimation of (log-)linear models with proper robust standard errors accounting for possible model mispecification; ii) appropriate functions for testing linearity when a parameter is either identifiable of non-identifiable, under then null hypothesis, by providing appropriate p-value bounds and bootstrap approximations; iii) new algorithms for generating (log-)linear and nonlinear network count time series models.
The paper is organized as follows.The next section introduces linear and log-linear network time series autoregressive models for count data.Details about inference for unknown model parameters are provided.An application to estimation of weekly number of influenza A & B cases from two Southern German states is given and some further model aspects are discussed.Then we focus on non-linear models and associated testing theory.Results concerning score tests for testing linearity in NAR models are discussed and applied to influenza data.We also address the issue of computational speed.A short section discussing simulation of network count time series shows the usefulness of this methodology.The paper concludes with a short discussion.

Poisson network models
Consider a known network with N nodes, indexed by i = 1, . . .N .The neighborhood structure of such a network is completely described by its adjacency matrix, say A = (a ij ) ∈ Rx N ×N where a ij = 1, if there is a directed edge from i to j, say i → j, and 0 otherwise.Undirected graphs are allowed (A = A ′ ), which means that the edge between two nodes, i and j, has no specific direction (say i ∼ j).This is common in geographical and epidemic networks (e.g.district i shares a border with district j, patient i has a contact with patient j).Self-relationships are excluded i.e. a ii = 0 for any i = 1, . . ., N .
Let Y be a count variable measured on each node of the network (i = 1, . . ., N ), over a window of time (t = 1, . . ., T ).The data is a N -dimensional vector of time series which is observed over the domain t = 1, 2 . . ., T ; in this way, a univariate time series is observed for each node, say Y i,t , with corresponding conditional expectation λ i,t .Denote by λ t = E(Y t |F t−1 ) with λ t = (λ 1,t , . . ., λ i,t , . . ., λ N,t ) ′ the conditional expectation vector of the counts with respect to their past history F t−1 .The following linear autoregressive network model takes into account the known relations between nodes where n i = j̸ =i a ij is the total number of connections starting from the node i, such that i → j; called out-degree.We call (1) linear Poisson Network Autoregression of order 1, abbreviated by PNAR(1); (Armillotta and Fokianos, 2022a).From the left hand side equation of (1), we observe that the process Y i,t is assumed to be marginally Poisson but the joint process depends upon a copula function described in simulations at the end of the paper.Note that β 0 , β 1 , β 2 > 0 since the conditional mean of the Poisson is positive.Model (1) postulates that, for every single node i, the marginal conditional mean of the process is regressed on: • the average count of the other nodes j ̸ = i which have a connection with i; the parameter β 1 is called network effect, as it measures the average impact of node i's connections; • the past count of the variable itself for i; the coefficient β 2 is called autoregressive effect because it provides an estimator for the impact of past count Y i,t−1 .
Model (1) implies that only nodes directly followed by the focal node i (i.e.i → j), possibly, have an impact on its mean process of counts.It is a reasonable assumption in many applications; for example, in a social network the activity of node k, which satisfies a ik = 0, does not affect node i.
Hence, (1) measures the effect of a network to the observed multivariate count time series.Moreover, the model accommodates different types of network connectivity i.e. a i,j does not necessarily take the values 1-0 (connected-not connected).For example, a i,j = 1/d i,j where d i,j is some measure of distance between node i and node j and a i,i = 0.In this way the network effect becomes a spatial network component; see the last paragraph of Knight et al. (2020, p.3) for a discussion about a similar set of weights.

Inference
Model (2), or (3), depends on the m-dimensional vector of unknown parameters θ = (β 0 , β 11 , . . ., β 1p , β 21 , . . ., β 2p , δ 1 , . . ., δ q ) ′ , with m = 1 + 2p + q.We use of quasi-maximum likelihood methodology for estimation of θ; see Wedderburn (1974) and Gourieroux et al. (1984).The Quasi Maximum Likelihood Estimator (QMLE) is the vector of parameters θ maximizing the function which is the so called pooled Poisson log-likelihood (up to a constant).Note that (4) is not necessarily the true log-likelihood of the process but it serves as an approximation.In particular, ( 4) is the log-likelihood function that would have been obtained if all time series were contemporaneously independent.However, the QMLE is not computed under the assumption of independence because ( 4) is simply a working log-likelihood function.The choice of maximizing ( 4) is justified for several reasons: i) full likelihood based on the joint process is complex (see the last section); ii) the optimization of (4) guarantees consistency and asymptotic normality of QMLE for the true parameter vector θ 0 ; iii) the QMLE is asymptotically equivalent to the MLE if the true probability mass function belongs to the linear exponential family (Gourieroux et al., 1984); iv) simplified computations entailing increased speed for estimation.Robustness of the QMLE in finite samples has been verified by Armillotta and Fokianos (2022a) through extensive simulation studies.
When considering the linear model ( 2), the score function is is the error sequence.Then, the empirical Hessian and conditional information matrices are given, respectively, by is the conditional covariance matrix evaluated at θ.Under suitable assumptions, Armillotta and Fokianos (2022a) when N → ∞ and T → ∞, where H and B are the theoretical limiting Hessian and information matrices, respectively, evaluated at the true value θ = θ 0 .Then, a suitable estimator for the standard errors of θ is the square-rooted main diagonal of the empirical "sandwich" covariance matrix, . Closely related works to ours have employed (4) for inference ; Paul et al. (2008) and Paul and Held (2011), among others.However, in such works (4) is viewed as the true log-likelihood of the model and standard errors are computed by using the which underestimates the real source of variation of the parameters when cross-section dependence among counts is present; see (6) which depends on the conditional covariance matrix of the process Y t .The package PNAR returns robust standard errors as independence among counts is not assumed for their calculation.Similar theory holds for the log-linear model (3); details can be found in the aforementioned works.

Influenza data
To illustrate the use of PNAR we apply the methodology to the dataset fluBYBW from the surveillance package (Meyer et al., 2017).This dataset includes information about the weekly number of in- Any type of non-negative matrix with zero main diagonal can be used as a valid adjacency matrix W. For instance consider weighted networks or inverse distance matrices, etc; see the last section for some alternatives when generating data.Optimization of ( 4) is implemented under the non-negativity constraint of coefficients satisfying the stationary condition.This is a nonlinear constrained optimization problem solved by means of a Sequential Quadratic programming (SQP) gradient-based algorithm (Kraft, 1994) of package nloptr (Ypma and Johnson, 2022).By default, the optimization is constrained in the stationary region; this can be removed by setting the option uncons = TRUE although this is not suggested because large sample properties of the estimators have been developed within the stationary region.The function lin estimnarpq() has three additional features: • maxeval: the maximum number of iterations for the optimization (the default is 100); • xtol rel: relative tolerance for the optimization termination condition (the default is 1e-8); • init: starting value of the optimization (the default is NULL).
When init = NULL starting values are computed internally by ordinary LS using a lower barrier value of 0.01 because the regression coefficients cannot assume negative values.However, users can provide their own initial values through the argument init.
The function lin estimnarpq returns as output a list consisting of the estimated coefficients, their associated standard errors, a z-test statistics with p-values, the score function evaluated at the optimum, the maximized log-likelihood, and the usual Akaike and Bayesian Information Criteria (AIC, BIC) accompanied by the Quasi IC (QIC) (Pan, 2001) which takes into account the fact that the log-likelihood (4) is a quasi log-likelihood; see Table 1 for the results, which are obtained in less than a second (see also Table 2).The score computed at the optimum values is of order 1e-5, on average, indicating successful convergence of the algorithm.
All the estimated coefficients are positive and significantly different from 0. The autoregressive effect β 2,h shows higher magnitude with respect to the network effects β 1,h since past counts of the same district are, in general, more informative than the neighboring cases.Both network and autoregressive parameters when p = 1 have a larger magnitude when compared to the corresponding coefficients at p = 2.This can be explained since influenza has an incubation period of only 1-4 days (with an average of 2 days) and a patient is still contagious for no more than 5-7 days after becoming sick1 so the case counts at first lag are more informative than the ones at second lag which are still important.The population covariate is significant with positive effect.Standard errors are computed by using the sandwich estimator.
To select the model order p, for the PNAR model ( 2), PNAR includes a function for estimating model parameters for a range of lag values.By default p ∈ {1, 2, . . ., 10}.This function returns the scatter plot of any IC (default is QIC) versus the lag order, for example lin_ic_plot(y = flu, W = A_flu, p = 1:10, Z = pop, ic = "AIC") Figure 2 shows the output of lin ic plot() for the case of AIC.Plots for the cases of BIC and QIC are similar and not shown.All information criteria point to the model with p = 9.However, from the corresponding estimation results reported in the Appendix, almost all β coefficients, which correspond to lags p ≥ 3, are close to zero and non significant.Therefore, we decide to retain p = 2, for parsimony.
We compare estimation of standard errors after fitting the linear PNAR model to influenza data using PNAR and surveillance packages (covariate Z is excluded).The latter follows the standard error estimation according to the approach described in Paul et al. (2008) and Paul and Held (2011).Comparing the above outputs, note that estimated coefficients are identical (subject to rounding errors) but standard errors, obtained by surveillance package, underestimate their value as it was explained at the end of the inference section.In addition, PNAR is more user-friendly for fitting (2) and it is almost twice as fast.Indeed, the average estimation time over 10 calls of the same function is 0.42 seconds in comparison to 0.22 by the PNAR package (see Table 2).
The usefulness of the PNAR package is not limited to epidemiological data.For example, the dataset crime, which is already built in the package, contains monthly number of burglaries within the census blocks on the south side of Chicago during 2010-2015 and includes a network matrix, called crime W, connecting two blocks sharing a border.The documentation of the functions lin estimnarpq() and log lin estimnarpq() provides an example of PNAR models applied to crime data.

Extending linearity
In this section, we give some motivating examples of nonlinear models for time series of networks as introduced by Armillotta and Fokianos (2022b).In addition, we provide testing procedures for testing the linearity assumption.For ease of presentation, denote by X i,t = n −1 i N j=1 a ij Y j,t the average neighbor network mean.When a drift in the intercept term of (2) is introduced, the following nonlinear Intercept Drift model, ID-PNAR(p, q), is obtained where γ ≥ 0. Model (7) shares similar to a linear model when the parameter γ takes small values and for γ = 0 reduces to (2).Instead, when γ takes values away from zero, model ( 7) introduces a perturbation, deviating from the linear model, which depends on the network mean at lag (t − d), where d is an additional delay parameter such that d = 1, 2, . . ., p. Model ( 7) is directly applicable when the baseline effect, β 0 , varies over time as a function of the network.
Another interesting class of nonlinear models are regime switching models, i.e. models allowing for the dynamics of the count process to depend on different regimes (e.g.exponentially expanding pandemic/ dying out pandemic).We give two such examples of models whose specification is based either on a smooth distortion or an abrupt transition.With the same notation as before, the Smooth Transition PNAR model, ST-PNAR(p, q) assumes a smooth transition between two regimes and it is defined by where γ ≥ 0 and α h ≥ 0, for h = 1, . . ., p.This models introduces a smooth regime switching behavior of the network effect making it possible to vary smoothly from β 1,h to β 1,h + α h , as γ varies from large to small values.The additional delay parameter d determines the time of nonlinear transition that can be chosen.When α h = 0, for h = 1, . . ., p in ( 8), the linear PNAR model ( 2) is recovered.In some applications the transition between regimes may be abrupt (e.g.financial market crashes).For this reason, we consider the Threshold PNAR model, T-PNAR(p, q), which is defined by where I(•) is the indicator function and γ ≥ 0 is now threshold parameter.Moreover, α 0 , α 1,h , α 2,h ≥ 0, for h = 1, . . ., p.When α 0 = α 11 = ... = α 2p = 0, model (9) reduces to (2).
For all nonlinear models, estimation of the unknown parameters is based on QMLE, following the discussion in the inference section.Therefore, analogous conclusions about estimation of regression parameters, their standard errors and model selection apply to the case of nonlinear models (7)-( 9).
Further details can be found in Armillotta and Fokianos (2022b).

Standard implementation of testing linearity
Testing linearity against several specific alternatives offers guidance about the type of nonlinear model to be fitted.Moreover, in certain cases where the linear model is nested within a nonlinear model, some nonlinear parameters may be inconsistently estimated (see the case ST-PNAR and T-PNAR models below) so testing linearity prevents incorrect estimation.
Then the restricted estimator is denoted by θ = ( β0 , β11 , . . ., β2p , δ1 , . . ., δq ) ′ and it is usually simpler to compute.The quasi score test statistic is given by (2) m 2 distribution, where m 2 is the number of nonlinear parameter tested (Armillotta and Fokianos, 2022b).Then, in case of model ( 7), m 2 = 1 and we can compute the p-value of the test statistic ( 10) by p = P(χ 2 1 ≥ M ), where M is the observed value of the test statistic LM T .

Non-standard implementation of testing linearity
Suppose now we wish to test linearity against (8).Then, considering the hypotheses H 0 : α 1 = • • • = α p = 0 versus H 1 : α h ̸ = 0 for some h = 1, . . ., p, we note that this testing problem is nonstandard because it is not possible to estimate the value of γ under H 0 .Note that the parameter γ exists in the partition of the score function (5) related to nonlinear parameters, in particular in ∂λ i,t (θ)/∂α h , and its associated covariance matrix.Hence, all relevant quantities for computing (10) are functions of γ; that is S (2) T ( θ, γ), Σ T ( θ, γ) and LM T (γ).The model is then subject to non-identifiable parameter γ under the null.So if the true model is the linear PNAR model but a ST-PNAR is estimated instead the smoothing parameter γ will not be consistently estimated.
Analogous conclusions hold for testing linearity against the T-PNAR model ( 9), where the threshold parameter γ is not identifiable under the null.When this issue arises, the standard theory does not apply and a chi-square type test is not suitable any more; see Davies (1987), Hansen (1996) and Armillotta and Fokianos (2022b), among others.It is clear that the value of the test changes by varying γ ∈ Γ, where Γ is some domain.A summary function of the test, computed under different values of γ, is then routinely employed in applications; a typical choice is g T = sup γ∈Γ LM T (γ); see Armillotta and Fokianos (2022b) who established the convergence of g T to g = sup γ∈Γ LM (γ), where g is a function of a chi-square process, LM (γ).The values of the latter asymptotic distribution cannot be tabulated, as they depends on unknown values of γ.Hence, we give methodology for computing p-values of such sup-type test statistic since they cannot be obtained otherwise.
Davies' bound.Since the space Γ = [γ L , γ U ] is usually assumed to be a closed interval, in practice, we take Γ F = (γ L , γ 1 , . . ., γ l , γ U ) i.e. a grid of values for the non-identifiable parameter γ, and g T is obtained as the maximum of the tests LM T (γ) computed over Γ F .Davies (1987) showed that P sup where M is the value of the test statistic g T computed in the available sample, Γ(•) is the gamma function, and V is the approximated total variation Equation ( 11) shows how to approximate the p-values of the sup-type test in a straightforward way.Indeed, by adding to the tail probability of a chi-square distribution a correction term, which depends on the total variation of the process, we obtain the desired bound.This method is attractive for its simplicity and speed even when the dimension N of the network is large.However, the method approximate p-values with their bound (11) leading to a conservative test.In addition, (11) cannot be applied to the T-PNAR model ( 9), because the total variation requires differentiability of the asymptotic distribution LM (γ) under the null hypothesis (Davies, 1987, p. 36), a condition that is not met for the case of T-PNAR models.
Bootstrapping the test statistic.Based on the previous arguments, we suggest an alternative p-value approximation of the test statistic employing stochastic permutations (Hansen, 1996;Armillotta and Fokianos, 2022b)-see Algorithm 1.
Algorithm 1 Score bootstrap 1: Obtain the constrained QMLE of the linear model (2), say θ 2: for j = 1, . . ., J do 3: Generate ν t,j ∼ N (0, 1) end for Compute the test LM Optimize LM ν j T (γ) for γ and take g j T = sup γ∈Γ LM ν j T (γ) 9: end for An approximation of the p-values is obtained from step 10 of Algorithm 1, where g T is the value of the test statistic computed on the available sample.When the number of bootstrap replications J is large enough, p J T provides a good approximation to the unknown p-values of the test.Then, the null hypothesis H 0 is rejected if p J T is smaller than a given significance level.

Revisiting the influenza data
We now apply the testing methodology described in the previous sections to influenza data.Con- Next consider testing linearity against the ST-PNAR(2) alternative (8).In this case, the test behaves in a non-standard way so we will be using the Davies' bound p-value (DV) for the sup- Extreme values for the range of γ ∈ Γ F are computed internally in such a way that γ U and γ L are those values of γ where, on average, the smoothing function exp(−γX 2 i,t−d ) is equal to 0.1 and 0.9, respectively.In this way, during the optimization procedure, the extremes of the function domain are excluded.For more details see the PNAR manual (Tsagris et al., 2023).The user can specify different values for γ L and γ U using the arguments gama L, gama U and the number of grid values len (the default is 100).Results suggest again a deviation from linearity of the model, i.e.

dv2.z
Test for linearity of PNAR(p) versus the non-linear ST-PNAR(p) data: coefficients of the PNAR(p, q)/time series data/order/lag/covariates /lower gamma/upper gamma/length chi-square-test statistic = 35.074,df = 2, p-value = 9.076e-08 alternative hypothesis: At least one coefficient of the non-linear component is not zero Next, we apply Algorithm 1 to compute the bootstrap p-values for the sup-type test statistic.In this case, the observed value of sup γ∈Γ F LM (γ) has to be computed.Initially, perform a global optimization of the LM T (γ) for the ST-PNAR(p) model, with respect to the nuisance scale parameter γ by using Brent's algorithm (Brent, 1973) in the interval [gama L to gama U], (see previous discussion for their computation).To ensure global optimality, the optimization is performed on runs at len-1 consecutive equidistant sub-intervals and the global optimum is determined by the maximum over those sub-intervals.The default value for len is 10.Then using the function global optimise LM stnarpq with the same arguments we obtain the optimal γ value and the corresponding value of the test statistic: The above output gives (among other information) p J T of step 10 of Algorithm 1 and an alternative corrected unbiased estimator for the p-value which is cp J T = (J + 1) −1 J j=1 I(g j T ≥ g T ) + 1 .Like in the case of DV p-value, linearity is rejected.
We work analogously for testing the PNAR(2) model versus the T-PNAR(2) model ( 9).Note that optimization of LM T (γ), in this case, is based on gama L and gama U which are obtained (by default) as the mean over i = 1, . . ., N of 20% and 80% quantiles of the empirical distribution of the network mean X i,t for t = 1, . . ., T .In this way, during the optimization process, the indicator function I(X i,t−d ≤ γ) avoids values close to 0 or 1.Alternatively, their value can be supplied by the user.The functions used are analogous to the functions used for the ST-PNAR(2).The test does not reject the null hypothesis of linearity, in the case of a threshold model.Overall, the analysis shows that the linear PNAR model may not be a suitable model to fit such epidemic data and nonlinear alternatives should be considered.In particular, evidence of a nonlinear drift in the intercept and of a regime switching mechanism is detected.In addition, it appears that a smooth regime switching mechanism might be more appropriate for the data.

Computational speed
Package PNAR is quite efficient in terms of computational speed, especially for estimation problems.
We have employed the Rfast and Rfast2 packages (Papadakis et al., 2023a,b) wherever possible to ensure computational speed.We run each function 10 times and compute the average time required to be executed.We use a laptop computer equipped with Intel Core i7 processor (3.00GHz) and 16 GB of RAM.Results are given in Table 2   robust standard errors, test statics and computational algorithms for testing model linearity and new simulation methodology for generating (log-)linear and nonlinear network count time series models.We showed that all these tasks are executed with a minimal computational effort.
There are a number of possible developments for PNAR.One possibility is to include several other nonlinear models and develop related linearity tests.In addition, developing a negative binomial quasi-likelihood estimation method offers more flexibility to model fitting.Alternative ways to compute p-values, for example employing different bootstrap approximation, may also be considered.
Further simulation methods for generating network count time series are easily accommodated by suitable modification of the copula Poisson algorithm.All these extensions provide users with new set of tools for inference in the broad framework of multivariate discrete-valued time series models.

Appendix
Output from estimation of PNAR model (2) with lag p = 9 and population covariate: Figure 1: Quarterly flu cases in two Southern German states Bavaria and Baden-Wuerttemberg for 2007.

Figure 2 :
Figure 2: Scatter plot of AIC for PNAR(p) model versus p.
J′ , where J = (O m 2 ×m 1 , I m 2 ), I s is a s × s identity matrix and O a×b is a a × b matrix of zeros.Σ T ( θ) is a the estimator for the unknown covariance matrix Σ = Var[S (2) T ( θ)].It can be proved that the quasi score test (10) converges, asymptotically, to a χ 2 testing linearity of the PNAR(2) model against the nonlinear ID-PNAR(2) (7) with d = 1.Analogous results have been obtained for d = 2 and therefore are omitted.The quasi score test (10) is computed by id2.z <-score_test_nonlinpq_h0(b = est2.z$coefs[,1], y = flu, W = A_flu, p = 2, d = 1, Z = pop) id2.zLinearity test against non-linear ID-PNAR(p) model data: coefficients of the PNAR(p, q)/time series data/order/lag/covariates/ chi-square-test statistic = 7.2318, df = 1, p-value = 0.007162 alternative hypothesis: True gamma parameter is greater than 0 where the first argument requires estimates under the null hypothesis H 0 : γ = 0 which have been obtained already.The rest of arguments follow previous syntax.For this testing problem the test is asymptotically chi-square distributed with 1 degree of freedom.The output lists the test statistic value and its corresponding p-value.There is strong indication to reject the linear model in favour of model (7).

Table 1 :
Estimation of linear PNAR model (2) for p = 1, 2 and Z = pop.Standard errors of coefficients are given in parentheses.

Table 2 :
Average computation times (in seconds) for the functions called in the text.

Table 3 :
A brief overview of the main functions in PNAR.Score test of linear PNAR versus ST-PNAR model by (11).globaloptimise LM stnarpq() Maximize test statistic of ST-PNAR for nuisance parameters.score test stnarpq j() Bootstrap score test of linear PNAR versus ST-PNAR model.
global optimise LM tnarpq() Maximize test statistic of T-PNAR for nuisance parameters.score test tnarpq j() Bootstrap score test of linear PNAR versus T-PNAR model.