which of ANN (nnetar) and ARIMA are better forecasting methods?

which of ANN (nnetar) and ARIMA are better forecasting methods?

**I conducted a forecasting study on water influent characteristics to
water treatment plant with to well-known ANN (nnetar) and arima from
forecast package. i write below code in r for doing cross validation on
the fotecasted values by medels. i modefied the example code for
crass-validation by Prof. Hyndman, i am not sure this modified code be
free of mistakes. but it provide interesting results. in many literature i
found that ANN models have better performances in forecasting because
these models can find non linearity behavior of the natural processes. but
after running my code MAE for cross-validation analysis and
Diebold-Mariano test says that ARIMA model have more accuracy than
nnetar!! please let me know am i in mistakes regard my code or this result
can be approved.
Var<-Alk " Alkalinity of water
library(forecast)
tsVar<- ts(Alk, start=1,end=length(Var),frequency=1)
library(Hmisc)
label(tsVar) <- "mgCaCO3/li"
label(Label)<-"Alk"
length(Var)
h= 10 # the horizon for forecasting
#End point of modeling
End=(length(Var)-h)
library(AID)
library(forecast)
L<-boxcoxnc(na.omit(Var), method = "sw", lam = seq(-5,5,0.01), plotit =
TRUE, rep = 30, p.method = "BY")
Lambda<- L$result[1,]
k <- End
n <- length(tsVar)
m<- n-k
For1<- matrix(NA,n-k,h)
For2<- matrix(NA,n-k,h)
Obs<- matrix(NA,n-k,h)
mae1<- matrix(NA,n-k,h)
mae2 <- matrix(NA,n-k,h)
error1<- matrix(NA,n-k,h)
error2 <- matrix(NA,n-k,h)
for(i in 1:(n-k))
{
yshort <- window(tsVar, end =(End-1)+i)
ynext <- window(tsVar, start= End+i, end=length(tsVar) )
fit1 <- nnetar(yshort,8, , repeats=20 ,lambda=Lambda)
fcast1 <- forecast(fit1, h=(h+1)-i)
fit2 <- Arima(yshort, order=c(1,1,1), seasonal=list(order=c(0,0,0)),
lambda= Lambda)
fcast2 <- forecast(fit2, h=(h+1)-i)
Obs[i, (m-length(ynext))+1:length(ynext)] <- ynext
For1[i, (m-length(ynext))+1:length(ynext)] <- abs(fcast1[['mean']])
For2[i, (m-length(ynext))+1:length(ynext)] <- abs(fcast2[['mean']])
mae1[i, (m-length(ynext))+1:length(ynext)] <- abs(fcast1[['mean']]-ynext)
mae2[i, (m-length(ynext))+1:length(ynext)] <- abs(fcast2[['mean']]-ynext)
error1[i, (m-length(ynext))+1:length(ynext)] <- (ynext-fcast1[['mean']])
error2[i, (m-length(ynext))+1:length(ynext)] <- (ynext-fcast2[['mean']])
}
# The Diebold-Mariano test compares the forecast accuracy of two forecast
# methods.
dm.test(colMeans(error1, na.rm=TRUE), colMeans(error2,na.rm=TRUE),
alternative=c("greater"),h=10, power=2)
# H0= two methods have the same forecast accuracy
#if "greater" selected H1= method 2 (ARIMA) is more accurate than method 1
mean(mae1, na.rm=TRUE)#MAE for NNAR
mean(mae2, na.rm=TRUE)#MAE for ARIMA
dev.new()
par(mfrow = c(2,2))
plot(1:h, colMeans(mae1,na.rm=TRUE), type="l", col=2, xlab="horizon",
ylab="MAE",ylim=c(min (colMeans(mae1,na.rm=TRUE)), max
(colMeans(mae1,na.rm=TRUE))))
lines(1:h, colMeans(mae2,na.rm=TRUE), type="l",col=3)
legend("topleft",legend=c("NNAR","ARIMA"),col=2:4,lty=1)
plot(1:h, colMeans(For1,na.rm=TRUE), type="l", col=2, xlab="horizon",
ylab="Forecasted",ylim=c(min (colMeans(Obs,na.rm=TRUE)), max
(colMeans(Obs,na.rm=TRUE))))
lines(1:h, colMeans(For2,na.rm=TRUE), type="l",col=3)
lines(1:h, colMeans(Obs,na.rm=TRUE), type="l",col=4)
legend("topleft",legend=c("NNAR","ARIMA","Obs"),col=2:4,lty=1)