############################################################
#
#  Fitting ARIMA and ARFIMA models
#
############################################################
require(stats)
require(graphics)
require(forecast)
library("arfima")

# ACF and PACF of ARFIMA data
# ---------------------------
lag = 200
k   = 2:lag
ds  = c(0.2,0.49)
par(mfrow=c(2,2))
for (d in ds){
  acf = c(d/(1-d),factorial(-d)/factorial(d-1)*k^(2*d-1))
  pacf = c(d/(1-d),d/(k-d))
  plot(c(1,k),acf,type="h",xlab="LAG")
  title(paste("ACF for d=",d,sep=""))
  plot(c(1,k),pacf,type="h",xlab="LAG",main="PACF")
}

# Simulating ARFIMA(1,2.3,0)
# ---------------------------
set.seed(6533)
sim <- arfima.sim(1000,model=list(phi=.2,dfrac=.3,dint=2))
par(mfrow=c(1,1))
plot(sim,type="b")
fit <- arfima(sim, order = c(1, 2, 0))
summary(fit)


# Revisiting the famous Annual Flow of Nile River
# -----------------------------------------------
?Nile
# Measurements of the annual flow of the river 
# Nile at Aswan (formerly Assuan), 1871–1970, 
# in 10^8 m^3, “with apparent changepoint near 
# 1898” (Cobb(1978), Table 1, p.249).

Nile
#Time Series:
#  Start = 1871 
#End = 1970 
#Frequency = 1 
#[1] 1120 1160  963 1210 1160 1160  813 1230 1370 1140  995  935
#[13] 1110  994 1020  960 1180  799  958 1140 1100 1210 1150 1250

par(mfrow=c(1,1))
plot(Nile)
abline(v=1898,col=2)

par(mfrow=c(1,2))
acf(Nile)
pacf(Nile)

fit1 = auto.arima(Nile)
fit2 = arfima(Nile) 

pred.arima=predict(fit1, 20)
pred.arfima=predict(fit2, 20)

par(mfrow=c(1,1))
plot(pred.arfima)
lines(pred.arima$pre,col=2,lwd=2)
lines(pred.arima$pre-2*pred.arima$se,lty=2,col=2,lwd=2)
lines(pred.arima$pre+2*pred.arima$se,lty=2,col=2,lwd=2)
lines(pred.arima$pre-pred.arima$se,lty=3,col=2,lwd=2)
lines(pred.arima$pre+pred.arima$se,lty=3,col=2,lwd=2)
legend("bottomleft",legend=c("ARIMA(1,1,1)","ARFIMA (d=0.364)"),
       col=c(2,4),lwd=2,bty="n")