MinimalQuadraticHaganWest.java

package org.drip.spline.pchip;

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

/*!
 * Copyright (C) 2020 Lakshmi Krishnamurthy
 * Copyright (C) 2019 Lakshmi Krishnamurthy
 * Copyright (C) 2018 Lakshmi Krishnamurthy
 * Copyright (C) 2017 Lakshmi Krishnamurthy
 * Copyright (C) 2016 Lakshmi Krishnamurthy
 * Copyright (C) 2015 Lakshmi Krishnamurthy
 * Copyright (C) 2014 Lakshmi Krishnamurthy
 * Copyright (C) 2013 Lakshmi Krishnamurthy
 * 
 *  This file is part of DROP, an open-source library targeting analytics/risk, transaction cost analytics,
 *  	asset liability management analytics, capital, exposure, and margin analytics, valuation adjustment
 *  	analytics, and portfolio construction analytics within and across fixed income, credit, commodity,
 *  	equity, FX, and structured products. It also includes auxiliary libraries for algorithm support,
 *  	numerical analysis, numerical optimization, spline builder, model validation, statistical learning,
 *  	and computational support.
 *  
 *  	https://lakshmidrip.github.io/DROP/
 *  
 *  DROP is composed of three modules:
 *  
 *  - DROP Product Core - https://lakshmidrip.github.io/DROP-Product-Core/
 *  - DROP Portfolio Core - https://lakshmidrip.github.io/DROP-Portfolio-Core/
 *  - DROP Computational Core - https://lakshmidrip.github.io/DROP-Computational-Core/
 * 
 * 	DROP Product Core implements libraries for the following:
 * 	- Fixed Income Analytics
 * 	- Loan Analytics
 * 	- Transaction Cost Analytics
 * 
 * 	DROP Portfolio Core implements libraries for the following:
 * 	- Asset Allocation Analytics
 *  - Asset Liability Management Analytics
 * 	- Capital Estimation Analytics
 * 	- Exposure Analytics
 * 	- Margin Analytics
 * 	- XVA Analytics
 * 
 * 	DROP Computational Core implements libraries for the following:
 * 	- Algorithm Support
 * 	- Computation Support
 * 	- Function Analysis
 *  - Model Validation
 * 	- Numerical Analysis
 * 	- Numerical Optimizer
 * 	- Spline Builder
 *  - Statistical Learning
 * 
 * 	Documentation for DROP is Spread Over:
 * 
 * 	- Main                     => https://lakshmidrip.github.io/DROP/
 * 	- Wiki                     => https://github.com/lakshmiDRIP/DROP/wiki
 * 	- GitHub                   => https://github.com/lakshmiDRIP/DROP
 * 	- Repo Layout Taxonomy     => https://github.com/lakshmiDRIP/DROP/blob/master/Taxonomy.md
 * 	- Javadoc                  => https://lakshmidrip.github.io/DROP/Javadoc/index.html
 * 	- Technical Specifications => https://github.com/lakshmiDRIP/DROP/tree/master/Docs/Internal
 * 	- Release Versions         => https://lakshmidrip.github.io/DROP/version.html
 * 	- Community Credits        => https://lakshmidrip.github.io/DROP/credits.html
 * 	- Issues Catalog           => https://github.com/lakshmiDRIP/DROP/issues
 * 	- JUnit                    => https://lakshmidrip.github.io/DROP/junit/index.html
 * 	- Jacoco                   => https://lakshmidrip.github.io/DROP/jacoco/index.html
 * 
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *   	you may not use this file except in compliance with the License.
 *   
 *  You may obtain a copy of the License at
 *  	http://www.apache.org/licenses/LICENSE-2.0
 *  
 *  Unless required by applicable law or agreed to in writing, software
 *  	distributed under the License is distributed on an "AS IS" BASIS,
 *  	WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  
 *  See the License for the specific language governing permissions and
 *  	limitations under the License.
 */

/**
 * <i>MinimalQuadraticHaganWest</i> implements the regime using the Hagan and West (2006) Minimal Quadratic
 * Estimator.
 *
 * <br><br>
 *  <ul>
 *		<li><b>Module </b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/ComputationalCore.md">Computational Core Module</a></li>
 *		<li><b>Library</b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/SplineBuilderLibrary.md">Spline Builder Library</a></li>
 *		<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>
 *		<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>
 *  </ul>
 * <br><br>
 *
 * @author Lakshmi Krishnamurthy
 */

public class MinimalQuadraticHaganWest {
	private double[] _adblA = null;
	private double[] _adblB = null;
	private double[] _adblC = null;
	private double[] _adblObservation = null;
	private double[] _adblPredictorOrdinate = null;
	private double _dblWeight = java.lang.Double.NaN;

	/**
	 * Create an instance of MinimalQuadraticHaganWest
	 * 
	 * @param adblPredictorOrdinate Array of Predictor Ordinates
	 * @param adblObservation Array of Observations
	 * @param dblWeight Relative Weights applied across the first and the second derivatives
	 * 
	 * @return Instance of MinimalQuadraticHaganWest
	 */

	public static final MinimalQuadraticHaganWest Create (
		final double[] adblPredictorOrdinate,
		final double[] adblObservation,
		final double dblWeight)
	{
		MinimalQuadraticHaganWest mchw = null;

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

			return null;
		}

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

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

		int iNumObservation = _adblObservation.length;

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

	private boolean setupCoefficients()
	{
		int iNumObservation = _adblObservation.length;
		_adblA = new double[iNumObservation];
		_adblB = new double[iNumObservation];
		_adblC = new double[iNumObservation];
		double[] adblH = new double[iNumObservation];
		double[] adblRHS = new double[3 * iNumObservation];
		double[][] aadblCoeffMatrix = new double[3 * iNumObservation][3 * iNumObservation];

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

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

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

		/*
		 * Setting up the coefficient linear constraint equation set
		 * 
		 * 	- Left index => Equation Index
		 * 	- Right Index => Coefficient Index
		 */

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

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

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

		for (int iEq = iNumObservation; iEq < 2 * iNumObservation - 1; ++iEq) {
			adblRHS[iEq] = 0.;
			int iSegmentIndex = iEq - iNumObservation;
			aadblCoeffMatrix[iEq][3 * iSegmentIndex] = 1.; // A_i
			aadblCoeffMatrix[iEq][3 * iSegmentIndex + 1] = adblH[iSegmentIndex]; // B_i
			aadblCoeffMatrix[iEq][3 * iSegmentIndex + 2] = adblH[iSegmentIndex] * adblH[iSegmentIndex]; // C_i
			aadblCoeffMatrix[iEq][3 * iSegmentIndex + 3] = -1.; // A_i+1
		}

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

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

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

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

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

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

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

		if (null == lssGaussianElimination) return false;

		double[] adblCoeff = lssGaussianElimination.getTransformedRHS();

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

		int iSegment = 0;

		for (int i = 0; i < 3 * iNumObservation; ++i) {
			if (0 == i % 3)
				_adblA[iSegment] = adblCoeff[i];
			else if (1 == i % 3)
				_adblB[iSegment] = adblCoeff[i];
			else if (2 == i % 3) {
				_adblC[iSegment] = adblCoeff[i];
				++iSegment;
			}
		}

		return true;
	}

	private int containingIndex (
		final double dblPredictorOrdinate,
		final boolean bIncludeLeft,
		final boolean bIncludeRight)
		throws java.lang.Exception
	{
		int iNumSegment = _adblA.length;

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

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

			if (bLeftValid && bRightValid) return i;
		}

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

	/**
	 * Calculate the Response Value given the Predictor Ordinate
	 * 
	 * @param dblPredictorOrdinate The Predictor Ordinate
	 * 
	 * @return The Response Value
	 * 
	 * @throws java.lang.Exception Thrown if the input is invalid
	 */

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

		return _adblA[i] + _adblB[i] * (dblPredictorOrdinate - _adblPredictorOrdinate[i]) + _adblC[i] *
			(dblPredictorOrdinate - _adblPredictorOrdinate[i]) * (dblPredictorOrdinate -
				_adblPredictorOrdinate[i]);
	}

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

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

		return adblConservedConstraint;
	}

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

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

		double[] adblConservedConstraint = mqhw.calcConservedConstraint();

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

		for (double dblTime = adblTime[0]; dblTime <= adblTime[adblTime.length - 1]; dblTime += 0.25)
			System.out.println ("Response[" + org.drip.numerical.common.FormatUtil.FormatDouble (dblTime, 2, 2,
				1.) + "] = " + org.drip.numerical.common.FormatUtil.FormatDouble (mqhw.responseValue (dblTime), 1,
					6, 1.));
	}
}