From 9c201fc0546417d5512de9cf3c87a90ff049f0f4 Mon Sep 17 00:00:00 2001 From: Alessandro Gnoatto Date: Wed, 24 Dec 2025 11:40:13 +0100 Subject: [PATCH 1/7] Introduce tree framework --- .../net/finmath/tree/AbstractTreeProduct.java | 100 +++++++++++++ .../OneDimensionalRiskFactorTreeModel.java | 86 +++++++++++ src/main/java/net/finmath/tree/TreeModel.java | 59 ++++++++ .../java/net/finmath/tree/TreeProduct.java | 23 +++ .../AbstractRecombiningTreeModel.java | 100 +++++++++++++ .../OneDimensionalAssetTreeModel.java | 29 ++++ .../models/BoyleTrinomial.java | 139 ++++++++++++++++++ .../models/CoxRossRubinsteinModel.java | 111 ++++++++++++++ .../models/JarrowRuddModel.java | 102 +++++++++++++ .../AbstractNonPathDependentProduct.java | 66 +++++++++ .../AbstractPathDependentProduct.java | 87 +++++++++++ .../products/AsianOption.java | 95 ++++++++++++ .../products/EuNonPathDependent.java | 50 +++++++ .../products/UsNonPathDependent.java | 65 ++++++++ .../EuropeanAndAmericanOptionPricingTest.java | 102 +++++++++++++ 15 files changed, 1214 insertions(+) create mode 100644 src/main/java/net/finmath/tree/AbstractTreeProduct.java create mode 100644 src/main/java/net/finmath/tree/OneDimensionalRiskFactorTreeModel.java create mode 100644 src/main/java/net/finmath/tree/TreeModel.java create mode 100644 src/main/java/net/finmath/tree/TreeProduct.java create mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/AbstractRecombiningTreeModel.java create mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/OneDimensionalAssetTreeModel.java create mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/models/BoyleTrinomial.java create mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/models/CoxRossRubinsteinModel.java create mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/models/JarrowRuddModel.java create mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractNonPathDependentProduct.java create mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractPathDependentProduct.java create mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/products/AsianOption.java create mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/products/EuNonPathDependent.java create mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/products/UsNonPathDependent.java create mode 100644 src/test/java/net/finmath/tree/EuropeanAndAmericanOptionPricingTest.java diff --git a/src/main/java/net/finmath/tree/AbstractTreeProduct.java b/src/main/java/net/finmath/tree/AbstractTreeProduct.java new file mode 100644 index 0000000000..c6a48317c4 --- /dev/null +++ b/src/main/java/net/finmath/tree/AbstractTreeProduct.java @@ -0,0 +1,100 @@ +package net.finmath.tree; + +import net.finmath.stochastic.RandomVariable; + +/** + * Base class for products that can be priced on a (recombining) with TreeModel. + * This abstraction provides a small template: + * getValue(TreeModel) returns the scalar time–0 price. + * getValue(double, TreeModel) returns the value as a RandomVariable at a given evaluation time. + * getValues(double, TreeModel) must be implemented by subclasses and + * should return the full vector of values (per state) at each time level, + * produced by a backward induction. By convention the element at index + * 0 corresponds to time 0. + * Subclasses (e.g. European, American, Asian ) implement their + * payoff logic and early–exercise features within getValues(double, TreeModel). + */ + +public abstract class AbstractTreeProduct implements TreeProduct { + + /** Contract maturity (in model time units). */ + private final double maturity; + + /** + * Creates a tree-priced product with the given maturity. + * + * @param maturity The contract maturity + * must be compatible with the model time grid. + */ + public AbstractTreeProduct(double maturity){ + this.maturity = maturity; + } + + /** + * Returns the time–0 present value under the given model + * This is a convenience method delegating to getValue(double, TreeModel)}ù + * with evaluationTime = 0.0 and extracting the scalar value from + * the returned RandomVariable. + * + * @param model The tree model to use (providing discounting and conditional expectations). + * @return The time–0 price. + * @throws IllegalArgumentException If model is null or the evaluation fails. + */ + @Override + public double getValue(TreeModel model ){ + RandomVariable v0 = getValue(0.0,model); + return v0.get(0); + } + /** + * Returns the product value at a given evaluation time as a RandomVariable. + * The default implementation calls getValues(double, TreeModel) and + * returns the first component (time–level 0). Subclasses should ensure + * their {@link #getValues(double, TreeModel)} respects this convention. + * + * @param evaluationTime The time at which the value is requested (must be >= 0). + * @param model The tree model to use. + * @return The value at evalutationTime as a random variable on the tree. + * @throws IllegalArgumentException If the inputs are invalid (see validate(double, TreeModel)). + */ + public RandomVariable getValue(double evaluationTime ,TreeModel model) { + RandomVariable[] levels = getValues(evaluationTime,model); + return levels[0]; + } + + /** + * Computes the vector of values at each time/state on the tree (via backward induction) + * Implementations should: + * Validate inputs via validate(double, TreeModel). + * Return an array whose element at index k is a RandomVariable + * representing the value at time level k (with length equal to the number of states at that level). + * + * @param evaluationTime The time at which valuation is performed + * @param model The tree model (spot lattice, discounting, conditional expectation). + * @return An array of value random variables, one per time level, indexed from 0 to maturity level. + */ + public abstract RandomVariable[] getValues(double evaluationTime, TreeModel model); + + /** + * Basic input validation helper used by implementations. + * + * @param t The evaluation time (must be >= 0). + * @param m The model (must not be null). + * @throws IllegalArgumentException If {@code m == null} or {@code t < 0}. + */ + protected void validate(double t, TreeModel m) { + if(m == null) throw new IllegalArgumentException("model null"); + if(t < 0) throw new IllegalArgumentException("evaluationTime < 0"); + } + + + /** + * Returns the contract maturity. + * + * @return The maturity T. + */ + protected double getMaturity(){ + return maturity; + } + +} + diff --git a/src/main/java/net/finmath/tree/OneDimensionalRiskFactorTreeModel.java b/src/main/java/net/finmath/tree/OneDimensionalRiskFactorTreeModel.java new file mode 100644 index 0000000000..70a3ded904 --- /dev/null +++ b/src/main/java/net/finmath/tree/OneDimensionalRiskFactorTreeModel.java @@ -0,0 +1,86 @@ +package net.finmath.tree; + +import java.util.ArrayList; +import java.util.List; +import java.util.function.DoubleUnaryOperator; +import net.finmath.montecarlo.RandomVariableFromDoubleArray; +import net.finmath.stochastic.RandomVariable; + +/** + * This abstract class encapsulates the logic of one-dimensional trees for the simulation of a + * risk factor. At this level, the class could be used for equities or interest rates (e.g. + * the Hull-White short rate model). + * + * @author Carlo Andrea Tramentozzi, Alessandro Gnoatto + */ +public abstract class OneDimensionalRiskFactorTreeModel implements TreeModel { + + /** Cache of level spot S_k */ + private final List riskFactorLevels = new ArrayList<>(); + + /** + * Return number of states at level k (binomial: k+1; trinomial: 2k+1) + * @param k + * @return number of states at level k + */ + public abstract int statesAt(int k); + + /** + * Builds all the realizations X_k[i] + * @param k + * @return + */ + protected abstract RandomVariable buildSpotLevel(int k); + + + /** + * Discounted Conditional expectation: from v_{k+1} (array) to v_k (array) + * @param vNext + * @param k + * @return the conditonal expectation + */ + protected abstract RandomVariable conditionalExpectation(RandomVariable vNext, int k); + + + /** + * Builds (if missing) and return X_k as array, using cache. + * @param k + * @return the risk factor at level k + */ + protected final RandomVariable ensureRiskFactorLevelArray(int k) { + while(riskFactorLevels.size() <= k) { + int next = riskFactorLevels.size(); + riskFactorLevels.add(buildSpotLevel(next)); + } + return riskFactorLevels.get(k); + } + + @Override + public RandomVariable getTransformedValuesAtGivenTimeRV(double time, DoubleUnaryOperator transformFunction) { + int k = (int) Math.round(time / getTimeStep()); + // Checking param to avoid out of bound exception + if (k < 0) k = 0; + if (k > getNumberOfTimes() - 1) k = getNumberOfTimes() - 1; + + RandomVariable sRV = ensureRiskFactorLevelArray(k); + double[] s = sRV.getRealizations(); + double[] v = new double[s.length]; + for (int i = 0; i < s.length; i++) { + v[i] = transformFunction.applyAsDouble(s[i]); + } + return new RandomVariableFromDoubleArray(k * getTimeStep(), v); + } + + @Override + public RandomVariable getConditionalExpectationRV(RandomVariable vNext, int k) { + RandomVariable vK = conditionalExpectation(vNext, k); // hook + return new RandomVariableFromDoubleArray(k * getTimeStep(), vK.getRealizations()); + } + + @Override + public RandomVariable getSpotAtGivenTimeIndexRV(int k) { + RandomVariable sRV = ensureRiskFactorLevelArray(k); + return new RandomVariableFromDoubleArray(k * getTimeStep(), sRV.getRealizations()); + } + +} diff --git a/src/main/java/net/finmath/tree/TreeModel.java b/src/main/java/net/finmath/tree/TreeModel.java new file mode 100644 index 0000000000..a01723a3da --- /dev/null +++ b/src/main/java/net/finmath/tree/TreeModel.java @@ -0,0 +1,59 @@ +package net.finmath.tree; + +import net.finmath.modelling.Model; +import net.finmath.stochastic.RandomVariable; +import java.util.function.DoubleUnaryOperator; + +/** + * General interface rappresenting a tree model with all the common methods + * + * @author Carlo Andrea Tramentozzi + * @version 1.0 + */ + +public interface TreeModel extends Model { + + double getTimeStep(); + int getNumberOfTimes(); + + + /** + * Returns the spot process at a given time transformed pointwise by a function. + * Equivalent to getSpotAtGivenTimeIndexRV(k).apply(transformFunction), but accepts a physical time. + * + * @param timeindex Physical time t; the model maps it to the nearest lattice index k (e.g. round(t/dt)). + * @param transformFunction Function f(s) applied element-wise to S_k. + * @return A RandomVariable at time k with length = number of states at level k + * (binomial: k+1; trinomial: 2k+1), containing f(S_k) pathwise. + */ + RandomVariable getTransformedValuesAtGivenTimeRV(double timeindex, DoubleUnaryOperator transformFunction); + + /** + * One-step risk-neutral discounted conditional expectation operator. + * Takes values defined at level k+1 and projects them back to level k: + * V_k = E^Q[ V_{k+1} | F_k ] * DF (DF is the one-step discount factor). + * + * @param optionValues RandomVariable living at level k+1. + * @param timeindex Level k to which the conditional expectation is taken (0 ≤ k < n). + * @return A RandomVariable at level k with length matching the number of states at k. + */ + RandomVariable getConditionalExpectationRV(RandomVariable optionValues, int timeindex); + + /** + * Returns the spot vector S_k as a RandomVariable at lattice level k. + * + * @param timeindex Level index k (0…n). + * @return A RandomVariable carrying S_k pathwise; length = number of states at level k + * (binomial: k+1; trinomial: 2k+1). The time stamp is typically t_k = k*dt. + */ + RandomVariable getSpotAtGivenTimeIndexRV(int timeindex); + + + /** Getter */ + //double getInitialPrice(); + //double getRiskFreeRate(); + //double getVolatility(); + double getLastTime(); + + +} diff --git a/src/main/java/net/finmath/tree/TreeProduct.java b/src/main/java/net/finmath/tree/TreeProduct.java new file mode 100644 index 0000000000..c0c5d7f1ab --- /dev/null +++ b/src/main/java/net/finmath/tree/TreeProduct.java @@ -0,0 +1,23 @@ +package net.finmath.tree; + +/** + * Contract for a product that can be priced on a (recombining) tree model. + * Implementations encapsulate the payoff and exercise style + * (e.g. European, American, path-dependent) and use the TreeModel's + * risk-neutral dynamics and discounting to compute a present value. + * */ + +public interface TreeProduct { + + /** + * Prices this product under the given tree model. + * + * @param model + * The tree model providing spot levels, discounting and + * conditional expectations (e.g. CRR, Jarrow–Rudd, Leisen–Reimer, + * or a trinomial model); must not be {@code null}. + * @return The present value (time 0) of the product under {@code model}, + * in the same currency units as the model’s numéraire. + */ + double getValue(TreeModel model ); +} diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/AbstractRecombiningTreeModel.java b/src/main/java/net/finmath/tree/assetderivativevaluation/AbstractRecombiningTreeModel.java new file mode 100644 index 0000000000..a1bcde7015 --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/AbstractRecombiningTreeModel.java @@ -0,0 +1,100 @@ +package net.finmath.tree.assetderivativevaluation; + +import net.finmath.montecarlo.RandomVariableFromDoubleArray; +import net.finmath.stochastic.RandomVariable; +import net.finmath.tree.OneDimensionalRiskFactorTreeModel; +import java.util.ArrayList; +import java.util.List; + +/** + * + * This class represent the recombining tree model, a general abstract class which collect the common properties + * about recombining trees. + * + * @author Carlo Andrea Tramentozzi + * @version 1.0 + */ + +public abstract class AbstractRecombiningTreeModel extends OneDimensionalRiskFactorTreeModel implements OneDimensionalAssetTreeModel { + + /** Common parameters of the model */ + private final double initialPrice; + private final double riskFreeRate; + private final double volatility; + + /** Temporal discretization */ + private final double timeStep; + private final double lastTime; + private final int numberOfTimes; + /** Cache of level spot S_k */ + private final List spotLevels = new ArrayList<>(); + + /** + * Constructs a recombining tree model on an equidistant time grid by + * specifying the time step. + * The constructor sets the common model parameters, computes + * the number of time points and initializes the level-0 cache of + * spot values with S_0. + * @param initialPrice The initial asset price S0 + * @param riskFreeRate The continuously compounded risk-free rate r used consistently with df(). + * @param volatility The volatility. + * @param lastTime The final time T of the grid. + * @param timeStep The time step Δt between consecutive grid points + */ + public AbstractRecombiningTreeModel(double initialPrice, double riskFreeRate, double volatility, double lastTime, double timeStep) { + this.initialPrice = initialPrice; + this.riskFreeRate = riskFreeRate; + this.volatility = volatility; + this.lastTime = lastTime; + this.timeStep = timeStep; + int steps = (int)Math.round(lastTime / timeStep); + this.numberOfTimes = steps + 1; + spotLevels.add(new RandomVariableFromDoubleArray(0.0, new double[] { initialPrice })); + } + + /** + * Constructs a recombining tree model on an equidistant time grid by + * specifying the time step. + * The constructor sets the common model parameters , computes + * the number of time points and initializes the level-0 cache of + * spot values with S_0. + * @param initialPrice The initial asset price S₀. + * @param riskFreeRate The continuously compounded risk-free rate r used consistently with df(). + * @param volatility The (log) volatility σ. + * @param lastTime The final time T of the grid. + * @param numberOfTimes The total number of grid points N (must be ≥ 2). The number of steps is N−1. + */ + public AbstractRecombiningTreeModel(double initialPrice, double riskFreeRate, double volatility, double lastTime, int numberOfTimes) { + this.initialPrice = initialPrice; + this.riskFreeRate = riskFreeRate; + this.volatility = volatility; + this.lastTime = lastTime; + this.numberOfTimes = numberOfTimes; + this.timeStep = lastTime / (numberOfTimes - 1); + + spotLevels.add(new RandomVariableFromDoubleArray(0.0, new double[] { initialPrice })); + } + + + /** One step discount factor(according to q): continuous */ + protected final double df() { + return Math.exp(-riskFreeRate * timeStep); + } + + + + /** Getter */ + @Override + public double getInitialPrice() { return initialPrice; } + @Override + public double getRiskFreeRate() { return riskFreeRate; } + @Override + public double getVolatility() { return volatility; } + + public double getTimeStep() { return timeStep; } + public double getLastTime() { return lastTime; } + public int getNumberOfTimes() { return numberOfTimes; } + public double getSpot() {return getInitialPrice(); } //Redefinition + + +} diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/OneDimensionalAssetTreeModel.java b/src/main/java/net/finmath/tree/assetderivativevaluation/OneDimensionalAssetTreeModel.java new file mode 100644 index 0000000000..0c1c67bb65 --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/OneDimensionalAssetTreeModel.java @@ -0,0 +1,29 @@ +package net.finmath.tree.assetderivativevaluation; + +/** + * From this level of abstraction, we mean that the tree should represent a stock price + * and not, for example, an interest rate. + * + * @author Alessandro Gnoatto + */ +public interface OneDimensionalAssetTreeModel { + + /** + * The risk factor is a stock + * @return the initial stock price + */ + public double getInitialPrice(); + + /** + * The risk free rate + * @return the risk free rate + */ + public double getRiskFreeRate(); + + /** + * Returns the volatility + * @return the volatility + */ + public double getVolatility(); + +} diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/models/BoyleTrinomial.java b/src/main/java/net/finmath/tree/assetderivativevaluation/models/BoyleTrinomial.java new file mode 100644 index 0000000000..a6dbe5e873 --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/models/BoyleTrinomial.java @@ -0,0 +1,139 @@ +package net.finmath.tree.assetderivativevaluation.models; + +import net.finmath.montecarlo.RandomVariableFromDoubleArray; +import net.finmath.stochastic.RandomVariable; +import net.finmath.tree.assetderivativevaluation.AbstractRecombiningTreeModel; + +/** + * This class is used in order to construct a trinomial model. The trinomial model is a discrete model + * for a stochastic process S, such that every time n we have + * S(n+1)=S(n)*M(n), + * where M(n) can take values u, 1 or 1/u with probabilities q_u, q_m and q_d, respectively. + * This is done under a martingale measure: it can be seen (see notes) that there exists infinitely many + * martingale measures once one fixed q_u. + * For any time index i, all possible value of the process are computed starting from the biggest one. + * + * @author Andrea Mazzon + * @version 1.0 + */ + +public class BoyleTrinomial extends AbstractRecombiningTreeModel { + + /** Specific parameters */ + private final double dt; + private final double r; + private final double u; + private final double d; + private final double pu, pm, pd; + + /** + * It constructs an object which represents the approximation of a Black-Scholes model via the Jarrow-Rudd model. + * + * @param spotPrice, the spot price of the asset modeled by the process + * @param riskFreeRate, the number r such that the value of a risk-free bond at time T is e^(rT) + * @param volatility, the log-volatility of the Black-Scholes model + * @param lastTime, the last time T in the time discretization 0=t_0= euTRI, "US(Call) TRI must be >= EU(Call) TRI"); + + assertEquals(euCRR, euTRI, 2*tol, "EU(Call) CRR vs TRI difference beyond tolerance"); + + + assertEquals(euJR, euCRR, 2*tol, "EU(Call) JR vs CRR differenza oltre tolleranza"); + assertEquals(euJR, euTRI, 2*tol, "EU(Call) JR vs TRI differenza oltre tolleranza"); + + assertEquals(usJR, usCRR, 2*tol, "US(Call) JR vs CRR differenza oltre tolleranza"); + assertEquals(usJR, usTRI, 2*tol, "US(Call) JR vs TRI differenza oltre tolleranza"); + + assertTrue(usCRR + 1e-12 >= euCRR, "US(Call) CRR must be >= EU(Call) CRR"); + } + + @Test + void testNonPathDep_Put() { + double spot = 100.0; + double rate = 0.03; + double vol = 0.25; + double maturity = 2.0; + int steps = 300; + double strike = 110.0; + double tol = 2e-2; + + CoxRossRubinsteinModel crr = new CoxRossRubinsteinModel(spot, rate, vol, maturity, steps); + JarrowRuddModel jr = new JarrowRuddModel(spot,rate,vol,maturity,steps); + BoyleTrinomial tri = new BoyleTrinomial(spot, rate, vol, maturity, steps); + + + + EuNonPathDependent euPut = new EuNonPathDependent(maturity, s -> Math.max(strike - s, 0.0)); + UsNonPathDependent usPut = new UsNonPathDependent(maturity, s -> Math.max(strike - s, 0.0)); + + + double euCRR = euPut.getValue(0.0, crr).getAverage(); + double euJR = euPut.getValue(0.0, jr).getAverage(); + double euTRI = euPut.getValue(0.0, tri).getAverage(); + + + double usCRR = usPut.getValue(0.0, crr).getAverage(); + double usJR = usPut.getValue(0.0,jr).getAverage(); + double usTRI = usPut.getValue(0.0, tri).getAverage(); + + assertTrue(usCRR + 1e-12 >= euCRR, "US(Put) CRR must be >= EU(Put) CRR"); + + assertEquals(euJR, euCRR, 2*tol, "EU(Put) JR vs CRR differenza oltre tolleranza"); + assertEquals(euJR, euTRI, 2*tol, "EU(Put) JR vs TRI differenza oltre tolleranza"); + + assertEquals(usJR, usCRR, 2*tol, "US(Put) JR vs CRR differenza oltre tolleranza"); + assertEquals(usJR, usTRI, 2*tol, "US(Put) JR vs TRI differenza oltre tolleranza"); + + assertTrue(usTRI + 1e-12 >= euTRI, "US(Put) TRI must be >= EU(Put) TRI"); + assertEquals(euCRR, euTRI, 2*tol, "EU(Put) CRR vs TRI difference beyond tolerance"); + } +} \ No newline at end of file From 07d22f18a453b9c5bf40c0d3279f10ee082238b2 Mon Sep 17 00:00:00 2001 From: Alessandro Gnoatto Date: Wed, 24 Dec 2025 11:46:15 +0100 Subject: [PATCH 2/7] Delete Asian Option --- .../products/AsianOption.java | 95 ------------------- 1 file changed, 95 deletions(-) delete mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/products/AsianOption.java diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/products/AsianOption.java b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AsianOption.java deleted file mode 100644 index e02102d1d2..0000000000 --- a/src/main/java/net/finmath/tree/assetderivativevaluation/products/AsianOption.java +++ /dev/null @@ -1,95 +0,0 @@ -package net.finmath.tree.assetderivativevaluation.products; - -import net.finmath.stochastic.RandomVariable; -import net.finmath.time.TimeDiscretization; -import net.finmath.tree.TreeModel; - -import java.util.HashSet; -import java.util.Set; -import java.util.function.DoubleUnaryOperator; - -/** - * Fixed-strike Asian option priced on a recombining tree. - * The averaging times are a subset of the tree time grid. - * - * The running average includes S0 if 0.0 is contained in the averaging times. - */ -public class AsianOption extends AbstractPathDependentProduct { - - private final TimeDiscretization averagingTimes; - private final Set averagingIndices; - private final int numberOfAveragingTimes; - - /** - * @param maturity Option maturity T - * @param averagingTimes Discrete averaging times (subset of model grid) - * @param payOffFunction Payoff applied to the arithmetic average - */ - public AsianOption( - double maturity, - TimeDiscretization averagingTimes, - DoubleUnaryOperator payOffFunction) { - super(maturity, payOffFunction); - this.averagingTimes = averagingTimes; - this.averagingIndices = new HashSet<>(); - this.numberOfAveragingTimes = averagingTimes.getNumberOfTimes(); - } - - @Override - public RandomVariable[] getValues(double evaluationTime, TreeModel model) { - - validate(evaluationTime, model); - - int maturityIndex = timeToIndex(getMaturity(), model); - int evaluationIndex = timeToIndex(evaluationTime, model); - - // Map averaging times → tree indices - averagingIndices.clear(); - for (int i = 0; i < averagingTimes.getNumberOfTimes(); i++) { - int idx = timeToIndex(averagingTimes.getTime(i), model); - averagingIndices.add(idx); - } - - // ---------- Forward construction of running sum ---------- - RandomVariable[] runningSum = new RandomVariable[maturityIndex + 1]; - - for (int k = 0; k <= maturityIndex; k++) { - RandomVariable Sk = model.getSpotAtGivenTimeIndexRV(k); - - if (k == 0) { - runningSum[k] = averagingIndices.contains(0) - ? Sk - : Sk.mult(0.0); - } else { - runningSum[k] = runningSum[k - 1]; - if (averagingIndices.contains(k)) { - runningSum[k] = runningSum[k].add(Sk); - } - } - } - - // ---------- Payoff at maturity ---------- - RandomVariable average = - runningSum[maturityIndex].div(numberOfAveragingTimes); - - RandomVariable value = - model.getTransformedValuesAtGivenTimeRV( - getMaturity(), - x -> getPayOffFunction().applyAsDouble(x) - ); - - // Replace payoff argument with average - value = average.apply(getPayOffFunction()); - - // ---------- Backward induction ---------- - RandomVariable[] values = new RandomVariable[maturityIndex - evaluationIndex + 1]; - values[values.length - 1] = value; - - for (int k = maturityIndex - 1; k >= evaluationIndex; k--) { - value = model.getConditionalExpectationRV(value, k); - values[k - evaluationIndex] = value; - } - - return values; - } -} From d279f952886aad240d9b37282b7209baf97040e1 Mon Sep 17 00:00:00 2001 From: Alessandro Gnoatto Date: Wed, 24 Dec 2025 12:41:28 +0100 Subject: [PATCH 3/7] Fix Junit test --- .../EuropeanAndAmericanOptionPricingTest.java | 63 ++++++++++--------- 1 file changed, 33 insertions(+), 30 deletions(-) diff --git a/src/test/java/net/finmath/tree/EuropeanAndAmericanOptionPricingTest.java b/src/test/java/net/finmath/tree/EuropeanAndAmericanOptionPricingTest.java index 28e5d9214b..fa45459538 100644 --- a/src/test/java/net/finmath/tree/EuropeanAndAmericanOptionPricingTest.java +++ b/src/test/java/net/finmath/tree/EuropeanAndAmericanOptionPricingTest.java @@ -1,7 +1,8 @@ package net.finmath.tree; -import org.junit.jupiter.api.Test; +import org.junit.Assert; +import org.junit.Test; import net.finmath.tree.assetderivativevaluation.models.BoyleTrinomial; import net.finmath.tree.assetderivativevaluation.models.CoxRossRubinsteinModel; @@ -9,13 +10,10 @@ import net.finmath.tree.assetderivativevaluation.products.EuNonPathDependent; import net.finmath.tree.assetderivativevaluation.products.UsNonPathDependent; -import static org.junit.jupiter.api.Assertions.assertEquals; -import static org.junit.jupiter.api.Assertions.assertTrue; - public class EuropeanAndAmericanOptionPricingTest { @Test - void testNonPathDep_Call() { + public void testNonPathDep_Call() { double spot = 100.0; double rate = 0.05; @@ -44,23 +42,25 @@ void testNonPathDep_Call() { double usTRI = usCall.getValue(0.0,tri).getAverage(); System.out.println(euCRR + " " + euJR + " " + euTRI); - - assertTrue(usTRI + 1e-12 >= euTRI, "US(Call) TRI must be >= EU(Call) TRI"); - - assertEquals(euCRR, euTRI, 2*tol, "EU(Call) CRR vs TRI difference beyond tolerance"); - - - assertEquals(euJR, euCRR, 2*tol, "EU(Call) JR vs CRR differenza oltre tolleranza"); - assertEquals(euJR, euTRI, 2*tol, "EU(Call) JR vs TRI differenza oltre tolleranza"); - - assertEquals(usJR, usCRR, 2*tol, "US(Call) JR vs CRR differenza oltre tolleranza"); - assertEquals(usJR, usTRI, 2*tol, "US(Call) JR vs TRI differenza oltre tolleranza"); - - assertTrue(usCRR + 1e-12 >= euCRR, "US(Call) CRR must be >= EU(Call) CRR"); + //"US(Call) TRI must be >= EU(Call) TRI" + Assert.assertTrue(usTRI + 1e-12 >= euTRI); + //"EU(Call) CRR vs TRI difference beyond tolerance" + Assert.assertEquals(euCRR, euTRI, 2*tol); + + //"EU(Call) JR vs CRR difference beyond tol" + Assert.assertEquals(euJR, euCRR, 2*tol); + //"EU(Call) JR vs TRI difference beyond tol" + Assert.assertEquals(euJR, euTRI, 2*tol); + //"US(Call) JR vs CRR difference beyond tol" + Assert.assertEquals(usJR, usCRR, 2*tol); + //"US(Call) JR vs TRI difference beyond tol" + Assert.assertEquals(usJR, usTRI, 2*tol); + //"US(Call) CRR must be >= EU(Call) CRR" + Assert.assertTrue(usCRR + 1e-12 >= euCRR); } @Test - void testNonPathDep_Put() { + public void testNonPathDep_Put() { double spot = 100.0; double rate = 0.03; double vol = 0.25; @@ -87,16 +87,19 @@ void testNonPathDep_Put() { double usCRR = usPut.getValue(0.0, crr).getAverage(); double usJR = usPut.getValue(0.0,jr).getAverage(); double usTRI = usPut.getValue(0.0, tri).getAverage(); - - assertTrue(usCRR + 1e-12 >= euCRR, "US(Put) CRR must be >= EU(Put) CRR"); - - assertEquals(euJR, euCRR, 2*tol, "EU(Put) JR vs CRR differenza oltre tolleranza"); - assertEquals(euJR, euTRI, 2*tol, "EU(Put) JR vs TRI differenza oltre tolleranza"); - - assertEquals(usJR, usCRR, 2*tol, "US(Put) JR vs CRR differenza oltre tolleranza"); - assertEquals(usJR, usTRI, 2*tol, "US(Put) JR vs TRI differenza oltre tolleranza"); - - assertTrue(usTRI + 1e-12 >= euTRI, "US(Put) TRI must be >= EU(Put) TRI"); - assertEquals(euCRR, euTRI, 2*tol, "EU(Put) CRR vs TRI difference beyond tolerance"); + //"US(Put) CRR must be >= EU(Put) CRR" + Assert.assertTrue(usCRR + 1e-12 >= euCRR); + //"EU(Put) JR vs CRR difference beyond tol" + Assert.assertEquals(euJR, euCRR, 2*tol); + //"EU(Put) JR vs TRI difference beyond tol" + Assert.assertEquals(euJR, euTRI, 2*tol); + //"US(Put) JR vs CRR difference beyond tol" + Assert.assertEquals(usJR, usCRR, 2*tol); + //"US(Put) JR vs TRI difference beyond tol" + Assert.assertEquals(usJR, usTRI, 2*tol); + //"US(Put) TRI must be >= EU(Put) TRI" + Assert.assertTrue(usTRI + 1e-12 >= euTRI); + //"EU(Put) CRR vs TRI difference beyond tolerance" + Assert.assertEquals(euCRR, euTRI, 2*tol); } } \ No newline at end of file From a3b8c1c8c5aca30ad8720edb0bea5119b9a575cb Mon Sep 17 00:00:00 2001 From: Alessandro Gnoatto Date: Thu, 5 Feb 2026 12:19:34 +0100 Subject: [PATCH 4/7] Fix path dependent products. Add Asian Options. Add test. --- .../net/finmath/tree/AbstractTreeProduct.java | 17 +- src/main/java/net/finmath/tree/TreeModel.java | 70 ++- .../AbstractRecombiningTreeModel.java | 3 +- .../models/BoyleTrinomial.java | 23 + .../models/CoxRossRubinsteinModel.java | 23 +- .../models/JarrowRuddModel.java | 23 + .../AbstractNonPathDependentProduct.java | 3 +- .../AbstractPathDependentProduct.java | 400 ++++++++++++++---- .../products/AsianOption.java | 74 ++++ .../products/FloatingStrikeAsianOption.java | 76 ++++ .../tree/PathDependentOptionTreeTest.java | 130 ++++++ 11 files changed, 752 insertions(+), 90 deletions(-) create mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/products/AsianOption.java create mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/products/FloatingStrikeAsianOption.java create mode 100644 src/test/java/net/finmath/tree/PathDependentOptionTreeTest.java diff --git a/src/main/java/net/finmath/tree/AbstractTreeProduct.java b/src/main/java/net/finmath/tree/AbstractTreeProduct.java index c6a48317c4..f459bfbbda 100644 --- a/src/main/java/net/finmath/tree/AbstractTreeProduct.java +++ b/src/main/java/net/finmath/tree/AbstractTreeProduct.java @@ -13,8 +13,11 @@ * 0 corresponds to time 0. * Subclasses (e.g. European, American, Asian ) implement their * payoff logic and early–exercise features within getValues(double, TreeModel). + * + * @author Carlo Andrea Tramentozzi + * @author Alessandro Gnoatto + * @author Andrea Mazzon */ - public abstract class AbstractTreeProduct implements TreeProduct { /** Contract maturity (in model time units). */ @@ -50,15 +53,21 @@ public double getValue(TreeModel model ){ * The default implementation calls getValues(double, TreeModel) and * returns the first component (time–level 0). Subclasses should ensure * their {@link #getValues(double, TreeModel)} respects this convention. + * Internally, the method creates all values from time zero and returns the level + * of interest. While this is inefficient for non-path-dependent products, it ensures + * a correct implementation for path dependent products. The alternative would be to provide as + * input the current value of the path-dependent functional, which would break the interface. * * @param evaluationTime The time at which the value is requested (must be >= 0). * @param model The tree model to use. * @return The value at evalutationTime as a random variable on the tree. * @throws IllegalArgumentException If the inputs are invalid (see validate(double, TreeModel)). */ - public RandomVariable getValue(double evaluationTime ,TreeModel model) { - RandomVariable[] levels = getValues(evaluationTime,model); - return levels[0]; + public RandomVariable getValue(double evaluationTime, TreeModel model) { + RandomVariable[] levels = getValues(0.0,model); + //This is dangerous as it assumes a uniform time discretization + int index = (int) Math.round(evaluationTime / model.getTimeStep()); + return levels[index]; } /** diff --git a/src/main/java/net/finmath/tree/TreeModel.java b/src/main/java/net/finmath/tree/TreeModel.java index a01723a3da..4ef01499a6 100644 --- a/src/main/java/net/finmath/tree/TreeModel.java +++ b/src/main/java/net/finmath/tree/TreeModel.java @@ -8,14 +8,15 @@ * General interface rappresenting a tree model with all the common methods * * @author Carlo Andrea Tramentozzi + * @author Alessandro Gnoatto + * @author Andrea Mazzon * @version 1.0 */ - public interface TreeModel extends Model { double getTimeStep(); int getNumberOfTimes(); - + double getLastTime(); /** * Returns the spot process at a given time transformed pointwise by a function. @@ -48,12 +49,67 @@ public interface TreeModel extends Model { */ RandomVariable getSpotAtGivenTimeIndexRV(int timeindex); + /** + * Returns the number of outgoing branches from a node at time index k. + *

+ * This method is intentionally defined per node to allow state-dependent branching + * (e.g. for trinomial with pruning, multinomial with barriers, etc.). For standard + * binomial / trinomial recombining trees this is constant (2 / 3) and independent + * of {@code stateIndex}. + * + * @param timeIndex The time index k (0 ≤ k < n). + * @param stateIndex The node / state index at time k (model-defined). + * @return The number of branches from that node to time k+1. + */ + default int getNumberOfBranches(int timeIndex, int stateIndex) { + throw new UnsupportedOperationException("Transition API not implemented by this TreeModel."); + } - /** Getter */ - //double getInitialPrice(); - //double getRiskFreeRate(); - //double getVolatility(); - double getLastTime(); + /** + * Returns the risk-neutral transition probability for a branch. + * Branch indices are model-defined but should be consistent with how the model + * evolves the risk factor (e.g. for binomial: 0=up,1=down; for trinomial: 0=up,1=mid,2=down). + * + * @param timeIndex The time index k (0 ≤ k < n). + * @param stateIndex The node / state index at time k (model-defined). + * @param branchIndex The branch index in {@code [0, getNumberOfBranches(k,i))}. + * @return Transition probability to the corresponding child on level k+1. + */ + default double getTransitionProbability(int timeIndex, int stateIndex, int branchIndex) { + throw new UnsupportedOperationException("Transition API not implemented by this TreeModel."); + } + + /** + * One-step discount factor from time k to k+1 (applied when taking the conditional expectation). + * For example, for a constant short rate r and step dt, DF = exp(-r*dt). + * + * @param timeIndex The time index k (0 ≤ k < n). + * @return One-step discount factor DF(k,k+1). + */ + default double getOneStepDiscountFactor(int timeIndex) { + throw new UnsupportedOperationException("Transition API not implemented by this TreeModel."); + } + /** + * Returns the model's child-index shift convention for recombining state indices. + * + * Convention: childIndex = parentIndex + shift[branchIndex]. + * + * Examples: + * - Binomial CRR/JR: {0, 1} interpreted as {up, down} + * - Trinomial Boyle: {0, 1, 2} interpreted as {up, mid, down} + * + * Path-dependent products that build a full non-recombining tree (exponential growth) + * can use this to map their parent recombining state index to the child's recombining + * state index when querying model spots. + * + * Default throws: models supporting path-dependent products should override this. + * + * @return shift array of length = number of branches. + */ + default int[] getChildStateIndexShift() { + throw new UnsupportedOperationException("Child state index shift not implemented by this TreeModel."); + } + } diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/AbstractRecombiningTreeModel.java b/src/main/java/net/finmath/tree/assetderivativevaluation/AbstractRecombiningTreeModel.java index a1bcde7015..f00539ee3c 100644 --- a/src/main/java/net/finmath/tree/assetderivativevaluation/AbstractRecombiningTreeModel.java +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/AbstractRecombiningTreeModel.java @@ -77,7 +77,8 @@ public AbstractRecombiningTreeModel(double initialPrice, double riskFreeRate, do /** One step discount factor(according to q): continuous */ - protected final double df() { + @Override + public double getOneStepDiscountFactor(int timeIndex) { return Math.exp(-riskFreeRate * timeStep); } diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/models/BoyleTrinomial.java b/src/main/java/net/finmath/tree/assetderivativevaluation/models/BoyleTrinomial.java index a6dbe5e873..beee80ce76 100644 --- a/src/main/java/net/finmath/tree/assetderivativevaluation/models/BoyleTrinomial.java +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/models/BoyleTrinomial.java @@ -136,4 +136,27 @@ protected RandomVariable conditionalExpectation(RandomVariable vNext, int k) { double time = k * dt; return new RandomVariableFromDoubleArray(time, now); } + + @Override + public int getNumberOfBranches(int timeIndex, int stateIndex) { + return 3; + } + + @Override + public double getTransitionProbability(int timeIndex, int stateIndex, int branchIndex) { + // Convention: 0 = up, 1 = middle, 2 = down + switch(branchIndex) { + case 0: return pu; + case 1: return pm; + case 2: return pd; + default: throw new IllegalArgumentException("Invalid branchIndex " + branchIndex + " for trinomial model."); + } + } + + @Override + public int[] getChildStateIndexShift() { + // Convention: childIndex = parentIndex + shift[branchIndex] + // children at next level correspond to indices stateIndex + {0,1,2} for {up, mid, down}. + return new int[] { 0, 1, 2}; + } } diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/models/CoxRossRubinsteinModel.java b/src/main/java/net/finmath/tree/assetderivativevaluation/models/CoxRossRubinsteinModel.java index f058afad26..242cbed19c 100644 --- a/src/main/java/net/finmath/tree/assetderivativevaluation/models/CoxRossRubinsteinModel.java +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/models/CoxRossRubinsteinModel.java @@ -94,7 +94,7 @@ protected RandomVariable buildSpotLevel(int k) { protected RandomVariable conditionalExpectation(RandomVariable vNext, int k) { double[] next = vNext.getRealizations(); double[] vK = new double[statesAt(k)]; - double disc = df(); + double disc = getOneStepDiscountFactor(k); for (int i = 0; i <= k; i++) { vK[i] = disc * (q * next[i] + (1.0 - q) * next[i + 1]); @@ -102,6 +102,27 @@ protected RandomVariable conditionalExpectation(RandomVariable vNext, int k) { return new RandomVariableFromDoubleArray(k * getTimeStep(), vK); } + @Override + public int getNumberOfBranches(int timeIndex, int stateIndex) { + return 2; + } + + @Override + public double getTransitionProbability(int timeIndex, int stateIndex, int branchIndex) { + // Convention: 0 = up, 1 = down + switch(branchIndex) { + case 0: return q; + case 1: return 1.0 - q; + default: throw new IllegalArgumentException("Invalid branchIndex " + branchIndex + " for binomial model."); + } + } + + @Override + public int[] getChildStateIndexShift() { + // Convention: childIndex = parentIndex + shift[branchIndex] + // For binomial recombining trees: up keeps index, down increases index by 1. + return new int[] { 0, 1 }; + } /** Getter */ public double getU() { return u; } diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/models/JarrowRuddModel.java b/src/main/java/net/finmath/tree/assetderivativevaluation/models/JarrowRuddModel.java index a56e0a5ccc..ec34f25ff3 100644 --- a/src/main/java/net/finmath/tree/assetderivativevaluation/models/JarrowRuddModel.java +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/models/JarrowRuddModel.java @@ -99,4 +99,27 @@ protected RandomVariable conditionalExpectation(RandomVariable vNext, int k) { return new RandomVariableFromDoubleArray(k * dt, res); } + + @Override + public int getNumberOfBranches(int timeIndex, int stateIndex) { + return 2; + } + + @Override + public double getTransitionProbability(int timeIndex, int stateIndex, int branchIndex) { + // Convention: 0 = up, 1 = down (Jarrow-Rudd uses p=0.5) + switch(branchIndex) { + case 0: return 0.5; + case 1: return 0.5; + default: throw new IllegalArgumentException("Invalid branchIndex " + branchIndex + " for binomial model."); + } + } + + @Override + public int[] getChildStateIndexShift() { + // Convention: childIndex = parentIndex + shift[branchIndex] + // For binomial recombining trees: up keeps index, down increases index by 1. + return new int[] { 0, 1 }; + } + } diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractNonPathDependentProduct.java b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractNonPathDependentProduct.java index 598263998b..1a8ecd1d27 100644 --- a/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractNonPathDependentProduct.java +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractNonPathDependentProduct.java @@ -13,8 +13,9 @@ * at maturity (e.g. s -> Math.max(s-K,0) for a call). Early–exercise features * (if any) and the backward induction logic are delegated to subclasses via * getValues(double, TreeModel). + * + * @author Carlo Andrea Tramentozzi */ - public abstract class AbstractNonPathDependentProduct extends AbstractTreeProduct { /** Payoff function f(S) applied at maturity (or when exercised). */ diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractPathDependentProduct.java b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractPathDependentProduct.java index a492755287..5c431302af 100644 --- a/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractPathDependentProduct.java +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractPathDependentProduct.java @@ -1,87 +1,335 @@ package net.finmath.tree.assetderivativevaluation.products; +import java.util.Arrays; + +import net.finmath.stochastic.RandomVariable; +import net.finmath.montecarlo.RandomVariableFromDoubleArray; import net.finmath.tree.AbstractTreeProduct; import net.finmath.tree.TreeModel; -import net.finmath.stochastic.RandomVariable; - -import java.util.function.DoubleUnaryOperator; /** - * Base class for path–dependent tree products (e.g. Asian options). - * It defines the common contract and utilities to: - * Store a payoff function f(A) or f(S) depending on the product - * Map continuous times to the model’s discrete time index - * Expose convenience getters for values at a given evaluation time. + * Base class for full (non-recombining) path-dependent products on a tree. + * + * Key property: + * - The number of nodes grows exponentially with the number of time steps: + * level j has branchingFactor^j nodes. + * - NO recombination / compression / bucketing. + * + * This class builds a "full path tree" of: + * - an underlying recombining state index (for reading spot from the model's spot level), + * - plus product-specific path state variables (e.g. running sum/count for Asian). We build the + * full (non-recombining) path tree (size branchingFactor^j after j steps) + * + * After the forward loop, it then prices by backward induction using TreeModel's transition API: + * - p = model.getTransitionProbability(timeIndex, recombStateIndex, branchIndex) + * - df = model.getOneStepDiscountFactor(timeIndex) + * + * Design choices: + * - We keep state variables in primitive arrays for speed and to avoid per-node objects. + * - We treat the model's spot as a function of (timeIndex, recombStateIndex) so we don't + * need the model to expose explicit "spot multipliers". + * + * IMPORTANT: + * - To evolve the recombining state index along a branch we require a "child state shift" + * convention. For the existing models (CRR/JR/Boyle) this matches their internal layout: + * Binomial: childShift = {0, 1} (up keeps index, down increases index) + * Trinomial (Boyle): childShift = {0, 1, 2} (up, mid, down) + * For a general multinomial model, one has to provide childShift explicitly. + * + * @author Alessandro Gnoatto */ - public abstract class AbstractPathDependentProduct extends AbstractTreeProduct { - /** Payoff function applied where appropriate (e.g. on a running average A). */ - private final DoubleUnaryOperator payOffFunction; - - /** - * Creates a path–dependent product. - * - * @param maturity Contract maturity in model time units. - * @param payOffFunction Payoff function (e.g., a -> Math.max(a-K, 0.0) for an Asian call). - */ - public AbstractPathDependentProduct(double maturity, DoubleUnaryOperator payOffFunction) { - super(maturity); - this.payOffFunction = payOffFunction; - } - - /** Getter for the payoff function used by concrete implementations. */ - protected DoubleUnaryOperator getPayOffFunction(){ - return payOffFunction; - } - - /** - * Maps a (continuous) time to the closest model time index using the model time step. - * - * @param maturity A time on the model grid (e.g., evaluation time or maturity). - * @param model The tree model providing the time step. - * @return The integer time index k = maturity / dt. - */ - protected final int timeToIndex(double maturity, TreeModel model) { - return (int) Math.round(maturity / model.getTimeStep()); - } - - /** - * Convenience scalar price at t = 0. - * - * @param model The tree model. - * @return The (scalar) price at time 0 obtained from getValue(double, TreeModel). - */ - public double getValue(TreeModel model ){ - RandomVariable v0 = getValue(0.0,model); - return v0.get(0); - } - - - /** - * Convenience accessor returning the value vector at a given evaluation time. - * This delegates to getValues(double, TreeModel) and returns its first element, - * which – by convention – corresponds to the evaluation time. - * - * @param evaluationTime Evaluation time (usually 0.0). - * @param model The tree model. - * @return The RandomVariable holding the option value across states at evalutationTime. - */ - public RandomVariable getValue(double evaluationTime ,TreeModel model) { - RandomVariable[] levels = getValues(evaluationTime,model); - return levels[0]; - } - - /** - * Computes the value process on the lattice from evaluationTime to maturity. - * Implementations : - * Build the path–dependent state (running average) where needed, Apply the payoff at maturity, - * Propagate values backward using the model’s discounted conditional expectation, possibly - * with approximations/projections specific to the product. - * @param evaluationTime The time at which the value is requested. - * @param model The tree model. - * @return An array of RandomVariable from evaluationTime (index 0) to maturity (last index). - */ - public abstract RandomVariable[] getValues(double evaluationTime, TreeModel model); + /** + * Fixing time indices on the model grid, e.g. {1,2,3,4,5}. + * These determine when path-state should be updated with spot. + */ + private final int[] fixingTimeIndices; + + + protected AbstractPathDependentProduct( + final double maturity, + final int[] fixingTimeIndices + ) { + super(maturity); + this.fixingTimeIndices = fixingTimeIndices != null ? fixingTimeIndices.clone() : new int[0]; + Arrays.sort(this.fixingTimeIndices); + } + + /** + * Number of (double-valued) path-state variables required by the product. + * Examples: + * - Fixed-strike Asian: 2 variables (sum, count) + * - Lookback (min): 1 variable (runningMin) + * - Lookback (min/max): 2 variables + */ + protected abstract int getNumberOfStateVariables(); + + /** + * Initialize product path-state at the root node. + * + * @param spotAtNode spot at (timeIndex, recombStateIndex) + * @param timeIndex time index for root (typically k0) + * @param isFixing whether this time index is part of averaging/fixing set + * @param stateOut output array length = getNumberOfStateVariables() + */ + protected abstract void initializeState( + double spotAtNode, + int timeIndex, + boolean isFixing, + double[] stateOut + ); + + /** + * Evolve product path-state from parent to a child node. + * + * @param parentState state at parent node + * @param spotAtChild spot at child node + * @param timeIndexChild time index of child + * @param isFixingChild whether child time is fixing + * @param childStateOut output state array (length = getNumberOfStateVariables()) + */ + protected abstract void evolveState( + double[] parentState, + double spotAtChild, + int timeIndexChild, + boolean isFixingChild, + double[] childStateOut + ); + + /** + * Terminal payoff at maturity. + * + * @param terminalState state variables at maturity + * @param spotAtMaturity spot at maturity node (often needed e.g. floating-strike Asian, lookback) + * @return payoff + */ + protected abstract double payoff( + double[] terminalState, + double spotAtMaturity + ); + + /** + * Converts a double time to a time index by rounding. + */ + protected final int timeToIndex(double time, TreeModel model) { + return (int) Math.round(time / model.getTimeStep()); + } + + /** + * Returns whether a given time index is a fixing time. + */ + protected final boolean isFixingTimeIndex(int timeIndex) { + return Arrays.binarySearch(fixingTimeIndices, timeIndex) >= 0; + } + + /** + * Determine branching factor and childStateShift convention. + */ + private int[] resolveChildShift(TreeModel model, int timeIndex, int stateIndex) { + int branchingFactor = model.getNumberOfBranches(timeIndex, stateIndex); + int[] childStateShift = model.getChildStateIndexShift(); + if(childStateShift != null) { + if(childStateShift.length != branchingFactor) { + throw new IllegalArgumentException( + "childStateShift length (" + childStateShift.length + ") " + + "does not match branching factor (" + branchingFactor + ")." + ); + } + return childStateShift; + } + + throw new IllegalArgumentException( + "No childStateShift provided for branching factor " + branchingFactor + + "." + ); + } + + /** + * Full non-recombining valuation. + * + * Notes about evaluationTime: + * - Path dependent contracts generally require the past path-state at evaluationTime. + * - This implementation assumes the path-state STARTS at evaluationTime (i.e. no past). + * This is exactly what you want when evaluationTime = 0. + * + * IMPORTANT WARNING: + * + * In the current structure of the library this method will always be called with evaluationTime = 0. + * When evaluationTime > 0. the logic contained in AbstractTreeProduct guarantees that the whole tree of the + * path dependent functional is reconstructed starting from time zero. After this the time instants before the + * requested time are discarded and the RandomVariable at the selected time is returned. Remember that we enjoy the + * Markov property only by extending the state space, i.e. with respect to the couple (stock, pathDependentFunctional) + */ + @Override + public final RandomVariable[] getValues(final double evaluationTime, final TreeModel model) { + + // Ensure transition API is supported + // (if not, TreeModel default methods will throw). + final int k0 = timeToIndex(evaluationTime, model); + final int n = timeToIndex(getMaturity(), model); + if(n < k0) throw new IllegalArgumentException("Maturity is before evaluation time."); + + // Root uses recombining state index 0 (consistent with European products). + final int rootRecombState = 0; + + // Determine branching and child shift convention at the root step. + // We assume constant branching factor. + final int[] childShift = resolveChildShift(model, k0, rootRecombState); + final int B = childShift.length; + + final int steps = n - k0; + + // levels[j] corresponds to timeIndex = k0 + j, and has B^j nodes. + final RandomVariable[] levels = new RandomVariable[steps + 1]; + + // Store recombining index per full node, and product state variables per node. + // recombIndex[j][node] = recombining state index at model time k0+j. + final int[][] recombIndex = new int[steps + 1][]; + final double[][] state = new double[steps + 1][]; // flattened: node-major, then var-major + + final int m = getNumberOfStateVariables(); + + // Allocate level 0 + recombIndex[0] = new int[] { rootRecombState }; + state[0] = new double[m]; // 1 node * m vars + + // Spot at root + final double spot0 = model.getSpotAtGivenTimeIndexRV(k0).get(rootRecombState); + + initializeState(spot0, k0, isFixingTimeIndex(k0), state[0]); + + // Forward expansion of path-state (full tree) + for(int j = 1; j <= steps; j++) { + final int timeIndex = k0 + j; + + final int nodes = powInt(B, j); + final int parentNodes = powInt(B, j - 1); + + recombIndex[j] = new int[nodes]; + state[j] = new double[nodes * m]; + + // Pre-fetch the recombining spot vector at this time index. + // We will index into it using recombIndex[j][node]. + final var spotRV = model.getSpotAtGivenTimeIndexRV(timeIndex); + + for(int parent = 0; parent < parentNodes; parent++) { + + final int parentRecomb = recombIndex[j - 1][parent]; + + // parent state slice + final int parentStateOffset = parent * m; + + for(int b = 0; b < B; b++) { + + final int child = parent * B + b; + + final int childRecomb = parentRecomb + childShift[b]; + recombIndex[j][child] = childRecomb; + + final double spotChild = spotRV.get(childRecomb); + + // child state slice + final int childStateOffset = child * m; + + // evolve state + evolveState( + state[j - 1], + parentStateOffset, + spotChild, + timeIndex, + isFixingTimeIndex(timeIndex), + state[j], + childStateOffset + ); + } + } + } + + // Terminal payoff values at maturity level + final int terminalNodes = powInt(B, steps); + final double[] vTerminal = new double[terminalNodes]; + + final var spotN = model.getSpotAtGivenTimeIndexRV(n); + + for(int node = 0; node < terminalNodes; node++) { + final int recomb = recombIndex[steps][node]; + final double spotAtMaturity = spotN.get(recomb); + + final double[] tmp = new double[m]; + System.arraycopy(state[steps], node * m, tmp, 0, m); + + vTerminal[node] = payoff(tmp, spotAtMaturity); + } + + levels[steps] = new RandomVariableFromDoubleArray(n * model.getTimeStep(), vTerminal); + + // Backward induction on full node set + double[] vNext = vTerminal; + + for(int j = steps - 1; j >= 0; j--) { + + final int timeIndex = k0 + j; + final int nodesHere = powInt(B, j); + + final double[] vHere = new double[nodesHere]; + + final double df = model.getOneStepDiscountFactor(timeIndex); + + for(int node = 0; node < nodesHere; node++) { + + final int parentRecomb = recombIndex[j][node]; + + double sum = 0.0; + + for(int b = 0; b < B; b++) { + + final int child = node * B + b; + + final double p = model.getTransitionProbability(timeIndex, parentRecomb, b); + + sum += p * vNext[child]; + } + + vHere[node] = df * sum; + } + + levels[j] = new RandomVariableFromDoubleArray(timeIndex * model.getTimeStep(), vHere); + vNext = vHere; + } + + return levels; + } + + /** + * Helper: evolveState where parent/child arrays are flattened. + * This avoids allocating per-node double[] objects. + */ + private void evolveState( + final double[] parentStateFlat, + final int parentOffset, + final double spotChild, + final int timeIndexChild, + final boolean isFixingChild, + final double[] childStateFlat, + final int childOffset + ) { + final int m = getNumberOfStateVariables(); + final double[] parent = new double[m]; + System.arraycopy(parentStateFlat, parentOffset, parent, 0, m); + + final double[] child = new double[m]; + evolveState(parent, spotChild, timeIndexChild, isFixingChild, child); + + System.arraycopy(child, 0, childStateFlat, childOffset, m); + } + private static int powInt(int base, int exp) { + if(exp < 0) throw new IllegalArgumentException("exp < 0"); + int r = 1; + for(int i = 0; i < exp; i++) { + r = Math.multiplyExact(r, base); // throws if overflow + } + return r; + } } diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/products/AsianOption.java b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AsianOption.java new file mode 100644 index 0000000000..cacaebf221 --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AsianOption.java @@ -0,0 +1,74 @@ +package net.finmath.tree.assetderivativevaluation.products; + +/** + * Fixed-strike Asian option using FULL non-recombining tree expansion. + * + * State variables: + * - state[0] = runningSum of fixing spots + * - state[1] = fixingCount + * + * Payoff at maturity: + * - max( (sum/count) - K, 0 ) + * + * The averaging can be performed on any subset of the time grid by specifying fixingTimeIndices. + * + * No path compression: number of nodes at level j is B^j. + */ +public class AsianOption extends AbstractPathDependentProduct { + + private final double strike; + + /** + * @param maturity maturity in model time units + * @param strike strike K + * @param fixingTimeIndices subset of model time indices included in the average + * @param childStateShift optional mapping for recombining index evolution; + * pass null for defaults (binomial: {0,1}, trinomial: {0,1,2}) + */ + public AsianOption( + final double maturity, + final double strike, + final int[] fixingTimeIndices + ) { + super(maturity, fixingTimeIndices); + this.strike = strike; + } + + @Override + protected int getNumberOfStateVariables() { + return 2; // sum, count + } + + @Override + protected void initializeState( + final double spotAtNode, + final int timeIndex, + final boolean isFixing, + final double[] stateOut + ) { + stateOut[0] = isFixing ? spotAtNode : 0.0; // sum + stateOut[1] = isFixing ? 1.0 : 0.0; // count stored as double for simplicity + } + + @Override + protected void evolveState( + final double[] parentState, + final double spotAtChild, + final int timeIndexChild, + final boolean isFixingChild, + final double[] childStateOut + ) { + childStateOut[0] = parentState[0] + (isFixingChild ? spotAtChild : 0.0); + childStateOut[1] = parentState[1] + (isFixingChild ? 1.0 : 0.0); + } + + @Override + protected double payoff(final double[] terminalState, final double spotAtMaturity) { + final double sum = terminalState[0]; + final double cnt = terminalState[1]; + + final double avg = (cnt > 0.0) ? (sum / cnt) : 0.0; + + return Math.max(avg - strike, 0.0); + } +} diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/products/FloatingStrikeAsianOption.java b/src/main/java/net/finmath/tree/assetderivativevaluation/products/FloatingStrikeAsianOption.java new file mode 100644 index 0000000000..4df0e8782d --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/FloatingStrikeAsianOption.java @@ -0,0 +1,76 @@ +package net.finmath.tree.assetderivativevaluation.products; + +/** + * Floating-strike Asian option (CALL) on a full non-recombining tree (exponential growth). + * + * Payoff at maturity: + * + * max( S(T) - A , 0 ) + * + * where A is the arithmetic average of the spot over the fixingTimeIndices. + * + * State variables: + * - state[0] = running sum of fixing spots + * - state[1] = fixing count + * + * No recombination / compression: + * - number of nodes at step j is B^j (B = branching factor of the model). + * + * @author Alessandro Gnoatto + */ +public class FloatingStrikeAsianOption extends AbstractPathDependentProduct { + + /** + * @param maturity maturity in model time units + * @param fixingTimeIndices subset of model time indices included in the average + */ + public FloatingStrikeAsianOption( + final double maturity, + final int[] fixingTimeIndices + ) { + super(maturity, fixingTimeIndices); + } + + @Override + protected int getNumberOfStateVariables() { + return 2; // running sum, fixing count + } + + @Override + protected void initializeState( + final double spotAtNode, + final int timeIndex, + final boolean isFixing, + final double[] stateOut + ) { + stateOut[0] = isFixing ? spotAtNode : 0.0; // sum + stateOut[1] = isFixing ? 1.0 : 0.0; // count (stored as double) + } + + @Override + protected void evolveState( + final double[] parentState, + final double spotAtChild, + final int timeIndexChild, + final boolean isFixingChild, + final double[] childStateOut + ) { + childStateOut[0] = parentState[0] + (isFixingChild ? spotAtChild : 0.0); + childStateOut[1] = parentState[1] + (isFixingChild ? 1.0 : 0.0); + } + + @Override + protected double payoff( + final double[] terminalState, + final double spotAtMaturity + ) { + final double sum = terminalState[0]; + final double cnt = terminalState[1]; + + // Defensive guard: empty fixing set + final double average = (cnt > 0.0) ? (sum / cnt) : 0.0; + + // Floating-strike CALL payoff + return Math.max(spotAtMaturity - average, 0.0); + } +} diff --git a/src/test/java/net/finmath/tree/PathDependentOptionTreeTest.java b/src/test/java/net/finmath/tree/PathDependentOptionTreeTest.java new file mode 100644 index 0000000000..a101995dd5 --- /dev/null +++ b/src/test/java/net/finmath/tree/PathDependentOptionTreeTest.java @@ -0,0 +1,130 @@ +package net.finmath.tree; + +import java.util.stream.IntStream; + +import org.junit.Assert; +import org.junit.Test; + +import net.finmath.exception.CalculationException; +import net.finmath.montecarlo.BrownianMotionFromMersenneRandomNumbers; +import net.finmath.montecarlo.assetderivativevaluation.AssetModelMonteCarloSimulationModel; +import net.finmath.montecarlo.assetderivativevaluation.MonteCarloAssetModel; +import net.finmath.montecarlo.assetderivativevaluation.models.BlackScholesModel; +import net.finmath.montecarlo.assetderivativevaluation.products.AsianOption; // Monte Carlo Asian +import net.finmath.montecarlo.process.EulerSchemeFromProcessModel; +import net.finmath.montecarlo.process.MonteCarloProcessFromProcessModel; +import net.finmath.stochastic.RandomVariable; +import net.finmath.time.TimeDiscretization; +import net.finmath.time.TimeDiscretizationFromArray; + +// Tree side: +import net.finmath.tree.assetderivativevaluation.models.CoxRossRubinsteinModel; + + +/** + * Compares a fixed-strike Asian option price between + * - Monte Carlo (Black-Scholes via Euler scheme) + * - Binomial tree (CRR) with a FULL non-recombining Asian product (exponential node growth). + * + * The tree AsianOption is assumed to be your class: + * net.finmath.tree.assetderivativevaluation.products.AsianOption + * + * and it is assumed to support: + * - fixed strike call/put + * - a subset of fixing time indices (or times -> indices). + */ +public class PathDependentOptionTreeTest { + + @Test + public void testAsianOptionTreeVsMonteCarlo() throws CalculationException { + + // ------------------------- + // Common contract parameters + // ------------------------- + final double S0 = 100.0; + final double r = 0.04; + final double sigma = 0.25; + + final double maturity = 1.0; + final double strike = 90.0; + + // ------------------------- + // Discretization + // ------------------------- + // Keep time steps modest because your tree Asian is exponential: 2^N nodes at maturity. + final int numberOfTimeSteps = 12; // => 2^12 = 4096 terminal nodes + final double dt = maturity / numberOfTimeSteps; + + // Averaging schedule (subset of the grid): + // Example: fixings from time index 1..N inclusive (exclude t=0). + final int[] fixingTimeIndices = IntStream.rangeClosed(0, numberOfTimeSteps).toArray(); + + // Convert fixing indices to fixing times for the MC product + final double[] fixingTimes = IntStream.of(fixingTimeIndices) + .mapToDouble(k -> k * dt) + .toArray(); + + // ------------------------- + // Monte Carlo pricing + // ------------------------- + final int numberOfPaths = 200_000; + final int seed = 31415; + + final TimeDiscretization timeDiscretization = + new TimeDiscretizationFromArray(0.0, numberOfTimeSteps, dt); + + final BlackScholesModel bsModel = new BlackScholesModel(S0, r, sigma); + + final MonteCarloProcessFromProcessModel process = + new EulerSchemeFromProcessModel( + bsModel, + new BrownianMotionFromMersenneRandomNumbers( + timeDiscretization, + 1, // numberOfFactors + numberOfPaths, + seed + ) + ); + + final AssetModelMonteCarloSimulationModel mcSimulation = + new MonteCarloAssetModel(bsModel, process); + + final AsianOption mcAsian = + new AsianOption( + maturity, + strike, + new TimeDiscretizationFromArray(fixingTimes) + ); + + final RandomVariable mcValueRV = mcAsian.getValue(0.0, mcSimulation); + final double mcValue = mcValueRV.getAverage(); + + // ------------------------- + // Tree pricing (CRR + AsianOption) + // ------------------------- + final TreeModel treeModel = + new CoxRossRubinsteinModel(S0, r, sigma, maturity, dt); + + final net.finmath.tree.assetderivativevaluation.products.AsianOption treeAsian = + new net.finmath.tree.assetderivativevaluation.products.AsianOption( + maturity, + strike, + fixingTimeIndices + ); + + final double treeValue = treeAsian.getValue(treeModel); + + // ------------------------- + // Comparison + // ------------------------- + System.out.println("Asian option (MC BS) : " + mcValue); + System.out.println("Asian option (Tree CRR): " + treeValue); + System.out.println("Abs difference : " + Math.abs(mcValue - treeValue)); + + // Tolerance: MC error ~ O(1/sqrt(paths)). With 200k paths, this is typically a few 1e-3..1e-2. + // Tree has discretization error too (CRR + dt). + final double tolerance = 1.0e-1; + + Assert.assertEquals("Tree vs Monte Carlo price mismatch", mcValue, treeValue, tolerance); + } +} From a6f823efa4f7447b44f75a6a2593a04620d1077f Mon Sep 17 00:00:00 2001 From: Alessandro Gnoatto Date: Wed, 18 Feb 2026 23:04:42 +0100 Subject: [PATCH 5/7] Fix around 1300 checkstyle issues --- .../marketdata/AffineDividendStream.java | 172 ++-- .../equities/marketdata/FlatYieldCurve.java | 35 +- .../equities/marketdata/YieldCurve.java | 142 ++-- .../BuehlerDividendForwardStructure.java | 297 +++---- .../models/EquityForwardStructure.java | 49 +- .../pricer/AnalyticOptionValuation.java | 255 +++--- .../equities/pricer/PdeOptionValuation.java | 765 ++++++++++-------- .../net/finmath/functions/HestonModel.java | 578 ++++++------- .../net/finmath/tree/AbstractTreeProduct.java | 86 +- .../OneDimensionalRiskFactorTreeModel.java | 81 +- .../AbstractRecombiningTreeModel.java | 4 +- .../OneDimensionalAssetTreeModel.java | 24 +- .../models/CoxRossRubinsteinModel.java | 2 +- .../models/package-info.java | 6 + .../package-info.java | 6 + .../AbstractPathDependentProduct.java | 507 ++++++------ .../products/AsianOption.java | 103 +-- .../products/FloatingStrikeAsianOption.java | 113 +-- .../products/package-info.java | 6 + .../java/net/finmath/tree/package-info.java | 6 + 20 files changed, 1699 insertions(+), 1538 deletions(-) create mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/models/package-info.java create mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/package-info.java create mode 100644 src/main/java/net/finmath/tree/assetderivativevaluation/products/package-info.java create mode 100644 src/main/java/net/finmath/tree/package-info.java diff --git a/src/main/java/net/finmath/equities/marketdata/AffineDividendStream.java b/src/main/java/net/finmath/equities/marketdata/AffineDividendStream.java index d265a4079a..142306130b 100644 --- a/src/main/java/net/finmath/equities/marketdata/AffineDividendStream.java +++ b/src/main/java/net/finmath/equities/marketdata/AffineDividendStream.java @@ -11,86 +11,98 @@ * * @author Andreas Grotz */ - public class AffineDividendStream { - private final AffineDividend[] dividendStream; - - public AffineDividendStream(final AffineDividend[] dividendStream) { - final var diviList = Arrays.asList(dividendStream); - diviList.sort(Comparator.comparing(pt -> pt.getDate())); - this.dividendStream = diviList.toArray(new AffineDividend[0]); - } - - public ArrayList getDividendDates() { - final var dates = new ArrayList(); - for (final AffineDividend divi : dividendStream) { - dates.add(divi.getDate()); - } - return dates; - } - - public double getDividend(final LocalDate date, final double stockPrice) { - for (final AffineDividend divi : dividendStream) { - if (divi.getDate() == date) { - return divi.getDividend(stockPrice); - } - } - return 0.0; - } - - public double getProportionalDividendFactor(final LocalDate date) { - for (final AffineDividend divi : dividendStream) { - if (divi.getDate() == date) { - return divi.getProportionalDividendFactor(); - } - } - return 1.0; - } - - public double getCashDividend(final LocalDate date) { - for (final AffineDividend divi : dividendStream) { - if (divi.getDate() == date) { - return divi.getCashDividend(); - } - } - return 0.0; - } - - public static AffineDividendStream getAffineDividendsFromCashDividends(AffineDividendStream cashDividends, - HashMap transformationFactors, LocalDate valDate, double spot, YieldCurve repoCurve) { - // This method takes a stream of cash dividends and converts them to affine dividends, - // by transforming a part of each cash dividend to a proportional dividend. - // The percentage of each cash dividend to be transformed to a proportional dividend - // is specified in the member propDividendFactor of the dividend. - // The transformation is done in an arbitrage-free way, i.e. the forward structure is preserved. - // This method is usefull in practice, where traders use dividend futures as input, and transform - // a part to a proportional dividend (the further away the dividend, the higher the proportional part - // and the lower the cash part. - - final var dates = cashDividends.getDividendDates(); - - final var affineDividends = new ArrayList(); - - for (final var date : dates) { - if (date.isBefore(valDate)) { - continue; - } - assert cashDividends.getProportionalDividendFactor( - date) == 0.0 : "Proportional dividend different from zero for date " + date; - final var cashDividend = cashDividends.getCashDividend(date); - var fwd = spot; - for (final var otherDate : dates) { - if (otherDate.isBefore(date) && !otherDate.isBefore(valDate)) { - fwd -= cashDividends.getCashDividend(otherDate) - * repoCurve.getForwardDiscountFactor(valDate, otherDate); - } - } - final var q = transformationFactors.get(date) * cashDividend - * repoCurve.getForwardDiscountFactor(valDate, date) / fwd; - affineDividends.add(new AffineDividend(date, (1.0 - transformationFactors.get(date)) * cashDividend, q)); - } - - return new AffineDividendStream(affineDividends.toArray(new AffineDividend[0])); - } + private final AffineDividend[] dividendStream; + + public AffineDividendStream(final AffineDividend[] dividendStream) { + final var diviList = Arrays.asList(dividendStream); + diviList.sort(Comparator.comparing(pt -> pt.getDate())); + this.dividendStream = diviList.toArray(new AffineDividend[0]); + } + + public ArrayList getDividendDates() { + final var dates = new ArrayList(); + for(final AffineDividend divi : dividendStream) { + dates.add(divi.getDate()); + } + return dates; + } + + public double getDividend(final LocalDate date, final double stockPrice) { + for(final AffineDividend divi : dividendStream) { + if(divi.getDate() == date) { + return divi.getDividend(stockPrice); + } + } + return 0.0; + } + + public double getProportionalDividendFactor(final LocalDate date) { + for(final AffineDividend divi : dividendStream) { + if(divi.getDate() == date) { + return divi.getProportionalDividendFactor(); + } + } + return 1.0; + } + + public double getCashDividend(final LocalDate date) { + for(final AffineDividend divi : dividendStream) { + if(divi.getDate() == date) { + return divi.getCashDividend(); + } + } + return 0.0; + } + + public static AffineDividendStream getAffineDividendsFromCashDividends( + AffineDividendStream cashDividends, + HashMap transformationFactors, + LocalDate valDate, + double spot, + YieldCurve repoCurve) { + + // This method takes a stream of cash dividends and converts them to affine dividends, + // by transforming a part of each cash dividend to a proportional dividend. + // The percentage of each cash dividend to be transformed to a proportional dividend + // is specified in the member propDividendFactor of the dividend. + // The transformation is done in an arbitrage-free way, i.e. the forward structure is preserved. + // This method is usefull in practice, where traders use dividend futures as input, and transform + // a part to a proportional dividend (the further away the dividend, the higher the proportional part + // and the lower the cash part. + + final var dates = cashDividends.getDividendDates(); + + final var affineDividends = new ArrayList(); + + for(final var date : dates) { + if(date.isBefore(valDate)) { + continue; + } + assert cashDividends.getProportionalDividendFactor(date) == 0.0 + : "Proportional dividend different from zero for date " + date; + + final var cashDividend = cashDividends.getCashDividend(date); + + var fwd = spot; + for(final var otherDate : dates) { + if(otherDate.isBefore(date) && !otherDate.isBefore(valDate)) { + fwd -= cashDividends.getCashDividend(otherDate) + * repoCurve.getForwardDiscountFactor(valDate, otherDate); + } + } + + final var q = transformationFactors.get(date) * cashDividend + * repoCurve.getForwardDiscountFactor(valDate, date) / fwd; + + affineDividends.add( + new AffineDividend( + date, + (1.0 - transformationFactors.get(date)) * cashDividend, + q)); + } + + return new AffineDividendStream(affineDividends.toArray(new AffineDividend[0])); + } } diff --git a/src/main/java/net/finmath/equities/marketdata/FlatYieldCurve.java b/src/main/java/net/finmath/equities/marketdata/FlatYieldCurve.java index 18a516768f..873d389c96 100644 --- a/src/main/java/net/finmath/equities/marketdata/FlatYieldCurve.java +++ b/src/main/java/net/finmath/equities/marketdata/FlatYieldCurve.java @@ -1,6 +1,7 @@ package net.finmath.equities.marketdata; import java.time.LocalDate; + import net.finmath.time.daycount.DayCountConvention; /** @@ -8,19 +9,33 @@ * * @author Andreas Grotz */ - public class FlatYieldCurve extends YieldCurve { - private final static int longTime = 100; + private static final int LONG_TIME = 100; - public FlatYieldCurve(final LocalDate curveDate, final double rate, final DayCountConvention dayCounter) { - super("NONE", curveDate, dayCounter, new LocalDate[] { curveDate.plusYears(longTime) }, new double[] { - Math.exp(-rate * dayCounter.getDaycountFraction(curveDate, curveDate.plusYears(longTime))) }); - } + public FlatYieldCurve( + final LocalDate curveDate, + final double rate, + final DayCountConvention dayCounter) { - public FlatYieldCurve rollToDate(LocalDate date) { - assert date.isAfter(baseCurve.getReferenceDate()) : "can only roll to future dates"; - return new FlatYieldCurve(date, getRate(baseCurve.getReferenceDate().plusYears(longTime)), dayCounter); - } + super( + "NONE", + curveDate, + dayCounter, + new LocalDate[] { curveDate.plusYears(LONG_TIME) }, + new double[] { + Math.exp( + -rate * dayCounter.getDaycountFraction( + curveDate, + curveDate.plusYears(LONG_TIME))) + }); + } + public FlatYieldCurve rollToDate(LocalDate date) { + assert date.isAfter(baseCurve.getReferenceDate()) : "can only roll to future dates"; + return new FlatYieldCurve( + date, + getRate(baseCurve.getReferenceDate().plusYears(LONG_TIME)), + dayCounter); + } } diff --git a/src/main/java/net/finmath/equities/marketdata/YieldCurve.java b/src/main/java/net/finmath/equities/marketdata/YieldCurve.java index 91cb81d91a..9d51df4449 100644 --- a/src/main/java/net/finmath/equities/marketdata/YieldCurve.java +++ b/src/main/java/net/finmath/equities/marketdata/YieldCurve.java @@ -14,69 +14,85 @@ * * @author Andreas Grotz */ - public class YieldCurve { - protected final LocalDate referenceDate; - protected final LocalDate[] discountDates; - protected final DayCountConvention dayCounter; - protected final DiscountCurveInterpolation baseCurve; - - public YieldCurve(final String name, final LocalDate referenceDate, final DayCountConvention dayCounter, - final LocalDate[] discountDates, final double[] discountFactors) { - this.dayCounter = dayCounter; - this.discountDates = discountDates; - double[] times = new double[discountDates.length]; - boolean[] isParameter = new boolean[discountDates.length]; - for (int i = 0; i < times.length; i++) { - times[i] = dayCounter.getDaycountFraction(referenceDate, discountDates[i]); - } - baseCurve = DiscountCurveInterpolation.createDiscountCurveFromDiscountFactors(name, referenceDate, times, - discountFactors, isParameter, InterpolationMethod.LINEAR, ExtrapolationMethod.CONSTANT, - InterpolationEntity.LOG_OF_VALUE_PER_TIME); - - this.referenceDate = referenceDate; - } - - public YieldCurve rollToDate(LocalDate date) { - assert date.isAfter(referenceDate) : "can only roll to future dates"; - LocalDate[] rolledDiscountDates = Arrays.stream(discountDates).filter(p -> p.isAfter(date)) - .toArray(LocalDate[]::new); - double[] rolledDiscountFactors = new double[rolledDiscountDates.length]; - for (int i = 0; i < rolledDiscountDates.length; i++) { - rolledDiscountFactors[i] = getForwardDiscountFactor(date, rolledDiscountDates[i]); - } - - return new YieldCurve(baseCurve.getName(), date, dayCounter, rolledDiscountDates, rolledDiscountFactors); - } - - public double getRate(double maturity) { - assert maturity >= 0.0 : "maturity must be positive"; - return baseCurve.getZeroRate(maturity); - } - - public double getRate(LocalDate date) { - return baseCurve.getZeroRate(dayCounter.getDaycountFraction(referenceDate, date)); - } - - public double getDiscountFactor(double maturity) { - assert maturity >= 0.0 : "maturity must be positive"; - return baseCurve.getDiscountFactor(maturity); - } - - public double getForwardDiscountFactor(double start, double expiry) { - assert start >= 0.0 : "start must be positive"; - assert expiry >= start : "start must be before expiry"; - return getDiscountFactor(expiry) / getDiscountFactor(start); - } - - public double getDiscountFactor(LocalDate date) { - return baseCurve.getDiscountFactor(dayCounter.getDaycountFraction(referenceDate, date)); - } - - public double getForwardDiscountFactor(LocalDate startDate, LocalDate endDate) { - assert !startDate.isBefore(referenceDate) : "start date must be after curve date"; - assert !endDate.isBefore(startDate) : "end date must be after start date"; - return getDiscountFactor(endDate) / getDiscountFactor(startDate); - } + protected final LocalDate referenceDate; + protected final LocalDate[] discountDates; + protected final DayCountConvention dayCounter; + protected final DiscountCurveInterpolation baseCurve; + + public YieldCurve( + final String name, + final LocalDate referenceDate, + final DayCountConvention dayCounter, + final LocalDate[] discountDates, + final double[] discountFactors) { + + this.dayCounter = dayCounter; + this.discountDates = discountDates; + + double[] times = new double[discountDates.length]; + boolean[] isParameter = new boolean[discountDates.length]; + + for(int i = 0; i < times.length; i++) { + times[i] = dayCounter.getDaycountFraction(referenceDate, discountDates[i]); + } + + baseCurve = DiscountCurveInterpolation.createDiscountCurveFromDiscountFactors( + name, + referenceDate, + times, + discountFactors, + isParameter, + InterpolationMethod.LINEAR, + ExtrapolationMethod.CONSTANT, + InterpolationEntity.LOG_OF_VALUE_PER_TIME); + + this.referenceDate = referenceDate; + } + + public YieldCurve rollToDate(LocalDate date) { + assert date.isAfter(referenceDate) : "can only roll to future dates"; + + LocalDate[] rolledDiscountDates = Arrays.stream(discountDates) + .filter(p -> p.isAfter(date)) + .toArray(LocalDate[]::new); + + double[] rolledDiscountFactors = new double[rolledDiscountDates.length]; + for(int i = 0; i < rolledDiscountDates.length; i++) { + rolledDiscountFactors[i] = getForwardDiscountFactor(date, rolledDiscountDates[i]); + } + + return new YieldCurve(baseCurve.getName(), date, dayCounter, rolledDiscountDates, rolledDiscountFactors); + } + + public double getRate(double maturity) { + assert maturity >= 0.0 : "maturity must be positive"; + return baseCurve.getZeroRate(maturity); + } + + public double getRate(LocalDate date) { + return baseCurve.getZeroRate(dayCounter.getDaycountFraction(referenceDate, date)); + } + + public double getDiscountFactor(double maturity) { + assert maturity >= 0.0 : "maturity must be positive"; + return baseCurve.getDiscountFactor(maturity); + } + + public double getForwardDiscountFactor(double start, double expiry) { + assert start >= 0.0 : "start must be positive"; + assert expiry >= start : "start must be before expiry"; + return getDiscountFactor(expiry) / getDiscountFactor(start); + } + + public double getDiscountFactor(LocalDate date) { + return baseCurve.getDiscountFactor(dayCounter.getDaycountFraction(referenceDate, date)); + } + + public double getForwardDiscountFactor(LocalDate startDate, LocalDate endDate) { + assert !startDate.isBefore(referenceDate) : "start date must be after curve date"; + assert !endDate.isBefore(startDate) : "end date must be after start date"; + return getDiscountFactor(endDate) / getDiscountFactor(startDate); + } } diff --git a/src/main/java/net/finmath/equities/models/BuehlerDividendForwardStructure.java b/src/main/java/net/finmath/equities/models/BuehlerDividendForwardStructure.java index c9b88af5a1..64b99aa438 100644 --- a/src/main/java/net/finmath/equities/models/BuehlerDividendForwardStructure.java +++ b/src/main/java/net/finmath/equities/models/BuehlerDividendForwardStructure.java @@ -2,6 +2,7 @@ import java.time.LocalDate; import java.util.HashMap; + import net.finmath.equities.marketdata.AffineDividendStream; import net.finmath.equities.marketdata.YieldCurve; import net.finmath.time.daycount.DayCountConvention; @@ -12,146 +13,162 @@ * * @author Andreas Grotz */ - public class BuehlerDividendForwardStructure implements EquityForwardStructure { - private final LocalDate valuationDate; - private final double spot; - private final YieldCurve repoCurve; - private final AffineDividendStream dividendStream; - private final DayCountConvention dayCounter; - private final HashMap dividendTimes; - - public BuehlerDividendForwardStructure(final LocalDate valuationDate, final double spot, final YieldCurve repoCurve, - final AffineDividendStream dividendStream, final DayCountConvention dayCounter) { - this.valuationDate = valuationDate; - this.spot = spot; - this.repoCurve = repoCurve; - this.dividendStream = dividendStream; - this.dayCounter = dayCounter; - dividendTimes = new HashMap(); - for (final var date : dividendStream.getDividendDates()) { - dividendTimes.put(date, dayCounter.getDaycountFraction(valuationDate, date)); - } - validate(); - } - - public void validate() { - assert getFutureDividendFactor(valuationDate) <= spot : "PV of future dividends is larger than spot."; - } - - @Override - public BuehlerDividendForwardStructure cloneWithNewSpot(double newSpot) { - return new BuehlerDividendForwardStructure(this.valuationDate, newSpot, this.repoCurve, this.dividendStream, - this.dayCounter); - } - - @Override - public BuehlerDividendForwardStructure cloneWithNewDate(LocalDate newDate) { - return new BuehlerDividendForwardStructure(newDate, this.spot, this.repoCurve.rollToDate(newDate), - this.dividendStream, this.dayCounter); - } - - @Override - public DividendModelType getDividendModel() { - return DividendModelType.Buehler; - } - - @Override - public LocalDate getValuationDate() { - return valuationDate; - } - - @Override - public double getSpot() { - return spot; - } - - @Override - public YieldCurve getRepoCurve() { - return repoCurve; - } - - @Override - public AffineDividendStream getDividendStream() { - return dividendStream; - } - - @Override - public double getGrowthDiscountFactor(double startTime, double endTime) { - var df = 1.0; - for (final var date : dividendStream.getDividendDates()) { - final var dividendTime = dividendTimes.get(date); - if (dividendTime > startTime && dividendTime <= endTime) { - df *= (1.0 - dividendStream.getProportionalDividendFactor(date)); - } - } - - return df / repoCurve.getForwardDiscountFactor(startTime, endTime); - } - - @Override - public double getGrowthDiscountFactor(LocalDate startDate, LocalDate endDate) { - final var startTime = dayCounter.getDaycountFraction(valuationDate, startDate); - final var endTime = dayCounter.getDaycountFraction(valuationDate, endDate); - return getGrowthDiscountFactor(startTime, endTime); - } - - @Override - public double getFutureDividendFactor(double valTime) { - var df = 0.0; - for (final var date : dividendStream.getDividendDates()) { - final var dividendTime = dividendTimes.get(date); - if (dividendTime > valTime) { - df += dividendStream.getCashDividend(date) / getGrowthDiscountFactor(valTime, dividendTime); - } - } - return df; - } - - @Override - public double getFutureDividendFactor(LocalDate valDate) { - final var valTime = dayCounter.getDaycountFraction(valuationDate, valDate); - return getFutureDividendFactor(valTime); - } - - @Override - public double getForward(double expiryTime) { - var forward = spot * getGrowthDiscountFactor(0.0, expiryTime); - for (final var date : dividendStream.getDividendDates()) { - final var dividendTime = dividendTimes.get(date); - if (dividendTime <= expiryTime) { - forward -= dividendStream.getCashDividend(date) * getGrowthDiscountFactor(dividendTime, expiryTime); - } - } - return forward; - } - - @Override - public double getForward(LocalDate expiryDate) { - final var expiryTime = dayCounter.getDaycountFraction(valuationDate, expiryDate); - return getForward(expiryTime); - } - - @Override - public double getDividendAdjustedStrike(double strike, double expiryTime) { - return strike - getFutureDividendFactor(expiryTime); - } - - @Override - public double getDividendAdjustedStrike(double strike, LocalDate expiryDate) { - return strike - getFutureDividendFactor(expiryDate); - } - - @Override - public double getLogMoneyness(double strike, double expiryTime) { - return Math.log(getDividendAdjustedStrike(strike, expiryTime) - / getDividendAdjustedStrike(getForward(expiryTime), expiryTime)); - } - - @Override - public double getLogMoneyness(double strike, LocalDate expiryDate) { - return Math.log(getDividendAdjustedStrike(strike, expiryDate) - / getDividendAdjustedStrike(getForward(expiryDate), expiryDate)); - } + private final LocalDate valuationDate; + private final double spot; + private final YieldCurve repoCurve; + private final AffineDividendStream dividendStream; + private final DayCountConvention dayCounter; + private final HashMap dividendTimes; + + public BuehlerDividendForwardStructure( + final LocalDate valuationDate, + final double spot, + final YieldCurve repoCurve, + final AffineDividendStream dividendStream, + final DayCountConvention dayCounter) { + + this.valuationDate = valuationDate; + this.spot = spot; + this.repoCurve = repoCurve; + this.dividendStream = dividendStream; + this.dayCounter = dayCounter; + + dividendTimes = new HashMap(); + for(final var date : dividendStream.getDividendDates()) { + dividendTimes.put(date, dayCounter.getDaycountFraction(valuationDate, date)); + } + + validate(); + } + + public void validate() { + assert getFutureDividendFactor(valuationDate) <= spot : "PV of future dividends is larger than spot."; + } + + @Override + public BuehlerDividendForwardStructure cloneWithNewSpot(double newSpot) { + return new BuehlerDividendForwardStructure( + this.valuationDate, + newSpot, + this.repoCurve, + this.dividendStream, + this.dayCounter); + } + + @Override + public BuehlerDividendForwardStructure cloneWithNewDate(LocalDate newDate) { + return new BuehlerDividendForwardStructure( + newDate, + this.spot, + this.repoCurve.rollToDate(newDate), + this.dividendStream, + this.dayCounter); + } + + @Override + public DividendModelType getDividendModel() { + return DividendModelType.Buehler; + } + + @Override + public LocalDate getValuationDate() { + return valuationDate; + } + + @Override + public double getSpot() { + return spot; + } + + @Override + public YieldCurve getRepoCurve() { + return repoCurve; + } + + @Override + public AffineDividendStream getDividendStream() { + return dividendStream; + } + + @Override + public double getGrowthDiscountFactor(double startTime, double endTime) { + var df = 1.0; + for(final var date : dividendStream.getDividendDates()) { + final var dividendTime = dividendTimes.get(date); + if(dividendTime > startTime && dividendTime <= endTime) { + df *= (1.0 - dividendStream.getProportionalDividendFactor(date)); + } + } + + return df / repoCurve.getForwardDiscountFactor(startTime, endTime); + } + + @Override + public double getGrowthDiscountFactor(LocalDate startDate, LocalDate endDate) { + final var startTime = dayCounter.getDaycountFraction(valuationDate, startDate); + final var endTime = dayCounter.getDaycountFraction(valuationDate, endDate); + return getGrowthDiscountFactor(startTime, endTime); + } + + @Override + public double getFutureDividendFactor(double valTime) { + var df = 0.0; + for(final var date : dividendStream.getDividendDates()) { + final var dividendTime = dividendTimes.get(date); + if(dividendTime > valTime) { + df += dividendStream.getCashDividend(date) / getGrowthDiscountFactor(valTime, dividendTime); + } + } + return df; + } + + @Override + public double getFutureDividendFactor(LocalDate valDate) { + final var valTime = dayCounter.getDaycountFraction(valuationDate, valDate); + return getFutureDividendFactor(valTime); + } + + @Override + public double getForward(double expiryTime) { + var forward = spot * getGrowthDiscountFactor(0.0, expiryTime); + for(final var date : dividendStream.getDividendDates()) { + final var dividendTime = dividendTimes.get(date); + if(dividendTime <= expiryTime) { + forward -= dividendStream.getCashDividend(date) * getGrowthDiscountFactor(dividendTime, expiryTime); + } + } + return forward; + } + + @Override + public double getForward(LocalDate expiryDate) { + final var expiryTime = dayCounter.getDaycountFraction(valuationDate, expiryDate); + return getForward(expiryTime); + } + + @Override + public double getDividendAdjustedStrike(double strike, double expiryTime) { + return strike - getFutureDividendFactor(expiryTime); + } + + @Override + public double getDividendAdjustedStrike(double strike, LocalDate expiryDate) { + return strike - getFutureDividendFactor(expiryDate); + } + + @Override + public double getLogMoneyness(double strike, double expiryTime) { + return Math.log( + getDividendAdjustedStrike(strike, expiryTime) + / getDividendAdjustedStrike(getForward(expiryTime), expiryTime)); + } + + @Override + public double getLogMoneyness(double strike, LocalDate expiryDate) { + return Math.log( + getDividendAdjustedStrike(strike, expiryDate) + / getDividendAdjustedStrike(getForward(expiryDate), expiryDate)); + } } diff --git a/src/main/java/net/finmath/equities/models/EquityForwardStructure.java b/src/main/java/net/finmath/equities/models/EquityForwardStructure.java index 493bba891a..d977d9b9c5 100644 --- a/src/main/java/net/finmath/equities/models/EquityForwardStructure.java +++ b/src/main/java/net/finmath/equities/models/EquityForwardStructure.java @@ -6,11 +6,12 @@ package net.finmath.equities.models; import java.time.LocalDate; + import net.finmath.equities.marketdata.AffineDividendStream; import net.finmath.equities.marketdata.YieldCurve; /** - * I to cover the forward structure of a stock, i.e. spot, repo curve and dividends. + * Interface to cover the forward structure of a stock, i.e. spot, repo curve and dividends. * Currently implemented is the Buehler dividend model, where the volatile part of the stock * excludes any future dividends. * @@ -26,44 +27,44 @@ */ public interface EquityForwardStructure { - enum DividendModelType { - None, - Buehler, - Escrowed, - HaugHaugLewis, - } + enum DividendModelType { + None, + Buehler, + Escrowed, + HaugHaugLewis, + } - DividendModelType getDividendModel(); + DividendModelType getDividendModel(); - LocalDate getValuationDate(); + LocalDate getValuationDate(); - double getSpot(); + double getSpot(); - YieldCurve getRepoCurve(); + YieldCurve getRepoCurve(); - AffineDividendStream getDividendStream(); + AffineDividendStream getDividendStream(); - EquityForwardStructure cloneWithNewSpot(double newSpot); + EquityForwardStructure cloneWithNewSpot(double newSpot); - EquityForwardStructure cloneWithNewDate(LocalDate newDate); + EquityForwardStructure cloneWithNewDate(LocalDate newDate); - double getGrowthDiscountFactor(double startTime, double endTime); + double getGrowthDiscountFactor(double startTime, double endTime); - double getGrowthDiscountFactor(LocalDate startDate, LocalDate endDate); + double getGrowthDiscountFactor(LocalDate startDate, LocalDate endDate); - double getFutureDividendFactor(double valTime); + double getFutureDividendFactor(double valTime); - double getFutureDividendFactor(LocalDate valDate); + double getFutureDividendFactor(LocalDate valDate); - double getForward(double expiryTime); + double getForward(double expiryTime); - double getForward(LocalDate expiryDate); + double getForward(LocalDate expiryDate); - double getDividendAdjustedStrike(double strike, double expiryTime); + double getDividendAdjustedStrike(double strike, double expiryTime); - double getDividendAdjustedStrike(double strike, LocalDate expiryDate); + double getDividendAdjustedStrike(double strike, LocalDate expiryDate); - double getLogMoneyness(double strike, double expiryTime); + double getLogMoneyness(double strike, double expiryTime); - double getLogMoneyness(double strike, LocalDate expiryDate); + double getLogMoneyness(double strike, LocalDate expiryDate); } diff --git a/src/main/java/net/finmath/equities/pricer/AnalyticOptionValuation.java b/src/main/java/net/finmath/equities/pricer/AnalyticOptionValuation.java index 8c76ccf185..9b993be292 100644 --- a/src/main/java/net/finmath/equities/pricer/AnalyticOptionValuation.java +++ b/src/main/java/net/finmath/equities/pricer/AnalyticOptionValuation.java @@ -1,7 +1,9 @@ package net.finmath.equities.pricer; import java.util.HashMap; + import org.apache.commons.lang3.NotImplementedException; + import net.finmath.equities.marketdata.YieldCurve; import net.finmath.equities.models.Black76Model; import net.finmath.equities.models.EquityForwardStructure; @@ -16,100 +18,165 @@ * * @author Andreas Grotz */ - public class AnalyticOptionValuation implements OptionValuation { - private final DayCountConvention dcc; - - public AnalyticOptionValuation(DayCountConvention dcc) { - this.dcc = dcc; - } - - @Override - public EquityValuationResult calculate(EquityValuationRequest request, EquityForwardStructure forwardStructure, - YieldCurve discountCurve, VolatilitySurface volaSurface) { - final var results = new HashMap(); - for (final var calcType : request.getCalcsRequested()) { - results.put(calcType, - calculate(request.getOption(), forwardStructure, discountCurve, volaSurface, calcType)); - } - - return new EquityValuationResult(request, results); - } - - public double calculate(Option option, EquityForwardStructure forwardStructure, YieldCurve discountCurve, - VolatilitySurface volaSurface, CalculationRequestType calcType) { - assert !option.isAmericanOption() : "Analytic pricer cannot handle American options."; - final var valDate = forwardStructure.getValuationDate(); - final var expiryDate = option.getExpiryDate(); - final var ttm = dcc.getDaycountFraction(forwardStructure.getValuationDate(), expiryDate); - final var forward = forwardStructure.getForward(expiryDate); - final var discountFactor = discountCurve.getDiscountFactor(expiryDate); - final var discountRate = discountCurve.getRate(expiryDate); - final var adjustedForward = forwardStructure.getDividendAdjustedStrike(forward, expiryDate); - final var adjustedStrike = forwardStructure.getDividendAdjustedStrike(option.getStrike(), expiryDate); - final var volatility = volaSurface.getVolatility(option.getStrike(), option.getExpiryDate(), forwardStructure); - - switch (calcType) { - case Price: - return Black76Model.optionPrice(1.0, adjustedStrike / adjustedForward, ttm, volatility, - option.isCallOption(), discountFactor * adjustedForward); - case EqDelta: - final var dFdS = forwardStructure.getGrowthDiscountFactor(valDate, expiryDate); - return dFdS * Black76Model.optionDelta(1.0, adjustedStrike / adjustedForward, ttm, volatility, - option.isCallOption(), discountFactor); - case EqGamma: - final var dFdS2 = Math.pow(forwardStructure.getGrowthDiscountFactor(valDate, expiryDate), 2); - return dFdS2 * Black76Model.optionGamma(1.0, adjustedStrike / adjustedForward, ttm, volatility, - option.isCallOption(), discountFactor / adjustedForward); - case EqVega: - return Black76Model.optionVega(1.0, adjustedStrike / adjustedForward, ttm, volatility, - option.isCallOption(), discountFactor * adjustedForward); - case Theta: - return Black76Model.optionTheta(1.0, adjustedStrike / adjustedForward, ttm, volatility, - option.isCallOption(), discountFactor * adjustedForward, discountRate); - default: - throw new NotImplementedException("Calculation for " + calcType + " not implemented yet."); - } - } - - public double getPrice(Option option, EquityForwardStructure forwardStructure, YieldCurve discountCurve, - VolatilitySurface volaSurface) { - return calculate(option, forwardStructure, discountCurve, volaSurface, CalculationRequestType.Price); - } - - public double getDelta(Option option, EquityForwardStructure forwardStructure, YieldCurve discountCurve, - VolatilitySurface volaSurface) { - return calculate(option, forwardStructure, discountCurve, volaSurface, CalculationRequestType.EqDelta); - } - - public double getGamma(Option option, EquityForwardStructure forwardStructure, YieldCurve discountCurve, - VolatilitySurface volaSurface) { - return calculate(option, forwardStructure, discountCurve, volaSurface, CalculationRequestType.EqGamma); - } - - public double getVega(Option option, EquityForwardStructure forwardStructure, YieldCurve discountCurve, - VolatilitySurface volaSurface) { - return calculate(option, forwardStructure, discountCurve, volaSurface, CalculationRequestType.EqVega); - } - - public double getTheta(Option option, EquityForwardStructure forwardStructure, YieldCurve discountCurve, - VolatilitySurface volaSurface) { - return calculate(option, forwardStructure, discountCurve, volaSurface, CalculationRequestType.Theta); - } - - public double getImpliedVolatility(Option option, EquityForwardStructure forwardStructure, YieldCurve discountCurve, - double price) { - assert !option.isAmericanOption() : "Analytic pricer cannot handle American options."; - final var expiryDate = option.getExpiryDate(); - final var ttm = dcc.getDaycountFraction(forwardStructure.getValuationDate(), expiryDate); - final var forward = forwardStructure.getForward(expiryDate); - final var discount = discountCurve.getDiscountFactor(expiryDate); - final var adjustedForward = forwardStructure.getDividendAdjustedStrike(forward, expiryDate); - final var adjustedStrike = forwardStructure.getDividendAdjustedStrike(option.getStrike(), expiryDate); - final var undiscountedPrice = price / discount / adjustedForward; - - return Black76Model.optionImpliedVolatility(1.0, adjustedStrike / adjustedForward, ttm, undiscountedPrice, - option.isCallOption()); - } + private final DayCountConvention dcc; + + public AnalyticOptionValuation(DayCountConvention dcc) { + this.dcc = dcc; + } + + @Override + public EquityValuationResult calculate( + EquityValuationRequest request, + EquityForwardStructure forwardStructure, + YieldCurve discountCurve, + VolatilitySurface volaSurface) { + + final var results = new HashMap(); + for(final var calcType : request.getCalcsRequested()) { + results.put( + calcType, + calculate(request.getOption(), forwardStructure, discountCurve, volaSurface, calcType)); + } + + return new EquityValuationResult(request, results); + } + + public double calculate( + Option option, + EquityForwardStructure forwardStructure, + YieldCurve discountCurve, + VolatilitySurface volaSurface, + CalculationRequestType calcType) { + + assert !option.isAmericanOption() : "Analytic pricer cannot handle American options."; + + final var valDate = forwardStructure.getValuationDate(); + final var expiryDate = option.getExpiryDate(); + final var ttm = dcc.getDaycountFraction(forwardStructure.getValuationDate(), expiryDate); + final var forward = forwardStructure.getForward(expiryDate); + final var discountFactor = discountCurve.getDiscountFactor(expiryDate); + final var discountRate = discountCurve.getRate(expiryDate); + final var adjustedForward = forwardStructure.getDividendAdjustedStrike(forward, expiryDate); + final var adjustedStrike = forwardStructure.getDividendAdjustedStrike(option.getStrike(), expiryDate); + final var volatility = volaSurface.getVolatility(option.getStrike(), option.getExpiryDate(), forwardStructure); + + switch(calcType) { + case Price: + return Black76Model.optionPrice( + 1.0, + adjustedStrike / adjustedForward, + ttm, + volatility, + option.isCallOption(), + discountFactor * adjustedForward); + + case EqDelta: + final var dFdS = forwardStructure.getGrowthDiscountFactor(valDate, expiryDate); + return dFdS * Black76Model.optionDelta( + 1.0, + adjustedStrike / adjustedForward, + ttm, + volatility, + option.isCallOption(), + discountFactor); + + case EqGamma: + final var dFdS2 = Math.pow(forwardStructure.getGrowthDiscountFactor(valDate, expiryDate), 2); + return dFdS2 * Black76Model.optionGamma( + 1.0, + adjustedStrike / adjustedForward, + ttm, + volatility, + option.isCallOption(), + discountFactor / adjustedForward); + + case EqVega: + return Black76Model.optionVega( + 1.0, + adjustedStrike / adjustedForward, + ttm, + volatility, + option.isCallOption(), + discountFactor * adjustedForward); + + case Theta: + return Black76Model.optionTheta( + 1.0, + adjustedStrike / adjustedForward, + ttm, + volatility, + option.isCallOption(), + discountFactor * adjustedForward, + discountRate); + + default: + throw new NotImplementedException("Calculation for " + calcType + " not implemented yet."); + } + } + + public double getPrice( + Option option, + EquityForwardStructure forwardStructure, + YieldCurve discountCurve, + VolatilitySurface volaSurface) { + return calculate(option, forwardStructure, discountCurve, volaSurface, CalculationRequestType.Price); + } + + public double getDelta( + Option option, + EquityForwardStructure forwardStructure, + YieldCurve discountCurve, + VolatilitySurface volaSurface) { + return calculate(option, forwardStructure, discountCurve, volaSurface, CalculationRequestType.EqDelta); + } + + public double getGamma( + Option option, + EquityForwardStructure forwardStructure, + YieldCurve discountCurve, + VolatilitySurface volaSurface) { + return calculate(option, forwardStructure, discountCurve, volaSurface, CalculationRequestType.EqGamma); + } + + public double getVega( + Option option, + EquityForwardStructure forwardStructure, + YieldCurve discountCurve, + VolatilitySurface volaSurface) { + return calculate(option, forwardStructure, discountCurve, volaSurface, CalculationRequestType.EqVega); + } + + public double getTheta( + Option option, + EquityForwardStructure forwardStructure, + YieldCurve discountCurve, + VolatilitySurface volaSurface) { + return calculate(option, forwardStructure, discountCurve, volaSurface, CalculationRequestType.Theta); + } + + public double getImpliedVolatility( + Option option, + EquityForwardStructure forwardStructure, + YieldCurve discountCurve, + double price) { + + assert !option.isAmericanOption() : "Analytic pricer cannot handle American options."; + + final var expiryDate = option.getExpiryDate(); + final var ttm = dcc.getDaycountFraction(forwardStructure.getValuationDate(), expiryDate); + final var forward = forwardStructure.getForward(expiryDate); + final var discount = discountCurve.getDiscountFactor(expiryDate); + final var adjustedForward = forwardStructure.getDividendAdjustedStrike(forward, expiryDate); + final var adjustedStrike = forwardStructure.getDividendAdjustedStrike(option.getStrike(), expiryDate); + final var undiscountedPrice = price / discount / adjustedForward; + + return Black76Model.optionImpliedVolatility( + 1.0, + adjustedStrike / adjustedForward, + ttm, + undiscountedPrice, + option.isCallOption()); + } } diff --git a/src/main/java/net/finmath/equities/pricer/PdeOptionValuation.java b/src/main/java/net/finmath/equities/pricer/PdeOptionValuation.java index e15a2a1adb..93144bbc66 100644 --- a/src/main/java/net/finmath/equities/pricer/PdeOptionValuation.java +++ b/src/main/java/net/finmath/equities/pricer/PdeOptionValuation.java @@ -4,10 +4,12 @@ import java.util.Collections; import java.util.Comparator; import java.util.HashMap; + import org.apache.commons.math3.linear.DecompositionSolver; import org.apache.commons.math3.linear.LUDecomposition; import org.apache.commons.math3.linear.MatrixUtils; import org.apache.commons.math3.linear.RealMatrix; + import net.finmath.equities.marketdata.YieldCurve; import net.finmath.equities.models.EquityForwardStructure; import net.finmath.equities.models.FlatVolatilitySurface; @@ -36,337 +38,434 @@ */ public class PdeOptionValuation implements OptionValuation { - private final int timeStepsPerYear; - private final double spaceMinForwardMultiple; - private final double spaceMaxForwardMultiple; - private final int spaceNbOfSteps; - private final double spaceStepSize; - private final ArrayList spots; - private final int spotIndex; - private final DayCountConvention dayCounter; - private final boolean isLvPricer; - private final boolean includeDividendDatesInGrid; - - public PdeOptionValuation(double spaceMinForwardMultiple, double spaceMaxForwardMultiple, int spaceNbPoints, - final int timeStepsPerYear, DayCountConvention dcc, final boolean isLvPricer, - final boolean includeDividendDatesInGrid) { - assert spaceMinForwardMultiple < 1.0 : "min multiple of forward must be below 1.0"; - assert spaceMaxForwardMultiple > 1.0 : "max multiple of forward must be below 1.0"; - - this.timeStepsPerYear = timeStepsPerYear; - this.dayCounter = dcc; - this.isLvPricer = isLvPricer; - this.includeDividendDatesInGrid = includeDividendDatesInGrid; - - // Set up the space grid for the pure volatility process - var tmpSpaceStepSize = (spaceMaxForwardMultiple - spaceMinForwardMultiple) / spaceNbPoints; - var tmpSpaceNbPoints = spaceNbPoints; - var tmpSpots = new ArrayList(); - for (int i = 0; i < tmpSpaceNbPoints; i++) { - tmpSpots.add(spaceMinForwardMultiple + tmpSpaceStepSize * i); - } - // The space grid needs to include the forward level 1.0 for the pure volatility process - // Hence if necessary, we increase the step size slightly to include it - final var lowerBound = Math.abs(Collections.binarySearch(tmpSpots, 1.0)) - 2; - if (!(tmpSpots.get(lowerBound) == 1.0)) { - tmpSpaceStepSize += (1.0 - tmpSpots.get(lowerBound)) / lowerBound; - tmpSpots = new ArrayList(); - tmpSpaceNbPoints = 0; - var tmpSpot = 0.0; - while (tmpSpot < spaceMaxForwardMultiple) { - tmpSpot = spaceMinForwardMultiple + tmpSpaceStepSize * tmpSpaceNbPoints; - tmpSpots.add(tmpSpot); - tmpSpaceNbPoints++; - } - } - - this.spaceMinForwardMultiple = spaceMinForwardMultiple; - this.spaceMaxForwardMultiple = tmpSpots.get(tmpSpots.size() - 1); - this.spaceNbOfSteps = tmpSpaceNbPoints; - spots = tmpSpots; - spaceStepSize = tmpSpaceStepSize; - spotIndex = lowerBound; - } - - @Override - public EquityValuationResult calculate(EquityValuationRequest request, EquityForwardStructure forwardStructure, - YieldCurve discountCurve, VolatilitySurface volaSurface) { - final var results = new HashMap(); - if (request.getCalcsRequested().isEmpty()) { - return new EquityValuationResult(request, results); - } - - double price = 0.0; - if (request.getCalcsRequested().contains(CalculationRequestType.EqDelta) - || request.getCalcsRequested().contains(CalculationRequestType.EqGamma)) { - final var spotSensis = getPdeSensis(request.getOption(), forwardStructure, discountCurve, volaSurface); - price = spotSensis[0]; - if (request.getCalcsRequested().contains(CalculationRequestType.EqDelta)) { - results.put(CalculationRequestType.EqDelta, spotSensis[1]); - } - if (request.getCalcsRequested().contains(CalculationRequestType.EqGamma)) { - results.put(CalculationRequestType.EqGamma, spotSensis[2]); - } - } else { - price = getPrice(request.getOption(), forwardStructure, discountCurve, volaSurface); - } - - if (request.getCalcsRequested().contains(CalculationRequestType.Price)) { - results.put(CalculationRequestType.Price, price); - } - - if (request.getCalcsRequested().contains(CalculationRequestType.EqVega)) { - final var volShift = 0.0001; // TODO Make part of class members - final var priceShifted = getPrice(request.getOption(), forwardStructure, discountCurve, - volaSurface.getShiftedSurface(volShift)); - results.put(CalculationRequestType.EqVega, (priceShifted - price) / volShift); - } - - return new EquityValuationResult(request, results); - } - - public double getPrice(Option option, EquityForwardStructure forwardStructure, YieldCurve discountCurve, - VolatilitySurface volSurface) { - return evolvePde(option, forwardStructure, discountCurve, volSurface, false)[0]; - } - - public double[] getPdeSensis(Option option, EquityForwardStructure forwardStructure, YieldCurve discountCurve, - VolatilitySurface volSurface) { - return evolvePde(option, forwardStructure, discountCurve, volSurface, true); - } - - public double getVega(Option option, EquityForwardStructure forwardStructure, YieldCurve discountCurve, - VolatilitySurface volSurface, double basePrice, double volShift) { - final var shiftedPrice = getPrice(option, forwardStructure, discountCurve, - volSurface.getShiftedSurface(volShift)); - return (shiftedPrice - basePrice) / volShift; - } - - public double getTheta(Option option, EquityForwardStructure forwardStructure, YieldCurve discountCurve, - VolatilitySurface volSurface, double basePrice) { - final var valDate = forwardStructure.getValuationDate(); - final var thetaDate = valDate.plusDays(1); - final var thetaSpot = forwardStructure.getForward(thetaDate); - final var shiftedFwdStructure = forwardStructure.cloneWithNewSpot(thetaSpot).cloneWithNewDate(thetaDate); - final var shiftedPrice = getPrice(option, shiftedFwdStructure, discountCurve, volSurface); - return (shiftedPrice - basePrice) / dayCounter.getDaycountFraction(valDate, thetaDate); - } - - private double[] evolvePde(Option option, EquityForwardStructure forwardStructure, YieldCurve discountCurve, - VolatilitySurface volSurface, boolean calculateSensis) { - // Get data - final var valDate = forwardStructure.getValuationDate(); - final var expiryDate = option.getExpiryDate(); - final var expiryTime = dayCounter.getDaycountFraction(valDate, expiryDate); - assert !forwardStructure.getValuationDate() - .isAfter(expiryDate) : "Valuation date must not be after option expiry"; - final var impliedVol = volSurface.getVolatility(option.getStrike(), expiryDate, forwardStructure); - var forward = forwardStructure.getForward(expiryDate); - var fdf = forwardStructure.getFutureDividendFactor(expiryDate); - - // Build matrices - final RealMatrix idMatrix = MatrixUtils.createRealIdentityMatrix(spaceNbOfSteps); - final RealMatrix tridiagMatrix = MatrixUtils.createRealMatrix(spaceNbOfSteps, spaceNbOfSteps); - final double spaceStepSq = spaceStepSize * spaceStepSize; - for (int i = 0; i < spaceNbOfSteps; i++) { - for (int j = 0; j < spaceNbOfSteps; j++) { - if (i == j) { - tridiagMatrix.setEntry(i, j, Math.pow(spots.get(i), 2) / spaceStepSq); - } else if (i == j - 1 || i == j + 1) { - tridiagMatrix.setEntry(i, j, -0.5 * Math.pow(spots.get(i), 2) / spaceStepSq); - } else { - tridiagMatrix.setEntry(i, j, 0); - } - } - } - - // Set initial values - var prices = MatrixUtils.createRealVector(new double[spaceNbOfSteps]); - for (int i = 0; i < spaceNbOfSteps; i++) { - prices.setEntry(i, option.getPayoff((forward - fdf) * spots.get(i) + fdf)); - } - - // Set time intervals to evolve the PDE (i.e. from dividend to dividend) - final var diviDates = forwardStructure.getDividendStream().getDividendDates(); - final var anchorTimes = new ArrayList(); - anchorTimes.add(0.0); - if (includeDividendDatesInGrid) { - for (final var date : diviDates) { - if (date.isAfter(valDate) && date.isBefore(expiryDate)) { - anchorTimes.add(dayCounter.getDaycountFraction(valDate, date)); - } - } - } - anchorTimes.add(expiryTime); - anchorTimes.sort(Comparator.comparing(pt -> pt)); - var lastAtmPrice = 0.0; - var dt = 0.0; - - // Evolve PDE - for (int a = anchorTimes.size() - 1; a > 0; a--) { - // Set time steps - final var timeInterval = anchorTimes.get(a) - anchorTimes.get(a - 1); - int timeNbOfSteps; - double timeStepSize; - if (timeStepsPerYear == 0) // Use optimal ratio of time and space step size - { - timeNbOfSteps = (int) Math.ceil(2 * impliedVol * Math.pow(timeInterval, 1.5) / spaceStepSize); - timeStepSize = timeInterval / timeNbOfSteps; - } else // Use time step size provided externally - { - timeNbOfSteps = (int) Math.floor(timeInterval * timeStepsPerYear); - timeStepSize = timeInterval / timeNbOfSteps; - } - - final var times = new ArrayList(); - for (int i = 0; i <= 4; i++) { - times.add(anchorTimes.get(a) - i * 0.25 * timeStepSize); - } - for (int i = timeNbOfSteps - 2; i >= 0; i--) { - times.add(anchorTimes.get(a - 1) + i * timeStepSize); - } - - // Evolve PDE in current time interval - for (int i = 1; i < times.size(); i++) { - lastAtmPrice = prices.getEntry(spotIndex); - dt = times.get(i - 1) - times.get(i); - double theta = 0.5; - if (i <= 4) { - theta = 1.0; - } - final var theta1 = 1.0 - theta; - final var volSq = impliedVol * impliedVol; - - RealMatrix implicitMatrix, explicitMatrix; - if (isLvPricer) { - implicitMatrix = tridiagMatrix.scalarMultiply(theta * dt); - explicitMatrix = tridiagMatrix.scalarMultiply(-theta1 * dt); - final var localVol = new double[spaceNbOfSteps]; - for (int s = 0; s < spaceNbOfSteps; s++) { - final var lv = volSurface.getLocalVolatility(Math.log(spots.get(s)), times.get(i - 1), - forwardStructure, spaceStepSize, dt); - localVol[s] = lv * lv; - } - - final var volaMatrix = MatrixUtils.createRealDiagonalMatrix(localVol); - implicitMatrix = volaMatrix.multiply(implicitMatrix); - explicitMatrix = volaMatrix.multiply(explicitMatrix); - } else { - implicitMatrix = tridiagMatrix.scalarMultiply(theta * dt * volSq); - explicitMatrix = tridiagMatrix.scalarMultiply(-theta1 * dt * volSq); - } - - implicitMatrix = idMatrix.add(implicitMatrix); - explicitMatrix = idMatrix.add(explicitMatrix); - - if (option.isAmericanOption()) { - // Use the penalty algorithm from Forsyth's 2001 paper to solve the - // linear complementary problem for the American exercise feature. - final var penaltyFactor = 1 / Math.min(timeStepSize * timeStepSize, spaceStepSize * spaceStepSize); - forward = forwardStructure.getForward(times.get(i)); - fdf = forwardStructure.getFutureDividendFactor(times.get(i)); - final var discountFactor = discountCurve.getForwardDiscountFactor(times.get(i), expiryTime); - final var payoffs = MatrixUtils.createRealVector(new double[spaceNbOfSteps]); - final var penaltyMatrix = MatrixUtils.createRealMatrix(spaceNbOfSteps, spaceNbOfSteps); - for (int j = 1; j < spaceNbOfSteps - 1; j++) { - final var payoff = option.getPayoff((forward - fdf) * spots.get(j) + fdf) / discountFactor; - payoffs.setEntry(j, payoff); - penaltyMatrix.setEntry(j, j, prices.getEntry(j) < payoff ? penaltyFactor : 0); - } - - final var b = explicitMatrix.operate(prices); - var oldPrices = prices.copy(); - final var oldPenaltyMatrix = penaltyMatrix.copy(); - final var tol = 1 / penaltyFactor; - int iterations = 0; - while (true) { - assert iterations++ < 100 : "Penalty algorithm for american exercise did not converge in 100 steps"; - final var c = b.add(penaltyMatrix.operate(payoffs)); - final var A = implicitMatrix.add(penaltyMatrix); - final DecompositionSolver solver = new LUDecomposition(A).getSolver(); - prices = solver.solve(c); - for (int j = 1; j < spaceNbOfSteps - 1; j++) { - penaltyMatrix.setEntry(j, j, prices.getEntry(j) < payoffs.getEntry(j) ? penaltyFactor : 0); - } - - if (penaltyMatrix.equals(oldPenaltyMatrix) || (prices.subtract(oldPrices).getLInfNorm()) - / Math.max(oldPrices.getLInfNorm(), 1.0) < tol) { - break; - } - oldPrices = prices.copy(); - } - - } else { - // Solve the PDE step directly - prices = explicitMatrix.operate(prices); - final DecompositionSolver solver = new LUDecomposition(implicitMatrix).getSolver(); - prices = solver.solve(prices); - } - - // Set boundary conditions - prices.setEntry(0, option.getPayoff((forward - fdf) * spaceMinForwardMultiple + fdf)); - prices.setEntry(spaceNbOfSteps - 1, option.getPayoff((forward - fdf) * spaceMaxForwardMultiple + fdf)); - } - } - - final var discountFactor = discountCurve.getDiscountFactor(expiryDate); - final var price = discountFactor * prices.getEntry(spotIndex); - - if (calculateSensis) { - final var dFdX = forwardStructure.getDividendAdjustedStrike(forwardStructure.getForward(expiryDate), - expiryDate); - final var dFdS = forwardStructure.getGrowthDiscountFactor(valDate, expiryDate); - final var delta = discountFactor * 0.5 * (prices.getEntry(spotIndex + 1) - prices.getEntry(spotIndex - 1)) - / spaceStepSize * dFdS / dFdX; - final var gamma = discountFactor - * (prices.getEntry(spotIndex + 1) + prices.getEntry(spotIndex - 1) - 2 * prices.getEntry(spotIndex)) - / spaceStepSq * dFdS * dFdS / dFdX / dFdX; - final var discountFactorTheta = discountCurve.getDiscountFactor(expiryTime - dt); - final var theta = (discountFactorTheta * lastAtmPrice - price) / dt; - return new double[] { price, delta, gamma, theta }; - } else { - return new double[] { price, Double.NaN, Double.NaN, Double.NaN }; - } - } - - public double getImpliedVolatility(Option option, EquityForwardStructure forwardStructure, YieldCurve discountCurve, - double price) { - double initialGuess = 0.25; - final var forward = forwardStructure.getForward(option.getExpiryDate()); - // Use analytic pricer as initial guess for Europeans and OTM Americans - // Use two bisection steps for ITM Americans - if (option.isAmericanOption() && option.getPayoff(forward) > 0.0) { - final var bisectionSolver = new BisectionSearch(0.00001, 1.0); - for (int i = 0; i < 3; i++) { - final double currentVol = bisectionSolver.getNextPoint(); - final double currentPrice = getPrice(option, forwardStructure, discountCurve, - new FlatVolatilitySurface(currentVol)); - - bisectionSolver.setValue(currentPrice - price); - } - initialGuess = bisectionSolver.getBestPoint(); - } else { - final var anaPricer = new AnalyticOptionValuation(dayCounter); - Option testOption; - if (option.isAmericanOption()) { - testOption = new EuropeanOption(option.getExpiryDate(), option.getStrike(), option.isCallOption()); - } else { - testOption = option; - } - initialGuess = anaPricer.getImpliedVolatility(testOption, forwardStructure, discountCurve, price); - - } - - // Solve for implied vol - final var solver = new SecantMethod(initialGuess, initialGuess * 1.01); - while (solver.getAccuracy() / price > 1e-3 && !solver.isDone()) { - final double currentVol = solver.getNextPoint(); - final double currentPrice = getPrice(option, forwardStructure, discountCurve, - new FlatVolatilitySurface(currentVol)); - - solver.setValue(currentPrice - price); - } - - return Math.abs(solver.getBestPoint()); // Note that the PDE only uses sigma^2 - } + private final int timeStepsPerYear; + private final double spaceMinForwardMultiple; + private final double spaceMaxForwardMultiple; + private final int spaceNbOfSteps; + private final double spaceStepSize; + private final ArrayList spots; + private final int spotIndex; + private final DayCountConvention dayCounter; + private final boolean isLvPricer; + private final boolean includeDividendDatesInGrid; + + public PdeOptionValuation( + double spaceMinForwardMultiple, + double spaceMaxForwardMultiple, + int spaceNbPoints, + final int timeStepsPerYear, + DayCountConvention dcc, + final boolean isLvPricer, + final boolean includeDividendDatesInGrid) { + + assert spaceMinForwardMultiple < 1.0 : "min multiple of forward must be below 1.0"; + assert spaceMaxForwardMultiple > 1.0 : "max multiple of forward must be below 1.0"; + + this.timeStepsPerYear = timeStepsPerYear; + this.dayCounter = dcc; + this.isLvPricer = isLvPricer; + this.includeDividendDatesInGrid = includeDividendDatesInGrid; + + // Set up the space grid for the pure volatility process + var tmpSpaceStepSize = (spaceMaxForwardMultiple - spaceMinForwardMultiple) / spaceNbPoints; + var tmpSpaceNbPoints = spaceNbPoints; + var tmpSpots = new ArrayList(); + for(int i = 0; i < tmpSpaceNbPoints; i++) { + tmpSpots.add(spaceMinForwardMultiple + tmpSpaceStepSize * i); + } + + // The space grid needs to include the forward level 1.0 for the pure volatility process. + // Hence if necessary, we increase the step size slightly to include it. + final var lowerBound = Math.abs(Collections.binarySearch(tmpSpots, 1.0)) - 2; + if(!(tmpSpots.get(lowerBound) == 1.0)) { + tmpSpaceStepSize += (1.0 - tmpSpots.get(lowerBound)) / lowerBound; + tmpSpots = new ArrayList(); + tmpSpaceNbPoints = 0; + var tmpSpot = 0.0; + while(tmpSpot < spaceMaxForwardMultiple) { + tmpSpot = spaceMinForwardMultiple + tmpSpaceStepSize * tmpSpaceNbPoints; + tmpSpots.add(tmpSpot); + tmpSpaceNbPoints++; + } + } + + this.spaceMinForwardMultiple = spaceMinForwardMultiple; + this.spaceMaxForwardMultiple = tmpSpots.get(tmpSpots.size() - 1); + this.spaceNbOfSteps = tmpSpaceNbPoints; + spots = tmpSpots; + spaceStepSize = tmpSpaceStepSize; + spotIndex = lowerBound; + } + + @Override + public EquityValuationResult calculate( + EquityValuationRequest request, + EquityForwardStructure forwardStructure, + YieldCurve discountCurve, + VolatilitySurface volaSurface) { + + final var results = new HashMap(); + if(request.getCalcsRequested().isEmpty()) { + return new EquityValuationResult(request, results); + } + + double price = 0.0; + if(request.getCalcsRequested().contains(CalculationRequestType.EqDelta) + || request.getCalcsRequested().contains(CalculationRequestType.EqGamma)) { + + final var spotSensis = getPdeSensis(request.getOption(), forwardStructure, discountCurve, volaSurface); + price = spotSensis[0]; + + if(request.getCalcsRequested().contains(CalculationRequestType.EqDelta)) { + results.put(CalculationRequestType.EqDelta, spotSensis[1]); + } + if(request.getCalcsRequested().contains(CalculationRequestType.EqGamma)) { + results.put(CalculationRequestType.EqGamma, spotSensis[2]); + } + } + else { + price = getPrice(request.getOption(), forwardStructure, discountCurve, volaSurface); + } + + if(request.getCalcsRequested().contains(CalculationRequestType.Price)) { + results.put(CalculationRequestType.Price, price); + } + + if(request.getCalcsRequested().contains(CalculationRequestType.EqVega)) { + final var volShift = 0.0001; // TODO Make part of class members + final var priceShifted = getPrice( + request.getOption(), + forwardStructure, + discountCurve, + volaSurface.getShiftedSurface(volShift)); + results.put(CalculationRequestType.EqVega, (priceShifted - price) / volShift); + } + + return new EquityValuationResult(request, results); + } + + public double getPrice( + Option option, + EquityForwardStructure forwardStructure, + YieldCurve discountCurve, + VolatilitySurface volSurface) { + return evolvePde(option, forwardStructure, discountCurve, volSurface, false)[0]; + } + + public double[] getPdeSensis( + Option option, + EquityForwardStructure forwardStructure, + YieldCurve discountCurve, + VolatilitySurface volSurface) { + return evolvePde(option, forwardStructure, discountCurve, volSurface, true); + } + + public double getVega( + Option option, + EquityForwardStructure forwardStructure, + YieldCurve discountCurve, + VolatilitySurface volSurface, + double basePrice, + double volShift) { + + final var shiftedPrice = getPrice( + option, + forwardStructure, + discountCurve, + volSurface.getShiftedSurface(volShift)); + return (shiftedPrice - basePrice) / volShift; + } + + public double getTheta( + Option option, + EquityForwardStructure forwardStructure, + YieldCurve discountCurve, + VolatilitySurface volSurface, + double basePrice) { + + final var valDate = forwardStructure.getValuationDate(); + final var thetaDate = valDate.plusDays(1); + final var thetaSpot = forwardStructure.getForward(thetaDate); + final var shiftedFwdStructure = forwardStructure.cloneWithNewSpot(thetaSpot).cloneWithNewDate(thetaDate); + final var shiftedPrice = getPrice(option, shiftedFwdStructure, discountCurve, volSurface); + return (shiftedPrice - basePrice) / dayCounter.getDaycountFraction(valDate, thetaDate); + } + + private double[] evolvePde( + Option option, + EquityForwardStructure forwardStructure, + YieldCurve discountCurve, + VolatilitySurface volSurface, + boolean calculateSensis) { + + // Get data + final var valDate = forwardStructure.getValuationDate(); + final var expiryDate = option.getExpiryDate(); + final var expiryTime = dayCounter.getDaycountFraction(valDate, expiryDate); + + assert !forwardStructure.getValuationDate().isAfter(expiryDate) + : "Valuation date must not be after option expiry"; + + final var impliedVol = volSurface.getVolatility(option.getStrike(), expiryDate, forwardStructure); + var forward = forwardStructure.getForward(expiryDate); + var fdf = forwardStructure.getFutureDividendFactor(expiryDate); + + // Build matrices + final RealMatrix idMatrix = MatrixUtils.createRealIdentityMatrix(spaceNbOfSteps); + final RealMatrix tridiagMatrix = MatrixUtils.createRealMatrix(spaceNbOfSteps, spaceNbOfSteps); + final double spaceStepSq = spaceStepSize * spaceStepSize; + + for(int i = 0; i < spaceNbOfSteps; i++) { + for(int j = 0; j < spaceNbOfSteps; j++) { + if(i == j) { + tridiagMatrix.setEntry(i, j, Math.pow(spots.get(i), 2) / spaceStepSq); + } + else if(i == j - 1 || i == j + 1) { + tridiagMatrix.setEntry(i, j, -0.5 * Math.pow(spots.get(i), 2) / spaceStepSq); + } + else { + tridiagMatrix.setEntry(i, j, 0); + } + } + } + + // Set initial values + var prices = MatrixUtils.createRealVector(new double[spaceNbOfSteps]); + for(int i = 0; i < spaceNbOfSteps; i++) { + prices.setEntry(i, option.getPayoff((forward - fdf) * spots.get(i) + fdf)); + } + + // Set time intervals to evolve the PDE (i.e. from dividend to dividend) + final var diviDates = forwardStructure.getDividendStream().getDividendDates(); + final var anchorTimes = new ArrayList(); + anchorTimes.add(0.0); + + if(includeDividendDatesInGrid) { + for(final var date : diviDates) { + if(date.isAfter(valDate) && date.isBefore(expiryDate)) { + anchorTimes.add(dayCounter.getDaycountFraction(valDate, date)); + } + } + } + + anchorTimes.add(expiryTime); + anchorTimes.sort(Comparator.comparing(pt -> pt)); + + var lastAtmPrice = 0.0; + var dt = 0.0; + + // Evolve PDE + for(int a = anchorTimes.size() - 1; a > 0; a--) { + + // Set time steps + final var timeInterval = anchorTimes.get(a) - anchorTimes.get(a - 1); + int timeNbOfSteps; + double timeStepSize; + + if(timeStepsPerYear == 0) { + // Use optimal ratio of time and space step size + timeNbOfSteps = (int) Math.ceil(2 * impliedVol * Math.pow(timeInterval, 1.5) / spaceStepSize); + timeStepSize = timeInterval / timeNbOfSteps; + } + else { + // Use time step size provided externally + timeNbOfSteps = (int) Math.floor(timeInterval * timeStepsPerYear); + timeStepSize = timeInterval / timeNbOfSteps; + } + + final var times = new ArrayList(); + for(int i = 0; i <= 4; i++) { + times.add(anchorTimes.get(a) - i * 0.25 * timeStepSize); + } + for(int i = timeNbOfSteps - 2; i >= 0; i--) { + times.add(anchorTimes.get(a - 1) + i * timeStepSize); + } + + // Evolve PDE in current time interval + for(int i = 1; i < times.size(); i++) { + + lastAtmPrice = prices.getEntry(spotIndex); + dt = times.get(i - 1) - times.get(i); + + double theta = 0.5; + if(i <= 4) { + theta = 1.0; + } + final var theta1 = 1.0 - theta; + final var volSq = impliedVol * impliedVol; + + RealMatrix implicitMatrix; + RealMatrix explicitMatrix; + + if(isLvPricer) { + implicitMatrix = tridiagMatrix.scalarMultiply(theta * dt); + explicitMatrix = tridiagMatrix.scalarMultiply(-theta1 * dt); + + final var localVol = new double[spaceNbOfSteps]; + for(int s = 0; s < spaceNbOfSteps; s++) { + final var lv = volSurface.getLocalVolatility( + Math.log(spots.get(s)), + times.get(i - 1), + forwardStructure, + spaceStepSize, + dt); + localVol[s] = lv * lv; + } + + final var volaMatrix = MatrixUtils.createRealDiagonalMatrix(localVol); + implicitMatrix = volaMatrix.multiply(implicitMatrix); + explicitMatrix = volaMatrix.multiply(explicitMatrix); + } + else { + implicitMatrix = tridiagMatrix.scalarMultiply(theta * dt * volSq); + explicitMatrix = tridiagMatrix.scalarMultiply(-theta1 * dt * volSq); + } + + implicitMatrix = idMatrix.add(implicitMatrix); + explicitMatrix = idMatrix.add(explicitMatrix); + + if(option.isAmericanOption()) { + // Use the penalty algorithm from Forsyth's 2001 paper to solve the + // linear complementary problem for the American exercise feature. + final var penaltyFactor = 1 / Math.min(timeStepSize * timeStepSize, spaceStepSize * spaceStepSize); + + forward = forwardStructure.getForward(times.get(i)); + fdf = forwardStructure.getFutureDividendFactor(times.get(i)); + final var discountFactor = discountCurve.getForwardDiscountFactor(times.get(i), expiryTime); + + final var payoffs = MatrixUtils.createRealVector(new double[spaceNbOfSteps]); + final var penaltyMatrix = MatrixUtils.createRealMatrix(spaceNbOfSteps, spaceNbOfSteps); + + for(int j = 1; j < spaceNbOfSteps - 1; j++) { + final var payoff = option.getPayoff((forward - fdf) * spots.get(j) + fdf) / discountFactor; + payoffs.setEntry(j, payoff); + penaltyMatrix.setEntry(j, j, prices.getEntry(j) < payoff ? penaltyFactor : 0); + } + + final var b = explicitMatrix.operate(prices); + var oldPrices = prices.copy(); + final var oldPenaltyMatrix = penaltyMatrix.copy(); + final var tol = 1 / penaltyFactor; + + int iterations = 0; + while(true) { + + assert iterations++ < 100 : "Penalty algorithm for american exercise did not converge in 100 steps"; + + final var c = b.add(penaltyMatrix.operate(payoffs)); + final var A = implicitMatrix.add(penaltyMatrix); + final DecompositionSolver solver = new LUDecomposition(A).getSolver(); + + prices = solver.solve(c); + + for(int j = 1; j < spaceNbOfSteps - 1; j++) { + penaltyMatrix.setEntry(j, j, prices.getEntry(j) < payoffs.getEntry(j) ? penaltyFactor : 0); + } + + if(penaltyMatrix.equals(oldPenaltyMatrix) + || (prices.subtract(oldPrices).getLInfNorm()) + / Math.max(oldPrices.getLInfNorm(), 1.0) < tol) { + break; + } + + oldPrices = prices.copy(); + } + } + else { + // Solve the PDE step directly + prices = explicitMatrix.operate(prices); + final DecompositionSolver solver = new LUDecomposition(implicitMatrix).getSolver(); + prices = solver.solve(prices); + } + + // Set boundary conditions + prices.setEntry(0, option.getPayoff((forward - fdf) * spaceMinForwardMultiple + fdf)); + prices.setEntry(spaceNbOfSteps - 1, option.getPayoff((forward - fdf) * spaceMaxForwardMultiple + fdf)); + } + } + + final var discountFactor = discountCurve.getDiscountFactor(expiryDate); + final var price = discountFactor * prices.getEntry(spotIndex); + + if(calculateSensis) { + + final var dFdX = forwardStructure.getDividendAdjustedStrike( + forwardStructure.getForward(expiryDate), + expiryDate); + final var dFdS = forwardStructure.getGrowthDiscountFactor(valDate, expiryDate); + + final var delta = discountFactor + * 0.5 * (prices.getEntry(spotIndex + 1) - prices.getEntry(spotIndex - 1)) + / spaceStepSize * dFdS / dFdX; + + final var gamma = discountFactor + * (prices.getEntry(spotIndex + 1) + prices.getEntry(spotIndex - 1) - 2 * prices.getEntry(spotIndex)) + / spaceStepSq * dFdS * dFdS / dFdX / dFdX; + + final var discountFactorTheta = discountCurve.getDiscountFactor(expiryTime - dt); + final var theta = (discountFactorTheta * lastAtmPrice - price) / dt; + + return new double[] { price, delta, gamma, theta }; + } + else { + return new double[] { price, Double.NaN, Double.NaN, Double.NaN }; + } + } + + public double getImpliedVolatility( + Option option, + EquityForwardStructure forwardStructure, + YieldCurve discountCurve, + double price) { + + double initialGuess = 0.25; + final var forward = forwardStructure.getForward(option.getExpiryDate()); + + // Use analytic pricer as initial guess for Europeans and OTM Americans + // Use two bisection steps for ITM Americans + if(option.isAmericanOption() && option.getPayoff(forward) > 0.0) { + + final var bisectionSolver = new BisectionSearch(0.00001, 1.0); + for(int i = 0; i < 3; i++) { + final double currentVol = bisectionSolver.getNextPoint(); + final double currentPrice = getPrice( + option, + forwardStructure, + discountCurve, + new FlatVolatilitySurface(currentVol)); + + bisectionSolver.setValue(currentPrice - price); + } + initialGuess = bisectionSolver.getBestPoint(); + } + else { + final var anaPricer = new AnalyticOptionValuation(dayCounter); + Option testOption; + + if(option.isAmericanOption()) { + testOption = new EuropeanOption(option.getExpiryDate(), option.getStrike(), option.isCallOption()); + } + else { + testOption = option; + } + + initialGuess = anaPricer.getImpliedVolatility(testOption, forwardStructure, discountCurve, price); + } + + // Solve for implied vol + final var solver = new SecantMethod(initialGuess, initialGuess * 1.01); + while(solver.getAccuracy() / price > 1e-3 && !solver.isDone()) { + + final double currentVol = solver.getNextPoint(); + final double currentPrice = getPrice( + option, + forwardStructure, + discountCurve, + new FlatVolatilitySurface(currentVol)); + + solver.setValue(currentPrice - price); + } + + return Math.abs(solver.getBestPoint()); // Note that the PDE only uses sigma^2 + } } diff --git a/src/main/java/net/finmath/functions/HestonModel.java b/src/main/java/net/finmath/functions/HestonModel.java index 02f74270d9..2c8e150821 100644 --- a/src/main/java/net/finmath/functions/HestonModel.java +++ b/src/main/java/net/finmath/functions/HestonModel.java @@ -13,361 +13,347 @@ /** * This class implements some functions as static class methods related to the Heston model. - * The calculation is performed by means of the FFT algorithm of Carr Madan applied to the gradient of the characteristic funtion. + * The calculation is performed by means of the FFT algorithm of Carr Madan applied to the gradient of the + * characteristic funtion. * * @author Alessandro Gnoatto */ public class HestonModel { /* - * Parameters for the FFT calculator + * Parameters for the FFT calculator. */ - final static InterpolationMethod intMethod =InterpolationMethod.LINEAR; - final static ExtrapolationMethod extMethod = ExtrapolationMethod.CONSTANT; - final static int numberOfPoints = 4096*2; - final static double gridSpacing = 0.4; + private static final InterpolationMethod INT_METHOD = InterpolationMethod.LINEAR; + private static final ExtrapolationMethod EXT_METHOD = ExtrapolationMethod.CONSTANT; + private static final int NUMBER_OF_POINTS = 4096 * 2; + private static final double GRID_SPACING = 0.4; - enum HestonGreek {DELTA, GAMMA, THETA, RHO, VEGA1, VANNA, VOLGA}; + enum HestonGreek {DELTA, GAMMA, THETA, RHO, VEGA1, VANNA, VOLGA} /** * Calculates the delta of a call option under a Heston model. - * + * * @param initialStockValue Initital value of the stock. * @param riskFreeRate The risk free rate. * @param dividendYield The dividend yield. - * @param kappa the speed of mean reversion. - * @param theta the long run mean of the volatility. - * @param sigma the volatility of variance. - * @param v0 the initial instantaneous variance - * @param rho correlation between the two Brownian motions - * @param optionMaturity the maturity of the option - * @param optionStrike the strike of the option. + * @param kappa The speed of mean reversion. + * @param theta The long run mean of the volatility. + * @param sigma The volatility of variance. + * @param v0 The initial instantaneous variance. + * @param rho Correlation between the two Brownian motions. + * @param optionMaturity The maturity of the option. + * @param optionStrike The strike of the option. * @return The delta of the option. */ public static double hestonOptionDelta( final double initialStockValue, final double riskFreeRate, final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double kappa, + final double theta, + final double sigma, + final double v0, final double rho, final double optionMaturity, - final double optionStrike) - { + final double optionStrike) { - HestonGreek whatToCompute = HestonGreek.DELTA; + final HestonGreek whatToCompute = HestonGreek.DELTA; - return hestonGreekCalculator(initialStockValue, + return hestonGreekCalculator( + initialStockValue, riskFreeRate, dividendYield, - kappa, - theta, - sigma, - v0, + kappa, + theta, + sigma, + v0, rho, optionMaturity, optionStrike, whatToCompute, - numberOfPoints, - gridSpacing); + NUMBER_OF_POINTS, + GRID_SPACING); } /** - * Calculates the gamma of a call option under a Heston model - * + * Calculates the gamma of a call option under a Heston model. + * * @param initialStockValue Initital value of the stock. * @param riskFreeRate The risk free rate. * @param dividendYield The dividend yield. - * @param kappa the speed of mean reversion. - * @param theta the long run mean of the volatility. - * @param sigma the volatility of variance. - * @param v0 the initial instantaneous variance - * @param rho correlation between the two Brownian motions - * @param optionMaturity the maturity of the option - * @param optionStrike the strike of the option. - * @return The gamma of the option + * @param kappa The speed of mean reversion. + * @param theta The long run mean of the volatility. + * @param sigma The volatility of variance. + * @param v0 The initial instantaneous variance. + * @param rho Correlation between the two Brownian motions. + * @param optionMaturity The maturity of the option. + * @param optionStrike The strike of the option. + * @return The gamma of the option. */ public static double hestonOptionGamma( final double initialStockValue, final double riskFreeRate, final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double kappa, + final double theta, + final double sigma, + final double v0, final double rho, final double optionMaturity, - final double optionStrike) - { - HestonGreek whatToCompute = HestonGreek.GAMMA; + final double optionStrike) { - return hestonGreekCalculator(initialStockValue, + final HestonGreek whatToCompute = HestonGreek.GAMMA; + + return hestonGreekCalculator( + initialStockValue, riskFreeRate, dividendYield, - kappa, - theta, - sigma, - v0, + kappa, + theta, + sigma, + v0, rho, optionMaturity, optionStrike, whatToCompute, - numberOfPoints, - gridSpacing); + NUMBER_OF_POINTS, + GRID_SPACING); } /** - * Calculates the theta of a call option under a Heston model - * + * Calculates the theta of a call option under a Heston model. + * * @param initialStockValue Initital value of the stock. * @param riskFreeRate The risk free rate. * @param dividendYield The dividend yield. - * @param kappa the speed of mean reversion. - * @param theta the long run mean of the volatility. - * @param sigma the volatility of variance. - * @param v0 the initial instantaneous variance - * @param rho correlation between the two Brownian motions - * @param optionMaturity the maturity of the option - * @param optionStrike the strike of the option. - * @return The theta of the option + * @param kappa The speed of mean reversion. + * @param theta The long run mean of the volatility. + * @param sigma The volatility of variance. + * @param v0 The initial instantaneous variance. + * @param rho Correlation between the two Brownian motions. + * @param optionMaturity The maturity of the option. + * @param optionStrike The strike of the option. + * @return The theta of the option. */ public static double hestonOptionTheta( final double initialStockValue, final double riskFreeRate, final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double kappa, + final double theta, + final double sigma, + final double v0, final double rho, final double optionMaturity, - final double optionStrike) - { - HestonGreek whatToCompute = HestonGreek.THETA; + final double optionStrike) { + + final HestonGreek whatToCompute = HestonGreek.THETA; - return hestonGreekCalculator(initialStockValue, + return hestonGreekCalculator( + initialStockValue, riskFreeRate, dividendYield, - kappa, - theta, - sigma, - v0, + kappa, + theta, + sigma, + v0, rho, optionMaturity, optionStrike, whatToCompute, - numberOfPoints, - gridSpacing); + NUMBER_OF_POINTS, + GRID_SPACING); } /** - * Calculates the rho of a call option under a Heston model - * + * Calculates the rho of a call option under a Heston model. + * * @param initialStockValue Initital value of the stock. * @param riskFreeRate The risk free rate. * @param dividendYield The dividend yield. - * @param kappa the speed of mean reversion. - * @param theta the long run mean of the volatility. - * @param sigma the volatility of variance. - * @param v0 the initial instantaneous variance - * @param rho correlation between the two Brownian motions - * @param optionMaturity the maturity of the option - * @param optionStrike the strike of the option. - * @return The rho of the option + * @param kappa The speed of mean reversion. + * @param theta The long run mean of the volatility. + * @param sigma The volatility of variance. + * @param v0 The initial instantaneous variance. + * @param rho Correlation between the two Brownian motions. + * @param optionMaturity The maturity of the option. + * @param optionStrike The strike of the option. + * @return The rho of the option. */ public static double hestonOptionRho( final double initialStockValue, final double riskFreeRate, final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double kappa, + final double theta, + final double sigma, + final double v0, final double rho, final double optionMaturity, - final double optionStrike) - { - HestonGreek whatToCompute = HestonGreek.RHO; + final double optionStrike) { + + final HestonGreek whatToCompute = HestonGreek.RHO; - return hestonGreekCalculator(initialStockValue, + return hestonGreekCalculator( + initialStockValue, riskFreeRate, dividendYield, - kappa, - theta, - sigma, - v0, + kappa, + theta, + sigma, + v0, rho, optionMaturity, optionStrike, whatToCompute, - numberOfPoints, - gridSpacing); + NUMBER_OF_POINTS, + GRID_SPACING); } /** - * Calculates the vega1 of a call option under a Heston model - * + * Calculates the vega1 of a call option under a Heston model. + * * @param initialStockValue Initital value of the stock. * @param riskFreeRate The risk free rate. * @param dividendYield The dividend yield. - * @param kappa the speed of mean reversion. - * @param theta the long run mean of the volatility. - * @param sigma the volatility of variance. - * @param v0 the initial instantaneous variance - * @param rho correlation between the two Brownian motions - * @param optionMaturity the maturity of the option - * @param optionStrike the strike of the option. - * @return The vega1 of the option + * @param kappa The speed of mean reversion. + * @param theta The long run mean of the volatility. + * @param sigma The volatility of variance. + * @param v0 The initial instantaneous variance. + * @param rho Correlation between the two Brownian motions. + * @param optionMaturity The maturity of the option. + * @param optionStrike The strike of the option. + * @return The vega1 of the option. */ public static double hestonOptionVega1( final double initialStockValue, final double riskFreeRate, final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double kappa, + final double theta, + final double sigma, + final double v0, final double rho, final double optionMaturity, - final double optionStrike) - { - HestonGreek whatToCompute = HestonGreek.VEGA1; + final double optionStrike) { - return hestonGreekCalculator(initialStockValue, + final HestonGreek whatToCompute = HestonGreek.VEGA1; + + return hestonGreekCalculator( + initialStockValue, riskFreeRate, dividendYield, - kappa, - theta, - sigma, - v0, + kappa, + theta, + sigma, + v0, rho, optionMaturity, optionStrike, whatToCompute, - numberOfPoints, - gridSpacing); + NUMBER_OF_POINTS, + GRID_SPACING); } /** - * Calculates the vanna of a call option under a Heston model - * + * Calculates the vanna of a call option under a Heston model. + * * @param initialStockValue Initital value of the stock. * @param riskFreeRate The risk free rate. * @param dividendYield The dividend yield. - * @param kappa the speed of mean reversion. - * @param theta the long run mean of the volatility. - * @param sigma the volatility of variance. - * @param v0 the initial instantaneous variance - * @param rho correlation between the two Brownian motions - * @param optionMaturity the maturity of the option - * @param optionStrike the strike of the option. - * @return The vanna of the option + * @param kappa The speed of mean reversion. + * @param theta The long run mean of the volatility. + * @param sigma The volatility of variance. + * @param v0 The initial instantaneous variance. + * @param rho Correlation between the two Brownian motions. + * @param optionMaturity The maturity of the option. + * @param optionStrike The strike of the option. + * @return The vanna of the option. */ public static double hestonOptionVanna( final double initialStockValue, final double riskFreeRate, final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double kappa, + final double theta, + final double sigma, + final double v0, final double rho, final double optionMaturity, - final double optionStrike) - { - HestonGreek whatToCompute = HestonGreek.VANNA; + final double optionStrike) { + + final HestonGreek whatToCompute = HestonGreek.VANNA; - return hestonGreekCalculator(initialStockValue, + return hestonGreekCalculator( + initialStockValue, riskFreeRate, dividendYield, - kappa, - theta, - sigma, - v0, + kappa, + theta, + sigma, + v0, rho, optionMaturity, optionStrike, whatToCompute, - numberOfPoints, - gridSpacing); + NUMBER_OF_POINTS, + GRID_SPACING); } /** - * Calculates the volga of a call option under a Heston model - * + * Calculates the volga of a call option under a Heston model. + * * @param initialStockValue Initital value of the stock. * @param riskFreeRate The risk free rate. * @param dividendYield The dividend yield. - * @param kappa the speed of mean reversion. - * @param theta the long run mean of the volatility. - * @param sigma the volatility of variance. - * @param v0 the initial instantaneous variance - * @param rho correlation between the two Brownian motions - * @param optionMaturity the maturity of the option - * @param optionStrike the strike of the option. - * @return The volga of the option + * @param kappa The speed of mean reversion. + * @param theta The long run mean of the volatility. + * @param sigma The volatility of variance. + * @param v0 The initial instantaneous variance. + * @param rho Correlation between the two Brownian motions. + * @param optionMaturity The maturity of the option. + * @param optionStrike The strike of the option. + * @return The volga of the option. */ public static double hestonOptionVolga( final double initialStockValue, final double riskFreeRate, final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double kappa, + final double theta, + final double sigma, + final double v0, final double rho, final double optionMaturity, - final double optionStrike) - { - HestonGreek whatToCompute = HestonGreek.VOLGA; + final double optionStrike) { + + final HestonGreek whatToCompute = HestonGreek.VOLGA; - return hestonGreekCalculator(initialStockValue, + return hestonGreekCalculator( + initialStockValue, riskFreeRate, dividendYield, - kappa, - theta, - sigma, - v0, + kappa, + theta, + sigma, + v0, rho, optionMaturity, optionStrike, whatToCompute, - numberOfPoints, - gridSpacing); + NUMBER_OF_POINTS, + GRID_SPACING); } - - /** - * Computes the gradient of the (discounted) characteristic function of the Heston model. - * More sensitivies can be added by introducing more cases in the switch statement. - * - * @param zeta - * @param initialStockValue - * @param riskFreeRate - * @param dividendYield - * @param kappa - * @param theta - * @param sigma - * @param v0 - * @param rho - * @param optionMaturity - * @param optionStrike - * @param whichGreek - * @param numberOfPoints - * @param gridSpacing - * @return - */ private static Complex hestonCharacteristicFunctionGradient( final Complex zeta, final double initialStockValue, final double riskFreeRate, final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double kappa, + final double theta, + final double sigma, + final double v0, final double rho, final double optionMaturity, final double optionStrike, @@ -375,151 +361,115 @@ private static Complex hestonCharacteristicFunctionGradient( final int numberOfPoints, final double gridSpacing) { - final double lambda = 2*Math.PI/(numberOfPoints*gridSpacing); + final double lambda = 2 * Math.PI / (numberOfPoints * gridSpacing); - double x = Math.log(initialStockValue); - double a = kappa * theta; + final double x = Math.log(initialStockValue); + final double a = kappa * theta; - Complex b, u, d, g, c, D, G, C; + Complex b; + Complex u; + Complex d; + Complex g; + Complex c; + Complex D; + Complex G; + Complex C; - // Initialize u = new Complex(-0.5, 0.0); b = new Complex(kappa + lambda, 0.0); - // Compute d - Complex term1 = (Complex.I.multiply(rho * sigma).multiply(zeta)).subtract(b); // (ρσiφ - b) - Complex powPart = term1.pow(2.0); // (...) - Complex term2 = new Complex(sigma * sigma, 0.0) + final Complex term1 = (Complex.I.multiply(rho * sigma).multiply(zeta)).subtract(b); + final Complex powPart = term1.pow(2.0); + final Complex term2 = new Complex(sigma * sigma, 0.0) .multiply((u.multiply(2.0).multiply(Complex.I).multiply(zeta)) - .subtract(zeta.multiply(zeta))); // σ²(2u i φ - φ²) + .subtract(zeta.multiply(zeta))); d = powPart.subtract(term2).sqrt(); - // Compute g - Complex numerator = (b.subtract(Complex.I.multiply(rho * sigma).multiply(zeta))).add(d); - Complex denominator = (b.subtract(Complex.I.multiply(rho * sigma).multiply(zeta))).subtract(d); + final Complex numerator = (b.subtract(Complex.I.multiply(rho * sigma).multiply(zeta))).add(d); + final Complex denominator = (b.subtract(Complex.I.multiply(rho * sigma).multiply(zeta))).subtract(d); g = numerator.divide(denominator); - // Initialize remaining - c = Complex.ZERO; - D = Complex.ZERO; - G = Complex.ZERO; - C = Complex.ZERO; - - // "Little Heston Trap" formulation - // c = 1.0 / g c = Complex.ONE.divide(g); - // Precompute some recurring parts - Complex exp_dT = d.multiply(-optionMaturity).exp(); // exp(-d*T) - Complex bMinusRhoSigmaIphiMinusD = b.subtract(Complex.I.multiply(rho * sigma).multiply(zeta)).subtract(d); + final Complex exp_dT = d.multiply(-optionMaturity).exp(); + final Complex bMinusRhoSigmaIphiMinusD = b.subtract(Complex.I.multiply(rho * sigma).multiply(zeta)).subtract(d); - // D = (b - rho*sigma*i*phi - d) / (sigma^2) * ((1 - exp(-d*T)) / (1 - c * exp(-d*T))) D = bMinusRhoSigmaIphiMinusD .divide(sigma * sigma) .multiply( (Complex.ONE.subtract(exp_dT)) - .divide(Complex.ONE.subtract(c.multiply(exp_dT))) - ); + .divide(Complex.ONE.subtract(c.multiply(exp_dT)))); - // G = (1 - c * exp(-d*T)) / (1 - c) G = (Complex.ONE.subtract(c.multiply(exp_dT))) .divide(Complex.ONE.subtract(c)); - // C = (r - q) * i * phi * T + a / sigma^2 * ((b - rho*sigma*i*phi - d) * T - 2.0 * Complex.Log(G)) - Complex firstTerm = Complex.I.multiply(zeta).multiply((riskFreeRate - dividendYield) * optionMaturity); - Complex secondTerm = new Complex(a / (sigma * sigma), 0.0) + final Complex firstTerm = Complex.I.multiply(zeta).multiply((riskFreeRate - dividendYield) * optionMaturity); + final Complex secondTerm = new Complex(a / (sigma * sigma), 0.0) .multiply( (bMinusRhoSigmaIphiMinusD.multiply(optionMaturity)) - .subtract(Complex.valueOf(2.0).multiply(G.log())) - ); + .subtract(Complex.valueOf(2.0).multiply(G.log()))); C = firstTerm.add(secondTerm); - // The characteristic function - Complex f = (C.add(D.multiply(v0)).add(Complex.I.multiply(zeta).multiply(x))).exp(); - + final Complex f = (C.add(D.multiply(v0)).add(Complex.I.multiply(zeta).multiply(x))).exp(); - // Return depending on the requested Greek - switch (whichGreek) { + switch(whichGreek) { case DELTA: - return f.multiply(Complex.I).multiply(zeta).divide(initialStockValue); + return f.multiply(Complex.I).multiply(zeta).divide(initialStockValue); case GAMMA: return Complex.I.multiply(zeta) .multiply(f) .multiply(Complex.I.multiply(zeta).subtract(Complex.ONE)) - .divide(initialStockValue * initialStockValue); + .divide(initialStockValue * initialStockValue); case THETA: - Complex numerator_dDdT = d.multiply(exp_dT) - .multiply(b.subtract(Complex.I.multiply(rho * sigma).multiply(zeta)).add(d)) - .multiply(g.subtract(Complex.ONE)); - Complex denominator_dDdT = Complex.valueOf(sigma * sigma) + final Complex numerator_dDdT = d.multiply(exp_dT) + .multiply(b.subtract(Complex.I.multiply(rho * sigma).multiply(zeta)).add(d)) + .multiply(g.subtract(Complex.ONE)); + final Complex denominator_dDdT = Complex.valueOf(sigma * sigma) .multiply(Complex.ONE.subtract(g.multiply(exp_dT)).pow(2.0)); - Complex dDdT = numerator_dDdT.divide(denominator_dDdT); + final Complex dDdT = numerator_dDdT.divide(denominator_dDdT); - Complex innerTerm = (b.subtract(Complex.I.multiply(rho * sigma).multiply(zeta)).add(d)) + final Complex innerTerm = (b.subtract(Complex.I.multiply(rho * sigma).multiply(zeta)).add(d)) .add( Complex.valueOf(2.0) - .multiply(g) - .multiply(d) - .multiply(exp_dT) - .divide(Complex.ONE.subtract(g.multiply(exp_dT))) - ); + .multiply(g) + .multiply(d) + .multiply(exp_dT) + .divide(Complex.ONE.subtract(g.multiply(exp_dT)))); - Complex dCdT = Complex.I.multiply(zeta).multiply(riskFreeRate) + final Complex dCdT = Complex.I.multiply(zeta).multiply(riskFreeRate) .add( Complex.valueOf(kappa * theta / (sigma * sigma)) - .multiply(innerTerm) - ); + .multiply(innerTerm)); return Complex.valueOf(riskFreeRate) .multiply(f) - .subtract( - f.multiply(dCdT.add(dDdT.multiply(v0))) - ); + .subtract(f.multiply(dCdT.add(dDdT.multiply(v0)))); case RHO: - Complex dCdr = Complex.I.multiply(zeta).multiply(optionMaturity); - return f.multiply(dCdr).subtract(f.multiply(optionMaturity)); + final Complex dCdr = Complex.I.multiply(zeta).multiply(optionMaturity); + return f.multiply(dCdr).subtract(f.multiply(optionMaturity)); case VEGA1: - return f.multiply(D); + return f.multiply(D); case VANNA: - return f.multiply(Complex.I).multiply(zeta).multiply(D).divide(initialStockValue); + return f.multiply(Complex.I).multiply(zeta).multiply(D).divide(initialStockValue); case VOLGA: return Complex.valueOf(2.0) .multiply(D) .multiply(f) - .multiply( - D.multiply(2.0 * v0).add(Complex.ONE) - ); + .multiply(D.multiply(2.0 * v0).add(Complex.ONE)); default: throw new IllegalArgumentException("Unknown Greek: " + whichGreek); } } - /** - * Service method that performs the Fourier inversion according to the FFT algorithm of Carr and Madan - * - * @param initialStockValue Initital value of the stock. - * @param riskFreeRate The risk free rate. - * @param dividendYield The dividend yield. - * @param kappa the speed of mean reversion. - * @param theta the long run mean of the volatility. - * @param sigma the volatility of variance. - * @param v0 the initial instantaneous variance - * @param rho correlation between the two Brownian motions - * @param optionMaturity the maturity of the option - * @param optionStrike the strike of the option. - * @param whichGreek - * @param numberOfPoints - * @param gridSpacing - * @return the requested calculation - */ - private static double hestonGreekCalculator(final double initialStockValue, + private static double hestonGreekCalculator( + final double initialStockValue, final double riskFreeRate, final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double kappa, + final double theta, + final double sigma, + final double v0, final double rho, final double optionMaturity, final double optionStrike, @@ -529,21 +479,25 @@ private static double hestonGreekCalculator(final double initialStockValue, final double lineOfIntegration = 1.2; - final double lambda = 2*Math.PI/(numberOfPoints*gridSpacing); //Equation 23 Carr and Madan - final double upperBound = (numberOfPoints * lambda)/2.0; //Equation 20 Carr and Madan + final double lambda = 2 * Math.PI / (numberOfPoints * gridSpacing); + final double upperBound = (numberOfPoints * lambda) / 2.0; final Complex[] integrandEvaluations = new Complex[numberOfPoints]; - for(int i = 0; i strikeToValue = new Function(){ + final RationalFunctionInterpolation interpolation = new RationalFunctionInterpolation( + strikeVector, + valuesVector, + INT_METHOD, + EXT_METHOD); + final Function strikeToValue = new Function() { @Override public Double apply(final Double t) { - return interpolation.getValue(t); + return interpolation.getValue(t); } - }; return strikeToValue.apply(optionStrike); diff --git a/src/main/java/net/finmath/tree/AbstractTreeProduct.java b/src/main/java/net/finmath/tree/AbstractTreeProduct.java index f459bfbbda..49ab67d558 100644 --- a/src/main/java/net/finmath/tree/AbstractTreeProduct.java +++ b/src/main/java/net/finmath/tree/AbstractTreeProduct.java @@ -5,15 +5,16 @@ /** * Base class for products that can be priced on a (recombining) with TreeModel. * This abstraction provides a small template: - * getValue(TreeModel) returns the scalar time–0 price. - * getValue(double, TreeModel) returns the value as a RandomVariable at a given evaluation time. - * getValues(double, TreeModel) must be implemented by subclasses and - * should return the full vector of values (per state) at each time level, - * produced by a backward induction. By convention the element at index - * 0 corresponds to time 0. - * Subclasses (e.g. European, American, Asian ) implement their - * payoff logic and early–exercise features within getValues(double, TreeModel). - * + *

    + *
  • getValue(TreeModel) returns the scalar time–0 price.
    • + *
  • getValue(double, TreeModel) returns the value as a RandomVariable at a given evaluation time.
    • + *
  • getValues(double, TreeModel) must be implemented by subclasses and should return the full vector of values + * (per state) at each time level, produced by a backward induction. By convention the element at index 0 + * corresponds to time 0.
    • + *
+ * Subclasses (e.g. European, American, Asian) implement their payoff logic and early–exercise features within + * getValues(double, TreeModel). + * * @author Carlo Andrea Tramentozzi * @author Alessandro Gnoatto * @author Andrea Mazzon @@ -26,58 +27,58 @@ public abstract class AbstractTreeProduct implements TreeProduct { /** * Creates a tree-priced product with the given maturity. * - * @param maturity The contract maturity - * must be compatible with the model time grid. + * @param maturity The contract maturity must be compatible with the model time grid. */ - public AbstractTreeProduct(double maturity){ + public AbstractTreeProduct(final double maturity) { this.maturity = maturity; } /** - * Returns the time–0 present value under the given model - * This is a convenience method delegating to getValue(double, TreeModel)}ù - * with evaluationTime = 0.0 and extracting the scalar value from - * the returned RandomVariable. + * Returns the time–0 present value under the given model. + * This is a convenience method delegating to getValue(double, TreeModel) + * with evaluationTime = 0.0 and extracting the scalar value from the returned RandomVariable. * * @param model The tree model to use (providing discounting and conditional expectations). * @return The time–0 price. * @throws IllegalArgumentException If model is null or the evaluation fails. */ @Override - public double getValue(TreeModel model ){ - RandomVariable v0 = getValue(0.0,model); + public double getValue(final TreeModel model) { + final RandomVariable v0 = getValue(0.0, model); return v0.get(0); } + /** * Returns the product value at a given evaluation time as a RandomVariable. - * The default implementation calls getValues(double, TreeModel) and - * returns the first component (time–level 0). Subclasses should ensure - * their {@link #getValues(double, TreeModel)} respects this convention. - * Internally, the method creates all values from time zero and returns the level - * of interest. While this is inefficient for non-path-dependent products, it ensures - * a correct implementation for path dependent products. The alternative would be to provide as - * input the current value of the path-dependent functional, which would break the interface. + * The default implementation calls getValues(double, TreeModel) and returns the corresponding component. + * Subclasses should ensure their {@link #getValues(double, TreeModel)} respects this convention. + * Internally, the method creates all values from time zero and returns the level of interest. While this is + * inefficient for non-path-dependent products, it ensures a correct implementation for path dependent products. + * The alternative would be to provide as input the current value of the path-dependent functional, which would + * break the interface. * * @param evaluationTime The time at which the value is requested (must be >= 0). * @param model The tree model to use. - * @return The value at evalutationTime as a random variable on the tree. + * @return The value at evaluationTime as a random variable on the tree. * @throws IllegalArgumentException If the inputs are invalid (see validate(double, TreeModel)). */ - public RandomVariable getValue(double evaluationTime, TreeModel model) { - RandomVariable[] levels = getValues(0.0,model); - //This is dangerous as it assumes a uniform time discretization - int index = (int) Math.round(evaluationTime / model.getTimeStep()); + public RandomVariable getValue(final double evaluationTime, final TreeModel model) { + final RandomVariable[] levels = getValues(0.0, model); + // This is dangerous as it assumes a uniform time discretization. + final int index = (int) Math.round(evaluationTime / model.getTimeStep()); return levels[index]; } /** - * Computes the vector of values at each time/state on the tree (via backward induction) + * Computes the vector of values at each time/state on the tree (via backward induction). * Implementations should: - * Validate inputs via validate(double, TreeModel). - * Return an array whose element at index k is a RandomVariable - * representing the value at time level k (with length equal to the number of states at that level). + *
    + *
  • Validate inputs via validate(double, TreeModel).
    • + *
  • Return an array whose element at index k is a RandomVariable representing the value at time level k + * (with length equal to the number of states at that level).
    • + *
* - * @param evaluationTime The time at which valuation is performed + * @param evaluationTime The time at which valuation is performed. * @param model The tree model (spot lattice, discounting, conditional expectation). * @return An array of value random variables, one per time level, indexed from 0 to maturity level. */ @@ -90,20 +91,21 @@ public RandomVariable getValue(double evaluationTime, TreeModel model) { * @param m The model (must not be null). * @throws IllegalArgumentException If {@code m == null} or {@code t < 0}. */ - protected void validate(double t, TreeModel m) { - if(m == null) throw new IllegalArgumentException("model null"); - if(t < 0) throw new IllegalArgumentException("evaluationTime < 0"); + protected void validate(final double t, final TreeModel m) { + if(m == null) { + throw new IllegalArgumentException("model null"); + } + if(t < 0) { + throw new IllegalArgumentException("evaluationTime < 0"); + } } - /** * Returns the contract maturity. * * @return The maturity T. */ - protected double getMaturity(){ + protected double getMaturity() { return maturity; } - } - diff --git a/src/main/java/net/finmath/tree/OneDimensionalRiskFactorTreeModel.java b/src/main/java/net/finmath/tree/OneDimensionalRiskFactorTreeModel.java index 70a3ded904..434dcf82b0 100644 --- a/src/main/java/net/finmath/tree/OneDimensionalRiskFactorTreeModel.java +++ b/src/main/java/net/finmath/tree/OneDimensionalRiskFactorTreeModel.java @@ -3,84 +3,91 @@ import java.util.ArrayList; import java.util.List; import java.util.function.DoubleUnaryOperator; + import net.finmath.montecarlo.RandomVariableFromDoubleArray; import net.finmath.stochastic.RandomVariable; /** - * This abstract class encapsulates the logic of one-dimensional trees for the simulation of a - * risk factor. At this level, the class could be used for equities or interest rates (e.g. - * the Hull-White short rate model). - * - * @author Carlo Andrea Tramentozzi, Alessandro Gnoatto + * This abstract class encapsulates the logic of one-dimensional trees for the simulation of a risk factor. + * At this level, the class could be used for equities or interest rates (e.g. the Hull-White short rate model). + * + * @author Carlo Andrea Tramentozzi + * @author Alessandro Gnoatto */ public abstract class OneDimensionalRiskFactorTreeModel implements TreeModel { - - /** Cache of level spot S_k */ + + /** Cache of level spot S_k. */ private final List riskFactorLevels = new ArrayList<>(); /** - * Return number of states at level k (binomial: k+1; trinomial: 2k+1) - * @param k - * @return number of states at level k + * Returns the number of states at level {@code k} (binomial: {@code k + 1}; trinomial: {@code 2k + 1}). + * + * @param k The level index. + * @return The number of states at level {@code k}. */ public abstract int statesAt(int k); /** - * Builds all the realizations X_k[i] - * @param k - * @return + * Builds all the realizations {@code X_k[i]}. + * + * @param k The level index. + * @return The spot (risk factor) values at level {@code k}. */ protected abstract RandomVariable buildSpotLevel(int k); - /** - * Discounted Conditional expectation: from v_{k+1} (array) to v_k (array) - * @param vNext - * @param k - * @return the conditonal expectation + * Discounted conditional expectation: from {@code v_{k+1}} to {@code v_k}. + * + * @param vNext The values at the next level. + * @param k The current level index. + * @return The conditional expectation. */ protected abstract RandomVariable conditionalExpectation(RandomVariable vNext, int k); - /** - * Builds (if missing) and return X_k as array, using cache. - * @param k - * @return the risk factor at level k + * Builds (if missing) and returns {@code X_k} using the cache. + * + * @param k The level index. + * @return The risk factor at level {@code k}. */ protected final RandomVariable ensureRiskFactorLevelArray(int k) { while(riskFactorLevels.size() <= k) { - int next = riskFactorLevels.size(); + final int next = riskFactorLevels.size(); riskFactorLevels.add(buildSpotLevel(next)); } return riskFactorLevels.get(k); } - + @Override - public RandomVariable getTransformedValuesAtGivenTimeRV(double time, DoubleUnaryOperator transformFunction) { + public RandomVariable getTransformedValuesAtGivenTimeRV(final double time, final DoubleUnaryOperator transformFunction) { int k = (int) Math.round(time / getTimeStep()); - // Checking param to avoid out of bound exception - if (k < 0) k = 0; - if (k > getNumberOfTimes() - 1) k = getNumberOfTimes() - 1; - RandomVariable sRV = ensureRiskFactorLevelArray(k); - double[] s = sRV.getRealizations(); - double[] v = new double[s.length]; - for (int i = 0; i < s.length; i++) { + // Checking param to avoid out of bound exception. + if(k < 0) { + k = 0; + } + if(k > getNumberOfTimes() - 1) { + k = getNumberOfTimes() - 1; + } + + final RandomVariable sRV = ensureRiskFactorLevelArray(k); + final double[] s = sRV.getRealizations(); + final double[] v = new double[s.length]; + for(int i = 0; i < s.length; i++) { v[i] = transformFunction.applyAsDouble(s[i]); } return new RandomVariableFromDoubleArray(k * getTimeStep(), v); } @Override - public RandomVariable getConditionalExpectationRV(RandomVariable vNext, int k) { - RandomVariable vK = conditionalExpectation(vNext, k); // hook + public RandomVariable getConditionalExpectationRV(final RandomVariable vNext, final int k) { + final RandomVariable vK = conditionalExpectation(vNext, k); // hook return new RandomVariableFromDoubleArray(k * getTimeStep(), vK.getRealizations()); } @Override - public RandomVariable getSpotAtGivenTimeIndexRV(int k) { - RandomVariable sRV = ensureRiskFactorLevelArray(k); + public RandomVariable getSpotAtGivenTimeIndexRV(final int k) { + final RandomVariable sRV = ensureRiskFactorLevelArray(k); return new RandomVariableFromDoubleArray(k * getTimeStep(), sRV.getRealizations()); } - } diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/AbstractRecombiningTreeModel.java b/src/main/java/net/finmath/tree/assetderivativevaluation/AbstractRecombiningTreeModel.java index f00539ee3c..da547a3898 100644 --- a/src/main/java/net/finmath/tree/assetderivativevaluation/AbstractRecombiningTreeModel.java +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/AbstractRecombiningTreeModel.java @@ -82,7 +82,7 @@ public double getOneStepDiscountFactor(int timeIndex) { return Math.exp(-riskFreeRate * timeStep); } - + /** Getter */ @Override @@ -91,7 +91,7 @@ public double getOneStepDiscountFactor(int timeIndex) { public double getRiskFreeRate() { return riskFreeRate; } @Override public double getVolatility() { return volatility; } - + public double getTimeStep() { return timeStep; } public double getLastTime() { return lastTime; } public int getNumberOfTimes() { return numberOfTimes; } diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/OneDimensionalAssetTreeModel.java b/src/main/java/net/finmath/tree/assetderivativevaluation/OneDimensionalAssetTreeModel.java index 0c1c67bb65..44dc629156 100644 --- a/src/main/java/net/finmath/tree/assetderivativevaluation/OneDimensionalAssetTreeModel.java +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/OneDimensionalAssetTreeModel.java @@ -3,27 +3,29 @@ /** * From this level of abstraction, we mean that the tree should represent a stock price * and not, for example, an interest rate. - * + * * @author Alessandro Gnoatto */ public interface OneDimensionalAssetTreeModel { /** - * The risk factor is a stock - * @return the initial stock price + * The risk factor is a stock. + * + * @return The initial stock price. */ - public double getInitialPrice(); + double getInitialPrice(); /** - * The risk free rate - * @return the risk free rate + * The risk free rate. + * + * @return The risk free rate. */ - public double getRiskFreeRate(); + double getRiskFreeRate(); /** - * Returns the volatility - * @return the volatility + * Returns the volatility. + * + * @return The volatility. */ - public double getVolatility(); - + double getVolatility(); } diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/models/CoxRossRubinsteinModel.java b/src/main/java/net/finmath/tree/assetderivativevaluation/models/CoxRossRubinsteinModel.java index 242cbed19c..4de8966886 100644 --- a/src/main/java/net/finmath/tree/assetderivativevaluation/models/CoxRossRubinsteinModel.java +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/models/CoxRossRubinsteinModel.java @@ -17,7 +17,7 @@ public class CoxRossRubinsteinModel extends AbstractRecombiningTreeModel { private final double u; private final double d; private final double q; - + /** * It constructs an object which represents the approximation of a Black-Scholes model via the Cox Ross * Rubinstein model. diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/models/package-info.java b/src/main/java/net/finmath/tree/assetderivativevaluation/models/package-info.java new file mode 100644 index 0000000000..f88a571e16 --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/models/package-info.java @@ -0,0 +1,6 @@ +/** + * Tree models for the valuation of equity derivatives. + * + * @author Alessandro Gnoatto + */ +package net.finmath.tree.assetderivativevaluation.models; diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/package-info.java b/src/main/java/net/finmath/tree/assetderivativevaluation/package-info.java new file mode 100644 index 0000000000..abb835eef0 --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/package-info.java @@ -0,0 +1,6 @@ +/** + * Provides algorithms related to equity derivative valuation via multinomial trees. + * + * @author Alessandro Gnoatto + */ +package net.finmath.tree.assetderivativevaluation; diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractPathDependentProduct.java b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractPathDependentProduct.java index 5c431302af..8df66a4f92 100644 --- a/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractPathDependentProduct.java +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractPathDependentProduct.java @@ -2,8 +2,8 @@ import java.util.Arrays; -import net.finmath.stochastic.RandomVariable; import net.finmath.montecarlo.RandomVariableFromDoubleArray; +import net.finmath.stochastic.RandomVariable; import net.finmath.tree.AbstractTreeProduct; import net.finmath.tree.TreeModel; @@ -11,325 +11,272 @@ * Base class for full (non-recombining) path-dependent products on a tree. * * Key property: - * - The number of nodes grows exponentially with the number of time steps: - * level j has branchingFactor^j nodes. - * - NO recombination / compression / bucketing. + *
    + *
  • The number of nodes grows exponentially with the number of time steps: + * level j has branchingFactor^j nodes.
    • + *
  • No recombination / compression / bucketing.
    • + *
* - * This class builds a "full path tree" of: - * - an underlying recombining state index (for reading spot from the model's spot level), - * - plus product-specific path state variables (e.g. running sum/count for Asian). We build the - * full (non-recombining) path tree (size branchingFactor^j after j steps) + * This class builds a full path tree of: + *
    + *
  • An underlying recombining state index (for reading spot from the model's spot level).
    • + *
  • Product-specific path state variables (e.g. running sum/count for Asian).
    • + *
* - * After the forward loop, it then prices by backward induction using TreeModel's transition API: - * - p = model.getTransitionProbability(timeIndex, recombStateIndex, branchIndex) - * - df = model.getOneStepDiscountFactor(timeIndex) + * After the forward loop, it prices by backward induction using the {@link TreeModel}'s transition API: + *
    + *
  • p = model.getTransitionProbability(timeIndex, recombStateIndex, branchIndex)
    • + *
  • df = model.getOneStepDiscountFactor(timeIndex)
    • + *
* * Design choices: - * - We keep state variables in primitive arrays for speed and to avoid per-node objects. - * - We treat the model's spot as a function of (timeIndex, recombStateIndex) so we don't - * need the model to expose explicit "spot multipliers". + *
    + *
  • State variables are kept in primitive arrays for speed and to avoid per-node objects.
    • + *
  • The model's spot is treated as a function of (timeIndex, recombStateIndex), so no explicit + * spot multipliers are required.
    • + *
* * IMPORTANT: - * - To evolve the recombining state index along a branch we require a "child state shift" - * convention. For the existing models (CRR/JR/Boyle) this matches their internal layout: - * Binomial: childShift = {0, 1} (up keeps index, down increases index) - * Trinomial (Boyle): childShift = {0, 1, 2} (up, mid, down) - * For a general multinomial model, one has to provide childShift explicitly. - * - * @author Alessandro Gnoatto + *
    + *
  • To evolve the recombining state index along a branch, we require a "child state shift" convention. + * For the existing models (CRR/JR/Boyle) this matches their internal layout: + *
      + *
    • Binomial: childShift = {0, 1}
      • + *
    • Trinomial (Boyle): childShift = {0, 1, 2}
      • + *
    + * For a general multinomial model, childShift must be provided explicitly.
    • + *
+ * + * @author Alessandro Gnoatto */ public abstract class AbstractPathDependentProduct extends AbstractTreeProduct { - /** - * Fixing time indices on the model grid, e.g. {1,2,3,4,5}. - * These determine when path-state should be updated with spot. - */ - private final int[] fixingTimeIndices; - - - protected AbstractPathDependentProduct( - final double maturity, - final int[] fixingTimeIndices - ) { - super(maturity); - this.fixingTimeIndices = fixingTimeIndices != null ? fixingTimeIndices.clone() : new int[0]; - Arrays.sort(this.fixingTimeIndices); - } - - /** - * Number of (double-valued) path-state variables required by the product. - * Examples: - * - Fixed-strike Asian: 2 variables (sum, count) - * - Lookback (min): 1 variable (runningMin) - * - Lookback (min/max): 2 variables - */ - protected abstract int getNumberOfStateVariables(); - - /** - * Initialize product path-state at the root node. - * - * @param spotAtNode spot at (timeIndex, recombStateIndex) - * @param timeIndex time index for root (typically k0) - * @param isFixing whether this time index is part of averaging/fixing set - * @param stateOut output array length = getNumberOfStateVariables() - */ - protected abstract void initializeState( - double spotAtNode, - int timeIndex, - boolean isFixing, - double[] stateOut - ); - - /** - * Evolve product path-state from parent to a child node. - * - * @param parentState state at parent node - * @param spotAtChild spot at child node - * @param timeIndexChild time index of child - * @param isFixingChild whether child time is fixing - * @param childStateOut output state array (length = getNumberOfStateVariables()) - */ - protected abstract void evolveState( - double[] parentState, - double spotAtChild, - int timeIndexChild, - boolean isFixingChild, - double[] childStateOut - ); - - /** - * Terminal payoff at maturity. - * - * @param terminalState state variables at maturity - * @param spotAtMaturity spot at maturity node (often needed e.g. floating-strike Asian, lookback) - * @return payoff - */ - protected abstract double payoff( - double[] terminalState, - double spotAtMaturity - ); - - /** - * Converts a double time to a time index by rounding. - */ - protected final int timeToIndex(double time, TreeModel model) { - return (int) Math.round(time / model.getTimeStep()); - } - - /** - * Returns whether a given time index is a fixing time. - */ - protected final boolean isFixingTimeIndex(int timeIndex) { - return Arrays.binarySearch(fixingTimeIndices, timeIndex) >= 0; - } - - /** - * Determine branching factor and childStateShift convention. - */ - private int[] resolveChildShift(TreeModel model, int timeIndex, int stateIndex) { - int branchingFactor = model.getNumberOfBranches(timeIndex, stateIndex); - int[] childStateShift = model.getChildStateIndexShift(); - if(childStateShift != null) { - if(childStateShift.length != branchingFactor) { - throw new IllegalArgumentException( - "childStateShift length (" + childStateShift.length + ") " - + "does not match branching factor (" + branchingFactor + ")." - ); - } - return childStateShift; - } - - throw new IllegalArgumentException( - "No childStateShift provided for branching factor " + branchingFactor - + "." - ); - } - - /** - * Full non-recombining valuation. - * - * Notes about evaluationTime: - * - Path dependent contracts generally require the past path-state at evaluationTime. - * - This implementation assumes the path-state STARTS at evaluationTime (i.e. no past). - * This is exactly what you want when evaluationTime = 0. - * - * IMPORTANT WARNING: - * - * In the current structure of the library this method will always be called with evaluationTime = 0. - * When evaluationTime > 0. the logic contained in AbstractTreeProduct guarantees that the whole tree of the - * path dependent functional is reconstructed starting from time zero. After this the time instants before the - * requested time are discarded and the RandomVariable at the selected time is returned. Remember that we enjoy the - * Markov property only by extending the state space, i.e. with respect to the couple (stock, pathDependentFunctional) - */ - @Override - public final RandomVariable[] getValues(final double evaluationTime, final TreeModel model) { - - // Ensure transition API is supported - // (if not, TreeModel default methods will throw). - final int k0 = timeToIndex(evaluationTime, model); - final int n = timeToIndex(getMaturity(), model); - if(n < k0) throw new IllegalArgumentException("Maturity is before evaluation time."); - - // Root uses recombining state index 0 (consistent with European products). - final int rootRecombState = 0; - - // Determine branching and child shift convention at the root step. - // We assume constant branching factor. - final int[] childShift = resolveChildShift(model, k0, rootRecombState); - final int B = childShift.length; - - final int steps = n - k0; - - // levels[j] corresponds to timeIndex = k0 + j, and has B^j nodes. - final RandomVariable[] levels = new RandomVariable[steps + 1]; - - // Store recombining index per full node, and product state variables per node. - // recombIndex[j][node] = recombining state index at model time k0+j. - final int[][] recombIndex = new int[steps + 1][]; - final double[][] state = new double[steps + 1][]; // flattened: node-major, then var-major - - final int m = getNumberOfStateVariables(); - - // Allocate level 0 - recombIndex[0] = new int[] { rootRecombState }; - state[0] = new double[m]; // 1 node * m vars - - // Spot at root - final double spot0 = model.getSpotAtGivenTimeIndexRV(k0).get(rootRecombState); - - initializeState(spot0, k0, isFixingTimeIndex(k0), state[0]); - - // Forward expansion of path-state (full tree) - for(int j = 1; j <= steps; j++) { - final int timeIndex = k0 + j; + /** + * Fixing time indices on the model grid, e.g. {1,2,3,4,5}. + * These determine when path-state should be updated with spot. + */ + private final int[] fixingTimeIndices; + + protected AbstractPathDependentProduct( + final double maturity, + final int[] fixingTimeIndices) { + super(maturity); + this.fixingTimeIndices = fixingTimeIndices != null ? fixingTimeIndices.clone() : new int[0]; + Arrays.sort(this.fixingTimeIndices); + } + + /** + * Number of (double-valued) path-state variables required by the product. + */ + protected abstract int getNumberOfStateVariables(); + + /** + * Initialize product path-state at the root node. + */ + protected abstract void initializeState( + double spotAtNode, + int timeIndex, + boolean isFixing, + double[] stateOut); + + /** + * Evolve product path-state from parent to a child node. + */ + protected abstract void evolveState( + double[] parentState, + double spotAtChild, + int timeIndexChild, + boolean isFixingChild, + double[] childStateOut); + + /** + * Terminal payoff at maturity. + */ + protected abstract double payoff( + double[] terminalState, + double spotAtMaturity); + + /** + * Converts a double time to a time index by rounding. + */ + protected final int timeToIndex(final double time, final TreeModel model) { + return (int) Math.round(time / model.getTimeStep()); + } + + /** + * Returns whether a given time index is a fixing time. + */ + protected final boolean isFixingTimeIndex(final int timeIndex) { + return Arrays.binarySearch(fixingTimeIndices, timeIndex) >= 0; + } + + /** + * Determine branching factor and childStateShift convention. + */ + private int[] resolveChildShift(final TreeModel model, final int timeIndex, final int stateIndex) { + final int branchingFactor = model.getNumberOfBranches(timeIndex, stateIndex); + final int[] childStateShift = model.getChildStateIndexShift(); + + if(childStateShift != null) { + if(childStateShift.length != branchingFactor) { + throw new IllegalArgumentException( + "childStateShift length (" + childStateShift.length + ") " + + "does not match branching factor (" + branchingFactor + ")."); + } + return childStateShift; + } + + throw new IllegalArgumentException( + "No childStateShift provided for branching factor " + branchingFactor + "."); + } + + @Override + public final RandomVariable[] getValues(final double evaluationTime, final TreeModel model) { + + final int k0 = timeToIndex(evaluationTime, model); + final int n = timeToIndex(getMaturity(), model); + if(n < k0) { + throw new IllegalArgumentException("Maturity is before evaluation time."); + } + + final int rootRecombState = 0; + + final int[] childShift = resolveChildShift(model, k0, rootRecombState); + final int B = childShift.length; + + final int steps = n - k0; + + final RandomVariable[] levels = new RandomVariable[steps + 1]; + + final int[][] recombIndex = new int[steps + 1][]; + final double[][] state = new double[steps + 1][]; + + final int m = getNumberOfStateVariables(); + + recombIndex[0] = new int[] { rootRecombState }; + state[0] = new double[m]; + + final double spot0 = model.getSpotAtGivenTimeIndexRV(k0).get(rootRecombState); + initializeState(spot0, k0, isFixingTimeIndex(k0), state[0]); + + for(int j = 1; j <= steps; j++) { + + final int timeIndex = k0 + j; + + final int nodes = powInt(B, j); + final int parentNodes = powInt(B, j - 1); + + recombIndex[j] = new int[nodes]; + state[j] = new double[nodes * m]; - final int nodes = powInt(B, j); - final int parentNodes = powInt(B, j - 1); + final RandomVariable spotRV = model.getSpotAtGivenTimeIndexRV(timeIndex); - recombIndex[j] = new int[nodes]; - state[j] = new double[nodes * m]; + for(int parent = 0; parent < parentNodes; parent++) { - // Pre-fetch the recombining spot vector at this time index. - // We will index into it using recombIndex[j][node]. - final var spotRV = model.getSpotAtGivenTimeIndexRV(timeIndex); - - for(int parent = 0; parent < parentNodes; parent++) { - - final int parentRecomb = recombIndex[j - 1][parent]; - - // parent state slice - final int parentStateOffset = parent * m; - - for(int b = 0; b < B; b++) { + final int parentRecomb = recombIndex[j - 1][parent]; + final int parentStateOffset = parent * m; - final int child = parent * B + b; + for(int b = 0; b < B; b++) { - final int childRecomb = parentRecomb + childShift[b]; - recombIndex[j][child] = childRecomb; + final int child = parent * B + b; + final int childRecomb = parentRecomb + childShift[b]; - final double spotChild = spotRV.get(childRecomb); + recombIndex[j][child] = childRecomb; - // child state slice - final int childStateOffset = child * m; + final double spotChild = spotRV.get(childRecomb); + final int childStateOffset = child * m; - // evolve state - evolveState( - state[j - 1], - parentStateOffset, - spotChild, - timeIndex, - isFixingTimeIndex(timeIndex), - state[j], - childStateOffset - ); - } - } - } + evolveState( + state[j - 1], + parentStateOffset, + spotChild, + timeIndex, + isFixingTimeIndex(timeIndex), + state[j], + childStateOffset); + } + } + } - // Terminal payoff values at maturity level - final int terminalNodes = powInt(B, steps); - final double[] vTerminal = new double[terminalNodes]; + final int terminalNodes = powInt(B, steps); + final double[] vTerminal = new double[terminalNodes]; - final var spotN = model.getSpotAtGivenTimeIndexRV(n); + final RandomVariable spotN = model.getSpotAtGivenTimeIndexRV(n); - for(int node = 0; node < terminalNodes; node++) { - final int recomb = recombIndex[steps][node]; - final double spotAtMaturity = spotN.get(recomb); + for(int node = 0; node < terminalNodes; node++) { - final double[] tmp = new double[m]; - System.arraycopy(state[steps], node * m, tmp, 0, m); + final int recomb = recombIndex[steps][node]; + final double spotAtMaturity = spotN.get(recomb); - vTerminal[node] = payoff(tmp, spotAtMaturity); - } + final double[] tmp = new double[m]; + System.arraycopy(state[steps], node * m, tmp, 0, m); - levels[steps] = new RandomVariableFromDoubleArray(n * model.getTimeStep(), vTerminal); + vTerminal[node] = payoff(tmp, spotAtMaturity); + } - // Backward induction on full node set - double[] vNext = vTerminal; + levels[steps] = new RandomVariableFromDoubleArray(n * model.getTimeStep(), vTerminal); - for(int j = steps - 1; j >= 0; j--) { + double[] vNext = vTerminal; - final int timeIndex = k0 + j; - final int nodesHere = powInt(B, j); + for(int j = steps - 1; j >= 0; j--) { - final double[] vHere = new double[nodesHere]; + final int timeIndex = k0 + j; + final int nodesHere = powInt(B, j); - final double df = model.getOneStepDiscountFactor(timeIndex); + final double[] vHere = new double[nodesHere]; + final double df = model.getOneStepDiscountFactor(timeIndex); - for(int node = 0; node < nodesHere; node++) { + for(int node = 0; node < nodesHere; node++) { - final int parentRecomb = recombIndex[j][node]; + final int parentRecomb = recombIndex[j][node]; + double sum = 0.0; - double sum = 0.0; + for(int b = 0; b < B; b++) { - for(int b = 0; b < B; b++) { + final int child = node * B + b; + final double p = model.getTransitionProbability(timeIndex, parentRecomb, b); - final int child = node * B + b; + sum += p * vNext[child]; + } - final double p = model.getTransitionProbability(timeIndex, parentRecomb, b); + vHere[node] = df * sum; + } - sum += p * vNext[child]; - } + levels[j] = new RandomVariableFromDoubleArray(timeIndex * model.getTimeStep(), vHere); + vNext = vHere; + } - vHere[node] = df * sum; - } + return levels; + } - levels[j] = new RandomVariableFromDoubleArray(timeIndex * model.getTimeStep(), vHere); - vNext = vHere; - } + private void evolveState( + final double[] parentStateFlat, + final int parentOffset, + final double spotChild, + final int timeIndexChild, + final boolean isFixingChild, + final double[] childStateFlat, + final int childOffset) { - return levels; - } + final int m = getNumberOfStateVariables(); - /** - * Helper: evolveState where parent/child arrays are flattened. - * This avoids allocating per-node double[] objects. - */ - private void evolveState( - final double[] parentStateFlat, - final int parentOffset, - final double spotChild, - final int timeIndexChild, - final boolean isFixingChild, - final double[] childStateFlat, - final int childOffset - ) { - final int m = getNumberOfStateVariables(); - final double[] parent = new double[m]; - System.arraycopy(parentStateFlat, parentOffset, parent, 0, m); + final double[] parent = new double[m]; + System.arraycopy(parentStateFlat, parentOffset, parent, 0, m); - final double[] child = new double[m]; - evolveState(parent, spotChild, timeIndexChild, isFixingChild, child); + final double[] child = new double[m]; + evolveState(parent, spotChild, timeIndexChild, isFixingChild, child); - System.arraycopy(child, 0, childStateFlat, childOffset, m); - } + System.arraycopy(child, 0, childStateFlat, childOffset, m); + } - private static int powInt(int base, int exp) { - if(exp < 0) throw new IllegalArgumentException("exp < 0"); - int r = 1; - for(int i = 0; i < exp; i++) { - r = Math.multiplyExact(r, base); // throws if overflow - } - return r; - } + private static int powInt(final int base, final int exp) { + if(exp < 0) { + throw new IllegalArgumentException("exp < 0"); + } + int r = 1; + for(int i = 0; i < exp; i++) { + r = Math.multiplyExact(r, base); + } + return r; + } } diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/products/AsianOption.java b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AsianOption.java index cacaebf221..06f93423be 100644 --- a/src/main/java/net/finmath/tree/assetderivativevaluation/products/AsianOption.java +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AsianOption.java @@ -4,11 +4,15 @@ * Fixed-strike Asian option using FULL non-recombining tree expansion. * * State variables: - * - state[0] = runningSum of fixing spots - * - state[1] = fixingCount + *
    + *
  • state[0] = running sum of fixing spots
    • + *
  • state[1] = fixing count
    • + *
* * Payoff at maturity: - * - max( (sum/count) - K, 0 ) + *
    + *
  • max((sum / count) - K, 0)
    • + *
* * The averaging can be performed on any subset of the time grid by specifying fixingTimeIndices. * @@ -16,59 +20,56 @@ */ public class AsianOption extends AbstractPathDependentProduct { - private final double strike; + private final double strike; - /** - * @param maturity maturity in model time units - * @param strike strike K - * @param fixingTimeIndices subset of model time indices included in the average - * @param childStateShift optional mapping for recombining index evolution; - * pass null for defaults (binomial: {0,1}, trinomial: {0,1,2}) - */ - public AsianOption( - final double maturity, - final double strike, - final int[] fixingTimeIndices - ) { - super(maturity, fixingTimeIndices); - this.strike = strike; - } + /** + * @param maturity Maturity in model time units. + * @param strike Strike K. + * @param fixingTimeIndices Subset of model time indices included in the average. + */ + public AsianOption( + final double maturity, + final double strike, + final int[] fixingTimeIndices) { + super(maturity, fixingTimeIndices); + this.strike = strike; + } - @Override - protected int getNumberOfStateVariables() { - return 2; // sum, count - } + @Override + protected int getNumberOfStateVariables() { + return 2; // sum, count + } - @Override - protected void initializeState( - final double spotAtNode, - final int timeIndex, - final boolean isFixing, - final double[] stateOut - ) { - stateOut[0] = isFixing ? spotAtNode : 0.0; // sum - stateOut[1] = isFixing ? 1.0 : 0.0; // count stored as double for simplicity - } + @Override + protected void initializeState( + final double spotAtNode, + final int timeIndex, + final boolean isFixing, + final double[] stateOut) { - @Override - protected void evolveState( - final double[] parentState, - final double spotAtChild, - final int timeIndexChild, - final boolean isFixingChild, - final double[] childStateOut - ) { - childStateOut[0] = parentState[0] + (isFixingChild ? spotAtChild : 0.0); - childStateOut[1] = parentState[1] + (isFixingChild ? 1.0 : 0.0); - } + stateOut[0] = isFixing ? spotAtNode : 0.0; // sum + stateOut[1] = isFixing ? 1.0 : 0.0; // count stored as double for simplicity + } - @Override - protected double payoff(final double[] terminalState, final double spotAtMaturity) { - final double sum = terminalState[0]; - final double cnt = terminalState[1]; + @Override + protected void evolveState( + final double[] parentState, + final double spotAtChild, + final int timeIndexChild, + final boolean isFixingChild, + final double[] childStateOut) { - final double avg = (cnt > 0.0) ? (sum / cnt) : 0.0; + childStateOut[0] = parentState[0] + (isFixingChild ? spotAtChild : 0.0); + childStateOut[1] = parentState[1] + (isFixingChild ? 1.0 : 0.0); + } - return Math.max(avg - strike, 0.0); - } + @Override + protected double payoff(final double[] terminalState, final double spotAtMaturity) { + final double sum = terminalState[0]; + final double cnt = terminalState[1]; + + final double avg = (cnt > 0.0) ? (sum / cnt) : 0.0; + + return Math.max(avg - strike, 0.0); + } } diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/products/FloatingStrikeAsianOption.java b/src/main/java/net/finmath/tree/assetderivativevaluation/products/FloatingStrikeAsianOption.java index 4df0e8782d..a35fc1679f 100644 --- a/src/main/java/net/finmath/tree/assetderivativevaluation/products/FloatingStrikeAsianOption.java +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/FloatingStrikeAsianOption.java @@ -4,73 +4,76 @@ * Floating-strike Asian option (CALL) on a full non-recombining tree (exponential growth). * * Payoff at maturity: - * - * max( S(T) - A , 0 ) - * + * 
+ *     max(S(T) - A, 0)
+ *
* where A is the arithmetic average of the spot over the fixingTimeIndices. * * State variables: - * - state[0] = running sum of fixing spots - * - state[1] = fixing count + *
    + *
  • state[0] = running sum of fixing spots
    • + *
  • state[1] = fixing count
    • + *
* * No recombination / compression: - * - number of nodes at step j is B^j (B = branching factor of the model). - * - * @author Alessandro Gnoatto + *
    + *
  • Number of nodes at step j is B^j (B = branching factor of the model).
    • + *
+ * + * @author Alessandro Gnoatto */ public class FloatingStrikeAsianOption extends AbstractPathDependentProduct { - /** - * @param maturity maturity in model time units - * @param fixingTimeIndices subset of model time indices included in the average - */ - public FloatingStrikeAsianOption( - final double maturity, - final int[] fixingTimeIndices - ) { - super(maturity, fixingTimeIndices); - } + /** + * @param maturity Maturity in model time units. + * @param fixingTimeIndices Subset of model time indices included in the average. + */ + public FloatingStrikeAsianOption( + final double maturity, + final int[] fixingTimeIndices) { + super(maturity, fixingTimeIndices); + } + + @Override + protected int getNumberOfStateVariables() { + return 2; // running sum, fixing count + } + + @Override + protected void initializeState( + final double spotAtNode, + final int timeIndex, + final boolean isFixing, + final double[] stateOut) { + + stateOut[0] = isFixing ? spotAtNode : 0.0; // sum + stateOut[1] = isFixing ? 1.0 : 0.0; // count (stored as double) + } - @Override - protected int getNumberOfStateVariables() { - return 2; // running sum, fixing count - } + @Override + protected void evolveState( + final double[] parentState, + final double spotAtChild, + final int timeIndexChild, + final boolean isFixingChild, + final double[] childStateOut) { - @Override - protected void initializeState( - final double spotAtNode, - final int timeIndex, - final boolean isFixing, - final double[] stateOut - ) { - stateOut[0] = isFixing ? spotAtNode : 0.0; // sum - stateOut[1] = isFixing ? 1.0 : 0.0; // count (stored as double) - } + childStateOut[0] = parentState[0] + (isFixingChild ? spotAtChild : 0.0); + childStateOut[1] = parentState[1] + (isFixingChild ? 1.0 : 0.0); + } - @Override - protected void evolveState( - final double[] parentState, - final double spotAtChild, - final int timeIndexChild, - final boolean isFixingChild, - final double[] childStateOut - ) { - childStateOut[0] = parentState[0] + (isFixingChild ? spotAtChild : 0.0); - childStateOut[1] = parentState[1] + (isFixingChild ? 1.0 : 0.0); - } + @Override + protected double payoff( + final double[] terminalState, + final double spotAtMaturity) { - @Override - protected double payoff( - final double[] terminalState, - final double spotAtMaturity - ) { - final double sum = terminalState[0]; - final double cnt = terminalState[1]; + final double sum = terminalState[0]; + final double cnt = terminalState[1]; - // Defensive guard: empty fixing set - final double average = (cnt > 0.0) ? (sum / cnt) : 0.0; + // Defensive guard: empty fixing set. + final double average = (cnt > 0.0) ? (sum / cnt) : 0.0; - // Floating-strike CALL payoff - return Math.max(spotAtMaturity - average, 0.0); - } + // Floating-strike CALL payoff. + return Math.max(spotAtMaturity - average, 0.0); + } } diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/products/package-info.java b/src/main/java/net/finmath/tree/assetderivativevaluation/products/package-info.java new file mode 100644 index 0000000000..51919101e8 --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/package-info.java @@ -0,0 +1,6 @@ +/** + * Product classes related to equity derivative valuation via multinomial trees. + * + * @author Alessandro Gnoatto + */ +package net.finmath.tree.assetderivativevaluation.products; diff --git a/src/main/java/net/finmath/tree/package-info.java b/src/main/java/net/finmath/tree/package-info.java new file mode 100644 index 0000000000..b210d1e61c --- /dev/null +++ b/src/main/java/net/finmath/tree/package-info.java @@ -0,0 +1,6 @@ +/** + * Provides algorithms related to derivative valuation via multinomial trees. + * + * @author Alessandro Gnoatto + */ +package net.finmath.tree; From c66e5ace413ef903cae1912741cd074db7d630ac Mon Sep 17 00:00:00 2001 From: Alessandro Gnoatto Date: Wed, 18 Feb 2026 23:08:48 +0100 Subject: [PATCH 6/7] revert to original Heston class --- .../net/finmath/functions/HestonModel.java | 704 ++++++++++-------- 1 file changed, 409 insertions(+), 295 deletions(-) diff --git a/src/main/java/net/finmath/functions/HestonModel.java b/src/main/java/net/finmath/functions/HestonModel.java index 2c8e150821..20579077fe 100644 --- a/src/main/java/net/finmath/functions/HestonModel.java +++ b/src/main/java/net/finmath/functions/HestonModel.java @@ -13,347 +13,421 @@ /** * This class implements some functions as static class methods related to the Heston model. - * The calculation is performed by means of the FFT algorithm of Carr Madan applied to the gradient of the - * characteristic funtion. + * The calculation is performed by means of the FFT algorithm of Carr Madan applied to the gradient of the characteristic funtion. + * + * The model is + * \[ + * dS(t) = r^{\text{c}}(t) S(t) dt + \sqrt{V(t)} S(t) dW_{1}(t), \quad S(0) = S_{0}, + * \] + * \[ + * dV(t) = \kappa ( \theta - V(t) ) dt + \xi \sqrt{V(t)} dW_{2}(t), \quad V(0) = \sigma^2, + * \] + * \[ + * dW_{1} dW_{2} = \rho dt + * \] + * \[ + * dN(t) = r^{\text{d}}(t) N(t) dt, \quad N(0) = N_{0}, + * \] + * where \( W \) is a Brownian motion. + * + * The free parameters of this model are: + *
+ *
\( S_{0} \)
spot - initial value of S
+ *
\( r^{\text{c}} \)
the drift of the stock, equals difference of repo rate and divident yield
+ *
\( \sigma \)
the initial volatility level (initial variance \( V_{0} = \sigma^{2} \))
+ *
\( r^{\text{d}} \)
the discount rate
+ *
\( \xi \)
the volatility of volatility
+ *
\( \theta \)
the mean reversion level of the stochastic volatility
+ *
\( \kappa \)
the mean reversion speed of the stochastic volatility
+ *
\( \rho \)
the correlation of the Brownian drivers
+ *
* * @author Alessandro Gnoatto */ public class HestonModel { /* - * Parameters for the FFT calculator. + * Parameters for the FFT calculator */ - private static final InterpolationMethod INT_METHOD = InterpolationMethod.LINEAR; - private static final ExtrapolationMethod EXT_METHOD = ExtrapolationMethod.CONSTANT; - private static final int NUMBER_OF_POINTS = 4096 * 2; - private static final double GRID_SPACING = 0.4; + final static InterpolationMethod intMethod =InterpolationMethod.LINEAR; + final static ExtrapolationMethod extMethod = ExtrapolationMethod.CONSTANT; + final static int numberOfPoints = 4096*2; + final static double gridSpacing = 0.4; - enum HestonGreek {DELTA, GAMMA, THETA, RHO, VEGA1, VANNA, VOLGA} + private enum HestonGreek {DELTA, GAMMA, THETA, RHO, VEGA1, VANNA, VOLGA}; /** * Calculates the delta of a call option under a Heston model. - * + * * @param initialStockValue Initital value of the stock. - * @param riskFreeRate The risk free rate. - * @param dividendYield The dividend yield. - * @param kappa The speed of mean reversion. - * @param theta The long run mean of the volatility. - * @param sigma The volatility of variance. - * @param v0 The initial instantaneous variance. - * @param rho Correlation between the two Brownian motions. - * @param optionMaturity The maturity of the option. - * @param optionStrike The strike of the option. + * @param riskFreeRate The risk free rate (the drift of the forward, repoRate - dividendYield = r-q). + * @param discountRate The discount rate. + * @param sigma the square root of the initial instantaneous variance (\( V_0 = sigma^2 \)) + * @param theta the long run mean of the volatility. + * @param kappa the speed of mean reversion. + * @param xi the volatility of variance. + * @param rho correlation between the two Brownian motions + * @param optionMaturity the maturity of the option + * @param optionStrike the strike of the option. * @return The delta of the option. */ public static double hestonOptionDelta( final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double discountRate, + final double sigma, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, - final double optionStrike) { + final double optionStrike) + { - final HestonGreek whatToCompute = HestonGreek.DELTA; + HestonGreek whatToCompute = HestonGreek.DELTA; - return hestonGreekCalculator( - initialStockValue, + return hestonGreekCalculator(initialStockValue, riskFreeRate, - dividendYield, - kappa, - theta, - sigma, - v0, + discountRate, + sigma, + theta, + kappa, + xi, rho, optionMaturity, optionStrike, whatToCompute, - NUMBER_OF_POINTS, - GRID_SPACING); + numberOfPoints, + gridSpacing); } /** - * Calculates the gamma of a call option under a Heston model. - * + * Calculates the gamma of a call option under a Heston model + * * @param initialStockValue Initital value of the stock. - * @param riskFreeRate The risk free rate. - * @param dividendYield The dividend yield. - * @param kappa The speed of mean reversion. - * @param theta The long run mean of the volatility. - * @param sigma The volatility of variance. - * @param v0 The initial instantaneous variance. - * @param rho Correlation between the two Brownian motions. - * @param optionMaturity The maturity of the option. - * @param optionStrike The strike of the option. - * @return The gamma of the option. + * @param riskFreeRate The risk free rate (the drift of the forward, repoRate - dividendYield = r-q). + * @param discountRate The discount rate. + * @param sigma the square root of the initial instantaneous variance (\( V_0 = sigma^2 \)) + * @param theta the long run mean of the volatility. + * @param kappa the speed of mean reversion. + * @param xi the volatility of variance. + * @param rho correlation between the two Brownian motions + * @param optionMaturity the maturity of the option + * @param optionStrike the strike of the option. + * @return The gamma of the option */ public static double hestonOptionGamma( final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double discountRate, + final double sigma, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, - final double optionStrike) { - - final HestonGreek whatToCompute = HestonGreek.GAMMA; + final double optionStrike) + { + HestonGreek whatToCompute = HestonGreek.GAMMA; - return hestonGreekCalculator( - initialStockValue, + return hestonGreekCalculator(initialStockValue, riskFreeRate, - dividendYield, - kappa, - theta, - sigma, - v0, + discountRate, + sigma, + theta, + kappa, + xi, rho, optionMaturity, optionStrike, whatToCompute, - NUMBER_OF_POINTS, - GRID_SPACING); + numberOfPoints, + gridSpacing); } /** - * Calculates the theta of a call option under a Heston model. - * + * Calculates the theta of a call option under a Heston model + * * @param initialStockValue Initital value of the stock. - * @param riskFreeRate The risk free rate. - * @param dividendYield The dividend yield. - * @param kappa The speed of mean reversion. - * @param theta The long run mean of the volatility. - * @param sigma The volatility of variance. - * @param v0 The initial instantaneous variance. - * @param rho Correlation between the two Brownian motions. - * @param optionMaturity The maturity of the option. - * @param optionStrike The strike of the option. - * @return The theta of the option. + * @param riskFreeRate The risk free rate (the drift of the forward, repoRate - dividendYield = r-q). + * @param discountRate The discount rate. + * @param sigma the square root of the initial instantaneous variance (\( V_0 = sigma^2 \)) + * @param theta the long run mean of the volatility. + * @param kappa the speed of mean reversion. + * @param xi the volatility of variance. + * @param rho correlation between the two Brownian motions + * @param optionMaturity the maturity of the option + * @param optionStrike the strike of the option. + * @return The theta of the option */ public static double hestonOptionTheta( final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double discountRate, + final double sigma, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, - final double optionStrike) { + final double optionStrike) + { + HestonGreek whatToCompute = HestonGreek.THETA; - final HestonGreek whatToCompute = HestonGreek.THETA; - - return hestonGreekCalculator( - initialStockValue, + return hestonGreekCalculator(initialStockValue, riskFreeRate, - dividendYield, - kappa, - theta, - sigma, - v0, + discountRate, + sigma, + theta, + kappa, + xi, rho, optionMaturity, optionStrike, whatToCompute, - NUMBER_OF_POINTS, - GRID_SPACING); + numberOfPoints, + gridSpacing); } /** - * Calculates the rho of a call option under a Heston model. - * + * Calculates the rho of a call option under a Heston model + * * @param initialStockValue Initital value of the stock. - * @param riskFreeRate The risk free rate. - * @param dividendYield The dividend yield. - * @param kappa The speed of mean reversion. - * @param theta The long run mean of the volatility. - * @param sigma The volatility of variance. - * @param v0 The initial instantaneous variance. - * @param rho Correlation between the two Brownian motions. - * @param optionMaturity The maturity of the option. - * @param optionStrike The strike of the option. - * @return The rho of the option. + * @param riskFreeRate The risk free rate (the drift of the forward, repoRate - dividendYield = r-q). + * @param discountRate The discount rate. + * @param sigma the square root of the initial instantaneous variance (\( V_0 = sigma^2 \)) + * @param theta the long run mean of the volatility. + * @param kappa the speed of mean reversion. + * @param xi the volatility of variance. + * @param rho correlation between the two Brownian motions + * @param optionMaturity the maturity of the option + * @param optionStrike the strike of the option. + * @return The rho of the option */ public static double hestonOptionRho( final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double discountRate, + final double sigma, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, - final double optionStrike) { - - final HestonGreek whatToCompute = HestonGreek.RHO; + final double optionStrike) + { + HestonGreek whatToCompute = HestonGreek.RHO; - return hestonGreekCalculator( - initialStockValue, + return hestonGreekCalculator(initialStockValue, riskFreeRate, - dividendYield, - kappa, - theta, - sigma, - v0, + discountRate, + sigma, + theta, + kappa, + xi, rho, optionMaturity, optionStrike, whatToCompute, - NUMBER_OF_POINTS, - GRID_SPACING); + numberOfPoints, + gridSpacing); } /** - * Calculates the vega1 of a call option under a Heston model. - * + * Calculates the vega of a call option under a Heston model, that is + * \( d/d \sigma \), where \( \sigma = \sqrt{v_0} \) is the square root of the initial variance of the model, + * i.e., the initial volatility. + * * @param initialStockValue Initital value of the stock. - * @param riskFreeRate The risk free rate. - * @param dividendYield The dividend yield. - * @param kappa The speed of mean reversion. - * @param theta The long run mean of the volatility. - * @param sigma The volatility of variance. - * @param v0 The initial instantaneous variance. - * @param rho Correlation between the two Brownian motions. - * @param optionMaturity The maturity of the option. - * @param optionStrike The strike of the option. - * @return The vega1 of the option. + * @param riskFreeRate The risk free rate (the drift of the forward, repoRate - dividendYield = r-q). + * @param discountRate The discount rate. + * @param sigma the square root of the initial instantaneous variance (\( V_0 = sigma^2 \)) + * @param theta the long run mean of the volatility. + * @param kappa the speed of mean reversion. + * @param xi the volatility of variance. + * @param rho correlation between the two Brownian motions + * @param optionMaturity the maturity of the option + * @param optionStrike the strike of the option. + * @return The vega of the option */ - public static double hestonOptionVega1( + public static double hestonOptionVega( final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double discountRate, + final double sigma, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, - final double optionStrike) { + final double optionStrike) + { + return hestonOptionVega1(initialStockValue, riskFreeRate, discountRate, sigma, theta, kappa, xi, rho, optionMaturity, optionStrike) * 2 * sigma; + } - final HestonGreek whatToCompute = HestonGreek.VEGA1; + /** + * Calculates the vega1 of a call option under a Heston model, that is + * \( d/d v_0 \), where \( v_0 \) is the initial variance of the model. + * + * @param initialStockValue Initital value of the stock. + * @param riskFreeRate The risk free rate (the drift of the forward, repoRate - dividendYield = r-q). + * @param discountRate The discount rate. + * @param sigma the square root of the initial instantaneous variance (\( V_0 = sigma^2 \)) + * @param theta the long run mean of the volatility. + * @param kappa the speed of mean reversion. + * @param xi the volatility of variance. + * @param rho correlation between the two Brownian motions + * @param optionMaturity the maturity of the option + * @param optionStrike the strike of the option. + * @return The vega1 of the option + */ + public static double hestonOptionVega1( + final double initialStockValue, + final double riskFreeRate, + final double discountRate, + final double sigma, + final double theta, + final double kappa, + final double xi, + final double rho, + final double optionMaturity, + final double optionStrike) + { + HestonGreek whatToCompute = HestonGreek.VEGA1; - return hestonGreekCalculator( - initialStockValue, + return hestonGreekCalculator(initialStockValue, riskFreeRate, - dividendYield, - kappa, - theta, - sigma, - v0, + discountRate, + sigma, + theta, + kappa, + xi, rho, optionMaturity, optionStrike, whatToCompute, - NUMBER_OF_POINTS, - GRID_SPACING); + numberOfPoints, + gridSpacing); } /** - * Calculates the vanna of a call option under a Heston model. - * + * Calculates the vanna of a call option under a Heston model + * * @param initialStockValue Initital value of the stock. - * @param riskFreeRate The risk free rate. - * @param dividendYield The dividend yield. - * @param kappa The speed of mean reversion. - * @param theta The long run mean of the volatility. - * @param sigma The volatility of variance. - * @param v0 The initial instantaneous variance. - * @param rho Correlation between the two Brownian motions. - * @param optionMaturity The maturity of the option. - * @param optionStrike The strike of the option. - * @return The vanna of the option. + * @param riskFreeRate The risk free rate (the drift of the forward, repoRate - dividendYield = r-q). + * @param discountRate The discount rate. + * @param sigma the square root of the initial instantaneous variance (\( V_0 = sigma^2 \)) + * @param theta the long run mean of the volatility. + * @param kappa the speed of mean reversion. + * @param xi the volatility of variance. + * @param rho correlation between the two Brownian motions + * @param optionMaturity the maturity of the option + * @param optionStrike the strike of the option. + * @return The vanna of the option */ public static double hestonOptionVanna( final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double discountRate, + final double sigma, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, - final double optionStrike) { + final double optionStrike) + { + HestonGreek whatToCompute = HestonGreek.VANNA; - final HestonGreek whatToCompute = HestonGreek.VANNA; - - return hestonGreekCalculator( - initialStockValue, + return hestonGreekCalculator(initialStockValue, riskFreeRate, - dividendYield, - kappa, - theta, - sigma, - v0, + discountRate, + sigma, + theta, + kappa, + xi, rho, optionMaturity, optionStrike, whatToCompute, - NUMBER_OF_POINTS, - GRID_SPACING); + numberOfPoints, + gridSpacing); } /** - * Calculates the volga of a call option under a Heston model. - * + * Calculates the volga of a call option under a Heston model + * * @param initialStockValue Initital value of the stock. - * @param riskFreeRate The risk free rate. - * @param dividendYield The dividend yield. - * @param kappa The speed of mean reversion. - * @param theta The long run mean of the volatility. - * @param sigma The volatility of variance. - * @param v0 The initial instantaneous variance. - * @param rho Correlation between the two Brownian motions. - * @param optionMaturity The maturity of the option. - * @param optionStrike The strike of the option. - * @return The volga of the option. + * @param riskFreeRate The risk free rate (the drift of the forward, repoRate - dividendYield = r-q). + * @param discountRate The discount rate. + * @param sigma the square root of the initial instantaneous variance (\( V_0 = sigma^2 \)) + * @param theta the long run mean of the volatility. + * @param kappa the speed of mean reversion. + * @param xi the volatility of variance. + * @param rho correlation between the two Brownian motions + * @param optionMaturity the maturity of the option + * @param optionStrike the strike of the option. + * @return The volga of the option */ public static double hestonOptionVolga( final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double discountRate, + final double sigma, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, - final double optionStrike) { - - final HestonGreek whatToCompute = HestonGreek.VOLGA; + final double optionStrike) + { + HestonGreek whatToCompute = HestonGreek.VOLGA; - return hestonGreekCalculator( - initialStockValue, + return hestonGreekCalculator(initialStockValue, riskFreeRate, - dividendYield, - kappa, - theta, - sigma, - v0, + discountRate, + sigma, + theta, + kappa, + xi, rho, optionMaturity, optionStrike, whatToCompute, - NUMBER_OF_POINTS, - GRID_SPACING); + numberOfPoints, + gridSpacing); } + + /** + * Computes the gradient of the (discounted) characteristic function of the Heston model. + * More sensitivies can be added by introducing more cases in the switch statement. + * + * @param zeta The argument of the characteristic function gradient + * @param initialStockValue Initital value of the stock. + * @param riskFreeRate The risk free rate (the drift of the forward, repoRate - dividendYield = r-q). + * @param discountRate The discount rate. + * @param sigma the square root of the initial instantaneous variance (\( V_0 = sigma^2 \)) + * @param theta the long run mean of the volatility. + * @param kappa the speed of mean reversion. + * @param xi the volatility of variance. + * @param rho correlation between the two Brownian motions + * @param optionMaturity the maturity of the option + * @param optionStrike the strike of the option. + * @param whichGreek + * @param numberOfPoints + * @param gridSpacing + * @return + */ private static Complex hestonCharacteristicFunctionGradient( final Complex zeta, final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double discountRate, + final double sigma, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, final double optionStrike, @@ -361,115 +435,151 @@ private static Complex hestonCharacteristicFunctionGradient( final int numberOfPoints, final double gridSpacing) { - final double lambda = 2 * Math.PI / (numberOfPoints * gridSpacing); + final double lambda = 2*Math.PI/(numberOfPoints*gridSpacing); + final double v0 = sigma*sigma; - final double x = Math.log(initialStockValue); - final double a = kappa * theta; + double x = Math.log(initialStockValue); + double a = kappa * theta; - Complex b; - Complex u; - Complex d; - Complex g; - Complex c; - Complex D; - Complex G; - Complex C; + Complex b, u, d, g, c, D, G, C; + // Initialize u = new Complex(-0.5, 0.0); b = new Complex(kappa + lambda, 0.0); - final Complex term1 = (Complex.I.multiply(rho * sigma).multiply(zeta)).subtract(b); - final Complex powPart = term1.pow(2.0); - final Complex term2 = new Complex(sigma * sigma, 0.0) + // Compute d + Complex term1 = (Complex.I.multiply(rho * xi).multiply(zeta)).subtract(b); // (ρσiφ - b) + Complex powPart = term1.pow(2.0); // (...) + Complex term2 = new Complex(xi * xi, 0.0) .multiply((u.multiply(2.0).multiply(Complex.I).multiply(zeta)) - .subtract(zeta.multiply(zeta))); + .subtract(zeta.multiply(zeta))); // σ²(2u i φ - φ²) d = powPart.subtract(term2).sqrt(); - final Complex numerator = (b.subtract(Complex.I.multiply(rho * sigma).multiply(zeta))).add(d); - final Complex denominator = (b.subtract(Complex.I.multiply(rho * sigma).multiply(zeta))).subtract(d); + // Compute g + Complex numerator = (b.subtract(Complex.I.multiply(rho * xi).multiply(zeta))).add(d); + Complex denominator = (b.subtract(Complex.I.multiply(rho * xi).multiply(zeta))).subtract(d); g = numerator.divide(denominator); + // Initialize remaining + c = Complex.ZERO; + D = Complex.ZERO; + G = Complex.ZERO; + C = Complex.ZERO; + + // "Little Heston Trap" formulation + // c = 1.0 / g c = Complex.ONE.divide(g); - final Complex exp_dT = d.multiply(-optionMaturity).exp(); - final Complex bMinusRhoSigmaIphiMinusD = b.subtract(Complex.I.multiply(rho * sigma).multiply(zeta)).subtract(d); + // Precompute some recurring parts + Complex exp_dT = d.multiply(-optionMaturity).exp(); // exp(-d*T) + Complex bMinusRhoSigmaIphiMinusD = b.subtract(Complex.I.multiply(rho * xi).multiply(zeta)).subtract(d); + // D = (b - rho*sigma*i*phi - d) / (sigma^2) * ((1 - exp(-d*T)) / (1 - c * exp(-d*T))) D = bMinusRhoSigmaIphiMinusD - .divide(sigma * sigma) + .divide(xi * xi) .multiply( (Complex.ONE.subtract(exp_dT)) - .divide(Complex.ONE.subtract(c.multiply(exp_dT)))); + .divide(Complex.ONE.subtract(c.multiply(exp_dT))) + ); + // G = (1 - c * exp(-d*T)) / (1 - c) G = (Complex.ONE.subtract(c.multiply(exp_dT))) .divide(Complex.ONE.subtract(c)); - final Complex firstTerm = Complex.I.multiply(zeta).multiply((riskFreeRate - dividendYield) * optionMaturity); - final Complex secondTerm = new Complex(a / (sigma * sigma), 0.0) + // C = (r - q) * i * phi * T + a / sigma^2 * ((b - rho*sigma*i*phi - d) * T - 2.0 * Complex.Log(G)) + Complex firstTerm = Complex.I.multiply(zeta).multiply(riskFreeRate * optionMaturity); + Complex secondTerm = new Complex(a / (xi * xi), 0.0) .multiply( (bMinusRhoSigmaIphiMinusD.multiply(optionMaturity)) - .subtract(Complex.valueOf(2.0).multiply(G.log()))); + .subtract(Complex.valueOf(2.0).multiply(G.log())) + ); C = firstTerm.add(secondTerm); - final Complex f = (C.add(D.multiply(v0)).add(Complex.I.multiply(zeta).multiply(x))).exp(); + // The characteristic function (discounted) + Complex f = (C.add(D.multiply(v0)).add(Complex.I.multiply(zeta).multiply(x))).add(-discountRate * optionMaturity).exp(); - switch(whichGreek) { + // Return depending on the requested Greek + switch (whichGreek) { case DELTA: - return f.multiply(Complex.I).multiply(zeta).divide(initialStockValue); + return f.multiply(Complex.I).multiply(zeta).divide(initialStockValue); case GAMMA: return Complex.I.multiply(zeta) .multiply(f) .multiply(Complex.I.multiply(zeta).subtract(Complex.ONE)) - .divide(initialStockValue * initialStockValue); + .divide(initialStockValue * initialStockValue); case THETA: - final Complex numerator_dDdT = d.multiply(exp_dT) - .multiply(b.subtract(Complex.I.multiply(rho * sigma).multiply(zeta)).add(d)) - .multiply(g.subtract(Complex.ONE)); - final Complex denominator_dDdT = Complex.valueOf(sigma * sigma) + Complex numerator_dDdT = d.multiply(exp_dT) + .multiply(b.subtract(Complex.I.multiply(rho * xi).multiply(zeta)).add(d)) + .multiply(g.subtract(Complex.ONE)); + Complex denominator_dDdT = Complex.valueOf(xi * xi) .multiply(Complex.ONE.subtract(g.multiply(exp_dT)).pow(2.0)); - final Complex dDdT = numerator_dDdT.divide(denominator_dDdT); + Complex dDdT = numerator_dDdT.divide(denominator_dDdT); - final Complex innerTerm = (b.subtract(Complex.I.multiply(rho * sigma).multiply(zeta)).add(d)) + Complex innerTerm = (b.subtract(Complex.I.multiply(rho * xi).multiply(zeta)).add(d)) .add( Complex.valueOf(2.0) - .multiply(g) - .multiply(d) - .multiply(exp_dT) - .divide(Complex.ONE.subtract(g.multiply(exp_dT)))); + .multiply(g) + .multiply(d) + .multiply(exp_dT) + .divide(Complex.ONE.subtract(g.multiply(exp_dT))) + ); - final Complex dCdT = Complex.I.multiply(zeta).multiply(riskFreeRate) + Complex dCdT = Complex.I.multiply(zeta).multiply(riskFreeRate) .add( - Complex.valueOf(kappa * theta / (sigma * sigma)) - .multiply(innerTerm)); + Complex.valueOf(kappa * theta / (xi * xi)) + .multiply(innerTerm) + ); return Complex.valueOf(riskFreeRate) .multiply(f) - .subtract(f.multiply(dCdT.add(dDdT.multiply(v0)))); + .subtract( + f.multiply(dCdT.add(dDdT.multiply(v0))) + ); case RHO: - final Complex dCdr = Complex.I.multiply(zeta).multiply(optionMaturity); - return f.multiply(dCdr).subtract(f.multiply(optionMaturity)); + Complex dCdr = Complex.I.multiply(zeta).multiply(optionMaturity); + return f.multiply(dCdr).subtract(f.multiply(optionMaturity)); case VEGA1: - return f.multiply(D); + return f.multiply(D); case VANNA: - return f.multiply(Complex.I).multiply(zeta).multiply(D).divide(initialStockValue); + return f.multiply(Complex.I).multiply(zeta).multiply(D).divide(initialStockValue); case VOLGA: return Complex.valueOf(2.0) .multiply(D) .multiply(f) - .multiply(D.multiply(2.0 * v0).add(Complex.ONE)); + .multiply( + D.multiply(2.0 * v0).add(Complex.ONE) + ); default: throw new IllegalArgumentException("Unknown Greek: " + whichGreek); } } - private static double hestonGreekCalculator( - final double initialStockValue, + /** + * Service method that performs the Fourier inversion according to the FFT algorithm of Carr and Madan + * + * @param initialStockValue Initital value of the stock. + * @param riskFreeRate The risk free rate (the drift of the forward, repoRate - dividendYield = r-q). + * @param discountRate The discount rate. + * @param sigma the square root of the initial instantaneous variance (\( V_0 = sigma^2 \)) + * @param theta the long run mean of the volatility. + * @param kappa the speed of mean reversion. + * @param xi the volatility of variance. + * @param rho correlation between the two Brownian motions + * @param optionMaturity the maturity of the option + * @param optionStrike the strike of the option. + * @param whichGreek + * @param numberOfPoints + * @param gridSpacing + * @return the requested calculation + */ + private static double hestonGreekCalculator(final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, - final double sigma, - final double v0, + final double discountRate, + final double sigma, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, final double optionStrike, @@ -479,25 +589,23 @@ private static double hestonGreekCalculator( final double lineOfIntegration = 1.2; - final double lambda = 2 * Math.PI / (numberOfPoints * gridSpacing); - final double upperBound = (numberOfPoints * lambda) / 2.0; + final double lambda = 2*Math.PI/(numberOfPoints*gridSpacing); //Equation 23 Carr and Madan + final double upperBound = (numberOfPoints * lambda)/2.0; //Equation 20 Carr and Madan + + final double v0 = sigma*sigma; final Complex[] integrandEvaluations = new Complex[numberOfPoints]; - for(int i = 0; i < numberOfPoints; i++) { + for(int i = 0; i strikeToValue = new Function(){ - final Function strikeToValue = new Function() { @Override public Double apply(final Double t) { - return interpolation.getValue(t); + return interpolation.getValue(t); } + }; return strikeToValue.apply(optionStrike); From fd72e7d9efcc6e85bee5721b031057d86ade6d4c Mon Sep 17 00:00:00 2001 From: Alessandro Gnoatto Date: Thu, 19 Feb 2026 10:53:26 +0100 Subject: [PATCH 7/7] Change class names to enforce library naming conventions. --- ...hDependent.java => AmericanNonPathDependent.java} | 5 ++--- ...hDependent.java => EuropeanNonPathDependent.java} | 5 ++--- .../tree/EuropeanAndAmericanOptionPricingTest.java | 12 ++++++------ 3 files changed, 10 insertions(+), 12 deletions(-) rename src/main/java/net/finmath/tree/assetderivativevaluation/products/{UsNonPathDependent.java => AmericanNonPathDependent.java} (93%) rename src/main/java/net/finmath/tree/assetderivativevaluation/products/{EuNonPathDependent.java => EuropeanNonPathDependent.java} (91%) diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/products/UsNonPathDependent.java b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AmericanNonPathDependent.java similarity index 93% rename from src/main/java/net/finmath/tree/assetderivativevaluation/products/UsNonPathDependent.java rename to src/main/java/net/finmath/tree/assetderivativevaluation/products/AmericanNonPathDependent.java index 2916fa0860..08e2d7cefb 100644 --- a/src/main/java/net/finmath/tree/assetderivativevaluation/products/UsNonPathDependent.java +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AmericanNonPathDependent.java @@ -15,8 +15,7 @@ * Immediate exercise value: f(S_k). * The option value is the pointwise maximum of the two (early exercise rule). */ - -public class UsNonPathDependent extends AbstractNonPathDependentProduct { +public class AmericanNonPathDependent extends AbstractNonPathDependentProduct { /** * Creates an American option with given maturity and payoff. @@ -24,7 +23,7 @@ public class UsNonPathDependent extends AbstractNonPathDependentProduct { * @param maturity Contract maturity (model time units). * @param payOffFunction Payoff function f(S) (e.g.,s -> Math.max(K - s, 0.0) for a put). */ - public UsNonPathDependent(double maturity, DoubleUnaryOperator payOffFunction) { + public AmericanNonPathDependent(double maturity, DoubleUnaryOperator payOffFunction) { super(maturity, payOffFunction); } diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/products/EuNonPathDependent.java b/src/main/java/net/finmath/tree/assetderivativevaluation/products/EuropeanNonPathDependent.java similarity index 91% rename from src/main/java/net/finmath/tree/assetderivativevaluation/products/EuNonPathDependent.java rename to src/main/java/net/finmath/tree/assetderivativevaluation/products/EuropeanNonPathDependent.java index 859c7f5787..69cb014d12 100644 --- a/src/main/java/net/finmath/tree/assetderivativevaluation/products/EuNonPathDependent.java +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/EuropeanNonPathDependent.java @@ -11,8 +11,7 @@ * by backward induction using the model’s discounted conditional expectation. * pricing uses evaluationTime = 0.0 => k0 = 0. */ - -public class EuNonPathDependent extends AbstractNonPathDependentProduct { +public class EuropeanNonPathDependent extends AbstractNonPathDependentProduct { /** * Creates a European option with given maturity and payoff. @@ -21,7 +20,7 @@ public class EuNonPathDependent extends AbstractNonPathDependentProduct { * @param payOffFunction Payoff function f(S) applied at T * (e.g., s -> Math.max(s-K,0.0) for a call). */ - public EuNonPathDependent(double maturity, DoubleUnaryOperator payOffFunction){ + public EuropeanNonPathDependent(double maturity, DoubleUnaryOperator payOffFunction){ super(maturity,payOffFunction); } diff --git a/src/test/java/net/finmath/tree/EuropeanAndAmericanOptionPricingTest.java b/src/test/java/net/finmath/tree/EuropeanAndAmericanOptionPricingTest.java index fa45459538..5baf489d8d 100644 --- a/src/test/java/net/finmath/tree/EuropeanAndAmericanOptionPricingTest.java +++ b/src/test/java/net/finmath/tree/EuropeanAndAmericanOptionPricingTest.java @@ -7,8 +7,8 @@ import net.finmath.tree.assetderivativevaluation.models.BoyleTrinomial; import net.finmath.tree.assetderivativevaluation.models.CoxRossRubinsteinModel; import net.finmath.tree.assetderivativevaluation.models.JarrowRuddModel; -import net.finmath.tree.assetderivativevaluation.products.EuNonPathDependent; -import net.finmath.tree.assetderivativevaluation.products.UsNonPathDependent; +import net.finmath.tree.assetderivativevaluation.products.EuropeanNonPathDependent; +import net.finmath.tree.assetderivativevaluation.products.AmericanNonPathDependent; public class EuropeanAndAmericanOptionPricingTest { @@ -29,8 +29,8 @@ public void testNonPathDep_Call() { - EuNonPathDependent euCall = new EuNonPathDependent(maturity, s -> Math.max(s - strike, 0.0)); - UsNonPathDependent usCall = new UsNonPathDependent(maturity, s -> Math.max(s - strike, 0.0)); + EuropeanNonPathDependent euCall = new EuropeanNonPathDependent(maturity, s -> Math.max(s - strike, 0.0)); + AmericanNonPathDependent usCall = new AmericanNonPathDependent(maturity, s -> Math.max(s - strike, 0.0)); double euCRR = euCall.getValue(crr); @@ -75,8 +75,8 @@ public void testNonPathDep_Put() { - EuNonPathDependent euPut = new EuNonPathDependent(maturity, s -> Math.max(strike - s, 0.0)); - UsNonPathDependent usPut = new UsNonPathDependent(maturity, s -> Math.max(strike - s, 0.0)); + EuropeanNonPathDependent euPut = new EuropeanNonPathDependent(maturity, s -> Math.max(strike - s, 0.0)); + AmericanNonPathDependent usPut = new AmericanNonPathDependent(maturity, s -> Math.max(strike - s, 0.0)); double euCRR = euPut.getValue(0.0, crr).getAverage();