####################################################################################
### sarfima.css: Conduct CSS estimation of the SARFIMA(p,d0,q)(ps,ds,qs)s model, ###
###  x={x_t}. This model is given by 			       			 ###
### (1-L)^m0 (1-L^s)^ms x_t - xreg = y_t,					 ###
### (1-L)^d0 (1-L^s)^ds y_t = u_t,						 ###
### u_t: the SARMA(p0,q0)(ps,qs) model.	                                         ###
### In this model, s is an even integer, m0 (ms) is a positive integer, and s,   ###
### m0, ms are known. xreg is an mean of the (seasonally) differenced process.   ###
### u_t is defined by the polynominals of backshift operator of the 	 	 ###
### SARIMA model in S-PLUS function when differencing paramters set to be zero.  ###
### See Section Univariate ARIMA modeling in Chapter 27, Guide to Statistics.    ###
### d0 and ds are real number and |d0+ds|,|ds| < 0.5 for stationary and          ###
### invertibility.                                                               ###
####################################################################################

### conv #########################################################################
# convolution function by M.Shibuya and R.Shibata (1992, p.78)
##################################################################################
 
conv <- function(a, b)
{
	m <- length(a)
	n <- length(b)
	aa <- c(rep(c(a, rep(0, n)), n - 1), a)
	z <- numeric(m + n - 1)
	.Fortran("dmatp",
		as.double(aa),
		as.integer(c(m + n - 1, n)),
		as.double(b),
		as.integer(c(n, 1)),
		z = z)$z
}



### pij.FId.gen ##################################################################
# returns n AR(infty) coefficients of an ARFIMA(0,d0,0) process 
# with a fractional filter (1-L)^d0. 
##################################################################################

pij.FId.gen <- function(d, n)
{
	result <- vector(length = n)
	result1 <- 1
	j <- 2:n
	result2 <- cumprod((j - 2 - d)/(j - 1))
	c(result1, result2)
}


### pij.SFId.gen ###################################################################
# returns n AR(infty) coefficients of an SARFIMA(0,0,0)(0,ds,0)ss process 
# with a fractional filter (1-L^ss)^ds. 
####################################################################################

pij.SFId.gen <- function(ds, ss, n){
ssn <- ceiling(n/ss)
nn <- ssn * ss 
vec0 <- rep(0, nn)
vec0[c(0:(ssn-1))*ss + 1] <- pij.FId.gen(ds, ssn)
vec0[1:n]
}

# or 
# pij.SFId.gen <- function(ds, ss, n)as.vector(matrix(c(pij.FId.gen(ds,n),rep(0,(ss-1)*n)), ss, n, byrow=T) )[1:n]

### pij.S.gen ######################################################################
# returns n AR(infty) coefficients of the SARFIMA(0,d0,0)(0,ds,0)ss process 
# with a fractional filter (1-L)^d0 (1-L^ss)^ds. 
####################################################################################

pij.S.gen <- function(d0,ds,ss,n)
(conv(pij.FId.gen(d0,n), as.vector(matrix(c(pij.FId.gen(ds,n),rep(0,(ss-1)*n)), ss, n, byrow=T) )[1:n]))[1:n]

# another diffinition of pij.S.gen
# pij.S.gen2 <- function(d0,ds,ss,n)conv(pij.FId.gen(d0,n), pij.SFId.gen(ds,ss,n))[1:n]


### frac.diff #########################################################################
# returns a vector of fractional differencing the data x. 
#######################################################################################

frac.diff <- function(x, d, lag = 1)
{
	n <- length(x)
	if(lag == 1)
		(conv(pij.FId.gen(d, n), x))[1:n]
	else (conv(as.vector(matrix(c(pij.FId.gen(d, n), rep(0, (lag - 1) * n)), lag, n, byrow = T))[1:n], x))[1:n]
}



### pij.SARMA.gen ######################################################################
# returns SARMA model's AR(infty) coefficietns.
# ARcoefs:  1 - ar1 * L - ar2 * L^2 - ... - arp * L^p, pvec <- c(ar1,ar2,...,arp)
# MAcoefs:  1 - ma1 * L - ma2 * L^2 - ... - maq * L^q, qvec <- c(ma1,ma2,...,maq)
# SARcoefs: 1 - sar1 * L - sar2 * L^2 - ... - sarP * L^P, Pvec <- c(sar1,sar2,...,sarP)
# SMAcoefs: 1 - sma1 * L - sma2 * L^2 - ... - smaQ * L^Q, Qvec <- c(sma1,sma2,...,smaQ)
#########################################################################################


ARcoefs.gen <- function(pvec, n)c(1, -pvec, rep(0, n))[1:n]


SARcoefs.gen <- function(Pvec, ss, n)
{
	P <- length(Pvec)
	c(as.vector(matrix(c(1,  - Pvec, rep(0, (ss - 1) * (P + 1))), ss, P + 1, byrow = T)), rep(0, n))[1:n]
}


MAcoefs.gen <- function(qvec, n)c(1, -qvec, rep(0, n))[1:n]


SMAcoefs.gen <- function(Qvec, ss, n)
{
	Q <- length(Qvec)
	c(as.vector(matrix(c(1, -Qvec, rep(0, (ss - 1) * (Q + 1))), ss, Q + 1, byrow = T)), rep(0, n))[1:n]
}


pij.SARMA.gen <- function(pvec, qvec, Pvec, Qvec, ss, n)
{
	arseq <- conv(ARcoefs.gen(pvec, n), SARcoefs.gen(Pvec, ss, n))[1:n]
	maseq <- conv(MAcoefs.gen(qvec, n), SMAcoefs.gen(Qvec, ss, n))[1:n]
	pij <- vector(length = n)
	pij[1] <- 1
	for(j in 1:(n - 1)) {
		pij[j + 1] <- arseq[j + 1] - t(maseq[2:(j + 1)]) %*% pij[j:1]
	}
	drop(pij)
}



### pij.SARFIMA.gen ##############################################################
# returns a vector of AR(infty) coefficients of the SARFIMA(p,d0,q)(ps,ds,qs)s
# model.
##################################################################################

pij.SARFIMA.gen <- function(d0, ds, pvec, qvec, Pvec, Qvec, ss, n)
conv( pij.S.gen(d0,ds,ss,n), pij.SARMA.gen(pvec, qvec, Pvec, Qvec, ss, n) )[1:n]



###################################################################################
# resid.SARFIMA.gen : returns a residual sequences of the SARFIMA model
###################################################################################

resid.SARFIMA.gen <- function(d0, ds, pvec, qvec, Pvec, Qvec, ss, n, x)
(conv(pij.SARFIMA.gen(d0, ds, pvec, qvec, Pvec, Qvec, ss, n), x[1:n]))[1:n]



### initial.d0ds.nlminb #################################################################
# calculate initial values of differencing parameters of the SARFIMA models. 
# ARGUMENTS: 
# d0: initial value if d0 is estimeted or preassigned value if d0 is not estimated
# ds: initial value if ds is estimeted or preassigned value if ds is not estimated
# 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.  
# ss: period
# n: length(x)
# x : a mean zero stationary process
# scale, cotrol: optional arguments of nlminb
#########################################################################################

# obj.fun.SARFIMA: objective function of the SSR of the SARFIMA model

obj.fun.SARFIMA <- function(d0, ds, pvec, qvec, Pvec, Qvec, ss, n, x)
sum(resid.SARFIMA.gen(d0, ds, pvec, qvec, Pvec, Qvec, ss, n, x)^2)

# derivertive matrix w.r.t d0 and ds 

matrix.deriv.e.d <- function(ss, e, d0.est = T, ds.est = T)
{
	n <- length(e)
	d0vec <- if(d0.est) (conv(c(0, -1/(1:(n-1))), e))[1:n] else NULL
	if(ds.est) {
		h <- floor(n/ss) + 1
		sk <- rep(0, h*ss)
		sk[ss * (1:h)] <- -1/(1:h)
		dsvec <- (conv(c(0, sk[1:(n-1)]), e))[1:n]
	}
	else dsvec <- NULL
	matrix(c(d0vec, dsvec), nrow=n)
}


obj.fun.d0ds <- function(d0, ds, ss, n, x)sum(resid.d0ds.gen(d0, ds, ss, n, x)^2)
resid.d0ds.gen <- function(d0, ds, ss, n, x)(conv(pij.S.gen(d0, ds, ss, n), x[1:n]))[1:n]

initial.d0ds.nlminb <- function(d0, ds, d0.est, ds.est, ss, n, x, scale = 1, control=NULL)
{
	assign("d0", d0, frame = 1)
	assign("ds", ds, frame = 1)
	assign("d0.est", d0.est, frame = 1)
	assign("ds.est", ds.est, frame = 1)
	assign("ss", ss, frame = 1)
	assign("n", n, frame = 1)
	assign("x", x, frame = 1)
	obj.in <- function(eta)
	{
		d0eta <- if(d0.est) eta[1] else d0
		dseta <- if(ds.est) eta[(d0.est + 1)] else ds
		obj.fun.d0ds(d0eta, dseta, ss, n, x) * 0.5 /n
	}
	gradhess.in <- function(eta)
	{
		d0eta <- if(d0.est) eta[1] else d0
		dseta <- if(ds.est) eta[d0.est + 1] else ds
		e <- resid.d0ds.gen(d0eta, dseta, ss, n, x)
		deriv.e <- matrix.deriv.e.d(ss, e, d0.est = d0.est, ds.est = ds.est)
		hess <- t(deriv.e) %*% deriv.e/n
		list(gradient = t(deriv.e) %*% e/n, 
		hessian = switch(d0.est+ds.est+1, NULL, hess[1,1], c(hess[1,1], hess[2,1], hess[2,2])))
	}
	if(all(!d0.est, !ds.est)) result <-  c(d0, ds)
	else {
	dstart <- switch(d0.est+ds.est+1, NULL, d0.est*d0 + ds.est*ds, c(d0, ds))
	dupper <- switch(d0.est+ds.est+1, NULL, d0.est + ds.est*0.5, c(1.0, 0.5))
	etaest <- nlminb(start = dstart, objective = obj.in, gradient = gradhess.in, hessian = T, scale = scale, control = control, lower=-dupper, upper = dupper)$parameters
	result <- switch(d0.est+2*ds.est+1, NULL, c(etaest, ds), c(d0, etaest), etaest)
	}
	result
}


### example ######################################################################
# data <- rnorm(100)
### estimate d0 & ds.
# initial.d0ds.nlminb(0, 0, d0.est=T, ds.est=T, 12, 100, data) 
### estimate d0 and ds is preassigned value (0).
# initial.d0ds.nlminb(0, 0, d0.est=T, ds.est=F, 12, 100, data)
### estimate ds and d0 is preassigned value (0). 
# initial.d0ds.nlminb(0, 0, d0.est=F, ds.est=T, 12, 100, data)
### d0 and ds are preassigned value, c(0,0). 
# initial.d0ds.nlminb(0, 0, d0.est=F, ds.est=F, 12, 100, data)
###################################################################################


### initial.sarma ######################################################################
# initial.sarma: calculate initial values of SARMA parameters of the SARFIMA models 
# by arima.mle. 
# ARGUMENTS: 
# d0: initial value of d0. d0 and ds are used for fractional differencing. 
# ds: initial value if ds. d0 and ds are used for fractional differencing. 
# p0, q0, ps, qs: the number of SARMA parameters. Namely, the number of AR, MA, SAR, SMA
# 	parameters, respectively. 
# ss: period
# n: length(x)
# x : a mean zero stationary process
# sarma.start: intial values of SARMA parameters such that 
#	sarma.start <- list(ar=c(...), ma=c(...), sar=c(...), sma=c(...)). 
# 	If sarma.start is not NULL, SARMA prameters are not estimated and returns values 
# 	of sarma.start.
#########################################################################################


initial.sarma <- function(d0, ds, p0, q0, ps, qs, ss, n, x, sarma.start = NULL)
{
	if(is.null(sarma.start) == T) {
		# series x is seasonally and fractionally differenced
		diffx <- frac.diff(frac.diff(x, d0), ds, lag = ss)	
		# model of arima.mle
		if((p0 == 0) && (q0 == 0) && (ps == 0) && (qs == 0))
			list(order = c(p0, q0, ps, qs), ar = NULL, ma = NULL, sar = NULL, sma = NULL)
		else if((p0 == 0) && (q0 == 0)) {
			fit <- arima.mle(diffx, list(order = c(ps, 0, qs), period = ss))
			init.ps <- if(ps > 0) fit$model$ar else NULL
			init.qs <- if(qs > 0) fit$model$ma else NULL
			list(order = c(p0, q0, ps, qs), ar = NULL, ma = NULL, sar = c(init.ps), sma = (init.qs))
		}
		else if((ps == 0) && (qs == 0)) {
			fit <- arima.mle(diffx, list(order = c(p0, 0, q0)))
			init.p0 <- if(p0 > 0) fit$model$ar else NULL
			init.q0 <- if(q0 > 0) fit$model$ma else NULL
			list(order = c(p0, q0, ps, qs), ar = init.p0, ma = init.q0, sar = NULL, sma = NULL)
		}
		else {
			fit <- arima.mle(diffx, list(list(order = c(p0, 0, q0)), list(order = c(ps, 0, qs), period = ss)))
			init.p0 <- if(p0 > 0) fit$model[[1]]$ar else NULL
			init.q0 <- if(q0 > 0) fit$model[[1]]$ma else NULL
			init.ps <- if(ps > 0) fit$model[[2]]$ar else NULL
			init.qs <- if(qs > 0) fit$model[[2]]$ma else NULL
			list(order = c(p0, q0, ps, qs), ar = c(init.p0), ma = c(init.q0), sar = c(init.ps), sma = c(init.qs))
		}
	}
	else list(order = c(p0, q0, ps, qs), ar = c(sarma.start$ar), ma = c(sarma.start$ma), sar = sarma.start$sar, sma = sarma.start$sma)
}
 

### example ######################################################################
# data <-rnorm(100)
# sarma.start <- list(ar=0,ma=0)
# initial.sarma(0,0,1,1,0,0,12,100,data)
# arima.mle(data,list(order=c(1,0,1)))$model
# initial.sarma(0,0,1,1,0,0,12,100,data,sarma.start)
##################################################################################


### matrix.deriv.e ###############################################################
# derivertive matrix w.r.t d0, ds, and SARMA parameters 
# ARGUMENTS: 
# p0vec, q0vec, psvec, qsvec: vectors of AR, MA, SAR, SMA parameters, respectively. 
# order: a list of the order of SARMA parameters 
# ss: period
# e: residual vector, n==length(e)
# 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.  
##################################################################################


matrix.deriv.e <- function(p0vec, q0vec, psvec, qsvec, order, ss, e, d0.est = T, ds.est = T)
{
	n <- length(e)
	p0 <- order[1]
	q0 <- order[2]
	ps <- order[3]
	qs <- order[4]
	col.mat <- d0.est + ds.est + p0 + q0 + ps + qs
	d0vec <- if(d0.est) (conv(c(0, -1/(1:(n-1))), e))[1:n] else NULL
	if(ds.est) {
		h <- floor(n/ss) + 1
		sk <- rep(0, h*ss)
		sk[ss * (1:h)] <- -1/(1:h)
		dsvec <- (conv(c(0, sk[1:(n-1)]), e))[1:n]
	}
	else dsvec <- NULL
	if(p0 > 0) {
		ARmat <- matrix(rep(c((conv(c(0, -inv.coef.gen(ARcoefs.gen(p0vec, n-1))), e))[1:n], 0), p0)[1:(n * p0)], n, p0)
		ARmat[col(ARmat) - row(ARmat) > 0] <- 0
	}
	else ARmat <- NULL
	if(q0 > 0) {
		MAmat <- matrix(rep(c((conv(c(0, inv.coef.gen(MAcoefs.gen(q0vec, n-1))), e))[1:n], 0), q0)[1:(n * q0)], n, q0)
		MAmat[col(MAmat) - row(MAmat) > 0] <- 0
	}
	else MAmat <- NULL
	if(ps > 0) {
		SARmat <- matrix(rep(c((conv(c(rep(0, ss), -inv.coef.gen(SARcoefs.gen(psvec, ss, n-ss))), e))[1:n], rep(0, ss)), ps)[1:(n * ps)], n, ps)
		SARmat[ss*col(SARmat) - row(SARmat) > 0] <- 0
	}
	else SARmat <- NULL
	if(qs > 0) {
		SMAmat <- matrix(rep(c((conv(c(rep(0, ss), inv.coef.gen(SMAcoefs.gen(qsvec, ss, n-ss))), e))[1:n], rep(0, ss)), qs)[1:(n * qs)], n, qs)
		SMAmat[ss*col(SMAmat) - row(SMAmat) > 0] <- 0
	}
	else SMAmat <- NULL
	matrix(c(d0vec, dsvec, as.vector(ARmat), as.vector(MAmat), as.vector(SARmat), as.vector(SMAmat)), n, col.mat)
}



### sarfima.nlse ###############################################################################
# calculate NLSEs of eta, where eta is a vector of d0, ds, SARMA parameters
# ARGUMENTS: 
# init.d0ds: initial values of d0 and ds
# init.SARMA: initial values of SARMA parameters
# ss: period
# 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.  
# scale, cotrol, lower, upper: optional arguments of nlminb
#########################################################################################

### convert a Hessian matrix to a nlminb's Hessian sequence 

nlminb.hessian.seq <- function(matx)
{
n <- nrow(matx)
matx[row(matx)- col(matx) < 0 ]  <- NA
as.vector(na.omit(as.vector(t(matx))))
}

### calculate NLSEs by nlminb 

sarfima.nlse.nlminb <- function(x, init.d0ds, init.SARMA, ss, d0.est = T, ds.est = T, scale=1, control=NULL, lower = -Inf, upper = Inf)
{
	assign("x", x, frame = 1)
	assign("n", length(x), frame = 1)
	assign("init.d0", init.d0ds[1], frame = 1)
	assign("init.ds", init.d0ds[2], frame = 1)
	init.d0.n <- if(d0.est) init.d0 else NULL
	init.ds.n <- if(ds.est) init.ds else NULL
	sarma.order <- c(init.SARMA$order)
	init.sarma <- c(c(init.SARMA$ar), c(init.SARMA$ma), c(init.SARMA$sar), c(init.SARMA$sma))
	assign("p0", sarma.order[1], frame = 1)
	assign("q0", sarma.order[2], frame = 1)
	assign("ps", sarma.order[3], frame = 1)
	assign("qs", sarma.order[4], frame = 1)
	assign("ss", ss, frame = 1)
	assign("d0.est", d0.est, frame = 1)
	assign("ds.est", ds.est, frame = 1)
	conv.test <- function(etafit)
	{
		conv.type <- etafit$message
		if(any(conv.type == "X CONVERGENCE", 
		 conv.type == "RELATIVE FUNCTION CONVERGENCE", 
		 conv.type == "BOTH X AND RELATIVE FUNCTION CONVERGENCE", 
		 conv.type == "ABSOLUTE FUNCTION CONVERGENCE")) TRUE
		else FALSE
	}
	obj.in <- function(eta)
	{
		d0eta <- if(d0.est) eta[1] else init.d0
		dseta <- if(ds.est) eta[(d0.est + 1)] else init.ds
		p0vec <- if(p0 > 0) eta[(d0.est + ds.est + 1):(d0.est + ds.est + p0)] else 0
		q0vec <- if(q0 > 0) eta[(d0.est + ds.est + p0 + 1):(d0.est + ds.est + p0 + q0)] else 0
		psvec <- if(ps > 0) eta[(d0.est + ds.est + p0 + q0 + 1):(d0.est + ds.est + p0 + q0 + ps)] else 0
		qsvec <- if(qs > 0) eta[(d0.est + ds.est + p0 + q0 + ps + 1):(d0.est + ds.est + p0 + q0 + ps + qs)] else 0
		obj.fun.SARFIMA(d0eta, dseta, p0vec, q0vec, psvec, qsvec, ss, n, x) * 0.5 /n
	}
	grad.in <- function(eta)
	{
		d0eta <- if(d0.est) eta[1] else init.d0
		dseta <- if(ds.est) eta[d0.est + 1] else init.ds
		p0vec <- if(p0 > 0) eta[(d0.est + ds.est + 1):(d0.est + ds.est + p0)] else 0
		q0vec <- if(q0 > 0) eta[(d0.est + ds.est + p0 + 1):(d0.est + ds.est + p0 + q0)] else 0
		psvec <- if(ps > 0) eta[(d0.est + ds.est + p0 + q0 + 1):(d0.est + ds.est + p0 + q0 + ps)] else 0
		qsvec <- if(qs > 0) eta[(d0.est + ds.est + p0 + q0 + ps + 1):(d0.est + ds.est + p0 + q0 + ps + qs)] else 0
		e <- resid.SARFIMA.gen(d0eta, dseta, p0vec, q0vec, psvec, qsvec, ss, n, x)
		deriv.e <- matrix.deriv.e(p0vec, q0vec, psvec, qsvec, c(p0, q0, ps, qs), ss, e, d0.est = d0.est, ds.est = ds.est)
		t(deriv.e) %*% e/n
	}
	gradhess.in <- function(eta)
	{
		d0eta <- if(d0.est) eta[1] else init.d0
		dseta <- if(ds.est) eta[d0.est + 1] else init.ds
		p0vec <- if(p0 > 0) eta[(d0.est + ds.est + 1):(d0.est + ds.est + p0)] else 0
		q0vec <- if(q0 > 0) eta[(d0.est + ds.est + p0 + 1):(d0.est + ds.est + p0 + q0)] else 0
		psvec <- if(ps > 0) eta[(d0.est + ds.est + p0 + q0 + 1):(d0.est + ds.est + p0 + q0 + ps)] else 0
		qsvec <- if(qs > 0) eta[(d0.est + ds.est + p0 + q0 + ps + 1):(d0.est + ds.est + p0 + q0 + ps + qs)] else 0
		e <- resid.SARFIMA.gen(d0eta, dseta, p0vec, q0vec, psvec, qsvec, ss, n, x)
		deriv.e <- matrix.deriv.e(p0vec, q0vec, psvec, qsvec, c(p0, q0, ps, qs), ss, e, d0.est = d0.est, ds.est = ds.est)
		list(gradient = t(deriv.e) %*% e/n, hessian = nlminb.hessian.seq(t(deriv.e) %*% deriv.e/n))
	}
	startvec <- c(init.d0.n, init.ds.n, init.sarma)
	etafit <- nlminb(start = startvec, objective = obj.in, gradient = gradhess.in, hessian = T, scale = scale, control = control, lower = lower, upper = upper)
	if(!conv.test(etafit)) etafit <- nlminb(start = startvec, objective = obj.in, gradient = grad.in, scale = scale, control = control, lower = lower, upper = upper)
	if(!conv.test(etafit)) etafit <- nlminb(start = startvec, objective = obj.in, scale = scale, control = control, lower = lower, upper = upper)
	etaest <- etafit$parameters
	d0eta <- if(d0.est) etaest[1] else NULL
	dseta <- if(ds.est) etaest[d0.est + 1] else NULL
	p0vec <- if(p0 > 0) etaest[(d0.est + ds.est + 1):(d0.est + ds.est + p0)] else NULL
	q0vec <- if(q0 > 0) etaest[(d0.est + ds.est + p0 + 1):(d0.est + ds.est + p0 + q0)] else NULL
	psvec <- if(ps > 0) etaest[(d0.est + ds.est + p0 + q0 + 1):(d0.est + ds.est + p0 + q0 + ps)] else NULL
	qsvec <- if(qs > 0) etaest[(d0.est + ds.est + p0 + q0 + ps + 1):(d0.est + ds.est + p0 + q0 + ps + qs)] else NULL
	estlist <- list(d0 = d0eta, ds = dseta, ar = p0vec, ma = q0vec, sar = psvec, sma = qsvec)
	list(estimates = estlist, converged = conv.test(etafit), conv.type = etafit$message, hessian = etafit$hessian)
}
 


##### MAIN FUNCTION: sarfima.css #################################################################################
# USAGE: 
# sarfima.css(x, model, xreg=0, d0.start=0, ds.start=0, sarma.start=NULL, output.summary=T, 
#	scale=1, control=NULL, lower = -Inf, upper = Inf)
#
# REQUIRED ARGUMENTS: 
# x 		a univariate time series or a vector. 
# model		a list specifying a SARFIMA model. The model should be a list with 
#		SARFIMA(p0,d0,q0)(ps,ds,qs) where the model is specified as follows:
# model <- list(list(order=c(p0,m0,q0)), list(order=c(ps,ms,qs), period=ss), list(d0.est=T, ds.est=T))
# (p0,q0,ps,qs)	order of SARMA parameter vector.
# m0		the number of differences.
# ms		the number of seasonal differences.
# period	an even integer (e.g., ss=4,12).
# 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.  
#
# OPTIONAL ARGUMENTS: 
# xreg		an estimates of mean of x. default is 0 (e.g., xreg=mean(x)).
# d0.start	initial value or preassigned value of d0. default is 0. 
# ds.start	initial value or preassigned value of ds. default is 0. 
# sarma.start   initial values of SARMA parameters.  default is NULL.  
# scale		optional arguments of nlminb
# cotrol	optional arguments of nlminb
# lower, upper	optional arguments of nlminb
#
# VALUE: 
# a list containing the following elements:
# model 	same as input.
# xreg 	same as input.
# estimates	estimates of SARFIMA parameters (d0, ds, SARMA parametes, and an innovation's variance)
# se.coef	standard error, upper 95% confidence limits,
#		and 95% confidence limits of estimates of SARFIMA parameters.
# var.coef 	covariance matrix of the parameter estimates. 
# convergence$converged
#	 	if the optimizer, nlminb, has apparently successfully converged to a minimum, 
#		then converged is TRUE; otherwise converged is FALSE. 
# convergence$conv.type
#	 	a character string describing the type of convergence by nlminb.
# loglik$loglik	value of loglikelihood, namely CSS function. 
# loglik$AIC 	AIC by CSS estimates. AIC is given by -2 * loglik + 2 * np. 
# loglik$BIC	BIC by CSS estimates. BIC is given by -2 * loglik + log(n.used) * np. 
# stinv$d0ds	If estimates d0 and ds satisfy stationary and invertible conditions, then stinv$d0ds is TRUE; 
#		otherwise stinv$d0ds is FALSE.
# stinv$sarma	If estimates of SARMA prameters satisfy stationary and invertible conditions, 
#		then stinv$sarma is TRUE; otherwise stinv$sarma is FALSE.
# n.used        the number of observations used to compute the CSS function. 
# np		the number of parameters. 
#
# WARNING: 	The mean of the series will not be estimated by sarfima.css unless you use the 
#		mean argument. sarfima.css assumes a zero mean series. 
# 		If stinv$sarma is FALSE, se.coef = NA and var.coef = NA.
##################################################################################################################


nullto0 <- function(a, b = 0)
{
	a[is.null(a)] <- b
	a
}


ev.SARMA <- function(pvec, period=NULL)
{
	p <- length(pvec)
	if(p == 1)
		m <- pvec
	else {
		if(!is.null(period)) {
		pvec <- SARcoefs.gen(-pvec, period, p*period +1)[-1]
		p <- p*period
		}
		matA <- matrix(0, p, p)
		matA[1,  ] <- pvec
		matA[row(matA) - col(matA) == 1] <- 1
		m <- eigen(matA)$values
	}
	m
}


### TRUE: there are equal elements between xvec and yvec.
### FALSE: there is no equal elements between xvec and yvec.

equal.elemet.test <- function(xvec, yvec)
{
	equal.test <- length(unlist(lapply(yvec, function(y, xvec)
	xvec[xvec == y], xvec = xvec)))
	if(is.zero(equal.test))
		FALSE  
	else TRUE
}




SARMA.test <- function(p0vec, q0vec, psvec, qsvec, p0, q0, ps, qs, ss)
{
	if(all(is.zero(p0), is.zero(q0), is.zero(ps), is.zero(qs)))
		TRUE
	else {
		e.p0vec <- if(p0 > 0) ev.SARMA(p0vec) else NULL
		e.q0vec <- if(q0 > 0) ev.SARMA(q0vec) else NULL
		e.psvec <- if(ps > 0) ev.SARMA(psvec, period = ss) else NULL
		e.qsvec <- if(qs > 0) ev.SARMA(qsvec, period = ss) else NULL
		maxev.test <- (max(Mod(c(e.p0vec, e.q0vec, e.psvec, e.qsvec))) < 1)
		arma.test <- if(all(p0 > 0, q0 > 0)) (!equal.elemet.test(e.p0vec, e.q0vec)) else TRUE
		sarsma.test <- if(all(ps > 0, qs > 0)) (!equal.elemet.test(e.psvec, e.qsvec)) else TRUE
		all(maxev.test, arma.test, sarsma.test)
	}
}




### extract non-zero estimates ###

sarfima.est.extract <- function(estvec, order, other.name=NULL, d0.est = T, ds.est = T)
{
	p0 <- order[1]
	q0 <- order[2]
	ps <- order[3]
	qs <- order[4]
	np <- d0.est + ds.est + p0 + q0 + ps + qs
	d0name <- if(d0.est) "d0" else  NULL
	dsname <- if(ds.est) "ds" else NULL
	ARname <- if(p0 > 0) paste("ar", 1:p0, sep="") else NULL
	MAname <- if(q0 > 0) paste("ma", 1:q0, sep="") else NULL
	SARname <- if(ps > 0) paste("sar", 1:ps, sep="") else NULL
	SMAname <- if(qs > 0) paste("sma", 1:qs, sep="") else NULL
	resultmat <- matrix(estvec, 1, np+2)
	dimnames(resultmat) <- list("estimates", c(d0name, dsname, ARname, MAname, SARname, SMAname, other.name))
	resultmat
}


### cat function ### 

sarfima.css.summary <- function(z)
{
sarimaorder <- data.frame(list(
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], 
period=z$model[[2]]$period, 
d0.est=z$model[[3]]$d0.est, ds.est=z$model[[3]]$ds.est, xreg0=is.zero(z$xreg), n.used=z$n.used, np=z$np), row.names="")
cat("\n","# sarfima model #","\n")
print(sarimaorder)
cat("\n","# estimates #","\n")
print(
sarfima.est.extract(unlist(z$estimates), 
c(z$model[[1]]$order[1], 
z$model[[1]]$order[3], 
z$model[[2]]$order[1], 
z$model[[2]]$order[3]), 
other.name=c("sigma2", "sigma"), 
d0.est=z$model[[3]]$d0.est, 
ds.est=z$model[[3]]$ds.est))
cat("\n","# stationary and invertible #","\n")
cat(" d0 & ds:", z$stinv$d0ds,  "   sarma:", z$stinv$sarma,"\n")
cat("\n","# 95% confidence limits and standard error of estimates #","\n")
print(z$se.coef)
cat("\n","# variance-covariance matrix of estimates #","\n")
print(z$var.coef)
cat("\n","# convergence #","\n")
cat(" converged:",z$convergence$converged,  "   type of convergence:",z$convergence$conv.type, "\n")
cat("\n","# information criterion #","\n")
cat(" AIC:", round(z$loglik$AIC, digits=5), "   BIC:", round(z$loglik$BIC, digits=5), "\n","\n")
invisible(z)
}



sarfima.css <- function(x, model, xreg=0, d0.start=0, ds.start=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
	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
	m0 <- arima.order[2]
	ms <- seasonal.order[2]
	assign("d0", d0.start, frame = 1)
	assign("ds", ds.start, frame = 1)
	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 # mean zero process
	assign("x", x, frame = 1)
	assign("n", length(x), frame = 1)
	np <- p0 + q0 + ps + qs + d0.est + ds.est + 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	### racf and Chi-square ###
	### outputs of no estimator case ###
		resultlist <- list(model = model, xreg = xreg, estimates = list(d0 = NULL, ds = NULL, ar = NULL, ma = NULL, sar = NULL, sma = NULL, sigma2 = sigma2, sigma = sqrt(sigma2)), var.coef=NULL, converged = T, conv.type = NULL, n.used = n, loglik = list(loglik = loglik, AIC = aic, BIC = bic), initial.value=list(d0ds=c(d0.start,ds.start), sarma=sarma.start), np=np )
	} else {
	### calculate initial value ###
	init.d0ds <- initial.d0ds.nlminb(d0, ds, 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, sarma.start=sarma.start)
	### sarfima.nlse ###
	fit <- sarfima.nlse.nlminb(x, init.d0ds, init.sarma, ss, d0.est=d0.est, ds.est=ds.est, scale =scale, control=control, lower = lower, upper = upper)
	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, b=d0)
	n.ds <- nullto0(resds, b=ds)
	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)
	estimates <- list(d0 = resd0, ds = resds, ar = resar, ma = resma, sar = ressar, sma = ressma, sigma2 = sigma2, sigma=sigma)
	### variane matrix and 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) {
	var.coef <- sarfima.css.var(resar, resma, ressar, ressma, c(p0,q0,ps,qs), ss, n, d0.est = d0.est, ds.est = ds.est)
	estvec <- unlist(list(estimates$d0,estimates$ds,estimates$ar,estimates$ma,estimates$sar, estimates$sma))
	se.coefvec <- sqrt(diag(var.coef))
	upp95se <- estvec + qnorm(0.975)*se.coefvec
	low95se <- estvec - qnorm(0.975)*se.coefvec
	se.coef <- matrix(c(upp95se, low95se, se.coefvec), nrow=3, byrow=T)
	dimnames(se.coef) <-  list(c("upper95","lower95","se.coef"), dimnames(var.coef)[[1]])}
	else {
	warning("Warning: estimated SARMA coefficients may be non-stationary or non-invertible")
	var.coef <- NA
	se.coef <- NA
	}
	testd0ds <- abs(n.d0 + n.ds)
	testds <- abs(n.ds)
	d0ds.sta.inv <- ((testd0ds > .5)|(n.ds > .5))
	if(d0ds.sta.inv) 
		warning("Warning: estimated d0 and ds are non-stationary or non-invertible")
	### Statistics ###
	loglik <-  - n * 0.5 * log(2 * pi * sigma2 + 1)
	aic <- -2 * loglik + 2 * np
	bic <- -2 * loglik + log(n) * np
	### outputs ###
	resultlist <- list(model=model, xreg=xreg, estimates = estimates, 
	se.coef=se.coef, var.coef=var.coef, 
	convergence =list(converged = fit$converged, conv.type = fit$conv.type), 
	loglik = list(loglik = loglik, AIC = aic, BIC = bic), 
	initial.value=list(d0ds=c(d0.start,ds.start), sarma=sarma.start), 
	stinv = list(d0ds = (!d0ds.sta.inv), sarma = sarma.sta.inv),
 	n.used=n, np=np   )
	}
	if(output.summary){
	sarfima.css.summary(resultlist)
	invisible(resultlist)
	}
	else resultlist
}





### example 1 ################################################################################
### ARIMA(1,1,1) with ar=-0.5 and ma=-0.5
# data <- arima.sim(100,model=list(ar=-.5, ndiff=1, ma=-.5)) 
### model 1 
### Fitting a SARFIMA model, SARFIMA(1,d0,1)(1,ds,1)s with m0=1, ms=0
# model1 <- list(list(order=c(1,1,1)), list(order=c(1,0,1), period=12), list(d0.est=T, ds.est=T))
# fit1 <- sarfima.css(data, model1)
# fit1$estimates
# sarfima.css.summary(fit1)
### model 2 
### Fitting ARIMA(1,1,1) by arima.mle
# fit2 <- arima.mle(data, list(list(order=c(1,1,1)), list(order=c(1,0,1), period=12)))$model
# fit2
### model 3 
### Fitting a SARFIMA model, SARFIMA(1,d0,1)(0,0,0)s with m0=1, ms=0
# model3 <- list(list(order=c(1,1,1)), list(order=c(0,0,0), period=12), list(d0.est=T, ds.est=F))
# fit3 <- sarfima.css(data, model3)
### model 4 
### Fitting a SARFIMA model, SARFIMA(1,0,1)(0,0,0)s with m0=1, ms=0, namely, ARIMA(1,1,1)
# model4 <- list(list(order=c(1,1,1)), list(order=c(0,0,0), period=12), list(d0.est=F, ds.est=F))
# fit4 <- sarfima.css(data, model4)
##############################################################################################



### example 2 ################################################################################
### SARIMA(1,1,1)(1,1,1)s with ar=0.3, ma=0.5, sar=-0.3, sma=-0.5, and period=12
# data <- arima.sim(500, model=list(list(ar=.3, ndiff=1, ma=.5), list(ar=-.3, ndiff=1, ma=-.5, period=12)) )[101:500]
### model 1
### Fitting a SARFIMA model, SARFIMA(2,d0,1)(1,ds,1)s with m0=1, ms=1
# model1 <- list(list(order=c(2,1,1)), list(order=c(1,1,1), period=12), list(d0.est=T, ds.est=T))
# fit1 <- sarfima.css(data, model1)
### model 2
### Fitting a SARFIMA model, SARFIMA(2,0,1)(1,0,1)s with m0=1, ms=1, namely, SARIMA(2,1,1)(1,1,1)s
# model2 <- list(list(order=c(2,1,1)), list(order=c(1,1,1), period=12), list(d0.est=F, ds.est=F))
# fit2<- sarfima.css(data, model2)
### model 3
### Fitting a SARFIMA model, SARFIMA(1,0,1)(1,0,1)s with m0=1, ms=1, namely, SARIMA(1,1,1)(1,1,1)s
# model3 <- list(list(order=c(1,1,1)), list(order=c(1,1,1), period=12), list(d0.est=F, ds.est=F))
# fit3<- sarfima.css(data, model3)
### model 4
### Fitting SARIMA(2,1,1)(1,1,1)s by arima.mle
# arima.mle(data, list(list(order=c(1,1,1)), list(order=c(1,1,1), period=12)))$model
##############################################################################################