WMA（加重移動平均）計算を高速化する

Question

私は15日間のバーで指数移動平均を計算しようとしていますが、各日（終わり）のバーのEMAの「進化」を見たいと思っています。つまり、私は15日間のバーを持っているということです。新しいデータが毎日入ってくると、私は新しい情報を使ってEMAを再計算したいと思います。実際には私は15日間のバーを持っていて、毎日の後に新しい15日間のバーが成長し始め、新しいバーはEMA計算に使用される予定です。

2012-01-01（この例では各カレンダー日のデータがあります）から始まり、2012-01-15の最後に15日の最初の完全なバーがあります。 2012-03-01に4つの完全な15日間のバーが完成したら、4バーEMA（EMA（x、n = 4））の計算を開始できます。 2012-03-02の終わりに、この瞬間までの情報を使用し、2012-03-02のEMAを計算します。2012-03-02のOHLCが進行中の15日のバーです。そこで、2012-03-02に4つのバーとバーを取り、EMA（x、n = 4）を計算します。その後、別の日を待って、進行中の新しい15日間のバーで何が起こったかを見てください（詳細については、to.period.cumulativeを参照してください）。そしてEMAの新しい価値を計算してください...そして今後15日間は...ファンクションEMAの詳細は下記をご覧ください。

以下、私が今までに何を思いついたのかを見つけてください。パフォーマンスは私には受け入れられないし、限られたR知識ではそれをもっと速くすることはできない。私のシステムでは、それは

user system elapsed 
    4.708 0.000 4.410 

を取る

library(quantmod) 

do.call.rbind <- function(lst) { 
    while(length(lst) > 1) { 
     idxlst <- seq(from=1, to=length(lst), by=2) 

     lst <- lapply(idxlst, function(i) { 
        if(i==length(lst)) { return(lst[[i]]) } 

        return(rbind(lst[[i]], lst[[i+1]])) 
       }) 
    } 
    lst[[1]] 
} 

to.period.cumulative <- function(x, name=NULL, period="days", numPeriods=15) { 
    if(is.null(name)) 
     name <- deparse(substitute(x)) 

    cnames <- c("Open", "High", "Low", "Close") 
    if (has.Vo(x)) 
     cnames <- c(cnames, "Volume") 

    cnames <- paste(name, cnames, sep=".") 

    if (quantmod:::is.OHLCV(x)) { 
     x <- OHLCV(x) 
     out <- do.call.rbind( 
       lapply(split(x, f=period, k=numPeriods), 
         function(x) cbind(rep(first(x[,1]), NROW(x[,1])), 
           cummax(x[,2]), cummin(x[,3]), x[,4], cumsum(x[,5])))) 
    } else if (quantmod:::is.OHLC(x)) { 
     x <- OHLC(x) 
     out <- do.call.rbind( 
       lapply(split(x, f=period, k=numPeriods), 
         function(x) cbind(rep(first(x[,1]), NROW(x[,1])), 
           cummax(x[,2]), cummin(x[,3]), x[,4]))) 
    } else { 
     stop("Object does not have OHLC(V).") 
    } 

    colnames(out) <- cnames 

    return(out) 
} 

EMA.cumulative<-function(cumulativeBars, nEMA = 4, period="days", numPeriods=15) { 
    barsEndptCl <- Cl(cumulativeBars[endpoints(cumulativeBars, on=period,  k=numPeriods)]) 

    # TODO: This is sloooooooooooooooooow... 
    outEMA <- do.call.rbind(
      lapply(split(Cl(cumulativeBars), period), 
        function(x) { 
         previousFullBars <- barsEndptCl[index(barsEndptCl) < last(index(x)), ] 
         if (NROW(previousFullBars) >= (nEMA - 1)) { 
           last(EMA(last(rbind(previousFullBars, x), n=(nEMA + 1)), n=nEMA)) 
         } else { 
          xts(NA, order.by=index(x)) 
         } 
        })) 

    colnames(outEMA) <- paste("EMA", nEMA, sep="") 

    return(outEMA) 
} 

getSymbols("SPY", from="2010-01-01") 

SPY.cumulative <- to.period.cumulative(SPY, , name="SPY") 

system.time(
     SPY.EMA <- EMA.cumulative(SPY.cumulative) 
) 

許容実行時間が1秒未満になります...それは純粋なRを使用してこれを達成することは可能ですか？

この投稿はOptimize moving averages calculation - is it possible?にリンクしていますが、回答がありませんでした。私は今、スピードアップしたいものの詳細な説明と再現可能な例を作成することができました。私は今質問がより意味をなさないことを願っています。

これをスピードアップする方法に関するアイデアは高く評価されています。

Answer 1

私はRを使用して私の質問に満足のいく解決策を見つけることができませんでした。だから私は古いツール、C言語を取りました。驚くべきことに、Rcpp、インラインなどのこの素晴らしいツールを使用して私を「押し込む」ように感謝します。私は将来、パフォーマンス要件があり、Rを使用して満たすことができないときはいつでも、CにRを追加し、パフォーマンスがそこにあると思います。だから、私のコードとパフォーマンスの問題の解決策を見てください。

# How to speedup cumulative EMA calculation 
# 
############################################################################### 

library(quantmod) 
library(Rcpp) 
library(inline) 
library(rbenchmark) 

do.call.rbind <- function(lst) { 
    while(length(lst) > 1) { 
     idxlst <- seq(from=1, to=length(lst), by=2) 

     lst <- lapply(idxlst, function(i) { 
        if(i==length(lst)) { return(lst[[i]]) } 

        return(rbind(lst[[i]], lst[[i+1]])) 
       }) 
    } 
    lst[[1]] 
} 

to.period.cumulative <- function(x, name=NULL, period="days", numPeriods=15) { 
    if(is.null(name)) 
     name <- deparse(substitute(x)) 

    cnames <- c("Open", "High", "Low", "Close") 
    if (has.Vo(x)) 
     cnames <- c(cnames, "Volume") 

    cnames <- paste(name, cnames, sep=".") 

    if (quantmod:::is.OHLCV(x)) { 
     x <- quantmod:::OHLCV(x) 
     out <- do.call.rbind( 
       lapply(split(x, f=period, k=numPeriods), 
         function(x) cbind(rep(first(x[,1]), NROW(x[,1])), 
           cummax(x[,2]), cummin(x[,3]), x[,4], cumsum(x[,5])))) 
    } else if (quantmod:::is.OHLC(x)) { 
     x <- OHLC(x) 
     out <- do.call.rbind( 
       lapply(split(x, f=period, k=numPeriods), 
         function(x) cbind(rep(first(x[,1]), NROW(x[,1])), 
           cummax(x[,2]), cummin(x[,3]), x[,4]))) 
    } else { 
     stop("Object does not have OHLC(V).") 
    } 

    colnames(out) <- cnames 

    return(out) 
} 

EMA.cumulative<-function(cumulativeBars, nEMA = 4, period="days", numPeriods=15) { 
    barsEndptCl <- Cl(cumulativeBars[endpoints(cumulativeBars, on=period, k=numPeriods)]) 

    # TODO: This is sloooooooooooooooooow... 
    outEMA <- do.call.rbind(
      lapply(split(Cl(cumulativeBars), period), 
        function(x) { 
         previousFullBars <- barsEndptCl[index(barsEndptCl) < last(index(x)), ] 
         if (NROW(previousFullBars) >= (nEMA - 1)) { 
           last(EMA(last(rbind(previousFullBars, x), n=(nEMA + 1)), n=nEMA)) 
         } else { 
          xts(NA, order.by=index(x)) 
         } 
        })) 

    colnames(outEMA) <- paste("EMA", nEMA, sep="") 

    return(outEMA) 
} 

EMA.c.c.code <- ' 
    /* Initalize loop and PROTECT counters */ 
    int i, P=0; 

    /* ensure that cumbars and fullbarsrep is double */ 
    if(TYPEOF(cumbars) != REALSXP) { 
     PROTECT(cumbars = coerceVector(cumbars, REALSXP)); P++; 
    } 

    /* Pointers to function arguments */ 
    double *d_cumbars = REAL(cumbars); 
    int i_nper = asInteger(nperiod); 
    int i_n = asInteger(n); 
    double d_ratio = asReal(ratio); 

    /* Input object length */ 
    int nr = nrows(cumbars); 

    /* Initalize result R object */ 
    SEXP result; 
    PROTECT(result = allocVector(REALSXP,nr)); P++; 
    double *d_result = REAL(result); 

    /* Find first non-NA input value */ 
    int beg = i_n*i_nper - 1; 
    d_result[beg] = 0; 
    for(i = 0; i <= beg; i++) { 
     /* Account for leading NAs in input */ 
     if(ISNA(d_cumbars[i])) { 
      d_result[i] = NA_REAL; 
      beg++; 
      d_result[beg] = 0; 
      continue; 
     } 
     /* Set leading NAs in output */ 
     if(i < beg) { 
      d_result[i] = NA_REAL; 
     } 
     /* Raw mean to start EMA - but only on full bars*/ 
     if ((i != 0) && (i%i_nper == (i_nper - 1))) { 
      d_result[beg] += d_cumbars[i]/i_n; 
     } 
    } 

    /* Loop over non-NA input values */ 
    int i_lookback = 0; 
    for(i = beg+1; i < nr; i++) { 
     i_lookback = i%i_nper; 

     if (i_lookback == 0) { 
      i_lookback = 1; 
     } 
     /*Previous result should be based only on full bars*/ 
     d_result[i] = d_cumbars[i] * d_ratio + d_result[i-i_lookback] * (1-d_ratio); 
    } 

    /* UNPROTECT R objects and return result */ 
    UNPROTECT(P); 
    return(result); 
' 

EMA.c.c <- cfunction(signature(cumbars="numeric", nperiod="numeric", n="numeric",  ratio="numeric"), EMA.c.c.code) 

EMA.cumulative.c<-function(cumulativeBars, nEMA = 4, period="days", numPeriods=15) { 
    ratio <- 2/(nEMA+1) 

    outEMA <- EMA.c.c(cumbars=Cl(cumulativeBars), nperiod=numPeriods, n=nEMA, ratio=ratio) 

    outEMA <- reclass(outEMA, Cl(cumulativeBars)) 

    colnames(outEMA) <- paste("EMA", nEMA, sep="") 

    return(outEMA) 
} 

getSymbols("SPY", from="2010-01-01") 

SPY.cumulative <- to.period.cumulative(SPY, name="SPY") 

system.time(
     SPY.EMA <- EMA.cumulative(SPY.cumulative) 
) 

system.time(
     SPY.EMA.c <- EMA.cumulative.c(SPY.cumulative) 
) 


res <- benchmark(EMA.cumulative(SPY.cumulative), EMA.cumulative.c(SPY.cumulative), 
     columns=c("test", "replications", "elapsed", "relative", "user.self", "sys.self"), 
     order="relative", 
     replications=10) 

print(res) 

EDIT：

> print(res) 
           test replications elapsed relative user.self 
2 EMA.cumulative.c(SPY.cumulative)   10 0.026 1.000  0.024 
1 EMA.cumulative(SPY.cumulative)   10 57.732 2220.462 56.755 

：私の面倒を超える性能向上の指示を与えるためにRは、ここにプリントアウトである（私は実際に私がループのダブル作成しておりますので、それは、良好なものとすることができると確信しています）

私の基準では、改善のSFタイプ...