// Licensed under the Apache License, Version 2.0
// © mihakralj
//@version=6
indicator("Jarque-Bera Test (JB)", "JB", overlay=false, precision=4)

// Helper function to get a window of series data into an array
_get_window_array(series float source, simple int length) =>
    float[] arr = array.new_float(length)
    for i = 0 to length - 1
        array.set(arr, i, source[length - 1 - i]) // Oldest to newest
    arr

// Helper function to calculate the k-th central moment
// m_k = sum((x_i - mean)^k) / n
_central_moment(float[] arr, int moment_order, float mean_val) =>
    n = array.size(arr)
    if n == 0
        float(na)
    else
        sum_pow_diff = 0.0
        for i = 0 to n - 1
            sum_pow_diff += math.pow(array.get(arr, i) - mean_val, moment_order)
        sum_pow_diff / n

//@function Calculates the Jarque-Bera statistic.
//@param source series float The input series.
//@param length simple int The lookback period (sample size). Min 10.
//@returns series float The Jarque-Bera statistic. Higher values suggest deviation from normality.
jb_stat(series float source, simple int length) => // Renamed function for clarity jb_stat -> jb
    if length < 10 // Need sufficient sample size
        float(na)
    else
        float[] window_arr = _get_window_array(source, length)
        
        float mean_val = array.avg(window_arr) // Pine Script built-in array average
        
        if na(mean_val)
            float(na)
        else
            // Calculate variance (2nd central moment)
            float m2 = _central_moment(window_arr, 2, mean_val)
            
            if na(m2) or m2 < 1e-10 // Avoid division by zero or near-zero std dev
                float(na) // Or 0 if data is truly constant, but JB is ill-defined
            else
                float stddev_val = math.sqrt(m2)
                
                // Calculate 3rd central moment for skewness
                float m3 = _central_moment(window_arr, 3, mean_val)
                float s = m3 / math.pow(stddev_val, 3) // Skewness
                
                // Calculate 4th central moment for kurtosis
                float m4 = _central_moment(window_arr, 4, mean_val)
                float k_raw = m4 / math.pow(stddev_val, 4) // Raw Kurtosis
                float ek = k_raw - 3.0 // Excess Kurtosis
                
                if na(s) or na(ek)
                    float(na)
                else
                    // Jarque-Bera statistic formula: JB = (n/6) * (S^2 + (EK^2)/4)
                    float jb_value_calc = (length / 6.0) * (s * s + (ek * ek) / 4.0) // Renamed internal variable
                    jb_value_calc

// Inputs
i_source = input.source(close, title="Source")
i_length = input.int(20, title="Lookback Period (Sample Size)", minval=10, maxval=200, tooltip="Number of bars for calculation. Min 10. Max 200 due to array processing. Higher values provide more stable estimates but lag more.")

// Calculation
jb_value = jb_stat(i_source, i_length) // Call renamed function

// Plot
plot(jb_value, "Jarque-Bera Statistic", color=color.teal, linewidth=2)

// Critical values for Chi-squared distribution with 2 degrees of freedom (approximate):
// Significance Level | Critical Value
// 10% (0.10)         | 4.605
// 5%  (0.05)         | 5.991
// 1%  (0.01)         | 9.210
// These lines can help interpret the JB statistic. If JB > critical value, reject null hypothesis of normality.
hline(4.605, "Critical Value (10%)", color.gray, linestyle=hline.style_dotted, linewidth=1)
hline(5.991, "Critical Value (5%)", color.orange, linestyle=hline.style_dashed, linewidth=1)
hline(9.210, "Critical Value (1%)", color.red, linestyle=hline.style_solid, linewidth=1)