Merge pull request #93 from maxsing/master

Bugfix DICEModel and submodels, added external forcing

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

@@ -10,33 +10,22 @@
 import java.util.logging.Logger;
 
 import net.finmath.climate.models.ClimateModel;
-import net.finmath.climate.models.dice.submodels.AbatementCostFunction;
-import net.finmath.climate.models.dice.submodels.CarbonConcentration3DScalar;
-import net.finmath.climate.models.dice.submodels.DamageFromTemperature;
-import net.finmath.climate.models.dice.submodels.EmissionExternalFunction;
-import net.finmath.climate.models.dice.submodels.EmissionIndustrialIntensityFunction;
-import net.finmath.climate.models.dice.submodels.EvolutionOfCapital;
-import net.finmath.climate.models.dice.submodels.EvolutionOfCarbonConcentration;
-import net.finmath.climate.models.dice.submodels.EvolutionOfPopulation;
-import net.finmath.climate.models.dice.submodels.EvolutionOfProductivity;
-import net.finmath.climate.models.dice.submodels.EvolutionOfTemperature;
-import net.finmath.climate.models.dice.submodels.ForcingFunction;
-import net.finmath.climate.models.dice.submodels.Temperature2DScalar;
+import net.finmath.climate.models.dice.submodels.*;
 import net.finmath.stochastic.RandomVariable;
 import net.finmath.stochastic.Scalar;
 import net.finmath.time.TimeDiscretization;
 
 /**
  * A simulation of a simplified DICE model (see <code>net.finmath.climate.models.dice.DICEModelTest</code> in src/test/java) for an example usage.
- * 
+ *
  * The model just composes the sub-models (evolution equations and functions) from the package {@link net.finmath.climate.models.dice.submodels}.
- * 
- * Note: The code makes some small simplification: it uses a constant savings rate and a constant external forcings.
+ *
+ * Note: The code uses exponential discounting.
  */
 public class DICEModel implements ClimateModel {
 
 	private static Logger logger = Logger.getLogger("net.finmath.climate");
-	
+
 	/*
 	 * Input to this class
 	 */
@@ -61,7 +50,7 @@ public class DICEModel implements ClimateModel {
 	private double[] value;
 
 	/**
-	 * 
+	 *
 	 * @param timeDiscretization
 	 * @param abatementFunction
 	 * @param savingsRateFunction
@@ -105,7 +94,7 @@ 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);
 
 		/*
@@ -129,14 +118,14 @@ private void init(Map<String, Object> modelProperties) {
 		// Model that describes the damage on the GBP as a function of the temperature-above-normal
 		final DoubleUnaryOperator damageFunction = new DamageFromTemperature();
 
-		final EmissionIndustrialIntensityFunction emissionIndustrialIntensityFunction = new EmissionIndustrialIntensityFunction();
+		final EmissionIndustrialIntensityFunction emissionIndustrialIntensityFunction = new EmissionIndustrialIntensityFunction(timeStep);
 
 		final Function<Double, Double> emissionExternalFunction = new EmissionExternalFunction();
 
 		final EvolutionOfCarbonConcentration evolutionOfCarbonConcentration = new EvolutionOfCarbonConcentration(timeDiscretization);
 
 		final ForcingFunction forcingFunction = new ForcingFunction();
-		final double forcingExternal = 1.0;
+		final ForcingExternalFunction forcingExternalFunction = new ForcingExternalFunction();
 
 		final EvolutionOfTemperature evolutionOfTemperature = new EvolutionOfTemperature(timeDiscretization);
 
@@ -171,6 +160,8 @@ private void init(Map<String, Object> modelProperties) {
 		population[0] = L0;
 		productivity[0] = A0;
 		double utilityDiscountedSum = 0;
+		//Emission intensity initial value, sigma(0) = e0/q0
+		double emissionIntensity = 35.85/105.5;
 
 		/*
 		 * Evolve
@@ -182,27 +173,24 @@ private void init(Map<String, Object> modelProperties) {
 			 * Evolve geo-physical quantities i -> i+1 (as a function of gdp[i])
 			 */
 
-			// Temperature
-
-			final double forcing = forcingFunction.apply(carbonConcentration[timeIndex], forcingExternal);
-			temperature[timeIndex+1] = evolutionOfTemperature.apply(timeIndex, temperature[timeIndex], forcing);
-
-
 			// Abatement
-
 			abatement[timeIndex] = abatementFunction.apply(timeDiscretization.getTime(timeIndex));
 
 			// Carbon
-
-			double emissionIndustrial = emissionIndustrialIntensityFunction.apply(time) * gdp[timeIndex];
+			double emissionIndustrial = emissionIntensity/(1-abatement[0]) * gdp[timeIndex];
 			double emissionExternal = emissionExternalFunction.apply(time);
-			emission[timeIndex] = (1 - abatement[timeIndex])/(1-abatement[0]) * emissionIndustrial + emissionExternal;
+			emission[timeIndex] = (1 - abatement[timeIndex]) * 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]);
 
+			// Temperature
+			double forcingExternal = forcingExternalFunction.apply(time+timeStep);
+			final double forcing = forcingFunction.apply(carbonConcentration[timeIndex+1], forcingExternal);
+			temperature[timeIndex+1] = evolutionOfTemperature.apply(timeIndex, temperature[timeIndex], forcing);
+
 			/*
 			 * Cost
 			 */
@@ -215,27 +203,30 @@ private void init(Map<String, Object> modelProperties) {
 			/*
 			 * Evolve economy i -> i+1 (as a function of temperature[i])
 			 */
-			
+
 			// Remaining gdp
 			double gdpNet = gdp[timeIndex] - damageCostAbsolute - abatementCostAbsolute;
 
 			/*
 			 * Equivalent (alternative way) to calculate the abatement
 			 */
-			double abatementCost = abatementCostFunction.apply(time, abatement[timeIndex]) * emissionIndustrialIntensityFunction.apply(time);
+			double abatementCost = abatementCostFunction.apply(time, abatement[timeIndex]) * emissionIntensity/(1-abatement[0]);
 			double gdpNet2 = gdp[timeIndex] * (1-damage[timeIndex] - abatementCost);
 			if(Math.abs(gdpNet2-gdpNet)/(1+Math.abs(gdpNet)) > 1E-10) logger.warning("Calculation of relative and absolute net GDP does not match.");
 
-
+			// Evolve emission intensity
+			emissionIntensity = emissionIndustrialIntensityFunction.apply(time, emissionIntensity);
+
+
 			// 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);
 
 			/*
@@ -253,7 +244,7 @@ private void init(Map<String, Object> modelProperties) {
 			 */
 			double alpha = 1.45;           // Elasticity of marginal utility of consumption (GAMS elasmu)
 			double C = consumption;
-			double utility = L*(Math.pow(C / (L/1000),1-alpha)-1)/(1-alpha);
+			double utility = population[timeIndex]*( (Math.pow(1000*C/(population[timeIndex]),1-alpha)-1) /(1-alpha)-1 );
 
 			/*
 			 * Discounted utility

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

@@ -17,8 +17,7 @@ public class DamageFromTemperature implements DoubleUnaryOperator {
 
 	/**
 	 * Create the damage function
-	 * \( T \mapsto \Omega(T) = D(T) / (1+D(T)) \), with
-	 * \( D(T) = a_{0} + a_{1} T + a_{2} T^{2} \), with
+	 * \( T \mapsto \Omega(T) = a_{0} + a_{1} T + a_{2} T^{2} \), with
 	 * \( T \) being temperature above pre-industrial.
 	 *
 	 * @param tempToDamage0 The constant term.
@@ -33,7 +32,7 @@ public DamageFromTemperature(double tempToDamage0, double tempToDamage1, double
 	}
 
 	/**
-	 * Create the damage function \( T \mapsto (a_{0} + a_{1} T + a_{2} T^{2})/(1+a_{0} + a_{1} T + a_{2} T^{2}) \), with \( T \) being temperature above pre-industrial,
+	 * Create the damage function \( T \mapsto (a_{0} + a_{1} T + a_{2} T^{2}) \), with \( T \) being temperature above pre-industrial,
 	 * using the default DICE (2016) parameters.
 	 */
 	public DamageFromTemperature() {
@@ -49,7 +48,7 @@ public DamageFromTemperature() {
 	public double applyAsDouble(double temperature) {
 		final double damage = tempToDamage0 + tempToDamage1 * temperature + tempToDamage2 * temperature * temperature;
 
-		return damage / (1+damage);
+		return damage;
 	}
 
 }

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

@@ -1,10 +1,10 @@
 package net.finmath.climate.models.dice.submodels;
 
-import java.util.function.Function;
+import java.util.function.BiFunction;
 
 /**
  * The function that maps time to emission intensity \( \sigma(t) \) (in kgCO2 / USD = GtCO2 / (10^12 USD)).
- * 
+ *
  * The emission intensity is the factor that is applied to the GDP to get the corresponding emissions.
  *
  * The function is modelled as an exponential decay, where the decay rate decays exponentially (double exponential).
@@ -13,7 +13,7 @@
  *
  * @author Christian Fries
  */
-public class EmissionIndustrialIntensityFunction implements Function<Double, Double> {
+public class EmissionIndustrialIntensityFunction implements BiFunction<Double, Double, Double> {
 
 	private static double e0 = 35.85;					// Initial emissions
 	private static double q0 = 105.5;					// Initial global output
@@ -22,33 +22,62 @@ public class EmissionIndustrialIntensityFunction implements Function<Double, Dou
 	//	private static double mu0 = 0.03;					// Initial mitigation rate
 	//	private static double sigma0 = e0/(q0*(1-mu0));		// Calculated initial emissions intensity
 
+	private final double timeStep;
 	private final double emissionIntensityInitial;		// sigma0;
 	private final double emissionIntensityRateInitial;	// = 0.0152;		// -g	// per year
 	private final double emissionIntensityRateDecay;	// exp decay rate corresponding to annual -0.001;			// -d	// per year
 
 	/**
-	 * 
-	 * @param emissionIntensityInitial The initial emission intensity. Unit: GtCO2 / (10^12 USD)
+	 * The evolution of the emission intensity
+	 * @param timeStep The size of one timeStep.
 	 * @param emissionIntensityRateInitial
 	 * @param emissionIntensityRateDecay
 	 */
-	public EmissionIndustrialIntensityFunction(double emissionIntensityInitial, double emissionIntensityRateInitial, double emissionIntensityRateDecay) {
+	public EmissionIndustrialIntensityFunction(double timeStep, double emissionIntensityInitial,
+	                                           double emissionIntensityRateInitial, double emissionIntensityRateDecay) {
 		super();
+		this.timeStep = timeStep;
 		this.emissionIntensityInitial = emissionIntensityInitial;
 		this.emissionIntensityRateInitial = emissionIntensityRateInitial;
 		this.emissionIntensityRateDecay = emissionIntensityRateDecay;
 	}
 
+	/**
+	 * @deprecated
+	 */
+	public EmissionIndustrialIntensityFunction(double emissionIntensityInitial, double emissionIntensityRateInitial,
+	                                           double emissionIntensityRateDecay) {
+		this(5, emissionIntensityInitial, emissionIntensityRateInitial, emissionIntensityRateDecay);
+	}
+
+	/**
+	 * @deprecated
+	 */
 	public EmissionIndustrialIntensityFunction() {
+		// These parameters give a good approximation of the original sigma(t), when using the apply(time) function
+		this(5, sigma0, 0.0152321293038455, 0.000491139147827816);
+	}
+
+	public EmissionIndustrialIntensityFunction(double timeStep) {
 		// Parameters from original model
-		this(sigma0, 0.0152, -Math.log(1-0.001));
+		this(timeStep, sigma0, 0.0152, -Math.log(1-0.001));
 	}
 
-	@Override
+	/**
+	 * @deprecated
+	 */
 	public Double apply(Double time) {
 		final double emissionIntensityRate = emissionIntensityRateInitial * Math.exp(-emissionIntensityRateDecay * time);
 		final double emissionIntensity = emissionIntensityInitial * Math.exp(-emissionIntensityRate * time);
 
 		return emissionIntensity;
 	}
+
+	@Override
+	public Double apply(Double time, Double emissionIntensity) {
+		final double emissionIntensityRate = emissionIntensityRateInitial * Math.exp(-emissionIntensityRateDecay * time);
+
+		return emissionIntensity*Math.exp(-emissionIntensityRate*timeStep);
+	}
+
 }

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

@@ -25,15 +25,13 @@ public EvolutionOfCapital(TimeDiscretization timeDiscretization, double capitalD
 	}
 
 	public EvolutionOfCapital(TimeDiscretization timeDiscretization) {
-		// capital deprecation per 1 year: 5th root of (1-0.1) per 5 years
-		this(timeDiscretization, -Math.log(1-0.1)/5.0);
+		// capital deprecation per 1 year
+		this(timeDiscretization, -Math.log(1-0.1));
 	}
 
 	@Override
 	public BiFunction<Double, Double, Double> apply(Integer timeIndex) {
 		double timeStep = timeDiscretization.getTimeStep(timeIndex);
-		return (Double capital, Double investment) -> {
-			return capital * Math.exp(-capitalDeprecation * timeStep) + investment * timeStep;
-		};
+		return (Double capital, Double investment) -> capital * Math.exp(-capitalDeprecation * timeStep) + investment * timeStep;
 	}
 }

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

@@ -9,9 +9,9 @@
 
 /**
  * The evolution of the carbon concentration M with a given emission E \( \mathrm{d}M(t) = \left( \Gamma_{M} M(t) + E(t) \right) \mathrm{d}t \).
- * 
+ *
  * The unit of \( M \) is GtC (Gigatons of Carbon).
- * 
+ *
  * The evolution is modelled as \( \mathrm{d}M(t) = \left( \Gamma_{M} M(t) + E(t) \right) \mathrm{d}t \right).
  * With the given {@link TimeDiscretization} it is approximated via an Euler-step
  * \(
@@ -30,7 +30,7 @@
  */
 public class EvolutionOfCarbonConcentration implements TriFunction<Integer, CarbonConcentration3DScalar, Double, CarbonConcentration3DScalar> {
 
-	private static double conversionGtCperGtCO2 = 3.0/11.0;
+	private static double conversionGtCperGtCO2 = 1/3.666;
 
 	private static double[][] transitionMatrix5YDefault;
 	// Original transition matrix is a 5Y transition matrix
@@ -47,7 +47,7 @@ public class EvolutionOfCarbonConcentration implements TriFunction<Integer, Carb
 		final double zeta22 = 1 - zeta12 - b23;
 		final double zeta32 = b23;
 		final double zeta23 = b23*(mueq/mleq);
-		final double zeta33 = 1 - b23;
+		final double zeta33 = 1 - zeta23;
 
 		transitionMatrix5YDefault = new double[][] { new double[] { zeta11, zeta12, 0.0 }, new double[] { zeta21, zeta22, zeta23 }, new double[] { 0.0, zeta32, zeta33 } };
 	}

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

@@ -21,9 +21,9 @@ public class EvolutionOfProductivity implements Function<Integer, Function<Doubl
 	/**
 	 * The evolution of the productivity (economy)
 	 * \(
-	 * 	A(t_{i+1}) = A(t_{i}) / (1 - ga * \exp(- deltaA * t))
+	 * 	A(t_{i+1}) = A(t_{i}) / (1 - ga * \exp(- deltaA * t))^{\frac{\delta t}{5}}
 	 * \)
-	 * 
+	 *
 	 * @param timeDiscretization The time discretization.
 	 * @param productivityGrowthRateInitial The initial productivity growth rate ga per 1Y.
 	 * @param productivityGrowthRateDecayRate The productivity growth decay rate per 1Y.
@@ -36,8 +36,8 @@ public EvolutionOfProductivity(TimeDiscretization timeDiscretization, double pro
 	}
 
 	public EvolutionOfProductivity(TimeDiscretization timeDiscretization) {
-		// Parameters from original model: initial productivity growth 0.076 per 5 years, decaying with 0.005 per 1 year.
-		this(timeDiscretization, 1-Math.pow(1-0.076,1.0/5.0), 0.005);
+		// Parameters from original model: initial productivity growth 0.076 per 1 year, decaying with 0.005 per 1 year.
+		this(timeDiscretization, 0.076, 0.005);
 	}
 
 	@Override
@@ -46,7 +46,7 @@ public Function<Double, Double> apply(Integer timeIndex) {
 		double timeStep = timeDiscretization.getTimeStep(timeIndex);
 		return (Double productivity) -> {
 			double productivityGrowthRate = productivityGrowthRateInitial * Math.exp(-productivityGrowthRateDecayRate * time);
-			return productivity * Math.pow(1 - productivityGrowthRate, -timeStep);
+			return productivity / (Math.exp(Math.log(1 - (productivityGrowthRate))*timeStep/5.0));
 		};
 	}
 }

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

@@ -0,0 +1,30 @@
+package net.finmath.climate.models.dice.submodels;
+
+import java.util.function.Function;
+
+/**
+ * The function models the external forcing as a linear function capped at 1.0
+ *
+ * @author Maximilian Singhof
+ */
+public class ForcingExternalFunction  implements Function<Double, Double> {
+
+	private final double intercept;
+	private final double linear;
+
+	public ForcingExternalFunction(double intercept, double linear) {
+		this.intercept = intercept;
+		this.linear = linear;
+	}
+
+	public ForcingExternalFunction() {
+		// Parameters from original model
+		this(0.5, 0.5/(5*17));
+	}
+
+	@Override
+	public Double apply(Double time) {
+		return Math.min(intercept + linear*time,1.0);
+	}
+}
+