hq (empty) → 0.1.0.0
raw patch · 46 files changed
+3673/−0 lines, 46 filesdep +basedep +bytestringdep +cassavasetup-changed
Dependencies added: base, bytestring, cassava, containers, conversion, data-default-class, erf, hmatrix, hmatrix-gsl, hmatrix-gsl-stats, hq, hspec, hspec-expectations, ieee754, math-functions, mersenne-random-pure64, monad-loops, mtl, random, random-fu, random-source, rvar, sorted-list, statistics, stm, text, time, vector, vector-algorithms
Files
- ChangeLog.md +3/−0
- LICENSE +30/−0
- README.md +1/−0
- Setup.hs +2/−0
- external/src/erf_cody.cpp +455/−0
- external/src/lets_be_rational.cpp +633/−0
- external/src/normaldistribution.cpp +147/−0
- external/src/rationalcubic.cpp +115/−0
- hq.cabal +186/−0
- src/Q/ContingentClaim.hs +68/−0
- src/Q/ContingentClaim/Options.hs +48/−0
- src/Q/Currencies/America.hs +25/−0
- src/Q/Currencies/Asia.hs +16/−0
- src/Q/Currencies/Europe.hs +34/−0
- src/Q/Currency.hs +20/−0
- src/Q/Greeks.hs +60/−0
- src/Q/Interpolation.hs +22/−0
- src/Q/MonteCarlo.hs +87/−0
- src/Q/Options.hs +35/−0
- src/Q/Options/Bachelier.hs +56/−0
- src/Q/Options/Black76.hs +60/−0
- src/Q/Options/BlackScholes.hs +85/−0
- src/Q/Options/ImpliedVol.hs +70/−0
- src/Q/Options/ImpliedVol/InterpolatingSmile.hs +27/−0
- src/Q/Options/ImpliedVol/LetsBeRational.hs +22/−0
- src/Q/Options/ImpliedVol/Normal.hs +127/−0
- src/Q/Options/ImpliedVol/SVI.hs +35/−0
- src/Q/Options/ImpliedVol/StrikeInterpolation.hs +31/−0
- src/Q/Options/ImpliedVol/Surface.hs +188/−0
- src/Q/Options/ImpliedVol/TimeInterpolation.hs +7/−0
- src/Q/Options/ImpliedVol/TimeSlice.hs +30/−0
- src/Q/Payoff.hs +47/−0
- src/Q/Plotting.hs +10/−0
- src/Q/SortedVector.hs +26/−0
- src/Q/Stats/Arima.hs +55/−0
- src/Q/Stats/TimeSeries.hs +82/−0
- src/Q/Stochastic.hs +5/−0
- src/Q/Stochastic/Discretize.hs +37/−0
- src/Q/Stochastic/Process.hs +97/−0
- src/Q/Time.hs +53/−0
- src/Q/Time/Date.hs +60/−0
- src/Q/Time/DayCounter.hs +53/−0
- src/Q/Types.hs +277/−0
- src/Q/Util/File.hs +32/−0
- test/bachelier/Spec.hs +62/−0
- test/normalimpliedvol/Spec.hs +52/−0
+ ChangeLog.md view
@@ -0,0 +1,3 @@+# Changelog for hq++## Unreleased changes
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Author name here (c) 2020++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Author name here nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,1 @@+# hq
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ external/src/erf_cody.cpp view
@@ -0,0 +1,455 @@+//+// Original Fortran code taken from http://www.netlib.org/specfun/erf, compiled with f2c, and adapted by hand.+//+// Created with command line f2c -C++ -c -a -krd -r8 cody_erf.f+//+// Translated by f2c (version 20100827).+//++//+// This source code resides at www.jaeckel.org/LetsBeRational.7z .+//+// ======================================================================================+// WARRANTY DISCLAIMER+// The Software is provided "as is" without warranty of any kind, either express or implied,+// including without limitation any implied warranties of condition, uninterrupted use,+// merchantability, fitness for a particular purpose, or non-infringement.+// ======================================================================================+//++#if defined( _DEBUG ) || defined( BOUNDS_CHECK_STL_ARRAYS )+#define _SECURE_SCL 1+#define _SECURE_SCL_THROWS 1+#define _SCL_SECURE_NO_WARNINGS+#define _HAS_ITERATOR_DEBUGGING 0+#else+#define _SECURE_SCL 0+#endif+#if defined(_MSC_VER)+# define NOMINMAX // to suppress MSVC's definitions of min() and max()+// These four pragmas are the equivalent to /fp:fast.+# pragma float_control( except, off )+# pragma float_control( precise, off )+# pragma fp_contract( on )+# pragma fenv_access( off )+#endif++#include "normaldistribution.h"+#include <math.h>+#include <float.h>++namespace {+ inline double d_int(const double x){ return( (x>0) ? floor(x) : -floor(-x) ); }+}++/*< SUBROUTINE CALERF(ARG,RESULT,JINT) >*/+double calerf(double x, const int jint) {++ static const double a[5] = { 3.1611237438705656,113.864154151050156,377.485237685302021,3209.37758913846947,.185777706184603153 };+ static const double b[4] = { 23.6012909523441209,244.024637934444173,1282.61652607737228,2844.23683343917062 };+ static const double c__[9] = { .564188496988670089,8.88314979438837594,66.1191906371416295,298.635138197400131,881.95222124176909,1712.04761263407058,2051.07837782607147,1230.33935479799725,2.15311535474403846e-8 };+ static const double d__[8] = { 15.7449261107098347,117.693950891312499,537.181101862009858,1621.38957456669019,3290.79923573345963,4362.61909014324716,3439.36767414372164,1230.33935480374942 };+ static const double p[6] = { .305326634961232344,.360344899949804439,.125781726111229246,.0160837851487422766,6.58749161529837803e-4,.0163153871373020978 };+ static const double q[5] = { 2.56852019228982242,1.87295284992346047,.527905102951428412,.0605183413124413191,.00233520497626869185 };++ static const double zero = 0.;+ static const double half = .5;+ static const double one = 1.;+ static const double two = 2.;+ static const double four = 4.;+ static const double sqrpi = 0.56418958354775628695;+ static const double thresh = .46875;+ static const double sixten = 16.;++ double y, del, ysq, xden, xnum, result;++ /* ------------------------------------------------------------------ */+ /* This packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x) */+ /* for a real argument x. It contains three FUNCTION type */+ /* subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX), */+ /* and one SUBROUTINE type subprogram, CALERF. The calling */+ /* statements for the primary entries are: */+ /* Y=ERF(X) (or Y=DERF(X)), */+ /* Y=ERFC(X) (or Y=DERFC(X)), */+ /* and */+ /* Y=ERFCX(X) (or Y=DERFCX(X)). */+ /* The routine CALERF is intended for internal packet use only, */+ /* all computations within the packet being concentrated in this */+ /* routine. The function subprograms invoke CALERF with the */+ /* statement */+ /* CALL CALERF(ARG,RESULT,JINT) */+ /* where the parameter usage is as follows */+ /* Function Parameters for CALERF */+ /* call ARG Result JINT */+ /* ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0 */+ /* ERFC(ARG) ABS(ARG) .LT. XBIG ERFC(ARG) 1 */+ /* ERFCX(ARG) XNEG .LT. ARG .LT. XMAX ERFCX(ARG) 2 */+ /* The main computation evaluates near-minimax approximations */+ /* from "Rational Chebyshev approximations for the error function" */+ /* by W. J. Cody, Math. Comp., 1969, PP. 631-638. This */+ /* transportable program uses rational functions that theoretically */+ /* approximate erf(x) and erfc(x) to at least 18 significant */+ /* decimal digits. The accuracy achieved depends on the arithmetic */+ /* system, the compiler, the intrinsic functions, and proper */+ /* selection of the machine-dependent constants. */+ /* ******************************************************************* */+ /* ******************************************************************* */+ /* Explanation of machine-dependent constants */+ /* XMIN = the smallest positive floating-point number. */+ /* XINF = the largest positive finite floating-point number. */+ /* XNEG = the largest negative argument acceptable to ERFCX; */+ /* the negative of the solution to the equation */+ /* 2*exp(x*x) = XINF. */+ /* XSMALL = argument below which erf(x) may be represented by */+ /* 2*x/sqrt(pi) and above which x*x will not underflow. */+ /* A conservative value is the largest machine number X */+ /* such that 1.0 + X = 1.0 to machine precision. */+ /* XBIG = largest argument acceptable to ERFC; solution to */+ /* the equation: W(x) * (1-0.5/x**2) = XMIN, where */+ /* W(x) = exp(-x*x)/[x*sqrt(pi)]. */+ /* XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to */+ /* machine precision. A conservative value is */+ /* 1/[2*sqrt(XSMALL)] */+ /* XMAX = largest acceptable argument to ERFCX; the minimum */+ /* of XINF and 1/[sqrt(pi)*XMIN]. */+ // The numbers below were preselected for IEEE .+ static const double xinf = 1.79e308;+ static const double xneg = -26.628;+ static const double xsmall = 1.11e-16;+ static const double xbig = 26.543;+ static const double xhuge = 6.71e7;+ static const double xmax = 2.53e307;+ /* Approximate values for some important machines are: */+ /* XMIN XINF XNEG XSMALL */+ /* CDC 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15 */+ /* CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15 */+ /* IEEE (IBM/XT, */+ /* SUN, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8 */+ /* IEEE (IBM/XT, */+ /* SUN, etc.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16 */+ /* IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17 */+ /* UNIVAC 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18 */+ /* VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17 */+ /* VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16 */+ /* XBIG XHUGE XMAX */+ /* CDC 7600 (S.P.) 25.922 8.39E+6 1.80X+293 */+ /* CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465 */+ /* IEEE (IBM/XT, */+ /* SUN, etc.) (S.P.) 9.194 2.90E+3 4.79E+37 */+ /* IEEE (IBM/XT, */+ /* SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307 */+ /* IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75 */+ /* UNIVAC 1108 (D.P.) 26.582 5.37D+8 8.98D+307 */+ /* VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38 */+ /* VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307 */+ /* ******************************************************************* */+ /* ******************************************************************* */+ /* Error returns */+ /* The program returns ERFC = 0 for ARG .GE. XBIG; */+ /* ERFCX = XINF for ARG .LT. XNEG; */+ /* and */+ /* ERFCX = 0 for ARG .GE. XMAX. */+ /* Intrinsic functions required are: */+ /* ABS, AINT, EXP */+ /* Author: W. J. Cody */+ /* Mathematics and Computer Science Division */+ /* Argonne National Laboratory */+ /* Argonne, IL 60439 */+ /* Latest modification: March 19, 1990 */+ /* ------------------------------------------------------------------ */+ /*< INTEGER I,JINT >*/+ /* S REAL */+ /*< >*/+ /*< DIMENSION A(5),B(4),C(9),D(8),P(6),Q(5) >*/+ /* ------------------------------------------------------------------ */+ /* Mathematical constants */+ /* ------------------------------------------------------------------ */+ /* S DATA FOUR,ONE,HALF,TWO,ZERO/4.0E0,1.0E0,0.5E0,2.0E0,0.0E0/, */+ /* S 1 SQRPI/5.6418958354775628695E-1/,THRESH/0.46875E0/, */+ /* S 2 SIXTEN/16.0E0/ */+ /*< >*/+ /* ------------------------------------------------------------------ */+ /* Machine-dependent constants */+ /* ------------------------------------------------------------------ */+ /* S DATA XINF,XNEG,XSMALL/3.40E+38,-9.382E0,5.96E-8/, */+ /* S 1 XBIG,XHUGE,XMAX/9.194E0,2.90E3,4.79E37/ */+ /*< >*/+ /* ------------------------------------------------------------------ */+ /* Coefficients for approximation to erf in first interval */+ /* ------------------------------------------------------------------ */+ /* S DATA A/3.16112374387056560E00,1.13864154151050156E02, */+ /* S 1 3.77485237685302021E02,3.20937758913846947E03, */+ /* S 2 1.85777706184603153E-1/ */+ /* S DATA B/2.36012909523441209E01,2.44024637934444173E02, */+ /* S 1 1.28261652607737228E03,2.84423683343917062E03/ */+ /*< >*/+ /*< >*/+ /* ------------------------------------------------------------------ */+ /* Coefficients for approximation to erfc in second interval */+ /* ------------------------------------------------------------------ */+ /* S DATA C/5.64188496988670089E-1,8.88314979438837594E0, */+ /* S 1 6.61191906371416295E01,2.98635138197400131E02, */+ /* S 2 8.81952221241769090E02,1.71204761263407058E03, */+ /* S 3 2.05107837782607147E03,1.23033935479799725E03, */+ /* S 4 2.15311535474403846E-8/ */+ /* S DATA D/1.57449261107098347E01,1.17693950891312499E02, */+ /* S 1 5.37181101862009858E02,1.62138957456669019E03, */+ /* S 2 3.29079923573345963E03,4.36261909014324716E03, */+ /* S 3 3.43936767414372164E03,1.23033935480374942E03/ */+ /*< >*/+ /*< >*/+ /* ------------------------------------------------------------------ */+ /* Coefficients for approximation to erfc in third interval */+ /* ------------------------------------------------------------------ */+ /* S DATA P/3.05326634961232344E-1,3.60344899949804439E-1, */+ /* S 1 1.25781726111229246E-1,1.60837851487422766E-2, */+ /* S 2 6.58749161529837803E-4,1.63153871373020978E-2/ */+ /* S DATA Q/2.56852019228982242E00,1.87295284992346047E00, */+ /* S 1 5.27905102951428412E-1,6.05183413124413191E-2, */+ /* S 2 2.33520497626869185E-3/ */+ /*< >*/+ /*< >*/+ /* ------------------------------------------------------------------ */+ /*< X = ARG >*/+ // x = *arg;+ /*< Y = ABS(X) >*/+ y = fabs(x);+ /*< IF (Y .LE. THRESH) THEN >*/+ if (y <= thresh) {+ /* ------------------------------------------------------------------ */+ /* Evaluate erf for |X| <= 0.46875 */+ /* ------------------------------------------------------------------ */+ /*< YSQ = ZERO >*/+ ysq = zero;+ /*< IF (Y .GT. XSMALL) YSQ = Y * Y >*/+ if (y > xsmall) {+ ysq = y * y;+ }+ /*< XNUM = A(5)*YSQ >*/+ xnum = a[4] * ysq;+ /*< XDEN = YSQ >*/+ xden = ysq;+ /*< DO 20 I = 1, 3 >*/+ for (int i__ = 1; i__ <= 3; ++i__) {+ /*< XNUM = (XNUM + A(I)) * YSQ >*/+ xnum = (xnum + a[i__ - 1]) * ysq;+ /*< XDEN = (XDEN + B(I)) * YSQ >*/+ xden = (xden + b[i__ - 1]) * ysq;+ /*< 20 CONTINUE >*/+ /* L20: */+ }+ /*< RESULT = X * (XNUM + A(4)) / (XDEN + B(4)) >*/+ result = x * (xnum + a[3]) / (xden + b[3]);+ /*< IF (JINT .NE. 0) RESULT = ONE - RESULT >*/+ if (jint != 0) {+ result = one - result;+ }+ /*< IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT >*/+ if (jint == 2) {+ result = exp(ysq) * result;+ }+ /*< GO TO 800 >*/+ goto L800;+ /* ------------------------------------------------------------------ */+ /* Evaluate erfc for 0.46875 <= |X| <= 4.0 */+ /* ------------------------------------------------------------------ */+ /*< ELSE IF (Y .LE. FOUR) THEN >*/+ } else if (y <= four) {+ /*< XNUM = C(9)*Y >*/+ xnum = c__[8] * y;+ /*< XDEN = Y >*/+ xden = y;+ /*< DO 120 I = 1, 7 >*/+ for (int i__ = 1; i__ <= 7; ++i__) {+ /*< XNUM = (XNUM + C(I)) * Y >*/+ xnum = (xnum + c__[i__ - 1]) * y;+ /*< XDEN = (XDEN + D(I)) * Y >*/+ xden = (xden + d__[i__ - 1]) * y;+ /*< 120 CONTINUE >*/+ /* L120: */+ }+ /*< RESULT = (XNUM + C(8)) / (XDEN + D(8)) >*/+ result = (xnum + c__[7]) / (xden + d__[7]);+ /*< IF (JINT .NE. 2) THEN >*/+ if (jint != 2) {+ /*< YSQ = AINT(Y*SIXTEN)/SIXTEN >*/+ double d__1 = y * sixten;+ ysq = d_int(d__1) / sixten;+ /*< DEL = (Y-YSQ)*(Y+YSQ) >*/+ del = (y - ysq) * (y + ysq);+ /*< RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT >*/+ d__1 = exp(-ysq * ysq) * exp(-del);+ result = d__1 * result;+ /*< END IF >*/+ }+ /* ------------------------------------------------------------------ */+ /* Evaluate erfc for |X| > 4.0 */+ /* ------------------------------------------------------------------ */+ /*< ELSE >*/+ } else {+ /*< RESULT = ZERO >*/+ result = zero;+ /*< IF (Y .GE. XBIG) THEN >*/+ if (y >= xbig) {+ /*< IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 300 >*/+ if (jint != 2 || y >= xmax) {+ goto L300;+ }+ /*< IF (Y .GE. XHUGE) THEN >*/+ if (y >= xhuge) {+ /*< RESULT = SQRPI / Y >*/+ result = sqrpi / y;+ /*< GO TO 300 >*/+ goto L300;+ /*< END IF >*/+ }+ /*< END IF >*/+ }+ /*< YSQ = ONE / (Y * Y) >*/+ ysq = one / (y * y);+ /*< XNUM = P(6)*YSQ >*/+ xnum = p[5] * ysq;+ /*< XDEN = YSQ >*/+ xden = ysq;+ /*< DO 240 I = 1, 4 >*/+ for (int i__ = 1; i__ <= 4; ++i__) {+ /*< XNUM = (XNUM + P(I)) * YSQ >*/+ xnum = (xnum + p[i__ - 1]) * ysq;+ /*< XDEN = (XDEN + Q(I)) * YSQ >*/+ xden = (xden + q[i__ - 1]) * ysq;+ /*< 240 CONTINUE >*/+ /* L240: */+ }+ /*< RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5)) >*/+ result = ysq * (xnum + p[4]) / (xden + q[4]);+ /*< RESULT = (SQRPI - RESULT) / Y >*/+ result = (sqrpi - result) / y;+ /*< IF (JINT .NE. 2) THEN >*/+ if (jint != 2) {+ /*< YSQ = AINT(Y*SIXTEN)/SIXTEN >*/+ double d__1 = y * sixten;+ ysq = d_int(d__1) / sixten;+ /*< DEL = (Y-YSQ)*(Y+YSQ) >*/+ del = (y - ysq) * (y + ysq);+ /*< RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT >*/+ d__1 = exp(-ysq * ysq) * exp(-del);+ result = d__1 * result;+ /*< END IF >*/+ }+ /*< END IF >*/+ }+ /* ------------------------------------------------------------------ */+ /* Fix up for negative argument, erf, etc. */+ /* ------------------------------------------------------------------ */+ /*< 300 IF (JINT .EQ. 0) THEN >*/+L300:+ if (jint == 0) {+ /*< RESULT = (HALF - RESULT) + HALF >*/+ result = (half - result) + half;+ /*< IF (X .LT. ZERO) RESULT = -RESULT >*/+ if (x < zero) {+ result = -(result);+ }+ /*< ELSE IF (JINT .EQ. 1) THEN >*/+ } else if (jint == 1) {+ /*< IF (X .LT. ZERO) RESULT = TWO - RESULT >*/+ if (x < zero) {+ result = two - result;+ }+ /*< ELSE >*/+ } else {+ /*< IF (X .LT. ZERO) THEN >*/+ if (x < zero) {+ /*< IF (X .LT. XNEG) THEN >*/+ if (x < xneg) {+ /*< RESULT = XINF >*/+ result = xinf;+ /*< ELSE >*/+ } else {+ /*< YSQ = AINT(X*SIXTEN)/SIXTEN >*/+ double d__1 = x * sixten;+ ysq = d_int(d__1) / sixten;+ /*< DEL = (X-YSQ)*(X+YSQ) >*/+ del = (x - ysq) * (x + ysq);+ /*< Y = EXP(YSQ*YSQ) * EXP(DEL) >*/+ y = exp(ysq * ysq) * exp(del);+ /*< RESULT = (Y+Y) - RESULT >*/+ result = y + y - result;+ /*< END IF >*/+ }+ /*< END IF >*/+ }+ /*< END IF >*/+ }+ /*< 800 RETURN >*/+L800:+ return result;+ /* ---------- Last card of CALERF ---------- */+ /*< END >*/+} /* calerf_ */++/* S REAL FUNCTION ERF(X) */+/*< DOUBLE PRECISION FUNCTION DERF(X) >*/+double erf_cody(double x){+ /* -------------------------------------------------------------------- */+ /* This subprogram computes approximate values for erf(x). */+ /* (see comments heading CALERF). */+ /* Author/date: W. J. Cody, January 8, 1985 */+ /* -------------------------------------------------------------------- */+ /*< INTEGER JINT >*/+ /* S REAL X, RESULT */+ /*< DOUBLE PRECISION X, RESULT >*/+ /* ------------------------------------------------------------------ */+ /*< JINT = 0 >*/+ /*< CALL CALERF(X,RESULT,JINT) >*/+ return calerf(x, 0);+ /* S ERF = RESULT */+ /*< DERF = RESULT >*/+ /*< RETURN >*/+ /* ---------- Last card of DERF ---------- */+ /*< END >*/+} /* derf_ */++/* S REAL FUNCTION ERFC(X) */+/*< DOUBLE PRECISION FUNCTION DERFC(X) >*/+double erfc_cody(double x) {+ /* -------------------------------------------------------------------- */+ /* This subprogram computes approximate values for erfc(x). */+ /* (see comments heading CALERF). */+ /* Author/date: W. J. Cody, January 8, 1985 */+ /* -------------------------------------------------------------------- */+ /*< INTEGER JINT >*/+ /* S REAL X, RESULT */+ /*< DOUBLE PRECISION X, RESULT >*/+ /* ------------------------------------------------------------------ */+ /*< JINT = 1 >*/+ /*< CALL CALERF(X,RESULT,JINT) >*/+ return calerf(x, 1);+ /* S ERFC = RESULT */+ /*< DERFC = RESULT >*/+ /*< RETURN >*/+ /* ---------- Last card of DERFC ---------- */+ /*< END >*/+} /* derfc_ */++/* S REAL FUNCTION ERFCX(X) */+/*< DOUBLE PRECISION FUNCTION DERFCX(X) >*/+double erfcx_cody(double x) {+ /* ------------------------------------------------------------------ */+ /* This subprogram computes approximate values for exp(x*x) * erfc(x). */+ /* (see comments heading CALERF). */+ /* Author/date: W. J. Cody, March 30, 1987 */+ /* ------------------------------------------------------------------ */+ /*< INTEGER JINT >*/+ /* S REAL X, RESULT */+ /*< DOUBLE PRECISION X, RESULT >*/+ /* ------------------------------------------------------------------ */+ /*< JINT = 2 >*/+ /*< CALL CALERF(X,RESULT,JINT) >*/+ return calerf(x, 2);+ /* S ERFCX = RESULT */+ /*< DERFCX = RESULT >*/+ /*< RETURN >*/+ /* ---------- Last card of DERFCX ---------- */+ /*< END >*/+} /* derfcx_ */
+ external/src/lets_be_rational.cpp view
@@ -0,0 +1,633 @@+//+// This source code resides at www.jaeckel.org/LetsBeRational.7z .+//+// ======================================================================================+// Copyright © 2013-2017 Peter Jäckel.+// +// Permission to use, copy, modify, and distribute this software is freely granted,+// provided that this notice is preserved.+//+// WARRANTY DISCLAIMER+// The Software is provided "as is" without warranty of any kind, either express or implied,+// including without limitation any implied warranties of condition, uninterrupted use,+// merchantability, fitness for a particular purpose, or non-infringement.+// ======================================================================================+//++#include "lets_be_rational.h"+// To cross-compile on a command line, you could just use something like+//+// i686-w64-mingw32-g++ -w -fpermissive -shared -DNDEBUG -O3 erf_cody.cpp rationalcubic.cpp normaldistribution.cpp lets_be_rational.cpp xlcall.cpp excel_registration.cpp xlcall32.lib -o lets_be_rational.xll -static-libstdc++ -static-libgcc -s+//+// To compile into a shared library on non-Windows systems, you could try+//+// g++ -fPIC -shared -DNDEBUG -Ofast erf_cody.cpp rationalcubic.cpp normaldistribution.cpp lets_be_rational.cpp -o lets_be_rational.so+//++#if defined(_MSC_VER)+# define NOMINMAX // to suppress MSVC's definitions of min() and max()+// These four pragmas are the equivalent to /fp:fast.+# pragma float_control( except, off )+# pragma float_control( precise, off )+# pragma fp_contract( on )+# pragma fenv_access( off )+#endif++#include "normaldistribution.h"+#include "rationalcubic.h"+#include <float.h>+#include <cmath>+#include <algorithm>+#if defined(_WIN32) || defined(_WIN64)+# include <windows.h>+#endif++#define TWO_PI 6.283185307179586476925286766559005768394338798750+#define SQRT_PI_OVER_TWO 1.253314137315500251207882642405522626503493370305 // sqrt(pi/2) to avoid misinterpretation.+#define SQRT_THREE 1.732050807568877293527446341505872366942805253810+#define SQRT_ONE_OVER_THREE 0.577350269189625764509148780501957455647601751270+#define TWO_PI_OVER_SQRT_TWENTY_SEVEN 1.209199576156145233729385505094770488189377498728 // 2*pi/sqrt(27)+#define PI_OVER_SIX 0.523598775598298873077107230546583814032861566563++namespace {+ static const double SQRT_DBL_EPSILON = sqrt(DBL_EPSILON);+ static const double FOURTH_ROOT_DBL_EPSILON = sqrt(SQRT_DBL_EPSILON);+ static const double EIGHTH_ROOT_DBL_EPSILON = sqrt(FOURTH_ROOT_DBL_EPSILON);+ static const double SIXTEENTH_ROOT_DBL_EPSILON = sqrt(EIGHTH_ROOT_DBL_EPSILON);+ static const double SQRT_DBL_MIN = sqrt(DBL_MIN);+ static const double SQRT_DBL_MAX = sqrt(DBL_MAX);++ // Set this to 0 if you want positive results for (positive) denormalised inputs, else to DBL_MIN.+ // Note that you cannot achieve full machine accuracy from denormalised inputs!+ static const double DENORMALISATION_CUTOFF = 0; ++ static const double VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC = -DBL_MAX;+ static const double VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM = DBL_MAX;++ inline bool is_below_horizon(double x){ return fabs(x) < DENORMALISATION_CUTOFF; } // This weeds out denormalised (a.k.a. 'subnormal') numbers.++ // See https://www.kernel.org/doc/Documentation/atomic_ops.txt for further details on this simplistic implementation of an atomic flag that is *not* volatile.+ typedef struct { +#if defined(_MSC_VER) || defined(_WIN32) || defined(_WIN64)+ long data;+#else+ int data;+#endif+ } atomic_t;++ static atomic_t implied_volatility_maximum_iterations = { 2 }; // (DBL_DIG*20)/3 ≈ 100 . Only needed when the iteration effectively alternates Householder/Halley/Newton steps and binary nesting due to roundoff truncation.++#ifdef ENABLE_SWITCHING_THE_OUTPUT_TO_ITERATION_COUNT+ static atomic_t implied_volatility_output_type = { 0 };+ inline double implied_volatility_output(int count, double volatility){ return implied_volatility_output_type.data>0 ? count : volatility; }+#else+ inline double implied_volatility_output(int count, double volatility){ return volatility; }+#endif++#ifdef ENABLE_CHANGING_THE_HOUSEHOLDER_METHOD_ORDER+ static atomic_t implied_volatility_householder_method_order = { 4 };+ inline double householder_factor(double newton, double halley, double hh3){+ return implied_volatility_householder_method_order.data > 3 ? (1+0.5*halley*newton)/(1+newton*(halley+hh3*newton/6)) : ( implied_volatility_householder_method_order.data > 2 ? 1/(1+0.5*halley*newton) : 1 );+ }+#else+ inline double householder_factor(double newton, double halley, double hh3){ return (1+0.5*halley*newton)/(1+newton*(halley+hh3*newton/6)); }+#endif++}++EXPORT_EXTERN_C double set_implied_volatility_maximum_iterations(double t){+ int i = (int)t;+ if (i>=0) {+#if defined(_MSC_VER) || defined(_WIN32) || defined(_WIN64)+ InterlockedExchange(&(implied_volatility_maximum_iterations.data),i);+#elif defined( __x86__ ) || defined( __x86_64__ )+ implied_volatility_maximum_iterations.data = i;+#else+# error Atomic operations not implemented for this platform.+#endif+ }+ return implied_volatility_maximum_iterations.data;+}++#ifdef ENABLE_SWITCHING_THE_OUTPUT_TO_ITERATION_COUNT+EXPORT_EXTERN_C double set_implied_volatility_output_type(double t){+ int i = (int)t;+#if defined(_MSC_VER) || defined(_WIN32) || defined(_WIN64)+ InterlockedExchange(&(implied_volatility_output_type.data),i);+#elif defined( __x86__ ) || defined( __x86_64__ )+ implied_volatility_output_type.data = i;+#else+# error Atomic operations not implemented for this platform.+#endif+ return implied_volatility_output_type.data;+}+#endif ++#ifdef ENABLE_CHANGING_THE_HOUSEHOLDER_METHOD_ORDER+EXPORT_EXTERN_C double set_implied_volatility_householder_method_order(double t){+ int i = (int)t;+ if (i>=0) {+#if defined(_MSC_VER) || defined(_WIN32) || defined(_WIN64)+ InterlockedExchange(&(implied_volatility_householder_method_order.data),i);+#elif defined( __x86__ ) || defined( __x86_64__ )+ implied_volatility_householder_method_order.data = i;+#else+# error Atomic operations not implemented for this platform.+#endif+ }+ return implied_volatility_householder_method_order.data;+}+#endif ++double normalised_intrinsic(double x, double q /* q=±1 */){+ if (q*x<=0)+ return 0;+ const double x2=x*x;+ if (x2<98*FOURTH_ROOT_DBL_EPSILON ) // The factor 98 is computed from last coefficient: √√92897280 = 98.1749+ return fabs( std::max( (q<0?-1:1)*x*(1+x2*((1.0/24.0)+x2*((1.0/1920.0)+x2*((1.0/322560.0)+(1.0/92897280.0)*x2)))) , 0.0 ) );+ const double b_max = exp(0.5*x), one_over_b_max = 1 / b_max;+ return fabs(std::max((q<0?-1:1)*(b_max-one_over_b_max),0.));+}++double normalised_intrinsic_call(double x){ return normalised_intrinsic(x,1); }++// Asymptotic expansion of+//+// b = Φ(h+t)·exp(x/2) - Φ(h-t)·exp(-x/2)+// with+// h = x/s and t = s/2+// which makes+// b = Φ(h+t)·exp(h·t) - Φ(h-t)·exp(-h·t)+//+// exp(-(h²+t²)/2)+// = --------------- · [ Y(h+t) - Y(h-t) ]+// √(2π)+// with+// Y(z) := Φ(z)/φ(z)+//+// for large negative (t-|h|) by the aid of Abramowitz & Stegun (26.2.12) where Φ(z) = φ(z)/|z|·[1-1/z^2+...].+// We define+// r+// A(h,t) := --- · [ Y(h+t) - Y(h-t) ]+// t+//+// with r := (h+t)·(h-t) and give an expansion for A(h,t) in q:=(h/r)² expressed in terms of e:=(t/h)² .+double asymptotic_expansion_of_normalised_black_call(double h, double t){+ const double e=(t/h)*(t/h), r=((h+t)*(h-t)), q=(h/r)*(h/r);+ // 17th order asymptotic expansion of A(h,t) in q, sufficient for Φ(h) [and thus y(h)] to have relative accuracy of 1.64E-16 for h <= η with η:=-10.+ const double asymptotic_expansion_sum = (2.0+q*(-6.0E0-2.0*e+3.0*q*(1.0E1+e*(2.0E1+2.0*e)+5.0*q*(-1.4E1+e*(-7.0E1+e*(-4.2E1-2.0*e))+7.0*q*(1.8E1+e*(1.68E2+e*(2.52E2+e*(7.2E1+2.0*e)))+9.0*q*(-2.2E1+e*(-3.3E2+e*(-9.24E2+e*(-6.6E2+e*(-1.1E2-2.0*e))))+1.1E1*q*(2.6E1+e*(5.72E2+e*(2.574E3+e*(3.432E3+e*(1.43E3+e*(1.56E2+2.0*e)))))+1.3E1*q*(-3.0E1+e*(-9.1E2+e*(-6.006E3+e*(-1.287E4+e*(-1.001E4+e*(-2.73E3+e*(-2.1E2-2.0*e))))))+1.5E1*q*(3.4E1+e*(1.36E3+e*(1.2376E4+e*(3.8896E4+e*(4.862E4+e*(2.4752E4+e*(4.76E3+e*(2.72E2+2.0*e)))))))+1.7E1*q*(-3.8E1+e*(-1.938E3+e*(-2.3256E4+e*(-1.00776E5+e*(-1.84756E5+e*(-1.51164E5+e*(-5.4264E4+e*(-7.752E3+e*(-3.42E2-2.0*e))))))))+1.9E1*q*(4.2E1+e*(2.66E3+e*(4.0698E4+e*(2.3256E5+e*(5.8786E5+e*(7.05432E5+e*(4.0698E5+e*(1.08528E5+e*(1.197E4+e*(4.2E2+2.0*e)))))))))+2.1E1*q*(-4.6E1+e*(-3.542E3+e*(-6.7298E4+e*(-4.90314E5+e*(-1.63438E6+e*(-2.704156E6+e*(-2.288132E6+e*(-9.80628E5+e*(-2.01894E5+e*(-1.771E4+e*(-5.06E2-2.0*e))))))))))+2.3E1*q*(5.0E1+e*(4.6E3+e*(1.0626E5+e*(9.614E5+e*(4.08595E6+e*(8.9148E6+e*(1.04006E7+e*(6.53752E6+e*(2.16315E6+e*(3.542E5+e*(2.53E4+e*(6.0E2+2.0*e)))))))))))+2.5E1*q*(-5.4E1+e*(-5.85E3+e*(-1.6146E5+e*(-1.77606E6+e*(-9.37365E6+e*(-2.607579E7+e*(-4.01166E7+e*(-3.476772E7+e*(-1.687257E7+e*(-4.44015E6+e*(-5.9202E5+e*(-3.51E4+e*(-7.02E2-2.0*e))))))))))))+2.7E1*q*(5.8E1+e*(7.308E3+e*(2.3751E5+e*(3.12156E6+e*(2.003001E7+e*(6.919458E7+e*(1.3572783E8+e*(1.5511752E8+e*(1.0379187E8+e*(4.006002E7+e*(8.58429E6+e*(9.5004E5+e*(4.7502E4+e*(8.12E2+2.0*e)))))))))))))+2.9E1*q*(-6.2E1+e*(-8.99E3+e*(-3.39822E5+e*(-5.25915E6+e*(-4.032015E7+e*(-1.6934463E8+e*(-4.1250615E8+e*(-6.0108039E8+e*(-5.3036505E8+e*(-2.8224105E8+e*(-8.870433E7+e*(-1.577745E7+e*(-1.472562E6+e*(-6.293E4+e*(-9.3E2-2.0*e))))))))))))))+3.1E1*q*(6.6E1+e*(1.0912E4+e*(4.74672E5+e*(8.544096E6+e*(7.71342E7+e*(3.8707344E8+e*(1.14633288E9+e*(2.07431664E9+e*(2.33360622E9+e*(1.6376184E9+e*(7.0963464E8+e*(1.8512208E8+e*(2.7768312E7+e*(2.215136E6+e*(8.184E4+e*(1.056E3+2.0*e)))))))))))))))+3.3E1*(-7.0E1+e*(-1.309E4+e*(-6.49264E5+e*(-1.344904E7+e*(-1.4121492E8+e*(-8.344518E8+e*(-2.9526756E9+e*(-6.49588632E9+e*(-9.0751353E9+e*(-8.1198579E9+e*(-4.6399188E9+e*(-1.6689036E9+e*(-3.67158792E8+e*(-4.707164E7+e*(-3.24632E6+e*(-1.0472E5+e*(-1.19E3-2.0*e)))))))))))))))))*q)))))))))))))))));+ const double b = ONE_OVER_SQRT_TWO_PI*exp((-0.5*(h*h+t*t)))*(t/r)*asymptotic_expansion_sum;+ return fabs(std::max(b , 0.));+}++namespace { /* η */ static const double asymptotic_expansion_accuracy_threshold = -10; }++double normalised_black_call_using_erfcx(double h, double t) {+ // Given h = x/s and t = s/2, the normalised Black function can be written as+ //+ // b(x,s) = Φ(x/s+s/2)·exp(x/2) - Φ(x/s-s/2)·exp(-x/2)+ // = Φ(h+t)·exp(h·t) - Φ(h-t)·exp(-h·t) . (*)+ //+ // It is mentioned in section 4 (and discussion of figures 2 and 3) of George Marsaglia's article "Evaluating the+ // Normal Distribution" (available at http://www.jstatsoft.org/v11/a05/paper) that the error of any cumulative normal+ // function Φ(z) is dominated by the hardware (or compiler implementation) accuracy of exp(-z²/2) which is not+ // reliably more than 14 digits when z is large. The accuracy of Φ(z) typically starts coming down to 14 digits when+ // z is around -8. For the (normalised) Black function, as above in (*), this means that we are subtracting two terms+ // that are each products of terms with about 14 digits of accuracy. The net result, in each of the products, is even+ // less accuracy, and then we are taking the difference of these terms, resulting in even less accuracy. When we are+ // using the asymptotic expansion asymptotic_expansion_of_normalised_black_call() invoked in the second branch at the+ // beginning of this function, we are using only *one* exponential instead of 4, and this improves accuracy. It+ // actually improves it a bit more than you would expect from the above logic, namely, almost the full two missing+ // digits (in 64 bit IEEE floating point). Unfortunately, going higher order in the asymptotic expansion will not+ // enable us to gain more accuracy (by extending the range in which we could use the expansion) since the asymptotic+ // expansion, being a divergent series, can never gain 16 digits of accuracy for z=-8 or just below. The best you can+ // get is about 15 digits (just), for about 35 terms in the series (26.2.12), which would result in an prohibitively+ // long expression in function asymptotic expansion asymptotic_expansion_of_normalised_black_call(). In this last branch,+ // here, we therefore take a different tack as follows.+ // The "scaled complementary error function" is defined as erfcx(z) = exp(z²)·erfc(z). Cody's implementation of this+ // function as published in "Rational Chebyshev approximations for the error function", W. J. Cody, Math. Comp., 1969, pp.+ // 631-638, uses rational functions that theoretically approximates erfcx(x) to at least 18 significant decimal digits,+ // *without* the use of the exponential function when x>4, which translates to about z<-5.66 in Φ(z). To make use of it,+ // we write+ // Φ(z) = exp(-z²/2)·erfcx(-z/√2)/2+ //+ // to transform the normalised black function to+ //+ // b = ½ · exp(-½(h²+t²)) · [ erfcx(-(h+t)/√2) - erfcx(-(h-t)/√2) ]+ //+ // which now involves only one exponential, instead of three, when |h|+|t| > 5.66 , and the difference inside the+ // square bracket is between the evaluation of two rational functions, which, typically, according to Marsaglia,+ // retains the full 16 digits of accuracy (or just a little less than that).+ //+ const double b = 0.5 * exp(-0.5*(h*h+t*t)) * ( erfcx_cody(-ONE_OVER_SQRT_TWO*(h+t)) - erfcx_cody(-ONE_OVER_SQRT_TWO*(h-t)) );+ return fabs(std::max(b,0.0));+}++// Calculation of+//+// b = Φ(h+t)·exp(h·t) - Φ(h-t)·exp(-h·t)+//+// exp(-(h²+t²)/2)+// = --------------- · [ Y(h+t) - Y(h-t) ]+// √(2π)+// with+// Y(z) := Φ(z)/φ(z)+//+// using an expansion of Y(h+t)-Y(h-t) for small t to twelvth order in t.+// Theoretically accurate to (better than) precision ε = 2.23E-16 when h<=0 and t < τ with τ := 2·ε^(1/16) ≈ 0.21.+// The main bottleneck for precision is the coefficient a:=1+h·Y(h) when |h|>1 .+double small_t_expansion_of_normalised_black_call(double h, double t){+ // Y(h) := Φ(h)/φ(h) = √(π/2)·erfcx(-h/√2)+ // a := 1+h·Y(h) --- Note that due to h<0, and h·Y(h) -> -1 (from above) as h -> -∞, we also have that a>0 and a -> 0 as h -> -∞+ // w := t² , h2 := h²+ const double a = 1+h*(0.5*SQRT_TWO_PI)*erfcx_cody(-ONE_OVER_SQRT_TWO*h), w=t*t, h2=h*h;+ const double expansion = 2*t*(a+w*((-1+3*a+a*h2)/6+w*((-7+15*a+h2*(-1+10*a+a*h2))/120+w*((-57+105*a+h2*(-18+105*a+h2*(-1+21*a+a*h2)))/5040+w*((-561+945*a+h2*(-285+1260*a+h2*(-33+378*a+h2*(-1+36*a+a*h2))))/362880+w*((-6555+10395*a+h2*(-4680+17325*a+h2*(-840+6930*a+h2*(-52+990*a+h2*(-1+55*a+a*h2)))))/39916800+((-89055+135135*a+h2*(-82845+270270*a+h2*(-20370+135135*a+h2*(-1926+25740*a+h2*(-75+2145*a+h2*(-1+78*a+a*h2))))))*w)/6227020800.0))))));+ const double b = ONE_OVER_SQRT_TWO_PI*exp((-0.5*(h*h+t*t)))*expansion;+ return fabs(std::max(b,0.0));+}++namespace { /* τ */ static const double small_t_expansion_of_normalised_black_threshold = 2*SIXTEENTH_ROOT_DBL_EPSILON; }++// b(x,s) = Φ(x/s+s/2)·exp(x/2) - Φ(x/s-s/2)·exp(-x/2)+// = Φ(h+t)·exp(x/2) - Φ(h-t)·exp(-x/2)+// with+// h = x/s and t = s/2+double normalised_black_call_using_norm_cdf(double x, double s){+ const double h = x/s, t = 0.5*s, b_max = exp(0.5*x), b = norm_cdf(h + t) * b_max - norm_cdf(h - t) / b_max;+ return fabs(std::max(b,0.0));+}++//+// Introduced on 2017-02-18+//+// b(x,s) = Φ(x/s+s/2)·exp(x/2) - Φ(x/s-s/2)·exp(-x/2)+// = Φ(h+t)·exp(x/2) - Φ(h-t)·exp(-x/2)+// = ½ · exp(-u²-v²) · [ erfcx(u-v) - erfcx(u+v) ]+// = ½ · [ exp(x/2)·erfc(u-v) - exp(-x/2)·erfc(u+v) ]+// = ½ · [ exp(x/2)·erfc(u-v) - exp(-u²-v²)·erfcx(u+v) ]+// = ½ · [ exp(-u²-v²)·erfcx(u-v) - exp(-x/2)·erfc(u+v) ]+// with+// h = x/s , t = s/2 ,+// and+// u = -h/√2 and v = t/√2 .+//+// Cody's erfc() and erfcx() functions each, for some values of their argument, involve the evaluation+// of the exponential function exp(). The normalised Black function requires additional evaluation(s)+// of the exponential function irrespective of which of the above formulations is used. However, the total+// number of exponential function evaluations can be minimised by a judicious choice of one of the above+// formulations depending on the input values and the branch logic in Cody's erfc() and erfcx().+//+double normalised_black_call_with_optimal_use_of_codys_functions(double x, double s){+ const double codys_threshold = 0.46875, h = x/s, t = 0.5*s, q1 = -ONE_OVER_SQRT_TWO*(h+t), q2 = -ONE_OVER_SQRT_TWO*(h-t);+ double two_b;+ if ( q1 < codys_threshold )+ if ( q2 < codys_threshold )+ two_b = exp(0.5*x)*erfc_cody(q1) - exp(-0.5*x)*erfc_cody(q2);+ else+ two_b = exp(0.5*x)*erfc_cody(q1) - exp(-0.5*(h*h+t*t))*erfcx_cody(q2);+ else+ if ( q2 < codys_threshold )+ two_b = exp(-0.5*(h*h+t*t))*erfcx_cody(q1) - exp(-0.5*x)*erfc_cody(q2);+ else+ two_b = exp(-0.5*(h*h+t*t)) * ( erfcx_cody(q1) - erfcx_cody(q2) );+ return fabs(std::max(0.5*two_b,0.0));+}++EXPORT_EXTERN_C double normalised_black_call(double x, double s) {+ if (x>0)+ return normalised_intrinsic_call(x)+normalised_black_call(-x,s); // In the money.+ if (s<=fabs(x)*DENORMALISATION_CUTOFF)+ return normalised_intrinsic_call(x); // sigma=0 -> intrinsic value.+ // Denote h := x/s and t := s/2.+ // We evaluate the condition |h|>|η|, i.e., h<η && t < τ+|h|-|η| avoiding any divisions by s , where η = asymptotic_expansion_accuracy_threshold and τ = small_t_expansion_of_normalised_black_threshold .+ if ( x < s*asymptotic_expansion_accuracy_threshold && 0.5*s*s+x < s*(small_t_expansion_of_normalised_black_threshold+asymptotic_expansion_accuracy_threshold) )+ return asymptotic_expansion_of_normalised_black_call(x/s,0.5*s);+ if ( 0.5*s < small_t_expansion_of_normalised_black_threshold )+ return small_t_expansion_of_normalised_black_call(x/s,0.5*s);+#ifdef DO_NOT_OPTIMISE_NORMALISED_BLACK_IN_REGIONS_3_AND_4_FOR_CODYS_FUNCTIONS+ // When b is more than, say, about 85% of b_max=exp(x/2), then b is dominated by the first of the two terms in the Black formula, and we retain more accuracy by not attempting to combine the two terms in any way.+ // We evaluate the condition h+t>0.85 avoiding any divisions by s.+ if ( x+0.5*s*s > s*0.85 )+ return normalised_black_call_using_norm_cdf(x,s);+ return normalised_black_call_using_erfcx(x/s,0.5*s);+#else+ return normalised_black_call_with_optimal_use_of_codys_functions(x,s);+#endif+}++inline double square(double x){ return x*x; }++EXPORT_EXTERN_C double normalised_vega(double x, double s) {+ const double ax = fabs(x);+ return (ax<=0) ? ONE_OVER_SQRT_TWO_PI*exp(-0.125*s*s) : ( (s<=0 || s<=ax*SQRT_DBL_MIN) ? 0 : ONE_OVER_SQRT_TWO_PI*exp(-0.5*(square(x/s)+square(0.5*s))) );+}++EXPORT_EXTERN_C double normalised_black(double x, double s, double q /* q=±1 */) { return normalised_black_call(q<0?-x:x,s); /* Reciprocal-strike call-put equivalence */ }++EXPORT_EXTERN_C double black(double F, double K, double sigma, double T, double q /* q=±1 */) {+ const double intrinsic = fabs(std::max((q<0?K-F:F-K),0.0));+ // Map in-the-money to out-of-the-money+ if (q*(F-K)>0)+ return intrinsic + black(F,K,sigma,T,-q);+ return std::max(intrinsic,(sqrt(F)*sqrt(K))*normalised_black(log(F/K),sigma*sqrt(T),q));+}++#ifdef COMPUTE_LOWER_MAP_DERIVATIVES_INDIVIDUALLY+double f_lower_map(const double x,const double s){ + if (is_below_horizon(x))+ return 0;+ if (is_below_horizon(s))+ return 0;+ const double z=SQRT_ONE_OVER_THREE*fabs(x)/s, Phi=norm_cdf(-z);+ return TWO_PI_OVER_SQRT_TWENTY_SEVEN*fabs(x)*(Phi*Phi*Phi);+}+double d_f_lower_map_d_beta(const double x,const double s){+ if (is_below_horizon(s))+ return 1;+ const double z=SQRT_ONE_OVER_THREE*fabs(x)/s, y = z*z, Phi=norm_cdf(-z);+ return TWO_PI*y*(Phi*Phi) * exp(y+0.125*s*s);+}+double d2_f_lower_map_d_beta2(const double x,const double s){+ const double ax=fabs(x), z=SQRT_ONE_OVER_THREE*ax/s, y = z*z, s2=s*s, Phi=norm_cdf(-z), phi=norm_pdf(z);+ return PI_OVER_SIX * y/(s2*s) * Phi * ( 8*SQRT_THREE*s*ax + (3*s2*(s2-8)-8*x*x)*Phi/phi ) * exp(2*y+0.25*s2);+}+void compute_f_lower_map_and_first_two_derivatives(const double x,const double s,double &f,double &fp,double &fpp){+ f = f_lower_map(x,s);+ fp = d_f_lower_map_d_beta(x,s);+ fpp = d2_f_lower_map_d_beta2(x,s);+}+#else+void compute_f_lower_map_and_first_two_derivatives(const double x,const double s,double &f,double &fp,double &fpp){+ const double ax=fabs(x), z=SQRT_ONE_OVER_THREE*ax/s, y = z*z, s2=s*s, Phi=norm_cdf(-z), phi=norm_pdf(z);+ fpp = PI_OVER_SIX * y/(s2*s) * Phi * ( 8*SQRT_THREE*s*ax + (3*s2*(s2-8)-8*x*x)*Phi/phi ) * exp(2*y+0.25*s2);+ if (is_below_horizon(s)) {+ fp = 1;+ f = 0;+ } else {+ const double Phi2=Phi*Phi;+ fp = TWO_PI*y*Phi2*exp(y+0.125*s*s);+ if (is_below_horizon(x))+ f = 0;+ else+ f = TWO_PI_OVER_SQRT_TWENTY_SEVEN*ax*(Phi2*Phi);+ }+}+#endif++double inverse_f_lower_map(const double x,const double f){+ return is_below_horizon(f) ? 0 : fabs(x/(SQRT_THREE*inverse_norm_cdf( std::pow( f/(TWO_PI_OVER_SQRT_TWENTY_SEVEN*fabs(x)) , 1./3.) ))); +}++#ifdef COMPUTE_UPPER_MAP_DERIVATIVES_INDIVIDUALLY+double f_upper_map(const double s){+ return norm_cdf(-0.5*s);+}+double d_f_upper_map_d_beta(const double x,const double s){+ return is_below_horizon(x) ? -0.5 : -0.5*exp(0.5*square(x/s));+}+double d2_f_upper_map_d_beta2(const double x,const double s){+ if (is_below_horizon(x))+ return 0;+ const double w = square(x/s);+ return SQRT_PI_OVER_TWO*exp(w+0.125*s*s)*w/s;+}+void compute_f_upper_map_and_first_two_derivatives(const double x,const double s,double &f,double &fp,double &fpp){+ f = f_upper_map(s);+ fp = d_f_upper_map_d_beta(x,s);+ fpp = d2_f_upper_map_d_beta2(x,s);+}+#else+void compute_f_upper_map_and_first_two_derivatives(const double x,const double s,double &f,double &fp,double &fpp){+ f = norm_cdf(-0.5*s);+ if (is_below_horizon(x)) {+ fp = -0.5;+ fpp = 0;+ } else {+ const double w = square(x/s);+ fp = -0.5*exp(0.5*w);+ fpp = SQRT_PI_OVER_TWO*exp(w+0.125*s*s)*w/s;+ }+}+#endif++double inverse_f_upper_map(double f){+ return -2.*inverse_norm_cdf(f);+}++// See http://en.wikipedia.org/wiki/Householder%27s_method for a detailed explanation of the third order Householder iteration.+//+// Given the objective function g(s) whose root x such that 0 = g(s) we seek, iterate+//+// s_n+1 = s_n - (g/g') · [ 1 - (g''/g')·(g/g') ] / [ 1 - (g/g')·( (g''/g') - (g'''/g')·(g/g')/6 ) ]+//+// Denoting newton:=-(g/g'), halley:=(g''/g'), and hh3:=(g'''/g'), this reads+//+// s_n+1 = s_n + newton · [ 1 + halley·newton/2 ] / [ 1 + newton·( halley + hh3·newton/6 ) ]+//+//+// NOTE that this function returns 0 when beta<intrinsic without any safety checks.+//+double unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_with_limited_iterations(double beta, double x, double q /* q=±1 */, int N){+ // Subtract intrinsic.+ if (q*x>0) {+ beta = fabs(std::max(beta-normalised_intrinsic(x, q),0.));+ q = -q;+ }+ // Map puts to calls+ if (q<0){+ x = -x;+ q = -q;+ }+ if (beta<=0) // For negative or zero prices we return 0.+ return implied_volatility_output(0,0);+ if (beta<DENORMALISATION_CUTOFF) // For positive but denormalised (a.k.a. 'subnormal') prices, we return 0 since it would be impossible to converge to full machine accuracy anyway.+ return implied_volatility_output(0,0);+ const double b_max = exp(0.5*x);+ if (beta>=b_max)+ return implied_volatility_output(0,VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM);+ int iterations=0, direction_reversal_count = 0;+ double f=-DBL_MAX, s=-DBL_MAX, ds=s, ds_previous=0, s_left=DBL_MIN, s_right=DBL_MAX;+ // The temptation is great to use the optimised form b_c = exp(x/2)/2-exp(-x/2)·Phi(sqrt(-2·x)) but that would require implementing all of the above types of round-off and over/underflow handling for this expression, too.+ const double s_c=sqrt(fabs(2*x)), b_c = normalised_black_call(x,s_c), v_c = normalised_vega(x, s_c);+ // Four branches.+ if ( beta<b_c ) {+ const double s_l = s_c - b_c/v_c, b_l = normalised_black_call(x,s_l);+ if (beta<b_l){+ double f_lower_map_l, d_f_lower_map_l_d_beta, d2_f_lower_map_l_d_beta2;+ compute_f_lower_map_and_first_two_derivatives(x,s_l,f_lower_map_l,d_f_lower_map_l_d_beta,d2_f_lower_map_l_d_beta2);+ const double r_ll=convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(0.,b_l,0.,f_lower_map_l,1.,d_f_lower_map_l_d_beta,d2_f_lower_map_l_d_beta2,true);+ f = rational_cubic_interpolation(beta,0.,b_l,0.,f_lower_map_l,1.,d_f_lower_map_l_d_beta,r_ll);+ if (!(f>0)) { // This can happen due to roundoff truncation for extreme values such as |x|>500.+ // We switch to quadratic interpolation using f(0)≡0, f(b_l), and f'(0)≡1 to specify the quadratic.+ const double t = beta/b_l;+ f = (f_lower_map_l*t + b_l*(1-t)) * t;+ }+ s = inverse_f_lower_map(x,f);+ s_right = s_l;+ //+ // In this branch, which comprises the lowest segment, the objective function is+ // g(s) = 1/ln(b(x,s)) - 1/ln(beta)+ // ≡ 1/ln(b(s)) - 1/ln(beta)+ // This makes+ // g' = -b'/(b·ln(b)²)+ // newton = -g/g' = (ln(beta)-ln(b))·ln(b)/ln(beta)·b/b'+ // halley = g''/g' = b''/b' - b'/b·(1+2/ln(b))+ // hh3 = g'''/g' = b'''/b' + 2(b'/b)²·(1+3/ln(b)·(1+1/ln(b))) - 3(b''/b)·(1+2/ln(b))+ //+ // The Householder(3) iteration is+ // s_n+1 = s_n + newton · [ 1 + halley·newton/2 ] / [ 1 + newton·( halley + hh3·newton/6 ) ]+ //+ for (; iterations<N && fabs(ds)>DBL_EPSILON*s; ++iterations){+ if (ds*ds_previous<0)+ ++direction_reversal_count;+ if ( iterations>0 && ( 3==direction_reversal_count || !(s>s_left && s<s_right) ) ) {+ // If looping inefficently, or the forecast step takes us outside the bracket, or onto its edges, switch to binary nesting.+ // NOTE that this can only really happen for very extreme values of |x|, such as |x| = |ln(F/K)| > 500.+ s = 0.5*(s_left+s_right);+ if (s_right-s_left<=DBL_EPSILON*s) break;+ direction_reversal_count = 0;+ ds = 0;+ }+ ds_previous=ds;+ const double b = normalised_black_call(x,s), bp = normalised_vega(x, s);+ if ( b>beta && s<s_right ) s_right=s; else if ( b<beta && s>s_left ) s_left=s; // Tighten the bracket if applicable.+ if (b<=0||bp<=0) // Numerical underflow. Switch to binary nesting for this iteration.+ ds = 0.5*(s_left+s_right)-s;+ else {+ const double ln_b=log(b), ln_beta=log(beta), bpob=bp/b, h=x/s, b_halley = h*h/s-s/4, newton = (ln_beta-ln_b)*ln_b/ln_beta/bpob, halley = b_halley-bpob*(1+2/ln_b);+ const double b_hh3 = b_halley*b_halley-3*square(h/s)-0.25, hh3 = b_hh3+2*square(bpob)*(1+3/ln_b*(1+1/ln_b))-3*b_halley*bpob*(1+2/ln_b);+ ds = newton * householder_factor(newton,halley,hh3);+ }+ s += ds = std::max(-0.5*s , ds );+ }+ return implied_volatility_output(iterations,s);+ } else {+ const double v_l = normalised_vega(x, s_l), r_lm = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(b_l,b_c,s_l,s_c,1/v_l,1/v_c,0.0,false);+ s = rational_cubic_interpolation(beta,b_l,b_c,s_l,s_c,1/v_l,1/v_c,r_lm);+ s_left = s_l;+ s_right = s_c;+ }+ } else {+ const double s_h = v_c>DBL_MIN ? s_c+(b_max-b_c)/v_c : s_c, b_h = normalised_black_call(x,s_h);+ if(beta<=b_h){+ const double v_h = normalised_vega(x, s_h), r_hm = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(b_c,b_h,s_c,s_h,1/v_c,1/v_h,0.0,false);+ s = rational_cubic_interpolation(beta,b_c,b_h,s_c,s_h,1/v_c,1/v_h,r_hm);+ s_left = s_c;+ s_right = s_h;+ } else {+ double f_upper_map_h, d_f_upper_map_h_d_beta, d2_f_upper_map_h_d_beta2;+ compute_f_upper_map_and_first_two_derivatives(x,s_h,f_upper_map_h,d_f_upper_map_h_d_beta,d2_f_upper_map_h_d_beta2);+ if ( d2_f_upper_map_h_d_beta2>-SQRT_DBL_MAX && d2_f_upper_map_h_d_beta2<SQRT_DBL_MAX ){+ const double r_hh = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(b_h,b_max,f_upper_map_h,0.,d_f_upper_map_h_d_beta,-0.5,d2_f_upper_map_h_d_beta2,true);+ f = rational_cubic_interpolation(beta,b_h,b_max,f_upper_map_h,0.,d_f_upper_map_h_d_beta,-0.5,r_hh);+ }+ if (f<=0) {+ const double h=b_max-b_h, t=(beta-b_h)/h;+ f = (f_upper_map_h*(1-t) + 0.5*h*t) * (1-t); // We switch to quadratic interpolation using f(b_h), f(b_max)≡0, and f'(b_max)≡-1/2 to specify the quadratic.+ }+ s = inverse_f_upper_map(f);+ s_left = s_h;+ if (beta>0.5*b_max) { // Else we better drop through and let the objective function be g(s) = b(x,s)-beta. + //+ // In this branch, which comprises the upper segment, the objective function is+ // g(s) = ln(b_max-beta)-ln(b_max-b(x,s))+ // ≡ ln((b_max-beta)/(b_max-b(s)))+ // This makes+ // g' = b'/(b_max-b)+ // newton = -g/g' = ln((b_max-b)/(b_max-beta))·(b_max-b)/b'+ // halley = g''/g' = b''/b' + b'/(b_max-b)+ // hh3 = g'''/g' = b'''/b' + g'·(2g'+3b''/b')+ // and the iteration is+ // s_n+1 = s_n + newton · [ 1 + halley·newton/2 ] / [ 1 + newton·( halley + hh3·newton/6 ) ].+ //+ for (; iterations<N && fabs(ds)>DBL_EPSILON*s; ++iterations){+ if (ds*ds_previous<0)+ ++direction_reversal_count;+ if ( iterations>0 && ( 3==direction_reversal_count || !(s>s_left && s<s_right) ) ) {+ // If looping inefficently, or the forecast step takes us outside the bracket, or onto its edges, switch to binary nesting.+ // NOTE that this can only really happen for very extreme values of |x|, such as |x| = |ln(F/K)| > 500.+ s = 0.5*(s_left+s_right);+ if (s_right-s_left<=DBL_EPSILON*s) break;+ direction_reversal_count = 0;+ ds = 0;+ }+ ds_previous=ds;+ const double b = normalised_black_call(x,s), bp = normalised_vega(x, s);+ if ( b>beta && s<s_right ) s_right=s; else if ( b<beta && s>s_left ) s_left=s; // Tighten the bracket if applicable.+ if (b>=b_max||bp<=DBL_MIN) // Numerical underflow. Switch to binary nesting for this iteration.+ ds = 0.5*(s_left+s_right)-s;+ else {+ const double b_max_minus_b = b_max-b, g = log((b_max-beta)/b_max_minus_b), gp = bp/b_max_minus_b;+ const double b_halley = square(x/s)/s-s/4, b_hh3 = b_halley*b_halley-3*square(x/(s*s))-0.25;+ const double newton = -g/gp, halley = b_halley+gp, hh3 = b_hh3+gp*(2*gp+3*b_halley);+ ds = newton * householder_factor(newton,halley,hh3);+ }+ s += ds = std::max(-0.5*s , ds );+ }+ return implied_volatility_output(iterations,s);+ }+ }+ }+ // In this branch, which comprises the two middle segments, the objective function is g(s) = b(x,s)-beta, or g(s) = b(s) - beta, for short.+ // This makes+ // newton = -g/g' = -(b-beta)/b'+ // halley = g''/g' = b''/b' = x²/s³-s/4+ // hh3 = g'''/g' = b'''/b' = halley² - 3·(x/s²)² - 1/4+ // and the iteration is+ // s_n+1 = s_n + newton · [ 1 + halley·newton/2 ] / [ 1 + newton·( halley + hh3·newton/6 ) ].+ //+ for (; iterations<N && fabs(ds)>DBL_EPSILON*s; ++iterations){+ if (ds*ds_previous<0)+ ++direction_reversal_count;+ if ( iterations>0 && ( 3==direction_reversal_count || !(s>s_left && s<s_right) ) ) {+ // If looping inefficently, or the forecast step takes us outside the bracket, or onto its edges, switch to binary nesting.+ // NOTE that this can only really happen for very extreme values of |x|, such as |x| = |ln(F/K)| > 500.+ s = 0.5*(s_left+s_right);+ if (s_right-s_left<=DBL_EPSILON*s) break;+ direction_reversal_count = 0;+ ds = 0;+ }+ ds_previous=ds;+ const double b = normalised_black_call(x,s), bp = normalised_vega(x, s);+ if ( b>beta && s<s_right ) s_right=s; else if ( b<beta && s>s_left ) s_left=s; // Tighten the bracket if applicable.+ const double newton = (beta-b)/bp, halley = square(x/s)/s-s/4, hh3 = halley*halley-3*square(x/(s*s))-0.25;+ s += ds = std::max(-0.5*s , newton * householder_factor(newton,halley,hh3) );+ }+ return implied_volatility_output(iterations,s);+}++EXPORT_EXTERN_C double implied_volatility_from_a_transformed_rational_guess_with_limited_iterations(double price, double F, double K, double T, double q /* q=±1 */, int N){+ const double intrinsic = fabs(std::max((q<0?K-F:F-K),0.0));+ if (price<intrinsic)+ return implied_volatility_output(0,VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC);+ const double max_price = (q<0?K:F);+ if (price>=max_price)+ return implied_volatility_output(0,VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM);+ const double x = log(F/K);+ // Map in-the-money to out-of-the-money+ if (q*x>0) {+ price = fabs(std::max(price-intrinsic,0.0));+ q = -q;+ }+ return unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_with_limited_iterations(price/(sqrt(F)*sqrt(K)), x, q, N)/sqrt(T);+}++EXPORT_EXTERN_C double implied_volatility_from_a_transformed_rational_guess(double price, double F, double K, double T, double q /* q=±1 */){+ return implied_volatility_from_a_transformed_rational_guess_with_limited_iterations(price,F,K,T,q,implied_volatility_maximum_iterations.data);+}++EXPORT_EXTERN_C double normalised_implied_volatility_from_a_transformed_rational_guess_with_limited_iterations(double beta, double x, double q /* q=±1 */, int N){+ // Map in-the-money to out-of-the-money+ if (q*x>0) {+ beta -= normalised_intrinsic(x, q);+ q = -q;+ }+ if (beta<0)+ return implied_volatility_output(0,VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC);+ return unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_with_limited_iterations(beta, x, q, N);+}++EXPORT_EXTERN_C double normalised_implied_volatility_from_a_transformed_rational_guess(double beta, double x, double q /* q=±1 */){+ return normalised_implied_volatility_from_a_transformed_rational_guess_with_limited_iterations(beta,x,q,implied_volatility_maximum_iterations.data);+}+
+ external/src/normaldistribution.cpp view
@@ -0,0 +1,147 @@+//+// normaldistribution.cpp+//++#if defined(_MSC_VER)+# define NOMINMAX // to suppress MSVC's definitions of min() and max()+// These four pragmas are the equivalent to /fp:fast.+# pragma float_control( except, off )+# pragma float_control( precise, off )+# pragma fp_contract( on )+# pragma fenv_access( off )+#endif++#include "normaldistribution.h"+#include <float.h>++namespace {+ // The asymptotic expansion Φ(z) = φ(z)/|z|·[1-1/z^2+...], Abramowitz & Stegun (26.2.12), suffices for Φ(z) to have+ // relative accuracy of 1.64E-16 for z<=-10 with 17 terms inside the square brackets (not counting the leading 1).+ // This translates to a maximum of about 9 iterations below, which is competitive with a call to erfc() and never+ // less accurate when z<=-10. Note that, as mentioned in section 4 (and discussion of figures 2 and 3) of George+ // Marsaglia's article "Evaluating the Normal Distribution" (available at http://www.jstatsoft.org/v11/a05/paper),+ // for values of x approaching -8 and below, the error of any cumulative normal function is actually dominated by+ // the hardware (or compiler implementation) accuracy of exp(-x²/2) which is not reliably more than 14 digits when+ // x becomes large. Still, we should switch to the asymptotic only when it is beneficial to do so.+ const double norm_cdf_asymptotic_expansion_first_threshold = -10.0;+ const double norm_cdf_asymptotic_expansion_second_threshold = -1/sqrt(DBL_EPSILON);+}++double norm_cdf(double z){+ if (z <= norm_cdf_asymptotic_expansion_first_threshold) {+ // Asymptotic expansion for very negative z following (26.2.12) on page 408+ // in M. Abramowitz and A. Stegun, Pocketbook of Mathematical Functions, ISBN 3-87144818-4.+ double sum = 1;+ if (z >= norm_cdf_asymptotic_expansion_second_threshold) {+ double zsqr = z * z, i = 1, g = 1, x, y, a = DBL_MAX, lasta;+ do {+ lasta = a;+ x = (4 * i - 3) / zsqr;+ y = x * ((4 * i - 1) / zsqr);+ a = g * (x - y);+ sum -= a;+ g *= y;+ ++i;+ a = fabs(a);+ } while (lasta > a && a >= fabs(sum * DBL_EPSILON));+ }+ return -norm_pdf(z) * sum / z;+ }+ return 0.5*erfc_cody( -z*ONE_OVER_SQRT_TWO );+}++double inverse_norm_cdf(double u){+ //+ // ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3+ //+ // Produces the normal deviate Z corresponding to a given lower+ // tail area of u; Z is accurate to about 1 part in 10**16.+ // see http://lib.stat.cmu.edu/apstat/241+ //+ const double split1 = 0.425;+ const double split2 = 5.0;+ const double const1 = 0.180625;+ const double const2 = 1.6;++ // Coefficients for P close to 0.5+ const double A0 = 3.3871328727963666080E0;+ const double A1 = 1.3314166789178437745E+2;+ const double A2 = 1.9715909503065514427E+3;+ const double A3 = 1.3731693765509461125E+4;+ const double A4 = 4.5921953931549871457E+4;+ const double A5 = 6.7265770927008700853E+4;+ const double A6 = 3.3430575583588128105E+4;+ const double A7 = 2.5090809287301226727E+3;+ const double B1 = 4.2313330701600911252E+1;+ const double B2 = 6.8718700749205790830E+2;+ const double B3 = 5.3941960214247511077E+3;+ const double B4 = 2.1213794301586595867E+4;+ const double B5 = 3.9307895800092710610E+4;+ const double B6 = 2.8729085735721942674E+4;+ const double B7 = 5.2264952788528545610E+3;+ // Coefficients for P not close to 0, 0.5 or 1.+ const double C0 = 1.42343711074968357734E0;+ const double C1 = 4.63033784615654529590E0;+ const double C2 = 5.76949722146069140550E0;+ const double C3 = 3.64784832476320460504E0;+ const double C4 = 1.27045825245236838258E0;+ const double C5 = 2.41780725177450611770E-1;+ const double C6 = 2.27238449892691845833E-2;+ const double C7 = 7.74545014278341407640E-4;+ const double D1 = 2.05319162663775882187E0;+ const double D2 = 1.67638483018380384940E0;+ const double D3 = 6.89767334985100004550E-1;+ const double D4 = 1.48103976427480074590E-1;+ const double D5 = 1.51986665636164571966E-2;+ const double D6 = 5.47593808499534494600E-4;+ const double D7 = 1.05075007164441684324E-9;+ // Coefficients for P very close to 0 or 1+ const double E0 = 6.65790464350110377720E0;+ const double E1 = 5.46378491116411436990E0;+ const double E2 = 1.78482653991729133580E0;+ const double E3 = 2.96560571828504891230E-1;+ const double E4 = 2.65321895265761230930E-2;+ const double E5 = 1.24266094738807843860E-3;+ const double E6 = 2.71155556874348757815E-5;+ const double E7 = 2.01033439929228813265E-7;+ const double F1 = 5.99832206555887937690E-1;+ const double F2 = 1.36929880922735805310E-1;+ const double F3 = 1.48753612908506148525E-2;+ const double F4 = 7.86869131145613259100E-4;+ const double F5 = 1.84631831751005468180E-5;+ const double F6 = 1.42151175831644588870E-7;+ const double F7 = 2.04426310338993978564E-15;++ if (u<=0)+ return log(u);+ if (u>=1)+ return log(1-u);++ const double q = u-0.5;+ if (fabs(q) <= split1)+ {+ const double r = const1 - q*q;+ return q * (((((((A7 * r + A6) * r + A5) * r + A4) * r + A3) * r + A2) * r + A1) * r + A0) /+ (((((((B7 * r + B6) * r + B5) * r + B4) * r + B3) * r + B2) * r + B1) * r + 1.0);+ }+ else+ {+ double r = q<0.0 ? u : 1.0-u;+ r = sqrt(-log(r));+ double ret;+ if (r < split2)+ {+ r = r - const2;+ ret = (((((((C7 * r + C6) * r + C5) * r + C4) * r + C3) * r + C2) * r + C1) * r + C0) /+ (((((((D7 * r + D6) * r + D5) * r + D4) * r + D3) * r + D2) * r + D1) * r + 1.0);+ }+ else+ {+ r = r - split2;+ ret = (((((((E7 * r + E6) * r + E5) * r + E4) * r + E3) * r + E2) * r + E1) * r + E0) /+ (((((((F7 * r + F6) * r + F5) * r + F4) * r + F3) * r + F2) * r + F1) * r + 1.0);+ }+ return q<0.0 ? -ret : ret;+ }+}+
+ external/src/rationalcubic.cpp view
@@ -0,0 +1,115 @@+//+// This source code resides at www.jaeckel.org/LetsBeRational.7z .+//+// ======================================================================================+// Copyright © 2013-2014 Peter Jäckel.+// +// Permission to use, copy, modify, and distribute this software is freely granted,+// provided that this notice is preserved.+//+// WARRANTY DISCLAIMER+// The Software is provided "as is" without warranty of any kind, either express or implied,+// including without limitation any implied warranties of condition, uninterrupted use,+// merchantability, fitness for a particular purpose, or non-infringement.+// ======================================================================================+//++#include "rationalcubic.h"++#if defined(_MSC_VER)+# define NOMINMAX // to suppress MSVC's definitions of min() and max()+// These four pragmas are the equivalent to /fp:fast.+// YOU NEED THESE FOR THE SAKE OF *ACCURACY* WHEN |x| IS LARGE, say, |x|>50.+// This is because they effectively enable the evaluation of certain+// expressions in 80 bit registers without loss of intermediate accuracy.+# pragma float_control( except, off )+# pragma float_control( precise, off )+# pragma fp_contract( on )+# pragma fenv_access( off )+#endif++#include <float.h>+#include <cmath>+#include <algorithm>++// Based on+//+// “Shape preserving piecewise rational interpolation”, R. Delbourgo, J.A. Gregory - SIAM journal on scientific and statistical computing, 1985 - SIAM.+// http://dspace.brunel.ac.uk/bitstream/2438/2200/1/TR_10_83.pdf [caveat emptor: there are some typographical errors in that draft version]+//++namespace {+ const double minimum_rational_cubic_control_parameter_value = -(1 - sqrt(DBL_EPSILON));+ const double maximum_rational_cubic_control_parameter_value = 2 / (DBL_EPSILON * DBL_EPSILON);+ inline bool is_zero(double x){ return fabs(x) < DBL_MIN; }+}++double rational_cubic_interpolation(double x, double x_l, double x_r, double y_l, double y_r, double d_l, double d_r, double r) {+ const double h = (x_r - x_l);+ if (fabs(h)<=0)+ return 0.5 * (y_l + y_r);+ // r should be greater than -1. We do not use assert(r > -1) here in order to allow values such as NaN to be propagated as they should.+ const double t = (x - x_l) / h;+ if ( ! (r >= maximum_rational_cubic_control_parameter_value) ) {+ const double t = (x - x_l) / h, omt = 1 - t, t2 = t * t, omt2 = omt * omt;+ // Formula (2.4) divided by formula (2.5)+ return (y_r * t2 * t + (r * y_r - h * d_r) * t2 * omt + (r * y_l + h * d_l) * t * omt2 + y_l * omt2 * omt) / (1 + (r - 3) * t * omt);+ }+ // Linear interpolation without over-or underflow.+ return y_r * t + y_l * (1 - t);+}++double rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(double x_l, double x_r, double y_l, double y_r, double d_l, double d_r, double second_derivative_l) {+ const double h = (x_r-x_l), numerator = 0.5*h*second_derivative_l+(d_r-d_l);+ if (is_zero(numerator))+ return 0;+ const double denominator = (y_r-y_l)/h-d_l;+ if (is_zero(denominator))+ return numerator>0 ? maximum_rational_cubic_control_parameter_value : minimum_rational_cubic_control_parameter_value;+ return numerator/denominator;+}++double rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(double x_l, double x_r, double y_l, double y_r, double d_l, double d_r, double second_derivative_r) {+ const double h = (x_r-x_l), numerator = 0.5*h*second_derivative_r+(d_r-d_l);+ if (is_zero(numerator))+ return 0;+ const double denominator = d_r-(y_r-y_l)/h;+ if (is_zero(denominator))+ return numerator>0 ? maximum_rational_cubic_control_parameter_value : minimum_rational_cubic_control_parameter_value;+ return numerator/denominator;+}++double minimum_rational_cubic_control_parameter(double d_l, double d_r, double s, bool preferShapePreservationOverSmoothness) {+ const bool monotonic = d_l * s >= 0 && d_r * s >= 0, convex = d_l <= s && s <= d_r, concave = d_l >= s && s >= d_r;+ if (!monotonic && !convex && !concave) // If 3==r_non_shape_preserving_target, this means revert to standard cubic.+ return minimum_rational_cubic_control_parameter_value;+ const double d_r_m_d_l = d_r - d_l, d_r_m_s = d_r - s, s_m_d_l = s - d_l;+ double r1 = -DBL_MAX, r2 = r1;+ // If monotonicity on this interval is possible, set r1 to satisfy the monotonicity condition (3.8).+ if (monotonic){+ if (!is_zero(s)) // (3.8), avoiding division by zero.+ r1 = (d_r + d_l) / s; // (3.8)+ else if (preferShapePreservationOverSmoothness) // If division by zero would occur, and shape preservation is preferred, set value to enforce linear interpolation.+ r1 = maximum_rational_cubic_control_parameter_value; // This value enforces linear interpolation.+ }+ if (convex || concave) {+ if (!(is_zero(s_m_d_l) || is_zero(d_r_m_s))) // (3.18), avoiding division by zero.+ r2 = std::max(fabs(d_r_m_d_l / d_r_m_s), fabs(d_r_m_d_l / s_m_d_l));+ else if (preferShapePreservationOverSmoothness)+ r2 = maximum_rational_cubic_control_parameter_value; // This value enforces linear interpolation.+ } else if (monotonic && preferShapePreservationOverSmoothness)+ r2 = maximum_rational_cubic_control_parameter_value; // This enforces linear interpolation along segments that are inconsistent with the slopes on the boundaries, e.g., a perfectly horizontal segment that has negative slopes on either edge.+ return std::max(minimum_rational_cubic_control_parameter_value, std::max(r1, r2));+}++double convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(double x_l, double x_r, double y_l, double y_r, double d_l, double d_r, double second_derivative_l, bool preferShapePreservationOverSmoothness) {+ const double r = rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(x_l, x_r, y_l, y_r, d_l, d_r, second_derivative_l);+ const double r_min = minimum_rational_cubic_control_parameter(d_l, d_r, (y_r-y_l)/(x_r-x_l), preferShapePreservationOverSmoothness);+ return std::max(r,r_min);+}++double convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(double x_l, double x_r, double y_l, double y_r, double d_l, double d_r, double second_derivative_r, bool preferShapePreservationOverSmoothness) {+ const double r = rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(x_l, x_r, y_l, y_r, d_l, d_r, second_derivative_r);+ const double r_min = minimum_rational_cubic_control_parameter(d_l, d_r, (y_r-y_l)/(x_r-x_l), preferShapePreservationOverSmoothness);+ return std::max(r,r_min);+}
+ hq.cabal view
@@ -0,0 +1,186 @@+cabal-version: 3.0++-- This file has been generated from package.yaml by hpack version 0.33.0.+--+-- see: https://github.com/sol/hpack+--+-- hash: ccfecaab3d1004f03e5a9e851269c0209561044c1b9aa90d587dc93a936394f0++name: hq+version: 0.1.0.0+synopsis: Quantitative Library+category: Finance+description: Please see the README on GitHub at <https://github.com/ghais/hq#readme>+homepage: https://github.com/github.com/ghais#readme+bug-reports: https://github.com/github.com/ghais/issues+author: Ghais+maintainer: 0x47@0x49.dev+copyright: 2020 Ghais Issa+license: BSD-3-Clause+license-file: LICENSE+build-type: Simple+extra-source-files:+ README.md+ ChangeLog.md++source-repository head+ type: git+ location: https://github.com/github.com/ghais/hq++library+ default-extensions:+ ConstraintKinds+ DeriveGeneric+ DerivingStrategies+ GeneralizedNewtypeDeriving+ InstanceSigs+ KindSignatures+ LambdaCase+ OverloadedStrings+ RecordWildCards+ ScopedTypeVariables+ StandaloneDeriving+ TupleSections+ TypeApplications+ ViewPatterns+ FlexibleContexts+ autogen-modules:+ Paths_hq+ exposed-modules:+ Q.SortedVector+ Q.ContingentClaim+ Q.ContingentClaim.Options+ Q.Currencies.America+ Q.Currencies.Asia+ Q.Currencies.Europe+ Q.Currency+ Q.Greeks+ Q.Interpolation+ Q.MonteCarlo+ Q.Options+ Q.Options.Bachelier+ Q.Options.Black76+ Q.Options.BlackScholes+ Q.Options.ImpliedVol+ Q.Options.ImpliedVol.TimeInterpolation+ Q.Options.ImpliedVol.InterpolatingSmile+ Q.Options.ImpliedVol.StrikeInterpolation+ Q.Options.ImpliedVol.LetsBeRational+ Q.Options.ImpliedVol.Normal+ Q.Options.ImpliedVol.Surface+ Q.Options.ImpliedVol.SVI+ Q.Options.ImpliedVol.TimeSlice+ Q.Payoff+ Q.Plotting+ Q.Stats.Arima+ Q.Stats.TimeSeries+ Q.Stochastic+ Q.Stochastic.Discretize+ Q.Stochastic.Process+ Q.Time+ Q.Time.Date+ Q.Time.DayCounter+ Q.Types+ Q.Util.File+ other-modules:+ Paths_hq+ hs-source-dirs:+ src+ include-dirs:+ external/src+ cxx-sources:+ external/src/lets_be_rational.cpp+ external/src/normaldistribution.cpp+ external/src/rationalcubic.cpp+ external/src/erf_cody.cpp+ build-depends:+ , base >=4.7 && <5+ , bytestring >=0.10 && <0.11+ , cassava >=0.5+ , containers >= 0.6.2 && <0.7+ , conversion >= 1.2 && <2+ , data-default-class >= 0.1 && <0.2+ , erf >= 2 && <3+ , hmatrix >= 0.18 && <0.3+ , hmatrix-gsl >= 0.19 && <0.20+ , hmatrix-gsl-stats >= 0.4.1+ , ieee754 >= 0.8 && <0.9+ , math-functions >= 0.3.4 && <0.4+ , mersenne-random-pure64 >= 0.2.2+ , monad-loops >= 0.4.3 && < 0.5+ , mtl >=2.2 && < 3+ , stm >= 2.5 && < 3+ , random >= 1.1 && < 2+ , random-fu >= 0.2 && < 0.3+ , random-source >= 0.3.0.11 && < 0.4+ , rvar >= 0.2 && < 0.3+ , sorted-list >= 0.2.1.0 && < 0.3+ , statistics >= 0.15.2 && < 0.16+ , text >= 1.2.4 && < 1.3+ , time >= 1.9 && < 2+ , vector >= 0.12.1 && < 0.13+ , vector-algorithms >= 0.8 && < 0.9+ default-language: Haskell2010++test-suite bachelier-test+ default-extensions:+ ConstraintKinds+ DeriveGeneric+ DerivingStrategies+ GeneralizedNewtypeDeriving+ InstanceSigs+ KindSignatures+ LambdaCase+ OverloadedStrings+ RecordWildCards+ ScopedTypeVariables+ StandaloneDeriving+ TupleSections+ TypeApplications+ ViewPatterns+ FlexibleContexts+ + type: exitcode-stdio-1.0+ main-is: Spec.hs+ other-modules:+ Paths_hq+ hs-source-dirs:+ test/bachelier+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ build-depends:+ , base >=4.7 && <5+ , hq+ , hspec >= 2.7+ , hspec-expectations >= 0.8+ default-language: Haskell2010++test-suite normalimpliedvol-test+ default-extensions:+ ConstraintKinds+ DeriveGeneric+ DerivingStrategies+ GeneralizedNewtypeDeriving+ InstanceSigs+ KindSignatures+ LambdaCase+ OverloadedStrings+ RecordWildCards+ ScopedTypeVariables+ StandaloneDeriving+ TupleSections+ TypeApplications+ ViewPatterns+ FlexibleContexts+ type: exitcode-stdio-1.0+ main-is: Spec.hs+ other-modules:+ Paths_hq+ hs-source-dirs:+ test/normalimpliedvol+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ build-depends:+ , base >=4.7 && <5+ , hq+ , hspec >= 2.7+ , hspec-expectations >= 0.8+ default-language: Haskell2010
+ src/Q/ContingentClaim.hs view
@@ -0,0 +1,68 @@+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE RankNTypes, ApplicativeDo #-}++module Q.ContingentClaim where++import Control.Monad.Reader+import Control.Monad.Writer.Strict+import Q.Types+import Data.Time+import qualified Data.Map as M++-- | A cash flow is a time and amount.+data CashFlow = CashFlow {+ cfTime :: LocalTime -- ^ The cash flow time.+ , cfAmount :: Double -- ^ The cash flow amount.+} deriving (Show, Eq)++-- | Stop at time t and potentially apply n payouts up to the monitoring time.+data CCProcessor a = CCProcessor {+ monitorTime :: LocalTime -- ^ Stopping time.+ , payouts :: [M.Map LocalTime a -> CashFlow] -- ^ list of payout functions at the stopping time.+}++-- | A claim contingent on some observable a.+newtype ContingentClaim a = ContingentClaim { unCC :: [CCProcessor a] }+-- ^ An example of an observable is a spot driven asset, such as a stock.++instance Monoid (ContingentClaim a) where+ mempty = ContingentClaim []++-- | multipley a contingent claim by its notional.+multiplier :: Double -> ContingentClaim a -> ContingentClaim a+multiplier notional (ContingentClaim ccProcessors) = ContingentClaim $ map scale ccProcessors where+ scale CCProcessor{ .. } = CCProcessor monitorTime (map scaledPayout payouts)+ scaledPayout payout = fmap (\ (CashFlow t v) -> CashFlow t (notional * v)) payout++-- | Change direction of the portfolio+short :: ContingentClaim a -> ContingentClaim a+short = multiplier (-1)+++instance Semigroup (ContingentClaim a) where+ c1 <> c2 = ContingentClaim $ combine (unCC c1) (unCC c2)+ where combine (cc1:ccs1) (cc2:ccs2)+ | monitorTime cc1 == monitorTime cc2 = let+ CCProcessor t mf = cc1+ CCProcessor _ mf' = cc2 in+ CCProcessor t (mf++mf') : combine ccs1 ccs2+ | monitorTime cc1 > monitorTime cc2 = cc2 : combine (cc1:ccs1) ccs2+ | otherwise = cc1 : combine ccs1 (cc2:ccs2)+ combine cs1 cs2 = cs1 ++ cs2++type CCBuilder w r a = WriterT w (Reader r) a++-- | Monitor an observable at the given time t.+monitor :: LocalTime -> CCBuilder (ContingentClaim a) (M.Map LocalTime a) a+monitor t = do+ tell $ ContingentClaim [CCProcessor t []] -- This step maintains the monitoring times.+ m <- ask -- This step gets the market data+ return $ m M.! t -- This step evaluate the market data at time t.++-- | Pay an amount at a given time+pay :: forall a. LocalTime -> CCBuilder (ContingentClaim a) (M.Map LocalTime a) CashFlow -> ContingentClaim a+pay t x = stoppingTimes <> ContingentClaim [CCProcessor t [payout]] where+ stoppingTimes = runReader (execWriterT x) M.empty+ payout = let r = fst <$> runWriterT x+ in runReader r+
+ src/Q/ContingentClaim/Options.hs view
@@ -0,0 +1,48 @@+module Q.ContingentClaim.Options where++import Data.Time+import Q.ContingentClaim+import Q.Types++vanillaPayout :: OptionType -- ^ Put or call+ -> Double -- ^ strike+ -> Double -- ^ Observable level+ -> Double -- ^ Payout+vanillaPayout Call k s = max (s - k) 0+vanillaPayout Put k s = max (k - s) 0+++spreadPayout :: OptionType -- ^ Put or call+ -> Double -- ^ Low strike+ -> Double -- ^ High strike+ -> Double -- ^ Observable level+ -> Double -- ^ Payout++straddlePayout :: Double -- ^ Strike+ -> Double -- ^ Observable+ -> Double -- ^ Payout+straddlePayout k s = (vanillaPayout Call k s) + (vanillaPayout Put k s)++spreadPayout Call lowStrike highStrike s = (vanillaPayout Call lowStrike s) - (vanillaPayout Call highStrike s)+spreadPayout Put lowStrike highStrike s = (vanillaPayout Put highStrike s) - (vanillaPayout Put lowStrike s)++vanillaOption :: OptionType -- ^ Option type+ -> Double -- ^ Strike+ -> LocalTime -- ^ Expiry+ -> ContingentClaim Double -- ^ Contingent claim+vanillaOption cp k t = pay t $ do+ s <- monitor t+ return $ CashFlow t $ vanillaPayout cp k s++callOption = vanillaOption Call+putOption = vanillaOption Put++-- | A call spread is a portfolio: \(C(K1, T) - C(K2 T) \) s.t. \( K1 < K2 \)+callSpread k1 k2 t = (vanillaOption Call k1 t) <> (short $ vanillaOption Call k2 t)++-- | A put spread is a portfolio: \(P(K2, T) - P(K1 T) \) s.t. \( K1 < K2 \)+putSpread k1 k2 t = (vanillaOption Put k2 t) <> (short $ vanillaOption Put k1 t)++-- | A straddle is a a portfolio :\(C(K, T) + Put(K, T)\)+straddle :: Double -> LocalTime -> ContingentClaim Double+straddle strike t = vanillaOption Put strike t <> vanillaOption Call strike t
+ src/Q/Currencies/America.hs view
@@ -0,0 +1,25 @@+module Q.Currencies.America+ (+ module Q.Currencies.America+ )+where++import Q.Currency++-- | Canadian dollar+cad :: Currency+cad = Currency {+ cName = "Canadian dollar"+ , cCode = "CAD"+ , cIsoCode = 124+ , cFracsPerUnit = 100+}++-- | U.S. dollar+usd :: Currency+usd = Currency {+ cName = "U.S. dollar"+ , cCode = "USD"+ , cIsoCode = 840+ , cFracsPerUnit = 100+}
+ src/Q/Currencies/Asia.hs view
@@ -0,0 +1,16 @@+module Q.Currencies.Asia+ (+ module Q.Currencies.Asia+ )+where++import Q.Currency++-- | Syrian Pounds+syp :: Currency+syp = Currency {+ cName = "Syrian pounds"+ , cCode = "SYP"+ , cIsoCode = 4217+ , cFracsPerUnit = 100+}
+ src/Q/Currencies/Europe.hs view
@@ -0,0 +1,34 @@+module Q.Currencies.Europe+ (+ module Q.Currencies.Europe+ )+where++import Q.Currency++-- | Swiss france+chf :: Currency+chf = Currency {+ cName = "Swiss franc"+ , cCode = "CHF"+ , cIsoCode = 756+ , cFracsPerUnit = 100+ }++-- | European Euro+eur :: Currency+eur = Currency {+ cName = "European Euro"+ , cCode = "EUR"+ , cIsoCode = 978+ , cFracsPerUnit = 100+ }++-- | British pound sterling+gbp :: Currency+gbp = Currency {+ cName = "British pound sterling"+ , cCode = "GBP"+ , cIsoCode = 826+ , cFracsPerUnit = 100+ }
+ src/Q/Currency.hs view
@@ -0,0 +1,20 @@+module Q.Currency+ (+ module Q.Currency+ )+where++-- | Currency specification+data Currency = Currency {+ -- | currency name, e.g. "U.S. dollar"+ cName :: String+ -- | ISO 4217 three-letter code, e.g. "USD"+ , cCode :: String+ -- | ISO 4217 numeric code, e.g. 840+ , cIsoCode :: Integer+ -- | number of fractionary parts in a unit+ , cFracsPerUnit :: Integer+ } deriving (Eq)++instance Show Currency where+ showsPrec _ x s = cCode x ++ s
+ src/Q/Greeks.hs view
@@ -0,0 +1,60 @@+{-# LANGUAGE MonoLocalBinds #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+module Q.Greeks+ (+ module Q.Types+ , module Q.Options+ , Bump (..)+ , DiffMethod(..)+ , Bumpable(..)+ , firstOrder+ ) where++import Q.Types+import Q.Options+import Data.Coerce++-- | A relative or an absolute bump. Used with numerical Greeks.+data Bump = Abs Double+ | Rel Double++data DiffMethod = ForwardDiff+ | BackwardDiff+ | CenteralDiff++class Bumpable a where+ bumpUp :: a -> Bump -> a+ bumpDown :: a -> Bump -> a+ stepSize :: a -> Bump -> Double++-- | Things we can bump to calculate Greeks such as 'Spot', 'Rate'..etc'+instance (Coercible a Double) => Bumpable a where+ bumpUp a (Abs bump) = coerce $ coerce a + bump+ bumpUp a (Rel bump) = coerce $ coerce a * (1 + bump)++ bumpDown a (Abs bump) = coerce $ coerce a - bump+ bumpDown a (Rel bump) = coerce $ coerce a * (1 - bump)++ stepSize _ (Abs bump) = bump+ stepSize s (Rel bump) = coerce s * bump++++firstOrder :: (Bumpable a) => DiffMethod -> Bump -> (a -> Double) -> a -> Double+firstOrder ForwardDiff b f a = df / dx+ where df = f a' - f a+ a' = bumpUp a b+ dx = stepSize a b :: Double++firstOrder BackwardDiff d f a = df / dx+ where df = f a - f a'+ a' = bumpDown a d+ dx = negate (stepSize a d )++firstOrder CenteralDiff b f a = df / dx+ where df = f u - f d+ u = bumpUp a b+ d = bumpDown a b+ dx = 2 * stepSize a b+
+ src/Q/Interpolation.hs view
@@ -0,0 +1,22 @@+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE QuantifiedConstraints #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+module Q.Interpolation where+import qualified Q.SortedVector as SV+import Numeric.GSL.Interpolation+import qualified Numeric.LinearAlgebra as V (Vector, fromList)+import Foreign (Storable)+import Data.List+class (Ord k, Storable k, Storable v) => Interpolator a k v where+ interpolate :: a -> [(k, v)] -> k -> v++class (Ord k, Storable k, Storable v) => InterpolatorV a k v where+ interpolateV :: a -> SV.SortedVector k -> V.Vector v -> k -> v++instance (Ord k, Storable k, Storable v, InterpolatorV a k v) => Interpolator a k v where+ interpolate a pts = interpolateV a xs' ys' where+ (xs, ys) = (unzip . sortOn fst) pts+ xs' = SV.fromSortedList xs+ ys' = V.fromList ys
+ src/Q/MonteCarlo.hs view
@@ -0,0 +1,87 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE InstanceSigs #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE NamedFieldPuns #-}+{-# LANGUAGE QuantifiedConstraints #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE ScopedTypeVariables #-}++module Q.MonteCarlo where+import Control.Monad.State+import Data.RVar+import Q.Stochastic.Discretize+import Q.Stochastic.Process+import Control.Monad+import Q.ContingentClaim+import Data.Random+import Q.Time+import Data.Time+import Statistics.Distribution (cumulative)+import Statistics.Distribution.Normal (standard)+import Q.ContingentClaim.Options+import Q.Types++type Path b = [(Time, b)]++-- |Summary type class aggregates all priced values of paths+class (PathPricer p) => Summary m p | m->p where+ -- | Updates summary with given priced pathes+ sSummarize :: m -> [p] -> m++ -- | Defines a metric, i.e. calculate distance between 2 summaries+ sNorm :: m -> m -> Double++-- | Path generator is a stochastic path generator+class PathGenerator m where+ pgMkNew :: m->IO m+ pgGenerate :: Integer -> m -> Path b++-- | Path pricer provides a price for given path+class PathPricer m where+ ppPrice :: m -> Path b -> m+++type MonteCarlo s a = StateT [(Time, s)] RVar a+++-- | Generate a single trajectory stopping at each provided time.+trajectory :: forall a b d. (StochasticProcess a b, Discretize d b) =>+ d -- ^ Discretization scheme+ -> a -- ^ The stochastic process+ -> b -- ^ \(S(0)\)+ -> [Time] -- ^ Stopping points \(\{t_i\}_i^n \) where \(t_i > 0\)+ -> [RVar b] -- ^ \(dW\)s. One for each stopping point.+ -> RVar [b] -- ^ \(S(0) \cup \{S(t_i)\}_i^n \) +trajectory disc p s0 times dws = reverse <$> evalStateT (onePath times dws) initState' where+ initState' :: [(Time, b)]+ initState' = [(0, s0)]++ onePath :: [Time] -> [RVar b] -> MonteCarlo b [b]+ onePath [] _ = do+ s <- get+ return $ map snd s+ onePath (t1:tn) (dw1:dws) = do+ s <- get+ let t0 = head s+ b <- lift $ pEvolve p disc t0 t1 dw1+ put $ (t1, b) : s+ onePath tn dws++-- | Generate multiple trajectories. See 'trajectory'+trajectories:: forall a b d. (StochasticProcess a b, Discretize d b) =>+ Int -- ^Num of trajectories+ -> d -- ^Discretization scheme+ -> a -- ^The stochastic process+ -> b -- ^\(S(0)\)+ -> [Time] -- ^Stopping points \(\{t_i\}_i^n \) where \(t_i > 0\)+ -> [RVar b] -- ^\(dW\)s. One for each stopping point.+ -> RVar [[b]] -- ^\(S(0) \cup \{S(t_i)\}_i^n \) +trajectories n disc p initState times dws = replicateM n $ trajectory disc p initState times dws++observationTimes :: ContingentClaim a -> [Day]+observationTimes = undefined++class Model a b | a -> b where+ discountFactor :: a -> YearFrac -> YearFrac -> RVar Rate+ evolve :: a -> YearFrac -> StateT (YearFrac, b) RVar Double
+ src/Q/Options.hs view
@@ -0,0 +1,35 @@++module Q.Options (+ Valuation(..)+ , intrinsinc+ , hasTimeValue+ , module Q.Types) where++import Numeric.IEEE+import Q.Types+++data Valuation = Valuation {+ vPremium :: Premium+ , vDelta :: Delta+ , vVega :: Vega+ , vGamma :: Gamma+} deriving (Show)+++-- | intrinsinc value of an option.+intrinsinc :: OptionType -> Forward -> Strike -> DF -> Double+intrinsinc Call (Forward f) (Strike k) (DF df) = max (f - k) 0+intrinsinc Put (Forward f) (Strike k) (DF df) = max (k - f) 0++-- | returns True if the undiscounted option premium is greater than the 'intrinsinc'+hasTimeValue ::+ OptionType+ -> Forward+ -> Strike+ -> Premium+ -> DF+ -> Bool+hasTimeValue cp f k p df = df `undiscount` p' - (intrinsinc cp f k df) > epsilon+ where (Premium p') = p+
+ src/Q/Options/Bachelier.hs view
@@ -0,0 +1,56 @@+{-# LANGUAGE MultiParamTypeClasses #-}++module Q.Options.Bachelier (+ Bachelier(..)+ , euOption+ , eucall+ , euput+ , module Q.Options+) where+import Data.Time ()+import Q.Stochastic.Discretize ()+import Q.Stochastic.Process ()+import Q.Time ()+import Statistics.Distribution (cumulative, density)+import Statistics.Distribution.Normal (standard)+import Control.Monad.State+import Data.Random (RVar, stdNormal)+import Q.MonteCarlo+import Q.Options+import Q.Types+++data Bachelier = Bachelier Forward Rate Vol deriving Show++-- | European option valuation with bachelier model.+euOption :: Bachelier -> YearFrac -> OptionType -> Strike -> Valuation+euOption (Bachelier (Forward f) (Rate r) (Vol sigma)) (YearFrac t) cp (Strike k)+ = Valuation premium delta vega gamma where+ premium = Premium $ df * (q*(f - k)*n(q*d1) + sigma*sqrt(t)/sqrt2Pi * (exp(-0.5 *d1 * d1)))+ delta = Delta $ df * n (q * d1)+ vega = Vega $ df * (sqrt t) / sqrt2Pi * (exp (-0.5 * d1 * d1))+ gamma = Gamma $ (df/(sigma * (sqrt t)))*(recip sqrt2Pi)*(exp(-0.5 *d1 * d1))+ d1 = (f - k) / (sigma * sqrt(t))+ q = cpi cp+ sqrt2Pi = sqrt (2*pi)+ df = exp $ (-r) * t+ n = cumulative standard++-- | see 'euOption'+euput b t = euOption b t Put++-- | see 'euOption'+eucall b t = euOption b t Call+++instance Model Bachelier Double where+ discountFactor (Bachelier _ r _) t1 t2 = return $ exp (scale dt r)+ where dt = t2 - t1++ evolve (Bachelier (Forward f) (Rate r) (Vol sigma)) (YearFrac t) = do+ (YearFrac t0, f0) <- get+ let dt = t - t0+ dW <- (lift stdNormal)::StateT (YearFrac, Double) RVar Double+ let ft = f0 * exp (r * dt) + sqrt(sigma*sigma/2*r * ((exp (2 * r * dt)) - 1)) * dW+ put (YearFrac t, ft)+ return ft
+ src/Q/Options/Black76.hs view
@@ -0,0 +1,60 @@+{-# LANGUAGE RecordWildCards #-}+module Q.Options.Black76+ (+ module Q.Options+ , Black76(..)+ , atmf+ , euOption+ , eucall+ , euput+ )+ where++import Q.Options+import Q.Types+import Statistics.Distribution (cumulative, density)+import Statistics.Distribution.Normal (standard)++data Black76 = Black76 {+ b76F :: Forward+ , b76DF :: DF+ , b76T :: YearFrac+ , b76Vol :: Vol+}++-- | At the money forward strike.+atmf :: Black76 -> Strike+atmf Black76{..} = Strike f+ where (Forward f) = b76F++-- | European option valuation with black 76+euOption :: Black76 -> OptionType -> Strike -> Valuation+euOption b76@Black76{..} cp k = Valuation premium delta vega gamma where+ (Forward f) = b76F+ n = cumulative standard+ (Vol sigmaSqt) = scale b76T b76Vol+ d1 = (dPlus b76F b76Vol k b76T)+ d2 = (dMinus b76F b76Vol k b76T)+ nd1 = n d1+ nd2 = n d2+ callDelta = b76DF `discount` nd1+ putDelta = b76DF `discount` (- (n (-d1)))+ vega = Vega $ b76DF `discount` (density standard d1 ) * f * sigmaSqt+ gamma = Gamma $ b76DF `discount` (density standard d1) / (f * sigmaSqt)+ premium = Premium $ case cp of+ Call -> b76DF `discount` (f * nd1 - nd2 * k')+ Put -> b76DF `discount` (n (-d2) * k' - n (-d1) * f)+ where (Strike k') = k+ delta | cp == Call = Delta $ callDelta+ | cp == Put = Delta $ putDelta++-- | see 'euOption'+euput b76 = euOption b76 Put++-- | see 'euOption'+eucall b76 = euOption b76 Call++dPlus (Forward f) (Vol sigma) (Strike k) (YearFrac t) =+ recip (sigma * sqrt t) * (log (f/k) + (0.5 * sigma * sigma) * t)+dMinus (Forward f) (Vol sigma) (Strike k) (YearFrac t) =+ recip (sigma * sqrt t) * (log (f/k) - (0.5 * sigma * sigma) * t)
+ src/Q/Options/BlackScholes.hs view
@@ -0,0 +1,85 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE RecordWildCards #-}++module Q.Options.BlackScholes (+ BlackScholes(..)+ , atmf+ , euOption+ , eucall+ , euput+ , module Q.Options+) where++import Control.Monad.State+import Data.Random hiding (Gamma)+import Data.Time+import Numeric.RootFinding+import Q.ContingentClaim.Options+import Q.MonteCarlo+import Q.Options+import Q.Stochastic.Discretize+import Q.Stochastic.Process+import Q.Time+import Q.Types+import Statistics.Distribution (cumulative, density)+import Statistics.Distribution.Normal (standard)+import qualified Q.Options.Black76 as B76++dcf = dcYearFraction ThirtyUSA++-- | Parameters for a simplified black scholes equation.+data BlackScholes = BlackScholes {+ bsSpot :: Spot -- ^ The asset's spot on the valuation date.+ , bsRate :: Rate -- ^ Risk free rate.+ , bsVol :: Vol -- ^ Volatility.+} deriving Show++++instance Model BlackScholes Double where+ discountFactor BlackScholes{..} t1 t2 = return $ exp (scale dt bsRate)+ where dt = t2 - t1++ evolve (BlackScholes spot (Rate r) (Vol sigma)) (YearFrac t) = do+ (YearFrac t0, s0) <- get+ let dt = t - t0+ dw <- (lift stdNormal)::StateT (YearFrac, Double) RVar Double+ let st = s0 * exp ((r - 0.5 * sigma * sigma) * dt + sigma * dw * sqrt dt)+ put (YearFrac t, st)+ return st++atmf :: BlackScholes -> YearFrac -> Strike+atmf BlackScholes{..} t = Strike $ s / d where+ (Rate d) = exp (scale t (-bsRate))+ (Spot s) = bsSpot++++-- | European option valuation with black scholes.+euOption :: BlackScholes -> YearFrac -> OptionType -> Strike ->Valuation+euOption bs@BlackScholes{..} t cp k =+ let b76 = B76.Black76 {+ b76F = forward bs t+ , b76DF = Q.Types.discountFactor t bsRate+ , b76T = t+ , b76Vol = bsVol+ }+ in B76.euOption b76 cp k++-- | see 'euOption'+euput bs t = euOption bs t Put+ +-- | see 'euOption'+eucall bs t = euOption bs t Call++forward BlackScholes{..} (YearFrac t) = Forward $ s * exp (r * t)+ where (Spot s) = bsSpot+ (Rate r) = bsRate++corradoMillerIniitalGuess bs@BlackScholes{..} cp (Strike k) (YearFrac t) (Premium premium) =+ (recip $ sqrt t) * ((sqrt (2 * pi)/ (s + discountedStrike)) + (premium - (s - discountedStrike)/2) + sqrt ((premium - (s - discountedStrike)/2)**2 - ((s - discountedStrike)**2/pi))) where+ discountedStrike = k * (exp $ (-r) * t)+ (Rate r) = bsRate+ (Spot s) = bsSpot++
+ src/Q/Options/ImpliedVol.hs view
@@ -0,0 +1,70 @@+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DuplicateRecordFields #-}+{-# LANGUAGE AllowAmbiguousTypes#-}++module Q.Options.ImpliedVol+ (+ module Q.Types+ , module Q.Options+ , LogRelStrike(..)+ , AbsRelStrike(..)+ , MoneynessForwardStrike(..)+ , LogMoneynessForwardStrike(..)+ , MoneynessSpotStrike(..)+ , LogMoneynessSpotStrike(..)+ , VolShift(..)+ , VolType(..)+ , euImpliedVol+ )+ where++import Q.Types+import Q.Options+import Q.Options.BlackScholes+import GHC.Generics (Generic)+import Data.Vector.Storable (Storable)+import qualified Q.Options.ImpliedVol.Normal as Bacherlier+import qualified Q.Options.ImpliedVol.LetsBeRational as B76+++newtype AbsRelStrike = AbsRel Double+ deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)++newtype LogRelStrike = LogRel Double+ deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)++newtype MoneynessForwardStrike = MoneynessForward Double+ deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)++newtype LogMoneynessForwardStrike = LogMoneynessForward Double+ deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)++newtype MoneynessSpotStrike = MoneynessSpot Double+ deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)++newtype LogMoneynessSpotStrike = LogMoneynessSpot Double+ deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)+++newtype VolShift = VolShift Double+ deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)++data VolType = Normal+ | LogNormal+ | ShiftedLogNormal VolShift+ deriving (Generic, Eq, Show, Read)+++euImpliedVol :: VolType -> OptionType -> Forward -> Strike -> YearFrac -> DF -> Premium -> Vol+euImpliedVol Normal cp f k t df premium =+ let r = rateFromDiscount t df+ in Bacherlier.euImpliedVol cp f k t r premium+euImpliedVol (ShiftedLogNormal (VolShift shift)) cp f k t df premium =+ let r = rateFromDiscount t df+ in B76.euImpliedVol cp (f + Forward shift) (k + Strike shift) t r premium+euImpliedVol LogNormal cp f k t df p = euImpliedVol (ShiftedLogNormal 0) cp f k t df p+++
+ src/Q/Options/ImpliedVol/InterpolatingSmile.hs view
@@ -0,0 +1,27 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE RecordWildCards #-}+module Q.Options.ImpliedVol.InterpolatingSmile where++import Q.Types+import Q.SortedVector (SortedVector)+import Numeric.LinearAlgebra (Vector)+import Q.Options.ImpliedVol.TimeSlice+import Q.Options.ImpliedVol.StrikeInterpolation+import Q.Interpolation+data InterpolatingSmile = StrikeSmile+ {+ smileForward :: Forward+ , smileTenor :: YearFrac+ , smileStrikes :: SortedVector Strike+ , smileVols :: Vector Vol+ , smileInterpolation :: StrikeInterpolation+ , smileExtrapolation :: StrikeExtrapolation+ , smileMinStrike :: Strike+ , smileMaxStrike :: Strike+ }++instance TimeSlice InterpolatingSmile Strike where+ totalVar smile@StrikeSmile{..} k = TotalVar $ scale smileTenor (sigma * sigma) where+ (Vol sigma) = impliedVol smile k++impliedVol StrikeSmile{..} = interpolateV smileInterpolation smileStrikes smileVols
+ src/Q/Options/ImpliedVol/LetsBeRational.hs view
@@ -0,0 +1,22 @@+{-# LANGUAGE ForeignFunctionInterface #-}+module Q.Options.ImpliedVol.LetsBeRational (+ euImpliedVol+) where++import Data.Coerce (coerce)+import Data.Number.Erf+import Foreign.C.Types+import Numeric.IEEE (epsilon, maxFinite, minNormal)+import Q.Options.BlackScholes+import Q.Options+import Q.Types+import Statistics.Distribution (cumulative, density, quantile)+import Statistics.Distribution.Normal (standard)++foreign import ccall+ "lets_be_rational.h implied_volatility_from_a_transformed_rational_guess" c_lbr ::+ CDouble -> CDouble -> CDouble -> CDouble -> CDouble -> CDouble++euImpliedVol :: OptionType -> Forward -> Strike -> YearFrac -> Rate -> Premium -> Vol+euImpliedVol cp (Forward f) (Strike k) (YearFrac t) (Rate r) (Premium p) =+ coerce $ c_lbr (CDouble p) (CDouble f) (CDouble k) (CDouble t) (CDouble (cpi cp))
+ src/Q/Options/ImpliedVol/Normal.hs view
@@ -0,0 +1,127 @@+{-# LANGUAGE RecordWildCards #-}+module Q.Options.ImpliedVol.Normal where+import Data.Default.Class+import Numeric.IEEE (epsilon, maxFinite, minNormal)+import Numeric.Natural+import Numeric.RootFinding+import Q.Options.Bachelier+import Q.Types+import Statistics.Distribution (cumulative, density)+import Statistics.Distribution.Normal (standard)+++-- | Method to use to calculate the normal implied vol.+data Method =+ Jackel -- ^ Jackel analytical formula approximation.+ | ChoKimKwak -- ^ J. Choi, K kim, and M. Kwak (2009)+ -- | Numerical root finding. Currently Ridders is used.+ | RootFinding {+ maxIter :: Natural -- ^ Maximum number of iterations.+ , tol :: Tolerance -- ^ Tolerance (relative or absolute)+ , bracket :: (Double, Double, Double) -- ^ Triple of @(low bound, initial+ -- guess, upper bound)@. If initial+ -- guess if out of bracket middle+ -- of bracket is taken as.+ }++instance Default Method where+ def = Jackel++-- | Default method implementation of 'euImpliedVolWith' using 'Jackel'.+euImpliedVol = euImpliedVolWith def++-- | Calcualte the bachelier option implied vol of a european option.+--+-- If the options premium does not have time value @'hasTimeValue'@ return 0.+euImpliedVolWith :: Method -> OptionType -> Forward -> Strike -> YearFrac -> Rate -> Premium -> Vol+euImpliedVolWith m cp f k t r p+ | hasTimeValue cp f k p df = euImpliedVolWith' m cp f k t r p+ | otherwise = Vol $ 0+ where df = discountFactor t r++euImpliedVolWith' Jackel cp (Forward f) (Strike k) (YearFrac t) (Rate r) (Premium p)+ -- Case where interest rate is not 0 we need undiscount. The paper is written+ -- for the undiscounted Bachelier option prices.+ | r /= 0+ = euImpliedVol cp (Forward f) (Strike k) (YearFrac t) (Rate 0) (Premium (p/df))+ -- Case of ATM. Solve directly.+ | abs (k - f) <= epsilon = Vol $ p * sqrt2Pi / (sqrt t)+ -- Case of ITM option. Calcualte vol of the out of the money option with Put-Call-Parity.+ | phiStarTilde >= 0+ = euImpliedVol (succ cp) (Forward f) (Strike k) (YearFrac t) (Rate r) (Premium p')+ -- General case for an out of the money option.+ | otherwise = let+ ẋ = if phiStarTilde < c then+ let g = 1 / (phiStarTilde - 0.5)+ ξ = (0.032114372355 - (g**2)*(0.016969777977 - (g**2)*(2.6207332461E-3-(9.6066952861E-5)*g**2)))+ /+ (1-(g**2)*(0.6635646938 - (g**2)*(0.14528712196 - 0.010472855461*g**2)))+ in g * (1 / sqrt2Pi + ξ*g**2)+ else+ let h = sqrt $ (-log (-phiStarTilde))+ in (9.4883409779-h*(9.6320903635-h*(0.58556997323 + 2.1464093351*h)))+ /+ (1-h*(0.65174820867 + h*(1.5120247828 + 6.6437847132E-5*h)))+ c = (-0.001882039271)+ x = ẋ + (3*q * ẋ * ẋ * (2 - q * ẋ * (2 + ẋ*ẋ)))+ /+ (6 + q*ẋ * ((-12) + ẋ *(6*q + ẋ * ((-6)*q*ẋ*(3+ẋ*ẋ)))))+ phiXBarTilde = (cumulative standard ẋ) + (density standard ẋ)/ẋ+ q = (phiXBarTilde-phiStarTilde)/ (density standard ẋ)+ in Vol $ (abs (k - f)) / (abs (x * sqrt t))+ where phiStarTilde = negate $ (abs (p - (max (theta * (f - k)) 0))) / (abs (k - f))+ theta = if cp == Call then 1 else -1+ phiTilde = (-theta) * p / (k - f)+ p' = cpi * df * (f - k) + p+ cpi = fromIntegral $ fromEnum cp --call put indicartor.+ df = exp $ (-r) * t+ sqrt2Pi = 2.506628274631000+++euImpliedVolWith' ChoKimKwak cp (Forward f) (Strike k) (YearFrac t) (Rate r) (Premium p) =+ let df = exp $ (-r) * t+ forwardPremium = p / df+ straddlePremium = case cp of Call -> 2 * forwardPremium - (f - k)+ Put -> 2 * forwardPremium + (f - k)+ nu' = (f - k) / straddlePremium+ nu = max (-1 + epsilon) (min nu' (1 - epsilon))+ eta | abs nu < sqrtEpsilon = 1+ | otherwise = nu / (atanh nu)+ heta = h eta+ in Vol $ sqrt (pi / (2 * t)) * straddlePremium * heta+++euImpliedVolWith' RootFinding{..} cp (Forward forward) k t r (Premium p) =+ let f vol = p' - p where+ (Premium p') = vPremium $ euOption b t cp k+ b = Bachelier (Forward forward) r (Vol vol)+ (lb, _, ub) = bracket+ root = ridders (RiddersParam (fromEnum maxIter) tol) (lb, ub) f+ in case root of (Root vol) -> Vol vol+ NotBracketed -> error "not bracketed"+ SearchFailed -> error "search failed"++sqrtEpsilon = sqrt epsilon+h eta = sqrt(eta) * (num / den) where+ a0 = 3.994961687345134e-1+ a1 = 2.100960795068497e+1+ a2 = 4.980340217855084e+1+ a3 = 5.988761102690991e+2+ a4 = 1.848489695437094e+3+ a5 = 6.106322407867059e+3+ a6 = 2.493415285349361e+4+ a7 = 1.266458051348246e+4++ b0 = 1.000000000000000e+0+ b1 = 4.990534153589422e+1+ b2 = 3.093573936743112e+1+ b3 = 1.495105008310999e+3+ b4 = 1.323614537899738e+3+ b5 = 1.598919697679745e+4+ b6 = 2.392008891720782e+4+ b7 = 3.608817108375034e+3+ b8 = -2.067719486400926e+2+ b9 = 1.174240599306013e+1++ num = a0 + eta * (a1 + eta * (a2 + eta * (a3 + eta * (a4 + eta * (a5 + eta * (a6 + eta * a7))))))+ den = b0 + eta * (b1 + eta * (b2 + eta * (b3 + eta * (b4 + eta * (b5 + eta * (b6 + eta * (b7 + eta * (b8 + eta * b9))))))))
+ src/Q/Options/ImpliedVol/SVI.hs view
@@ -0,0 +1,35 @@+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE MultiParamTypeClasses #-}+module Q.Options.ImpliedVol.SVI where+import Q.Types (Forward (..), Strike (..), Vol (..),+ YearFrac (..))+import Q.Options.ImpliedVol.TimeSlice+import GHC.Generics (Generic)+import Q.Greeks (Bump, Bumpable(..))++newtype Alpha = Alpha Double deriving (Generic, Eq, Show, Ord, Num, Fractional, Floating)+newtype Beta = Beta Double deriving (Generic, Eq, Show, Ord, Num, Fractional, Floating)+newtype Rho = Rho Double deriving (Generic, Eq, Show, Ord, Num, Fractional, Floating)+newtype M = M Double deriving (Generic, Eq, Show, Ord, Num, Fractional, Floating)+newtype Sigma = Sigma Double deriving (Generic, Eq, Show, Ord, Num, Fractional, Floating)++-- | Stochastic volatility inspired parameterization of the vol surface.+data SVI = RSVI -- ^ The original raw SVI representation from Gatheral+ Alpha -- ^ Corresponds to a vertical translation of the smile.+ Beta -- ^ Slope of call and put wings.+ Rho -- ^ A counter clock wise rotation of the smile.+ M -- ^ translate the smile to the right+ Sigma -- ^ ATM curviture of the smile.++instance TimeSlice SVI LogRelStrike where+ totalVar (RSVI (Alpha 𝜶) (Beta 𝜷) (Rho 𝛒) (M 𝐦) (Sigma 𝛔)) (LogRel 𝐤) =+ TotalVar $ 𝜶 + 𝜷 * (𝛒 * (𝐤 - 𝐦) + sqrt ((𝐤 - 𝐦) ** 2 + 𝛔 * 𝛔))+++isValidSVI (RSVI (Alpha 𝜶) (Beta 𝜷) (Rho 𝛒) (M 𝐦) (Sigma 𝛔)) =+ 𝜷 >= 0+ && abs 𝛒 < 1+ && 𝛔 > 0+ && 𝜶 + 𝜷 * 𝛔 * sqrt (1 -𝛒*𝛒) >= 0
+ src/Q/Options/ImpliedVol/StrikeInterpolation.hs view
@@ -0,0 +1,31 @@+{-# LANGUAGE MultiParamTypeClasses #-}+module Q.Options.ImpliedVol.StrikeInterpolation where++import Data.Coerce+import qualified Numeric.GSL.Interpolation as GSL+import Q.Interpolation+import Q.SortedVector+import Q.Types++data StrikeInterpolation = Linear+ | CubicNatural+ | CubicAkima+ | CubicMonotone++data StrikeExtrapolation = Constant+ | ConstantGradient+ | ConstantCurvature++instance InterpolatorV StrikeInterpolation Strike Vol where+ interpolateV Linear (SortedVector strikes) vols (Strike k) =+ Vol $ GSL.evaluateV GSL.Linear (coerce strikes) (coerce vols) k++ interpolateV CubicNatural (SortedVector strikes) vols (Strike k) =+ Vol $ GSL.evaluateV GSL.CSpline (coerce strikes) (coerce vols) k++ interpolateV CubicAkima (SortedVector strikes) vols (Strike k) =+ Vol $ GSL.evaluateV GSL.Akima (coerce strikes) (coerce vols) k++ interpolateV CubicMonotone (SortedVector strikes) vols (Strike k) =+ Vol $ GSL.evaluateV GSL.Steffen (coerce strikes) (coerce vols) k+
+ src/Q/Options/ImpliedVol/Surface.hs view
@@ -0,0 +1,188 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE QuantifiedConstraints #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE ScopedTypeVariables #-}+module Q.Options.ImpliedVol.Surface+ (+ Surface(..)+ , totalVarKT+ , fwdTotalVarKT+ , volKT)++where++import qualified Data.Map as M+import qualified Data.Map.Strict as M+import Data.Maybe (fromJust)+import Numeric.LinearAlgebra hiding (maxElement,+ minElement)+import qualified Q.Options.Bachelier as Bacherlier+import qualified Q.Options.Black76 as B76+import Q.Options.ImpliedVol+import Q.Options.ImpliedVol.InterpolatingSmile+import Q.Options.ImpliedVol.StrikeInterpolation+import Q.Options.ImpliedVol.TimeInterpolation+import Q.Options.ImpliedVol.TimeSlice+import Q.SortedVector+import Q.Types++-- | Implied volatility surface where the strikes are in the space of 'k' and+-- implied volatility time slice is 'v'.+data Surface v k = Surface+ {+ surfaceSpot :: Spot -- ^ Spot.+ , surfaceTenors :: SortedVector YearFrac -- ^ Ordered list of tenors.+ , surfaceForwardCurve :: YearFrac -> Forward -- ^ The forward curve.+ , surfaceDiscountCurve :: YearFrac -> DF -- ^ The discount curve.+ , surfaceAtmTotalVar :: YearFrac -> TotalVar -- ^ A spline of the at the money total variance.+ , surfaceVols :: M.Map YearFrac v -- ^ Map from tenor to 'TimeSlice'+ , surfaceTimeInterpolation :: TimeInterpolation -- ^ Method of interpolation between tenors.+ , surfaceType :: VolType -- ^ The type of surface.+ }++totalVarKT :: (StrikeSpace k, TimeSlice v k) => Surface v k -> Strike -> YearFrac -> TotalVar+totalVarKT surface@Surface{..} k t | t <= minElement surfaceTenors =+ extrapolateTotalVarFrom (minElement surfaceTenors) surface k t+ | t >= maxElement surfaceTenors =+ extrapolateTotalVarFrom (maxElement surfaceTenors) surface k t+ | otherwise =+ timeInterpolate surfaceTimeInterpolation surface k t+++volKT :: (StrikeSpace k, TimeSlice v k) => Surface v k -> Strike -> YearFrac -> Vol+volKT surface k t = totalVarToVol (totalVarKT surface k t) t++fwdTotalVarKT :: ( StrikeSpace k, TimeSlice v k) => Surface v k -> Strike -> YearFrac -> Strike -> YearFrac -> TotalVar+fwdTotalVarKT surface@Surface{..} k1 t1 k2 t2 = TotalVar $ (totalVarKT2 - totalVarKT1) / (unYearFrac t2 - unYearFrac t1)+ where (TotalVar totalVarKT1) = totalVarKT surface k1 t1+ (TotalVar totalVarKT2) = totalVarKT surface k2 t2+++class StrikeSpace k where+ strikeSpaceToCash :: k -> YearFrac -> Spot -> Forward -> Vol -> VolShift -> Strike+ cashToStrikeSpace :: Strike -> YearFrac -> Spot -> Forward -> Vol -> VolShift -> k++instance StrikeSpace Strike where+ strikeSpaceToCash x _ _ _ _ _ = x+ cashToStrikeSpace k _ _ _ _ _ = k+++instance StrikeSpace AbsRelStrike where+ strikeSpaceToCash (AbsRel x) _ _ (Forward f) _ _ = Strike $ x + f+ cashToStrikeSpace (Strike k) _ _ (Forward f) _ _ = AbsRel $ k - f++instance StrikeSpace LogRelStrike where+ strikeSpaceToCash (LogRel x) _ _ (Forward f) _ _ = Strike $ f * exp x+ cashToStrikeSpace (Strike k) _ _ (Forward f) _ _ = LogRel $ log $ k - f++instance StrikeSpace MoneynessForwardStrike where+ strikeSpaceToCash (MoneynessForward x) (YearFrac t) _ (Forward f) (Vol atmVol) _ =+ Strike $ x * sqrt t * atmVol + f++ cashToStrikeSpace (Strike k) (YearFrac t) _ (Forward f) (Vol atmVol) _ =+ MoneynessForward $ (k - f) / (atmVol * sqrt t)++instance StrikeSpace LogMoneynessForwardStrike where+ strikeSpaceToCash (LogMoneynessForward x) (YearFrac t) _ (Forward f) (Vol atmVol) (VolShift slnShift) =+ Strike $ (f + slnShift) * exp (x * (sqrt t) * atmVol) - slnShift++ cashToStrikeSpace (Strike k) (YearFrac t) _ (Forward f) (Vol atmVol) (VolShift slnShift) =+ LogMoneynessForward $ (log ((k - slnShift) / (f + slnShift))) / (atmVol * sqrt t)++instance StrikeSpace LogMoneynessSpotStrike where+ strikeSpaceToCash (LogMoneynessSpot x) (YearFrac t) (Spot s) _ (Vol atmVol) (VolShift slnShift) =+ Strike $ (s + slnShift) * exp (x * (sqrt t) * atmVol) - slnShift++ cashToStrikeSpace (Strike k) (YearFrac t) (Spot s) _ (Vol atmVol) (VolShift slnShift) =+ LogMoneynessSpot $ (log ((k - slnShift) / (s + slnShift))) / (atmVol * sqrt t)++instance StrikeSpace MoneynessSpotStrike where+ strikeSpaceToCash (MoneynessSpot x) (YearFrac t) (Spot s) _ (Vol atmVol) _ =+ Strike $ x * (sqrt t) * atmVol + s++ cashToStrikeSpace (Strike k) (YearFrac t) (Spot s) _ (Vol atmVol) _ =+ MoneynessSpot $ (k - s) / (atmVol * sqrt t)+++slnShift Surface{..} = case surfaceType of+ (ShiftedLogNormal shift) -> shift+ _ -> VolShift 0+extrapolateTotalVarFrom :: forall v k. (StrikeSpace k, TimeSlice v k) => YearFrac -> Surface v k -> Strike -> YearFrac -> TotalVar+extrapolateTotalVarFrom t0 surface@Surface{..} k t = let+ f0 = surfaceForwardCurve t0+ atmVol0 = totalVarToVol (surfaceAtmTotalVar t0) t0+ f = surfaceForwardCurve t+ spot = surfaceSpot+ atmTotalVar = surfaceAtmTotalVar t+ atmVol = totalVarToVol atmTotalVar t+ x = cashToStrikeSpace k t spot f atmVol (slnShift surface)::k+ k' = strikeSpaceToCash x t0 spot f0 atmVol0 (slnShift surface)+ x' = cashToStrikeSpace k' t0 spot f0 atmVol (slnShift surface)::k+ in totalVar (surfaceVols M.! t0) x'+++timeInterpolate :: forall v k. (StrikeSpace k, TimeSlice v k) => TimeInterpolation -> Surface v k -> Strike -> YearFrac -> TotalVar+timeInterpolate Gatheral surface@Surface{..} k11 t =+ let (t1, smile1) = fromJust $ M.lookupLE t surfaceVols+ (t2, smile2) = fromJust $ M.lookupGE t surfaceVols+ (TotalVar thetaT) = surfaceAtmTotalVar t+ (TotalVar thetaT1) = surfaceAtmTotalVar t1+ (TotalVar thetaT2) = surfaceAtmTotalVar t2+ alphaT = (sqrt thetaT2 - sqrt thetaT1 ) / (sqrt thetaT1 - sqrt thetaT)+ atmVol = totalVarToVol (TotalVar thetaT) t+ (Forward f) = surfaceForwardCurve t+ (Forward f1) = surfaceForwardCurve t1+ (Forward f2) = surfaceForwardCurve t2+ df = surfaceDiscountCurve t+ df1 = surfaceDiscountCurve t1+ df2 = surfaceDiscountCurve t2+ k1 = k11 $*$ (f1 / f)+ k2 = k11 $*$ (f2 / f)+ x1 = cashToStrikeSpace k1 t1 surfaceSpot (Forward f1) atmVol (slnShift surface)::k+ x2 = cashToStrikeSpace k2 t2 surfaceSpot (Forward f2) atmVol (slnShift surface)::k+ vol1 = totalVarToVol (totalVar smile1 x1) t1+ vol2 = totalVarToVol (totalVar smile2 x2) t2+ (Premium premium1) = vPremium $ euOption surfaceType k1 f1 df1 t1 vol1+ (Premium premium2) = vPremium $ euOption surfaceType k1 f2 df2 t1 vol1+ premiumT = Premium $ (alphaT * premium1 $/$ k1 + (1- alphaT)* premium2 $/$ k2) $*$ k11+ x1 :: k+ x2 :: k+ in if surfaceType == Normal then+ volToTotalVar (euImpliedVol Normal Call (Forward f) k11 t df premiumT) t+ else+ volToTotalVar (euImpliedVol LogNormal Call (Forward f) k11 t df premiumT) t++timeInterpolate LinearInVol surface@Surface{..} k t =+ let f = surfaceForwardCurve t+ atmVol = totalVarToVol (surfaceAtmTotalVar t) t+ x = cashToStrikeSpace k t surfaceSpot f atmVol (slnShift surface)::k+ (t1, smile1) = fromJust $ M.lookupLE t surfaceVols+ (t2, smile2) = fromJust $ M.lookupGE t surfaceVols+ (Vol vol1) = totalVarToVol (totalVar smile1 x) t1+ (Vol vol2) = totalVarToVol (totalVar smile2 x) t2+ in volToTotalVar (Vol $ linearInterpolate (t1, vol1) (t2, vol2) t) t++timeInterpolate LinearInTotalVar surface@Surface{..} k t =+ let f = surfaceForwardCurve t+ atmVol = totalVarToVol (surfaceAtmTotalVar t) t+ x = cashToStrikeSpace k t surfaceSpot f atmVol (slnShift surface)::k+ (t1, smile1) = fromJust $ M.lookupLE t surfaceVols+ (t2, smile2) = fromJust $ M.lookupGE t surfaceVols+ (TotalVar tv1) = totalVar smile1 x+ (TotalVar tv2) = totalVar smile2 x+ in TotalVar $ linearInterpolate (t1, tv1) (t2, tv2) t++euOption Normal k f (DF df) (YearFrac t) vol =+ let r = Rate $ -(log df) / t+ bacherlier = Bacherlier.Bachelier (Forward f) r vol+ in Bacherlier.eucall bacherlier (YearFrac t) k++euOption _ k f df t vol =+ let b76 = B76.Black76 (Forward f) df t vol+ in B76.eucall b76 k++linearInterpolate (YearFrac t1, v1) (YearFrac t2, v2) (YearFrac t) =+ v1 + (v2 - v1)*(t - t1) / (t2 - t1)
+ src/Q/Options/ImpliedVol/TimeInterpolation.hs view
@@ -0,0 +1,7 @@+module Q.Options.ImpliedVol.TimeInterpolation where++import Q.Options.ImpliedVol.TimeSlice+data TimeInterpolation = LinearInVol | LinearInTotalVar | Gatheral+data TimeExtrapolation = TerminalMoneyness++
+ src/Q/Options/ImpliedVol/TimeSlice.hs view
@@ -0,0 +1,30 @@+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE DuplicateRecordFields #-}+{-# LANGUAGE AllowAmbiguousTypes#-}+{-# LANGUAGE FunctionalDependencies #-}+module Q.Options.ImpliedVol.TimeSlice+ (+ module Q.Types+ , module Q.Options.ImpliedVol+ , TimeSlice(..)+ )+where++import Q.Types+import Q.Options.ImpliedVol+import Q.Options.Black76++class TimeSlice v k where+ totalVar :: v -> k -> TotalVar++instance TimeSlice (k -> TotalVar) k where+ totalVar f = f++instance TimeSlice Black76 k where+ totalVar Black76{..} _ = TotalVar $ vol * vol * t+ where (Vol vol) = b76Vol+ (YearFrac t) = b76T+
+ src/Q/Payoff.hs view
@@ -0,0 +1,47 @@+module Q.Payoff where++import Q.Types+import Q.Time+++class Payoff a where+ payoff :: (Obs1 b) => a -- ^ The instrument.+ -> b -- ^ The observable at the payoff time.+ -> Cash -- ^ Payoff amount.++-- | Path independent payoffs based on a fixed strike.+data StrikedPayoff =+ -- | Vanilla option payoff \(max (s - k, 0)\)+ -- for call and \(max (k - s, 0)\) for put+ PlainVanillaPayoff+ OptionType -- ^ Call/Put indicator+ Strike -- ^ Strike \(k\)+ -- | Payoff with strike expressed as percentage+ | PercentagePayoff+ OptionType -- ^ Call/Put indicator+ Strike -- ^ Strike in percentage.+ -- | Binary asset or nothing payoff.+ | AssetOrNothingPayoff+ OptionType -- ^ Call/Put indicator+ Strike -- ^ Strike \(k\)+ -- | Binary cash or nothing payoff.+ | CashOrNothingPayoff+ OptionType -- ^ Call/Put indicator+ Strike -- ^ Strike \(k\)+ Cash -- ^ Cash amount.+ ++instance Payoff StrikedPayoff where+ payoff (PlainVanillaPayoff cp (Strike k)) obs = Cash $ max ((cpi cp) * (s - k)) 0 where+ s = get1 obs++ payoff (PercentagePayoff cp (Strike k)) _ = Cash $ max ((cpi cp) * (1 - k)) 0++ payoff (AssetOrNothingPayoff cp (Strike k)) obs+ | (cpi cp) * (s - k) > 0 = Cash $ s+ | otherwise = Cash $ 0+ where s = get1 obs+ payoff (CashOrNothingPayoff cp (Strike k) (Cash amount)) obs+ | (cpi cp) * (s - k) > 0 = Cash $ amount+ | otherwise = 0+ where s = get1 obs
+ src/Q/Plotting.hs view
@@ -0,0 +1,10 @@+{-|+Module : Q.Plotting+Description : A collection of plotting tools i found useful.+-}+{-# LANGUAGE OverloadedStrings #-}+module Q.Plotting where+import qualified Data.Text as T++colorPairs :: [(T.Text, T.Text)]+colorPairs = cycle [("#001f3f", "#FF851B"), ("#0074D9", "#FF4136"),("#7FDBFF", "#85144b"), ("#3D9970", "#B10DC9")]
+ src/Q/SortedVector.hs view
@@ -0,0 +1,26 @@+{-# LANGUAGE FlexibleContexts #-}+module Q.SortedVector+ (+ fromList+ , fromVector+ , fromSortedList+ , SortedVector(..)+ , minElement+ , maxElement+ ) where++import qualified Data.Vector.Algorithms.Merge as V (sort)++import Data.Vector.Storable (Storable)+import qualified Data.Vector.Storable as V (Vector (..), fromList, length, head, last, modify)+import Q.Types++newtype SortedVector a = SortedVector (V.Vector a)++fromList as = SortedVector (V.modify V.sort $ V.fromList as)+fromVector v = SortedVector (V.modify V.sort v)+fromSortedList xs = SortedVector $ V.fromList xs+++minElement (SortedVector v) = V.head v+maxElement (SortedVector v) = V.last v
+ src/Q/Stats/Arima.hs view
@@ -0,0 +1,55 @@+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE InstanceSigs #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE QuantifiedConstraints #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE UndecidableInstances #-}++module Q.Stats.Arima where+import Control.Monad.State+import Data.Foldable+import Data.Functor.Identity+import Data.Random+import Data.Random.Source+import Data.RVar+import Data.Time+import Numeric.LinearAlgebra+import Q.Stats.TimeSeries+import System.Random.Mersenne.Pure64+import Data.Random.Distribution+import Data.Random.Distribution.Poisson+import Data.Random.Distribution.T+import Data.RVar+import Statistics.Sample++data Ewma d = Ewma Double d++--ll :: (Ewma d) -> [DataPoint Double] -> (Double -> Double)+ll (Ewma lambda d) datapoints = mapM ll_ datapoints where+ ll_ :: DataPoint LocalTime Double -> State Double Double+ ll_ x@(DataPoint _ v) = do+ vart <- get+ let vart2 = lambda * vart + (1 - lambda) * v * v+ put vart2+ return $ logPdf d (sqrt (v * v / vart))+++--forecast :: (Distribution d Double) => (Ewma d) -> Int ->+forecast :: forall d. (Distribution d Double) => Ewma (d Double) -> StateT Double RVar Double+forecast (Ewma lambda d) = do+ y <- lift $ rvar d+ vart <- get+ let vart2 = lambda * vart + (1 - lambda) * y * y+ put vart2+ return (y * sqrt vart)+++--forecastN :: Distribution d Double => Ewma (d Double) -> Int -> Double -> RVar ([Double], Double)+forecastN ewma var0 n = sample $ runStateT (replicateM n (forecast ewma)) var0+
+ src/Q/Stats/TimeSeries.hs view
@@ -0,0 +1,82 @@+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE RankNTypes #-}++module Q.Stats.TimeSeries where+import qualified Data.ByteString.Lazy as B+import Data.Csv ((.:))+import qualified Data.Csv as Csv+import qualified Data.Map as M+import Data.Maybe (fromJust)+import qualified Data.Text as T+import Data.Time (Day, LocalTime (LocalTime), midnight)+import Data.Time.Format ()+import Data.Time.Format.ISO8601 (FormatExtension (BasicFormat),+ calendarFormat, formatParseM,+ formatShow, localTimeFormat,+ timeOfDayFormat)+import Data.Vector (Vector, toList)+import GHC.Generics (Generic)+-- A single data point with a time and value.+data DataPoint a b = DataPoint {+ dpT :: a -- ^Time+ , dpV :: b -- ^Value+ } deriving (Generic, Show, Eq, Ord)++{-|+Read a a csv row with 2 columns: `date,value` where `date` is+in shortened iso format. (with our without time)+-}+instance Csv.FromNamedRecord (DataPoint LocalTime Double) where+ parseNamedRecord m = DataPoint+ <$> fmap (fromJust . parseDateTime) (m .: "date")+ <*> (m .: "value")++{-|+Read a a csv row with 2 columns: `date,value` where `date` is+in year fractions.+-}+instance Csv.FromNamedRecord (DataPoint Double Double) where+ parseNamedRecord m = DataPoint+ <$> (m .: "date")+ <*> (m .: "value")+++parseDateTime :: String -> Maybe LocalTime+parseDateTime iso_datetime =+ if length iso_datetime == 8 then+ parseDay iso_datetime+ else+ formatParseM localTimeFormat' iso_datetime++localTimeFormat' = localTimeFormat (calendarFormat BasicFormat) (timeOfDayFormat BasicFormat)+dayFormat' = calendarFormat BasicFormat++parseTime :: String -> Maybe LocalTime+parseTime = formatParseM localTimeFormat'++parseDay :: String -> Maybe LocalTime+parseDay iso_date = do+ day <- formatParseM dayFormat' iso_date+ return $ LocalTime day midnight++dayToString :: Day -> T.Text+dayToString = T.pack . formatShow dayFormat'++dateToString :: LocalTime -> String+dateToString = formatShow (localTimeFormat (calendarFormat BasicFormat) (timeOfDayFormat BasicFormat))++read :: forall a. (Csv.FromNamedRecord a) => FilePath -> IO [a]+read f = do+ s <- B.readFile f+ let records = Csv.decodeByName s+ case records of (Left s) -> fail s+ (Right (header, rows)) -> return $ toList rows++++valuesOnly :: [DataPoint a b] -> [b]+valuesOnly = fmap dpV++toPair (DataPoint d v) = (d, v)
+ src/Q/Stochastic.hs view
@@ -0,0 +1,5 @@+module Q.Stochastic ( module Q ) where++import Q.Stochastic.Process as Q+import Q.Stochastic.Discretize as Q+
+ src/Q/Stochastic/Discretize.hs view
@@ -0,0 +1,37 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE QuantifiedConstraints #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UndecidableInstances #-}+module Q.Stochastic.Discretize+ where++import Data.Functor+import Data.RVar+import Numeric.LinearAlgebra+import Q.Stochastic.Process+-- |Euler discretization of stochastic processes+newtype Euler = Euler { eDt :: Double }+ deriving (Show, Eq)++-- | Euler end-point discretization of stochastic processes+newtype EndEuler = EndEuler { eeDt :: Double }+ deriving (Show, Eq)+++instance Discretize Euler Double where+ dDrift p Euler{..} s0 = pDrift p s0 <&> (* eDt)+ dDiff p Euler{..} b = (pDiff p b) <&> (* (sqrt eDt))+ dDt _ Euler{..} _ = eDt++instance Discretize Euler (Vector Double) where+ dDrift p Euler{..} s0 = pDrift p s0 <&> (scale eDt)+ dDiff p Euler{..} b = (pDiff p b) <&> (scale (sqrt eDt))+ dDt _ Euler{..} _ = eDt++instance (forall a b. StochasticProcess a Double) => Discretize EndEuler Double where+ dDrift p EndEuler{..} s0@(t0, x0) = pDrift p (t0 + eeDt, x0) <&> (* eeDt)+ dDiff p EndEuler{..} s0@(t0, x0) = pDiff p (t0 + eeDt, x0) <&> (* (sqrt eeDt))+ dDt _ e _ = eeDt e
+ src/Q/Stochastic/Process.hs view
@@ -0,0 +1,97 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE QuantifiedConstraints #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TupleSections #-}+module Q.Stochastic.Process+ where+import Control.Monad+import Control.Monad.State+import Data.List (foldl')+import Data.RVar+import Data.Random+import Numeric.LinearAlgebra++rwalkState :: RVarT (State Double) Double+rwalkState = do+ prev <- lift get+ change <- rvarT StdNormal++ let new = prev + change+ lift (put new)+ return new++type Time = Double++-- Dont know why this wasn't done.+-- Is there an easier way to do this where we either lift or return?+instance (Num a) => Num (RVarT m a) where+ (+) = liftM2 (+)+ (-) = liftM2 (-)+ (*) = liftM2 (*)+ abs = liftM abs+ signum = liftM signum+ fromInteger x = return $ fromInteger x++++-- |Discretization of stochastic process over given interval+class (Num b) => Discretize d b where+ -- |Discretization of the drift process.+ dDrift :: (StochasticProcess a b) => a -> d -> (Time, b) -> RVar b+ -- |Discretization of the diffusion process.+ dDiff :: (StochasticProcess a b) => a -> d -> (Time, b) -> RVar b+ -- |dt used.+ dDt :: (StochasticProcess a b) => a -> d -> (Time, b) -> Time+++-- |A stochastic process of the form \(dX_t = \mu(X_t, t)dt + \sigma(S_t, t)dB_t \)+class (Num b) => StochasticProcess a b where+ -- |The process drift.+ pDrift :: a -> (Time, b) -> RVar b+ -- |The process diffusion.+ pDiff :: a -> (Time, b) -> RVar b++ -- |Evolve a process from a given state to a given time.+ pEvolve :: (Discretize d b) => a -- ^The process+ -> d -- ^Discretization scheme+ -> (Time, b) -- ^Initial state+ -> Time -- ^Target time t.+ -> RVar b -- ^\(dB_i\).+ -> RVar b -- ^\(X(t)\).+ pEvolve p disc s0@(t0, x0) t dw = do+ if t0 >= t then return x0 else do+ s'@(t', b') <- pEvolve' p disc s0 dw+ if t' >= t then return b' else pEvolve p disc s' t dw++ -- |Similar to evolve, but evolves the process with the discretization scheme \(dt\).+ pEvolve' :: (Discretize d b, Num b) => a -> d -> (Time, b) -> RVar b -> RVar (Time, b)+ pEvolve' process discr s@(t, b) dw = do+ let !newT = t + dDt process discr s+ !newX = do+ drift <- dDrift process discr s+ diff <- dDiff process discr s+ dw' <- dw+ return $ b + drift + diff * dw'+ newX :: RVar b++ (newT,) <$> newX++-- |Geometric Brownian motion+data GeometricBrownian = GeometricBrownian {+ gbDrift :: Double -- ^Drift+ , gbDiff :: Double -- ^Vol+} deriving (Show)+++instance StochasticProcess GeometricBrownian Double where+-- pDrift :: GeometricBrownian -> (Time, Double) -> RVar Double+ pDrift p (_, x) = return $ gbDrift p * x -- drift is prpotional to the spot.+ pDiff p (_, x) = return $ gbDiff p * x -- diffisuion is also prportional to the spot.+++-- | Ito process+data ItoProcess = ItoProcess {+ ipDrift :: (Time, Double) -> Double,+ ipDiff :: (Time, Double) -> Double+}
+ src/Q/Time.hs view
@@ -0,0 +1,53 @@+{-# LANGUAGE OverloadedStrings #-}++module Q.Time+ ( module Q.Time.Date+ , module Q.Time.DayCounter+ , parseDay+ , parseLocalTime+ ) where++import qualified Data.ByteString as B+import Data.ByteString.Char8 (unpack)+import Data.Csv (FromField (..), ToField (..), record,+ toField, (.!), (.:))+import Data.Maybe (fromJust)+import Data.Time+import Data.Time.Format+import Data.Time.Format.ISO8601+import Data.Vector (Vector, toList)+import Q.Time.Date+import Q.Time.DayCounter++-- | Converts a shortened ISO08601 date string, or datetime to a 'LocalTime'.+-- If the string doesn't contain time then 'midnight' is used.+parseLocalTime :: String -> Maybe LocalTime+parseLocalTime iso_datetime =+ if length iso_datetime == 8 then do+ day <- formatParseM dayFormat' iso_datetime+ return $ LocalTime day midnight+ else+ formatParseM localTimeFormat' iso_datetime++-- | Converts a shortned ISO08601 date to a 'Day'+parseDay :: String -> Maybe Day+parseDay = formatParseM dayFormat'+++-- | basic ISO08601 date/time format.+localTimeFormat' = localTimeFormat dayFormat' timeFormat'+-- | basic ISO08601 time format.+timeFormat' = timeOfDayFormat BasicFormat+-- | basic ISO08601 day format.+dayFormat' = calendarFormat BasicFormat++-- | Format a date as an basic ISO08601 format.+dateToString :: LocalTime -> String+dateToString date = formatShow localTimeFormat' date+++instance ToField Day where+ toField d = toField $ formatShow dayFormat' d+instance FromField Day where+ parseField s = pure $ fromJust (parseDay (unpack s))+
+ src/Q/Time/Date.hs view
@@ -0,0 +1,60 @@+{-# LANGUAGE DeriveGeneric #-}+module Q.Time.Date (Calendar(..)) where++import Data.Time+import GHC.Generics++{- |Business Day conventions+ - These conventions specify the algorithm used to adjust a date in case it is not a valid business day.+ -}+data BusinessDayConvention =+ Following -- ^Choose the first business day after the holiday + | ModifiedFollowing {- ^Choose the first business day after+ the given holiday unless it belongs+ to a different month, in which case+ choose the first business day before+ the holiday -} + | Preceding -- ^Choose the first business day before the holiday+ | ModifiedPreceding {- ^Choose the first business day before+ the given holiday unless it belongs+ to a different month, in which case+ choose the first business day after+ the holiday. -}+ | Unadjusted -- ^Do not adjust+ deriving (Generic, Show, Eq, Enum)++-- | Defines a holidays for given calendar. Corresponds to calendar class in QuantLib+class Calendar m where+ isHoliday :: m -> (Integer, Int, Int) -> Bool+ isWeekend :: m -> Day -> Bool++ isBusinessDay :: m -> Day -> Bool+ isBusinessDay m d = not (isHoliday m $ toGregorian d)++ hBusinessDayBetween :: m -> (Day, Day) -> Int+ hBusinessDayBetween m (fd, td) = foldl countDays 0 listOfDates+ where countDays counter x = counter + fromEnum (isBusinessDay m x)+ listOfDates = getDaysBetween (fd, td)++ hNextBusinessDay :: m -> Day -> Day+ hNextBusinessDay m d | isBusinessDay m nextDay = nextDay+ | otherwise = getNextBusinessDay m nextDay+ where nextDay = addDays 1 d++++-- | Generate a list of all dates inbetween+getDaysBetween :: (Day, Day) -> [Day]+getDaysBetween (fd, td) = reverse $ generator fd []+ where generator date x+ | date < td = generator nextDate (nextDate : x)+ | otherwise = x+ where nextDate = addDays 1 date++-- | Gets the next working day+getNextBusinessDay :: Calendar a => a -> Day -> Day+getNextBusinessDay m d+ | isBusinessDay m nextDay = nextDay+ | otherwise = getNextBusinessDay m nextDay+ where nextDay = addDays 1 d+
+ src/Q/Time/DayCounter.hs view
@@ -0,0 +1,53 @@+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE UndecidableInstances #-}+module Q.Time.DayCounter (+ DayCounter(..),+ Thirty360(..)+ ) where++import Data.Time.Calendar+import GHC.Generics+-- |Day counter type class+class DayCounter m where+ dcName :: m -> String -- ^Name of the day counter.+ dcCount :: m -> Day -> Day -> Int -- ^Number of business days inbetween+ dcYearFraction :: m -> Day -> Day -> Double -- ^Year fraction between 2 dates.+++-- | Thirty day counters as in QuantLib+data Thirty360 = ThirtyUSA | ThirtyEuropean | ThirtyItalian+ deriving (Generic, Eq, Show, Read)++instance DayCounter Thirty360 where+ dcName ThirtyUSA = "Thirty USA"+ dcName ThirtyEuropean = "Thirty Euro"+ dcName ThirtyItalian = "Thirty Italian"++ dcYearFraction dc fromDate toDate = fromIntegral (dcCount dc fromDate toDate) / 360.0++ dcCount ThirtyUSA fd td = 360*(yy2-yy1) + 30*(mm2-mm1-1) + max 0 (30-dd1) + min 30 dd2+ where (yy1, mm1, dd1) = intGregorian fd+ (yy2, m2, d2) = intGregorian td+ (dd2, mm2) = adjust dd1 d2 m2+ adjust x1 x2 z2+ | x2 == 31 && x1 < 30 = (1, z2+1)+ | otherwise = (x2, z2)+++ dcCount ThirtyEuropean fd td = 360*(yy2-yy1) + 30*(m2-m1-1) + max 0 (30-d1) + min 30 d2+ where (yy1, m1, d1) = intGregorian fd+ (yy2, m2, d2) = intGregorian td++ dcCount ThirtyItalian fd td = 360*(yy2-yy1) + 30*(mm2-mm1-1) + max 0 (30-dd1) + min 30 dd2+ where (yy1, mm1, d1) = intGregorian fd+ (yy2, mm2, d2) = intGregorian td+ dd1 = adjust d1 mm1+ dd2 = adjust d2 mm2+ adjust x1 z1+ | z1 == 2 && x1 > 27 = 30+ | otherwise = x1++intGregorian :: Day -> (Int, Int, Int)+intGregorian date = (fromIntegral y, m, d)+ where (y, m, d) = toGregorian date
+ src/Q/Types.hs view
@@ -0,0 +1,277 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE OverloadedStrings #-}++module Q.Types (+ Observables1(..)+ , Observables2(..)+ , Observables3(..)+ , Observables4(..)+ , Observables5(..)+ , OptionType(..)+ , Cash(..)+ , Spot(..)+ , Obs1(..)+ , Obs2(..)+ , Obs3(..)+ , Obs4(..)+ , Obs5(..)+ , Strike(..)+ , Forward(..)+ , Premium(..)+ , Delta(..)+ , Vega(..)+ , Gamma(..)+ , Expiry(..)+ , YearFrac(..)+ , Rate(..)+ , DF(..)+ , Vol(..)+ , TotalVar(..)+ , TimeScaleable(..)+ , cpi+ , discountFactor+ , discount+ , undiscount+ , rateFromDiscount+ , totalVarToVol+ , volToTotalVar+ , ($*$)+ , ($/$)+ , ($+$)+ ) where++import qualified Data.ByteString as B+import Data.Csv (FromField (..), ToField (..))+import Data.Time+import GHC.Generics (Generic)+import Q.Time+import Q.Time.Date+import Foreign (Storable)+import Numeric.LinearAlgebra (Element(..))+import Data.Coerce+-- | Type for Put or Calls+data OptionType = Put | Call deriving (Generic, Eq, Show, Read, Bounded)+instance Enum OptionType where+ succ Call = Put+ succ Put = Call++ pred = succ+ toEnum x = if signum x == 1 then Call else Put+ fromEnum Call = 1+ fromEnum Put = -1+++cpi Call = 1+cpi Put = -1++newtype Cash = Cash Double deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable) ++newtype Spot = Spot Double deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)+newtype Forward = Forward Double deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)+newtype Strike = Strike Double deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)++($*$) :: (Coercible a Double, Coercible b Double) => a -> b -> a+x1 $*$ x2 = coerce $ (coerce x1::Double) * (coerce x2::Double)++($/$) :: (Coercible a Double, Coercible b Double) => a -> b -> a+x1 $/$ x2 = coerce $ (coerce x1::Double) / (coerce x2::Double)++($+$) :: (Coercible a Double, Coercible b Double) => a -> b -> a+x1 $+$ x2 = coerce $ (coerce x1::Double) + (coerce x2::Double)++++-- Later on i should add roll.+newtype Expiry = Expiry Day deriving (Generic, Eq, Show, Read, Ord)++newtype Premium = Premium Double deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)+newtype Delta = Delta Double deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)+newtype Vega = Vega Double deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)+newtype Gamma = Gamma Double deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)++newtype YearFrac = YearFrac {unYearFrac:: Double} deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)++newtype Rate = Rate Double deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)+newtype DF = DF Double deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)++discountFactor (YearFrac t) (Rate r) = DF $ exp ((-r) * t)+discount (DF df) p = p * df+undiscount (DF df) p = p / df++rateFromDiscount (YearFrac t) (DF df) = Rate $ - (log df) / t++newtype Vol = Vol Double deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)+-- | (\w(S_0, K, T) = \sigma_{BS}(S_0, K, T)T \)+newtype TotalVar = TotalVar Double deriving (Generic, Eq, Show, Read, Ord, Num, Fractional, Real, RealFrac, RealFloat, Floating, Storable)++totalVarToVol (TotalVar v) (YearFrac t) = Vol $ sqrt (v / t)+volToTotalVar (Vol sigma) (YearFrac t) = TotalVar $ sigma * sigma * t++instance FromField OptionType where+ parseField s | (s == "C" || s == "c") = pure Call+ | (s == "P" || s == "p") = pure Put+instance ToField OptionType where+ toField Call = toField ("C"::B.ByteString)+ toField Put = toField ("P"::B.ByteString)++instance FromField Spot where+ parseField s = Spot <$> parseField s+instance ToField Spot where+ toField (Spot k) = toField k++instance FromField Cash where+ parseField s = Cash <$> parseField s+instance ToField Cash where+ toField (Cash k) = toField k+++instance FromField Strike where+ parseField s = Strike <$> parseField s+instance ToField Strike where+ toField (Strike k) = toField k++instance FromField Expiry where+ parseField s = Expiry <$> parseField s+instance ToField Expiry where+ toField (Expiry k) = toField k++instance FromField Premium where+ parseField s = Premium <$> parseField s+instance ToField Premium where+ toField (Premium k) = toField k++instance FromField Delta where+ parseField s = Delta <$> parseField s+instance ToField Delta where+ toField (Delta k) = toField k++instance FromField Vega where+ parseField s = Vega <$> parseField s+instance ToField Vega where+ toField (Vega k) = toField k++instance FromField Gamma where+ parseField s = Gamma <$> parseField s+instance ToField Gamma where+ toField (Gamma k) = toField k++instance FromField YearFrac where+ parseField s = YearFrac <$> parseField s+instance ToField YearFrac where+ toField (YearFrac k) = toField k++instance FromField Rate where+ parseField s = Rate <$> parseField s+instance ToField Rate where+ toField (Rate k) = toField k+++instance FromField Vol where+ parseField s = Vol <$> parseField s+instance ToField Vol where+ toField (Vol k) = toField k++-- | Represents concepts that scale as a function of time such as 'Vol'+class TimeScaleable a where+ scale :: YearFrac -> a -> a++instance TimeScaleable Double where+ scale (YearFrac t) y = y * t+ +instance TimeScaleable Rate where+ scale (YearFrac t) (Rate r) = Rate $ r * t+instance TimeScaleable Vol where+ scale (YearFrac t) (Vol sigma) = Vol $ sigma * sqrt t+++-- | Single-observable container.+data Observables1 = Observables1 {-# UNPACK #-} !Double+-- | Two observable container.+data Observables2 = Observables2 {-# UNPACK #-} !Double {-# UNPACK #-} !Double+-- | Three observable container.+data Observables3 = Observables3 {-# UNPACK #-} !Double {-# UNPACK #-} !Double+ {-# UNPACK #-} !Double+-- | Four observable container.+data Observables4 = Observables4 {-# UNPACK #-} !Double {-# UNPACK #-} !Double+ {-# UNPACK #-} !Double {-# UNPACK #-} !Double+-- | Five observable container.+data Observables5 = Observables5 {-# UNPACK #-} !Double {-# UNPACK #-} !Double+ {-# UNPACK #-} !Double {-# UNPACK #-} !Double+ {-# UNPACK #-} !Double++class Obs1 a where+ get1 :: a -> Double++class (Obs1 a) => Obs2 a where+ get2 :: a -> Double++class (Obs2 a) => Obs3 a where+ get3 :: a -> Double++class (Obs3 a) => Obs4 a where+ get4 :: a -> Double++class (Obs4 a) => Obs5 a where+ get5 :: a -> Double++instance Obs1 Observables1 where+ get1 (Observables1 x) = x+ {-# INLINE get1 #-}++instance Obs1 Observables2 where+ get1 (Observables2 x _) = x+ {-# INLINE get1 #-}++instance Obs1 Observables3 where+ get1 (Observables3 x _ _) = x+ {-# INLINE get1 #-}++instance Obs1 Observables4 where+ get1 (Observables4 x _ _ _) = x+ {-# INLINE get1 #-}++instance Obs1 Observables5 where+ get1 (Observables5 x _ _ _ _) = x+ {-# INLINE get1 #-}++instance Obs2 Observables2 where+ get2 (Observables2 _ x) = x+ {-# INLINE get2 #-}++instance Obs2 Observables3 where+ get2 (Observables3 _ x _) = x+ {-# INLINE get2 #-}++instance Obs2 Observables4 where+ get2 (Observables4 _ x _ _) = x+ {-# INLINE get2 #-}++instance Obs2 Observables5 where+ get2 (Observables5 _ x _ _ _) = x+ {-# INLINE get2 #-}++instance Obs3 Observables3 where+ get3 (Observables3 _ _ x) = x+ {-# INLINE get3 #-}++instance Obs3 Observables4 where+ get3 (Observables4 _ _ x _) = x+ {-# INLINE get3 #-}++instance Obs3 Observables5 where+ get3 (Observables5 _ _ x _ _) = x+ {-# INLINE get3 #-}++instance Obs4 Observables4 where+ get4 (Observables4 _ _ _ x) = x+ {-# INLINE get4 #-}++instance Obs4 Observables5 where+ get4 (Observables5 _ _ _ x _) = x+ {-# INLINE get4 #-}++instance Obs5 Observables5 where+ get5 (Observables5 _ _ _ _ x) = x+ {-# INLINE get5 #-}
+ src/Q/Util/File.hs view
@@ -0,0 +1,32 @@+{-# LANGUAGE OverloadedStrings #-}+-- |+-- Module : Data.Random.Distribution.MultivariateNormal+-- Copyright : (c) 2016 FP Complete Corporation+-- License : MIT (see LICENSE)+-- Maintainer : dominic@steinitz.org+module Q.Util.File (write)+ where++import Numeric.LinearAlgebra+import Control.Monad+import qualified Data.ByteString.Char8 as C+import Data.Csv+import Data.Char (ord)+import qualified Data.ByteString.Lazy as B+import Data.Random+++rowToRecord :: (Show t) => [t] -> Record+rowToRecord x = record $ map (C.pack . show) x++write :: (Show t) => [[t]] -> [String] -> FilePath -> IO ()+write m header path = do+ let out = (encodeWith opt s) where+ opt = defaultEncodeOptions { encDelimiter = fromIntegral (ord ','), encQuoting = QuoteNone }+ rows :: [Record]+ rows = map rowToRecord $ m+ header_ = record $ map C.pack header+ s = if null header_ then rows else header_:rows+ B.writeFile path out++
+ test/bachelier/Spec.hs view
@@ -0,0 +1,62 @@+{-# LANGUAGE OverloadedStrings #-}+module Main where+import Test.Hspec hiding (shouldBe)+import Q.Options.Bachelier+import Q.Types+import Test.Hspec.Expectations+import Control.Monad (unless)+import Q.SortedVector+closeTo x y = compareWith (\x y -> (abs $ (x - y)) <= 1e-7) errorMessage x y where+ errorMessage = "Is not close to"+ compareWith :: (HasCallStack, Show a) => (a -> a -> Bool) -> String -> a -> a -> Expectation+ compareWith comparator errorDesc result expected = expectTrue errorMsg (comparator expected result)+ where errorMsg = show result ++ " " ++ errorDesc ++ " " ++ show expected+ expectTrue msg b = unless b (expectationFailure msg)++testOptionValuation b k t v expected = do+ let p = vPremium expected+ delta = vDelta expected+ vega = vVega expected+ gamma = vGamma expected+ it ("is priced at " ++ (show p)) $ do+ vPremium v `closeTo` p+ it ("has a " ++ (show delta)) $ do+ vDelta v `closeTo` delta+ it ("has a " ++ (show vega)) $ do+ vVega v `closeTo` vega+ it ("has a " ++ (show gamma)) $ do+ vGamma v `closeTo` gamma+++main :: IO ()+main = hspec $ do+ describe "bachelier" $ do+ context "When asset price is positive ($100)" $ do+ let f = Forward 100+ context "When interest rate is zero (0%)" $ do+ let r = Rate 0+ context "When volatility is $20" $ do+ let vol = Vol 20+ context "1Y 'Call' option atm strike ($100)" $ do+ let k = Strike 100+ t = YearFrac 1+ b = Bachelier f r vol+ v = eucall b t k+ let expected = Valuation+ (Premium 7.9788456)+ (Delta 0.5)+ (Vega 0.3989422)+ (Gamma 0.01994711)+ testOptionValuation b k t v expected+ context "1Y 'Put' option atm strike ($100)" $ do+ let k = Strike 100+ t = YearFrac 1+ b = Bachelier f r vol+ v = euput b t k+ let expected = Valuation+ (Premium 7.9788456)+ (Delta 0.5)+ (Vega 0.3989422)+ (Gamma 0.01994711)+ testOptionValuation b k t v expected+
+ test/normalimpliedvol/Spec.hs view
@@ -0,0 +1,52 @@+module Main where+import Control.Monad (guard, unless, when)+import Data.List (intercalate)+import Q.Options.Bachelier+import Q.Options.ImpliedVol.Normal+import Q.Options+import Q.Types+import Test.Hspec hiding (shouldBe)+import Test.Hspec.Expectations++closeTo x y = compareWith (\x y -> (abs $ (x - y)) / (max x y) <= 1e-2) errorMessage x y where+ errorMessage = "Is not close to"+ compareWith :: (HasCallStack, Show a) => (a -> a -> Bool) -> String -> a -> a -> Expectation+ compareWith comparator errorDesc result expected = expectTrue errorMsg (comparator expected result)+ where errorMsg = show result ++ " " ++ errorDesc ++ " " ++ show expected+ expectTrue msg b = unless b (expectationFailure msg)++test cp t (k, b@(Bachelier f r sigma))= do+ let v = euOption b t cp k+ p = vPremium v+ sigma' = euImpliedVolWith Jackel cp f k t r p+ df = discountFactor t r+ when (hasTimeValue cp f k p df) $+ it (intercalate ", " [show t, show f, show df, show cp, show k, show p, show sigma]) $ do+ sigma' `closeTo` sigma+++runTests f r strikes vols t = do+ let bs = [Bachelier f r sigma | sigma <- vols]+ testCases = [(k, b) | k <- strikes, b <- bs]+ context "Call Option" $ do+ mapM_ (test Call t) testCases+ context "Put Option" $ do+ mapM_ (test Put t) testCases++main = hspec $ do+ describe "bachelier european implied vol" $ do+ context "When asset price is positive ($100)" $ do+ let strikes = [Strike k | k <- [80,81..120]]+ vols = [Vol sigma | sigma <- [1,2..200]]+ let f = Forward 100+ context "1Y option" $ do+ let t = YearFrac 1+ context "When interest rate is zero (0%)" $ do+ let r = Rate 0+ runTests f r strikes vols t+ context "When interest rate is slightly positive (1%)" $ do+ let r = Rate 0.01+ runTests f r strikes vols t+ context "When interest rate is slightly negative (-1%)" $ do+ let r = Rate (-0.01)+ runTests f r strikes vols t