garsim: Simulate a Generalized Autoregressive Time Series
In gsarima: Two Functions for Generalized SARIMA Time Series Simulation
Description
Usage
Arguments
Details
Value
Note
Author(s)
References
See Also
Examples
Description
Simulate a time series using a general autoregressive model.
Usage
1
2
3
Arguments
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").
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
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))
Example output
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
gsarima documentation built on Sept. 4, 2020, 1:08 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Note Author(s) References See Also Examples
Description
Simulate a time series using a general autoregressive model.
Usage
1 2 3 |
Arguments
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"). |
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
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))
|
Example output
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
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.