From 1d7f4ed538efcbafe5527d5e26635f160d3a003a Mon Sep 17 00:00:00 2001 From: mrFendel Date: Thu, 6 Apr 2023 03:17:01 +0300 Subject: [PATCH 01/15] shiftOp and diff in SeriesAlgebra --- .../kscience/kmath/series/SeriesAlgebra.kt | 36 ++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) 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 3cd2212f6..bdd74c4fe 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 @@ -3,16 +3,21 @@ package space.kscience.kmath.series import space.kscience.kmath.operations.BufferAlgebra import space.kscience.kmath.operations.Ring import space.kscience.kmath.operations.RingOps +import space.kscience.kmath.operations.reduce import space.kscience.kmath.stat.StatisticalAlgebra import space.kscience.kmath.structures.Buffer +import space.kscience.kmath.structures.BufferFactory import space.kscience.kmath.structures.BufferView import kotlin.math.max import kotlin.math.min + +// TODO: check if ranges are intersected @PublishedApi internal fun IntRange.intersect(other: IntRange): IntRange = max(first, other.first)..min(last, other.last) + @PublishedApi internal val IntRange.size: Int get() = last - first + 1 @@ -33,11 +38,14 @@ public interface Series : Buffer { public val position: Int } + + public val Series.absoluteIndices: IntRange get() = position until position + size /** * A [BufferView] with index offset (both positive and negative) and possible size change */ + private class SeriesImpl( override val origin: Buffer, override val position: Int, @@ -107,6 +115,8 @@ public class SeriesAlgebra, out BA : BufferAlgebra, L>( * Get a label buffer for given buffer. */ public val Buffer.labels: List get() = indices.map(labelResolver) + // TODO: there can be troubles with label consistency after moving position argument + // TODO: so offset should be reflected in the labelResolver also /** @@ -145,6 +155,17 @@ public class SeriesAlgebra, out BA : BufferAlgebra, L>( return accumulator } + // TODO: add fold with recorded accumulation +// public inline fun Buffer.traceFold(initial: R, operation: A.(acc: R, T) -> R): Buffer { +// var tempBuffer = elementAlgebra.bufferFactory(this.size) {i -> getAbsolute(i)} +// var accumulator = initial +// for (index in this.indices) { +// accumulator = elementAlgebra.operation(accumulator, getAbsolute(index)) +// tempBuffer.set(index, accumulator) +// } +// return elementAlgebra.bufferFactory(this.size) {i -> tempBuffer.getAbsolute(i)} +// } + public inline fun Buffer.foldWithLabel(initial: R, operation: A.(acc: R, arg: T, label: L) -> R): R { val labels = labels var accumulator = initial @@ -154,7 +175,7 @@ public class SeriesAlgebra, out BA : BufferAlgebra, L>( } /** - * Zip two buffers in the range whe they overlap + * Zip two buffers in the range where they overlap */ public inline fun Buffer.zip( other: Buffer, @@ -169,11 +190,24 @@ public class SeriesAlgebra, out BA : BufferAlgebra, L>( }.moveTo(newRange.first) } + /** + * Zip buffer with itself, but shifted + * */ + public inline fun Buffer.shiftOp( + shift: Int = 1, + crossinline operation: A.(left: T, right: T) -> T + ): Buffer { + val shifted = this.moveTo(this.offset+shift) + return zip(shifted, operation) + } + override fun Buffer.unaryMinus(): Buffer = map { -it } override fun add(left: Buffer, right: Buffer): Series = left.zip(right) { l, r -> l + r } override fun multiply(left: Buffer, right: Buffer): Buffer = left.zip(right) { l, r -> l * r } + + public inline fun Buffer.diff(): Buffer = this.shiftOp {l, r -> r - l} } public fun , BA : BufferAlgebra, L> BA.seriesAlgebra(labels: Iterable): SeriesAlgebra { From ba26c7020e76a5b2384f37d515471f287a37617c Mon Sep 17 00:00:00 2001 From: mrFendel Date: Thu, 6 Apr 2023 17:58:29 +0300 Subject: [PATCH 02/15] started TimeSeriesAlgebra --- .../kscience/kmath/series/SeriesAlgebra.kt | 11 ------ .../kmath/series/TimeSeriesAlgebra.kt | 37 +++++++++++++++++++ 2 files changed, 37 insertions(+), 11 deletions(-) create mode 100644 kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/TimeSeriesAlgebra.kt 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 bdd74c4fe..b6d952955 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 @@ -155,17 +155,6 @@ public class SeriesAlgebra, out BA : BufferAlgebra, L>( return accumulator } - // TODO: add fold with recorded accumulation -// public inline fun Buffer.traceFold(initial: R, operation: A.(acc: R, T) -> R): Buffer { -// var tempBuffer = elementAlgebra.bufferFactory(this.size) {i -> getAbsolute(i)} -// var accumulator = initial -// for (index in this.indices) { -// accumulator = elementAlgebra.operation(accumulator, getAbsolute(index)) -// tempBuffer.set(index, accumulator) -// } -// return elementAlgebra.bufferFactory(this.size) {i -> tempBuffer.getAbsolute(i)} -// } - public inline fun Buffer.foldWithLabel(initial: R, operation: A.(acc: R, arg: T, label: L) -> R): R { val labels = labels var accumulator = initial diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/TimeSeriesAlgebra.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/TimeSeriesAlgebra.kt new file mode 100644 index 000000000..37abe710e --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/TimeSeriesAlgebra.kt @@ -0,0 +1,37 @@ +/* + * 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.structures.Buffer + + + +public interface TimeSeries: Series { + public override val origin: Buffer + public override val position: Int + + // TODO: Specify the types of DateTime that can be contained in timeStamp Buffer + public val timeStamp: Buffer // Buffer of some datetime instances like: Instant, LocalDate, LocalTime... +} + +private class TimeSeriesImpl( + override val origin: Buffer, + override val timeStamp: Buffer, + override val position: Int, + override val size: Int = origin.size, +) : TimeSeries by origin { // TODO: manage with delegation + + init { + require(size > 0) { "Size must be positive" } + require(size <= origin.size) { "Slice size is larger than the original buffer" } + require(size <= timeStamp.size) { "Slice size is larger than the timeStamp buffer" } + } + +// override fun toString(): String = "$origin-->${position}" +} + + +// TODO: add conversion to Buffer of Pairs ? From 4db091d8984d269cc8abb08dc369277ef78bb942 Mon Sep 17 00:00:00 2001 From: mrFendel Date: Fri, 7 Apr 2023 12:39:30 +0300 Subject: [PATCH 03/15] deleted TimeSeriesAlgebra --- .../kmath/series/TimeSeriesAlgebra.kt | 37 ------------------- 1 file changed, 37 deletions(-) delete mode 100644 kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/TimeSeriesAlgebra.kt diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/TimeSeriesAlgebra.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/TimeSeriesAlgebra.kt deleted file mode 100644 index 37abe710e..000000000 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/TimeSeriesAlgebra.kt +++ /dev/null @@ -1,37 +0,0 @@ -/* - * 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.structures.Buffer - - - -public interface TimeSeries: Series { - public override val origin: Buffer - public override val position: Int - - // TODO: Specify the types of DateTime that can be contained in timeStamp Buffer - public val timeStamp: Buffer // Buffer of some datetime instances like: Instant, LocalDate, LocalTime... -} - -private class TimeSeriesImpl( - override val origin: Buffer, - override val timeStamp: Buffer, - override val position: Int, - override val size: Int = origin.size, -) : TimeSeries by origin { // TODO: manage with delegation - - init { - require(size > 0) { "Size must be positive" } - require(size <= origin.size) { "Slice size is larger than the original buffer" } - require(size <= timeStamp.size) { "Slice size is larger than the timeStamp buffer" } - } - -// override fun toString(): String = "$origin-->${position}" -} - - -// TODO: add conversion to Buffer of Pairs ? From 31d1cc774af7ce56a8d401691ecb08e9a24a2b1c Mon Sep 17 00:00:00 2001 From: mrFendel Date: Tue, 11 Apr 2023 20:31:04 +0300 Subject: [PATCH 04/15] added shiftOperartion and diff --- .../space/kscience/kmath/series/SeriesAlgebra.kt | 13 +++++++++++++ 1 file changed, 13 insertions(+) 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..4b7f8b83a 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 @@ -199,12 +199,25 @@ public open class SeriesAlgebra, out BA : BufferAlgebra } } + /** + * Zip buffer with itself, but shifted + * */ + public inline fun Buffer.shiftOp( + shift: Int = 1, + crossinline operation: A.(left: T, right: T) -> T + ): Buffer { + val shifted = this.moveTo(this.startOffset+shift) + return zip(shifted, operation) + } + override fun Buffer.unaryMinus(): Buffer = map { -it } override fun add(left: Buffer, right: Buffer): Series = left.zip(right) { l, r -> l + r } override fun multiply(left: Buffer, right: Buffer): Buffer = left.zip(right) { l, r -> l * r } + public inline fun Buffer.diff(): Buffer = this.shiftOp {l, r -> r - l} + public companion object } From 2b83560da8ff62b58d0df8bb334a76f3db3ccf09 Mon Sep 17 00:00:00 2001 From: mrFendel Date: Wed, 12 Apr 2023 22:24:48 +0300 Subject: [PATCH 05/15] Variance Ratio function --- .../kscience/kmath/series/VarianceRatio.kt | 30 +++++++++++++++++++ 1 file changed, 30 insertions(+) create mode 100644 kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatio.kt diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatio.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatio.kt new file mode 100644 index 000000000..aa616487b --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatio.kt @@ -0,0 +1,30 @@ +/* + * 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.DoubleBufferOps.Companion.map +import space.kscience.kmath.operations.algebra +import space.kscience.kmath.operations.bufferAlgebra +import space.kscience.kmath.operations.fold + + +fun varianceRatio(series: Series, shift: Int): Double { + val mean = series.fold(0.0) {acc, value -> acc + value} / series.size + val demeanedSquares = series.map { power(it - mean, 2) } + val variance = demeanedSquares.fold(0.0) {acc, value -> acc + value} + + with(Double.algebra.bufferAlgebra.seriesAlgebra()) { + val seriesAgg = series + for (i in -1..-shift + 1) { + seriesAgg.shiftOp(i) { v1, v2 -> v1 + v2 } + } + + val demeanedSquaresAgg = seriesAgg.map { power(it - shift * mean, 2) } + val varianceAgg = demeanedSquaresAgg.fold(0.0) { acc, value -> acc + value } + + return varianceAgg * (series.size - 1) / variance / (series.size - shift + 1) / (1 - shift / series.size) + } +} From a68ebef26dda77e015c1d8144999c40083ee303f Mon Sep 17 00:00:00 2001 From: mrFendel Date: Thu, 13 Apr 2023 03:38:10 +0300 Subject: [PATCH 06/15] zScore for variance Ratio Test --- .../kscience/kmath/series/VarianceRatio.kt | 42 ++++++++++++++----- 1 file changed, 32 insertions(+), 10 deletions(-) diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatio.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatio.kt index aa616487b..2dfcd6b00 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatio.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatio.kt @@ -5,26 +5,48 @@ package space.kscience.kmath.series +import space.kscience.kmath.distributions.NormalDistribution 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 -fun varianceRatio(series: Series, shift: Int): Double { - val mean = series.fold(0.0) {acc, value -> acc + value} / series.size +// TODO: add p-value +public data class VarianceRatioTestResult(val varianceRatio: Double, val zScore: Double) + +public fun varianceRatioTest(series: Series, shift: Int, homoscedastic: Boolean): VarianceRatioTestResult { + + val sum = { x: Double, y: Double -> x + y } + + val mean = series.fold(0.0, sum) / series.size val demeanedSquares = series.map { power(it - mean, 2) } - val variance = demeanedSquares.fold(0.0) {acc, value -> acc + value} + val variance = demeanedSquares.fold(0.0, sum) with(Double.algebra.bufferAlgebra.seriesAlgebra()) { - val seriesAgg = series - for (i in -1..-shift + 1) { - seriesAgg.shiftOp(i) { v1, v2 -> v1 + v2 } + for (i in -1..-shift + 1) { series.shiftOp(i) { v1, v2 -> v1 + v2 } } + val demeanedSquaresAgg = series.map { power(it - shift * mean, 2) } + val varianceAgg = demeanedSquaresAgg.fold(0.0, sum) + + val varianceRatio = + varianceAgg * (series.size - 1) / variance / (series.size - shift + 1) / (1 - shift / series.size) + + + // calculating asymptotic variance + var phi: Double + if (homoscedastic) { // under homoscedastic null hypothesis + phi = 2 * (2 * shift - 1.0) * (shift - 1.0) / (3 * shift * series.size) + } else { // under homoscedastic null hypothesis + phi = 0.0 + for (j in 1.. v1 * v2 } + val delta = series.size * shiftedProd.fold(0.0, sum) / variance.pow(2) + phi += delta * 4 * (shift - j) * (shift - j) / shift / shift // TODO: refactor with square + } } - val demeanedSquaresAgg = seriesAgg.map { power(it - shift * mean, 2) } - val varianceAgg = demeanedSquaresAgg.fold(0.0) { acc, value -> acc + value } - - return varianceAgg * (series.size - 1) / variance / (series.size - shift + 1) / (1 - shift / series.size) + val zScore = (varianceRatio - 1) / phi.pow(0.5) + return VarianceRatioTestResult(varianceRatio, zScore) } } From 0ce1861ab4724912e796ffa92aa940fdaa0735a0 Mon Sep 17 00:00:00 2001 From: mrFendel Date: Thu, 13 Apr 2023 03:47:36 +0300 Subject: [PATCH 07/15] refactoring --- .../series/{VarianceRatio.kt => VarianceRatioTest.kt} | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) rename kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/{VarianceRatio.kt => VarianceRatioTest.kt} (90%) diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatio.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatioTest.kt similarity index 90% rename from kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatio.kt rename to kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatioTest.kt index 2dfcd6b00..77163f638 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatio.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatioTest.kt @@ -5,7 +5,6 @@ package space.kscience.kmath.series -import space.kscience.kmath.distributions.NormalDistribution import space.kscience.kmath.operations.DoubleBufferOps.Companion.map import space.kscience.kmath.operations.DoubleField.pow import space.kscience.kmath.operations.algebra @@ -18,6 +17,12 @@ public data class VarianceRatioTestResult(val varianceRatio: Double, val zScore: public fun varianceRatioTest(series: Series, shift: Int, homoscedastic: Boolean): VarianceRatioTestResult { + /** + * Calculate 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 + * **/ + val sum = { x: Double, y: Double -> x + y } val mean = series.fold(0.0, sum) / series.size From a91b43a52db60388a4947c1a277e6b26bd3e6997 Mon Sep 17 00:00:00 2001 From: mrFendel Date: Thu, 13 Apr 2023 17:52:14 +0300 Subject: [PATCH 08/15] tests for varianceRatio --- .../kmath/series/VarianceRatioTest.kt | 4 +- .../kmath/series/TestVarianceRatioTest.kt | 51 +++++++++++++++++++ 2 files changed, 53 insertions(+), 2 deletions(-) create mode 100644 kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestVarianceRatioTest.kt 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 77163f638..b769d78a3 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 @@ -24,10 +24,10 @@ public fun varianceRatioTest(series: Series, shift: Int, homoscedastic: * **/ val sum = { x: Double, y: Double -> x + y } - + //TODO: catch if shift is too large val mean = series.fold(0.0, sum) / series.size val demeanedSquares = series.map { power(it - mean, 2) } - val variance = demeanedSquares.fold(0.0, sum) + val variance = demeanedSquares.fold(0.0, sum) // TODO: catch if variance is zero with(Double.algebra.bufferAlgebra.seriesAlgebra()) { for (i in -1..-shift + 1) { series.shiftOp(i) { v1, v2 -> v1 + v2 } } 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..7c31663bc --- /dev/null +++ b/kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestVarianceRatioTest.kt @@ -0,0 +1,51 @@ +/* + * 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 volatileData() { + with(Double.algebra.bufferAlgebra.seriesAlgebra()) { + val volatileData = series(10) { sin(PI * it + PI/2) + 1.0} + val resultHomo = varianceRatioTest(volatileData, 2, homoscedastic = true) + 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) + } + } + + @Test + fun negativeData() { + with(Double.algebra.bufferAlgebra.seriesAlgebra()) { + val volatileData = series(10) { sin(PI * it)} + val resultHomo = varianceRatioTest(volatileData, 2, homoscedastic = true) + assertEquals(1.142857, resultHomo.varianceRatio, 1e-6) + // homoscedastic zScore + assertEquals(0.451753, resultHomo.zScore, 1e-6) + val resultHetero = varianceRatioTest(volatileData, 2, homoscedastic = false) + // heteroscedastic zScore + assertEquals(2.462591, resultHetero.zScore, 1e-6) + } + } + +// @Test +// fun zeroVolatility() { +// with(Double.algebra.bufferAlgebra.seriesAlgebra()) { +// val volatileData = series(10) { 1.0 } +// val result = varianceRatioTest(volatileData, 2, homoscedastic = true) +// } +// } +} \ No newline at end of file From 5b95923bb9072c3437e8a3f0b65449aa90b896ef Mon Sep 17 00:00:00 2001 From: mrFendel Date: Fri, 14 Apr 2023 06:36:20 +0300 Subject: [PATCH 09/15] fixed zip in SereiesAlgebra + tests for VarianceRatio --- .../kscience/kmath/series/SeriesAlgebra.kt | 2 +- .../kmath/series/VarianceRatioTest.kt | 22 ++++++---- .../kmath/series/TestVarianceRatioTest.kt | 40 +++++++++++++------ 3 files changed, 43 insertions(+), 21 deletions(-) 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 4b7f8b83a..9efbd629c 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, out BA : BufferAlgebra crossinline operation: A.(left: T, right: T) -> T, ): Series { 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) 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 b769d78a3..9a00b1be2 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 @@ -25,17 +25,22 @@ public fun varianceRatioTest(series: Series, shift: Int, homoscedastic: val sum = { x: Double, y: Double -> x + y } //TODO: catch if shift is too large - val mean = series.fold(0.0, sum) / series.size - val demeanedSquares = series.map { power(it - mean, 2) } - val variance = demeanedSquares.fold(0.0, sum) // TODO: catch if variance is zero - with(Double.algebra.bufferAlgebra.seriesAlgebra()) { - for (i in -1..-shift + 1) { series.shiftOp(i) { v1, v2 -> v1 + v2 } } - val demeanedSquaresAgg = series.map { power(it - shift * mean, 2) } + val mean = series.fold(0.0, sum) / series.size + val demeanedSquares = series.map { power(it - mean, 2) } + val variance = demeanedSquares.fold(0.0, sum) // TODO: catch if variance is zero + + + var seriesAgg = series + for (i in 1.. v1 + v2 } + } + + val demeanedSquaresAgg = seriesAgg.map { power(it - shift * mean, 2) } val varianceAgg = demeanedSquaresAgg.fold(0.0, sum) val varianceRatio = - varianceAgg * (series.size - 1) / variance / (series.size - shift + 1) / (1 - shift / series.size) + varianceAgg * (series.size.toDouble() - 1) / variance / (series.size.toDouble() - shift.toDouble() + 1) / (1 - shift.toDouble()/series.size.toDouble()) / shift.toDouble() // calculating asymptotic variance @@ -44,8 +49,9 @@ public fun varianceRatioTest(series: Series, shift: Int, homoscedastic: phi = 2 * (2 * shift - 1.0) * (shift - 1.0) / (3 * shift * series.size) } else { // under homoscedastic null hypothesis phi = 0.0 + var shiftedProd = demeanedSquares for (j in 1.. v1 * v2 } + shiftedProd = shiftedProd.zip(demeanedSquares.moveTo(j)) { v1, v2 -> v1 * v2 } val delta = series.size * shiftedProd.fold(0.0, sum) / variance.pow(2) phi += delta * 4 * (shift - j) * (shift - j) / shift / shift // TODO: refactor with square } 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 7c31663bc..8dcdb3a14 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,6 +13,21 @@ import kotlin.test.assertEquals class TestVarianceRatioTest { + // TODO: refactor Heteroscedastic zScore + @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) +// val resultHetero = varianceRatioTest(monotonicData, 2, homoscedastic = false) +// // heteroscedastic zScore +// assertEquals(3.253248, resultHetero.zScore, 1e-6) + } + } + @Test fun volatileData() { with(Double.algebra.bufferAlgebra.seriesAlgebra()) { @@ -21,31 +36,32 @@ 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(-3.535533, resultHetero.zScore, 1e-6) } } @Test fun negativeData() { with(Double.algebra.bufferAlgebra.seriesAlgebra()) { - val volatileData = series(10) { sin(PI * it)} - val resultHomo = varianceRatioTest(volatileData, 2, homoscedastic = true) - assertEquals(1.142857, resultHomo.varianceRatio, 1e-6) + val negativeData = series(10) { sin(it * 1.2)} + val resultHomo = varianceRatioTest(negativeData, 3, homoscedastic = true) + assertEquals(1.240031, resultHomo.varianceRatio, 1e-6) // homoscedastic zScore - assertEquals(0.451753, resultHomo.zScore, 1e-6) - val resultHetero = varianceRatioTest(volatileData, 2, homoscedastic = false) - // heteroscedastic zScore - assertEquals(2.462591, resultHetero.zScore, 1e-6) + assertEquals(0.509183, resultHomo.zScore, 1e-6) +// val resultHetero = varianceRatioTest(negativeData, 3, homoscedastic = false) +// // heteroscedastic zScore +// assertEquals(0.661798, resultHetero.zScore, 1e-6) } } + //TODO: add zero volatility Test, logReturns test, big shift Test // @Test // fun zeroVolatility() { // with(Double.algebra.bufferAlgebra.seriesAlgebra()) { -// val volatileData = series(10) { 1.0 } -// val result = varianceRatioTest(volatileData, 2, homoscedastic = true) +// val zeroVolData = series(10) { 1.0 } +// val result = varianceRatioTest(zeroVolData, 2, homoscedastic = true) // } // } } \ No newline at end of file From e6da61c52a7c7b38cabe5621014993eef5bdf2d2 Mon Sep 17 00:00:00 2001 From: mrFendel Date: Tue, 18 Apr 2023 01:53:07 +0300 Subject: [PATCH 10/15] refactoring --- .../kscience/kmath/series/SeriesAlgebra.kt | 6 +++--- .../kscience/kmath/series/VarianceRatioTest.kt | 18 +++++++++++------- 2 files changed, 14 insertions(+), 10 deletions(-) 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 9efbd629c..72b9ca7a2 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 @@ -202,11 +202,11 @@ public open class SeriesAlgebra, out BA : BufferAlgebra /** * Zip buffer with itself, but shifted * */ - public inline fun Buffer.shiftOp( + public inline fun Buffer.zipWithShift( shift: Int = 1, crossinline operation: A.(left: T, right: T) -> T ): Buffer { - val shifted = this.moveTo(this.startOffset+shift) + val shifted = this.moveBy(shift) return zip(shifted, operation) } @@ -216,7 +216,7 @@ public open class SeriesAlgebra, out BA : BufferAlgebra override fun multiply(left: Buffer, right: Buffer): Buffer = left.zip(right) { l, r -> l * r } - public inline fun Buffer.diff(): Buffer = this.shiftOp {l, r -> r - l} + public fun Buffer.diff(shift: Int=1): Buffer = 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 9a00b1be2..45bc836fe 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 @@ -12,10 +12,14 @@ import space.kscience.kmath.operations.bufferAlgebra import space.kscience.kmath.operations.fold -// TODO: add p-value +// TODO: add p-value with formula: 2*(1 - cdf(|zScore|)) public data class VarianceRatioTestResult(val varianceRatio: Double, val zScore: Double) + /** + * Container class for Variance Ratio Test result: + * ratio itself, corresponding Z-score, also it's p-value + * **/ -public fun varianceRatioTest(series: Series, shift: Int, homoscedastic: Boolean): VarianceRatioTestResult { +public fun varianceRatioTest(series: Series, shift: Int, homoscedastic: Boolean=true): VarianceRatioTestResult { /** * Calculate the Z statistic and the p-value for the Lo and MacKinlay's Variance Ratio test (1987) @@ -44,17 +48,17 @@ public fun varianceRatioTest(series: Series, shift: Int, homoscedastic: // calculating asymptotic variance - var phi: Double - if (homoscedastic) { // under homoscedastic null hypothesis - phi = 2 * (2 * shift - 1.0) * (shift - 1.0) / (3 * shift * series.size) + 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 - phi = 0.0 + var accumulator = 0.0 var shiftedProd = demeanedSquares for (j in 1.. v1 * v2 } val delta = series.size * shiftedProd.fold(0.0, sum) / variance.pow(2) - phi += delta * 4 * (shift - j) * (shift - j) / shift / shift // TODO: refactor with square + accumulator += delta * 4 * (shift - j) * (shift - j) / shift / shift // TODO: refactor with square } + accumulator } val zScore = (varianceRatio - 1) / phi.pow(0.5) From 98781c83adcbb62ea0a83f5c905b6bf6a85ea59c Mon Sep 17 00:00:00 2001 From: mrFendel Date: Tue, 18 Apr 2023 19:16:10 +0300 Subject: [PATCH 11/15] fixed bug with heteroscedastic z-score in Variance Ratio Test --- .../kmath/series/VarianceRatioTest.kt | 13 ++++++------- .../kmath/series/TestVarianceRatioTest.kt | 19 +++++++++---------- 2 files changed, 15 insertions(+), 17 deletions(-) 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, 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, 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.. 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) } } From 0193349f9413a79fb1b06e05a437b8e609e0a6c3 Mon Sep 17 00:00:00 2001 From: mrFendel Date: Wed, 19 Apr 2023 01:36:54 +0300 Subject: [PATCH 12/15] requirements, default parameters, new Test for varianceRatioTest --- .../space/kscience/kmath/ejml/_generated.kt | 2 +- .../kmath/series/VarianceRatioTest.kt | 12 ++++++----- .../kmath/series/TestVarianceRatioTest.kt | 21 ++++++++++--------- 3 files changed, 19 insertions(+), 16 deletions(-) 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/series/VarianceRatioTest.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/VarianceRatioTest.kt index 8afc01c81..d81a6c317 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 @@ -9,11 +9,10 @@ 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|)) -public data class VarianceRatioTestResult(val varianceRatio: Double, val zScore: Double) +public data class VarianceRatioTestResult(val varianceRatio: Double=1.0, val zScore: Double=0.0) /** * Container class for Variance Ratio Test result: * ratio itself, corresponding Z-score, also it's p-value @@ -27,12 +26,15 @@ public fun varianceRatioTest(series: Series, shift: Int, homoscedastic: * https://ssrn.com/abstract=346975 * **/ + 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 } - //TODO: catch if shift is too large + 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) // TODO: catch if variance is zero + val variance = demeanedSquares.fold(0.0, sum) + if (variance == 0.0) return VarianceRatioTestResult() var seriesAgg = series @@ -53,7 +55,7 @@ public fun varianceRatioTest(series: Series, shift: Int, homoscedastic: } else { // under heteroscedastic null hypothesis var accumulator = 0.0 for (j in 1.. v1 * v2 }.fold(0.0, sum) / variance.pow(2) accumulator += delta * 4 * (shift - j).toDouble().pow(2) / shift.toDouble().pow(2) } 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 0ece143eb..8e559276c 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 @@ -31,7 +31,7 @@ class TestVarianceRatioTest { fun volatileData() { with(Double.algebra.bufferAlgebra.seriesAlgebra()) { val volatileData = series(10) { sin(PI * it + PI/2) + 1.0} - val resultHomo = varianceRatioTest(volatileData, 2, homoscedastic = true) + val resultHomo = varianceRatioTest(volatileData, 2) assertEquals(0.0, resultHomo.varianceRatio, 1e-6) // homoscedastic zScore assertEquals(-3.162277, resultHomo.zScore, 1e-6) @@ -45,7 +45,7 @@ class TestVarianceRatioTest { fun negativeData() { with(Double.algebra.bufferAlgebra.seriesAlgebra()) { val negativeData = series(10) { sin(it * 1.2)} - val resultHomo = varianceRatioTest(negativeData, 3, homoscedastic = true) + val resultHomo = varianceRatioTest(negativeData, 3) assertEquals(1.240031, resultHomo.varianceRatio, 1e-6) // homoscedastic zScore assertEquals(0.509183, resultHomo.zScore, 1e-6) @@ -55,12 +55,13 @@ class TestVarianceRatioTest { } } - //TODO: add zero volatility Test, logReturns test, big shift Test -// @Test -// fun zeroVolatility() { -// with(Double.algebra.bufferAlgebra.seriesAlgebra()) { -// val zeroVolData = series(10) { 1.0 } -// val result = varianceRatioTest(zeroVolData, 2, homoscedastic = true) -// } -// } + @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) + } + } } \ No newline at end of file From 1e27af9cf577478a2b6b582e9882022f91e36610 Mon Sep 17 00:00:00 2001 From: mrFendel Date: Wed, 19 Apr 2023 17:13:47 +0300 Subject: [PATCH 13/15] - Zelen-Severo CDF aproximation - p-value for varianceRatioTest --- .../kmath/series/VarianceRatioTest.kt | 11 +++++++-- .../kscience/kmath/series/zSNormalCdf.kt | 24 +++++++++++++++++++ .../kmath/series/TestVarianceRatioTest.kt | 5 ++++ 3 files changed, 38 insertions(+), 2 deletions(-) create mode 100644 kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/zSNormalCdf.kt 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 d81a6c317..280a9db64 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 @@ -9,10 +9,11 @@ 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 kotlin.math.absoluteValue // TODO: add p-value with formula: 2*(1 - cdf(|zScore|)) -public data class VarianceRatioTestResult(val varianceRatio: Double=1.0, val zScore: Double=0.0) +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 @@ -23,6 +24,7 @@ public fun varianceRatioTest(series: Series, shift: Int, homoscedastic: /** * 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 * **/ @@ -63,6 +65,11 @@ public fun varianceRatioTest(series: Series, shift: Int, homoscedastic: } val zScore = (varianceRatio - 1) / phi.pow(0.5) - return VarianceRatioTestResult(varianceRatio, zScore) + val pValue = 2*(1 - zSNormalCDF(zScore.absoluteValue)) + return VarianceRatioTestResult(varianceRatio, zScore, pValue) } } + + + + diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/zSNormalCdf.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/zSNormalCdf.kt new file mode 100644 index 000000000..1eb733ea9 --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/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.series + +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. + * */ + + 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/commonTest/kotlin/space/kscience/kmath/series/TestVarianceRatioTest.kt b/kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestVarianceRatioTest.kt index 8e559276c..afc0d541d 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 @@ -21,9 +21,11 @@ class TestVarianceRatioTest { 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) } } @@ -35,9 +37,11 @@ class TestVarianceRatioTest { 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) } } @@ -62,6 +66,7 @@ class TestVarianceRatioTest { 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 From 16385b5f4e8e15fcb27774ff41b5443be0cd9642 Mon Sep 17 00:00:00 2001 From: mrFendel Date: Fri, 5 May 2023 18:45:54 +0300 Subject: [PATCH 14/15] -- refactoring --- .../{series => distributions}/zSNormalCdf.kt | 12 +-- .../kscience/kmath/series/SeriesAlgebra.kt | 2 +- .../kmath/series/VarianceRatioTest.kt | 89 +++++++++---------- 3 files changed, 51 insertions(+), 52 deletions(-) rename kmath-stat/src/commonMain/kotlin/space/kscience/kmath/{series => distributions}/zSNormalCdf.kt (79%) 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, out BA : BufferAlgebra override fun multiply(left: Buffer, right: Buffer): Buffer = left.zip(right) { l, r -> l * r } - public fun Buffer.diff(shift: Int=1): Buffer = this.zipWithShift(shift) {l, r -> r - l} + public fun Buffer.difference(shift: Int=1): Buffer = 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, 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.varianceRatioTest(series: Series, 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.. 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.. 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.. 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.. 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) } From f44039e30915530e260c1de6fb9d955c3b76d73b Mon Sep 17 00:00:00 2001 From: mrFendel Date: Fri, 5 May 2023 18:50:10 +0300 Subject: [PATCH 15/15] -- refactoring --- .../kotlin/space/kscience/kmath/series/VarianceRatioTest.kt | 2 -- 1 file changed, 2 deletions(-) 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 3c6747ed2..6a2f5e426 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 @@ -7,8 +7,6 @@ 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 import space.kscience.kmath.operations.fold import kotlin.math.absoluteValue