BeginPackage["Options`HestonVanilla`"]

(************************************************************************************
author: Uwe Wystup, wystup@mathfinance.de
date  : November 1999
*************************************************************************************)
AuxFunc::usage = "computing the ingredients for HestonVanilla"
HestonVanilla::usage = "Steven Heston's Stochastic Volatility Model\n
to price European put and call options\n
HestonVanilla[k,o,sigma,rho,lambda,r,rf,v,S,K,tau,cp,G]\n
The input parameters are:\n
k:      MeanReversion\n
o:      LongRunVariance\n
sigma:  VolaVolatility\n
rho:    Correlation\n
lambda: VolaRiskPremium\n
r:      domestic RiskFreeRate\n
rf:     foreign RiskFreeRate\n
v:      CurrentVariance\n
S:      AssetPrice\n
K:      ExercisePrice\n
tau:    ExpirationTime\n
cp :    1 for call, -1 for put\n
G:      Greek\n
		0 : value\n
		1 : spot delta\n
		2 : spot gamma\n
		3 : theta (in years)\n
		4 : vega (wrt v)\n
		5 : domestic rho\n
		6 : foreign rho\n
		7 : vomma (wrt v)\n
		21: dual delta (wrt K)\n
		22: dual gamma (wrt K)"


Begin["`Private`"]

AuxFunc[k_,o_,sigma_,rho_,lambda_,
	r_,rf_,v_,S_,K_,tau_,fi_] := 
  Block[{u,a,b,x,rsf,d,g,cl,dl,f},
  u = {0.5,-0.5};
  a = k*o;
  b = {k+lambda-rho*sigma, k+lambda};
  x = Log[S];
  rsf = rho*sigma*fi;

  d = Table [
	Sqrt[(I*rsf - b[[j]])^2 - 
	sigma^2*(2*I*u[[j]]*fi - fi^2)],
	{j,2}];

  g = Table [
 	(b[[j]] - I*rsf + d[[j]]) /
	(b[[j]] - I*rsf - d[[j]]),
	{j,2}];

  cl = Table [
	I*(r-rf)*fi*tau + 
	(a/sigma^2)*((b[[j]]-I*rsf+d[[j]])*tau - 
	2*Log[(1-g[[j]]*Exp[d[[j]]*tau])/(1-g[[j]])]),
	{j,2}];

  dl = Table [
	((b[[j]] - I*rsf + d[[j]])/sigma^2) * 
	((1-Exp[d[[j]]*tau])/
	 (1-g[[j]]*Exp[d[[j]]*tau])),
	{j,2}];

  f = Table [
	Exp[cl[[j]] + dl[[j]]*v + I*fi*x], {j,2}];

  Table [{
    Re[(Exp[-I*fi*Log[K]] * f[[j]])/(I*fi)],
    Re[(Exp[-I*fi*Log[K]] * f[[j]])],
    Re[(dl[[j]]*Exp[-I*fi*Log[K]] * f[[j]])/(I*fi)],
    Re[(dl[[j]]*dl[[j]]*Exp[-I*fi*Log[K]] * f[[j]])/(I*fi)],
    d[[j]],
    f[[j]]
    },{j,2}]];


HestonVanilla[k_,o_,sigma_,rho_,lambda_,
	r_,rf_,v_,S_,K_,tau_,cp_,G_] :=   
Switch[G,
0, (*value*)
    	Block[{p,j},
		p = Table[
		0.5+(1/N[Pi])*NIntegrate[
		AuxFunc[k,o,sigma,rho,lambda,Log[1+r],Log[1+rf],v,S,K,tau,
		fi][[j,1]], {fi,0.,100}],
		{j,2}];
	S*Exp[-Log[1+rf]*tau]*(p[[1]]-(1-cp)/2) 
	- K*Exp[-Log[1+r]*tau]*(p[[2]]-(1-cp)/2)],
1, (*delta*)
    	Block[{p,j},
		p = Table[
		0.5+(1/N[Pi])*NIntegrate[
		AuxFunc[k,o,sigma,rho,lambda,Log[1+r],Log[1+rf],v,S,K,tau,
		fi][[j,1]], {fi,0.,100}],
		{j,1}];
	Exp[-Log[1+rf]*tau]*(p[[1]]-(1-cp)/2)],
2, (*gamma*)
	Block[{dp,j},
    		dp=Table[
		(1/N[Pi])*NIntegrate[
		AuxFunc[k,o,sigma,rho,lambda,Log[1+r],Log[1+rf],v,S,K,tau,
		fi][[j,2]], {fi,0.,100}],
		{j,1}];
	Exp[-Log[1+rf]*tau]/S*dp[[1]]],
3, (*theta*)
	50000*
	(HestonVanilla[k,o,sigma,rho,lambda,r,rf,v,S,K,tau-0.00001,cp,0]
	-HestonVanilla[k,o,sigma,rho,lambda,r,rf,v,S,K,tau+0.00001,cp,0]),
4, (*vega*)
	Block[{dpv,j},
	 	dpv=Table[
		(1/N[Pi])*NIntegrate[
		AuxFunc[k,o,sigma,rho,lambda,Log[1+r],Log[1+rf],v,S,K,tau,
		fi][[j,3]], {fi,0.,100}],
		{j,2}];
	S*Exp[-Log[1+rf]*tau]*dpv[[1]] 
	- K*Exp[-Log[1+r]*tau]*dpv[[2]]],
5, (*rho*)
	Block[{p,j},
   		p = Table[
		0.5+(1/N[Pi])*NIntegrate[
		AuxFunc[k,o,sigma,rho,lambda,Log[1+r],Log[1+rf],v,S,K,tau,
		fi][[j,1]], {fi,0.,100}],
		{j,2,2}];
    	K*tau*Exp[-Log[1+r]*tau]*(p[[1]]-(1-cp)/2)/(1+r)],
6, (*rhof*)
	Block[{p,j},
   		p = Table[
		0.5+(1/N[Pi])*NIntegrate[
		AuxFunc[k,o,sigma,rho,lambda,Log[1+r],Log[1+rf],v,S,K,tau,
		fi][[j,1]], {fi,0.,100}],
		{j,1}];
    	S*tau*Exp[-Log[1+rf]*tau]*((1-cp)/2-p[[1]])/(1+rf)],
7, (*vomma*)
	Block[{d2pv,j},
		d2pv=Table[
		(1/N[Pi])*NIntegrate[
		AuxFunc[k,o,sigma,rho,lambda,Log[1+r],Log[1+rf],v,S,K,tau,
		fi][[j,4]], {fi,0.,100}],
		{j,2}];
	S*Exp[-Log[1+rf]*tau]*d2pv[[1]] 
	- K*Exp[-Log[1+r]*tau]*d2pv[[2]]],
21, (*dual delta*)
	Block[{p,j},
   		p = Table[
		0.5+(1/N[Pi])*NIntegrate[
		AuxFunc[k,o,sigma,rho,lambda,Log[1+r],Log[1+rf],v,S,K,tau,
		fi][[j,1]], {fi,0.,100}],
		{j,2,2}];
	Exp[-Log[1+r]*tau]*((1-cp)/2-p[[1]])],
22, (*dual gamma*)
	Block[{dp,j},
    		dp=Table[
		(1/N[Pi])*NIntegrate[
		AuxFunc[k,o,sigma,rho,lambda,Log[1+r],Log[1+rf],v,S,K,tau,
		fi][[j,2]], {fi,0.,100}],
		{j,2,2}];
	Exp[-Log[1+r]*tau]/K*dp[[1]]],
101, (* Real part *)
	Block[{dp,j},
    		dp=Table[
		AuxFunc[k,o,sigma,rho,lambda,Log[1+r],Log[1+rf],v,S,K,tau,
		2][[j,1]],
		{j,2}];
	dp[[2]]],
111, (* d *)
	Block[{dp,j},
    		dp=Table[
		AuxFunc[k,o,sigma,rho,lambda,Log[1+r],Log[1+rf],v,S,K,tau,
		2][[j,5]],
		{j,2}];
	dp[[2]]],
112, (*probability factors*)
    	Block[{p,j},
		p = Table[
		0.5+(1/N[Pi])*NIntegrate[
		AuxFunc[k,o,sigma,rho,lambda,Log[1+r],Log[1+rf],v,S,K,tau,
		fi][[j,1]], {fi,0.,100}],
		{j,2}];
	p[[1]]],
30, (*the Heston Integrand*)
	AuxFunc[k,o,sigma,rho,lambda,Log[1+r],Log[1+rf],v,S,K,tau,cp][[1,1]],
_,0];

End[]

EndPackage[]