inst/unitTests/runit.likTests.R In glarma: Generalized Linear Autoregressive Moving Average Models

```### Unit tests of function likTests

### Functions with name test.* are run by R CMD check or by make if
### LEVEL=1 in call to make
### Functions with name levelntest.* are run by make if
### LEVEL=n in call to make
### Functions with name graphicstest.* are run by make if
### LEVEL=graphics in call to make

test.likTests <- function(){
## Poisson Distribution
data(Polio)
y <- Polio[, 2]
X <- as.matrix(Polio[, 3:8])
glarmamod <- glarma(y, X, thetaLags = c(1,2,5), type = "Poi", method = "FS",
residuals = "Pearson", maxit = 100, grad = 2.22e-16)

lik.tests1 <- likTests(glarmamod)

## Test for LR test
standard.LRTest1 <- c(27.1926023919697, 5.36460455780041e-06)
checkTrue(abs(lik.tests1[1, 1] - standard.LRTest1[1]) < 10^(-7))
checkTrue(abs(lik.tests1[1, 2] - standard.LRTest1[2]) < 10^(-7))

## Test for Wald test
standard.WaldTest1 <- c(38.1193252957133, 2.66674532456435e-08)
checkTrue(abs(lik.tests1[2, 1] - standard.WaldTest1[1]) < 10^(-7))
checkTrue(abs(lik.tests1[2, 2] - standard.WaldTest1[2]) < 10^(-7))

## Binomial Distribution
data(RobberyConvict)
datalen <- dim(RobberyConvict)[1]
monthmat <- matrix(0, nrow = datalen, ncol = 12)
dimnames(monthmat) <- list(NULL, c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
"Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
months <- unique(months(strptime(RobberyConvict\$Date, "%m/%d/%Y"),
abbreviate = TRUE))

for (j in 1:12) {
monthmat[months(strptime(RobberyConvict\$Date, "%m/%d/%Y"),
abbreviate = TRUE) == months[j], j] <-1
}

RobberyConvict <- cbind(rep(1, datalen), RobberyConvict, monthmat)
rm(monthmat)

## LOWER COURT ROBBERY
y1 <- RobberyConvict\$LC.Y
n1 <- RobberyConvict\$LC.N

Y <- cbind(y1, n1 - y1)

glm.LCRobbery <- glm(Y ~ -1 + Incpt + Step.2001 +
I(Feb + Mar + Apr + May + Jun + Jul) +
I(Aug + Sep + Oct + Nov + Dec),
data = RobberyConvict, family = binomial(link = logit),
na.action = na.omit,x = TRUE)

X <- glm.LCRobbery\$x

## Newton Raphson
glarmamod <- glarma(Y, X, phiLags = c(1), type = "Bin", method = "NR",
residuals = "Pearson", maxit = 100, grad = 1e-6)

lik.tests2 <- likTests(glarmamod)
standard.LRTest2 <- c(6.11042760345322, 0.0134386612016913)
checkTrue(abs(lik.tests2[1, 1] - standard.LRTest2[1]) < 10^(-7))
checkTrue(abs(lik.tests2[1, 2] - standard.LRTest2[2]) < 10^(-7))
standard.WaldTest2 <- c(6.14430921515927, 0.0131835654128953)
checkTrue(abs(lik.tests2[2, 1] - standard.WaldTest2[1]) < 10^(-7))
checkTrue(abs(lik.tests2[2, 2] - standard.WaldTest2[2]) < 10^(-7))

## Negative Binomial Distribution
data(Asthma)
y <- Asthma[, 1]
X <- as.matrix(Asthma[, 2:16])
glarmamod <- glarma(y, X, thetaLags = 7, type = "NegBin", method = "NR",
residuals = "Pearson", maxit = 100, grad = 1e-6)

lik.tests3 <- likTests(glarmamod)
standard.LRTest3 <- c(7.04717827753939, 0.00793902172903205)
checkTrue(abs(lik.tests3[1, 1] - standard.LRTest3[1]) < 10^(-7))
checkTrue(abs(lik.tests3[1, 2] - standard.LRTest3[2]) < 10^(-7))
standard.WaldTest3 <- c(5.14686152737022, 0.0232884286514632)
checkTrue(abs(lik.tests2[2, 1] - standard.WaldTest2[1]) < 10^(-7))
checkTrue(abs(lik.tests2[2, 2] - standard.WaldTest2[2]) < 10^(-7))
}
```

Try the glarma package in your browser

Any scripts or data that you put into this service are public.

glarma documentation built on Feb. 9, 2018, 6:08 a.m.