diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/zSNormalCdf.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/zSNormalCdf.kt
similarity index 79%
rename from kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/zSNormalCdf.kt
rename to kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/zSNormalCdf.kt
index 1eb733ea9..adfbbc9cd 100644
--- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/zSNormalCdf.kt
+++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/zSNormalCdf.kt
@@ -3,19 +3,19 @@
  * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
  */
 
-package space.kscience.kmath.series
+package space.kscience.kmath.distributions
 
 import space.kscience.kmath.operations.DoubleField.pow
 import kotlin.math.PI
 import kotlin.math.absoluteValue
 import kotlin.math.exp
 
-public fun zSNormalCDF(x: Double): Double {
 
-    /**
-     * Zelen & Severo approximation for the standard normal CDF.
-     * The error is bounded by 7.5 * 10e-8.
-     * */
+/**
+ * Zelen & Severo approximation for the standard normal CDF.
+ * The error is bounded by 7.5 * 10e-8.
+ * */
+public fun zSNormalCDF(x: Double): Double {
 
     val t = 1 / (1 + 0.2316419 * x.absoluteValue)
     val summ = 0.319381530*t - 0.356563782*t.pow(2) + 1.781477937*t.pow(3) - 1.821255978*t.pow(4) + 1.330274429*t.pow(5)
diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt
index 72b9ca7a2..cabff25e6 100644
--- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt
+++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt
@@ -216,7 +216,7 @@ public open class SeriesAlgebra<T, out A : Ring<T>, out BA : BufferAlgebra<T, A>
 
     override fun multiply(left: Buffer<T>, right: Buffer<T>): Buffer<T> = left.zip(right) { l, r -> l * r }
 
-    public fun Buffer<T>.diff(shift: Int=1): Buffer<T> = this.zipWithShift(shift) {l, r -> r - l}
+    public fun Buffer<T>.difference(shift: Int=1): Buffer<T> = this.zipWithShift(shift) {l, r -> r - l}
 
     public companion object
 }
diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatioTest.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatioTest.kt
index 280a9db64..3c6747ed2 100644
--- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatioTest.kt
+++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatioTest.kt
@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.series
 
+import space.kscience.kmath.distributions.zSNormalCDF
 import space.kscience.kmath.operations.DoubleField.pow
 import space.kscience.kmath.operations.algebra
 import space.kscience.kmath.operations.bufferAlgebra
@@ -12,62 +13,60 @@ import space.kscience.kmath.operations.fold
 import kotlin.math.absoluteValue
 
 
-// TODO: add p-value with formula: 2*(1 - cdf(|zScore|))
+/**
+ * Container class for Variance Ratio Test result:
+ * ratio itself, corresponding Z-score, also it's p-value
+ * **/
 public data class VarianceRatioTestResult(val varianceRatio: Double=1.0, val zScore: Double=0.0, val pValue: Double=0.5)
-    /**
-     * Container class for Variance Ratio Test result:
-     * ratio itself, corresponding Z-score, also it's p-value
-     * **/
 
-public fun varianceRatioTest(series: Series<Double>, shift: Int, homoscedastic: Boolean=true): VarianceRatioTestResult {
 
-    /**
-     * Calculates the Z-statistic and the p-value for the Lo and MacKinlay's Variance Ratio test (1987)
-     * under Homoscedastic or Heteroscedstic assumptions
-     * with two-sided p-value test
-     * 	https://ssrn.com/abstract=346975
-     * **/
+/**
+ * Calculates the Z-statistic and the p-value for the Lo and MacKinlay's Variance Ratio test (1987)
+ * under Homoscedastic or Heteroscedstic assumptions
+ * with two-sided p-value test
+ * 	https://ssrn.com/abstract=346975
+ * **/
+public fun SeriesAlgebra<Double, *, *, *>.varianceRatioTest(series: Series<Double>, shift: Int, homoscedastic: Boolean=true): VarianceRatioTestResult {
 
     require(shift > 1) {"Shift must be greater than one"}
     require(shift < series.size) {"Shift must be smaller than sample size"}
     val sum = { x: Double, y: Double -> x + y }
 
-    with(Double.algebra.bufferAlgebra.seriesAlgebra()) {
-        val mean = series.fold(0.0, sum) / series.size
-        val demeanedSquares = series.map { power(it - mean, 2) }
-        val variance = demeanedSquares.fold(0.0, sum)
-        if (variance == 0.0) return VarianceRatioTestResult()
+
+    val mean = series.fold(0.0, sum) / series.size
+    val demeanedSquares = series.map { (it - mean).pow(2) }
+    val variance = demeanedSquares.fold(0.0, sum)
+    if (variance == 0.0) return VarianceRatioTestResult()
 
 
-        var seriesAgg = series
-        for (i in 1..<shift) {
-            seriesAgg = seriesAgg.zip(series.moveTo(i)) { v1, v2 -> v1 + v2 }
-        }
-
-        val demeanedSquaresAgg = seriesAgg.map { power(it - shift * mean, 2) }
-        val varianceAgg = demeanedSquaresAgg.fold(0.0, sum)
-
-        val varianceRatio =
-            varianceAgg * (series.size.toDouble() - 1) / variance / (series.size.toDouble() - shift.toDouble() + 1) / (1 - shift.toDouble()/series.size.toDouble()) / shift.toDouble()
-
-
-        // calculating asymptotic variance
-        val phi = if (homoscedastic) {  // under homoscedastic null hypothesis
-            2 * (2 * shift - 1.0) * (shift - 1.0) / (3 * shift * series.size)
-        } else { // under heteroscedastic null hypothesis
-            var accumulator = 0.0
-            for (j in 1..<shift) {
-                val temp = demeanedSquares
-                val delta = series.size * temp.zipWithShift(j) { v1, v2 -> v1 * v2 }.fold(0.0, sum) / variance.pow(2)
-                accumulator += delta * 4 * (shift - j).toDouble().pow(2) / shift.toDouble().pow(2)
-            }
-            accumulator
-        }
-
-        val zScore = (varianceRatio - 1) / phi.pow(0.5)
-        val pValue = 2*(1 - zSNormalCDF(zScore.absoluteValue))
-        return VarianceRatioTestResult(varianceRatio, zScore, pValue)
+    var seriesAgg = series
+    for (i in 1..<shift) {
+        seriesAgg = seriesAgg.zip(series.moveTo(i)) { v1, v2 -> v1 + v2 }
     }
+
+    val demeanedSquaresAgg = seriesAgg.map { (it - shift * mean).pow(2) }
+    val varianceAgg = demeanedSquaresAgg.fold(0.0, sum)
+
+    val varianceRatio =
+        varianceAgg * (series.size.toDouble() - 1) / variance / (series.size.toDouble() - shift.toDouble() + 1) / (1 - shift.toDouble()/series.size.toDouble()) / shift.toDouble()
+
+
+    // calculating asymptotic variance
+    val phi = if (homoscedastic) {  // under homoscedastic null hypothesis
+        2 * (2 * shift - 1.0) * (shift - 1.0) / (3 * shift * series.size)
+    } else { // under heteroscedastic null hypothesis
+        var accumulator = 0.0
+        for (j in 1..<shift) {
+            val temp = demeanedSquares
+            val delta = series.size * temp.zipWithShift(j) { v1, v2 -> v1 * v2 }.fold(0.0, sum) / variance.pow(2)
+            accumulator += delta * 4 * (shift - j).toDouble().pow(2) / shift.toDouble().pow(2)
+        }
+        accumulator
+    }
+
+    val zScore = (varianceRatio - 1) / phi.pow(0.5)
+    val pValue = 2*(1 - zSNormalCDF(zScore.absoluteValue))
+    return VarianceRatioTestResult(varianceRatio, zScore, pValue)
 }