NonLinearLeastSquareWithPenalty { public LeastSquareWithPenaltyResults solve(final DoubleMatrix1D observedValues, final Function1D<DoubleMatrix1D, DoubleMatrix1D> func, final DoubleMatrix1D startPos, final DoubleMatrix2D penalty) { final int n = observedValues.getNumberOfElements(); final VectorFieldFirstOrderDifferentiator jac = new VectorFieldFirstOrderDifferentiator(); return solve(observedValues, new DoubleMatrix1D(n, 1.0), func, jac.differentiate(func), startPos, penalty); } NonLinearLeastSquareWithPenalty();  NonLinearLeastSquareWithPenalty(final Decomposition<?> decomposition);  NonLinearLeastSquareWithPenalty(final double eps);  NonLinearLeastSquareWithPenalty(final Decomposition<?> decomposition, final double eps);  NonLinearLeastSquareWithPenalty(final Decomposition<?> decomposition, final MatrixAlgebra algebra, final double eps); LeastSquareWithPenaltyResults solve(final DoubleMatrix1D observedValues, final Function1D<DoubleMatrix1D, DoubleMatrix1D> func, final DoubleMatrix1D startPos, final DoubleMatrix2D penalty); LeastSquareWithPenaltyResults solve(final DoubleMatrix1D observedValues, final DoubleMatrix1D sigma, final Function1D<DoubleMatrix1D, DoubleMatrix1D> func, final DoubleMatrix1D startPos,
      final DoubleMatrix2D penalty); LeastSquareWithPenaltyResults solve(final DoubleMatrix1D observedValues, final DoubleMatrix1D sigma, final Function1D<DoubleMatrix1D, DoubleMatrix1D> func, final DoubleMatrix1D startPos,
      final DoubleMatrix2D penalty, final Function1D<DoubleMatrix1D, Boolean> allowedValue); LeastSquareWithPenaltyResults solve(final DoubleMatrix1D observedValues, final DoubleMatrix1D sigma, final Function1D<DoubleMatrix1D, DoubleMatrix1D> func,
      final Function1D<DoubleMatrix1D, DoubleMatrix2D> jac, final DoubleMatrix1D startPos, final DoubleMatrix2D penalty); LeastSquareWithPenaltyResults solve(final DoubleMatrix1D observedValues, final DoubleMatrix1D sigma, final Function1D<DoubleMatrix1D, DoubleMatrix1D> func,
      final Function1D<DoubleMatrix1D, DoubleMatrix2D> jac, final DoubleMatrix1D startPos, final DoubleMatrix2D penalty, final Function1D<DoubleMatrix1D, Boolean> allowedValue); static final Function1D<DoubleMatrix1D, Boolean> UNCONSTRAINED; }
@Test public void linearTest() { final boolean print = false; if (print) { System.out.println("NonLinearLeastSquareWithPenaltyTest.linearTest"); } final int nWeights = 20; final int diffOrder = 2; final double lambda = 100.0; final DoubleMatrix2D penalty = (DoubleMatrix2D) MA.scale(getPenaltyMatrix(nWeights, diffOrder), lambda); final int[] onIndex = new int[] {1, 4, 11, 12, 15, 17 }; final double[] obs = new double[] {0, 1.0, 1.0, 1.0, 0.0, 0.0 }; final int n = onIndex.length; final Function1D<DoubleMatrix1D, DoubleMatrix1D> func = new Function1D<DoubleMatrix1D, DoubleMatrix1D>() { @Override public DoubleMatrix1D evaluate(final DoubleMatrix1D x) { final double[] temp = new double[n]; for (int i = 0; i < n; i++) { temp[i] = x.getEntry(onIndex[i]); } return new DoubleMatrix1D(temp); } }; final Function1D<DoubleMatrix1D, DoubleMatrix2D> jac = new Function1D<DoubleMatrix1D, DoubleMatrix2D>() { @Override public DoubleMatrix2D evaluate(final DoubleMatrix1D x) { final DoubleMatrix2D res = new DoubleMatrix2D(n, nWeights); for (int i = 0; i < n; i++) { res.getData()[i][onIndex[i]] = 1.0; } return res; } }; final RandomEngine ran = new MersenneTwister64(MersenneTwister.DEFAULT_SEED); final double[] temp = new double[nWeights]; for (int i = 0; i < nWeights; i++) { temp[i] = ran.nextDouble(); } final DoubleMatrix1D start = new DoubleMatrix1D(temp); final LeastSquareWithPenaltyResults lsRes = NLLSWP.solve(new DoubleMatrix1D(obs), new DoubleMatrix1D(n, 0.01), func, jac, start, penalty); if (print) { System.out.println("chi2: " + lsRes.getChiSq()); System.out.println(lsRes.getFitParameters()); } for (int i = 0; i < n; i++) { assertEquals(obs[i], lsRes.getFitParameters().getEntry(onIndex[i]), 0.01); } double expPen = 20.87912357454752; assertEquals(expPen, lsRes.getPenalty(), 1e-9); }
@Test public void volTermStructureTest() { double lowStrike = 0.001; double highStrike = 0.05; final double fwd = 0.023; final double t = 2.0; double alpha = 0.2; double beta = 0.9; double rho = -0.2; double nu = 0.5; SABRFormulaData sabr = new SABRFormulaData(alpha, beta, rho, nu); final double[] strikes = new double[] {0.002, 0.004, 0.006, 0.01, 0.02, 0.03, 0.05 }; SABRHaganVolatilityFunction volFunc = new SABRHaganVolatilityFunction(); double[] vols = volFunc.getVolatilityFunction(fwd, strikes, t).evaluate(sabr); final int nStrikes = strikes.length; double[] prices = new double[nStrikes]; double[] vega = new double[nStrikes]; for (int i = 0; i < nStrikes; i++) { prices[i] = BlackFormulaRepository.price(fwd, strikes[i], t, vols[i], true); vega[i] = BlackFormulaRepository.vega(fwd, strikes[i], t, vols[i]); } BasisFunctionGenerator gen = new BasisFunctionGenerator(); final List<Function1D<Double, Double>> basisFuncs = gen.generateSet(lowStrike, highStrike, 12, 2); final int nWeights = basisFuncs.size(); double lambda = 10000.0; DoubleMatrix2D p = (DoubleMatrix2D) MA.scale(getPenaltyMatrix(nWeights, 2), lambda); Function1D<DoubleMatrix1D, DoubleMatrix1D> priceFunc = new Function1D<DoubleMatrix1D, DoubleMatrix1D>() { @Override public DoubleMatrix1D evaluate(DoubleMatrix1D w) { BasisFunctionAggregation<Double> agg = new BasisFunctionAggregation<>(basisFuncs, w.getData()); DoubleMatrix1D res = new DoubleMatrix1D(nStrikes); double[] data = res.getData(); for (int i = 0; i < nStrikes; i++) { data[i] = BlackFormulaRepository.price(fwd, strikes[i], t, agg.evaluate(strikes[i]), true); } return res; } }; Function1D<DoubleMatrix1D, DoubleMatrix2D> priceJac = new Function1D<DoubleMatrix1D, DoubleMatrix2D>() { @Override public DoubleMatrix2D evaluate(DoubleMatrix1D w) { BasisFunctionAggregation<Double> agg = new BasisFunctionAggregation<>(basisFuncs, w.getData()); DoubleMatrix2D res = new DoubleMatrix2D(nStrikes, nWeights); double[][] data = res.getData(); for (int i = 0; i < nStrikes; i++) { double vega = BlackFormulaRepository.vega(fwd, strikes[i], t, agg.evaluate(strikes[i])); DoubleMatrix1D dSigmaDW = agg.weightSensitivity(strikes[i]); double[] temp = dSigmaDW.toArray(); for (int j = 0; j < nWeights; j++) { temp[j] *= vega; } data[i] = temp; } return res; } }; double oneBP = 1e-4; DoubleMatrix1D sigma = new DoubleMatrix1D(nStrikes); for (int i = 0; i < nStrikes; i++) { sigma.getData()[i] = vega[i] * 0.1 * oneBP; } DoubleMatrix1D start = new DoubleMatrix1D(nWeights, 0.3); LeastSquareWithPenaltyResults res = NLLSWP.solve(new DoubleMatrix1D(prices), sigma, priceFunc, priceJac, start, p); DoubleMatrix1D weights = res.getFitParameters(); BasisFunctionAggregation<Double> agg = new BasisFunctionAggregation<>(basisFuncs, weights.getData()); for (int i = 0; i < nStrikes; i++) { assertEquals(vols[i], agg.evaluate(strikes[i]), 1e-5); } double expChi2 = 0.06671526449702014; assertEquals(expChi2, res.getChiSq(), 1e-8); NonLinearLeastSquareWithPenalty nllswp = new NonLinearLeastSquareWithPenalty(1e-6); res = nllswp.solve(new DoubleMatrix1D(prices), sigma, priceFunc, priceJac, start, p); assertEquals(expChi2, res.getChiSq(), 1e-8); nllswp = new NonLinearLeastSquareWithPenalty(DecompositionFactory.LU_COMMONS); res = nllswp.solve(new DoubleMatrix1D(prices), sigma, priceFunc, priceJac, start, p); assertEquals(expChi2, res.getChiSq(), 2e-8); nllswp = new NonLinearLeastSquareWithPenalty(DecompositionFactory.QR_COMMONS); res = nllswp.solve(new DoubleMatrix1D(prices), sigma, priceFunc, priceJac, start, p); assertEquals(expChi2, res.getChiSq(), 2e-8); nllswp = new NonLinearLeastSquareWithPenalty(new CholeskyDecompositionCommons()); res = nllswp.solve(new DoubleMatrix1D(prices), sigma, priceFunc, priceJac, start, p); assertEquals(expChi2, res.getChiSq(), 2e-8); nllswp = new NonLinearLeastSquareWithPenalty(new CholeskyDecompositionCommons(), 1e-4); assertEquals(expChi2, res.getChiSq(), 2e-8); }
@Test public void test() { final boolean print = false; if (print) { System.out.println("NonLinearLeastSquareWithPenaltyTest"); } final int nWeights = B_SPLINES.size(); final LeastSquareResults res = NLLSWP.solve(new DoubleMatrix1D(RATES), new DoubleMatrix1D(RATES.length, 1e-4), WEIGHTS_TO_SWAP_FUNC, new DoubleMatrix1D(nWeights, 0.03), PENALTY_MAT); if (print) { System.out.println("chi2: " + res.getChiSq()); System.out.println(); } final DoubleMatrix1D fittedSwaps = WEIGHTS_TO_SWAP_FUNC.evaluate(res.getFitParameters()); for (int i = 0; i < N_SWAPS; i++) { if (print) { System.out.println("swap rates: " + RATES[i] + "\t" + fittedSwaps.getEntry(i)); } assertEquals(RATES[i], fittedSwaps.getEntry(i), 1e-4); } if (print) { final Function1D<Double, Double> func = new BasisFunctionAggregation<>(B_SPLINES, res.getFitParameters().getData()); System.out.println(); System.out.println("t\t yield"); for (int i = 0; i < 101; i++) { final double t = i * 20.0 / 100; System.out.println(t + "\t" + func.evaluate(t)); } } }