#
#  SPLUS VERSION
#
arfima.mle<-function(series, m = 30)
{
	#
	# Computes MLE for an ARFIMA(0,d,0) model, missing values allowed,
	# (the "m" first values must be observed)
	#
	#  output:
	#                   d             : estimate of "d"
	#	std          : standard deviation of the estimate of d
	#	student   : t-test value for d
	#	missing   : number of missing values in the series
	#	noise.std : white noise standard deviation
	#
	aux1 <- nlminb(start = 0.25, arfima.loglik, lo = 0., up = 0.5, series = series, m = m)$parameters
	aux2 <- hess(arfima.loglik, aux1, series = series)
	aux2 <- sqrt(2./aux2)
	aux3 <- aux1/aux2
	aux4 <- length(series[series == "NA"])
	order <- c(m, 0., 0.)
	ar <-  - psi1( - aux1, m)
	model <- list(order, ar, 0.)
	names(model) <- c("order", "ar", "ndiff")
	aux5 <- arima.filt(series - mean(series, na.rm = T), model = model)$sigma2
	aux5 <- sqrt(aux5)
	aux <- list(aux1, aux2, aux3, aux4, aux5)
	names(aux) <- c("d", "std", "student", "missing", "noise.std")
	return(aux)
}
#
# AUXILIARY FUNCTIONS
#
arfima.loglik<-function(x, series, m = 30)
{
	order <- c(m, 0., 0.)
	ar <-  - psi1( - x, m)
	model <- list(order, ar, 0.)
	names(model) <- c("order", "ar", "ndiff")
	y <- arima.filt(series - mean(series, na.rm = T), model = model)$loglik
	return(y)
}
psi1<-function(d, n)
{
	psi <- c(1:n)
	psi[1] <- d
	if(n > 1) {
		for(i in 2:n)
			psi[i] <- ((i - 1 + d) * psi[i - 1])/i
	}
	return(psi)
}
#
#Available from the website: http://www.public.iastate.edu/~pdixon/stat534/Splus/hess.do
#
hess<-function(f, x0, ...)
{
	# compute numerical approx. to Hessian of f, evaluated at x0
	#  usually need to pass additional parameters (e.g. data)
	# N.B. this uses no numerical sophistication,
	#   so expect 3-4 figure accuracy, at best
	# input: f: name of function that defines log likelihood (or negative of it)
	#        x0: scalar or vector of parameters that give the point at which
	#             you want the hessian estimated (usually will be the mle )
	# output: matrix of 2nd derivatives.  Size is n x n, where n = length of x0
	n <- length(x0)
	grad <- rep(0., n)
	mdelta <- hess <- matrix(0., nrow = n, ncol = n)
	delta <- 0.0001 * (sign(x0) + (x0 == 0.)) * pmax(abs(x0), 0.01)
	diag(mdelta) <- delta
	f0 <- f(x0, ...)
	for(i in 1.:n) {
		grad[i] <- f(x0 + mdelta[, i], ...)
	}
	for(i in 1.:n) {
		for(j in i:n) {
			hess[i, j] <- f(x0 + mdelta[, i] + mdelta[, j], ...) - grad[i] - grad[j]
			hess[j, i] <- hess[i, j]
		}
	}
	(hess + f0)/outer(delta, delta, "*")
}