###########################################################
###   Testing for the integration order, d0 and ds      ###
###########################################################

### MAIN PROGRAM: godfrey.d0ds.test ##################################################
# Testing d0 and ds by Godfrey's Auxiliary regression method, e.g.,
# H0: d0=d0' and ds=ds' vs H1: d0 != d0' or ds != ds'
# returns p-value of the test statistics.
# USAGE: godfrey.d0ds.test(x, model, xreg = 0, H0.d = c(0,0), sarma.start=NULL, 
#		 output.summary=T, scale=1, control=NULL, lower = -Inf, upper = Inf)
# REQUIRED ARGUMENTS: 
# x 		a univariate time serires or a vector
# model		specify the order of the SARMA model like sarfima.css
# OPTIONAL ARGUMENTS: 
# xreg 		estimator of the mean of x. default is NULL
# H0.d	 	a null hypothesis for d0 and ds. default is d0=0, ds=0.
# sarma.start	initial value of SARMA parameters
# output.summary
#		If TRUE, summary of outputs is displayed. 
# scale, control, lower, upper
#	 	optional arguements for nlminb
# VALUE: 
# teststat$H0.d	same as input.
# teststat$chisq, teststat$pvalue
#		value and p-value of chisquared test statistics.
# model, xreg	same as input.
# estimates$ar, estimates$ma, estimates$sar, estimates$sma, estimates$sigma2
# 		estimates of SARFIMA parameters.
# converged 	if the optimizer, nlminb, has apparently successfully converged to a minimum, 
#		then converged is TRUE; otherwise converged is FALSE. 
# conv.type 	a character string describing the type of convergence by nlminb.
# stinv		If estimates of SARMA prameters satisfy stationary and invertible conditions, 
#		then stinv is TRUE; otherwise stinv is FALSE. 
# n.used        the number of observations used to compute the CSS function. 
# np		the number of parameters. 
######################################################################################


godfrey.d0ds.test <- function(x, model, xreg = 0, H0.d = c(0,0), sarma.start=NULL, output.summary=T, scale=1, control=NULL, lower = -Inf, upper = Inf){
	### assignement ###
	arima.order <- model[[1]]$order
	seasonal.order <- model[[2]]$order
	p0 <- arima.order[1]
	q0 <- arima.order[3]
	ps <- seasonal.order[1]
	qs <- seasonal.order[3]
	ss <- model[[2]]$period
	d0 <- H0.d[1]
	ds <- H0.d[2]
	m0 <- arima.order[2]
	ms <- seasonal.order[2]
	x <- as.vector(x)
	x <- if(m0>0) diff(x, lag=1, differences=m0) else x
	x <- if(ms>0) diff(x, lag=ss, differences=ms) else x
	x <- x - xreg
	n <- length(x)
	np <- p0 + ps + q0 + qs +1
	if(is.one(np)) {
	est0list <- list(d0 = NULL, ds = NULL, ar = NULL, ma = NULL, sar = NULL, sma = NULL)
	fit <- list(estimates = est0list, converged = TRUE, conv.type = NULL, hessian = NULL)
	} else {
	### estimate SARMA parameters ###
	init.sarma <- initial.sarma(d0, ds, p0, q0, ps, qs, ss, n, x, sarma.start=sarma.start)
	fit <- sarfima.nlse.nlminb(x, H0.d, init.sarma, ss, d0.est = F, ds.est = F, scale=scale, control=control, lower = lower, upper = upper)
	}
	resar <- fit$estimates$ar
	resma <- fit$estimates$ma
	ressar <- fit$estimates$sar
	ressma <- fit$estimates$sma
	n.ar <- nullto0(resar)
	n.ma <- nullto0(resma)
	n.sar <- nullto0(ressar)
	n.sma <- nullto0(ressma)
	sigma2 <- obj.fun.SARFIMA(d0, ds, n.ar, n.ma, n.sar, n.sma, ss, n, x)/n
	### check the stationarity/invertibily of the fitted coeffs ###
	sarma.sta.inv <- SARMA.test(n.ar, n.ma, n.sar, n.sma, p0, q0, ps, qs, ss)
	if(!sarma.sta.inv)
	warning("Warning: estimated SARMA coefficients are non-stationary or non-invertible")
	### calculate Godfrey's test statistics ###
	e <- resid.SARFIMA.gen(d0, ds, n.ar, n.ma, n.sar, n.sma, ss, n, x)
	matX <- matrix.deriv.e(n.ar, n.ma, n.sar, n.sma, c(p0, q0, ps, qs), ss, e, d0.est = T, ds.est = T)
	Xe <- t(matX) %*% e
	invXX <- solve(t(matX) %*% matX)
	chisq <- c((t(Xe) %*% invXX %*% Xe)/sigma2)
	pvalue <- 1 - pchisq(chisq, 2)
	### output ###
	estlist <- list(ar = resar, ma = resma, sar = ressar, sma = ressma, sigma2 = sigma2)
	result <- list(teststat = list(H0.d = H0.d, chisq = chisq, pvalue = pvalue), model=model, xreg=xreg, estimates = estlist, converged = fit$converged, conv.type = fit$conv.type, stinv = (sarma.sta.inv),  n.used= n, np=np )
	if(output.summary){
	godfrey.d0ds.summary(result)
	invisible(result)
	}
	else result
}



### cat function 

godfrey.d0ds.summary <- function(z)
{
p0 <- z$model[[1]]$order[1]
m0 <- z$model[[1]]$order[2]
q0 <- z$model[[1]]$order[3]
ps <- z$model[[2]]$order[1]
ms <- z$model[[2]]$order[2]
qs <- z$model[[2]]$order[3]
ss <- z$model[[2]]$period
sarimaorder <- data.frame(list(p0=p0, m0=m0, q0=q0, ps=ps, ms=ms, qs=qs, period=ss, 
xreg0=is.zero( sum((z$xreg)^2) ), n.used=z$n.used, np=z$np, convd =z$converged, stinv=z$stinv), row.names="")
cat("\n","# hypothesis #","\n")
cat(" H0: ( d0 , ds ) = (", z$teststat$H0.d[1], ",", z$teststat$H0.d[2],")", 
"  vs  " ,
"H1: ( d0 , ds ) != (", z$teststat$H0.d[1], ",", z$teststat$H0.d[2],")", "\n")
cat("\n","# test statistics #","\n")
print(data.frame(list(ChiSq=z$teststat$chisq, pvalue = z$teststat$pvalue, df=2), row.names=""))
cat("\n","# estimated sarima model #","\n")
print(sarimaorder)
cat("\n","# estimates of nuisance parameters #", "\n")
print(unlist(z$estimates))
cat("\n")
invisible(z)
}








### MAIN PROGRAM: sarfima.d0.test ######################################################
# Testing d0 by Score test statistics, e.g.,
# H0: d0=d0'vs H1: d0 != d0'
# returns p-value of the test statistics.
# USAGE: sarfima.d0.test(x, model, xreg = 0, H0.d = 0, ds.est=F, sarma.start=NULL, 
#		output.summary=T, scale=1, control=NULL, lower = -Inf, upper = Inf)
# REQUIRED ARGUMENTS: 
# x 		a univariate time serires or a vector
# model		specify the order of the SARMA model like sarfima.css
# OPTIONAL ARGUMENTS: 
# xreg 		estimator of the mean of x. default is NULL
# H0.d	 	a null hypothesis for d0. default is d0=0.
# ds.est 	ds is estimeted if TRUE and ds is not estimeted if FALSE.  
# sarma.start	initial value of SARMA parameters
# output.summary
#		If TRUE, summary of outputs is displayed. 
# scale, control, lower, upper
#	 	optional arguements for nlminb
# VALUE: 
# teststat$H0.d	same as input.
# teststat$teststat
#		value of test statistics.
# teststat$pvaluetwo, teststat$pvalueright, teststat$pvalueleft
#		p-value of two-sided, right-sided, and left-sided of tests. 
# model, ds.est, xreg	
#		same as input.
# estimates	estimates of SARFIMA parameters.
# converged	if the optimizer, nlminb, has apparently successfully converged to a minimum, 
#		then converged is TRUE; otherwise converged is FALSE. 
# conv.type 	a character string describing the type of convergence by nlminb.
# stinv		If estimates of SARMA prameters satisfy stationary and invertible conditions, 
#		then stinv is TRUE; otherwise stinv is FALSE. 
# n.used        the number of observations used to compute the CSS function. 
# np		the number of parameters. 
######################################################################################


sarfima.d0.test <- function(x, model, xreg = 0, H0.d = 0, ds.est=F, sarma.start=NULL, output.summary=T, scale=1, control=NULL, lower = -Inf, upper = Inf){
	### assignement ###
	arima.order <- model[[1]]$order
	seasonal.order <- model[[2]]$order
	p0 <- arima.order[1]
	q0 <- arima.order[3]
	ps <- seasonal.order[1]
	qs <- seasonal.order[3]
	ss <- model[[2]]$period
	m0 <- arima.order[2]
	ms <- seasonal.order[2]
	x <- as.vector(x)
	x <- if(m0>0) diff(x, lag=1, differences=m0) else x
	x <- if(ms>0) diff(x, lag=ss, differences=ms) else x
	x <- x - xreg
	n <- length(x)
	np <- p0 + ps + q0 + qs + ds.est + 1
	if(is.one(np)) {
	est0list <- list(d0 = NULL, ds = NULL, ar = NULL, ma = NULL, sar = NULL, sma = NULL)
	fit <- list(estimates = est0list, converged = TRUE, conv.type = NULL, hessian = NULL)
	} else {
	### estimate SARMA parameters ###
	init.sarma <- initial.sarma(H0.d, 0, p0, q0, ps, qs, ss, n, x, sarma.start=sarma.start)
	fit <- sarfima.nlse.nlminb(x, c(H0.d, 0), init.sarma, ss, d0.est = F, ds.est = ds.est, scale=scale, control=control, lower = lower, upper = upper)
	}
	resds <- fit$estimates$ds
	resar <- fit$estimates$ar
	resma <- fit$estimates$ma
	ressar <- fit$estimates$sar
	ressma <- fit$estimates$sma
	n.ds <- nullto0(resds)
	n.ar <- nullto0(resar)
	n.ma <- nullto0(resma)
	n.sar <- nullto0(ressar)
	n.sma <- nullto0(ressma)
	sigma2 <- obj.fun.SARFIMA(H0.d, n.ds, n.ar, n.ma, n.sar, n.sma, ss, n, x)/n
	### check the stationarity/invertibily of the fitted coeffs ###
	if(abs(n.ds)>0.5) 
	warning("Warning: estimated ds is non-stationary or non-invertible")
	sarma.sta.inv <- SARMA.test(n.ar, n.ma, n.sar, n.sma, p0, q0, ps, qs, ss)
	if(!sarma.sta.inv)
	warning("Warning: estimated SARMA coefficients are non-stationary or non-invertible")
	### calculate Tanaka's test statistics ###
	e <- resid.SARFIMA.gen(H0.d, n.ds, n.ar, n.ma, n.sar, n.sma, ss, n, x)
	racf <- meanzero.acf(e, n-1)[-1]
	if(is.one(np)) std.err <- sqrt(pi^2/6) 
	else {
	var.coef <- sarfima.css.var(resar, resma, ressar, ressma, c(p0,q0,ps,qs), ss, 1, d0.est = F, ds.est = ds.est)
	std.err <- 1/sqrt(var.coef[1,1])
	}
	teststat <- sqrt(n*(n+2))*sum( 	racf/( (1:(n-1)) * sqrt((n-1):1) ) )/std.err
	pvaluetwo <- 1 - pchisq(teststat^2, 1)
	pvalueright <- 1 - pnorm(teststat)
	pvalueleft <- pnorm(teststat)
	### output ###
	estlist <- list(ds = resds, ar = resar, ma = resma, sar = ressar, sma = ressma, sigma2 = sigma2)
	result <- list(teststat = list(H0.d = H0.d, 
	teststat = teststat, pvaluetwo = pvaluetwo, pvalueright = pvalueright, pvalueleft = pvalueleft),
	 model=model, ds.est=ds.est, xreg=xreg, estimates = estlist, converged = fit$converged, 
	conv.type = fit$conv.type, stinv = (sarma.sta.inv),  n.used= n, np=np )
	if(output.summary){
	sarfima.d0.summary(result)
	invisible(result)
	}
	else result
}


### cat function 

sarfima.d0.summary <- function(z)
{
p0 <- z$model[[1]]$order[1]
m0 <- z$model[[1]]$order[2]
q0 <- z$model[[1]]$order[3]
ps <- z$model[[2]]$order[1]
ms <- z$model[[2]]$order[2]
qs <- z$model[[2]]$order[3]
ss <- z$model[[2]]$period
ds.est <- z$ds.est
sarimaorder <- data.frame(list(p0=p0, m0=m0, q0=q0, ps=ps, ms=ms, qs=qs, period=ss,
xreg0=is.zero( sum((z$xreg)^2) ), n.used=z$n.used, np=z$np, convd =z$converged, stinv=z$stinv), row.names="")
cat("\n","# hypothesis #","\n")
cat(" H0: d0 = ", z$teststat$H0.d, "  vs  " , " H0: d0 != ", z$teststat$H0.d, "\n")
cat("\n","# test statistics #","\n")
print(data.frame(list(teststat=z$teststat$teststat, pvaluetwo = z$teststat$pvaluetwo, pvalueleft = z$teststat$pvalueleft, pvalueright = z$teststat$pvalueright), row.names=""))
cat("\n", "# estimated model #",  "ds.est=", ds.est, "\n")
print(sarimaorder)
cat("\n","# estimates of nuisance parameters #", "\n")
print(unlist(z$estimates))
cat("\n")
invisible(z)
}






### MAIN PROGRAM: sarfima.ds.test ######################################################
# Testing ds by Score test statistics, e.g.,
# H0: ds=ds'vs H1: ds != ds'
# returns p-value of the test statistics.
# USAGE: sarfima.ds.test(x, model, xreg = 0, H0.d = 0, d0.est=F, sarma.start=NULL, 
# 		output.summary=T, scale=1, control=NULL, lower = -Inf, upper = Inf)
# REQUIRED ARGUMENTS: 
# x 		a univariate time serires or a vector
# model		specify the order of the SARMA model like sarfima.css
# OPTIONAL ARGUMENTS: 
# xreg 		estimator of the mean of x. default is NULL
# H0.d	 	a null hypothesis for ds. default is ds=0.
# d0.est 	d0 is estimeted if TRUE and d0 is not estimeted if FALSE.  
# sarma.start	initial value of SARMA parameters
# output.summary
#		If TRUE, summary of outputs is displayed. 
# scale, control, lower, upper
#	 	optional arguements for nlminb
# VALUE: 
# teststat$H0.d	same as input.
# teststat$teststat
#		value of test statistics.
# teststat$pvaluetwo, teststat$pvalueright, teststat$pvalueleft
#		p-value of two-sided, right-sided, and left-sided of tests. 
# model, d0.est, xreg	
#		same as input.
# estimates	estimates of SARFIMA parameters.
# converged	if the optimizer, nlminb, has apparently successfully converged to a minimum, 
#		then converged is TRUE; otherwise converged is FALSE. 
# conv.type 	a character string describing the type of convergence by nlminb.
# stinv		If estimates of SARMA prameters satisfy stationary and invertible conditions, 
#		then stinv is TRUE; otherwise stinv is FALSE. 
# n.used        the number of observations used to compute the CSS function. 
# np		the number of parameters. 
######################################################################################


sarfima.ds.test <- function(x, model, xreg = 0, H0.d = 0, d0.est=F, sarma.start=NULL, output.summary=T, scale=1, control=NULL, lower = -Inf, upper = Inf){
	### assignement ###
	arima.order <- model[[1]]$order
	seasonal.order <- model[[2]]$order
	p0 <- arima.order[1]
	q0 <- arima.order[3]
	ps <- seasonal.order[1]
	qs <- seasonal.order[3]
	ss <- model[[2]]$period
	m0 <- arima.order[2]
	ms <- seasonal.order[2]
	x <- as.vector(x)
	x <- if(m0>0) diff(x, lag=1, differences=m0) else x
	x <- if(ms>0) diff(x, lag=ss, differences=ms) else x
	x <- x - xreg
	n <- length(x)
	np <- p0 + ps + q0 + qs + d0.est + 1
	if(is.one(np)) {
	est0list <- list(d0 = NULL, ds = NULL, ar = NULL, ma = NULL, sar = NULL, sma = NULL)
	fit <- list(estimates = est0list, converged = TRUE, conv.type = NULL, hessian = NULL)
	} else {
	### estimate SARMA parameters ###
	init.sarma <- initial.sarma(0, H0.d, p0, q0, ps, qs, ss, n, x, sarma.start=sarma.start)
	fit <- sarfima.nlse.nlminb(x, c(0, H0.d), init.sarma, ss, d0.est = d0.est, ds.est = F, scale=scale, control=control, lower = lower, upper = upper)
	}
	resd0 <- fit$estimates$d0
	resar <- fit$estimates$ar
	resma <- fit$estimates$ma
	ressar <- fit$estimates$sar
	ressma <- fit$estimates$sma
	n.d0 <- nullto0(resd0)
	n.ar <- nullto0(resar)
	n.ma <- nullto0(resma)
	n.sar <- nullto0(ressar)
	n.sma <- nullto0(ressma)
	sigma2 <- obj.fun.SARFIMA(n.d0, H0.d, n.ar, n.ma, n.sar, n.sma, ss, n, x)/n
	### check the stationarity/invertibily of the fitted coeffs ###
	if(abs(n.d0)>0.5) 
	warning("Warning: estimated d0 is non-stationary or non-invertible")
	sarma.sta.inv <- SARMA.test(n.ar, n.ma, n.sar, n.sma, p0, q0, ps, qs, ss)
	if(!sarma.sta.inv)
	warning("Warning: estimated SARMA coefficients are non-stationary or non-invertible")
	### calculate Tanaka's test statistics ###
	e <- resid.SARFIMA.gen(n.d0, H0.d, n.ar, n.ma, n.sar, n.sma, ss, n, x)
	racf <- meanzero.acf(e, n-1)[-1]
	if(is.one(np)) std.err <- sqrt(pi^2/6) 
	else {
	var.coef <- sarfima.css.var(resar, resma, ressar, ressma, c(p0,q0,ps,qs), ss, 1, d0.est = d0.est, d0.est = F)
	std.err <- 1/sqrt(var.coef[2,2])
	}
	nn <- floor((n-1)/ss)
	teststat <- sqrt(n*(n+2))*sum( 	racf[c(1:nn)*ss]/( (1:nn) * sqrt(n-c(1:nn)*ss) )  )/std.err
	pvaluetwo <- 1 - pchisq(teststat^2, 1)
	pvalueright <- 1 - pnorm(teststat)
	pvalueleft <- pnorm(teststat)
	### output ###
	estlist <- list(d0 = resd0, ar = resar, ma = resma, sar = ressar, sma = ressma, sigma2 = sigma2)
	result <- list(teststat = list(H0.d = H0.d, 
	teststat = teststat, pvaluetwo = pvaluetwo, pvalueright = pvalueright, pvalueleft = pvalueleft),
	 model=model, d0.est=d0.est, xreg=xreg, estimates = estlist, converged = fit$converged, 
	conv.type = fit$conv.type, stinv = (sarma.sta.inv),  n.used= n, np=np )
	if(output.summary){
	sarfima.ds.summary(result)
	invisible(result)
	}
	else result
}


### cat function

sarfima.ds.summary <- function(z)
{
p0 <- z$model[[1]]$order[1]
m0 <- z$model[[1]]$order[2]
q0 <- z$model[[1]]$order[3]
ps <- z$model[[2]]$order[1]
ms <- z$model[[2]]$order[2]
qs <- z$model[[2]]$order[3]
ss <- z$model[[2]]$period
d0.est <- z$d0.est
sarimaorder <- data.frame(list(p0=p0, m0=m0, q0=q0, ps=ps, ms=ms, qs=qs, period=ss,
xreg0=is.zero( sum((z$xreg)^2) ), n.used=z$n.used, np=z$np, convd =z$converged, stinv=z$stinv), row.names="")
cat("\n","# hypothesis #","\n")
cat(" H0: ds = ", z$teststat$H0.d, "  vs  " , " H0: ds != ", z$teststat$H0.d, "\n")
cat("\n","# test statistics #","\n")
print(data.frame(list(teststat=z$teststat$teststat, pvaluetwo = z$teststat$pvaluetwo, pvalueleft = z$teststat$pvalueleft, pvalueright = z$teststat$pvalueright), row.names=""))
cat("\n", "# estimated model #",  "d0.est=", d0.est, "\n")
print(sarimaorder)
cat("\n","# estimates of nuisance parameters #", "\n")
print(unlist(z$estimates))
cat("\n")
invisible(z)
}







### example 1 ################################################################################
# data <- arima.fracdiff.sim(model = list(d = 0.3, ar = -0.2, ma = 0.4), n = 200) 
# model1 <- list(list(order=c(1,0,1)), list(order=c(0,0,0), period=12))
# godfrey.d0ds.test(data, model1)
# result <- godfrey.d0ds.test(data, model=model1, H0.d = c(0.3,0) )
# godfrey.d0ds.summary(result)
# sarfima.d0.test(data, model1)
# sarfima.ds.test(data, model1)
# sarfima.ds.test(data, model1, d0.est=T)
##############################################################################################

### example 2 ################################################################################
# data <- cumsum( arima.fracdiff.sim(model = list(d = 0.3, ar = -0.2, ma = 0.4), n = 500) + 1)[301:500]
# model2 <- list(list(order=c(1,1,1)), list(order=c(0,0,0), period=12))
# godfrey.d0ds.test(data, model2, xreg=mean(diff(data)))
# godfrey.d0ds.test(data, model2, xreg=mean(diff(data)),  H0.d = c(0.3,0) )
# sarfima.d0.test(data, model2, xreg=mean(diff(data)))
# sarfima.ds.test(data, model2, xreg=mean(diff(data)))
# sarfima.ds.test(data, model2, d0.est=T, xreg=mean(diff(data)))
##############################################################################################