HestonStochasticVolatilityAlgorithm.java

  1. package org.drip.pricer.option;

  3. /*
  4.  * -*- mode: java; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
  5.  */

  7. /*!
  8.  * Copyright (C) 2020 Lakshmi Krishnamurthy
  9.  * Copyright (C) 2019 Lakshmi Krishnamurthy
  10.  * Copyright (C) 2018 Lakshmi Krishnamurthy
  11.  * Copyright (C) 2017 Lakshmi Krishnamurthy
  12.  * Copyright (C) 2016 Lakshmi Krishnamurthy
  13.  * Copyright (C) 2015 Lakshmi Krishnamurthy
  14.  * Copyright (C) 2014 Lakshmi Krishnamurthy
  15.  * 
  16.  *  This file is part of DROP, an open-source library targeting analytics/risk, transaction cost analytics,
  17.  *      asset liability management analytics, capital, exposure, and margin analytics, valuation adjustment
  18.  *      analytics, and portfolio construction analytics within and across fixed income, credit, commodity,
  19.  *      equity, FX, and structured products. It also includes auxiliary libraries for algorithm support,
  20.  *      numerical analysis, numerical optimization, spline builder, model validation, statistical learning,
  21.  *      and computational support.
  22.  *  
  23.  *      https://lakshmidrip.github.io/DROP/
  24.  *  
  25.  *  DROP is composed of three modules:
  26.  *  
  27.  *  - DROP Product Core - https://lakshmidrip.github.io/DROP-Product-Core/
  28.  *  - DROP Portfolio Core - https://lakshmidrip.github.io/DROP-Portfolio-Core/
  29.  *  - DROP Computational Core - https://lakshmidrip.github.io/DROP-Computational-Core/
  30.  * 
  31.  *  DROP Product Core implements libraries for the following:
  32.  *  - Fixed Income Analytics
  33.  *  - Loan Analytics
  34.  *  - Transaction Cost Analytics
  35.  * 
  36.  *  DROP Portfolio Core implements libraries for the following:
  37.  *  - Asset Allocation Analytics
  38.  *  - Asset Liability Management Analytics
  39.  *  - Capital Estimation Analytics
  40.  *  - Exposure Analytics
  41.  *  - Margin Analytics
  42.  *  - XVA Analytics
  43.  * 
  44.  *  DROP Computational Core implements libraries for the following:
  45.  *  - Algorithm Support
  46.  *  - Computation Support
  47.  *  - Function Analysis
  48.  *  - Model Validation
  49.  *  - Numerical Analysis
  50.  *  - Numerical Optimizer
  51.  *  - Spline Builder
  52.  *  - Statistical Learning
  53.  * 
  54.  *  Documentation for DROP is Spread Over:
  55.  * 
  56.  *  - Main                     => https://lakshmidrip.github.io/DROP/
  57.  *  - Wiki                     => https://github.com/lakshmiDRIP/DROP/wiki
  58.  *  - GitHub                   => https://github.com/lakshmiDRIP/DROP
  59.  *  - Repo Layout Taxonomy     => https://github.com/lakshmiDRIP/DROP/blob/master/Taxonomy.md
  60.  *  - Javadoc                  => https://lakshmidrip.github.io/DROP/Javadoc/index.html
  61.  *  - Technical Specifications => https://github.com/lakshmiDRIP/DROP/tree/master/Docs/Internal
  62.  *  - Release Versions         => https://lakshmidrip.github.io/DROP/version.html
  63.  *  - Community Credits        => https://lakshmidrip.github.io/DROP/credits.html
  64.  *  - Issues Catalog           => https://github.com/lakshmiDRIP/DROP/issues
  65.  *  - JUnit                    => https://lakshmidrip.github.io/DROP/junit/index.html
  66.  *  - Jacoco                   => https://lakshmidrip.github.io/DROP/jacoco/index.html
  67.  * 
  68.  *  Licensed under the Apache License, Version 2.0 (the "License");
  69.  *      you may not use this file except in compliance with the License.
  70.  *   
  71.  *  You may obtain a copy of the License at
  72.  *      http://www.apache.org/licenses/LICENSE-2.0
  73.  *  
  74.  *  Unless required by applicable law or agreed to in writing, software
  75.  *      distributed under the License is distributed on an "AS IS" BASIS,
  76.  *      WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  77.  *  
  78.  *  See the License for the specific language governing permissions and
  79.  *      limitations under the License.
  80.  */

  82. /**
  83.  * <i>HestonStochasticVolatilityAlgorithm</i> implements the Heston 1993 Stochastic Volatility European Call
  84.  * and Put Options Pricer.
  85.  *
  86.  *  <br><br>
  87.  *  <ul>
  88.  *      <li><b>Module </b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/ProductCore.md">Product Core Module</a></li>
  89.  *      <li><b>Library</b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/FixedIncomeAnalyticsLibrary.md">Fixed Income Analytics</a></li>
  90.  *      <li><b>Project</b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/src/main/java/org/drip/pricer/README.md">Custom Pricing Algorithms and the Derivative Fokker Planck Trajectory Generators</a></li>
  91.  *      <li><b>Package</b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/src/main/java/org/drip/pricer/option/README.md">Deterministic/Stochastic Volatility Settings/Greeks</a></li>
  92.  *  </ul>
  93.  * <br><br>
  94.  * 
  95.  * @author Lakshmi Krishnamurthy
  96.  */

  98. public class HestonStochasticVolatilityAlgorithm extends org.drip.pricer.option.FokkerPlanckGenerator {

  100.     /**
  101.      * Payoff Transformation Type - The Original Heston 1993 Scheme
  102.      */

  104.     public static final int PAYOFF_TRANSFORM_SCHEME_HESTON_1993 = 1;

  106.     /**
  107.      * Payoff Transformation Type - The Albrecher, Mayer, Schoutens, and Tistaert Scheme
  108.      */

  110.     public static final int PAYOFF_TRANSFORM_SCHEME_AMST_2007 = 1;

  112.     private static final double FOURIER_FREQ_INIT = 0.01;
  113.     private static final double FOURIER_FREQ_INCREMENT = 0.1;
  114.     private static final double FOURIER_FREQ_FINAL = 25.;

  116.     private org.drip.param.pricer.HestonOptionPricerParams _fphp = null;

  118.     class PhaseCorrectedF {
  119.         double _dblCorrectedPhase = java.lang.Double.NaN;
  120.         org.drip.function.definition.CartesianComplexNumber _cnF = null;

  122.         PhaseCorrectedF (
  123.             final double dblCorrectedPhase,
  124.             final org.drip.function.definition.CartesianComplexNumber cnF)
  125.         {
  126.             _cnF = cnF;
  127.             _dblCorrectedPhase = dblCorrectedPhase;
  128.         }
  129.     }

  131.     private PhaseCorrectedF fourierTransformHeston93 (
  132.         final double dblStrike,
  133.         final double dbTimeToExpiry,
  134.         final double dblRiskFreeRate,
  135.         final double dblSpot,
  136.         final double dblInitialVolatility,
  137.         final double dblA,
  138.         final double dblFreq,
  139.         final double dblB,
  140.         final double dblU,
  141.         final org.drip.numerical.fourier.RotationCountPhaseTracker rcpt)
  142.     {
  143.         try {
  144.             org.drip.function.definition.CartesianComplexNumber cnSmallDLHS = new org.drip.function.definition.CartesianComplexNumber (dblB,
  145.                 -1. * _fphp.rho() * _fphp.sigma() * dblFreq);

  147.             org.drip.function.definition.CartesianComplexNumber cnSmallD = org.drip.function.definition.CartesianComplexNumber.Square
  148.                 (cnSmallDLHS);

  150.             if (null == cnSmallD) return null;

  152.             double dblSigmaScaler = _fphp.sigma() * _fphp.sigma();

  154.             if (null == (cnSmallD = org.drip.function.definition.CartesianComplexNumber.Add (cnSmallD, new
  155.                 org.drip.function.definition.CartesianComplexNumber (dblSigmaScaler * dblFreq * dblFreq, -2. * dblSigmaScaler
  156.                     * dblFreq * dblU))))
  157.                 return null;

  159.             if (null == (cnSmallD = org.drip.function.definition.CartesianComplexNumber.SquareRoot (cnSmallD))) return null;

  161.             org.drip.function.definition.CartesianComplexNumber cnGNumerator = org.drip.function.definition.CartesianComplexNumber.Subtract
  162.                 (cnSmallDLHS, cnSmallD);

  164.             if (null == cnGNumerator) return null;

  166.             org.drip.function.definition.CartesianComplexNumber cnG = org.drip.function.definition.CartesianComplexNumber.Add (cnSmallDLHS,
  167.                 cnSmallD);

  169.             if (null == cnG) return null;

  171.             if (null == (cnG = org.drip.function.definition.CartesianComplexNumber.Divide (cnGNumerator, cnG))) return null;

  173.             int iM = 0;
  174.             int iN = 0;

  176.             if (org.drip.numerical.fourier.PhaseAdjuster.MULTI_VALUE_BRANCH_POWER_PHASE_TRACKER_KAHL_JACKEL ==
  177.                 _fphp.phaseTrackerType()) {
  178.                 iM = (int) ((cnG.argument() + java.lang.Math.PI) / (2. * java.lang.Math.PI));

  180.                 iN = (int) ((cnG.argument() + (dbTimeToExpiry * cnSmallD.argument()) + java.lang.Math.PI) /
  181.                     (2. * java.lang.Math.PI));
  182.             }

  184.             org.drip.function.definition.CartesianComplexNumber cnExpTTEScaledSmallD =
  185.                 org.drip.function.definition.CartesianComplexNumber.Scale (cnSmallD, -1. * dbTimeToExpiry);

  187.             if (null == cnExpTTEScaledSmallD) return null;

  189.             if (null == (cnExpTTEScaledSmallD = org.drip.function.definition.CartesianComplexNumber.Exponentiate
  190.                 (cnExpTTEScaledSmallD)))
  191.                 return null;

  193.             org.drip.function.definition.CartesianComplexNumber cnD = new org.drip.function.definition.CartesianComplexNumber (1. -
  194.                 cnExpTTEScaledSmallD.real(), -1. * cnExpTTEScaledSmallD.imaginary());

  196.             org.drip.function.definition.CartesianComplexNumber cnInvGExpTTEScaledSmallD =
  197.                 org.drip.function.definition.CartesianComplexNumber.Multiply (cnExpTTEScaledSmallD, cnG);

  199.             if (null == cnInvGExpTTEScaledSmallD) return null;

  201.             cnInvGExpTTEScaledSmallD = new org.drip.function.definition.CartesianComplexNumber (1. -
  202.                 cnInvGExpTTEScaledSmallD.real(), -1. * cnInvGExpTTEScaledSmallD.imaginary());

  204.             if (null == (cnD = org.drip.function.definition.CartesianComplexNumber.Divide (cnD, cnInvGExpTTEScaledSmallD)))
  205.                 return null;

  207.             if (null == (cnD = org.drip.function.definition.CartesianComplexNumber.Multiply (cnGNumerator, cnD)))
  208.                 return null;

  210.             dblSigmaScaler = 1. / dblSigmaScaler;

  212.             if (null == (cnD = org.drip.function.definition.CartesianComplexNumber.Scale (cnD, dblSigmaScaler))) return null;

  214.             org.drip.function.definition.CartesianComplexNumber cnC = new org.drip.function.definition.CartesianComplexNumber (1. -
  215.                 cnG.real(), -1. * cnG.imaginary());

  217.             if (org.drip.numerical.fourier.PhaseAdjuster.MULTI_VALUE_BRANCH_POWER_PHASE_TRACKER_KAHL_JACKEL ==
  218.                 _fphp.phaseTrackerType()) {
  219.                 if (null == (cnC = org.drip.numerical.fourier.PhaseAdjuster.PowerLogPhaseTracker
  220.                     (cnInvGExpTTEScaledSmallD, cnC, iN, iM)))
  221.                     return null;
  222.             } else if (org.drip.numerical.fourier.PhaseAdjuster.MULTI_VALUE_BRANCH_PHASE_TRACKER_ROTATION_COUNT
  223.                 == _fphp.phaseTrackerType()) {
  224.                 if (null == (cnC = org.drip.function.definition.CartesianComplexNumber.Logarithm (cnC))) return null;

  226.                 cnC = new org.drip.function.definition.CartesianComplexNumber (cnC.real(), rcpt.updateAndApply
  227.                     (cnC.argument(), true));
  228.             }

  230.             double dblCorrectedPhase = cnC.argument();

  232.             if (null == (cnC = org.drip.function.definition.CartesianComplexNumber.Scale (cnC, -2.))) return null;

  234.             org.drip.function.definition.CartesianComplexNumber cnTTEScaledGNumerator =
  235.                 org.drip.function.definition.CartesianComplexNumber.Scale (cnGNumerator, dbTimeToExpiry);

  237.             if (null == cnTTEScaledGNumerator) return null;

  239.             if (null == (cnC = org.drip.function.definition.CartesianComplexNumber.Add (cnTTEScaledGNumerator, cnC)))
  240.                 return null;

  242.             if (null == (cnC = org.drip.function.definition.CartesianComplexNumber.Scale (cnC, dblA * dblSigmaScaler)))
  243.                 return null;

  245.             if (null == (cnC = org.drip.function.definition.CartesianComplexNumber.Add (new
  246.                 org.drip.function.definition.CartesianComplexNumber (0., dblRiskFreeRate * dbTimeToExpiry * dblFreq),
  247.                     cnC)))
  248.                 return null;

  250.             org.drip.function.definition.CartesianComplexNumber cnF = org.drip.function.definition.CartesianComplexNumber.Scale (cnD,
  251.                 dblInitialVolatility);

  253.             if (null == cnF) return null;

  255.             if (null == (cnF = org.drip.function.definition.CartesianComplexNumber.Add (cnF, new
  256.                 org.drip.function.definition.CartesianComplexNumber (0., java.lang.Math.log (dblSpot) * dblFreq))))
  257.                 return null;

  259.             if (null == (cnF = org.drip.function.definition.CartesianComplexNumber.Add (cnF, cnC))) return null;

  261.             if (null == (cnF = org.drip.function.definition.CartesianComplexNumber.Add (cnF, new
  262.                 org.drip.function.definition.CartesianComplexNumber (0., -1. * java.lang.Math.log (dblStrike) * dblFreq))))
  263.                 return null;

  265.             if (null == (cnF = org.drip.function.definition.CartesianComplexNumber.Exponentiate (cnF))) return null;

  267.             if (null == (cnF = org.drip.function.definition.CartesianComplexNumber.Divide (cnF, new
  268.                 org.drip.function.definition.CartesianComplexNumber (0., dblFreq))))
  269.                 return null;

  271.             return new PhaseCorrectedF (dblCorrectedPhase, cnF);
  272.         } catch (java.lang.Exception e) {
  273.             e.printStackTrace();
  274.         }

  276.         return null;
  277.     }

  279.     private PhaseCorrectedF fourierTransformAMST07 (
  280.         final double dblStrike,
  281.         final double dbTimeToExpiry,
  282.         final double dblRiskFreeRate,
  283.         final double dblSpot,
  284.         final double dblInitialVolatility,
  285.         final double dblA,
  286.         final double dblFreq,
  287.         final double dblB,
  288.         final double dblU,
  289.         final org.drip.numerical.fourier.RotationCountPhaseTracker rcpt)
  290.     {
  291.         try {
  292.             org.drip.function.definition.CartesianComplexNumber cnSmallDLHS = new org.drip.function.definition.CartesianComplexNumber (dblB,
  293.                 -1. * _fphp.rho() * _fphp.sigma() * dblFreq);

  295.             org.drip.function.definition.CartesianComplexNumber cnSmallD = org.drip.function.definition.CartesianComplexNumber.Square
  296.                 (cnSmallDLHS);

  298.             if (null == cnSmallD) return null;

  300.             double dblSigmaScaler = _fphp.sigma() * _fphp.sigma();

  302.             if (null == (cnSmallD = org.drip.function.definition.CartesianComplexNumber.Add (cnSmallD, new
  303.                 org.drip.function.definition.CartesianComplexNumber (dblSigmaScaler * dblFreq * dblFreq, -2. * dblSigmaScaler
  304.                     * dblFreq * dblU))))
  305.                 return null;

  307.             if (null == (cnSmallD = org.drip.function.definition.CartesianComplexNumber.SquareRoot (cnSmallD))) return null;

  309.             org.drip.function.definition.CartesianComplexNumber cnGNumerator = org.drip.function.definition.CartesianComplexNumber.Add
  310.                 (cnSmallDLHS, cnSmallD);

  312.             if (null == cnGNumerator) return null;

  314.             org.drip.function.definition.CartesianComplexNumber cnG = org.drip.function.definition.CartesianComplexNumber.Subtract
  315.                 (cnSmallDLHS, cnSmallD);

  317.             if (null == cnG) return null;

  319.             if (null == (cnG = org.drip.function.definition.CartesianComplexNumber.Divide (cnGNumerator, cnG))) return null;

  321.             int iM = 0;
  322.             int iN = 0;

  324.             if (org.drip.numerical.fourier.PhaseAdjuster.MULTI_VALUE_BRANCH_POWER_PHASE_TRACKER_KAHL_JACKEL ==
  325.                 _fphp.phaseTrackerType()) {
  326.                 iM = (int) ((cnG.argument() + java.lang.Math.PI) / (2. * java.lang.Math.PI));

  328.                 iN = (int) ((cnG.argument() + (dbTimeToExpiry * cnSmallD.argument()) + java.lang.Math.PI) /
  329.                     (2. * java.lang.Math.PI));
  330.             }

  332.             org.drip.function.definition.CartesianComplexNumber cnExpTTEScaledSmallD =
  333.                 org.drip.function.definition.CartesianComplexNumber.Scale (cnSmallD, dbTimeToExpiry);

  335.             if (null == cnExpTTEScaledSmallD) return null;

  337.             if (null == (cnExpTTEScaledSmallD = org.drip.function.definition.CartesianComplexNumber.Exponentiate
  338.                 (cnExpTTEScaledSmallD)))
  339.                 return null;

  341.             org.drip.function.definition.CartesianComplexNumber cnD = new org.drip.function.definition.CartesianComplexNumber (1. -
  342.                 cnExpTTEScaledSmallD.real(), -1. * cnExpTTEScaledSmallD.imaginary());

  344.             org.drip.function.definition.CartesianComplexNumber cnInvGExpTTEScaledSmallD =
  345.                 org.drip.function.definition.CartesianComplexNumber.Multiply (cnG, cnExpTTEScaledSmallD);

  347.             if (null == cnInvGExpTTEScaledSmallD) return null;

  349.             cnInvGExpTTEScaledSmallD = new org.drip.function.definition.CartesianComplexNumber (1. -
  350.                 cnInvGExpTTEScaledSmallD.real(), -1. * cnInvGExpTTEScaledSmallD.imaginary());

  352.             if (null == (cnD = org.drip.function.definition.CartesianComplexNumber.Divide (cnD, cnInvGExpTTEScaledSmallD)))
  353.                 return null;

  355.             if (null == (cnD = org.drip.function.definition.CartesianComplexNumber.Multiply (cnGNumerator, cnD)))
  356.                 return null;

  358.             dblSigmaScaler = 1. / dblSigmaScaler;

  360.             if (null == (cnD = org.drip.function.definition.CartesianComplexNumber.Scale (cnD, dblSigmaScaler))) return null;

  362.             org.drip.function.definition.CartesianComplexNumber cnC = new org.drip.function.definition.CartesianComplexNumber (1. -
  363.                 cnG.real(), -1. * cnG.imaginary());

  365.             if (org.drip.numerical.fourier.PhaseAdjuster.MULTI_VALUE_BRANCH_POWER_PHASE_TRACKER_KAHL_JACKEL ==
  366.                 _fphp.phaseTrackerType()) {
  367.                 if (null == (cnC = org.drip.numerical.fourier.PhaseAdjuster.PowerLogPhaseTracker
  368.                     (cnInvGExpTTEScaledSmallD, cnC, iN, iM)))
  369.                     return null;
  370.             } else if (org.drip.numerical.fourier.PhaseAdjuster.MULTI_VALUE_BRANCH_PHASE_TRACKER_ROTATION_COUNT
  371.                 == _fphp.phaseTrackerType()) {
  372.                 if (null == (cnC = org.drip.function.definition.CartesianComplexNumber.Logarithm (cnC))) return null;

  374.                 cnC = new org.drip.function.definition.CartesianComplexNumber (cnC.real(), rcpt.updateAndApply
  375.                     (cnC.argument(), true));
  376.             }

  378.             double dblCorrectedPhase = cnC.argument();

  380.             if (null == (cnC = org.drip.function.definition.CartesianComplexNumber.Scale (cnC, -2.))) return null;

  382.             org.drip.function.definition.CartesianComplexNumber cnTTEScaledGNumerator =
  383.                 org.drip.function.definition.CartesianComplexNumber.Scale (cnGNumerator, dbTimeToExpiry);

  385.             if (null == cnTTEScaledGNumerator) return null;

  387.             if (null == (cnC = org.drip.function.definition.CartesianComplexNumber.Add (cnTTEScaledGNumerator, cnC)))
  388.                 return null;

  390.             if (null == (cnC = org.drip.function.definition.CartesianComplexNumber.Scale (cnC, dblA * dblSigmaScaler)))
  391.                 return null;

  393.             if (null == (cnC = org.drip.function.definition.CartesianComplexNumber.Add (new
  394.                 org.drip.function.definition.CartesianComplexNumber (0., dblRiskFreeRate * dbTimeToExpiry * dblFreq),
  395.                     cnC)))
  396.                 return null;

  398.             org.drip.function.definition.CartesianComplexNumber cnF = org.drip.function.definition.CartesianComplexNumber.Scale (cnD,
  399.                 dblInitialVolatility);

  401.             if (null == cnF) return null;

  403.             if (null == (cnF = org.drip.function.definition.CartesianComplexNumber.Add (cnF, new
  404.                 org.drip.function.definition.CartesianComplexNumber (0., java.lang.Math.log (dblSpot) * dblFreq))))
  405.                 return null;

  407.             if (null == (cnF = org.drip.function.definition.CartesianComplexNumber.Add (cnF, cnC))) return null;

  409.             if (null == (cnF = org.drip.function.definition.CartesianComplexNumber.Add (cnF, new
  410.                 org.drip.function.definition.CartesianComplexNumber (0., -1. * java.lang.Math.log (dblStrike) * dblFreq))))
  411.                 return null;

  413.             if (null == (cnF = org.drip.function.definition.CartesianComplexNumber.Exponentiate (cnF))) return null;

  415.             if (null == (cnF = org.drip.function.definition.CartesianComplexNumber.Divide (cnF, new
  416.                 org.drip.function.definition.CartesianComplexNumber (0., dblFreq))))
  417.                 return null;

  419.             return new PhaseCorrectedF (dblCorrectedPhase, cnF);
  420.         } catch (java.lang.Exception e) {
  421.             e.printStackTrace();
  422.         }

  424.         return null;
  425.     }

  427.     private PhaseCorrectedF payoffTransform (
  428.         final double dblStrike,
  429.         final double dbTimeToExpiry,
  430.         final double dblRiskFreeRate,
  431.         final double dblSpot,
  432.         final double dblInitialVolatility,
  433.         final double dblA,
  434.         final double dblFreq,
  435.         final double dblB,
  436.         final double dblU,
  437.         final org.drip.numerical.fourier.RotationCountPhaseTracker rcpt)
  438.     {
  439.         if (org.drip.numerical.fourier.PhaseAdjuster.MULTI_VALUE_BRANCH_PHASE_TRACKER_ROTATION_COUNT ==
  440.             _fphp.phaseTrackerType() && null == rcpt)
  441.         {
  442.             return null;
  443.         }

  445.         if (PAYOFF_TRANSFORM_SCHEME_HESTON_1993 == _fphp.payoffTransformScheme())
  446.             return fourierTransformHeston93 (dblStrike, dbTimeToExpiry, dblRiskFreeRate, dblSpot,
  447.                 dblInitialVolatility, dblA, dblFreq, dblB, dblU, rcpt);

  449.         if (PAYOFF_TRANSFORM_SCHEME_AMST_2007 == _fphp.payoffTransformScheme())
  450.             return fourierTransformAMST07 (dblStrike, dbTimeToExpiry, dblRiskFreeRate, dblSpot,
  451.                 dblInitialVolatility, dblA, dblFreq, dblB, dblU, rcpt);

  453.         return null;
  454.     }

  456.     /**
  457.      * HestonStochasticVolatilityAlgorithm constructor
  458.      * 
  459.      * @param fphp The Heston Algorithm Parameters
  460.      * 
  461.      * @throws java.lang.Exception Thrown if the Inputs are Invalid
  462.      */

  464.     public HestonStochasticVolatilityAlgorithm (
  465.         final org.drip.param.pricer.HestonOptionPricerParams fphp)
  466.         throws java.lang.Exception
  467.     {
  468.         if (null == (_fphp = fphp))
  469.             throw new java.lang.Exception ("HestonStochasticVolatilityAlgorithm ctr: Invalid Inputs");
  470.     }

  472.     /**
  473.      * Record the Details of a Single Phase Adjustment Run
  474.      * 
  475.      * @param dblStrike Strike
  476.      * @param dbTimeToExpiry TTE
  477.      * @param dblRiskFreeRate Risk Free Rate
  478.      * @param dblSpot Spot
  479.      * @param dblInitialVolatility Initial Volatility
  480.      * @param bLeft TRUE - Phase Correction applied to Left
  481.      * 
  482.      * @return Map of the Phase Correction Record
  483.      */

  485.     public java.util.Map<java.lang.Double, java.lang.Double> recordPhase (
  486.         final double dblStrike,
  487.         final double dbTimeToExpiry,
  488.         final double dblRiskFreeRate,
  489.         final double dblSpot,
  490.         final double dblInitialVolatility,
  491.         final boolean bLeft)
  492.     {
  493.         if (!org.drip.numerical.common.NumberUtil.IsValid (dblStrike) ||!org.drip.numerical.common.NumberUtil.IsValid
  494.             (dblSpot) ||!org.drip.numerical.common.NumberUtil.IsValid (dblInitialVolatility) ||
  495.                 !org.drip.numerical.common.NumberUtil.IsValid (dbTimeToExpiry) ||
  496.                     !org.drip.numerical.common.NumberUtil.IsValid (dblRiskFreeRate))
  497.             return null;

  499.         int i = 0;
  500.         double dblU1 = 0.5;
  501.         double dblU2 = -0.5;
  502.         double dblPreviousPhase = 0.;

  504.         org.drip.numerical.fourier.RotationCountPhaseTracker rcpt =
  505.             org.drip.numerical.fourier.PhaseAdjuster.MULTI_VALUE_BRANCH_PHASE_TRACKER_ROTATION_COUNT ==
  506.                 _fphp.phaseTrackerType() ? new org.drip.numerical.fourier.RotationCountPhaseTracker() : null;

  508.         double dblA = _fphp.kappa() * _fphp.theta();

  510.         double dblB2 = _fphp.kappa() + _fphp.lambda();

  512.         double dblB1 = dblB2 - _fphp.rho() * _fphp.sigma();

  514.         java.util.Map<java.lang.Double, java.lang.Double> mapPhaseRun = new
  515.             java.util.TreeMap<java.lang.Double, java.lang.Double>();

  517.         for (double dblFreq = FOURIER_FREQ_INIT; dblFreq <= FOURIER_FREQ_FINAL; dblFreq +=
  518.             FOURIER_FREQ_INCREMENT, ++i) {
  519.             PhaseCorrectedF pcf = bLeft ? payoffTransform (dblStrike, dbTimeToExpiry, dblRiskFreeRate,
  520.                 dblSpot, dblInitialVolatility, dblA, dblFreq, dblB1, dblU1, rcpt) : payoffTransform
  521.                     (dblStrike, dbTimeToExpiry, dblRiskFreeRate, dblSpot, dblInitialVolatility, dblA,
  522.                         dblFreq, dblB2, dblU2, rcpt);

  524.             if (null != rcpt) {
  525.                 if (0 == i)
  526.                     dblPreviousPhase = rcpt.getPreviousPhase();
  527.                 else if (1 == i) {
  528.                     double dblCurrentPhase = rcpt.getPreviousPhase();

  530.                     if (dblCurrentPhase < dblPreviousPhase) {
  531.                         if (!rcpt.setDirection
  532.                             (org.drip.numerical.fourier.RotationCountPhaseTracker.APPLY_BACKWARD))
  533.                             return null;
  534.                     } else if (dblCurrentPhase > dblPreviousPhase) {
  535.                         if (!rcpt.setDirection
  536.                             (org.drip.numerical.fourier.RotationCountPhaseTracker.APPLY_FORWARD))
  537.                             return null;
  538.                     } else
  539.                         return null;
  540.                 }
  541.             }

  543.             mapPhaseRun.put (dblFreq, pcf._dblCorrectedPhase);
  544.         }

  546.         return mapPhaseRun;
  547.     }

  549.     @Override public double payoff (
  550.         final double dblStrike,
  551.         final double dblTimeToExpiry,
  552.         final double dblRiskFreeRate,
  553.         final double dblUnderlier,
  554.         final boolean bIsPut,
  555.         final boolean bIsForward,
  556.         final double dblInitialVolatility,
  557.         final boolean bAsPrice)
  558.         throws java.lang.Exception
  559.     {
  560.         if (!org.drip.numerical.common.NumberUtil.IsValid (dblStrike) ||
  561.             !org.drip.numerical.common.NumberUtil.IsValid (dblUnderlier) ||
  562.                 !org.drip.numerical.common.NumberUtil.IsValid (dblInitialVolatility) ||
  563.                     !org.drip.numerical.common.NumberUtil.IsValid (dblTimeToExpiry) ||
  564.                         !org.drip.numerical.common.NumberUtil.IsValid (dblRiskFreeRate))
  565.             throw new java.lang.Exception ("HestonStochasticVolatilityAlgorithm::payoff => Invalid Inputs");

  567.         org.drip.numerical.fourier.RotationCountPhaseTracker rcpt1 =
  568.             org.drip.numerical.fourier.PhaseAdjuster.MULTI_VALUE_BRANCH_PHASE_TRACKER_ROTATION_COUNT ==
  569.                 _fphp.phaseTrackerType() ? new org.drip.numerical.fourier.RotationCountPhaseTracker() : null;

  571.         org.drip.numerical.fourier.RotationCountPhaseTracker rcpt2 =
  572.             org.drip.numerical.fourier.PhaseAdjuster.MULTI_VALUE_BRANCH_PHASE_TRACKER_ROTATION_COUNT ==
  573.                 _fphp.phaseTrackerType() ? new org.drip.numerical.fourier.RotationCountPhaseTracker() : null;

  575.         double dblA = _fphp.kappa() * _fphp.theta();

  577.         double dblB2 = _fphp.kappa() + _fphp.lambda();

  579.         double dblB1 = dblB2 - _fphp.rho() * _fphp.sigma();

  581.         double dblDF = java.lang.Math.exp (-1. * dblRiskFreeRate * dblTimeToExpiry);

  583.         int i = 0;
  584.         double dblU1 = 0.5;
  585.         double dblU2 = -0.5;
  586.         double dblCallProb1 = 0.;
  587.         double dblCallProb2 = 0.;
  588.         double dblPreviousPhase = 0.;
  589.         double dblSpot = bIsForward ? dblUnderlier * dblDF : dblUnderlier;

  591.         for (double dblFreq = FOURIER_FREQ_INIT; dblFreq <= FOURIER_FREQ_FINAL; dblFreq +=
  592.             FOURIER_FREQ_INCREMENT, ++i) {
  593.             PhaseCorrectedF pcf1 = payoffTransform (dblStrike, dblTimeToExpiry, dblRiskFreeRate, dblSpot,
  594.                 dblInitialVolatility, dblA, dblFreq, dblB1, dblU1, rcpt1);

  596.             if (null != rcpt1) {
  597.                 if (0 == i)
  598.                     dblPreviousPhase = rcpt1.getPreviousPhase();
  599.                 else if (1 == i) {
  600.                     double dblCurrentPhase = rcpt1.getPreviousPhase();

  602.                     if (dblCurrentPhase < dblPreviousPhase) {
  603.                         if (!rcpt1.setDirection
  604.                             (org.drip.numerical.fourier.RotationCountPhaseTracker.APPLY_BACKWARD))
  605.                             throw new java.lang.Exception
  606.                                 ("HestonStochasticVolatilityAlgorithm::payoff => Cannot compute payoff");
  607.                     } else if (dblCurrentPhase > dblPreviousPhase) {
  608.                         if (!rcpt1.setDirection
  609.                             (org.drip.numerical.fourier.RotationCountPhaseTracker.APPLY_FORWARD))
  610.                             throw new java.lang.Exception
  611.                                 ("HestonStochasticVolatilityAlgorithm::payoff => Cannot compute payoff");
  612.                     } else
  613.                         throw new java.lang.Exception
  614.                             ("HestonStochasticVolatilityAlgorithm::payoff => Cannot compute payoff");
  615.                 }
  616.             }

  618.             PhaseCorrectedF pcf2 = payoffTransform (dblStrike, dblTimeToExpiry, dblRiskFreeRate, dblSpot,
  619.                 dblInitialVolatility, dblA, dblFreq, dblB2, dblU2, rcpt2);

  621.             if (null != rcpt2) {
  622.                 if (0 == i)
  623.                     dblPreviousPhase = rcpt2.getPreviousPhase();
  624.                 else if (1 == i) {
  625.                     double dblCurrentPhase = rcpt2.getPreviousPhase();

  627.                     if (dblCurrentPhase < dblPreviousPhase) {
  628.                         if (!rcpt2.setDirection
  629.                             (org.drip.numerical.fourier.RotationCountPhaseTracker.APPLY_BACKWARD))
  630.                             throw new java.lang.Exception
  631.                                 ("HestonStochasticVolatilityAlgorithm::payoff => Cannot compute payoff");
  632.                     } else if (dblCurrentPhase > dblPreviousPhase) {
  633.                         if (!rcpt2.setDirection
  634.                             (org.drip.numerical.fourier.RotationCountPhaseTracker.APPLY_FORWARD))
  635.                             throw new java.lang.Exception
  636.                                 ("HestonStochasticVolatilityAlgorithm::payoff => Cannot compute payoff");
  637.                     } else
  638.                         throw new java.lang.Exception
  639.                             ("HestonStochasticVolatilityAlgorithm::payoff => Cannot compute payoff");
  640.                 }
  641.             }

  643.             dblCallProb1 += pcf1._cnF.real() * FOURIER_FREQ_INCREMENT;

  645.             dblCallProb2 += pcf2._cnF.real() * FOURIER_FREQ_INCREMENT;
  646.         }

  648.         double dblForward = dblSpot / dblDF;
  649.         double dblPIScaler = 1. / java.lang.Math.PI;
  650.         double dblCallPayoff = dblForward * (0.5 + dblCallProb1 * dblPIScaler) - dblStrike * (0.5 +
  651.             dblCallProb2 * dblPIScaler);

  653.         if (!bAsPrice) return bIsPut ? dblCallPayoff + dblStrike - dblForward : dblCallPayoff;

  655.         return bIsPut ? dblDF * (dblCallPayoff + dblStrike - dblForward) : dblDF * dblCallPayoff;
  656.     }

  658.     @Override public org.drip.pricer.option.Greeks greeks (
  659.         final double dblStrike,
  660.         final double dblTimeToExpiry,
  661.         final double dblRiskFreeRate,
  662.         final double dblUnderlier,
  663.         final boolean bIsPut,
  664.         final boolean bIsForward,
  665.         final double dblInitialVolatility)
  666.     {
  667.         if (!org.drip.numerical.common.NumberUtil.IsValid (dblStrike) ||
  668.             !org.drip.numerical.common.NumberUtil.IsValid (dblUnderlier) ||
  669.                 !org.drip.numerical.common.NumberUtil.IsValid (dblInitialVolatility) ||
  670.                     !org.drip.numerical.common.NumberUtil.IsValid (dblTimeToExpiry) ||
  671.                         !org.drip.numerical.common.NumberUtil.IsValid (dblRiskFreeRate))
  672.             return null;

  674.         org.drip.numerical.fourier.RotationCountPhaseTracker rcpt1 =
  675.             org.drip.numerical.fourier.PhaseAdjuster.MULTI_VALUE_BRANCH_PHASE_TRACKER_ROTATION_COUNT ==
  676.                 _fphp.phaseTrackerType() ? new org.drip.numerical.fourier.RotationCountPhaseTracker() : null;

  678.         org.drip.numerical.fourier.RotationCountPhaseTracker rcpt2 =
  679.             org.drip.numerical.fourier.PhaseAdjuster.MULTI_VALUE_BRANCH_PHASE_TRACKER_ROTATION_COUNT ==
  680.                 _fphp.phaseTrackerType() ? new org.drip.numerical.fourier.RotationCountPhaseTracker() : null;

  682.         double dblA = _fphp.kappa() * _fphp.theta();

  684.         double dblB2 = _fphp.kappa() + _fphp.lambda();

  686.         double dblB1 = dblB2 - _fphp.rho() * _fphp.sigma();

  688.         double dblDF = java.lang.Math.exp (-1. * dblRiskFreeRate * dblTimeToExpiry);

  690.         int i = 0;
  691.         double dblU1 = 0.5;
  692.         double dblU2 = -0.5;
  693.         double dblCallProb1 = 0.;
  694.         double dblCallProb2 = 0.;
  695.         double dblPreviousPhase = 0.;
  696.         double dblSpot = bIsForward ? dblUnderlier * dblDF : dblUnderlier;

  698.         for (double dblFreq = FOURIER_FREQ_INIT; dblFreq <= FOURIER_FREQ_FINAL; dblFreq +=
  699.             FOURIER_FREQ_INCREMENT, ++i) {
  700.             PhaseCorrectedF pcf1 = payoffTransform (dblStrike, dblTimeToExpiry, dblRiskFreeRate, dblSpot,
  701.                 dblInitialVolatility, dblA, dblFreq, dblB1, dblU1, rcpt1);

  703.             if (null != rcpt1) {
  704.                 if (0 == i)
  705.                     dblPreviousPhase = rcpt1.getPreviousPhase();
  706.                 else if (1 == i) {
  707.                     double dblCurrentPhase = rcpt1.getPreviousPhase();

  709.                     if (dblCurrentPhase < dblPreviousPhase) {
  710.                         if (!rcpt1.setDirection
  711.                             (org.drip.numerical.fourier.RotationCountPhaseTracker.APPLY_BACKWARD))
  712.                             return null;
  713.                     } else if (dblCurrentPhase > dblPreviousPhase) {
  714.                         if (!rcpt1.setDirection
  715.                             (org.drip.numerical.fourier.RotationCountPhaseTracker.APPLY_FORWARD))
  716.                             return null;
  717.                     } else
  718.                         return null;
  719.                 }
  720.             }

  722.             PhaseCorrectedF pcf2 = payoffTransform (dblStrike, dblTimeToExpiry, dblRiskFreeRate, dblSpot,
  723.                 dblInitialVolatility, dblA, dblFreq, dblB2, dblU2, rcpt2);

  725.             if (null != rcpt2) {
  726.                 if (0 == i)
  727.                     dblPreviousPhase = rcpt2.getPreviousPhase();
  728.                 else if (1 == i) {
  729.                     double dblCurrentPhase = rcpt2.getPreviousPhase();

  731.                     if (dblCurrentPhase < dblPreviousPhase) {
  732.                         if (!rcpt2.setDirection
  733.                             (org.drip.numerical.fourier.RotationCountPhaseTracker.APPLY_BACKWARD))
  734.                             return null;
  735.                     } else if (dblCurrentPhase > dblPreviousPhase) {
  736.                         if (!rcpt2.setDirection
  737.                             (org.drip.numerical.fourier.RotationCountPhaseTracker.APPLY_FORWARD))
  738.                             return null;
  739.                     } else
  740.                         return null;
  741.                 }
  742.             }

  744.             dblCallProb1 += pcf1._cnF.real() * FOURIER_FREQ_INCREMENT;

  746.             dblCallProb2 += pcf2._cnF.real() * FOURIER_FREQ_INCREMENT;
  747.         }

  749.         double dblForward = dblSpot / dblDF;
  750.         double dblPIScaler = 1. / java.lang.Math.PI;
  751.         dblCallProb1 = 0.5 + dblCallProb1 * dblPIScaler;
  752.         dblCallProb2 = 0.5 + dblCallProb2 * dblPIScaler;
  753.         double dblATMCallPayoff = dblForward * (dblCallProb1 - dblCallProb2);
  754.         double dblCallPrice = dblSpot * dblCallProb1 - dblStrike * dblDF * dblCallProb2;
  755.         double dblExpectedCallPayoff = dblForward * dblCallProb1 - dblStrike * dblDF * dblCallProb2;

  757.         try {
  758.             if (!bIsPut)
  759.                 return new org.drip.pricer.option.Greeks (
  760.                     dblDF,
  761.                     dblInitialVolatility,
  762.                     dblExpectedCallPayoff,
  763.                     dblATMCallPayoff,
  764.                     dblCallPrice,
  765.                     dblCallProb1,
  766.                     dblCallProb2,
  767.                     dblCallProb1,
  768.                     java.lang.Double.NaN,
  769.                     java.lang.Double.NaN,
  770.                     java.lang.Double.NaN,
  771.                     java.lang.Double.NaN,
  772.                     java.lang.Double.NaN,
  773.                     java.lang.Double.NaN,
  774.                     java.lang.Double.NaN,
  775.                     java.lang.Double.NaN,
  776.                     java.lang.Double.NaN,
  777.                     java.lang.Double.NaN,
  778.                     java.lang.Double.NaN
  779.                 );

  781.             double dblPutPriceFromParity = dblCallPrice + dblStrike * dblDF - dblSpot;

  783.             return new org.drip.pricer.option.PutGreeks (
  784.                 dblDF,
  785.                 dblInitialVolatility,
  786.                 dblExpectedCallPayoff + dblStrike - dblSpot,
  787.                 dblATMCallPayoff,
  788.                 dblPutPriceFromParity,
  789.                 dblPutPriceFromParity,
  790.                 java.lang.Double.NaN,
  791.                 java.lang.Double.NaN,
  792.                 java.lang.Double.NaN,
  793.                 java.lang.Double.NaN,
  794.                 java.lang.Double.NaN,
  795.                 java.lang.Double.NaN,
  796.                 java.lang.Double.NaN,
  797.                 java.lang.Double.NaN,
  798.                 java.lang.Double.NaN,
  799.                 java.lang.Double.NaN,
  800.                 java.lang.Double.NaN,
  801.                 java.lang.Double.NaN,
  802.                 java.lang.Double.NaN,
  803.                 java.lang.Double.NaN
  804.             );
  805.         } catch (java.lang.Exception e) {
  806.             e.printStackTrace();
  807.         }

  809.         return null;
  810.     }
  811. }
Created with JaCoCo 0.7.9.201702052155