# fonction d'autocorrelation
gamma<-function(x,h){ n<-length(x); h<-abs(h);x<-x-mean(x)
gamma<-sum(x[1:(n-h)]*x[(h+1):n])/n }
rho<-function(x,h) rho<-gamma(x,h)/gamma(x,0)
# fonction acf avec bandes de significativité d'un bruit blanc semi-fort
################################
# acf with significance limits for uncorrelated but nonindependent white noises 
# NP is the method in which the generalized Bartlett formula is  estimated nonparametrically (Figure 2)
# PAR is the method in which the generalized Bartlett formula assuming a GARCH(1,1) model (Figure 3)
################################

nl.acf<-function(x,main=NULL,method='NP',xlab="Lag",ylab=NULL){
n<-length(x)
nlag<-as.integer(min(10*log10(n),n-1))
acf.val<-sapply(c(1:nlag),function(h) rho(x,h))

if(method=='NP'){
x2<-x^2
var<-1+(sapply(c(1:nlag),function(h) gamma(x2,h)))/gamma(x,0)^2
band<-sqrt(var/n)
minval<-1.2*min(acf.val,-1.96*band,-1.96/sqrt(n))
maxval<-1.2*max(acf.val,1.96*band,1.96/sqrt(n))
acf(x,ylim=c(minval,maxval),main=main,xlab=xlab,ylab=ylab)
lines(c(1:nlag),-1.96*band,lty=1,col='red')
lines(c(1:nlag),1.96*band,lty=1,col='red')}

if(method=='PAR'){
alpha0<-0.09
beta0<-0.89
omega0<-var(x)*(1-alpha0-beta0)
est<-estimgarch11(omega0,alpha0,beta0,x)
omega<-est$coef[1]
alpha<-est$coef[2]
beta<-est$coef[3]
den<-1-3*alpha^2-beta^2-2*alpha*beta
if(den<0)return(NA)
E.nu.2<-2*omega^2*(1+alpha+beta)/((1-alpha-beta)*den)
rho2<-NULL
rho2[1]<-alpha*(1-beta^2-alpha*beta)/(1-beta^2-2*alpha*beta)
for(h in (2:nlag))rho2[h]<-(alpha+beta)*rho2[h-1]
gam2<-NULL
gam2[1]<-E.nu.2*(1+alpha^2/(1-(alpha+beta)^2))
for(h in (2:(nlag+1)))gam2[h]<-rho2[h-1]*gam2[1]
gam.0<-omega/(1-alpha-beta)
var<-1+gam2[2:(nlag+1)]/gam.0^2
band<-sqrt(var/n)
minval<-1.2*min(acf.val,-1.96*band,-1.96/sqrt(n))
maxval<-1.2*max(acf.val,1.96*band,1.96/sqrt(n))
acf(x,ylim=c(minval,maxval),main=main,xlab=xlab,ylab=ylab)
lines(c(1:nlag),-1.96*band,lty=1,col='red')
lines(c(1:nlag),1.96*band,lty=1,col='red')}
}

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

indices<-c("yen.csv", "dollar.csv","Livre.csv","SWF.csv")
name.indices<-c("Yen","Dollar","Livre","SWF")
nbind<-length(indices)

op <- par(mfrow = c(3, 2),cex.main=0.8) # 2 x 2 pictures on one plot


for(ind in (1:nbind)){
data <- read.table(indices[ind],header=TRUE,sep=";")
taux<-data$ser[!is.na(data$ser)]
ntot<-length(taux)
rend<-rep(0,ntot); rend[2:ntot]<-log(taux[2:ntot]/taux[1:(ntot-1)])*100 
ser<-rend[2:ntot]-mean(rend[2:ntot])
nl.acf(ser,main=name.indices[ind])
}
par(op)