R-beta: ar

[Paul Gilbert]

I have been trying to get a working version of ar, since I have used it
in several calculations in the test suite for my time series library.
The following limited version (order.max must be specified and other
short comings) works more or less, but the results differ by more than
I would expect from those given by Splus. I have tried several
variations with no success. If anyone can see a reason I would much
appreciate help.

Further below is a limited version of acf, which reproduces Splus
results a little better.

Paul Gilbert

____________

ar.test <-function(x, aic=T, order.max=NULL, method="yule-walker")
   {if(method=="burg") stop("burg method for ar not yet implemented.")
    warning(" ar function not complete and not checked.")   
    if (is.vector(x))x <- matrix(x, length(x),1)
    sampleT <- nrow(x)
    if (is.null(order.max)) order.max <- round(10*log10(sampleT/ncol(x)))
    AC <- array(NA, c(order.max+1, ncol(x), ncol(x)))
    # x <- x - t(array(apply(z,2,sum)/nrow(x), rev(dim(x))))
    for (i in 0:order.max)
	  {#Om<- (t(x[(i+1):sampleT,,drop=F]) %*% 
           #        x[1:(sampleT-i),,drop=F])/(sampleT-i)
	   Om<- cor(x[(i+1):sampleT,,drop=F],
                    x[1:(sampleT-i),,drop=F]) # cor seems better than var
           if(i==0) Om0 <- solve(Om)
           Om <- Om0 %*% Om
           AC[i+1,,] <-  Om
          }

   # now solve yule walker eqns.
   n <- ncol(x)
   a <- matrix(NA, n*(order.max), n*(order.max) )
   # using AC[1,,] in place of I
   # there must be a better way (with outer?)
    for (i in 0:(order.max-1))
       for (j in 0:(order.max-1))
	  a[(1+i*n):(n+i*n),(1+j*n):(n+j*n)] <- AC[abs(i-j)+1,,] 
   AR <-solve(a, matrix(aperm(AC[2:(order.max+1),,,drop=F],c(2,1,3)), n*order.max,n))
   AR <- aperm(array(AR,c(dim(AC)-c(1,0,0))), c(2,1,3)) 
   order <- order.max # need aic here
   list(ar=AR, order=order, order.max=order.max, method=method)
   }


acf <-function (residual, plot = F, type = "correlation") 
       {if (plot) warning(" acf plot not yet supported.")
        if(0==charmatch(type,c("covariance","correlation","partial"),nomatch=0))
             stop("type not allowed in acf")
        if (is.vector(residual))residual <- matrix(residual, length(residual),1)
	sampleT <- nrow(residual)
        N <- round(10*(log10(sampleT)-log10(ncol(residual)))) 
        acf <- array(NA, c(N, ncol(residual), ncol(residual)))
	for (i in 0:(N-1))
	  {Om<- var(residual[(i+1):sampleT,,drop=F], 
                    residual[1:(sampleT-i),,drop=F])
           if (type=="correlation")
                {if(i==0) Om0 <- diag(1/sqrt(diag(Om)),nrow=nrow(Om))
                            # nrow above for univariate case
                 Om <- Om0 %*% Om %*% Om0
                }
           acf[i+1,,] <-  Om 
          }
        if (type=="partial") 
          {warning("acf type partial not yet supported. 0 value being returned")
           acf <- array(0, dim(acf))
          }
	list(acf=acf, type=type )
       }



=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
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
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=