Moments
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © dontspamnolan

//@version=5

// @description: Calculates the moments of a source based on Moments (Mean,Variance,Skewness,Kurtosis) [pig]
library("Moments")

// --- Constants
centerLine = 0
zScoreSmallSample = 1.96 // Absolute zscore for N < 50 95% confidence
zScoreLargeSample = 3.29 // Absolute zscore for N > 50 95% confidence
skewReferenceValue = 2.0 // Skewness reference value  for N > 300
kurtosisReferenceValue = 10.0 // Excess Kurtosis reference value  for N > 300 (value for kurtosis is 7, excess kurtosis = 7 + 3)

// @function Calculates log returns of a series (e.g log percentage change)
// @param src Source to use for the returns calculation (e.g. close).
// @returns Log percentage returns of a series
export logReturns(float src = close) =>
    returns = src / src[1]
    math.log(returns)
 
// @function Calculates the mean of a series using ta.sma
// @param src Source to use for the mean calculation (e.g. close).
// @param length Length to use mean calculation (e.g. 14).
// @returns The sma of the source over the length provided.
export mean(float src, int length) =>
    ta.sma(src, length)

// @function Calculates the variance of a series
// @param src Source to use for the variance calculation (e.g. close).
// @param length Length to use for the variance calculation (e.g. 14).
// @returns The variance of the source over the length provided.
export variance(float src, int length) =>
    sample = array.new_float(0)
    sampleMean = mean(src, length)
    for i=0 to length-1
        array.push(sample, math.pow(src[i] - sampleMean, 2))
    sum = array.sum(sample)
    sum / (length-1)

// @function Calculates the standard deviation of a series
// @param src Source to use for the standard deviation calculation (e.g. close).
// @param length Length to use for the standard deviation calculation (e.g. 14).
// @returns The standard deviation of the source over the length provided.
export standardDeviation(float src, int length) =>
    math.sqrt(variance(src, length))

// @function Calculates the skewness of a series
// @param src Source to use for the skewness calculation (e.g. close).
// @param length Length to use for the skewness calculation (e.g. 14).
// @returns The skewness of the source over the length provided.
export skewness(float src, int length) =>
    stdev = standardDeviation(src, length)
    sampleMean = mean(src, length)
    sample = array.new_float(0)
    for i=0 to length-1
        array.push(sample, math.pow(src[i] - sampleMean, 3))
    sum = array.sum(sample)
    a = sum * length
    b = (length-1) * (length-2) * math.pow(stdev, 3)
    a / b

// @function Calculates the kurtosis of a series
// @param src Source to use for the kurtosis calculation (e.g. close).
// @param length Length to use for the kurtosis calculation (e.g. 14).
// @returns The kurtosis of the source over the length provided.
export kurtosis(float src, int length) =>
    stdev = standardDeviation(src, length)
    sampleMean = mean(src, length)
    sample = array.new_float(0)
    for i=0 to length-1
        array.push(sample, math.pow(src[i] - sampleMean, 4))
    sum = array.sum(sample)
    a = length * (length+1) * sum
    b = (length-1) * (length-2) * (length-3) * math.pow(stdev, 4)
    (a / b) - (3 * math.pow((length-1), 2)) / ((length-2) * (length-3))    

// @function Estimates the standard error of skewness based on sample size
// @param sampleSize The number of samples used for calculating standard error.
// @returns The standard error estimate for skewness based on the sample size provided.
export skewnessStandardError(int sampleSize) =>
    math.sqrt((6 * sampleSize * (sampleSize - 1)) / ((sampleSize - 2) * (sampleSize + 1) * (sampleSize + 3)))

// @function Estimates the standard error of kurtosis based on sample size
// @param sampleSize The number of samples used for calculating standard error.
// @returns The standard error estimate for kurtosis based on the sample size provided.    
export kurtosisStandardError(int sampleSize) =>
    skewnessError = skewnessStandardError(sampleSize)
    2 * skewnessError * math.sqrt((math.pow(sampleSize, 2) - 1) / ((sampleSize - 3) * (sampleSize + 5)))

// @function Estimates the critical value of skewness based on sample size
// @param sampleSize The number of samples used for calculating critical value.
// @returns The critical value estimate for skewness based on the sample size provided.     
export skewnessCriticalValue(int sampleSize) =>
    if sampleSize > 300
        skewReferenceValue
    else
        if sampleSize > 50
            zScoreLargeSample * skewnessStandardError(sampleSize)
        else
            zScoreSmallSample * skewnessStandardError(sampleSize)

// @function Estimates the critical value of kurtosis based on sample size
// @param sampleSize The number of samples used for calculating critical value.
// @returns The critical value estimate for kurtosis based on the sample size provided. 
export kurtosisCriticalValue(int sampleSize) =>
    if sampleSize > 300
        kurtosisReferenceValue
    else
        if sampleSize > 50
            zScoreLargeSample * kurtosisStandardError(sampleSize)
        else
            zScoreSmallSample * kurtosisStandardError(sampleSize)
            


// Example
returns = logReturns(close) * 100
length = 14
plot(returns, "Log Returns", color=color.green)
plot(variance(returns, length), "Variance", color=color.blue)
plot(standardDeviation(returns, length), "Stdev", color=color.purple)
plot(skewness(returns, length), "Skewness", color=color.lime)
plot(kurtosis(returns, length), "Kurtosis", color=color.red)
plot(skewnessCriticalValue(length), "Skewness Critical Value", color=color.white)
plot(kurtosisCriticalValue(length), "Kurtosis Critical Value", color=color.black)