diff --git a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt
index c56583fa8..8ad7f7293 100644
--- a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt
+++ b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt
@@ -19,9 +19,9 @@ import org.ejml.sparse.csc.factory.DecompositionFactory_DSCC
 import org.ejml.sparse.csc.factory.DecompositionFactory_FSCC
 import org.ejml.sparse.csc.factory.LinearSolverFactory_DSCC
 import org.ejml.sparse.csc.factory.LinearSolverFactory_FSCC
-import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.linear.*
 import space.kscience.kmath.linear.Matrix
+import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.nd.StructureFeature
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.FloatField
diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/zSNormalCdf.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/zSNormalCdf.kt
new file mode 100644
index 000000000..adfbbc9cd
--- /dev/null
+++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/zSNormalCdf.kt
@@ -0,0 +1,24 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * 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.distributions
+
+import space.kscience.kmath.operations.DoubleField.pow
+import kotlin.math.PI
+import kotlin.math.absoluteValue
+import kotlin.math.exp
+
+
+/**
+ * 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)
+    val temp = summ * exp(-x.absoluteValue.pow(2) / 2) / (2 * PI).pow(0.5)
+    return if (x >= 0)  1 - temp else temp
+}
\ No newline at end of file
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 3ccbca5a1..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
@@ -191,7 +191,7 @@ public open class SeriesAlgebra<T, out A : Ring<T>, out BA : BufferAlgebra<T, A>
         crossinline operation: A.(left: T, right: T) -> T,
     ): Series<T> {
         val newRange = offsetIndices.intersect(other.offsetIndices)
-        return seriesByOffset(startOffset = newRange.first, size = newRange.last - newRange.first) { offset ->
+        return seriesByOffset(startOffset = newRange.first, size = newRange.last + 1 - newRange.first) { offset ->
             elementAlgebra.operation(
                 getByOffset(offset),
                 other.getByOffset(offset)
@@ -199,12 +199,25 @@ public open class SeriesAlgebra<T, out A : Ring<T>, out BA : BufferAlgebra<T, A>
         }
     }
 
+    /**
+     * Zip buffer with itself, but shifted
+     * */
+    public inline fun Buffer<T>.zipWithShift(
+        shift: Int = 1,
+        crossinline operation: A.(left: T, right: T) -> T
+    ): Buffer<T> {
+        val shifted = this.moveBy(shift)
+        return zip(shifted, operation)
+    }
+
     override fun Buffer<T>.unaryMinus(): Buffer<T> = map { -it }
 
     override fun add(left: Buffer<T>, right: Buffer<T>): Series<T> = left.zip(right) { l, r -> l + r }
 
     override fun multiply(left: Buffer<T>, right: Buffer<T>): Buffer<T> = left.zip(right) { l, r -> l * r }
 
+    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
new file mode 100644
index 000000000..6a2f5e426
--- /dev/null
+++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatioTest.kt
@@ -0,0 +1,72 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * 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
+
+import space.kscience.kmath.distributions.zSNormalCDF
+import space.kscience.kmath.operations.DoubleField.pow
+import space.kscience.kmath.operations.fold
+import kotlin.math.absoluteValue
+
+
+/**
+ * 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)
+
+
+/**
+ * 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 }
+
+
+    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 { (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)
+}
+
+
+
+
diff --git a/kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestVarianceRatioTest.kt b/kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestVarianceRatioTest.kt
new file mode 100644
index 000000000..afc0d541d
--- /dev/null
+++ b/kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestVarianceRatioTest.kt
@@ -0,0 +1,72 @@
+/*
+ * Copyright 2018-2023 KMath contributors.
+ * 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
+
+import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.operations.bufferAlgebra
+import kotlin.math.PI
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class TestVarianceRatioTest {
+
+    @Test
+    fun monotonicData() {
+        with(Double.algebra.bufferAlgebra.seriesAlgebra()) {
+            val monotonicData = series(10) { it * 1.0 }
+            val resultHomo = varianceRatioTest(monotonicData, 2, homoscedastic = true)
+            assertEquals(1.818181, resultHomo.varianceRatio, 1e-6)
+            // homoscedastic zScore
+            assertEquals(2.587318, resultHomo.zScore, 1e-6)
+            assertEquals(.0096, resultHomo.pValue, 1e-4)
+            val resultHetero = varianceRatioTest(monotonicData, 2, homoscedastic = false)
+            // heteroscedastic zScore
+            assertEquals(0.819424, resultHetero.zScore, 1e-6)
+            assertEquals(.4125, resultHetero.pValue, 1e-4)
+        }
+    }
+
+    @Test
+    fun volatileData() {
+        with(Double.algebra.bufferAlgebra.seriesAlgebra()) {
+            val volatileData = series(10) { sin(PI * it + PI/2) + 1.0}
+            val resultHomo = varianceRatioTest(volatileData, 2)
+            assertEquals(0.0, resultHomo.varianceRatio, 1e-6)
+            // homoscedastic zScore
+            assertEquals(-3.162277, resultHomo.zScore, 1e-6)
+            assertEquals(.0015, resultHomo.pValue, 1e-4)
+            val resultHetero = varianceRatioTest(volatileData, 2, homoscedastic = false)
+            // heteroscedastic zScore
+            assertEquals(-1.0540925, resultHetero.zScore, 1e-6)
+            assertEquals(.2918, resultHetero.pValue, 1e-4)
+        }
+    }
+
+    @Test
+    fun negativeData() {
+        with(Double.algebra.bufferAlgebra.seriesAlgebra()) {
+            val negativeData = series(10) { sin(it * 1.2)}
+            val resultHomo = varianceRatioTest(negativeData, 3)
+            assertEquals(1.240031, resultHomo.varianceRatio, 1e-6)
+            // homoscedastic zScore
+            assertEquals(0.509183, resultHomo.zScore, 1e-6)
+            val resultHetero = varianceRatioTest(negativeData, 3, homoscedastic = false)
+            // heteroscedastic zScore
+            assertEquals(0.209202, resultHetero.zScore, 1e-6)
+        }
+    }
+
+    @Test
+    fun zeroVolatility() {
+        with(Double.algebra.bufferAlgebra.seriesAlgebra()) {
+            val zeroVolData = series(10) { 0.0 }
+            val result = varianceRatioTest(zeroVolData, 4)
+            assertEquals(1.0, result.varianceRatio, 1e-6)
+            assertEquals(0.0, result.zScore, 1e-6)
+            assertEquals(0.5, result.pValue, 1e-4)
+        }
+    }
+}
\ No newline at end of file