Barndorff-Nielsen and Shephard Jump Statistic [Loxx]
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © loxx

// The following comments and descriptions are from from "Problems in the Application of Jump Detection Tests to Stock Price Data" by 
// Michael William Schwert Professor George Tauchen, Faculty Advisor"
//
// This Indicator applies several jump detection tests to intraday stock price data sampled at various 
// frequencies. It finds that the choice of sampling frequency has an effect on both the amount
// of jumps detected by these tests, as well as the timing of those jumps. Furthermore, although 
// these tests are designed to identify the same phenomenon, they find different amounts and timing 
// of jumps when performed on the same data. These results suggest that these jump detection tests 
// are probably identifying different types of jump behavior in stock price data, so they are not 
// really substitutes for one another.

//@version=5
indicator("Barndorff-Nielsen and Shephard Jump Statistic [Loxx]",
     shorttitle = "BNSJS [Loxx]", 
     overlay = false, 
     timeframe="", 
     timeframe_gaps = true)
    
import loxx/loxxexpandedsourcetypes/4
import RicardoSantos/MathSpecialFunctionsGamma/1 as gmm

greencolor = #2DD204
redcolor = #D2042D 

//
//     Inverse of Normal distribution function
//
//
//
// SYNOPSIS:
//
// double x, y, ndtri();
//
// x = ndtri( y );
//
//
//
// DESCRIPTION:
//
// Returns the argument, x, for which the area under the
// Gaussian probability density function (integrated from
// minus infinity to x) is equal to y.
//
//
// For small arguments 0 < y < exp(-2), the program computes
// z = sqrt( -2.0 // log(y) );  then the approximation is
// x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).
// There are two rational functions P/Q, one for 0 < y < exp(-32)
// and the other for y up to exp(-2).  For larger arguments,
// w = y - 0.5, and  x/sqrt(2pi) = w + w////3 R(w////2)/S(w////2)).
//
//
// ACCURACY:
//
//                      Relative error:
// arithmetic   domain        # trials      peak         rms
//    IEEE     0.125, 1        20000       7.2e-16     1.3e-16
//    IEEE     3e-308, 0.135   50000       4.6e-16     9.8e-17
//
//
// ERROR MESSAGES:
//
//   message         condition    value returned
// ndtri domain       x <= 0        -1E10
// ndtri domain       x >= 1         1E10
//
 
// Evaluates polynomial of degree N
polevl(float x, float[] coef, int n)=>
    float ans = array.get(coef, 0)
    for i = 1 to n
        ans := ans * x + array.get(coef, i)
    ans

// Evaluates polynomial of degree N with assumption that coef[N] = 1.0
p1evl(float x, float[] coef, int n)=>
    float ans = x + array.get(coef, 0)
    for i = 1 to n - 1
        ans := ans * x + array.get(coef, i)
    ans

// Cephes Math Library Release 2.1:  January, 1989
// Copyright 1984, 1987, 1989 by Stephen L. Moshier
// Direct inquiries to 30 Frost Street, Cambridge, MA 02140

// sqrt(2pi) 
float s2pi = 2.50662827463100050242E0

// approximation for 0 <= |y - 0.5| <= 3/8 
var array<float> P0 = array.new_float(5, 0)
array.set(P0, 0, -5.99633501014107895267E1)
array.set(P0, 1, 9.80010754185999661536E1)
array.set(P0, 2, -5.66762857469070293439E1)
array.set(P0, 3, 1.39312609387279679503E1)
array.set(P0, 4, -1.23916583867381258016E0)

var array<float> Q0 = array.new_float(8, 0)
array.set(Q0, 0, 1.95448858338141759834E0)
array.set(Q0, 1, 4.67627912898881538453E0)
array.set(Q0, 2, 8.63602421390890590575E1)
array.set(Q0, 3, -2.25462687854119370527E2)
array.set(Q0, 4, 2.00260212380060660359E2)
array.set(Q0, 5, -8.20372256168333339912E1)
array.set(Q0, 6, 1.59056225126211695515E1)
array.set(Q0, 7, -1.18331621121330003142E0)

// Approximation for interval z = sqrt(-2 log y ) between 2 and 8
// i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
var array<float> P1 = array.new_float(9, 0)
array.set(P1, 0, 4.05544892305962419923E0)
array.set(P1, 1, 3.15251094599893866154E1)
array.set(P1, 2, 5.71628192246421288162E1)
array.set(P1, 3, 4.40805073893200834700E1)
array.set(P1, 4, 1.46849561928858024014E1)
array.set(P1, 5, 2.18663306850790267539E0)
array.set(P1, 6, -1.40256079171354495875E-1)
array.set(P1, 7, -3.50424626827848203418E-2)
array.set(P1, 8, -8.57456785154685413611E-4)

var array<float> Q1 = array.new_float(8, 0)
array.set(Q1, 0, 1.57799883256466749731E1)
array.set(Q1, 1, 4.53907635128879210584E1)
array.set(Q1, 2, 4.13172038254672030440E1)
array.set(Q1, 3, 1.50425385692907503408E1)
array.set(Q1, 4, 2.50464946208309415979E0)
array.set(Q1, 5, -1.42182922854787788574E-1)
array.set(Q1, 6, -3.80806407691578277194E-2)
array.set(Q1, 7, -9.33259480895457427372E-4)

// Approximation for interval z = sqrt(-2 log y ) between 8 and 64
// i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
var array<float> P2 = array.new_float(9, 0)
array.set(P2, 0, 3.23774891776946035970E0)
array.set(P2, 1, 6.91522889068984211695E0)
array.set(P2, 2, 3.93881025292474443415E0)
array.set(P2, 3, 1.33303460815807542389E0)
array.set(P2, 4, 2.01485389549179081538E-1)
array.set(P2, 5, 1.23716634817820021358E-2)
array.set(P2, 6, 3.01581553508235416007E-4)
array.set(P2, 7, 2.65806974686737550832E-6)
array.set(P2, 8, 6.23974539184983293730E-9)

var array<float> Q2 = array.new_float(8, 0)
array.set(Q2, 0, 6.02427039364742014255E0)
array.set(Q2, 1, 3.67983563856160859403E0)
array.set(Q2, 2, 1.37702099489081330271E0)
array.set(Q2, 3, 2.16236993594496635890E-1)
array.set(Q2, 4, 1.34204006088543189037E-2)
array.set(Q2, 5, 3.28014464682127739104E-4)
array.set(Q2, 6, 2.89247864745380683936E-6)
array.set(Q2, 7, 6.79019408009981274425E-9)

ndtri(float y0)=>
    float x = 0
    float y = 0
    float z = 0
    float y2 = 0
    float x0 = 0
    float x1 = 0
    int code = 0
    if not(y0 <= 0. or y0 >= 1.)
        code := 1
        y := y0
        // 0.135... = exp(-2) 
        if y > (1.0 - 0.13533528323661269189)
            y := 1.0 - y
            code := 0
        if y > 0.13533528323661269189
            y := y - 0.5
            y2 := y * y
            x := y + y * (y2 * polevl(y2, P0, 4) / p1evl(y2, Q0, 8))
            x := x * s2pi
        else 
            x := math.sqrt(-2.0 *  math.log(y))
            x0 := x - math.log(x) / x
            z := 1.0 / x
            // y > exp(-32) = 1.2664165549e-14 
            if x < 8.
                x1 := z * polevl(z, P1, 8) / p1evl(z, Q1, 8)
            else
                x1 := z * polevl(z, P2, 8) / p1evl(z, Q2, 8)
            x := x0 - x1
            if (code != 0)
                x := -x
    else 
        if y0 <= 0.
            x := -1E10
        if y0 >= 1
            x := 1E10

    x

// Barndorff-Nielsen and Shephard (2004, 2006) developed a test that 
// uses high-frequency price data to determine whether there is a jump over 
// the course of a day. Their test compares two measures of variance: Realized
// Variance, which converges to the integrated variance plus a jump component 
// as the time between observations approaches zero; and Bipower Variation,
// which converges to the integrated variance as the time between observations 
// approaches zero, and is robust to jumps in the price path, an important 
// fact for this application. The integrated variance of a price process is the 
// integral of the square of the variance taken over the course of a day. 
// Since prices cannot be observed continuously, one cannot calculate integrated 
// variance exactly, and must estimate it instead.
//
// For our purposes here, this is caculated as variance = log(p[ti]/p[ti-1]). 
// This the geometric return from time ti-1 to time ti. Then, Realized Variance
// and Bipower Variation are described by the following functions:
//
// realizedVariance(float src, int per) 
//
// and
//
// bipowerVariance(float src, int per)

realizedVariance(float src, int per)=>
    float realizedvariance = 0
    for i = 0 to per 
        realizedvariance += math.pow(src[i], 2)
    realizedvariance

bipowerVariance(float src, int per)=>
    float tempvar = 0
    for i = 1 to per 
        tempvar += src[i] * src[i - 1]
    float bipowervariance = (math.pi / 2.0) * (per / (per - 1.0)) * tempvar
    bipowervariance

// To develop a statistical test to determine whether there is a significant difference between
// RV and BV, one needs an estimate of integrated quarticity. Andersen, Bollerslev, and Diebold
// (2004) recommend using a jump-robust realized Tri-Power Quarticity,
// where Barndorff-Nielsen and Shephard recommend the realized Quad-Power Quarticity. Quad-Power 
// Quarticity won't be covered here in this indicator

tripowerQuarticity(float src, int per)=>
    float mu = math.pow(math.pow(2, 2 / 3) * gmm.GammaLn(7 / 6) * math.pow(gmm.GammaLn(1 / 2), -1), -3)
    float tripower = 0
    for i = 2 to per 
        tripower += math.pow(src[i - 2], 4 / 3) * math.pow(src[i - 1], 4 / 3) * math.pow(src[i], 4 / 3)
    float tpquarticity = per * mu * (per / (per - 2)) * tripower
    tpquarticity

barndorffNielsenStatistic(float src, int per)=>
    float logreturns = math.log(src / src[1])
    float abslogreturns = math.abs(logreturns)
    float realized_variance = realizedVariance(logreturns, per)
    float bipower_variance =  bipowerVariance(abslogreturns, per)
    float tripower = tripowerQuarticity(abslogreturns, per)

    // Huang and Tauchen (2005) also consider Relative Jump, 
    // a measure that approximates the percentage of total 
    // variance attributable to jumps.This statistic approximates the 
    // ratio of the sum of squared jumps to the total variance and 
    // is useful because it scales  out long-term trends in volatility
    // so one can compare the  relative contribution of jumps to
    // the variance of two price  series with different volatilities.

    float relative_jump = (realized_variance - bipower_variance) / realized_variance

    // There are numerous potential test statistics that use these 
    // estimates of integrated variance and integrated quarticity, 
    // but Huang and Tauchen find that these two asymptotically standard 
    // normal test statistics perform best in simulations:

    float statistic = relative_jump / math.sqrt((math.pow(math.pi / 2, 2) + math.pi - 5) * (1.0 / per) * math.max(1, tripower / (math.pow(bipower_variance, 2))))
    statistic

smthtype = input.string("Kaufman", "Heikin-Ashi Better Caculation Type", options = ["AMA", "T3", "Kaufman"], group = "Source Settings")
srcin = input.string("Close", "Source", group= "Source Settings", 
     options = 
     ["Close", "Open", "High", "Low", "Median", "Typical", "Weighted", "Average", "Average Median Body", "Trend Biased", "Trend Biased (Extreme)", 
     "HA Close", "HA Open", "HA High", "HA Low", "HA Median", "HA Typical", "HA Weighted", "HA Average", "HA Average Median Body", "HA Trend Biased", "HA Trend Biased (Extreme)",
     "HAB Close", "HAB Open", "HAB High", "HAB Low", "HAB Median", "HAB Typical", "HAB Weighted", "HAB Average", "HAB Average Median Body", "HAB Trend Biased", "HAB Trend Biased (Extreme)"])

per = input.int(16, "Jump Window", group = "Basic Settings", tooltip = "For SPY benchmark: 7, 16, 78, 110, 156, and 270 returns for sampling intervals of 1 week, 1 day, 1 hour, 30 minutes, 15 minutes, and 5 minutes, respectively.")
alpha = input.float(0.1, "Percent-Point Function (PPF) Alpha", group = "Basic Settings", minval = 0, maxval = 1, step = 0.01)

colorbars = input.bool(true, "Color bars?", group = "UI Options")
showSigs = input.bool(true, "Show signals?", group = "UI Options")

kfl=input.float(0.666, title="* Kaufman's Adaptive MA (KAMA) Only - Fast End", group = "Moving Average Inputs")
ksl=input.float(0.0645, title="* Kaufman's Adaptive MA (KAMA) Only - Slow End", group = "Moving Average Inputs")
amafl = input.int(2, title="* Adaptive Moving Average (AMA) Only - Fast", group = "Moving Average Inputs")
amasl = input.int(30, title="* Adaptive Moving Average (AMA) Only - Slow", group = "Moving Average Inputs")

haclose = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, close)
haopen = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, open)
hahigh = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, high)
halow = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, low)
hamedian = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, hl2)
hatypical = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, hlc3)
haweighted = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, hlcc4)
haaverage = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, ohlc4)

src = switch srcin
	"Close" => loxxexpandedsourcetypes.rclose()
	"Open" => loxxexpandedsourcetypes.ropen()
	"High" => loxxexpandedsourcetypes.rhigh()
	"Low" => loxxexpandedsourcetypes.rlow()
	"Median" => loxxexpandedsourcetypes.rmedian()
	"Typical" => loxxexpandedsourcetypes.rtypical()
	"Weighted" => loxxexpandedsourcetypes.rweighted()
	"Average" => loxxexpandedsourcetypes.raverage()
    "Average Median Body" => loxxexpandedsourcetypes.ravemedbody()
	"Trend Biased" => loxxexpandedsourcetypes.rtrendb()
	"Trend Biased (Extreme)" => loxxexpandedsourcetypes.rtrendbext()
	"HA Close" => loxxexpandedsourcetypes.haclose(haclose)
	"HA Open" => loxxexpandedsourcetypes.haopen(haopen)
	"HA High" => loxxexpandedsourcetypes.hahigh(hahigh)
	"HA Low" => loxxexpandedsourcetypes.halow(halow)
	"HA Median" => loxxexpandedsourcetypes.hamedian(hamedian)
	"HA Typical" => loxxexpandedsourcetypes.hatypical(hatypical)
	"HA Weighted" => loxxexpandedsourcetypes.haweighted(haweighted)
	"HA Average" => loxxexpandedsourcetypes.haaverage(haaverage)
    "HA Average Median Body" => loxxexpandedsourcetypes.haavemedbody(haclose, haopen)
	"HA Trend Biased" => loxxexpandedsourcetypes.hatrendb(haclose, haopen, hahigh, halow)
	"HA Trend Biased (Extreme)" => loxxexpandedsourcetypes.hatrendbext(haclose, haopen, hahigh, halow)
	"HAB Close" => loxxexpandedsourcetypes.habclose(smthtype, amafl, amasl, kfl, ksl)
	"HAB Open" => loxxexpandedsourcetypes.habopen(smthtype, amafl, amasl, kfl, ksl)
	"HAB High" => loxxexpandedsourcetypes.habhigh(smthtype, amafl, amasl, kfl, ksl)
	"HAB Low" => loxxexpandedsourcetypes.hablow(smthtype, amafl, amasl, kfl, ksl)
	"HAB Median" => loxxexpandedsourcetypes.habmedian(smthtype, amafl, amasl, kfl, ksl)
	"HAB Typical" => loxxexpandedsourcetypes.habtypical(smthtype, amafl, amasl, kfl, ksl)
	"HAB Weighted" => loxxexpandedsourcetypes.habweighted(smthtype, amafl, amasl, kfl, ksl)
	"HAB Average" => loxxexpandedsourcetypes.habaverage(smthtype, amafl, amasl, kfl, ksl)
    "HAB Average Median Body" => loxxexpandedsourcetypes.habavemedbody(smthtype, amafl, amasl, kfl, ksl)
	"HAB Trend Biased" => loxxexpandedsourcetypes.habtrendb(smthtype, amafl, amasl, kfl, ksl)
	"HAB Trend Biased (Extreme)" => loxxexpandedsourcetypes.habtrendbext(smthtype, amafl, amasl, kfl, ksl)
	=> haclose

statistic = barndorffNielsenStatistic(src, per)
quantile = ndtri(1 - alpha)

colorout = statistic < quantile ? greencolor : color.yellow
plot(statistic, "BNSJS", style = plot.style_histogram, linewidth = 2, color = colorout)

plot(quantile, "Cutoff", color = redcolor)

colorout2 = statistic < quantile ? na : color.yellow

barcolor(colorbars ? colorout2 : na)

plot(3, "Empirical Z-Score Maximum", color = bar_index % 2 ? color.gray : na)

gojump = ta.crossover(statistic, quantile)

plotshape(showSigs and gojump, title = "Jump Detected", color = color.yellow, textcolor = color.yellow, text = "J", style = shape.cross, location = location.bottom, size = size.auto)
alertcondition(gojump, title = "Jump Detected", message = "Barndorff-Nielsen and Shephard Jump Statistic [Loxx]: Long\nSymbol: {{ticker}}\nPrice: {{close}}")