#############################################
##################################
#####   This file contains the functions:
#####   1.cepstral.coef
#####   2.lik3.samples
#####   3.der.lik3.samples
#####   4.Fishinfmatr
#####   5.st.dev
#####   6.correct4zeros
#####
#####   These functions are used to estimate the cepstral coefficients
#######################################################

#############################################
### The function cepstral.coef calculates the maximum likelihood
### estimator of theta, a vector that corresponds to the cepstral
### coefficients. It also calculates an estimate for the standard deviation
### of these coefficients, a value for the likelihood ratio test, the p-value
### for the test, the value of the Akaike's Information Criterion and an indicator of
### whether the method has successfully converged.
###
### Required input:
### 1.data should be a GxN matrix, G being the number of time series and N
###   the number of observations for each time series.
### 2. p should be the order of the model used to fit the data.
### 3. nspans should be the width of the modified Daniell smoother to
###    be used to smooth the periodogram.
###
#############################################
cepstral.coef<-function(data,p,nspans)
{
	# Number of time series
	G <- dim(data)[1]

	# Number of observations for each time series
	N <- dim(data)[2]

	# Vector for te cepstral coefficients
	coef <- rep(NA, (G-1)*(p+1))

	# The value of the likelihood ratio test
	dev.new <- rep(NA,1)

	# Indicator for convergence. 0 indicates successful convergence
	message <- rep(NA,1)

	# The number of values disregarded in the periodogram. This should be zero!!!
	ntaper <- 0

	# The matrix for the estimates of the spectral density
	ffts <-matrix(NA, nrow=G, ncol=floor(N/2)-1)

	for(i in 1:G)
	{
		ffts[i,] <- spectrum(data[i,],spans=nspans,taper=ntaper,plot=F, fast=F)$spec[-floor(N/2)]
		ffts[i,] <- correct4zeros(ffts[i,])
	}

	# The matrix for the values of the periodogram at the Fourier frequencies
	# for each times series
	X<-matrix(NA, nrow=p, ncol=floor(N/2)-1)

	for (i in 1:p)
	{
		X[i,] <-2*cos((2*i)*pi*spectrum(data[1,], plot=F, fast=F)$freq[-floor(N/2)])
	}

	# The first (G-1) columns of the matrix "stuff" contain the estimates for the ratios of the spectral density
	# of each of the first (G-1) time series to the spectral density of the  last time series
	stuff <-matrix(NA, nrow=floor(N/2)-1, ncol=G-1+p)
	for (j in 1:G-1)
	{
		stuff[,j] <-log(ffts[j,]/ffts[G,])
	}

	# The last p columns of the matrix "stuff" contain the the values
	# of the periodogram at the Fourier frequencies for each times series
	for(k in G:(G-1+p))
	{
		stuff[,k]<-X[k-G+1,]
	}

	# The vector for the initial values of the cepstral coefficients to be used in the function optim
	a<-rep(1,(G-1)*(p+1))

	outputs     <- optim(fn=lik3.samples, par=c(a), data=stuff, numG=G, numN =N, gr=der.lik3.samples, method="BFGS" )
	outputsnew  <- lik3.samples(theta=rep(0, (G-1)*(p+1)),data=stuff,numG=G, numN =N)
	# The estimate for the cepstral coefficients
	coef		  <- outputs$par
	# Indicator for successful convergence
	message	  <- outputs$convergence
	# The value of the likelihood ratio test
	dev.new     <- 2*(outputsnew-outputs$value)
	# The p-value
	pvalue		  <- 1 - pchisq(dev.new,(G-1)*(p+1))
	# The value of the AIC
	AIC 		  <- 2*lik3.samples(coef,data=stuff,numG=G, numN =N)+ 2*(G-1)*(p-1)


	Cepstral_Coef <- matrix(coef,(p+1),(G-1))
	Column.Names <- paste("theta", 1:(G-1), sep = "")
	colnames(Cepstral_Coef)<- Column.Names


	Ceps_Coef.StD <- matrix(st.dev(G,N,p),(p+1),(G-1))
	colnames(Ceps_Coef.StD)<- Column.Names

	return(list(Time_Series_Num=G,Time_Series_Length=N,Order_Of_Model=p,Cepstral_Coefficients=Cepstral_Coef,Ceps_Coef.Standard_Deviation=Ceps_Coef.StD,
				Likelihood_Ratio_Test=dev.new,P_value=pvalue,Akaikes_Information_Criterion=AIC,Convergence=message))

}
#############################################
##################################
#################################
#########Log likelihood function
#########for the general model as presented in the
############paper in Technometrics (Fokianos & Savvides)
############# with  covariates
#######################################################



lik3.samples<-function(theta, data,numG,numN)
{
	G <- numG
	N <- numN
	p <- dim(data)[2]-G+1
	s<-rep(NA,G-1)
	u<-matrix(NA, nrow=floor(N/2)-1, ncol=G-1)
	a<-rep(NA,floor(N/2)-1)
	linear<-matrix(NA, nrow=floor(N/2)-1, ncol=G-1)

	for (i in 1:(G-1))
	{
		d<-matrix(NA, nrow=floor(N/2)-1, ncol=p)

		for (m in 1:p)
		{
			d[,m]<- theta[(i-1)*p+i+m]*data[,(G-1)+m]
		}

		l<-rep(NA,floor(N/2)-1)

		for(k in 1:floor(N/2)-1)
		{
			l[k]<- sum(d[k,])
		}

		linear[,i]<- data[,i]-theta[(i-1)*p+i]-l
		s[i]<--sum(linear[,i])
		u[,i]<- exp(linear[,i])
	}

	for (j in 1:(floor(N/2)-1))
	{
		a[j]<-sum(u[j,])
	}

	ll <-  sum(s)+G*sum(log(1+a))

	return( ll )
}


########################################
########################################
#########Score function
#########for the general model as presented in the
############paper in Technometrics (Fokianos & Savvides)
##############with  covariates
#######################################################



der.lik3.samples <- function(theta, data,numG,numN)
{
	G <- numG
	N <- numN
	p <- dim(data)[2]-G+1
	s<-rep(NA,G-1)
	u<-matrix(NA, nrow=floor(N/2)-1, ncol=G-1)
	a<-rep(NA,floor(N/2)-1)
	linear<-matrix(NA, nrow=floor(N/2)-1, ncol=G-1)

	for (i in 1:(G-1))
	{
		d<-matrix(NA, nrow=floor(N/2)-1, ncol=p)

		for (m in 1:p)
		{
			d[,m]<- theta[(i-1)*p+i+m]*data[,(G-1)+m]
		}

		l<-rep(NA,floor(N/2)-1)

		for(k in 1:floor(N/2)-1)
		{
			l[k]<- sum(d[k,])
		}

		linear[,i]<- data[,i]-theta[(i-1)*p+i]-l

		s[i]<--sum(linear[,i])
		u[,i]<- exp(linear[,i])
	}

	for (j in 1:(floor(N/2)-1))
	{
		a[j]<-sum(u[j,])
	}

	b<-rep(NA,G-1)
	v<-matrix(NA, nrow=G-1, ncol=p)

	for (i in 1:(G-1))
	{
		b[i]<-dim(data)[1] -G*sum( exp(linear[,i])/(1+a))

		for(h in 1:p)
		{
			v[i,h]<-sum(data[,G-1+h])-G*sum(exp(linear[,i])*data[,G-1+h]/(1+a))
		}
	}

	w<-matrix(NA, nrow=G-1, ncol=p+1)

	for(m in 1:(G-1))
	{
		w[m,]<-c(b[m],v[m,])
	}

	A<-t(w)

	return( c(as.vector(A)) )
}

##########################################################
##########################################################
##########################################################

Fishinfmatr<-function(G,N,p)
{
	A<-rep(NA,(p+1)*(p+1))
	dim(A)<-c((p+1),(p+1))
	B<-rep(NA,(p+1)*(p+1))
	dim(B)<-c((p+1),(p+1))
	d<-floor(N/2)
	x<-rep(NA,(p+1)*d)
	dim(x)<-c(p+1,d)
	for(j in 0:p)
	{
		x[j+1,]<-cos(j*2*pi*(1:d)/N)
	}
	a<-(G-1)/(G)
	b<--1/(G)
	A[1,1]<-a*sum(x[1,]*x[1,])
	B[1,1]<-b*sum(x[1,]*x[1,])
	for(m in 2:(p+1))
	{
		A[1,m]<-2*a*sum(x[1,]*x[m,])
		A[m,1]<-2*a*sum(x[m,]*x[1,])
		B[1,m]<-2*b*sum(x[1,]*x[m,])
		B[m,1]<-2*b*sum(x[m,]*x[1,])
	}
	for(l in 2:(p+1))
	{
		for(k in 2:(p+1))
		{
			A[l,k]<-4*a*sum(x[l,]*x[k,])
			B[l,k]<-4*b*sum(x[l,]*x[k,])
		}
	}
	return(list(x,A,B))
}
##################################################################


st.dev <-function(G,N,p)
{
	k <- Fishinfmatr(G,N,p)
	A <- k[[2]]
	B <- k[[3]]

	J <- matrix(1, nrow=G-1, ncol=G-1)
	InfoMatr <- kronecker(diag(1,G-1), A)+kronecker(J-diag(1,G-1), B)
	###################  is the fisher information matrix.
	######################################################################

	return(sqrt(diag(solve(InfoMatr))))
}

#################################################################

correct4zeros<-function(mat)
{
    min2<-0
    smat<-sort(mat)
    if( max(mat) == 0)
    {
        print ("ERROR ... Maximum value equals to zero")
        return(mat)
    }
    for(i in 1:length(smat))
    {
        if( smat[i]>0 )
        {
            min2=smat[i]
            break
        }
    }
    for(i in 1:length(mat))
    {
        if(mat[i]==0) mat[i] = min2
    }
    return(mat)
}

###################################################################
###################################################################