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