# Fast Robust Moments

Fast, numerically robust computation of weighted moments via 'Rcpp'. Supports computation on vectors and matrices, and Monoidal append of moments. Moments and cumulants over running fixed length windows can be computed, as well as over time-based windows. Moment computations are via a generalization of Welford's method, as described by Bennett et. (2009) .

Fast Robust Moments -- Pick Three!

Fast, numerically robust, higher order moments in R, computed via Rcpp, mostly as an exercise to learn Rcpp. Supports computation on vectors and matrices, and Monoidal append (and unappend) of moments. Computations are via the Welford-Terriberry algorithm, as described by Bennett et al.

-- Steven E. Pav, [email protected]

## Installation

This package can be installed from CRAN, via drat, or from github:

# Basic Usage

Currently the package functionality can be divided into the following:

• Functions which reduce a vector to an array of moments.
• Functions which take a vector to a matrix of the running moments.
• Functions which transform a vector to some normalized form, like a centered, rescaled, z-scored sample, or a summarized form, like the running Sharpe or t-stat.
• Functions for computing the covariance of a vector robustly.
• Object representations of moments with join and unjoin methods.

## Summary moments

A function which computes, say, the kurtosis, typically also computes the mean and standard deviation, and has performed enough computation to easily return the skew. However, the default functions in R for higher order moments discard these lower order moments. So, for example, if you wish to compute Merten's form for the standard error of the Sharpe ratio, you have to call separate functions to compute the kurtosis, skew, standard deviation, and mean.

The summary functions in fromo return all the moments up to some order, namely the functions `sd3`, `skew4`, and `kurt5`. The latter of these, `kurt5` returns an array of length 5 containing the excess kurtosis, the skewness, the standard deviation, the mean, and the observation count. (The number in the function name denotes the length of the output.) Along the same lines, there are summarizing functions that compute centered moments, standardized moments, and 'raw' cumulants:

• `cent_moments`: return a `k+1`-vector of the `k`th centered moment, the `k-1`th, all the way down to the 2nd (the variance), then the mean and the observation count.
• `std_moments`: return a `k+1`-vector of the `k`th standardized moment, the `k-1`th, all the way down to the 3rd, then the standard deviation, the mean, and the observation count.
• `cent_cumulants`: computes the centered cumulants (yes, this is redundant, but they are not standardized). return a `k+1`-vector of the `k`th raw cumulant, the `k-1`th, all the way down to the second, then the mean, and the observation count.
• `std_cumulants`: computes the standardized (and, of course, centered) cumulants. return a `k+1`-vector of the `k`th standardized cumulant, all the way down to the third, then the variance, the mean, and the observation count.
``````## [1]   47.276   -0.047    3.986   10.092 1000.000
``````
``````## [1]  3.0e+00 -5.9e-03  2.0e+00  1.0e+01  1.0e+03
``````
``````## [1]   -0.388   -0.047    3.986   10.092 1000.000
``````
``````## [1] -2.4e-02 -5.9e-03  4.0e+00  1.0e+01  1.0e+03
``````

### Speed

In theory these operations should be just as fast as the default functions, but faster than calling multiple default functions. Here is a speed comparison of the basic moment computations:

``````## Error in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels) else paste0(labels, : factor level [5] is duplicated
``````
``````## Unit: milliseconds
##         expr  min   lq mean median   uq  max neval    cld
##     kurt5(x) 1408 1409 1414   1413 1418 1424    10     e
##     skew4(x)  801  802  813    804  815  871    10    d
##       sd3(x)   66   66   67     66   67   69    10  b
##     dumbk(x) 1647 1650 1674   1655 1712 1734    10      f
##  kurtosis(x)  802  804  813    806  811  869    10    d
##  skewness(x)  777  780  789    781  787  822    10   c
##        sd(x)   45   45   46     46   46   47    10  b
##      mean(x)   16   17   17     17   17   18    10 a
``````

## Weight! Weight!

Many of the methods now support the computation of weighted moments. There are a few options around weights: whether to check them for negative values, whether to normalize them to unit mean.

``````## [1] 2.9e+00 1.2e-02 1.0e+00 1.0e-02 1.0e+03
``````
``````## [1] 3.1e+00 4.1e-02 1.0e+00 1.3e-02 1.0e+03
``````
``````## [1]   3.072   0.041   1.001   0.013 493.941
``````
``````## Unit: milliseconds
##             expr min  lq mean median  uq max neval cld
##  sd3(x, wts = w)  80  81   83     83  84  87   100  a
##    slow_sd(x, w) 234 243  253    247 254 382   100   b
``````

## Monoid mumbo-jumbo

The `as.centsums` object performs the summary (centralized) moment computation, and stores the centralized sums. There is a print method that shows raw, centralized, and standardized moments of the ingested data. This object supports concatenation and unconcatenation. These should satisfy 'monoidal homomorphism', meaning that concatenation and taking moments commute with each other. So if you have two vectors, `x1` and `x2`, the following should be equal: `c(as.centsums(x1,4),as.centsums(x2,4))` and `as.centsums(c(x1,x2),4)`. Moreover, the following should also be equal: `as.centsums(c(x1,x2),4) %-% as.centsums(x2,4))` and `as.centsums(x1,4)`. This is a small step of the way towards fast machine learning methods (along the lines of Mike Izbicki's Hlearn library).

Some demo code:

``````##           class: centsums
##     raw moments: 100 0.0051 0.09 -0.00092 0.014 -0.00043 0.0027
## central moments: 0 0.09 -0.0023 0.014 -0.00079 0.0027
##     std moments: 0 1 -0.086 1.8 -0.33 3.8
``````

We also have 'raw' join and unjoin methods, not nicely wrapped:

### For multivariate input

There is also code for computing co-sums and co-moments, though as of this writing only up to order 2. Some demo code for the monoidal stuff here:

``````## An object of class "centcosums"
## Slot "cosums":
##          [,1]    [,2]   [,3]     [,4]    [,5]
## [1,] 100.0000  -0.093  0.045  -0.0046   0.046
## [2,]  -0.0934 111.012  4.941 -16.4822   6.660
## [3,]   0.0450   4.941 71.230   0.8505   5.501
## [4,]  -0.0046 -16.482  0.850 117.3456  13.738
## [5,]   0.0463   6.660  5.501  13.7379 100.781
##
## Slot "order":
## [1] 2
``````

## Running moments

Since an online algorithm is used, we can compute cumulative running moments. Moreover, we can remove observations, and thus compute moments over a fixed length lookback window. The code checks for negative even moments caused by roundoff, and restarts the computation to correct; periodic recomputation can be forced by an input parameter.

A demonstration:

``````##       excess_kurtosis  skew stdev   mean nobs
##  [1,]             NaN   NaN   NaN -1.207    1
##  [2,]             NaN   NaN  1.05 -0.465    2
##  [3,]             NaN -0.34  1.16  0.052    3
##  [4,]          -1.520 -0.13  1.53 -0.548    4
##  [5,]          -1.254 -0.50  1.39 -0.352    5
##  [6,]          -0.860 -0.79  1.30 -0.209    6
##  [7,]          -0.714 -0.70  1.19 -0.261    7
##  [8,]          -0.525 -0.64  1.11 -0.297    8
##  [9,]          -0.331 -0.58  1.04 -0.327    9
## [10,]          -0.331 -0.42  1.00 -0.383   10
## [11,]           0.262 -0.65  0.95 -0.310   10
## [12,]           0.017 -0.30  0.95 -0.438   10
## [13,]           0.699 -0.61  0.79 -0.624   10
## [14,]          -0.939  0.69  0.53 -0.383   10
## [15,]          -0.296  0.99  0.64 -0.330   10
## [16,]           1.078  1.33  0.57 -0.391   10
## [17,]           1.069  1.32  0.57 -0.385   10
## [18,]           0.868  1.29  0.60 -0.421   10
## [19,]           0.799  1.31  0.61 -0.449   10
## [20,]           1.193  1.50  1.07 -0.118   10
``````
``````##         alt5
##  [1,]    NaN    NaN
##  [2,] -2.000    NaN
##  [3,] -1.500    NaN
##  [4,] -1.520 -1.520
##  [5,] -1.254 -1.254
##  [6,] -0.860 -0.860
##  [7,] -0.714 -0.714
##  [8,] -0.525 -0.525
##  [9,] -0.331 -0.331
## [10,] -0.331 -0.331
## [11,]  0.262  0.262
## [12,]  0.017  0.017
## [13,]  0.699  0.699
## [14,] -0.939 -0.939
## [15,] -0.296 -0.296
## [16,]  1.078  1.078
## [17,]  1.069  1.069
## [18,]  0.868  0.868
## [19,]  0.799  0.799
## [20,]  1.193  1.193
``````

If you like rolling computations, do also check out the following packages (I believe they are all on CRAN):

Of these three, it seems that `RollingWindow` implements the optimal algorithm of reusing computations, while the other two packages gain efficiency from parallelization and implementation in C++.

## Running 'scale' operations

Through template magic, the same code was modified to perform running centering, scaling, z-scoring and so on:

``````##              altz
##  [1,]   NaN    NA
##  [2,]  0.71  0.71
##  [3,]  0.89  0.89
##  [4,] -1.18 -1.18
##  [5,]  0.56  0.56
##  [6,]  0.55  0.55
##  [7,] -0.26 -0.26
##  [8,] -0.23 -0.23
##  [9,] -0.23 -0.23
## [10,] -0.51 -0.51
## [11,] -0.17 -0.17
## [12,] -0.59 -0.59
## [13,] -0.19 -0.19
## [14,]  0.84  0.84
## [15,]  2.02  2.02
## [16,]  0.49  0.49
## [17,] -0.22 -0.22
## [18,] -0.82 -0.82
## [19,] -0.64 -0.64
## [20,]  2.37  2.37
``````

A list of the available running functions:

• `running_centered` : from the current value, subtract the mean over the trailing window.
• `running_scaled`: divide the current value by the standard deviation over the trailing window.
• `running_zscored`: from the current value, subtract the mean then divide by the standard deviation over the trailing window.
• `running_sharpe`: divide the mean by the standard deviation over the trailing window. There is a boolean flag to also compute and return the Mertens' form of the standard error of the Sharpe ratio over the trailing window in the second column.
• `running_tstat`: compute the t-stat over the trailing window.
• `running_cumulants`: computes cumulants over the trailing window.
• `running_apx_quantiles`: computes approximate quantiles over the trailing window based on the cumulants and the Cornish-Fisher approximation.
• `running_apx_median`: uses `running_apx_quantiles` to give the approximate median over the trailing window.

The functions `running_centered`, `running_scaled` and `running_zscored` take an optional `lookahead` parameter that allows you to peek ahead (or behind if negative) to the computed moments for comparing against the current value. These are not supported for `running_sharpe` or `running_tstat` because they do not have an idea of the 'current value'.

Here is an example of using the lookahead to z-score some data, compared to a purely time-safe lookback. Around a timestamp of 1000, you can see the difference in outcomes from the two methods:

### Time-Based Running Computations

The standard running moments computations listed above work on a running window of a fixed number of observations. However, sometimes one needs to compute running moments over a different kind of window. The most common form of this is over time-based windows. For example, the following computations:

• Compute the total sales over the past six months, as of every day.
• Compute the volatility of an asset's daily returns, over a yearly window, computed at the end of every trading month.

These are now supported in `fromo` via the `t_running` class of functions, which are like the `running` functions, but accept also the 'times' at which the input are marked, and optionally also the times at which one will 'look back' to perform the computations. The times can be computed implicitly as the cumulative sum of given (non-negative) time deltas.

Here is an example of computing the volatility of daily 'returns' of the Fama French Market factor, based on a one year window, computed at month ends:

date sd mean num_days
2016-07-31 1.09 0.02 251
2016-08-31 0.98 0.05 253
2016-09-30 0.94 0.06 253
2016-10-31 0.91 0.02 252
2016-11-30 0.91 0.04 253
2016-12-31 0.86 0.05 252

And the plot of the time series:

## Efficiency

We make every attempt to balance numerical robustness, computational efficiency and memory usage. As a bit of strawman-bashing, here we microbenchmark the running Z-score computation against the naive algorithm:

``````## Unit: microseconds
##                     expr    min     lq   mean median     uq    max neval cld
##  running_zscored(x, 250)    314    330    352    343    362    449   100  a
##      dumb_zscore(x, 250) 215985 228979 244549 240660 257504 345485   100   b
``````

### Timing against the roll package

More seriously, here we compare the `running_sd3` function, which computes the standard deviation, mean and number of elements with the `roll_sd` and `roll_mean` functions from the roll package.

``````## Unit: microseconds
##                      expr  min   lq mean median   uq  max neval cld
##      running_sd3(xm, 250) 2990 3007 3246   3322 3408 4335   100  b
##  roll::roll_mean(xm, 250)  731  739  849    883  911 1382   100 a
##    roll::roll_sd(xm, 250) 3452 3641 3764   3778 3814 5187   100   c
``````

OK, that's not a fair comparison: `roll_mean` is optimized to work columwise on a matrix. Let's unbash this strawman. I create a function using `fromo::running_sd3` to compute a running mean or running standard deviation columnwise on a matrix, then compare that to `roll_mean` and `roll_sd`:

``````## Error in roll::roll_mean(xm, wins, parallel_for = "cols"): unused argument (parallel_for = "cols")
``````

This was somewhat unexpected. I did not think the `fromo` functions would be faster than `roll_sd`, much less `roll_mean`. From these benchmarks, I suspected that my computer was not taking advantage of parallelization, since the different calls to `roll_mean` had the same timings. To check, I timed `roll_mean` for the same data size, but for larger rolling windows. During the following timing, I confirmed that all eight cores on my laptop were at 100% utilization. I believe the problem is that `roll_mean` is literally recomputing the moments over the entire window for every cell of the output, instead of reusing computations, which `fromo` mostly does:

``````## Unit: microseconds
##                                     expr   min    lq  mean median    uq   max neval   cld
##     roll::roll_mean(xm, 10, min_obs = 3)   630   731   842    786   831  5080   100 a
##    roll::roll_mean(xm, 100, min_obs = 3)   645   731   809    786   835  2352   100 a
##   roll::roll_mean(xm, 1000, min_obs = 3)   642   723   791    765   814  2548   100 a
##  roll::roll_mean(xm, 10000, min_obs = 3)   619   718   810    760   818  2493   100 a
##             fromo_mu(xm, 10, min_df = 3)  5707  5898  6201   6097  6330  8060   100  b
##            fromo_mu(xm, 100, min_df = 3)  6888  7063  7485   7395  7648  9460   100   c
##           fromo_mu(xm, 1000, min_df = 3) 17937 18193 19319  18981 19947 23671   100    d
##          fromo_mu(xm, 10000, min_df = 3) 70518 71236 73727  72855 74553 92000   100     e
``````

The runtime for operations from `roll` grow with the window size. The equivalent operations from `fromo` appear to also consume more time. In theory they would be invariant with respect to window size, but I coded them to 'restart' the computation periodically for improved accuracy. The user has control over how often this happens, in order to balance speed and accuracy. Here I set that parameter very large to show that runtimes need not grow with window size:

``````## Unit: milliseconds
##                                                  expr min  lq mean median  uq  max neval cld
##     fromo_mu(xm, 10, min_df = 3, restart_period = rp) 6.1 6.5  6.8    6.5 6.8  9.7   100   a
##    fromo_mu(xm, 100, min_df = 3, restart_period = rp) 6.2 6.5  6.8    6.5 6.7 10.6   100   a
##   fromo_mu(xm, 1000, min_df = 3, restart_period = rp) 6.2 6.5  6.9    6.6 6.8  9.7   100   a
##  fromo_mu(xm, 10000, min_df = 3, restart_period = rp) 6.5 6.6  7.0    6.7 7.0 10.4   100   a
``````

Here are some more benchmarks, also against the `rollingWindow` package, for running sums:

``````## Unit: microseconds
##                                                      expr min  lq mean median  uq max neval   cld
##                       running_sum(x, wins, na_rm = FALSE)  67  72   74     73  75  84   100 a
##                        RollingWindow::RollingSum(x, wins) 104 108  110    110 111 134   100    d
##                        running_sum(x, wins, na_rm = TRUE)  97 103  106    105 107 174   100   c
##  RollingWindow::RollingSum(x, wins, na_method = "ignore") 268 278  282    281 286 306   100     e
##                                  roll::roll_sum(xm, wins)  92  95   98     96  99 155   100  b
``````

And running means:

``````## Unit: microseconds
##                                                          expr min  lq mean median  uq max neval   cld
##  running_mean(x, wins, na_rm = FALSE, restart_period = 1e+05)  66  71   74     74  75 102   100 a
##                           RollingWindow::RollingMean(x, wins) 129 135  138    138 139 185   100    d
##   running_mean(x, wins, na_rm = TRUE, restart_period = 1e+05)  95 102  105    105 106 138   100  b
##     RollingWindow::RollingMean(x, wins, na_method = "ignore") 288 302  308    309 313 348   100     e
##                                     roll::roll_mean(xm, wins) 103 109  112    111 113 165   100   c
``````

# Reference manual

install.packages("fromo")

0.2.1 by Steven E. Pav, 2 years ago

https://github.com/shabbychef/fromo

Report a bug at https://github.com/shabbychef/fromo/issues

Browse source code at https://github.com/cran/fromo

Authors: Steven E. Pav [aut, cre]

Documentation:   PDF Manual