#Script for comparing FGN/ARMA forecast performance NileMin
#Also compare forecasts for TrenchForecast and predict.Arima
data(NileMin)
z<-NileMin
K<-100 #number of out-of-sample data values
z1<-z[1:(length(z)-K)] #training data
z2<-z[-(1:(length(z)-K))] #testing data
sdz<-sd(z2)
#
err1<-err2<-err3<-numeric(K+1)
for (i in 1:K) {
    if (i > 1)
        zz<-c(z1,z2[1:i])
    else 
        zz<-z1
    outzz<-FitFGN(zz)
    H<-outzz$H
    mu<-outzz$muHat
    rFGN<-var(zz)*FGNAcf(0:(length(zz)-1), H)
    znew<-z2[i]
    zf<-TrenchForecast(c(zz,znew), r, zm, length(zz), maxLead=3)$Forecasts
    err1[i]<-z2[i]-zf[1]
    if (i < K)   err2[i]<-z2[i+1]-zf[2]
    if (i < K-1) err3[i]<-z2[i+2]-zf[3]
    }
err2<-err2[-K]
err3<-err3[-c(K-1,K)]
rmse1<-sqrt(mean(err1^2))
rmse2<-sqrt(mean(err2^2))
rmse3<-sqrt(mean(err3^2))
FGNrmse<-c(rmse1,rmse2,rmse3)
#
#ARMA(p,q) fit recursively, use TrenchForecast
p<-2
q<-1
err1<-err2<-err3<-numeric(K+1)
for (i in 1:K) {
    if (i > 1)
        zz<-c(z1,z2[1:i])
    else 
        zz<-z1
    anszz<-arima(zz, c(p,0,q))
phi<-coef(anszz)[1:p]
theta<-coef(anszz)[(p+1):(p+q)]
zm<-coef(anszz)[p+q+1]
r<-var(zz)*ARMAacf(ar=phi, ma=theta, lag.max=n1-1)
znew<-z2[i]
    zf<-TrenchForecast(c(zz,znew), r, zm, length(zz), maxLead=3)$Forecasts
    err1[i]<-z2[i]-zf[1]
    if (i < K)   err2[i]<-z2[i+1]-zf[2]
    if (i < K-1) err3[i]<-z2[i+2]-zf[3]
    }
err2<-err2[-K]
err3<-err3[-c(K-1,K)]
rmse1<-sqrt(mean(err1^2))
rmse2<-sqrt(mean(err2^2))
rmse3<-sqrt(mean(err3^2))
ARMArmse<-c(rmse1,rmse2,rmse3)
#
#
#ARMA(p,q) fit recursively, use predict.ARIMA
err1<-err2<-err3<-numeric(K+1)
for (i in 1:K) {
    if (i > 1)
        zz<-c(z1,z2[1:i])
    else 
        zz<-z1
    zf<-predict(arima(zz, order=c(2,0,1)), n.ahead=3)$pred
    err1[i]<-z2[i]-zf[1]
    if (i < K)   err2[i]<-z2[i+1]-zf[2]
    if (i < K-1) err3[i]<-z2[i+2]-zf[3]
    }
err2<-err2[-K]
err3<-err3[-c(K-1,K)]
rmse1<-sqrt(mean(err1^2))
rmse2<-sqrt(mean(err2^2))
rmse3<-sqrt(mean(err3^2))
ARMAPrmse<-c(rmse1,rmse2,rmse3)
#
#tabulate result
tb<-matrix(c(FGNrmse,sdz,ARMArmse,sdz,ARMAPrmse,sdz),ncol=3)
dimnames(tb)<-list(c("lead1","lead2","lead3","infty"),c("FGN","ARMA/Trench","ARMA/predict"))