Title:  Poisson Network Autoregressive Models 

Description:  Quasi likelihoodbased methods for estimating linear and loglinear Poisson Network Autoregression models with p lags and covariates. Tools for testing the linearity versus several nonlinear alternatives. Tools for simulation of multivariate count distributions, from linear and nonlinear PNAR models, by using a specific copula construction. References include: Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. <arXiv:2104.06296>. Armillotta, M. and K. Fokianos (2022b). Testing linearity for network autoregressive models. <arXiv:2202.03852>. Armillotta, M., Tsagris, M. and Fokianos, K. (2022c). The Rpackage PNAR for modelling count network time series. <arXiv:2211.02582>. 
Authors:  Michail Tsagris [aut, cre], Mirko Armillotta [aut, cph], Konstantinos Fokianos [aut] 
Maintainer:  Michail Tsagris <mtsagris@uoc.gr> 
License:  GPL (>= 2) 
Version:  1.6 
Built:  20240206 07:49:39 UTC 
Source:  CRAN 
Quasi likelihoodbased methods for estimating linear and loglinear Poisson Network Autoregression models with p lags and covariates. Tools for testing the linearity versus several nonlinear alternatives. Tools for simulation of multivariate count distributions, from linear and nonlinear PNAR models, by using a specific copula construction. References include: Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296. Armillotta, M. and K. Fokianos (2022b). Testing linearity for network autoregressive models. https://arxiv.org/abs/2202.03852. Armillotta, M., Tsagris, M. and Fokianos, K. (2022c). The Rpackage PNAR for modelling count network time series. https://arxiv.org/abs/2211.02582.
Package:  PNAR 
Type:  Package 
Version:  1.6 
Date:  20231009 
License:  GPL(>=2) 
Disclaimer: Dr Mirko Armillotta and Konstantinos Fokianos wrote the initial functions. Dr Tsagris modified them, created the package and he is the maintainer.
We would to like to acknowledge Manos Papadakis for his help with the "htest" class object and S3 methods (print() and summary() functions).
Michail Tsagris, Mirko Armillotta and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Armillotta, M. and K. Fokianos (2022b). Testing linearity for network autoregressive models. https://arxiv.org/abs/2202.03852
This function generates a network from the Stochastic Block Model with $K$
blocks.
adja(N, K, alpha, directed = FALSE)
N 
The number of nodes on the network. 
K 
The number of blocks. Each block has dimension 
alpha 
The network density. A value in 
directed 
Logical scalar, whether to generate a directed network or not. If TRUE a directed network is generated. 
For each pair of nodes it performs a Bernoulli trial with values 1 "draw an edge", 0 "otherwise".
The probabilities of these trials are bigger if the two nodes are in the same block, lower otherwise, and they are specified based on the number of nodes on the network $N$
and network density $alpha$
:
Probability to draw an edge for a pair of nodes in the same block: $\alpha*N^{0.3}$
.
Probability to draw an edge for a pair of nodes in different blocks: $\alpha*N^{1}$
.
A rownormalized nonnegative matrix describing the network. The main diagonal entries of the matrix are zeros, all the other entries are nonnegative and the sum of elements over the rows equals one.
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Faust, K. and S. Wasserman (1992). Blockmodels: Interpretation and evaluation. Social Networks, 14, 561.
W < adja(N = 20, K = 5, alpha = 0.1)
This function generates a network from the ErdosRenyi model.
adja_gnp(N, alpha, directed = FALSE)
N 
The number of nodes on the network. 
alpha 
The network density. A value in 
directed 
Logical scalar, whether to generate a directed network. If TRUE a directed network is generated. 
For each pair of nodes it performs a Bernoulli trial with values 1 "draw an edge", 0 "otherwise".
Each trial has the same probability of having an edge; this is equal to $\alpha*N^{0.3}$
, specified based on the number of nodes on the network $N$
and the network density alpha
.
A rownormalized nonnegative matrix describing the network. The main diagonal entries of the matrix are zeros, all the other entries are nonnegative and the maximum sum of elements over the rows equals one.
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Erdos, P. and A. Renyi (1959). On random graphs. Publicationes Mathematicae, 6, 290297.
W < adja_gnp(N = 20, alpha= 0.1)
Monthly number of burglaries on the south side of Chicago (552 blocks) during 20102015 (72 temporal observations).
crime
A time series object ("ts" class) with multivariate time series, a matrix with 72 rows and 552 columns.
Clark and Dixon (2021), available at https://github.com/nick3703/ChicagoData.
Clark, N. J. and P. M. Dixon (2021). A class of spatially correlated selfexciting statistical models. Spatial Statistics, 43, 118.
crime_W, lin_estimnarpq, log_lin_estimnarpq
data(crime)
data(crime_W)
mod1 < lin_estimnarpq( crime, crime_W, p = 1)
mod2 < log_lin_estimnarpq( crime, crime_W, p = 1)
Nonnegative rownormized adjacency matrix describing the network structure between Chicago census blocks.
crime_W
A matrix with 552 rows and 552 columns.
Clark and Dixon (2021), available at https://github.com/nick3703/ChicagoData.
Clark, N. J. and P. M. Dixon (2021). A class of spatially correlated selfexciting statistical models. Spatial Statistics, 43, 118.
crime, lin_estimnarpq, log_lin_estimnarpq
data(crime)
data(crime_W)
mod1 < lin_estimnarpq(crime, crime_W, p = 1)
mod2 < log_lin_estimnarpq(crime, crime_W, p = 1)
This function counts the number of events within a specified time.
getN(x, tt = 1)
x 
A matrix of (positive) interevent times. 
tt 
A positive time. 
The number of events within time tt
(possibly 0), for each column of
x
.
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
rcopula, poisson.MODpq, poisson.MODpq.log
x < rcopula(n = 100, N = 50, rho = 0.3)
getN(x)
Global optimization of the linearity test statistic for the Smooth Transition
Poisson Network Autoregressive model of order $p$
with $q$
covariates
(STPNAR($p$
)) with respect to the nuisance scale parameter $\gamma$
.
global_optimise_LM_stnarpq(gama_L = NULL, gama_U = NULL, len = 10, b, y, W,
p, d, Z = NULL, tol = 1e9)
gama_L 
The lower value of the 
gama_U 
The upper value of the 
len 
The number of increments to consider for the 
b 
The estimated parameters from the linear model, in the following order:
(intercept, network parameters, autoregressive parameters, covariates).
The length of the vector should be 
y 
A 
W 
The 
p 
The number of lags in the model. 
d 
The lag parameter of nonlinear variable (should be between 1 and 
Z 
An 
tol 
Tolerance level for the optimizer. 
The function optimizes the quasi score test statistic, under the null assumption of linearity, for testing linearity of Poisson Network Autoregressive model of order $p$
against the following STPNAR($p$
) model, with respect to the unknown nuisance parameter ($\gamma$
). For each node of the network $i=1,...,N$
over the time sample $t=1,...,TT$
$\lambda_{i,t}=\beta_{0}+\sum_{h=1}^{p}(\beta_{1h}X_{i,th}+\beta_{2h}Y_{i,th}+\alpha_{h}e^{\gamma X_{i,td}^{2}}X_{i,th})+\sum_{l=1}^{q}\delta_{l}Z_{i,l}$
where $X_{i,t}=\sum_{j=1}^{N}W_{ij}Y_{j,t}$
is the network effect, i.e. the weighted average impact of node $i$
connections, with the weights of the mean being $W_{ij}$
, the single element of the network matrix $W$
. The sequence $\lambda_{i,t}$
is the expectation of $Y_{i,t}$
, conditional to its past values.
The null hypothesis of the test is defined as $H_{0}: \alpha_{1}=...=\alpha_{p}=0$
, versus the alternative that at least one among $\alpha_{h}$
is not 0. The test statistic has the form
$LM(\gamma)=S^{'}(\hat{\theta},\gamma)\Sigma^{1}(\hat{\theta},\gamma)S(\hat{\theta},\gamma)$
where
$S(\hat{\theta},\gamma)=\sum_{t=1}^{TT}\sum_{i=1}^{N}\left(\frac{Y_{i,t}}{\lambda_{i,t}(\hat{\theta},\gamma)}1\right) \frac{\partial\lambda_{i,t}(\hat{\theta},\gamma)}{\partial\alpha}$
is the partition of the quasi score related to the vector of nonlinear parameters $\alpha=(\alpha_{1},...,\alpha_{p})$
, evaluated at the estimated parameters $\hat{\theta}$
under the null assumption $H_{0}$
(linear model) and $\Sigma(\hat{\theta},\gamma)$
is the variance of $S(\hat{\theta},\gamma)$
.
The optimization employes the Brent algorithm (Brent, 1973) applied in the interval from gama_L
to gama_U
. To be sure that the global optimum is found, the optimization is performed at (len
1) consecutive equidistant subintervals and then the maximum over them is taken as global optimum.
The values of gama_L
and gama_U
are computed internally as gama_L
$=\log(0.9)/X^{2}$
and gama_U
$=\log(0.1)/X^{2}$
, where $X$
is the overall mean of $X_{i,t}$
over the nodes $i=1,...,N$
and times $t=1,...,TT$
. Since the nonlinear function $e^{\gamma X_{i,td}^{2}}$
ranges between 0 and 1, by considering $X$
to be a representative value for the network mean, gama_U
and gama_L
would be the values of $\gamma$
leading the nonlinear switching function to be 0.1 and 0.9, respectively, so that in the optimization procedure the extremes of the function domain are excluded. Alternatively, their value can be supplied by the user. For details see Armillotta and Fokianos (2022b, Sec. 45).
A list including:
gama 
The optimum value of the 
supLM 
The value of the objective function at the optimum. 
int 
A vector with the extremes points of subintervals. 
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Armillotta, M. and K. Fokianos (2022b). Testing linearity for network autoregressive models. https://arxiv.org/abs/2202.03852
Armillotta, M., Tsagris, M. and Fokianos, K. (2022c). The Rpackage PNAR for modelling count network time series. https://arxiv.org/abs/2211.02582
Brent, R. (1973) Algorithms for Minimization without Derivatives. PrenticeHall, Englewood Cliffs N.J.
score_test_stnarpq_j, global_optimise_LM_tnarpq,
score_test_tnarpq_j
data(crime)
data(crime_W)
mod1 < lin_estimnarpq(crime, crime_W, p = 1)
b < mod1$coefs[, 1]
global_optimise_LM_stnarpq(b = b, y = crime, W = crime_W, p = 1, d = 1)
Global optimization of the linearity test statistic for the Threshold Poisson
Network Autoregressive model of order $p$
with $q$
covariates
(TPNAR($p$
)) with respect to the nuisance threshold parameter $\gamma$
.
global_optimise_LM_tnarpq(gama_L = NULL, gama_U = NULL, len = 10, b, y, W,
p, d, Z = NULL, tol = 1e9)
gama_L 
The lower value of the 
gama_U 
The upper value of the 
len 
The number of increments to consider for the 
b 
The estimated parameters from the linear model, in the following order:
(intercept, network parameters, autoregressive parameters, covariates).
The dimension of the vector should be 
y 
A 
W 
The 
p 
The number of lags in the model. 
d 
The lag parameter of nonlinear variable (should be between 1 and 
Z 
An 
tol 
Tolerance level for the optimizer. 
The function optimizes the quasi score test statistic, under the null assumption of linearity,
for testing linearity of Poisson Network Autoregressive model of order $p$
against the following TPNAR($p$
) model, with respect to the unknown nuisance parameter ($\gamma$
). For each node of the network $i=1,...,N$
over
the time sample $t=1,...,TT$
$\lambda_{i,t}=\beta_{0}+\sum_{h=1}^{p}\left[\beta_{1h}X_{i,th}+\beta_{2h}Y_{i,th}+(\alpha_{0}+\alpha_{1h}X_{i,th} +\alpha_{2h}Y_{i,th})I(X_{i,td}\leq\gamma)\right]+\sum_{l=1}^{q}\delta_{l}Z_{i,l}$
where $X_{i,t}=\sum_{j=1}^{N} W_{ij}Y_{j,t}$
is the network effect, i.e. the weighted average impact of node $i$
connections, with the weights of the mean being $W_{ij}$
, the single element of the network matrix $W$
, and $I()$
is the indicator function. The sequence $\lambda_{i,t}$
is the expectation of $Y_{i,t}$
, conditional to its past values.
The null hypothesis of the test is defined as $H_{0}: \alpha_{0}=\alpha_{11}=...=\alpha_{2p}=0$
, versus the alternative that at least one among $\alpha_{s,h}$
is not $0$
, for $s=0,1,2$
. The test statistic has the form
$LM(\gamma)=S^{'}(\hat{\theta},\gamma)\Sigma^{1}(\hat{\theta},\gamma)S(\hat{\theta},\gamma)$
where
$S(\hat{\theta},\gamma)=\sum_{t=1}^{TT}\sum_{i=1}^{N}\left(\frac{Y_{i,t}}{\lambda_{i,t}(\hat{\theta},\gamma)}1\right)\frac{\partial\lambda_{i,t}(\hat{\theta},\gamma)}{\partial\alpha}$
is the partition of the quasi score related to the vector of nonlinear parameters $\alpha=(\alpha_{0},...,\alpha_{2p})$
, evaluated at the estimated parameters $\hat{\theta}$
under the null assumption $H_{0}$
(linear model) and $\Sigma(\hat{\theta},\gamma)$
is the variance of $S(\hat{\theta},\gamma)$
.
The optimization employes the Brent algorithm (Brent, 1973) applied in the interval from gama_L
to gama_U
. To be sure that the global optimum is found, the optimization is performed at (len
1) consecutive equidistant subintervals and then the maximum over them is taken as global optimum.
The values of gama_L
and gama_U
are computed internally as the mean over $i=1,...,N$
of $20\%$
and $80\%$
quantile of the empirical distribution of the network mean $X_{i,t}$
for $t=1,...,TT$
. In this way the optimization is performed for values of $\gamma$
such that the indicator function $I(X_{i,td}\leq\gamma)$
is not always close to 0 or 1. Alternatively, their value can be supplied by the user. For details see Armillotta and Fokianos (2022b, Sec. 45).
A list including:
gama 
The optimum value of the 
supLM 
The value of the objective function at the optimum. 
int 
A vector with the extremes points of subintervals. 
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Armillotta, M. and K. Fokianos (2022b). Testing linearity for network autoregressive models. https://arxiv.org/abs/2202.03852
Armillotta, M., Tsagris, M. and Fokianos, K. (2022c). The Rpackage PNAR for modelling count network time series. https://arxiv.org/abs/2211.02582
Brent, R. (1973) Algorithms for Minimization without Derivatives. PrenticeHall, Englewood Cliffs N.J.
score_test_tnarpq_j, global_optimise_LM_stnarpq,
score_test_stnarpq_j
data(crime)
data(crime_W)
mod1 < lin_estimnarpq(crime, crime_W, p = 2)
b < mod1$coefs[, 1]
global_optimise_LM_tnarpq(b = b, y = crime, W = crime_W, p = 2, d = 1)
Estimation of the linear Poisson Network Autoregressive model of order $p$
with $q$
covariates (PNAR($p$
)).
lin_estimnarpq(y, W, p, Z = NULL, uncons = FALSE, init = NULL,
xtol_rel = 1e8, maxeval = 100)
y 
A 
W 
The 
p 
The number of lags in the model. 
Z 
An 
uncons 
logical, if TRUE an unconstrained optimization is run (default is FALSE). 
init 
A vector of starting values for the optimization algorithm. If this is NULL, the function computes them internally. 
xtol_rel 
The stopping tolerance of the optimization algorithm. 
maxeval 
The maximum number of evalutions the optimization algorithm will perform. 
This function performs constrained estimation of the linear Poisson NAR($p$
) model with $q$
nonnegative valued covariates, for each node of the network $i=1,...,N$
over the time sample $t=1,...,TT$
, defined as
$\lambda_{i,t}=\beta_{0}+\sum_{h=1}^{p}(\beta_{1h}X_{i,th}+\beta_{2h}Y_{i,th})+\sum_{l=1}^{q}\delta_{l}Z_{i,l},$
where $X_{i,t}=\sum_{j=1}^{N}W_{ij}Y_{j,t}$
is the network effect, i.e. the weighted average impact of node $i$
connections, with the weights of the mean being $W_{ij}$
, the single element of the network matrix $W$
. The sequence $\lambda_{i,t}$
is the expecation of $Y_{i,t}$
, conditional to its past values. The parameter $\beta_{0}$
is the intercept of the model, $\beta_{1h}$
are the network coefficients, $\beta_{2h}$
are the autoregressive parameters, and $\delta_{l}$
are the coefficients assocciated to the covariates $Z_{i,l}$
.
The estimation of the parameters of the model is performed by Quasi Maximum Likelihood Estimation (QMLE), maximizing the following quasi loglikelihood
$l(\theta)=\sum_{t=1}^{TT}\sum_{i=1}^{N}\left[Y_{i,t}\log\lambda_{i,t}(\theta)\lambda_{i,t}(\theta)\right]$
with respect to the vector of unknown parameters $\theta$
described above. The coefficients are defined only in the nonnegative real line.
By default, the optimization is constrained in the stationary region where $\sum_{h=1}^{p}(\beta_{1h}+\beta_{2h})<1$
; this can be removed by setting uncons = TRUE
. However, the model estimates might be inconsistent if the estimated parameters lie outside the stationary region.
The ordinary least squares estimates are employed as starting values of the optimization procedure. Robust standard errors and ztests are also returned.
A list with attribute class "PNAR" including:
coefs 
A matrix with the estimated QMLE coefficients, their standard errors their Ztest statistics and the relevant pvalues computed via the standard normal approximation. 
score 
The value of the quasi score function at the optimization point. It should be close to 0 if the optimization is successful. 
loglik 
The value of the maximized quasi loglikelihood. 
ic 
A vector with the Akaike information criterion (AIC), the Bayesian information criterion (BIC) and the Quasi information criterion (QIC). 
Alternatively, these can be printed via the function summary.PNAR
.
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Armillotta, M., Tsagris, M. and Fokianos, K. (2022c). The Rpackage PNAR for modelling count network time series. https://arxiv.org/abs/2211.02582
data(crime)
data(crime_W)
mod1 < lin_estimnarpq(crime, crime_W, p = 2)
summary(mod1)
Scatter plot of information criteria versus the number of lags in the linear Poisson Network Autoregressive model of order $p$
with $q$
covariates (PNAR($p$
)).
lin_ic_plot(y, W, p = 1:10, Z = NULL, uncons = FALSE, ic = "QIC")
y 
A 
W 
The 
p 
A vector with integer numbers, the range of lags in the model, for which the AIC, BIC and QIC will be computed. 
Z 
An 
uncons 
Logical, if TRUE an unconstrained optimization without stationarity constraints is performed (default is FALSE). 
ic 
The information criterion you want to plot, "QIC" (default value), "AIC" or "BIC". 
The function computes the AIC, BIC or QIC for a range of lag orders of the
linear Poisson Network Autoregressive model of order $p$
with $q$
covariates (PNAR($p$
)).
A scatter plot with the lag order versus either QIC (default), AIC or BIC, and a vector with their values, for each lag order.
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022). Poisson network autoregression. https://arxiv.org/abs/2104.06296
lin_estimnarpq, log_lin_ic_plot
data(crime)
data(crime_W)
lin_ic_plot(crime, crime_W, p = 1:3)
Starting values for the linear Poisson Network Autoregressive model of order
$p$
with $q$
covariates (PNAR($p$
)).
lin_narpq_init(y, W, p, Z = NULL)
y 
A 
W 
The 
p 
The number of lags in the model. 
Z 
An 
The function computes starting values to be used in the function lin_estimnarpq
.
These are simply the ordinary least squares estimators with a correction.
If any of the the resulting coefficients is negative they become equal to 0.001
A vector with the initial values.
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022). Poisson network autoregression. https://arxiv.org/abs/2104.06296
data(crime)
data(crime_W)
x0 < lin_narpq_init(crime, crime_W, p = 2)
Estimation of the loglinear Poisson Network Autoregressive model of order
$p$
with $q$
covariates (logPNAR($p$
)).
log_lin_estimnarpq(y, W, p, Z = NULL, uncons = FALSE, init = NULL,
xtol_rel = 1e8, maxeval = 100)
y 
A 
W 
The 
p 
The number of lags in the model. 
Z 
An 
uncons 
logical, if TRUE an unconstrained optimization is performed (default is FALSE). 
init 
A vector of starting values for the optimization algorithm. If this is NULL, the function computes them internally. 
xtol_rel 
The stopping tolerance of the optimization algorithm. 
maxeval 
The maximum number of evalutions the optimization algorithm will perform. 
This function performs a constrained estimation of the linear Poisson NAR($p$
) model with $q$
nonnegative valued covariates, for each node of the network $i=1,...,N$
over the time sample $t=1,...,TT$
, defined as
$\nu_{i,t}=\beta_{0}+\sum_{h=1}^{p}(\beta_{1h}X_{i,th}+\beta_{2h}Y_{i,th})+\sum_{l=1}^{q}\delta_{l}Z_{i,l},$
where $X_{i,t}=\sum_{j=1}^{N}W_{ij}Y_{j,t}$
is the network effect, i.e. the weighted average impact of node $i$
connections, with the weights of the mean being $W_{ij}$
, the single element of the network matrix $W$
. The sequence $\nu_{i,t}$
is the log of the expectation of $Y_{i,t}$
, conditional to its past values. The parameter $\beta_{0}$
is the intercept of the model, $\beta_{1h}$
are the network coefficients, $\beta_{2h}$
are the autoregressive parameters, and $\delta_{l}$
are the coefficients assocciated to the covariates $Z_{i,l}$
.
The estimation of the parameters of the model is performed by Quasi Maximum Likelihood Estimation (QMLE), maximizing the following quasi loglikelihood
$l(\theta)=\sum_{t=1}^{TT}\sum_{i=1}^{N}\left[Y_{i,t}\nu_{i,t}(\theta)e^{\nu_{i,t}(\theta)}\right]$
with respect to the vector of unknown parameters $\theta$
described above.
By default, the optimization is constrained in the stationary region where $\sum_{h=1}^{p}(\beta_{1h}+\beta_{2h})<1$
; this can be removed by setting uncons = TRUE
. However, the model estimates might be inconsistent if the estimated parameters lie outside the stationary region.
The ordinary least squares estimates are employed as starting values of the optimization procedure. Robust standard errors and ztests are also returned.
A list with attribute class "PNAR" including:
coefs 
A matrix with the estimated QMLE coefficients, their standard errors, their Ztest statistics and the relevant pvalues computed via the standard normal approximation. 
score 
The value of the quasi score function at the optimization point. It should be close to 0 if the optimization is successful. 
loglik 
The value of the maximized quasi loglikelihood. 
ic 
A vector with the Akaike information criterion (AIC), the Bayesian information criterion (BIC) and the Quasi information criterion (QIC). 
Alternatively, these can be printed via the function summary.PNAR
.
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Armillotta, M., Tsagris, M. and Fokianos, K. (2022c). The Rpackage PNAR for modelling count network time series. https://arxiv.org/abs/2211.02582
data(crime)
data(crime_W)
mod1 < log_lin_estimnarpq(crime, crime_W, p = 2)
summary(mod1)
Scatter plot of information criteria versus the number of lags in loglinear Poisson Network Autoregressive model of order
$p$
with $q$
covariates (logPNAR($p$
)).
log_lin_ic_plot(y, W, p = 1:10, Z = NULL, uncons = FALSE, ic = "QIC")
y 
A 
W 
The 
p 
A vector with integer numbers, the range of lags in the model, for which the AIC, BIC and QIC will be computed. 
Z 
An 
uncons 
Logical, if TRUE an unconstrained optimization without stationarity constraints is performed (default is FALSE). 
ic 
The information criterion you want to plot, "QIC" (default value), "AIC" or "BIC". 
The function computes the AIC, BIC or QIC for a range of lag orders of the
loglinear Poisson Network Autoregressive model of order $p$
with $q$
covariates (PNAR($p$
)).
A scatter plot with the lag order versus either QIC (default), AIC or BIC, and a vector with their values, for each lag order.
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022). Poisson network autoregression. https://arxiv.org/abs/2104.06296
log_lin_estimnarpq, lin_ic_plot
data(crime)
data(crime_W)
log_lin_ic_plot(crime, crime_W, p = 1:3)
Starting values for the loglinear Poisson Network Autoregressive model of order
$p$
with $q$
covariates (logPNAR($p$
)).
log_lin_narpq_init(y, W, p, Z = NULL)
y 
A 
W 
The 
p 
The number of lags in the model. 
Z 
An 
This function computes initial values for the loglinear Poisson Network
Autoregressive model of order $p$
with $q$
covariates (logPNAR($p$
))
with stationarity conditions. These initial values are simply the ordinary least
squares estimators with a correction.
A vector with the initial values.
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022). Poisson network autoregression. https://arxiv.org/abs/2104.06296
data(crime)
data(crime_W)
mod1 < log_lin_narpq_init(crime, crime_W, p = 2)
Generation of multivariate count time series from a linear Poisson Network Autoregressive
model of order $p$
with $q$
covariates (PNAR($p$
)).
poisson.MODpq(b, W, p, Z = NULL, TT, N, copula = "gaussian",
corrtype = "equicorrelation", rho, dof = 1)
b 
The coefficients of the model, in the following order: (intercept, network
parameters, autoregressive parameters, covariates). The dimension of the vector
should be 
W 
The 
p 
The number of lags in the model. 
Z 
An 
TT 
The temporal sample size. 
N 
The number of nodes on the network. 
copula 
Which copula function to use? The choices are "gaussian", "t", or "clayton". 
rho 
The value of the copula parameter ( 
corrtype 
Used only for elliptical copulas. The type of correlation matrix employed for
the copula; it will either be the "equicorrelation" or "toeplitz". The
"equicorrelation" option generates a correlation matrix where all the offdiagonal
entries equal 
dof 
The degrees of freedom for Student's t copula. 
This function generates counts from a linear Poisson NAR($p$
) model, where $q$
non timevarying
covariates are allowed as well. The counts are simulated from $Y_{t}=N_{t}(\lambda_{t})$
, where
$N_{t}$
is a sequence of $N$
dimensional IID Poisson count processes, with intensity 1, and
whose structure of dependence is modelled through a copula construction $C(\rho)$
on their associated exponential waiting times random variables. For details see Armillotta and Fokianos (2022, Sec. 2.12.2).
The sequence $\lambda_{i,t}$
is the expectation of $Y_{i,t}$
, conditional to its past values and it is generated by means of the following PNAR($p$
) model. For each node of the network $i=1,...,N$
over the time sample $t=1,...,TT$
$\lambda_{i,t}=\beta_{0}+\sum_{h=1}^{p}(\beta_{1h}X_{i,th}+\beta_{2h}Y_{i,th})+\sum_{l=1}^{q}\delta_{l}Z_{i,l}$
where $X_{i,t}=\sum_{j=1}^{N}W_{ij}Y_{j,t}$
is the network effect, i.e. the weighted average impact of node $i$
connections, with the weights of the mean being $W_{ij}$
, the single element of the network matrix $W$
. The parameter $\beta_{0}$
is the intercept of the model, $\beta_{1h}$
are the network coefficients, $\beta_{2h}$
are the autoregressive parameters, and $\delta_{l}$
are the coefficients assocciated to the covariates $Z_{i,l}$
.
A list including:
p2R 
The Toeplitz correlation matrix, if employed in the copula or NULL else. 
lambda 
A 
y 
A 
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Fokianos, K., Stove, B., Tjostheim, D., and P. Doukhan (2020). Multivariate count autoregression. Bernoulli, 26(1), 471499.
poisson.MODpq.log, poisson.MODpq.nonlin,
poisson.MODpq.stnar, poisson.MODpq.tnar
W < adja( N = 20, K = 5, alpha= 0.5)
y < poisson.MODpq( b = c(0.5, 0.3, 0.2), W = W, p = 1, Z = NULL,
TT = 1000, N = 20, copula = "gaussian",
corrtype = "equicorrelation", rho = 0.5)$y
Generation of counts from a loglinear Poisson Network Autoregressive model of
order $p$
with $q$
covariates (logPNAR($p$
)).
poisson.MODpq.log(b, W, p, Z = NULL, TT, N, copula = "gaussian",
corrtype = "equicorrelation", rho, dof = 1)
b 
The coefficients of the model, in the following order: (intercept, network
parameters, autoregressive parameters, covariates). The dimension of the vector
should be 
W 
The 
p 
The number of lags in the model. 
Z 
An 
TT 
The temporal sample size. 
N 
The number of nodes on the network. 
copula 
Which copula function to use? The "gaussian", "t", or "clayton". 
rho 
The the value of the copula parameter ( 
corrtype 
Used only for elliptical copulas. The type of correlation matrix employed for
the copula; it will either be the "equicorrelation" or "toeplitz". The
"equicorrelation" option generates a correlation matrix where all the offdiagonal
entries equal 
dof 
The degrees of freedom for Student's t copula. 
This function generates counts from a loglinear Poisson NAR($p$
) model, where $q$
non timevarying covariates are allowed as well. The counts are simulated from $Y_{t}=N_{t}(e^{\nu_{t}})$
, where $N_{t}$
is a sequence of $N$
dimensional IID Poisson count processes, with intensity 1, and whose structure of dependence is modelled through a copula construction $C(\rho)$
on their associated exponential waiting times random variables. For details see Armillotta and Fokianos (2022, Sec. 2.12.2).
The sequence $\nu_{t}$
is the log of the expecation of $Y_{t}$
, conditional to its past values and it is generated by means of the following logPNAR($p$
) model. For each node of the network $i=1,...,N$
over the time sample $t=1,...,TT$
$\nu_{i,t}=\beta_{0}+\sum_{h=1}^{p}(\beta_{1h}X_{i,th}+\beta_{2h}Y_{i,th})+\sum_{l=1}^{q}\delta_{l}Z_{i,l}$
where $X_{i,t}=\sum_{j=1}^{N}W_{ij}Y_{j,t}$
is the network effect, i.e. the weighted average impact of node $i$
connections, with the weights of the mean being $W_{ij}$
, the single element of the network matrix $W$
. The parameter $\beta_{0}$
is the intercept of the model, $\beta_{1h}$
are the network coefficients, $\beta_{2h}$
are the autoregressive parameters, and $\delta_{l}$
are the coefficients assocciated to the covariates $Z_{i,l}$
.
A list including:
p2R 
The Toeplitz correlation matrix, if employed in the copula or NULL else. 
log_lambda 
A 
y 
A 
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Fokianos, K., Stove, B., Tjostheim, D., and P. Doukhan (2020). Multivariate count autoregression. Bernoulli, 26(1), 471499.
poisson.MODpq, poisson.MODpq.nonlin,
poisson.MODpq.stnar, poisson.MODpq.tnar
W < adja( N = 20, K = 5, alpha= 0.5)
y < poisson.MODpq.log( b = c(0.5, 0.3, 0.2), W = W, p = 1,
Z = NULL, TT = 1000, N = 20, copula = "gaussian",
corrtype = "equicorrelation", rho = 0.5)$y
Generation of counts from a nonlinear Intercept Drift Poisson Network
Autoregressive model of order $p$
with $q$
covariates (IDPNAR($p$
)).
poisson.MODpq.nonlin(b, W, gama, p, d, Z = NULL, TT, N, copula = "gaussian",
corrtype = "equicorrelation", rho, dof = 1)
b 
The linear coefficients of the model, in the following order: (intercept,
network parameters, autoregressive parameters, covariates). The dimension of
the vector should be 
W 
The 
gama 
A scalar nonlinear intercept drift parameter. 
p 
The number of lags in the model. 
d 
The lag parameter of nonlinear variable (should be between 1 and 
Z 
An 
TT 
The temporal sample size. 
N 
The number of nodes on the network. 
copula 
Which copula function to use? The "gaussian", "t", or "clayton". 
rho 
The value of the copula parameter ( 
corrtype 
Used only for elliptical copulas. The type of correlation matrix employed for
the copula; it will either be the "equicorrelation" or "toeplitz". The
"equicorrelation" option generates a correlation matrix where all the offdiagonal
entries equal 
dof 
The degrees of freedom for Student's t copula. 
This function generates counts from a nonlinear Intercept Drift Poisson NAR($p$
) model, where $q$
non timevarying covariates are allowed as well. The counts are simulated from $Y_{t}=N_{t}(\lambda_{t})$
, where $N_{t}$
is a sequence of $N$
dimensional IID Poisson count processes, with intensity 1, and whose structure of dependence is modelled through a copula construction $C(\rho)$
on their associated exponential waiting times random variables. For details see Armillotta and Fokianos (2022a, Sec. 2.12.2). The sequence $\lambda_{i,t}$
is the expecation of $Y_{i,t}$
, conditional to its past values and it is generated by means of the following IDPNAR($p$
) model. For each node of the network $i=1,...,N$
over the time sample $t=1,...,TT$
$\lambda_{i,t}=\frac{\beta_{0}}{(1+X_{i,td})^{\gamma}}+\sum_{h=1}^{p}(\beta_{1h}X_{i,th}+\beta_{2h}Y_{i,th})+\sum_{l=1}^{q}\delta_{l}Z_{i,l}$
where $X_{i,t}=\sum_{j=1}^{N}W_{ij}Y_{j,t}$
is the network effect, i.e. the weighted average impact of node $i$
connections, with the weights of the mean being $W_{ij}$
, the single element of the network matrix $W$
.
The parameter $\beta_{0}$
is the intercept of the model, $\beta_{1h}$
are the network coefficients, $\beta_{2h}$
are the autoregressive parameters, $\gamma$
is the nonlinear coefficient associated with the intercept drift, and $\delta_{l}$
are the coefficients assocciated with the covariates $Z_{i,l}$
. The coefficient $d$
is considered as an extra parameter defining the lag of the network effect in the nonlinear part of the model and is left to be set by the user. For details on IDPNAR models see Armillotta and Fokianos (2022b, Sec. 2).
A list including:
p2R 
The Toeplitz correlation matrix, if employed in the copula or NULL else. 
lambda 
A 
y 
A 
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Armillotta, M. and K. Fokianos (2022b). Testing linearity for network autoregressive models. https://arxiv.org/abs/2202.03852
poisson.MODpq, poisson.MODpq.log,
poisson.MODpq.stnar, poisson.MODpq.tnar
W < adja( N = 20, K = 5, alpha= 0.5)
y < poisson.MODpq.nonlin( b = c(0.5, 0.3, 0.2), W = W, gama = 1, p = 1,
d = 1, Z = NULL, TT = 1000, N = 20, copula = "gaussian",
corrtype = "equicorrelation", rho = 0.5)$y
Generation of multivariate count time series from a nonlinear Smooth Transition Poisson Network
Autoregressive model of order $p$
with $q$
covariates (STPNAR($p$
)).
poisson.MODpq.stnar(b, W, gama, a, p, d, Z = NULL, TT, N, copula = "gaussian",
corrtype = "equicorrelation", rho, dof = 1)
b 
The linear coefficients of the model, in the following order: (intercept, network
parameters, autoregressive parameters, covariates). The dimension of the vector
should be 
W 
The 
gama 
The scalar nuisance smoothing parameter. 
a 
Vector of nonlinear parameters. The dimension of the vector should be 
p 
The number of lags in the model. 
d 
The lag parameter of nonlinear variable (should be between 1 and 
Z 
An 
TT 
The temporal sample size. 
N 
The number of nodes on the network. 
copula 
Which copula function to use? The choices are "gaussian", "t", or "clayton". 
rho 
The value of the copula parameter ( 
corrtype 
Used only for elliptical copulas. The type of correlation matrix employed for
the copula; it will either be the "equicorrelation" or "toeplitz". The
"equicorrelation" option generates a correlation matrix where all the offdiagonal
entries equal 
dof 
The degrees of freedom for Student's t copula. 
This function generates counts from a nonlinear Smooth Transition Poisson NAR($p$
) model, where
$q$
non timevarying covariates are allowed as well. The counts are simulated from $Y_{t}=N_{t}(\lambda_{t})$
, where $N_{t}$
is a sequence of $N$
dimensional IID Poisson count processes, with intensity 1, and whose structure of dependence is modelled through a copula construction $C(\rho)$
on their associated exponential waiting times random variables. For details see Armillotta and Fokianos (2022a, Sec. 2.12.2).
The sequence $\lambda_{i,t}$
is the expecation of $Y_{i,t}$
, conditional to its past values and it is generated by means of the following STPNAR($p$
) model. For each node of the network $i=1,...,N$
over the time sample $t=1,...,TT$
$\lambda_{i,t}=\beta_{0}+\sum_{h=1}^{p}(\beta_{1h}X_{i,th}+\beta_{2h}Y_{i,th}+\alpha_{h}e^{\gamma X_{i,td}^{2}}X_{i,th})+\sum_{l=1}^{q}\delta_{l}Z_{i,l}$
where $X_{i,t}=\sum_{j=1}^{N}W_{ij}Y_{j,t}$
is the network effect, i.e. the weighted average impact of node $i$
connections, with the weights of the mean being $W_{ij}$
, the single element of the network matrix $W$
.
The parameter $\beta_{0}$
is the intercept of the model, $\beta_{1h}$
are the network coefficients, $\beta_{2h}$
are the autoregressive parameters, $\alpha_{h}$
are the nonlinear smooth transition parameters, $\gamma$
is the nuisance smoothing parameter, and $\delta_{l}$
are the coefficients assocciated to the covariates $Z_{i,l}$
. The coefficient $d$
is considered as an extra parameter defining the lag of the network effect in the nonlinear part of the model and is left to be set by the user. For details on STPNAR models see Armillotta and Fokianos (2022b, Sec. 2).
A list including:
p2R 
The Toeplitz correlation matrix, if employed in the copula or NULL else. 
lambda 
A 
y 
A 
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Armillotta, M. and K. Fokianos (2022b). Testing linearity for network autoregressive models. https://arxiv.org/abs/2202.03852
poisson.MODpq, poisson.MODpq.log,
poisson.MODpq.nonlin, poisson.MODpq.tnar
W < adja( N = 20, K = 5, alpha= 0.5)
y < poisson.MODpq.stnar( b = c(0.5, 0.3, 0.2), W = W, gama = 0.2, a = 0.4,
p = 1, d = 1, Z = NULL, TT = 1000, N = 20, copula = "gaussian",
corrtype = "equicorrelation", rho = 0.5)$y
Generation of multivariate count time series from a nonlinear Threshold Poisson network Autoregressive
model of order $p$
with $q$
covariates (TPNAR($p$
)).
poisson.MODpq.tnar(b, W, gama, a, p, d, Z = NULL, TT, N, copula = "gaussian",
corrtype = "equicorrelation", rho, dof = 1)
b 
The linear coefficients of the model, in the following order: (intercept, network
parameters, autoregressive parameters, covariates). The dimension of the vector
should be 
W 
The 
gama 
The scalar nuisance threshold parameter. 
a 
Vector of nonlinear parameters. The dimension of the vector should be

p 
The number of lags in the model. 
d 
The lag parameter of nonlinear variable (should be between 1 and 
Z 
An 
TT 
The temporal sample size. 
N 
The number of nodes on the network. 
copula 
Which copula function to use? The "gaussian", "t", or "clayton". 
rho 
The value of the copula parameter ( 
corrtype 
Used only for elliptical copulas. The type of correlation matrix employed for
the copula; it will either be the "equicorrelation" or "toeplitz". The
"equicorrelation" option generates a correlation matrix where all the offdiagonal
entries equal 
dof 
The degrees of freedom for Student's t copula. 
This function generates counts from a nonlinear Threshold Poisson NAR($p$
) model, where $q$
non timevarying covariates are allowed as well. The counts are simulated from $Y_{t}=N_{t}(\lambda_{t})$
, where $N_{t}$
is a sequence of $N$
dimensional IID Poisson count processes, with intensity 1, and whose structure of dependence is modelled through a copula construction $C(\rho)$
on their associated exponential waiting times random variables. For details see Armillotta and Fokianos (2022a, Sec. 2.12.2).
The sequence $\lambda_{i,t}$
is the expecation of $Y_{i,t}$
, conditional to its past values and it is generated by means of the following TPNAR($p$
) model. For each node of the network $i=1,...,N$
over the time sample $t=1,...,TT$
$\lambda_{i,t}=\beta_{0}+\sum_{h=1}^{p}\left[\beta_{1h}X_{i,th}+\beta_{2h}Y_{i,th}+(\alpha_{0}+\alpha_{1h}X_{i,th}+\alpha_{2h}Y_{i,th})I(X_{i,td}\leq\gamma)\right]+\sum_{l=1}^{q}\delta_{l}Z_{i,l}$
where $X_{i,t}=\sum_{j=1}^{N}W_{ij}Y_{j,t}$
is the network effect, i.e. the weighted average impact of node $i$
connections, with the weights of the mean being $W_{ij}$
, the single element of the network matrix $W$
, and $I()$
is the indicator function.
The parameter $\beta_{0}$
is the intercept of the model, $\beta_{1h}$
are the network coefficients, $\beta_{2h}$
are the autoregressive parameters, the $\alpha$
vector of nonlinear parameters is divided as follows: $\alpha_{0}$
is the intercept, $\alpha_{1h}$
are the network coefficients, $\alpha_{2h}$
are the autoregressive parameters; $\gamma$
is the nuisance threshold parameter, and $\delta_{l}$
are the coefficients assocciated to the covariates $Z_{i,l}$
. The coefficient $d$
is considered as an extra parameter defining the lag of the network effect in the nonlinear part of the model and is left to be set by the user. For details on TPNAR models see Armillotta and Fokianos (2022b, Sec. 2).
A list including:
p2R 
The Toeplitz correlation matrix, if employed in the copula or NULL else. 
lambda 
A 
y 
A 
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Armillotta, M. and K. Fokianos (2022b). Testing linearity for network autoregressive models. https://arxiv.org/abs/2202.03852
poisson.MODpq, poisson.MODpq.log,
poisson.MODpq.nonlin, poisson.MODpq.stnar
W < adja( N = 20, K = 5, alpha= 0.5)
y < poisson.MODpq.tnar( b = c(0.5, 0.3, 0.2), W = W, gama = 1,
a = c(0.2, 0.2, 0.2), p = 1, d = 1, Z = NULL, TT = 1000, N = 20,
copula = "gaussian", corrtype = "equicorrelation", rho = 0.5)$y
Random number generation of copula functions.
rcopula(n, N, copula = "gaussian", corrtype = "equicorrelation",
rho, dof = 1, cholR = NULL)
n 
The number of random values to generate. 
N 
The number of variables for which random valeus will be generated. 
copula 
Which copula function to use? The "gaussian", "t", or "clayton". 
rho 
The the value of the copula parameter ( 
corrtype 
Used only for elliptical copulas. The type of correlation matrix employed for
the copula; it will either be the "equicorrelation" or "toeplitz". The
"equicorrelation" option generates a correlation matrix where all the offdiagonal
entries equal 
dof 
The degrees of freedom for Student's t copula. 
cholR 
An alternative input for elliptic copulas, providing directly the Cholesky decomposition for a specific correlation matrix to be passed, otherwise leave it NULL. 
This function generates random copula values from Gaussian, Student's t, or Clayton copulas based on a single copula paremeter and different correlation structures.
An $n \times N$
matrix with the simulated copula values.
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Nelsen, Roger B. (1999). An Introduction to Copulas, Springer.
getN, poisson.MODpq, poisson.MODpq.log
u < rcopula(n = 100, N = 50, rho = 0.3)
Quasi score test for testing linearity of Poisson Network Autoregressive model
of order $p$
against the nonlinear Intercep Drift (ID) version
(IDPNAR($p$
)).
score_test_nonlinpq_h0(b, y, W, p, d, Z = NULL)
b 
The estimated parameters from the linear PNAR model, in the following order:
(intercept, network parameters, autoregressive parameters, covariates).
The dimension of the vector should be 
y 
A 
W 
The 
p 
The number of lags in the model. 
d 
The lag parameter of nonlinear variable (should be between 1 and 
Z 
An 
The function computes the quasi score test for testing linearity of Poisson Network Autoregressive model of order $p$
against the following IDPNAR($p$
) model. For each node of the network $i=1,...,N$
over the time sample $t=1,...,TT$
$\lambda_{i,t}=\frac{\beta_{0}}{(1+X_{i,td})^{\gamma}}+\sum_{h=1}^{p}(\beta_{1h}X_{i,th}+\beta_{2h}Y_{i,th})+\sum_{l=1}^{q}\delta_{l}Z_{i,l}$
where $X_{i,t}=\sum_{j=1}^{N}W_{ij}Y_{j,t}$
is the network effect, i.e. the weighted average impact of node $i$
connections, with the weights of the mean being $W_{ij}$
, the single element of the network matrix $W$
. The sequence $\lambda_{i,t}$
is the expectation of $Y_{i,t}$
conditional to its past values.
The null hypothesis of the test is defined as $H_{0}: \gamma=0$
, versus the alternative $H_{1}: \gamma >0$
. The test statistic has the form
$LM=S^{'}(\hat{\theta})\Sigma^{1}(\hat{\theta})S(\hat{\theta}),$
where
$S(\hat{\theta})=\sum_{t=1}^{TT}\sum_{i=1}^{N}\left(\frac{Y_{i,t}}{\lambda_{i,t}(\hat{\theta})}1\right)\frac{\partial\lambda_{i,t}(\hat{\theta})}{\partial\gamma}$
is the partition of the quasi score related to the nonlinear parameter $\gamma$
, evaluated at the estimated parameters $\hat{\theta}$
under the null assumption $H_{0}$
(linear model), and $\Sigma(\hat{\theta})$
is the variance of $S(\hat{\theta})$
. Under $H_{0}$
, the test asymptotically follows the $\chi^2$
distribution with 1 degree of freedom. For details see Armillotta and Fokianos (2022b, Sec. 4).
A list with attribute class "htest" including:
statistic 
The value of the 
parameter 
The degrees of freedom of the 
p.value 
The pvalue of the 
null.value 
The value of the 
alternative 
The alternative hypothesis, 
method 
The name of the test. 
data.name 
Information on the arguments used. 
Alternatively, these can be printed via the function summary.nonlin
.
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Armillotta, M. and K. Fokianos (2022b). Testing linearity for network autoregressive models. https://arxiv.org/abs/2202.03852
Armillotta, M., Tsagris, M. and Fokianos, K. (2022c). The Rpackage PNAR for modelling count network time series. https://arxiv.org/abs/2211.02582
score_test_stnarpq_j, score_test_tnarpq_j,
lin_estimnarpq
data(crime)
data(crime_W)
mod1 < lin_estimnarpq(crime, crime_W, p = 2)
ca < mod1$coefs[, 1]
score_test_nonlinpq_h0(ca, crime, crime_W, p = 2, d = 1)
Computation of Davies bound pvalue for the suptype test for testing linearity
of Poisson Network Autoregressive model of order $p$
(PNAR($p$
)) versus
the nonlinear Smooth Transition alternative (STPNAR($p$
)).
score_test_stnarpq_DV(b, y, W, p, d, Z = NULL, gama_L = NULL,
gama_U = NULL, len = 100)
b 
The estimated parameters from the linear model, in the following order:
(intercept, network parameters, autoregressive parameters, covariates).
The dimension of the vector should be 
y 
A 
W 
The 
p 
The number of lags in the model. 
d 
The lag parameter of nonlinear variable (should be between 1 and 
Z 
An 
gama_L 
The lower value of the nuisance parameter 
gama_U 
The upper value of the nuisance parameter 
len 
The length of the grid of values of 
The function computes an upperbound for the pvalue of the suptype test for testing linearity of Poisson Network
Autoregressive model of order $p$
(PNAR($p$
)) versus the following Smooth Transition alternative (STPNAR($p$
)).
For each node of the network $i=1,...,N$
over the time sample $t=1,...,TT$
$\lambda_{i,t}=\beta_{0}+\sum_{h=1}^{p}(\beta_{1h}X_{i,th}+\beta_{2h}Y_{i,th}+\alpha_{h}e^{\gamma X_{i,td}^{2}}X_{i,th})+\sum_{l=1}^{q}\delta_{l}Z_{i,l}$
where $X_{i,t}=\sum_{j=1}^{N}W_{ij}Y_{j,t}$
is the network effect, i.e. the weighted average impact of node $i$
connections, with the weights of the mean being $W_{ij}$
, the single element of the network matrix $W$
. The sequence $\lambda_{i,t}$
is the expectation of $Y_{i,t}$
, conditional to its past values.
The null hypothesis of the test is defined as $H_{0}: \alpha_{1}=...=\alpha_{p}=0$
, versus the alternative that at least one among $\alpha_{h}$
is not 0. The test statistic has the form
$LM(\gamma)=S^{'}(\hat{\theta},\gamma)\Sigma^{1}(\hat{\theta},\gamma)S(\hat{\theta},\gamma),$
where
$S(\hat{\theta},\gamma)=\sum_{t=1}^{TT}\sum_{i=1}^{N}\left(\frac{Y_{i,t}}{\lambda_{i,t}(\hat{\theta},\gamma)}1\right)\frac{\partial\lambda_{i,t}(\hat{\theta},\gamma)}{\partial\alpha}$
is the partition of the quasi score related to the vector of nonlinear parameters $\alpha=(\alpha_{1},...,\alpha_{p})$
, evaluated at the estimated parameters $\hat{\theta}$
under the null assumption $H_{0}$
(linear model), and $\Sigma(\hat{\theta},\gamma)$
is the variance of $S(\hat{\theta},\gamma)$
. Since the test statistic depends on an unknown nuisance parameter ($\gamma$
), the supremum of the statistic is considered in the test, $\sup_{\gamma}LM(\gamma)$
. The function computes the bound of the pvalue, suggested by Davies (1987), for the test statistic $\sup_{\gamma}LM(\gamma)$
, with scalar nuisance parameter $\gamma$
, as follows.
$P(\chi^{2}_{k} \geq M)+V M^{1/2(k1)}\frac{e^{M/2}2^{k/2}}{\Gamma(k/2)}$
where $M$
is the maximum of the test statistic $LM(\gamma)$
, computed by the available sample, over a grid of values for the nuisance parameter $\gamma_{F}=(\gamma_{L},\gamma_{1},...,\gamma_{l},\gamma_{U})$
; $k$
is the number of nonlinear parameters tested. So the first summand of the bound is just the pvalue of a chisquare test with $k$
degrees of freedom. The second summand is a correction term depending on $V$
, which is the approximated total variation computed as
$V=LM^{1/2}(\gamma_{1})LM^{1/2}(\gamma_{L})+LM^{1/2}(\gamma_{2})LM^{1/2}(\gamma_{1})+...+LM^{1/2}(\gamma_{U})LM^{1/2}(\gamma_{l}).$
The feasible bound allows to approximate the pvalues of the suptype test in a straightforward way, by adding to the tail probability of a chisquare distribution a correction term which depends on the total variation of the process. For details see Armillotta and Fokianos (2022b, Sec. 5).
The values of gama_L
and gama_U
are computed internally as gama_L
$=\log(0.9)/X^{2}$
and gama_U
$=\log(0.1)/X^{2}$
, where $X$
is the overall mean of $X_{i,t}$
over the nodes $i=1,...,N$
and times $t=1,...,TT$
. Since the nonlinear function $e^{\gamma X_{i,td}^{2}}$
ranges between 0 and 1, by considering $X$
to be a representative value for the network mean, gama_U
and gama_L
would be the values of $\gamma$
leading the nonlinear switching function to be 0.1 and 0.9, respectively, so that in the optimization procedure the extremes of the function domain are excluded. Alternatively, their values can be supplied by the user.
A list including:
DV 
The Davies bound of pvalues for sup test. 
supLM 
The value of the sup test statistic in the sample 
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Armillotta, M. and K. Fokianos (2022b). Testing linearity for network autoregressive models. https://arxiv.org/abs/2202.03852
Davies, R. B. (1987). Hypothesis testing when a nuisance parameter is present only under the alternative. Biometrika 74, 3343.
Armillotta, M., Tsagris, M. and Fokianos, K. (2022c). The Rpackage PNAR for modelling count network time series. https://arxiv.org/abs/2211.02582
score_test_stnarpq_j, global_optimise_LM_stnarpq
data(crime)
data(crime_W)
mod1 < lin_estimnarpq(crime, crime_W, p = 1)
ca < mod1$coefs[, 1]
score_test_stnarpq_DV(ca, crime, crime_W, p = 1, d = 1)
Computation of bootstrap pvalue for the suptype test for testing linearity of
Poisson Network Autoregressive model of order $p$
(PNAR($p$
)) versus the
nonlinear Smooth Transition alternative (STPNAR($p$
)).
score_test_stnarpq_j(supLM, b, y, W, p, d, Z = NULL, J = 499,
gama_L = NULL, gama_U = NULL, tol = 1e9, ncores = 1, seed = NULL)
supLM 
The optimized value of the test statistic. See the function

b 
The estimated parameters from the linear model, in the following order:
(intercept, network parameters, autoregressive parameters, covariates).
The dimension of the vector should be 
y 
A 
W 
The 
p 
The number of lags in the model. 
d 
The lag parameter of nonlinear variable (should be between 1 and 
Z 
An 
J 
The number of bootstrap samples to draw. 
gama_L 
The lower value of the nuisance parameter 
gama_U 
The upper value of the nuisance parameter 
tol 
Tolerance level for the optimizer. 
ncores 
Number of cores to use for parallel computing. By default the number of cores is set to 1 (no parallel computing). Note: If for some reason the parallel does not work then load the doParallel package yourseleves. 
seed 
To replicate the results use a seed for the generator, an integer number. 
The function computes a bootstrap pvalue for the suptype test for testing linearity of Poisson Network Autoregressive model of order $p$
(PNAR($p$
)) versus the following Smooth Transition alternative (STPNAR($p$
)). For each node of the network $i=1,...,N$
over the time sample $t=1,...,TT$
$\lambda_{i,t}=\beta_{0}+\sum_{h=1}^{p}(\beta_{1h}X_{i,th}+\beta_{2h}Y_{i,th}+\alpha_{h}e^{\gamma X_{i,td}^{2}}X_{i,th})+\sum_{l=1}^{q}\delta_{l}Z_{i,l},$
where $X_{i,t}=\sum_{j=1}^{N}W_{ij}Y_{j,t}$
is the network effect, i.e. the weighted average impact of node $i$
connections, with the weights of the mean being $W_{ij}$
, the single element of the network matrix $W$
. The sequence $\lambda_{i,t}$
is the expectation of $Y_{i,t}$
, conditional to its past values.
The null hypothesis of the test is defined as $H_{0}: \alpha_{1}=...=\alpha_{p}=0$
, versus the alternative that at least one among $\alpha_{h}$
is not $0$
. The test statistic has the form
$LM(\gamma)=S^{'}(\hat{\theta},\gamma)\Sigma^{1}(\hat{\theta},\gamma)S(\hat{\theta},\gamma)$
where
$S(\hat{\theta},\gamma)=\sum_{t=1}^{TT}\sum_{i=1}^{N}\left(\frac{Y_{i,t}}{\lambda_{i,t}(\hat{\theta},\gamma)}1\right)\frac{\partial\lambda_{i,t}(\hat{\theta},\gamma)}{\partial\alpha}$
is the partition of the quasi score related to the vector of nonlinear parameters $\alpha=(\alpha_{1},...,\alpha_{p})$
, evaluated at the estimated parameters $\hat{\theta}$
under the null assumption $H_{0}$
(linear model), and $\Sigma(\hat{\theta},\gamma)$
is the variance of $S(\hat{\theta},\gamma)$
.
Since the test statistic depends on an unknown nuisance parameter ($\gamma$
), the supremum of the statistic is considered in the test, $\sup_{\gamma}LM(\gamma)$
. This value can be computed for the available sample by using the function global_optimise_LM_stnarpq
and should be supplied here as an input supLM
.
The function performs the bootstrap resampling of the test statistic $\sup_{\gamma}LM(\gamma)$
by employing Gaussian perturbations of the score $S(\hat{\theta},\gamma)$
. For details see Armillotta and Fokianos (2022b, Sec. 5).
The values of gama_L
and gama_U
are computed internally as gama_L
$=\log(0.9)/X^{2}$
and gama_U
$=\log(0.1)/X^{2}$
, where $X$
is the overall mean of $X_{i,t}$
over the nodes $i=1,...,N$
and times $t=1,...,TT$
. Since the nonlinear function $e^{\gamma X_{i,td}^{2}}$
ranges between 0 and 1, by considering $X$
to be a representative value for the network mean, gama_U
and gama_L
would be the values of $\gamma$
leading the nonlinear switching function to be 0.1 and 0.9, respectively, so that in the optimization procedure the extremes of the function domain are excluded. Alternatively, their value can be supplied by the user.
Note: For large datasets the function may require few minutes to run. Parallel computing is suggested to speed up the computations.
A list including:
pJ 
The bootstrap pvalue of the sup test. 
cpJ 
The adjusted version of bootstrap pvalue of the sup test. 
gamaj 
The optimal values of the 
supLMj 
The values of perturbed test statistic at the optimum point 
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Armillotta, M. and K. Fokianos (2022b). Testing linearity for network autoregressive models. https://arxiv.org/abs/2202.03852
Armillotta, M., Tsagris, M. and Fokianos, K. (2022c). The Rpackage PNAR for modelling count network time series. https://arxiv.org/abs/2211.02582
score_test_stnarpq_DV, global_optimise_LM_stnarpq,
score_test_tnarpq_j
# load data
data(crime)
data(crime_W)
#estimate linear PNAR model
mod1 < lin_estimnarpq(crime, crime_W, p = 2)
b < mod1$coefs[, 1]
g < global_optimise_LM_stnarpq(b = b, y = crime, W = crime_W, p = 2, d = 1)
supg < g$supLM
score_test_stnarpq_j(supLM = supg, b = b, y = crime, W = crime_W, p = 2, d = 1, J = 5)
Computation of bootstrap pvalue for the suptype test for testing linearity of
Poisson Network Autoregressive model of order $p$
(PNAR($p$
)) versus the
nonlinear Threshold alternative (TPNAR($p$
)).
score_test_tnarpq_j(supLM, b, y, W, p, d, Z = NULL, J = 499,
gama_L = NULL, gama_U = NULL, tol = 1e9, ncores = 1, seed = NULL)
supLM 
The optimized value of the test statistic. See the function

b 
The estimated parameters from the linear model, in the following order:
(intercept, network parameters, autoregressive parameters, covariates).
The dimension of the vector should be 
y 
A 
W 
The 
p 
The number of lags in the model. 
d 
The lag parameter of nonlinear variable (should be between 1 and 
Z 
An 
J 
The number of bootstrap samples to draw. 
gama_L 
The lower value of the nuisance parameter 
gama_U 
The upper value of the nuisance parameter 
tol 
Tolerance level for the optimizer. 
ncores 
Number of cores to use for parallel computing. By default the number of cores is set to 1 (no parallel computing). Note: If for some reason the parallel does not work then load the doParallel package yourseleves. 
seed 
To replicate the results use a seed for the generator, an integer number. 
The function computes a bootstrap pvalue for the suptype test for testing linearity of Poisson Network Autoregressive model of order $p$
(PNAR($p$
)) versus the following Threshold alternative (TPNAR($p$
)). For each node of the network $i=1,...,N$
over the time sample $t=1,...,TT$
$\lambda_{i,t}=\beta_{0}+\sum_{h=1}^{p}\left[\beta_{1h}X_{i,th}+\beta_{2h}Y_{i,th}+(\alpha_{0}+\alpha_{1h}X_{i,th}+\alpha_{2h}Y_{i,th})I(X_{i,td}\leq\gamma)\right]+\sum_{l=1}^{q}\delta_{l}Z_{i,l}$
where $X_{i,t}=\sum_{j=1}^{N}W_{ij}Y_{j,t}$
is the network effect, i.e. the weighted average impact of node $i$
connections, with the weights of the mean being $W_{ij}$
, the single element of the network matrix $W$
. The sequence $\lambda_{i,t}$
is the expectation of $Y_{i,t}$
, conditional to its past values.
The null hypothesis of the test is defined as $H_{0}: \alpha_{0}=\alpha_{11}=...=\alpha_{2p}=0$
, versus the alternative that at least one among $\alpha_{s,h}$
is not $0$
, for $s=0,1,2$
. The test statistic has the form
$LM(\gamma)=S^{'}(\hat{\theta},\gamma)\Sigma^{1}(\hat{\theta},\gamma)S(\hat{\theta},\gamma)$
where
$S(\hat{\theta}, \gamma)=\sum_{t=1}^{TT}\sum_{i=1}^{N}\left(\frac{Y_{i,t}}{\lambda_{i,t}(\hat{\theta}, \gamma)}1\right)\frac{\partial\lambda_{i,t}(\hat{\theta},\gamma)}{\partial\alpha}$
is the partition of the quasi score related to the vector of nonlinear parameters $\alpha=(\alpha_{0},...,\alpha_{2p})$
, evaluated at the estimated parameters $\hat{\theta}$
under the null assumption $H_{0}$
(linear model), and $\Sigma(\hat{\theta},\gamma)$
is the variance of $S(\hat{\theta},\gamma)$
.
Since the test statistic depends on an unknown nuisance parameter ($\gamma$
), the supremum of the statistic is considered in the test, $\sup_{\gamma}LM(\gamma)$
. This value can be computed for the available sample by using the function global_optimise_LM_tnarpq
and should be supplied here as an input supLM
.
The function performs the bootstrap resampling of the test statistic $\sup_{\gamma}LM(\gamma)$
by employing Gaussian perturbations of the score $S(\hat{\theta},\gamma)$
. For details see Armillotta and Fokianos (2022b, Sec. 5).
The values of gama_L
and gama_U
are computed internally 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,...,TT$
. In this way the optimization is performed for values of $\gamma$
such that the indicator function $I(X_{i,td}\leq\gamma)$
is not always close to 0 or 1. Alternatively, their value can be supplied by the user. For details see Armillotta and Fokianos (2022b, Sec. 45).
Note: For large datasets the function may require few minutes to run. Parallel computing is suggested to speed up the computations.
A list including:
pJ 
The bootstrap pvalue of the sup test. 
cpJ 
The adjusted version of bootstrap pvalue of the sup test. 
gamaj 
The optimal values of the 
supLMj 
The values of perturbed test statistic at the optimum point 
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Armillotta, M. and K. Fokianos (2022b). Testing linearity for network autoregressive models. https://arxiv.org/abs/2202.03852
Armillotta, M., Tsagris, M. and Fokianos, K. (2022c). The Rpackage PNAR for modelling count network time series. https://arxiv.org/abs/2211.02582
global_optimise_LM_tnarpq,
global_optimise_LM_stnarpq, score_test_stnarpq_j
# load data
data(crime)
data(crime_W)
#estimate linear PNAR model
mod1 < lin_estimnarpq(crime, crime_W, p = 2)
b < mod1$coefs[, 1]
g < global_optimise_LM_tnarpq(b = b, y = crime, W = crime_W, p = 2, d = 1)
supg < g$supLM
score_test_tnarpq_j(supLM = supg, b = b, y = crime, W = crime_W, p = 2, d = 1, J = 5)
S3 methods for extracting the results of the bound pvalue for testing for smooth transition effects on PNAR($p$
) model.
## S3 method for class 'DV'
summary(object, ...)
## S3 method for class 'summary.DV'
print(x, ...)
## S3 method for class 'DV'
print(x, ...)
object 
An object containing the results of the function 
x 
An object containing the results of the function 
... 
Extra arguments the user can pass. 
The functions print the output of the bound pvalue for testing for smooth transition effects on PNAR($p$
) model.
The functions print the results of the function score_test_stnarpq_DV
.
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Armillotta, M. and K. Fokianos (2022b). Testing linearity for network autoregressive models. https://arxiv.org/abs/2202.03852
Davies, R. B. (1987). Hypothesis testing when a nuisance parameter is present only under the alternative. Biometrika 74, 3343.
Armillotta, M., Tsagris, M. and Fokianos, K. (2022c). The Rpackage PNAR for modelling count network time series. https://arxiv.org/abs/2211.02582
data(crime)
data(crime_W)
mod1 < lin_estimnarpq(crime, crime_W, p = 1)
ca < mod1$coefs[, 1]
a < score_test_stnarpq_DV(ca, crime, crime_W, p = 1, d = 1)
print(a)
summary(a)
S3 methods for extracting the results of the nonlinear hypothesis test.
## S3 method for class 'nonlin'
summary(object, ...)
## S3 method for class 'summary.nonlin'
print(x, ...)
## S3 method for class 'nonlin'
print(x, ...)
object 
An object containing the results of the function 
x 
An object containing the results of the function 
... 
Extra arguments the user can pass. 
The functions print the output of the nonlinear hypothesis test.
The functions print the results of the function score_test_nonlinpq_h0
.
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Armillotta, M. and K. Fokianos (2022b). Testing linearity for network autoregressive models. https://arxiv.org/abs/2202.03852
Armillotta, M., Tsagris, M. and Fokianos, K. (2022c). The Rpackage PNAR for modelling count network time series. https://arxiv.org/abs/2211.02582
data(crime)
data(crime_W)
mod1 < lin_estimnarpq(crime, crime_W, p = 2)
ca < mod1$coefs[, 1]
a < score_test_nonlinpq_h0(ca, crime, crime_W, p = 2, d = 1)
print(a)
summary(a)
S3 methods for extracting the results of the estimation functions.
## S3 method for class 'PNAR'
summary(object, ...)
## S3 method for class 'summary.PNAR'
print(x, ...)
## S3 method for class 'PNAR'
print(x, ...)
object 
An object containing the results of the estimation function 
x 
An object containing the results of the estimation function 
... 
Extra arguments the user can pass. 
These functions print the output of the estimation functions.
The print.PNAR() function prints the coefficients of the model. The summary.PNAR() function prints the output in the lm() style.
Mirko Armillotta, Michail Tsagris and Konstantinos Fokianos.
Armillotta, M. and K. Fokianos (2022a). Poisson network autoregression. https://arxiv.org/abs/2104.06296
Armillotta, M., Tsagris, M. and Fokianos, K. (2022c). The Rpackage PNAR for modelling count network time series. https://arxiv.org/abs/2211.02582
data(crime)
data(crime_W)
mod1 < lin_estimnarpq(crime, crime_W, p = 2)
mod1
print(mod1)
summary(mod1)