Two Functions for Generalized SARIMA Time Series Simulation

Description

Write SARIMA models in (finite) AR representation and simulate generalized multiplicative seasonal autoregressive moving average (time) series The methodology of this method is described in Briet OJT, Amerasinghe PH, and Vounatsou P (2013) <doi:10.1371/journal.pone.0065761>.

Details

Use arrep() for converting the SARIMA function into AR representation, and use garsim() to simulate.

Author(s)

Olivier Briet <o.briet@gmail.com>

Maintainer: Olivier Briet <o.briet@gmail.com>

References

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.

Compute the Autoregressive Representation of a SARIMA Model

Description

Invert (invertible) SARIMA(p, d, q, P, D, Q) models to ar representation.

Usage

arrep(notation = "arima", phi = c(rep(0, 10)), d = 0, theta = c(rep(0, 10)),
 Phi = c(rep(0, 10)), D = 0, Theta = c(rep(0, 10)), frequency = 1)

Arguments

Details

For input, positive values of phi, theta, Phi and Theta indicate positive dependence. Implemented for p,q,P,Q element of c(0,1,2,3,4,5,6,7,8,9,10). The ar representation is truncated at coefficients less than 1.0e-10. Values of theta, Theta near non invertibility (-1 or 1) will not be practical and will cause near infinite lags, especially for Theta and large frequency.

Value

A vector containing a truncated autoregressive representation of a SARIMA model. This can be used as input for the function gar.sim.

Author(s)

Olivier Briet o.briet@gmail.com

References

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.

See Also

'garsim'.

Examples

phi<-c(0.5, 0.3, 0.1)
theta<-c(0.6, 0.2, 0.2)
ar<-arrep(phi=phi, theta=theta, frequency=12)
check<-(acf2AR(ARMAacf(ar=phi, ma=theta, lag.max = 100, pacf = FALSE))[100,1:length(ar)])
as.data.frame(cbind(ar,check))

phi<-c(0.2,0.5)
theta<-c(0.4)
Phi<-c(0.6)
Theta<-c(0.3)
d<-2
D<-1
frequency<-12
ar<-arrep(phi=phi, theta=theta, Phi=Phi, Theta=Theta, frequency= frequency, d=d, D=D)
N<-500
intercept<-10
data.sim <- garsim(n=(N+length(ar)),phi=ar, X=matrix(rep(intercept,(N+ length(ar)))),
beta=1, sd=1) 
y<-data.sim[1+length(ar): (N+length(ar))]
tsy<-ts(y, freq= frequency)
plot(tsy)
arima(tsy, order=c(2,2,1), seasonal=list(order=c(1,1,1)))

Simulate a Generalized Autoregressive Time Series

Description

Simulate a time series using a general autoregressive model.

Usage

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)

Arguments

Details

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

Value

An autoregressive series of length n. Note that the first p data do not have autoregressive structure.

Note

Version 0.1-2: bug corrected and garmaFit example given, Version 0.1-5: garmaFit example removed due to archiving of package gamlss.util

Author(s)

Olivier Briet o.briet@gmail.com

References

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.

See Also

'rnegbin(MASS)', 'arrep'.

Examples

  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))