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..20579077fe 100644 --- a/src/main/java/net/finmath/functions/HestonModel.java +++ b/src/main/java/net/finmath/functions/HestonModel.java @@ -14,6 +14,33 @@ /** * 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 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 */ @@ -27,18 +54,18 @@ public class HestonModel { 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 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 sigma the volatility of variance. - * @param v0 the initial instantaneous variance + * @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. @@ -47,11 +74,11 @@ enum HestonGreek {DELTA, GAMMA, THETA, RHO, VEGA1, VANNA, VOLGA}; public static double hestonOptionDelta( final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, + final double discountRate, final double sigma, - final double v0, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, final double optionStrike) @@ -61,11 +88,11 @@ public static double hestonOptionDelta( return hestonGreekCalculator(initialStockValue, riskFreeRate, - dividendYield, - kappa, - theta, + discountRate, sigma, - v0, + theta, + kappa, + xi, rho, optionMaturity, optionStrike, @@ -78,12 +105,12 @@ public static double hestonOptionDelta( * 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 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 sigma the volatility of variance. - * @param v0 the initial instantaneous variance + * @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. @@ -92,11 +119,11 @@ public static double hestonOptionDelta( public static double hestonOptionGamma( final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, + final double discountRate, final double sigma, - final double v0, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, final double optionStrike) @@ -105,11 +132,11 @@ public static double hestonOptionGamma( return hestonGreekCalculator(initialStockValue, riskFreeRate, - dividendYield, - kappa, - theta, + discountRate, sigma, - v0, + theta, + kappa, + xi, rho, optionMaturity, optionStrike, @@ -122,12 +149,12 @@ public static double hestonOptionGamma( * 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 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 sigma the volatility of variance. - * @param v0 the initial instantaneous variance + * @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. @@ -136,11 +163,11 @@ public static double hestonOptionGamma( public static double hestonOptionTheta( final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, + final double discountRate, final double sigma, - final double v0, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, final double optionStrike) @@ -149,11 +176,11 @@ public static double hestonOptionTheta( return hestonGreekCalculator(initialStockValue, riskFreeRate, - dividendYield, - kappa, - theta, + discountRate, sigma, - v0, + theta, + kappa, + xi, rho, optionMaturity, optionStrike, @@ -166,12 +193,12 @@ public static double hestonOptionTheta( * 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 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 sigma the volatility of variance. - * @param v0 the initial instantaneous variance + * @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. @@ -180,11 +207,11 @@ public static double hestonOptionTheta( public static double hestonOptionRho( final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, + final double discountRate, final double sigma, - final double v0, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, final double optionStrike) @@ -193,11 +220,11 @@ public static double hestonOptionRho( return hestonGreekCalculator(initialStockValue, riskFreeRate, - dividendYield, - kappa, - theta, + discountRate, sigma, - v0, + theta, + kappa, + xi, rho, optionMaturity, optionStrike, @@ -207,15 +234,48 @@ public static double hestonOptionRho( } /** - * 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 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 hestonOptionVega( + 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) + { + return hestonOptionVega1(initialStockValue, riskFreeRate, discountRate, sigma, theta, kappa, xi, rho, optionMaturity, optionStrike) * 2 * sigma; + } + + /** + * 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 sigma the volatility of variance. - * @param v0 the initial instantaneous variance + * @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. @@ -224,11 +284,11 @@ public static double hestonOptionRho( public static double hestonOptionVega1( final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, + final double discountRate, final double sigma, - final double v0, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, final double optionStrike) @@ -237,11 +297,11 @@ public static double hestonOptionVega1( return hestonGreekCalculator(initialStockValue, riskFreeRate, - dividendYield, - kappa, - theta, + discountRate, sigma, - v0, + theta, + kappa, + xi, rho, optionMaturity, optionStrike, @@ -254,12 +314,12 @@ public static double hestonOptionVega1( * 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 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 sigma the volatility of variance. - * @param v0 the initial instantaneous variance + * @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. @@ -268,11 +328,11 @@ public static double hestonOptionVega1( public static double hestonOptionVanna( final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, + final double discountRate, final double sigma, - final double v0, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, final double optionStrike) @@ -281,11 +341,11 @@ public static double hestonOptionVanna( return hestonGreekCalculator(initialStockValue, riskFreeRate, - dividendYield, - kappa, - theta, + discountRate, sigma, - v0, + theta, + kappa, + xi, rho, optionMaturity, optionStrike, @@ -298,12 +358,12 @@ public static double hestonOptionVanna( * 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 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 sigma the volatility of variance. - * @param v0 the initial instantaneous variance + * @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. @@ -312,11 +372,11 @@ public static double hestonOptionVanna( public static double hestonOptionVolga( final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, + final double discountRate, final double sigma, - final double v0, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, final double optionStrike) @@ -325,11 +385,11 @@ public static double hestonOptionVolga( return hestonGreekCalculator(initialStockValue, riskFreeRate, - dividendYield, - kappa, - theta, + discountRate, sigma, - v0, + theta, + kappa, + xi, rho, optionMaturity, optionStrike, @@ -343,17 +403,17 @@ public static double hestonOptionVolga( * 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 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 @@ -363,11 +423,11 @@ private static Complex hestonCharacteristicFunctionGradient( final Complex zeta, final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, + final double discountRate, final double sigma, - final double v0, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, final double optionStrike, @@ -376,6 +436,7 @@ private static Complex hestonCharacteristicFunctionGradient( final double gridSpacing) { final double lambda = 2*Math.PI/(numberOfPoints*gridSpacing); + final double v0 = sigma*sigma; double x = Math.log(initialStockValue); double a = kappa * theta; @@ -387,16 +448,16 @@ private static Complex hestonCharacteristicFunctionGradient( b = new Complex(kappa + lambda, 0.0); // Compute d - Complex term1 = (Complex.I.multiply(rho * sigma).multiply(zeta)).subtract(b); // (ρσiφ - b) + Complex term1 = (Complex.I.multiply(rho * xi).multiply(zeta)).subtract(b); // (ρσiφ - b) Complex powPart = term1.pow(2.0); // (...) - Complex term2 = new Complex(sigma * sigma, 0.0) + Complex term2 = new Complex(xi * xi, 0.0) .multiply((u.multiply(2.0).multiply(Complex.I).multiply(zeta)) .subtract(zeta.multiply(zeta))); // σ²(2u i φ - φ²) 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); + 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 @@ -411,11 +472,11 @@ private static Complex hestonCharacteristicFunctionGradient( // 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); + 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))) @@ -426,8 +487,8 @@ private static Complex hestonCharacteristicFunctionGradient( .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) + 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())) @@ -435,9 +496,8 @@ private static Complex hestonCharacteristicFunctionGradient( C = firstTerm.add(secondTerm); - // The characteristic function - 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(); // Return depending on the requested Greek switch (whichGreek) { @@ -450,13 +510,13 @@ private static Complex hestonCharacteristicFunctionGradient( .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(b.subtract(Complex.I.multiply(rho * xi).multiply(zeta)).add(d)) .multiply(g.subtract(Complex.ONE)); - Complex denominator_dDdT = Complex.valueOf(sigma * sigma) + Complex denominator_dDdT = Complex.valueOf(xi * xi) .multiply(Complex.ONE.subtract(g.multiply(exp_dT)).pow(2.0)); Complex dDdT = numerator_dDdT.divide(denominator_dDdT); - 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) @@ -467,7 +527,7 @@ private static Complex hestonCharacteristicFunctionGradient( Complex dCdT = Complex.I.multiply(zeta).multiply(riskFreeRate) .add( - Complex.valueOf(kappa * theta / (sigma * sigma)) + Complex.valueOf(kappa * theta / (xi * xi)) .multiply(innerTerm) ); @@ -499,12 +559,12 @@ private static Complex hestonCharacteristicFunctionGradient( * 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 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 sigma the volatility of variance. - * @param v0 the initial instantaneous variance + * @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. @@ -515,11 +575,11 @@ private static Complex hestonCharacteristicFunctionGradient( */ private static double hestonGreekCalculator(final double initialStockValue, final double riskFreeRate, - final double dividendYield, - final double kappa, - final double theta, + final double discountRate, final double sigma, - final double v0, + final double theta, + final double kappa, + final double xi, final double rho, final double optionMaturity, final double optionStrike, @@ -532,6 +592,8 @@ private static double hestonGreekCalculator(final double initialStockValue, 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 + *
  • 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 + */ +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(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. + * + * @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(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 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 evaluationTime as a random variable on the tree. + * @throws IllegalArgumentException If the inputs are invalid (see validate(double, TreeModel)). + */ + 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). + * 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(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() { + 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..434dcf82b0 --- /dev/null +++ b/src/main/java/net/finmath/tree/OneDimensionalRiskFactorTreeModel.java @@ -0,0 +1,93 @@ +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 + * @author Alessandro Gnoatto + */ +public abstract class OneDimensionalRiskFactorTreeModel implements TreeModel { + + /** Cache of level spot S_k. */ + private final List riskFactorLevels = new ArrayList<>(); + + /** + * 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 {@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 {@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 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) { + final int next = riskFactorLevels.size(); + riskFactorLevels.add(buildSpotLevel(next)); + } + return riskFactorLevels.get(k); + } + + @Override + 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; + } + + 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(final RandomVariable vNext, final int k) { + final RandomVariable vK = conditionalExpectation(vNext, k); // hook + return new RandomVariableFromDoubleArray(k * getTimeStep(), vK.getRealizations()); + } + + @Override + 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/TreeModel.java b/src/main/java/net/finmath/tree/TreeModel.java new file mode 100644 index 0000000000..4ef01499a6 --- /dev/null +++ b/src/main/java/net/finmath/tree/TreeModel.java @@ -0,0 +1,115 @@ +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 + * @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. + * 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); + + /** + * 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."); + } + + /** + * 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/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..da547a3898 --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/AbstractRecombiningTreeModel.java @@ -0,0 +1,101 @@ +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 */ + @Override + public double getOneStepDiscountFactor(int timeIndex) { + 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..44dc629156 --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/OneDimensionalAssetTreeModel.java @@ -0,0 +1,31 @@ +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. + */ + double getInitialPrice(); + + /** + * The risk free rate. + * + * @return The risk free rate. + */ + double getRiskFreeRate(); + + /** + * Returns the volatility. + * + * @return The volatility. + */ + 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..beee80ce76 --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/models/BoyleTrinomial.java @@ -0,0 +1,162 @@ +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 S0 * u^(k-i) * d^i. + * @param k level of depth */ + @Override + protected RandomVariable buildSpotLevel(int k) { + double[] s = new double[statesAt(k)]; + double s0 = getInitialPrice(); + int down = 0; + int up = 0; + for (int i = 0; i <= k; i++) { + up = k - i; + down = i; + s[i] = s0 * Math.pow(u, up) * Math.pow(d, down); + } + return new RandomVariableFromDoubleArray(k*getTimeStep(),s); + } + + /** + * Discounted conditional expectation : V_k[i] = df() * ( q * V_{k+1}[i] + (1-q) * V_{k+1}[i+1] ). + */ + protected RandomVariable conditionalExpectation(RandomVariable vNext, int k) { + double[] next = vNext.getRealizations(); + double[] vK = new double[statesAt(k)]; + double disc = getOneStepDiscountFactor(k); + + for (int i = 0; i <= k; i++) { + vK[i] = disc * (q * next[i] + (1.0 - q) * next[i + 1]); + } + 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; } + public double getD() { return d; } + public double getQ() { return q; } +} + diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/models/JarrowRuddModel.java b/src/main/java/net/finmath/tree/assetderivativevaluation/models/JarrowRuddModel.java new file mode 100644 index 0000000000..ec34f25ff3 --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/models/JarrowRuddModel.java @@ -0,0 +1,125 @@ +package net.finmath.tree.assetderivativevaluation.models; + +import net.finmath.montecarlo.RandomVariableFromDoubleArray; +import net.finmath.stochastic.RandomVariable; +import net.finmath.tree.assetderivativevaluation.AbstractRecombiningTreeModel; + +/** + * This class represents the approximation of a Black-Scholes model via the Jarrow-Rudd model. + * It extends ApproximatingBinomialModel. The only method that is implemented here computes the values + * of the up and down movements of the Binomial model. + * + * @author Andrea Mazzon + * + */ + +public class JarrowRuddModel extends AbstractRecombiningTreeModel { + + /** + * Constructs a Jarrow-Rudd model using an equidistant time step. + * + * @param spotPrice the initial 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<... S0 * u^(k-i) * d^i. + * @param k level of depth */ + @Override + protected RandomVariable buildSpotLevel(int k) { + double s0 = getInitialPrice(); + double dt = getTimeStep(); + double sigma = getVolatility(); + double r = getRiskFreeRate(); + + double muStar = (r - 0.5*sigma*sigma) * dt; + double nu = sigma * Math.sqrt(dt); + + double u = Math.exp(muStar + nu); + double d = Math.exp(muStar - nu); + + double[] level = new double[k + 1]; + for (int i = 0; i <= k; i++) { + int ups = k - i; + int downs = i; + level[i] = s0 * Math.pow(u, ups) * Math.pow(d, downs); + } + return new RandomVariableFromDoubleArray(k * dt, level); + } + + /** + * Discounted conditional expectation : V_k[i] = df() * ( q * V_{k+1}[i] + (1-q) * V_{k+1}[i+1] ). + */ + @Override + protected RandomVariable conditionalExpectation(RandomVariable vNext, int k) { + final double r = getRiskFreeRate(); + final double dt = getTimeStep(); + + final double disc = Math.exp(-r * dt); + final double p = 0.5; + + final double[] next = vNext.getRealizations(); + final int expected = k + 2; + if (next == null || next.length != expected) { + throw new IllegalArgumentException("vNext length " + (next == null ? "null" : next.length) + + " != expected " + expected + " for time index k=" + k); + } + + final int nHere = statesAt(k); + final double[] res = new double[nHere]; + for (int i = 0; i < nHere; i++) { + res[i] = disc * (p * next[i] + (1.0 - p) * next[i + 1]); + } + + 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/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/AbstractNonPathDependentProduct.java b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractNonPathDependentProduct.java new file mode 100644 index 0000000000..1a8ecd1d27 --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractNonPathDependentProduct.java @@ -0,0 +1,67 @@ +package net.finmath.tree.assetderivativevaluation.products; + +import net.finmath.stochastic.RandomVariable; +import net.finmath.tree.AbstractTreeProduct; +import net.finmath.tree.TreeModel; + +import java.util.function.DoubleUnaryOperator; + + +/** + * Base class for non–path-dependent options priced on a TreeModel. + * The payoff is fully specified by a DoubleUnaryOperator acting on the spot + * 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). */ + private final DoubleUnaryOperator payOffFunction; + + /** + * Creates a non–path-dependent option with a given maturity and payoff. + * + * @param maturity + * Contract maturity (model time units), typically T. + * @param payOffFunction + * Payoff function f(S) (e.g. call/put); must not be null. + */ + public AbstractNonPathDependentProduct(double maturity, DoubleUnaryOperator payOffFunction ){ + super(maturity); + this.payOffFunction = payOffFunction; + } + + /** + * Returns the payoff function f(S). + * + * @return The payoff function. + */ + protected DoubleUnaryOperator getPayOffFunction(){ + return payOffFunction; + } + + /** + * Converts a model time to the corresponding lattice time index by rounding + * to the nearest step: index = round(time / model.getTimeStep()). + * This helper is useful to map maturity (or any evaluation time aligned to the + * model discretization) to an integer level of the tree. + * @param maturity A model time (e.g. T) to map to a lattice index. + * @param model The tree model providing the time step. + * @return The nearest lattice index for the given time. + */ + protected final int timeToIndex(double maturity, TreeModel model){ + return (int)Math.round(maturity/model.getTimeStep()); + } + + /** + * Computes the vector of option values on the lattice (via backward induction). + * @param evaluationTime The time at which valuation is performed (e.g. 0.0). + * @param model The tree model to use; must not be null. + * @return The array of value random variables per time level. + */ + @Override + public abstract RandomVariable[] getValues(double evaluationTime, TreeModel model); +} diff --git a/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractPathDependentProduct.java b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractPathDependentProduct.java new file mode 100644 index 0000000000..8df66a4f92 --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AbstractPathDependentProduct.java @@ -0,0 +1,282 @@ +package net.finmath.tree.assetderivativevaluation.products; + +import java.util.Arrays; + +import net.finmath.montecarlo.RandomVariableFromDoubleArray; +import net.finmath.stochastic.RandomVariable; +import net.finmath.tree.AbstractTreeProduct; +import net.finmath.tree.TreeModel; + +/** + * 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).
      • + *
    • Product-specific path state variables (e.g. running sum/count for Asian).
      • + *
    + * + * 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: + *
      + *
    • 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}
        • + *
      • 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. + */ + 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 RandomVariable spotRV = model.getSpotAtGivenTimeIndexRV(timeIndex); + + for(int parent = 0; parent < parentNodes; parent++) { + + final int parentRecomb = recombIndex[j - 1][parent]; + 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); + final int childStateOffset = child * m; + + evolveState( + state[j - 1], + parentStateOffset, + spotChild, + timeIndex, + isFixingTimeIndex(timeIndex), + state[j], + childStateOffset); + } + } + } + + final int terminalNodes = powInt(B, steps); + final double[] vTerminal = new double[terminalNodes]; + + 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); + + 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); + + 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; + } + + 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(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/AmericanNonPathDependent.java b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AmericanNonPathDependent.java new file mode 100644 index 0000000000..08e2d7cefb --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AmericanNonPathDependent.java @@ -0,0 +1,64 @@ +package net.finmath.tree.assetderivativevaluation.products; + + +import net.finmath.stochastic.RandomVariable; +import net.finmath.tree.TreeModel; + +import java.util.function.DoubleUnaryOperator; + +/** + * American (non–path-dependent) option priced on a recombining TreeModel. + * The product allows early exercise at every lattice time prior to maturity. + * At each step k it compares: + * Continuation value: discounted conditional expectation + * E^Q[V_{k+1} | F_k] provided by the model, and + * Immediate exercise value: f(S_k). + * The option value is the pointwise maximum of the two (early exercise rule). + */ +public class AmericanNonPathDependent extends AbstractNonPathDependentProduct { + + /** + * Creates an American option with given maturity and payoff. + * + * @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 AmericanNonPathDependent(double maturity, DoubleUnaryOperator payOffFunction) { + super(maturity, payOffFunction); + } + + /** + * Performs backward induction with early exercise: + * Set terminal values at maturity: levels[n-k0] = f(S_T). + * For k = n-1,...,k0: + * Continuation: continuation = DF * E^Q[ levels[k+1-k0] | F_k ] via getConditionalExpectationRV(RandomVariable, int) + * Exercise: exercise = f(S_k) via getSpotAtGivenTimeIndexRV(int). + * American rule: levels[k-k0] = max(continuation, exercise). + * + * @param evaluationTime The time at which the value is requested (typically 0.0). + * @param model The tree model (spot lattice, discounting, conditional expectations). + * @return Array of value random variables from index k0 to n; + */ + @Override + public RandomVariable[] getValues(double evaluationTime, TreeModel model) { + final int k0 = timeToIndex(evaluationTime, model); + final int n = model.getNumberOfTimes() - 1; + final RandomVariable[] levels = new RandomVariable[n - k0 + 1]; + levels[n - k0] = model.getTransformedValuesAtGivenTimeRV(model.getLastTime(), getPayOffFunction()); + + for (int k = n - 1; k >= k0; --k) { + + final RandomVariable continuation = model.getConditionalExpectationRV(levels[(k + 1) - k0], k); + + final RandomVariable spotK = model.getSpotAtGivenTimeIndexRV(k); + final RandomVariable exercise = spotK.apply(getPayOffFunction()); + + levels[k -k0] = continuation.floor(exercise); + + } + + return levels; + + } +} + 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..06f93423be --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/AsianOption.java @@ -0,0 +1,75 @@ +package net.finmath.tree.assetderivativevaluation.products; + +/** + * Fixed-strike Asian option using FULL non-recombining tree expansion. + * + * State variables: + *
      + *
    • state[0] = running sum of fixing spots
      • + *
    • state[1] = fixing count
      • + *
    + * + * 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. + */ + 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/EuropeanNonPathDependent.java b/src/main/java/net/finmath/tree/assetderivativevaluation/products/EuropeanNonPathDependent.java new file mode 100644 index 0000000000..69cb014d12 --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/EuropeanNonPathDependent.java @@ -0,0 +1,49 @@ +package net.finmath.tree.assetderivativevaluation.products; + +import net.finmath.stochastic.RandomVariable; +import net.finmath.tree.TreeModel; + +import java.util.function.DoubleUnaryOperator; + +/** + * European (non–path-dependent) option priced on a recombining TreeModel. + * The payoff f(S_T) is applied only at maturity and the price is obtained + * by backward induction using the model’s discounted conditional expectation. + * pricing uses evaluationTime = 0.0 => k0 = 0. + */ +public class EuropeanNonPathDependent extends AbstractNonPathDependentProduct { + + /** + * Creates a European option with given maturity and payoff. + * + * @param maturity Contract maturity (model time units). + * @param payOffFunction Payoff function f(S) applied at T + * (e.g., s -> Math.max(s-K,0.0) for a call). + */ + public EuropeanNonPathDependent(double maturity, DoubleUnaryOperator payOffFunction){ + super(maturity,payOffFunction); + } + + /** + * Performs backward induction under the given tree model: + * @param evaluationTime The time at which the value is requested (usually 0.0). + * @param model The tree model providing spot lattice and CE/discounting. + * @return An array of RandomVariable with values from time index k0 to n. + * By convention, values[0] corresponds to k0 and + * values[n-k0] corresponds to maturity. + */ + @Override + public RandomVariable[] getValues(double evaluationTime, TreeModel model) { + final int k0 = timeToIndex(evaluationTime,model); + final int n = model.getNumberOfTimes()-1; + final RandomVariable[] values = new RandomVariable[n -k0 +1]; + values[n-k0] = model.getTransformedValuesAtGivenTimeRV(model.getLastTime(),getPayOffFunction()); + + for (int timeIndex = n - 1; timeIndex >= k0; timeIndex--) { + values[timeIndex - k0] = model.getConditionalExpectationRV(values[(timeIndex+1)], timeIndex); + } + return values; + } + + } + 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..a35fc1679f --- /dev/null +++ b/src/main/java/net/finmath/tree/assetderivativevaluation/products/FloatingStrikeAsianOption.java @@ -0,0 +1,79 @@ +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/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; diff --git a/src/test/java/net/finmath/tree/EuropeanAndAmericanOptionPricingTest.java b/src/test/java/net/finmath/tree/EuropeanAndAmericanOptionPricingTest.java new file mode 100644 index 0000000000..5baf489d8d --- /dev/null +++ b/src/test/java/net/finmath/tree/EuropeanAndAmericanOptionPricingTest.java @@ -0,0 +1,105 @@ +package net.finmath.tree; + + +import org.junit.Assert; +import org.junit.Test; + +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.EuropeanNonPathDependent; +import net.finmath.tree.assetderivativevaluation.products.AmericanNonPathDependent; + +public class EuropeanAndAmericanOptionPricingTest { + + @Test + public void testNonPathDep_Call() { + + double spot = 100.0; + double rate = 0.05; + double vol = 0.2; + double maturity = 1.0; + int steps = 100; + double strike = 110.0; + double tol = 1e-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); + + + + 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); + double euJR = euCall.getValue(0.0,jr).getAverage(); + double euTRI = euCall.getValue(0.0,tri).getAverage(); + + double usCRR = usCall.getValue(0.0, crr).getAverage(); + double usJR = usCall.getValue(0.0,jr).getAverage(); + double usTRI = usCall.getValue(0.0,tri).getAverage(); + System.out.println(euCRR + " " + euJR + " " + euTRI); + + //"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 + public 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); + + + + 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(); + 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(); + //"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 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); + } +}