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 45bc836fe..8afc01c81 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,11 +5,11 @@
 
 package space.kscience.kmath.series
 
-import space.kscience.kmath.operations.DoubleBufferOps.Companion.map
 import space.kscience.kmath.operations.DoubleField.pow
 import space.kscience.kmath.operations.algebra
 import space.kscience.kmath.operations.bufferAlgebra
 import space.kscience.kmath.operations.fold
+import space.kscience.kmath.structures.slice
 
 
 // TODO: add p-value with formula: 2*(1 - cdf(|zScore|))
@@ -22,7 +22,7 @@ public data class VarianceRatioTestResult(val varianceRatio: Double, val zScore:
 public fun varianceRatioTest(series: Series<Double>, shift: Int, homoscedastic: Boolean=true): VarianceRatioTestResult {
 
     /**
-     * Calculate the Z statistic and the p-value for the Lo and MacKinlay's Variance Ratio test (1987)
+     * Calculates the Z-statistic and the p-value for the Lo and MacKinlay's Variance Ratio test (1987)
      * under Homoscedastic or Heteroscedstic assumptions
      * 	https://ssrn.com/abstract=346975
      * **/
@@ -50,13 +50,12 @@ public fun varianceRatioTest(series: Series<Double>, shift: Int, homoscedastic:
         // 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 homoscedastic null hypothesis
+        } else { // under heteroscedastic null hypothesis
             var accumulator = 0.0
-            var shiftedProd = demeanedSquares
             for (j in 1..<shift) {
-                shiftedProd = shiftedProd.zip(demeanedSquares.moveTo(j)) { v1, v2 -> v1 * v2 }
-                val delta = series.size * shiftedProd.fold(0.0, sum) / variance.pow(2)
-                accumulator += delta * 4 * (shift - j) * (shift - j) / shift / shift // TODO: refactor with square
+                var 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
         }
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
index 8dcdb3a14..0ece143eb 100644
--- a/kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestVarianceRatioTest.kt
+++ b/kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestVarianceRatioTest.kt
@@ -13,7 +13,6 @@ import kotlin.test.assertEquals
 
 class TestVarianceRatioTest {
 
-    // TODO: refactor Heteroscedastic zScore
     @Test
     fun monotonicData() {
         with(Double.algebra.bufferAlgebra.seriesAlgebra()) {
@@ -22,9 +21,9 @@ class TestVarianceRatioTest {
             assertEquals(1.818181, resultHomo.varianceRatio, 1e-6)
             // homoscedastic zScore
             assertEquals(2.587318, resultHomo.zScore, 1e-6)
-//            val resultHetero = varianceRatioTest(monotonicData, 2, homoscedastic = false)
-//            // heteroscedastic zScore
-//            assertEquals(3.253248, resultHetero.zScore, 1e-6)
+            val resultHetero = varianceRatioTest(monotonicData, 2, homoscedastic = false)
+            // heteroscedastic zScore
+            assertEquals(0.819424, resultHetero.zScore, 1e-6)
         }
     }
 
@@ -36,9 +35,9 @@ class TestVarianceRatioTest {
             assertEquals(0.0, resultHomo.varianceRatio, 1e-6)
             // homoscedastic zScore
             assertEquals(-3.162277, resultHomo.zScore, 1e-6)
-//            val resultHetero = varianceRatioTest(volatileData, 2, homoscedastic = false)
-//            // heteroscedastic zScore
-//            assertEquals(-3.535533, resultHetero.zScore, 1e-6)
+            val resultHetero = varianceRatioTest(volatileData, 2, homoscedastic = false)
+            // heteroscedastic zScore
+            assertEquals(-1.0540925, resultHetero.zScore, 1e-6)
         }
     }
 
@@ -50,9 +49,9 @@ class TestVarianceRatioTest {
             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.661798, resultHetero.zScore, 1e-6)
+            val resultHetero = varianceRatioTest(negativeData, 3, homoscedastic = false)
+            // heteroscedastic zScore
+            assertEquals(0.209202, resultHetero.zScore, 1e-6)
         }
     }