GrangerCoh=function(x, y, filter=T, frequency.range=F,plot=T)
{
# Computes the Pierce version of the Granger Coefficient of Coherency (GrangerCoh)
# as discussed in Lemmens, A., Croux, C., and Dekimpe, M.G. (2008)
#  ``Measuring and Testing Granger Causality over the Spectrum: 
#   An application to European Production Expectation Surveys,'' International Journal of Forecasting.
#
# Input: 
# x and y:  two univariate time series of the same length. 
#           They are supposed to be filtered.
# filter: if filter==F, the series are not filtered yet and a filtering is done within this program with an AR-model
# frequency range: vector of frequencies for which the GrangerCoh needs to be computed
# plot: if plot==T, a plot of the Granger Coefficient of Coherency over the spectrum,
#       together with the critical value for a test of the zero Grancer causality, is provided
#
# 
# Output:
# omega:  vector of frequencies for which the GrangerCoh is computed
# GCcoefcoh: vector of values for the  GrangerCoh 
# cv: critical value for the test of zero GrangerCoh
#
#Example 1: x=arima.sim(n=100,list(ar=c(0.5,0.3)));
#           y=arima.sim(n=100,list(ar=c(0.2,0.3,0.2,0.1)));
#           GrangerCoh(x, y, filter=T);
#Example 2: x=rnorm(100);y=rnorm(100);
#           GrangerCoh(x, y)


Tlength=length(x);
M=floor(sqrt(Tlength));

if (filter==F)
{
xfilter=ar(x,order.max=M)
yfilter=ar(y,order.max=M)
nmissing=max(xfilter$order,yfilter$order)
x=as.vector(xfilter$res);y=as.vector(yfilter$res)
x=x[(nmissing+1):(Tlength)];
y=y[(nmissing+1):(Tlength)];
}

Tlength=length(x);
M=floor(sqrt(Tlength));

gamma=acf(cbind(y,x),plot=F,type="covariance");
lag=gamma$lag;

gammax=gamma$acf[,2,2]
allgammax=c(rev(gammax[-1]),gammax);
gammay=gamma$acf[,1,1]
allgammay=c(rev(gammay[-1]),gammay);

posgamma=gamma$acf[-1,1,2];
nulgamma=gamma$acf[1,1,2];
neggamma=gamma$acf[-1,2,1];
allgamma=c(rev(neggamma),nulgamma,posgamma);
        
k=c(rev(lag[-1,2,1]),lag[,1,1]);
posk=lag[-1,1,1];
lambda=pmax(1-abs(k)/M,0);
poslambda=pmax(1-abs(posk)/M,0);

omega=frequency.range
if (frequency.range==F) omega=seq(0,pi,length=64)
iunit=complex(re=0,im=1);

outermatrix=outer(omega,k, FUN="*");
result2=exp(-outermatrix*iunit);
crossspectrum=(result2%*%as.vector(lambda*allgamma))/(2*pi);

posoutermatrix=outer(omega,posk, FUN="*");
posresult2=exp(-posoutermatrix*iunit);
poscrossspectrum=(posresult2%*%as.vector(poslambda*posgamma))/(2*pi);
        

spectrumx=Re((result2%*%as.vector(lambda*allgammax))/(2*pi));
spectrumy=Re((result2%*%as.vector(lambda*allgammay))/(2*pi));

GCcoefcoh=sqrt((Re(poscrossspectrum)^2+Im(poscrossspectrum)^2)/(spectrumx*spectrumy));
nprime=Tlength/sum(poslambda^2);   
cv=sqrt(qchisq(0.95,2)/(2*nprime-2));

if (plot==T)
{
ymax=max(c(max(max(GCcoefcoh),cv)*1.1,1));
plot(omega,GCcoefcoh,xlab="frequency",ylab="Granger  Coherence",
     main="Granger Coefficient of Coherence",type='l',ylim=c(0,ymax))
abline(h=cv);
}

list(omega=omega,GCcoefcoh=GCcoefcoh,cv=cv);   
}