Simulate a Generalized Autoregressive Time Series
Simulate a time series using a general autoregressive model.
garsim(n, phi, X = matrix(0, nrow = n), beta = as.matrix(0), sd = 1, family = "gaussian", transform.Xbeta = "identity", link = "identity", minimum = 0, zero.correction = "zq1", c = 1, theta = 0)
n: The number of simulated values.phi: A vector of autoregressive parameters of length p.X: An n by m optional matrix of external covariates, optionally including an intercept (recommended for family = "poisson").beta: An m vector of coefficients.sd: Standard deviation for Gaussian family.family: Distribution family, defaults to "gaussian".transform.Xbeta: Optional transformation for the product of covariates and coefficients, see Details.link: The link function, defaults to "identity".minimum: A minimum value for the mean parameter of the Poisson and Negative Binomialdistributions (only applicable for link= "identity" and family = c("poisson","negative.binomial")). Defaults to 0. A small positive value will allow non-stationary series to "grow" after encountering a simulated value of 0.zero.correction: Method for transformation for dealing with zero values (only applicable when link = "log"), see Details.c: The constant used for transformation before taking the logarithm (only applicable when link = "log"). A value between 0 and 1 is recommended.theta: Parameter theta (for family = "negative.binomial").Implemented are the following models: 1) family = "gaussian", link = "identity" 2) family = "poisson", link = "identity" 3) family = "poisson", link = "identity", transform.Xbeta = "exponential" 4) family = "poisson", link = "log", zero.correction = "zq1" 5) family = "poisson", link = "log", zero.correction = "zq2" 6) family = "negative.binomial", link = "identity" 7) family = "negative.binomial", link = "identity", transform.Xbeta = "exponential" 8) family = "negative.binomial", link = "log", zero.correction = "zq1" 9) family = "negative.binomial", link = "log", zero.correction = "zq2"
Models 1 to 4 are within the family of GARMA models of Benjamin and colleagues 2003 Model 2 is the extension to higher order p of a Poisson CLAR(1) model proposed by Grunwald and colleagues (2000). Model 3 is a modification of the PAR(p) data generating process (https://personal.utdallas.edu/~pxb054000/code/pests.r) of Brandt and Williams (2001). Note that for psi = 0, the model reduces to a standard Poisson model with log-link function. For a model without external variables (only an intercept), the transformation of Xbeta has no consequence and then model 3 is the same as model 2. Model 4 corresponds to model 2.2 of Zeger and Qaqish (1988). The value c is only added to values of zero prior to taking the log. Models 6 to 9 are similar but with negative binomial distribution
An autoregressive series of length n. Note that the first p data do not have autoregressive structure.
Briet OJT, Amerasinghe PH, Vounatsou P: Generalized seasonal autoregressive integrated moving average models for count data with application to malaria time series with low case numbers. PLoS ONE, 2013, 8(6): e65761. doi:10.1371/journal.pone.0065761 https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0065761 If you use the gsarima package, please cite the above reference.
Brandt PT, Williams JT: A linear Poisson autoregressive model: The PAR(p). Political Analysis 2001, 9.
Benjamin MA, Rigby RA, Stasinopoulos DM: Generalized Autoregressive Moving Average Models. Journal of the American Statistical Association 2003, 98:214-223.
Zeger SL, Qaqish B: Markov regression models for time series: a quasi-likelihood approach. Biometrics 1988, 44:1019-1031
Grunwald G, Hyndman R, Tedesco L, Tweedie R: Non-Gaussian conditional linear AR(1) models. Australian & New Zealand Journal of Statistics 2000, 42:479-495.
Olivier Briet o.briet@gmail.com
Version 0.1-2: bug corrected and garmaFit example given, Version 0.1-5: garmaFit example removed due to archiving of package gamlss.util
'rnegbin(MASS)', 'arrep'.
N<-1000 ar<-c(0.8) intercept<-2 frequency<-1 x<- rnorm(N) beta.x<-0.7 #Gaussian simulation with covariate X=matrix(c(rep(intercept, N+length(ar)), rep(0, length(ar)), x), ncol=2) y.sim <- garsim(n=(N+length(ar)),phi=ar, X=X, beta=c(1,beta.x), sd=sqrt(1)) y<-y.sim[(1+length(ar)):(N+length(ar))] tsy<-ts(y, freq=frequency) plot(tsy) arima(tsy, order=c(1,0,0), xreg=x) #Gaussian simulation with covariate and deterministic seasonality through first order harmonic ar<-c(1.4,-0.4) frequency<-12 beta.x<-c(0.7,4,4) X<-matrix(nrow= (N+ length(ar)), ncol=3) for (t in 1: length(ar)){ X[t,1]<-0 X[t,2]<-sin(2*pi*(t- length(ar))/frequency) X[t,3]<- cos(2*pi*(t- length(ar))/frequency) } for (t in (1+ length(ar)): (N+ length(ar))){ X[t,1]<-x[t- length(ar)] X[t,2]<-sin(2*pi*(t- length(ar))/frequency) X[t,3]<- cos(2*pi*(t- length(ar))/frequency) } y.sim <- garsim(n=(N+length(ar)),phi=ar, X=X, beta= beta.x, sd=sqrt(1)) y<-y.sim[(1+length(ar)):(N+length(ar))] tsy<-ts(y, freq=frequency) plot(tsy) Xreg<-matrix(nrow= N, ncol=3) for (t in 1: N){ Xreg[t,1]<-x[t] Xreg[t,2]<-sin(2*pi*t/frequency) Xreg[t,3]<- cos(2*pi*t/frequency) } arimares<-arima(tsy, order=c(1,1,0), xreg=Xreg) tsdiag(arimares) arimares #Negative binomial simulation with covariate ar<-c(0.8) frequency<-1 beta.x<-0.7 X=matrix(c(rep(log(intercept), N+length(ar)), rep(0, length(ar)), x), ncol=2) y.sim <- garsim(n=(N+length(ar)), phi=ar, beta=c(1,beta.x), link= "log", family= "negative.binomial", zero.correction = "zq1", c=1, theta=5, X=X) y<-y.sim[(1+length(ar)):(N+length(ar))] tsy<-ts(y, freq=frequency) plot(tsy) #Poisson ARMA(1,1) with identity link and negative auto correlation N<-500 phi<-c(-0.8) theta<-c(0.6) ar<-arrep(phi=phi, theta=theta) check<-(acf2AR(ARMAacf(ar=phi, ma=theta, lag.max = 100, pacf = FALSE))[100,1:length(ar)]) as.data.frame(cbind(ar,check)) intercept<-100 frequency<-1 X=matrix(c(rep(intercept, N+length(ar))), ncol=1) y.sim <- garsim(n=(N+length(ar)), phi=ar, beta=c(1), link= "identity", family= "poisson", minimum = -100, X=X) y<-y.sim[(1+length(ar)):(N+length(ar))] tsy<-ts(y, freq=frequency) plot(tsy) #Poisson AR(1) with identity link and negative auto correlation N<-1000 ar<-c(-0.8) intercept<-100 frequency<-1 X=matrix(c(rep(intercept, N+length(ar))), ncol=1) y.sim <- garsim(n=(N+length(ar)), phi=ar, beta=c(1), link= "identity", family= "poisson", minimum = -100, X=X) y<-y.sim[(1+length(ar)):(N+length(ar))] tsy<-ts(y, freq=frequency) plot(tsy) #Example of negative binomial GSARIMA(2,1,0,0,0,1)x phi<-c(0.5,0.2) theta<-c(0) Theta<-c(0.5) Phi<-c(0) d<-c(1) D<-c(0) frequency<-12 ar<-arrep(phi=phi, theta=theta, Phi=Phi, Theta=Theta, frequency= frequency, d=d, D=D) N<-c(1000) intercept<-c(10) x<- rnorm(N) beta.x<-c(0.7) X<-matrix(c(rep(log(intercept), N+length(ar)), rep(0, length(ar)), x), ncol=2) c<-c(1) y.sim <- garsim(n=(N+length(ar)), phi=ar, beta=c(1,beta.x), link= "log", family= "negative.binomial", zero.correction = "zq1", c=c, theta=5, X=X) y<-y.sim[(1+length(ar)):(N+length(ar))] tsy<-ts(y, freq=frequency) plot(tsy) plot(log(tsy))