data(NiagaraToxic)
head(NiagaraToxic)
#Example from thesis
## Not run:
# #Diagnostic checks and bootstrap confidence intervals
# Zdf <- NiagaraTmeans <- apply(A, MARGIN=2, fun=mean)oxic
# z <- log(Zdf$toxic)
# iz <- c("o", "L")[1+Zdf$cQ]
# #
# #CENARMA(1,1)
# cenarma(z, iz, p=1, q=1)
# #fit CENAR(1)
# cenarma(z, iz, p=1)
# #
# #diagnostic checks########
# #test CENARMA(1,1)
# SimModel <- function(OUTCENARMA){
# outSim <- boot.Cenarma(OUTCENARMA)
# }
# FitModel <- function(outSim){
# z <- outSim$y
# iz <- outSim$iy
# ans <- cenarma(z, iz, p=1, q=1)
# res <- resid(ans$outarima)
# list(res=res)
# }
# OUTCENARMA <- cenarma(y=NiagaraToxic$toxic, iy=c("o", "L")[1+NiagaraToxic$cQ], p=1, q=1)
# func <- list(SimModel=SimModel, FitModel=FitModel)
# start.time <- proc.time()[3]
# outp <- portest(OUTCENARMA$outarima, lags=5:25, nslaves=8, NREP=10^3, func=func, test="LjungBox")
# total.time <- proc.time()[3]-start.time
# total.time
# plot(outp[,1], outp[,4], xlab="lag", ylab="P-Value", cex=1.5, col="blue", pch=18, ylim=c(0,1))
# abline(col="red", h=0.05)
# #
# #test CENAR(1)
# SimModel <- function(OUTCENARMA){
# boot.Cenarma(OUTCENARMA)
# }
# FitModel <- function(outSim){
# z <- outSim$y
# iz <- outSim$iy
# ans <- cenarma(z, iz, p=1)
# res <- resid(ans$outarima)
# list(res=res)
# }
# OUTCENARMA <- cenarma(y=log(NiagaraToxic$toxic), iy=c("o", "L")[1+NiagaraToxic$cQ], p=1)
# func <- list(SimModel=SimModel, FitModel=FitModel)
# start.time <- proc.time()[3]
# outp <- portest(OUTCENARMA$outarima, lags=5:25, nslaves=8, NREP=10^3, func=func, test="LjungBox")
# total.time <- proc.time()[3]-start.time
# total.time
# plot(outp[,1], outp[,4], xlab="lag", ylab="P-Value", cex=1.5, col="blue", pch=18, ylim=c(0,1))
# abline(col="red", h=0.05)
# #
# #bootstrap confidence intervals
# #CENARMA(1,1)
# OUTCENARMA <- cenarma(y=log(NiagaraToxic$toxic), iy=c("o", "L")[1+NiagaraToxic$cQ], p=1, q=1)
# nboot <- 1000
# A <- matrix(numeric(nboot*3), ncol=3)
# start.time <- proc.time()[3]
# for (iboot in 1:nboot){
# out <- boot.Cenarma(OUTCENARMA)
# A[iboot, ] <- coef(cenarma(y=out$y, iy=out$iy, p=1, q=1)$outarima)
# }
# total.time <- proc.time()[3]-start.time
# total.time
# means <- apply(A, MARGIN=2, FUN=mean)
# means
# LO <- apply(A, MARGIN=2, function(x) quantile(x, 0.025))
# HI <- apply(A, MARGIN=2, function(x) quantile(x, 0.975))
# ansARMA11 <- matrix(c(LO,HI), nrow=3, dimnames=list(c("phi","theta","mu"), c("lo", "hi")))
# #CEAR(1)
# OUTCENARMA <- cenarma(y=log(NiagaraToxic$toxic), iy=c("o", "L")[1+NiagaraToxic$cQ], p=1)
# nboot <- 1000
# A <- matrix(numeric(nboot*2), ncol=2)
# start.time <- proc.time()[3]
# for (iboot in 1:nboot){
# out <- boot.Cenarma(OUTCENARMA)
# A[iboot, ] <- coef(cenarma(y=out$y, iy=out$iy, p=1)$outarima)
# }
# total.time <- proc.time()[3]-start.time
# total.time
# means <- apply(A, MARGIN=2, FUN=mean)
# means
# LO <- apply(A, MARGIN=2, function(x) quantile(x, 0.025))
# HI <- apply(A, MARGIN=2, function(x) quantile(x, 0.975))
# ansAR2 <- matrix(c(LO,HI), nrow=2, dimnames=list(c("phi","mu"), c("lo", "hi")))
# #summary
# ansARMA11
# ansAR2
# ## End(Not run)
Run the code above in your browser using DataLab