{-# LANGUAGE RecordWildCards #-}
module Q.Options.ImpliedVol.Normal where
import           Data.Default.Class
import           Numeric.IEEE                   (epsilon, maxFinite, minNormal)
import           Numeric.Natural
import           Numeric.RootFinding
import           Q.Options.Bachelier
import           Q.Types
import           Statistics.Distribution        (cumulative, density)
import           Statistics.Distribution.Normal (standard)


-- | Method to use to calculate the normal implied vol.
data Method =
    Jackel        -- ^ Jackel analytical formula approximation.
  | ChoKimKwak    -- ^ J. Choi, K kim, and M. Kwak (2009)
  -- | Numerical root finding. Currently Ridders is used.
  | RootFinding {
        Method -> Natural
maxIter ::  Natural                 -- ^ Maximum number of iterations.
      , Method -> Tolerance
tol     ::  Tolerance               -- ^ Tolerance (relative or absolute)
      , Method -> (Double, Double, Double)
bracket :: (Double, Double, Double) -- ^ Triple of @(low bound, initial
                                            --   guess, upper bound)@. If initial
                                            --   guess if out of bracket middle
                                            --   of bracket is taken as.
        }

instance Default Method where
  def :: Method
def = Method
Jackel

-- | Default method implementation of 'euImpliedVolWith' using 'Jackel'.
euImpliedVol :: OptionType
-> Forward -> Strike -> YearFrac -> Rate -> Premium -> Vol
euImpliedVol = Method
-> OptionType
-> Forward
-> Strike
-> YearFrac
-> Rate
-> Premium
-> Vol
euImpliedVolWith Method
forall a. Default a => a
def

-- | Calcualte the bachelier option implied vol of a european option.
--
-- If the options premium does not have time value @'hasTimeValue'@ return 0.
euImpliedVolWith :: Method -> OptionType -> Forward -> Strike -> YearFrac -> Rate -> Premium -> Vol
euImpliedVolWith :: Method
-> OptionType
-> Forward
-> Strike
-> YearFrac
-> Rate
-> Premium
-> Vol
euImpliedVolWith Method
m OptionType
cp Forward
f Strike
k YearFrac
t Rate
r Premium
p
  | OptionType -> Forward -> Strike -> Premium -> DF -> Bool
hasTimeValue OptionType
cp Forward
f Strike
k Premium
p DF
df = Method
-> OptionType
-> Forward
-> Strike
-> YearFrac
-> Rate
-> Premium
-> Vol
euImpliedVolWith' Method
m OptionType
cp Forward
f Strike
k YearFrac
t Rate
r Premium
p
  | Bool
otherwise                = Double -> Vol
Vol (Double -> Vol) -> Double -> Vol
forall a b. (a -> b) -> a -> b
$ Double
0
  where df :: DF
df = YearFrac -> Rate -> DF
discountFactor YearFrac
t Rate
r

euImpliedVolWith' :: Method
-> OptionType
-> Forward
-> Strike
-> YearFrac
-> Rate
-> Premium
-> Vol
euImpliedVolWith' Method
Jackel OptionType
cp (Forward Double
f) (Strike Double
k) (YearFrac Double
t) (Rate Double
r) (Premium Double
p)
  -- Case where interest rate is not 0 we need undiscount. The paper is written
  -- for the undiscounted Bachelier option prices.
  | Double
r Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0
    = OptionType
-> Forward -> Strike -> YearFrac -> Rate -> Premium -> Vol
euImpliedVol OptionType
cp (Double -> Forward
Forward Double
f) (Double -> Strike
Strike Double
k) (Double -> YearFrac
YearFrac Double
t) (Double -> Rate
Rate Double
0) (Double -> Premium
Premium (Double
pDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
df))
  -- Case of ATM. Solve directly.
  | Double -> Double
forall a. Num a => a -> a
abs (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
f) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
forall a. IEEE a => a
epsilon = Double -> Vol
Vol (Double -> Vol) -> Double -> Vol
forall a b. (a -> b) -> a -> b
$ Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sqrt2Pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double -> Double
forall a. Floating a => a -> a
sqrt Double
t)
  -- Case of ITM option. Calcualte vol of the out of the money option with Put-Call-Parity.
  | Double
phiStarTilde Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0
    = OptionType
-> Forward -> Strike -> YearFrac -> Rate -> Premium -> Vol
euImpliedVol (OptionType -> OptionType
forall a. Enum a => a -> a
succ OptionType
cp) (Double -> Forward
Forward Double
f) (Double -> Strike
Strike Double
k) (Double -> YearFrac
YearFrac Double
t) (Double -> Rate
Rate Double
r) (Double -> Premium
Premium Double
p')
  -- General case for an out of the money option.
  | Bool
otherwise  = let
      ẋ :: Double
ẋ      = if Double
phiStarTilde Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
c then
                 let g :: Double
g = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
phiStarTilde Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5)
                     ξ :: Double
ξ = (Double
0.032114372355 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
gDouble -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
2)Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
0.016969777977 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
gDouble -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
2)Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
2.6207332461E-3Double -> Double -> Double
forall a. Num a => a -> a -> a
-(Double
9.6066952861E-5)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
gDouble -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
2)))
                         Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/
                         (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-(Double
gDouble -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
2)Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
0.6635646938 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
gDouble -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
2)Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
0.14528712196 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.010472855461Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
gDouble -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
2)))
                 in Double
g Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sqrt2Pi Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ξDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
gDouble -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
2)
               else
                 let h :: Double
h = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ (-Double -> Double
forall a. Floating a => a -> a
log (-Double
phiStarTilde))
                 in (Double
9.4883409779Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
hDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
9.6320903635Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
hDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
0.58556997323 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2.1464093351Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
h)))
                    Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/
                    (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
hDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
0.65174820867 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
hDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
1.5120247828 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
6.6437847132E-5Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
h)))
      c :: Double
c       = (-Double
0.001882039271)
      x :: Double
x       = Double
ẋ Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
3Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
ẋ Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
ẋ Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
ẋ Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ẋDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ẋ)))
                    Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/
                    (Double
6 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
qDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ẋ Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((-Double
12) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ẋ Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
6Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ẋ Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((-Double
6)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
qDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ẋDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
3Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
ẋDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ẋ)))))
      phiXBarTilde :: Double
phiXBarTilde = (NormalDistribution -> Double -> Double
forall d. Distribution d => d -> Double -> Double
cumulative NormalDistribution
standard Double
ẋ) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (NormalDistribution -> Double -> Double
forall d. ContDistr d => d -> Double -> Double
density NormalDistribution
standard Double
ẋ)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
ẋ
      q :: Double
q       =  (Double
phiXBarTildeDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
phiStarTilde)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (NormalDistribution -> Double -> Double
forall d. ContDistr d => d -> Double -> Double
density NormalDistribution
standard Double
ẋ)
    in Double -> Vol
Vol (Double -> Vol) -> Double -> Vol
forall a b. (a -> b) -> a -> b
$ (Double -> Double
forall a. Num a => a -> a
abs (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
f)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double -> Double
forall a. Num a => a -> a
abs (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt Double
t))
  where phiStarTilde :: Double
phiStarTilde = Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double
forall a. Num a => a -> a
abs (Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (Double
theta Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
k)) Double
0))) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double -> Double
forall a. Num a => a -> a
abs (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
f))
        theta :: Double
theta        = if OptionType
cp OptionType -> OptionType -> Bool
forall a. Eq a => a -> a -> Bool
== OptionType
Call then Double
1 else -Double
1
        phiTilde :: Double
phiTilde     = (-Double
theta) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
p Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
f)
        p' :: Double
p'           = Double
cpi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
df Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
k) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
p
        cpi :: Double
cpi          = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Double) -> Int -> Double
forall a b. (a -> b) -> a -> b
$ OptionType -> Int
forall a. Enum a => a -> Int
fromEnum OptionType
cp --call put indicartor.
        df :: Double
df           = Double -> Double
forall a. Floating a => a -> a
exp (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ (-Double
r) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t
        sqrt2Pi :: Double
sqrt2Pi      = Double
2.506628274631000


euImpliedVolWith' Method
ChoKimKwak OptionType
cp (Forward Double
f) (Strike Double
k) (YearFrac Double
t) (Rate Double
r) (Premium Double
p) =
  let df :: Double
df              = Double -> Double
forall a. Floating a => a -> a
exp (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ (-Double
r) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t
      forwardPremium :: Double
forwardPremium  = Double
p Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
df
      straddlePremium :: Double
straddlePremium = case OptionType
cp of OptionType
Call -> Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forwardPremium Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
k)
                                   OptionType
Put  -> Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forwardPremium Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
k)
      nu' :: Double
nu'             = (Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
k) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
straddlePremium
      nu :: Double
nu              = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (-Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
forall a. IEEE a => a
epsilon) (Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
nu' (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
forall a. IEEE a => a
epsilon))
      eta :: Double
eta             | Double -> Double
forall a. Num a => a -> a
abs Double
nu Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
sqrtEpsilon = Double
1
                      | Bool
otherwise            = Double
nu Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double -> Double
forall a. Floating a => a -> a
atanh Double
nu)
      heta :: Double
heta            = Double -> Double
forall a. Floating a => a -> a
h Double
eta
  in Double -> Vol
Vol (Double -> Vol) -> Double -> Vol
forall a b. (a -> b) -> a -> b
$ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t)) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
straddlePremium Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
heta


euImpliedVolWith' RootFinding{Natural
(Double, Double, Double)
Tolerance
bracket :: (Double, Double, Double)
tol :: Tolerance
maxIter :: Natural
bracket :: Method -> (Double, Double, Double)
tol :: Method -> Tolerance
maxIter :: Method -> Natural
..}  OptionType
cp (Forward Double
forward) Strike
k YearFrac
t Rate
r (Premium Double
p) =
  let f :: Double -> Double
f Double
vol        = Double
p' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p  where
        (Premium Double
p') = Valuation -> Premium
vPremium (Valuation -> Premium) -> Valuation -> Premium
forall a b. (a -> b) -> a -> b
$ Bachelier -> YearFrac -> OptionType -> Strike -> Valuation
euOption Bachelier
b YearFrac
t OptionType
cp Strike
k
        b :: Bachelier
b            = Forward -> Rate -> Vol -> Bachelier
Bachelier (Double -> Forward
Forward Double
forward) Rate
r (Double -> Vol
Vol Double
vol)
      (Double
lb, Double
_, Double
ub) = (Double, Double, Double)
bracket
      root :: Root Double
root = RiddersParam
-> (Double, Double) -> (Double -> Double) -> Root Double
ridders (Int -> Tolerance -> RiddersParam
RiddersParam (Natural -> Int
forall a. Enum a => a -> Int
fromEnum Natural
maxIter) Tolerance
tol) (Double
lb, Double
ub) Double -> Double
f
  in case Root Double
root of (Root Double
vol)   -> Double -> Vol
Vol Double
vol
                  Root Double
NotBracketed -> [Char] -> Vol
forall a. HasCallStack => [Char] -> a
error [Char]
"not bracketed"
                  Root Double
SearchFailed -> [Char] -> Vol
forall a. HasCallStack => [Char] -> a
error [Char]
"search failed"

sqrtEpsilon :: Double
sqrtEpsilon = Double -> Double
forall a. Floating a => a -> a
sqrt Double
forall a. IEEE a => a
epsilon
h :: p -> p
h p
eta = p -> p
forall a. Floating a => a -> a
sqrt(p
eta) p -> p -> p
forall a. Num a => a -> a -> a
* (p
num p -> p -> p
forall a. Fractional a => a -> a -> a
/ p
den) where
  a0 :: p
a0          = p
3.994961687345134e-1
  a1 :: p
a1          = p
2.100960795068497e+1
  a2 :: p
a2          = p
4.980340217855084e+1
  a3 :: p
a3          = p
5.988761102690991e+2
  a4 :: p
a4          = p
1.848489695437094e+3
  a5 :: p
a5          = p
6.106322407867059e+3
  a6 :: p
a6          = p
2.493415285349361e+4
  a7 :: p
a7          = p
1.266458051348246e+4

  b0 :: p
b0          = p
1.000000000000000e+0
  b1 :: p
b1          = p
4.990534153589422e+1
  b2 :: p
b2          = p
3.093573936743112e+1
  b3 :: p
b3          = p
1.495105008310999e+3
  b4 :: p
b4          = p
1.323614537899738e+3
  b5 :: p
b5          = p
1.598919697679745e+4
  b6 :: p
b6          = p
2.392008891720782e+4
  b7 :: p
b7          = p
3.608817108375034e+3
  b8 :: p
b8          = -p
2.067719486400926e+2
  b9 :: p
b9          = p
1.174240599306013e+1

  num :: p
num = p
a0 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* (p
a1 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* (p
a2 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* (p
a3 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* (p
a4 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* (p
a5 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* (p
a6 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* p
a7))))))
  den :: p
den = p
b0 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* (p
b1 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* (p
b2 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* (p
b3 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* (p
b4 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* (p
b5 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* (p
b6 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* (p
b7 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* (p
b8 p -> p -> p
forall a. Num a => a -> a -> a
+ p
eta p -> p -> p
forall a. Num a => a -> a -> a
* p
b9))))))))