# GARCH(1,1) estimation objf <- function(x, eps){ n <- length(eps) omega <- x[1] alpha <- x[2] beta <- x[3] if((alpha+beta>0.9999)|(alpha+beta=='NaN'))return(Inf) sigma2<-rep(0,n) sigma2[1]<-omega for (t in 2:n) sigma2[t]<-omega+alpha*eps[t-1]**2+beta*sigma2[t-1] qml <- mean(log(sigma2[2:n])+eps[2:n]**2/sigma2[2:n]) qml } estimgarch11<- function(omega0,alpha0,beta0,eps0) { n<-length(eps0) norm<-sd(eps0) eps<-(eps0-mean(eps0))/norm valinit<-c(omega0,alpha0,beta0) res <- nlminb(valinit,objf, lower=c(0.0000001,0,0), upper=c(Inf,0.9999,0.9999), eps=eps) omega<-res$par[1]*norm**2 alpha<-res$par[2] beta<-res$par[3] list(coef=c(omega,alpha,beta), minimum=res$objective) } # autocorrelation function 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) ################################ # 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("cac.csv","dax.csv","ftse.csv","nikkei.csv","smi.csv","sp500.csv") name.indices<-c("CAC","DAX","FTSE","Nikkei","SMI","SP500") 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=",") close<-rev(data$Close) date<-rev(data$Date) ntot<-length(close) rend<-rep(0,ntot); rend[2:ntot]<-log(close[2:ntot]/close[1:(ntot-1)])*100 ser<-rend[2:ntot]-mean(rend[2:ntot]) if(name.indices[ind]=="Nasdaq")ser<-ser[-which.min(ser)] nl.acf(ser,main=name.indices[ind],method='NP') } par(op) 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=",") close<-rev(data$Close) date<-rev(data$Date) ntot<-length(close) rend<-rep(0,ntot); rend[2:ntot]<-log(close[2:ntot]/close[1:(ntot-1)])*100 ser<-rend[2:ntot]-mean(rend[2:ntot]) if(name.indices[ind]=="Nasdaq")ser<-ser[-which.min(ser)] nl.acf(ser,main=name.indices[ind],method='PAR') } par(op)