hasquant-0.7.0.0: test/hspec/QuantLib/Spec/PricingEngine.hs
-- Coverage for 'QuantLib.PricingEngine' entry points ported from proven-correct
-- @test/smoke/*.hs@ scripts that were never wired into @stack test --coverage@ (see
-- CLAUDE.md: coverage is only measured over @test\/hspec\/**@ + @test\/example\/**@).
-- Each block below preserves the source smoke script's own reasoning in its header
-- comment; only the assertion style changed (SmokeCheck's checkClose/checkEq/checkWith,
-- which 'error' on failure, become hspec 'shouldSatisfy'/'shouldBe').
module QuantLib.Spec.PricingEngine (spec) where
import Control.Monad(forM, forM_, when)
import Data.Time.Calendar(addDays, addGregorianYearsClip)
import Test.Hspec
import qualified Data.Vector.Storable as V
import Data.List(unzip6)
import Data.List.NonEmpty(fromList, NonEmpty((:|)))
import qualified QuantLib.Context as Context
import QuantLib.Time.Date
import QuantLib.Time.Calendar
import QuantLib.Time.Schedule
import QuantLib.InterestRate(Compounding(..), VolatilityType(..))
import QuantLib.Quote hiding(value)
import QuantLib.TermStructure.Yield
import QuantLib.TermStructure.Volatility hiding(atmStrike, atmVol)
import QuantLib.CashFlow(fixedDividend)
import qualified QuantLib.CashFlow as CF
import qualified QuantLib.Index.InterestRate as IR
import QuantLib.Instrument
import QuantLib.Instrument.Option hiding(deltaForward, vega, rho, dividendRho, strikeSensitivity, itmCashProbability)
import qualified QuantLib.Instrument.Option as Opt(rho, vega, dividendRho)
import QuantLib.Instrument.CapFloor(cap)
import QuantLib.Instrument.Swap(varianceSwap, vanillaSwap, floatingLeg, SwapType(..))
import QuantLib.Process hiding(thetaAt)
import QuantLib.Model hiding(setPricingEngine, value, discount)
import QuantLib.Math(RngTrait(..), StatisticsTrait(..), PolynomialType(..), BinomialTree(..), FdmScheme(..), boxedRealMatrix, ComplexLogFormula(..)
,SobolDirectionIntegers(..), realMatrixRows, realMatrixColumns, realMatrixData)
import QuantLib.Method(fdmBlackScholesMesher, fdmMesherComposite, fdmMesherLocations)
import QuantLib.PricingEngine hiding(alpha, delta, gamma, theta)
import qualified QuantLib.PricingEngine as Calc
import QuantLib.Spec.Helpers(closePrec)
normalPdf :: Double -> Double
normalPdf x = exp (-x * x / 2) / sqrt (2 * pi)
allFinite :: [Double] -> Bool
allFinite = all (\v -> not (isNaN v) && not (isInfinite v))
allTrees :: [BinomialTree]
allTrees =
[ JarrowRudd, CoxRossRubinstein, AdditiveEQPBinomialTree, Trigeorgis, Tian
, LeisenReimer, Joshi4, ExtendedJarrowRudd, ExtendedCoxRossRubinstein
, ExtendedAdditiveEQPBinomialTree, ExtendedTrigeorgis, ExtendedTian
, ExtendedLeisenReimer, ExtendedJoshi4
]
spec :: Spec
spec = do
-- Leans on invariants that hold for any
-- correct implementation of the Black-76/Bachelier formula (value = discount*(forward*alpha
-- + x*beta), with alpha/beta's *derivatives* independent of Call-vs-Put) rather than
-- hand-derived reference numbers: ctor-pair agreement, put-call parity, Call/Put agreement
-- on second-order greeks, agreement with the independently-implemented blackFormula/
-- bachelierBlackFormula free functions, and closed-form vanna/volga checks from
-- blackcalculator.cpp's own documented formula.
-- Coverage for the free functions bound from ql/pricingengines/blackformula.hpp -- both the
-- long-standing ones (blackFormula, blackCashItmProbability, blackImpliedStdDev*,
-- blackStdDevDerivative, blackVolDerivative, bachelierBlackFormula), which had none, and the
-- ones added alongside this block. Upstream's test-suite/blackformula.cpp asserts properties,
-- not cached numbers -- round-tripping an implied vol back through the pricing formula, and a
-- mean-value-theorem bracket on the forward derivative -- so those are ported as properties
-- rather than invented golden values. Where upstream has no test at all
-- (blackFormulaAssetItmProbability, blackFormulaStdDevSecondDerivative, and the Bachelier
-- derivative/probability pair) the check is a closed-form identity against a function bound
-- independently of it: for a Black call the premium decomposes as
-- discount*(F*N(d1) - K*N(d2)) with N(d1)/N(d2) exactly the asset/cash ITM probabilities and
-- dPremium/dF = discount*N(d1); the Bachelier analogues follow the same shape with a single N(d).
describe "blackFormula free functions" $ do
let fwd = 100.0; tte = 1.7; r = 0.1; df = exp (-r * tte)
vol = 0.3; sd = vol * sqrt tte
strikes = [50, 60, 70, 80, 90, 100, 110, 125, 150, 200, 300] :: [Double]
types = [Call, Put]
-- test-suite/blackformula.cpp: testRadoicicStefanicaImpliedVol. Same fixture (T=1.7, r=0.1,
-- forward=100, vol=0.3, the same 11 strikes) and the same 0.02 vol tolerance; the RS
-- approximation is closed-form, so this is an accuracy bound, not a solver convergence check.
it "blackImpliedStdDevApproximationRs recovers the generating vol to upstream's 0.02 tolerance" $
forM_ strikes $ \k -> forM_ types $ \t -> do
mv <- blackFormula t k fwd sd df 0.0
estSd <- blackImpliedStdDevApproximationRs t k fwd mv df 0.0
(estSd / sqrt tte) `shouldSatisfy` closePrec vol 0.02
-- test-suite/blackformula.cpp: testRadoicicStefanicaLowerBound, the figure-3.1 sweep from
-- "Tighter Bounds for Implied Volatility". Two separate claims: the approximation is within
-- 0.05 of the true stdDev, and (for a non-negligible premium) it is a *lower* bound.
it "blackImpliedStdDevApproximationRs stays a lower bound within 0.05 across the\
\ Gatheral-Matic-Radoicic-Stefanica sweep" $ do
let k = 1.2
strike = exp k * 1.0
forM_ [0.17, 0.18 .. 2.89 :: Double] $ \s -> do
c <- blackFormula Call strike 1.0 s 1.0 0.0
est <- blackImpliedStdDevApproximationRs Call strike 1.0 c 1.0 0.0
est `shouldSatisfy` (not . isNaN)
(s - est) `shouldSatisfy` (\e -> abs e <= 0.05)
when (c > 1e-6) $ (s - est) `shouldSatisfy` (>= 0.0)
-- test-suite/blackformula.cpp: testImpliedVolAdaptiveSuccessiveOverRelaxation. 'Nothing' for
-- the guess is upstream's Null<Real>(), i.e. "start from the RS approximation"; upstream
-- allows 10x the requested solver accuracy as the assertion tolerance.
it "blackImpliedStdDevLiRs inverts blackFormula to 10x its requested accuracy, over\
\ displacements" $ do
let tol = 1e-8
forM_ strikes $ \k -> forM_ types $ \t -> forM_ [0.0, 0.01, 0.05 :: Double] $ \displacement -> do
mv <- blackFormula t k fwd sd df displacement
impl <- blackImpliedStdDevLiRs t k fwd mv df displacement Nothing 1.0 tol 100
impl `shouldSatisfy` closePrec sd (10 * tol)
-- test-suite/blackformula.cpp: testChambersImpliedVol. Chambers-Nawalkha needs the ATM
-- premium as a second input; upstream measures its error moneyness-weighted and one-sided
-- (the approximation may undershoot freely, but must not overshoot by more than 5e-4).
it "blackImpliedStdDevChambers does not overshoot the true stdDev beyond upstream's\
\ moneyness-weighted 5e-4" $ do
let tol = 5.0e-4
displacements = [0.0, 0.001, 0.005, 0.01, 0.02] :: [Double]
fwds = [-0.001, 0.0, 0.005, 0.01, 0.02, 0.05] :: [Double]
ks = [-0.01, -0.005, -0.001, 0.0, 0.001, 0.005, 0.01, 0.02, 0.05, 0.1] :: [Double]
sds = [0.1, 0.15, 0.2, 0.3, 0.5, 0.6, 0.7, 0.8, 1.0, 1.5, 2.0] :: [Double]
discounts = [1.0, 0.95, 0.8, 1.1] :: [Double]
forM_ types $ \t -> forM_ displacements $ \displacement -> forM_ fwds $ \f ->
forM_ ks $ \k -> forM_ sds $ \s -> forM_ discounts $ \disc ->
when (f + displacement > 0.0 && k + displacement > 0.0) $ do
premium <- blackFormula t k f s disc displacement
atmPremium <- blackFormula t f f s disc displacement
iSd <- blackImpliedStdDevChambers t k f premium atmPremium disc displacement
let moneyness0 = (k + displacement) / (f + displacement)
moneyness = if moneyness0 > 1.0 then 1.0 / moneyness0 else moneyness0
((iSd - s) / s * moneyness) `shouldSatisfy` (<= tol)
-- test-suite/blackformula.cpp: testBachelierImpliedVol, including its two very tight
-- tolerances (1e-12 for the Choi approximation, 1e-15 for the Jaeckel formula) and its
-- zero-, deep-ITM- and deep-OTM-strike sweep in units of the standard deviation. Note both
-- take the time to expiry and return a *volatility*, unlike every Black function here.
it "bachelierImpliedVol and bachelierImpliedVolChoi invert bachelierBlackFormula to\
\ upstream's 1e-15 / 1e-12 tolerances" $ do
let bfwd = 1.0; bpvol = 0.01; btte = 10.0
bsd = bpvol * sqrt btte
bdisc = 0.95
forM_ [-3, -2, -1, -0.5, 0, 0.5, 1, 2, 3 :: Double] $ \i -> do
let k = bfwd - i * bpvol * sqrt btte
prem <- bachelierBlackFormula Call k bfwd bsd bdisc
choi <- bachelierImpliedVolChoi Call k bfwd btte prem bdisc
choi `shouldSatisfy` closePrec bpvol 1.0e-12
exact <- bachelierImpliedVol Call k bfwd btte prem bdisc
exact `shouldSatisfy` closePrec bpvol 1.0e-15
-- test-suite/blackformula.cpp: assertBlackFormulaForwardDerivative /
-- assertBachelierBlackFormulaForwardDerivative, including the zero-strike and zero-vol edge
-- cases each drives. The mean value theorem puts the bumped difference quotient between the
-- analytic derivative at the two ends of the bump, for any function monotonic across it.
it "blackForwardDerivative and bachelierForwardDerivative bracket their own bumped\
\ difference quotient, at zero strike and zero vol too" $ do
let bfwd = 1.0; btte = 10.0; bdisc = 0.95; displacement = 0.01
bump = 0.0001; eps = 1e-10
brackets d bd approx = max d bd + eps > approx && approx > min d bd - eps
forM_ [(0.01 :: Double, [0.1, 0.5, 1.0, 1.5, 2.0 :: Double]), (0.01, [0.0]), (0.0, [0.1, 1.0, 2.0])] $
\(bpvol, ks) -> do
let bsd = bpvol * sqrt btte
forM_ ks $ \k -> forM_ types $ \t -> do
d <- blackForwardDerivative t k bfwd bsd bdisc displacement
bd <- blackForwardDerivative t k (bfwd + bump) bsd bdisc displacement
p0 <- blackFormula t k bfwd bsd bdisc displacement
p1 <- blackFormula t k (bfwd + bump) bsd bdisc displacement
brackets d bd ((p1 - p0) / bump) `shouldBe` True
bad <- bachelierForwardDerivative t k bfwd bsd bdisc
bbd <- bachelierForwardDerivative t k (bfwd + bump) bsd bdisc
bp0 <- bachelierBlackFormula t k bfwd bsd bdisc
bp1 <- bachelierBlackFormula t k (bfwd + bump) bsd bdisc
brackets bad bbd ((bp1 - bp0) / bump) `shouldBe` True
-- No upstream test covers these, so each is pinned to a closed-form identity involving a
-- separately bound function. For a Black call, premium = discount*(F*N(d1) - K*N(d2)) with
-- N(d1) = blackAssetItmProbability and N(d2) = blackCashItmProbability, and
-- dPremium/dF = discount*N(d1) = blackForwardDerivative.
it "blackAssetItmProbability satisfies the Black decomposition and equals\
\ blackForwardDerivative/discount" $
forM_ strikes $ \k -> forM_ types $ \t -> do
pa <- blackAssetItmProbability t k fwd sd 0.0
pc <- blackCashItmProbability t k fwd sd 0.0
premium <- blackFormula t k fwd sd df 0.0
let sign = if t == Call then 1.0 else -1.0
(sign * df * (fwd * pa - k * pc)) `shouldSatisfy` closePrec premium (1e-12 * max 1 premium)
fd <- blackForwardDerivative t k fwd sd df 0.0
(fd / df) `shouldSatisfy` closePrec (sign * pa) 1e-12
-- Bachelier analogue of the identity above: premium = discount*((F-K)*N(d) + s*phi(d)) with
-- N(d) = bachelierAssetItmProbability and dPremium/dF = discount*N(d).
it "bachelierAssetItmProbability matches the Bachelier decomposition and\
\ bachelierForwardDerivative/discount" $ do
let bfwd = 1.0; bsd = 0.01 * sqrt 10.0; bdisc = 0.95
forM_ [0.9, 0.95, 1.0, 1.05, 1.1 :: Double] $ \k -> forM_ types $ \t -> do
let sign = if t == Call then 1.0 else -1.0
d = sign * (bfwd - k) / bsd
pa <- bachelierAssetItmProbability t k bfwd bsd
prem <- bachelierBlackFormula t k bfwd bsd bdisc
(bdisc * (sign * (bfwd - k) * pa + bsd * normalPdf d)) `shouldSatisfy` closePrec prem 1e-14
fd <- bachelierForwardDerivative t k bfwd bsd bdisc
(fd / bdisc) `shouldSatisfy` closePrec (sign * pa) 1e-13
-- Both StdDev derivative families are checked against a central difference of the function
-- they differentiate: blackStdDevDerivative against blackFormula, blackStdDevSecondDerivative
-- against blackStdDevDerivative, and bachelierStdDevDerivative against bachelierBlackFormula.
-- The central difference is O(h^2) accurate, hence the 1e-6-relative tolerance rather than a
-- tighter one; this also covers blackVolDerivative, which is the stdDev derivative scaled by
-- 1/sqrt(T) less the discount-rate term, i.e. sqrt(T)*blackStdDevDerivative here.
it "blackStdDevDerivative, blackStdDevSecondDerivative, blackVolDerivative and\
\ bachelierStdDevDerivative match central differences of what they differentiate" $ do
let h = 1e-5
forM_ strikes $ \k -> do
d1v <- blackStdDevDerivative k fwd sd df 0.0
pUp <- blackFormula Call k fwd (sd + h) df 0.0
pDn <- blackFormula Call k fwd (sd - h) df 0.0
d1v `shouldSatisfy` closePrec ((pUp - pDn) / (2 * h)) (1e-6 * max 1 (abs d1v))
d2v <- blackStdDevSecondDerivative k fwd sd df 0.0
dUp <- blackStdDevDerivative k fwd (sd + h) df 0.0
dDn <- blackStdDevDerivative k fwd (sd - h) df 0.0
d2v `shouldSatisfy` closePrec ((dUp - dDn) / (2 * h)) (1e-5 * max 1 (abs d2v))
volD <- blackVolDerivative k fwd sd tte df 0.0
volD `shouldSatisfy` closePrec (sqrt tte * d1v) (1e-9 * max 1 (abs volD))
bd <- bachelierStdDevDerivative k fwd sd df
bUp <- bachelierBlackFormula Call k fwd (sd + h) df
bDn <- bachelierBlackFormula Call k fwd (sd - h) df
bd `shouldSatisfy` closePrec ((bUp - bDn) / (2 * h)) (1e-6 * max 1 (abs bd))
-- The pre-existing iterative solver, which had no coverage either. Its cheap closed-form
-- sibling is the Brenner-Subrahmanyan\/Feinstein ATM formula extended by Corrado-Miller, so
-- it is only accurate near the money: measured across this fixture it is within 0.005 vol
-- over strikes 80..125 and degrades to ~0.18 at strike 300. The assertion is split
-- accordingly rather than pinned to one invented tolerance -- in the wings it only has to
-- stay a usable finite seed for the exact solver, which is all upstream uses it for.
it "blackImpliedStdDev inverts blackFormula, and blackImpliedStdDevApproximation is accurate\
\ near the money and a finite seed in the wings" $
forM_ strikes $ \k -> forM_ types $ \t -> do
mv <- blackFormula t k fwd sd df 0.0
impl <- blackImpliedStdDev t k fwd mv df 0.0 sd 1e-10 100
impl `shouldSatisfy` closePrec sd 1e-8
appr <- blackImpliedStdDevApproximation t k fwd mv df 0.0
appr `shouldSatisfy` (\v -> v > 0 && not (isNaN v) && not (isInfinite v))
when (k >= 80 && k <= 125) $ (appr / sqrt tte) `shouldSatisfy` closePrec vol 0.01
describe "BlackCalculator / BlackScholesCalculator / BachelierCalculator" $ do
let strike = 100.0; forward = 105.0; stdDev = 0.25; disc = 0.97
spot = 103.0; maturity = 2.0
it "BlackCalculator: ctors agree, put-call parity holds, Call/Put share second-order greeks,\
\ value matches blackFormula, vanna/volga match their closed forms" $ do
callBC <- blackCalculator (Strike Call strike) forward stdDev disc
putBC <- blackCalculator (Strike Put strike) forward stdDev disc
callBC2 <- blackCalculator (StrikePayoff (PlainVanilla (PlainVanillaPayoff Call strike))) forward stdDev disc
putBC2 <- blackCalculator (StrikePayoff (PlainVanilla (PlainVanillaPayoff Put strike))) forward stdDev disc
callVal <- value callBC
callVal2 <- value callBC2
callVal2 `shouldBe` callVal
putVal <- value putBC
putVal2 <- value putBC2
putVal2 `shouldBe` putVal
(callVal - putVal) `shouldSatisfy` closePrec (disc * (forward - strike)) 1e-10
callDF <- deltaForward callBC
putDF <- deltaForward putBC
(callDF - putDF) `shouldSatisfy` closePrec disc 1e-10
forM_ [gammaForward, (`vega` maturity), strikeGamma] $ \f -> do
c <- f callBC
p <- f putBC
p `shouldBe` c
refCallVal <- blackFormula Call strike forward stdDev disc 0.0
callVal `shouldSatisfy` closePrec refCallVal 1e-12
refPutVal <- blackFormula Put strike forward stdDev disc 0.0
putVal `shouldSatisfy` closePrec refPutVal 1e-12
let d1 = log (forward / strike) / stdDev + 0.5 * stdDev
d2 = d1 - stdDev
callVega <- vega callBC maturity
callVanna <- blackVanna callBC spot maturity
callVanna `shouldSatisfy` closePrec (-d2 / (spot * stdDev) * callVega) 1e-10
callVolga <- volga callBC maturity
callVolga `shouldSatisfy` closePrec (callVega * d1 * d2 / stdDev) 1e-10
it "BlackScholesCalculator: ctors agree, inherited GenBlackCalculator methods match the\
\ equivalent BlackCalculator exactly, and its own no-spot overrides match BlackCalculator's\
\ spot-taking versions at its stored spot" $ do
callBC <- blackCalculator (Strike Call strike) forward stdDev disc
callVal <- value callBC
let growth = 1.0
bscSpot = forward * disc / growth
callBSC <- blackScholesCalculator (Strike Call strike) bscSpot growth stdDev disc
callBSC2 <- blackScholesCalculator (StrikePayoff (PlainVanilla (PlainVanillaPayoff Call strike))) bscSpot growth stdDev disc
callBSCVal <- Calc.value callBSC
callBSCVal2 <- Calc.value callBSC2
callBSCVal2 `shouldBe` callBSCVal
callBSCVal `shouldBe` callVal
let checkInherited fBSC fBC = do
fromBSC <- fBSC callBSC
fromBC <- fBC callBC
fromBSC `shouldBe` fromBC
checkInherited Calc.deltaForward Calc.deltaForward
checkInherited (`Calc.rho` maturity) (`Calc.rho` maturity)
checkInherited (`Calc.dividendRho` maturity) (`Calc.dividendRho` maturity)
checkInherited Calc.strikeSensitivity Calc.strikeSensitivity
checkInherited Calc.strikeGamma Calc.strikeGamma
checkInherited (`Calc.vega` maturity) (`Calc.vega` maturity)
checkInherited (`Calc.volga` maturity) (`Calc.volga` maturity)
checkInherited Calc.itmAssetProbability Calc.itmAssetProbability
checkInherited Calc.itmCashProbability Calc.itmCashProbability
checkInherited Calc.alpha Calc.alpha
checkInherited Calc.beta Calc.beta
bscVanna <- blackVanna callBSC bscSpot maturity
bcVannaAtBscSpot <- blackVanna callBC bscSpot maturity
bscVanna `shouldBe` bcVannaAtBscSpot
bscDelta <- blackScholesDelta callBSC
bcDeltaAtSpot <- Calc.delta callBC bscSpot
bscDelta `shouldBe` bcDeltaAtSpot
bscElasticity <- blackScholesElasticity callBSC
bcElasticityAtSpot <- Calc.elasticity callBC bscSpot
bscElasticity `shouldBe` bcElasticityAtSpot
bscGamma <- blackScholesGamma callBSC
bcGammaAtSpot <- Calc.gamma callBC bscSpot
bscGamma `shouldBe` bcGammaAtSpot
bscTheta <- blackScholesTheta callBSC maturity
bcThetaAtSpot <- Calc.theta callBC bscSpot maturity
bscTheta `shouldBe` bcThetaAtSpot
bscThetaPerDay <- blackScholesThetaPerDay callBSC maturity
bcThetaPerDayAtSpot <- Calc.thetaPerDay callBC bscSpot maturity
bscThetaPerDay `shouldBe` bcThetaPerDayAtSpot
it "BachelierCalculator: ctors agree, put-call parity holds, Call/Put share second-order\
\ greeks, value matches bachelierBlackFormula, vanna/volga/vega match their closed forms" $ do
-- BachelierCalculator's stdDev is an *absolute* normal-model volatility (e.g. rate
-- points), not a Black-style relative one -- a stdDev of realistic magnitude for the
-- forward/strike spread is needed or every second-order greek degenerates to ~1e-87
-- in the tail, passing every check without exercising the formula.
let bachelierStdDev = 8.0
callNC <- bachelierCalculator (Strike Call strike) forward bachelierStdDev disc
putNC <- bachelierCalculator (Strike Put strike) forward bachelierStdDev disc
callNC2 <- bachelierCalculator (StrikePayoff (PlainVanilla (PlainVanillaPayoff Call strike))) forward bachelierStdDev disc
putNC2 <- bachelierCalculator (StrikePayoff (PlainVanilla (PlainVanillaPayoff Put strike))) forward bachelierStdDev disc
callNVal <- Calc.value callNC
callNVal2 <- Calc.value callNC2
callNVal2 `shouldBe` callNVal
putNVal <- Calc.value putNC
putNVal2 <- Calc.value putNC2
putNVal2 `shouldBe` putNVal
(callNVal - putNVal) `shouldSatisfy` closePrec (disc * (forward - strike)) 1e-10
callNDF <- Calc.deltaForward callNC
putNDF <- Calc.deltaForward putNC
(callNDF - putNDF) `shouldSatisfy` closePrec disc 1e-10
forM_ [ Calc.gammaForward, (`Calc.vega` maturity)
, Calc.strikeGamma, (`bachelierVanna` maturity), (`Calc.volga` maturity)
] $ \f -> do
c <- f callNC
p <- f putNC
p `shouldBe` c
refCallNVal <- bachelierBlackFormula Call strike forward bachelierStdDev disc
callNVal `shouldSatisfy` closePrec refCallNVal 1e-12
refPutNVal <- bachelierBlackFormula Put strike forward bachelierStdDev disc
putNVal `shouldSatisfy` closePrec refPutNVal 1e-12
let d = (forward - strike) / bachelierStdDev
nd = normalPdf d
callNVega <- Calc.vega callNC maturity
callNVanna <- bachelierVanna callNC maturity
callNVanna `shouldSatisfy` closePrec (-d * nd * sqrt maturity / bachelierStdDev) 1e-9
callNVolga <- Calc.volga callNC maturity
callNVolga `shouldSatisfy` closePrec (d * d / bachelierStdDev * callNVega) 1e-9
callNVega `shouldSatisfy` closePrec (disc * sqrt maturity * nd) 1e-9
-- volatility(strike)/variance(strike) must
-- exactly match the already-bound unsafeShiftedSabrVolatility formula (enum-dispatched
-- through a C-side VolatilityType cast, same class of bug as the CPIInterpolationType
-- incident); the Date- and Time-based ctors of SabrSmileSection/NoArbSabrSmileSection must
-- agree with each other and (for NoArb) differ from the plain SABR smile; and
-- SabrInterpolatedSmileSection must calibrate back to the SABR parameters that generated its
-- input vols.
describe "SabrSmileSection / NoArbSabrSmileSection / SabrInterpolatedSmileSection" $ do
let forward = 0.03; expiry = 5.0; alpha_ = 0.04; beta_ = 0.5; nu = 0.4; rho_ = -0.2; shift = 0.0
it "matches unsafeShiftedSabrVolatility for both VolatilityType cases, and Normal /= ShiftedLognormal" $
Context.keepingSettingsGc $ do
forM_ [ShiftedLognormal, Normal] $ \volType -> do
section <- sabrSmileSection expiry forward alpha_ beta_ nu rho_ shift volType
forM_ [0.01, 0.02, 0.03, 0.04, 0.05 :: Double] $ \strike -> do
got <- smileSectionVolatility section strike
expected <- unsafeShiftedSabrVolatility strike forward expiry alpha_ beta_ nu rho_ shift volType
got `shouldBe` expected
var <- smileSectionVariance section strike
var `shouldSatisfy` closePrec (expected * expected * expiry) 1e-12
volShiftedLognormal <- sabrSmileSection expiry forward alpha_ beta_ nu rho_ shift ShiftedLognormal
volAtAtm1 <- smileSectionVolatility volShiftedLognormal forward
volNormal <- sabrSmileSection expiry forward alpha_ beta_ nu rho_ shift Normal
volAtAtm2 <- smileSectionVolatility volNormal forward
volAtAtm1 `shouldNotBe` volAtAtm2
it "SabrSmileSection/NoArbSabrSmileSection Date- and Time-based ctors agree, and NoArb\
\ differs from the plain SabrSmileSection" $
Context.keepingSettingsGc $ do
refDate <- today
Context.setEvaluationDate (Just refDate)
let expiryDays = 1826 :: Int -- ~5y in actual days
expiryFromDays = fromIntegral expiryDays / 365.0 :: Double
optionDate <- addPeriod refDate (expiryDays, Days)
act365 <- dayCounter Actual365FixedStandard
sectionByTime <- sabrSmileSection expiryFromDays forward alpha_ beta_ nu rho_ shift ShiftedLognormal
sectionByDate <- sabrSmileSectionAtDate optionDate act365 (Just refDate) forward alpha_ beta_ nu rho_ shift ShiftedLognormal
forM_ [0.01, 0.02, 0.03, 0.04, 0.05 :: Double] $ \strike -> do
volT <- smileSectionVolatility sectionByTime strike
volD <- smileSectionVolatility sectionByDate strike
volD `shouldSatisfy` closePrec volT 1e-12
noArbByTime <- noArbSabrSmileSection (RateAtTime expiryFromDays) forward alpha_ beta_ nu rho_ shift ShiftedLognormal
noArbByDate <- noArbSabrSmileSection (RateAtDate optionDate act365) forward alpha_ beta_ nu rho_ shift ShiftedLognormal
differences <- forM [0.01, 0.02, 0.03, 0.04, 0.05 :: Double] $ \strike -> do
volT <- smileSectionVolatility noArbByTime strike
volD <- smileSectionVolatility noArbByDate strike
volD `shouldSatisfy` closePrec volT 1e-12
plainVol <- smileSectionVolatility sectionByTime strike
return (volT /= plainVol)
or differences `shouldBe` True
it "SabrInterpolatedSmileSection calibrates back to the generating SABR parameters,\
\ including through the AsSmileSection upcast" $
Context.keepingSettingsGc $ do
let strikes = [0.01, 0.02, 0.03, 0.04, 0.05]
refVols <- mapM (\k -> unsafeShiftedSabrVolatility k forward expiry alpha_ beta_ nu rho_ shift ShiftedLognormal) strikes
atmVol <- unsafeShiftedSabrVolatility forward forward expiry alpha_ beta_ nu rho_ shift ShiftedLognormal
now <- today
Context.setEvaluationDate (Just now)
optionDate <- addPeriod now (round (expiry * 365) :: Int, Days)
forwardQuote <- simpleQuote forward
atmVolQ <- simpleQuote atmVol
refVolQuotes <- mapM simpleQuote refVols
interp <- sabrInterpolatedSmileSection optionDate forwardQuote (fromList $ zip strikes refVolQuotes) False atmVolQ
alpha_ beta_ nu rho_ defaultSabrInterpolatedSmileSectionOpts
rms <- sabrRmsError interp
maxErr <- sabrMaxError interp
rms `shouldSatisfy` (< 1e-6)
maxErr `shouldSatisfy` (< 1e-6)
-- the upcast escape hatch: volatility through the generic SmileSection interface (only
-- reachable this way now that the concrete type has no smileSectionVolatility of its own).
generic <- sabrAsSmileSection interp
forM_ (zip strikes refVols) $ \(k, expected) -> do
got <- smileSectionVolatility generic k
got `shouldSatisfy` closePrec expected 1e-6
-- SviInterpolatedSmileSection has no upstream test-suite fixture of its own (only
-- SviSmileSection's direct-parameter construction is exercised in test-suite/svivolatility.cpp),
-- so this is a round-trip against the already-bound sviSmileSection: generate vols at known SVI
-- parameters, feed them back in as quotes, and require calibration to reproduce them closely.
-- Same shape as the SabrInterpolatedSmileSection check above, plus one asymmetric-fixed-flag
-- case (only mIsFixed set) to pin the argument order across the six adjacent Bools
-- (aIsFixed/bIsFixed/sigmaIsFixed/rhoIsFixed/mIsFixed/vegaWeighted) -- an all-False round-trip
-- alone wouldn't catch a transposition there.
describe "SviInterpolatedSmileSection" $ do
let a_ = -0.0666; b_ = 0.229; sigma_ = 0.337; rho_ = 0.439; m_ = 0.193
strikes = [0.5, 0.8, 1.0, 1.2, 1.5, 2.0]
forward = 1.0
it "calibrates back to the generating SVI parameters, including through the\
\ AsSmileSection upcast" $
Context.keepingSettingsGc $ do
refDate <- today
Context.setEvaluationDate (Just refDate)
optionDate <- addPeriod refDate (180, Days)
act365 <- dayCounter Actual365FixedStandard
generating <- sviSmileSection optionDate forward a_ b_ sigma_ rho_ m_ act365
refVols <- mapM (smileSectionVolatility generating) strikes
atmVol <- smileSectionVolatility generating forward
forwardQuote <- simpleQuote forward
atmVolQ <- simpleQuote atmVol
refVolQuotes <- mapM simpleQuote refVols
interp <- sviInterpolatedSmileSection optionDate forwardQuote (fromList $ zip strikes refVolQuotes)
False atmVolQ a_ b_ sigma_ rho_ m_ False False False False False True Nothing Nothing act365
rms <- sviRmsError interp
maxErr <- sviMaxError interp
rms `shouldSatisfy` (< 1e-6)
maxErr `shouldSatisfy` (< 1e-6)
-- the upcast escape hatch: volatility through the generic SmileSection interface (only
-- reachable this way now that the concrete type has no smileSectionVolatility of its own).
genericSection <- sviAsSmileSection interp
forM_ (zip strikes refVols) $ \(k, expected) -> do
got <- smileSectionVolatility genericSection k
got `shouldSatisfy` closePrec expected 1e-6
it "mIsFixed pins m at the seed value while the other four parameters still calibrate" $
Context.keepingSettingsGc $ do
refDate <- today
Context.setEvaluationDate (Just refDate)
optionDate <- addPeriod refDate (180, Days)
act365 <- dayCounter Actual365FixedStandard
generating <- sviSmileSection optionDate forward a_ b_ sigma_ rho_ m_ act365
refVols <- mapM (smileSectionVolatility generating) strikes
atmVol <- smileSectionVolatility generating forward
forwardQuote <- simpleQuote forward
atmVolQ <- simpleQuote atmVol
refVolQuotes <- mapM simpleQuote refVols
interp <- sviInterpolatedSmileSection optionDate forwardQuote (fromList $ zip strikes refVolQuotes)
False atmVolQ a_ b_ sigma_ rho_ m_ False False False False True True Nothing Nothing act365
calibratedM <- sviM interp
calibratedM `shouldBe` m_
rms <- sviRmsError interp
rms `shouldSatisfy` (< 1e-6)
-- NoArbSabrInterpolatedSmileSection has no upstream test-suite fixture of its own either (only
-- NoArbSabrInterpolation's direct construction is exercised in test-suite/interpolations.cpp),
-- so this is the same round-trip shape as SviInterpolatedSmileSection above: generate vols at
-- known no-arb SABR parameters via the already-bound noArbSabrSmileSection, feed them back in
-- as quotes, and require calibration to reproduce them closely.
describe "NoArbSabrInterpolatedSmileSection" $
it "calibrates back to the generating no-arb SABR parameters, including through the\
\ AsSmileSection upcast" $
Context.keepingSettingsGc $ do
let forward = 0.03; alpha_ = 0.04; beta_ = 0.5; nu = 0.4; rho_ = -0.2
strikes = [0.01, 0.02, 0.03, 0.04, 0.05]
refDate <- today
Context.setEvaluationDate (Just refDate)
optionDate <- addPeriod refDate (round (5.0 * 365 :: Double) :: Int, Days)
act365 <- dayCounter Actual365FixedStandard
generating <- noArbSabrSmileSection (RateAtDate optionDate act365) forward alpha_ beta_ nu rho_ 0 ShiftedLognormal
refVols <- mapM (smileSectionVolatility generating) strikes
atmVol <- smileSectionVolatility generating forward
forwardQuote <- simpleQuote forward
atmVolQ <- simpleQuote atmVol
refVolQuotes <- mapM simpleQuote refVols
interp <- noArbSabrInterpolatedSmileSection optionDate forwardQuote (fromList $ zip strikes refVolQuotes)
False atmVolQ alpha_ beta_ nu rho_ False False False False True Nothing Nothing act365
rms <- noArbSabrRmsError interp
maxErr <- noArbSabrMaxError interp
rms `shouldSatisfy` (< 1e-6)
maxErr `shouldSatisfy` (< 1e-6)
-- the upcast escape hatch: volatility through the generic SmileSection interface (only
-- reachable this way now that the concrete type has no smileSectionVolatility of its own).
genericSection <- noArbSabrAsSmileSection interp
forM_ (zip strikes refVols) $ \(k, expected) -> do
got <- smileSectionVolatility genericSection k
got `shouldSatisfy` closePrec expected 1e-6
-- ported from test-suite/zabr.cpp::testConsistency: at gamma=1, ZabrSmileSection<Evaluation>
-- must (nearly) coincide with the Hagan 2002 SABR closed form (already bound as
-- sabrSmileSection), across all four ZabrEvaluation modes. fdRefinement=2 (not upstream's
-- default of 5) for the ZabrFullFd case matches upstream's own speed-up for this test; the
-- strike grid here is coarser than upstream's 7000-point sweep -- upstream's own tolerance
-- (1e-4 absolute on price) doesn't need that density to be meaningfully checked.
describe "ZabrSmileSection vs. SabrSmileSection (gamma=1) consistency" $
it "optionPrice agrees closely across a strike sweep, for every ZabrEvaluation mode" $ do
let tau = 5.0; forward = 0.03; alpha_ = 0.08; beta_ = 0.70; nu = 0.20; rho_ = -0.30
gamma_ = 1.0
strikes = [0.0001, 0.0071 .. 0.70] :: [Double]
tol = 1e-4
sabr <- sabrSmileSection tau forward alpha_ beta_ nu rho_ 0 ShiftedLognormal
forM_ [(ZabrShortMaturityLognormal, 5), (ZabrShortMaturityNormal, 5),
(ZabrLocalVolatility, 5), (ZabrFullFd, 2)] $ \(evaluation, fdRefinement) -> do
zabr <- zabrSmileSection evaluation (RateAtTime tau) forward alpha_ beta_ nu rho_ gamma_ [] fdRefinement
forM_ strikes $ \k -> do
c0 <- smileSectionOptionPrice sabr k Call 1.0
z <- smileSectionOptionPrice zabr k Call 1.0
z `shouldSatisfy` closePrec c0 tol
-- Verify that the StatisticsTrait axis reaches each Monte Carlo engine's
-- second template parameter. Nothing in the type system catches a StatisticsTrait value
-- being silently ignored (a copy-paste slip could alias all four cases to the same
-- instantiation), so each case is constructed and priced under a fixed nonzero seed.
describe "MC engine StatisticsTrait dispatch" $ do
let refDate = 15 `january` 2024
maturity = 15 `january` 2025
stats = [Statistics, GaussianStatistics, GeneralStatistics, IncrementalStatistics]
flatProc = do
dc <- dayCounter Actual365FixedStandard
cal <- calendar Null
spot <- simpleQuote 100.0
rfQ <- simpleQuote 0.03
rf <- flatForward (ReferenceDate refDate) rfQ dc Continuous Annual
volQ <- simpleQuote 0.20
vol <- blackConstantVol (CalendarReferenceDate refDate) cal volQ dc
blackScholesProcess spot rf vol EulerDiscretization False
europeanNpvUnder stat = do
proc' <- flatProc
let payoff = PlainVanilla (PlainVanillaPayoff Call 100)
exercise = European (EuropeanExercise maturity)
opt <- vanillaOption payoff exercise
eng <- mcEuropeanEngine PseudoRandom stat proc' (Just 1) Nothing False False (Just 1000) Nothing Nothing 42
setPricingEngine opt eng
npv opt
americanNpvUnder stat = do
proc' <- flatProc
let payoff = PlainVanilla (PlainVanillaPayoff Put 100)
exercise = American Nothing maturity False
opt <- vanillaOption payoff exercise
eng <- mcAmericanEngine PseudoRandom stat proc' (Just 50) Nothing True False Nothing (Just 0.02) Nothing 42
2 Monomial (Just 512) Nothing Nothing
setPricingEngine opt eng
npv opt
varianceSwapNpvUnder stat = do
proc' <- flatProc
sw <- varianceSwap Long 0.04 10000 refDate maturity
eng <- mcVarianceSwapEngine PseudoRandom stat proc'
(Just 52) Nothing False False (Just 1000) Nothing Nothing 42
setPricingEngine sw eng
npv sw
it "mcEuropeanEngine: every StatisticsTrait case constructs and prices finitely" $
Context.keepingSettingsGc $ do
Context.setEvaluationDate (Just refDate)
results <- mapM europeanNpvUnder stats
allFinite results `shouldBe` True
it "mcAmericanEngine: every StatisticsTrait case constructs and prices finitely" $
Context.keepingSettingsGc $ do
Context.setEvaluationDate (Just refDate)
results <- mapM americanNpvUnder stats
allFinite results `shouldBe` True
it "mcVarianceSwapEngine: every StatisticsTrait case constructs, prices finitely, and\
\ (same seed/process/timesteps) agrees closely across accumulators" $
Context.keepingSettingsGc $ do
Context.setEvaluationDate (Just refDate)
results <- mapM varianceSwapNpvUnder stats
allFinite results `shouldBe` True
let mn = minimum results; mx = maximum results
((mx - mn) / abs mn) `shouldSatisfy` (< 1e-6)
-- Check BlackDeltaCalculator against cached rows
-- from blackdeltacalculator.cpp's testDeltaValues, VannaVolgaBarrierEngine and
-- AnalyticDoubleBarrierEngine against cached NPVs from barrieroption.cpp/
-- doublebarrieroption.cpp, AnalyticPartialTimeBarrierOptionEngine against
-- partialtimebarrieroption.cpp, and every remaining barrier/double-barrier engine
-- (AnalyticBinaryBarrierEngine, FdBlackScholesBarrierEngine, BinomialBarrierEngine over
-- every BinomialTree, VannaVolgaDoubleBarrierEngine, BinomialDoubleBarrierEngine over every
-- BinomialTree, MCDoubleBarrierEngine over every RngTrait) constructed and priced at least
-- once, guarding against the binomial-tree/RNG switch tables being incomplete or
-- mis-ordered (the CPIInterpolationType-style gotcha in CLAUDE.md).
describe "Barrier / DoubleBarrier engines" $ do
it "BlackDeltaCalculator reproduces cached deltas per DeltaType, and AtmType/DeltaType do not alias" $ do
let checkOne ot dt spot dDf fDf stdDev strike expected = do
calc <- blackDeltaCalculator ot dt spot dDf fDf stdDev
delta_ <- deltaFromStrike calc strike
delta_ `shouldSatisfy` closePrec expected 0.15
checkOne Call Spot 1.421 0.997306 0.992266 0.1180654 1.608080 0.15
checkOne Call PaSpot 1.421 0.997306 0.992266 0.1180654 1.600545 0.15
checkOne Call Fwd 1.421 0.997306 0.992266 0.1180654 1.609029 0.15
checkOne Call PaFwd 1.421 0.997306 0.992266 0.1180654 1.601550 0.15
calc <- blackDeltaCalculator Call Spot 1.30265 0.99979 0.98508 0.11638
strikes <- forM [AtmSpot, AtmFwd, AtmDeltaNeutral] $ atmStrike calc
let dedup [] = []
dedup (x:xs) = x : dedup (filter (\y -> abs (y - x) > 1.0e-6) xs)
length (dedup strikes) `shouldBe` 3
fwdCalc <- blackDeltaCalculator Call Fwd 1.30265 0.99979 0.98508 0.11638
fwdStrike <- atmStrike fwdCalc AtmFwd
putCall50Strike <- atmStrike fwdCalc AtmPutCall50
abs (fwdStrike - putCall50Strike) `shouldSatisfy` (> 1.0e-6)
it "VannaVolgaBarrierEngine reproduces the cached UpOut EUR call value from barrieroption.cpp" $
Context.keepingSettingsGc $ do
let today' = 5 `march` 2013
barrier = 1.5; strike = 1.13321; s = 1.30265; q = 0.0003541; r = 0.0033871; t = 1 :: Double
vol25Put = 0.10087; volAtm = 0.08925; vol25Call = 0.08463; vol = 0.11638
Context.setEvaluationDate (Just today')
dc <- dayCounter Actual365FixedStandard
spotQ <- simpleQuote s
qQ <- simpleQuote q
rQ <- simpleQuote r
qTS <- flatForward (ReferenceDate today') qQ dc Continuous Annual
rTS <- flatForward (ReferenceDate today') rQ dc Continuous Annual
vol25PutQ <- simpleQuote vol25Put
volAtmQ <- simpleQuote volAtm
vol25CallQ <- simpleQuote vol25Call
qDisc <- discount qTS (TimePoint t) False
rDisc <- discount rTS (TimePoint t) False
let forward = s * qDisc / rDisc
bsVanillaPrice <- blackFormula Call strike forward (vol * sqrt t) rDisc 0.0
volAtmQuote <- atmVolQuote volAtmQ Fwd t AtmDeltaNeutral
vol25PutQuote <- deltaVolQuote (-0.25) vol25PutQ t Fwd
vol25CallQuote <- deltaVolQuote 0.25 vol25CallQ t Fwd
let payoff = PlainVanilla (PlainVanillaPayoff Call strike)
exercise = European (EuropeanExercise (addDays (round (t * 365)) today'))
opt <- barrierOption UpOut barrier 0 payoff exercise
optInst <- asOneAssetOption opt
engine <- vannaVolgaBarrierEngine volAtmQuote vol25PutQuote vol25CallQuote spotQ rTS qTS True bsVanillaPrice
setPricingEngine optInst engine
price <- npv optInst
price `shouldSatisfy` closePrec 0.148127 2.0e-3
it "AnalyticDoubleBarrierEngine reproduces the cached KnockOut call value from doublebarrieroption.cpp" $
Context.keepingSettingsGc $ do
let today' = 1 `january` 2020
barrierLo = 50.0; barrierHi = 150.0; strike = 100.0
s = 100.0; q = 0.0; r = 0.1; t = 0.25 :: Double; vol = 0.15
Context.setEvaluationDate (Just today')
dc <- dayCounter (Actual360 False)
spotQ <- simpleQuote s
qQ <- simpleQuote q
rQ <- simpleQuote r
volQ <- simpleQuote vol
qTS <- flatForward (ReferenceDate today') qQ dc Continuous Annual
rTS <- flatForward (ReferenceDate today') rQ dc Continuous Annual
tgt <- calendar TARGET
volTS <- blackConstantVol (CalendarReferenceDate today') tgt volQ dc
proc <- blackScholesMertonProcess spotQ qTS rTS volTS EulerDiscretization False
let payoff = PlainVanilla (PlainVanillaPayoff Call strike)
exercise = European (EuropeanExercise (addDays (round (t * 360)) today'))
opt <- doubleBarrierOption KnockOut barrierLo barrierHi 0 payoff exercise
optInst <- asOneAssetOption opt
engine <- analyticDoubleBarrierEngine proc 5
setPricingEngine optInst engine
price <- npv optInst
price `shouldSatisfy` closePrec 4.3515 1.0e-3
it "AnalyticPartialTimeBarrierOptionEngine reproduces the cached DownOut/EndB1 value from\
\ partialtimebarrieroption.cpp" $
Context.keepingSettingsGc $ do
let today' = 1 `january` 2020
barrier = 100.0; rebate = 0.0; strike = 90.0
s = 95.0; q = 0.0; r = 0.1 :: Double; vol = 0.25
Context.setEvaluationDate (Just today')
dc <- dayCounter (Actual360 False)
spotQ <- simpleQuote s
qQ <- simpleQuote q
rQ <- simpleQuote r
volQ <- simpleQuote vol
qTS <- flatForward (ReferenceDate today') qQ dc Continuous Annual
rTS <- flatForward (ReferenceDate today') rQ dc Continuous Annual
tgt <- calendar TARGET
volTS <- blackConstantVol (CalendarReferenceDate today') tgt volQ dc
proc <- blackScholesMertonProcess spotQ qTS rTS volTS EulerDiscretization False
let payoff = PlainVanilla (PlainVanillaPayoff Call strike)
exercise = European (EuropeanExercise (addDays 360 today'))
coverEventDate = addDays 1 today'
opt <- partialTimeBarrierOption DownOut EndB1 barrier rebate coverEventDate payoff exercise
engine <- analyticPartialTimeBarrierOptionEngine proc
setPricingEngine opt engine
price <- npv opt
price `shouldSatisfy` closePrec 0.0393 1.0e-4
it "every remaining barrier/double-barrier engine constructs and prices finitely, over every\
\ BinomialTree/RngTrait case" $
Context.keepingSettingsGc $ do
let today' = 1 `january` 2020 :: Day
Context.setEvaluationDate (Just today')
dc <- dayCounter Actual365FixedStandard
spotQ <- simpleQuote 100
qQ <- simpleQuote 0.01
rQ <- simpleQuote 0.02
volQ <- simpleQuote 0.2
qTS <- flatForward (ReferenceDate today') qQ dc Continuous Annual
rTS <- flatForward (ReferenceDate today') rQ dc Continuous Annual
tgt <- calendar TARGET
volTS <- blackConstantVol (CalendarReferenceDate today') tgt volQ dc
proc <- blackScholesMertonProcess spotQ qTS rTS volTS EulerDiscretization False
let payoff = PlainVanilla (PlainVanillaPayoff Call 100)
exercise = European (EuropeanExercise (addDays 180 today'))
binaryPayoff = CashOrNothing Call 100 10
americanExercise = American Nothing (addDays 180 today') True
binOpt <- barrierOption UpOut 130 0 binaryPayoff americanExercise >>= asOneAssetOption
analyticBinaryBarrierEngine proc >>= setPricingEngine binOpt
binNpv <- npv binOpt
binNpv `shouldSatisfy` (\v -> not (isNaN v) && not (isInfinite v))
barOpt <- barrierOption UpOut 130 0 payoff exercise >>= asOneAssetOption
fdBlackScholesBarrierEngine proc 100 100 0 Douglas False 0.0 >>= setPricingEngine barOpt
fdNpv <- npv barOpt
fdNpv `shouldSatisfy` (\v -> not (isNaN v) && not (isInfinite v))
binomialNpvs <- forM allTrees $ \ty -> do
barOpt' <- barrierOption UpOut 130 0 payoff exercise >>= asOneAssetOption
binomialBarrierEngine ty proc 200 0 >>= setPricingEngine barOpt'
npv barOpt'
allFinite binomialNpvs `shouldBe` True
dblOpt <- doubleBarrierOption KnockOut 70 130 0 payoff exercise >>= asOneAssetOption
let dblT = 180 / 365 :: Double
qDisc <- discount qTS (TimePoint dblT) False
rDisc <- discount rTS (TimePoint dblT) False
bsVanillaPrice <- blackFormula Call 100 (100 * qDisc / rDisc) (0.2 * sqrt dblT) rDisc 0.0
volAtmQ <- simpleQuote 0.2
vol25PutQ <- simpleQuote 0.22
vol25CallQ <- simpleQuote 0.18
volAtmQuote <- atmVolQuote volAtmQ Fwd dblT AtmDeltaNeutral
vol25PutQuote <- deltaVolQuote (-0.25) vol25PutQ dblT Fwd
vol25CallQuote <- deltaVolQuote 0.25 vol25CallQ dblT Fwd
vannaVolgaDoubleBarrierEngine volAtmQuote vol25PutQuote vol25CallQuote spotQ rTS qTS True bsVanillaPrice 5
>>= setPricingEngine dblOpt
vvdbNpv <- npv dblOpt
vvdbNpv `shouldSatisfy` (\v -> not (isNaN v) && not (isInfinite v))
binomialDblNpvs <- forM allTrees $ \ty -> do
dblOpt' <- doubleBarrierOption KnockOut 70 130 0 payoff exercise >>= asOneAssetOption
binomialDoubleBarrierEngine ty proc 200 >>= setPricingEngine dblOpt'
npv dblOpt'
allFinite binomialDblNpvs `shouldBe` True
mcDblNpvs <- forM [PseudoRandom, LowDiscrepancy] $ \rng -> do
dblOpt' <- doubleBarrierOption KnockOut 70 130 0 payoff exercise >>= asOneAssetOption
mcDoubleBarrierEngine rng Statistics proc (Just 20) Nothing False False (Just 1024) Nothing Nothing 42
>>= setPricingEngine dblOpt'
npv dblOpt'
allFinite mcDblNpvs `shouldBe` True
-- fdHestonBarrierEngine/fdHestonDoubleBarrierEngine are otherwise untested: no upstream
-- golden fixture was found reachable through hasquant's current Heston-model bindings,
-- so (matching the "constructs and prices finitely" precedent just above) this just
-- confirms both actually price rather than crash/NaN.
hp <- hestonProcess rTS (Just qTS) spotQ 0.04 2.0 0.04 0.5 (-0.5) QuadraticExponentialMartingale
hm <- hestonModel hp
fdHBarOpt <- barrierOption UpOut 130 0 payoff exercise >>= asOneAssetOption
fdHestonBarrierEngine hm [] 20 100 20 0 Douglas Nothing 1.0 >>= setPricingEngine fdHBarOpt
fdHBarNpv <- npv fdHBarOpt
fdHBarNpv `shouldSatisfy` (\v -> not (isNaN v) && not (isInfinite v))
fdHDblOpt <- doubleBarrierOption KnockOut 70 130 0 payoff exercise >>= asOneAssetOption
fdHestonDoubleBarrierEngine hm 20 100 20 0 Douglas Nothing 1.0 >>= setPricingEngine fdHDblOpt
fdHDblNpv <- npv fdHDblOpt
fdHDblNpv `shouldSatisfy` (\v -> not (isNaN v) && not (isInfinite v))
it "mcBarrierEngine reproduces the cached DownIn call value from barrieroption.cpp" $
-- cached reference from testHaugValues's Barrier::DownIn row (barrier=90, vol=0.10,
-- expected 0.07187): analytic (AnalyticBarrierEngine) and MC (MakeMCBarrierEngine) are
-- both checked against the same literal there, so the MC engine is checked here too,
-- at upstream's own relative tolerance.
Context.keepingSettingsGc $ do
let today' = 1 `january` 2020
expected = 0.07187 :: Double
Context.setEvaluationDate (Just today')
dc <- dayCounter (Actual360 False)
spotQ <- simpleQuote 100.0
qQ <- simpleQuote 0.02
rQ <- simpleQuote 0.05
volQ <- simpleQuote 0.10
qTS <- flatForward (ReferenceDate today') qQ dc Continuous Annual
rTS <- flatForward (ReferenceDate today') rQ dc Continuous Annual
tgt <- calendar TARGET
volTS <- blackConstantVol (CalendarReferenceDate today') tgt volQ dc
proc <- blackScholesMertonProcess spotQ qTS rTS volTS EulerDiscretization False
let payoff = PlainVanilla (PlainVanillaPayoff Call 100.0)
exercise = European (EuropeanExercise (addDays 360 today'))
opt <- barrierOption DownIn 90.0 0 payoff exercise >>= asOneAssetOption
eng <- mcBarrierEngine LowDiscrepancy Statistics proc Nothing (Just 1) True False (Just 131071) Nothing (Just 1048575) False 5
setPricingEngine opt eng
v <- npv opt
v `shouldSatisfy` closePrec expected (2.0e-2 * expected)
-- Golden values come from
-- digitaloption.cpp's testCashAtExpiryOrNothingAmericanValues and
-- testAssetAtExpiryOrNothingAmericanValues, on a payoff-at-expiry American exercise --
-- knockin picks AnalyticDigitalAmericanEngine vs. AnalyticDigitalAmericanKOEngine.
describe "AnalyticDigitalAmericanEngine / AnalyticDigitalAmericanKoEngine" $ do
let today' = 28 `august` 2026
priceCase mkPayoff (ty, strike, spot, q, r, tDays, vol, knockIn, _expected) =
Context.keepingSettingsGc $ do
Context.setEvaluationDate (Just today')
dc <- dayCounter (Actual360 False)
cal <- calendar Null
spotQ <- simpleQuote spot
qQ <- simpleQuote q
rQ <- simpleQuote r
volQ <- simpleQuote vol
qTS <- flatForward (SettlementDays 0 cal) qQ dc Continuous Annual
rTS <- flatForward (SettlementDays 0 cal) rQ dc Continuous Annual
volTS <- blackConstantVol (CalendarSettlementDays 0) cal volQ dc
process <- blackScholesMertonProcess spotQ qTS rTS volTS EulerDiscretization False
engine <- if knockIn then analyticDigitalAmericanEngine process
else analyticDigitalAmericanKoEngine process
let exDate = addDays tDays today'
opt <- europeanOption (mkPayoff ty strike) (American Nothing exDate True)
>>= asOneAssetOption >>= asOption >>= asInstrument
setPricingEngine opt engine
npv opt
cashOrNothingCases =
[ (Put, 100, 105, 0.00, 0.10, 180, 0.20, True, 9.3604 :: Double)
, (Call, 100, 95, 0.00, 0.10, 180, 0.20, True, 11.2223)
, (Put, 100, 105, 0.00, 0.10, 180, 0.20, False, 4.9081)
, (Call, 100, 95, 0.00, 0.10, 180, 0.20, False, 3.0461)
]
assetOrNothingCases =
[ (Put, 100, 105, 0.00, 0.10, 180, 0.20, True, 64.8426 :: Double)
, (Call, 100, 95, 0.00, 0.10, 180, 0.20, True, 77.7017)
, (Put, 100, 105, 0.00, 0.10, 180, 0.20, False, 40.1574)
, (Call, 100, 95, 0.00, 0.10, 180, 0.20, False, 17.2983)
]
it "cash-(at-expiry)-or-nothing American digital reproduces digitaloption.cpp" $
forM_ cashOrNothingCases $ \c@(_,_,_,_,_,_,_,_,expected) -> do
v <- priceCase (\ty s -> CashOrNothing ty s 15.0) c
v `shouldSatisfy` closePrec expected 1e-4
it "asset-(at-expiry)-or-nothing American digital reproduces digitaloption.cpp" $
forM_ assetOrNothingCases $ \c@(_,_,_,_,_,_,_,_,expected) -> do
v <- priceCase AssetOrNothing c
v `shouldSatisfy` closePrec expected 1e-4
-- cached references from digitaloption.cpp::testMCCashAtHit: cash-(at-hit)-or-nothing
-- American digital, priced via MakeMCDigitalEngine (default payoffAtExpiry=False, i.e. the
-- cash is paid at the moment the strike is hit, not at exercise).
it "mcDigitalEngine reproduces digitaloption.cpp's cash-at-hit values" $
Context.keepingSettingsGc $ do
let evalDate' = 1 `january` 2020
cases = [ (Put, 100.0, 105.0, 0.20, 0.10, 0.5 :: Double, 0.20, 12.2715 :: Double)
, (Call, 100.0, 95.0, 0.20, 0.10, 0.5, 0.20, 8.9109)
]
Context.setEvaluationDate (Just evalDate')
dc <- dayCounter (Actual360 False)
forM_ cases $ \(ty, strike, spot, q, r, t, vol, expected) -> do
spotQ <- simpleQuote spot
qQ <- simpleQuote q
rQ <- simpleQuote r
volQ <- simpleQuote vol
qTS <- flatForward (ReferenceDate evalDate') qQ dc Continuous Annual
rTS <- flatForward (ReferenceDate evalDate') rQ dc Continuous Annual
tgt <- calendar TARGET
volTS <- blackConstantVol (CalendarReferenceDate evalDate') tgt volQ dc
proc <- blackScholesMertonProcess spotQ qTS rTS volTS EulerDiscretization False
let exDate = addDays (round (t * 360 :: Double)) evalDate'
opt <- europeanOption (CashOrNothing ty strike 15.0) (American (Just evalDate') exDate False)
>>= asOneAssetOption >>= asOption >>= asInstrument
eng <- mcDigitalEngine LowDiscrepancy Statistics proc Nothing (Just 90) True False (Just 16383) Nothing (Just 1000000) 1
setPricingEngine opt eng
v <- npv opt
v `shouldSatisfy` closePrec expected 1.0e-2
-- Ported from test-suite/americanoption.cpp::testQdAmericanEngines' "standard put option" /
-- "call-put parity on standard option" cached edge cases -- a 10-year put and its
-- put-call-symmetric call (r and q swapped) share the same American value. QdPlus is only an
-- initial-guess-quality approximation (see its class haddock), and indeed doesn't converge to
-- the true American price here: an independent LeisenReimer(20001) binomial reprices this same
-- contract at ~23.0000, ~0.026 away from QdPlus's cached 22.9738 -- so qdPlusAmericanEngine is
-- checked against upstream's own cached (regression, not "true price") value, while
-- qdFpAmericanEngine -- the higher-precision fixed-point refinement -- is cross-checked against
-- that binomial reference instead, across all three schemes and all three FixedPointEquations.
-- testAndersenLakeHighPrecisionExample/testBulkQdFpAmericanEngine aren't ported directly: both
-- need QdFpLegendreScheme/QdFpGaussLobattoScheme, deliberately unbound (QdFpScheme only exposes
-- the three built-in static-factory schemes).
describe "QD+ / QD-FP American engines" $ do
let today' = 1 `june` 2022
qdPlusExpected = 22.97383256003585 :: Double
putCase = (Put, 100.0, 120.0, 0.10, 0.03)
cases = [ putCase
, (Call, 120.0, 100.0, 0.03, 0.10)
]
mkOption ty spot strike r q = do
spotQ <- simpleQuote spot
qQ <- simpleQuote q
rQ <- simpleQuote r
volQ <- simpleQuote 0.25
dc <- dayCounter Actual365FixedStandard
qTS <- flatForward (ReferenceDate today') qQ dc Continuous Annual
rTS <- flatForward (ReferenceDate today') rQ dc Continuous Annual
tgt <- calendar TARGET
volTS <- blackConstantVol (CalendarReferenceDate today') tgt volQ dc
process <- blackScholesMertonProcess spotQ qTS rTS volTS EulerDiscretization False
let exDate = addDays 3650 today'
opt <- vanillaOption (PlainVanilla (PlainVanillaPayoff ty strike)) (American Nothing exDate False)
pure (process, opt)
it "qdPlusAmericanEngine reproduces americanoption.cpp's standard put/call cached values" $
Context.keepingSettingsGc $ do
Context.setEvaluationDate (Just today')
forM_ cases $ \(ty, spot, strike, r, q) -> do
(process, opt) <- mkOption ty spot strike r q
eng <- qdPlusAmericanEngine process 8 Halley 1e-10 Nothing
setPricingEngine opt eng
v <- npv opt
v `shouldSatisfy` closePrec qdPlusExpected 1e-8
it "qdFpAmericanEngine agrees with a converged binomial price across every scheme/equation" $
Context.keepingSettingsGc $ do
Context.setEvaluationDate (Just today')
let (ty, spot, strike, r, q) = putCase
(biProcess, biOpt) <- mkOption ty spot strike r q
biEng <- binomialVanillaEngine LeisenReimer biProcess 20001
setPricingEngine biOpt biEng
binomialPrice <- npv biOpt
forM_ [FastScheme, AccurateScheme, HighPrecisionScheme] $ \scheme ->
forM_ [FP_A, FP_B, Auto] $ \fpEquation -> do
(process, opt) <- mkOption ty spot strike r q
eng <- qdFpAmericanEngine process scheme fpEquation
setPricingEngine opt eng
v <- npv opt
v `shouldSatisfy` closePrec binomialPrice 5e-3
-- Ported from test-suite/quantooption.cpp. Each case wires a QuantoEngine<Instr,Engine>
-- instantiation (a GeneralizedBlackScholesProcess plus a foreign risk-free curve, an
-- exchange-rate vol surface, and a correlation quote) around the matching base engine, and
-- checks the resulting quanto-adjusted NPV against Haug's cached literature values.
describe "Quanto engines" $ do
let today' = 1 `january` 2020
setupFlat dc tgt s q r vol fxr fxv corr = do
spotQ <- simpleQuote s
qQ <- simpleQuote q
rQ <- simpleQuote r
volQ <- simpleQuote vol
fxrQ <- simpleQuote fxr
fxvQ <- simpleQuote fxv
corrQ <- simpleQuote corr
qTS <- flatForward (ReferenceDate today') qQ dc Continuous Annual
rTS <- flatForward (ReferenceDate today') rQ dc Continuous Annual
volTS <- blackConstantVol (CalendarReferenceDate today') tgt volQ dc
fxrTS <- flatForward (ReferenceDate today') fxrQ dc Continuous Annual
fxVolTS <- blackConstantVol (CalendarReferenceDate today') tgt fxvQ dc
proc <- blackScholesMertonProcess spotQ qTS rTS volTS EulerDiscretization False
return (proc, fxrTS, fxVolTS, corrQ)
it "QuantoEngine<VanillaOption,AnalyticEuropeanEngine> reproduces testValues" $
Context.keepingSettingsGc $ do
Context.setEvaluationDate (Just today')
dc <- dayCounter (Actual360 False)
tgt <- calendar TARGET
let cases = [ (Call, 105.0, 100.0, 0.04, 0.08, 0.5 :: Double, 0.2, 0.05, 0.10, 0.3, 5.3280 / 1.5 :: Double)
, (Put, 105.0, 100.0, 0.04, 0.08, 0.5, 0.2, 0.05, 0.10, 0.3, 8.1636)
]
forM_ cases $ \(ty, strike, s, q, r, t, vol, fxr, fxv, corr, expected) -> do
(proc, fxrTS, fxVolTS, corrQ) <- setupFlat dc tgt s q r vol fxr fxv corr
engine <- quantoEuropeanEngine proc fxrTS fxVolTS corrQ
let exDate = addDays (round (t * 360 :: Double)) today'
opt <- quantoVanillaOption (PlainVanilla (PlainVanillaPayoff ty strike)) (European (EuropeanExercise exDate))
optInst <- asOneAssetOption opt
setPricingEngine optInst engine
v <- npv optInst
v `shouldSatisfy` closePrec expected 1.0e-4
-- Ported from quantooption.cpp's testGreeks: for every (type, strike) pair, bump each of
-- spot/rRate/qRate/vol/fxRate/fxVol/correlation (and the evaluation date, for theta) by a
-- relative 1e-4 and compare the resulting central-difference estimate against the engine's
-- own analytic greek, at every combination of the other quotes' levels.
-- Every greek but theta holds upstream's own 1e-5 tolerance. theta's central difference
-- bumps the evaluation date by a full 2 days -- a far larger relative step (~0.3% of the
-- 2-year maturity) than the other greeks' 1e-4-relative quote bumps -- so its own
-- discretization error dominates at the vol=1.20 (120% annual) extreme; measured worst case
-- across the whole sweep is ~3.1e-4 (Put, strike=150, vol=1.2, fxVol=1.2, corr=0.9), so 5e-4
-- keeps this a real check while accommodating it, per CLAUDE.md's numeric-tolerance rule.
it "QuantoEngine<VanillaOption,AnalyticEuropeanEngine> reproduces testGreeks" $
Context.keepingSettingsGc $ do
Context.setEvaluationDate (Just today')
dc <- dayCounter (Actual360 False)
tgt <- calendar TARGET
spotQ <- simpleQuote (0.0 :: Double)
qRateQ <- simpleQuote (0.0 :: Double)
rRateQ <- simpleQuote (0.0 :: Double)
volQ <- simpleQuote (0.0 :: Double)
fxRateQ <- simpleQuote (0.0 :: Double)
fxVolQ <- simpleQuote (0.0 :: Double)
corrQ <- simpleQuote (0.0 :: Double)
qTS <- flatForward (SettlementDays 0 tgt) qRateQ dc Continuous Annual
rTS <- flatForward (SettlementDays 0 tgt) rRateQ dc Continuous Annual
volTS <- blackConstantVol (CalendarSettlementDays 0) tgt volQ dc
fxrTS <- flatForward (SettlementDays 0 tgt) fxRateQ dc Continuous Annual
fxVolTS <- blackConstantVol (CalendarSettlementDays 0) tgt fxVolQ dc
proc <- blackScholesMertonProcess spotQ qTS rTS volTS EulerDiscretization False
engine <- quantoEuropeanEngine proc fxrTS fxVolTS corrQ
let u = 100.0 :: Double
qRates = [0.04, 0.05 :: Double]
rRates = [0.01, 0.05, 0.15 :: Double]
vols = [0.11, 1.20 :: Double]
corrs = [0.10, 0.90 :: Double]
innerCases = [ (q, r, v, fxr, fxv, corr)
| q <- qRates, r <- rRates, v <- vols
, fxr <- rRates, fxv <- vols, corr <- corrs
]
forM_ [Call, Put] $ \ty -> forM_ [50.0, 99.5, 100.0, 100.5, 150.0 :: Double] $ \strike -> do
let exDate = addGregorianYearsClip 2 today'
opt <- quantoVanillaOption (PlainVanilla (PlainVanillaPayoff ty strike)) (European (EuropeanExercise exDate))
optInst <- asOneAssetOption opt
setPricingEngine optInst engine
forM_ innerCases $ \(q, r, v, fxr, fxv, corr) -> do
_ <- setValue spotQ u
_ <- setValue qRateQ q
_ <- setValue rRateQ r
_ <- setValue volQ v
_ <- setValue fxRateQ fxr
_ <- setValue fxVolQ fxv
_ <- setValue corrQ corr
val <- npv optInst
when (val > u * 1.0e-5) $ do
calcDelta <- delta optInst
calcGamma <- gamma optInst
calcTheta <- theta optInst
calcRho <- Opt.rho optInst
calcDivRho <- Opt.dividendRho optInst
calcVega <- Opt.vega optInst
calcQrho <- qrho opt
calcQvega <- qvega opt
calcQlambda <- qlambda opt
let du = u * 1.0e-4
_ <- setValue spotQ (u + du); valueP <- npv optInst; deltaP <- delta optInst
_ <- setValue spotQ (u - du); valueM <- npv optInst; deltaM <- delta optInst
_ <- setValue spotQ u
let expDelta = (valueP - valueM) / (2 * du)
expGamma = (deltaP - deltaM) / (2 * du)
let dr = r * 1.0e-4
_ <- setValue rRateQ (r + dr); rValueP <- npv optInst
_ <- setValue rRateQ (r - dr); rValueM <- npv optInst
_ <- setValue rRateQ r
let expRho = (rValueP - rValueM) / (2 * dr)
let dq = q * 1.0e-4
_ <- setValue qRateQ (q + dq); qValueP <- npv optInst
_ <- setValue qRateQ (q - dq); qValueM <- npv optInst
_ <- setValue qRateQ q
let expDivRho = (qValueP - qValueM) / (2 * dq)
let dv = v * 1.0e-4
_ <- setValue volQ (v + dv); vValueP <- npv optInst
_ <- setValue volQ (v - dv); vValueM <- npv optInst
_ <- setValue volQ v
let expVega = (vValueP - vValueM) / (2 * dv)
let dfxr = fxr * 1.0e-4
_ <- setValue fxRateQ (fxr + dfxr); fxrValueP <- npv optInst
_ <- setValue fxRateQ (fxr - dfxr); fxrValueM <- npv optInst
_ <- setValue fxRateQ fxr
let expQrho = (fxrValueP - fxrValueM) / (2 * dfxr)
let dfxv = fxv * 1.0e-4
_ <- setValue fxVolQ (fxv + dfxv); fxvValueP <- npv optInst
_ <- setValue fxVolQ (fxv - dfxv); fxvValueM <- npv optInst
_ <- setValue fxVolQ fxv
let expQvega = (fxvValueP - fxvValueM) / (2 * dfxv)
let dcorr = corr * 1.0e-4
_ <- setValue corrQ (corr + dcorr); corrValueP <- npv optInst
_ <- setValue corrQ (corr - dcorr); corrValueM <- npv optInst
_ <- setValue corrQ corr
let expQlambda = (corrValueP - corrValueM) / (2 * dcorr)
dTyears <- yearFraction dc (addDays (-1) today') (addDays 1 today') Nothing Nothing
Context.setEvaluationDate (Just (addDays (-1) today'))
thetaValueM <- npv optInst
Context.setEvaluationDate (Just (addDays 1 today'))
thetaValueP <- npv optInst
Context.setEvaluationDate (Just today')
let expTheta = (thetaValueP - thetaValueM) / dTyears
relErr expctd calcd = abs (expctd - calcd) / u
relErr expDelta calcDelta `shouldSatisfy` (< 1.0e-5)
relErr expGamma calcGamma `shouldSatisfy` (< 1.0e-5)
relErr expTheta calcTheta `shouldSatisfy` (< 5.0e-4)
relErr expRho calcRho `shouldSatisfy` (< 1.0e-5)
relErr expDivRho calcDivRho `shouldSatisfy` (< 1.0e-5)
relErr expVega calcVega `shouldSatisfy` (< 1.0e-5)
relErr expQrho calcQrho `shouldSatisfy` (< 1.0e-5)
relErr expQvega calcQvega `shouldSatisfy` (< 1.0e-5)
relErr expQlambda calcQlambda `shouldSatisfy` (< 1.0e-5)
it "QuantoEngine<ForwardVanillaOption,ForwardVanillaEngine<AnalyticEuropeanEngine>> reproduces testForwardValues" $
Context.keepingSettingsGc $ do
Context.setEvaluationDate (Just today')
dc <- dayCounter (Actual360 False)
tgt <- calendar TARGET
let cases = [ (Call, 1.05 :: Double, 100.0, 0.04, 0.08, 0.00 :: Double, 0.5 :: Double, 0.20, 0.05, 0.10, 0.3, 5.3280 / 1.5 :: Double)
, (Put, 1.05, 100.0, 0.04, 0.08, 0.00, 0.5, 0.20, 0.05, 0.10, 0.3, 8.1636)
, (Call, 1.05, 100.0, 0.04, 0.08, 0.25, 0.5, 0.20, 0.05, 0.10, 0.3, 2.0171)
, (Put, 1.05, 100.0, 0.04, 0.08, 0.25, 0.5, 0.20, 0.05, 0.10, 0.3, 6.7296)
]
forM_ cases $ \(ty, moneyness, s, q, r, start, t, vol, fxr, fxv, corr, expected) -> do
(proc, fxrTS, fxVolTS, corrQ) <- setupFlat dc tgt s q r vol fxr fxv corr
engine <- quantoForwardEuropeanEngine proc fxrTS fxVolTS corrQ
let exDate = addDays (round (t * 360 :: Double)) today'
resetDate = addDays (round (start * 360 :: Double)) today'
opt <- quantoForwardVanillaOption moneyness resetDate (PlainVanilla (PlainVanillaPayoff ty 0.0)) (European (EuropeanExercise exDate))
optInst <- asOneAssetOption opt
setPricingEngine optInst engine
v <- npv optInst
v `shouldSatisfy` closePrec expected 1.0e-4
-- Ported from quantooption.cpp's testForwardGreeks, same bump-and-revalue shape as testGreeks
-- above but over 'QuantoForwardVanillaOption' (type, moneyness, resetMonths) combinations.
-- theta's tolerance is loosened further than testGreeks' (see the comment there): the
-- forward-starting payoff's moneyness-relative strike makes its date sensitivity more
-- nonlinear, and the measured worst case across the sweep is ~2.5e-3 (Put, moneyness=1.1,
-- reset=6m, vol=1.2, fxVol=1.2, corr=0.9); every other greek still holds 1e-5.
it "QuantoEngine<ForwardVanillaOption,ForwardVanillaEngine<AnalyticEuropeanEngine>> reproduces testForwardGreeks" $
Context.keepingSettingsGc $ do
Context.setEvaluationDate (Just today')
dc <- dayCounter (Actual360 False)
tgt <- calendar TARGET
spotQ <- simpleQuote (0.0 :: Double)
qRateQ <- simpleQuote (0.0 :: Double)
rRateQ <- simpleQuote (0.0 :: Double)
volQ <- simpleQuote (0.0 :: Double)
fxRateQ <- simpleQuote (0.0 :: Double)
fxVolQ <- simpleQuote (0.0 :: Double)
corrQ <- simpleQuote (0.0 :: Double)
qTS <- flatForward (SettlementDays 0 tgt) qRateQ dc Continuous Annual
rTS <- flatForward (SettlementDays 0 tgt) rRateQ dc Continuous Annual
volTS <- blackConstantVol (CalendarSettlementDays 0) tgt volQ dc
fxrTS <- flatForward (SettlementDays 0 tgt) fxRateQ dc Continuous Annual
fxVolTS <- blackConstantVol (CalendarSettlementDays 0) tgt fxVolQ dc
proc <- blackScholesMertonProcess spotQ qTS rTS volTS EulerDiscretization False
engine <- quantoForwardEuropeanEngine proc fxrTS fxVolTS corrQ
let u = 100.0 :: Double
qRates = [0.04, 0.05 :: Double]
rRates = [0.01, 0.05, 0.15 :: Double]
vols = [0.11, 1.20 :: Double]
corrs = [0.10, 0.90 :: Double]
innerCases = [ (q, r, v, fxr, fxv, corr)
| q <- qRates, r <- rRates, v <- vols
, fxr <- rRates, fxv <- vols, corr <- corrs
]
forM_ [Call, Put] $ \ty -> forM_ [0.9, 1.0, 1.1 :: Double] $ \moneyness ->
forM_ [6, 9 :: Integer] $ \startMonth -> do
let exDate = addGregorianYearsClip 2 today'
resetDate <- addPeriod today' (fromInteger startMonth, Months)
opt <- quantoForwardVanillaOption moneyness resetDate (PlainVanilla (PlainVanillaPayoff ty 0.0)) (European (EuropeanExercise exDate))
optInst <- asOneAssetOption opt
setPricingEngine optInst engine
forM_ innerCases $ \(q, r, v, fxr, fxv, corr) -> do
_ <- setValue spotQ u
_ <- setValue qRateQ q
_ <- setValue rRateQ r
_ <- setValue volQ v
_ <- setValue fxRateQ fxr
_ <- setValue fxVolQ fxv
_ <- setValue corrQ corr
val <- npv optInst
when (val > u * 1.0e-5) $ do
calcDelta <- delta optInst
calcGamma <- gamma optInst
calcTheta <- theta optInst
calcRho <- Opt.rho optInst
calcDivRho <- Opt.dividendRho optInst
calcVega <- Opt.vega optInst
calcQrho <- qrho opt
calcQvega <- qvega opt
calcQlambda <- qlambda opt
let du = u * 1.0e-4
_ <- setValue spotQ (u + du); valueP <- npv optInst; deltaP <- delta optInst
_ <- setValue spotQ (u - du); valueM <- npv optInst; deltaM <- delta optInst
_ <- setValue spotQ u
let expDelta = (valueP - valueM) / (2 * du)
expGamma = (deltaP - deltaM) / (2 * du)
let dr = r * 1.0e-4
_ <- setValue rRateQ (r + dr); rValueP <- npv optInst
_ <- setValue rRateQ (r - dr); rValueM <- npv optInst
_ <- setValue rRateQ r
let expRho = (rValueP - rValueM) / (2 * dr)
let dq = q * 1.0e-4
_ <- setValue qRateQ (q + dq); qValueP <- npv optInst
_ <- setValue qRateQ (q - dq); qValueM <- npv optInst
_ <- setValue qRateQ q
let expDivRho = (qValueP - qValueM) / (2 * dq)
let dv = v * 1.0e-4
_ <- setValue volQ (v + dv); vValueP <- npv optInst
_ <- setValue volQ (v - dv); vValueM <- npv optInst
_ <- setValue volQ v
let expVega = (vValueP - vValueM) / (2 * dv)
let dfxr = fxr * 1.0e-4
_ <- setValue fxRateQ (fxr + dfxr); fxrValueP <- npv optInst
_ <- setValue fxRateQ (fxr - dfxr); fxrValueM <- npv optInst
_ <- setValue fxRateQ fxr
let expQrho = (fxrValueP - fxrValueM) / (2 * dfxr)
let dfxv = fxv * 1.0e-4
_ <- setValue fxVolQ (fxv + dfxv); fxvValueP <- npv optInst
_ <- setValue fxVolQ (fxv - dfxv); fxvValueM <- npv optInst
_ <- setValue fxVolQ fxv
let expQvega = (fxvValueP - fxvValueM) / (2 * dfxv)
let dcorr = corr * 1.0e-4
_ <- setValue corrQ (corr + dcorr); corrValueP <- npv optInst
_ <- setValue corrQ (corr - dcorr); corrValueM <- npv optInst
_ <- setValue corrQ corr
let expQlambda = (corrValueP - corrValueM) / (2 * dcorr)
dTyears <- yearFraction dc (addDays (-1) today') (addDays 1 today') Nothing Nothing
Context.setEvaluationDate (Just (addDays (-1) today'))
thetaValueM <- npv optInst
Context.setEvaluationDate (Just (addDays 1 today'))
thetaValueP <- npv optInst
Context.setEvaluationDate (Just today')
let expTheta = (thetaValueP - thetaValueM) / dTyears
relErr expctd calcd = abs (expctd - calcd) / u
relErr expDelta calcDelta `shouldSatisfy` (< 1.0e-5)
relErr expGamma calcGamma `shouldSatisfy` (< 1.0e-5)
relErr expTheta calcTheta `shouldSatisfy` (< 3.0e-3)
relErr expRho calcRho `shouldSatisfy` (< 1.0e-5)
relErr expDivRho calcDivRho `shouldSatisfy` (< 1.0e-5)
relErr expVega calcVega `shouldSatisfy` (< 1.0e-5)
relErr expQrho calcQrho `shouldSatisfy` (< 1.0e-5)
relErr expQvega calcQvega `shouldSatisfy` (< 1.0e-5)
relErr expQlambda calcQlambda `shouldSatisfy` (< 1.0e-5)
it "QuantoEngine<ForwardVanillaOption,ForwardPerformanceVanillaEngine<AnalyticEuropeanEngine>> reproduces testForwardPerformanceValues" $
Context.keepingSettingsGc $ do
Context.setEvaluationDate (Just today')
dc <- dayCounter (Actual360 False)
tgt <- calendar TARGET
let cases = [ (Call, 1.05 :: Double, 100.0, 0.04, 0.08, 0.00 :: Double, 0.5 :: Double, 0.20, 0.05, 0.10, 0.3, 5.3280 / 150 :: Double)
, (Put, 1.05, 100.0, 0.04, 0.08, 0.00, 0.5, 0.20, 0.05, 0.10, 0.3, 0.0816)
, (Call, 1.05, 100.0, 0.04, 0.08, 0.25, 0.5, 0.20, 0.05, 0.10, 0.3, 0.0201)
, (Put, 1.05, 100.0, 0.04, 0.08, 0.25, 0.5, 0.20, 0.05, 0.10, 0.3, 0.0672)
]
forM_ cases $ \(ty, moneyness, s, q, r, start, t, vol, fxr, fxv, corr, expected) -> do
(proc, fxrTS, fxVolTS, corrQ) <- setupFlat dc tgt s q r vol fxr fxv corr
engine <- quantoForwardPerformanceEuropeanEngine proc fxrTS fxVolTS corrQ
let exDate = addDays (round (t * 360 :: Double)) today'
resetDate = addDays (round (start * 360 :: Double)) today'
opt <- quantoForwardVanillaOption moneyness resetDate (PlainVanilla (PlainVanillaPayoff ty 0.0)) (European (EuropeanExercise exDate))
optInst <- asOneAssetOption opt
setPricingEngine optInst engine
v <- npv optInst
v `shouldSatisfy` closePrec expected 1.0e-4
it "QuantoEngine<BarrierOption,AnalyticBarrierEngine> reproduces testBarrierValues" $
Context.keepingSettingsGc $ do
Context.setEvaluationDate (Just today')
dc <- dayCounter (Actual360 False)
tgt <- calendar TARGET
let cases = [ (DownOut, 95.0 :: Double, 3.0 :: Double, Call, 100.0 :: Double, 90.0 :: Double, 0.04, 0.0212, 0.50 :: Double, 0.25, 0.05, 0.2, 0.3, 8.247 :: Double, 0.5 :: Double)
, (DownOut, 95.0, 3.0, Put, 100.0, 90.0, 0.04, 0.0212, 0.50, 0.25, 0.05, 0.2, 0.3, 2.274, 0.5)
, (DownIn, 95.0, 0.0, Put, 100.0, 90.0, 0.04, 0.0212, 0.50, 0.25, 0.05, 0.2, 0.3, 2.85, 0.5)
]
forM_ cases $ \(barType, barrier, rebate, ty, s, strike, q, r, t, vol, fxr, fxv, corr, expected, tol) -> do
(proc, fxrTS, fxVolTS, corrQ) <- setupFlat dc tgt s q r vol fxr fxv corr
engine <- quantoBarrierEngine proc fxrTS fxVolTS corrQ
let exDate = addDays (round (t * 360 :: Double)) today'
opt <- quantoBarrierOption barType barrier rebate (PlainVanilla (PlainVanillaPayoff ty strike)) (European (EuropeanExercise exDate))
optInst <- asOneAssetOption opt
setPricingEngine optInst engine
v <- npv optInst
v `shouldSatisfy` closePrec expected tol
it "QuantoEngine<DoubleBarrierOption,AnalyticDoubleBarrierEngine> reproduces testDoubleBarrierValues" $
Context.keepingSettingsGc $ do
Context.setEvaluationDate (Just today')
dc <- dayCounter (Actual360 False)
tgt <- calendar TARGET
let cases = [ (KnockOut, 50.0 :: Double, 150.0 :: Double, 0.0 :: Double, Call, 100.0 :: Double, 100.0 :: Double, 0.00 :: Double, 0.1 :: Double, 0.25 :: Double, 0.15, 0.05, 0.2, 0.3, 3.4623 :: Double)
, (KnockOut, 90.0, 110.0, 0.0, Call, 100.0, 100.0, 0.00, 0.1, 0.50, 0.15, 0.05, 0.2, 0.3, 0.5236)
, (KnockOut, 90.0, 110.0, 0.0, Put, 100.0, 100.0, 0.00, 0.1, 0.25, 0.15, 0.05, 0.2, 0.3, 1.1320)
, (KnockIn, 80.0, 120.0, 0.0, Call, 100.0, 102.0, 0.00, 0.1, 0.25, 0.25, 0.05, 0.2, 0.3, 2.6313)
, (KnockIn, 80.0, 120.0, 0.0, Call, 100.0, 102.0, 0.00, 0.1, 0.50, 0.15, 0.05, 0.2, 0.3, 1.9305)
]
forM_ cases $ \(barType, barLo, barHi, rebate, ty, s, strike, q, r, t, vol, fxr, fxv, corr, expected) -> do
(proc, fxrTS, fxVolTS, corrQ) <- setupFlat dc tgt s q r vol fxr fxv corr
engine <- quantoDoubleBarrierEngine proc fxrTS fxVolTS corrQ
let exDate = addDays (round (t * 360 :: Double)) today'
opt <- quantoDoubleBarrierOption barType barLo barHi rebate (PlainVanilla (PlainVanillaPayoff ty strike)) (European (EuropeanExercise exDate))
optInst <- asOneAssetOption opt
setPricingEngine optInst engine
v <- npv optInst
v `shouldSatisfy` closePrec expected 1.0e-4
-- Ported from quantooption.cpp's testFDMQuantoHelper, testPDEOptionValues, and
-- testAmericanQuantoOption: the FDM-side building blocks (FdmQuantoHelper's own quanto drift
-- adjustment, and the FdmBlackScholesMesher grid it feeds into) and the FD-vs-analytic /
-- FD-vs-FD cross-checks that exercise the empty and non-empty dividend dispatch paths of
-- fdBlackScholesVanillaEngineQuanto and fdHestonVanillaEngineQuanto.
describe "FdmQuantoHelper / FD quanto engines" $ do
it "FdmQuantoHelper.quantoAdjustment and FdmBlackScholesMesher grid bounds reproduce testFDMQuantoHelper" $
Context.keepingSettingsGc $ do
let today' = 22 `april` 2019
s = 100.0 :: Double
domesticR = 0.1 :: Double
foreignR = 0.2 :: Double
q = 0.3 :: Double
vol = 0.3 :: Double
fxVol = 0.2 :: Double
exchRateATMlevel = 1.0 :: Double
equityFxCorrelation = -0.75 :: Double
Context.setEvaluationDate (Just today')
dc <- dayCounter (Actual360 False)
tgt <- calendar TARGET
domesticRQ <- simpleQuote domesticR
domesticTS <- flatForward (ReferenceDate today') domesticRQ dc Continuous Annual
divQ <- simpleQuote q
divTS <- flatForward (ReferenceDate today') divQ dc Continuous Annual
volQ <- simpleQuote vol
volTS <- blackConstantVol (CalendarReferenceDate today') tgt volQ dc
spotQ <- simpleQuote s
bsmProcess <- blackScholesMertonProcess spotQ divTS domesticTS volTS EulerDiscretization False
foreignRQ <- simpleQuote foreignR
foreignTS <- flatForward (ReferenceDate today') foreignRQ dc Continuous Annual
fxVolQ <- simpleQuote fxVol
fxVolTS <- blackConstantVol (CalendarReferenceDate today') tgt fxVolQ dc
fdmHelper <- fdmQuantoHelper domesticTS foreignTS fxVolTS equityFxCorrelation exchRateATMlevel
calculatedQuantoAdj <- quantoAdjustment fdmHelper vol 0.0 1.0
let expectedQuantoAdj = domesticR - foreignR + equityFxCorrelation * vol * fxVol
calculatedQuantoAdj `shouldSatisfy` closePrec expectedQuantoAdj 1.0e-10
maturityDate <- addPeriod today' (6, Months)
maturityTime <- yearFraction dc today' maturityDate Nothing Nothing
let eps = 0.0002 :: Double
scalingFactor = 1.25 :: Double
mesher1d <- fdmBlackScholesMesher 3 bsmProcess maturityTime s Nothing Nothing eps scalingFactor Nothing Nothing [] (Just fdmHelper) 0.0
fdmMesher <- fdmMesherComposite [mesher1d]
loc <- fdmMesherLocations fdmMesher 0
let loc0 = loc V.! 0
-- InverseCumulativeNormal()(1 - eps), eps = 0.0002 -- a fixed literal (like the other
-- ported golden values) rather than re-deriving boost's inverse normal CDF here.
let normInvEps = 3.5400837992061738 :: Double
sigmaSqrtT = vol * sqrt maturityTime
qQuanto = q + expectedQuantoAdj
expectedDriftRate = domesticR - qQuanto
logFwd = log s + expectedDriftRate * maturityTime
xMin = logFwd - sigmaSqrtT * normInvEps * scalingFactor
xMax = log s + sigmaSqrtT * normInvEps * scalingFactor
loc0 `shouldSatisfy` closePrec xMin 1.0e-6
V.last loc `shouldSatisfy` closePrec xMax 1.0e-6
it "fdBlackScholesVanillaEngineQuanto reproduces QuantoEngine<VanillaOption,AnalyticEuropeanEngine> (testPDEOptionValues)" $ do
let today' = 21 `april` 2019
cases = [ (Call, 105.0 :: Double, 100.0 :: Double, 0.04 :: Double, 0.08 :: Double, 0.5 :: Double, 0.2 :: Double, 0.05 :: Double, 0.10 :: Double, 0.3 :: Double)
, (Call, 100.0, 100.0, 0.16, 0.08, 0.25, 0.15, 0.05, 0.20, -0.3)
, (Call, 105.0, 100.0, 0.04, 0.08, 0.5, 0.2, 0.05, 0.10, 0.3)
, (Put, 105.0, 100.0, 0.04, 0.08, 0.5, 0.2, 0.05, 0.10, 0.3)
, (Call, 0.0, 100.0, 0.04, 0.08, 0.3, 0.3, 0.05, 0.10, 0.75)
]
forM_ cases $ \(ty, strike, s, q, r, t, vol, fxr, fxv, corr) ->
Context.keepingSettingsGc $ do
Context.setEvaluationDate (Just today')
dc <- dayCounter (Actual360 False)
tgt <- calendar TARGET
spotQ <- simpleQuote s
rQ <- simpleQuote r
domesticTS <- flatForward (ReferenceDate today') rQ dc Continuous Annual
divQ <- simpleQuote q
divTS <- flatForward (ReferenceDate today') divQ dc Continuous Annual
volQ <- simpleQuote vol
volTS <- blackConstantVol (CalendarReferenceDate today') tgt volQ dc
bsmProcess <- blackScholesMertonProcess spotQ divTS domesticTS volTS EulerDiscretization False
foreignRQ <- simpleQuote fxr
foreignTS <- flatForward (ReferenceDate today') foreignRQ dc Continuous Annual
fxVolQ <- simpleQuote fxv
fxVolTS <- blackConstantVol (CalendarReferenceDate today') tgt fxVolQ dc
quantoHelper <- fdmQuantoHelper domesticTS foreignTS fxVolTS corr 1.0
let exDate = addDays (round (t * 360 :: Double)) today'
opt <- vanillaOption (PlainVanilla (PlainVanillaPayoff ty strike)) (European (EuropeanExercise exDate))
pdeEngine <- fdBlackScholesVanillaEngineQuanto bsmProcess [] (Just quantoHelper) (round (t * 200 :: Double)) 500 1
Douglas False 0.0 CashDividendSpot
setPricingEngine opt pdeEngine
optInst <- asOneAssetOption opt
calcNpv <- npv optInst
calcDelta <- delta optInst
corrQ <- simpleQuote corr
analyticEngine <- quantoEuropeanEngine bsmProcess foreignTS fxVolTS corrQ
setPricingEngine opt analyticEngine
expNpv <- npv optInst
expDelta <- delta optInst
closePrec expNpv 2.0e-4 calcNpv `shouldBe` True
closePrec expDelta 1.0e-4 calcDelta `shouldBe` True
it "quanto FD engines with dividends reproduce testAmericanQuantoOption" $
Context.keepingSettingsGc $ do
let today' = 21 `april` 2019
domesticR = 0.025 :: Double
foreignR = 0.075 :: Double
q = 0.03 :: Double
vol = 0.3 :: Double
fxVol = 0.15 :: Double
equityFxCorrelation = -0.75 :: Double
strike = 105.0 :: Double
expected = 8.90611734 :: Double
tol = 1.0e-4 :: Double
Context.setEvaluationDate (Just today')
dc <- dayCounter Actual365FixedStandard
maturity <- addPeriod today' (9, Months)
domesticRQ <- simpleQuote domesticR
domesticTS <- flatForward (ReferenceDate today') domesticRQ dc Continuous Annual
divQ <- simpleQuote q
divTS <- flatForward (ReferenceDate today') divQ dc Continuous Annual
volQ <- simpleQuote vol
tgt <- calendar TARGET
volTS <- blackConstantVol (CalendarReferenceDate today') tgt volQ dc
spotQ <- simpleQuote (100.0 :: Double)
bsmProcess <- blackScholesMertonProcess spotQ divTS domesticTS volTS EulerDiscretization False
foreignRQ <- simpleQuote foreignR
foreignTS <- flatForward (ReferenceDate today') foreignRQ dc Continuous Annual
fxVolQ <- simpleQuote fxVol
fxVolTS <- blackConstantVol (CalendarReferenceDate today') tgt fxVolQ dc
quantoHelper <- fdmQuantoHelper domesticTS foreignTS fxVolTS equityFxCorrelation 1.0
divDate <- addPeriod today' (6, Months)
dividends <- sequence [fixedDividend 8.0 divDate]
opt <- vanillaOption (PlainVanilla (PlainVanillaPayoff Call strike)) (American Nothing maturity False)
optInst <- asOneAssetOption opt
bsEngine <- fdBlackScholesVanillaEngineQuanto bsmProcess dividends (Just quantoHelper) 100 400 1
Douglas False 0.0 CashDividendSpot
setPricingEngine opt bsEngine
bsCalculated <- npv optInst
closePrec expected tol bsCalculated `shouldBe` True
localVolEngine <- fdBlackScholesVanillaEngineQuanto bsmProcess dividends (Just quantoHelper) 100 400 1
Douglas False 0.0 CashDividendSpot
setPricingEngine opt localVolEngine
localVolCalculated <- npv optInst
closePrec expected tol localVolCalculated `shouldBe` True
closePrec bsCalculated 1.0e-6 localVolCalculated `shouldBe` True
divOpt <- vanillaOption (PlainVanilla (PlainVanillaPayoff Call strike)) (American Nothing maturity False)
divOptInst <- asOneAssetOption divOpt
let v0 = vol * vol
kappa = 1.0 :: Double
theta0 = v0
sigma = 1.0e-4 :: Double
hestonRho = 0.0 :: Double
hp <- hestonProcess domesticTS (Just divTS) spotQ v0 kappa theta0 sigma hestonRho QuadraticExponentialMartingale
hm <- hestonModel hp
hestonEngine <- fdHestonVanillaEngineQuanto hm dividends (Just quantoHelper) 100 400 3 1 Hundsdorfer Nothing 1.0
setPricingEngine divOpt hestonEngine
hestonCalculated <- npv divOptInst
closePrec expected tol hestonCalculated `shouldBe` True
constVolQ <- simpleQuote (2.0 :: Double)
localConstVol <- localConstantVol (ReferenceDate today') constVolQ dc
hp05 <- hestonProcess domesticTS (Just divTS) spotQ (0.25 * v0) kappa (0.25 * theta0) sigma hestonRho QuadraticExponentialMartingale
hm05 <- hestonModel hp05
hestonSlvEngine <- fdHestonVanillaEngineQuanto hm05 dividends (Just quantoHelper) 100 400 3 1 Hundsdorfer (Just localConstVol) 1.0
setPricingEngine divOpt hestonSlvEngine
hestonSlvCalculated <- npv divOptInst
closePrec expected tol hestonSlvCalculated `shouldBe` True
-- Ported from test-suite/twoassetbarrieroption.cpp's testHaugValues: a barrier option on two
-- correlated assets, where the first asset's value is compared to the strike and the second's
-- is monitored against the barrier (Heynen and Kat's formulas via AnalyticTwoAssetBarrierEngine).
describe "Two-asset barrier engine" $
it "AnalyticTwoAssetBarrierEngine reproduces twoassetbarrieroption.cpp's testHaugValues" $
Context.keepingSettingsGc $ do
let today' = 1 `january` 2020
cases = [ (DownOut, Call, 95.0 :: Double, 90.0 :: Double, 0.5 :: Double, 0.08 :: Double, 6.6592 :: Double)
, (UpOut, Call, 105.0, 90.0, -0.5, 0.08, 4.6670)
, (DownOut, Put, 95.0, 90.0, -0.5, 0.08, 0.6184)
, (UpOut, Put, 105.0, 100.0, 0.0, 0.08, 0.8246)
]
Context.setEvaluationDate (Just today')
dc <- dayCounter (Actual360 False)
tgt <- calendar TARGET
rQ <- simpleQuote (0.0 :: Double)
rTS <- flatForward (ReferenceDate today') rQ dc Continuous Annual
forM_ cases $ \(barType, ty, barrier, strike, corr, r, expected) -> do
_ <- setValue rQ r
s1Q <- simpleQuote (100.0 :: Double)
q1Q <- simpleQuote (0.0 :: Double)
v1Q <- simpleQuote (0.2 :: Double)
s2Q <- simpleQuote (100.0 :: Double)
q2Q <- simpleQuote (0.0 :: Double)
v2Q <- simpleQuote (0.2 :: Double)
rhoQ <- simpleQuote corr
q1TS <- flatForward (ReferenceDate today') q1Q dc Continuous Annual
q2TS <- flatForward (ReferenceDate today') q2Q dc Continuous Annual
vol1TS <- blackConstantVol (CalendarReferenceDate today') tgt v1Q dc
vol2TS <- blackConstantVol (CalendarReferenceDate today') tgt v2Q dc
proc1 <- blackScholesMertonProcess s1Q q1TS rTS vol1TS EulerDiscretization False
proc2 <- blackScholesMertonProcess s2Q q2TS rTS vol2TS EulerDiscretization False
engine <- analyticTwoAssetBarrierEngine proc1 proc2 rhoQ
let exDate = addDays 180 today'
inst <- twoAssetBarrierOption barType barrier (PlainVanilla (PlainVanillaPayoff ty strike)) (European (EuropeanExercise exDate))
setPricingEngine inst engine
v <- npv inst
v `shouldSatisfy` closePrec expected 4.0e-3
-- The Gaussian1dModels example already covers the model calibration and swaption paths.
-- These are its distinct cap/floor-engine scenarios, including the caplet-smile constructor
-- used by test-suite/markovfunctional.cpp's testVanillaEngines.
describe "Gaussian1d cap/floor engine" $
it "prices a cap for GSR and MarkovFunctional, and the caplet-calibrated Markov model agrees with Black" $
Context.keepingSettingsGc $ do
cal <- calendar TARGET
originalEvalDate <- Context.evaluationDate
evalDate <- adjust cal originalEvalDate Following
Context.setEvaluationDate (Just evalDate)
settlement <- advance cal evalDate (2, Days) Following False
dc365 <- dayCounter Actual365FixedStandard
thirty360bb <- dayCounter Thirty360BondBasis
act360 <- dayCounter (Actual360 False)
flatQ <- simpleQuote 0.03
ts <- flatForward (ReferenceDate settlement) flatQ dc365 Continuous Annual
euribor6m <- IR.iborIndex IR.Euribor6M (Just ts)
start <- advance cal settlement (1, Years) ModifiedFollowing False
maturity <- advance cal start (10, Years) ModifiedFollowing False
fixedSchedule <- schedule (Just start) maturity (1, Years) cal ModifiedFollowing ModifiedFollowing Forward False Nothing Nothing
floatSchedule <- schedule (Just start) maturity (6, Months) cal ModifiedFollowing ModifiedFollowing Forward False Nothing Nothing
swp <- vanillaSwap Payer 1.0 fixedSchedule 0.03 thirty360bb floatSchedule euribor6m 0.0 act360 (Just ModifiedFollowing) Nothing
floatLeg <- floatingLeg swp
capfl <- cap floatLeg (fromList [0.03])
stepDates <- mapM (\n -> advance cal evalDate (n, Years) Following False) [1, 2 :: Int]
gsrInitialVolQuote <- simpleQuote 0.01
gsrStepVolQuotes <- mapM simpleQuote [0.01, 0.01]
gsrReversionQuote <- simpleQuote 0.01
gsrModel <- gsr ts gsrInitialVolQuote (zip stepDates gsrStepVolQuotes) gsrReversionQuote 60.0
gsrVols <- volatilities gsrModel
gsrVols `shouldSatisfy` (\xs -> length xs == 3 && all (closePrec 0.01 1.0e-12) xs)
gsrModel' <- asGaussian1dModel gsrModel
gsrEngine <- gaussian1dCapFloorEngine gsrModel' 64 7.0 True False (Just ts)
setPricingEngine capfl gsrEngine
gsrNpv <- npv capfl
gsrNpv `shouldSatisfy` (\x -> not (isNaN x || isInfinite x) && x >= 0)
swapBase <- IR.liborSwapIndex IR.EuriborSwapIsdaFixA (10, Years) (Just ts) (Just ts)
swaptionVolQ <- simpleQuote 0.20
swaptionVol <- constantSwaptionVolatility (CalendarReferenceDate evalDate) cal ModifiedFollowing swaptionVolQ dc365 ShiftedLognormal 0.0
cmsExpiries <- mapM (\n -> advance cal evalDate (n, Years) Following False) [1, 2, 3 :: Int]
markov <- markovFunctional ts 0.01 0.01 [] swaptionVol (fromList $ zip cmsExpiries $ replicate 3 (10, Years)) swapBase 16
markovVols <- volatilities markov
markovVols `shouldSatisfy` (\xs -> not (null xs) && all (\x -> not (isNaN x || isInfinite x)) xs)
markovModel <- asGaussian1dModel markov
markovEngine <- gaussian1dCapFloorEngine markovModel 8 5.0 True False (Just ts)
setPricingEngine capfl markovEngine
markovNpv <- npv capfl
markovNpv `shouldSatisfy` (\x -> not (isNaN x || isInfinite x) && x >= 0)
capletExpiries <- CF.toCouponLeg floatLeg >>= CF.couponAccrualStartDates
capletVolQ <- simpleQuote 0.20
capletVol <- constantOptionletVolatility (CalendarSettlementDays 0) cal ModifiedFollowing capletVolQ dc365 ShiftedLognormal 0.0
markovCaplet <- markovFunctionalCaplet ts 0.01 0.01 [] capletVol (fromList capletExpiries) euribor6m 16
markovCapletModel <- asGaussian1dModel markovCaplet
blackEngine <- blackCapFloorEngineFromVolatilityStructure ts capletVol
setPricingEngine capfl blackEngine
blackNpv <- npv capfl
markovCapletEngine <- gaussian1dCapFloorEngine markovCapletModel 64 7.0 True False (Just ts)
setPricingEngine capfl markovCapletEngine
markovCapletNpv <- npv capfl
markovCapletNpv `shouldSatisfy` closePrec blackNpv 1.0e-4
-- Basket/spread pricing engines added alongside test-suite/basketoption.cpp's already-bound
-- core (BasketOption, BasketPayoff, StulzEngine, KirkEngine, MCEuropeanBasketEngine,
-- MCAmericanBasketEngine). One 'it' per upstream test function, per this file's usual
-- convention (see the quanto-engine blocks above) -- ported from
-- test/smoke/CheckBasketSpreadEngines.hs, which stays as the standalone smoke version.
describe "Basket and spread pricing engines" $ do
it "testEuroTwoValues: StulzEngine/KirkEngine vs. Fd2dBlackScholesVanillaEngine/MCEuropeanBasketEngine on a representative row subset" $
Context.keepingSettingsGc $ do
let evalDate = 1 `march` 2024
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter (Actual360 False)
cal <- calendar TARGET
-- basketType: 0=Min, 1=Max, 2=Spread
let rows :: [(Int, OptionType, Double, Double, Double, Double, Double, Double, Double, Double, Double, Double, Double, Double)]
rows =
-- basketType, type, strike, s1, s2, q1, q2, r, t, v1, v2, corr, result, tol
[ (0, Call, 100.0, 100.0, 100.0, 0.00, 0.00, 0.05, 1.00, 0.30, 0.30, 0.90, 10.898, 1.0e-3)
, (0, Call, 100.0, 100.0, 100.0, 0.00, 0.00, 0.05, 1.00, 0.30, 0.30, 0.10, 4.413, 1.0e-3)
, (1, Call, 100.0, 100.0, 100.0, 0.00, 0.00, 0.05, 1.00, 0.30, 0.30, 0.90, 17.565, 1.0e-3)
, (1, Call, 100.0, 80.0, 120.0, 0.00, 0.00, 0.05, 1.00, 0.30, 0.30, 0.30, 30.141, 1.0e-3)
, (0, Put, 100.0, 100.0, 100.0, 0.00, 0.00, 0.05, 1.00, 0.30, 0.30, 0.50, 13.890, 1.0e-3)
, (1, Put, 100.0, 100.0, 100.0, 0.00, 0.00, 0.05, 1.00, 0.30, 0.30, 0.30, 3.967, 1.1e-3)
, (0, Call, 98.0, 100.0, 105.0, 0.06, 0.09, 0.05, 0.50, 0.11, 0.16, 0.63, 2.9340, 1.0e-4)
, (1, Call, 98.0, 100.0, 105.0, 0.06, 0.09, 0.05, 0.50, 0.11, 0.16, 0.63, 8.0701, 1.0e-4)
, (0, Put, 98.0, 100.0, 105.0, 0.06, 0.09, 0.05, 0.50, 0.11, 0.16, 0.63, 3.5224, 1.0e-4)
, (2, Call, 3.0, 122.0, 120.0, 0.00, 0.00, 0.10, 0.1, 0.20, 0.20, -0.5, 4.7530, 1.0e-3)
, (2, Call, 3.0, 122.0, 120.0, 0.00, 0.00, 0.10, 0.5, 0.25, 0.20, 0.5, 7.0067, 1.0e-3)
, (2, Call, 3.0, 122.0, 120.0, 0.00, 0.00, 0.10, 0.5, 0.20, 0.25, -0.5, 12.1483, 1.0e-3)
]
forM_ rows $ \(basketType, ty, strike, s1, s2, q1, q2, r, t, v1, v2, corr, result, tol) -> do
let exDate = addDays (round (t * 360 :: Double)) evalDate
spot1 <- simpleQuote s1
spot2 <- simpleQuote s2
qQ1 <- simpleQuote q1
qQ2 <- simpleQuote q2
qTS1 <- flatForward (ReferenceDate evalDate) qQ1 dc Continuous Annual
qTS2 <- flatForward (ReferenceDate evalDate) qQ2 dc Continuous Annual
rQ <- simpleQuote r
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
vQ1 <- simpleQuote v1
vQ2 <- simpleQuote v2
volTS1 <- blackConstantVol (CalendarReferenceDate evalDate) cal vQ1 dc
volTS2 <- blackConstantVol (CalendarReferenceDate evalDate) cal vQ2 dc
(analyticEngine, p1, p2) <- case basketType of
2 -> do
bp1 <- blackProcess spot1 rTS volTS1 EulerDiscretization False
bp2 <- blackProcess spot2 rTS volTS2 EulerDiscretization False
gp1 <- asGeneralizedBlackScholesProcess bp1
gp2 <- asGeneralizedBlackScholesProcess bp2
kirk <- kirkEngine bp1 bp2 corr
pure (kirk, gp1, gp2)
_ -> do
gp1 <- blackScholesMertonProcess spot1 qTS1 rTS volTS1 EulerDiscretization False
gp2 <- blackScholesMertonProcess spot2 qTS2 rTS volTS2 EulerDiscretization False
stulz <- stulzEngine gp1 gp2 corr
pure (stulz, gp1, gp2)
let payoff = PlainVanillaPayoff ty strike
basket = case basketType of
0 -> Min (plainVanillaPayoff payoff)
1 -> Max (plainVanillaPayoff payoff)
_ -> Spread (plainVanillaPayoff payoff)
opt <- basketOption basket (European (EuropeanExercise exDate))
setPricingEngine opt analyticEngine
calculated <- npv opt
calculated `shouldSatisfy` closePrec result tol
rhoMatrix <- either error pure (boxedRealMatrix 2 2 [1, corr, corr, 1])
fd2d <- fd2dBlackScholesVanillaEngine p1 p2 corr 50 50 15 0 Hundsdorfer False (-1.0e10)
setPricingEngine opt fd2d
fdCalculated <- npv opt
fdCalculated `shouldSatisfy` closePrec result (0.01 * result)
procArr <- stochasticProcessArray (p1 :| [p2]) rhoMatrix
mc <- mcEuropeanBasketEngine PseudoRandom Statistics procArr Nothing (Just 1) False False (Just 10000) Nothing Nothing 42
setPricingEngine opt mc
mcCalculated <- npv opt
mcCalculated `shouldSatisfy` closePrec result (0.01 * s1)
it "testBarraquandThreeValues: MCEuropeanBasketEngine/MCAmericanBasketEngine reproduce Barraquand-Martineau Table 3" $
Context.keepingSettingsGc $ do
evalDate <- Context.evaluationDate
dc <- dayCounter (Actual360 False)
cal <- calendar TARGET
let rows :: [(OptionType, Double, Double, Double, Double, Double)]
-- optionType=Put, basketType=Max always here (the only live Table-3 rows upstream
-- leaves un-commented); strike, t (months, 30 days/month), corr, euroValue, amValue
rows =
[ (Put, 35.0, 1.0, 0.0, 0.00, 0.00)
, (Put, 40.0, 1.0, 0.0, 0.13, 0.23)
, (Put, 45.0, 1.0, 0.0, 2.26, 5.00)
, (Put, 40.0, 4.0, 0.0, 0.25, 0.44)
, (Put, 45.0, 4.0, 0.0, 1.55, 5.00)
, (Put, 45.0, 7.0, 0.0, 1.41, 5.00)
, (Put, 40.0, 7.0, 0.5, 0.91, 1.19)
]
spot1 <- simpleQuote 40.0
spot2 <- simpleQuote 40.0
spot3 <- simpleQuote 40.0
qQ <- simpleQuote 0.0
qTS <- flatForward (ReferenceDate evalDate) qQ dc Continuous Annual
rQ <- simpleQuote 0.05
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
vQ1 <- simpleQuote 0.20
vQ2 <- simpleQuote 0.30
vQ3 <- simpleQuote 0.50
volTS1 <- blackConstantVol (CalendarReferenceDate evalDate) cal vQ1 dc
volTS2 <- blackConstantVol (CalendarReferenceDate evalDate) cal vQ2 dc
volTS3 <- blackConstantVol (CalendarReferenceDate evalDate) cal vQ3 dc
p1 <- blackScholesMertonProcess spot1 qTS rTS volTS1 EulerDiscretization False
p2 <- blackScholesMertonProcess spot2 qTS rTS volTS2 EulerDiscretization False
p3 <- blackScholesMertonProcess spot3 qTS rTS volTS3 EulerDiscretization False
forM_ rows $ \(ty, strike, t, corr, euroValue, amValue) -> do
let exDate = addDays (round (t * 30 :: Double)) evalDate
payoff = plainVanillaPayoff (PlainVanillaPayoff ty strike)
basket = Max payoff
rhoMatrix <- either error pure (boxedRealMatrix 3 3 [1, corr, corr, corr, 1, corr, corr, corr, 1])
procArr <- stochasticProcessArray (p1 :| [p2, p3]) rhoMatrix
euroOpt <- basketOption basket (European (EuropeanExercise exDate))
mcQuasi <- mcEuropeanBasketEngine LowDiscrepancy Statistics procArr Nothing (Just 1) False False (Just 8091) Nothing Nothing 42
setPricingEngine euroOpt mcQuasi
euroCalculated <- npv euroOpt
euroCalculated `shouldSatisfy` closePrec euroValue (0.01 * 40.0)
amOpt <- basketOption basket (American (Just evalDate) exDate False)
mcLsmc <- mcAmericanBasketEngine PseudoRandom procArr (Just 500) Nothing False True (Just 1000) Nothing Nothing 1 (Just 250) 2 Monomial
setPricingEngine amOpt mcLsmc
amCalculated <- npv amOpt
amCalculated `shouldSatisfy` closePrec amValue (0.01 * 40.0)
it "testTavellaValues: MCAmericanBasketEngine reproduces Tavella's cached three-asset American call value" $
Context.keepingSettingsGc $ do
evalDate <- Context.evaluationDate
dc <- dayCounter (Actual360 False)
cal <- calendar TARGET
spot1 <- simpleQuote 100.0
spot2 <- simpleQuote 100.0
spot3 <- simpleQuote 100.0
qQ <- simpleQuote 0.1
qTS <- flatForward (ReferenceDate evalDate) qQ dc Continuous Annual
rQ <- simpleQuote 0.05
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
vQ1 <- simpleQuote 0.20
vQ2 <- simpleQuote 0.20
vQ3 <- simpleQuote 0.20
volTS1 <- blackConstantVol (CalendarReferenceDate evalDate) cal vQ1 dc
volTS2 <- blackConstantVol (CalendarReferenceDate evalDate) cal vQ2 dc
volTS3 <- blackConstantVol (CalendarReferenceDate evalDate) cal vQ3 dc
p1 <- blackScholesMertonProcess spot1 qTS rTS volTS1 EulerDiscretization False
p2 <- blackScholesMertonProcess spot2 qTS rTS volTS2 EulerDiscretization False
p3 <- blackScholesMertonProcess spot3 qTS rTS volTS3 EulerDiscretization False
rhoMatrix <- either error pure (boxedRealMatrix 3 3 [1, -0.25, 0.25, -0.25, 1, 0.3, 0.25, 0.3, 1])
procArr <- stochasticProcessArray (p1 :| [p2, p3]) rhoMatrix
let exDate = addDays (round (3.0 * 360 :: Double)) evalDate
payoff = plainVanillaPayoff (PlainVanillaPayoff Call 100.0)
basket = Max payoff
opt <- basketOption basket (American (Just evalDate) exDate False)
eng <- mcAmericanBasketEngine PseudoRandom procArr (Just 20) Nothing False True (Just 10000) Nothing Nothing 1 (Just 2500) 2 Monomial
setPricingEngine opt eng
calculated <- npv opt
est <- errorEstimate opt
calculated `shouldSatisfy` closePrec 18.082 (0.01 * 100.0)
est `shouldSatisfy` (\x -> not (isNaN x || isInfinite x))
it "testOneDAmericanValues: single-asset MaxBasketPayoff American reduces to the 1-D put table (sliceOne)" $
Context.keepingSettingsGc $ do
evalDate <- Context.evaluationDate
dc <- dayCounter (Actual360 False)
cal <- calendar TARGET
let rows :: [(Double, Double, Double)]
-- strike=100, r=0.06, t=0.5, vol=0.4 fixed across sliceOne; spot, expected, tol
rows =
[ (80.00, 21.6059, 1.0e-2)
, (85.00, 18.0374, 1.0e-2)
, (90.00, 14.9187, 1.0e-2)
, (95.00, 12.2314, 1.0e-2)
, (100.00, 9.9458, 1.0e-2)
]
spot1 <- simpleQuote 0.0
qQ <- simpleQuote 0.0
qTS <- flatForward (ReferenceDate evalDate) qQ dc Continuous Annual
rQ <- simpleQuote 0.06
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
vQ1 <- simpleQuote 0.4
volTS1 <- blackConstantVol (CalendarReferenceDate evalDate) cal vQ1 dc
p1 <- blackScholesMertonProcess spot1 qTS rTS volTS1 EulerDiscretization False
rhoMatrix <- either error pure (boxedRealMatrix 1 1 [1])
procArr <- stochasticProcessArray (p1 :| []) rhoMatrix
let exDate = addDays (round (0.5 * 360 :: Double)) evalDate
eng <- mcAmericanBasketEngine PseudoRandom procArr (Just 52) Nothing False True (Just 10000) Nothing Nothing 1 (Just 2500) 2 Monomial
forM_ rows $ \(s, expected, tol) -> do
_ <- setValue spot1 s
let payoff = plainVanillaPayoff (PlainVanillaPayoff Put 100.0)
basket = Max payoff
opt <- basketOption basket (American (Just evalDate) exDate False)
setPricingEngine opt eng
calculated <- npv opt
calculated `shouldSatisfy` closePrec expected (tol * s)
it "testOddSamples: MCAmericanBasketEngine survives an odd required-sample count (antithetic off-by-one regression)" $
Context.keepingSettingsGc $ do
evalDate <- Context.evaluationDate
dc <- dayCounter (Actual360 False)
cal <- calendar TARGET
spot1 <- simpleQuote 80.0
qQ <- simpleQuote 0.0
qTS <- flatForward (ReferenceDate evalDate) qQ dc Continuous Annual
rQ <- simpleQuote 0.06
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
vQ1 <- simpleQuote 0.4
volTS1 <- blackConstantVol (CalendarReferenceDate evalDate) cal vQ1 dc
p1 <- blackScholesMertonProcess spot1 qTS rTS volTS1 EulerDiscretization False
rhoMatrix <- either error pure (boxedRealMatrix 1 1 [1])
procArr <- stochasticProcessArray (p1 :| []) rhoMatrix
let exDate = addDays (round (0.5 * 360 :: Double)) evalDate
payoff = plainVanillaPayoff (PlainVanillaPayoff Put 100.0)
basket = Max payoff
opt <- basketOption basket (American (Just evalDate) exDate False)
eng <- mcAmericanBasketEngine PseudoRandom procArr (Just 53) Nothing False True (Just 10001) Nothing Nothing 1 (Just 2500) 2 Monomial
setPricingEngine opt eng
calculated <- npv opt
calculated `shouldSatisfy` closePrec 21.6059 (1.0e-2 * 80.0)
it "testLocalVolatilitySpreadOption: Fd2dBlackScholesVanillaEngine on two Heston-implied local-vol surfaces" $
Context.keepingSettingsGc $ do
let evalDate = 21 `september` 2017
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter (Actual360 False)
maturity <- addPeriod evalDate (3, Months)
rQ <- simpleQuote 0.07
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
qQ <- simpleQuote 0.03
qTS <- flatForward (ReferenceDate evalDate) qQ dc Continuous Annual
let s1Value = 100.0 :: Double
s2Value = 110.0
s1 <- simpleQuote s1Value
s2 <- simpleQuote s2Value
hp1 <- hestonProcess rTS (Just qTS) s1 0.09 1.0 0.06 0.6 (-0.75) QuadraticExponentialMartingale
hm1 <- hestonModel hp1
hp2 <- hestonProcess rTS (Just qTS) s2 0.1 2.0 0.07 0.8 0.85 QuadraticExponentialMartingale
hm2 <- hestonModel hp2
vol1 <- hestonBlackVolSurface hm1 AngledContour 160
vol2 <- hestonBlackVolSurface hm2 AngledContour 160
let rho_ = -0.6
spreadStrike = s2Value - s1Value
bs1 <- blackScholesMertonProcess s1 qTS rTS vol1 EulerDiscretization False
bs2 <- blackScholesMertonProcess s2 qTS rTS vol2 EulerDiscretization False
opt <- basketOption (Spread (plainVanillaPayoff (PlainVanillaPayoff Call spreadStrike))) (European (EuropeanExercise maturity))
eng <- fd2dBlackScholesVanillaEngine bs1 bs2 rho_ 11 11 6 0 Hundsdorfer True 0.25
setPricingEngine opt eng
calculated <- npv opt
calculated `shouldSatisfy` closePrec 2.561 0.01
it "test2DPDEGreeks: Fd2dBlackScholesVanillaEngine's delta/gamma vs. KirkEngine bump-and-reprice" $
Context.keepingSettingsGc $ do
evalDate <- Context.evaluationDate
dc <- dayCounter Actual365FixedStandard
let maturity = addDays 1095 evalDate
let s1 = 100.0 :: Double
s2 = 100.0
rho_ = 0.5
strike = s1 - s2
cal <- calendar TARGET
spot1 <- simpleQuote s1
spot2 <- simpleQuote s2
rQ <- simpleQuote 0.013
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
vQ <- simpleQuote 0.2
volTS <- blackConstantVol (CalendarReferenceDate evalDate) cal vQ dc
p1 <- blackProcess spot1 rTS volTS EulerDiscretization False
p2 <- blackProcess spot2 rTS volTS EulerDiscretization False
gp1 <- asGeneralizedBlackScholesProcess p1
gp2 <- asGeneralizedBlackScholesProcess p2
opt <- basketOption (Spread (plainVanillaPayoff (PlainVanillaPayoff Call strike))) (European (EuropeanExercise maturity))
fd2d <- fd2dBlackScholesVanillaEngine gp1 gp2 rho_ 100 100 50 0 Hundsdorfer False (-1.0e10)
setPricingEngine opt fd2d
calculatedDelta <- delta opt
calculatedGamma <- gamma opt
kirk <- kirkEngine p1 p2 rho_
setPricingEngine opt kirk
npv0 <- npv opt
let eps = 1.0
_ <- setValue spot1 (s1 + eps)
_ <- setValue spot2 (s2 + eps)
npvUp <- npv opt
_ <- setValue spot1 (s1 - eps)
_ <- setValue spot2 (s2 - eps)
npvDown <- npv opt
let expectedDelta = (npvUp - npvDown) / (2 * eps)
expectedGamma = (npvUp + npvDown - 2 * npv0) / (eps * eps)
tol = 0.0005
calculatedDelta `shouldSatisfy` closePrec expectedDelta tol
calculatedGamma `shouldSatisfy` closePrec expectedGamma tol
it "testBjerksundStenslandSpreadEngine: reproduces the cached put value and call-put parity" $
Context.keepingSettingsGc $ do
let evalDate = 1 `march` 2024
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter Actual365FixedStandard
cal <- calendar TARGET
maturity <- addPeriod evalDate (12, Months)
let f1 = 100 :: Double
f2 = 110 :: Double
rho_ = 0.75 :: Double
spreadStrike = 5 :: Double
rQ <- simpleQuote 0.05
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
v1Q <- simpleQuote 0.25
v2Q <- simpleQuote 0.35
vol1TS <- blackConstantVol (CalendarReferenceDate evalDate) cal v1Q dc
vol2TS <- blackConstantVol (CalendarReferenceDate evalDate) cal v2Q dc
s1 <- simpleQuote f1
s2 <- simpleQuote f2
p1 <- blackScholesMertonProcess s1 rTS rTS vol1TS EulerDiscretization False
p2 <- blackScholesMertonProcess s2 rTS rTS vol2TS EulerDiscretization False
bjEngine <- bjerksundStenslandSpreadEngine p1 p2 rho_
putOpt <- basketOption (Spread (plainVanillaPayoff (PlainVanillaPayoff Put spreadStrike))) (European (EuropeanExercise maturity))
setPricingEngine putOpt bjEngine
putNPV <- npv putOpt
putNPV `shouldSatisfy` closePrec 17.850835947276213 1.0e-6
callOpt <- basketOption (Spread (plainVanillaPayoff (PlainVanillaPayoff Call spreadStrike))) (European (EuropeanExercise maturity))
setPricingEngine callOpt bjEngine
callNPV <- npv callOpt
df <- discount rTS (DatePoint maturity) False
((callNPV - putNPV) / df) `shouldSatisfy` closePrec (f1 - f2 - spreadStrike) 1.0e-3
it "testOperatorSplittingSpreadEngine: reproduces the full Kirk-vs-Strang(First/Second) rho table" $
Context.keepingSettingsGc $ do
let evalDate = 1 `march` 2024
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter Actual365FixedStandard
cal <- calendar TARGET
maturity <- addPeriod evalDate (12, Months)
rQ <- simpleQuote 0.05
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
-- forward-adjusted BlackProcess inputs (f1=110*dq1/df, f2=90*dq2/df)
dq1Q <- simpleQuote 0.03
dq2Q <- simpleQuote 0.02
dq1TS <- flatForward (ReferenceDate evalDate) dq1Q dc Continuous Annual
dq2TS <- flatForward (ReferenceDate evalDate) dq2Q dc Continuous Annual
dfR <- discount rTS (DatePoint maturity) False
dq1 <- discount dq1TS (DatePoint maturity) False
dq2 <- discount dq2TS (DatePoint maturity) False
let f1' = 110 * dq1 / dfR
f2' = 90 * dq2 / dfR
v1Q' <- simpleQuote 0.3
v2Q' <- simpleQuote 0.2
vol1TS' <- blackConstantVol (CalendarReferenceDate evalDate) cal v1Q' dc
vol2TS' <- blackConstantVol (CalendarReferenceDate evalDate) cal v2Q' dc
f1Q <- simpleQuote f1'
f2Q <- simpleQuote f2'
bp1' <- blackProcess f1Q rTS vol1TS' EulerDiscretization False
bp2' <- blackProcess f2Q rTS vol2TS' EulerDiscretization False
p1' <- asGeneralizedBlackScholesProcess bp1'
p2' <- asGeneralizedBlackScholesProcess bp2'
let opsRows =
[ (-0.9 :: Double, 18.9323 :: Double, 18.9361 :: Double)
, (-0.7, 18.0092, 18.012)
, (-0.5, 17.0325, 17.0344)
, (-0.4, 16.5211, 16.5227)
, (-0.3, 15.9925, 15.9937)
, (-0.2, 15.4449, 15.4458)
, (-0.1, 14.8762, 14.8768)
, (0.0, 14.284, 14.2843)
, (0.1, 13.6651, 13.6654)
, (0.2, 13.016, 13.0161)
, (0.3, 12.3319, 12.3319)
, (0.4, 11.6067, 11.6067)
, (0.5, 10.8323, 10.8323)
, (0.7, 9.0863, 9.0862)
, (0.9, 6.9148, 6.9134)
]
opsOpt <- basketOption (Spread (plainVanillaPayoff (PlainVanillaPayoff Call 20.0))) (European (EuropeanExercise maturity))
forM_ opsRows $ \(r, exp1, exp2) -> do
e1 <- operatorSplittingSpreadEngine p1' p2' r First
setPricingEngine opsOpt e1
v1 <- npv opsOpt
v1 `shouldSatisfy` closePrec exp1 1.0e-3
e2 <- operatorSplittingSpreadEngine p1' p2' r Second
setPricingEngine opsOpt e2
v2 <- npv opsOpt
v2 `shouldSatisfy` closePrec exp2 5.0e-3
it "testStrangSplittingSpreadEngineVsMathematica: Kirk/OperatorSplitting(First/Second) reproduce cached Mathematica values" $
Context.keepingSettingsGc $ do
let evalDate = 27 `may` 2024
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter Actual365FixedStandard
cal <- calendar TARGET
rTS <- simpleQuote 0.05 >>= \rQ -> flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
vol2TS <- simpleQuote 0.2 >>= \vQ -> blackConstantVol (CalendarReferenceDate evalDate) cal vQ dc
let s1 = 110.0 :: Double
s2 = 90.0 :: Double
-- T, K, vol1, rho_, kirkNPV, strang1, strang2
rows =
[ (5.0, 20, 0.1, 0.6, 15.39520956886349, 15.39641179190707, 15.41992212706643)
, (10.0, 20, 0.1, 0.6, 22.91537136258191, 22.89480115264337, 22.95919510928365)
, (20.0, 20, 0.1, 0.6, 33.69859018569740, 33.59697949481467, 33.73582501903848)
, (1.0, 20, 0.3, 0.6, 10.9751711157804, 10.97662152028116, 10.97661321814579)
, (2.0, 20, 0.3, 0.6, 15.68896063758723, 15.69277461480688, 15.69275497617036)
, (1.0, 10, 0.3, 0.6, 16.10447007803242, 16.10494344785443, 16.10494658134660)
, (1.0, 40, 0.3, 0.6, 4.657519189575983, 4.657079657030094, 4.656973008981588)
, (1.0, 60, 0.3, 0.6, 1.837359067901817, 1.831230481909945, 1.831241843743509)
, (1.0, 20, 0.5, 0.6, 18.79838447214884, 18.79674735337080, 18.79654551825391)
, (1.0, 20, 0.3, -0.9, 20.17112122874686, 20.14780367419582, 20.15151348149147)
, (1.0, 20, 0.3, 0.0, 15.38036208157481, 15.37697052349819, 15.37728179978961)
, (2.0, 20, 0.3, -0.5, 25.80847626931109, 25.77323435009942, 25.77810550213640)
] ::
[(Double, Double, Double, Double, Double, Double, Double)]
forM_ rows $ \(t, strike, vol1, rho_, kirkNPV, strang1, strang2) -> do
let maturityDate = addDays (round (t * 365 :: Double)) evalDate
dr <- discount rTS (DatePoint maturityDate) False
let f1 = s1 / dr
f2 = s2 / dr
vol1TS <- simpleQuote vol1 >>= \vQ -> blackConstantVol (CalendarReferenceDate evalDate) cal vQ dc
f1Q <- simpleQuote f1
f2Q <- simpleQuote f2
p1 <- blackProcess f1Q rTS vol1TS EulerDiscretization False
p2 <- blackProcess f2Q rTS vol2TS EulerDiscretization False
gp1 <- asGeneralizedBlackScholesProcess p1
gp2 <- asGeneralizedBlackScholesProcess p2
opt <- basketOption (Spread (plainVanillaPayoff (PlainVanillaPayoff Call strike))) (European (EuropeanExercise maturityDate))
kirk <- kirkEngine p1 p2 rho_
setPricingEngine opt kirk
kirkCalc <- npv opt
kirkCalc `shouldSatisfy` closePrec kirkNPV (1.0e-4 * abs kirkNPV)
os1 <- operatorSplittingSpreadEngine gp1 gp2 rho_ First
setPricingEngine opt os1
strang1Calc <- npv opt
strang1Calc `shouldSatisfy` closePrec strang1 (1.0e-4 * abs strang1)
os2 <- operatorSplittingSpreadEngine gp1 gp2 rho_ Second
setPricingEngine opt os2
strang2Calc <- npv opt
strang2Calc `shouldSatisfy` closePrec strang2 (1.0e-4 * abs strang2)
it "testPDEvsApproximations: Kirk/BjerksundStensland/OperatorSplitting/Pearson/GaussianCopula track Fd2d across type/rho_/rate/spot" $
Context.keepingSettingsGc $ do
let evalDate = 5 `february` 2024
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter Actual365FixedStandard
cal <- calendar TARGET
maturity <- addPeriod evalDate (6, Months)
let strike = 5.0 :: Double
s1Q <- simpleQuote 100.0
s2Q <- simpleQuote 100.0
rQ <- simpleQuote 0.05
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
v1Q <- simpleQuote 0.25
v2Q <- simpleQuote 0.4
vol1TS <- blackConstantVol (CalendarReferenceDate evalDate) cal v1Q dc
vol2TS <- blackConstantVol (CalendarReferenceDate evalDate) cal v2Q dc
bp1 <- blackProcess s1Q rTS vol1TS EulerDiscretization False
bp2 <- blackProcess s2Q rTS vol2TS EulerDiscretization False
p1 <- asGeneralizedBlackScholesProcess bp1
p2 <- asGeneralizedBlackScholesProcess bp2
diffs <- fmap concat $
forM [Call, Put] $ \ty -> do
opt <- basketOption (Spread (plainVanillaPayoff (PlainVanillaPayoff ty strike))) (European (EuropeanExercise maturity))
fmap concat $
forM [-0.75, 0.0, 0.9] $ \rho_ -> do
kirk <- kirkEngine bp1 bp2 rho_
bs2014 <- bjerksundStenslandSpreadEngine p1 p2 rho_
os1 <- operatorSplittingSpreadEngine p1 p2 rho_ First
os2 <- operatorSplittingSpreadEngine p1 p2 rho_ Second
pearson <- pearsonSpreadEngine p1 p2 rho_ 1.0e-10 10000 8.0
gauss <- gaussianCopulaSpreadEngine p1 p2 rho_ 64
fd2d <- fd2dBlackScholesVanillaEngine p1 p2 rho_ 50 50 15 0 Hundsdorfer False (-1.0e10)
fmap concat $
forM [0.0, 0.05, 0.2] $ \rate -> do
_ <- setValue rQ rate
forM [75.0, 90.0, 100.0, 105.0, 175.0] $ \spot -> do
_ <- setValue s2Q spot
setPricingEngine opt fd2d
fdNPV <- npv opt
setPricingEngine opt kirk
kirkNPV <- npv opt
setPricingEngine opt bs2014
bs2014NPV <- npv opt
setPricingEngine opt os1
os1NPV <- npv opt
setPricingEngine opt os2
os2NPV <- npv opt
setPricingEngine opt pearson
pearsonNPV <- npv opt
setPricingEngine opt gauss
gaussNPV <- npv opt
pure
( kirkNPV - fdNPV
, bs2014NPV - fdNPV
, os1NPV - fdNPV
, os2NPV - fdNPV
, pearsonNPV - fdNPV
, gaussNPV - fdNPV
)
let stdDev :: [Double] -> Double
stdDev xs =
let n = fromIntegral (length xs) :: Double
m = sum xs / n
in sqrt (sum [(x - m) ^ (2 :: Int) | x <- xs] / (n - 1))
(kirkDiffs, bs2014Diffs, os1Diffs, os2Diffs, pearsonDiffs, gaussDiffs) =
unzip6 diffs
stdDev kirkDiffs `shouldSatisfy` (< 0.03)
stdDev bs2014Diffs `shouldSatisfy` (< 0.02)
stdDev os1Diffs `shouldSatisfy` (< 0.02)
stdDev os2Diffs `shouldSatisfy` (< 0.02)
stdDev pearsonDiffs `shouldSatisfy` (< 0.02)
stdDev gaussDiffs `shouldSatisfy` (< 0.02)
it "ChoiBasketEngine/DengLiZhouBasketEngine/SingleFactorBsmBasketEngine self-consistency vs. MCEuropeanBasketEngine" $
Context.keepingSettingsGc $ do
let evalDate = 1 `march` 2024
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter Actual365FixedStandard
cal <- calendar TARGET
maturity <- addPeriod evalDate (12, Months)
-- Choi/DengLiZhou (accept Average or Spread) and SingleFactorBsm (Average only)
-- self-consistency against the already-bound MCEuropeanBasketEngine
let brho = 0.3 :: Double
bq1 <- simpleQuote 0.0
bq2 <- simpleQuote 0.0
bqTS1 <- flatForward (ReferenceDate evalDate) bq1 dc Continuous Annual
bqTS2 <- flatForward (ReferenceDate evalDate) bq2 dc Continuous Annual
brQ <- simpleQuote 0.05
brTS <- flatForward (ReferenceDate evalDate) brQ dc Continuous Annual
bv1 <- simpleQuote 0.3
bv2 <- simpleQuote 0.3
bvolTS1 <- blackConstantVol (CalendarReferenceDate evalDate) cal bv1 dc
bvolTS2 <- blackConstantVol (CalendarReferenceDate evalDate) cal bv2 dc
bs1Q <- simpleQuote 100
bs2Q <- simpleQuote 100
bp1 <- blackScholesMertonProcess bs1Q bqTS1 brTS bvolTS1 EulerDiscretization False
bp2 <- blackScholesMertonProcess bs2Q bqTS2 brTS bvolTS2 EulerDiscretization False
rhoMatrix <- either error pure (boxedRealMatrix 2 2 [1, brho, brho, 1])
procArr <- stochasticProcessArray (bp1 :| [bp2]) rhoMatrix
mc <- mcEuropeanBasketEngine PseudoRandom Statistics procArr (Just 1) Nothing False False (Just 20000) Nothing Nothing 42
choi <- choiBasketEngine (bp1 :| [bp2]) rhoMatrix 10.0 100000 False False
dlz <- dengLiZhouBasketEngine (bp1 :| [bp2]) rhoMatrix
sfb <- singleFactorBsmBasketEngine (bp1 :| [bp2]) 1.0e-8
spreadOpt <- basketOption (Spread (plainVanillaPayoff (PlainVanillaPayoff Call 0.0))) (European (EuropeanExercise maturity))
setPricingEngine spreadOpt mc
mcV <- npv spreadOpt
setPricingEngine spreadOpt choi
choiV <- npv spreadOpt
setPricingEngine spreadOpt dlz
dlzV <- npv spreadOpt
choiV `shouldSatisfy` closePrec mcV (0.02 * mcV)
dlzV `shouldSatisfy` closePrec mcV (0.05 * mcV)
-- SingleFactorBsmBasketEngine assumes every underlying is driven by one common factor, so
-- it is only verified where that assumption actually holds (rho_=1.0)
rhoMatrix1 <- either error pure (boxedRealMatrix 2 2 [1, 1, 1, 1])
procArr1 <- stochasticProcessArray (bp1 :| [bp2]) rhoMatrix1
mc1 <- mcEuropeanBasketEngine PseudoRandom Statistics procArr1 (Just 1) Nothing False False (Just 20000) Nothing Nothing 42
avgOpt <- basketOption (Average (plainVanillaPayoff (PlainVanillaPayoff Call 100.0)) 2) (European (EuropeanExercise maturity))
setPricingEngine avgOpt mc1
mcAvgV1 <- npv avgOpt
setPricingEngine avgOpt sfb
sfbV1 <- npv avgOpt
sfbV1 `shouldSatisfy` closePrec mcAvgV1 (0.02 * mcAvgV1)
it "FdndimBlackScholesVanillaEngine (both grid forms) vs. Fd2dBlackScholesVanillaEngine" $
Context.keepingSettingsGc $ do
let evalDate = 1 `march` 2024
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter Actual365FixedStandard
cal <- calendar TARGET
maturity <- addPeriod evalDate (12, Months)
let rho_ = 0.75 :: Double
rQ <- simpleQuote 0.05
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
dq1Q <- simpleQuote 0.03
dq2Q <- simpleQuote 0.02
dq1TS <- flatForward (ReferenceDate evalDate) dq1Q dc Continuous Annual
dq2TS <- flatForward (ReferenceDate evalDate) dq2Q dc Continuous Annual
dfR <- discount rTS (DatePoint maturity) False
dq1 <- discount dq1TS (DatePoint maturity) False
dq2 <- discount dq2TS (DatePoint maturity) False
let f1' = 110 * dq1 / dfR
f2' = 90 * dq2 / dfR
v1Q' <- simpleQuote 0.3
v2Q' <- simpleQuote 0.2
vol1TS' <- blackConstantVol (CalendarReferenceDate evalDate) cal v1Q' dc
vol2TS' <- blackConstantVol (CalendarReferenceDate evalDate) cal v2Q' dc
f1Q <- simpleQuote f1'
f2Q <- simpleQuote f2'
bp1' <- blackProcess f1Q rTS vol1TS' EulerDiscretization False
bp2' <- blackProcess f2Q rTS vol2TS' EulerDiscretization False
p1' <- asGeneralizedBlackScholesProcess bp1'
p2' <- asGeneralizedBlackScholesProcess bp2'
fd2d <- fd2dBlackScholesVanillaEngine p1' p2' rho_ 100 100 50 0 Hundsdorfer False (-1.0e10)
crossOpt <- basketOption (Spread (plainVanillaPayoff (PlainVanillaPayoff Call 20.0))) (European (EuropeanExercise maturity))
setPricingEngine crossOpt fd2d
fd2dV <- npv crossOpt
rhoMatrix2 <- either error pure (boxedRealMatrix 2 2 [1, rho_, rho_, 1])
fdndim1 <- fdndimBlackScholesVanillaEngine (p1' :| [p2']) rhoMatrix2 (AxisGrids (50 :| [50])) 50 0 Douglas
fdndim2 <- fdndimBlackScholesVanillaEngine (p1' :| [p2']) rhoMatrix2 (UniformGrid 100) 50 0 Douglas
setPricingEngine crossOpt fdndim1
fdndim1V <- npv crossOpt
setPricingEngine crossOpt fdndim2
fdndim2V <- npv crossOpt
fdndim1V `shouldSatisfy` closePrec fd2dV 0.1
fdndim2V `shouldSatisfy` closePrec fd2dV 0.1
-- Ported from hestonmodel.cpp's testAlanLewisReferencePrices (COS-engine case only; the
-- other six engines checked against the same table upstream aren't re-derived here) and
-- testCosHestonEngineTruncation. Alan Lewis's posted reference prices
-- (http://wilmott.com/messageview.cfm?catid=34&threadid=90957) are checked at upstream's own
-- 1e-12 relative tolerance -- the tightest golden-value check in this file -- so 'closePrec'
-- (an absolute-tolerance helper) isn't reused; tolerance is scaled to each expected value
-- inline instead, per CLAUDE.md's "scale a numeric tolerance to the magnitude of the value"
-- rule.
describe "COS Heston engine" $ do
let closeRel expected relTol actual = abs (actual - expected) < relTol * abs expected
it "reproduces hestonmodel.cpp's testAlanLewisReferencePrices" $
Context.keepingSettingsGc $ do
let evalDate = 5 `july` 2002
maturity = 5 `july` 2003
v0 = 0.04; kappa = 4.0; theta_ = 0.25; sigma = 1.0; rho_ = -0.5 :: Double
cases =
[ (80.0 :: Double, 7.958878113256768285213263077598987193482161301733 :: Double, 26.774758743998854221382195325726949201687074848341 :: Double)
, (90.0, 12.017966707346304987709573290236471654992071308187, 20.933349000596710388139445766564068085476194042256)
, (100.0, 17.055270961270109413522653999411000974895436309183, 16.070154917028834278213466703938231827658768230714)
, (110.0, 23.017825898442800538908781834822560777763225722188, 12.132211516709844867860534767549426052805766831181)
, (120.0, 29.811026202682471843340682293165857439167301370697, 9.024913483457835636553375454092357136489051667150)
]
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter Actual365FixedStandard
rQ <- simpleQuote 0.01
qQ <- simpleQuote 0.02
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
qTS <- flatForward (ReferenceDate evalDate) qQ dc Continuous Annual
s0 <- simpleQuote 100.0
hp <- hestonProcess rTS (Just qTS) s0 v0 kappa theta_ sigma rho_ QuadraticExponentialMartingale
hm <- hestonModel hp
engine <- cosHestonEngine hm 20.0 400
forM_ cases $ \(strike, expectedPut, expectedCall) -> do
putOpt <- vanillaOption (PlainVanilla (PlainVanillaPayoff Put strike)) (European (EuropeanExercise maturity))
setPricingEngine putOpt engine
putInst <- asOneAssetOption putOpt
putV <- npv putInst
putV `shouldSatisfy` closeRel expectedPut 1.0e-12
callOpt <- vanillaOption (PlainVanilla (PlainVanillaPayoff Call strike)) (European (EuropeanExercise maturity))
setPricingEngine callOpt engine
callInst <- asOneAssetOption callOpt
callV <- npv callInst
callV `shouldSatisfy` closeRel expectedCall 1.0e-12
it "reproduces hestonmodel.cpp's testCosHestonEngineTruncation (near-zero deep OTM price)" $
Context.keepingSettingsGc $ do
let evalDate = 22 `august` 2022
maturity = 23 `august` 2022
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter Actual365FixedStandard
rQ <- simpleQuote 0.0
qQ <- simpleQuote 0.0
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
qTS <- flatForward (ReferenceDate evalDate) qQ dc Continuous Annual
s0 <- simpleQuote 100.0
hp <- hestonProcess rTS (Just qTS) s0 0.007 0.8 0.007 0.1 (-0.2) QuadraticExponentialMartingale
hm <- hestonModel hp
-- upstream calls COSHestonEngine(model) with no explicit L/N, taking its defaults (16, 200)
engine <- cosHestonEngine hm 16.0 200
opt <- vanillaOption (PlainVanilla (PlainVanillaPayoff Call 200.0)) (European (EuropeanExercise maturity))
setPricingEngine opt engine
optInst <- asOneAssetOption opt
v <- npv optInst
v `shouldSatisfy` closePrec 0.0 1.0e-7
-- Ported from hestonmodel.cpp's testAnalyticPDFHestonEngine (plain-vanilla case only; the
-- digital-via-call-spread case in the same upstream test isn't re-derived here). Self-
-- consistency: the transition-density integration engine must reprice a plain vanilla call
-- to within upstream's tolerance of the semi-analytic AnalyticHestonEngine.
describe "Analytic PDF Heston engine" $
it "reproduces hestonmodel.cpp's testAnalyticPDFHestonEngine plain-vanilla case" $
Context.keepingSettingsGc $ do
let evalDate = 5 `january` 2014
maturity = 5 `july` 2014
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter Actual365FixedStandard
rQ <- simpleQuote 0.07
qQ <- simpleQuote 0.185
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
qTS <- flatForward (ReferenceDate evalDate) qQ dc Continuous Annual
s0 <- simpleQuote 100.0
hp <- hestonProcess rTS (Just qTS) s0 0.1 4.0 0.05 1.0 (-0.5) QuadraticExponentialMartingale
hm <- hestonModel hp
pdfEngine <- analyticPdfHestonEngine hm 1.0e-6 10000
analyticEngine <- analyticHestonEngine hm (IntegrationOrder 178)
toleranceEngine <- analyticHestonEngine hm (IntegrationTolerance 1.0e-8 10000)
forM_ [40.0, 60.0 .. 180.0 :: Double] $ \strike -> do
opt <- vanillaOption (PlainVanilla (PlainVanillaPayoff Call strike)) (European (EuropeanExercise maturity))
optInst <- asOneAssetOption opt
setPricingEngine opt analyticEngine
expected <- npv optInst
setPricingEngine opt toleranceEngine
toleranceValue <- npv optInst
toleranceValue `shouldSatisfy` closePrec expected 1.0e-7
setPricingEngine opt pdfEngine
calculated <- npv optInst
calculated `shouldSatisfy` closePrec expected 3.0e-6
-- Ported from batesmodel.cpp's testAnalyticVsMCPricing (FD-vs-analytic case only; the
-- Monte-Carlo-vs-analytic comparison in the same upstream test isn't re-derived here).
-- Self-consistency: the partial-integro finite-difference engine must reprice within 0.2
-- (upstream's own absolute tolerance) of the semi-analytic BatesEngine, across upstream's
-- four named model fixtures.
describe "FD Bates vanilla engine" $
it "reproduces batesmodel.cpp's testAnalyticVsMCPricing FD-vs-analytic case" $
Context.keepingSettingsGc $ do
let evalDate = 30 `march` 2007
maturity = 30 `march` 2012
strike = 100.0 :: Double
lambda = 2.0; nu = -0.2; delta_ = 0.1 :: Double
cases =
[ ("t'Hout case 1" :: String, 0.04 :: Double, 1.5 :: Double, 0.04 :: Double, 0.3 :: Double, -0.9 :: Double, 0.025 :: Double, 0.0 :: Double)
, ("Ikonen-Toivanen", 0.0625, 5.0, 0.16, 0.9, 0.1, 0.1, 0.0)
, ("Kahl-Jaeckel", 0.16, 1.0, 0.16, 2.0, -0.8, 0.0, 0.0)
, ("Equity case", 0.07, 2.0, 0.04, 0.55, -0.8, 0.03, 0.035)
]
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter ActualActualISDA
forM_ cases $ \(name, v0, kappa, theta_, sigma, rho_, r, q) -> do
rQ <- simpleQuote r
qQ <- simpleQuote q
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
qTS <- flatForward (ReferenceDate evalDate) qQ dc Continuous Annual
s0 <- simpleQuote 100.0
bp <- batesProcess rTS qTS s0 v0 kappa theta_ sigma rho_ lambda nu delta_ QuadraticExponentialMartingale
bm <- batesModel bp
fdEngine <- fdBatesVanillaEngine bm [] 50 100 30 0 Hundsdorfer
analyticEngine <- batesEngine bm (IntegrationOrder 160)
opt <- vanillaOption (PlainVanilla (PlainVanillaPayoff Put strike)) (European (EuropeanExercise maturity))
optInst <- asOneAssetOption opt
setPricingEngine opt analyticEngine
expected <- npv optInst
setPricingEngine opt fdEngine
fdV <- npv optInst
(name, closePrec expected 0.2 fdV) `shouldBe` (name, True)
-- Ported from forwardoption.cpp's testHestonMCPrices, "Test 1": a near-zero-vol-of-vol Heston
-- process (kappa/sigma = 1e-8) is observationally a flat Black-Scholes process, so its
-- MC-priced forward-starting option must reprice within upstream's own per-moneyness
-- tolerance of the closed-form ForwardVanillaEngine<AnalyticEuropeanEngine> price. Only this
-- self-contained sub-case is ported (not the file's second sub-case, which additionally
-- exercises AnalyticHestonForwardEuropeanEngine/AnalyticHestonEngine consistency at reset=0 —
-- out of scope for "MC forward Heston engine" coverage specifically).
describe "MC forward Heston engine" $
it "reproduces forwardoption.cpp's testHestonMCPrices flat-Heston-vs-analytic-BS case" $
Context.keepingSettingsGc $ do
let evalDate = 2 `january` 2024
maturity = addGregorianYearsClip 1 evalDate
reset = addDays 182 evalDate
q = 0.04; r = 0.01; sigmaBs = 0.245; s = 100.0 :: Double
v0 = sigmaBs * sigmaBs; kappa = 1.0e-8; theta_ = sigmaBs * sigmaBs; sigma = 1.0e-8; rho_ = -0.93 :: Double
moneyness = [0.8, 0.9, 1.0, 1.1, 1.2 :: Double]
tolCall = [7.0e-4, 8.0e-4, 6.0e-4, 5.0e-4, 5.0e-4]
tolPut = [6.0e-4, 5.0e-4, 6.0e-4, 1.0e-3, 1.0e-3]
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter (Actual360 False)
rQ <- simpleQuote r
qQ <- simpleQuote q
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
qTS <- flatForward (ReferenceDate evalDate) qQ dc Continuous Annual
volQ <- simpleQuote sigmaBs
cal <- calendar Null
volTS <- blackConstantVol (CalendarReferenceDate evalDate) cal volQ dc
spotQ <- simpleQuote s
bsProcess <- blackScholesMertonProcess spotQ qTS rTS volTS EulerDiscretization False
analyticEngine <- forwardEuropeanEngine bsProcess
hp <- hestonProcess rTS (Just qTS) spotQ v0 kappa theta_ sigma rho_ QuadraticExponentialMartingale
mcEngine <- mcForwardEuropeanHestonEngine LowDiscrepancy Statistics hp (Just 50) Nothing False (Just 4095) Nothing Nothing 42 False
forM_ [(Call, tolCall), (Put, tolPut)] $ \(optType, tols) ->
forM_ (zip moneyness tols) $ \(mny, tol) -> do
let payoff = PlainVanilla (PlainVanillaPayoff optType 0.0)
exercise = European (EuropeanExercise maturity)
opt <- forwardVanillaOption mny reset payoff exercise
setPricingEngine opt analyticEngine
analyticPrice <- npv opt
setPricingEngine opt mcEngine
mcPrice <- npv opt
let relErr = abs (analyticPrice - mcPrice) / s
(mny, relErr <= tol) `shouldBe` (mny, True)
-- Ported from americanoption.cpp's testFDShoutNPV (golden-value table), testZeroVolFDShoutNPV
-- and testLargeDividendShoutNPV (both self-consistency).
describe "FdBlackScholesShoutEngine" $ do
it "reproduces americanoption.cpp's testFDShoutNPV" $
Context.keepingSettingsGc $ do
let evalDate = 4 `february` 2021
cases =
[ (Put, 105.0 :: Double, 19.136 :: Double)
, (Call, 105.0, 28.211)
, (Put, 120.0, 28.02)
, (Call, 80.0, 40.785)
]
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter Actual365FixedStandard
s0 <- simpleQuote 100.0
qQ <- simpleQuote 0.03
rQ <- simpleQuote 0.06
volQ <- simpleQuote 0.25
qTS <- flatForward (ReferenceDate evalDate) qQ dc Continuous Annual
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
tgt <- calendar TARGET
volTS <- blackConstantVol (CalendarReferenceDate evalDate) tgt volQ dc
process <- blackScholesMertonProcess s0 qTS rTS volTS EulerDiscretization False
maturity <- addPeriod evalDate (5, Years)
engine <- fdBlackScholesShoutEngine process [] 400 200 0 Hundsdorfer
forM_ cases $ \(ty, strike, expected) -> do
opt <- vanillaOption (PlainVanilla (PlainVanillaPayoff ty strike)) (American Nothing maturity False)
optInst <- asOneAssetOption opt
setPricingEngine opt engine
v <- npv optInst
v `shouldSatisfy` closePrec expected 2.0e-2
it "reproduces americanoption.cpp's testZeroVolFDShoutNPV (shout with a discrete dividend matches the American NPV once undiscounted through the ex-date)" $
Context.keepingSettingsGc $ do
let evalDate = 14 `february` 2021
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter Actual365FixedStandard
s0 <- simpleQuote 100.0
qQ <- simpleQuote 0.03
rQ <- simpleQuote 0.07
volQ <- simpleQuote 1.0e-6
qTS <- flatForward (ReferenceDate evalDate) qQ dc Continuous Annual
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
tgt <- calendar TARGET
volTS <- blackConstantVol (CalendarReferenceDate evalDate) tgt volQ dc
process <- blackScholesMertonProcess s0 qTS rTS volTS EulerDiscretization False
maturity <- addPeriod evalDate (1, Years)
divDate <- addPeriod evalDate (3, Months)
dividends <- sequence [fixedDividend 10.0 divDate]
americanOpt <- vanillaOption (PlainVanilla (PlainVanillaPayoff Put 100.0)) (American (Just evalDate) maturity False)
americanInst <- asOneAssetOption americanOpt
americanEngine <- fdBlackScholesVanillaEngine process dividends 50 50 1 Douglas False 0.0 CashDividendSpot
setPricingEngine americanOpt americanEngine
americanNPV <- npv americanInst
shoutOpt <- vanillaOption (PlainVanilla (PlainVanillaPayoff Put 100.0)) (American (Just evalDate) maturity False)
shoutInst <- asOneAssetOption shoutOpt
shoutEngine <- fdBlackScholesShoutEngine process dividends 50 50 0 Hundsdorfer
setPricingEngine shoutOpt shoutEngine
shoutNPV <- npv shoutInst
rMaturityDf <- discount rTS (DatePoint maturity) True
rDivDateDf <- discount rTS (DatePoint divDate) True
let df = rMaturityDf / rDivDateDf
(shoutNPV / df) `shouldSatisfy` closePrec americanNPV 1.0e-3
it "reproduces americanoption.cpp's testLargeDividendShoutNPV" $
Context.keepingSettingsGc $ do
let evalDate = 21 `february` 2021
strike = 80.0 :: Double
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter Actual365FixedStandard
s0 <- simpleQuote 100.0
qQ <- simpleQuote 0.0
rQ <- simpleQuote 0.0
volQ <- simpleQuote 0.25
qTS <- flatForward (ReferenceDate evalDate) qQ dc Continuous Annual
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
tgt <- calendar TARGET
volTS <- blackConstantVol (CalendarReferenceDate evalDate) tgt volQ dc
process <- blackScholesMertonProcess s0 qTS rTS volTS EulerDiscretization False
maturity <- addPeriod evalDate (6, Months)
divDate <- addPeriod evalDate (3, Months)
dividends <- sequence [fixedDividend 30.0 divDate]
opt <- vanillaOption (PlainVanilla (PlainVanillaPayoff Call strike)) (American (Just evalDate) maturity False)
optInst <- asOneAssetOption opt
engine <- fdBlackScholesShoutEngine process dividends 100 400 0 Hundsdorfer
setPricingEngine opt engine
calculated <- npv optInst
refOpt <- vanillaOption (PlainVanilla (PlainVanillaPayoff Call strike)) (American (Just evalDate) divDate False)
refInst <- asOneAssetOption refOpt
refEngine <- fdBlackScholesShoutEngine process [] 100 400 0 Hundsdorfer
setPricingEngine refOpt refEngine
refNPV <- npv refInst
rMaturityDf <- discount rTS (DatePoint maturity) True
rDivDateDf <- discount rTS (DatePoint divDate) True
let expected = refNPV * rMaturityDf / rDivDateDf
calculated `shouldSatisfy` closePrec expected 5.0e-2
-- Uses hestonslvmodel.cpp's model-construction fixture and FDM parameters, but
-- a flat local vol and a plain Heston parameter set in place of upstream's
-- createSmoothImpliedVol surface and its implied-calibrated parameters. The density-grid
-- shape check pins the 'RealMatrix' layout (rows = varianceGrid length, cols = spotGrid length).
-- The @logging = False@ case also pins empty-log handling; see 'hestonSlvFdmLogEntries'.
describe "HestonSLV model" $ do
it "builds MC/FDM Heston-SLV models with a consistent density-grid layout (LONG)" $
Context.keepingSettingsGc $ do
let evalDate = 5 `march` 2016
Context.setEvaluationDate (Just evalDate)
dc <- dayCounter Actual365FixedStandard
rQ <- simpleQuote 0.01
qQ <- simpleQuote 0.02
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
qTS <- flatForward (ReferenceDate evalDate) qQ dc Continuous Annual
s0 <- simpleQuote 100.0
localVolQ <- simpleQuote 0.3
hp <- hestonProcess rTS (Just qTS) s0 0.09 1.0 0.06 0.4 (-0.75) HestonFullTruncation
hm <- hestonModel hp
end <- addPeriod evalDate (1, Years)
localVolTS <- localConstantVol (ReferenceDate evalDate) localVolQ dc
factory <- sobolBrownianGeneratorFactory Diagonal 1234 JoeKuoD7
mc <- hestonSlvMcModel localVolTS hm factory end 91 201 32768 [] 1.0
mcLeverage <- leverageFunction mc
slv <- hestonSlvProcess hp mcLeverage 1.0
n <- factors slv
n `shouldBe` 2
let fdmParams = HestonSLVFokkerPlanckFdmParams
51 151 500 50 100.0 5 2 0.1 1.0e-4 10000
1.0e-5 1.0e-5 2.5e-6 1.0 0.1 0.9 1.0e-5
ZeroCorrelation Log ModifiedCraigSneyd
fdm <- hestonSlvFdmModel localVolTS hm end fdmParams True [] 1.0
fdmLeverage <- leverageFunction fdm
fdmVol <- localVol fdmLeverage end 100 True
fdmVol `shouldSatisfy` (\v -> v > 0 && not (isNaN v || isInfinite v))
logs <- hestonSlvFdmLogEntries fdm
case logs of
[] -> expectationFailure "FDM logging produced no diagnostic snapshots"
entry : _ -> do
let density = hestonSlvLogDensity entry
nVar = V.length (hestonSlvLogVarianceCoordinates entry)
nSpot = V.length (hestonSlvLogSpotCoordinates entry)
realMatrixRows density `shouldBe` fromIntegral nVar
realMatrixColumns density `shouldBe` fromIntegral nSpot
V.length (realMatrixData density) `shouldBe` nVar * nSpot
fdmNoLog <- hestonSlvFdmModel localVolTS hm end fdmParams False [] 1.0
noLogs <- hestonSlvFdmLogEntries fdmNoLog
noLogs `shouldBe` []
-- Ported from hestonslvmodel.cpp's testMonteCarloVsFdmPricing: the FD Heston-SLV engine's
-- price must be unaffected by a "mixing factor" applied to a differently-parameterized
-- Heston model paired with the same leverage function (mixingFactor scales the mixing
-- model's contribution toward zero, so the two engines are constructed to price
-- identically by upstream's own design). The Monte-Carlo leg of the same upstream test
-- (MCEuropeanHestonEngine<..., HestonSLVProcess>) isn't ported: 'HestonSLVProcess' is a
-- 'GenStochasticProcess' leaf outside the 'GenHestonProcess' family hasquant's
-- 'mcEuropeanHestonEngine' requires, so that specific engine/process combination isn't
-- constructible from hasquant's current bindings.
it "reproduces hestonslvmodel.cpp's testMonteCarloVsFdmPricing mixing-factor FDM consistency (LONG)" $
Context.keepingSettingsGc $ do
let evalDate = 5 `december` 2015
v0 = 0.19; kappa = 2.0; theta_ = 0.18; sigma = 0.8; rho_ = -0.75 :: Double
strikes = [100.0, 110.0 :: Double]
Context.setEvaluationDate (Just evalDate)
maturity <- addPeriod evalDate (1, Years)
dc <- dayCounter ActualActualISDA
s0 <- simpleQuote 100.0
rQ <- simpleQuote 0.05
qQ <- simpleQuote 0.02
rTS <- flatForward (ReferenceDate evalDate) rQ dc Continuous Annual
qTS <- flatForward (ReferenceDate evalDate) qQ dc Continuous Annual
hp <- hestonProcess rTS (Just qTS) s0 v0 kappa theta_ sigma rho_ QuadraticExponentialMartingale
hm <- hestonModel hp
leverageQ <- simpleQuote 0.25
leverageFct <- localConstantVol (ReferenceDate evalDate) leverageQ dc
fdEngine <- fdHestonVanillaEngine hm [] 51 401 101 0 ModifiedCraigSneyd (Just leverageFct) 1.0
mixHp <- hestonProcess rTS (Just qTS) s0 v0 kappa theta_ (sigma * 10) rho_ QuadraticExponentialMartingale
mixHm <- hestonModel mixHp
fdEngineMix <- fdHestonVanillaEngine mixHm [] 51 401 101 0 ModifiedCraigSneyd (Just leverageFct) 0.1
forM_ strikes $ \strike -> do
opt <- vanillaOption (PlainVanilla (PlainVanillaPayoff Call strike)) (European (EuropeanExercise maturity))
optInst <- asOneAssetOption opt
setPricingEngine opt fdEngine
priceFDM <- npv optInst
setPricingEngine opt fdEngineMix
priceFDMWithMix <- npv optInst
(strike, priceFDMWithMix) `shouldBe` (strike, priceFDM)