time series
Paul Gilbert
pgilbert@bank-banque-canada.ca
Tue, 16 Dec 1997 15:08:38 -0500
>I'm very interested in time series anlaysis (acf, pacf, arima modelling, and
>all other) is a possible enhancement for the WINDOWS 95 version of R!(?)
I am still hoping to make a Unix version of my DSE library available this week.
If I don't make it then it will be after the New Year as I'm away for two weeks.
I will make all the source code available, so the only real sense in which this
is a Unix version is that the installation procedures are set up for Unix.
The library does a lot of linear, time-invariant, vector ARIMA and state space
stuff, but does not do ar, acf, pacf or arima.mle as supplied in Splus. Below is
a rough start at ar and acf, but they are not yet working right. If anyone is
interested in finishing them I would be very happy. I haven't tried pacf and
I've never had any interest in arima.mle as it did not work for vector series.
The library should work (with the exception of a few non-critical functions) in
the existing Windows version of R, as I had most of it working in an older
version under Linux. I was going to wait until Robert has the Window version up
to 0.60, but the main problem is that I have little familiarity with Windows. So
if anyone wants to work on that, again, I would be more than happy.
There is a stand alone subset of the library (I call tframe) which is a kernel
for my other code. I'm fairly sure this work in the Windows version of R. This
provides a "classes and methods" way of handling the time dimension of an
object, which I think is very useful for making code that can fairly easily
handle new time representations. It also provides a replacement for tsmatrix and
some plotting routines which are useful for time series.
Paul Gilbert
___________
ar.test <-function(x, aic=T, order.max=2, method="yule-walker")
{# this seems to work, in the sense of producing estimates which are
# asymtotically as good as (or as bad as) Splus ar, at least for
# stationary data. But estimates are not the same.
# For non-stationary data Splus ar does better, but neither method is really
valid.
if(method=="burg") stop("burg method for ar not yet implemented.")
warning(" ar function not complete and not checked & default order.max=2.")
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 <- unclass(x)
tsp(x) <-NULL
x <- x - t(array(apply(x,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)
# # nrow above for univariate case
# Om <- Om0 %*% Om #Yule-Walker eqn. should solve without this
AC[i+1,,] <- Om
}
# now solve yule walker eqns.
n <- ncol(x)
a <- matrix(NA, n*(order.max), n*(order.max) )
b <- matrix(NA, n*(order.max), n )
# 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,,]
for (i in 1:order.max)
b[(1+(i-1)*n):(i*n),] <- AC[i+1,,]
AR <-solve(a, b)
if (n==1) AR <- matrix(AR, length(AR), 1)
#AR <- solve(t(a),b) # the off-diag values may require this??
# ar <- aperm(array(AR,c(dim(AC)-c(1,0,0))), c(2,1,3))
ar <- array(NA, c(order.max, ncol(x), ncol(x)))
for (i in 1:order.max)
ar[i,,] <- AR[(1+(i-1)*n):(i*n),]
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<- cov(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-devel 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-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._