R1NonCentral.java

  1. package org.drip.measure.chisquare;

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

  77. /**
  78.  * <i>R1NonCentral</i> implements the Distribution Table for the R<sup>1</sup> Non-central Chi-Square
  79.  *  Distribution. The References are:
  80.  * 
  81.  * <br><br>
  82.  *  <ul>
  83.  *      <li>
  84.  *          Johnson, N. L., S. Kotz, and N. Balakrishnan (1995): <i>Continuous Univariate Distributions
  85.  *              2<sup>nd</sup> Edition</i> <b>John Wiley and Sons</b>
  86.  *      </li>
  87.  *      <li>
  88.  *          Muirhead, R. (2005): <i>Aspects of Multivariate Statistical Theory 2<sup>nd</sup> Edition</i>
  89.  *              <b>Wiley</b>
  90.  *      </li>
  91.  *      <li>
  92.  *          Non-central Chi-Squared Distribution (2019): Chi-Squared Function
  93.  *              https://en.wikipedia.org/wiki/Noncentral_chi-squared_distribution
  94.  *      </li>
  95.  *      <li>
  96.  *          Sankaran, M. (1963): Approximations to the Non-Central Chi-Square Distribution <i>Biometrika</i>
  97.  *              <b>50 (1-2)</b> 199-204
  98.  *      </li>
  99.  *      <li>
  100.  *          Young, D. S. (2010): tolerance: An R Package for Estimating Tolerance Intervals <i>Journal of
  101.  *              Statistical Software</i> <b>36 (5)</b> 1-39
  102.  *      </li>
  103.  *  </ul>
  104.  *
  105.  *  <br><br>
  106.  *  <ul>
  107.  *      <li><b>Module </b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/ComputationalCore.md">Computational Core Module</a></li>
  108.  *      <li><b>Library</b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/NumericalAnalysisLibrary.md">Numerical Analysis Library</a></li>
  109.  *      <li><b>Project</b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/src/main/java/org/drip/measure/README.md">R<sup>d</sup> Continuous/Discrete Probability Measures</a></li>
  110.  *      <li><b>Package</b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/src/main/java/org/drip/measure/chisquare/README.md">Chi-Square Distribution Implementation/Properties</a></li>
  111.  *  </ul>
  112.  *
  113.  * @author Lakshmi Krishnamurthy
  114.  */

  116. public class R1NonCentral
  117.     extends org.drip.measure.continuous.R1Univariate
  118. {
  119.     private double _cdfScaler = java.lang.Double.NaN;
  120.     private org.drip.function.definition.R1ToR1 _gammaEstimator = null;
  121.     private org.drip.function.definition.R1ToR1 _digammaEstimator = null;
  122.     private org.drip.function.definition.R2ToR1 _lowerIncompleteGammaEstimator = null;
  123.     private org.drip.measure.chisquare.R1NonCentralParameters _r1NonCentralParameters = null;
  124.     private org.drip.specialfunction.definition.ModifiedBesselFirstKindEstimator
  125.         _modifiedBesselFirstKindEstimator = null;

  127.     private double nonCentralMoment (
  128.         final int n,
  129.         final double[] fourLeadingRawMoments)
  130.         throws java.lang.Exception
  131.     {
  132.         if (n <= 3)
  133.         {
  134.             return fourLeadingRawMoments[n - 1];
  135.         }

  137.         double nonCentralMoment = cumulant (
  138.             n
  139.         );

  141.         double degreesOfFreedom = _r1NonCentralParameters.degreesOfFreedom();

  143.         double nonCentralityParameter = _r1NonCentralParameters.nonCentralityParameter();

  145.         for (int j = 1;
  146.             j < n;
  147.             ++j)
  148.         {
  149.             nonCentralMoment = nonCentralMoment +
  150.                 _gammaEstimator.evaluate (
  151.                     n
  152.                 ) * java.lang.Math.pow (
  153.                     2.,
  154.                     j - 1.
  155.                 ) * (
  156.                     degreesOfFreedom + j * nonCentralityParameter
  157.                 ) * nonCentralMoment (
  158.                     n - j,
  159.                     fourLeadingRawMoments
  160.                 ) / _gammaEstimator.evaluate (
  161.                     n - j + 1
  162.                 );
  163.         }

  165.         return nonCentralMoment;
  166.     }

  168.     /**
  169.      * Construct the Standard Instance of R1NonCentral
  170.      * 
  171.      * @param degreesOfFreedom Degrees of Freedom
  172.      * @param nonCentralityParameter Non-centrality Parameter
  173.      * @param gammaEstimator Gamma Estimator
  174.      * @param digammaEstimator Digamma Estimator
  175.      * @param lowerIncompleteGammaEstimator Lower Incomplete Gamma Estimator
  176.      * @param modifiedBesselFirstKindEstimator Modified Bessel First Kind Estimator
  177.      * 
  178.      * @return The Standard Instance of R1NonCentral
  179.      */

  181.     public static final R1NonCentral Standard (
  182.         final double degreesOfFreedom,
  183.         final double nonCentralityParameter,
  184.         final org.drip.function.definition.R1ToR1 gammaEstimator,
  185.         final org.drip.function.definition.R1ToR1 digammaEstimator,
  186.         final org.drip.function.definition.R2ToR1 lowerIncompleteGammaEstimator,
  187.         final org.drip.specialfunction.definition.ModifiedBesselFirstKindEstimator
  188.             modifiedBesselFirstKindEstimator)
  189.     {
  190.         try
  191.         {
  192.             return new R1NonCentral (
  193.                 new org.drip.measure.chisquare.R1NonCentralParameters (
  194.                     degreesOfFreedom,
  195.                     nonCentralityParameter
  196.                 ),
  197.                 gammaEstimator,
  198.                 digammaEstimator,
  199.                 lowerIncompleteGammaEstimator,
  200.                 modifiedBesselFirstKindEstimator
  201.             );
  202.         }
  203.         catch (java.lang.Exception e)
  204.         {
  205.             e.printStackTrace();
  206.         }

  208.         return null;
  209.     }

  211.     /**
  212.      * R1NonCentral Constructor
  213.      * 
  214.      * @param r1NonCentralParameters R<sup>1</sup> Non-central Parameters
  215.      * @param gammaEstimator Gamma Estimator
  216.      * @param digammaEstimator Digamma Estimator
  217.      * @param lowerIncompleteGammaEstimator Lower Incomplete Gamma Estimator
  218.      * @param modifiedBesselFirstKindEstimator Modified Bessel First Kind Estimator
  219.      * 
  220.      * @throws java.lang.Exception Thrown if the Inputs are Invalid
  221.      */

  223.     public R1NonCentral (
  224.         final org.drip.measure.chisquare.R1NonCentralParameters r1NonCentralParameters,
  225.         final org.drip.function.definition.R1ToR1 gammaEstimator,
  226.         final org.drip.function.definition.R1ToR1 digammaEstimator,
  227.         final org.drip.function.definition.R2ToR1 lowerIncompleteGammaEstimator,
  228.         final org.drip.specialfunction.definition.ModifiedBesselFirstKindEstimator
  229.             modifiedBesselFirstKindEstimator)
  230.         throws java.lang.Exception
  231.     {
  232.         if (null == (_r1NonCentralParameters = r1NonCentralParameters) ||
  233.             null == (_gammaEstimator = gammaEstimator) ||
  234.             null == (_digammaEstimator = digammaEstimator) ||
  235.             null == (_lowerIncompleteGammaEstimator = lowerIncompleteGammaEstimator) ||
  236.             null == (_modifiedBesselFirstKindEstimator = modifiedBesselFirstKindEstimator)
  237.         )
  238.         {
  239.             throw new java.lang.Exception (
  240.                 "R1NonCentral Constructor => Invalid Inputs"
  241.             );
  242.         }

  244.         _cdfScaler = java.lang.Math.exp (
  245.             -0.5 * _r1NonCentralParameters.nonCentralityParameter()
  246.         );
  247.     }

  249.     /**
  250.      * Retrieve the R<sup>1</sup> Non-Central Parameters
  251.      * 
  252.      * @return The R<sup>1</sup> Non-Central Parameters
  253.      */

  255.     public org.drip.measure.chisquare.R1NonCentralParameters parameters()
  256.     {
  257.         return _r1NonCentralParameters;
  258.     }

  260.     /**
  261.      * Retrieve the Gamma Estimator
  262.      * 
  263.      * @return Gamma Estimator
  264.      */

  266.     public org.drip.function.definition.R1ToR1 gammaEstimator()
  267.     {
  268.         return _gammaEstimator;
  269.     }

  271.     /**
  272.      * Retrieve the Digamma Estimator
  273.      * 
  274.      * @return Digamma Estimator
  275.      */

  277.     public org.drip.function.definition.R1ToR1 digammaEstimator()
  278.     {
  279.         return _digammaEstimator;
  280.     }

  282.     /**
  283.      * Retrieve the Lower Incomplete Gamma Estimator
  284.      * 
  285.      * @return Lower Incomplete Gamma Estimator
  286.      */

  288.     public org.drip.function.definition.R2ToR1 lowerIncompleteGammaEstimator()
  289.     {
  290.         return _lowerIncompleteGammaEstimator;
  291.     }

  293.     /**
  294.      * Retrieve the Modified Bessel First Kind Estimator
  295.      * 
  296.      * @return Modified Bessel First Kind Estimator
  297.      */

  299.     public org.drip.specialfunction.definition.ModifiedBesselFirstKindEstimator
  300.         modifiedBesselFirstKindEstimator()
  301.     {
  302.         return _modifiedBesselFirstKindEstimator;
  303.     }

  305.     @Override public double[] support()
  306.     {
  307.         return new double[]
  308.         {
  309.             0.,
  310.             java.lang.Double.POSITIVE_INFINITY
  311.         };
  312.     }

  314.     @Override public double density (
  315.         final double x)
  316.         throws java.lang.Exception
  317.     {
  318.         if (!supported (
  319.             x
  320.         ))
  321.         {
  322.             throw new java.lang.Exception (
  323.                 "R1NonCentral::density => Variate not in Range"
  324.             );
  325.         }

  327.         double degreesOfFreedom = _r1NonCentralParameters.degreesOfFreedom();

  329.         double nonCentralityParameter = _r1NonCentralParameters.nonCentralityParameter();

  331.         return 0.5 * java.lang.Math.pow (
  332.             x / nonCentralityParameter,
  333.             0.25 * degreesOfFreedom - 0.5
  334.         ) * java.lang.Math.exp (
  335.             -0.5 * (x + nonCentralityParameter)
  336.         ) * _modifiedBesselFirstKindEstimator.bigI (
  337.             0.5 * degreesOfFreedom - 1.,
  338.             java.lang.Math.sqrt (
  339.                 x * nonCentralityParameter
  340.             )
  341.         );
  342.     }

  344.     @Override public double cumulative (
  345.         final double t)
  346.         throws java.lang.Exception
  347.     {
  348.         if (!supported (
  349.             t
  350.         ))
  351.         {
  352.             throw new java.lang.Exception (
  353.                 "R1NonCentral::cumulative => Invalid Inputs"
  354.             );
  355.         }

  357.         int termCount = 10;
  358.         double cumulative = 0.;

  360.         double degreesOfFreedom = _r1NonCentralParameters.degreesOfFreedom();

  362.         double nonCentralityParameter = _r1NonCentralParameters.nonCentralityParameter();

  364.         for (int termIndex = 0;
  365.             termIndex < termCount;
  366.             ++termIndex)
  367.         {
  368.             cumulative = cumulative +
  369.                 java.lang.Math.pow (
  370.                     0.5 * nonCentralityParameter,
  371.                     termIndex
  372.                 ) * new org.drip.measure.chisquare.R1Central (
  373.                     degreesOfFreedom + 2 * termIndex,
  374.                     _gammaEstimator,
  375.                     _digammaEstimator,
  376.                     _lowerIncompleteGammaEstimator
  377.                 ).cumulative (
  378.                     t
  379.                 ) / org.drip.numerical.common.NumberUtil.Factorial (
  380.                     termIndex
  381.                 );
  382.         }

  384.         return _cdfScaler * cumulative;
  385.     }

  387.     @Override public double mean()
  388.         throws java.lang.Exception
  389.     {
  390.         return _r1NonCentralParameters.degreesOfFreedom() + _r1NonCentralParameters.nonCentralityParameter();
  391.     }

  393.     @Override public double variance()
  394.         throws java.lang.Exception
  395.     {
  396.         return 2. * _r1NonCentralParameters.degreesOfFreedom() +
  397.             4. * _r1NonCentralParameters.nonCentralityParameter();
  398.     }

  400.     @Override public double skewness()
  401.         throws java.lang.Exception
  402.     {
  403.         double degreesOfFreedom = _r1NonCentralParameters.degreesOfFreedom();

  405.         double nonCentralityParameter = _r1NonCentralParameters.nonCentralityParameter();

  407.         return (degreesOfFreedom + 3. * nonCentralityParameter) * java.lang.Math.pow (
  408.             2.,
  409.             degreesOfFreedom + 2. * nonCentralityParameter
  410.         );
  411.     }

  413.     @Override public double excessKurtosis()
  414.         throws java.lang.Exception
  415.     {
  416.         double degreesOfFreedom = _r1NonCentralParameters.degreesOfFreedom();

  418.         double nonCentralityParameter = _r1NonCentralParameters.nonCentralityParameter();

  420.         double kPlusTwoLambda = degreesOfFreedom + 2. * nonCentralityParameter;
  421.         return 12. * (degreesOfFreedom + 4. * nonCentralityParameter) /
  422.             (kPlusTwoLambda* kPlusTwoLambda);
  423.     }

  425.     @Override public org.drip.function.definition.R1ToR1 momentGeneratingFunction()
  426.     {
  427.         return new org.drip.function.definition.R1ToR1 (
  428.             null
  429.         )
  430.         {
  431.             @Override public double evaluate (
  432.                 final double t)
  433.                 throws java.lang.Exception
  434.             {
  435.                 if (!org.drip.numerical.common.NumberUtil.IsValid (
  436.                         t
  437.                     ) || t > 0.5
  438.                 )
  439.                 {
  440.                     throw new java.lang.Exception (
  441.                         "R1NonCentral::momentGeneratingFunction::evaluate => Invalid Input"
  442.                     );
  443.                 }

  445.                 double oneMinusTwoT = 1. - 2. * t;

  447.                 return java.lang.Math.exp (
  448.                     _r1NonCentralParameters.nonCentralityParameter() * t / oneMinusTwoT
  449.                 ) / java.lang.Math.pow (
  450.                     oneMinusTwoT,
  451.                     0.5 * _r1NonCentralParameters.degreesOfFreedom()
  452.                 );
  453.             }
  454.         };
  455.     }

  457.     /**
  458.      * Compute the Cumulant
  459.      * 
  460.      * @param n Cumulant Index
  461.      * 
  462.      * @return The Cumulant
  463.      * 
  464.      * @throws java.lang.Exception Thrown if the Inputs are Invalid
  465.      */

  467.     public double cumulant (
  468.         final int n)
  469.         throws java.lang.Exception
  470.     {
  471.         if (0 > n)
  472.         {
  473.             throw new java.lang.Exception (
  474.                 "R1NonCentral::cumulant => Invalid Inputs"
  475.             );
  476.         }

  478.         return (
  479.             _r1NonCentralParameters.degreesOfFreedom() + n * _r1NonCentralParameters.nonCentralityParameter()
  480.         ) * java.lang.Math.pow (
  481.             2.,
  482.             n - 1.
  483.         ) * _gammaEstimator.evaluate (
  484.             n
  485.         );
  486.     }

  488.     /**
  489.      * Compute the Leading Non-central Moments
  490.      * 
  491.      * @return Leading Non-central Moments
  492.      */

  494.     public double[] leadingRawMoments()
  495.     {
  496.         double degreesOfFreedom = _r1NonCentralParameters.degreesOfFreedom();

  498.         double nonCentralityParameter = _r1NonCentralParameters.nonCentralityParameter();

  500.         double[] fourLeadingRawMomentArray = new double[4];
  501.         fourLeadingRawMomentArray[0] = nonCentralityParameter + degreesOfFreedom;
  502.         fourLeadingRawMomentArray[1] =
  503.             fourLeadingRawMomentArray[0] * fourLeadingRawMomentArray[0] +
  504.             2. * (nonCentralityParameter + 2. * degreesOfFreedom);
  505.         fourLeadingRawMomentArray[2] =
  506.             fourLeadingRawMomentArray[0] * fourLeadingRawMomentArray[0] * fourLeadingRawMomentArray[0] +
  507.             6. * fourLeadingRawMomentArray[0] * (nonCentralityParameter + 2. * degreesOfFreedom) +
  508.             8. * (nonCentralityParameter + 3. * degreesOfFreedom);
  509.         fourLeadingRawMomentArray[3] =
  510.             fourLeadingRawMomentArray[0] * fourLeadingRawMomentArray[0] * fourLeadingRawMomentArray[0] *
  511.                 fourLeadingRawMomentArray[0] +
  512.             12. * fourLeadingRawMomentArray[0] * fourLeadingRawMomentArray[0] *
  513.                 (nonCentralityParameter + 2. * degreesOfFreedom) +
  514.             4. * (
  515.                 11. * degreesOfFreedom * degreesOfFreedom +
  516.                 44. * degreesOfFreedom * nonCentralityParameter +
  517.                 36. * nonCentralityParameter * nonCentralityParameter
  518.             ) +
  519.             48. * (nonCentralityParameter + 4. * degreesOfFreedom);
  520.         return fourLeadingRawMomentArray;
  521.     }

  523.     /**
  524.      * Compute the Leading central Moments
  525.      * 
  526.      * @return Leading central Moments
  527.      */

  529.     public double[] leadingCentralMoments()
  530.     {
  531.         double degreesOfFreedom = _r1NonCentralParameters.degreesOfFreedom();

  533.         double nonCentralityParameter = _r1NonCentralParameters.nonCentralityParameter();

  535.         double[] fourLeadingCentralMomentArray = new double[4];
  536.         fourLeadingCentralMomentArray[0] = 0.;
  537.         fourLeadingCentralMomentArray[1] =
  538.             2. * (nonCentralityParameter + 2. * degreesOfFreedom);
  539.         fourLeadingCentralMomentArray[2] =
  540.             8. * (nonCentralityParameter + 3. * degreesOfFreedom);
  541.         fourLeadingCentralMomentArray[3] =
  542.             12. * (nonCentralityParameter + 2. * degreesOfFreedom) *
  543.                 (nonCentralityParameter + 2. * degreesOfFreedom) +
  544.             48. * (nonCentralityParameter + 4. * degreesOfFreedom);
  545.         return fourLeadingCentralMomentArray;
  546.     }

  548.     /**
  549.      * Compute the Non-central Moment
  550.      * 
  551.      * @param n Moment Index
  552.      * 
  553.      * @return The Non-central Moment
  554.      * 
  555.      * @throws java.lang.Exception Thrown if the Inputs are Invalid
  556.      */

  558.     public double nonCentralMoment (
  559.         final int n)
  560.         throws java.lang.Exception
  561.     {
  562.         if (0 >= n)
  563.         {
  564.             throw new java.lang.Exception (
  565.                 "R1NonCentral::nonCentralMoment => Invalid Inputs"
  566.             );
  567.         }

  569.         return nonCentralMoment (
  570.             n,
  571.             leadingRawMoments()
  572.         );
  573.     }
  574. }
Created with JaCoCo 0.7.9.201702052155