######################################################################
#INPUT
# evec  : vector of (e_1, e_2, ... , e_k). upper limit of range of k-multiple integration
#	e_1 =< e_2 =< ... =< e_k 	 
# dfvec : vector of (d_1, d_2, ... , d_k). degrees of freedom of portmanteau statisitcs 
# 	d_1 =< d_2 =< ... =< d_k and even integers
# 	k = length(evec) = length(dfvec) > 1 
#######################################################################


# convolution function 
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
}



# subprograms
Bpnp <- function(mjvec, djvec)
{
	p <- length(djvec)
	n1 <- sum(mjvec + 1)
	njvec <- n1 - c(0, cumsum(mjvec + 1)[1:p - 1])
	B <- matrix(0, p, n1 + 1)
	n2vec <- c(0:njvec[2])
	b2i <- (factorial(mjvec[1]) * djvec[1]^(mjvec[1] + 1 + n2vec))/factorial(mjvec[1] + 1 + n2vec)
	b2ji <- (djvec[2] - djvec[1])^n2vec/factorial(n2vec)
	B[2, n2vec + 1] <- conv(b2i, b2ji)[n2vec + 1] * factorial(n2vec)
	if(p == 2)
		B
	else {
		for(k in 3:p) {
			nkvec <- c(0:njvec[k])
			bki <- (factorial(mjvec[k - 1]) * B[k - 1, mjvec[k - 1] + 1 + nkvec + 1])/factorial(mjvec[k - 1] + 1 + nkvec)
			bkji <- (djvec[k] - djvec[k - 1])^nkvec/factorial(nkvec)
			B[k, nkvec + 1] <- conv(bki, bkji)[nkvec + 1] * factorial(nkvec)
		}
		B
	}
}


Apnp <- function(mjvec, djvec)
{
	p <- length(djvec)
	n1 <- sum(mjvec + 1)
	njvec <- n1 - c(0, cumsum(mjvec + 1)[1:p - 1])
	A <- matrix(0, p, n1 + 1)
	B <- Bpnp(mjvec, djvec)
	d1 <- djvec[1]
	A[1, 1] <- -2 * (exp(-0.5 * d1) - 1)
	for(j in 1:n1) {
		A[1, j + 1] <- -2 * exp(-0.5 * d1) * d1^j + 2 * j * A[1, j]
	}
	for(k in 2:p) {
		A[k, 1] <- -2 * exp(-0.5 * djvec[k]) * B[k, 1] + 2 * A[k - 1, mjvec[k - 1] + 1]
		nk <- njvec[k]
		for(j in 1:nk) {
			A[k, j + 1] <- -2 * exp(-0.5 * djvec[k]) * B[k, j + 1] + 2 * j * A[k, j]
		}
	}
	A
}



############################################################################
# MAIN PROGRAM 
# p	: positve integer larger than 2. 
# dfvec : p-vector of dfs of portmanteau statistics where dfs are even integers 
#	and dfvec[1]<dfvec[2]<...<dfvec[p]
# evec  : p-vector of critical values of the portmanteau statistics where 
#	0<=evec[1]<=evec[2]<=...<=evec[p]
############################################################################

############################################################################
# jcdfport : function to calculate joint cumulative distribution function 
#	 (j.c.d.f.) of the  portmanteau statistics, i.e. 
#        Pr( chi^2(dfvec[1])<=evec[1], ... , chi^2(dfvec[p])<=evec[p] )
############################################################################

jcdfport <- function(dfvec, evec)
{
	mjvec <- c(dfvec[1], diff(dfvec))/2 - 1
	p <- length(dfvec)
	if(p == 1)
		pchisq(evec, df = dfvec)
	else Apnp(mjvec, evec)[p, mjvec[p] + 1]/prod((2^(mjvec + 1)) * factorial(mjvec))
}


############################################################################
# cjcdfport : function to calculate conditional j.c.d.f of the 
#	 portmanteau statistics, i.e. 
# Pr( chi2(dfvec[1])<=evec[1], ... , chi2(dfvec[p-1])<=evec[p-1] | chi2(dfvec[p])==evec[p] )
############################################################################

cjcdfport <- function(dfvec, evec)
{
	mjvec <- c(dfvec[1], diff(dfvec))/2 - 1
	p <- length(dfvec)
	hvec <- rep(evec[p], p)
	Bpnp(mjvec, evec)[p, mjvec[p] + 1]/Bpnp(mjvec, hvec)[p, mjvec[p] + 1]
}


############################################################################
# betaport : function to calculate rejection value and beta of the 
#	 portmanteau statistics given alpha, i.e. 
# INPUT: 
# alpha : significance level of multiple portamnateaus tests
# dfvec : p-vector of dfs of portmanteau statistics where dfs are even integers 
#	and dfvec[1]<dfvec[2]<...<dfvec[p]
# OUTPUT: beta and djbeta
############################################################################

derivbetaport <- function(beta, dfvec)
{
	djbeta <- qchisq(1 - beta, dfvec)
	p <- length(dfvec)
	probvec <- vector(length = p)
	dfdiff <- c(dfvec[1], diff(dfvec))
	if(p == 2) {
		probvec[1] <- pchisq(djbeta[p] - djbeta[1], df = dfdiff[p])
		probvec[p] <- cjcdfport(dfvec, djbeta) 
	}
	else {
		djmd1 <- c(djbeta - djbeta[1])
		probvec[1] <- jcdfport(cumsum(dfdiff[-1]), djmd1[-1])
		probvec[p] <- cjcdfport(dfvec, djbeta)
		for(k in 2:(p - 1)) {
			probvec[k] <- cjcdfport(dfvec[1:k], djbeta[1:k]) * jcdfport(
				cumsum(dfdiff[(k + 1):p]), djmd1[(k + 1):p])
		}
	}
 - sum(probvec)
}

# test 
derivbetaport(0.3, c(2,4,8))


betaport <- function(alpha, dfvec, conv=0.1^5, maxit=100)
{
	mjvec <- c(dfvec[1], diff(dfvec))/2 - 1
	p <- length(dfvec)
	beta <- alpha
	djbeta <- qchisq(1 - beta, df=dfvec)
	numer <- jcdfport(dfvec, djbeta) - (1 - alpha)
	denom <- derivbetaport(beta, dfvec)
	bias <- numer/denom
	it <- 0
	while(max(abs(bias)) > conv && (it <- it + 1) < maxit) {
		numer <- jcdfport(dfvec, djbeta) - (1 - alpha)
		denom <- derivbetaport(beta, dfvec)
		bias <- numer/denom
		beta <- beta - bias
		djbeta <- qchisq(1 - beta, df=dfvec)
	}
	list(beta = c(min(max(0, beta), 1)), evec = c(djbeta))
}
 

### EXAMPLE: Table 1 DF(2,4,6) #####
dfvec <- c(2,6,10) 
betaport(0.01, dfvec) # alpha =0.01 is given.
betaport(0.05, dfvec) # alpha =0.05 is given.
betaport(0.10, dfvec) # alpha =0.10 is given.

beta <- 0.01
divector01 <- qchisq(1-beta, df=dfvec)
list(alpha=1 - jcdfport(dfvec, divector01), beta=beta, divec=divector01)
beta <- 0.05
divector05 <- qchisq(1-beta, df=dfvec)
list(alpha=1 - jcdfport(dfvec, divector05), beta=beta, divec=divector05)
beta <- 0.10
divector10 <- qchisq(1-beta, df=dfvec)
list(alpha=1 - jcdfport(dfvec, divector10), beta=beta, divec=divector10)