Statistics

algorithms
r
Author

Jason Foster

Published

April 6, 2025

Usage

library(roll)
library(microbenchmark)
options(microbenchmark.unit = "us")
n_vars <- 10
n_obs <- 1000
weights <- 0.9 ^ (n_obs:1)

x <- matrix(rnorm(n_obs * n_vars), nrow = n_obs, ncol = n_vars)
y <- matrix(rnorm(n_obs), nrow = n_obs, ncol = 1)
x_lgl <- x < 0

Rolling any

result <- microbenchmark("125" = roll_any(x_lgl, width = 125, min_obs = 1),
                         "250" = roll_any(x_lgl, width = 250, min_obs = 1),
                         "500" = roll_any(x_lgl, width = 500, min_obs = 1),
                         "1000" = roll_any(x_lgl, width = 1000, min_obs = 1))
print(result)
Unit: microseconds
 expr   min     lq    mean median     uq    max neval
  125  99.2 189.30 216.681 207.40 238.05  380.9   100
  250 158.2 192.05 222.938 208.00 247.55  348.1   100
  500 116.1 184.05 227.864 207.90 231.20 1386.2   100
 1000  88.0 174.65 198.288 193.15 217.25  324.6   100

Rolling all

result <- microbenchmark("125" = roll_all(x_lgl, width = 125, min_obs = 1),
                         "250" = roll_all(x_lgl, width = 250, min_obs = 1),
                         "500" = roll_all(x_lgl, width = 500, min_obs = 1),
                         "1000" = roll_all(x_lgl, width = 1000, min_obs = 1))
print(result)
Unit: microseconds
 expr   min     lq    mean median     uq   max neval
  125 122.5 130.15 157.580 150.30 177.10 296.9   100
  250 121.4 128.55 154.670 146.75 171.75 373.4   100
  500 116.0 124.65 145.850 137.30 165.45 234.9   100
 1000 102.3 116.85 138.581 133.65 161.00 221.3   100

Rolling sums

\[ \begin{aligned} &\text{Expanding window} \\ &\bullet\text{sum}_{x}\leftarrow\lambda\times\text{sum}_{x}+\text{w}_{new}\times\text{x}_{new}\\ &\text{Rolling window}\\ &\bullet\text{sum}_{x}\leftarrow\lambda\times\text{sum}_{x}+\text{w}_{new}\times\text{x}_{new}-\lambda\times\text{w}_{old}\times\text{x}_{old} \end{aligned} \]

result <- microbenchmark("125" = roll_sum(x, width = 125, min_obs = 1, weights = weights),
                         "250" = roll_sum(x, width = 250, min_obs = 1, weights = weights),
                         "500" = roll_sum(x, width = 500, min_obs = 1, weights = weights),
                         "1000" = roll_sum(x, width = 1000, min_obs = 1, weights = weights))
print(result)
Unit: microseconds
 expr   min    lq    mean median    uq   max neval
  125 132.9 149.3 176.270 165.45 187.8 300.0   100
  250 137.5 149.8 175.095 170.75 189.6 299.6   100
  500 131.2 145.2 175.581 157.70 193.8 475.1   100
 1000 126.4 142.7 171.267 163.65 189.7 338.8   100

Rolling products

\[ \begin{aligned} &\text{Expanding window}\\ &\bullet\text{prod}_{w}\leftarrow\text{prod}_{w}\times\text{w}_{new}\\ &\bullet\text{prod}_{x}\leftarrow\text{prod}_{x}\times\text{x}_{new}\\ &\text{Rolling window}\\ &\bullet\text{prod}_{x}\leftarrow\text{prod}_{x}\times\text{x}_{new}/\text{x}_{old} \end{aligned} \]

result <- microbenchmark("125" = roll_prod(x, width = 125, min_obs = 1, weights = weights),
                         "250" = roll_prod(x, width = 250, min_obs = 1, weights = weights),
                         "500" = roll_prod(x, width = 500, min_obs = 1, weights = weights),
                         "1000" = roll_prod(x, width = 1000, min_obs = 1, weights = weights))
print(result)
Unit: microseconds
 expr   min     lq    mean median     uq   max neval
  125 313.9 344.75 374.579 364.00 392.05 688.0   100
  250 314.2 338.10 377.010 367.25 399.80 817.1   100
  500 221.0 251.35 280.739 279.25 305.20 384.3   100
 1000 223.3 243.90 273.549 268.35 295.10 445.2   100

Rolling means

\[ \begin{aligned} &\text{Expanding window}\\ &\bullet\text{sum}_{w}\leftarrow\text{sum}_{w}+\text{w}_{new}\\ &\bullet\text{sum}_{x}\leftarrow\lambda\times\text{sum}_{x}+\text{w}_{new}\times\text{x}_{new}\\ &\text{Rolling window}\\ &\bullet\text{sum}_{x}\leftarrow\lambda\times\text{sum}_{x}+\text{w}_{new}\times\text{x}_{new}-\lambda\times\text{w}_{old}\times \text{x}_{old} \end{aligned} \]

result <- microbenchmark("125" = roll_mean(x, width = 125, min_obs = 1, weights = weights),
                         "250" = roll_mean(x, width = 250, min_obs = 1, weights = weights),
                         "500" = roll_mean(x, width = 500, min_obs = 1, weights = weights),
                         "1000" = roll_mean(x, width = 1000, min_obs = 1, weights = weights))
print(result)
Unit: microseconds
 expr   min     lq    mean median     uq   max neval
  125 134.7 147.05 166.045 156.45 176.25 249.5   100
  250 137.0 147.20 169.917 162.65 178.50 299.1   100
  500 134.5 144.35 165.778 154.55 177.95 242.1   100
 1000 133.5 140.75 163.923 152.25 171.05 374.2   100

Rolling minimums

result <- microbenchmark("125" = roll_min(x, width = 125, min_obs = 1, weights = weights),
                         "250" = roll_min(x, width = 250, min_obs = 1, weights = weights),
                         "500" = roll_min(x, width = 500, min_obs = 1, weights = weights),
                         "1000" = roll_min(x, width = 1000, min_obs = 1, weights = weights))
print(result)
Unit: microseconds
 expr   min     lq    mean median     uq    max neval
  125 151.8 162.20 187.623 181.25 205.05  314.4   100
  250 151.1 161.85 191.762 181.90 202.35  374.3   100
  500 150.1 164.45 220.960 184.25 201.35 3229.3   100
 1000 148.4 169.25 197.918 186.40 213.65  467.3   100

Rolling maximums

result <- microbenchmark("125" = roll_max(x, width = 125, min_obs = 1, weights = weights),
                         "250" = roll_max(x, width = 250, min_obs = 1, weights = weights),
                         "500" = roll_max(x, width = 500, min_obs = 1, weights = weights),
                         "1000" = roll_max(x, width = 1000, min_obs = 1, weights = weights))
print(result)
Unit: microseconds
 expr   min     lq    mean median    uq   max neval
  125 151.0 171.05 198.998 185.45 223.5 445.0   100
  250 149.5 166.00 199.048 185.30 224.1 383.7   100
  500 147.9 174.75 203.926 198.40 217.7 448.3   100
 1000 148.7 171.70 200.275 192.05 224.3 345.6   100

Rolling index of minimums

result <- microbenchmark("125" = roll_idxmin(x, width = 125, min_obs = 1, weights = weights),
                         "250" = roll_idxmin(x, width = 250, min_obs = 1, weights = weights),
                         "500" = roll_idxmin(x, width = 500, min_obs = 1, weights = weights),
                         "1000" = roll_idxmin(x, width = 1000, min_obs = 1, weights = weights))
print(result)
Unit: microseconds
 expr   min    lq    mean median    uq   max neval
  125 146.9 165.0 196.658 185.15 214.2 479.4   100
  250 141.5 163.7 192.498 186.30 210.0 308.5   100
  500 142.0 169.0 195.321 185.15 211.3 298.0   100
 1000 140.8 162.2 189.930 178.25 208.2 326.9   100

Rolling index of maximums

result <- microbenchmark("125" = roll_idxmax(x, width = 125, min_obs = 1, weights = weights),
                         "250" = roll_idxmax(x, width = 250, min_obs = 1, weights = weights),
                         "500" = roll_idxmax(x, width = 500, min_obs = 1, weights = weights),
                         "1000" = roll_idxmax(x, width = 1000, min_obs = 1, weights = weights))
print(result)
Unit: microseconds
 expr   min     lq    mean median     uq   max neval
  125 146.2 168.95 194.865 181.15 213.75 455.3   100
  250 145.8 169.50 189.566 181.20 203.75 278.6   100
  500 141.8 167.20 192.937 186.45 211.25 283.6   100
 1000 145.9 165.95 190.519 180.20 200.60 366.0   100

Rolling medians

# "'online' is only supported for equal 'weights'"
result <- microbenchmark("125" = roll_median(x, width = 125, min_obs = 1),
                         "250" = roll_median(x, width = 250, min_obs = 1),
                         "500" = roll_median(x, width = 500, min_obs = 1),
                         "1000" = roll_median(x, width = 1000, min_obs = 1))
print(result)
Unit: microseconds
 expr    min      lq     mean  median      uq    max neval
  125 1070.3 1136.40 1158.972 1146.60 1162.50 1991.3   100
  250 1037.9 1127.15 1186.371 1140.35 1152.55 5956.6   100
  500  980.6 1044.60 1058.530 1058.00 1073.25 1122.0   100
 1000  742.0  816.90  887.779  828.90  844.90 5909.2   100

Rolling quantiles

# "'online' is only supported for equal 'weights'"
result <- microbenchmark("125" = roll_quantile(x, width = 125, min_obs = 1),
                         "250" = roll_quantile(x, width = 250, min_obs = 1),
                         "500" = roll_quantile(x, width = 500, min_obs = 1),
                         "1000" = roll_quantile(x, width = 1000, min_obs = 1))
print(result)
Unit: microseconds
 expr    min      lq     mean  median      uq    max neval
  125 1103.1 1138.55 1162.070 1157.90 1180.30 1259.1   100
  250 1083.7 1134.95 1160.616 1146.65 1175.00 1746.7   100
  500 1025.8 1053.30 1076.989 1067.40 1092.05 1240.8   100
 1000  786.5  818.15  853.537  840.00  863.20 1782.5   100

Rolling variances

\[ \begin{aligned} &\text{Expanding window}\\ &\bullet\text{sum}_{w}\leftarrow\text{sum}_{w}+\text{w}_{new}\\ &\bullet\text{sumsq}_{w}\leftarrow\text{sumsq}_{w}+\text{w}_{new}^{2}\\ &\bullet\text{sumsq}_{x}\leftarrow\lambda\times\text{sumsq}_{x}+\text{w}_{new}\times (\text{x}_{new}-\text{mean}_{x})(\text{x}_{new}-\text{mean}_{prev_x})\\ &\text{Rolling window}\\ &\bullet\text{sumsq}_{x}\leftarrow\lambda\times\text{sumsq}_{x}+\text{w}_{new}\times (\text{x}_{new}-\text{mean}_{x})(\text{x}_{new}-\text{mean}_{prev_x})-\\ &\lambda\times\text{w}_{old}\times (\text{x}_{old}-\text{mean}_{x})(\text{x}_{old}-\text{mean}_{prev_x}) \end{aligned} \]

result <- microbenchmark("125" = roll_var(x, width = 125, min_obs = 1, weights = weights),
                         "250" = roll_var(x, width = 250, min_obs = 1, weights = weights),
                         "500" = roll_var(x, width = 500, min_obs = 1, weights = weights),
                         "1000" = roll_var(x, width = 1000, min_obs = 1, weights = weights))
print(result)
Unit: microseconds
 expr   min     lq    mean median     uq   max neval
  125 164.3 176.00 199.081 188.70 208.15 350.4   100
  250 162.4 171.45 193.467 180.50 202.30 432.2   100
  500 158.6 167.65 187.993 179.15 201.20 289.8   100
 1000 150.3 156.75 177.466 163.45 189.30 344.1   100

Rolling standard deviations

result <- microbenchmark("125" = roll_sd(x, width = 125, min_obs = 1, weights = weights),
                         "250" = roll_sd(x, width = 250, min_obs = 1, weights = weights),
                         "500" = roll_sd(x, width = 500, min_obs = 1, weights = weights),
                         "1000" = roll_sd(x, width = 1000, min_obs = 1, weights = weights))
print(result)
Unit: microseconds
 expr   min     lq    mean median     uq   max neval
  125 172.5 178.30 207.702 196.95 223.30 482.6   100
  250 165.3 186.05 214.730 200.40 226.80 434.2   100
  500 160.8 173.55 199.249 186.35 217.20 346.4   100
 1000 154.4 166.05 192.072 183.85 208.05 320.1   100

Rolling scaling and centering

result <- microbenchmark("125" = roll_scale(x, width = 125, min_obs = 1, weights = weights),
                         "250" = roll_scale(x, width = 250, min_obs = 1, weights = weights),
                         "500" = roll_scale(x, width = 500, min_obs = 1, weights = weights),
                         "1000" = roll_scale(x, width = 1000, min_obs = 1, weights = weights))
print(result)
Unit: microseconds
 expr   min     lq    mean median    uq   max neval
  125 185.4 200.45 229.177 220.65 242.0 413.9   100
  250 184.7 199.60 229.178 212.45 239.3 369.3   100
  500 175.6 190.15 214.534 205.95 228.2 382.4   100
 1000 171.4 186.35 207.516 198.25 216.1 501.7   100

Rolling covariances

\[ \begin{aligned} &\text{Expanding window}\\ &\bullet\text{sum}_{w}\leftarrow\text{sum}_{w}+\text{w}_{new}\\ &\bullet\text{sumsq}_{w}\leftarrow\text{sumsq}_{w}+\text{w}_{new}^{2}\\ &\bullet\text{sumsq}_{xy}\leftarrow\lambda\times\text{sumsq}_{xy}+\text{w}_{new}\times (\text{x}_{new}-\text{mean}_{x})(\text{y}_{new}-\text{mean}_{prev_y})\\ &\text{Rolling window}\\ &\bullet\text{sumsq}_{xy}\leftarrow\lambda\times\text{sumsq}_{xy}+\text{w}_{new}\times (\text{x}_{new}-\text{mean}_{x})(\text{y}_{new}-\text{mean}_{prev_y})-\\ &\lambda\times\text{w}_{old}\times (\text{x}_{old}-\text{mean}_{x})(\text{y}_{old}-\text{mean}_{prev_y}) \end{aligned} \]

result <- microbenchmark("125" = roll_cov(x, width = 125, min_obs = 1, weights = weights),
                         "250" = roll_cov(x, width = 250, min_obs = 1, weights = weights),
                         "500" = roll_cov(x, width = 500, min_obs = 1, weights = weights),
                         "1000" = roll_cov(x, width = 1000, min_obs = 1, weights = weights))
print(result)
Unit: microseconds
 expr   min     lq     mean  median      uq    max neval
  125 983.5 1210.2 1373.523 1301.85 1391.15 7948.9   100
  250 974.8 1155.4 1263.895 1271.35 1353.90 1749.1   100
  500 953.9 1088.2 1319.174 1189.90 1315.00 5374.9   100
 1000 832.2  960.9 1113.695 1053.00 1119.55 8046.0   100

Rolling correlations

result <- microbenchmark("125" = roll_cor(x, width = 125, min_obs = 1, weights = weights),
                         "250" = roll_cor(x, width = 250, min_obs = 1, weights = weights),
                         "500" = roll_cor(x, width = 500, min_obs = 1, weights = weights),
                         "1000" = roll_cor(x, width = 1000, min_obs = 1, weights = weights))
print(result)
Unit: microseconds
 expr    min      lq     mean  median      uq     max neval
  125 1298.0 1958.90 2832.886 2935.60 3325.65 12622.0   100
  250 1301.2 1987.00 2706.550 2790.90 3179.30 11272.3   100
  500 1139.0 2396.40 2736.490 2771.25 3112.25 10832.3   100
 1000  946.8 1506.05 2173.297 2238.15 2535.00 11392.7   100

Rolling crossproducts

result <- microbenchmark("125" = roll_crossprod(x, width = 125, min_obs = 1, weights = weights),
                         "250" = roll_crossprod(x, width = 250, min_obs = 1, weights = weights),
                         "500" = roll_crossprod(x, width = 500, min_obs = 1, weights = weights),
                         "1000" = roll_crossprod(x, width = 1000, min_obs = 1, weights = weights))
print(result)
Unit: microseconds
 expr   min      lq     mean  median      uq     max neval
  125 894.2 1908.10 2224.067 2007.55 2287.45 10764.4   100
  250 816.4 1816.60 1971.640 1940.10 2210.40 10537.5   100
  500 772.2 1786.10 1880.740 1879.30 2179.00  2773.9   100
 1000 735.1 1658.95 1909.115 1731.25 2016.80 10204.6   100

Rolling linear models

\[ \begin{aligned} &\text{coef}=\text{cov}_{xx}^{-1}\times\text{cov}_{xy}\\ &\text{intercept}=\text{mean}_{y}-\text{coef}\times\text{mean}_{x}\\ &\text{rsq}=\frac{\text{coef}^{T}\times\text{cov}_{xx}\times\text{coef}}{\text{var}_{y}}\\ &\text{var}_{resid}=\frac{(1-\text{rsq})(\text{var}_{y})(\text{sum}_{w}-\text{sumsq}_{w}/\text{sum}_{w})}{\text{n}_{rows}-\text{n}_{cols}}\\ &\text{xx}=\text{cov}_{xx}\times(\text{sum}_{w}-\text{sumsq}_{w}/\text{sum}_{w})\\ &\text{se}_{coef}=\sqrt{\text{var}_{resid}\times\text{diag}(\text{xx}^{-1})}\\ &\text{se}_{intercept}=\sqrt{\text{var}_{resid}\left(1/\text{sum}_{w}+\text{mean}_{x}^{T}\text{xx}^{-1}\text{mean}_{x}\right)} \end{aligned} \]

result <- microbenchmark("125" = roll_lm(x, y, width = 125, min_obs = 1, weights = weights),
                         "250" = roll_lm(x, y, width = 250, min_obs = 1, weights = weights),
                         "500" = roll_lm(x, y, width = 500, min_obs = 1, weights = weights),
                         "1000" = roll_lm(x, y, width = 1000, min_obs = 1, weights = weights))
print(result)
Unit: microseconds
 expr    min      lq     mean  median      uq     max neval
  125 2730.2 3040.25 4572.549 3290.40 6764.15 11336.2   100
  250 2796.3 3070.25 5095.439 3734.95 7294.90 20005.7   100
  500 2736.8 3000.55 4939.481 3221.80 6938.95 17028.6   100
 1000 2580.4 2898.30 4554.263 3227.85 6591.30 12103.0   100

References