time series

[Paul Gilbert]

>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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._