The Hagan Formula Technology

Rederivations and Extension

Olivier Croissant

Natixis Fixed Income

The equation to Solve

Feynman Kac Theorem : equivalence between

Computing the expectation

E[(Subscript[S, T])-K]  (Subscript[dS, t])=(Subscript[α, t])C[(Subscript[S, t])] (Subscript[dW, t ]); (Subscript[dα, t])=ν (Subscript[dW, t]) ; (Subscript[dα, t])(Subscript[dW, t])=ρ dt

Computing the solution of the PDE

∂P[t,S,α]+(1/2)(α^2)((C[S])^2)(((∂^2)P[t,S,α])/(∂(S^2)))+ρ ν (α^2)C[S](((∂^2)P[t,S,α])/(∂S∂α))+(1/2)(ν^2)(α^2)(((∂^2)P[t,S,α])/(∂(α^2)))=0

P[T,S,α]=(α^2)δ[S-K]

Financial  Interpretation

The price of a call on a stochastic volatility market

The price of the right to buy at  K  at the date T  something whose current price is (Subscript[S, 0])

Solution Strategy

No closed form, only an approximation based on differential geometry and a singular perturbation

Minackshisundaram-Pleijel Expansion

P[T,(Subscript[S, 0]),(Subscript[α, 0]),(Subscript[S, T]),(Subscript[α, T])]=(1/((4πT)))(e^(-(((d((Subscript[S, 0]),(Subscript[α, 0]),(Subscript[S, T]),(Subscript[α, T])))^2)/(2T))))P Δ(1+T (((R+Q))/6)) ;  P : work of the connection; Δ Van Vleck determiant;d geodesic distanceR scalar curvature;  
Q is the mass term;

We have to integrate with respect to α  because we do not observe the (Subscript[α, T])

The smart change of variable

To integrate, we use a Laplace Formula  : (∫_0^∞)(e^( λ f[x]))dx  when λ->∞ , ∃(Subscript[x, s])  /(f^(,))[(Subscript[x, s])]=0 and  (f^(,,))[(Subscript[x, s])]<0

(∫_0^∞)(e^( λ f[x]))dx≅ (Sqrt[((2π)/(λ|(f^(,,))[(Subscript[x, s])]|))])(e^(λf[(Subscript[x, s])]))

therefore

(∫_0^∞)(e^(-(((d[(Subscript[S, 0]),(Subscript[α, 0]),(Subscript[S, T]),(Subscript[α, T])])^2)/(2T)))) d(Subscript[α, T])≅(e^(-(((d[(Subscript[S, 0]),(Subscript[α, 0]),(Subscript[S, T]),(Subscript[α, H])])^2)/(2T))))   where (Subscript[α, H]) is the projection on the hyperplan  (Subscript[S, T])=(OverBar[(Subscript[S, T])])

So a smart  change of variable  is : {(Subscript[S, T]),(Subscript[α, T])}->{d[(Subscript[S, 0]),(Subscript[α, 0]),(Subscript[S, T]),(Subscript[α, H])],γ=(Subscript[α, T])}

What happens on the γ direction will be neglected

The Extraction of factor technique

Let Assume that we have to solve a generic  ((∂P)/(∂t))+(1/2)(α^2)((C[S])^2)(((∂^2)P)/(∂(S^2)))+ρ ν (α^2)C[S](((∂^2)P)/(∂S∂α))+(1/2)(ν^2)(α^2)(((∂^2)P)/(∂(α^2)))+A((∂P)/(∂S))+B((∂P)/(∂α))+M P=0

we assume P[t,S,α] =Q[t,S,α] (Subscript[P, 0])[t,S,α]

then the genreric PDE transforms into :

((∂P)/(∂t))+((1/2)(α^2)((C[S])^2)Q)(((∂^2) (Subscript[P, 0]))/(∂(S^2)))+ρ ν (α^2)C[S] Q[t,S,α](((∂^2)(Subscript[P, 0]))/(∂S∂α))+(1/2)(ν^2)(α^2)Q(((∂^2) (Subscript[P, 0]))/(∂(α^2)))+((α^2)((C[S])^2)((∂ Q)/(∂S))+ρ ν (α^2)C[S]((∂Q)/(∂α))+A Q[t,S,α])((∂(Subscript[P, 0]))/(∂S))+
((ν^2)(α^2)((∂Q)/(∂α))+ρ ν (α^2)C[S]((∂Q)/(∂S))+B Q)((∂(Subscript[P, 0]))/(∂α))+((1/2)(α^2)((C[S])^2)(((∂^2)Q)/(∂(S^2)))+ρ ν (α^2)C[S](((∂^2)Q)/(∂S∂α))+(1/2)(ν^2)(α^2)(((∂^2)Q)/(∂(α^2)))+M Q)(Subscript[P, 0])=0

If we want to remove A ≠0

we just have to find a Q[S]  such that  (α^2)((C[S])^2)((∂ Q)/(∂S))+A Q[t,S,α] =0

If we want to remove B ≠0

we just have to find a Q[α]  such that  (ν^2)(α^2)((∂Q)/(∂α))+B Q =0

Consequence for the mass term:

every factor removal implies a change of mass equal  ΔM=(((1/2)(α^2)((C[S])^2)(((∂^2)Q)/(∂(S^2)))+ρ ν (α^2)C[S](((∂^2)Q)/(∂S∂α))+(1/2)(ν^2)(α^2)(((∂^2)Q)/(∂(α^2))))/Q)

Primitives

Clear[a1,v,u,x,X,H,M,f,K,S,F,T,α,β,ρ,ν,λ,μ,φ,ω,g,B,z,x,b1,b2,Z,I0,I1,I2,favg,zavg,f1,H]

Jetmatrix[u_,v_,x_,y_]:={
{D[u,x],D[v,x],0,0,0},
{D[u,y],D[v,y],0,0,0},
{D[u,x,x],D[v,x,x],D[u,x]*D[u,x],2D[u,x]*D[v,x],D[v,x]*D[v,x]},
{D[u,x,y],D[v,x,y],D[u,x]*D[u,y],D[u,x]*D[v,y]+D[u,y]*D[v,x],D[v,x]*D[v,y]},
{D[u,y,y],D[v,y,y],D[u,y]*D[u,y],2D[u,y]*D[v,y],D[v,y]*D[v,y]}
}

ChangeOfVariable[supersextivector_,v1_,v2_,v1new_,v2new_,v1newDef_,v2newDef_]:=
Module[{jet,U ,rules,eq,eqnew,function,sextivector=supersextivector[[1]],initialcondition=supersextivector[[2]],newsextivector,newinitialcondition},
eq=sextivector[[1]] D[function[v1,v2],v1]+
sextivector[[2]] D[function[v1,v2],v2]+
sextivector[[3]] D[function[v1,v2],v1,v1]+
sextivector[[4]] D[function[v1,v2],v1,v2]+
sextivector[[5]] D[function[v1,v2],v2,v2]+
sextivector[[6]]function[v1,v2] ;
U={D[function[v1new,v2new],v1new],D[function[v1new,v2new],v2new],D[function[v1new,v2new],v1new,v1new],D[function[v1new,v2new],v1new,v2new],D[function[v1new,v2new],v2new,v2new]};
jet=Simplify[Jetmatrix[v1newDef[v1,v2],v2newDef[v1,v2],v1,v2].U] ;
rules= {D[function[v1,v2],v1]->jet[[1]],D[function[v1,v2],v2]->jet[[2]],D[function[v1,v2],v1,v1]->jet[[3]],D[function[v1,v2],v1,v2]->jet[[4]],D[function[v1,v2],v2,v2]->jet[[5]],function[v1,v2]->function[v1new,v2new]};
eqnew=eq /. rules;
newsextivector=Simplify[{
Coefficient[eqnew,U[[1]]],
Coefficient[eqnew,U[[2]]],
Coefficient[eqnew,U[[3]]],
Coefficient[eqnew,U[[4]]],
Coefficient[eqnew,U[[5]]],
Coefficient[eqnew,function[v1new,v2new]]}];
newinitialcondition=initialcondition;
newinitialcondition[[1]]=initialcondition[[1]] Derivative[0,1][v2newDef][v1,initialcondition[[2]] ] /. t->0;
newinitialcondition[[2]]=v2newDef[v1,initialcondition[[2]] ];
{newsextivector,newinitialcondition}
]

ExtractFactor[supersextivector_,factor_,functionname_,v1_,v2_,newfunctionname_]:=
Module[{eq,neweq,f0,newf0,factor0,sextivector=supersextivector[[1]],initialcondition=supersextivector[[2]],newsextivector,newinitialcondition},
factor0=factor /. {functionname->f0};
f0[v1_,v2_]:=factor0 newf0[v1,v2];
eq=sextivector[[1]] D[f0[v1,v2],v1]+
sextivector[[2]] D[f0[v1,v2],v2]+
sextivector[[3]] D[f0[v1,v2],v1,v1]+
sextivector[[4]] D[f0[v1,v2],v1,v2]+
sextivector[[5]] D[f0[v1,v2],v2,v2]+
sextivector[[6]]f0[v1,v2] ;
neweq=Simplify[eq];
newEq1=Collect[neweq,{D[newf0[v1,v2],v1],D[newf0[v1,v2],v2], D[newf0[v1,v2],v1,v1],D[newf0[v1,v2],v1,v2],D[newf0[v1,v2],v2,v2],newf0[v1,v2]}];
newsextivector=Simplify[({
Coefficient[neweq,D[newf0[v1,v2],v1]],
Coefficient[neweq,D[newf0[v1,v2],v2]],
Coefficient[neweq, D[newf0[v1,v2],v1,v1]],
Coefficient[neweq, D[newf0[v1,v2],v1,v2]],
Coefficient[neweq,D[newf0[v1,v2],v2,v2]],
Coefficient[neweq,newf0[v1,v2]]}/factor0) /. {newf0->newfunctionname}];
newinitialcondition=initialcondition;
newinitialcondition[[1]]=newinitialcondition[[1]] /  (factor /. {v2->initialcondition[[2]],t->0});
{newsextivector,newinitialcondition}
]

PruneEpsilonOrder[exp_,order_]:=Normal[Series[Numerator[exp],{ϵ,0,order}]]/Denominator[exp]

ExpPruneOrder[exp_,factor_,norder_]:=Module[{n},Sum[Coefficient[Numerator[exp],factor,n] (factor^n),{n,0,norder}]/Denominator[exp]]

EDPPruneOrder[supersextivector_,rank_,factor_,norder_]:=Module[{sextivector=supersextivector[[1]],initialcondition=supersextivector[[2]],newsextivector},
newsextivector=sextivector;
newsextivector[[rank]]=Simplify[ExpPruneOrder[newsextivector[[rank]],factor,norder]];
{newsextivector,initialcondition}
]

PrintEquaDiff[sextiv_]:=Module[{},
Print["Coef{((∂H)/(∂α))}=",sextiv[[1,1]]];
Print["Coef{((∂H)/(∂f))}=",sextiv[[1,2]]];
Print["Coef{(((∂^2)H)/(∂(α^2)))}=",sextiv[[1,3]]];
Print["Coef{(((∂^2)H)/(∂α ∂f))}=",sextiv[[1,4]]];
Print["Coef{(((∂^2)H)/(∂(f^2)))}=",sextiv[[1,5]]];
Print["Coef{H}=",sextiv[[1,6]]];
Print["Coef{delta}=",sextiv[[2,1]]];
Print["abcisse{delta}=",sextiv[[2,2]]];
]

PrintEquaDiffDeltaPart[sextiv_]:=Module[{},
Print["Coef{delta}=",sextiv[[2,1]]];
Print["abcisse{delta}=",sextiv[[2,2]]];
]

Mathematica representation of the Hagan PDE  with its boundary condition

vec0={{0 ,0,(1/2)   (ϵ^2) (ν^2) (α^2),(α^2) (ϵ^2) ν ρ C[f],(1/2) (α^2) (ϵ^2) ((C[f])^2),0},{(α^2),K}}

PrintEquaDiff[vec0]

First Change of Variable : 
  purpose : to remove the skew generating local vol part

v[a_,f_]:=1/(ϵ α) Integrate[1/C[s],{s,K,f}]

u[a_,f_]:=a

vec2=FullSimplify[ChangeOfVariable[vec0,α,f,α,z,u,v] /.  {(∫_K^f)(1/(C[s]))\[DifferentialD]s->ϵ α z,C[f]->B[ϵ α z],C'[f]->((B'[ϵ α z])/(B[ϵ α z]))}];

PrintEquaDiff[vec2]

Extraction of  (α/(ϵ C[K])) to make the inital condition be δ (z)

we see that the leading terms make the solution gaussian at first order  (equation 7.3 .15 of the paper )

vec21=FullSimplify[ExtractFactor[vec2,α/(ϵ C[K]),X,α,z,H]] ;vec21[[1]] // ColumnForm

PrintEquaDiffDeltaPart[vec21]

PrintEquaDiff[vec2]

Extraction of (Sqrt[B[α z ϵ]/B[0]]) to prepare for zeroing the first derivative with respect to z

vec22=FullSimplify[ExtractFactor[vec21, (Sqrt[B[α z ϵ]/B[0]]),X,α,z,H]/. {B'[ϵ z α]->b1 B[ϵ z α], B''[ϵ z α]->b2  B[ϵ z α]} ];

PrintEquaDiff[vec21]

We see the coefficient of z is ((α ϵ B'[z α ϵ])/(2 B[z α ϵ])) which negligeable at the (ϵ^(2 ))order

Removal of negligeable terms

we neglect terms of order (ϵ^2) and higher for the ∂α term

vec23=FullSimplify[EDPPruneOrder[vec22, 1,ϵ,0]];PrintEquaDiff[vec23]

we neglect terms of order (ϵ^2) and higher for the ∂αz term

vec24=FullSimplify[EDPPruneOrder[vec23, 4,ϵ,1]];vec24[[1]] // ColumnForm

we neglect terms of order (ϵ^2) and higher for the ∂αα term  (* eq 7.3.21 *)

vec25=FullSimplify[EDPPruneOrder[vec24, 3,ϵ,1]];vec25[[1]] // ColumnForm

we neglect terms of order ϵ and higher for the ∂αz term  (* eq 7.3 .24 *)

vec26=FullSimplify[EDPPruneOrder[vec25, 4,ϵ,0]];vec26[[1]] // ColumnForm

PrintEquaDiff[vec26]

To remove the derivatives z, first let determine the function to extract

let introduce the function to extract : f1

Clear[f1];vec27a=FullSimplify[ExtractFactor[vec26,f1[z],X,α,z,H]];vec27a[[1]] //ColumnForm

let solve the differential equation

SolForf=DSolve[vec27a[[1,2]]==0,f1,z]

Let affect f1 to the solution

Clear[f1];f1=f1 /. Simplify[SolForf[[1]] /. C[1]->1 ]

We can check that f[0] = 1  , such that the initial boundary will be perturbated

f1[0]

Let' s perform the extraction

vec27=ExtractFactor[vec26,f1[z]/f1[0],X,α,z,H];vec27[[1]] //ColumnForm

PrintEquaDiffDeltaPart[vec27]

We look for  a change of variable to completly standardize the second \
derivative with respect to z

Clear[xx];xx[α_,z_]:=f[z]

αα[α_,z_]:= α

vec3a=FullSimplify[ChangeOfVariable[vec27,α,z,α,x,αα,xx]] ;vec3a[[1]] // ColumnForm

chse1=DSolve[vec3a[[1,5]]==1/2,f,z]

xx=Function[{α,z},Evaluate[(f /. chse1[[2,1]])[z]-C[1]]]

Clear[xx1];xx1=Function[{α,z},Evaluate[Simplify[xx[α,z]-xx[α,0]]]]

We perform the change of variable  : we should get 1/2 on the 5 th line

vec3=FullSimplify[ChangeOfVariable[vec27,α,z,α,x,αα,xx1],α>0] ;vec3[[1]] // ColumnForm

vec3=Simplify[vec3, α>0];

PrintEquaDiffDeltaPart[vec3]

Simplification  used to prove the identities

We want to prove that vec3[[1,2]]=-((Z''[x])/(2 Z'[x])) for some function Z[x]

zz[x_,ν_]:=(Sinh[ν x]-ρ (Cosh[ν x] -1))/ν

Z1[z_,ν_]:=1/ν (Log[(Sqrt[1-2ρ ν z+ν^2 z^2 ])-ρ+ν z]-Log[1-ρ])

Computation of the Derivative ratio

vect3proxi=FullSimplify[FullSimplify[-Derivative[2,0][zz][x,ϵ ν]/(2Derivative[1,0][zz][x,ϵ ν])  /. {x->Z1[z,ϵ ν]}]//. {Log[q_]-Log[q1_]->Log[q/q1],Log[q_]+Log[q1_]->Log[q q1],(Sqrt[1+((q_)/((q1_)^2))])->((Sqrt[(q1^2)+q])/q1)}]

Simplify[vec3[[1,2]]-vect3proxi]

So it matches

We assume the element is of the shape precedingly checked

vec3[[1,2]]=-((Z''[x])/(2 Z'[x]))

vec3a=vec3 /. {z->Z[x]} ;vec3a[[1]] // ColumnForm

Clear[g];vec3b=ExtractFactor[vec3a,g[x],X,α,x,H];vec3b[[1]] // ColumnForm

Clear[g];DSolve[Simplify[(vec3b[[1,2]]/.{ }),Z'[x]>0]==0,g,x]

We perform the extraction

g[y_]:=(Sqrt[Z'[y]])

vec3c=ExtractFactor[vec3a,g[x],X,α,x,H];vec3c[[1]] // ColumnForm

We simplify

Clear[I0,I1,I2,b1,b2,z,x];vec3d=FullSimplify[vec3c /. {Z[x]->z,Z'[x]->I0,Z''[x]->d I1 I0,Derivative[3][Z][x]->(d^2)(I2 (I0^2)+(I1^(2 ))I0)}]; vec3d[[1]] //ColumnForm

We can keep only the leading terms here we : 
  only the (ϵ^2) order

Clear[I0,I1,I2,b1,b2,z,x];vec3e=PruneEpsilonOrder[vec3d,2]; vec3e[[1]] //ColumnForm

Solution of the partial differential equation

Clear[x,Q]

Q[x_,t_]:=(1/(Sqrt[2π t]))(E^(-((x^2)/(2t))+(∫_0^t)κ[s]\[DifferentialD]s))

FullSimplify[Derivative[0,1][Q][x,t]-(1/2)Derivative[2,0][Q][x,t]-κ[t]Q[x,t]]

Example of Implementation

SABRHaganVol[f_,K_,α_,β_,ρ_,ν_,T_]:=Module[{z,x},z=(((f^(1-β))-(K^(1-β)))/(α (1-β)));
x=((Log[((-ρ+ν z+(Sqrt[1-2 ν ρ z+(ν^2) (z^2)]))/(1-ρ))])/ν);
(α/(((f K))^((1-β)/(2 )))) (z/x)  (((1+((1/24) (α^2) (((1-β))^2) (((f K))^(-1+β))+(1/4) α β (((f K))^((1/2) (-1+β))) ν ρ+(1/24) (ν^2) (2-3 (ρ^2))) T))/((1+(1/24) (((1-β))^2) ((Log[(f/K)])^2)+(((((1-β))^4) ((Log[(f/K)])^4))/1920))))  ]

Module[{S=0.05,β=0.7,ρ1=-0.4,ρ2=-0.,ρ3=0.4,ν=0.2,α=0.11,T=10},
Plot[{ SABRHaganVol[S,K,α,β,ρ1,ν,T], SABRHaganVol[S,K,α,β,ρ2,ν,T], SABRHaganVol[S,K,α,β,ρ3,ν,T]},{K,0.01,0.2}]]

Extension

Drift of the volatility

Important feature that can be used to better match time structure of the observed volatilities

Mean reversion of the volatility

Important  feature  to model volatility  for long maturities (30 years).

BiSABR

Modeling of option on the spread between two indexes modeled as SABR: Huge positions todays.

Reference

MANAGING SMILE RISK (Working Paper) :
 PATRICK S.HAGAN∗, DEEP KUMAR†, ANDREW S.LESNIEWSKI‡, AND DIANA E.WOODWARD (1999)