A couple of questions:
- have you looked at the code/etc. available under acf() and plot.acf()
in the ts library? Does it do what you want? (My guess is that it
doesn't, that it provides confidence intervals on the null ACF under two
different hypotheses, not confidence intervals on the estimates
themselves.) Computing the ACF itself is certainly faster using the
built-in acf(), but that may not be what's taking your time. Without
looking *too* closely at the code, I wonder if you can use (1) fft() or
(2) filter() to speed things up.
Ben Bolker
On Tue, 15 Jan 2002, Brian Scholl wrote:
> 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
>
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>
--
318 Carr Hall bolker at zoo.ufl.edu
Zoology Department, University of Florida http://www.zoo.ufl.edu/bolker
Box 118525 (ph) 352-392-5697
Gainesville, FL 32611-8525 (fax) 352-392-3704
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._