####################################################################
###  sarfima.css.ms: Select best models in terms of AIC and BIC  ###
###  Note: this program works S-PLUS ver 6.2 because this        ###
###  uses new S-PLUS function "try".                            ###
####################################################################


### generate all lists of SARFIMA models ###

sarfima.allpatterns <- function( ss, diff.order=c(0,0), max.sarma=c(1,1,1,1), max.np.sarfima=NULL, d0.est=T, ds.est=T )
{
	assign("m0", diff.order[1], frame =1)
	assign("ms", diff.order[2], frame =1)
	assign("ss", ss, frame =1)
	modellist <- expand.grid(p0=c(0:max.sarma[1]), q0=c(0:max.sarma[2]), ps=c(0:max.sarma[3]), qs=c(0:max.sarma[4]), d0.est=c(0:d0.est), ds.est=c(0:ds.est))
	if(!is.null(max.np.sarfima)) {
	npvec <- as.matrix(modellist) %*% rep(1,6)
	modellist <- modellist[row(npvec)[npvec <= max.np.sarfima], ]
	}
	assign("sarfima.grid", data.frame(modellist, ID=c(1:nrow(modellist))), frame = 1)
	lappfun <- function(rowseq)
	{
	est.order <- sarfima.grid[rowseq,]
	list(
	list(order=c(est.order$p0, m0, est.order$q0)), 
	list(order=c(est.order$ps, ms, est.order$qs), period=ss), 
	list(d0.est=is.one(est.order$d0.est), ds.est=is.one(est.order$ds.est), 
	np=est.order$p0+est.order$q0+est.order$ps+est.order$qs+est.order$d0.est+est.order$ds.est+1, 
	ID=est.order$ID))
	}
	lapply(1:nrow(sarfima.grid), lappfun)
}


### cat function ###

sarfima.ms.summary <- function(z)
{
	assign("z", z, frame = 1)
	if(is.null(z$detail)) {
		cat("# short summary #", "\n")
		print(z$result)
		cat("\n", "# number of models #", z$n.model, "\n", "\n")
		invisible(z)
	}
	else {
		cat("# summary of AIC model selection  #", "\n")
		n <- nrow(z$result$AIC)
		zlist <- list(ID = z$result$AIC[, 1], AIC = z$result$AIC[, 2], BIC = z$result$AIC[, 3], 
		ChiSq = z$result$AIC[, 4], RACF = z$result$AIC[, 5], SARMA = z$result$AIC[, 6], 
		d0ds = z$result$AIC[, 7], convd = z$result$AIC[, 8], p0 = z$result$AIC[, 9], 
		m0 = z$result$AIC[, 10], q0 = z$result$AIC[, 11], ps = z$result$AIC[, 12], 
		ms = z$result$AIC[, 13], qs = z$result$AIC[, 14], d0.est = z$result$AIC[, 15], 
		ds.est = z$result$AIC[, 16], np = z$result$AIC[, 17])
		maxd0 <- if(length(zlist$d0.est[zlist$d0.est == T]) > 0) 1 else 0
		maxds <- if(length(zlist$ds.est[zlist$ds.est == T]) > 0) 1 else 0
		maxp0 <- max(zlist$p0)
		maxq0 <- max(zlist$q0)
		maxps <- max(zlist$ps)
		maxqs <- max(zlist$qs)
		maxnp <- max(zlist$np)
		if(is.one(maxd0)) {
			d0vec <- sapply(1:n, function(j)
			round(nullto0(z$detail$AIC[[j]]$estimates$d0, NA), 3))
			d0name <- "d0"
		}
		else {
			d0vec <- NULL
			d0name <- NULL
		}
		if(is.one(maxds)) {
			dsvec <- sapply(1:n, function(j)
			round(nullto0(z$detail$AIC[[j]]$estimates$ds, NA), 3))
			dsname <- "ds"
		}
		else {
			dsvec <- NULL
			dsname <- NULL
		}
		if(maxp0 > 0) {
			assign("seq.ar", expand.grid(1:n, 1:maxp0), frame = 1)
			seq.list <- sapply(1:nrow(seq.ar), function(j)
			list(seq.ar[j,  ]))
			arsapfun <- function(elm)
			round(nullto0(z$detail$AIC[[elm$Var1]]$estimates$ar[elm$Var2], NA), 3)
			ARmat <- matrix(sapply(seq.list, arsapfun), n, maxp0)
			ARname <- paste("ar", 1:maxp0, sep = "")
		}
		else {
			ARmat <- NULL
			ARname <- NULL
		}
		if(maxq0 > 0) {
			assign("seq.ma", expand.grid(1:n, 1:maxq0), frame = 1)
			seq.list <- sapply(1:nrow(seq.ma), function(j)
			list(seq.ma[j,  ]))
			masapfun <- function(elm)
			round(nullto0(z$detail$AIC[[elm$Var1]]$estimates$ma[elm$Var2], NA), 3)
			MAmat <- matrix(sapply(seq.list, masapfun), n, maxq0)
			MAname <- paste("ma", 1:maxq0, sep = "")
		}
		else {
			MAmat <- NULL
			MAname <- NULL
		}
		if(maxps > 0) {
			assign("seq.sar", expand.grid(1:n, 1:maxps), frame = 1)
			seq.list <- sapply(1:nrow(seq.sar), function(j)
			list(seq.sar[j,  ]))
			sarsapfun <- function(elm)
			round(nullto0(z$detail$AIC[[elm$Var1]]$estimates$sar[elm$Var2], NA), 3)
			SARmat <- matrix(sapply(seq.list, sarsapfun), n, maxps)
			SARname <- paste("sar", 1:maxps, sep = "")
		}
		else {
			SARmat <- NULL
			SARname <- NULL
		}
		if(maxqs > 0) {
			assign("seq.sma", expand.grid(1:n, 1:maxqs), frame = 1)
			seq.list <- sapply(1:nrow(seq.sma), function(j)
			list(seq.sma[j,  ]))
			smasapfun <- function(elm)
			round(nullto0(z$detail$AIC[[elm$Var1]]$estimates$sma[elm$Var2], NA), 3)
			SMAmat <- matrix(sapply(seq.list, smasapfun), n, maxqs)
			SMAname <- paste("sma", 1:maxqs, sep = "")
		}
		else {
			SMAmat <- NULL
			SMAname <- NULL
		}
		sigma2vec <- sapply(1:n, function(j)
		signif(nullto0(z$detail$AIC[[j]]$estimates$sigma2, NA), 5))
		estmat <- matrix(c(d0vec, dsvec, as.vector(ARmat), as.vector(MAmat), as.vector(SARmat), as.vector(SMAmat), sigma2vec), nrow = n)
		dimnames(estmat) <- list(NULL, c(d0name, dsname, ARname, MAname, SARname, SMAname, "sigma2"))
		resultAIC <- data.frame(data.frame(zlist), estmat)
		print(resultAIC)
		invisible(resultAIC)
		cat("\n", "# summary of BIC model selection  #", "\n")
		n <- nrow(z$result$BIC)
		zlist <- list(ID = z$result$BIC[, 1], BIC = z$result$BIC[, 2], AIC = z$result$BIC[, 3],
		ChiSq = z$result$BIC[, 4], RACF = z$result$BIC[, 5], SARMA = z$result$BIC[, 6], 
		d0ds = z$result$BIC[, 7], convd = z$result$BIC[, 8], p0 = z$result$BIC[, 9], 
		m0 = z$result$BIC[, 10], q0 = z$result$BIC[, 11], ps = z$result$BIC[, 12], 
		ms = z$result$BIC[, 13], qs = z$result$BIC[, 14], d0.est = z$result$BIC[, 15], 
		ds.est = z$result$BIC[, 16], np = z$result$BIC[, 17])
		maxd0 <- if(length(zlist$d0.est[zlist$d0.est == T]) > 0) 1 else 0
		maxds <- if(length(zlist$ds.est[zlist$ds.est == T]) > 0) 1 else 0
		maxp0 <- max(zlist$p0)
		maxq0 <- max(zlist$q0)
		maxps <- max(zlist$ps)
		maxqs <- max(zlist$qs)
		maxnp <- max(zlist$np)
		if(is.one(maxd0)) {
			d0vec <- sapply(1:n, function(j)
			round(nullto0(z$detail$BIC[[j]]$estimates$d0, NA), 3))
			d0name <- "d0"
		}
		else {
			d0vec <- NULL
			d0name <- NULL
		}
		if(is.one(maxds)) {
			dsvec <- sapply(1:n, function(j)
			round(nullto0(z$detail$BIC[[j]]$estimates$ds, NA), 3))
			dsname <- "ds"
		}
		else {
			dsvec <- NULL
			dsname <- NULL
		}
		if(maxp0 > 0) {
			assign("seq.ar", expand.grid(1:n, 1:maxp0), frame = 1)
			seq.list <- sapply(1:nrow(seq.ar), function(j)
			list(seq.ar[j,  ]))
			arsapfun <- function(elm)
			round(nullto0(z$detail$BIC[[elm$Var1]]$estimates$ar[elm$Var2], NA), 3)
			ARmat <- matrix(sapply(seq.list, arsapfun), n, maxp0)
			ARname <- paste("ar", 1:maxp0, sep = "")
		}
		else {
			ARmat <- NULL
			ARname <- NULL
		}
		if(maxq0 > 0) {
			assign("seq.ma", expand.grid(1:n, 1:maxq0), frame = 1)
			seq.list <- sapply(1:nrow(seq.ma), function(j)
			list(seq.ma[j,  ]))
			masapfun <- function(elm)
			round(nullto0(z$detail$BIC[[elm$Var1]]$estimates$ma[elm$Var2], NA), 3)
			MAmat <- matrix(sapply(seq.list, masapfun), n, maxq0)
			MAname <- paste("ma", 1:maxq0, sep = "")
		}
		else {
			MAmat <- NULL
			MAname <- NULL
		}
		if(maxps > 0) {
			assign("seq.sar", expand.grid(1:n, 1:maxps), frame = 1)
			seq.list <- sapply(1:nrow(seq.sar), function(j)
			list(seq.sar[j,  ]))
			sarsapfun <- function(elm)
			round(nullto0(z$detail$BIC[[elm$Var1]]$estimates$sar[elm$Var2], NA), 3)
			SARmat <- matrix(sapply(seq.list, sarsapfun), n, maxps)
			SARname <- paste("sar", 1:maxps, sep = "")
		}
		else {
			SARmat <- NULL
			SARname <- NULL
		}
		if(maxqs > 0) {
			assign("seq.sma", expand.grid(1:n, 1:maxqs), frame = 1)
			seq.list <- sapply(1:nrow(seq.sma), function(j)
			list(seq.sma[j,  ]))
			smasapfun <- function(elm)
			round(nullto0(z$detail$BIC[[elm$Var1]]$estimates$sma[elm$Var2], NA), 3)
			SMAmat <- matrix(sapply(seq.list, smasapfun), n, maxqs)
			SMAname <- paste("sma", 1:maxqs, sep = "")
		}
		else {
			SMAmat <- NULL
			SMAname <- NULL
		}
		sigma2vec <- sapply(1:n, function(j)
		signif(nullto0(z$detail$BIC[[j]]$estimates$sigma2, NA), 5))
		estmat <- matrix(c(d0vec, dsvec, as.vector(ARmat), as.vector(MAmat), as.vector(SARmat), as.vector(SMAmat), sigma2vec), nrow = n)
		dimnames(estmat) <- list(NULL, c(d0name, dsname, ARname, MAname, SARname, SMAname, "sigma2"))
		resultBIC <- data.frame(data.frame(zlist), estmat)
		print(resultBIC)
		if(!is.null(z$result$Error)) {
			cat("\n","# Error models #", "\n")
			print(z$result$Error)
			print(z$result$message)
		}
		cat("\n", "# number of models #", "    total:", z$n.model, "   error:", z$n.Error, "\n", "\n")
		invisible(z$n.model)
	}
}


### inner function

sarfima.css.ms.in <- function(model, x, n, xreg = 0, maxlag = 5, control = NULL)
{
	### assignement ###
	arima.order <- model[[1]]$order
	seasonal.order <- model[[2]]$order
	assign("p0", arima.order[1], frame = 1)
	assign("q0", arima.order[3], frame = 1)
	assign("ps", seasonal.order[1], frame = 1)
	assign("qs", seasonal.order[3], frame = 1)
	assign("ss", model[[2]]$period, frame = 1)
	d0.est <- model[[3]]$d0.est
	ds.est <- model[[3]]$ds.est
	np <- model[[3]]$np
	m <- max(min(maxlag * ss, floor(n/3)), ss + 1, np + 1)
	### no estimator case ###
	if(is.one(np)) {
		sigma2 <- sum(x^2)/n
		sigma <- sqrt(sigma2)
		loglik <-  - n * 0.5 * ( log(2 * pi * sigma2) + 1)
		aic <- -2 * loglik + 2 * np
		bic <- -2 * loglik + log(n) * np
		residual <- x
		racf <- meanzero.acf(residual, m)
		mod.Portmanteau <- n * (n + 2) * cumsum(racf[2:(m + 1)]^2/(n - c(1:m)))
		pvalue.chi <- 1 - pchisq(mod.Portmanteau[np:m], c(1:(m - np + 1)))
		chisq.rej.rate <- round((100 * length(pvalue.chi[pvalue.chi < 0.05]))/length(pvalue.chi), digits = 1)
		std.err.racf <- sqrt( (n-c(1:m))/(n*(n+2)) )
		pvalue.racf <- 2*(1 - pnorm( abs(racf[2:(m+1)])/std.err.racf ))
		racf.rej.rate <- round((100 * length(pvalue.racf[pvalue.racf < 0.05]))/length(pvalue.racf), digits = 1)
		list(model = model, xreg = xreg, estimates = list(d0 = NULL, ds = NULL, ar = NULL, ma = NULL, sar = NULL, sma = NULL, sigma2 = sigma2), converged = T, conv.type = NULL, n.used = n, loglik = list(loglik = loglik,
			AIC = aic, BIC = bic), stinv = list(SARMA = T, d0ds = T), chisq = chisq.rej.rate, racf = racf.rej.rate)
	}
	else {
		### calculate initial value ###	
		init.d0ds <- initial.d0ds.nlminb(0, 0, d0.est, ds.est, ss, n, x, control = control)
		init.sarma <- initial.sarma(init.d0ds[1], init.d0ds[2], p0, q0, ps, qs, ss, n, x)
		### sarfima.nlse ###
		fit <- sarfima.nlse.nlminb(x, init.d0ds, init.sarma, ss, d0.est = d0.est, ds.est = ds.est, control = control)
		resd0 <- fit$estimates$d0
		resds <- fit$estimates$ds
		resar <- fit$estimates$ar
		resma <- fit$estimates$ma
		ressar <- fit$estimates$sar
		ressma <- fit$estimates$sma
		n.d0 <- nullto0(resd0)
		n.ds <- nullto0(resds)
		n.ar <- nullto0(resar)
		n.ma <- nullto0(resma)
		n.sar <- nullto0(ressar)
		n.sma <- nullto0(ressma)
		sigma2 <- obj.fun.SARFIMA(n.d0, n.ds, n.ar, n.ma, n.sar, n.sma, ss, n, x)/n
		sigma <- sqrt(sigma2)
		### variane matrix and check the stationarity/invertibily of the fitted coeffs ###
		testSARMA <- SARMA.test(n.ar, n.ma, n.sar, n.sma, p0, q0, ps, qs, ss)
		testd0ds <- if((abs(n.d0 + n.ds) > 0.5) | (abs(n.ds) > 0.5)) F else T
		### Statistics ###
		loglik <-  - n * 0.5 * log(2 * pi * sigma2 + 1)
		aic <- -2 * loglik + 2 * np
		bic <- -2 * loglik + log(n) * np
		### racf and Chi-square ###
		residual <- resid.SARFIMA.gen(n.d0, n.ds, n.ar, n.ma, n.sar, n.sma, ss, n, x)
		racf <- meanzero.acf(residual, m)
		mod.Portmanteau <- n*(n+2)*cumsum(racf[2:(m+1)]^2/(n-c(1:m)))
		pvalue.chi <-  1-pchisq(mod.Portmanteau[(np+1):m], c(1:(m-np)) )
		chisq.rej.rate <- round((100 * length(pvalue.chi[pvalue.chi < 0.05]))/length(pvalue.chi), digits = 1)
		if(testSARMA) {
			var.coef <- sarfima.css.var(resar, resma, ressar, ressma, c(p0, q0, ps, qs), ss, n, d0.est = d0.est, ds.est = ds.est)
			Jm <- matrix.Jm(resar, resma, ressar, ressma, c(p0,q0,ps,qs), maxlag, ss, d0.est = d0.est, ds.est = ds.est)[1:m,]
			std.err.racf <- sqrt( diag( diag(m) - n*(Jm %*% var.coef %*% t(Jm)) ) * (n-c(1:m))/(n*(n+2))  ) 
			pvalue.racf <- 2*(1 - pnorm( abs(racf[2:(m+1)])/std.err.racf ))
			racf.rej.rate <- round((100 * length(pvalue.racf[pvalue.racf < 0.05]))/length(pvalue.racf), digits = 1)
		}
		else racf.rej.rate <- NA
		list(model = model, xreg = xreg, estimates = list(d0 = resd0, ds = resds, ar = resar, ma = resma, sar = ressar, sma = ressma, sigma2 = sigma2), converged = fit$converged, conv.type = fit$conv.type, n.used = n,
			loglik = list(loglik = loglik, AIC = aic, BIC = bic), stinv = list(SARMA = testSARMA, d0ds = testd0ds), chisq = chisq.rej.rate, racf = racf.rej.rate)
	}
}
 


##### MAIN FUNCTION: sarfima.css.ms ############################################################################
# USAGE: 
# sarfima.css.ms(x, period, xreg=0, diff.order=c(0,0), n.best=10, max.sarma=c(1,1,1,1), d0.est=T, ds.est=T, 
#	 max.np=NULL, maxlag = 5, summary.output=T, detail=F, control=NULL)
#
# REQUIRED ARGUMENTS: 
# x 		a univariate time series or a vector. 
# period	an even integer (e.g., ss=4,12).
# OPTIONAL ARGUMENTS: 
# xreg		an estimates of mean of x. default is 0 (e.g., xreg=mean(x)).
# diff.order	the number of differencing, m0 and the number of seasonal differencing, ms. 
# n.best 	upper bound of the number of the selected best models.
# max.sarma	maximum order of SARMA parameters.
# d0.est	if TRUE, d0 is estimeted and if FALSE, d0 is not estimeted.  
# ds.est 	if TRUE, ds is estimeted and if FALSE, ds is not estimeted.  
# max.np	maximum number of parameters.
# maxlag 	maxlag is used for computing maximum lag ACFs of residuals, PACFs of residuals, 
#		and Portmanteau goodness of fit statistics, where maximum lag, m, is given by 
#		m = max(min(maxlag * ss, floor(n/3)), ss + 1, np+1)
# summary.output
#		logical flag: if TRUE, it produces summary of the best models. 
# detail	logical flag: if TRUE, it produces detailed summary of the best models. 
# pass.only	logocal flag: if TRUE, the following models are selected; 
#		1) rate RACF test and Chisquare test rejected lags is less than 10 %, 
#		2) d0, ds, and SARMA parameters satisfy stationary and invertible conditions.
#		3) The estimated SARFIMA parameters all converge.
# cotrol	optional arguments of nlminb
#
# VALUE: 
# a list containing the following elements:
# result$AIC 	AIC by CSS estimates. AIC is given by -2 * loglik + 2 * np. 
# result$BIC	BIC by CSS estimates. BIC is given by -2 * loglik + log(n.used) * np. 
# n.model	the number of models of candidates. 
# n.Error	the number of error models. 
# detail	estimated results with model ID like results of sarfima.css. 
#		If optional arguements detail = FALSE, detail is NULL. 
#		If optional arguements detail = TRUE, detail follows the order of the ranking of AIC and BIC. 
##################################################################################################################

sarfima.css.ms <- function(x, period, xreg = 0, diff.order = c(0, 0), n.best = 10, max.sarma = c(1, 1, 1, 1), 
d0.est = T, ds.est = T, max.np = NULL, maxlag = 5, summary.output = T, detail = F, pass.only = F, control = NULL)
{
	x <- as.vector(x)
	m0 <- diff.order[1]
	ms <- diff.order[2]
	x <- if(m0 > 0) diff(x, lag = 1, differences = m0) else x
	x <- if(ms > 0) diff(x, lag = period, differences = ms) else x
	x <- x - xreg
	# mean zero process
	assign("x", x, frame = 1)
	assign("n", length(x), frame = 1)
	assign("xreg", xreg, frame = 1)
	assign("maxlag", maxlag, frame = 1)
	assign("control", control, frame = 1)
	max.np.sarfima <- if(!is.null(max.np)) max.np - 1
	assign("model.lists", sarfima.allpatterns(period, diff.order = diff.order, max.sarma = max.sarma, max.np.sarfima = max.np.sarfima, d0.est = d0.est, ds.est = 
		ds.est), frame = 1)
	modelnm <- length(model.lists)
	est.result.full <- lapply(1:modelnm, function(a)
	try(sarfima.css.ms.in(model.lists[[a]], x = x, n = n, xreg = xreg, maxlag = maxlag, control = control)))
	assign("est.result.F", model.lists[sapply(est.result.full, is, "Error")], frame = 1)
	assign("est.result.F.em", est.result.full[sapply(est.result.full, is, "Error")], frame = 1)
	### table of est.result.F ### 
	modelnm.F <- length(est.result.F)
	if(modelnm.F > 0) {
	result.df.F <- data.frame(list(ID = sapply(1:modelnm.F, 
	function(j)est.result.F[[j]][[3]]$ID), 
	p0 = sapply(1:modelnm.F, function(j)est.result.F[[j]][[1]]$order[1]), 
	m0 = sapply(1:modelnm.F, function(j)est.result.F[[j]][[1]]$order[2]), 
	q0 = sapply(1:modelnm.F, function(j)est.result.F[[j]][[1]]$order[3]), 
	ps = sapply(1:modelnm.F, function(j)est.result.F[[j]][[2]]$order[1]), 
	ms = sapply(1:modelnm.F, function(j)est.result.F[[j]][[2]]$order[2]), 
	qs = sapply(1:modelnm.F, function(j)est.result.F[[j]][[2]]$order[3]), 
	d0.est = sapply(1:modelnm.F, function(j)est.result.F[[j]][[3]]$d0.est), 
	ds.est = sapply(1:modelnm.F, function(j)est.result.F[[j]][[3]]$ds.est), 
	np = sapply(1:modelnm.F, function(j)est.result.F[[j]][[3]]$np)), 
	row.names = paste(" ", 1:modelnm.F))
		result.df.F.em <- data.frame(ErrorMessage = sapply(1:modelnm.F, 
		function(j)est.result.F.em[[j]][[1]]), row.names = paste(" ", 1:modelnm.F))
	}
	else {
		result.df.F <- NULL
		result.df.F.em <- NULL
	}
	### pass.only ###
	est.result.ne <- est.result.full[!sapply(est.result.full, is, "Error")]
	if(pass.only) {
		pass.only.ID <- lapply(1:modelnm, function(j, est.result)
		{
			est <- est.result[[j]]
			chisq <- est$chisq <= 10
			racf <- est$racf <= 10
			SARMA <- est$stinv$SARMA
			d0ds <- est$stinv$d0ds
			convd <- est$converged
			if(all(chisq, racf, SARMA, d0ds, convd))
				est$model[[3]]$ID
		}
		, est.result = est.result.ne)
		pass.only <- lapply(unlist(pass.only.ID), function(j, est.result)
		est.result[[j]], est.result = est.result.ne)
		if(is.zero(length(pass.only)))
			warning("Warning: No results passed on the test given by pass.only assignment")
		else est.result.ne <- pass.only
	}
	assign("est.result", est.result.ne, frame = 1)
	###  table of est.result  ### 
	modelnm.T <- length(est.result)
	n.best <- min(modelnm.T, n.best)
	ms.mat <- sapply(1:modelnm.T, function(j)
	c(est.result[[j]]$loglik$AIC, est.result[[j]]$loglik$BIC))
	### sort by AIC ###
	assign("est.result.sort.A", lapply(order(ms.mat[1,  ])[1:n.best], function(j)
	est.result[[j]]), frame = 1)
	result.df.A <- data.frame(list(ID = sapply(1:n.best, function(j)
	est.result.sort.A[[j]]$model[[3]]$ID), AIC = sapply(1:n.best, function(j)
	signif(est.result.sort.A[[j]]$loglik$AIC, 5)), BIC = sapply(1:n.best, function(j)
	signif(est.result.sort.A[[j]]$loglik$BIC, 5)), ChiSq = sapply(1:n.best, function(j)
	est.result.sort.A[[j]]$chisq), RACF = sapply(1:n.best, function(j)
	est.result.sort.A[[j]]$racf), SARMA = sapply(1:n.best, function(j)
	est.result.sort.A[[j]]$stinv$SARMA), d0ds = sapply(1:n.best, function(j)
	est.result.sort.A[[j]]$stinv$d0ds), convd = sapply(1:n.best, function(j)
	est.result.sort.A[[j]]$converged), p0 = sapply(1:n.best, function(j)
	est.result.sort.A[[j]]$model[[1]]$order[1]), m0 = sapply(1:n.best, function(j)
	est.result.sort.A[[j]]$model[[1]]$order[2]), q0 = sapply(1:n.best, function(j)
	est.result.sort.A[[j]]$model[[1]]$order[3]), ps = sapply(1:n.best, function(j)
	est.result.sort.A[[j]]$model[[2]]$order[1]), ms = sapply(1:n.best, function(j)
	est.result.sort.A[[j]]$model[[2]]$order[2]), qs = sapply(1:n.best, function(j)
	est.result.sort.A[[j]]$model[[2]]$order[3]), d0.est = sapply(1:n.best, function(j)
	est.result.sort.A[[j]]$model[[3]]$d0.est), ds.est = sapply(1:n.best, function(j)
	est.result.sort.A[[j]]$model[[3]]$ds.est), np = sapply(1:n.best, function(j)
	est.result.sort.A[[j]]$model[[3]]$np)), row.names = paste(" ", 1:n.best))
	### sort by BIC ###
	assign("est.result.sort.B", lapply(order(ms.mat[2,  ])[1:n.best], function(j)
	est.result[[j]]), frame = 1)
	result.df.B <- data.frame(list(ID = sapply(1:n.best, function(j)
	est.result.sort.B[[j]]$model[[3]]$ID), BIC = sapply(1:n.best, function(j)
	signif(est.result.sort.B[[j]]$loglik$BIC, 5)), AIC = sapply(1:n.best, function(j)
	signif(est.result.sort.B[[j]]$loglik$AIC, 5)), ChiSq = sapply(1:n.best, function(j)
	est.result.sort.B[[j]]$chisq), RACF = sapply(1:n.best, function(j)
	est.result.sort.B[[j]]$racf), SARMA = sapply(1:n.best, function(j)
	est.result.sort.B[[j]]$stinv$SARMA), d0ds = sapply(1:n.best, function(j)
	est.result.sort.B[[j]]$stinv$d0ds), convd = sapply(1:n.best, function(j)
	est.result.sort.B[[j]]$converged), p0 = sapply(1:n.best, function(j)
	est.result.sort.B[[j]]$model[[1]]$order[1]), m0 = sapply(1:n.best, function(j)
	est.result.sort.B[[j]]$model[[1]]$order[2]), q0 = sapply(1:n.best, function(j)
	est.result.sort.B[[j]]$model[[1]]$order[3]), ps = sapply(1:n.best, function(j)
	est.result.sort.B[[j]]$model[[2]]$order[1]), ms = sapply(1:n.best, function(j)
	est.result.sort.B[[j]]$model[[2]]$order[2]), qs = sapply(1:n.best, function(j)
	est.result.sort.B[[j]]$model[[2]]$order[3]), d0.est = sapply(1:n.best, function(j)
	est.result.sort.B[[j]]$model[[3]]$d0.est), ds.est = sapply(1:n.best, function(j)
	est.result.sort.B[[j]]$model[[3]]$ds.est), np = sapply(1:n.best, function(j)
	est.result.sort.B[[j]]$model[[3]]$np)), row.names = paste(" ", 1:n.best))
	detail.list <- if(detail) list(AIC = est.result.sort.A, BIC = est.result.sort.B) else NULL
	resultlist <- list(result = list(AIC = result.df.A, BIC = result.df.B, Error = result.df.F, message = result.df.F.em), n.model = modelnm, n.Error = modelnm.F,
		detail = detail.list)
	if(summary.output) {
		sarfima.ms.summary(resultlist)
		invisible(resultlist)
	}
	else resultlist
}
 







### example 1 ################################################################################
data <- rts(arima.fracdiff.sim(model = list(d = 0.3, ar = -0.4, ma = 0.4), n = 200), start = c(1980,1), frequency = 12)
result11 <- sarfima.css.ms(data, period=12, max.sarma=c(1,1,1,1), n.best=3, max.np=4, detail=F)
result12 <- sarfima.css.ms(data, period=12, max.sarma=c(1,1,1,1), n.best=3, max.np=4, detail=T)
##############################################################################################

### example 2 ################################################################################
data <- arima.fracdiff.sim(model = list(d = 0.3, ar = -0.2, ma = 0.4), n = 200)
data <- cumsum(data+1)[101:200]
result21 <- sarfima.css.ms(data, period=12, xreg=mean(diff(data)), n.best=3, diff.order=c(1,0), max.sarma=c(2,2,0,0), detail=T)
result22 <- sarfima.css.ms(data, period=12, xreg=mean(diff(data)), n.best=3, diff.order=c(1,0), max.sarma=c(2,2,0,0), ds.est=F, detail=T)
result23 <- sarfima.css.ms(data, period=12, xreg=mean(diff(data)), n.best=3, diff.order=c(1,0), max.sarma=c(2,2,0,0), ds.est=F, max.np=4, detail=T, pass.only = T)
##############################################################################################

### summary of tables ##########################################################################
#  In short summary (detail=F), summary of AIC modelselection and  summary of BIC modelselection (detail=T),
# the numbers in column of AIC (BIC) denote the ranking of models in terms of AIC (BIC).
#  Note that if SARMA is FALSE, RACF is not calculated and set to be NA because it is impossible 
# to calculate a variance-covariance matrix of estimated parameters of the SARFIMA model.
# ID		the model identification within total models.
# AIC		AIC by CSS estimates. AIC is given by -2 * loglik + 2 * np. 
# BIC		BIC by CSS estimates. BIC is given by -2 * loglik + log(n.used) * np.
# ChiSq		100 * (rejected number of Modified portmateau tests) / (total number of Modified portmateau tests).
# RACF		100 * (rejected number of RACF tests) / (total number of RACF tests).
# SARMA		logical flag: if TRUE, SARMA parameters satisfy stationary and invertible conditions. 
# d0ds		logical flag: if TRUE, d0 and ds satisfy stationary and invertible conditions. 
# convd		logical flag: if TRUE, the estimated SARFIMA parameters all converge. 
# p0, m0, q0, ps, ms, qs, 
#		the order of the SARIMA(p0,m0,q0)(ps,ms,qs)s model. 
# d0.est	if TRUE, d0 is estimeted and if FALSE, d0 is not estimeted.  
# ds.est 	if TRUE, ds is estimeted and if FALSE, ds is not estimeted.  
# np		number of parameters defined by p0 + q0 + ps + qs + d0.est + ds.est + 1.
# d0, ds, ar1,..., sigma2, 
#		estimates. NA indicates the corresponding parameter is not estimated and set to be 0.   
################################################################################################


> data <- rts(arima.fracdiff.sim(model = list(d = 0.3, ar = -0.4, ma = 0.4), n = 200), start = c(1980,1), frequency = 12)
> result11 <- sarfima.css.ms(data, period=12, max.sarma=c(1,1,1,1), n.best=3, max.np=4, detail=F)
# short summary # 
$AIC:
    ID    AIC    BIC ChiSq RACF SARMA d0ds convd p0 m0 q0 ps ms qs d0.est ds.est np 
  1 19 371.29 384.48   0.0  3.3     T    T     T  1  0  1  0  0  0      T      F  4
  2 17 372.68 382.57   3.4  6.7     T    T     T  1  0  0  0  0  0      T      F  3
  3 39 374.63 387.82  12.3  6.7     T    T     T  1  0  0  0  0  0      T      T  4

$BIC:
    ID    BIC    AIC ChiSq RACF SARMA d0ds convd p0 m0 q0 ps ms qs d0.est ds.est np 
  1 17 382.57 372.68   3.4  6.7     T    T     T  1  0  0  0  0  0      T      F  3
  2  2 383.26 376.66  61.0  8.3     T    T     T  1  0  0  0  0  0      F      F  2
  3 19 384.48 371.29   0.0  3.3     T    T     T  1  0  1  0  0  0      T      F  4

$Error:
NULL

$message:
NULL


 # number of models # 42 
 
> result12 <- sarfima.css.ms(data, period=12, max.sarma=c(1,1,1,1), n.best=3, max.np=4, detail=T)
# summary of AIC model selection  # 
  ID    AIC    BIC ChiSq RACF SARMA d0ds convd p0 m0 q0 ps ms qs d0.est ds.est np    d0     ds    ar1   ma1  sigma2 
1 19 371.29 384.48   0.0  3.3     T    T     T  1  0  1  0  0  0      T      F  4 0.299     NA -0.416 0.324 0.81965
2 17 372.68 382.57   3.4  6.7     T    T     T  1  0  0  0  0  0      T      F  3 0.139     NA -0.534    NA 0.83636
3 39 374.63 387.82  12.3  6.7     T    T     T  1  0  0  0  0  0      T      T  4 0.141 -0.015 -0.537    NA 0.83612

 # summary of BIC model selection  # 
  ID    BIC    AIC ChiSq RACF SARMA d0ds convd p0 m0 q0 ps ms qs d0.est ds.est np    d0    ar1   ma1  sigma2 
1 17 382.57 372.68   3.4  6.7     T    T     T  1  0  0  0  0  0      T      F  3 0.139 -0.534    NA 0.83636
2  2 383.26 376.66  61.0  8.3     T    T     T  1  0  0  0  0  0      F      F  2    NA -0.429    NA 0.86660
3 19 384.48 371.29   0.0  3.3     T    T     T  1  0  1  0  0  0      T      F  4 0.299 -0.416 0.324 0.81965

 # number of models #     total: 42    error: 0 
 
> ##############################################################################################
> 
> 
> 
> ### example 2 ################################################################################
> data <- arima.fracdiff.sim(model = list(d = 0.3, ar = -0.2, ma = 0.4), n = 200)
> data <- cumsum(data+1)[101:200]
> result21 <- sarfima.css.ms(data, period=12, xreg=mean(diff(data)), n.best=3, diff.order=c(1,0), max.sarma=c(2,2,0,0), detail=T)
# summary of AIC model selection  # 
  ID    AIC    BIC ChiSq RACF SARMA d0ds convd p0 m0 q0 ps ms qs d0.est ds.est np     ds    ar1    ar2   ma1    ma2  sigma2 
1 27 208.32 223.89     0  3.0     T    T     T  2  1  2  0  0  0      F      T  6 -0.149  0.065 -0.685 0.340 -0.848 0.99704
2  9 208.37 221.35     0  0.0     T    T     T  2  1  2  0  0  0      F      F  5     NA  0.066 -0.571 0.355 -0.761 1.02130
3  2 209.19 214.38     0  9.1     T    T     T  1  1  0  0  0  0      F      F  2     NA -0.275     NA    NA     NA 1.10540

 # summary of BIC model selection  # 
  ID    BIC    AIC ChiSq RACF SARMA d0ds convd p0 m0 q0 ps ms qs d0.est ds.est np    ar1   ma1 sigma2 
1  2 214.38 209.19   0.0  9.1     T    T     T  1  1  0  0  0  0      F      F  2 -0.275    NA 1.1054
2  4 215.17 209.97   0.0  9.1     T    T     T  0  1  1  0  0  0      F      F  2     NA 0.238 1.1155
3  1 216.60 214.01  78.8  9.1     T    T     T  0  1  0  0  0  0      F      F  1     NA    NA 1.1956

 # number of models #     total: 36    error: 0 
 
> result22 <- sarfima.css.ms(data, period=12, xreg=mean(diff(data)), n.best=3, diff.order=c(1,0), max.sarma=c(2,2,0,0), ds.est=F, detail=T)
# summary of AIC model selection  # 
  ID    AIC    BIC ChiSq RACF SARMA d0ds convd p0 m0 q0 ps ms qs d0.est ds.est np    ar1    ar2   ma1    ma2 sigma2 
1  9 208.37 221.35     0  0.0     T    T     T  2  1  2  0  0  0      F      F  5  0.066 -0.571 0.355 -0.761 1.0213
2  2 209.19 214.38     0  9.1     T    T     T  1  1  0  0  0  0      F      F  2 -0.275     NA    NA     NA 1.1054
3  7 209.25 217.04     0  6.1     T    T     T  0  1  2  0  0  0      F      F  3     NA     NA 0.332 -0.230 1.0809

 # summary of BIC model selection  # 
  ID    BIC    AIC ChiSq RACF SARMA d0ds convd p0 m0 q0 ps ms qs d0.est ds.est np    ar1   ma1 sigma2 
1  2 214.38 209.19   0.0  9.1     T    T     T  1  1  0  0  0  0      F      F  2 -0.275    NA 1.1054
2  4 215.17 209.97   0.0  9.1     T    T     T  0  1  1  0  0  0      F      F  2     NA 0.238 1.1155
3  1 216.60 214.01  78.8  9.1     T    T     T  0  1  0  0  0  0      F      F  1     NA    NA 1.1956

 # number of models #     total: 18    error: 0 
 
> result23 <- sarfima.css.ms(data, period=12, xreg=mean(diff(data)), n.best=3, diff.order=c(1,0), max.sarma=c(2,2,0,0), ds.est=F, max.np=4, detail=T, pass.only = T)
# summary of AIC model selection  # 
  ID    AIC    BIC ChiSq RACF SARMA d0ds convd p0 m0 q0 ps ms qs d0.est ds.est np    ar1   ma1   ma2 sigma2 
1  2 209.19 214.38     0  9.1     T    T     T  1  1  0  0  0  0      F      F  2 -0.275    NA    NA 1.1054
2  7 209.25 217.04     0  6.1     T    T     T  0  1  2  0  0  0      F      F  3     NA 0.332 -0.23 1.0809
3  4 209.97 215.17     0  9.1     T    T     T  0  1  1  0  0  0      F      F  2     NA 0.238    NA 1.1155

 # summary of BIC model selection  # 
  ID    BIC    AIC ChiSq RACF SARMA d0ds convd p0 m0 q0 ps ms qs d0.est ds.est np    ar1   ma1   ma2 sigma2 
1  2 214.38 209.19     0  9.1     T    T     T  1  1  0  0  0  0      F      F  2 -0.275    NA    NA 1.1054
2  4 215.17 209.97     0  9.1     T    T     T  0  1  1  0  0  0      F      F  2     NA 0.238    NA 1.1155
3  7 217.04 209.25     0  6.1     T    T     T  0  1  2  0  0  0      F      F  3     NA 0.332 -0.23 1.0809

 # number of models #     total: 14    error: 0 
 
> ##############################################################################################
> 
>