inst/doc/v23i05_table15a.R

#Script for comparing FGN/ARMA forecast performance
#This script requires a Beowulf cluster computer with Rmpi library
#attach library
library(FGN)
#input data 
data(NileMin)
y<-NileMin
#ARMA(p,q) model order
p<-2
q<-1
#
nB<-10^4 #Number of bootstrap iterations
n<-length(y) #length of series
K<-100 #number of out-of-sample data values
n1<-n-K #length of training series
sdy<-sd(y) #sd of original series
#syd is the long-run prediction sd
outy<-FitFGN(y) #fix model for bootstrap
#
#combine all into together except outy,
inputdata=c(p,q,n,K)

onebootFGNARMA=function(outy,inputdata){
    MAXIT <- 10
    p <- inputdata[1]
    q <- inputdata[2]
    n <- inputdata[3]
    K <- inputdata[4]
    n1=n-K
#
#FGN fit to z1 and forecast using z2. 
#FGN and ARMA use independent bootstraps
    NotOK <- TRUE
    ITER<-0
    while (NotOK && ITER<MAXIT){
	z<-Boot(outy)
 	z1<-z[1:n1] #training data
   	z2<-z[-(1:n1)] #testing data
    	outz1<-tryCatch(FitFGN(z1), error = function(e) FALSE)
    	if (!(is.logical(outz1)))  NotOK <- FALSE
    	else ITER<-ITER+1
    }
    H<-outz1$H
    mu<-outz1$muHat
    rFGN<-var(z1)*FGNAcf(0:(n+3-1), H)
    F<-TrenchForecast(c(z1,z2), rFGN, mu, n1, maxLead=3)$Forecasts
    nF<-nrow(F)
    err1<-z2-F[,1][-nF]
    err2<-z2[-1]-F[,2][-c(nF,(nF-1))]
    err3<-z2[-c(1,2)]-F[,3][-c(nF,(nF-1),(nF-2))]
    rmse1<-sqrt(mean(err1^2))
    rmse2<-sqrt(mean(err2^2))
    rmse3<-sqrt(mean(err3^2))
    FGNrmse<-c(ITER,rmse1,rmse2,rmse3)
#
#ARMA(p,q) fit to z1 and forecast using z2
    NotOK <- TRUE
    ITER<-0
    while (NotOK && ITER<MAXIT){
	z<-Boot(outy)
 	z1<-z[1:n1] #training data
   	z2<-z[-(1:n1)] #testing data
    	outz1<-tryCatch(arima(z1,c(p,0,q)), error = function(e) FALSE)
    	if (!(is.logical(outz1)))  NotOK <- FALSE
    	else ITER<-ITER+1
    }
    z1<-z[1:n1] #training data
    z2<-z[-(1:n1)] #testing data
    phi<-theta<-numeric(0)
    if (p>0) phi<-coef(ansz1)[1:p]
    if (q>0) theta<-coef(ansz1)[(p+1):(p+q)]
    zm<-coef(ansz1)[p+q+1]
    sigma2<-ansz1$sigma2
    r<-var(z)*ARMAacf(ar=phi, ma=theta, lag.max=n+3-1)
    F<-TrenchForecast(c(z1,z2), r, zm, n1, maxLead=3)$Forecasts
    err1<-z2-F[,1][-nF]
    err2<-z2[-1]-F[,2][-c(nF,(nF-1))]
    err3<-z2[-c(1,2)]-F[,3][-c(nF,(nF-1),(nF-2))]
    rmse1<-sqrt(mean(err1^2))
    rmse2<-sqrt(mean(err2^2))
    rmse3<-sqrt(mean(err3^2))
    ARMArmse<-c(ITER, rmse1, rmse2, rmse3)
#Average mse
    c(FGNrmse, ARMArmse)
}

####################################################
#start Rmpi and setup R slaves
library(Rmpi)
#setup slaves
mpi.spawn.Rslaves()
#send data and function to all slaves
mpi.bcast.Robj2slave(inputdata)
mpi.bcast.Robj2slave(outy)
mpi.bcast.Robj2slave(onebootFGNARMA)

#tell slaves to load FGN
mpi.bcast.cmd(library(FGN))

#setup parallel RNG. seed can be specified.
mpi.setup.rngstream(19000713)

startTime<-proc.time()
#start parallel bootstrap
out <- mpi.parReplicate(nB, onebootFGNARMA(outy=outy,inputdata=inputdata))
failed<-apply(out[c(1,5),], 1, sum)
averRMSE=apply(out[-c(1,5),], 1, mean)
endTime<-proc.time()
totalTime<-endTime-startTime
totalTime<-(totalTime/3600)[3]
#
#tabulate result
tb<-matrix(c(averRMSE[1:3],sdy,failed[1],averRMSE[4:6],sdy,failed[2]),ncol=2)
dimnames(tb)<-list(c("lead1","lead2","lead3","infty","failed"),c("FGN","ARMA"))
#
#output and close R
save(tb, out, totalTime, file="tbBoot.Rdata")  #
write(totalTime, "totalTimeMPI.txt")
write(t(round(tb,3)), file="tbMPI.txt", ncolumns=2)
print(paste("Elapsed Time = ", totalTime , "Hours"))
print(paste("nB = ",nB))
print(tb)

#close all slaves and quit
mpi.close.Rslaves()
mpi.quit()

Try the ltsa package in your browser

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

ltsa documentation built on May 2, 2019, 4:01 a.m.