Efficient computation of "variable (number of points included)" moving average in R

Question

Hi everybody,

I'm trying to implement a variable exponential moving average on a time series of intraday data (i.e 10 seconds). By variable, I mean that the size of the window included in the moving average depends on another factor (i.e. volatility). I was thinking of the following:

MA(t)=alpha(t)*price(t) + (1-alpha(t))MA(t-1),

where alpha corresponds for example to a changing volatility index.

In a backtest on huge series (more than 100000) points, this computation causes me "troubles". I have the complete vectors alpha and price, but for the current values of MA I always need the value just calculated before. Thus, so far I do not see a vectorized solution????

Another idea, I had, was trying to directly apply the implemented EMA(..,n=f()) function to every data point, by always having a different value for f(). But I do not find a fast solution neither so far.

Would be very kind if somebody could help me with my problem??? Even other suggestions of how constructing a variable moving average would be great.

Thx a lot in advance Martin

Answer 1

For timeseries, see the function rollmean in the zoo package.

You actually don't calculate a moving average, but some kind of a weighted cumulative average. A (weighted) moving average would be something like :

price <- runif(100,10,1000)
alpha <- rbeta(100,1,0.5)

tp <- embed(price,2)
ta <- embed(alpha,2)

MA1 <- apply(cbind(tp,ta),1,function(x){
    weighted.mean(x[1:2],w=2*x[3:4]/sum(x))
})

Make sure you rescale the weights so they sum to the amount of observations.

For your own calculation, you could try something like :

MAt <- price*alpha

ma.MAt <- matrix(rep(MAt,each=n),nrow=n)
ma.MAt[upper.tri(ma.MAt)] <- 0

tt1 <- sapply(1:n,function(x){
  tmp <- rev(c(rep(0,n-x),1,cumprod(rev(alpha[1:(x-1)])))[1:n])
  sum(ma.MAt[i,]*tmp)
})

This calculates the averages as linear combinations of MAt, with weights defined by the cumulative product of alpha.

On a sidenote : I assumed the index to lie somewhere between 0 and 1.

Answer 2

A very efficient moving average operation is also possible via filter():

  ## create a weight vector -- this one has equal weights, other schemes possible
  weights <- rep(1/nobs, nobs)     

  ## and apply it as a one-sided moving average calculations, see help(filter)
  movavg <- as.vector(filter(somevector, weights, method="convolution", side=1)) 

That was left-sided only, other choices are possible.

Answer 3

I just added a VMA function to the TTR package to do this. For example:

library(quantmod)  # loads TTR
getSymbols("SPY")
SPY$absCMO <- abs(CMO(Cl(SPY),20))/100
SPY$vma <- VMA(Cl(SPY), SPY$absCMO)
chartSeries(SPY,TA="addTA(SPY$vma,on=1,col='blue')")

x <- xts(rnorm(1e6),Sys.time()-1e6:1)
y <- xts(runif(1e6),Sys.time()-1e6:1)
system.time(VMA(x,y))  # < 0.5s on a 2.2Ghz Centrino

A couple notes from the documentation:

‘VMA’ calculate a variable-length moving average based on the absolute value of ‘w’. Higher (lower) values of ‘w’ will cause ‘VMA’ to react faster (slower).

The pre-compiled binaries should be on R-forge within 24 hours.