MinimalQuadraticHaganWest.java

  1. package org.drip.spline.pchip;

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

  83. /**
  84.  * <i>MinimalQuadraticHaganWest</i> implements the regime using the Hagan and West (2006) Minimal Quadratic
  85.  * Estimator.
  86.  *
  87.  * <br><br>
  88.  *  <ul>
  89.  *      <li><b>Module </b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/ComputationalCore.md">Computational Core Module</a></li>
  90.  *      <li><b>Library</b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/SplineBuilderLibrary.md">Spline Builder Library</a></li>
  91.  *      <li><b>Project</b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/src/main/java/org/drip/spline/README.md">Basis Splines and Linear Compounders across a Broad Family of Spline Basis Functions</a></li>
  92.  *      <li><b>Package</b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/src/main/java/org/drip/spline/pchip/README.md">Monotone Convex Themed PCHIP Splines</a></li>
  93.  *  </ul>
  94.  * <br><br>
  95.  *
  96.  * @author Lakshmi Krishnamurthy
  97.  */

  99. public class MinimalQuadraticHaganWest {
  100.     private double[] _adblA = null;
  101.     private double[] _adblB = null;
  102.     private double[] _adblC = null;
  103.     private double[] _adblObservation = null;
  104.     private double[] _adblPredictorOrdinate = null;
  105.     private double _dblWeight = java.lang.Double.NaN;

  107.     /**
  108.      * Create an instance of MinimalQuadraticHaganWest
  109.      * 
  110.      * @param adblPredictorOrdinate Array of Predictor Ordinates
  111.      * @param adblObservation Array of Observations
  112.      * @param dblWeight Relative Weights applied across the first and the second derivatives
  113.      * 
  114.      * @return Instance of MinimalQuadraticHaganWest
  115.      */

  117.     public static final MinimalQuadraticHaganWest Create (
  118.         final double[] adblPredictorOrdinate,
  119.         final double[] adblObservation,
  120.         final double dblWeight)
  121.     {
  122.         MinimalQuadraticHaganWest mchw = null;

  124.         try {
  125.             mchw = new MinimalQuadraticHaganWest (adblPredictorOrdinate, adblObservation, dblWeight);
  126.         } catch (java.lang.Exception e) {
  127.             e.printStackTrace();

  129.             return null;
  130.         }

  132.         return mchw.setupCoefficients() ? mchw : null;
  133.     }

  135.     private MinimalQuadraticHaganWest (
  136.         final double[] adblPredictorOrdinate,
  137.         final double[] adblObservation,
  138.         final double dblWeight)
  139.         throws java.lang.Exception
  140.     {
  141.         if (null == (_adblObservation = adblObservation) || null == (_adblPredictorOrdinate =
  142.             adblPredictorOrdinate) || !org.drip.numerical.common.NumberUtil.IsValid (_dblWeight = dblWeight))
  143.             throw new java.lang.Exception ("MinimalQuadraticHaganWest ctr: Invalid Inputs!");

  145.         int iNumObservation = _adblObservation.length;

  147.         if (1 >= iNumObservation || iNumObservation + 1 != _adblPredictorOrdinate.length)
  148.             throw new java.lang.Exception ("MinimalQuadraticHaganWest ctr: Invalid Inputs!");
  149.     }

  151.     private boolean setupCoefficients()
  152.     {
  153.         int iNumObservation = _adblObservation.length;
  154.         _adblA = new double[iNumObservation];
  155.         _adblB = new double[iNumObservation];
  156.         _adblC = new double[iNumObservation];
  157.         double[] adblH = new double[iNumObservation];
  158.         double[] adblRHS = new double[3 * iNumObservation];
  159.         double[][] aadblCoeffMatrix = new double[3 * iNumObservation][3 * iNumObservation];

  161.         for (int i = 0; i < 3 * iNumObservation; ++i) {
  162.             adblRHS[i] = 0.;

  164.             for (int j = 0; j < 3 * iNumObservation; ++j)
  165.                 aadblCoeffMatrix[i][j] = 0.;
  166.         }

  168.         for (int i = 0; i < iNumObservation; ++i)
  169.             adblH[i] = _adblPredictorOrdinate[i + 1] - _adblPredictorOrdinate[i];

  171.         /*
  172.          * Setting up the coefficient linear constraint equation set
  173.          * 
  174.          *  - Left index => Equation Index
  175.          *  - Right Index => Coefficient Index
  176.          */

  178.         /*
  179.          * Set up the conserved quantities; Laid out as:
  180.          *      A_i + (H_i / 2.) * B_i + (H_i * H_i / 3.) * C_i = Observation_i
  181.          */

  183.         for (int iEq = 0; iEq < iNumObservation; ++iEq) {
  184.             int iSegmentIndex = iEq;
  185.             adblRHS[iEq] = _adblObservation[iEq]; // Z_i
  186.             aadblCoeffMatrix[iEq][3 * iSegmentIndex] = 1.; // A_i
  187.             aadblCoeffMatrix[iEq][3 * iSegmentIndex + 1] = 0.5 * adblH[iSegmentIndex]; // B_i
  188.             aadblCoeffMatrix[iEq][3 * iSegmentIndex + 2] = adblH[iSegmentIndex] * adblH[iSegmentIndex] / 3.; // C_i
  189.         }

  191.         /*
  192.          * Set up the continuity constraints; Laid out as:
  193.          *      A_i + H_i * B_i + (H_i * H_i) * C_i - A_i+1 = 0.
  194.          */

  196.         for (int iEq = iNumObservation; iEq < 2 * iNumObservation - 1; ++iEq) {
  197.             adblRHS[iEq] = 0.;
  198.             int iSegmentIndex = iEq - iNumObservation;
  199.             aadblCoeffMatrix[iEq][3 * iSegmentIndex] = 1.; // A_i
  200.             aadblCoeffMatrix[iEq][3 * iSegmentIndex + 1] = adblH[iSegmentIndex]; // B_i
  201.             aadblCoeffMatrix[iEq][3 * iSegmentIndex + 2] = adblH[iSegmentIndex] * adblH[iSegmentIndex]; // C_i
  202.             aadblCoeffMatrix[iEq][3 * iSegmentIndex + 3] = -1.; // A_i+1
  203.         }

  205.         /*
  206.          * Set up the derivative penalty minimizer; Laid out as:
  207.          *      w * B_i + (2. * H_i) * C_i - w * B_i+1 = 0.
  208.          */

  210.         for (int iEq = 2 * iNumObservation - 1; iEq < 3 * iNumObservation - 2; ++iEq) {
  211.             adblRHS[iEq] = 0.;
  212.             int iSegmentIndex = iEq - 2 * iNumObservation + 1;
  213.             aadblCoeffMatrix[iEq][3 * iSegmentIndex + 1] = _dblWeight; // B_i
  214.             aadblCoeffMatrix[iEq][3 * iSegmentIndex + 2] = 2. * adblH[iSegmentIndex]; // C_i
  215.             aadblCoeffMatrix[iEq][3 * iSegmentIndex + 4] = -1. * _dblWeight; // B_i+1
  216.         }

  218.         /*
  219.          * Left Boundary Condition: Starting Left Slope is zero, i.e., B_0 = 0.
  220.          */

  222.         adblRHS[3 * iNumObservation - 2] = 0.;
  223.         aadblCoeffMatrix[3 * iNumObservation - 2][1] = 1.;

  225.         /*
  226.          * Right Boundary Condition: Final First Derivative is zero, i.e., B_n-1 = 0.
  227.          */

  229.         adblRHS[3 * iNumObservation - 1] = 0.;
  230.         aadblCoeffMatrix[3 * iNumObservation - 1][3 * iNumObservation - 2] = 1.;

  232.         org.drip.numerical.linearalgebra.LinearizationOutput lssGaussianElimination =
  233.             org.drip.numerical.linearalgebra.LinearSystemSolver.SolveUsingGaussianElimination (aadblCoeffMatrix,
  234.                 adblRHS);

  236.         if (null == lssGaussianElimination) return false;

  238.         double[] adblCoeff = lssGaussianElimination.getTransformedRHS();

  240.         if (null == adblCoeff || 3 * iNumObservation != adblCoeff.length) return false;

  242.         int iSegment = 0;

  244.         for (int i = 0; i < 3 * iNumObservation; ++i) {
  245.             if (0 == i % 3)
  246.                 _adblA[iSegment] = adblCoeff[i];
  247.             else if (1 == i % 3)
  248.                 _adblB[iSegment] = adblCoeff[i];
  249.             else if (2 == i % 3) {
  250.                 _adblC[iSegment] = adblCoeff[i];
  251.                 ++iSegment;
  252.             }
  253.         }

  255.         return true;
  256.     }

  258.     private int containingIndex (
  259.         final double dblPredictorOrdinate,
  260.         final boolean bIncludeLeft,
  261.         final boolean bIncludeRight)
  262.         throws java.lang.Exception
  263.     {
  264.         int iNumSegment = _adblA.length;

  266.         for (int i = 0 ; i < iNumSegment; ++i) {
  267.             boolean bLeftValid = bIncludeLeft ? _adblPredictorOrdinate[i] <= dblPredictorOrdinate :
  268.                 _adblPredictorOrdinate[i] < dblPredictorOrdinate;

  270.             boolean bRightValid = bIncludeRight ? _adblPredictorOrdinate[i + 1] >= dblPredictorOrdinate :
  271.                 _adblPredictorOrdinate[i + 1] > dblPredictorOrdinate;

  273.             if (bLeftValid && bRightValid) return i;
  274.         }

  276.         throw new java.lang.Exception
  277.             ("MinimalQuadraticHaganWest::containingIndex => Cannot locate Containing Index");
  278.     }

  280.     /**
  281.      * Calculate the Response Value given the Predictor Ordinate
  282.      * 
  283.      * @param dblPredictorOrdinate The Predictor Ordinate
  284.      * 
  285.      * @return The Response Value
  286.      * 
  287.      * @throws java.lang.Exception Thrown if the input is invalid
  288.      */

  290.     public double responseValue (
  291.         final double dblPredictorOrdinate)
  292.         throws java.lang.Exception
  293.     {
  294.         int i = containingIndex (dblPredictorOrdinate, true, true);

  296.         return _adblA[i] + _adblB[i] * (dblPredictorOrdinate - _adblPredictorOrdinate[i]) + _adblC[i] *
  297.             (dblPredictorOrdinate - _adblPredictorOrdinate[i]) * (dblPredictorOrdinate -
  298.                 _adblPredictorOrdinate[i]);
  299.     }

  301.     public double[] calcConservedConstraint()
  302.     {
  303.         int iNumObservation = _adblObservation.length;
  304.         double[] adblConservedConstraint = new double[iNumObservation];

  306.         for (int i = 0; i < iNumObservation; ++i)
  307.             adblConservedConstraint[i] = _adblA[i] + _adblB[i] * 0.5 + _adblC[i] / 3.;

  309.         return adblConservedConstraint;
  310.     }

  312.     public static final void main (
  313.         final String[] astrArgs)
  314.         throws Exception
  315.     {
  316.         double[] adblTime = new double[] {0., 1.0, 2.0};
  317.         double[] adblForwardRate = new double[] {0.02, 0.026};

  319.         MinimalQuadraticHaganWest mqhw = MinimalQuadraticHaganWest.Create (adblTime, adblForwardRate, 0.5);

  321.         double[] adblConservedConstraint = mqhw.calcConservedConstraint();

  323.         for (int i = 0; i < adblConservedConstraint.length; ++i)
  324.             System.out.println ("Conserved Constraint[" + i + "] => " +
  325.                 org.drip.numerical.common.FormatUtil.FormatDouble (adblConservedConstraint[i], 1, 6, 1.));

  327.         for (double dblTime = adblTime[0]; dblTime <= adblTime[adblTime.length - 1]; dblTime += 0.25)
  328.             System.out.println ("Response[" + org.drip.numerical.common.FormatUtil.FormatDouble (dblTime, 2, 2,
  329.                 1.) + "] = " + org.drip.numerical.common.FormatUtil.FormatDouble (mqhw.responseValue (dblTime), 1,
  330.                     6, 1.));
  331.     }
  332. }
Created with JaCoCo 0.7.9.201702052155