On Fri, Nov 4, 2011 at 9:24 AM, Gabor Grothendieck
<ggrothendieck@gmail.com>wrote:

> On Fri, Nov 4, 2011 at 9:09 AM, Matthew Clegg <matthewcleggphd@gmail.com>
> wrote:
> > Hello R-Sig-Finance members:
> >
> > I was wondering if anyone has contributed functions that are similar
> > to the zoo roll* functions but which operate on fixed-length time
> > windows?  For example, suppose I have a zoo-based object consisting
> > of the daily closing prices of a stock, and I wish to know for each
> > date, what was the volatility over the succeeding 30 calendar days?
> > Probably many people would settle for something like:
> >  rollapply (log(lag(P))-log(P), 21, sd, align="left") * sqrt(252)
> > (where P is the price series).  However, this is an approximation.
> > Not all periods of 30 calendar days include precisely 21 trading days.
> >
> > This seems like an obvious enough question that I would think that it
> > has been asked (and answered) many times before, but I could not find
> > a reference to the recommended solution.
> >
> > If no one has tackled this problem before, I might try to put together
> > a small library of functions that are like roll* but which operate
> > on fixed time windows.  I am including an example of one such function
> > below.
> >
> > Matthew Clegg
> >
> > [snip]
> >
>
> Here is a one liner (two if you count making the result into a zoo object):
>
> > z <- zoo(1:25)
> > zz <- sapply(seq_along(z), function(i) sum(z[time(z) <= time(z)[i] &
> time(z) > time(z)[i] - 3]))
> > zoo(zz, time(z))
>  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
>  1  3  6  9 12 15 18 21 24 27 30 33 36 39 42 45 48 51 54 57 60 63 66 69 72
>
> --
> Statistics & Software Consulting
> GKX Group, GKX Associates Inc.
> tel: 1-877-GKX-GROUP
> email: ggrothendieck at gmail.com
>


Aha!  That's an elegant solution and another great illustration of the
power of vector processing in R.

I found that after tweaking my code, I could achieve a significant
improvement in running time over this sapply()-based one liner.  The
following table compares the running times for various lengths of the
underlying zoo vector:

> compare_times()
      sapply ztw_sum
100    0.058   0.003
200    0.108   0.004
400    0.222   0.008
600    0.354   0.011
800    0.499   0.016
1000   0.651   0.018
2000   1.713   0.036
4000   5.023   0.073
8000  18.143   0.148
16000 64.770   0.290

These times were obtained on a 2.66GHz Intel Mac running Mac OS 10.7.2
and R 2.14.0.  Similar results were obtained in the Amazon cloud,
using a 32-bit Small Linux instance.

One thing that I did not properly appreciate was the high cost of
indexing into zoo vectors.  In the original code that I presented,
the main loop contained a number of expressions that involved
indexing into a zoo vector.  For example, one such line read
    sum <- sum + as.numeric(X[j]);
Before the beginning of the loop, I added:
    Xc <- coredata(X)
Then, I changed each zoo access (X[i]) to an access of Xc[i], e.g.,
    sum <- sum + Xc[j];
The original code ran slower than the sapply() one liner, but making
this one change provided a dramatic speedup.

Source code follows ...

compare_times <- function (W=3, N=c(100, 200, 400, 600, 800, 1000, 2000,
4000, 8000, 16000),
                           index.base=as.Date("1970-01-01")) {
  funs <- c(sapply = function(z) { zoo(sapply(seq_along(z), function(i)
sum(z[time(z) <= time(z)[i] & time(z) >= time(z)[i] - W])), order.by=index(z))
},
#            rollapply = function(z) { rollapply(z, W+1, sum, partial=TRUE,
align="right") },
            ztw_sum = function(z) { ztw_sum(z, W, partial=TRUE) }
            )

  ftime <- function (f, N) { sapply(N, function(k) { system.time(f(zoo(1:k,
order.by=index.base+1:k)))["user.self"] }) }
  E <- sapply (funs, function(F) { ftime(F, N)})
  rownames(E) <- N;

  E;
}

ztw_sum <- function (X, delta, align="right", partial=FALSE) {
  # Zoo Time Window Sum
  #
  # On input, X is a zoo-based numeric vector and delta is a time
difference.
  # Constructs a zoo-based numeric vector of partial sums from X.  The
values
  # included in a partial sum are those whose associated timestamps are
  # within delta of the corresponding element from X.
  #
  # If align="right", then result[i] is a sum of those elements
  # X[j] such that
  #    0 <= timestamp[i] - timestamp[j] <= delta,
  # where timestamp[i] is the timestamp (index) associated with the
  # i-th element of X.  Conversely, if align="left", then result[i] is a
  # sum of those elements X[j] such that
  #    0 <= timestamp[j] - timestamp[i] <= delta.
  #
  # Parameters:
  # X:        A zoo-based numeric vector with a time-based index type.
  # delta:    An object of type difftime specifying the size of
  #           the time window.
  # align:    Specifies whether the sum for a given index should
  #           be computed using elements of lower timestamps ("right")
  #           or higher timestamps ("left").
  # partial:  If TRUE, then partial sums are computed for elements
  #           at the left (respectively, right) end of the vector.
  #
  # Returns a zoo-based numeric vector of partial sums.

  if (!inherits(X, "zoo")) {
    stop ("X must be a vector of type zoo");
  }

  Xc <- coredata(X)
  if (!inherits(Xc, "numeric") && !inherits(Xc, "integer")) {
    stop ("coredata(X) must be a numeric or integer vector");
  } else if (delta <= 0) {
    stop ("delta must be positive");
  } else if ((align != "left") && (align != "right")) {
    stop ("align must be from c('left', 'right')");
  }

  timestamp <- index(X)
  timestampd <- timestamp - delta

  if (class(timestamp)[1] %in% c("Date", "POSIXlt", "POSIXct")) {
    # This can achieve a dramatic speedup.  It seems that doing
    # comparisons of the internal date and time types is quite slow.
    timestamp <- as.numeric(timestamp);
    timestampd <- as.numeric(timestampd);
  }

  N <- length(X)
  R <- rep(NA, N)
  sum <- 0;  # The current partial sum

  if (align == "right") {
    # Invariants:
    #   (a) 0 < i <= j <= N
    #   (b) 0 <= timestamp(j) - timestamp(i) <= delta
    i <- 1;  # The left edge of the current window
    for (j in 1:N) {
      if (!is.na(Xc[j])) {
        sum <- sum + Xc[j];
      }
      while (timestampd[j] > timestamp[i]) {
        if (!is.na(Xc[i])) {
          sum <- sum - Xc[i];
        }
        i <- i+1;
      }
      if ((i > 1) || partial || (timestampd[j] == timestamp[i])) {
        R[j] <- sum;
      }
    }
  } else { # align == "left"
    # Invariants:
    #   (a) 0 < j <= i <= N
    #   (b) 0 <= timestamp(i) - timestamp(j) <= delta
    i <- N;  # The right edge of the current window
    for (j in N:1) {
      if (!is.na(Xc[j])) {
        sum <- sum + Xc[j];
      }
      while (timestampd[i] > timestamp[j]) {
        if (!is.na(Xc[i])) {
          sum <- sum - Xc[i];
        }
        i <- i-1;
      }
      if ((i < N) || partial || (timestamp[j] == timestampd[i])) {
        R[j] <- sum;
      }
    }
  }

  class(R) <- class(Xc)
  zoo(R, order.by=index(X))
}


-- 
Matthew Clegg
matthewcleggphd@gmail.com

	[[alternative HTML version deleted]]

