###################################################################
### sarfima.css.diag: compute diagnostics for the SARFIMA Model ###
### This program uses new S-PLUS function "stdev", which        ###
### calculates standard deviation.                              ###
###################################################################

### sample autocorrelation functions for mean zero series
meanzero.acf <- function(z,maxlag)
{
n <- length(z)
maxlag <- maxlag+1
m <- rep(c(z,rep(0,maxlag)), maxlag)[1:((n+maxlag-1)*maxlag)]
dim(m) <- c(n+maxlag-1, maxlag)
result <- z %*% m[1:n,  ]
c(result/result[1])
}




### inverted filter's coefficients  
### uses n series of the coefficients 

inv.coef.gen <- function(coef)
{
	n <- length(coef)
	inv <- vector(length = n)
	inv[1] <- 1
	if(n > 1) {
		for(j in 1:(n - 1)) {
			inv[j + 1] <-  - t(coef[2:(j + 1)]) %*% inv[j:1]
		}
	}
	drop(inv)
}





### matrix, Jm #########################################
# m = h * s : nubmer of row of matrix Jm
########################################################

matrix.Jm <- function(p0vec, q0vec, psvec, qsvec, order, h, ss, d0.est = T, ds.est = T)
{
	p0 <- order[1]
	q0 <- order[2]
	ps <- order[3]
	qs <- order[4]
	colJm <- d0.est + ds.est + p0 + q0 + ps + qs
	m <- h * ss
	if(d0.est == T) {
		d0vec <- 1/(1:m)
		d0name <- "d0"
	}
	else {
		d0vec <- NULL
		d0name <- NULL
	}
	if(ds.est == T) {
		dsvec <- rep(0, m)
		dsvec[ss * (1:h)] <- 1/(1:h)
		dsname <- "ds"
	}
	else {
		dsvec <- NULL
		dsname <- NULL
	}
	if(p0 > 0) {
		ARmat <- matrix(rep(c(inv.coef.gen(ARcoefs.gen(p0vec, m)), 0), p0)[1:(m * p0)], m, p0)
		ARmat[col(ARmat) - row(ARmat) > 0] <- 0
		ARname <- paste("ar",1:p0, sep="") 
	}
	else {
		ARmat <- NULL
		ARname <- NULL
	}	
	if(q0 > 0) {
		MAmat <-  - matrix(rep(c(inv.coef.gen(MAcoefs.gen(q0vec, m)), 0), q0)[1:(m * q0)], m, q0)
		MAmat[col(MAmat) - row(MAmat) > 0] <- 0
		MAname <- paste("ma",1:q0, sep="") 
	}
	else {
		MAmat <- NULL
		MAname <- NULL
	}
	if(ps > 0) {
		SARmat <- matrix(rep(c(rep(0, ss - 1), inv.coef.gen(SARcoefs.gen(psvec, ss, m - ss + 1)), rep(0, ss)), ps)[1:(m * ps)], m, ps)
		SARmat[ss * col(SARmat) - row(SARmat) > 0] <- 0
		SARname <- paste("sar",1:ps, sep="") 
	}
	else {
		SARmat <- NULL
		SARname <- NULL
	}
	if(qs > 0) {
		SMAmat <-  - matrix(rep(c(rep(0, ss - 1), inv.coef.gen(SMAcoefs.gen(qsvec, ss, m - ss + 1)), rep(0, ss)), qs)[1:(m * qs)], m, qs)
		SMAmat[ss * col(SMAmat) - row(SMAmat) > 0] <- 0
		SMAname <- paste("sma",1:qs, sep="") 
	}
	else {
		SMAmat <- NULL
		SMAname <- NULL
	}
	resultmat <- matrix(c(d0vec, dsvec, as.vector(ARmat), as.vector(MAmat), as.vector(SARmat), as.vector(SMAmat)), m, colJm)
	dimnames(resultmat) <- list(NULL, c(d0name, dsname, ARname, MAname, SARname, SMAname))
	resultmat
}


# example 
# matrix.Jm(0.1, 0 ,0, 0, c(1,0,0,0), 5, 12, d0.est=T, ds.est=F)


sarfima.css.Fisher <- function(p0vec, q0vec, psvec, qsvec, order, ss, d0.est = T, ds.est = T, M =50)
{
	if((d0.est == T) && (ds.est == T)) {
		Jm <- matrix.Jm(p0vec, q0vec, psvec, qsvec, order, M, ss, d0.est = T, ds.est = T)
		Fisher <- t(Jm) %*% Jm
		Fisher[1:2, 1:2] <- matrix(c(pi^2/6, pi^2/(6 * ss), pi^2/(6 * ss), pi^2/6), 2, 2)
	}
	else if((d0.est == T) && (ds.est == F)) {
		Jm <- matrix.Jm(p0vec, q0vec, psvec, qsvec, order, M, ss, d0.est = T, ds.est = F)
		Fisher <- t(Jm) %*% Jm
		Fisher[1, 1] <- pi^2/6
	}
	else if((d0.est == F) && (ds.est == T)) {
		Jm <- matrix.Jm(p0vec, q0vec, psvec, qsvec, order, M, ss, d0.est = F, ds.est = T)
		Fisher <- t(Jm) %*% Jm
		Fisher[1, 1] <- pi^2/6
	}
	else {
		Jm <- matrix.Jm(p0vec, q0vec, psvec, qsvec, order, M, ss, d0.est = F, ds.est = F)
		Fisher <- t(Jm) %*% Jm
	}
	Fisher
}


### variance-covariance matrix of CSS estimates

sarfima.css.var <- function(p0vec, q0vec, psvec, qsvec, order, ss, n, d0.est = T, ds.est = T, M =50)
{
	Fisher <- sarfima.css.Fisher(p0vec, q0vec, psvec, qsvec, order, ss, d0.est = d0.est, ds.est = ds.est, M =M)
	solve(Fisher)/n
}





### MAIN PROGRAM: sarfima.css.diag #############################################################
# compute diagnostics for the SARFIMA Model: 
# plot std residuals, residual autocorrelations tests, and modified portmateau tests. 
#                              
# USAGE: sarfima.css.diag(z, x, maxlag =5, diag.plot=T, invisible.output=T)
#
# REQUIRED ARGUMENTS: 
# z 		a list of the output from sarfima.css
# x 		a univariate time series or a vector. 
#
# OPTIONAL ARGUMENTS: 
# 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)
# diag.plot	logical flag: if TRUE, the diagnostics will be plotted. 
# invisible.output
#		logical flag: if TRUE, values will be hidden. 
#
# VALUE: 
# a list containing the following elements:
# resid		residuals
# std.resid	standardized residuals
# racf$lag 	lag of RACF
# racf$racf	RACF
# racf$pvale	p-value of test statistics by RACF, which test residuals are 
# 		approximately white noise sequences. 
# racf$rejlag	rejected lags of RACF tests with 5% siginificance level.
# pacf$lag 	lag of PACF of residuals
# pacf$pacf	PACF of residuals
# chisquare$df	degrees of freedom of the Chi-squared test statistics
# chisquare$chisquare
#		the Chi-squared test statistics
# chisquare$pvalue 
# 		p-value of the Chi-squared test statistics
# chisquare$rejdf 
# 		rejected dfs of Chi-squared tests with 5% siginificance level.
#################################################################################################

sarfima.css.diag <- function(z, x, maxlag =5, diag.plot=T, invisible.output=T)
{
	arima.order <- z$model[[1]]$order
	seasonal.order <- z$model[[2]]$order
	p0 <- arima.order[1]
	m0 <- arima.order[2]
	q0 <- arima.order[3]
	ps <- seasonal.order[1]
	ms <- seasonal.order[2]
	qs <- seasonal.order[3]
	ss <- z$model[[2]]$period
	d0.est <- z$model[[3]]$d0.est
	ds.est <- z$model[[3]]$ds.est
	d0 <- z$initial.value$d0ds[1]
	ds <- z$initial.value$d0ds[2]
	xreg <- z$xreg
	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 <- z$n.used
	np <- p0 + q0 + ps + qs + d0.est + ds.est
	m <- max(min(maxlag * ss, floor(n/3)), ss + 1, np+1)
	### estimates ###
	resd0 <- z$estimates$d0
	resds <- z$estimates$ds
	resar <- z$estimates$ar
	resma <- z$estimates$ma
	ressar <- z$estimates$sar
	ressma <- z$estimates$sma
	var.coef <- z$var.coef
	### residual ###
	residual <- resid.SARFIMA.gen(nullto0(resd0, b=d0), nullto0(resds, b=ds), nullto0(resar), nullto0(resma), nullto0(ressar), nullto0(ressma), ss, n, x)
	std.resid <- residual/stdev(residual)
	racf <- meanzero.acf(residual, m) 
	if(is.zero(np)) std.err.racf <- sqrt( (n-c(1:m))/(n*(n+2)) )
	else {
	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 ) 
	}
	pvalue.racf <- 2*(1 - pnorm( abs(racf[2:(m+1)])/std.err.racf ))
	upp95bound.racf <- std.err.racf * qnorm(0.975)
	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)) )
	pacf.resid <- c(acf(residual, type="partial", lag=m, plot=F)$acf)
	rej.lag.racf <- if(any(pvalue.racf<0.05)) c(1:m)[pvalue.racf<0.05] else NULL
	rej.lag.chi <- if(any(pvalue.chi<0.05)) c(1:m)[pvalue.chi<0.05] else NULL
	### plot (4 graphs) ###
	cexv <- 0.9 
	if(diag.plot){
	par( mfrow = c(4,1)  )
	### Plot of Standardized Residuals ###
	nn <- length(std.resid)
	plot(1:nn, std.resid, type='l', xlab="", ylab="", main="Plot of Standardized Residuals", cex=cexv  ) #alima.mle ‚Ĺ‚Ítype='h'
	lines(-1:(nn+1),rep(0,(nn+3)),type='l')
	### ACF Plot of Residuals ###
	ylimlow <- max( min(racf, -upp95bound.racf) - 0.1, -1.0 )
	if(is.null(rej.lag.racf)) xlabracf <- "lag" 
	else {
	rej.lag.racf.plot <- if(length(rej.lag.racf)>10) paste(c(paste(rej.lag.racf[1:10]), "..."), collapse=" ") else paste(rej.lag.racf)
	xlabracf <- paste(c("lag", "  ( rejected 5% :",  rej.lag.racf.plot, ")"), collapse=" ")
	}
	plot(0:m, racf,  type='h', xlab=xlabracf, ylab="", ylim=c(ylimlow,1), main="ACF Plot of Residuals", cex=cexv  )
	lines(-1:(m+1),rep(0,(m+3)),type='l')
	lines(1:m, upp95bound.racf, lty=2) # UCL 95%
	lines(1:m, -upp95bound.racf, lty=2) # LCL 95 %
	points(1:m, pvalue.racf, cex=cexv)
	### PACF Plot of Residuals ###
	plot( 1:m, pacf.resid, type='h', xlab="lag", ylab="", main="PACF Plot of Residuals" , cex=cexv )
	lines(-1:(m+1),rep(0,(m+3)),type='l')
	### Ljung-Box Chi-Squared Statistics ###
	if(is.null(rej.lag.chi)) xlabchi <- "df" 
	else {
	rej.lag.chi.plot <- if(length(rej.lag.chi)>10) paste(c(paste(rej.lag.chi[1:10]), "..."), collapse=" ") else paste(rej.lag.chi)
	xlabchi <- paste(c("df", "  ( rejected 5% :",  rej.lag.chi.plot, ")"), collapse=" ")
	}
	plot(1:(m-np), pvalue.chi, xlab=xlabchi, ylab="", main="P-values of Ljung-Box Chi-Squared Statistics", cex=cexv)
	}
	### list ###
	result <- list(resid=residual, std.resid= std.resid, racf=list(lag=c(1:m), racf=racf[2:(m+1)], pvalue=pvalue.racf, rejlag=rej.lag.racf), pacf=list(lag=c(1:m), pacf=pacf.resid), chisquare = list(df=c(1:(m-np)), chisquare=mod.Portmanteau[(np+1):m], pvalue=pvalue.chi, rejdf=rej.lag.chi) )
	if(invisible.output) invisible(result) else result
}




### 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), list(d0.est=T, ds.est=F))
 fit1 <- sarfima.css(data, model1)
 sarfima.css.diag(fit1, data)
 model2 <- list(list(order=c(0,0,0)), list(order=c(0,0,0), period=12), list(d0.est=F, ds.est=T))
 fit2 <- sarfima.css(data, model2)
 sarfima.css.diag(fit2, data)
##############################################################################################