Merge branch 'master' of https://github.com/finmath/finmath-lib.git

Author: cfries

src/main/java/net/finmath/climate/models/dice/DICEModel.java

@@ -1,7 +1,11 @@
 package net.finmath.climate.models.dice;
 
 import java.util.Arrays;
+import java.util.Map;
 import java.util.function.DoubleUnaryOperator;
+import java.util.function.Function;
+import java.util.function.IntToDoubleFunction;
+import java.util.function.Predicate;
 import java.util.function.UnaryOperator;
 import java.util.logging.Logger;
 
@@ -56,7 +60,15 @@ public class DICEModel implements ClimateModel {
 	private double[] welfare;
 	private double[] value;
 
-	public DICEModel(TimeDiscretization timeDiscretization, UnaryOperator<Double> abatementFunction, UnaryOperator<Double> savingsRateFunction, double discountRate) {
+	/**
+	 * 
+	 * @param timeDiscretization
+	 * @param abatementFunction
+	 * @param savingsRateFunction
+	 * @param discountRate
+	 * @param modelProperties A key value map of optional model properties or parameters.
+	 */
+	public DICEModel(TimeDiscretization timeDiscretization, UnaryOperator<Double> abatementFunction, UnaryOperator<Double> savingsRateFunction, double discountRate, Map<String, Object> modelProperties) {
 		super();
 		this.timeDiscretization = timeDiscretization;
 		this.abatementFunction = abatementFunction;
@@ -77,14 +89,23 @@ public DICEModel(TimeDiscretization timeDiscretization, UnaryOperator<Double> ab
 		welfare = new double[numberOfTimes];
 		value = new double[numberOfTimes];
 
-		this.init();
+		this.init(modelProperties);
+	}
+
+	public DICEModel(TimeDiscretization timeDiscretization, UnaryOperator<Double> abatementFunction, UnaryOperator<Double> savingsRateFunction, double discountRate) {
+		this(timeDiscretization, abatementFunction, savingsRateFunction, discountRate, Map.of());
 	}
 
 	public DICEModel(TimeDiscretization timeDiscretization, UnaryOperator<Double> abatementFunction) {
 		this(timeDiscretization, abatementFunction, t -> 0.259029014481802, 0.03);
 	}
 
-	private void init() {
+	private void init(Map<String, Object> modelProperties) {
+
+		Predicate<Integer> isTimeIndexToShift = (Predicate<Integer>) modelProperties.getOrDefault("isTimeIndexToShift", (Predicate<Integer>) i -> true);
+		double initialEmissionShift = (double) modelProperties.getOrDefault("initialEmissionShift", 0.0);
+		double initialConsumptionShift = (double) modelProperties.getOrDefault("initialConsumptionShift", 0.0);
+		
 		final double timeStep = timeDiscretization.getTimeStep(0);
 
 		/*
@@ -109,12 +130,13 @@ private void init() {
 		final DoubleUnaryOperator damageFunction = new DamageFromTemperature();
 
 		final EmissionIndustrialIntensityFunction emissionIndustrialIntensityFunction = new EmissionIndustrialIntensityFunction();
-		final EmissionExternalFunction emissionExternalFunction = new EmissionExternalFunction();
+
+		final Function<Double, Double> emissionExternalFunction = new EmissionExternalFunction();
 
 		final EvolutionOfCarbonConcentration evolutionOfCarbonConcentration = new EvolutionOfCarbonConcentration(timeDiscretization);
 
 		final ForcingFunction forcingFunction = new ForcingFunction();
-		final double forcingExternal = 1.0/5.0;		// per year
+		final double forcingExternal = 1.0;
 
 		final EvolutionOfTemperature evolutionOfTemperature = new EvolutionOfTemperature(timeDiscretization);
 
@@ -175,6 +197,10 @@ private void init() {
 			double emissionIndustrial = emissionIndustrialIntensityFunction.apply(time) * gdp[timeIndex];
 			double emissionExternal = emissionExternalFunction.apply(time);
 			emission[timeIndex] = (1 - abatement[timeIndex])/(1-abatement[0]) * emissionIndustrial + emissionExternal;
+
+			// Allow for an external shift to the emissions (e.g. to calculate SCC).
+			emission[timeIndex] += isTimeIndexToShift.test(timeIndex) ? initialEmissionShift : 0.0;
+
 			carbonConcentration[timeIndex+1] = evolutionOfCarbonConcentration.apply(timeIndex, carbonConcentration[timeIndex], emission[timeIndex]);
 
 			/*
@@ -204,9 +230,12 @@ private void init() {
 			// Constant from the original model - in the original model this is a time varying control variable.
 			double savingsRate = savingsRateFunction.apply(time);	//0.259029014481802;
 
-			double consumption = (1-savingsRate) * gdpNet;
+			double consumption = (1-savingsRate) * gdpNet;			
 			double investment = savingsRate * gdpNet;
 
+			// Allow for an external shift to the emissions (e.g. to calculate SCC).
+			consumption += isTimeIndexToShift.test(timeIndex) ? initialConsumptionShift : 0.0;
+					
 			capital[timeIndex+1] = evolutionOfCapital.apply(timeIndex).apply(capital[timeIndex], investment);
 
 			/*

src/main/java/net/finmath/climate/models/dice/submodels/AbatementCostFunction.java

@@ -32,7 +32,7 @@ public class AbatementCostFunction implements BiFunction<Double, Double, Double>
 	 *
 	 * @param backstopPriceInitial USD per ton CO2 to abate 100% per year.
 	 * @param backstopRate Annual rate by which the price declines.
-	 * @param abatementExponent Exponent for &mu;
+	 * @param abatementExponent Exponent for &mu; (sometimes called theta_2)
 	 */
 	public AbatementCostFunction(double backstopPriceInitial, double backstopRate, double abatementExponent) {
 		super();
@@ -49,7 +49,7 @@ public AbatementCostFunction() {
 	@Override
 	public Double apply(Double time, Double abatement) {
 		final double backstopPrice = backstopPriceInitial * Math.exp(-backstopRate * time);
-		final double abatementCost = backstopPrice * Math.pow(abatement , abatementExponent)/abatementExponent;
+		final double abatementCost = backstopPrice * Math.pow(abatement, abatementExponent)/abatementExponent;
 
 		return abatementCost;
 	}

src/main/java/net/finmath/climate/models/dice/submodels/CarbonConcentration3DScalar.java

@@ -1,11 +1,13 @@
 package net.finmath.climate.models.dice.submodels;
 
+import org.apache.commons.lang3.Validate;
+
 import net.finmath.climate.models.CarbonConcentration;
 import net.finmath.stochastic.RandomVariable;
 import net.finmath.stochastic.Scalar;
 
 /**
- * State vector representing carbon concentration.
+ * State vector representing carbon concentration in units of GtC.
  *
  * @author Christian Fries
  */
@@ -17,6 +19,9 @@ public class CarbonConcentration3DScalar implements CarbonConcentration {
 
 	public CarbonConcentration3DScalar(double carbonConcentrationInAtmosphere, double carbonConcentrationInShallowOcean, double carbonConcentrationInLowerOcean) {
 		super();
+		Validate.isTrue(carbonConcentrationInAtmosphere >= 0, "carbonConcentrationInAtmosphere must not be negative.", carbonConcentrationInAtmosphere);
+		Validate.isTrue(carbonConcentrationInShallowOcean >= 0, "carbonConcentrationInShallowOcean must not be negative.", carbonConcentrationInShallowOcean);
+		Validate.isTrue(carbonConcentrationInLowerOcean >= 0, "carbonConcentrationInLowerOcean must not be negative.", carbonConcentrationInLowerOcean);
 		this.carbonConcentrationInAtmosphere = carbonConcentrationInAtmosphere;
 		this.carbonConcentrationInShallowOcean = carbonConcentrationInShallowOcean;
 		this.carbonConcentrationInLowerOcean = carbonConcentrationInLowerOcean;
@@ -26,6 +31,9 @@ public CarbonConcentration3DScalar(double[] carbonConcentration) {
 		this(carbonConcentration[0], carbonConcentration[1], carbonConcentration[2]);
 	}
 
+	/**
+	 * Create a state vector with 851 GtC in atmosphere, 460 GtC in shallow ocean, 1740 in lower ocean.
+	 */
 	public CarbonConcentration3DScalar() {
 		this(851, 460, 1740);		// GtC
 	}

src/main/java/net/finmath/climate/models/dice/submodels/EmissionIndustrialIntensityFunction.java

@@ -24,7 +24,7 @@ public class EmissionIndustrialIntensityFunction implements Function<Double, Dou
 
 	private final double emissionIntensityInitial;		// sigma0;
 	private final double emissionIntensityRateInitial;	// = 0.0152;		// -g	// per year
-	private final double emissionIntensityRateDecay;	// = 0.001;			// -d	// per year
+	private final double emissionIntensityRateDecay;	// exp decay rate corresponding to annual -0.001;			// -d	// per year
 
 	/**
 	 * 

src/main/java/net/finmath/climate/models/dice/submodels/EvolutionOfTemperature.java

@@ -9,22 +9,24 @@
 
 /**
  *
- * The evolution of the temperature \( \mathrm{d}T(t) = \left( \Gamma_{T} T(t) + F(t) \right) \mathrm{d}t \).
+ * The evolution of the temperature \( \mathrm{d}T(t) = \left( \Gamma_{T} T(t) + \xi \cdot F(t) \right) \mathrm{d}t \).
  * 
  * The unit of \( T \) is K (Kelvin).
  * 
- * The evolution is modelled as \( \mathrm{d}T(t) = \left( \Gamma_{T} T(t) + F(t) \right) \mathrm{d}t \).
+ * This is a function of timeIndex, previous temperature and forcing.
+ * 
+ * The evolution is modelled as \( \mathrm{d}T(t) = \left( \Gamma_{T} T(t) + \xi \cdot F(t) \right) \mathrm{d}t \).
  * With the given {@link TimeDiscretization} it is approximated via an Euler-step
  * \(
- * 	T(t_{i+1}) = \Phi T(t_{i}) + (forcingToTemp \cdot (forcing, 0, 0) \cdot \Delta t_{i}
+ * 	T(t_{i+1}) = \Phi T(t_{i}) + (forcingToTemp \cdot (forcing, 0) \cdot \Delta t_{i}
  * \)
  * where \( \Phi = (1 + \Gamma_{T} \Delta t_{i}) \).
  *
  * @author Christian Fries
  */
 public class EvolutionOfTemperature implements TriFunction<Integer, Temperature2DScalar, Double, Temperature2DScalar> {
 
-	private static double forcingToTempDefault = 0.1005;		// sometimes called xi1 or c1		// TODO check role of time step
+	private static double forcingToTemp5YDefault = 0.1005;	// (climate sensitivity) sometimes called xi1 or c1 (original parameter was per 5 year)
 
 	private static double[][] transitionMatrix5YDefault;
 	static {
@@ -33,8 +35,8 @@ public class EvolutionOfTemperature implements TriFunction<Integer, Temperature2
 		final double c3 = 0.088;			// Transfer coefficient upper to lower stratum
 		final double c4 = 0.025;
 
-		final double phi11 = 1-forcingToTempDefault*((fco22x/t2xco2) + c3);
-		final double phi12 = forcingToTempDefault*c3;
+		final double phi11 = 1-forcingToTemp5YDefault*((fco22x/t2xco2) + c3);
+		final double phi12 = forcingToTemp5YDefault*c3;
 		final double phi21 = c4;
 		final double phi22 = 1-c4;
 
@@ -48,7 +50,7 @@ public class EvolutionOfTemperature implements TriFunction<Integer, Temperature2
 	/**
 	 * @param timeDiscretization The time discretization.
 	 * @param transitionMatrices Transition matrix \( \Phi \) for each time step.
-	 * @param forcingToTemp The scaling coefficient for the external forcing.
+	 * @param forcingToTemp The scaling coefficient for the external forcing to temperature change per year (this is per 1Y).
 	 */
 	public EvolutionOfTemperature(TimeDiscretization timeDiscretization, Function<Integer, double[][]> transitionMatrices, double forcingToTemp) {
 		super();
@@ -61,7 +63,7 @@ public EvolutionOfTemperature(TimeDiscretization timeDiscretization) {
 		Function<Integer, Double> timeSteps = ((Integer timeIndex) -> { return timeDiscretization.getTimeStep(timeIndex); });
 		this.timeDiscretization = timeDiscretization;
 		transitionMatrices = timeSteps.andThen(Cached.<Double, double[][]>of(timeStep -> timeStep == 5.0 ? transitionMatrix5YDefault : LinearAlgebra.matrixPow(transitionMatrix5YDefault, (Double)timeStep/5.0)));
-		this.forcingToTemp = forcingToTempDefault;
+		this.forcingToTemp = forcingToTemp5YDefault/5; // Rescale to per 1 Y
 	}
 
 	@Override

src/main/java/net/finmath/climate/models/dice/submodels/ForcingFunction.java

@@ -9,8 +9,9 @@
  */
 public class ForcingFunction implements BiFunction<CarbonConcentration3DScalar, Double, Double> {
 
+	// Parameters of the orignal model
 	private final double carbonConcentrationBase = 580;
-	private final double forcingPerCarbonDoubling = 3.6813/5;			// Parameter in original model was per 5 year
+	private final double forcingPerCarbonDoubling = 3.6813;
 
 	@Override
 	public Double apply(CarbonConcentration3DScalar carbonConcentration, Double forcingExternal) {

src/main/java/net/finmath/climate/models/dice/submodels/Temperature2DScalar.java

@@ -5,7 +5,7 @@
 import net.finmath.stochastic.Scalar;
 
 /**
- * State vector representing temperature above pre-industrial level.
+ * State vector representing temperature above pre-industrial level in Kelvin (K).
  *
  * @author Christian Fries
  */
@@ -14,6 +14,12 @@ public class Temperature2DScalar implements Temperature {
 	private final double temperatureOfAtmosphere;
 	private final double temperatureOfLandAndOcean;
 
+	/**
+	 * Create a temperature vector.
+	 * 
+	 * @param temperatureOfAtmosphere Temperature over pre-industrial of the atmosphere.
+	 * @param temperatureOfLandAndOcean Temperature over pre-industrial of land and ocean.
+	 */
 	public Temperature2DScalar(double temperatureOfAtmosphere, double temperatureOfLandAndOcean) {
 		super();
 		this.temperatureOfAtmosphere = temperatureOfAtmosphere;

src/main/java/net/finmath/equities/marketdata/FlatYieldCurve.java

@@ -25,6 +25,11 @@ public FlatYieldCurve(
 		this.rate = rate;
 		this.dayCounter = dayCounter;
 	}
+	
+	public FlatYieldCurve rollToDate(LocalDate date) 
+	{
+		return new FlatYieldCurve(date, rate, dayCounter);
+	}
 
 	public double getRate(double maturity)
 	{

src/main/java/net/finmath/equities/models/BuehlerDividendForwardStructure.java

@@ -65,7 +65,7 @@ public BuehlerDividendForwardStructure cloneWithNewDate(LocalDate newDate)
 		return new BuehlerDividendForwardStructure(
 				newDate,
 				this.spot,
-				this.repoCurve,
+				this.repoCurve.rollToDate(newDate),
 				this.dividendStream,
 				this.dayCounter);
 	}

src/main/java/net/finmath/equities/pricer/AnalyticOptionValuation.java

@@ -58,8 +58,7 @@ public double calculate(
 		final var ttm = dcc.getDaycountFraction(forwardStructure.getValuationDate(), expiryDate);
 		final var forward = forwardStructure.getForward(expiryDate);
 		final var discountFactor = discountCurve.getDiscountFactor(expiryDate);
-		//var repoRate = forwardStructure.getRepoCurve().getRate(expiryDate);
-		final var repoRate = discountCurve.getRate(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(
@@ -111,14 +110,7 @@ public double calculate(
 					volatility,
 					option.isCallOption(),
 					discountFactor * adjustedForward,
-					repoRate)
-					- repoRate * Black76Model.optionPrice(
-							1.0,
-							adjustedStrike / adjustedForward,
-							ttm,
-							volatility,
-							option.isCallOption(),
-							discountFactor * adjustedForward);
+					discountRate);
 		default:
 			throw new NotImplementedException("Calculation for " + calcType + " not implemented yet.");
 		}

src/main/java/net/finmath/equities/pricer/PdeOptionValuation.java

@@ -401,7 +401,8 @@ else if (i == j - 1 || i == j + 1)
 					* 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 theta = (discountFactor * lastAtmPrice - price) / dt;
+			final var discountFactorTheta = discountCurve.getDiscountFactor(expiryTime - dt);
+			final var theta = (discountFactorTheta * lastAtmPrice - price) / dt;
 			return new double[] {price, delta, gamma, theta};
 		}
 		else

src/main/java/net/finmath/finitedifference/models/FDMConstantElasticityOfVarianceModel.java

@@ -105,9 +105,6 @@ public double getNumStandardDeviations() {
 		return numStandardDeviations;
 	}
 
-	/* (non-Javadoc)
-	 * @see net.finmath.finitedifference.models.FiniteDifference1DModel#valueOptionWithThetaMethod(net.finmath.finitedifference.products.FDMEuropeanCallOption, double)
-	 */
 	@Override
 	public double[][] getValue(double evaluationTime, double time, DoubleUnaryOperator values, FiniteDifference1DBoundary boundary) {
 		final FDMThetaMethod solver = new FDMThetaMethod(this, boundary, time, center, theta);

src/main/java/net/finmath/time/TimeDiscretization.java

@@ -76,8 +76,28 @@ public interface TimeDiscretization extends Iterable<Double> {
 	 */
 	int getTimeIndexNearestGreaterOrEqual(double time);
 
+	/**
+	 * Returns the first time in the time discretization.
+	 * 
+	 * @return The first time in the time discretization.
+	 */
+	default double getFirstTime() {
+		return getTime(0);
+	}
+
+	/**
+	 * Returns the last time in the time discretization.
+	 * 
+	 * @return The last time in the time discretization.
+	 */
+	default double getLastTime() {
+		return getTime(getNumberOfTimes()-1);
+	}
+
+	
 	/**
 	 * Return a clone of this time discretization as <code>double[]</code>.
+	 * 
 	 * @return The time discretization as <code>double[]</code>
 	 */
 	double[] getAsDoubleArray();

src/test/java/net/finmath/equities/AnalyticOptionValuationTest.java

@@ -133,8 +133,9 @@ public void Test_sensis() throws CalculationException
 			final var anaTheta = pricer.getTheta(option, fwdStructure, discountCurve, flatVol);
 			final var thetaDate = valDate.plusDays(1);
 			final var thetaSpot = fwdStructure.getForward(thetaDate);
+			final var thetaCurve = discountCurve.rollToDate(thetaDate);
 			final var shiftedFwdStructure = fwdStructure.cloneWithNewSpot(thetaSpot).cloneWithNewDate(thetaDate);
-			final var priceTheta = pricer.getPrice(option, shiftedFwdStructure, discountCurve, flatVol);
+			final var priceTheta = pricer.getPrice(option, shiftedFwdStructure, thetaCurve, flatVol);
 			final var fdTheta = (priceTheta - price) / dcc.getDaycountFraction(valDate, thetaDate);
 
 

src/test/java/net/finmath/equities/PdeOptionPricerTest.java

@@ -345,7 +345,7 @@ public void Test_europeanSensis() throws CalculationException
 			assertEquals("Pde Vega deviates to much from analytic value for "  + (isCall ? "Call" : "Put"),
 					0.0, pdeVega/anaVega - 1.0, 0.02);
 			assertEquals("Pde Theta deviates to much from analytic value for "  + (isCall ? "Call" : "Put"),
-					0.0, pdeSensis[3]/anaTheta - 1.0, 0.02);
+					0.0, pdeSensis[3]/anaTheta - 1.0, 0.04);
 		}
 	}