# Remove missing values from end points
na.ends <- function(x)
{
    tspx <- tsp(x)
    # Strip initial and final missing values
    nonmiss <- (1:length(x))[!is.na(x)]
    if(length(nonmiss)==0)
        stop("No non-missing data")
    j <- nonmiss[1]
    k <- nonmiss[length(nonmiss)]
    x <- x[j:k]
    if(!is.null(tspx))
        x <- ts(x,start=tspx[1]+(j-1)/tspx[3],f=tspx[3])
    return(x)
}

# Add back missing values at ends
# x is original series. y is the series with NAs removed at ends.
# returns y with the nas put back at beginning but not end.
undo.na.ends <- function(x,y)
{
    n <- length(x)
    nonmiss <- (1:length(x))[!is.na(x)]
    j <- nonmiss[1]
    k <- nonmiss[length(nonmiss)]
    if(j>1)
        y <- c(rep(NA,j-1),y)
    if(k<n)
        y <- c(y,rep(NA,n-k))
    tspx <- tsp(x)
    if(!is.null(tspx))
        tsp(y) <- tsp(x)
    return(y)
}

## Undifference
unfracdiff <- function(x,y,n,h,d)
{
    bin.c <- (-1)^(0:(n+h)) * choose(d, (0:(n+h)))
    b <- numeric(n)
    xnew <- LHS <- numeric(h)
    RHS <- cumsum(y)
    bs <- cumsum(bin.c[1:h])
    b <- bin.c[(1:n) + 1]
    xnew[1] <- RHS[1] <- y[1] - sum(b * rev(x))
    if(h>1)
    {
        for (k in 2:h) 
        {
            b <- b + bin.c[(1:n) + k]
            RHS[k] <- RHS[k] - sum(b * rev(x))
            LHS[k] <- sum(rev(xnew[1:(k-1)]) * bs[2:k])
            xnew[k] <- RHS[k] - LHS[k]
        }
    }
	tspx <- tsp(x)
    if(is.null(tspx))
        tspx <- c(1,length(x),1)
	return(ts(xnew,f=tspx[3],s=tspx[2]+1/tspx[3]))    
}

## Automatic ARFIMA modelling
## Will return Arima object if d < 0.01 to prevent estimation problems
arfima <- function(x, drange = c(0, 0.5), estim = c("mle","ls"),...)
{
    estim <- match.arg(estim)
    
    require(fracdiff)
    
    # Strip initial and final missing values
    xx <- na.ends(x)
    
    # Remove mean
    meanx <- mean(xx)
    xx <- xx - meanx
    
    # Choose differencing parameter with AR(2) proxy to handle correlations
    warn <- options(warn=-1)$warn
    fit <- fracdiff(xx,nar=2)
    options(warn=warn)
   
    # Choose p and q
    d <- fit$d
    y <- diffseries(xx, d=d)
    fit <- auto.arima(y, max.P=0, max.Q=0, stationary=TRUE, ...)
    
    # Refit model using fracdiff
    warn <- options(warn=-1)$warn
    fit <- fracdiff(xx, nar=fit$arma[1], nma=fit$arma[2])
    options(warn=warn)
    
    # Refine parameters with MLE
    if(estim=="mle")
    {
        y <- diffseries(xx, d=fit$d)
        p <- length(fit$ar)
        q <- length(fit$ma)
        fit2 <- try(Arima(y,order=c(p,0,q),include.mean=FALSE))
        if(class(fit2) != "try-error")
        {
            if(p>0)
                fit$ar <- fit2$coef[1:p]
            if(q>0)
                fit$ma <- -fit2$coef[p+(1:q)]
            fit$residuals <- fit2$residuals
        }
        else
            warning("MLE estimation failed. Returning LS estimates")
    }
    
    # Add things to model that will be needed by forecast.fracdiff
    fit$x <- xx + meanx
    fit$residuals <- undo.na.ends(x,residuals(fit))
    fit$x <- x
    fit$fitted <- fit$x - fit$residuals
    
    return(fit)
}

# Forecast the output of fracdiff() or arfima()

forecast.fracdiff <- function(object, h=10, level=c(80,95), fan=FALSE, ...) 
{
    # Extract data
    if (is.element("x", names(object))) 
        x <- object$x
    else 
        x <- object$x <- eval.parent(parse(text=as.character(object$call)[2]))
    
    xx <- na.ends(x)
    n <- length(xx)
    
    meanx <- mean(xx)
    xx <- xx - meanx
    
    # Construct ARMA part of model and forecast with it
    y <- diffseries(xx, d=object$d)
    fit <- arima(y, order=c(length(object$ar),0,length(object$ma)), include.mean=FALSE, fixed=c(object$ar,-object$ma))
    fcast.y <- forecast(fit, h=h, level=level)

    # Undifference
    fcast.x <- unfracdiff(xx,fcast.y$mean,n,h,object$d)
    
    # Binomial coefficient for expansion of d
    bin.c <- (-1)^(0:(n+h)) * choose(object$d,(0:(n+h)))

    #Cumulative forecasts of y and forecast of y
    # b <- numeric(n)
    # fcast.x <- LHS <- numeric(h)
    # RHS <- cumsum(fcast.y$mean)
    # bs <- cumsum(bin.c[1:h])
    # b <- bin.c[(1:n)+1]
    # fcast.x[1] <- RHS[1] <- fcast.y$mean[1] - sum(b*rev(xx))
    # if(h>1)
    # {
        # for (k in 2:h)
        # {
            # b <- b + bin.c[(1:n)+k]
            # RHS[k] <- RHS[k] - sum(b*rev(xx))
            # LHS[k] <- sum(rev(fcast.x[1:(k-1)]) * bs[2:k])
            # fcast.x[k] <- RHS[k] - LHS[k]
        # }
    # }
    
    # Extract stuff from ARMA model
    p <- fit$arma[1]
    q <- fit$arma[2]
    phi <- theta <- numeric(h)
    if(p > 0)
        phi[1:p] <- fit$coef[1:p]
    if(q > 0)
        theta[1:q] <- fit$coef[p+(1:q)]

    # Calculate psi weights
    new.phi <- psi <- numeric(h)
    psi[1] <- new.phi[1] <- 1
    if(h>1)
    {
        new.phi[2:h] <- -bin.c[2:h]
        for (i in 2:h) 
        {
            if(p>0)
                new.phi[i] <- sum(phi[1:(i-1)] * bin.c[(i-1):1]) - bin.c[i]
            psi[i] <- sum(new.phi[2:i] * rev(psi[1:(i-1)])) + theta[i-1]
        }
    }
    
    # Compute forecast variances
    fse <- sqrt(cumsum(psi^2) * fit$sigma2)
    
    # Compute prediction intervals
    if (fan) 
        level <- seq(51, 99, by = 3)
    else 
    {
        if (min(level) > 0 & max(level) < 1) 
            level <- 100 * level
        else if (min(level) < 0 | max(level) > 99.99) 
            stop("Confidence limit out of range")
    }
    nint <- length(level)
    upper <- lower <- matrix(NA, ncol = nint, nrow=h)
    for (i in 1:nint) 
    {
        qq <- qnorm(0.5 * (1 + level[i]/100))
        lower[, i] <- fcast.x - qq * fse
        upper[, i] <- fcast.x + qq * fse
    }
    colnames(lower) = colnames(upper) = paste(level, "%", sep = "")

    res <- undo.na.ends(x,residuals(fit))
    data.tsp <- tsp(x)
    if(is.null(data.tsp))
        data.tsp <- c(1,length(x),1)
    mean.fcast <- ts(fcast.x+meanx, f=data.tsp[3], s=data.tsp[2] + 1/data.tsp[3])
    lower <- ts(lower+meanx, f=data.tsp[3], s=data.tsp[2] + 1/data.tsp[3])
    upper <- ts(upper+meanx, f=data.tsp[3], s=data.tsp[2] + 1/data.tsp[3])
    method <- paste("ARFIMA(",p,",",round(object$d,2),",",q,")")
    return(structure(list(x=x, mean=mean.fcast, upper=upper, lower=lower, 
        level=level, method=method, xname=deparse(substitute(x)), model=object, 
        residuals=res, fitted=x-res), class="forecast"))
}

# Residuals from arfima() or fracdiff()

residuals.fracdiff <- function(object, ...)
{
    require(fracdiff)

    if(!is.null(object$residuals))   # Object produced by arfima()
        return(object$residuals)
    else                             # Object produced by fracdiff()
    {
        if (is.element("x", names(object))) 
            x <- object$x
        else 
            x <- eval.parent(parse(text=as.character(object$call)[2]))
        y <- diffseries(x - mean(x), d=object$d)
        fit <- arima(y, order=c(length(object$ar),0,length(object$ma)), include.mean=FALSE, fixed=c(object$ar,object$ma))
        return(residuals(fit))
    }
}

# Fitted values from arfima() or fracdiff()

fitted.fracdiff <- function(object, ...)
{
    if(!is.null(object$fitted))      # Object produced by arfima()
        return(object$fitted)
    else
    {
        if (is.element("x", names(object))) 
            x <- object$x
        else 
            x <- eval.parent(parse(text=as.character(object$call)[2]))
        return(x-residuals(object))
    }
}