klong<-function(d,ar,ma,n){

# generate a fractionally-differenced ARIMA(1,d,1) model given initial values

	x <- rts(arima.fracdiff.sim(model = list(d=d, ar=ar, ma=ma), n = n))

	# estimate the parameters in an ARIMA(1,d,1) model for the simulated series

	xfit<-arima.fracdiff(x, model = list(ar = ar, ma = ma))

	split.screen(c(2, 1))	#split screen into two parts
	split.screen(c(1, 3), screen = 2)	#split the bottom graph into three parts
	screen(1)
	ts.plot(x)
	title(main = "Fractionally Differenced ARMA Process")
	screen(3)
	hist(x)
	screen(4)
	xacf <- acf(x)
	screen(5)
	acf(x,type="partial")
	close.screen(all = T)	#lag.plot(x,layout=c(2,5),lags=10)
	spectrum(x,method="ar")
	
	return(x,xfit)

}

klongfr<-function(x,fit,nforecast){
	
		ar<-fit$ar
		ma<-fit$ma
		d<-fit$d
		
		# This function is written for an FARIMA(1,d,q), AR(1) is hardwired but can be modified
	
		# Now, let's find the ARMA(infinity,q) representation so that we can use
		# the forecasting tools for ARMA models
	
		n<-length(x)		# Number of observations

		kj<-seq(1,100)
		pij<-gamma(kj-d)/(gamma(kj+1)*gamma(-d))   #Coefficients of binomial expansion of fractional
															# difference operator
		kar<-c(ar-pij[1],ar*pij[1:99]-pij[2:100])  # Multiplied by AR(1) operator

		# Look at coefficients to see when they are very close to zero
		# Suppose that point is lag 10
		
		approxmodel<-list(ar=kar[1:10],ma=ma)
	
		xpfore<-c(x,rep(NA,nforecast))
	 	forecast<-arima.filt(xpfore,model=approxmodel)
	    pred<-forecast$pred[(n+1):(n+nforecast)]
		
		# Plot forecast
		xpfore<-c(x,pred)
		xpflb<-c(x,pred-1.96*sqrt(forecast$var.pred[(n+1):(n+nforecast)]))
		xpfub<-c(x,pred+1.96*sqrt(forecast$var.pred[(n+1):(n+nforecast)]))
		nf<-n-50
		nl<-n+nforecast
		tsplot(xpfore[nf:nl],xpflb[nf:nl],xpfub[nf:nl])
		title(main="Last 50 points plus forecasted values with 95% PI")
		
		#Examining the diagnostics
	    resid<- x[11:n]-forecast$pred[11:n]  # First 10 are NA
		stresid<-resid/sqrt(forecast$var.pred[11:n])
		par(mfrow=c(2,2))
		tsplot(stresid)
		boxplot(stresid)
		acf(stresid)
		acf(stresid,type="partial")
		par(mfrow=c(1,1))
			
		return(forecast,resid,stresid)

}

klonga<-klong(.4,.8,.2,1000)
klongfra<-klongfr(klonga$x,klonga$xfit$model,10)

klongb<-klong(.1,.9,0,1000)
klongfrb<-klongfr(klongb$x,klongb$xfit$model,10)