### FIVAR_VARFI.R
## Copyright (c) 2008, Rebecca Sela
## All rights reserved.
##
## Redistribution and use in source and binary forms, with or without
## modification, are permitted provided that the following conditions are met:
##     * Redistributions of source code must retain the above copyright
##       notice, this list of conditions and the following disclaimer.
##     * Redistributions in binary form must reproduce the above copyright
##       notice, this list of conditions and the following disclaimer in the
##       documentation and/or other materials provided with the distribution.
##     * Neither the name of the <organization> nor the
##       names of its contributors may be used to endorse or promote products
##       derived from this software without specific prior written permission.
##
## THIS SOFTWARE IS PROVIDED BY Rebecca Sela ``AS IS'' AND ANY
## EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
## WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
## DISCLAIMED. IN NO EVENT SHALL Rebecca Sela BE LIABLE FOR ANY
## DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
## (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
## LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
## ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
## (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
## SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
##
## This code can be cited by citing the paper:
## "Computationally Efficient Methods for Two Multivariate Fractionally 
## Integrated Models" by Sela and Hurvich.
##
## Additional information (and this code) is available from:
## http://pages.stern.nyu.edu/~rsela/VARFI/code.html
##
## Please report bugs to Rebecca Sela <rsela@stern.nyu.edu>.
##


################ Functions for Estimation ##################
### FIVAR.GPH
## This function computes the GPH estimator of the d parameters in 
## a FIVAR model.  
## Inputs:
##	X - KTx1 vector of observations, grouped by series
##	K - number of series
##	numPoints [T^0.4] - number of points to use in the log periodogram regression
##	verbose [FALSE] - should the regression output be printed?
## Output: vector of d estimates
FIVAR.GPH <- function(X, K, numPoints=T^0.4, verbose=FALSE){
    # Compute T
    T <- length(X)/K

    # Compute cross-periodogram of the series
    pgram <- rep(0,T*K*K)
    dim(pgram)<-c(T,K,K)
    # FFT's of the series, one at a time
    xmatrix <- VectorToMatrix(X,K)
    xft<-xmatrix
    for(j in 1:K)
	xft[j,] <- fft(xmatrix[j,])
    # Compute the periodogram based on these data
    for(j in 1:K){
	for(k in 1:K){
	    pgram[,j,k]<- xft[j,]*Conj(xft[k,])/(2*pi*T)
	}
    }

    # Compute d using the log periodogram regression
    d <- rep(0,K)
    for(i in 1:K){
	# Regression
	result <- lsfit(log(1:numPoints),log(pgram[2:(numPoints+1),i,i]))
	if(verbose) ls.print(result)
	d[i] <- -result$coefficients[2]/2
    }

    return(d)
}






### Function: FIVAR.loglikelihood
## This function computes the log likelihood of the given
## FIVAR(1,d) model.
## Inputs:
##	X - KTx1 vector of observations, grouped by series
##	A1 - KxK matrix describing the AR(1)
##	Sigma - KxK covariance matrix
##	d - vector of length K of differencing parameters
##	Mean (0) - vector of length K of means
##	determinant [TRUE] - should the regression approximation to the determinant
## 		be used?  (otherwise, the naive approximation is used)
##	MinVal [-1E10] - For some extreme parameter values, the determinant
##		calculation method returns NaN or a very small number.
##		In that case, we return MinVal and print a warning.
FIVAR.loglikelihood <- function(X, A1, Sigma, d, Mean=rep(0,K), determinant=TRUE, MinVal=-1E10){
     K <- dim(A1)[1]
     T <- length(X)/K

     # Demean X
     X <- MultivariateDemean(X, Mean)

     # Compute the autocovariance sequence
     covs <- FIVAR.autoCov(A1, Sigma, d, T)

     # Compute the quadratic form
     qForm <- QuadraticFormWithPCG(covs, X)

     # Compute the determinant
     if(determinant){
     	logdet <- FIVAR.logdeterminant(n=T, Sigma=Sigma, F=A1, d=d)
     } else {
	logdet <- det(Sigma)^T
     }

     if(identical(-0.5*logdet - 0.5*qForm, NaN) || identical(-0.5*logdet - 0.5*qForm, -Inf)){
	warning("Infinite or non-numeric FIVAR loglikelihood; minimum value returned")   
	return(MinVal)
     }

     if(-0.5*logdet - 0.5*qForm<MinVal){
	warning("FIVAR log likelihood less than minimum value; actual value is ", -0.5*logdet - 0.5*qForm,
		" Minimum Value of  ", MinVal, " returned.")
	return(MinVal)
    }

     return(-0.5*logdet - 0.5*qForm)
}


### Function: FIVAR.MLE
## This function computes the maximum likelihood estimator
## for a FIVAR(1,d) model.
## Inputs:
##	X - KTx1 vector of observations, grouped by series
##	K - number of series
##	init [NULL] - initial values as a list with:
##		A1 - KxK matrix
##		Sigma - KxK matrix
##		d - K-vector of differencing parameters
##	determinant [TRUE] - should the regression approximation to the determinant
## 		be used?  (otherwise, the naive approximation is used)
##	demean [FALSE] - should the data be demeaned before estimation begins?
##	verbose [FALSE] - should the parameters be printed at each iteration?
##	boundaryDistance [0.01] - how close to the boundaries of the parameter
##		space is the search allowed to go?
##	computeHessian [FALSE] - should the Hessian be computed for the original parameterization?
## Output: list of parameters
## 3/23/08 - Now, instead of bounding parameter values by Inf or -Inf when
## this is allowed, we bound them by 1/boundaryDistance.
FIVAR.MLE <- function(X, K, init=NULL, determinant=TRUE, demean=FALSE, verbose=FALSE, boundaryDistance=0.01, computeHessian=FALSE){
    # Demean, if desired
    if(demean){
	T <- length(X)/K
	Mean <- vector(length=K)
	for(i in 1:K){
	    Mean[i] <- mean(X[((i-1)*T+1):(i*T)])
    	}
        X <- MultivariateDemean(X, Mean)
    } else {
	Mean <- rep(0,K)
    }

    # Set initial values, if necessary
    if(identical(init,NULL)){
	A1 <- 0.3*diag(K)
	Sigma <- diag(K)
	d <- rep(0.1,K)
    } else {
	if(identical(init$A1,NULL)){
	    A1 <- 0.3*diag(K)
	} else {
	    A1 <- init$A1
	}
	if(identical(init$Sigma,NULL)){
	    Sigma <- diag(K)
	} else {
	    Sigma <- init$Sigma
	}
	if(identical(init$d,NULL)){
	    d <- rep(0.1,K)
	} else {
	    d <- init$d
	}
    }

    initial <- CombineParameters(A1=A1, Sigma=Sigma, d=d, maxSVD=1-boundaryDistance)

    # Set up lower and upper bounds
    lowerSigma <- vector(length=K*(K+1)/2)

    for(i in 1:K){
	lowerSigma[(i*(i-1)/2+1):(i*(i+1)/2)] <- c(rep(-1/boundaryDistance,i-1),boundaryDistance)
    }
    lower <- c(rep(-1/boundaryDistance,K^2),
		lowerSigma,
		rep(-0.5+boundaryDistance,K))
    upper <- c(rep(Inf,K^2),
		rep(1/boundaryDistance,K*(K+1)/2),
		rep(0.5-boundaryDistance,K))


    # Create the log likelihood based on these X values
    llik <- function(params){
	# Separate out parameters
	paramlist <- ParseParameters(params, K, maxSVD=1-boundaryDistance)

	if(verbose){
	    print("Current Parameter Values and Likelihood")
	    print(paramlist)
	}
	
	loglik <- FIVAR.loglikelihood(X, A1=paramlist$A1, Sigma=paramlist$Sigma, d=paramlist$d, determinant=determinant)
	if(verbose){
	    print(loglik)
	}
	return(-1*loglik)
    }

    # Optimize with box constraints
    optim <- optim(par=initial,fn=llik,lower=lower, upper=upper, hessian=TRUE) 

    if(verbose){ 
	print(paste("Maximum Likelihood:", -1*optim$value))
	print(paste("Convergence Code (0 indicates successful convergence): ", optim$convergence))
	print(paste("Number of iterations: ", optim$iterations))
    }

    params <- ParseParameters(optim$par,K, maxSVD=1-boundaryDistance)

    # Compute the Hessian (using only the upper triangle of Sigma)
    lFunction <- function(params){
	A1 <- params[1:(K^2)]
	dim(A1) <- c(K,K)
	Sigma <- matrix(nrow=K, ncol=K)
	nextStart <- K^2+1
	for(i in 1:K){
	    Sigma[i,i:K] <- params[nextStart:(nextStart+K-i)]
	    Sigma[i:K,i] <- params[nextStart:(nextStart+K-i)]
	    nextStart <- nextStart+K-i+1
	}
	d <- params[(K^2+K*(K+1)/2+1):(K^2+K*(K+1)/2+K)]
	return(FIVAR.loglikelihood(X, A1=A1, Sigma=Sigma, d=d, determinant=determinant))
    }

    if(computeHessian){
    	hessian <- ComputeHessian(f=lFunction, params=c(params$A1, params$Sigma[!lower.tri(params$Sigma)], params$d), initialStep=boundaryDistance/2)
    } else {
	hessian <- 0
    }
 
    return(append(params, list(Hessian=hessian, originalHessian=optim$hessian, loglik=-1*optim$value, Mean=Mean)))
}




### Function: FIVAR.spectrum
## This function computes the spectrum of a FIVAR process at a given frequency.
## Inputs:
## nu - frequency for which the spectrum should be computed
## A1 - the VAR(1) term 
## Sigma - the innovation variance
## d - the differencing parameters
## Output: the value of the spectrum at that point
FIVAR.spectrum <- function(nu, A1, Sigma, d){
    K <- length(d)
    return(solve(diag((1-exp(-1i*nu))^d)) %*% 
	solve(diag(K) - A1*exp(-1i*nu)) %*%
	Sigma %*%
	t(solve(diag(K) - A1*exp(+1i*nu))) %*%
	t(solve(diag((1-exp(+1i*nu))^d))))
}

### Function: FIVAR.WhittleLikelihood
## This function takes in the cross-periodogram of the data
## and the parameters and returns the Whittle approximation 
## to the log likelihood.
## Inputs:
## 	pgram - cross-periodogram of the data (as a KxKxT array)
## 	A - the VAR(1) term 
## 	Sigma - the innovation variance
## 	d - the differencing parameters
##	MinVal [-1E10] - For some extreme parameter values, the determinant
##		calculation method returns NaN or a very small number.
##		In that case, we return MinVal and print a warning.
##	zeroCutoff [1E-12] - the largest value that we consider "0" when checking 
##		that the imaginary part is 0 (if the imaginary part exceeds zeroCutoff 
##		in absolute value, we print a warning and return MinVal)
## Output: the Whittle approximation to the log likelihood
FIVAR.WhittleLikelihood <- function(pgram, A, Sigma, d, MinVal=-1E10, zeroCutoff=1E-12){
    # Extract dimensions
    K <- dim(pgram)[3]
    T <- dim(pgram)[1]

    # Compute the sum of the trace and of the log determinants
    sum <- 0 
    for(t in 1:(T-1)){
    	# Compute the spectrum at this point
	f <- FIVAR.spectrum(2*pi*t/T, A, Sigma, d)

	# This uses Dunsmuir & Hannan's approximation
	sum <- sum - 0.5*(log(det(Sigma))+sum(diag(solve(f) %*% pgram[t,,])))
    }

    # If things have gone haywire, there will be a non-negligible imaginary part
    if(abs(Im(sum))>zeroCutoff){
	warning("Computed VARFI Whittle log likelihood is ", sum,
		" Minimum Value of  ", MinVal, " returned.")
	return(MinVal)   
    }

     if(Re(sum)<MinVal){
	warning("FIVAR Whittle log likelihood less than minimum value; actual value is ", sum,
		" Minimum Value of  ", MinVal, " returned.")
	return(MinVal)
    }

    return(Re(sum))
}


### Function: FIVAR.Whittle
## This function computes the Whittle estimator for a FIVAR(1,d) model.
## Inputs:
##	X - KTx1 vector of observations, grouped by series
##	K - number of series
##	init [NULL] - initial values as a list with:
##		A1 - KxK matrix
##		Sigma - KxK matrix
##		d - K-vector of differencing parameters
##	verbose [FALSE] - should the parameters be printed at each iteration?
##	boundaryDistance [0.01] - how close to the boundaries of the parameter
##		space is the search allowed to go?
##	demean [FALSE] - should the data be demeaned, with the means estimated by the sample mean?
##	computeHessian [FALSE] - should the Hessian be computed for the original parameterization?
## Output: list of parameters
FIVAR.Whittle <- function(X, K, init=NULL, demean=FALSE, verbose=FALSE, boundaryDistance=0.01, computeHessian=FALSE){
    # Demean, if desired (matters only in the estimation of the means)
    if(demean){
	T <- length(X)/K
	Mean <- vector(length=K)
	for(i in 1:K){
	    Mean[i] <- mean(X[((i-1)*T+1):(i*T)])
    	}
	X <- MultivariateDemean(X, Mean)
    } else {
	Mean <- rep(0,K)
    }

    T <- length(X)/K
    # Set initial values, if necessary
    if(identical(init,NULL)){
	A1 <- 0.3*diag(K)
	Sigma <- diag(K)
	d <- rep(0.1,K)
    } else {
	if(identical(init$A1,NULL)){
	    A1 <- 0.3*diag(K)
	} else {
	    A1 <- init$A1
	}
	if(identical(init$Sigma,NULL)){
	    Sigma <- diag(K)
	} else {
	    Sigma <- init$Sigma
	}
	if(identical(init$d,NULL)){
	    d <- rep(0.1,K)
	} else {
	    d <- init$d
	}
    }

    initial <- CombineParameters(A1=A1, Sigma=Sigma, d=d, maxSVD=1-boundaryDistance)

    # Set up lower and upper bounds
    lowerSigma <- vector(length=K*(K+1)/2)

    for(i in 1:K){
	lowerSigma[(i*(i-1)/2+1):(i*(i+1)/2)] <- c(rep(-1/boundaryDistance,i-1),boundaryDistance)
    }
    lower <- c(rep(-1/boundaryDistance,K^2),
		lowerSigma,
		rep(-0.5+boundaryDistance,K))
    upper <- c(rep(Inf,K^2),
		rep(1/boundaryDistance,K*(K+1)/2),
		rep(0.5-boundaryDistance,K))
		


    # Compute cross-periodogram of the series
    pgram <- rep(0,T*K*K)
    dim(pgram)<-c(T,K,K)
    # FFT's of the series, one at a time
    xmatrix <- VectorToMatrix(X,K)
    xft<-xmatrix
    for(j in 1:K)
	xft[j,] <- fft(xmatrix[j,])
    # Compute the periodogram based on these data
    for(j in 1:K){
	for(k in 1:K){
	    pgram[,j,k]<- xft[j,]*Conj(xft[k,])/(2*pi*T)
	}
    }

    # Create the log likelihood based on these X values
    llik <- function(params){
	# Separate out parameters
	paramlist <- ParseParameters(params, K, maxSVD=1-boundaryDistance)

	if(verbose){
	    print("Current Parameter Values and Likelihood")
	    print(paramlist)
	}
	
	loglik <- FIVAR.WhittleLikelihood(pgram, paramlist$A1, paramlist$Sigma, paramlist$d)
	if(verbose){
	    print(loglik)
	}
	return(-1*loglik)
    }

    # Optimize with box constraints
    optim <- optim(par=initial,fn=llik,lower=lower, upper=upper, hessian=TRUE) 
    #optim <- nlminb(start=initial,objective=llik,lower=lower, upper=upper) 

    if(verbose){
	print(paste("Maximum Likelihood:", -1*optim$value))
	print(paste("Convergence Code (0 or 1 indicates successful convergence): ", optim$convergence))
	print(paste("Number of iterations: ", optim$iterations))
    }
 
    params <- ParseParameters(optim$par,K, maxSVD=1-boundaryDistance)

    # Compute the Hessian (using only the upper triangle of Sigma)
    lFunction <- function(params){
	A1 <- params[1:(K^2)]
	dim(A1) <- c(K,K)
	Sigma <- matrix(nrow=K, ncol=K)
	nextStart <- K^2+1
	for(i in 1:K){
	    Sigma[i,i:K] <- params[nextStart:(nextStart+K-i)]
	    Sigma[i:K,i] <- params[nextStart:(nextStart+K-i)]
	    nextStart <- nextStart + K-i+1
	}
	d <- params[(K^2+K*(K+1)/2+1):(K^2+K*(K+1)/2+K)]

	return(FIVAR.WhittleLikelihood(pgram, A1, Sigma, d))
    }
    if(computeHessian){
    	hessian <- ComputeHessian(f=lFunction, params=c(params$A1, params$Sigma[!lower.tri(params$Sigma)], params$d), initialStep=boundaryDistance/2)
    } else {
	hessian <- 0
    }

    return(append(params, list(Hessian=hessian, originalHessian=optim$hessian, loglik=-1*optim$value, Mean=Mean)))
}

### FIVAR.SowellLogLikelihood
## This function computes the log likelihood of data, given the covariances and the observed values.
## Inputs:
##	X - vector of data
##	A1 - KxK matrix describing the AR(1)
##	Sigma - KxK covariance matrix
##	d - vector of length K of differencing parameters
##	determinant [TRUE] - should the determinant be included in 
##		the log likelihood?
FIVAR.SowellLogLikelihood <- function(X, A1, Sigma, d, determinant=TRUE){
     K <- dim(A1)[1]
     T <- length(X)/K

     X <- SortByTime(X, K, T)

     # Compute the autocovariance sequence
     covarray <- FIVAR.autoCov(A1, Sigma, d, T)

     # Compute the matrix decomposition
     sMatrices <- Sowell.AllMatrices(covarray)

     # Compute the quadratic form
     qForm <- Sowell.quadraticForm(X, sowellMatrices=sMatrices)

     # Compute the determinant
     if(determinant){
     	logdet <- Sowell.logdeterminant(sowellMatrices=sMatrices)
     } else {
	logdet <- det(Sigma)*T
     }

     return(-0.5*logdet - 0.5*qForm)
}

### Function: FIVAR.Sowell
## This function computes the maximum likelihood estimator
## for a FIVAR(1,d) model.
## Inputs:
##	X - KTx1 vector of observations, grouped by series
##	K - number of series
##	init [NULL] - initial values as a list with:
##		A1 - KxK matrix
##		Sigma - KxK matrix
##		d - K-vector of differencing parameters
##	determinant [FAlSE] - should the determinant be included in 
##		the log likelihood?
##	demean [FALSE] - should the data be demeaned before estimation begins?
##	verbose [FALSE] - should the parameters be printed at each iteration?
##	boundaryDistance [0.01] - how close to the boundaries of the parameter
##		space is the search allowed to go?
## Output: list of parameters
FIVAR.Sowell <- function(X, K, init=NULL,determinant=TRUE, demean=FALSE, verbose=FALSE, boundaryDistance=0.01){
    # Demean, if desired
    if(demean){
	T <- length(X)/K
	for(i in 1:K){
	    X[((i-1)*T+1):(i*T)] <- X[((i-1)*T+1):(i*T)] - mean(X[((i-1)*T+1):(i*T)])
    	}
    }

    # Set initial values, if necessary
    if(identical(init,NULL)){
	A1 <- 0.3*diag(K)
	Sigma <- diag(K)
	d <- rep(0.1,K)
    } else {
	if(identical(init$A1,NULL)){
	    A1 <- 0.3*diag(K)
	} else {
	    A1 <- init$A1
	}
	if(identical(init$Sigma,NULL)){
	    Sigma <- diag(K)
	} else {
	    Sigma <- init$Sigma
	}
	if(identical(init$d,NULL)){
	    d <- rep(0.1,K)
	} else {
	    d <- init$d
	}
    }

    initial <- CombineParameters(A1=A1, Sigma=Sigma, d=d, maxSVD=1-boundaryDistance)

    # Set up lower and upper bounds
    lowerSigma <- vector(length=K*(K+1)/2)

    for(i in 1:K){
	lowerSigma[(i*(i-1)/2+1):(i*(i+1)/2)] <- c(rep(-Inf,i-1),boundaryDistance)
    }
    lower <- c(rep(-Inf,K^2),
		lowerSigma,
		rep(-0.5+boundaryDistance,K))
    upper <- c(rep(Inf,K^2),
		rep(Inf, K*(K+1)/2),
		rep(0.5-boundaryDistance,K))
		

    # Create the log likelihood based on these X values
    llik <- function(params){
	# Separate out parameters
	paramlist <- ParseParameters(params, K, maxSVD=1-boundaryDistance)

	if(verbose){
	    print("Current Parameter Values and Likelihood")
	    print(paramlist)
	}
	
	loglik <- FIVAR.SowellLogLikelihood(X, paramlist$A1, paramlist$Sigma, paramlist$d, determinant=determinant)
	if(verbose){
	    print(loglik)
	}
	return(-1*loglik)
    }

    # Optimize with box constraints
    optim <- nlminb(initial,llik,lower=lower, upper=upper) 

    if(verbose){
	print(paste("Number of iterations:",optim$iterations))
	print(paste("Convergence Code (0 indicates successful convergence): ", optim$convergence))
	print(paste("Maximum Likelihood:", -1*optim$objective))
    }
 
    #return(ParseParameters(optim$par,K, maxSVD=1-boundaryDistance))
    return(append(ParseParameters(optim$par,K, maxSVD=1-boundaryDistance), list(loglik=-1*optim$objective)))

}



### Function: VARFI.loglikelihood
## This function computes the log likelihood of the given
## VARFI(1,d) model.
## Inputs:
##	X - KTx1 vector of observations, grouped by series
##	A1 - KxK matrix describing the AR(1)
##	Sigma - KxK covariance matrix
##	d - vector of length K of differencing parameters
##	Mean (0) - vector of length K containing the mean of each series
##	determinant [TRUE] - should the regression approximation to the determinant
## 		be used?  (otherwise, the naive approximation is used)
##	MinVal [-1E10] - For some extreme parameter values, the determinant
##		calculation method returns NaN or a very small number.
##		In that case, we return MinVal and print a warning.
VARFI.loglikelihood <- function(X, A1, Sigma, d, Mean=rep(0,K), determinant=TRUE, MinVal=-1E10){
     K <- dim(A1)[1]
     T <- length(X)/K

     # Demean X
     X <- MultivariateDemean(X, Mean)

     # Compute the autocovariance sequence
     covs <- VARFI.autoCov(A1, Sigma, d, T)

     # Compute the quadratic form
     qForm <- QuadraticFormWithPCG(covs, X)

     # Compute the determinant
     if(determinant){
     	logdet <- VARFI.logdeterminant(n=T, Sigma=Sigma, F=A1, d=d)
     } else {
	logdet <- det(Sigma)^T
     }

     if(is.nan(-0.5*logdet - 0.5*qForm)){
	message("VARFI log likelihood not a number,  Minimum Value of  ", MinVal, " returned.")
	return(MinVal)
    }

     if(-0.5*logdet - 0.5*qForm<MinVal){
	message("VARFI log likelihood less than minimum value; actual value is ", -0.5*logdet - 0.5*qForm,
		" Minimum Value of  ", MinVal, " returned.")
	return(MinVal)
    }
     return(-0.5*logdet - 0.5*qForm)
}


### Function: VARFI.MLE
## This function computes the maximum likelihood estimator
## for a VARFI(1,d) model.
## Inputs:
##	X - KTx1 vector of observations, grouped by series
##	K - number of series
##	init [NULL] - initial values as a list with:
##		A1 - KxK matrix
##		Sigma - KxK matrix
##		d - K-vector of differencing parameters
##	determinant [FAlSE] - should the determinant be included in 
##		the log likelihood?
##	verbose [FALSE] - should the parameters be printed at each iteration?
##	demean [FALSE] - should the data be demeaned before estimation begins?
##	boundaryDistance [0.01] - how close to the boundaries of the parameter
##		space is the search allowed to go?
##	computeHessian [FALSE] - should the Hessian be computed for the original parameterization?
## Output: list of parameters
VARFI.MLE <- function(X, K, init=NULL,determinant=TRUE, demean=FALSE, verbose=FALSE, boundaryDistance=0.01, computeHessian=FALSE){
    # Demean, if desired
    if(demean){
	T <- length(X)/K
	Mean <- vector(length=K)
	for(i in 1:K){
	    Mean[i] <- mean(X[((i-1)*T+1):(i*T)])
    	}
	X <- MultivariateDemean(X, Mean)
    } else {
	Mean <- rep(0,K)
    }

    # Set initial values, if necessary
    if(identical(init,NULL)){
	A1 <- 0.3*diag(K)+c(0.1,.1,rep(0,K^2-2))
	Sigma <- diag(K)
	d <- rep(0.1,K)
    } else {
	if(identical(init$A1,NULL)){
	    A1 <- 0.3*diag(K)
	} else {
	    A1 <- init$A1
	}
	if(identical(init$Sigma,NULL)){
	    Sigma <- diag(K)
	} else {
	    Sigma <- init$Sigma
	}
	if(identical(init$d,NULL)){
	    d <- rep(0.01,K)
	} else {
	    d <- init$d
	}
    }

    initial <- CombineParameters(A1=A1, Sigma=Sigma, d=d, maxSVD=1-boundaryDistance)

    # Set up lower and upper bounds
    lowerSigma <- vector(length=K*(K+1)/2)

    for(i in 1:K){
	lowerSigma[(i*(i-1)/2+1):(i*(i+1)/2)] <- c(rep(-1/boundaryDistance,i-1),boundaryDistance)
    }
    lower <- c(rep(-1/boundaryDistance,K^2),
		lowerSigma,
		rep(-0.5+boundaryDistance,K))
    upper <- c(rep(1/boundaryDistance,K^2),
		rep(1/boundaryDistance, K*(K+1)/2),
		rep(0.5-boundaryDistance,K))
		

    # Create the log likelihood based on these X values
    llik <- function(params){
	# Separate out parameters
	paramlist <- ParseParameters(params, K, maxSVD=1-boundaryDistance)

	if(verbose){
	    print("Current Parameter Values and Likelihood")
	    print(paramlist)
	}
	
	loglik <- VARFI.loglikelihood(X, paramlist$A1, paramlist$Sigma, paramlist$d, determinant=determinant)
	if(verbose){
	    print(loglik)
	}
	return(-1*loglik)
    }

    # Optimize with box constraints
    optim <- optim(par=initial,fn=llik,lower=lower, upper=upper, hessian=TRUE) 

    if(verbose){
	print(paste("Maximum Likelihood:", -1*optim$objective))
	print(paste("Convergence Code (0 or 1 indicates successful convergence): ", optim$convergence))
	print(paste("Number of iterations: ", optim$iterations))
    }

    params <- ParseParameters(optim$par,K, maxSVD=1-boundaryDistance)
 
    # Compute the Hessian (using only the upper triangle of Sigma)
    lFunction <- function(params){
	A1 <- params[1:(K^2)]
	dim(A1) <- c(K,K)
	Sigma <- matrix(nrow=K, ncol=K)
	nextStart <- K^2+1
	for(i in 1:K){
	    Sigma[i,i:K] <- params[nextStart:(nextStart+K-i)]
	    Sigma[i:K,i] <- params[nextStart:(nextStart+K-i)]
	    nextStart <- nextStart + K-i+1
	}
	d <- params[(K^2+K*(K+1)/2+1):(K^2+K*(K+1)/2+K)]
	return(VARFI.loglikelihood(X, A1, Sigma, d, determinant=determinant))
    }
    if(computeHessian){
    	hessian <- ComputeHessian(f=lFunction, params=c(params$A1, params$Sigma[!lower.tri(params$Sigma)], params$d), initialStep=boundaryDistance/2)
    } else {
	hessian <- 0
    }
 
    return(append(params, list(Hessian=hessian, originalHessian=optim$hessian, loglik=-1*optim$value, Mean=Mean)))
}



### Function: VARFI.spectrum
## This function computes the spectrum of a FIVAR
## process at a given frequency.
## Inputs:
## nu - frequency for which the spectrum should be computed
## A1 - the VAR(1) term 
## Sigma - the innovation variance
## d - the differencing parameters
## Output: the value of the spectrum at that point
VARFI.spectrum <- function(nu, A1, Sigma, d){
    K <- length(d)
    ans<- (solve(diag(K) - A1*exp(-1i*nu)) %*%
	solve(diag((1-exp(-1i*nu))^d)) %*% 
	Sigma %*%
	t(solve(diag((1-exp(+1i*nu))^d))) %*%
	t(solve(diag(K) - A1*exp(+1i*nu))))

    return(ans)
}

### Function:  VARFI.WhittleLikelihood
## This function takes in the cross-periodogram of the data
## and the parameters and returns the Whittle approximation 
## to the log likelihood.
## Inputs:
## pgram - cross-periodogram of the data (as a KxKxT array)
## A - the VAR(1) term 
## Sigma - the innovation variance
## d - the differencing parameters
## MinVal [-1E-10] - The minimum value which the likelihood will return
##	(if the computed likelihood is lower, then this will be returned, for
##	the sake of the optimizer
## zeroCutoff [1E-12] - The largest value which counts as "0" when we ensure
##	that the imaginary part of the likelihood is "0"
## Output: the Whittle approximation to the log likelihood
VARFI.WhittleLikelihood <- function(pgram, A, Sigma, d, MinVal=-1E10, zeroCutoff=1E-12){
    # Extract dimensions
    K <- dim(pgram)[3]
    T <- dim(pgram)[1]

    # Compute the sum of the trace and of the log determinants
    sum <- 0 
    for(t in 1:(T-1)){
    	# Compute the spectrum at this point
	f <- VARFI.spectrum(2*pi*t/T, A, Sigma, d)

	# This uses Dunsmuir & Hannan's approximation
	sum <- sum - 0.5*(log(det(Sigma))+sum(diag(solve(f) %*% pgram[t,,])))
    }

    # If things have gone haywire, there will be a non-negligible imaginary part
    if(abs(Im(sum))>zeroCutoff){
	warning("Computed VARFI Whittle log likelihood is ", sum,
		" Minimum Value of  ", MinVal, " returned.")
	return(MinVal)   
    }

     if(Re(sum)<MinVal){
	warning("VARFI Whittle log likelihood less than minimum value; actual value is ", sum,
		" Minimum Value of  ", MinVal, " returned.")
	return(MinVal)
   }

    return(Re(sum))
}



### Function: VARFI.Whittle
## This function computes the Whittle estimator
## for a FIVAR(1,d) model.
## Inputs:
##	X - KTx1 vector of observations, grouped by series
##	K - number of series
##	init [NULL] - initial values as a list with:
##		A1 - KxK matrix
##		Sigma - KxK matrix
##		d - K-vector of differencing parameters
##	verbose [FALSE] - should the parameters be printed at each iteration?
##	demean [FALSE] - should the data be demeaned before estimation begins?
##	boundaryDistance [0.01] - how close to the boundaries of the parameter
##		space is the search allowed to go?
##	computeHessian [FALSE] - should the Hessian be computed for the original parameterization?
## Output: list of parameters
VARFI.Whittle <- function(X, K, init=NULL, verbose=FALSE, demean=FALSE, boundaryDistance=0.01, computeHessian=FALSE){
    # Demean, if desired
    if(demean){
	T <- length(X)/K
	Mean <- vector(length=K)
	for(i in 1:K){
	    Mean[i] <- mean(X[((i-1)*T+1):(i*T)])
    	}
	X <- MultivariateDemean(X, Mean)
    } else {
	Mean <- rep(0,K)
    }


    T <- length(X)/K


    # Set initial values, if necessary
    if(identical(init,NULL)){
	A1 <- 0.3*diag(K)
	Sigma <- diag(K)
	d <- rep(0.1,K)
    } else {
	if(identical(init$A1,NULL)){
	    A1 <- 0.3*diag(K)
	} else {
	    A1 <- init$A1
	}
	if(identical(init$Sigma,NULL)){
	    Sigma <- diag(K)
	} else {
	    Sigma <- init$Sigma
	}
	if(identical(init$d,NULL)){
	    d <- rep(0.1,K)
	} else {
	    d <- init$d
	}
    }

    initial <- CombineParameters(A1=A1, Sigma=Sigma, d=d, maxSVD=1-boundaryDistance)

    # Set up lower and upper bounds
    lowerSigma <- vector(length=K*(K+1)/2)

    for(i in 1:K){
	lowerSigma[(i*(i-1)/2+1):(i*(i+1)/2)] <- c(rep(-1/boundaryDistance,i-1),boundaryDistance)
    }
    lower <- c(rep(-1/boundaryDistance,K^2),
		lowerSigma,
		rep(-0.5+boundaryDistance,K))
    upper <- c(rep(1/boundaryDistance,K^2),
		rep(1/boundaryDistance, K*(K+1)/2),
		rep(0.5-boundaryDistance,K))
		
		


    # Compute cross-periodogram of the series
    pgram <- rep(0,T*K*K)
    dim(pgram)<-c(T,K,K)
    # FFT's of the series, one at a time
    xmatrix <- VectorToMatrix(X,K)
    xft<-xmatrix
    for(j in 1:K)
	xft[j,] <- fft(xmatrix[j,])
    # Compute the periodogram based on these data
    for(j in 1:K){
	for(k in 1:K){
	    pgram[,j,k]<- xft[j,]*Conj(xft[k,])/(2*pi*T)
	}
    }

    # Create the log likelihood based on these X values
    llik <- function(params){
	# Separate out parameters
	paramlist <- ParseParameters(params, K, maxSVD=1-boundaryDistance)

	if(verbose){
	    print("Current Parameter Values and Likelihood")
	    print(paramlist)
	}
	
	loglik <- VARFI.WhittleLikelihood(pgram, paramlist$A1, paramlist$Sigma, paramlist$d)
	if(verbose){
	    print(loglik)
	}
	return(-1*loglik)
    }

    # Optimize with box constraints
    optim <- optim(par=initial,fn=llik,lower=lower, upper=upper, hessian=TRUE) 

    if(verbose){
	print(paste("Maximum Likelihood:", -1*optim$objective))
	print(paste("Convergence Code (0 or 1 indicates successful convergence): ", optim$convergence))
	print(paste("Number of iterations: ", optim$iterations))
    }
 
     params <- ParseParameters(optim$par,K, maxSVD=1-boundaryDistance)
 
    # Compute the Hessian (using only the upper triangle of Sigma)
    lFunction <- function(params){
	A1 <- params[1:(K^2)]
	dim(A1) <- c(2,2)
	Sigma <- matrix(nrow=K, ncol=K)
	nextStart <- K^2+1
	for(i in 1:K){
	    Sigma[i,i:K] <- params[nextStart:(nextStart+K-i)]
	    Sigma[i:K,i] <- params[nextStart:(nextStart+K-i)]
	    nextStart <- nextStart + K-i+1
	}
	d <- params[(K^2+K*(K+1)/2+1):(K^2+K*(K+1)/2+K)]
	return(VARFI.WhittleLikelihood(pgram, A1, Sigma, d))
    }
    if(computeHessian){
    	hessian <- ComputeHessian(f=lFunction, params=c(params$A1, params$Sigma[!lower.tri(params$Sigma)], params$d), initialStep=boundaryDistance/2)
    } else {
	hessian <- 0
    }
 
    return(append(params, list(Hessian=hessian, originalHessian=optim$hessian, loglik=-1*optim$value, Mean=Mean)))
}

### VARFI.SowellLogLikelihood
## This function computes the log likelihood of data,
## given the covariances and the observed values.
## Inputs:
##	X - vector of data
##	A1 - KxK matrix describing the AR(1)
##	Sigma - KxK covariance matrix
##	d - vector of length K of differencing parameters
##	determinant [TRUE] - should the determinant be included in 
##		the log likelihood?
VARFI.SowellLogLikelihood <- function(X, A1, Sigma, d, determinant=TRUE){
     K <- dim(A1)[1]
     T <- length(X)/K

     X <- SortByTime(X, K, T)

     # Compute the autocovariance sequence
     covarray <- VARFI.autoCov(A1, Sigma, d, T)

     # Compute the matrix decomposition
     sMatrices <- Sowell.AllMatrices(covarray)

     # Compute the quadratic form
     qForm <- Sowell.quadraticForm(X, sowellMatrices=sMatrices)

     # Compute the determinant
     if(determinant){
     	logdet <- Sowell.logdeterminant(sowellMatrices=sMatrices)
     } else {
	logdet <- det(Sigma)*T
     }

     return(-0.5*logdet - 0.5*qForm)
}

### Function: VARFI.Sowell
## This function computes the maximum likelihood estimator
## for a VARFI(1,d) model.
## Inputs:
##	X - KTx1 vector of observations, grouped by series
##	K - number of series
##	init [NULL] - initial values as a list with:
##		A1 - KxK matrix
##		Sigma - KxK matrix
##		d - K-vector of differencing parameters
##	demean [FALSE] - should the data be demeaned before estimation begins?
##	determinant [FAlSE] - should the determinant be included in 
##		the log likelihood?
##	verbose [FALSE] - should the parameters be printed at each iteration?
##	boundaryDistance [0.01] - how close to the boundaries of the parameter
##		space is the search allowed to go?
## Output: list of parameters
VARFI.Sowell <- function(X, K, init=NULL,determinant=TRUE, demean=FALSE, verbose=FALSE, boundaryDistance=0.01){
    # Demean, if desired
    if(demean){
	T <- length(X)/K
	for(i in 1:K){
	    X[((i-1)*T+1):(i*T)] <- X[((i-1)*T+1):(i*T)] - mean(X[((i-1)*T+1):(i*T)])
    	}
    }

    # Set initial values, if necessary
    if(identical(init,NULL)){
	A1 <- 0.3*diag(K)+c(0.1,.1,rep(0,K^2-2))
	Sigma <- diag(K)
	d <- rep(0.1,K)
    } else {
	if(identical(init$A1,NULL)){
	    A1 <- 0.3*diag(K)+c(0.1,.1,rep(0,K^2-2))
	} else {
	    A1 <- init$A1
	}
	if(identical(init$Sigma,NULL)){
	    Sigma <- diag(K)
	} else {
	    Sigma <- init$Sigma
	}
	if(identical(init$d,NULL)){
	    d <- rep(0.1,K)
	} else {
	    d <- init$d
	}
    }

    initial <- CombineParameters(A1=A1, Sigma=Sigma, d=d, maxSVD=1-boundaryDistance)

    # Set up lower and upper bounds
    lowerSigma <- vector(length=K*(K+1)/2)

    for(i in 1:K){
	lowerSigma[(i*(i-1)/2+1):(i*(i+1)/2)] <- c(rep(-Inf,i-1),boundaryDistance)
    }
    lower <- c(rep(-Inf,K^2),
		lowerSigma,
		rep(-0.5+boundaryDistance,K))
    upper <- c(rep(Inf,K^2),
		rep(Inf, K*(K+1)/2),
		rep(0.5-boundaryDistance,K))
		

    # Create the log likelihood based on these X values
    llik <- function(params){
	# Separate out parameters
	paramlist <- ParseParameters(params, K, maxSVD=1-boundaryDistance)

	if(verbose){
	    print("Current Parameter Values and Likelihood")
	    print(paramlist)
	}
	
	loglik <- VARFI.SowellLogLikelihood(X, paramlist$A1, paramlist$Sigma, paramlist$d, determinant=determinant)
	if(verbose){
	    print(loglik)
	}
	return(-1*loglik)
    }

    # Optimize with box constraints
    optim <- nlminb(initial,llik,lower=lower, upper=upper) 

    if(verbose){
	print(paste("Number of iterations:",optim$iterations))
	print(paste("Convergence Code (0 indicates successful convergence): ", optim$convergence))
	print(paste("Maximum Likelihood:", -1*optim$objective))
    }
 
    #return(ParseParameters(optim$par,K, maxSVD=1-boundaryDistance))
    return(append(ParseParameters(optim$par,K, maxSVD=1-boundaryDistance), list(loglik=-1*optim$objective)))

}



################ Functions for using the models #############


### Function: CrossPeriodogram
## This function takes in a vector of length KT
## and returns the cross-periodogram of the data.
## Inputs:
## X - KT vector of data OR a KxT matrix of data
## K [0] - number of time series in X if X is a vector; 0 if X is a matrix
## Output: The TxKxK cross-periodogram of the data
CrossPeriodogram <- function(X, K=0){
    if(K>0){
	T <- length(X)/K
	xmatrix <- VectorToMatrix(X,K)
    } else {
	xmatrix <- X
	T <- dim(xmatrix)[2]
	K <- dim(xmatrix)[1]
	# Warning if K > T (since in my applications, that
	# probably means that a matrix should be transposed)
	if(K > T) warning("More series than time periods in the cross-periodogram")
    }

    # Compute cross-periodogram of the series
    pgram <- rep(0,T*K*K)
    dim(pgram)<-c(T,K,K)
    # FFT's of the series, one at a time
#    xmatrix <- VectorToMatrix(X,K)
    xft<-xmatrix
    for(j in 1:K)
	xft[j,] <- fft(xmatrix[j,])
    # Compute the periodogram based on these data
    for(j in 1:K){
	for(k in 1:K){
	    pgram[,j,k]<- xft[j,]*Conj(xft[k,])/(2*pi*T)
	}
    }

    return(pgram)
}







### FIVAR.predict
## This function predicts h lags ahead for a FIVAR process,
## using the relationship:
## E(X_{T+h}|X) = Cov(X, X_{T+h}) Cov(X)^{-1}X
## where X =(X_1,...,X_T).
## Inputs:
## X - vector of length KT containing the observations
## lags [1] - a vector containing all of the number of steps ahead that the prediction should be computed
## paramlist - a list containing the parameters:
##		A1 - KxK VAR(1) matrix
##		Sigma - KxK innovation variance matrix
##		d - K-vector of differencing parameters
##		Mean - K-vector of means
## Output: a matrix of size K x length(lags), in which the predictions
##	are given in the same order as the lags were given
FIVAR.predict <- function(X, lags=1, paramlist){
    # Extract parameters
    A1 <- paramlist$A1
    Sigma <- paramlist$Sigma
    d <- paramlist$d
    Mean <- paramlist$Mean
    K <- length(d)
    T <- length(X)/K

    lags <- as.vector(lags)
    if(min(lags)<1){
	stop("Prediction is for future values (lags>0).")
    }
    maxlag <- max(lags)

    # Get the covariances needed for both the quadratic form and the prediction
    covs <- FIVAR.autoCov(A1, Sigma, d, T+maxlag)
    center <- T+maxlag

    # Demean X
    X <- MultivariateDemean(X, Mean)

    # Compute Cov(X)^{-1}X
    prod1 <- MultivariatePCG(covs[,,(center+T-1):(center-T+1)],X)

    # Compute the vector of predictions, in the same order as the list of
    # lags.
    result <- matrix(nrow=K, ncol=length(lags))
    # The matrix to hold Cov(X, X_T+h).
    firstfactor <- matrix(nrow=K, ncol=K*T)
    for(i in 1:length(lags)){
	h <- lags[i]

	# Set up Cov(X, X_T+h)
	for(k in 1:K){
	    for(t in 1:T){
		firstfactor[,(k-1)*T+t] <- covs[k,,center-T+t-h]
	    }
	}

	# Multiply
	result[,i] <- firstfactor %*% prod1
    }

    # Add the means back on
    for(k in 1:K){
	result[k,] <- result[k,] + Mean[k]
    }

    return(result)
}



### VARFI.predict
## This function predicts h lags ahead for a VARFI process,
## using the relationship:
## E(X_{T+h}|X) = Cov(X, X_{T+h}) Cov(X)^{-1}X
## where X =(X_1,...,X_T).
## Inputs:
## X - vector of length KT containing the observations
## lags [1] - a vector containing all of the number of steps ahead that the prediction should be computed
## paramlist - a list containing the parameters:
##		A1 - KxK VAR(1) matrix
##		Sigma - KxK innovation variance matrix
##		d - K-vector of differencing parameters
##		Mean - K-vector of means
## Output: a matrix of size K x length(lags), in which the predictions
##	are given in the same order as the lags were given
VARFI.predict <- function(X, lags=1, paramlist){
    # Extract parameters
    A1 <- paramlist$A1
    Sigma <- paramlist$Sigma
    d <- paramlist$d
    Mean <- paramlist$Mean
    K <- length(d)
    T <- length(X)/K

    lags <- as.vector(lags)
    if(min(lags)<1){
	stop("Prediction is for future values (lags>0).")
    }
    maxlag <- max(lags)

    # Get the covariances needed for both the quadratic form and the prediction
    covs <- VARFI.autoCov(A1, Sigma, d, T+maxlag)
    center <- T+maxlag

    # Demean X
    X <- MultivariateDemean(X, Mean)

    # Compute Cov(X)^{-1}X
    prod1 <- MultivariatePCG(covs[,,(center+T-1):(center-T+1)],X)

    # Compute the vector of predictions, in the same order as the list of
    # lags.
    result <- matrix(nrow=K, ncol=length(lags))
    # The matrix to hold Cov(X, X_T+h).
    firstfactor <- matrix(nrow=K, ncol=K*T)
    for(i in 1:length(lags)){
	h <- lags[i]

	# Set up Cov(X, X_T+h)
	for(k in 1:K){
	    for(t in 1:T){
		firstfactor[,(k-1)*T+t] <- covs[k,,center-T+t-h]
	    }
	}

	# Multiply
	result[,i] <- firstfactor %*% prod1
    }

    # Add the means back on
    for(k in 1:K){
	result[k,] <- result[k,] + Mean[k]
    }

    return(result)
}

########### Autocovariances


### Function: FIVAR.autoCov
## This function takes in the parameters of a FIVAR(1,d,0) process
## and returns the specified number of autocovariances
## Inputs:
##	F [NULL] - KxK matrix of coefficients of lagged variables (NULL implies 0)
##	Sigma - KxK symmetric and positive definite covariance matrix of the errors
##	d - vector of length K of differencing parameters
##	n - number of lags returned (from -(n-1) to n-1)
##	tol (1E-10) - maximum allowable sum of dropped VAR covariance terms
##	verbose [FALSE] - determines whether the number of lags of the VAR covariances used is printed
## Output: a K x K x n array, with the r-th autocovariance at [.,.,r-1]	
FIVAR.autoCov <- function(F=NULL,Sigma,d,n,tol=1E-10,verbose=FALSE){
    # Dealing quickly with the FIVAR(0,d,0) case
    if(identical(F,NULL)){
	return(FIVAR.autoCov0(Sigma,d,n,verbose))
    }
    if(max(abs(F))==0){
	return(FIVAR.autoCov0(Sigma,d,n,verbose))
    }

    K <- dim(F)[1]

    # Check that the d's imply stationarity
    if(max(d)>=0.5){
	stop(paste("Differencing parameters must be less than 0.5.  Current differencing parameters: ", d))
    }

    # Figure out M for the finite-sample approximation and check for stationarity
    # This is now based on using the largest singular value (the Euclidean norm)
    normF <- svd(F,nu=0,nv=0)$d[1]
    if(abs(normF)>=1){
	stop("VAR(1) matrix has singular values greater than one; autocovariances cannot be calculated.")
    }

    cov0 <- VAR.cov0(F,Sigma)
    norm0 <-  svd(cov0,nu=0,nv=0)$d[1]
    maxlag <- ceiling((log(1-normF)+log(tol)-log(norm0))/log(normF) + 1)
    if(verbose) print(c("Lags of VAR needed:",maxlag))

    # Get the VAR covariances
    VARcov <- array(dim=c(K,K,2*maxlag+1))
    # Store autocovariances 0 - maxlag
    VARcov[,,(maxlag+1):(2*maxlag+1)] <- VAR.autoCov(F,Sigma,maxlag,known0=cov0)
    # Transpose to get autocovariances at negative lags
    for(i in 1:maxlag)
	VARcov[,,i] <- t(VARcov[,,2*maxlag+2-i])

    # Get the ARFIMA covariances
    # We need n+max(maxlag, n+1) covariances of the ARFIMA process
    lagsARFIMA <- n+max(maxlag, n+1)
    ARFIMAcov <- array(dim=c(K,K,2*lagsARFIMA+1))
    for(i in 1:K){
	for(j in 1:K){
	    ARFIMAcov[i,j,] <- ARFIMA.crossCov(d[i],d[j],(-lagsARFIMA):lagsARFIMA)
	}
    }

    # Set up storage for the result
    result <- array(dim=c(K,K,2*(n-1)+1))

    for(i in 1:K){
	for(j in 1:K){
	   conv <- convolution(VARcov[i,j,],ARFIMAcov[i,j,])
    	   # Set up the range to extract from each convolution
	   # We first find the middle of the returned array (which has
	   # length equal to one less than the sum of the two input 
	   # arrays).
	   middle <- (length(conv)+1)/2
	   # Then, we extract (n-1) on either side of the middle.
	   result[i,j,] <- conv[(middle-(n-1)):(middle+(n-1))]
	}
     }

    return(result)
}


### Function: Sowell.autoCov
## This function takes in the parameters of a VARFIMA process
## and returns the specified number of autocovariances
## Inputs:
##	F - KxK matrix of coefficients of lagged variables
##	Sigma - KxK symmetric and positive definite covariance matrix of the errors
##	d - vector of length K of differencing parameters
##	n - number of lags returned (from -(n-1) to n-1)
## Output: a K x K x n array, with the r-th autocovariance at [.,.,r-1]
Sowell.autoCov <- function(F, Sigma, d, n, V=NULL){
    K <- dim(F)[1]
    if(identical(V, NULL)) V <- diag(K)

    # Compute the Wold decomposition of this process (Sowell, page 11)
    # A(L) is a scalar polynomial of degree K with roots equal to the 
    # eigenvalues of F; B(L) multiplies the eigenvectors by the matrix
    # of polynomials 
    M <- K-1
    H <- K
    # rho is the vector containing all the eigenvalues of F
    eigF <- eigen(F)
    rho <- eigF$values
    P <- eigF$vectors

    # If there are repeated eigenvalues, this method fails
    if(sum(duplicated(rho))>0)
	stop("The eigenvalues of F are not all unique; Sowell cannot be used")

    ### B[i,j,n] is B_{i,j}(n) from page 12 
    # First, compute Lambda, the diagonal matrix of polynomials
    # with Lambda(k) equal to the product of (1-rho_jL), where 
    # we omit j=k
    Lambda <- array(data=0,dim=c(K,K,M+1))
    for(i in 1:K){
	Lambda[i,i,] <- c(1,rep(0,M))
	for(j in 1:K){
	    if(i != j){
		Lambda[i,i,] <- Lambda[i,i,]-c(0, rho[j] * Lambda[i,i,1:M])
	    }
	}
    }
    # B is P*Lambda*P^{-1} 
    B <- array(data=0,dim=c(K,K,M+1))
    for(m in 1:(M+1)){
	B[,,m] <- (P) %*% Lambda[,,m] %*% solve(P)
    }

    # psi(i,j), based on the equation on page 12
    psi <- array(dim=c(K,K,2*M+1))
    for(i in 1:K){
	for(j in 1:K){
	    for(l in -M:M){
	    	sum <- 0
	    	    for(h in 1:K){
		    	for(t in 1:K){
			    # Sowell's summation endpoints are wrong
		    	    #for(s in max(0,l):min(M,M-l)){
		    	    for(s in max(0,l):min(M,M+l)){
			    	sum <- sum+Sigma[h,t]*B[i,h,s+1]*B[j,t,s-l+1]
		    	    }
		        }
	            }
	    	psi[i,j,l+M+1] <- sum
	    }
	}
    }

    # zeta_j (page 11)
    zeta <- vector(length=H)
    for(j in 1:H){
	zeta[j] <- 1/rho[j]
	for(i in 1:H){
	    zeta[j] <- zeta[j]/(1-rho[i]*rho[j])
	}
	for(m in 1:H){
	    if(m != j){
		zeta[j] <- zeta[j]/(rho[j]-rho[m])
	    }
	}
    }

    # Compute matrix of Sowell's C(di,dj, h, rho)
    C <- SowellC(d,rho,H+M+n,H)

    # Compute all the sums (equation 7, page 13)
    autoCov <- array(dim=c(K,K,2*(n-1)+1))
    for(i in 1:K){
	for(j in 1:K){
	    for(s in (-(n-1)):(n-1)){
		autoCov[i,j,s+n] <- 0
		for(l in (-1*M):M){
		    for(m in 1:H){
			for(nn in 1:K){
			    for(r in 1:K){
				autoCov[i,j,s+n] <- autoCov[i,j,s+n] + V[i,nn]*V[j,r]*psi[nn,r,(M+1)+l]*zeta[m]*C[nn,r,m,(H+M+n)+H+l-s]
				if(V[i,nn]*V[j,r]*psi[nn,r,l+M+1]*zeta[m]*C[i,j,m,(H+M+n)+H+l-s] !=0){
				}
			    }
			}
		    }
		}
	    }
	}
    }

    return(autoCov)

}

### Function: FIVAR.integralCov
## This function computes the autocovariance matrices of a FIVAR(1,d,0)
## process using numerical integration methods.
## Inputs:
##	F - KxK matrix of coefficients of lagged variables
##	Sigma - KxK symmetric and positive definite covariance matrix of the errors
##	d - vector of length K of differencing parameters
##	n - number of lags returned (from -(n-1) to n-1)
## Output: a K x K x n array, with the r-th autocovariance at [.,.,r-1]
FIVAR.integralCov <- function(F,Sigma,d,n){
    K <- dim(F)[1]

    # Set up storage for results
    autoCov <- array(dim=c(K,K,2*(n-1)+1))

    for(i in 1:K){
	for(j in 1:K){
	    print(paste("Current pair:", i,j))
	    # This function defines the [i,j] component of the spectral density;
	    # The argument MUST be scalar
   	    oneSpectrum <- function(lambda){
		# First compute the whole matrix product; then, we extract [i,j] element
		Dinv <- diag((1-exp(-1i*lambda))^-d)
		A <- diag(2) - F*exp(-1i*lambda)		
		Ainv <- solve(A)
		
		product <- (1/(2*pi)) * Dinv %*% Ainv %*% Sigma %*%
				t(Conj(Ainv)) %*% Conj(Dinv)
		answer <- product[i,j]	

		# If the sum is Inf or NaN, return 0 (this happens only on a 
		# set of measure 0
		if(identical(Re(answer),Inf) || identical(Re(answer),NaN))  answer <- 0
    		
		return(answer)
	    }
	    # This function allows oneSpectrum to be applied to a list
	    spectrum <- function(lambda) {
		sapply(lambda, oneSpectrum)
	    }
	    for(r in (-n+1):(n-1)){
	    	# Set up a function to define the integral
		integrand <- function(lambda){
		    Re(exp(1i*r*lambda)*spectrum(lambda))
		}

	        result <- integrate(integrand, lower=-pi,upper=pi)


	        autoCov[i,j,r+n] <- result$value
	    }
	}
    }

    return(autoCov)
} 
	   


### Function: VARFI.autoCov
## This function computes the autocovariances of a VARFI(1,d) process.
## Inputs:
##	F [NULL] - KxK matrix of coefficients of lagged variables (NULL implies 0)
##	Sigma - KxK symmetric and positive definite covariance matrix of the errors
##	d - vector of length K of differencing parameters
##	n - number of lags returned (from -(n-1) to n-1)
##	tol (1E-10) - maximum allowable sum of dropped VAR covariance terms
## 	verbose [FALSE] - determines whether the number of lags of the VAR covariances used is printed
## Output: a K x K x 2n-1 array, with the r-th autocovariance at [.,.,r-1]
VARFI.autoCov <- function(F=NULL,Sigma,d,n,tol=1E-10,verbose=FALSE){
    # Dealing quickly with the VARFIMA(0,d,0) case
    if(identical(F,NULL)){
	return(VARFI.autoCov0(Sigma,d,n,verbose))
    }
    if(max(abs(F))==0){
	return(VARFI.autoCov0(Sigma,d,n,verbose))
    }

    K <- dim(F)[1]

    # Eigenvalue decomposition of F
    eigF <- eigen(F)

    ### Compute maximum lag needed over all (i,j) pairs
    G <- max(VARFI.autoCov0(Sigma=Sigma, d=d, n=1)[,,1])
    maxLambda = max(abs(eigF$values))
    maxlag <- ceiling((log(tol)+2*log(1-maxLambda)+log(1+maxLambda))/log(maxLambda))
    if(verbose) print(paste("Number of lags used in the approximation: ",maxlag))

    # Compute L's up to this lag
    L <- array(dim=c(K,K,2*maxlag-1))
    for(i in 1:K){
	for(j in 1:K){
	    L[i,j,maxlag:(2*maxlag-1)] <- eigF$values[i]^(0:(maxlag-1))/(1-eigF$values[i]*Conj(eigF$values[j]))
	    L[i,j,(maxlag-1):1] <- Conj(eigF$values[j])^(1:(maxlag-1))/(1-eigF$values[i]*Conj(eigF$values[j]))
  	}
    }
  
    ### Compute H
    # The integral of the inner matrix of Q is just the lag 0 autocovariance when F=0
    # We need n+max(maxlag, n+1) covariances of the ARFIMA process
    lagsARFIMA <- n+max(maxlag, n+1)
    Hinner <- VARFI.autoCov0(Sigma=Sigma, d=d, n=n+max(maxlag,n+1))
    # Multiply by V's on either side at each lag
    Vleft <- solve(eigF$vectors)
    H <- Hinner
    for(r in 1:dim(Hinner)[3]){
    	H[,,r] <- Vleft %*% Hinner[,,r] %*% Conj(t(Vleft))
    }

    # Convolutions
    innerResult <- array(dim=c(K,K,2*n-1))
    for(i in 1:K){ 
	for(j in 1:K){
	    conv <- convolution(L[i,j,],H[i,j,])
	    # We first find the middle of the returned array (which has
	    # length equal to one less than the sum of the two input 
	    # arrays).
	    middle <- (length(conv)+1)/2
	    # Then, we extract (n-1) on either side of the middle.
	    innerResult[i,j,] <- conv[(middle-(n-1)):(middle+(n-1))]
  	}
    }

    # Multiply by the outer V's
    result <- innerResult
    for(r in 1:dim(innerResult)[3]){
	result[,,r] <- eigF$vectors %*% innerResult[,,r] %*% Conj(t(eigF$vectors))
    }

    return(Re(result))
}


### Function: VARFI.integralCov
## This function computes the autocovariance matrices of a VARFI(1,d,0)
## process using numerical integration methods.
## Inputs:
##	F - KxK matrix of coefficients of lagged variables
##	Sigma - KxK symmetric and positive definite covariance matrix of the errors
##	d - vector of length K of differencing parameters
##	n - number of lags returned (from -(n-1) to n-1)
## Output: a K x K x n array, with the r-th autocovariance at [.,.,r-1]
VARFI.integralCov <- function(F,Sigma,d,n){
    K <- dim(F)[1]

    # Set up storage for results
    autoCov <- array(dim=c(K,K,2*(n-1)+1))

    for(i in 1:K){
	for(j in 1:K){
	    # This function defines the [i,j] component of the spectral density;
	    # The argument MUST be scalar
   	    oneSpectrum <- function(lambda){
		# First compute the whole matrix product; then, we extract [i,j] element
		Dinv <- diag((1-exp(-1i*lambda))^-d)
		A <- diag(2) - F*exp(-1i*lambda)		
		Ainv <- solve(A)
		
		product <- (1/(2*pi)) * Ainv %*% Dinv %*% Sigma %*%
				Conj(Dinv) %*% t(Conj(Ainv)) 
		answer <- product[i,j]	

		# If the sum is Inf or NaN, return 0 (this happens only on a 
		# set of measure 0
		if(identical(Re(answer),Inf) || identical(Re(answer),NaN))  answer <- 0

		return(answer)
	    }
	    # This function allows oneSpectrum to be applied to a list
	    spectrum <- function(lambda) {
		sapply(lambda, oneSpectrum)
	    }
	    for(r in (-n+1):(n-1)){
	    	# Set up a function to define the integral
		integrand <- function(lambda){
		    Re(exp(1i*r*lambda)*spectrum(lambda))
		}

	        result <- integrate(integrand, lower=-pi,upper=pi)

	        autoCov[i,j,r+n] <- result$value
	    }
	}
    }

    return(autoCov)
} 





### Function: Cointegrate
## This function takes in an array of covariances 
## for an uncointegrated model and a cointegrating matrix 
## and returns the covariances of the cointegrated model.
## Inputs:
##	oldCovs - array of covariances of the uncointegrated model
##	coMatrix - cointegrating matrix
## Output: newCovs - array of covariances of the cointegrated model
Cointegrate <- function(oldCovs, coMatrix){
    T <- dim(oldCovs)[3]

    newCovs <- array(dim=dim(oldCovs))
    invCoMatrix <- solve(coMatrix)
    for(t in 1:T){
	newCovs[,,t] <- invCoMatrix %*% oldCovs[,,t] %*% t(invCoMatrix)
    }

    return(newCovs)
}


########## Simulation


### Function: FIVAR.simulationSetup
## This function sets up everything necessary to simulate from a 
## FIVAR process.
## Input:
##	T - number of periods to simulate
##	d - vector of differencing parameters
##	Sigma - innovation variance matrix
##	F [NULL] - VAR(1) parameter
##	verbose [FALSE] - should intermediate messages be printed?
##	retry [2] - how many times should the function get additional
##		covariances to make the circulant embedding positive
##		definite if it fails the first time
## Output: the array of blocks needed for FFT simulation, or FALSE
##	if a positive definite circulant embedding was not found
FIVAR.simulationSetup <- function(T, d, Sigma, F=NULL, verbose=FALSE, retry=2){
    # Pick a number of covariances, M, such that 2M-1
    # is easy to factor (a power of 3) and M>=T
    M <- (3^ceiling(log(2*T-1,base=3))+1)/2

    while(retry >=0){
    	# Get array of covariances
    	if(identical(F,NULL)){
	    covarray <- FIVAR.autoCov0(Sigma, d, M, verbose=verbose)
    	} else {
	    covarray <- FIVAR.autoCov(F,Sigma, d, M, verbose=verbose)
    	}

        blocks <- SetUpBlockSqrt(covarray, verbose=verbose)

  	if(!identical(blocks,FALSE)) break

	# If the blocks fail (because the circulant embedding matrix
	# was not positive definite), decrement the retry counter
	# and increase M to try again
	retry <- retry-1
	M <- 3*M-1
    }

    if(retry==-1) return(FALSE)
    if(verbose) print(paste("Number of lags used in simulation = ", M))

    return(blocks)
}

### Function: VARFI.simulationSetup
## This function sets up everything necessary to simulate from a 
## VARFI process.
## Input:
##	T - number of periods to simulate
##	d - vector of differencing parameters
##	Sigma - innovation variance matrix
##	F [NULL] - VAR(1) parameter
##	verbose [FALSE] - should intermediate messages be printed?
##	retry [2] - how many times should the function get additional
##		covariances to make the circulant embedding positive
##		definite if it fails the first time
## Output: the array of blocks needed for FFT simulation, or FALSE
##	if a positive definite circulant embedding was not found
VARFI.simulationSetup <- function(T, d, Sigma, F=NULL, verbose=FALSE, retry=2){
    # Pick a number of covariances, M, such that 2M-1
    # is easy to factor (a power of 3) and M>=T
    M <- (3^ceiling(log(2*T-1,base=3))+1)/2

    while(retry >=0){
    	# Get array of covariances
    	if(identical(F,NULL)){
	    covarray <- VARFI.autoCov0(Sigma, d, M, verbose=verbose)
    	} else {
	    covarray <- VARFI.autoCov(F,Sigma, d, M, verbose=verbose)
    	}

        blocks <- SetUpBlockSqrt(covarray, verbose=verbose)

  	if(!identical(blocks,FALSE)) break

	# If the blocks fail (because the circulant embedding matrix
	# was not positive definite), decrement the retry counter
	# and increase M to try again
	retry <- retry-1
	M <- 3*M-1
    }

    if(retry==-1) return(FALSE)
    if(verbose) print(paste("Number of lags used in simulation = ", M))

    return(blocks)
}


### Function: FFTsimulate
## This function simulates T periods of a k-variate process, given the 
## block matrix defining the block circulant embedding of the covariance
## matrix. 
## Input:
##	T - number of periods to simulate
##	SqrtBlocks - a KxKxM array containing the square roots of the 
##		blocks of the transformed block circulant embedding (M>2T)
## Output: a matrix of size KxT containing the simulated values
FFTsimulate <- function(T, SqrtBlocks){
    # Extract dimensions
    K <- dim(SqrtBlocks)[1]
    M <- dim(SqrtBlocks)[3]

    # Generate Y = (F x I)*U
    Y <- vector(length=K*M)
    for(i in 1:K){
	Y[((i-1)*M+1):(i*M)] <- SimulateFourierNormal(M)
    }

    # Multiply by P
    Y2 <- SortByTime(Y, K, M)


    # Multiply by B^(1/2)
    Y3 <- Y2
    for(i in 1:M){
	Y3[((i-1)*K + 1):(i*K)] <- SqrtBlocks[,,i] %*% Y2[((i-1)*K + 1):(i*K)]
    }

    # Multiply by P*
    Y4 <- SortBySeries(Y3, K, M)

    # Multiply by (F x I)
    Y5 <- Y4
    for(i in 1:K){
	Y5[((i-1)*M+1):(i*M)] <- fft(Y4[((i-1)*M+1):(i*M)], inverse=FALSE)
    }

    # Generate a subset of size T
    X <- vector(length=K*T)
    for(i in 1:K){
	X[((i-1)*T+1):(i*T)] <- Y5[((i-1)*M+1):((i-1)*M+T)]
    }
    
    return(Re(X))
}    


### Function: SetUpBlockSqrt
## This function takes in a block Toeplitz covariance
## matrix (in the form of a KxKx(2T-1) array) and
## returns a KxKx(2T) array in which the t-th level is the 
## square root matrix needed for simulation.
## It checks that the circulant embedding matrix implicit
## in the computation is positive definite and returns FALSE
## if it is not.
## Input: covarray - KxKx(2T-1) array of covariances
## Output: cholarray - KxKx(2T) array of cholesky decompositions
SetUpBlockSqrt <- function(covarray, verbose=FALSE){
    # Extract dimenstions
    T <- (dim(covarray)[3]+1)/2
    K <- dim(covarray)[1]

    # The circulant embedding
    circembed <- SmallBlockCircEmbed(covarray)

    # We pass through each block and
    # compute the eigenvalues of the circulant embedding
    Aeigen<- rep(0, K*K*(2*T-1))
    dim(Aeigen) <- c(K,K,2*T-1)
   
    for(i in 1:K){
	for(j in 1:K){
	    Aeigen[i,j,] <- fft(circembed[i,j,],inverse=TRUE)
	}
    }

    # Compute square roots
    sqrtarray <- rep(0, K*K*(2*T-1))
    dim(sqrtarray) <- c(K,K,2*T-1)

    for(t in 1:(2*T-1)){    
	eig <- eigen(Aeigen[,,t])
	# Check positivity constraint
	if(min(Re(eig$values))<=0){
	    if(verbose){ 
		print(paste("Positivity constraint failed at level ", t))
		print(Aeigen[,,t])
		print(eig)
	    }
	    return(FALSE)
	}
	# Print a warning if the eigenvalues were not real (THIS SHOULD NOT HAPPEN)
	if(sum(abs(eig$values-Re(eig$values)))>1E-8){
	    warning("Eigenvalues had non-negligible imaginary parts in SetUpBlockSqrt matrix at level ", t)
	}
	sqrtarray[,,t] <- eig$vectors %*% diag(sqrt(eig$values) ) %*% solve(eig$vectors)
    }

    return(sqrtarray)
}


############### Functions used by the code ###################

### Function: ReconstructFromCholesky
## This function takes in the Cholesky decomposition
## of a matrix, written as a vector by column, and
## reconstructrs the original matrix.
## Input: cM - vector of length K(K+1)/2
##	 K [0] - K (so that we don't need to figure it out)
## Output: M = cM %*% cM
ReconstructFromCholesky <- function(cM,K=0){
    if(K==0){
    	K <- floor(sqrt(2*length(cM)))
    }

    cholMatrix <- matrix(data=0,nrow=K,ncol=K)

    for(i in 1:K){
	cholMatrix[1:i,i] <- cM[(i*(i-1)/2+1):(i*(i+1)/2)]
    }

    return(t(cholMatrix) %*% cholMatrix)
}

### Function: ConvertToCholesky
## This function takes a (positive definite) matrix and returns the Cholesky decomposition
## of a matrix, written as a vector by row.
## Input: M - a square, positive definite matrix
ConvertToCholesky <- function(M){
    K <- dim(M)[2]

    cholMatrix <- chol(M)

    cM <- vector(length=K*(K+1)/2)

    for(i in 1:K){
	cM[(i*(i-1)/2+1):(i*(i+1)/2)] <- cholMatrix[1:i,i]
    }

    return(cM)
}

### Function: ConvertToAK
## This function takes in a matrix that necessarily has singular
## values less than 1 and returns a representation that
## can range over all possible matrices.
## From Ansley-Kohn (1986) lemma 2.2:
## B^-1 * B^(-1)' = I-PP'
## A = BP
## 1/27/08 (RS): To ensure that the maximum singular value doesn't
## get too close to one, we now rescale.  If the maximum singular
## value of P1 is greater than maxSVD, this will fail.
ConvertToAK <- function(P1, maxSVD=0.99){
    K <- dim(P1)[1]
    scaledP1 <- P1/maxSVD
    B <- solve(chol(diag(K)-scaledP1 %*% t(scaledP1)))
    A <- t(B) %*% scaledP1
    return(A)
}

### Function: ConstructFromAK
## This function takes the representation of a matrix that
## ensures that it will have SVD less than 1 and converts
## it to a matrix.
## From Ansley-Kohn (1986) lemma 2.2:
## BB'=I+AA'
## P=B^-1 A
## 1/27/08 (RS): To ensure that the maximum singular value doesn't
## get too close to one, we now rescale.
ConstructFromAK <- function(A1, maxSVD=0.99){
    K <- dim(A1)[1]
    B <- chol(diag(K)+A1 %*% t(A1))
    scaledP <- t(solve(B)) %*% A1
    P <- scaledP * maxSVD
    return(P)
}

### Function: ParseParameters
## This function takes in a vector containing the parameters 
## (A1's A-K representation, Sigma's Cholesky decomposition, d) and returns a 
## list of the parameters separately.
ParseParameters <- function(paramVec,K, maxSVD=0.99){
    result <- list(A1=NULL, Sigma=NULL, d=NULL)

    result$Sigma <- ReconstructFromCholesky(paramVec[(K^2+1):(K^2+K*(K+1)/2)])

    M <- matrix(data=paramVec[1:K^2],nrow=K,ncol=K)
    result$A1 <- ConstructFromAK(M, maxSVD=maxSVD)
    result$d <- paramVec[(1.5*K^2+K/2+1):(1.5*(K^2+K))]

    return(result)
}

### Function: CombineParameters
## This function takes in the parameters and returns a vector
## containing representations of each (the Ansley-Kohn representation
## for A1 and the Cholesky decomposition of Sigma).
CombineParameters <- function(A1, Sigma, d, maxSVD=0.99){
    K <- length(d)

    result <- vector(length=1.5*(K^2+K))

    result[1:K^2] <- ConvertToAK(A1, maxSVD=maxSVD)
    result[(K^2+1):(K^2+K*(K+1)/2)] <- ConvertToCholesky(Sigma)
    result[(1.5*K^2+K/2+1):(1.5*(K^2+K))] <- d

    return(result)
}

### Function: VectorToMatrix
## This function takes in a vector of T observations for K 
## series and returns a KxT matrix where each column
## is one series.
## Inputs:
## Xvector - vector of length KT, sorted by series
## K - number of series
## Output: KxT matrix
VectorToMatrix <- function(Xvector, K){
    T <- length(Xvector)/K

    Xmatrix <- matrix(nrow=K,ncol=T)

    for(i in 1:K){
	Xmatrix[i,] <- Xvector[((i-1)*T+1):(i*T)]
    }

    return(Xmatrix)
}


### Function: MultivariateDemean
## This function takes in a vector containing K time series and 
## a K-vector of means and returns the demeaned versions of the time
## series.
## Inputs:
##	X - KT vector of data
##	means - K vector of means
MultivariateDemean <- function(X, means){
    K <- length(means)
    T <- length(X)/K

    result <- vector(length=K*T)
    for(i in 1:K){
	result[((i-1)*T+1):(i*T)] <- X[((i-1)*T+1):(i*T)] - means[i]
    }

    return(result)
}

### Function: ExtractSubset
## This function takes in a vector of length KT containing
## K time series and returns the first S observations of
## each time series.
## Inputs:
##	fullData - vector of length KT
##	K - number of time series
##	S - the length of the subset desired
ExtractSubset <- function(fullData, K, S){
    T <- length(fullData)/K

    result <- vector(length=K*S)

    for(i in 1:K){
	result[((i-1)*S+1):(i*S)] <- fullData[((i-1)*T+1):((i-1)*T+S)]
    }

    return(result)
}

### Function: ComputeHessian
## This function takes in a vector of function values and a 
## function and computes the Hessian at that set of values.
## Inputs: 
##	f - Function for which the Hessian is to be evaluated
##	params - Point at which the Hessian is to be evaluated
##	initialStep [1E-2] - Initial stepsize used for evaluation
##	minStep [1E-5] - Minimum step size allowed for evaluation
##	relTol [1E-4] - Percent difference between computed Hessians at consecutive
##		step sizes allowed to be considered convergent
## 7/2/08 (RS): We now choose the step size dynamically for each pair, 
## reducing it until the estimates stabilize.
ComputeHessian <- function(f, params, initialStep=1E-2, relTol=1E-4, minStep=1E-5){
    n <- length(params)

    result <- matrix(nrow=n, ncol=n)

    # Compute the value at this point
    f0 <- f(params)

    # Compute the diagonal elements
    for(i in 1:n){
	stepsize <- initialStep
	nextstep <- rep(0,n)
	nextstep[i] <- stepsize
	result1 <- (f(params + 2*nextstep)-2*f(params + nextstep)+f0)/stepsize^2
	while(stepsize > minStep){
	    stepsize <- stepsize/2
	    nextstep <- rep(0,n)
	    nextstep[i] <- stepsize
	    result2 <- (f(params + 2*nextstep)-2*f(params + nextstep)+f0)/stepsize^2
	    if(abs(result1-result2)/abs(result1)<relTol) break
	    result1 <- result2
	}
	result[i,i] <- result2
    }

    # Compute the off-diagonal elements
    for(i in 2:n){
	for(j in 1:(i-1)){
	    stepsize <- initialStep
	    nextstep1 <- rep(0,n)
	    nextstep1[i] <- stepsize
	    nextstep2 <- rep(0,n)
	    nextstep2[j] <- stepsize
	    result1 <- (f(params+nextstep1+nextstep2)-
		f(params+nextstep1) - f(params+nextstep2) + f0)/stepsize^2
	    while(stepsize > minStep){
		stepsize <- stepsize/2
	    	nextstep1 <- rep(0,n)
	    	nextstep1[i] <- stepsize
	    	nextstep2 <- rep(0,n)
	    	nextstep2[j] <- stepsize
	    	result2 <- (f(params+nextstep1+nextstep2)-
			f(params+nextstep1) - f(params+nextstep2) + f0)/stepsize^2
		if(abs(result1-result2)/abs(result1)<relTol) break
		result1 <- result2
	    }
	    result[i,j] <- result2
	    result[j,i] <- result[i,j]
	}
    }
	    
    return(result)

}




### Function: HyperSum
## This function computes the sum of terms of the form:
## 	Gamma(n+a)Gamma(n+b)/Gamma(n+c)n! z^n
## times Gamma(c)/Gamma(a)Gamma(b), from n=0 to infinity.
## (To do this, we continue adding terms until z^n/(1-z) is less than the 
## convergence criterion -- this is an APPROXIMATION to the true error.)
## Inputs: a - first numerator parameter
##	  b - second numerator parameter
##	  c - denominator parameter
##	  z - function argument (being raised to n^th powers)
##	  maxError [1E-10] - maximum value of z^n/(1-z) that allows convergence
##	  maxSummands [1000] - maximum number of summands allowed (after this,
##		the sum will be returned with a warning
# Output: sum
HyperSum <- function(a,b,c,z, maxError=1E-10, maxSummands=1000){
    sum <- 0
    n <- 0
    error <- 1/(1-z)

    ConstantTerm <- exp(lgamma(c)-lgamma(a)-lgamma(b))

    # If any of a, b, or c is 0, compute the constant term directly
    if(identical(ConstantTerm,NaN)){
	ConstantTerm <- gamma(c)/(gamma(a)*gamma(b))
    }

    while(n<=maxSummands & error>maxError){
	logNextTerm <- lgamma(n+a)+lgamma(n+b)-lgamma(n+c)-lgamma(n+1)
	NextTerm <- ConstantTerm*exp(logNextTerm)*z^n
	sum <- sum+NextTerm
	error <- z*error
	n <- n+1
    } 
    if(n==maxSummands) warning(c("HyperSum: Maximum number of summands reached with error ", error))
    return(sum)
}


### Function: SowellC
## This function computes all of Sowell's C(di,dj, h, rho)
## for h between -(MaxH-1) and +(MaxH-1), given a vector of 
## d's and rho's.  
## Inputs:
##	d - vector of differencing parameters (length K)
##	rho - vector of length M
##	MaxH - maximum lag (h) needed to calculate
##	H - H (needed to compute rho^(2*H))
## Output: A KxKxMx(2*(MaxH-1)+1) array
SowellC <- function(d,rho, MaxH,H){
    K <- length(d)
    M <- length(rho)

    # Store the hypergeometric functions using the recursion suggested on page 26 (if d[i] and d[j] are not 0):
    # We store F(d[i]+h, 1, 1-d[j]+h,rho[m]); the other term is found by switching
    # i and j and reversing the sign of h
    if(max(abs(d))!=0){
    	allF <- array(dim=c(K,K,M,2*MaxH-1))
	allFcheck <- allF
    	for(i in 1:K){
	    for(j in 1:K){
	    	for(m in 1:M){
		    allF[i,j,m,MaxH] <- HyperSum(d[j],1,1-d[i],rho[m])

		    # If rho[m]==0, recursive methods don't work
		    if(rho[m]==0){
		     	for(h in 1:(MaxH-1)){
		    	    # Forward
			    allF[i,j,m,MaxH+h] <- HyperSum(d[j]+h, 1, 1-d[i]+h,rho[m])
		    	    # Backward
			    allF[i,j,m,MaxH-h] <- HyperSum(d[j]-h, 1, 1-d[i]-h,rho[m])
			}
		    } else {		
		    # Compute the other h terms recursively
		     	for(h in 1:(MaxH-1)){
		    	    # Forward
			    allF[i,j,m,MaxH+h] <- HyperSum(d[j]+h, 1, 1-d[i]+h,rho[m])
		    	    # Backward
			    allF[i,j,m,MaxH-h] <- HyperSum(d[j]-h, 1, 1-d[i]-h,rho[m])
			}

		    }
	    	}
	    }
    	}
    }

    # Storage for output
    C <- array(dim=c(K,K,M,2*MaxH-1))
    for(i in 1:K){
	for(j in 1:K){
	    for(m in 1:M){
		for(h in (-(MaxH-1)):(MaxH-1)){
		    if(d[i]!=0 && d[j] !=0){
			C[i,j,m,h+MaxH] <- gamma(1-d[i]-d[j])*gamma(d[j]+h)/(gamma(1-d[i]+h)*gamma(1-d[j])*
				gamma(d[j]))*(rho[m]^(2*H)*allF[i,j,m,MaxH+h] +	allF[j,i,m,MaxH-h]-1)
		    }
		    if(d[i]==0 && d[j]==0){
			if(h <=0){
			    C[i,j,m,h+MaxH] <- rho[m]^(2*K-h)
			}
			if(h>0){
			    C[i,j,m,h+MaxH] <- rho[m]^h
			}
		    }
		    if(d[i]==0 && d[j] !=0){
			if(h <=0){
			    C[i,j,m,h+MaxH] <- rho[m]^(2*K-h)*HyperSum(d[j],1,1,rho[m])
			}
			if(h>0){
			    sum <- 0
			    for(q in 1:h){
				sum <- sum + rho[m]^q*gamma(d[j]+h-q)*gamma(1+h)/gamma(d[j]+h)/gamma(1+h-q)
			    }
			    C[i,j,m,h+MaxH] <- gamma(d[j]+h)/gamma(d[j])/gamma(2+h)*
				(rho[m]^(2*K)*HyperSum(d[j]+h,1,1+h,rho[m])+sum)
			}
		    }
		    if(d[i] !=0 && d[j]==0){
			if(h <=0){
			    sum <- 0
			    for(l in 0:(-h)){
				sum <- sum + rho[m]^l*gamma(d[i]-h-l)*gamma(1-h)/gamma(d[i]-h)/gamma(1-h-l)
			    }
			    C[i,j,m,h+MaxH] <- gamma(d[i]-h)/gamma(d[i])/gamma(1-h)*
				(rho[m]^(2*K)*sum + HyperSum(d[i]-h,1,1-h,rho[m])-1)
			}
			if(h>0){
			    C[i,j,m,h+MaxH] <-rho[m]^h*HyperSum(d[i],1,1,rho)
			}
		    }
		}
	    }
	}
    }

    return(C)
}




### Function: ConditionNumber
## This function computes the condition number of a matrix,
## which is the ratio of the largest singular value to 
## the smallest singular value.
## Input: A - KxK matrix
ConditionNumber <- function(A){
    svds <- svd(A, nu=0, nv=0)$d
    return(svds[1]/svds[dim(A)[1]])
}



### Function: Sowell.AllMatrices
## This function computes and stores all the intermediate
## matrices from the Sowell recursions: A(n,k), v(n),
## D(n), Abar(n,k), vbar(n), and Dbar(n).  It returns all
## of these matrices as a list of lists of matrices.
## 5/31/08 (RS): This function now prints a warning if
## the condition number of any computed matrix exceeds a given
## value (1000 is the default).
## Inputs:
##	covarray - a KxKx(2n+1) array containing the
##		autocovariance matrices from lags -n to n
##	condCheck [1000] - the minimum condition number for which
##		this function should print a warning
## Output: a list containing
##	A - a list of n lists, with k matrices in the n-th list
##	v - a list of n+1 matrices
##	D - a list of n matrices
##	Abar - a list of n lists, with k matrices in the n-th list
##	vbar - a list of n+1 matrices
##	Dbar - a list of n matrices
Sowell.AllMatrices <- function(covarray, condCheck=1000){
    # Catching the case where K=1
    if(identical(NULL, dim(covarray))){
	dim(covarray) <- c(1,1,length(covarray))
    }

    # Get dimensions
    K <- dim(covarray)[1]
    n <- (dim(covarray)[3]-1)/2

    # Set up storage for matrix lists
    A <- vector("list",n)
    v <- vector("list",n+1)
    D <- vector("list",n)
    Abar <- vector("list",n)
    vbar <- vector("list",n+1)
    Dbar <- vector("list",n)


    # Initialize
    v[[1]] <- covarray[,,n+1]
    vbar[[1]] <- covarray[,,n+1]
    D[[1]] <- covarray[,,n+2]
    Dbar[[1]] <- covarray[,,n]
    A[[1]] <- vector("list",1)
    A[[1]][[1]] <- D[[1]] %*% solve(vbar[[1]])
    Abar[[1]] <- vector("list",1)
    Abar[[1]][[1]] <- Dbar[[1]] %*% solve(v[[1]])


    # Iterate
    for(i in 2:(n)){
   	# v(n)
	v[[i]] <- covarray[,,n+1]
	vbar[[i]] <- covarray[,,n+1]
	for(j in 1:(i-1)){
	    v[[i]] <- v[[i]] - A[[i-1]][[j]] %*% covarray[,,n+1-j]
	    vbar[[i]] <- vbar[[i]] - Abar[[i-1]][[j]] %*% covarray[,,n+1+j]
	}
	
	if(Im(ConditionNumber(v[[i]]))!= 0){
	    warning(paste("Condition number of v at location ", i, " equals ",ConditionNumber(v[[i]])))
	} else if(ConditionNumber(v[[i]])>condCheck){
	    warning(paste("Condition number of v at location ", i, " equals ",ConditionNumber(v[[i]])))
	} 

  	# D(n)
	D[[i]] <- covarray[,,n+1+i]
	Dbar[[i]] <- covarray[,,n+1-i]
	for(j in 1:(i-1)){
	    D[[i]] <- D[[i]] - A[[i-1]][[i-j]] %*% covarray[,,n+1+j]
	    Dbar[[i]] <- Dbar[[i]] - Abar[[i-1]][[i-j]] %*% covarray[,,n+1-j]
	}

	if(Im(ConditionNumber(D[[i]]))!= 0){
	    warning("Condition number of D at location ", i, " equals ",ConditionNumber(D[[i]]))
	} else if(ConditionNumber(D[[i]])>condCheck){
	    warning("Condition number of D at location ", i, " equals ",ConditionNumber(D[[i]]))
	} 


	# Set up storage list for A(n,k)
	A[[i]] <- vector("list",i)
	Abar[[i]] <- vector("list",i)


	# A(n,n)
	A[[i]][[i]] <- D[[i]] %*% solve(vbar[[i]])
	Abar[[i]][[i]] <- Dbar[[i]] %*% solve(v[[i]])

	# A(n,k)
	for(k in 1:(i-1)){
	    A[[i]][[k]] <- A[[i-1]][[k]] - A[[i]][[i]] %*% Abar[[i-1]][[i-k]]
	    Abar[[i]][[k]] <- Abar[[i-1]][[k]] - Abar[[i]][[i]] %*% A[[i-1]][[i-k]]
	}

	if(Im(ConditionNumber(A[[i]][[i]]))!= 0){
	    warning(paste("Condition number of A(i,i) at location ", i, " equals ",ConditionNumber(A[[i]][[i]])))
	} else if(ConditionNumber(A[[i]][[i]])>condCheck){
	    warning(paste("Condition number of A(i,i) at location ", i, " equals ",ConditionNumber(A[[i]][[i]])))
	} 
    }

    # The final v
    v[[n+1]] <- covarray[,,n+1]
    vbar[[n+1]] <- covarray[,,n+1]
    for(j in 1:(n+1-1)){
	v[[n+1]] <- v[[n+1]] - A[[n+1-1]][[j]] %*% covarray[,,n+1-j]
	vbar[[n+1]] <- vbar[[n+1]] - Abar[[n+1-1]][[j]] %*% covarray[,,n+1+j]
    }

    return(list(A=A,v=v,D=D,Abar=Abar,vbar=vbar,Dbar=Dbar))
}

### Function: Sowell.determinant
# This function computes the determinant of a covariance matrix
# using Sowell's method.  It can take in either an array of
# covariances (in which case it computes all the Sowell matrices
# first) or the list of Sowell matrices directly.
# Input: One of
#	covarray - a KxKx(2n+1) array containing the
#		autocovariance matrices from lags -n to n
# 	sowellMatrices - in the format given by
#		Sowell.allMatrices
# Output: the determinant of the corresponding covariance matrix
Sowell.determinant <- function(covarray=NULL, sowellMatrices=NULL, n=0){
    # Get the Sowell matrices if necessary
    if(identical(sowellMatrices, NULL))
	sowellMatrices <- Sowell.AllMatrices(covarray)

    # Compute the determinant as the product of the first n 
    # determinants of v 
    if(n==0)  n <- length(sowellMatrices$v) 
    det <- 1
    for(i in 1:n){
	det <- det*det(sowellMatrices$v[[i]])
    }

    return(det)
}

### Function: Sowell.logdeterminant
## This function computes the determinant of a covariance matrix
## using Sowell's method.  It can take in either an array of
## covariances (in which case it computes all the Sowell matrices
## first) or the list of Sowell matrices directly.
## Input: One of
##	covarray - a KxKx(2n+1) array containing the
##		autocovariance matrices from lags -n to n
##	sowellMatrices - in the format given by
##		Sowell.allMatrices
##	n [0] - the number of data points for which the determinant
##		should be calculated; if n=0, then this number is determined
##		from the dimension of the array of covariances of the length
##		of the lists in the Sowell matrix decomposition
## Output: the determinant of the corresponding covariance matrix
Sowell.logdeterminant <- function(covarray=NULL, sowellMatrices=NULL, n=0){
    # Get the Sowell matrices if necessary
    if(identical(sowellMatrices, NULL))
	sowellMatrices <- Sowell.AllMatrices(covarray)

    # Compute the determinant as the product of the first n 
    # determinants of v 
    if(n==0)  n <- length(sowellMatrices$v) 
    logdet <- 0
    for(i in 1:n){
	logdet <- logdet+log(det(sowellMatrices$v[[i]]))
    }

    return(logdet)
}

### Function: MatrixSqrt
## This function computes the square root of a matrix 
## (using the eigenvalue decomposition).
## Input: A - square matrix
## Output: A^(1/2)
MatrixSqrt <- function(A){
    eig <- eigen(A)

    return(eig$vectors %*% diag(sqrt(eig$values) ) %*% solve(eig$vectors))
}

### Function: MatrixInvSqrt
## This function computes the inverse of the square root of a matrix 
## (using the eigenvalue decomposition).
## Input: A - square matrix
## Output: A^(1/2)
MatrixInvSqrt <- function(A){
    eig <- eigen(A)

    return(eig$vectors %*% diag(1/sqrt(eig$values) ) %*% solve(eig$vectors))
}


### Function: Sowell.quadraticForm
## This function computes the quadratic form, X'Sigma^{-1}X, given 
## X and either the array of covariances that determine Sigma or 
## the Sowell matrices implied by Sigma
## Inputs:
##	X - vector of length nK, grouped by TIME (not series)
##	covarray [NULL] - array of covariances that describes Sigma
##	sowellMatrices [NULL] - list of list of matrices from the Sowell
##		algorithm
Sowell.quadraticForm <- function(X, covarray=NULL, sowellMatrices=NULL){
    # Get the Sowell matrices if necessary
    if(identical(sowellMatrices, NULL)){
	sowellMatrices <- Sowell.AllMatrices(covarray)
    }

    # Get dimensions
    K <- dim(sowellMatrices$v[[1]])[1]
    T <- length(X)/K

    # Using Result 2, we compute beta'X
    betaX <- vector(length=T*K)
    # Deal with i=0 separately
    betaX[1:K] <- MatrixInvSqrt(sowellMatrices$vbar[[1]]) %*% X[1:K]
    for(i in 2:T){
	betaX[((i-1)*K+1):(i*K)] <- X[((i-1)*K+1):(i*K)]
	for(j in 1:(i-1)){
	    betaX[((i-1)*K+1):(i*K)] <- betaX[((i-1)*K+1):(i*K)] - 
				sowellMatrices$Abar[[i-1]][[i-j]] %*% X[((j-1)*K+1):(j*K)] 
	}
	betaX[((i-1)*K+1):(i*K)] <- MatrixInvSqrt(sowellMatrices$vbar[[i]]) %*% betaX[((i-1)*K+1):(i*K)]
    }

    return(t(betaX) %*% betaX)
}


### Function: convolution
## This function takes in two sequences, possibly of different lengths,
## and computes the (non-circular) convolution, sum(x(j)y(m-j)), for all
## j where the two are defined and all m.  
## Inputs: x - vector of length J
##	  y - vector of length K
## Output: All convolutions in which at least one product, x(j)y(m-j), is defined.
convolution <- function(x,y){
    J <- length(x)
    K <- length(y)

    # Using direct sums if one series is much longer than the other
    if(K < log(J)){
	result <- rep(0, J+K-1)

	# Deal with initial values
	for(i in 1:(K-1)){
	    result[i] <- sum(x[1:i]*rev(y[1:i]))
 	}
	# Deal with the middle
	for(i in K:J){
	    result[i] <- sum(x[(i-K+1):i]*rev(y))
	}
	# Deal with final values
	for(i in 1:(K-1)){
	    result[J+i] <- sum(x[(J-i+1):J]*rev(y[(K-i+1):K]))
	}
    } else if(K < log(J)) {
	result <- rep(0, J+K-1)

	# Deal with initial values
	for(i in 1:(J-1)){
	    result[i] <- sum(x[1:i]*rev(y[1:i]))
 	}
	# Deal with the middle
	for(i in J:K){
	    result[i] <- sum(x*rev(y[(i-J+1):i]))
	}
	# Deal with final values
	for(i in 1:(J-1)){
	    result[J+i] <- sum(x[(J-i+1):J]*rev(y[(K-i+1):K]))
	}

    } else {
    	# If the series are of comparable lengths, we use FFT:
    	# Total length should be at least K+J to avoid overlap; we move up to 
    	# the next power of two, for optimizing speed
    	longLength <- 2^(floor(log(K+J)/log(2))+1)
    
    	longX <- c(x,rep(0,longLength-J))
    	longY <- c(y,rep(0,longLength-K))

    	# Compute the convolution
    	longResult <- fft(fft(longX)*fft(longY),inverse=TRUE)/longLength

    	# Extract the non-zero convolutions
    	# (there are J+K-1 of these; they are the first terms)
    	result <- Re(longResult[1:(J+K-1)])
    }

    return(result)
}


### Function: VAR.cov0
## This function computes the unconditional covariance (autocovariance at
## lag 0) of an VAR(1) process described by KxK lag coefficient matrix 
## F and KxK error covariance matrix Sigma.
## Inputs:
##	F - KxK matrix of coefficients on lagged variables
##	Sigma - KxK symmetric and positive definite covariance matrix of the errors
## Output: KxK matrix with the unconditional covariance.
VAR.cov0 <- function(F,Sigma){
    K <- dim(Sigma)[1]

    # Vectorize Sigma
    vecSigma <- Sigma
    dim(vecSigma)<- NULL

    result <- solve(diag(K*K)-kronecker(F,F),vecSigma)

    # Un-vectorize the result
    dim(result)<-c(K,K)

    return(result)
}    

### Function: VAR.autoCov
## This function computes the first M autocovariances of a 
## VAR(1) process described by KxK lag coefficient matrix 
## F and KxK error covariance matrix Sigma. 
## Inputs:
##	F - KxK matrix of coefficients of lagged variables
##	Sigma - KxK symmetric and positive definite covariance matrix of the errors
##	M - number of autocovariances returned (from 0 to M-1)
##	known0 [NULL] - autocovariance of lag 0, if already known
## Output: a KxKx(M+1) array, in which [.,.,r] is the r-th autocovariance matrix
VAR.autoCov <- function(F,Sigma,M, known0=NULL){
    # Create array to store results
    autocovs <- array(dim=c(dim(F),M+1))

    # Find the autocovariance at the 0'th lag
    if(identical(known0,NULL)){
    	autocovs[,,1] <- VAR.cov0(F,Sigma)
    } else {
	autocovs[,,1]<-known0
    }

    # Iteratively compute other autocovariances
    for(i in 1:M)
	autocovs[,,i+1]<-F %*% autocovs[,,i]

    return(autocovs)
}

### Function: VARtoVAR1
## This function takes in the coefficients and innovation variance
## of a VAR(p) and returns the coefficients and innvovation variance
## of the VAR(1) to which it corresponds.
## Inputs:
## 	Alist <- list of KxK matrices, where the first matrix
##		is the coefficient on lag t-1, and the 
##		last matrix is the coefficient on lag t-p
##	Sigma <- innovation variance of the original process
## Output: a list containing:
##	bigA <- Kp x Kp matrix which is the matrix on lag t-1 in the
##		augmented model
##	bigSigma <- Kp x Kp matrix which is the innovation variance of
##		the augmented model
VARtoVAR1 <- function(Alist, Sigma){
    # Determine p and K
    p <- length(Alist)
    K <- dim(Sigma)[1]

    # Set up results
    bigA <- matrix(data=0,nrow=K*p,ncol=K*p)
    bigSigma <- matrix(data=0,nrow=K*p,ncol=K*p)
    
    # Populate bigA
    bigA[(K+1):(p*K),1:(K*(p-1))] <- diag((p-1)*K)
    for(i in 1:p){
	bigA[1:K,((i-1)*K+1):(i*K)] <- Alist[[i]]
    }

    bigSigma[1:K,1:K] <- Sigma

    return(list(bigA=bigA,bigSigma=bigSigma))
}



### Function: ARFIMA.MAcoeff
## This function computes the j-th moving average coefficient
## in the MA infinity expansion of an ARFIMA(0,d,0) model.
## The formula is:
## Gamma(j-d)/Gamma(j+1)Gamma(d)
## Inputs:
##	j - which coefficient in the MA expansion
##	d - the differencing parameter of the ARFIMA process
ARFIMA.MAcoeff <- function(j,d){
    # If d=0, this is white noise
    if(d==0){
	if(length(j)==1){
	    if(j==0) return(1)
	    return(0)
	}
     	result <- rep(0,length(j))
	for(i in 1:length(j)){
	    if(j[i]==0) result[i]<-1
	}
	return(result)
    }


    # When j is a vector, we work one-at-a-time
    result <- rep(0,length(j))
    for(i in 1:length(j)){
	if(j[i]>=0){
		if(j[i]+d <0){
    	    		lresult <- lgamma(j[i]+1)
    	    		result[i] <- exp(lresult)/gamma(d)*gamma(j[i]+d)
		} else {
    	    		lresult <- lgamma(j[i]+d)-lgamma(j[i]+1)
    	    		result[i] <- exp(lresult)/gamma(d)
		}
	}
    }
    return(result)
}

### Function: ARFIMA.autoCov
## This function computes the lag-h autocovariance of an 
## ARFIMA(0,d,0) model.
## Inputs:
##	h - lag of the autocovariance
##	d - differencing parameter
##	sigmasq [1] - innovation variance
ARFIMA.autoCov <- function(h,d, sigmasq=1){
    # If d=0, this is white noise
    if(d==0){
	if(length(h)==1){
	    if(h==0) return(sigmasq)
	    return(0)
	}
     	result <- rep(0,length(h))
	for(i in 1:length(h)){
	    if(h[i]==0) result[i]<-sigmasq
	}
	return(result)
    }

    h <- abs(h)
    lresult <- lgamma(1-2*d)-lgamma(h-d+1)-lgamma(1-h-d)
    return(exp(lresult)*sigmasq)
}

### function: ARFIMA.crossCov
## This computes the cross-covariance of ARFIMA(0,d,0) processes driven 
## by the same white noise, using Sowell's (1989) expression.
## Inputs:
##	d1 - differencing parameter for the first process
##	d2 - differencing parameter for the second process
## 	m - the autocovariance lag to compute 
##	sigmasq [1] - the variance of the innovations
# 8/14/09 (RS): Using lgamma doesn't work when the gamma function should give a negative
# answer.  When d1 or d2 is negative, we handle those cases separately). (Also,
# m < 0 can cause problems, and we're fixed that, too.)
ARFIMA.crossCov <- function(d1, d2, m, sigmasq=1, verbose=FALSE){
    # Ensure that d1 and d2 imply a stationary model (NOT CHECKING INVERTIBILITY)
    if(max(d1,d2) >= 0.5) stop(paste("Cannot compute covariances for a model with differencing parameter of 0.5 or greater; value of ", max(d1,d2), " given"))


    # If d1 or d2 is 0, the result the the MA coefficients of the
    # other differencing parameter
    # To deal with cases when d1 or d2 is very close to 0, we also 
    # check whether lgamma(d1) or lgamma(d2) is Inf.
    if(d1==0 || identical(lgamma(1-d1-1),Inf)){
	return(ARFIMA.MAcoeff(m,d2))
    }
    if(d2==0 || identical(lgamma(1-d2-1),Inf)){
	return(ARFIMA.MAcoeff(-1*m,d1))
    }


    # Pass through the m's one at a time:
    # If m is positive, use the formula as is (using gamma instead of
    # lgamma if m+d1 < 0)
    # If m is negative, switch d1 and d2 and the sign of m.
    result <- vector(length=length(m))
    for(i in 1:length(m)){
	if(m[i] < 0){
		thism <- -1*m[i]
		thisd1 <- d2
		thisd2 <- d1
	} else {
	 	thism <- m[i]
		thisd1 <- d1
		thisd2 <- d2
	}

    	# When d1 is negative, we have to be more careful with the lgamma function
    	if(thisd1>0){
        	lresult <- log(sigmasq)+lgamma(1-thisd1-thisd2)+lgamma(thisd1+thism)-lgamma(thisd1)-lgamma(1-thisd1)-lgamma(1+thism-thisd2)
        	result[i] <- exp(lresult)
    	} else {
    		multiplier <- sigmasq
	    	logpart <- lgamma(1-thisd1-thisd2)-lgamma(1-thisd1)-lgamma(1+thism-thisd2)
    		if(thisd1+thism>0){
			multiplier <- multiplier/gamma(thisd1)
			logpart <- logpart +lgamma(thisd1+thism)
     		} else {
			multiplier <- multiplier/gamma(thisd1)*gamma(thisd1+thism)
     		}
		result[i] <- multiplier*exp(logpart)
	}
    }

    return(result)
}


### Function: FIVAR.autoCov0
## This function takes in the parameters of a FIVAR(0,d,0) process
## and returns the specified number of autocovariances
## Inputs:
##	Sigma - KxK symmetric and positive definite covariance matrix of the errors
##	d - vector of length K of differencing parameters
##	n - number of lags returned (from -(n-1) to n-1)
## Output: a K x K x n array, with the r-th autocovariance at [.,.,r-1]
FIVAR.autoCov0 <- function(Sigma,d,n,verbose=FALSE){
    K <- dim(Sigma)[1]

    # Get the ARFIMA covariances
    ARFIMAcov <- array(dim=c(K,K,2*n-1))
    for(i in 1:K){
	for(j in 1:K){
	    ARFIMAcov[i,j,] <- Sigma[i,j]*ARFIMA.crossCov(d[i],d[j],(-n+1):(n-1))
	}
    }
    return(ARFIMAcov)
}



### Function: VARFI.autoCov0
## This function takes in the parameters of a VARFIMA(0,d,0) process
## and returns the specified number of autocovariances. (Identical to 
## FIVAR.autoCov0, since they are the same when there is no AR part.)
## Inputs:
##	Sigma - KxK symmetric and positive definite covariance matrix of the errors
##	d - vector of length K of differencing parameters
##	n - number of lags returned (from -(n-1) to n-1)
## Output: a K x K x n array, with the r-th autocovariance at [.,.,r-1]
VARFI.autoCov0 <- function(Sigma,d,n,verbose=FALSE){
    K <- dim(Sigma)[1]

    # Get the ARFIMA covariances
    ARFIMAcov <- array(dim=c(K,K,2*n-1))
    for(i in 1:K){
	for(j in 1:K){
	    ARFIMAcov[i,j,] <- Sigma[i,j]*ARFIMA.crossCov(d[i],d[j],(-n+1):(n-1))
	}
    }
    return(ARFIMAcov)
}



### function: SmallBlockCircEmbed
# This function takes in a KxKx(2T-1) array describing a 
# block Toeplitz matrix and returns a KxKx(2T-1) array 
# describing the block circulant matrix in which it can be
# embedded.  Unlike circembed in the multipcg algorithm,
# this circulant embedding does NOT add an additional diagonal
# containing the value on the main diagonal.
SmallBlockCircEmbed <- function(toeplitzarray){
    oldsize <- (dim(toeplitzarray)[3]+1)/2

    result <- toeplitzarray

    result[,,1:oldsize] <- toeplitzarray[,,oldsize:(2*oldsize-1)]
    result[,,(oldsize+1):(2*oldsize-1)] <- toeplitzarray[,,1:(oldsize-1)]

    return(result)
}



### Function: SimulateFourierNormal
# This function simulates N normal random variables with covariance matrix
# equal to the Fourier matrix F*F.
# Input: n - number of random variables to simulate
# Output: a complex vector of length n whose components have covariance
#	matrix equal to FF*
SimulateFourierNormal <- function(n){
    U <- rnorm(ceiling(n/2))
    V <- rnorm(floor(n/2))

    result <- vector(length=n)
    result[1] <- U[ceiling(n/2)]/sqrt(n)
    if(n %% 2 ==0){
	result[n/2+1] <- V[n/2]
    }
    for(i in 1:(ceiling(n/2)-1)){
	result[i+1] <- (U[i]+1i*V[i])/sqrt(2*n)
	result[n-i+1] <- (U[i]-1i*V[i])/sqrt(2*n)
    }

    return(result)
}

### Function: SortByTime
# This function takes in a vector of length K*T, assumed to be grouped 
# by series, and groups it by time. 
# Inputs:
#	x - vector of length K*T
#	K - number of series
#	T - number of periods per series
# Output: the elements of x in a different order
SortByTime <- function(x, K, T){
    # Check dimensions
    if(length(x) != K*T) stop("SortByTime: Length of vector does match product of given dimensions")

    sortX <- x
    for(k in 1:K){
	for(t in 1:T){
	    sortX[(t-1)*K+k] <- x[(k-1)*T+t]
	}
    }

    return(sortX)
}

### Function: SortBySeries
# This function takes in a vector of length K*T, assumed to be grouped 
# by time, and groups it by series. 
# Inputs:
#	x - vector of length K*T
#	K - number of series
#	T - number of periods per series
# Output: the elements of x in a different order
SortBySeries <- function(x, K, T){
    # Check dimensions
    if(length(x) != K*T) stop("SortByTime: Length of vector does match product of given dimensions")

    sortX <- x
    for(k in 1:K){
	for(t in 1:T){
	    sortX[(k-1)*T+t] <- x[(t-1)*K+k]
	}
    }

    return(sortX)
}


## Function: ListOfDeterminants
# This function takes in a list of matrices and returns a vector containing their
# determinants.
ListOfDeterminants <- function(matlist){
    # Determinant doesn't work if the matrices are 1x1
    if(length(matlist[[1]])==1){
	result <- unlist(matlist)
    } else {
    	n <- length(matlist)
    	result <- vector(length=n)
	for(i in 1:n) 
	    result[i] <- det(matlist[[i]])
    }

    return(result)
}

### Function: PCGpredvar
# This function uses PCG to compute the prediction variance (v)
# for a particular n.  The covariances can be specified either by
# d and Sigma (the differencing parameters and innovation variance 
# of an ARFIMA(0,d,0)) or by the array of autocovariances from -n to n.
# Inputs:
# n - the lag of the prediction variance
# d [NULL] - the vector of differencing parameters
# Sigma [NULL] - the innovation variance
# covarray [NULL] - a KxKx(2n+1) array containing the autocovariances from
#	lags -n to n
# Output: The prediction variance matrix, v
PCGpredvar <- function(n, d=NULL, Sigma=NULL, F=NULL, covarray=NULL){
    # Create the array of covariances if necessary
    if(identical(covarray, NULL)){
	covarray <- FIVAR.autoCov(Sigma=Sigma, d=d, n=n+1, F=F)
    }

    # Catching the case where K=1
    if(identical(NULL, dim(covarray))){
	dim(covarray) <- c(1,1,length(covarray))
    }

    # Extract dimensions
    K <- dim(covarray)[1]

    # Set up Gamma matrices
    Gamma <- matrix(nrow=K*n,ncol=K)
    for(i in 1:K){
	for(j in 1:K){
	    Gamma[((i-1)*n+1):(i*n),j] <- covarray[i,j,(n):(1)]
	}
    }

    # Set up a representation for Sigma by deleting the first and
    # last covariances
    Sigmarowcol <- covarray[1:K,1:K,2:(2*n)]
    if(n==1) dim(Sigmarowcol) <- c(K,K,1)

    # Compute Sigma^(-1)*Gamma
    product <- matrix(nrow=K*n,ncol=K)
    for(i in 1:K){
	product[,i] <- MultivariatePCG(Sigmarowcol, Gamma[,i])
    }

    # Compute Gamma'*Sigma^(-1)*Gamma
    product2 <- t(Gamma) %*% product

    # v = v(0)-product2
    return(covarray[,,n+1] - product2)

}


### Function: RegressionApprox
# This function uses Sowell to compute some initial v-determinants and
# PCG to compute some later ones.  Then, it runs the regression of 
# sqrt(v(j))*j^2 on j and uses this function to interpolate the
# remaining determinants.
# It returns the entire sequence 
# of approximate determinants (which will
# presumably be added together to get the answer).  
# Inputs:
#	n - Maximum determinant to be calculated
# 	covarray - a KxKx(2n+1) array containing the autocovariances from
#		lags -n to n
#	initialSowell [32] - how many Sowell matrices should be computed for the initial part of the regression
#	PCGpoints [n-1] - the points at which PCG should be used to compute the exact value
# Output: logdets - a vector of prediction variance determinants from 0 to n-1
RegressionApprox <- function(n, covarray, verbose=FALSE, initialSowell=min(32,n-1), PCGpoints=c(n-1)){
    middle <- (dim(covarray)[3]+1)/2

    # Set up storage
    logdets <- vector(length=n)

    # Compute initial Sowell matrices until the errors are less 
    # than maxPredVarError
    sowellMatrices <- Sowell.AllMatrices(covarray[,,(middle-initialSowell):(middle+initialSowell)])

    # Copy the log determinants into the result
    logdets[1:(initialSowell+1)] <- log(ListOfDeterminants(sowellMatrices$v))

    # Return this if we have computed all the determinants
    if(initialSowell==n-1) return(logdets)

    # Storage for the points for the correction regression
    fitX <- 1:(initialSowell)
    originalY <- logdets[2:(initialSowell+1)]

    # PCG for final point and any intermediate points
    fitX <- c(fitX, PCGpoints)

    for(i in 1:length(PCGpoints)){
    	originalY <- c(originalY,
		log(det(PCGpredvar(PCGpoints[i],covarray=covarray[,,(middle-PCGpoints[i]):(middle+PCGpoints[i])]))))
    }

    # Fit line using this kappa
    fitY <- sqrt(exp(originalY)*fitX^2)
	
    # Create data frames for fitting and prediction
    frame <- data.frame(fitX,fitY)

    # Fit the model
    estmodel <- lm(fitY~fitX, data=frame)
    if(verbose){ 
    	print(paste("Estimated Line for n times the square root of the determinant"))
	print(summary(estmodel))
    }

    # Predict for all the other points
    logdets[(initialSowell+1):n] <- 2*(log(predict(estmodel,data.frame("fitX"=initialSowell:(n-1))))-log(initialSowell:(n-1)))

    return(logdets)
}

	
	
### Function: FIVAR.logdeterminant
# This function computes the log determinant of the covariance
# matrix of n realizations of a vector ARFIMA(0,d,0) process,
# to a given level of accuracy (in log terms).
# Inputs:
#	Sigma - KxK covariance of the white noise
#	d - vector of differencing parameters of length K
#	n - sample size
#	F [NULL] - KxK VAR(1) matrix (if any)
#	tol [0.001] - error tolerance in the determinant
# Output: the determinant of the Kn x Kn covariance matrix 
#	of all the observations
FIVAR.logdeterminant <- function(n, covarray=NULL, Sigma=NULL, d=NULL,  F=NULL, tol=0.001,verbose=FALSE){
    # Get all the autocovariances, if necessary (this requires n of them, not n-1)
    if(identical(covarray,NULL)){
    	covarray <- FIVAR.autoCov(F, Sigma, d, n+1)
    } 

    # Get the list of terms of the determinant
    predvars <- RegressionApprox(n,covarray,verbose=verbose)

    return(sum((predvars)))
}

### Function: VARFI.logdeterminant
# This function computes the log determinant of the covariance
# matrix of n realizations of a VARFI(1,d) process
# using the regression approximation.
# NOTE: FAILS IF n<4.  SHOULD PROBABLY FIX THIS.
# Inputs:
#	Sigma - KxK covariance of the white noise
#	d - vector of differencing parameters of length K
#	n - sample size
#	F [NULL] - KxK VAR(1) matrix (if any)
#	tol [0.001] - error tolerance in the determinant
# Output: the determinant of the Kn x Kn covariance matrix 
#	of all the observations
VARFI.logdeterminant <- function(n, covarray=NULL, Sigma=NULL, d=NULL,  F=NULL, tol=0.001,verbose=FALSE){
    # Get all the autocovariances, if necessary (this requires n of them, not n-1)
    if(identical(covarray,NULL)){
    	covarray <- VARFI.autoCov(F, Sigma, d, n+1)
    } 

    # Get the list of terms of the determinant
    predvars <- RegressionApprox(n,covarray,verbose=verbose)

    return(sum(predvars))
}

### Global variables to set
maxiters <- 100
maxerror <- .Machine$double.eps



### function: fftmultwitheigen
# This function does an FFT multiplication assuming that the eigenvalues
# of a circulant matrix are given.
# Inputs: eigen - Tx1 vector of eigenvalues of a circulant matrix
#	  x - Tx1 vector
fftmultwitheigen <- function(eigen, x){
    n <- length(eigen)
    product <- fft(eigen * fft(x),inverse=TRUE)/n
    return(Re(product))
}

### function: blockcircmult
# This function multiplies a block circulant matrix by a vector,
# given the eigenvalues of each block.
# Inputs: eigen - a K x K x T matrix, with the eigenvalues of block
#		  [i,j] in location [i,j,*]
#	  v - KT x 1 vector to be multiplied
blockcircmult <- function(eigen,v){
    # Extract sizes
    K <- dim(eigen)[2]
    T <- dim(eigen)[3]

    product <- rep(0,K*T)
    for(i in 1:K){
	for(j in 1:K){
	    product[((i-1)*T+1):(i*T)] <- product[((i-1)*T+1):(i*T)] +
			fftmultwitheigen(eigen[i,j,1:T],v[((j-1)*T+1):(j*T)])
	}
    }

    return(product)
}


### function: blocktoepmult
# This function multiplies a block Toeplitz matrix by a vector,
# given the eigenvalues of the circulant embedding of each block.
# Inputs: eigen - a K x K x 2T matrix, with the eigenvalues of block
#		  [i,j] in location [i,j,*]
#	  v - KT x 1 vector to be multiplied
# NO SIZE CHECKING
blocktoepmult <- function(eigen, v){
    K <- dim(eigen)[2]
    T <- (length(v))/K

    # Add zeroes to v to correspond to circulant embedding
    longv <- rep(0,2*K*T)
    for(i in 1:K){
	longv[((i-1)*2*T+1):((i-1)*2*T+T)]<-v[((i-1)*T+1):(i*T)]
    }

    longproduct <- blockcircmult(eigen,longv)

    # Extract the corresponding part of v
    product <- rep(0,K*T)
    for(i in 1:K){
	product[((i-1)*T+1):(i*T)]<-longproduct[((i-1)*2*T+1):((i-1)*2*T+T)]
    }
    
    return(product)
}


### function:circembedFromFirstRow
# This function takes the first row of a SYMMETRIC Toeplitz matrix
# and returns the first row
# of its embedding in a circulant.
# Input: firstA - the Tx1 vector with the first row of the Toeplitz matrix
circembedFromFirstRow <- function(firstA){
    n <- length(firstA)
    result <- rep(0, 2*n)
    result[1:n] <- firstA
    result[n+1] <- firstA[1]
    for(i in 2:n){
	result[(n+i)] <- firstA[n-abs(i-1)+1]
    }

    return(result)
}

### function:circembedFromFirstRowAndColumn
# This function takes the first row and column of a Toeplitz matrix
# and returns the first row
# of its embedding in a circulant.
# 10/23/07 (RS): Fixed in the special case of n=1.
# Input: firstRowA - the Tx1 vector with the first row of the Toeplitz matrix
# 	     firstColA - the Tx1 vector with the first column of the Toeplitz matrix
circembedFromFirstRowAndColumn <- function(firstRowA,firstColA){
    n <- length(firstRowA)
    # Special case when n=1
    if(n==1) return(c(firstRowA,firstColA))

    result <- rep(0, 2*n)
    result[1:n] <- firstRowA
    result[n+1] <- firstRowA[1]
    for(i in 2:n){
	result[(n+i)] <- firstColA[n-abs(i-1)+1]
    }

    return(result)
}

### function: toepmultFromFirstRow
# This function multipies a Toeplitz matrix by a vector.
# (This does the eigenvalue computation; eigenvalues do not need
# to be passed in.)
# Inputs: firstA -  the first row of a SYMMETRIC Toeplitz matrix
#	  v - N a 1 vector
# Output: A*v
toepmultFromFirstRow <- function(firstA,v){
    N <- length(v)
    circA <- circembedFromFirstRow(firstA)
    longv <- rep(0,2*N)
    longv[1:N] <- v
    longresult <- fftmultwitheigen(fft(circA),longv)
    return(longresult[1:N])
}


### function: circapproxFromFirstRow
# This function finds the first row of the circulant approximation
# of a SYMMETRIC Toeplitz matrix, in which only the first row is given.
# 10/23/07 (RS): Fixed in the special case of n=1.
# Input: firstRowA - the Tx1 vector with the first row of the Toeplitz matrix
circapproxFromFirstRow <- function(firstA){
    n <- length(firstA)
    if(n==1) return(firstA)

    c <- rep(0, n)
    c[1] <- firstA[1]
    for(i in 2:n){
	c[i] <- ((i-1)*firstA[n-i+2]+(n-i+1)*firstA[i])/n
    }
    return(c)

}

### function: circapproxFromFirstRowAndColumn
# This function finds the first row of the circulant approximation
# of a Toeplitz matrix, which is specified by the first row and the
# first column.
# 10/23/07 (RS): Fixed in the special case of n=1.
# Input: firstRowA - the Tx1 vector with the first row of the Toeplitz matrix
# 	     firstColA - the Tx1 vector with the first row of the Toeplitz matrix
circapproxFromFirstRowAndColumn <- function(firstRowA, firstColA){
    n <- length(firstRowA)
    if(n==1) return(firstRowA)

    c <- rep(0, n)
    c[1] <- firstRowA[1]
    for(i in 2:n){
	c[i] <- ((i-1)*firstColA[n-i+2]+(n-i+1)*firstRowA[i])/n
    }
    return(c)

}

### function: multicg
# This function performs the CG algorithm for any multivariate model. 
# It takes in the covariance matrix and inverted preconditioner, each already
# in eigenvalue form, and the vector of observations (Y), and returns an 
# estimate of the signal.
# Inputs:
# Aeigen - K x K x 2T matrix containing the eigenvalues of the circulant emedding
#	   of the covariance matrix blocks
# Y - the vector in the linear equation Ax=Y
# verbose - TRUE if number of iterations should be printed (FALSE by default)
# Output:  the TKx1 solution
multicg <- function(Aeigen, Y, verbose=FALSE){
    # Compute dimensions
    K <- dim(Aeigen)[1]/2
    T <- length(Y)/K

    # Size checking
	if(dim(Aeigen)[2] != 2*K | dim(Aeigen)[3] !=T){
		stop("Dimensions of eigenvalues of matrix and vector not compatible")
	}

    # Use CG (no preconditioning here!)
    x <- rep(0,K*T)
    r <- Y
    d <- r
    counter <- 0
    # We compute the error as r'*r
    error <- t(r) %*% r
    while(counter<maxiters && error>maxerror){
	# Compute A*d only once
	Ad <- blocktoepmult(Aeigen,d)
	alpha <- error/(t(d) %*% Ad)
	x <- x + alpha[1]*d
	rnew <- r - alpha[1]*Ad
	errornew <- t(rnew) %*% rnew
	beta <- errornew/error
	d <- rnew + beta[1]*d
	r <- rnew
	error <- errornew[1]
	counter <- counter+1
    }


    if(verbose) print(c("Number of iterations:",counter))
    return(x)   
}


### function: multipcg
# This function performs the PCG algorithm for any multivariate model 
# It takes in the covariance matrix and inverted preconditioner, each already
# in eigenvalue form, and the vector of observations (Y), and returns an 
# estimate of the signal.
# Inputs:
# Aeigen - K x K x 2T matrix containing the eigenvalues of the circulant emedding
#	   of the covariance matrix blocks
# lambdainv - K x K x T matrix containing the eigenvalues of the circulant
# 	      approximations to the covariance matrix blocks
# Y - the vector in the linear equation Ax=Y
# verbose - TRUE if number of iterations should be printed (FALSE by default)
# stopOnNegative [FALSE] - if the algorithm encounters a negative error, should it stop?
#	(if not, a warning is printed, and the result is returned)
# Output:  the TKx1 solution
multipcg <- function(Aeigen, lambdainv, Y, verbose=FALSE, stopOnNegative=FALSE){
    # Compute dimensions
    K <- dim(Aeigen)[1]
    T <- dim(lambdainv)[3]

    # Size checking
	if(dim(Aeigen)[2] != K | dim(Aeigen)[3] != 2*T){
		stop("Dimensions of eigenvalues of matrix and preconditioner not compatible")	}
	if(dim(lambdainv)[1] != K | dim(lambdainv)[2] != K){
		stop("Dimensions of eigenvalues of matrix and preconditioner not compatible")	}
	if(length(Y) != K*T){
		stop("Dimension of vector incompatible with dimensions of matrix")
	}

    # Use PCG with this preconditioner
    x <- rep(0,K*T)
    r <- Y
    d <- blockcircmult(lambdainv,r)
    counter <- 0
    # We compute the error as r'*Cinv*r
    error <- t(r) %*% d
    while(counter<maxiters && error>maxerror){
	# Compute A*d only once
	Ad <- blocktoepmult(Aeigen,d)
	alpha <- error/(t(d) %*% Ad)
	x <- x + alpha[1]*d
	rnew <- r - alpha[1]*Ad
	# Compute Cinv * r only once
	Cinvrnew <- blockcircmult(lambdainv,rnew)
	errornew <- t(rnew) %*% Cinvrnew
	beta <- errornew/error
	d <- Cinvrnew + beta[1]*d
	r <- rnew
	error <- errornew[1]
	if(error < 0){
	    if(stopOnNegative){
	    	stop("PCG algorithm has error < 0; preconditioner is not positive definite.")
	    } else {
	    	warning("PCG algorithm has error < 0; preconditioner is not positive definite.")
	    }
	}
	counter <- counter+1
    }


    if(verbose==TRUE) print(c("Number of iterations:",counter))
    return(x)   
}


### function: multipcgModel1
# This function performs the PCG algorithm for any signal+noise model 
# like the one specified in the paper.  It takes in the autocovariances 
# (autocov) of the signal, the covariance matrix (Tau) of the noise, the vector of
# factor loadings (C), and the vector of observations (Y), and returns an 
# estimate of the signal.
# Inputs:
# autocov - the vector of autocovariances of the signal (length T)
# Tau - covariance matrix of noise (K x K)
# C - factor loadings (K x 1)
# Y - observations (KT x 1, observations grouped by series)
# Output: X (T x 1)
# NO SIZE CHECKING 
multipcgModel1 <- function(autocov, Tau, C, Y){
    # Compute dimensions
    K <- dim(Tau)[1]
    T <- length(autocov)

    # We pass through each block and (based on the first row):
    # Compute the eigenvalues for circulant embedding
    # Compute the circulant approximation and its eigenvalues
    Aeigen<- rep(0, 2*K*T*K)
    dim(Aeigen) <- c(K,K,2*T)
    lambda<- rep(0, K*T*K)
    dim(lambda) <- c(K,K,T)
   
    # A vector with 1 in the first location and 0 elsewhere (for the
    # variance contribution of the noise)
    NoiseEffect <- rep(0,T)
    NoiseEffect[1] <- 1

    for(i in 1:K){
	for(j in 1:K){
     	    currblock <- C[i]*C[j]*autocov+Tau[i,j]*NoiseEffect
	    Aeigen[i,j,] <- fft(circembedFromFirstRow(currblock))
	    lambda[i,j,] <- fft(circapproxFromFirstRow(currblock))
	}
    }
    
    ### Finish preconditioner
    # The blocks along the diagonal from Chan & Olkin
    # are the submatrices [1:K,1:K,t].
    # Using built-in inversion for these.
    lambdainv <- lambda
    for(t in 1:T){
	lambdainv[,,t] <- solve(lambda[,,t])
    }

    return(multipcg(Aeigen,lambdainv,Y))
}

### function: multipcgModel2                      
# This function performs the PCG algorithm for any signal+noise model 
# like the one specified in the paper, in which the signal is differenced
# to make it stationary (and therefore the noise is non-invertible).  
# It takes in the covariance matrix 
# (Sigma) of the signal, the covariance matrix (Tau) of the (undifferenced) 
# noise, the vector of
# factor loadings (C), and the vector of observations (Y), and returns an 
# estimate of the signal.
# Inputs:
# autocov - the vector of autocovariances of the signal (length T)
# Tau - covariance matrix of (undifferenced) noise (K x K)
# C - factor loadings (K x 1)
# Y - observations (KT x 1, observations grouped by series)
# Output: X (T x 1)
# NO SIZE CHECKING 
multipcgModel2 <- function(autocov, Tau, C, Y){
    # Compute dimensions
    K <- dim(Tau)[1]
    T <- length(autocov)

    # We pass through each block and:
    # Compute the eigenvalues for circulant embedding
    # Compute the circulant approximation and its eigenvalues
    Aeigen<- rep(0, 2*K*T*K)
    dim(Aeigen) <- c(K,K,2*T)
    lambda<- rep(0, K*T*K)
    dim(lambda) <- c(K,K,T)
   
    # A vector with 2 in the first location, -1 in the second
    # location, and 0 elsewhere (for the
    # variance contribution of the differenced noise)
    NoiseEffect <- rep(0,T)
    NoiseEffect[1] <- 2
    NoiseEffect[2] <- -1

    for(i in 1:K){
	for(j in 1:K){
     	    currblock <- C[i]*C[j]*autocov+Tau[i,j]*NoiseEffect
	    Aeigen[i,j,] <- fft(circembedFromFirstRow(currblock))
	    lambda[i,j,] <- fft(circapproxFromFirstRow(currblock))
	}
    }


    
    ### Finish preconditioner
    # The blocks along the diagonal from Chan & Olkin
    # are the submatrices [1:K,1:K,t].
    # USING BUILT-IN INVERSION FOR THESE.
    lambdainv <- lambda
    for(t in 1:T){
	lambdainv[,,t] <- solve(lambda[,,t])
    }

    return(multipcg(Aeigen,lambdainv,Y))
}





############ Functions for General Use ########################

### Function: MultivariateCG
# This is the general wrapper function for the multivariate conjugate
# gradient algorithm.  This function solves the system Ax=b for x, when
# A is a block Toeplitz matrix (with K^2 blocks of size TxT).
# This takes in a K x K x T or K x K x (2T-1) matrix 
# containing the first row OR the first column and row of each block
# of the block Toeplitz matrix.  (If the first row is specified, then
# the matrix is assumed to be symmetric.  Otherwise, the first row and
# column must be specified, in this order:
# a(T,1), a(T-1,1), ..., a(1,1), a(1,2), ..., a(1,T)
# This corresponds to giving the covariances from lag -(T-1) to lag (T-1)
# in a time series context.)
# Inputs:
# firstrowcol - a KxKxT or KxKx(2T-1) matrix, in which the vector 
#			describing the first row (and possibly 
#			column) of the (i,j) block is in (i,j,:)
# Y - KT x 1 vector
# symmetric [FALSE] - indicator of whether the the Toeplitz matrix
#		      is being specified by only the first row, since
#		      it is symmetric
# verbose [FALSE] - should the number of iterations be printed?
# Output: X (T x 1)
MultivariateCG <- function(firstrowcol, Y, symmetric=FALSE,verbose=FALSE){
    # Compute dimensions
    K <- dim(firstrowcol)[1]
    T <- length(Y)/K

    # Size checking
	if(dim(firstrowcol)[2] !=K){
		stop("Number of blocks per row and per column are not equal.")
	}
	if(symmetric & dim(firstrowcol)[3] !=T){
		stop("Size of vector and size of blocks incompatible")
	}
	if(!symmetric & dim(firstrowcol)[3] !=(2*T-1)){
		stop("Size of vector and size of blocks incompatible")
	}
	
	# Separate rows and columns
	if(symmetric){
		firstrow <- firstrowcol
		firstcol <- firstrowcol
	} else {
		firstrow <- firstrowcol[1:K,1:K,T:(2*T-1)]
		firstcol <- firstrowcol[1:K,1:K,T:1]
		# Fix dimensions of firstrow and firstcol (in case K=1)
		dim(firstrow) <- c(K,K,T)
		dim(firstcol) <- c(K,K,T)
	}
	

    # We pass through each block and (based on the first row):
    # Compute the eigenvalues for circulant embedding
    # No need for the circulant preconditioner, though.
    Aeigen<- rep(0, 2*K*T*K)
    dim(Aeigen) <- c(K,K,2*T)
   
    for(i in 1:K){
	for(j in 1:K){
	    Aeigen[i,j,] <- fft(circembedFromFirstRowAndColumn(firstrow[i,j,],firstcol[i,j,]),inverse=TRUE)
	}
    }
    
    return(multicg(Aeigen,Y,verbose=verbose))
}


### Function: MultivariatePCG
# This is the general wrapper function for the multivariate preconditioned
# conjugate gradient algorithm.  This function solves the system Ax=b for x, # when A is a block Toeplitz matrix (with K^2 blocks of size TxT).
# This takes in a K x K x T or K x K x (2T-1) matrix 
# containing the first row OR the first column and row of each block
# of the block Toeplitz matrix.  (If the first row is specified, then
# the matrix is assumed to be symmetric.  Otherwise, the first row and
# column must be specified, in this order:
# a(T,1), a(T-1,1), ..., a(1,1), a(1,2), ..., a(1,T)
# This corresponds to giving the covariances from lag -(T-1) to lag (T-1)
# in a time series context.)
# Inputs:
# firstrowcol - a KxKxT or KxKx(2T-1) array, in which the vector 
#			describing the first row (and possibly 
#			column) of the (i,j) block is in (i,j,:)
# Y - KT x 1 vector
# symmetric [FALSE] - indicator of whether the the Toeplitz matrix
#		      is being specified by only the first row, since
#		      it is symmetric
# verbose [FALSE] - should the number of iterations be printed?
# Output: X (T x 1)
MultivariatePCG.old <- function(firstrowcol, Y, symmetric=FALSE, verbose=FALSE){
    # Dealing with the case where K=1
    if(identical(dim(firstrowcol), NULL)) dim(firstrowcol) <- c(1,1,length(firstrowcol))

    # Compute dimensions
    K <- dim(firstrowcol)[1]
    T <- length(Y)/K

    # Dealing with the case where T=1
    if(length(dim(firstrowcol))==2) dim(firstrowcol) <- c(K,K,T)

    # Size checking
	if(dim(firstrowcol)[2] !=K){
		stop("Number of blocks per row and per column are not equal.")
	}
	if(symmetric & dim(firstrowcol)[3] !=T){
		stop("Size of vector and size of blocks incompatible")
	}
	if(!symmetric & dim(firstrowcol)[3] !=(2*T-1)){
		stop("Size of vector and size of blocks incompatible")
	}
	
	# Separate rows and columns
	if(symmetric){
		firstrow <- firstrowcol
		firstcol <- firstrowcol
	} else {
		firstrow <- firstrowcol[1:K,1:K,T:(2*T-1)]
		firstcol <- firstrowcol[1:K,1:K,T:1]
		# Fix dimensions of firstrow and firstcol (in case K=1)
		dim(firstrow) <- c(K,K,T)
		dim(firstcol) <- c(K,K,T)
	}
	

    # We pass through each block and (based on the first row):
    # Compute the eigenvalues for circulant embedding
    # Compute the circulant approximation and its eigenvalues
    Aeigen<- rep(0, 2*K*T*K)
    dim(Aeigen) <- c(K,K,2*T)
    lambda<- rep(0, K*T*K)
    dim(lambda) <- c(K,K,T)
   
    for(i in 1:K){
	for(j in 1:K){
	    Aeigen[i,j,] <- fft(circembedFromFirstRowAndColumn(firstrow[i,j,],firstcol[i,j,]),inverse=TRUE)
	    lambda[i,j,] <- fft(circapproxFromFirstRowAndColumn(firstrow[i,j,],firstcol[i,j,]),inverse=TRUE)
	}
    }
    
    ### Finish preconditioner
    # The blocks along the diagonal from Chan & Olkin
    # are the submatrices [1:K,1:K,t].
    # USING BUILT-IN INVERSION FOR THESE.
    lambdainv <- lambda
    for(t in 1:T){
	lambdainv[,,t] <- solve(lambda[,,t])
    }

    return(multipcg(Aeigen,lambdainv,Y,verbose))
}


### Function: MultivariatePCG
# This is the general wrapper function for the multivariate preconditioned
# conjugate gradient algorithm.  This function solves the system Ax=b for x, # when A is a block Toeplitz matrix (with K^2 blocks of size TxT).
# This takes in a K x K x T or K x K x (2T-1) matrix 
# containing the first row OR the first column and row of each block
# of the block Toeplitz matrix.  (If the first row is specified, then
# the matrix is assumed to be symmetric.  Otherwise, the first row and
# column must be specified, in this order:
# a(T,1), a(T-1,1), ..., a(1,1), a(1,2), ..., a(1,T)
# This corresponds to giving the covariances from lag -(T-1) to lag (T-1)
# in a time series context.)
# Inputs:
# firstrowcol - a KxKxT or KxKx(2T-1) array, in which the vector 
#			describing the first row (and possibly 
#			column) of the (i,j) block is in (i,j,:)
# Y - KT x 1 vector
# symmetric [FALSE] - indicator of whether the the Toeplitz matrix
#		      is being specified by only the first row, since
#		      it is symmetric
# verbose [FALSE] - should the number of iterations be printed?
# preconditioner [c("Jackson", "ChanOlkin")] - which preconditioner should be used?
# Output: X (T x 1)
# 5/28/08 (RS): This now offers the option of using either the Chan and Olkin preconditioner
#	or the Jackson kernel preconditioner of R. Chan, Yip, Ng
MultivariatePCG <- function(firstrowcol, Y, symmetric=FALSE, verbose=FALSE, preconditioner="ChanOlkin"){
    # Dealing with the case where K=1
    if(identical(dim(firstrowcol), NULL)) dim(firstrowcol) <- c(1,1,length(firstrowcol))

    # Compute dimensions
    K <- dim(firstrowcol)[1]
    T <- length(Y)/K

    # Dealing with the case where T=1
    if(length(dim(firstrowcol))==2) dim(firstrowcol) <- c(K,K,T)

    # Size checking
	if(dim(firstrowcol)[2] !=K){
		stop("Number of blocks per row and per column are not equal.")
	}
	if(symmetric & dim(firstrowcol)[3] !=T){
		stop("Size of vector and size of blocks incompatible")
	}
	if(!symmetric & dim(firstrowcol)[3] !=(2*T-1)){
		stop("Size of vector and size of blocks incompatible")
	}
	
	# Separate rows and columns
	if(symmetric){
		firstrow <- firstrowcol
		firstcol <- firstrowcol
		firstrowcol <- array(dim=c(K,K,2*T-1))
		firstrowcol[1:K,1:K,T:(2*T-1)] <- firstrow
		firstrowcol[1:K,1:K,T:1] <- firstcol
	} else {
		firstrow <- firstrowcol[1:K,1:K,T:(2*T-1)]
		firstcol <- firstrowcol[1:K,1:K,T:1]
		# Fix dimensions of firstrow and firstcol (in case K=1)
		dim(firstrow) <- c(K,K,T)
		dim(firstcol) <- c(K,K,T)
	}
	

    # We pass through each block and (based on the first row):
    # Compute the eigenvalues for circulant embedding
    Aeigen<- rep(0, 2*K*T*K)
    dim(Aeigen) <- c(K,K,2*T)
    for(i in 1:K){
	for(j in 1:K){
	    Aeigen[i,j,] <- fft(circembedFromFirstRowAndColumn(firstrow[i,j,],firstcol[i,j,]),inverse=TRUE)
	}
    }


    # Computing the eigenvalues of the preconditioner
    lambda<- rep(0, K*T*K)
    dim(lambda) <- c(K,K,T)
   
    if(preconditioner=="ChanOlkin"){
    	for(i in 1:K){
	    for(j in 1:K){
	    	lambda[i,j,] <- fft(circapproxFromFirstRowAndColumn(firstrow[i,j,],firstcol[i,j,]),inverse=TRUE)
	    }
    	}
	#print(lambda)
    } else if(preconditioner=="Jackson"){
	# This code borrows heavily from CH's code
	ntilde <- ceiling((T)/2)
	FejerWindow <- c((1:ntilde)/ntilde, ((ntilde-1):1)/ntilde)
	if(ntilde == 1){
	    FejerWindow <- c(1)
	} else {
	    FejerWindow <- c((1:ntilde)/ntilde, ((ntilde-1):1)/ntilde)
	}

	# This convolution is twice as long (minus 1), which allows for asymmetry in the block
	JacksonWindow <- convolution(FejerWindow,FejerWindow)
	JacksonWindow <- JacksonWindow/max(JacksonWindow)
	# In some cases, we must pad the beginning and the end
	if(ntilde==(T)/2) JacksonWindow <- c(0,JacksonWindow,0)

    	for(i in 1:K){
	    for(j in 1:K){
		dk <- firstrowcol[i,j,]*JacksonWindow;
		# To sum dk*exp(ik*theta) for k=-n+1,...,n-1, we have to FFT the two parts
		# of dk separately (first, k<=0, then k>=0)
		if(T==1){
		    lambda[i,j,] <- fft(dk[1],inverse=TRUE)
		} else {
	    	    lambda[i,j,] <- fft(dk[T:(2*T-1)],inverse=TRUE) + fft(c(0, dk[(T-1):1]), inverse=FALSE)
		}
	    }
    	}
	#print(lambda)
    } else {
	stop("Unknown preconditioner given")
    }
    
    ### Finish preconditioner
    # The blocks along the diagonal from Chan & Olkin
    # are the submatrices [1:K,1:K,t].
    # USING BUILT-IN INVERSION FOR THESE.
    lambdainv <- lambda
    for(t in 1:T){
	lambdainv[,,t] <- solve(lambda[,,t])
    }

    return(multipcg(Aeigen,lambdainv,Y,verbose))
}

### function: SignalPlusNoiseExtraction
# This function extracts the signal from a signal+noise model.  To do this,
# we use the formula from the paper, using MultiPCG to solve Omega^{-1}x.
# Inputs:
# autocov - the vector of autocovariances of the signal (length T)
# Tau - covariance matrix of noise (K x K)
# C - factor loadings (K x 1)
# Y - observations (KT x 1, observations grouped by series)
# alpha (NULL) - means of X (they will be subtracted off first)
# Output: X (T x 1)
SignalPlusNoiseExtraction <- function(autocov, Tau, C, Y,alpha=NULL){
    # Compute dimensions
    K <- dim(Tau)[1]
    T <- length(autocov)

    # Demean Y if necessary
    if(!identical(alpha, NULL)){
	for(k in 1:K){
	    Y[((k-1)*T+1):(k*T)] <- Y[((k-1)*T+1):(k*T)]-alpha[k]
	}
    }

    OmegainvX <- multipcgModel1(autocov, Tau, C, Y)

    # Multiply by c*Sigma
    h <- rep(0,T)
    for(i in 1:K){
	h <- h + toepmultFromFirstRow(C[i]*autocov,OmegainvX[((i-1)*T+1):(i*T)])
    }

    return(h)
}
    
### function: SignalPlusDiffNoiseExtraction
# This function extracts the signal from a signal+noise model in 
# which the variables have been differenced to make the signal 
# stationary (and consequiently made the noise non-invertible).  To do this,
# we use the formula from the paper, using MultiPCG to solve Omega^{-1}x.
# Inputs:
# autocov - the vector of autocovariances of the signal (length T)
# Tau - covariance matrix of (undifferenced) noise (K x K)
# C - factor loadings (K x 1)
# Y - observations (KT x 1, observations grouped by series, ALREADY DIFFERENCED)
# alpha (NULL) - means of X (they will be subtracted off first)
# Output: X (T x 1)
# 5/9/07 (RS): Added in the possibility of alpha.
SignalPlusDiffNoiseExtraction <- function(Sigma, Tau, C, Y,alpha=NULL){
    # Compute dimensions
    K <- dim(Tau)[1]
    T <- length(autocov)

    # Demean Y if necessary
    if(!identical(alpha, NULL)){
	for(k in 1:K){
	    Y[((k-1)*T+1):(k*T)] <- Y[((k-1)*T+1):(k*T)]-alpha[k]
	}
    }

    OmegainvX <- multipcgModel2(autocov, Tau, C, Y)

    # Multiply by c*Sigma
    h <- rep(0,T)
    for(i in 1:K){
	h <- h + toepmultFromFirstRow(C[i]*autocov,OmegainvX[((i-1)*T+1):(i*T)])
    }

    return(h)
}
 

### Function: QuadraticFormWithPCG
# This function computes the quadratic form, x'Vx using the multivariate
# preconditioned gradient algorithm, when V is a block Toeplitz matrix.
# Inputs: 
# firstrowcol - a KxKxT or KxKx(2T-1) matrix, in which the vector 
#				  describing the first row (and possibly 
#				  column) of the (i,j) block is in (i,j,:)
# x - KT x 1 vector
# symmetric [FALSE] - indicator of whether the the Toeplitz matrix
#					     is being specified by only the first row, since
#						  it is symmetric
# Output: x'Vx (1 x 1)
QuadraticFormWithPCG <- function(firstrowcol, x, symmetric=FALSE){
    firstproduct <- MultivariatePCG(firstrowcol, x, symmetric)
    secondproduct <- x %*% firstproduct

    # Return the first element, since result is a matrix
    return(secondproduct[1,1])
}