BlackModel.clj

(use 'com.lambder.deriva.core)

(def N 
  '(/ 1
      (+ 1
         (exp (- 
                (* -0.07056 (pow x 3))
                (* -1.5976 x)))))

(def d1 '(/ (+ (/ F K) (* T (/ (sq sigma) 2)))))

(def d2 '(/ (- (/ F K) (* T (/ (sq sigma) 2)))))

(def call
    `(* 
        (exp (neg (* r T))
        (- 
          (* F ~(bind N x d1))
          (* K ~(bind N x d2)))))

(def put
    `(* 
        (exp (neg (* r T))
        (- 
          (* F ~(bind N x (neg d1)))
          (* K ~(bind N x (neg d2))))))

(defn black-expression [call?] call put)


;; usage

(def black-model-with-sensitivities (∂ (bind (black-expression true) T 0.523) F K r))

(black-model-with-sensitivities 12.3 14.3 0.03)
(black-model-with-sensitivities 12.3 11.0 0.03)
(black-model-with-sensitivities 12.3 11.0 0.02)

BlackModel.java

import static com.lambder.deriva.Deriva.*;

public class Formulas {

  public static Expression black(final boolean isCall) {

    // Logistic approximation of Cumulated Standard Normal Distribution
    // 1/( e^(-0.07056 * x^3 - 1.5976*x) + 1)
    Expression N =  div(1.0,
        add(
            exp(
                sub(
                    mul(
                        -0.07056,
                        pow('x', 3)),
                    mul(-1.5976, 'x'))),
            1.0));

    // ( F/K+T*σ^2/2 ) / σ*sqrt(T)
    Expression d1 = div(
        add(
            div('F', 'K'),
            mul(div(sq("sigma"), 2.0), 'T')),
        mul("sigma", sqrt('T')));

    // ( F/K-T*σ^2/2 ) / σ*sqrt(T)
    Expression d2 = div(
        sub(
            div('F', 'K'),
            mul(div(sq("sigma"), 2.0), 'T')),
        mul("sigma", sqrt('T')));

    // e^(-r*T) * ( F*N(d1)-K*N(d2) )
    Expression call = mul(
        exp(neg(mul('r', 'T'))),
        sub(
            mul('F', N.bind('x', d1)),
            mul('K', N.bind('x', d2))));

    // e^(-r*T) * ( F*N(-d2)-K*N(-d1) )
    Expression put = mul(
        exp(neg(mul('r', 'T'))),
        sub(
            mul('F', N.bind('x', neg(d2))),
            mul('K', N.bind('x', neg(d1)))));

    return isCall ? call : put;    
  }

  // usage
  public static void main(String[] args) {
    // lets fix timeToExpiry to 0.523 and get only strike, forward and lognormalVol sensitivities
    Expression blackModel = black(true).bind('T', 0.523);
    Function1D fun = d(blackModel, 'F', 'K', 'r').function('F', 'K', 'r');
    fun.execute(12.3, 14.3, 0.03);
    fun.execute(12.3, 11.0, 0.03);
    fun.execute(12.3, 11.0, 0.02);
  }
}

PartialDifferentiation.java

Expression expr = sin(mul(sq('x'), sq('y')));
Vector1D g_expr = vector(expr, d(expr, 'x'), d(expr, 'y'));
System.out.println(g_expr.describe());

Expression:
[(sin (mul (sq x) (sq y)))
 (d (sin (mul (sq x) (sq y))) x)
 (d (sin (mul (sq x) (sq y))) y)]

gets turned into:
(vector
 (sin (* (sq x) (sq y)))
 (* (cos (* (sq x) (sq y))) (* (* 2 x) (sq y)))
 (* (cos (* (sq x) (sq y))) (* (* 2 y) (sq x))))

 and into:
final double G__33 = x;  //  x
final double G__32 = sq( G__33 );  //  (sq x)
final double G__31 = y;  //  y
final double G__30 = 2;  //  2
final double G__29 = G__30 * G__31;  //  (* 2 y)
final double G__28 = G__29 * G__32;  //  (* (* 2 y) (sq x))
final double G__26 = sq( G__31 );  //  (sq y)
final double G__23 = G__32 * G__26;  //  (* (sq x) (sq y))
final double G__22 = cos( G__23 );  //  (cos (* (sq x) (sq y)))
final double G__21 = G__22 * G__28;  //  (* (cos (* (sq x) (sq y))) (* (* 2 y) (sq x)))
final double G__16 = G__30 * G__33;  //  (* 2 x)
final double G__15 = G__16 * G__26;  //  (* (* 2 x) (sq y))
final double G__8 = G__22 * G__15;  //  (* (cos (* (sq x) (sq y))) (* (* 2 x) (sq y)))
final double G__2 = sin( G__23 );  //  (sin (* (sq x) (sq y)))
final double G__1 = [ G__2, G__8, G__21 ];  //  (vector (sin (* (sq x) (sq y))) (* (cos (* (sq x) (sq y))) (* (* 2 x) (sq y))) (* (cos (* (sq x) (sq y))) (* (* 2 y) (sq x))))