From: Brian Scholl (brianscholl@yahoo.com)
Date: Wed 16 Jan 2002 - 08:25:02 EST
Message-id: <20020115212502.36233.qmail@web20308.mail.yahoo.com>
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@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
This archive was generated by hypermail 2.1.3 : Thu 17 Jan 2002 - 11:10:11 EST