#
#  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, "*")
}







#
#  R 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)
	#
	aux1 <- optim(par = 0.25, gr=NULL, arfima.loglik, method = "L-BFGS-B", lo = 0., up = 0.5, series = series, m = m)$par
	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)
}
#
#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, "*")
}












#
#SPLUS VERSION
#
arfima.whittle<-function(series)           
{    
 #    
 #   Computes Whittle MLE for ARFIMA(0,d,0) models    
 #    
  nlminb(start = 0.25, whittle.loglik, lo = 0., up = 0.499, series = series)$parameters  
}    


#
# R VERSION
#
arfima.whittle<-function(series)           
{    
 #    
 #   Computes Whittle MLE for ARFIMA(0,d,0) models    
 #    
 optim(par = 0.25, whittle.loglik,gr=NULL,method = "L-BFGS-B", lo = 0., up = 0.499, series = series)$par
}    



#
#SPLUS VERSION
#
arfima.whittle<-function(series)           
{    
 #    
 #   Computes Whittle MLE for ARFIMA(0,d,0) models    
 #    
  nlminb(start = 0.25, whittle.loglik, lo = 0., up = 0.499, series = series)$parameters  
}    

#
# AUXILIARY FUNCTIONS
# 

whittle.loglik<-function(x, series)  
{  
#
#   Whittle Loglikelihood
#
 series <- series - mean(series)  
 a <- fft(series)  
 a <- Mod(a)^2  
 n <- length(series)  
 a <- a/(2 * pi * n)  
 m <- n/2  
 #  
 #Fourier frequencies  
 #  
 w <- (2 * pi * (1:m))/n  
 #  
 #Spectral Density:  
 #  
 b <- fn.density(w, x)  
 #  
 #  Calculate sigma^2  
 #  
 sigma2 <- (2 * sum(a[1:m]/b))/n  
 #  
 #  Whittle Log-likelihood  
 #  
 loglik <- 2 * pi * (sum(log(b)) + sum(a[1:m]/b)/sigma2)  
 return(loglik/n + pi * log(sigma2))  
}  


fn.density<-function(x, d)
{
#
#  Spectral density of an ARFIMA(0,d,0) process
#
	a <- (2 * sin(x/2))^(-2 * d)
	a <- a/(2 * pi)
	return(a)
}