forked from kscience/kmath
-- refactoring
This commit is contained in:
parent
1e27af9cf5
commit
16385b5f4e
@ -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)
|
@ -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
|
||||
}
|
||||
|
@ -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)
|
||||
}
|
||||
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user