diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,11 @@ +# Changelog + +`HQu` uses [PVP Versioning][1]. +The changelog is available [on GitHub][2]. + +## 0.0.0.0 + +* Initially created. + +[1]: https://pvp.haskell.org +[2]: https://github.com/ghais/HQu/releases diff --git a/HQu.cabal b/HQu.cabal new file mode 100644 --- /dev/null +++ b/HQu.cabal @@ -0,0 +1,173 @@ +cabal-version: 2.4 +name: HQu +version: 0.0.0.0 +synopsis: quantitative finance library +description: General purpose quantitative finance library +homepage: https://github.com/ghais/HQu +bug-reports: https://github.com/ghais/HQu/issues +license: MIT +license-file: LICENSE +author: Ghais +maintainer: Ghais <0x47@0x49.dev> +copyright: 2021 Ghais +category: Finance,Math +build-type: Simple +extra-doc-files: README.md + CHANGELOG.md +tested-with: GHC == 8.10.4 + +source-repository head + type: git + location: https://github.com/ghais/HQu.git + +common common-options + build-depends: base >= 4.11.1.0 && < 5 + + + ghc-options: -Wall + -Wcompat + -Widentities + -Wincomplete-uni-patterns + -Wincomplete-record-updates + if impl(ghc >= 8.0) + ghc-options: -Wredundant-constraints + if impl(ghc >= 8.2) + ghc-options: -fhide-source-paths + if impl(ghc >= 8.4) + ghc-options: -Wmissing-export-lists + -Wpartial-fields + if impl(ghc >= 8.8) + ghc-options: -Wmissing-deriving-strategies + + default-language: Haskell2010 + default-extensions: ConstraintKinds + DeriveGeneric + DerivingStrategies + GeneralizedNewtypeDeriving + InstanceSigs + KindSignatures + LambdaCase + OverloadedStrings + RecordWildCards + ScopedTypeVariables + StandaloneDeriving + TupleSections + TypeApplications + ViewPatterns + +library + import: common-options + hs-source-dirs: src + exposed-modules: HQu + 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_HQu + autogen-modules: + Paths_HQu + include-dirs: + external/include + cxx-sources: + external/src/lets_be_rational.cpp + external/src/normaldistribution.cpp + external/src/rationalcubic.cpp + external/src/erf_cody.cpp + install-includes: + importexport.h + lets_be_rational.h + normaldistribution.h + rationalcubic.h + build-depends: + 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.30 + , 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 + +test-suite bachelier-test + import: common-options + type: exitcode-stdio-1.0 + hs-source-dirs: test/bachelier + main-is: Spec.hs + build-depends: HQu + , hspec >= 2.7 + , hspec-expectations >= 0.8 + ghc-options: -threaded + -rtsopts + -with-rtsopts=-N + +test-suite normalimpliedvol-test + import: common-options + type: exitcode-stdio-1.0 + hs-source-dirs: test/normalimpliedvol + main-is: Spec.hs + build-depends: HQu + , hspec >= 2.7 + , hspec-expectations >= 0.8 + ghc-options: -threaded + -rtsopts + -with-rtsopts=-N + +benchmark HQu-benchmark + import: common-options + type: exitcode-stdio-1.0 + hs-source-dirs: benchmark + main-is: Main.hs + build-depends: gauge + , HQu + ghc-options: -threaded + -rtsopts + -with-rtsopts=-N diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2021 Ghais + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,11 @@ +# HQu + +[![GitHub CI](https://github.com/ghais/HQu/workflows/CI/badge.svg)](https://github.com/ghais/HQu/actions) +[![Build status](https://img.shields.io/travis/ghais/HQu.svg?logo=travis)](https://travis-ci.org/ghais/HQu) +[![Windows build status](https://ci.appveyor.com/api/projects/status/github/ghais/HQu?branch=master&svg=true)](https://ci.appveyor.com/project/ghais/HQu) +[![Hackage](https://img.shields.io/hackage/v/HQu.svg?logo=haskell)](https://hackage.haskell.org/package/HQu) +[![Stackage Lts](http://stackage.org/package/HQu/badge/lts)](http://stackage.org/lts/package/HQu) +[![Stackage Nightly](http://stackage.org/package/HQu/badge/nightly)](http://stackage.org/nightly/package/HQu) +[![MIT license](https://img.shields.io/badge/license-MIT-blue.svg)](LICENSE) + +General purpose quantitative finance library diff --git a/benchmark/Main.hs b/benchmark/Main.hs new file mode 100644 --- /dev/null +++ b/benchmark/Main.hs @@ -0,0 +1,7 @@ +module Main (main) where + +import Gauge.Main + + +main :: IO () +main = defaultMain [bench "const" (whnf const ())] diff --git a/external/include/importexport.h b/external/include/importexport.h new file mode 100644 --- /dev/null +++ b/external/include/importexport.h @@ -0,0 +1,36 @@ +// +// 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. +// ====================================================================================== +// +#ifndef IMPORTEXPORT_H +#define IMPORTEXPORT_H + +#if defined(_WIN32) || defined(_WIN64) +# define EXPORT __declspec(dllexport) +# define IMPORT __declspec(dllimport) +# else +# define EXPORT +# define IMPORT +#endif + +#ifdef __cplusplus +# define EXTERN_C extern "C" +#else +# define EXTERN_C +#endif + +# define EXPORT_EXTERN_C EXTERN_C EXPORT +# define IMPORT_EXTERN_C EXTERN_C IMPORT + +#endif // IMPORTEXPORT_H diff --git a/external/include/lets_be_rational.h b/external/include/lets_be_rational.h new file mode 100644 --- /dev/null +++ b/external/include/lets_be_rational.h @@ -0,0 +1,36 @@ +// +// 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. +// ====================================================================================== +// +#ifndef LETS_BE_RATIONAL_H +#define LETS_BE_RATIONAL_H + +#include "importexport.h" + +#define ENABLE_SWITCHING_THE_OUTPUT_TO_ITERATION_COUNT +#define ENABLE_CHANGING_THE_HOUSEHOLDER_METHOD_ORDER + +EXPORT_EXTERN_C double set_implied_volatility_maximum_iterations(double n); +EXPORT_EXTERN_C double set_implied_volatility_output_type(double k); +EXPORT_EXTERN_C double set_implied_volatility_householder_method_order(double m); +EXPORT_EXTERN_C double normalised_black_call(double x, double s); +EXPORT_EXTERN_C double normalised_vega(double x, double s); +EXPORT_EXTERN_C double normalised_black(double x, double s, double q /* q=±1 */); +EXPORT_EXTERN_C double black(double F, double K, double sigma, double T, double q /* q=±1 */); +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); +EXPORT_EXTERN_C double normalised_implied_volatility_from_a_transformed_rational_guess(double beta, double x, double q /* q=±1 */); +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); +EXPORT_EXTERN_C double implied_volatility_from_a_transformed_rational_guess(double price, double F, double K, double T, double q /* q=±1 */); + +#endif // NORMAL_DISTRIBUTION_H diff --git a/external/include/normaldistribution.h b/external/include/normaldistribution.h new file mode 100644 --- /dev/null +++ b/external/include/normaldistribution.h @@ -0,0 +1,34 @@ +// +// 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. +// ====================================================================================== +// +#ifndef NORMAL_DISTRIBUTION_H +#define NORMAL_DISTRIBUTION_H + +#include +#include +#include "importexport.h" + +#define ONE_OVER_SQRT_TWO 0.7071067811865475244008443621048490392848359376887 +#define ONE_OVER_SQRT_TWO_PI 0.3989422804014326779399460599343818684758586311649 +#define SQRT_TWO_PI 2.506628274631000502415765284811045253006986740610 + +EXPORT_EXTERN_C double erf_cody(double z); +EXPORT_EXTERN_C double erfc_cody(double z); +EXPORT_EXTERN_C double erfcx_cody(double z); +EXPORT_EXTERN_C double norm_cdf(double z); +inline double norm_pdf(double x){ return ONE_OVER_SQRT_TWO_PI*exp(-.5*x*x); } +EXPORT_EXTERN_C double inverse_norm_cdf(double u); + +#endif // NORMAL_DISTRIBUTION_H diff --git a/external/include/rationalcubic.h b/external/include/rationalcubic.h new file mode 100644 --- /dev/null +++ b/external/include/rationalcubic.h @@ -0,0 +1,34 @@ +// +// 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. +// ====================================================================================== +// +#ifndef RATIONAL_CUBIC_H +#define RATIONAL_CUBIC_H + +// 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] +// + +#include "importexport.h" + +EXPORT_EXTERN_C 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); +EXPORT_EXTERN_C 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); +EXPORT_EXTERN_C 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); +EXPORT_EXTERN_C double minimum_rational_cubic_control_parameter(double d_l, double d_r, double s, bool preferShapePreservationOverSmoothness); +EXPORT_EXTERN_C 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); +EXPORT_EXTERN_C 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); + +#endif // RATIONAL_CUBIC_H diff --git a/external/src/erf_cody.cpp b/external/src/erf_cody.cpp new file mode 100644 --- /dev/null +++ b/external/src/erf_cody.cpp @@ -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 +#include + +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_ */ diff --git a/external/src/lets_be_rational.cpp b/external/src/lets_be_rational.cpp new file mode 100644 --- /dev/null +++ b/external/src/lets_be_rational.cpp @@ -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 +#include +#include +#if defined(_WIN32) || defined(_WIN64) +# include +#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 beta0) { + 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=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 ( beta0)) { // 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 (; iterationsDBL_EPSILON*s; ++iterations){ + if (ds*ds_previous<0) + ++direction_reversal_count; + if ( iterations>0 && ( 3==direction_reversal_count || !(s>s_left && s 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 && ss_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_beta20.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 (; iterationsDBL_EPSILON*s; ++iterations){ + if (ds*ds_previous<0) + ++direction_reversal_count; + if ( iterations>0 && ( 3==direction_reversal_count || !(s>s_left && s 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 && ss_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 (; iterationsDBL_EPSILON*s; ++iterations){ + if (ds*ds_previous<0) + ++direction_reversal_count; + if ( iterations>0 && ( 3==direction_reversal_count || !(s>s_left && s 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 && ss_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=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); +} + diff --git a/external/src/normaldistribution.cpp b/external/src/normaldistribution.cpp new file mode 100644 --- /dev/null +++ b/external/src/normaldistribution.cpp @@ -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 + +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; + } +} + diff --git a/external/src/rationalcubic.cpp b/external/src/rationalcubic.cpp new file mode 100644 --- /dev/null +++ b/external/src/rationalcubic.cpp @@ -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 +#include +#include + +// 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); +} diff --git a/src/HQu.hs b/src/HQu.hs new file mode 100644 --- /dev/null +++ b/src/HQu.hs @@ -0,0 +1,15 @@ +{- | +Copyright: (c) 2021 Ghais +SPDX-License-Identifier: MIT +Maintainer: Ghais <0x47@0x49.dev> + +General purpose quantitative finance library +-} + +module HQu + ( someFunc + ) where + + +someFunc :: IO () +someFunc = putStrLn ("someFunc" :: String) diff --git a/src/Q/ContingentClaim.hs b/src/Q/ContingentClaim.hs new file mode 100644 --- /dev/null +++ b/src/Q/ContingentClaim.hs @@ -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 + diff --git a/src/Q/ContingentClaim/Options.hs b/src/Q/ContingentClaim/Options.hs new file mode 100644 --- /dev/null +++ b/src/Q/ContingentClaim/Options.hs @@ -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 diff --git a/src/Q/Currencies/America.hs b/src/Q/Currencies/America.hs new file mode 100644 --- /dev/null +++ b/src/Q/Currencies/America.hs @@ -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 +} diff --git a/src/Q/Currencies/Asia.hs b/src/Q/Currencies/Asia.hs new file mode 100644 --- /dev/null +++ b/src/Q/Currencies/Asia.hs @@ -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 +} diff --git a/src/Q/Currencies/Europe.hs b/src/Q/Currencies/Europe.hs new file mode 100644 --- /dev/null +++ b/src/Q/Currencies/Europe.hs @@ -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 + } diff --git a/src/Q/Currency.hs b/src/Q/Currency.hs new file mode 100644 --- /dev/null +++ b/src/Q/Currency.hs @@ -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 diff --git a/src/Q/Greeks.hs b/src/Q/Greeks.hs new file mode 100644 --- /dev/null +++ b/src/Q/Greeks.hs @@ -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 + diff --git a/src/Q/Interpolation.hs b/src/Q/Interpolation.hs new file mode 100644 --- /dev/null +++ b/src/Q/Interpolation.hs @@ -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 diff --git a/src/Q/MonteCarlo.hs b/src/Q/MonteCarlo.hs new file mode 100644 --- /dev/null +++ b/src/Q/MonteCarlo.hs @@ -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 diff --git a/src/Q/Options.hs b/src/Q/Options.hs new file mode 100644 --- /dev/null +++ b/src/Q/Options.hs @@ -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 + diff --git a/src/Q/Options/Bachelier.hs b/src/Q/Options/Bachelier.hs new file mode 100644 --- /dev/null +++ b/src/Q/Options/Bachelier.hs @@ -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 diff --git a/src/Q/Options/Black76.hs b/src/Q/Options/Black76.hs new file mode 100644 --- /dev/null +++ b/src/Q/Options/Black76.hs @@ -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) diff --git a/src/Q/Options/BlackScholes.hs b/src/Q/Options/BlackScholes.hs new file mode 100644 --- /dev/null +++ b/src/Q/Options/BlackScholes.hs @@ -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 + + diff --git a/src/Q/Options/ImpliedVol.hs b/src/Q/Options/ImpliedVol.hs new file mode 100644 --- /dev/null +++ b/src/Q/Options/ImpliedVol.hs @@ -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 + + + diff --git a/src/Q/Options/ImpliedVol/InterpolatingSmile.hs b/src/Q/Options/ImpliedVol/InterpolatingSmile.hs new file mode 100644 --- /dev/null +++ b/src/Q/Options/ImpliedVol/InterpolatingSmile.hs @@ -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 diff --git a/src/Q/Options/ImpliedVol/LetsBeRational.hs b/src/Q/Options/ImpliedVol/LetsBeRational.hs new file mode 100644 --- /dev/null +++ b/src/Q/Options/ImpliedVol/LetsBeRational.hs @@ -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)) diff --git a/src/Q/Options/ImpliedVol/Normal.hs b/src/Q/Options/ImpliedVol/Normal.hs new file mode 100644 --- /dev/null +++ b/src/Q/Options/ImpliedVol/Normal.hs @@ -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)))))))) diff --git a/src/Q/Options/ImpliedVol/SVI.hs b/src/Q/Options/ImpliedVol/SVI.hs new file mode 100644 --- /dev/null +++ b/src/Q/Options/ImpliedVol/SVI.hs @@ -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 diff --git a/src/Q/Options/ImpliedVol/StrikeInterpolation.hs b/src/Q/Options/ImpliedVol/StrikeInterpolation.hs new file mode 100644 --- /dev/null +++ b/src/Q/Options/ImpliedVol/StrikeInterpolation.hs @@ -0,0 +1,34 @@ +{-# 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) = + -- The interpolation method should be Steffen but until the next + -- version of hmatrix-gsl is release i am using Akima. + Vol $ GSL.evaluateV GSL.Akima (coerce strikes) (coerce vols) k + diff --git a/src/Q/Options/ImpliedVol/Surface.hs b/src/Q/Options/ImpliedVol/Surface.hs new file mode 100644 --- /dev/null +++ b/src/Q/Options/ImpliedVol/Surface.hs @@ -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) diff --git a/src/Q/Options/ImpliedVol/TimeInterpolation.hs b/src/Q/Options/ImpliedVol/TimeInterpolation.hs new file mode 100644 --- /dev/null +++ b/src/Q/Options/ImpliedVol/TimeInterpolation.hs @@ -0,0 +1,7 @@ +module Q.Options.ImpliedVol.TimeInterpolation where + +import Q.Options.ImpliedVol.TimeSlice +data TimeInterpolation = LinearInVol | LinearInTotalVar | Gatheral +data TimeExtrapolation = TerminalMoneyness + + diff --git a/src/Q/Options/ImpliedVol/TimeSlice.hs b/src/Q/Options/ImpliedVol/TimeSlice.hs new file mode 100644 --- /dev/null +++ b/src/Q/Options/ImpliedVol/TimeSlice.hs @@ -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 + diff --git a/src/Q/Payoff.hs b/src/Q/Payoff.hs new file mode 100644 --- /dev/null +++ b/src/Q/Payoff.hs @@ -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 diff --git a/src/Q/Plotting.hs b/src/Q/Plotting.hs new file mode 100644 --- /dev/null +++ b/src/Q/Plotting.hs @@ -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")] diff --git a/src/Q/SortedVector.hs b/src/Q/SortedVector.hs new file mode 100644 --- /dev/null +++ b/src/Q/SortedVector.hs @@ -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 diff --git a/src/Q/Stats/Arima.hs b/src/Q/Stats/Arima.hs new file mode 100644 --- /dev/null +++ b/src/Q/Stats/Arima.hs @@ -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 + diff --git a/src/Q/Stats/TimeSeries.hs b/src/Q/Stats/TimeSeries.hs new file mode 100644 --- /dev/null +++ b/src/Q/Stats/TimeSeries.hs @@ -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) diff --git a/src/Q/Stochastic.hs b/src/Q/Stochastic.hs new file mode 100644 --- /dev/null +++ b/src/Q/Stochastic.hs @@ -0,0 +1,5 @@ +module Q.Stochastic ( module Q ) where + +import Q.Stochastic.Process as Q +import Q.Stochastic.Discretize as Q + diff --git a/src/Q/Stochastic/Discretize.hs b/src/Q/Stochastic/Discretize.hs new file mode 100644 --- /dev/null +++ b/src/Q/Stochastic/Discretize.hs @@ -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 diff --git a/src/Q/Stochastic/Process.hs b/src/Q/Stochastic/Process.hs new file mode 100644 --- /dev/null +++ b/src/Q/Stochastic/Process.hs @@ -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 +} diff --git a/src/Q/Time.hs b/src/Q/Time.hs new file mode 100644 --- /dev/null +++ b/src/Q/Time.hs @@ -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)) + diff --git a/src/Q/Time/Date.hs b/src/Q/Time/Date.hs new file mode 100644 --- /dev/null +++ b/src/Q/Time/Date.hs @@ -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 + diff --git a/src/Q/Time/DayCounter.hs b/src/Q/Time/DayCounter.hs new file mode 100644 --- /dev/null +++ b/src/Q/Time/DayCounter.hs @@ -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 diff --git a/src/Q/Types.hs b/src/Q/Types.hs new file mode 100644 --- /dev/null +++ b/src/Q/Types.hs @@ -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 #-} diff --git a/src/Q/Util/File.hs b/src/Q/Util/File.hs new file mode 100644 --- /dev/null +++ b/src/Q/Util/File.hs @@ -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 + + diff --git a/test/bachelier/Spec.hs b/test/bachelier/Spec.hs new file mode 100644 --- /dev/null +++ b/test/bachelier/Spec.hs @@ -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 + diff --git a/test/normalimpliedvol/Spec.hs b/test/normalimpliedvol/Spec.hs new file mode 100644 --- /dev/null +++ b/test/normalimpliedvol/Spec.hs @@ -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