From 0ce8640d4f03e547bac84ee521063e079bfa5eba Mon Sep 17 00:00:00 2001 From: Giovanni La Cagnina Date: Tue, 7 May 2024 17:13:25 +0200 Subject: [PATCH 1/2] [RD-PRICING] using hData.Strikes instead of Strike --- Heston/HestonOptimizationProblem.cs | 766 ++++++++++++++-------------- 1 file changed, 386 insertions(+), 380 deletions(-) diff --git a/Heston/HestonOptimizationProblem.cs b/Heston/HestonOptimizationProblem.cs index 0b341aa..24184c4 100644 --- a/Heston/HestonOptimizationProblem.cs +++ b/Heston/HestonOptimizationProblem.cs @@ -1,91 +1,91 @@ -/* Copyright (C) 2009-2012 Fairmat SRL (info@fairmat.com, http://www.fairmat.com/) - * Author(s): Enrico Degiuli (enrico.degiuli@fairmat.com) - * Matteo Tesser (matteo.tesser@fairmat.com) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU Lesser General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU Lesser General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . +/* Copyright (C) 2009-2012 Fairmat SRL (info@fairmat.com, http://www.fairmat.com/) + * Author(s): Enrico Degiuli (enrico.degiuli@fairmat.com) + * Matteo Tesser (matteo.tesser@fairmat.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . */ -/* - * Notes - * http://www.zeliade.com/whitepapers/zwp-0004.pdf +/* + * Notes + * http://www.zeliade.com/whitepapers/zwp-0004.pdf */ -using System; -using System.Collections.Generic; -using System.Threading.Tasks; +using System; +using System.Collections.Generic; +using System.Threading.Tasks; using DVPLI; -using DVPLDOM; -using DVPLI.TaskScheduler; -using Fairmat.MarketData; +using DVPLDOM; +using DVPLI.TaskScheduler; +using Fairmat.MarketData; using Fairmat.Optimization; - -namespace HestonEstimator -{ - - enum EWeighting + +namespace HestonEstimator +{ + + enum EWeighting { CONSTANT = 0, LINEAR = 1, LOG = 2 - } - - - /// - /// Describes the time-dependent Heston optimization problem. - /// - public class HestonCallOptimizationProblem : IOptimizationProblem + } + + + /// + /// Describes the time-dependent Heston optimization problem. + /// + public class HestonCallOptimizationProblem : IOptimizationProblem { - protected CallPriceMarketData cpmd; - /// - /// The call market price matrix. - /// - internal Matrix callMarketPrice; - - /// - /// The maturity vector relative to the callMarketPrice matrix. - /// - private Vector maturity; - - /// - /// The strike vector relative to callMarketPrice matrix. - /// - protected Vector strike; - - /// - /// The average rate vector up to a given maturity - /// - protected Vector rate; - - /// - /// The avereage dividend yield vector up to a given maturity - /// - protected Vector dividendYield; - - /// - /// Averegate drift up to a given maturity - /// - //private Vector drift; - - /// - /// Establish whether to use the boundary penalty function or not. - /// - public bool useBoundPenalty = false; - - /// - /// Establish whether to use the Feller penalty function or not. - /// Note: it's affect a lot the multi-level single linkage algorithm performance. + protected CallPriceMarketData cpmd; + /// + /// The call market price matrix. + /// + internal Matrix callMarketPrice; + + /// + /// The maturity vector relative to the callMarketPrice matrix. + /// + private Vector maturity; + + /// + /// The strike vector relative to callMarketPrice matrix. + /// + protected Vector strike; + + /// + /// The average rate vector up to a given maturity + /// + protected Vector rate; + + /// + /// The avereage dividend yield vector up to a given maturity + /// + protected Vector dividendYield; + + /// + /// Averegate drift up to a given maturity + /// + //private Vector drift; + + /// + /// Establish whether to use the boundary penalty function or not. + /// + public bool useBoundPenalty = false; + + /// + /// Establish whether to use the Feller penalty function or not. + /// Note: it's affect a lot the multi-level single linkage algorithm performance. /// public bool useFellerPenalty = false; @@ -104,134 +104,140 @@ public class HestonCallOptimizationProblem : IOptimizationProblem /// Builds objective function on relative pricing error. /// static bool optimizeRelativeError = false; - static double pricingMin = 0.0001; + static double pricingMin = 0.0001; static internal bool displayObjInfo = false; static internal double optionThreshold = 0.000; static internal EWeighting weighting = EWeighting.LOG; internal static bool calibrateOnCallOptions = true; internal static bool calibrateOnPutOptions = true; - - /// - /// Small value used in the boundary penalty function. - /// - protected double smallValue = 1e-4; - - /// - /// Value that weights the boundary penalty function. - /// - protected double k1 = 1e10; - - /// - /// Value that weights the Feller inequality penalty function. - /// - protected double k2 = 1e2; - - /// - /// The number of call option on which calibration is performed. - /// + + /// + /// Small value used in the boundary penalty function. + /// + protected double smallValue = 1e-4; + + /// + /// Value that weights the boundary penalty function. + /// + protected double k1 = 1e10; + + /// + /// Value that weights the Feller inequality penalty function. + /// + protected double k2 = 1e2; + + /// + /// The number of call option on which calibration is performed. + /// protected internal int numCall; protected internal int numPut; /// /// Total volume for calls and puts. /// - protected internal double totalVolume; - /// - /// Process starting value. - /// + protected internal double totalVolume; + /// + /// Process starting value. + /// protected double s0; protected Vector matBound; protected Vector strikeBound; - /// - /// Initializes a new instance of the HestonCallOptimizationProblem class using the - /// EquityCalibrationData data structure. - /// - /// - /// An EquityCalibrationData object containing market data for calibration. - /// - /// - /// A vector containing the minimum and maximum values - /// for maturities to be used in calibration. - /// - /// - /// A vector containing the minimum and maximum values - /// for strikes to be used in calibration. - /// - public HestonCallOptimizationProblem(EquityCalibrationData equityCalData, Vector matBound, Vector strikeBound) + /// + /// Initializes a new instance of the HestonCallOptimizationProblem class using the + /// EquityCalibrationData data structure. + /// + /// + /// An EquityCalibrationData object containing market data for calibration. + /// + /// + /// A vector containing the minimum and maximum values + /// for maturities to be used in calibration. + /// + /// + /// A vector containing the minimum and maximum values + /// for strikes to be used in calibration. + /// + public HestonCallOptimizationProblem(EquityCalibrationData equityCalData, Vector matBound, Vector strikeBound) { this.cpmd = equityCalData.Hdata; this.matBound = matBound; - this.strikeBound = strikeBound; - SetVariables(equityCalData.Hdata.CallPrice, equityCalData.Hdata.Maturity, - equityCalData.Hdata.Strike, equityCalData.CallMatrixRiskFreeRate, - equityCalData.CallMatrixDividendYield, equityCalData.Hdata.S0); + this.strikeBound = strikeBound; + var strike = equityCalData.Hdata.Strikes.GetRowReference(0).Clone(); + + SetVariables( + callMarketPrice: equityCalData.Hdata.CallPrice, + maturity: equityCalData.Hdata.Maturity, + strike: strike, + rate: equityCalData.CallMatrixRiskFreeRate, + s0: equityCalData.Hdata.S0, + dividendYield: equityCalData.CallMatrixDividendYield); - displayObjInfo = false; - } - - /// - /// Initializes a new instance of the HestonCallOptimizationProblem class. - /// - /// A matrix containing call option market prices. - /// - /// Vector of call option maturities relative to callMarketPrice matrix. - /// - /// - /// Vector of call option strikes relative to callMarketPrice matrix. - /// - /// - /// Vector of zero coupon bond rates calculated relative to maturity vector maturities. - /// - /// - /// Vector of dividend yield rates calculated relative to maturity vector maturities. - /// - /// Index/Equity value at the time of calibration. - /// - /// A vector containing the minimum and maximum values - /// for maturities to be used in calibration. - /// - /// - /// A vector containing the minimum and maximum values - /// for strikes to be used in calibration. - /// - [Obsolete] - HestonCallOptimizationProblem(Matrix callMarketPrice, Vector maturity, Vector strike, Vector rate, Vector dividendYield, double s0, Vector matBound, Vector strikeBound) - { + displayObjInfo = false; + } + + /// + /// Initializes a new instance of the HestonCallOptimizationProblem class. + /// + /// A matrix containing call option market prices. + /// + /// Vector of call option maturities relative to callMarketPrice matrix. + /// + /// + /// Vector of call option strikes relative to callMarketPrice matrix. + /// + /// + /// Vector of zero coupon bond rates calculated relative to maturity vector maturities. + /// + /// + /// Vector of dividend yield rates calculated relative to maturity vector maturities. + /// + /// Index/Equity value at the time of calibration. + /// + /// A vector containing the minimum and maximum values + /// for maturities to be used in calibration. + /// + /// + /// A vector containing the minimum and maximum values + /// for strikes to be used in calibration. + /// + [Obsolete] + HestonCallOptimizationProblem(Matrix callMarketPrice, Vector maturity, Vector strike, Vector rate, Vector dividendYield, double s0, Vector matBound, Vector strikeBound) + { this.matBound=matBound; - this.strikeBound = strikeBound; - SetVariables(callMarketPrice, maturity, strike, rate, - dividendYield, s0); - } - - /// - /// Sets several variables used to solve the optimization problem. - /// - /// A matrix containing call option market prices. - /// - /// Vector of call option maturities relative to callMarketPrice matrix. - /// - /// - /// Vector of call option strikes relative to callMarketPrice matrix. - /// - /// - /// Vector of zero coupon bond rates calculated relative to maturity vector maturities. - /// - /// - /// Vector of dividend yield rates calculated relative to maturity vector maturities. - /// - /// Index/Equity value at the time of calibration. - /// - /// A vector containing the minimum and maximum values - /// for maturities to be used in calibration. - /// - /// - /// A vector containing the minimum and maximum values - /// for strikes to be used in calibration. - private void SetVariables(Matrix callMarketPrice, Vector maturity, Vector strike, Vector rate, Vector dividendYield, double s0) - { + this.strikeBound = strikeBound; + SetVariables(callMarketPrice, maturity, strike, rate, + dividendYield, s0); + } + + /// + /// Sets several variables used to solve the optimization problem. + /// + /// A matrix containing call option market prices. + /// + /// Vector of call option maturities relative to callMarketPrice matrix. + /// + /// + /// Vector of call option strikes relative to callMarketPrice matrix. + /// + /// + /// Vector of zero coupon bond rates calculated relative to maturity vector maturities. + /// + /// + /// Vector of dividend yield rates calculated relative to maturity vector maturities. + /// + /// Index/Equity value at the time of calibration. + /// + /// A vector containing the minimum and maximum values + /// for maturities to be used in calibration. + /// + /// + /// A vector containing the minimum and maximum values + /// for strikes to be used in calibration. + private void SetVariables(Matrix callMarketPrice, Vector maturity, Vector strike, Vector rate, Vector dividendYield, double s0) + { this.s0 = s0; var dy = new PFunction(maturity, dividendYield); @@ -247,7 +253,7 @@ private void SetVariables(Matrix callMarketPrice, Vector maturity, Vector strike this.rate = rate; this.maturity = maturity; this.strike = strike; - this.callMarketPrice = callMarketPrice; + this.callMarketPrice = callMarketPrice; this.numCall = 0; callWeight = new Matrix(this.callMarketPrice.R, this.callMarketPrice.C); @@ -299,61 +305,61 @@ private void SetVariables(Matrix callMarketPrice, Vector maturity, Vector strike Console.WriteLine("Strike Bounds:\t" + strikeBound); Console.WriteLine("Maturity Bounds:\t" + matBound); Console.WriteLine("Lb:\t" + Bounds.Lb); - Console.WriteLine("Ub:\t" + Bounds.Ub); + Console.WriteLine("Ub:\t" + Bounds.Ub); if(Engine.Verbose>=2) - PutCallTest(); + PutCallTest(); + } + + + /// + /// Finds the couple of integer {i,j} so that the values + /// vector[i], vector[i+1],..., vector[j-1], vector[j] + /// are all included in the interval (bound[0], bound[1]) + /// while vector[i-1] and vector[j+1] are not. + /// + /// + /// Vector to be filtered. + /// + /// + /// Bounds to be used for the filtering. + /// + /// + /// Integer array representing the index couple. + /// + private int[] FindExtremes(Vector vector, Vector bounds) + { + int[] result = new int[2]; + int i = 0; + + while (vector[i] < bounds[0]) + i++; + result[0] = i; + i = vector.Length - 1; + while (vector[i] > bounds[1]) + i--; + result[1] = i; + return result; } - - /// - /// Finds the couple of integer {i,j} so that the values - /// vector[i], vector[i+1],..., vector[j-1], vector[j] - /// are all included in the interval (bound[0], bound[1]) - /// while vector[i-1] and vector[j+1] are not. - /// - /// - /// Vector to be filtered. - /// - /// - /// Bounds to be used for the filtering. - /// - /// - /// Integer array representing the index couple. - /// - private int[] FindExtremes(Vector vector, Vector bounds) - { - int[] result = new int[2]; - int i = 0; - - while (vector[i] < bounds[0]) - i++; - result[0] = i; - i = vector.Length - 1; - while (vector[i] > bounds[1]) - i--; - result[1] = i; - return result; - } - - #region IOptimizationProblem Members - + #region IOptimizationProblem Members + /// /// If volatility is observed, it may be useful to use this information - /// + /// protected double v0Min=0.001; protected double v0Max = 1; - protected double vLastMin = 0.001; - - /// - /// Gets the bounds for the optimization. - /// - public Bounds Bounds - { - get - { - Bounds bounds = new Bounds(); - - // The order of parameters is k, theta, sigma, rho, V0 (and optionally div.y) + protected double vLastMin = 0.001; + + /// + /// Gets the bounds for the optimization. + /// + public Bounds Bounds + { + get + { + Bounds bounds = new Bounds(); + + // The order of parameters is k, theta, sigma, rho, V0 (and optionally div.y) if (HestonConstantDriftEstimator.impliedDividends) { bounds.Lb = (Vector)new double[] { 0, v0Min, 0.001, -1, vLastMin, 0 }; @@ -363,74 +369,74 @@ public Bounds Bounds { bounds.Lb = (Vector)new double[] { 0, v0Min, 0.001, -1, vLastMin }; bounds.Ub = (Vector)new double[] { 15,v0Max, 2, 1, 1 }; - } - return bounds; - } - } - - /// - /// This method is unused but part of the interface. - /// - /// The parameter is not used. - /// Null as it's unused. - public DVPLI.Vector G(DVPLI.Vector x) - { - throw new NotImplementedException(); - } - - /// - /// This method is unused but part of the interface. - /// - /// The parameter is not used. - /// Nothing the function always throws a NotImplementedException. - public DVPLI.Vector Grad(DVPLI.Vector x) - { - throw new NotImplementedException(); - } - - /// - /// Gets a value indicating whether there are non linear constraints in this - /// optimization problem. In this case there are not. - /// - public bool HasNonLinearConstraints - { - get - { - return false; - } - } - - /// - /// Gets null as we have no linear constrains defined. - /// - public virtual LinearConstraints LinearIneqConstraints - { - get - { - return null; - } - } - - /// - /// Calibration objective function. - /// - /// - /// The vector of parameters. - /// - /// - /// Objective function value. - /// - public virtual double Obj(DVPLI.Vector x) - { + } + return bounds; + } + } + + /// + /// This method is unused but part of the interface. + /// + /// The parameter is not used. + /// Null as it's unused. + public DVPLI.Vector G(DVPLI.Vector x) + { + throw new NotImplementedException(); + } + + /// + /// This method is unused but part of the interface. + /// + /// The parameter is not used. + /// Nothing the function always throws a NotImplementedException. + public DVPLI.Vector Grad(DVPLI.Vector x) + { + throw new NotImplementedException(); + } + + /// + /// Gets a value indicating whether there are non linear constraints in this + /// optimization problem. In this case there are not. + /// + public bool HasNonLinearConstraints + { + get + { + return false; + } + } + + /// + /// Gets null as we have no linear constrains defined. + /// + public virtual LinearConstraints LinearIneqConstraints + { + get + { + return null; + } + } + + /// + /// Calibration objective function. + /// + /// + /// The vector of parameters. + /// + /// + /// Objective function value. + /// + public virtual double Obj(DVPLI.Vector x) + { double sum = 0; - if(Engine.MultiThread && !displayObjInfo) - { - // Instantiate parallel computation if enabled. - List tl = new List(); - - // Context contains both input parameters and outputs. - List context = new List(); - for (int r = 0; r < this.callMarketPrice.R; r++) + if(Engine.MultiThread && !displayObjInfo) + { + // Instantiate parallel computation if enabled. + List tl = new List(); + + // Context contains both input parameters and outputs. + List context = new List(); + for (int r = 0; r < this.callMarketPrice.R; r++) { if (this.maturity[r] >= this.matBound[0]) { @@ -445,18 +451,18 @@ public virtual double Obj(DVPLI.Vector x) hc.dividend = this.dividendYield[r]; hc.row = r; tl.Add(Task.Factory.StartNew(this.CalculateSingleRow, hc)); - } - } - - tsScheduler.WaitTaskList(tl); - for (int r = 0; r < tl.Count; r++) - sum += context[r].sum; - } - else - { - // Sequential version of the code, used when parallel computation is disabled. - HestonCall hc = new HestonCall(this, x, this.s0); - for (int r = 0; r < this.callMarketPrice.R; r++) + } + } + + tsScheduler.WaitTaskList(tl); + for (int r = 0; r < tl.Count; r++) + sum += context[r].sum; + } + else + { + // Sequential version of the code, used when parallel computation is disabled. + HestonCall hc = new HestonCall(this, x, this.s0); + for (int r = 0; r < this.callMarketPrice.R; r++) { if (this.maturity[r] >= this.matBound[0]) { @@ -470,7 +476,7 @@ public virtual double Obj(DVPLI.Vector x) hc.sum = 0; this.CalculateSingleRow(hc); sum += hc.sum; - } + } } var pricingErrors = hc.hestonCallPrice - this.callMarketPrice; @@ -494,37 +500,37 @@ public virtual double Obj(DVPLI.Vector x) Console.WriteLine(this.callMarketPrice[Range.New(0,RR),Range.New(0,CC)]); Console.WriteLine("Pricing Errors"); Console.WriteLine(pricingErrors[Range.New(0, RR), Range.New(0, CC)]); - } - objCount++; - } - - //Calculate average distance... - sum = Math.Sqrt( sum /this.totalVolume); - if (this.useBoundPenalty) - sum += this.BoundPenalty(x); - - if (this.useFellerPenalty) - sum +=this.FellerPenalty(x); - - return sum; + } + objCount++; + } + + //Calculate average distance... + sum = Math.Sqrt( sum /this.totalVolume); + if (this.useBoundPenalty) + sum += this.BoundPenalty(x); + + if (this.useFellerPenalty) + sum +=this.FellerPenalty(x); + + return sum; } static int objCount = 0; - - //Average pricing error - internal static double avgPricingError = 0; - - /// - /// Calculates a single row of the objective function. Basically - /// options with the same maturity and different strikes. - /// - /// - /// An object of type containing the context. - /// - private void CalculateSingleRow(object context) + + //Average pricing error + internal static double avgPricingError = 0; + + /// + /// Calculates a single row of the objective function. Basically + /// options with the same maturity and different strikes. + /// + /// + /// An object of type containing the context. + /// + private void CalculateSingleRow(object context) { // use the interpolation to calculate the prices CalculateSingleRowWithInterpolation(context); - return; + return; } @@ -623,52 +629,52 @@ private void CalculateSingleRowWithInterpolation(object context) } return; - } - - - static PFunction BuildInterpolationObject(List strike, List price, DVPLUtils.EInterpolationType interpolationType = DVPLUtils.EInterpolationType.SPLINE) + } + + + static PFunction BuildInterpolationObject(List strike, List price, DVPLUtils.EInterpolationType interpolationType = DVPLUtils.EInterpolationType.SPLINE) { PFunction fun = new PFunction((Vector)strike.ToArray(), (Vector)price.ToArray()); fun.m_Function.iType = interpolationType; return fun; - } - - #endregion - - /// - /// Penalty function relative to bounds. - /// - /// - /// Vector of parameters. - /// - /// - /// Penalty value. - /// - private double BoundPenalty(Vector x) - { - Vector t1, t2, t3; - t1 = Bounds.Lb + this.smallValue - x; - t2 = -Bounds.Ub + this.smallValue + x; - - t3 = t1 * t1 * (t1 > 0) + t2 * t2 * (t2 > 0); - return this.k1 * t3.Sum(); - } - - /// - /// Penalty function in order to satisfy Feller condition: - /// 2k theta >= sigma^2 - /// + } + + #endregion + + /// + /// Penalty function relative to bounds. + /// /// - /// Vector of parameters: x=[k, theta, sigma,...] - /// - /// - /// Penalty value. - /// - protected double FellerPenalty(Vector x) - { - double result; - result = Math.Max(0, x[2] * x[2] - 2 * x[0] * x[1]); - return this.s0*this.k2 * result * result; + /// Vector of parameters. + /// + /// + /// Penalty value. + /// + private double BoundPenalty(Vector x) + { + Vector t1, t2, t3; + t1 = Bounds.Lb + this.smallValue - x; + t2 = -Bounds.Ub + this.smallValue + x; + + t3 = t1 * t1 * (t1 > 0) + t2 * t2 * (t2 > 0); + return this.k1 * t3.Sum(); + } + + /// + /// Penalty function in order to satisfy Feller condition: + /// 2k theta >= sigma^2 + /// + /// + /// Vector of parameters: x=[k, theta, sigma,...] + /// + /// + /// Penalty value. + /// + protected double FellerPenalty(Vector x) + { + double result; + result = Math.Max(0, x[2] * x[2] - 2 * x[0] * x[1]); + return this.s0*this.k2 * result * result; } /// @@ -730,5 +736,5 @@ double CalculateWeight(double volume) } - } -} + } +} From e551e259b97b7cb054e7aad6394583dd9dcaad33 Mon Sep 17 00:00:00 2001 From: Giovanni La Cagnina Date: Tue, 7 May 2024 17:25:02 +0200 Subject: [PATCH 2/2] [RD-PRICING] refactoring --- Heston/HestonOptimizationProblem.cs | 758 ++++++++++++++-------------- 1 file changed, 379 insertions(+), 379 deletions(-) diff --git a/Heston/HestonOptimizationProblem.cs b/Heston/HestonOptimizationProblem.cs index 24184c4..2239f9e 100644 --- a/Heston/HestonOptimizationProblem.cs +++ b/Heston/HestonOptimizationProblem.cs @@ -1,91 +1,91 @@ -/* Copyright (C) 2009-2012 Fairmat SRL (info@fairmat.com, http://www.fairmat.com/) - * Author(s): Enrico Degiuli (enrico.degiuli@fairmat.com) - * Matteo Tesser (matteo.tesser@fairmat.com) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU Lesser General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU Lesser General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . +/* Copyright (C) 2009-2012 Fairmat SRL (info@fairmat.com, http://www.fairmat.com/) + * Author(s): Enrico Degiuli (enrico.degiuli@fairmat.com) + * Matteo Tesser (matteo.tesser@fairmat.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . */ -/* - * Notes - * http://www.zeliade.com/whitepapers/zwp-0004.pdf +/* + * Notes + * http://www.zeliade.com/whitepapers/zwp-0004.pdf */ -using System; -using System.Collections.Generic; -using System.Threading.Tasks; +using System; +using System.Collections.Generic; +using System.Threading.Tasks; using DVPLI; -using DVPLDOM; -using DVPLI.TaskScheduler; -using Fairmat.MarketData; +using DVPLDOM; +using DVPLI.TaskScheduler; +using Fairmat.MarketData; using Fairmat.Optimization; - -namespace HestonEstimator -{ - - enum EWeighting + +namespace HestonEstimator +{ + + enum EWeighting { CONSTANT = 0, LINEAR = 1, LOG = 2 - } - - - /// - /// Describes the time-dependent Heston optimization problem. - /// - public class HestonCallOptimizationProblem : IOptimizationProblem + } + + + /// + /// Describes the time-dependent Heston optimization problem. + /// + public class HestonCallOptimizationProblem : IOptimizationProblem { - protected CallPriceMarketData cpmd; - /// - /// The call market price matrix. - /// - internal Matrix callMarketPrice; - - /// - /// The maturity vector relative to the callMarketPrice matrix. - /// - private Vector maturity; - - /// - /// The strike vector relative to callMarketPrice matrix. - /// - protected Vector strike; - - /// - /// The average rate vector up to a given maturity - /// - protected Vector rate; - - /// - /// The avereage dividend yield vector up to a given maturity - /// - protected Vector dividendYield; - - /// - /// Averegate drift up to a given maturity - /// - //private Vector drift; - - /// - /// Establish whether to use the boundary penalty function or not. - /// - public bool useBoundPenalty = false; - - /// - /// Establish whether to use the Feller penalty function or not. - /// Note: it's affect a lot the multi-level single linkage algorithm performance. + protected CallPriceMarketData cpmd; + /// + /// The call market price matrix. + /// + internal Matrix callMarketPrice; + + /// + /// The maturity vector relative to the callMarketPrice matrix. + /// + private Vector maturity; + + /// + /// The strike vector relative to callMarketPrice matrix. + /// + protected Vector strike; + + /// + /// The average rate vector up to a given maturity + /// + protected Vector rate; + + /// + /// The avereage dividend yield vector up to a given maturity + /// + protected Vector dividendYield; + + /// + /// Averegate drift up to a given maturity + /// + //private Vector drift; + + /// + /// Establish whether to use the boundary penalty function or not. + /// + public bool useBoundPenalty = false; + + /// + /// Establish whether to use the Feller penalty function or not. + /// Note: it's affect a lot the multi-level single linkage algorithm performance. /// public bool useFellerPenalty = false; @@ -104,62 +104,62 @@ public class HestonCallOptimizationProblem : IOptimizationProblem /// Builds objective function on relative pricing error. /// static bool optimizeRelativeError = false; - static double pricingMin = 0.0001; + static double pricingMin = 0.0001; static internal bool displayObjInfo = false; static internal double optionThreshold = 0.000; static internal EWeighting weighting = EWeighting.LOG; internal static bool calibrateOnCallOptions = true; internal static bool calibrateOnPutOptions = true; - - /// - /// Small value used in the boundary penalty function. - /// - protected double smallValue = 1e-4; - - /// - /// Value that weights the boundary penalty function. - /// - protected double k1 = 1e10; - - /// - /// Value that weights the Feller inequality penalty function. - /// - protected double k2 = 1e2; - - /// - /// The number of call option on which calibration is performed. - /// + + /// + /// Small value used in the boundary penalty function. + /// + protected double smallValue = 1e-4; + + /// + /// Value that weights the boundary penalty function. + /// + protected double k1 = 1e10; + + /// + /// Value that weights the Feller inequality penalty function. + /// + protected double k2 = 1e2; + + /// + /// The number of call option on which calibration is performed. + /// protected internal int numCall; protected internal int numPut; /// /// Total volume for calls and puts. /// - protected internal double totalVolume; - /// - /// Process starting value. - /// + protected internal double totalVolume; + /// + /// Process starting value. + /// protected double s0; protected Vector matBound; protected Vector strikeBound; - /// - /// Initializes a new instance of the HestonCallOptimizationProblem class using the - /// EquityCalibrationData data structure. - /// - /// - /// An EquityCalibrationData object containing market data for calibration. - /// - /// - /// A vector containing the minimum and maximum values - /// for maturities to be used in calibration. - /// - /// - /// A vector containing the minimum and maximum values - /// for strikes to be used in calibration. - /// - public HestonCallOptimizationProblem(EquityCalibrationData equityCalData, Vector matBound, Vector strikeBound) + /// + /// Initializes a new instance of the HestonCallOptimizationProblem class using the + /// EquityCalibrationData data structure. + /// + /// + /// An EquityCalibrationData object containing market data for calibration. + /// + /// + /// A vector containing the minimum and maximum values + /// for maturities to be used in calibration. + /// + /// + /// A vector containing the minimum and maximum values + /// for strikes to be used in calibration. + /// + public HestonCallOptimizationProblem(EquityCalibrationData equityCalData, Vector matBound, Vector strikeBound) { this.cpmd = equityCalData.Hdata; this.matBound = matBound; @@ -172,72 +172,72 @@ public HestonCallOptimizationProblem(EquityCalibrationData equityCalData, Vector strike: strike, rate: equityCalData.CallMatrixRiskFreeRate, s0: equityCalData.Hdata.S0, - dividendYield: equityCalData.CallMatrixDividendYield); - - - displayObjInfo = false; - } - - /// - /// Initializes a new instance of the HestonCallOptimizationProblem class. - /// - /// A matrix containing call option market prices. - /// - /// Vector of call option maturities relative to callMarketPrice matrix. - /// - /// - /// Vector of call option strikes relative to callMarketPrice matrix. - /// - /// - /// Vector of zero coupon bond rates calculated relative to maturity vector maturities. - /// - /// - /// Vector of dividend yield rates calculated relative to maturity vector maturities. - /// - /// Index/Equity value at the time of calibration. - /// - /// A vector containing the minimum and maximum values - /// for maturities to be used in calibration. - /// - /// - /// A vector containing the minimum and maximum values - /// for strikes to be used in calibration. - /// - [Obsolete] - HestonCallOptimizationProblem(Matrix callMarketPrice, Vector maturity, Vector strike, Vector rate, Vector dividendYield, double s0, Vector matBound, Vector strikeBound) - { + dividendYield: equityCalData.CallMatrixDividendYield); + + + displayObjInfo = false; + } + + /// + /// Initializes a new instance of the HestonCallOptimizationProblem class. + /// + /// A matrix containing call option market prices. + /// + /// Vector of call option maturities relative to callMarketPrice matrix. + /// + /// + /// Vector of call option strikes relative to callMarketPrice matrix. + /// + /// + /// Vector of zero coupon bond rates calculated relative to maturity vector maturities. + /// + /// + /// Vector of dividend yield rates calculated relative to maturity vector maturities. + /// + /// Index/Equity value at the time of calibration. + /// + /// A vector containing the minimum and maximum values + /// for maturities to be used in calibration. + /// + /// + /// A vector containing the minimum and maximum values + /// for strikes to be used in calibration. + /// + [Obsolete] + HestonCallOptimizationProblem(Matrix callMarketPrice, Vector maturity, Vector strike, Vector rate, Vector dividendYield, double s0, Vector matBound, Vector strikeBound) + { this.matBound=matBound; - this.strikeBound = strikeBound; - SetVariables(callMarketPrice, maturity, strike, rate, - dividendYield, s0); - } - - /// - /// Sets several variables used to solve the optimization problem. - /// - /// A matrix containing call option market prices. - /// - /// Vector of call option maturities relative to callMarketPrice matrix. - /// - /// - /// Vector of call option strikes relative to callMarketPrice matrix. - /// - /// - /// Vector of zero coupon bond rates calculated relative to maturity vector maturities. - /// - /// - /// Vector of dividend yield rates calculated relative to maturity vector maturities. - /// - /// Index/Equity value at the time of calibration. - /// - /// A vector containing the minimum and maximum values - /// for maturities to be used in calibration. - /// - /// - /// A vector containing the minimum and maximum values - /// for strikes to be used in calibration. - private void SetVariables(Matrix callMarketPrice, Vector maturity, Vector strike, Vector rate, Vector dividendYield, double s0) - { + this.strikeBound = strikeBound; + SetVariables(callMarketPrice, maturity, strike, rate, + dividendYield, s0); + } + + /// + /// Sets several variables used to solve the optimization problem. + /// + /// A matrix containing call option market prices. + /// + /// Vector of call option maturities relative to callMarketPrice matrix. + /// + /// + /// Vector of call option strikes relative to callMarketPrice matrix. + /// + /// + /// Vector of zero coupon bond rates calculated relative to maturity vector maturities. + /// + /// + /// Vector of dividend yield rates calculated relative to maturity vector maturities. + /// + /// Index/Equity value at the time of calibration. + /// + /// A vector containing the minimum and maximum values + /// for maturities to be used in calibration. + /// + /// + /// A vector containing the minimum and maximum values + /// for strikes to be used in calibration. + private void SetVariables(Matrix callMarketPrice, Vector maturity, Vector strike, Vector rate, Vector dividendYield, double s0) + { this.s0 = s0; var dy = new PFunction(maturity, dividendYield); @@ -253,7 +253,7 @@ private void SetVariables(Matrix callMarketPrice, Vector maturity, Vector strike this.rate = rate; this.maturity = maturity; this.strike = strike; - this.callMarketPrice = callMarketPrice; + this.callMarketPrice = callMarketPrice; this.numCall = 0; callWeight = new Matrix(this.callMarketPrice.R, this.callMarketPrice.C); @@ -305,61 +305,61 @@ private void SetVariables(Matrix callMarketPrice, Vector maturity, Vector strike Console.WriteLine("Strike Bounds:\t" + strikeBound); Console.WriteLine("Maturity Bounds:\t" + matBound); Console.WriteLine("Lb:\t" + Bounds.Lb); - Console.WriteLine("Ub:\t" + Bounds.Ub); + Console.WriteLine("Ub:\t" + Bounds.Ub); if(Engine.Verbose>=2) - PutCallTest(); - } - - - /// - /// Finds the couple of integer {i,j} so that the values - /// vector[i], vector[i+1],..., vector[j-1], vector[j] - /// are all included in the interval (bound[0], bound[1]) - /// while vector[i-1] and vector[j+1] are not. - /// - /// - /// Vector to be filtered. - /// - /// - /// Bounds to be used for the filtering. - /// - /// - /// Integer array representing the index couple. - /// - private int[] FindExtremes(Vector vector, Vector bounds) - { - int[] result = new int[2]; - int i = 0; - - while (vector[i] < bounds[0]) - i++; - result[0] = i; - i = vector.Length - 1; - while (vector[i] > bounds[1]) - i--; - result[1] = i; - return result; + PutCallTest(); } - #region IOptimizationProblem Members - + + /// + /// Finds the couple of integer {i,j} so that the values + /// vector[i], vector[i+1],..., vector[j-1], vector[j] + /// are all included in the interval (bound[0], bound[1]) + /// while vector[i-1] and vector[j+1] are not. + /// + /// + /// Vector to be filtered. + /// + /// + /// Bounds to be used for the filtering. + /// + /// + /// Integer array representing the index couple. + /// + private int[] FindExtremes(Vector vector, Vector bounds) + { + int[] result = new int[2]; + int i = 0; + + while (vector[i] < bounds[0]) + i++; + result[0] = i; + i = vector.Length - 1; + while (vector[i] > bounds[1]) + i--; + result[1] = i; + return result; + } + + #region IOptimizationProblem Members + /// /// If volatility is observed, it may be useful to use this information - /// + /// protected double v0Min=0.001; protected double v0Max = 1; - protected double vLastMin = 0.001; - - /// - /// Gets the bounds for the optimization. - /// - public Bounds Bounds - { - get - { - Bounds bounds = new Bounds(); - - // The order of parameters is k, theta, sigma, rho, V0 (and optionally div.y) + protected double vLastMin = 0.001; + + /// + /// Gets the bounds for the optimization. + /// + public Bounds Bounds + { + get + { + Bounds bounds = new Bounds(); + + // The order of parameters is k, theta, sigma, rho, V0 (and optionally div.y) if (HestonConstantDriftEstimator.impliedDividends) { bounds.Lb = (Vector)new double[] { 0, v0Min, 0.001, -1, vLastMin, 0 }; @@ -369,74 +369,74 @@ public Bounds Bounds { bounds.Lb = (Vector)new double[] { 0, v0Min, 0.001, -1, vLastMin }; bounds.Ub = (Vector)new double[] { 15,v0Max, 2, 1, 1 }; - } - return bounds; - } - } - - /// - /// This method is unused but part of the interface. - /// - /// The parameter is not used. - /// Null as it's unused. - public DVPLI.Vector G(DVPLI.Vector x) - { - throw new NotImplementedException(); - } - - /// - /// This method is unused but part of the interface. - /// - /// The parameter is not used. - /// Nothing the function always throws a NotImplementedException. - public DVPLI.Vector Grad(DVPLI.Vector x) - { - throw new NotImplementedException(); - } - - /// - /// Gets a value indicating whether there are non linear constraints in this - /// optimization problem. In this case there are not. - /// - public bool HasNonLinearConstraints - { - get - { - return false; - } - } - - /// - /// Gets null as we have no linear constrains defined. - /// - public virtual LinearConstraints LinearIneqConstraints - { - get - { - return null; - } - } - - /// - /// Calibration objective function. - /// - /// - /// The vector of parameters. - /// - /// - /// Objective function value. - /// - public virtual double Obj(DVPLI.Vector x) - { + } + return bounds; + } + } + + /// + /// This method is unused but part of the interface. + /// + /// The parameter is not used. + /// Null as it's unused. + public DVPLI.Vector G(DVPLI.Vector x) + { + throw new NotImplementedException(); + } + + /// + /// This method is unused but part of the interface. + /// + /// The parameter is not used. + /// Nothing the function always throws a NotImplementedException. + public DVPLI.Vector Grad(DVPLI.Vector x) + { + throw new NotImplementedException(); + } + + /// + /// Gets a value indicating whether there are non linear constraints in this + /// optimization problem. In this case there are not. + /// + public bool HasNonLinearConstraints + { + get + { + return false; + } + } + + /// + /// Gets null as we have no linear constrains defined. + /// + public virtual LinearConstraints LinearIneqConstraints + { + get + { + return null; + } + } + + /// + /// Calibration objective function. + /// + /// + /// The vector of parameters. + /// + /// + /// Objective function value. + /// + public virtual double Obj(DVPLI.Vector x) + { double sum = 0; - if(Engine.MultiThread && !displayObjInfo) - { - // Instantiate parallel computation if enabled. - List tl = new List(); - - // Context contains both input parameters and outputs. - List context = new List(); - for (int r = 0; r < this.callMarketPrice.R; r++) + if(Engine.MultiThread && !displayObjInfo) + { + // Instantiate parallel computation if enabled. + List tl = new List(); + + // Context contains both input parameters and outputs. + List context = new List(); + for (int r = 0; r < this.callMarketPrice.R; r++) { if (this.maturity[r] >= this.matBound[0]) { @@ -451,18 +451,18 @@ public virtual double Obj(DVPLI.Vector x) hc.dividend = this.dividendYield[r]; hc.row = r; tl.Add(Task.Factory.StartNew(this.CalculateSingleRow, hc)); - } - } - - tsScheduler.WaitTaskList(tl); - for (int r = 0; r < tl.Count; r++) - sum += context[r].sum; - } - else - { - // Sequential version of the code, used when parallel computation is disabled. - HestonCall hc = new HestonCall(this, x, this.s0); - for (int r = 0; r < this.callMarketPrice.R; r++) + } + } + + tsScheduler.WaitTaskList(tl); + for (int r = 0; r < tl.Count; r++) + sum += context[r].sum; + } + else + { + // Sequential version of the code, used when parallel computation is disabled. + HestonCall hc = new HestonCall(this, x, this.s0); + for (int r = 0; r < this.callMarketPrice.R; r++) { if (this.maturity[r] >= this.matBound[0]) { @@ -476,7 +476,7 @@ public virtual double Obj(DVPLI.Vector x) hc.sum = 0; this.CalculateSingleRow(hc); sum += hc.sum; - } + } } var pricingErrors = hc.hestonCallPrice - this.callMarketPrice; @@ -500,37 +500,37 @@ public virtual double Obj(DVPLI.Vector x) Console.WriteLine(this.callMarketPrice[Range.New(0,RR),Range.New(0,CC)]); Console.WriteLine("Pricing Errors"); Console.WriteLine(pricingErrors[Range.New(0, RR), Range.New(0, CC)]); - } - objCount++; - } - - //Calculate average distance... - sum = Math.Sqrt( sum /this.totalVolume); - if (this.useBoundPenalty) - sum += this.BoundPenalty(x); - - if (this.useFellerPenalty) - sum +=this.FellerPenalty(x); - - return sum; + } + objCount++; + } + + //Calculate average distance... + sum = Math.Sqrt( sum /this.totalVolume); + if (this.useBoundPenalty) + sum += this.BoundPenalty(x); + + if (this.useFellerPenalty) + sum +=this.FellerPenalty(x); + + return sum; } static int objCount = 0; - - //Average pricing error - internal static double avgPricingError = 0; - - /// - /// Calculates a single row of the objective function. Basically - /// options with the same maturity and different strikes. - /// - /// - /// An object of type containing the context. - /// - private void CalculateSingleRow(object context) + + //Average pricing error + internal static double avgPricingError = 0; + + /// + /// Calculates a single row of the objective function. Basically + /// options with the same maturity and different strikes. + /// + /// + /// An object of type containing the context. + /// + private void CalculateSingleRow(object context) { // use the interpolation to calculate the prices CalculateSingleRowWithInterpolation(context); - return; + return; } @@ -629,52 +629,52 @@ private void CalculateSingleRowWithInterpolation(object context) } return; - } - - - static PFunction BuildInterpolationObject(List strike, List price, DVPLUtils.EInterpolationType interpolationType = DVPLUtils.EInterpolationType.SPLINE) + } + + + static PFunction BuildInterpolationObject(List strike, List price, DVPLUtils.EInterpolationType interpolationType = DVPLUtils.EInterpolationType.SPLINE) { PFunction fun = new PFunction((Vector)strike.ToArray(), (Vector)price.ToArray()); fun.m_Function.iType = interpolationType; return fun; - } - - #endregion - - /// - /// Penalty function relative to bounds. - /// + } + + #endregion + + /// + /// Penalty function relative to bounds. + /// + /// + /// Vector of parameters. + /// + /// + /// Penalty value. + /// + private double BoundPenalty(Vector x) + { + Vector t1, t2, t3; + t1 = Bounds.Lb + this.smallValue - x; + t2 = -Bounds.Ub + this.smallValue + x; + + t3 = t1 * t1 * (t1 > 0) + t2 * t2 * (t2 > 0); + return this.k1 * t3.Sum(); + } + + /// + /// Penalty function in order to satisfy Feller condition: + /// 2k theta >= sigma^2 + /// /// - /// Vector of parameters. - /// - /// - /// Penalty value. - /// - private double BoundPenalty(Vector x) - { - Vector t1, t2, t3; - t1 = Bounds.Lb + this.smallValue - x; - t2 = -Bounds.Ub + this.smallValue + x; - - t3 = t1 * t1 * (t1 > 0) + t2 * t2 * (t2 > 0); - return this.k1 * t3.Sum(); - } - - /// - /// Penalty function in order to satisfy Feller condition: - /// 2k theta >= sigma^2 - /// - /// - /// Vector of parameters: x=[k, theta, sigma,...] - /// - /// - /// Penalty value. - /// - protected double FellerPenalty(Vector x) - { - double result; - result = Math.Max(0, x[2] * x[2] - 2 * x[0] * x[1]); - return this.s0*this.k2 * result * result; + /// Vector of parameters: x=[k, theta, sigma,...] + /// + /// + /// Penalty value. + /// + protected double FellerPenalty(Vector x) + { + double result; + result = Math.Max(0, x[2] * x[2] - 2 * x[0] * x[1]); + return this.s0*this.k2 * result * result; } /// @@ -736,5 +736,5 @@ double CalculateWeight(double volume) } - } -} + } +}