Bat optimization algorithm with special levy jump implementation on

Searching an optimized solution of a function. on formula with 2 parameters.

The js fiddle

http://jsfiddle.net/sebadiaz/6YrQT/99/

The source code

function ln(val) {

    return Math.log(val) / Math.log(Math.E);

}

function pollarMethod() {

    do {

        U1 = Math.random();

        U2 = Math.random();

        V1 = 2 * U1 - 1.0;

        V2 = 2 * U2 - 1.0;

        S = V1 * V1 + V2 * V2;

    }

    while (S >= 1);

    X = V1 * Math.sqrt((-2 * ln(S)) / S);

    return X;

}

function poissonProcess(lambda, T) {

    var t = 0;

    var k = 0;

    var S = new Array();

    while (t <= T) {

        var r = Math.random();

        t = t - ln(r) / lambda;

        if (t > T) return S;

        S[k] = t;

        k = k + 1;

    }

    return S;

}

function levyJump(t, b, c, lambda) {

    var L = b * t;

    L = L + Math.sqrt(c) * pollarMethod();

    var nbLoop = poissonProcess(lambda, t);

    for (var i = 0; i < nbLoop; i++) {

        L += Math.random();

    }

    return L;

}

function fOptimize(x) {

    return x[0] * x[0] + x[1] * x[1]*x[1];

}

function bat(fOptimize, nbx, nbbat, maxiter, fmin1, fmax1, fmin2, fmax2) {

    var xxh = new Array();

    var vv = new Array();

    var xx = new Array();

    var ff1 = new Array();

    var ff2 = new Array();

    var rr = new Array();

    var rrOrig = new Array();

    var AA = new Array();

    var globalSolution = 99999999;

    var XSTAR = new Array();

    for (var j = 0; j < nbx; j++)

    XSTAR[j] = Math.random();

    for (var i = 0; i < nbbat; i++) {

        xx[i] = new Array();

        vv[i] = new Array();

        ff1[i] = new Array();

        ff2[i] = new Array();

        rr[i] = Math.random();

        rrOrig[i] = rr[i];

        AA[i] = Math.random();

        for (var j = 0; j < nbx; j++) {

            xx[i][j] = Math.random() ;

            vv[i][j] = Math.random() ;

            ff1[i][j] = Math.random() ;

            ff2[i][j] = Math.random() ;

        }

    }

    var t = 0;

    while (t < maxiter) {

        var actualSolution = 99999999;

        var xstar = new Array();

        for (var j = 0; j < nbx; j++)

        xstar[j] = Math.random();

        for (var i = 0; i < nbbat; i++) {

            for (var j = 0; j < nbx; j++) {

                ff1[i][j] = fmin1 + (fmax1 - fmin1) * Math.random()+levyJump(t, 0.5, 0.9, 2.0);

                vv[i][j] = vv[i][j] + xx[i][j] - XSTAR[j] * ff1[i][j];

                xx[i][j] = xx[i][j] + vv[i][j];

                /*

                ff1[i][j] = ((fmin1 - fmax1) * t / 2 + fmax1) * Math.random();

                ff2[i][j] = ((fmin2 - fmax2) * t / 2 + fmax2) * Math.random();

                var tempx = XSTAR[j] + ff1[i][j] * xx[i][j] + ff2[i][j] * xx[i][j];

                xx[i][j] = tempx + (Math.random() - 0.5) * levyJump(1, 0.9, 20.0, 1.0);

                */

            }

            var res = fOptimize(xx[i]);

            if (res < actualSolution) {

                xstar = xx[i];

                res = actualSolution;

            }

        }

        //Local solution search

        for (var i = 0; i < nbbat; i++) {

            if (Math.random() > rr[i]) {

                for (var j = 0; j < nbbat; j++) {

                    var newX = new Array();

                    for (var k = 0; k < nbx; k++) {

                        newX = xx[j][k] + (Math.random() * 2 - 1.0) * AA[i];

                    }

                    var res = fOptimize(newX);

                    if (res < actualSolution) {

                        xstar = newX;

                        res = actualSolution;

                    }

                }

            }

        }

        //Global solution search

        for (var i = 0; i < nbbat; i++) {

            if (Math.random() > rr[i]) {

                for (var j = 0; j < nbbat; j++) {

                    var newX = new Array();

                    for (var k = 0; k < nbx; k++) {

                        newX = xx[j][k] + (Math.random() * 2 - 1.0) * rr[i];

                    }

                    var res = fOptimize(newX);

                    if (res < actualSolution) {

                        xstar = newX;

                        res = actualSolution;

                    }

                }

            }

        }

        if (actualSolution < globalSolution) {

            for (var j = 0; j < nbbat; j++) {

                AA[j] = 0.7 * AA[j];

                rr[j] = rrOrig * (1 - Math.exp(-0.7 * t));

            }

            actualSolution=globalSolution;

            XSTAR=xstar;

        }

        xxh[t] = xx;

        t = t + 1;

        

    }

    xxh[xxh.length] = XSTAR;

    return xxh;

}

var str = '<svg xmlns="http://www.w3.org/2000/svg" version="1.1">';

var previousvalue = 0;

for (var i = 0; i < 5; i++) {

    val = levyJump(i, 0.9, 20.0, 1.0);

    str += '<line x1="' + i * 10 + '" y1="' + previousvalue + '" x2="' + (i + 1) * 10 + '" y2="' + val * 10 + '" style="stroke:rgb(255,0,0);stroke-width:2" />';

    previousvalue = val * 10;

}

var previousX = 0;

var previousY = 0;

var str = '<svg xmlns="http://www.w3.org/2000/svg" version="1.1">';

var X = bat(fOptimize, 2, 50, 30, 0, 5, 0, 5);

var maxX=0;

var minX=9999999;

var maxY=0;

var minY=9999999;

for (var i = 0; i < X.length-1; i++) {

    var x = X[i][0][0];

    var y = X[i][0][1];

    if (x>maxX)

        maxX=x;

    if (y<minY)

        minY=y;

    if (y>maxY)

        maxY=y;

    if (x<minX)

        minX=x;

}

var maxouX=1;

var maxouY=1;

if (Math.abs(maxX)>maxouX)

    maxouX=Math.abs(maxX)/15;

if (Math.abs(minX)>maxouX)

    maxouX=Math.abs(minX)/15;

if (Math.abs(maxY)>maxouY)

    maxouY=Math.abs(maxY)/15;

if (Math.abs(minY)>maxouY)

    maxouY=Math.abs(minY)/15;

for (var j = 0; j < X[0].length; j++) {

    var colorq=j*5;

    var colorp=j*5;

    var colorr=j*5;

    if(Math.random()>0.5)

        colorr=0;

    if(Math.random()>0.5)

        colorp=0;

    if(Math.random()>0.5)

        colorq=0;

for (var i = 0; i < X.length-1; i++) {

    var x = Math.abs(X[i][j][0])/maxouX;

    var y = Math.abs(X[i][j][1])/maxouY;

   

    str += '<line x1="' + previousX + '" y1="' + previousY + '" x2="' + x * 10 + '" y2="' + y * 10 + '" style="stroke:rgb('+colorr+','+colorq+','+colorp+');stroke-width:2" />';

    previousX = x * 10;

    previousY = y * 10;

}

}

$("#graph").append(str + '</svg>');

$("#graph").append("Resultat=" + X[X.length - 1][0] + " " + X[X.length - 1][1]);

Levy Jump

try jsfiddle

function ln(val){
    return Math.log(val) / Math.log(Math.E);
}

function pollarMethod()
{
do{
U1=Math.random();
U2=Math.random();
V1=2*U1-1.0;
V2=2*U2-1.0;
S=V1*V1+V2*V2;
}
while(S>=1);

X=V1*Math.sqrt((-2*ln(S))/S);
return X;
}

function poissonProcess(lambda,T)
{
    var t=0;
    var k=0;
    var S=new Array();
    while(t<=T)
    {
        var r=Math.random();
        t=t-ln(r)/lambda
        if (t>T) return S;
        S[k]=t;
        k=k+1;
    }
    return S;
}

function levyJump(t,b,c,lambda)
{
    var L=b*t;
    L=L+Math.sqrt(c)*pollarMethod();
    var nbLoop=poissonProcess(lambda,t);
    for (var i=0;i<nbLoop;i++)
    {
        L+=Math.random();
    }
    return L;
}
var str='<svg xmlns="http://www.w3.org/2000/svg" version="1.1">';
var previousvalue=0;

for (var i = 0; i < 50; i++) {
    val=levyJump(i,0.9,20.0,1.0);
    str+='<line x1="'+i*10+'" y1="'+previousvalue+'" x2="'+(i+1)*10+'" y2="'+val*10+'" style="stroke:rgb(255,0,0);stroke-width:2" />';
    previousvalue=val*10;
   
}
$("#graph").append(str+'</svg>' );

Non Homogeneous Poisson Process

Try jsfiddle

function ln(val){
    return Math.log(val) / Math.log(Math.E);
}


function nonHomogeneousPoissonProcess(lambdaTArray,T)
{
    var lambda=0.0;
    for (var i = 0; i < lambdaTArray.length; i++) {
     
        if (lambdaTArray[i]>lambda)
            lambda=lambdaTArray[i]+1;
    }
 
    var t=0;
    var k=0;
    var S=new Array();
    while(t<=T)
    {
     
        var r=Math.random();
        t=t-ln(r)/lambda;
        if (t>T) return S;
        var s=Math.random();
        if (s<=lambdaTArray[parseInt(t, 10)]/lambda)
        {
            S[k]=t;
            k=k+1;
        }
    }
    return S;
}
var str='<svg xmlns="http://www.w3.org/2000/svg" version="1.1">';
var lambdaTArray=new Array();
var T=50
for (var i = 0; i < T; i++)
    lambdaTArray[i]=Math.random();
var Liste=nonHomogeneousPoissonProcess(lambdaTArray,T);
var previousvalue=0.0;
for (var i = 0; i < Liste.length; i++) {
 
    str+='<line x1="'+i*10+'" y1="'+previousvalue+'" x2="'+(i+1)*10+'" y2="'+Liste[i]*10+'" style="stroke:rgb(255,0,0);stroke-width:2" />';
    previousvalue=Liste[i]*10;
 
}
$("#graph").append(str+'</svg>' );

Poisson process

Try my jsFiddle

function ln(val){
    return Math.log(val) / Math.log(Math.E);
}


function poissonProcess(lambda,T)
{
    var t=0;
    var k=0;
    var S=new Array();
    while(t<=T)
    {
        var r=Math.random();
        t=t-ln(r)/lambda
        if (t>T) return S;
        S[k]=t;
        k=k+1;
    }
    return S;
}

Spread Option

function CallSpreadOption(S1,S2,X,r,sigma1,sigma2,rho,T)
{
 sigmaC1=sigma1*sigma1;
 sigmaC2=sigma2*sigma2;
 Sa=S1/(S2+X);
 nSig=Math.sqrt(sigmaC1+Math.pow(sigma2*S2/(S2+X),2)-2*rho*sigma1*sigma2*S2/(S2+X));
 d1=(ln(Sa)+0.5*nSig*nSig*T)/Math.sqrt(T)/nSig;
 d2=d1-Math.sqrt(T)*nSig;
 
 
 return (S2+X)*Math.exp(-r*T)*((Sa*NCDF2(d1))-NCDF2(d2));
}

function PutSpreadOption(S1,S2,X,r,sigma1,sigma2,rho,T)
{
 sigmaC1=sigma1*sigma1;
 sigmaC2=sigma2*sigma2;
 Sa=S1/(S2+X);
 nSig=Math.sqrt(sigmaC1+Math.pow(sigma2*S2/(S2+X),2)-2*rho*sigma1*sigma2*S2/(S2+X));
 d1=(ln(Sa)+0.5*nSig*nSig*T)/Math.sqrt(T)/nSig;
 d2=d1-Math.sqrt(T)*nSig;
 
 
 return (S2+X)*Math.exp(-r*T)*((NCDF2(-d2)-Sa*NCDF2(-d1)));
}

Caplet/Floorlet option

function CapletOption(F,X,Notional,d,basis,r,sigma,T1,T2)
{


 sigma2=sigma*sigma;
 d1=(ln(F/X)+(sigma2/2/(T2-T1)))/Math.sqrt(T2-T1)/sigma;
 d2=d1-Math.sqrt(T2-T1)*sigma;
 
 return (Notional*d/basis)/(1+F*d/basis)*Math.exp(-r*(T2-T1))*(F*NCDF2(d1)-X*NCDF2(d2));
}


function FloorletOption(F,X,Notional,d,basis,r,sigma,T1,T2)
{


 sigma2=sigma*sigma;
 d1=(ln(F/X)+(sigma2/2/(T2-T1)))/Math.sqrt(T2-T1)/sigma;
 d2=d1-Math.sqrt(T2-T1)*sigma;
 
 return (Notional*d/basis)/(1+F*d/basis)*Math.exp(-r*(T2-T1))*(X*NCDF2(-d2)-F*NCDF2(-d1));
}

Compound Option Sharp Monte Carlo

function calculateMCCompoundOptionSSharp(oStart,oEnd,precision,BASE,X1,X2,r,sigma,q,T1,T2)
{
 xn=oStart;
 minValue=xn;
 minPrecision=888888888888;
 fxprev=0;
 counter=0;
 xn1=xn;
 xn2=0;
 basee=BASE*100;
 beginBase=0;
 xnOld=0;
 fxmOld=0
 while (counter<oEnd)
 {
  
  xn=Math.random()*basee+beginBase-0.5*basee;
  if (xn>0)
  {
   fxm1=CallBlackScholes(xn,X1,r,sigma,q,T2-T1,0)-X2;
   
   va=Math.abs(fxm1);
   if (va <= precision)
   {
    
    return xn;
   }
   if (minPrecision>va)
   {
    minPrecision=va;
    minValue=xn;
   }
   if (counter>0)
   {
    if (Math.abs(fxm1)<Math.abs(fxmOld))
    {
     beginBase=xn;
     basee=beginBase*0.1;
    } 
    else
    {
     beginBase=xnOld;
     basee=beginBase*0.1;
    }
    
   }
   xnOld=xn;
   fxmOld=fxm1;
  }
  counter++;
  
 }
 
 
 return minValue;
 
}

function calculateLinCompoundOptionSSharp(oStart,oEnd,precision,BASE,X1,X2,r,sigma,q,T1,T2)
{
 xn=oStart;
 minValue=BASE;
 
 fxprev=0;
 counterade=0;
 xn1=xn;
 xn2=0;
 basee=BASE*100;
 beginBase=0;
 xnOld=0;
 fxmOld=0;
 maxvalue=BASE;
 if (maxvalue<X1)
  maxvalue=X1;
 if (maxvalue<X2)
  maxvalue=X2;
 mratio=5*maxvalue/oEnd;
 minPrecision=maxvalue*100000;
 while (counterade<oEnd)
 {
  po=counterade;
  xn=po*mratio;
  if (xn>0)
  {
   fxm1=CallBlackScholes(xn,X1,r,sigma,q,T2-T1,0)-X2;
   if (isNaN(fxm1))
   {
    counterade=counterade+1;
    continue;
   }
   va=Math.abs(fxm1);
   if (va <= precision)
   {
    
    return xn;
   }
   if (!isNaN(va)&&minPrecision>va )
   {
    
    minPrecision=va;
    minValue=xn;
   }
   
  }
  counterade=counterade+1;
  
 }
 
 
 return minValue;
 
}

function e(val){
 return Math.exp(val);
}
function exp(val){
 return e(val);
}
function N2(a,b,r)
{
 return CumulativeBivariateNormal(a,b,r);
}
function N(a)
{
 return NCDF2(a);
}

function CallCallCompoundOption(S,X1,X2,r,sigma,q,T1,T2)
{
 sigma2=sigma*sigma;
 //search bossvalue
 Base=S;
 if (X1>Base)
 Base=X1;
 if (X2>Base)
 Base=X2;
 Sb=calculateLinCompoundOptionSSharp(0,9999,0.1,S,X1,X2,r,sigma,q,T1,T2);
 a1=(ln(S/Sb)+(r-q+0.5*sigma2)*T1)/(sigma*Math.sqrt(T1));
 a2=a1-sigma*Math.sqrt(T1);
 b1=(ln(S/X2)+(r-q+0.5*sigma2)*T2)/(sigma*Math.sqrt(T2));
 b2=b1-(sigma*Math.sqrt(T2));
 return S*e(-q*T2)*N2(a1,b1,Math.sqrt(T1/T2))-X2*e(-r*T2)*N2(a2,b2,Math.sqrt(T1/T2))-X1*e(-r*T1)*N(a2);
}

function CallPutCompoundOption(S,X1,X2,r,sigma,q,T1,T2)
{
 sigma2=sigma*sigma;
 //search bossvalue
 Base=S;
 if (X1>Base)
 Base=X1;
 if (X2>Base)
 Base=X2;
 Sb=calculateLinCompoundOptionSSharp(0,9999,0.1,S,X1,X2,r,sigma,q,T1,T2);
 
 a1=(ln(S/Sb)+(r-q+0.5*sigma2)*T1)/(sigma*Math.sqrt(T1));
 a2=a1-sigma*Math.sqrt(T1);
 b1=(ln(S/X2)+(r-q+0.5*sigma2)*T2)/(sigma*Math.sqrt(T2));
 b2=b1-(sigma*Math.sqrt(T2));
 return S*e(-q*T2)*N2(a1,-b1,-Math.sqrt(T1/T2))-X2*e(-r*T2)*N2(a2,-b2,-Math.sqrt(T1/T2))+X1*e(-r*T1)*N(a2);
}
function PutCallCompoundOption(S,X1,X2,r,sigma,q,T1,T2)
{
 sigma2=sigma*sigma;
 //search bossvalue
 Base=S;
 if (X1>Base)
 Base=X1;
 if (X2>Base)
 Base=X2;
 Sb=calculateLinCompoundOptionSSharp(0,9999,0.1,S,X1,X2,r,sigma,q,T1,T2);
 
 a1=(ln(S/Sb)+(r-q+0.5*sigma2)*T1)/(sigma*Math.sqrt(T1));
 a2=a1-sigma*Math.sqrt(T1);
 b1=(ln(S/X2)+(r-q+0.5*sigma2)*T2)/(sigma*Math.sqrt(T2));
 b2=b1-(sigma*Math.sqrt(T2));
 return X2*e(-r*T2)*N2(-a2,b2,-Math.sqrt(T1/T2))-S*e(-q*T2)*N2(-a1,b1,-Math.sqrt(T1/T2))+X1*e(-r*T1)*N(-a2);
}
function PutPutCompoundOption(S,X1,X2,r,sigma,q,T1,T2)
{
 sigma2=sigma*sigma;
 //search bossvalue
 Base=S;
 if (X1>Base)
 Base=X1;
 if (X2>Base)
 Base=X2;
 Sb=calculateLinCompoundOptionSSharp(0,9999,0.1,S,X1,X2,r,sigma,q,T1,T2);
 
 a1=(ln(S/Sb)+(r-q+0.5*sigma2)*T1)/(sigma*Math.sqrt(T1));
 a2=a1-sigma*Math.sqrt(T1);
 b1=(ln(S/X2)+(r-q+0.5*sigma2)*T2)/(sigma*Math.sqrt(T2));
 b2=b1-(sigma*Math.sqrt(T2));
 return X2*e(-r*T2)*N2(-a2,-b2,Math.sqrt(T1/T2))-S*e(-q*T2)*N2(-a1,-b1,Math.sqrt(T1/T2))-X1*e(-r*T1)*N(-a2);
}

Bivariate Normal cumulative distribution function

function BiVariateNCDF(a , b , rho )
{
 aarray=new Array;
 barray=new Array;
    rho1=0.0;
 rho2=0.0;
 delta=0.0;
    ap=0.0;
 bp=0.0;
    Sum=0;
    i=0;
 j=0;
 Pi=0;
    denum=0.0;
 aarray[0]=0.24840615;
 aarray[1]=0.39233107;
 aarray[2]=0.21141819;
 aarray[3]=0.03324666;
 aarray[4]=0.00082485334;
 barray[0]=0.10024215;
 barray[1]=0.48281397;
 barray[2]=1.0609498;
 barray[3]=1.7797294;
 barray[4]=2.6697604;
 Pi = 3.14159265358979;
    
    ap = a / Math.sqrt(2 * (1 - rho * rho));
    bp = b / Math.sqrt(2 * (1 - rho * rho ));
 BivN=0.0;
 if (a <= 0 && b <= 0 && rho <= 0)
 {
  Sum = 0;
        for (i = 0; i< 5;i++)
   for (j = 0; j< 5;j++)
                Sum = Sum + aarray[i] * aarray[j] * f(barray[i], barray[j], ap, bp, rho);
        BivN = Sum * Math.sqrt(1 - rho * rho) / Pi;
 }
 else if( a <= 0 && b >= 0 && rho >= 0)
  BivN = NCDF2(a) - BiVariateNCDF(a, -b, -rho);
 else if(a >= 0 && b <= 0 && rho >= 0 )
 {
  BivN = NCDF2(b) - BiVariateNCDF(-a, b, -rho);
 }
 else if(a >= 0 && b >= 0 && rho <= 0 )
  BivN = NCDF2(a) + NCDF2(b) - 1 + BiVariateNCDF(-a, -b, rho);
 else if(a * b * rho >= 0)
 {
  denum = Math.sqrt(Math.abs(a*a-2*rho*a*b+b*b));
        rho1 = (rho * a - b) * Sgn(a) / denum;
        rho2 = (rho * b - a) * Sgn(b) / denum;
        delta = (1 - Sgn(a) * Sgn(b)) / 4;
  
        BivN = BiVariateNCDF(a, 0, rho1) + BiVariateNCDF(b, 0, rho2) - delta;
 }
 
 return BivN;
}


function CumulativeBivariateNormal(a, b, rho)
 {
    aarray=new Array;
 aarray[0]= 0.24840615;
 aarray[1]=0.39233107;
 aarray[2]=0.21141819;
 aarray[3]=0.03324666;
 aarray[4]=0.00082485334;
    barray=new Array;
 barray[0]=0.10024215;
 barray[1]=0.48281397;
 barray[2]=1.0609498;
 barray[3]=1.7797294;
 barray[4]=2.6697604;
  //   static Real aarray[4]= {0.325303, 0.4211071, 0.1334425, 0.006374323};
  //static Real barray[4]={0.1337764, 0.6243247, 1.3425378, 2.2626645};
    rho1=0;  rho2=0;  delta =0;
 ap=0;  bp =0;
    Sum=0; Pi=0;
    sgnA=0;sgnB=0;
    i=0; j=0 ;
    BivN=0; denum=0;

    sgnA=Sgn(a);
    sgnB=Sgn(b);
    
    Pi = 3.14159265358979;
    if(Math.pow(rho,2)==1)
        rho=Sgn(rho)*0.99999;
    ap = a / Math.sqrt(2 * (1 - Math.pow(rho,2)));
    bp = b / Math.sqrt(2 * (1 - Math.pow(rho,2)));
     
    if( (a <= 0) && (b <= 0) && (rho <= 0) ){
        Sum = 0;
        for (i = 0 ;i< 5;i++)
           for (j = 0 ;j< 5;j++)
                Sum +=  + aarray[i] * aarray[j] * SubFunctionForBivariateNormal(barray[i], barray[j], ap, bp, rho);
        BivN = Sum * Math.sqrt(1 - Math.pow(rho,2)) / Pi;
    }  
    else if ((a <= 0 )&& (b >= 0) && (rho >= 0)) 
        BivN = NCDF2(a) - CumulativeBivariateNormal(a, -b, -rho); 
    else if ((a >= 0) && (b <= 0) && (rho >= 0)) 
        BivN = NCDF2(b) - CumulativeBivariateNormal(-a, b, -rho);
    else if ((a >= 0) && (b >= 0) && (rho <= 0))
        BivN = NCDF2(a) + NCDF2(b) - 1 + CumulativeBivariateNormal(-a, -b, rho);
    else if (a * b * rho >= 0) {
        denum = Math.sqrt(Math.pow(a,2) - 2 * rho * a * b + Math.pow(b,2));
        rho1 = (rho * a - b) * sgnA / denum;
        rho2 = (rho * b - a) * sgnB / denum;
        delta = (1 - sgnA*sgnB) / 4;
        BivN = CumulativeBivariateNormal(a, 0, rho1) + CumulativeBivariateNormal(b, 0, rho2) - delta;
    }       
    return BivN;
}
 
function SubFunctionForBivariateNormal( X, y, ap, bp, rho){
    r = ap * (2 * X - ap) + bp * (2 * y - bp) + 2 * rho * (X - ap) * (y - bp);
   return Math.exp(r);
}

Asian Option Arithmetic Levy Form (not stable)

function CallAsianOptionArithmeticLevyForm(S,X,r,sigma,q,T,t)
{
 sigmaC=sigma*sigma;
 M=(2*S*S/(r-q+sigmaC)*Math.exp((2*(r-q)+sigmaC)*(T-t)-1)/(2*(r-q)+sigmaC))-(Math.exp((r-q)*(T-t))-1)/(r-q);
 L=M/T/T;
 Sz=S/(r-q)/T*(Math.exp(-q*(T-t))-Math.exp(-r*(T-t)));
 K=ln(L)-2*(r*(T-t)+ln(Sz));
 Xz=X-S*((T-T-t)/T);
 d1=1/Math.sqrt(K)*(ln(L)*0.5-ln(Xz));
 d2=d1-Math.sqrt(K);
 
 return Sz*NCDF2(d1)-Xz*Math.exp(-r*(T-t))*NCDF2(d2);
}

function PutAsianOptionArithmeticLevyForm(S,X,r,sigma,q,T,t)
{
 sigmaC=sigma*sigma;
 M=(2*S*S/(r-q+sigmaC)*Math.exp((2*(r-q)+sigmaC)*(T-t)-1)/(2*(r-q)+sigmaC))-(Math.exp((r-q)*(T-t))-1)/(r-q);
 L=M/T/T;
 
 Sz=S/(r-q)/T*(Math.exp(-q*(T-t))-Math.exp(-r*(T-t)));

 K=ln(L)-2*(r*(T-t)+ln(Sz));
 
 Xz=X-S*((T-T-t)/T);
 
 d1=1/Math.sqrt(K)*(ln(L)/2-ln(Xz));
 d2=d1-Math.sqrt(K);
 
 return CallAsianOptionArithmeticLevyForm(S,X,r,sigma,q,T,t)-Sz+Xz*Math.exp(-r*(T-t));
}

Asian Option Arithmetic Turbunn Wakwman Form

function CallAsianOptionArithmeticTurbunnWakemanForm(S,K,r,sigma,q,T,t)
{
 sigmaC=sigma*sigma;
 M1=(Math.exp((r-q)*T)-Math.exp((r-q)*t))/((r-q)*(T-t));
 
 M2p1=2*Math.exp((2*(r-q)+sigmaC)*T)/((r-q+sigmaC)*(2*r-2*q+sigmaC)*(T-t)*(T-t));
 M2p2=2*Math.exp((2*(r-q)+sigmaC)*t)/((r-q)*(T-t)*(T-t))*((1/(2*(r-q)+sigmaC))-(Math.exp(((r-q))*T)/(r-q+sigmaC)));
 M2=M2p1+M2p2;
 b=ln(M1)/T;
 sigmaA=Math.sqrt((ln(M2)/T)-2*b);
 
 d1=(ln(S/K)+(b+0.5*sigmaC)*(T-t))/(sigmaA*Math.sqrt(T-t));
 d2=d1-(sigmaA*Math.sqrt(T-t));
 return S*Math.exp((b-r)*(T-t))*NCDF2(d1)-K*Math.exp(-r*(T-t))*NCDF2(d2);
}

function PutAsianOptionArithmeticTurbunnWakemanForm(S,K,r,sigma,q,T,t)
{
 sigmaC=sigma*sigma;
 M1=(Math.exp((r-q)*T)-Math.exp((r-q)*t))/((r-q)*(T-t));
 
 M2p1=2*Math.exp((2*(r-q)+sigmaC)*T)/((r-q+sigmaC)*(2*r-2*q+sigmaC)*(T-t)*(T-t));
 M2p2=(2*Math.exp((2*(r-q)+sigmaC)*t)/((r-q)*(T-t)*(T-t)))*((1/(2*(r-q)+sigmaC))-(Math.exp(((r-q))*T)/(r-q+sigmaC)));
 M2=M2p1+M2p2;
 b=ln(M1)/T;
 sigmaA=Math.sqrt((ln(M2)/T)-2*b);
 
 d1=(ln(S/K)+(b+0.5*sigmaC)*(T-t))/(sigmaA*Math.sqrt(T-t));
 d2=d1-(sigmaA*Math.sqrt(T-t));
 return K*Math.exp(-r*(T-t))*NCDF2(-d2)-S*Math.exp((b-r)*(T-t))*NCDF2(-d1);
}

Asian Option Geometric Closed Form

function CallAsianOptionGeometricClosedForm(S,K,r,sigma,q,T,t)
{
 sigmaA=sigma/Math.sqrt(3);
 b=0.5*(r-q-sigma*sigma/6);
 d1=(ln(S/K)+(b+0.5*sigmaA*sigmaA)*(T-t))/(sigmaA*Math.sqrt(T-t));
 d2=(ln(S/K)+(b-0.5*sigmaA*sigmaA)*(T-t))/(sigmaA*Math.sqrt(T-t));
 return S*Math.exp((b-r)*(T-t))*NCDF2(d1)-K*Math.exp(-r*(T-t))*NCDF2(d2);
}

function PutAsianOptionGeometricClosedForm(S,K,r,sigma,q,T,t)
{
 sigmaA=sigma/Math.sqrt(3);
 b=0.5*(r-q-sigma*sigma/6);
 d1=(ln(S/K)+(b+0.5*sigmaA*sigmaA)*(T-t))/(sigmaA*Math.sqrt(T-t));
 d2=(ln(S/K)+(b-0.5*sigmaA*sigmaA)*(T-t))/(sigmaA*Math.sqrt(T-t));
 return K*Math.exp(-r*(T-t))*NCDF2(-d2)-S*Math.exp((b-r)*(T-t))*NCDF2(-d1);
}

Barrier Down/UP Out/IN Call PUt

function DownAndOutCall(S,K,r,sigma,q,barrier,T,t)
{
 
 alpha=0.5*(1-(r-q)*0.5/Math.sqrt(sigma*sigma));
 
 
 if (barrier<=K)
 {
  value=CallBlackScholes(S,K,r,sigma,q,T,t)-Math.pow(S/barrier,2*alpha)*CallBlackScholes(barrier*barrier/S,K,r,sigma,q,T,t);
 }
 else
 {
  value=CallBlackScholes(S,barrier,r,sigma,q,T,t)+(barrier-K)*CallBlackScholesDigit(S,barrier,r,sigma,q,T,t)-Math.pow(S/barrier,2*alpha)*(CallBlackScholes(barrier*barrier/S,barrier,r,sigma,q,T,t)+(barrier-K)*CallBlackScholesDigit(barrier*barrier/S,barrier,r,sigma,q,T,t));
 }
 return value;
}

function DownAndOutPut(S,K,r,sigma,q,barrier,T,t)
{
 alpha=0.5*(1-(r-q)*0.5/Math.sqrt(sigma*sigma));
 
 value=PutBlackScholes(S,K,r,sigma,q,T,t)-PutBlackScholes(S,barrier,r,sigma,q,T,t)-(K-barrier)*PutBlackScholesDigit(S,barrier,r,sigma,q,T,t)-Math.pow(S/barrier,2*alpha)*(PutBlackScholes(barrier*barrier/S,K,r,sigma,q,T,t)-PutBlackScholes(barrier*barrier/S,barrier,r,sigma,q,T,t)+(K-barrier)*PutBlackScholesDigit(barrier*barrier/S,barrier,r,sigma,q,T,t));
 return value;
}

function DownAndInCall(S,K,r,sigma,q,barrier,T,t)
{
 return CallBlackScholes(S,K,r,sigma,q,T,t)-DownAndOutCall(S,K,r,sigma,q,barrier,T,t);
}

function DownAndInPut(S,K,r,sigma,q,barrier,T,t)
{
 
 return PutBlackScholes(S,K,r,sigma,q,T,t)-DownAndOutPut(S,K,r,sigma,q,barrier,T,t);
}

function UpAndOutCall(S,K,r,sigma,q,barrier,T,t)
{
 
 alpha=0.5*(1-(r-q)*0.5/Math.sqrt(sigma*sigma));
 
 value=CallBlackScholes(S,K,r,sigma,q,T,t)-CallBlackScholes(S,barrier,r,sigma,q,T,t)-(barrier-K)*CallBlackScholesDigit(S,barrier,r,sigma,q,T,t)-Math.pow(S/barrier,2*alpha)*(CallBlackScholes(barrier*barrier/S,K,r,sigma,q,T,t)-CallBlackScholes(barrier*barrier/S,barrier,r,sigma,q,T,t)+(barrier-K)*CallBlackScholesDigit(barrier*barrier/S,barrier,r,sigma,q,T,t));
 
 return value;
}

function UpAndOutPut(S,K,r,sigma,q,barrier,T,t)
{
 
 alpha=0.5*(1-(r-q)*0.5/Math.sqrt(sigma*sigma));
 if (barrier>K)
  value=PutBlackScholes(S,K,r,sigma,q,T,t)-Math.pow(S/barrier,2*alpha)*PutBlackScholes(barrier*barrier/S,K,r,sigma,q,T,t);
 else
  value=PutBlackScholes(S,K,r,sigma,q,T,t)+(K-barrier)*PutBlackScholesDigit(S,barrier,r,sigma,q,T,t)-Math.pow(S/barrier,2*alpha)*(PutBlackScholes(barrier*barrier/S,barrier,r,sigma,q,T,t)+(K-barrier)*PutBlackScholesDigit(barrier*barrier/S,barrier,r,sigma,q,T,t));
 
 return value;
}

function UpAndInCall(S,K,r,sigma,q,barrier,T,t)
{
 return CallBlackScholes(S,K,r,sigma,q,T,t)-UpAndOutCall(S,K,r,sigma,q,barrier,T,t);
}

function UpAndInPut(S,K,r,sigma,q,barrier,T,t)
{
 
 return PutBlackScholes(S,K,r,sigma,q,T,t)-UpAndOutPut(S,K,r,sigma,q,barrier,T,t);
}

Digital option with Black and Scholle

unction CallBlackScholesDigit(S,K,r,sigma,q,T,t)
{
 
 return (NCDF(D2BlackScholes(S,K,r,sigma,q,T,t),0,1)*K*Math.exp(r*(T-t)*-1.0));
}

function PutBlackScholesDigit(S,K,r,sigma,q,T,t)
{
 
 return (NCDF(-D2BlackScholes(S,K,r,sigma,q,T,t),0,1)*K*Math.exp(r*(T-t)*-1.0));
}

the Garman–Kohlhagen model

function D1GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t)
{
 sigma2=sigma*sigma;
 TT=T-t;
 F=S*Math.exp((rd-q)*TT);
 r1=ln(S/K)+TT*(rd-rf+sigma2*0.5);
 r2=sigma*Math.sqrt(TT);
 return r1/r2;
 //value=1/(sigma*Math.sqrt(T-t))*
}
function D2GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t)
{
 return D1GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t)-sigma*Math.sqrt(T-t);
}
function CallGarmanKohlhagen(S,K,rd,rf,sigma,q,T,t)
{
 
 return NCDF(D1GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t),0,1)*S*Math.exp(-rf*(T-t))-(NCDF(D2GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t),0,1)*K*Math.exp(rd*(T-t)*-1.0));
}

function PutGarmanKohlhagen(S,K,rd,rf,sigma,q,T,t)
{
 
 return (NCDF(-D2GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t),0,1)*K*Math.exp(rd*(T-t)*-1.0))-NCDF(-D1GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t),0,1)*S*Math.exp(-rf*(T-t));
}

function CallGarmanKohlhagenDelta(S,K,rd,rf,sigma,q,T,t)
{
 return NCDF(D1GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t),0,1);
}
function PutGarmanKohlhagenDelta(S,K,rd,rf,sigma,q,T,t)
{
 return NCDF(D1GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t),0,1)-1;
}
function CallGarmanKohlhagenGamma(S,K,rd,rf,sigma,q,T,t)
{
 return NormalDens(D1GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t))/(S*sigma*Math.sqrt(T-t));
}
function PutGarmanKohlhagenGamma(S,K,rd,rf,sigma,q,T,t)
{
 return CallGarmanKohlhagenGamma(S,K,rd,rf,sigma,q,T,t);
}
function CallGarmanKohlhagenVega(S,K,rd,rf,sigma,q,T,t)
{
 return NormalDens(D1GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t))*(S*Math.sqrt(T-t));
}
function PutGarmanKohlhagenVega(S,K,rd,rf,sigma,q,T,t)
{
 return CallGarmanKohlhagenVega(S,K,rd,rf,sigma,q,T,t);
}
function CallGarmanKohlhagenTheta(S,K,rd,rf,sigma,q,T,t)
{
 return (S*NormalDens(D1GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t))*sigma/(2*Math.sqrt(T-t)))-(rd*K*Math.exp(-rd*(T-t))*NCDF(D2GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t),0,1));
}
function PutGarmanKohlhagenTheta(S,K,rd,rf,sigma,q,T,t)
{
 return (S*NormalDens(D1GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t))*sigma/(2*Math.sqrt(T-t)))+(rd*K*Math.exp(-rd*(T-t))*NCDF(-D2GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t),0,1));
}
function CallGarmanKohlhagenRho(S,K,rd,rf,sigma,q,T,t)
{
 return K*(T-t)*Math.exp(-rd*(T-t))*NCDF(D2GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t),0,1);
}
function PutGarmanKohlhagenRho(S,K,rd,rf,sigma,q,T,t)
{
 return -K*(T-t)*Math.exp(-rd*(T-t))*NCDF(-D2GarmanKohlhagen(S,K,rd,rf,sigma,q,T,t),0,1);
}

Monte carlo option pricing

function pollarMethod()
{
 do{
  U1=Math.random();
  U2=Math.random();
  V1=2*U1-1.0;
  V2=2*U2-1.0;
  S=V1*V1+V2*V2;
 }
 while(S>=1);
 
 X=V1*Math.sqrt((-2*ln(S))/S);
 return X;
}

function callMontecarloEuropean(time,S,K,r,sigma,q,nbSim)
{
 
 var res=new Array;
 R=(r-q-0.5*sigma*sigma)*time;
 SD=sigma*Math.sqrt(time);
 sumPayof=0;
 for (i=0;i<nbSim;i++)
 {
  ST=S*Math.exp(R+SD*pollarMethod());
  res[i]=ST;
  if ((ST-K)>=0)
   sumPayof=sumPayof+(ST-K);
 }
 //Discounting
 res[nbSim]=Math.exp(-r*time)*(sumPayof/nbSim);
 return res;

}

function putMontecarloEuropean(time,S,K,r,sigma,q,nbSim)
{
 var res=new Array;
 R=(r-q-0.5*sigma*sigma)*time;
 SD=sigma*Math.sqrt(time);
 sumPayof=0;
 for (i=0;i<nbSim;i++)
 {
  ST=S*Math.exp(R+SD*pollarMethod());
  res[i]=ST;
  if ((K-ST)>=0)
   sumPayof=sumPayof+(K-ST);
 }
 //Discounting
 res[nbSim]=Math.exp(-r*time)*(sumPayof/nbSim);
 return res;

}

American trinomial pricing

function AmericanTrinomialPrice(T, stock, K, r, sigma, q, n,callput)
{
 N=n; 

 rate = r;
 
 
 div =q; 

 dt = T/N; 

 x = sigma*Math.sqrt(3.0*dt); 
 v = rate - div - sigma*sigma*0.5; 
 dis = Math.exp(-rate*dt);  
 edx = Math.exp(x);   


 pu = 0.5*(dt*(sigma*sigma + v*v*dt)/x/x + (v*dt/x));  //up probability
 pm = 1.0 - dt*(sigma*sigma + v*v*dt)/x/x;  //middle probability
 pd = 0.5*(dt*(sigma*sigma + v*v*dt)/x/x - (v*dt/x));  //down probability

 C=new Array;
 S=new Array;
 for(i=N;i>=0;i--)
 {
  S[i]=new Array;
 }
 for(i=N;i>=0;i--)
  {

   S[i][-i] = stock*Math.exp(-i*x);      for(j=-i+1;j<=i;j++)
    {
     S[i][j]=S[i][j-1]*edx;  
    }
  }

 for(i=N-1;i>=0;i--)
  C[i]=new Array;
 C[N]=new Array;
 for(j=-N;j<=N;j++)
 {
  if (callput=="c")
  {
   C[N][j] = max(S[N][j]-K,0.0);
  }
  else
  {
   C[N][j] = max(K-S[N][j],0.0);
  }
 }


 for(i=N-1;i>=0;i--)
  {
   for( j=-i;j<=i;j++)
    {
     C[i][j] = max(dis*(pu*C[i+1][j+1] + pm*C[i+1][j] + pd*C[i+1][j-1]),K-S[i][j]);
    }
  }

 return C;
}

European trinomial pricing

function EuropeanTrinomialPrice(T, stock, K, r, sigma, q, n,callput)
{
 N=n;  
 rate = r;
 div =q;
 dt = T/N;  

 x = sigma*Math.sqrt(3.0*dt); 
 v = rate - div - sigma*sigma*0.5;
 dis = Math.exp(-rate*dt); 
 edx = Math.exp(x);  


 pu = 0.5*(dt*(sigma*sigma + v*v*dt)/x/x + (v*dt/x));  //up probability
 pm = 1.0 - dt*(sigma*sigma + v*v*dt)/x/x;  //middle probability
 pd = 0.5*(dt*(sigma*sigma + v*v*dt)/x/x - (v*dt/x));  //down probability


 C=new Array;
 S=new Array;
 

 S[N]=new Array;

 for(i=-N;i<N+1;i++)
 {
  S[i]=new Array;
 }
 for(i=N;i>=0;i--)
  {

   S[i][-i] = stock*Math.exp(-i*x);   

   for(j=-i+1;j<=i;j++)
    {
     S[i][j]=S[i][j-1]*edx;  
    }
  }

 for(i=-N;i<N+1;i++)
 {
  C[i]=new Array;
 }
 C[N]=new Array;
 for(j=-N;j<=N;j++)
 {
  if (callput=="c")
  {
   C[N][j] = max(S[N][j]-K,0.0);
  }
  else
  {
   C[N][j] = max(K-S[N][j],0.0);
  }
 }

 for(i=N-1;i>=0;i--)
  {
   for( j=-i;j<=i;j++)
    {
     C[i][j] = dis*(pu*C[i+1][j+1] + pm*C[i+1][j] + pd*C[i+1][j-1]);
    }
  }


return C;
 

}

American binomial pricing

function AmericanBinomialPrice(T, stock, K, r, sigma, q, n,callput)
{
  N=n;  

 rate = r;
 
 dt = T/N;  
 u = Math.exp(sigma * Math.sqrt(dt));
 d = 1/u;
 dis = Math.exp(-rate*dt);  
 p = (Math.exp((rate-q)*dt)-d)/(u-d);  

 S=new Array;

 for(i=N;i>=0;i--)
  S[i]=new Array;
 for(i=N;i>=0;i--)
  {

   S[i][0] = stock*Math.pow(d,i); 

   for(j=1;j<=i;j++)
    {
     S[i][j]=S[i][j-1]*u/d;  
    }
  }



 C=new Array;
 for(j=0;j<=N;j++)
  C[j]=new Array;
 for(j=0;j<=N;j++)
 {
  if (callput=="c")
  {
   C[N][j] = max(S[N][j]-K,0.0);
  }else
  {
   C[N][j] = max(K-S[N][j],0.0);
  }
 }

 for(i=N-1;i>=0;i--)
  {
   for( j=0;j<=i;j++)
    {



     C[i][j] = max(dis*(p*C[i+1][j+1] + (1.0-p)*C[i+1][j]),K-S[i][j]);


    }
  }


 return C;
}

European binomial pricing

function EuropeanBinomialPrice(T, stock, K, r, sigma, q, n,callput)
{
 
 N=n;  

 rate = r;
 S =new Array;

 dt = T/N;  

 u = Math.exp(sigma * Math.sqrt(dt));
 d = 1/u;

 dis = Math.exp(-rate*dt);  

 p = (Math.exp((rate-q)*dt)-d)/(u-d);  

 C=new Array;

 st = 1.0;
 for ( i=0;i<N;i++)
  st = st*d;
 S[0] = stock*st; 
           
 udivd=u/d;
 for( i=1;i<N+1;i++)
 {
  S[i] = S[i-1]*udivd;
 }


 for( j=0;j<=N;j++)
 {
  C[j]=new Array;
 }
 for( j=0;j<=N;j++)
 {
  if (callput=="c")
  {
   C[N][j] = max(S[j]-K,0);
  }
  else
  {
   C[N][j] = max(K-S[j],0);
  }
 }


 for( i=N-1;i>=0;i--)
 {
  for( j=0;j<=i;j++)
  {
   C[i][j] = dis*(p*C[i+1][j+1] + (1.0-p)*C[i+1][j]);
  }
 }


return C;
 

}

Black and cholles

A list of javascript function to calculate Call and put, delta, gamma, vega ,Theta and rho

function D1BlackScholes(S,K,r,sigma,q,T,t)
{
 sigma2=sigma*sigma;
 TT=T-t;
 
 r1=ln(S/K)+TT*(r-q+sigma2*0.5);
 
 r2=sigma*Math.sqrt(TT);
 
 return r1/r2;
 //value=1/(sigma*Math.sqrt(T-t))*
}
function D2BlackScholes(S,K,r,sigma,q,T,t)
{
 return (ln(S/K)+(T-t)*(r-q-(sigma*sigma*0.5)))/(sigma*Math.sqrt(T-t));
}
function CallBlackScholes(S,K,r,sigma,q,T,t)
{
 
 return NCDF(D1BlackScholes(S,K,r,sigma,q,T,t),0,1)*S-(NCDF(D2BlackScholes(S,K,r,sigma,q,T,t),0,1)*K*Math.exp(-r*(T-t)));
}

function PutBlackScholes(S,K,r,sigma,q,T,t)
{
 
 return NCDF(-D2BlackScholes(S,K,r,sigma,q,T,t),0,1)*K*Math.exp(-r*(T-t))-NCDF(-D1BlackScholes(S,K,r,sigma,q,T,t),0,1)*S;
}

function CallBlackScholesDelta(S,K,r,sigma,q,T,t)
{
 return NCDF(D1BlackScholes(S,K,r,sigma,q,T,t),0,1);
}
function PutBlackScholesDelta(S,K,r,sigma,q,T,t)
{
 return NCDF(D1BlackScholes(S,K,r,sigma,q,T,t),0,1)-1;
}
function CallBlackScholesGamma(S,K,r,sigma,q,T,t)
{
 return NormalDens(D1BlackScholes(S,K,r,sigma,q,T,t))/(S*sigma*Math.sqrt(T-t));
}
function PutBlackScholesGamma(S,K,r,sigma,q,T,t)
{
 return CallBlackScholesGamma(S,K,r,sigma,q,T,t);
}
function CallBlackScholesVega(S,K,r,sigma,q,T,t)
{
 return NormalDens(D1BlackScholes(S,K,r,sigma,q,T,t))*(S*Math.sqrt(T-t));
}
function PutBlackScholesVega(S,K,r,sigma,q,T,t)
{
 return CallBlackScholesVega(S,K,r,sigma,q,T,t);
}
function CallBlackScholesTheta(S,K,r,sigma,q,T,t)
{
 return (S*NormalDens(D1BlackScholes(S,K,r,sigma,q,T,t))*sigma/(2*Math.sqrt(T-t)))-(r*K*Math.exp(-r*(T-t))*NCDF(D2BlackScholes(S,K,r,sigma,q,T,t),0,1));
}
function PutBlackScholesTheta(S,K,r,sigma,q,T,t)
{
 return (S*NormalDens(D1BlackScholes(S,K,r,sigma,q,T,t))*sigma/(2*Math.sqrt(T-t)))+(r*K*Math.exp(-r*(T-t))*NCDF(-D2BlackScholes(S,K,r,sigma,q,T,t),0,1));
}
function CallBlackScholesRho(S,K,r,sigma,q,T,t)
{
 return K*(T-t)*Math.exp(-r*(T-t))*NCDF(D2BlackScholes(S,K,r,sigma,q,T,t),0,1);
}
function PutBlackScholesRho(S,K,r,sigma,q,T,t)
{
 return -K*(T-t)*Math.exp(-r*(T-t))*NCDF(-D2BlackScholes(S,K,r,sigma,q,T,t),0,1);
}

Start on Normal cumulative distribution function

This a javascript code of this function

For resolve the non centered and reduced distribution we use a Gauss error function.
//Erf function based

function erfAppro6(x)
{
 returne=1.0;
  a1=0.0705230784; a2=0.0422820123; a3=0.0092705272; a4=0.0001520143; a5=0.0002765672; a6=0.0000430638;
 returne=returne-1/Math.pow((1+a1*x+0.230389*x*x+a2*x*x*x+a3*x*x*x*x+a4*x*x*x*x*x+a5*x*x*x*x*x*x+a6*x*x*x*x*x*x*x),16); 
 return returne;
}

//Normal centered and reduct mu==0 && sigma==1
function NCDF2(z)

{
 
 res=1/(1+Math.exp(-(0.07056*z*z*z+1.5976*z)));
 if (res>1.0) res=1.0;
 return res;
}

//Normal with mu & sigma
function NCDF(x,mu,sigma)
{
 if (mu==0&&sigma==1)
 return NCDF2(x);
 res=0.5*(1+erfAppro6(((x-mu)/(sigma*Math.sqrt(2)))));
 if (res>1.0) res=1.0;
 return res;
}