garmaFit: A function to fit a GARMA model

Description Usage Arguments Details Value Note Author(s) References See Also Examples

View source: R/garmaFit.R

Description

This function is for fitting a GARMA model, see Benjamin et al. (2003).

Usage

1
2
3
4
5
garmaFit(formula = formula(data), order = c(0, 0), 
         weights = NULL, data = sys.parent(), 
         family = NO(), alpha = 0.1, 
         phi.start = NULL, theta.start = NULL, 
         tail = max(order), control = list())

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 [email protected], Bob Rigby [email protected] 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

gamlss.family, gamlss

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 29, 2017, 8:58 p.m.