(* :Title: SVJD *)

(* :Context: SVJD` *)

(* :Author: Niels Elken Sønderby *)

(* :Mathematica Version: 4.0 *)

(* :Copyright: Copyright (C) 2003 Niels Elken Sønderby *)

(* Mathematica only implementation of the SVJD model *)

BeginPackage["SVJD`"]

(* Usage messages *)

Keyfigures::usage = 
"Keyfigures[instrument, model] gives all keyfigures (including Value)
for the given instrument under the given model."

Value::usage = 
"Value[instrument, model] is a short form 
for Value /. Keyfigures[instrument, model] 
and gives the value of the instrument under the given model."

CallOption::usage =
"CallOption is an option type where the buyer has the right to buy the 
underlying at maturity."

PutOption::usage =
"PutOption is an option type where the buyer has the right to sell the 
underlying at maturity."

EuropeanOption::usage = 
"EuropeanOption[strike, residualtime, type] represents a european
option of the given type (CallOption, PutOption or Straddle), exercise 
price (strike) and time to maturity (residualtime)."

BlackScholesModel::usage =
"BlackScholesModel[stockprice, riskfreerate, dividend, vol] represents a
model where the underlying asset is worth stockprice at time zero and the
risk free rate, dividend yield and volatility are as given."

SVJDModel::usage =
"SVJDModel[underlying, riskFreeRate, dividendYield, volatility,
volatilityOfVolatility, steadyStateVolatility,
meanReversionRate, correlationUnderlyingVolatility, 
jumpIntensity, jumpMean, jumpStandardDeviation] represents ..."

PDF::usage = "PDF[instr, model, x]"

integr::usage = "test"

H::usage = "test"

DH::usage = "test"

Variance::usage = "test"

XX::usage = "test"

YY::usage = "test"

CC::usage = "test"

Begin["`Private`"]

XX[j_,u_] := - ((alpha T) / sigmav^2) (rho sigmav u - beta[j] - gamma[j, u])

YY[j_,u_] := - ((2 alpha) / sigmav^2) * Log[1 + 0.5 (rho sigmav u - beta[j] - gamma[j, u]) * ((1 - Exp[gamma[j, u] T]) / gamma[j, u])]


TranslateSVJDModel[model_] := 
model /. 
  SVJDModel[underlying_, riskFreeRate_, dividendYield_, volatility_, 
      volatilityOfVolatility_, steadyStateVolatility_, meanReversionRate_, 
      correlationUnderlyingVolatility_,  jumpIntensity_, jumpMean_, 
      jumpStandardDeviation_] :>
    {S = underlying, 
     r = riskFreeRate, rf = dividendYield, 
     V = volatility^2, 
     sigmav = volatilityOfVolatility, 
     alpha = (steadyStateVolatility^2)*meanReversionRate, 
     betastar = meanReversionRate, 
     rho = correlationUnderlyingVolatility,  
     lambdastar = jumpIntensity, 
     kmeanstar = jumpMean, 
     delta = jumpStandardDeviation}

TranslateInstrument[instr_] := 
    instr /. 
    EuropeanOption[strike_, residualTime_, typ_] :> 
    {K = strike, T = residualTime, type = typ}

Keyfigures[EuropeanOption[instrpar__],
	   SVJDModel[modelpar__],
	   opts___?OptionQ] ^:= 
	   Module[{},
	   TranslateSVJDModel[SVJDModel[modelpar]];
	   TranslateInstrument[EuropeanOption[instrpar]];
	   Switch[type,
	   CallOption,
	       {Value -> 
	        S Exp[-rf T] P[1, Log[K/S], V, T, K] - Exp[-r T] K P[2, Log[K/S], V, T, K]},
	   PutOption,
	       {Value -> 
	        S Exp[-rf T] (P[1, Log[K/S], V, T, K] - 1) 
	        - Exp[-r T] K (P[2, Log[K/S], V, T, K] - 1)}
           ]
           ]                     
           
P[j_, x_, V_, T_, K_] := 0.5 + (1/Pi) NIntegrate[integr[j, u, x], {u, 0, 300}]

integr[j_, u_, x_] := Im[F[j, I u] Exp[-I u x]] / u

F[j_, u_] := Exp[CC[j, u] + DD[j, u] V + lambdastar T (1 + kmeanstar)^(my[j] + 0.5) *
                  ((((1 + kmeanstar)^u) Exp[delta^2 (my[j] u + u^2 / 2)]) - 1)]

CC[j_, u_] := ((r - rf - lambdastar kmeanstar) u T 
                 - ((alpha T) / sigmav^2) (rho sigmav u - beta[j] - gamma[j, u])
                 - ((2 alpha) / sigmav^2) * 
                   Log[1 + 0.5 (rho sigmav u - beta[j] - gamma[j, u]) *
                           ((1 - Exp[gamma[j, u] T]) / gamma[j, u])])

DD[j_, u_] :=  -2 (my[j] u + 0.5 u^2) / 
		  (rho sigmav u - beta[j] + gamma[j, u] *
		      (1 + Exp[gamma[j, u] T]) / (1 - Exp[gamma[j, u] T]))

gamma[j_, u_] := Sqrt[(rho sigmav u - beta[j])^2 - 2 sigmav^2 (my[j] u + 0.5 u^2)]
                
my[1] = 0.5; my[2] = -0.5;

beta[1] := betastar - rho sigmav; beta[2] := betastar;

(* Nielsen (1999) p. 49 *)
H[u_] := CC[2, u] + DD[2, u] V + lambdastar T (1 + kmeanstar)^(my[2] + 0.5) *
                  ((((1 + kmeanstar)^u) Exp[delta^2 (my[2] u + u^2 / 2)]) - 1)
                  
DH[u_] := D[H[t],t] /. t->u

Variance[EuropeanOption[instrpar__],
    SVJDModel[modelpar__], 
    opts___?OptionQ] ^:= 
	   Module[{},
	   TranslateSVJDModel[SVJDModel[modelpar]];
	   TranslateInstrument[EuropeanOption[instrpar]];
	   D[H[u], u, u] /. u -> 0	
	   ]

BlackScholesVolatility[EuropeanOption[instrpar__],
    SVJDModel[modelpar__], 
    opts___?OptionQ] ^:= 
	   Module[{},
	   TranslateSVJDModel[SVJDModel[modelpar]];
	   TranslateInstrument[EuropeanOption[instrpar]];
	   Sqrt[(D[H[u], u, u] /. u -> 0)/T]
	   ]

BlackScholesModel[EuropeanOption[instrpar__],
    model_SVJDModel, 
    opts___?OptionQ] ^:= BlackScholesModel[model[[1]], model[[2]],     
    		             model[[3]], 
    		             BlackScholesVolatility[EuropeanOption[instrpar],
    		                                    model]]

PDF[EuropeanOption[strike_, maturity_, type_],
    SVJDModel[underlying_, riskFreeRate_, 
              dividendYield_, volatility_, 
              volatilityOfVolatility_, steadyStateVolatility_,
              meanReversionRate_, correlationUnderlyingVolatility_, 
              jumpIntensity_, jumpMean_, jumpStandardDeviation_], 
    x_,
    opts___?OptionQ] ^:= 
	   Module[{},
	   TranslateSVJDModel[SVJDModel[underlying, riskFreeRate, 
              dividendYield, volatility, 
              volatilityOfVolatility, steadyStateVolatility,
              meanReversionRate, correlationUnderlyingVolatility, 
              jumpIntensity, jumpMean, jumpStandardDeviation]];
	   TranslateInstrument[EuropeanOption[strike, maturity, type]];
           1/Pi NIntegrate[Re[F[2, I u] Exp[-I u x]], {u, 0, 300}]
           ]                     

CDFN[x_] := (1 + Erf[x/Sqrt[2]])/2

PDFN[x_] := 1 / Sqrt[2 Pi] * Exp[-x^2 / 2]

PDFN[m_,s_,x_] := 1/(E^((-m + x)^2/(2*s^2))*Sqrt[2*Pi]*s)

PDF[EuropeanOption[strike_, maturity_, type_], 
    BlackScholesModel[underlying_, riskFreeRate_, dividendYield_, 
      volatility_], x_, opts___?OptionQ] ^:= 
  PDFN[(riskFreeRate - 
              dividendYield - (volatility^2) / 2)maturity, 
      volatility Sqrt[maturity], x]

Value[x___] := Value /. Keyfigures[x]

End[]

EndPackage[]