garmaFit: A function to fit a GARMA model
In gamlss.util: GAMLSS Utilities
Description
Usage
Arguments
Details
Value
Note
Author(s)
References
See Also
Examples
Description
This function is for fitting a GARMA model, see Benjamin et al. (2003).
Usage
1
2
3
4
5
Arguments
formula
A formula for linear terms i.e. like in lm()
order
order specify the order of the generalised arm model
weights
prior weighs, they are working like in gamlss
data
the relevant data.frame
family
A gamlss.family distribution
alpha
This parameter is used in the definition of the link function of the response variable i.e. \log(y^*) will be y^*=max(y, α)
phi.start
starting values for the AR parameters
theta.start
starting values for the MA part
tail
how many observation from the tall of the response variable should be suppressed
control
control for optim() or nlminb() function use for optimisation.
Details
The model is described in Benjamin et al. (2003). The implementation here is more general that it allows all the gamlss.family distributions to be fitted rather than only for the exponential family which was described in the original paper. Note that in this formulation only the mu can be modelled as ARMA.
Value
It returns a fitted garma model.
Note
There is no check done whether the fitted model is stationary.
Author(s)
Mikis Stasinopoulos d.stasinopoulos@londonmet.ac.uk, Bob Rigby r.rigby@londonmet.ac.uk and Vlasios Voudouris
References
Benjamin M. A., Rigby R. A. and Stasinopoulos D.M. (2003) Generalised Autoregressive Moving Average Models. J. Am. Statist. Ass., 98, 214-223.
Rigby, R. A. and Stasinopoulos D. M. (2005). Generalized additive models for location, scale and shape,(with discussion),
Appl. Statist., 54, part 3, pp 507-554.
Stasinopoulos D. M., Rigby R.A. and Akantziliotou C. (2006) Instructions on how to use the GAMLSS package in R.
Accompanying documentation in the current GAMLSS help files, (see also http://www.gamlss.org/).
Stasinopoulos D. M. Rigby R.A. (2007) Generalized additive models for location scale and shape (GAMLSS) in R.
Journal of Statistical Software, Vol. 23, Issue 7, Dec 2007, http://www.jstatsoft.org/v23/i07.
See Also
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
data(polio)
ti <- as.numeric(time(polio))
mo <- as.factor(cycle(polio))
x1 <- 0:167 #Index used in Tutz p197
x2 <- cos(2*pi*x1/12)
x3 <- sin(2*pi*x1/12)
x4 <- cos(2*pi*x1/6)
x5 <- sin(2*pi*x1/6)
# all the data here
da <-data.frame(polio,x1,x2,x3,x4,x5, ti, mo)
rm(ti,mo,x1,x2,x3,x4,x5)
#-------------------------------------------------------------------
# with linear trend
m00 <- garmaFit(polio~x1+x2+x3+x4+x5, data=da, order=c(0,0), family=NBI, tail=3) #
m10 <- garmaFit(polio~x1+x2+x3+x4+x5, data=da, order=c(1,0), family=NBI, tail=3) #
## Not run:
m01 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(0,1), data=da, family=NBI, tail=3)
m20 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(2,0), data=da, family=NBI, tail=3)
m11 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(1,1), data=da, family=NBI, tail=3)
m02 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(0,2), data=da, family=NBI, tail=3)
m30 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(3,0), data=da, family=NBI, tail=3)
m21 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(2,1), data=da, family=NBI, tail=3)
m12 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(1,2), data=da, family=NBI, tail=3)
m03 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(0,3), data=da, family=NBI, tail=3)
AIC(m00,m10,m01,m20,m11,m02,m30,m21,m12,m03 , k=0)
AIC(m00,m10,m01,m20,m11,m02,m30,m21,m12,m03 , k=log(168))
# without linear trend
n00 <- garmaFit(polio~x2+x3+x4+x5, data=da, order=c(0,0), family=NBI, tail=3) #
n10 <- garmaFit(polio~x2+x3+x4+x5, data=da, order=c(1,0), family=NBI, tail=3) # OK
n01 <- garmaFit(polio~x2+x3+x4+x5, order=c(0,1), data=da, family=NBI, tail=3)
n20 <- garmaFit(polio~x2+x3+x4+x5, order=c(2,0), data=da, family=NBI, tail=3)
n11 <- garmaFit(polio~x2+x3+x4+x5, order=c(1,1), data=da, family=NBI, tail=3)
n02 <- garmaFit(polio~x2+x3+x4+x5, order=c(0,2), data=da, family=NBI, tail=3)
n30 <- garmaFit(polio~x2+x3+x4+x5, order=c(3,0), data=da, family=NBI, tail=3)
n21 <- garmaFit(polio~x2+x3+x4+x5, order=c(2,1), data=da, family=NBI, tail=3)
n12 <- garmaFit(polio~x2+x3+x4+x5, order=c(1,2), data=da, family=NBI, tail=3)
n03 <- garmaFit(polio~x2+x3+x4+x5, order=c(0,3), data=da, family=NBI, tail=3)
AIC(m00,n10,n01,n20,n11,n02,n30,n21,n12, k=0)
AIC(m00,n10,n01,n20,n11,n02,n30,n21,n12, k=log(168))
## End(Not run)
Example output
Loading required package: gamlss.dist
Loading required package: MASS
Loading required package: gamlss
Loading required package: splines
Loading required package: gamlss.data
Loading required package: nlme
Loading required package: parallel
********** GAMLSS Version 5.0-2 **********
For more on GAMLSS look at http://www.gamlss.org/
Type gamlssNews() to see new features/changes/bug fixes.
Loading required package: zoo
Attaching package: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
deviance of linear model= 507.656
deviance of linear model= 507.656
deviance of garma model= 494.8478
deviance of linear model= 507.656
deviance of garma model= 496.5089
deviance of linear model= 507.656
deviance of garma model= 490.655
deviance of linear model= 507.656
deviance of garma model= 492.868
deviance of linear model= 507.656
deviance of garma model= 487.9239
deviance of linear model= 507.656
deviance of garma model= 486.3349
deviance of linear model= 507.656
deviance of garma model= 486.1983
deviance of linear model= 507.656
deviance of garma model= 484.5784
deviance of linear model= 507.656
deviance of garma model= 484.6816
df AIC
m12 10 484.5784
m03 10 484.6816
m21 10 486.1983
m30 10 486.3349
m02 9 487.9239
m20 9 490.6550
m11 9 492.8680
m10 8 494.8478
m01 8 496.5089
m00 7 507.6560
df AIC
m02 9 534.0396
m12 10 535.8180
m10 8 535.8395
m03 10 535.9212
m20 9 536.7707
m21 10 537.4380
m01 8 537.5006
m30 10 537.5746
m11 9 538.9836
m00 7 543.5237
deviance of linear model= 513.0696
deviance of linear model= 513.0696
deviance of garma model= 498.9811
deviance of linear model= 513.0696
deviance of garma model= 501.3473
deviance of linear model= 513.0696
deviance of garma model= 493.5762
deviance of linear model= 513.0696
deviance of garma model= 495.0234
deviance of linear model= 513.0696
deviance of garma model= 491.0733
deviance of linear model= 513.0696
deviance of garma model= 490.305
deviance of linear model= 513.0696
deviance of garma model= 489.7187
deviance of linear model= 513.0696
deviance of garma model= 488.452
deviance of linear model= 513.0696
Error in filter(some, -theta, "r", init = rep(0, order[2])) :
missing values in 'filter'
OOPS IT FAILED: let us try once more
deviance of linear model= 513.0696
deviance of garma model= 488.7644
df AIC
n12 9 488.4520
n21 9 489.7187
n30 9 490.3050
n02 8 491.0733
n20 8 493.5762
n11 8 495.0234
n10 7 498.9811
n01 7 501.3473
m00 7 507.6560
df AIC
n02 8 532.0650
n12 9 534.5676
n20 8 534.5679
n10 7 534.8489
n21 9 535.8344
n11 8 536.0151
n30 9 536.4207
n01 7 537.2151
m00 7 543.5237
gamlss.util documentation built on May 2, 2019, 7:10 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
This function is for fitting a GARMA model, see Benjamin et al. (2003).
Usage
1 2 3 4 5 |
Arguments
formula |
A formula for linear terms i.e. like in |
order |
|
weights |
prior weighs, they are working like in |
data |
the relevant |
family |
A |
alpha |
This parameter is used in the definition of the link function of the response variable i.e. \log(y^*) will be y^*=max(y, α) |
phi.start |
starting values for the AR parameters |
theta.start |
starting values for the MA part |
tail |
how many observation from the tall of the response variable should be suppressed |
control |
control for |
Details
The model is described in Benjamin et al. (2003). The implementation here is more general that it allows all the gamlss.family distributions to be fitted rather than only for the exponential family which was described in the original paper. Note that in this formulation only the mu can be modelled as ARMA.
Value
It returns a fitted garma model.
Note
There is no check done whether the fitted model is stationary.
Author(s)
Mikis Stasinopoulos d.stasinopoulos@londonmet.ac.uk, Bob Rigby r.rigby@londonmet.ac.uk and Vlasios Voudouris
References
Benjamin M. A., Rigby R. A. and Stasinopoulos D.M. (2003) Generalised Autoregressive Moving Average Models. J. Am. Statist. Ass., 98, 214-223.
Rigby, R. A. and Stasinopoulos D. M. (2005). Generalized additive models for location, scale and shape,(with discussion), Appl. Statist., 54, part 3, pp 507-554.
Stasinopoulos D. M., Rigby R.A. and Akantziliotou C. (2006) Instructions on how to use the GAMLSS package in R. Accompanying documentation in the current GAMLSS help files, (see also http://www.gamlss.org/).
Stasinopoulos D. M. Rigby R.A. (2007) Generalized additive models for location scale and shape (GAMLSS) in R. Journal of Statistical Software, Vol. 23, Issue 7, Dec 2007, http://www.jstatsoft.org/v23/i07.
See Also
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 | data(polio)
ti <- as.numeric(time(polio))
mo <- as.factor(cycle(polio))
x1 <- 0:167 #Index used in Tutz p197
x2 <- cos(2*pi*x1/12)
x3 <- sin(2*pi*x1/12)
x4 <- cos(2*pi*x1/6)
x5 <- sin(2*pi*x1/6)
# all the data here
da <-data.frame(polio,x1,x2,x3,x4,x5, ti, mo)
rm(ti,mo,x1,x2,x3,x4,x5)
#-------------------------------------------------------------------
# with linear trend
m00 <- garmaFit(polio~x1+x2+x3+x4+x5, data=da, order=c(0,0), family=NBI, tail=3) #
m10 <- garmaFit(polio~x1+x2+x3+x4+x5, data=da, order=c(1,0), family=NBI, tail=3) #
## Not run:
m01 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(0,1), data=da, family=NBI, tail=3)
m20 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(2,0), data=da, family=NBI, tail=3)
m11 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(1,1), data=da, family=NBI, tail=3)
m02 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(0,2), data=da, family=NBI, tail=3)
m30 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(3,0), data=da, family=NBI, tail=3)
m21 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(2,1), data=da, family=NBI, tail=3)
m12 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(1,2), data=da, family=NBI, tail=3)
m03 <- garmaFit(polio~x1+x2+x3+x4+x5, order=c(0,3), data=da, family=NBI, tail=3)
AIC(m00,m10,m01,m20,m11,m02,m30,m21,m12,m03 , k=0)
AIC(m00,m10,m01,m20,m11,m02,m30,m21,m12,m03 , k=log(168))
# without linear trend
n00 <- garmaFit(polio~x2+x3+x4+x5, data=da, order=c(0,0), family=NBI, tail=3) #
n10 <- garmaFit(polio~x2+x3+x4+x5, data=da, order=c(1,0), family=NBI, tail=3) # OK
n01 <- garmaFit(polio~x2+x3+x4+x5, order=c(0,1), data=da, family=NBI, tail=3)
n20 <- garmaFit(polio~x2+x3+x4+x5, order=c(2,0), data=da, family=NBI, tail=3)
n11 <- garmaFit(polio~x2+x3+x4+x5, order=c(1,1), data=da, family=NBI, tail=3)
n02 <- garmaFit(polio~x2+x3+x4+x5, order=c(0,2), data=da, family=NBI, tail=3)
n30 <- garmaFit(polio~x2+x3+x4+x5, order=c(3,0), data=da, family=NBI, tail=3)
n21 <- garmaFit(polio~x2+x3+x4+x5, order=c(2,1), data=da, family=NBI, tail=3)
n12 <- garmaFit(polio~x2+x3+x4+x5, order=c(1,2), data=da, family=NBI, tail=3)
n03 <- garmaFit(polio~x2+x3+x4+x5, order=c(0,3), data=da, family=NBI, tail=3)
AIC(m00,n10,n01,n20,n11,n02,n30,n21,n12, k=0)
AIC(m00,n10,n01,n20,n11,n02,n30,n21,n12, k=log(168))
## End(Not run)
|
Example output
Loading required package: gamlss.dist
Loading required package: MASS
Loading required package: gamlss
Loading required package: splines
Loading required package: gamlss.data
Loading required package: nlme
Loading required package: parallel
********** GAMLSS Version 5.0-2 **********
For more on GAMLSS look at http://www.gamlss.org/
Type gamlssNews() to see new features/changes/bug fixes.
Loading required package: zoo
Attaching package: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
deviance of linear model= 507.656
deviance of linear model= 507.656
deviance of garma model= 494.8478
deviance of linear model= 507.656
deviance of garma model= 496.5089
deviance of linear model= 507.656
deviance of garma model= 490.655
deviance of linear model= 507.656
deviance of garma model= 492.868
deviance of linear model= 507.656
deviance of garma model= 487.9239
deviance of linear model= 507.656
deviance of garma model= 486.3349
deviance of linear model= 507.656
deviance of garma model= 486.1983
deviance of linear model= 507.656
deviance of garma model= 484.5784
deviance of linear model= 507.656
deviance of garma model= 484.6816
df AIC
m12 10 484.5784
m03 10 484.6816
m21 10 486.1983
m30 10 486.3349
m02 9 487.9239
m20 9 490.6550
m11 9 492.8680
m10 8 494.8478
m01 8 496.5089
m00 7 507.6560
df AIC
m02 9 534.0396
m12 10 535.8180
m10 8 535.8395
m03 10 535.9212
m20 9 536.7707
m21 10 537.4380
m01 8 537.5006
m30 10 537.5746
m11 9 538.9836
m00 7 543.5237
deviance of linear model= 513.0696
deviance of linear model= 513.0696
deviance of garma model= 498.9811
deviance of linear model= 513.0696
deviance of garma model= 501.3473
deviance of linear model= 513.0696
deviance of garma model= 493.5762
deviance of linear model= 513.0696
deviance of garma model= 495.0234
deviance of linear model= 513.0696
deviance of garma model= 491.0733
deviance of linear model= 513.0696
deviance of garma model= 490.305
deviance of linear model= 513.0696
deviance of garma model= 489.7187
deviance of linear model= 513.0696
deviance of garma model= 488.452
deviance of linear model= 513.0696
Error in filter(some, -theta, "r", init = rep(0, order[2])) :
missing values in 'filter'
OOPS IT FAILED: let us try once more
deviance of linear model= 513.0696
deviance of garma model= 488.7644
df AIC
n12 9 488.4520
n21 9 489.7187
n30 9 490.3050
n02 8 491.0733
n20 8 493.5762
n11 8 495.0234
n10 7 498.9811
n01 7 501.3473
m00 7 507.6560
df AIC
n02 8 532.0650
n12 9 534.5676
n20 8 534.5679
n10 7 534.8489
n21 9 535.8344
n11 8 536.0151
n30 9 536.4207
n01 7 537.2151
m00 7 543.5237
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.