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

//@version=5
//@description double precision algorithm to compute the cumulative trivariate normal distribution 
//	found in A.Genz, Numerical computation of rectangular bivariate and trivariate normal 
//	and t probabilities”, Statistics and Computing, 14, (3), 2004. The cumulative trivariate 
//	normal is needed to price window barrier options, see G.F. Armstrong, Valuation formulae
// 	or window barrier options”, Applied Mathematical Finance, 8, 2001.
// references:
// 		https://link.springer.com/article/10.1023/B%3ASTCO.0000035304.20635.31
//      https://www.tandfonline.com/doi/abs/10.1080/13504860210124607
//		https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.353.1954&rep=rep1&type=pdf
//		The Complete Guide to Option Pricing Formulas, 2nd ed. (Espen Gaarder Haug) 
//
library("ctnd")

import loxx/cbnd/1

TRVFND(float T, float B2, float B3, float RHO21, float RHO31, float rho)=>
	TRVFND = math.exp(-T * T / 2) * cbnd.CBND3(-T * RHO21 + B2, -T * RHO31 + B3, rho)
	TRVFND

KRNRDD(float A, float b, float ABSERR, float b2p,float b3p, float RHO21, float RHO31, float rho)=>
//
//           The abscissae and weights are given for the interval (-1,1)
//           because of symmetry only the positive abscisae and their
//           corresponding weights are given.
//
//           XGK    - abscissae of the 2N+1-point Kronrod rule:
//                    XGK(2), XGK(4), ...  N-point Gauss rule abscissae;
//                    XGK(1), XGK(3), ...  abscissae optimally added
//                    to the N-point Gauss rule.
//
//           WGK    - weights of the 2N+1-point Kronrod rule.
//
//           WG     - weights of the N-point Gauss rule.

	float ABSCIS = 0
	float CENTER = 0
	float FC = 0
	float FUNSUM = 0
	float HFLGTH = 0
	float RESLTG = 0
	float RESLTK = 0

	var int N = 12
	//

	float KRNRDD = 0
	float[] WG = array.new<float>(int((N + 1) / 2), 0) 
	float[] WGK = array.new<float>(N, 0) 
	float[] XGK = array.new<float>(N, 0) 

	array.set(WG, 1, 0.0556685671161745)
	array.set(WG, 2, 0.125580369464905)
	array.set(WG, 3, 0.186290210927735)
	array.set(WG, 4, 0.233193764591991)
	array.set(WG, 5, 0.262804544510248)
	array.set(WG, 0, 0.272925086777901)
	//
	array.set(XGK, 1, 0.996369613889543)
	array.set(XGK, 2, 0.978228658146057)
	array.set(XGK, 3, 0.941677108578068)
	array.set(XGK, 4, 0.887062599768095)
	array.set(XGK, 5, 0.816057456656221)
	array.set(XGK, 6, 0.730152005574049)
	array.set(XGK, 7, 0.630599520161965)
	array.set(XGK, 8, 0.519096129206812)
	array.set(XGK, 9, 0.397944140952378)
	array.set(XGK, 10, 0.269543155952345)
	array.set(XGK, 11, 0.136113000799362)
	array.set(XGK, 0, 0)
	//
	array.set(WGK, 1, 0.00976544104596129)
	array.set(WGK, 2, 0.0271565546821044)
	array.set(WGK, 3, 0.0458293785644267)
	array.set(WGK, 4, 0.0630974247503748)
	array.set(WGK, 5, 0.0786645719322276)
	array.set(WGK, 6, 0.0929530985969007)
	array.set(WGK, 7, 0.105872074481389)
	array.set(WGK, 8, 0.116739502461047)
	array.set(WGK, 9, 0.125158799100319)
	array.set(WGK, 10, 0.131280684229806)
	array.set(WGK, 11, 0.135193572799885)
	array.set(WGK, 0, 0.136577794711118)

//
//
//           List of major variables
//
//           CENTER  - mid point of the interval
//           HFLGTH  - half-length of the interval
//           ABSCIS   - abscissae
//           RESLTG   - result of the N-point Gauss formula
//           RESLTK   - result of the 2N+1-point Kronrod formula
//
//

	HFLGTH := (b - A) / 2
	CENTER := (b + A) / 2
//
//           compute the 2N+1-point Kronrod approximation to
//           the integral, and estimate the absolute error.
//

	FC := TRVFND(CENTER, b2p, b3p, RHO21, RHO31, rho)
	RESLTG := FC * array.get(WG, 0)
	RESLTK := FC * array.get(WGK, 0)
	for j = 1 to N - 1
		ABSCIS := HFLGTH * array.get(XGK, j)
		FUNSUM := TRVFND(CENTER - ABSCIS, b2p, b3p, RHO21, RHO31, rho) 
			 + TRVFND(CENTER + ABSCIS, b2p, b3p, RHO21, RHO31, rho)

		RESLTK := RESLTK + array.get(WGK, j) * FUNSUM
		if j % 2 == 0 
			RESLTG := RESLTG + array.get(WG, int(j / 2)) * FUNSUM
	KRNRDD := RESLTK * HFLGTH
	KRNRDD
	//ABSERR := 3 * math.abs((RESLTK - RESLTG) * HFLGTH)

ADONED(float A, float b, float TOL, float b2p, float b3p, float RHO21, float RHO31, float rho)=>
//
//     One Dimensional Globally Adaptive Integration Function
//
	int IM = 1
	int IP = 1
	var int NL = 100
	float[] EI = array.new<float>(NL, 0) 
	float[] AI = array.new<float>(NL, 0) 
	float[] BI = array.new<float>(NL, 0) 
	float[] FI = array.new<float>(NL, 0) 
	float FIN = 0
	float ERR = 0
		
	array.set(AI, IP, A)
	array.set(BI, IP, b)
	array.set(FI, IP, KRNRDD(array.get(AI, IP), array.get(BI, IP), array.get(EI, IP), b2p, b3p, RHO21, RHO31, rho))

	while IM < NL - 1
		IM := IM + 1
		array.set(BI, IM, array.get(BI, IP))
		array.set(AI, IM, array.get(AI, IP) + array.get(BI, IP) / 2)
		array.set(BI, IP, array.get(AI, IM))
		FIN := array.get(FI, IP)
		array.set(FI, IP, KRNRDD(array.get(AI, IP), array.get(BI, IP), array.get(EI, IP), b2p, b3p, RHO21, RHO31, rho))
		array.set(FI, IM, KRNRDD(array.get(AI, IM), array.get(BI, IM), array.get(EI, IM), b2p, b3p, RHO21, RHO31, rho))
		ERR := math.abs(FIN - array.get(FI, IP) - array.get(FI, IM)) / 2
		array.set(EI, IP, array.get(EI, IP) + ERR)
		array.set(EI, IM, array.get(EI, IM) + ERR)
		IP := 1
		ERR := 0
		FIN := 0
		for i = 1 to IM
			if array.get(EI, i) > array.get(EI, IP) 
				IP := i
			FIN := FIN + array.get(FI, i)
			ERR := ERR + array.get(EI, i)
	ADONED = FIN
	ADONED


// @function Returns the Cumulative Trivariate Normal Distribution 
// @param LIMIT1 float, 
// @param LIMIT2 float, 
// @param LIMIT3 float, 
// @param SIGMA1 float, 
// @param SIGMA2 float, 
// @param SIGMA3 float, 
// @returns float.
export CTND(float LIMIT1, float LIMIT2, float LIMIT3, float SIGMA1, float SIGMA2, float SIGMA3)=>
//
//     A function for computing trivariate normal probabilities.
//     This function uses an algorithm given in the paper
//        "Numerical Computation of Bivariate and
//             Trivariate Normal Probabilities",
//     by
//       Alan Genz
//       Department of Mathematics
//       Washington State University
//       Pullman, WA 99164-3113
//       Email : alangenz@wsu.edu
//
// CTND calculates the probability that X(I) < LIMIT(I), I = 1, 2, 3.
//
// Parameters
//
//   LIMIT  DOUBLE PRECISION array of three upper integration limits.
//   SIGMA  DOUBLE PRECISION array of three correlation coefficients,
//          SIGMA should contain the lower left portion of the
//          correlation matrix R.
//          SIGMA(1) = R(2,1), SIGMA(2) = R(3,1), SIGMA(3) = R(3,2).
//
//    CTND cuts the outer integral over -infinity to B1 to
//      an integral from -8.5 to B1 and then uses an adaptive
//      integration method to compute the integral of a bivariate
//      normal distribution function.
//
	bool TAIL = false
//
//     Bivariate normal distribution function CBND is required.
//
	float SQ21 = 0
	float SQ31 = 0
	float rho = 0
	float b1 = 0
	float B2 = 0
	float B3 = 0
	float b2p = 0
	float b3p  = 0
	float RHO21 = 0
	float RHO31 = 0
	float RHO32  = 0
	var float SQTWPI = 2.506628274631
	var float XCUT = -8.5
	var float EPS = 5E-16
	float CTND = 0
	
	b1 := LIMIT1
	B2 := LIMIT2
	B3 := LIMIT3
	RHO21 := SIGMA1
	RHO31 := SIGMA2
	RHO32 := SIGMA3
	if math.abs(B2) > math.max(math.abs(b1), math.abs(B3))
		b1 := B2
		B2 := LIMIT1
		RHO31 := RHO32
		RHO32 := SIGMA2
	else if math.abs(B3) > math.max(math.abs(b1), math.abs(B2)) 
		b1 := B3
		B3 := LIMIT1
		RHO21 := RHO32
		RHO32 := SIGMA1
	
	TAIL := false
	
	if b1 > 0 
		TAIL := true
		b1 := -b1
		RHO21 := -RHO21
		RHO31 := -RHO31
	
	if b1 > XCUT 
		if 2 * math.abs(RHO21) < 1 
			SQ21 := math.sqrt(1 - RHO21 * RHO21)
		else
			SQ21 := math.sqrt((1 - RHO21) * (1 + RHO21))
		
		if 2 * math.abs(RHO31) < 1 
			SQ31 := math.sqrt(1 - RHO31 * RHO31)
		else
			SQ31 := math.sqrt((1 - RHO31) * (1 + RHO31))
		
		rho := (RHO32 - RHO21 * RHO31) / (SQ21 * SQ31)
		b2p := B2 / SQ21
		RHO21 := RHO21 / SQ21
		b3p := B3 / SQ31
		RHO31 := RHO31 / SQ31
		CTND := ADONED(XCUT, b1, EPS, b2p, b3p, RHO21, RHO31, rho) / SQTWPI
	else
		CTND := 0
	
	if TAIL == true 
		CTND := cbnd.CBND3(B2, B3, RHO32) - CTND
	CTND

// example usage
float out = 0
if barstate.islast
	out := CTND(0., 0., 0., 0., 0., 0.)
plot(out)
