Description Usage Arguments Details Value Note Author(s) References See Also Examples
Simulate a time series using a general autoregressive model.
1 2 3 |
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.
Version 0.1-2: bug corrected and garmaFit example given, Version 0.1-5: garmaFit example removed due to archiving of package gamlss.util
Olivier Briet o.briet@gmail.com
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.
'rnegbin(MASS)', 'arrep'.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 | 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))
|
Call:
arima(x = tsy, order = c(1, 0, 0), xreg = x)
Coefficients:
ar1 intercept x
0.8105 2.1212 0.6693
s.e. 0.0184 0.1643 0.0238
sigma^2 estimated as 0.9782: log likelihood = -1408.47, aic = 2824.95
Call:
arima(x = tsy, order = c(1, 1, 0), xreg = Xreg)
Coefficients:
ar1 Xreg1 Xreg2 Xreg3
0.3711 0.7010 4.0195 3.9885
s.e. 0.0294 0.0177 0.1217 0.1215
sigma^2 estimated as 0.9804: log likelihood = -1407.73, aic = 2825.45
ar check
1 -2.000000e-01 -2.000000e-01
2 1.200000e-01 1.200000e-01
3 -7.200000e-02 -7.200000e-02
4 4.320000e-02 4.320000e-02
5 -2.592000e-02 -2.592000e-02
6 1.555200e-02 1.555200e-02
7 -9.331200e-03 -9.331200e-03
8 5.598720e-03 5.598720e-03
9 -3.359232e-03 -3.359232e-03
10 2.015539e-03 2.015539e-03
11 -1.209324e-03 -1.209324e-03
12 7.255941e-04 7.255941e-04
13 -4.353565e-04 -4.353565e-04
14 2.612139e-04 2.612139e-04
15 -1.567283e-04 -1.567283e-04
16 9.403700e-05 9.403700e-05
17 -5.642220e-05 -5.642220e-05
18 3.385332e-05 3.385332e-05
19 -2.031199e-05 -2.031199e-05
20 1.218719e-05 1.218719e-05
21 -7.312317e-06 -7.312317e-06
22 4.387390e-06 4.387390e-06
23 -2.632434e-06 -2.632434e-06
24 1.579460e-06 1.579460e-06
25 -9.476763e-07 -9.476763e-07
26 5.686058e-07 5.686058e-07
27 -3.411635e-07 -3.411635e-07
28 2.046981e-07 2.046981e-07
29 -1.228188e-07 -1.228188e-07
30 7.369131e-08 7.369131e-08
31 -4.421478e-08 -4.421478e-08
32 2.652887e-08 2.652887e-08
33 -1.591732e-08 -1.591732e-08
34 9.550393e-09 9.550393e-09
35 -5.730236e-09 -5.730236e-09
36 3.438142e-09 3.438142e-09
37 -2.062885e-09 -2.062885e-09
38 1.237731e-09 1.237731e-09
39 -7.426386e-10 -7.426386e-10
40 4.455832e-10 4.455832e-10
41 -2.673499e-10 -2.673499e-10
42 1.604099e-10 1.604099e-10
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.