[R] acf conf intervals +speed

[Brian Scholl]

Hi,

I'm trying to obtain confidence intervals for auto and
cross correlation estimates.  I've adapted code made
available by Stock and Watson that uses the Bartlett
Kernel and the delta method.  In R it runs really,
really slow because of the loops it uses and I have 9
series that I'd like to examine (81 total
combinations).  It was easy enough to replace one of
the while loops with a vector operation, but the
others are tough.  Anyone have either alternative code
or a suggestion for how to replace the other loops
below? 

Thanks



initc<-function(x,y){
n<-min(length(x),length(y))

# -- Demean -- @
xm<-x-mean(x) 
ym<-y-mean(y) 

# -- Cross Products -- @

a<-xm*ym 
b<-xm*xm 
c<-ym*ym 

cx0<-mean(b) 
cy0<-mean(c) 
list(cx0=cx0,cy0=cy0)
}




corse<-function(x,y,cx0,cy0,mix=1){
# FROM code by stock and watson from HOM paper
#corse.prc
 #  10/28/97, mww
#   Compute Correlation between x and y 
 #  Also compute se of estimated cor using
  # HAC-Delta Method Calculation 
  
  n<-min(length(x),length(y))
# -- Demean -- @
xm<-x-mean(x) 
ym<-y-mean(y) 

# -- Cross Products -- @
a<-xm*ym 
b<-xm*xm 
c<-ym*ym 

ma<-mean(a) 
mb<-mean(b) 
mc<-mean(c) 
cor<-ma/(sqrt(cx0*cy0))
if (mix==1){
	cor<-ma/(sqrt(mb*mc)) 
}

# -- Compute SE of Cor using VAR Spectral Est. and
delta method -- @
am<-a-mean(a) 
bm<-b-mean(b) 
cm<-c-mean(c) 
d<-cbind(ts(am),ts(bm),ts(cm)) 

# Form Bartlett Kernel @
kern<-matrix(0,n+1,1) 
ii <- 0  
seq<-(0:n)
kern<-1-(seq/(n+1))
#while (ii <= n){
# kern[ii+1]<-(1-ii/(n+1)) 
# ii<-ii+1  
#}
lam<-(t(d) %*% d)/nrow(d) 
vd<-kern[1]*lam

ii<-1
#while (ii <=n){
while (ii < (n-1)){
 lam<-(t(d[(1+ii):nrow(d),]) %*% d[1:(nrow(d)-ii),]) /
nrow(d) 
 vd<-vd+kern[ii+1]*(lam+t(lam)) 
 ii<-ii+1  
}
lam<-(t(d[(1+ii):nrow(d),]) %o% d[1:(nrow(d)-ii),]) /
nrow(d) 
lam<-lam[1,,]
vd<-vd+kern[ii+1]*(lam+t(lam)) 

vd<-vd/nrow(d) 

# Delta Method @
g<-matrix(0,3,1) 
g[1]<-1/(sqrt(mb*mc))    # Grad wrt a @
g[2]<-(-0.5)*(cor/mb)    # Grad wrt b @
g[3]<-(-0.5)*(cor/mc)    # Grad wrt c @

vcor<-t(g)%*%vd%*%g
secor<-sqrt(vcor) 

list(cor=cor,secor=secor)
}

acf2<-function(x,y=x,lagmax=(length(x)-10),m="ACF",mix=0){
	L<-min(length(x),length(y))
	cormat<-matrix(0,lagmax,1)
	secormat<-matrix(0,lagmax,1)
   
   x<-x[1:L]
	y<-y[1:L]
	c0<-initc(x,y)
	cx0<-c0$cx0
	cy0<-c0$cy0
	for (i in 1:lagmax){
  	 x1<-x[1:(L-i)]
  	 y1<-y[(i+1):L]
  	 c<-corse(x1,y1,cx0,cy0)
  	 cormat[i]<-c$cor
  	 secormat[i]<-c$secor
	}

	stop<-lagmax-sum(is.na(cormat))

plot(1:stop,cormat[1:stop],type="l",main=m,xlab="lag",ylab="ACF")

lines(cormat[1:stop]+2*secormat[1:stop],col="red",lty=2)

lines(cormat[1:stop]-2*secormat[1:stop],col="red",lty=2)
	lines(matrix(0,lagmax,1),col="black")
	cv<-1.96/sqrt(L)
	lines(matrix(cv,lagmax,1),col="blue",lty=2)

	list(cormat=cormat,secormat=secormat)
}




__________________________________________________


http://promo.yahoo.com/videomail/
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._