From 3f278de42d9486d01faabe618724cadfab698634 Mon Sep 17 00:00:00 2001
From: Bart Arendarczyk <s1924442@ed.ac.uk>
Date: Tue, 2 May 2023 14:45:13 +0100
Subject: [PATCH] Stocks can no longer go negative. Improved price update.

---
 src/ac/ed/lurg/InternationalMarket.java    |  2 +-
 src/ac/ed/lurg/ModelConfig.java            |  2 +-
 src/ac/ed/lurg/country/CountryAgent.java   | 75 +++++++++-------------
 src/ac/ed/lurg/country/CountryManager.java |  4 ++
 src/ac/ed/lurg/country/GlobalPrice.java    | 71 ++++++++++++--------
 src/ac/ed/lurg/types/CropType.java         | 22 +++----
 6 files changed, 91 insertions(+), 85 deletions(-)

diff --git a/src/ac/ed/lurg/InternationalMarket.java b/src/ac/ed/lurg/InternationalMarket.java
index 36d9bd57..2a1436f7 100644
--- a/src/ac/ed/lurg/InternationalMarket.java
+++ b/src/ac/ed/lurg/InternationalMarket.java
@@ -143,7 +143,7 @@ public class InternationalMarket {
 			LogWriter.println(timestep.getYear() + " Updating " + crop.getGamsName() + " prices", 2);
 			GlobalPrice adjustedPrice = prevPrice.createWithUpdatedMarketPrices(imports, exportsBeforeTransportLosses,
 					timestep, totalProduction.get(crop), true, crop.getStockToUseRatio());
-			LogWriter.println( String.format("Price for %s updated from %s to %s \n", crop.getGamsName(), prevPrice, adjustedPrice), 2);
+			LogWriter.println( String.format("Price for %s updated from %s \nto %s \n", crop.getGamsName(), prevPrice, adjustedPrice), 2);
 			if (adjustedPrice.getStockLevel() < 0)
 				LogWriter.println("Global stocks are below zero" + crop.getGamsName() + ", " + timestep.getYear(), 2);
 
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 9de22fae..485f730a 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -389,7 +389,7 @@ public class ModelConfig {
 	public static final ModelFitType DEMAND_ANIMAL_PROD_FIT = ModelFitType.findByName(getProperty("DEMAND_ANIMAL_PROD_FIT", "loglinear"));
 	public static final ModelFitType DEMAND_NON_ANIMAL_PROD_FIT = ModelFitType.findByName(getProperty("DEMAND_NON_ANIMAL_PROD_FIT", "loglinear"));
 	public static final boolean LIMIT_DEMAND_FRACTION = getBooleanProperty("LIMIT_DEMAND_FRACTION", true);
-	public static final boolean DEMAND_FRACT_BY_COST = getBooleanProperty("DEMAND_FRACT_BY_COST", false);;
+	public static final boolean DEMAND_FRACT_BY_COST = getBooleanProperty("DEMAND_FRACT_BY_COST", true);;
 	public static final double MAX_INCOME_PROPORTION_FOOD_SPEND = getDoubleProperty("MAX_INCOME_PROPORTION_FOOD_SPEND", 0.7);
 	
 	public static final double PASTURE_HARVEST_FRACTION = getDoubleProperty("PASTURE_HARVEST_FRACTION", 0.5);
diff --git a/src/ac/ed/lurg/country/CountryAgent.java b/src/ac/ed/lurg/country/CountryAgent.java
index 75c4dfd5..5f557b17 100644
--- a/src/ac/ed/lurg/country/CountryAgent.java
+++ b/src/ac/ed/lurg/country/CountryAgent.java
@@ -164,12 +164,12 @@ public class CountryAgent extends AbstractCountryAgent {
 		for (Map.Entry<CropType, CropUsageData> entry : previousGamsRasterOutput.getCropUsageData().entrySet()) {
 			CropUsageData cropUsage = entry.getValue();
 			CropType crop = entry.getKey();
+			GlobalPrice cropPrice = currentWorldPrices.get(crop);
 			double baseTrade = cropUsage.getNetImportsExpected();
 			double changeUp, changeDown;
 
 			// max of supply overall, needed for when imports are supplying feed and change is zero to allow for production to replace imports
-			// minimum 1 so trade balance can still change when supply is 0
-			double maxOfProdOrSupply = Math.max(cropUsage.getProductionExpected() + Math.max(baseTrade, 0), 1);
+			double maxOfProdOrSupply = cropUsage.getProductionExpected() + Math.max(baseTrade, 0);
 
 			if (ModelConfig.IS_CALIBRATION_RUN && currentTimestep.getTimestep() <= ModelConfig.END_FIRST_STAGE_CALIBRATION) {
 				changeUp = changeDown = 0;
@@ -177,24 +177,18 @@ public class CountryAgent extends AbstractCountryAgent {
 				changeDown = changeUp = maxOfProdOrSupply * ModelConfig.MAX_IMPORT_CHANGE;
 
 				if (crop.isImportedCrop()) {
-					// Guide export level towards import level
-					if (baseTrade < 0) {
-						if (crop.equals(CropType.RUMINANTS)) {
-							int foo = 1;
-						}
-						double ratio = currentWorldPrices.get(crop).getExportAdjustment(crop.getStockToUseRatio());
-						baseTrade *= ratio;
-						cropUsage.updateNetImports(cropUsage.getNetImportsExpected() * ratio);
-					}
+					changeUp = cropPrice.getTradeChangeUp(maxOfProdOrSupply);
+					changeDown = cropPrice.getTradeChangeDown(maxOfProdOrSupply, crop.getStockToUseRatio());
 
-					// If stock about to go negative temporarily increase import price
-					if (currentWorldPrices.get(crop).isStockCritical()) {
-						currentCountryPrices.get(crop).adjustImportPrice(2.0);
+					if (changeUp < 0) { // need to reduce imports to avoid negative stocks
+						baseTrade += changeUp;
+						cropUsage.updateNetImports(baseTrade);
+						changeUp = 0;
 					}
 				}
 			}
 
-			if (CropType.ENERGY_CROPS.equals(crop) && baseTrade == 0) { // could apply this logic for all crops?
+			if (CropType.ENERGY_CROPS.equals(crop) && baseTrade == 0) { // need this in case we start with no global production
 				changeDown = changeUp = ModelConfig.MAX_IMPORT_CHANGE * currentGen2EcDemand;
 			}
 			
@@ -217,7 +211,8 @@ public class CountryAgent extends AbstractCountryAgent {
 
 			importConstraints.put(crop, tradeConstraint);
 		}
-		
+
+		// Carbon trade constraints
 		TradeConstraint carbonTradeConstraint;
 		CarbonUsageData carbonUsageData = previousGamsRasterOutput.getCarbonUsageData();
 		
@@ -228,19 +223,14 @@ public class CountryAgent extends AbstractCountryAgent {
 			if (Math.abs(baseTrade) < 1e-3) { // to set initial trade when starting from zero
 				changeUp = changeDown = getCurrentCarbonDemand();
 			} else {
-				double maxOfProdOrSupply = Math.max(carbonUsageData.getCarbonCredits() + Math.max(baseTrade, 0), 1);
-				changeUp = changeDown = maxOfProdOrSupply * ModelConfig.MAX_IMPORT_CHANGE;
-
-				// Guide export level towards import level
-				if (baseTrade < 0) {
-					double ratio = globalCarbonPrice.getExportAdjustment(ModelConfig.DEFAULT_STOCK_USE_RATIO);
-					baseTrade *= ratio;
-					carbonUsageData.updateNetImports(carbonUsageData.getNetCarbonImport() * ratio);
-				}
-
-				// If stock about to go negative temporarily increase import price
-				if (globalCarbonPrice.isStockCritical()) {
-					currentCarbonPrice.adjustImportPrice(2.0);
+				double maxOfProdOrSupply = carbonUsageData.getCarbonCredits() + Math.max(baseTrade, 0);
+				changeUp = globalCarbonPrice.getTradeChangeUp(maxOfProdOrSupply);
+				changeDown = globalCarbonPrice.getTradeChangeDown(maxOfProdOrSupply, ModelConfig.DEFAULT_STOCK_USE_RATIO);
+
+				if (changeUp < 0) { // need to reduce imports to avoid negative stocks
+					baseTrade += changeUp;
+					carbonUsageData.updateNetImports(baseTrade);
+					changeUp = 0;
 				}
 			}
 
@@ -250,9 +240,8 @@ public class CountryAgent extends AbstractCountryAgent {
 			carbonTradeConstraint = new TradeConstraint(0, 0);
 		}
 
-		Map<WoodType, TradeConstraint> woodTradeConstraints = new HashMap<WoodType, TradeConstraint>();
-
 		// Timber import/export constraints
+		Map<WoodType, TradeConstraint> woodTradeConstraints = new HashMap<WoodType, TradeConstraint>();
 		Map<WoodType, WoodUsageData> woodUsageData = previousGamsRasterOutput.getWoodUsageData();
 
 		for (Map.Entry<WoodType, WoodUsageData> entry : woodUsageData.entrySet()) {
@@ -262,27 +251,21 @@ public class CountryAgent extends AbstractCountryAgent {
 				continue;
 			}
 			WoodUsageData woodUsage = entry.getValue();
+			GlobalPrice woodPrice = globalWoodPrices.get(woodType);
 			double baseTrade = woodUsage.getNetImport();
-			double changeUp;
-			double changeDown;
-
-			double maxOfProdOrSupply = Math.max(woodUsage.getProduction() + Math.max(baseTrade, 0), 1);
+			double changeUp, changeDown;
+			double maxOfProdOrSupply = woodUsage.getProduction() + Math.max(baseTrade, 0);
 
 			if (ModelConfig.IS_CALIBRATION_RUN && currentTimestep.getTimestep() <= ModelConfig.END_FIRST_STAGE_CALIBRATION) {
 				changeUp = changeDown = 0.0;
 			} else {
-				changeUp = changeDown = maxOfProdOrSupply * ModelConfig.MAX_IMPORT_CHANGE;
-
-				// Guide export level towards import level
-				if (baseTrade < 0) {
-					double ratio = globalWoodPrices.get(woodType).getExportAdjustment(ModelConfig.DEFAULT_STOCK_USE_RATIO);
-					baseTrade *= ratio;
-					woodUsageData.get(woodType).updateNetImports(woodUsageData.get(woodType).getNetImport() * ratio);
-				}
+				changeUp = woodPrice.getTradeChangeUp(maxOfProdOrSupply);
+				changeDown = woodPrice.getTradeChangeDown(maxOfProdOrSupply, ModelConfig.DEFAULT_STOCK_USE_RATIO);
 
-				// If stock about to go negative temporarily increase import price
-				if (globalWoodPrices.get(woodType).isStockCritical()) {
-					currentWoodPrices.get(woodType).adjustImportPrice(2.0);
+				if (changeUp < 0) { // need to reduce imports to avoid negative stocks
+					baseTrade += changeUp;
+					woodUsage.updateNetImports(baseTrade);
+					changeUp = 0;
 				}
 			}
 
diff --git a/src/ac/ed/lurg/country/CountryManager.java b/src/ac/ed/lurg/country/CountryManager.java
index 14a5a8e6..2b537143 100644
--- a/src/ac/ed/lurg/country/CountryManager.java
+++ b/src/ac/ed/lurg/country/CountryManager.java
@@ -86,4 +86,8 @@ public class CountryManager {
 		return codeMap.containsKey(isoCode);
 	}
 
+	public int getCountryCount() {
+		return mapFromCompositeCountry.size();
+	}
+
 }
diff --git a/src/ac/ed/lurg/country/GlobalPrice.java b/src/ac/ed/lurg/country/GlobalPrice.java
index 6c47365e..a8719fd7 100644
--- a/src/ac/ed/lurg/country/GlobalPrice.java
+++ b/src/ac/ed/lurg/country/GlobalPrice.java
@@ -19,9 +19,11 @@ public class GlobalPrice implements Serializable {
 	private double stockUseRatio;
 	private double referencePrice;
 	private double prevStockChange;
+	private double production;
 
 	private GlobalPrice(Timestep timestep, double exportPrice, double stockLevel, double importAmount, double exportAmountBeforeLoss,
-						double transportLosses, double referencePrice, double stockUseRatio, double prevStockChange) {
+						double transportLosses, double referencePrice, double stockUseRatio, double prevStockChange,
+						double production) {
 		this.timestep = timestep;
 		this.exportPrice = exportPrice;
 		this.stockLevel = stockLevel;
@@ -31,12 +33,13 @@ public class GlobalPrice implements Serializable {
 		this.transportLosses = transportLosses;
 		this.referencePrice = referencePrice;
 		this.prevStockChange = prevStockChange;
+		this.production = production;
 		LogWriter.println("Created " + this, 3);
 	}
 	
 	public static GlobalPrice createInitial(double exportPrice, double intitalStock) {
 		return new GlobalPrice(new Timestep(ModelConfig.START_TIMESTEP-1), exportPrice, intitalStock, Double.NaN,
-				Double.NaN, Double.NaN, Double.NaN, Double.NaN, 0.0);
+				Double.NaN, Double.NaN, Double.NaN, Double.NaN, 0.0, intitalStock / 0.1);
 	}
 	
 	public double getExportPrice(){
@@ -53,6 +56,8 @@ public class GlobalPrice implements Serializable {
 		double tradeBarrierMultiplier = ModelConfig.updateParameterForShocks(currentYear, "TRADE_BARRIER_MULTIPLIER");
 		double transportLossesRate =ModelConfig.updateParameterForShocks(currentYear, "TRANSPORT_LOSSES");
 		double transportCosts = ModelConfig.updateParameterForShocks(currentYear, "TRANSPORT_COST");
+
+		//transportCosts = timestep.getYear() >= 2023 ? transportCosts * 1.5 : transportCosts;
 		
 		double importPrice = exportPrice * (1 + transportCosts) * (1.0 + countryTradeBarrier*tradeBarrierMultiplier) / (1.0 - transportLossesRate);
 		return importPrice;
@@ -89,9 +94,9 @@ public class GlobalPrice implements Serializable {
 			else
 				updatedStock = stockLevel + stockChange;
 
-			double targetStock = (production + newImports - newExportAmountBeforeLoss) * targetStockUseRatio;
+			double targetStock = production * targetStockUseRatio;
 			double updatedStockUseRatio = production > 0 ? updatedStock / production : 0;
-				
+
 			LogWriter.println(String.format("     imports %.2f, exports %.2f", newImports, newExportAmountBeforeLoss - newTransportLosses), 2);
 			LogWriter.println(String.format("     updatedStock %.2f, previous %.2f (last timestep %.2f), stockChange %.2f, oldDiff %.2f", updatedStock, stockLevel, stockLevel-oldDiff, stockChange, oldDiff), 2);
 			
@@ -101,14 +106,25 @@ public class GlobalPrice implements Serializable {
 				adjustment = 1;
 			}
 			else if (ModelConfig.PRICE_UPDATE_METHOD.equals("StockFlow")) {
-				// Sum of target stock and trade flow volume. Trying to avoid large swings in prices when stocks low
-				double stockAndFlow = stockLevel + Math.max(newImports, newExportAmountBeforeLoss);
+				// Market lambda lower when the rate of change in stock is decreasing
+				// Acts as a damping mechanism to reduce price fluctuation
+				double lambda = Math.atan2(Math.abs(stockChange), Math.abs(prevStockChange)) / (0.5 * Math.PI) * ModelConfig.MARKET_LAMBDA;
+
+				// Prices decrease more slowly than they increase (20% slower according to FAO food price index)
+				if (stockChange > 0)
+					lambda *= 0.8;
+
+				LogWriter.println("lambda="+lambda);
+
+				// Adjust price based on stock change
+				adjustment = 1 - lambda * stockChange / production;
 
-				// (1) stockChange analogous to velocity of change
-				// (2) (stockChange - prevStockChange) analogous to acceleration of change
-				// We need (2) to avoid large changes in prices when trade imbalance is already reducing
-				adjustment = 1 - ModelConfig.MARKET_LAMBDA * (stockChange + 0.5 * (stockChange - prevStockChange)) / stockAndFlow;
+				// If below target stock and stocks not already increasing, nudge prices to increase stock.
+				if (updatedStock < targetStock & stockChange < 0) {
+					adjustment *= 1 - 0.5 * ModelConfig.MARKET_LAMBDA * (updatedStock - targetStock) / production;
+				}
 			}
+
 			else if (ModelConfig.PRICE_UPDATE_METHOD.equals("MarketImbalance")) {
 				double ratio = stockChange/(production-stockChange);
 				adjustment =  Math.exp(-ratio * ModelConfig.MARKET_LAMBDA);
@@ -139,7 +155,7 @@ public class GlobalPrice implements Serializable {
 			double refPrice = ModelConfig.IS_CALIBRATION_RUN ? exportPrice : referencePrice;  // during calibration reference price isn't fixed, but it is after that
 
 			return new GlobalPrice(thisTimeStep, newPrice, updatedStock, newImports, newExportAmountBeforeLoss,
-					newTransportLosses, refPrice, updatedStockUseRatio, stockChange);
+					newTransportLosses, refPrice, updatedStockUseRatio, stockChange, production);
 		}
 		else {
 			LogWriter.printlnWarning(String.format("Price for not updated (still %s), as no imports and no exports for it", exportPrice));
@@ -150,7 +166,7 @@ public class GlobalPrice implements Serializable {
 	@Override
 	public String toString() {
 		return "GlobalPrice [timestep=" + timestep.getTimestep() + ", exportPrice=" + exportPrice + ", importAmount=" + importAmount
-				+ ", stockLevel=" + stockLevel + ", referencePrice=" + referencePrice + "]";
+				+ ", stockLevel=" + stockLevel + ", referencePrice=" + referencePrice + ", production=" + production + "]";
 	}
 
 	public double getReferencePrice() {
@@ -160,31 +176,34 @@ public class GlobalPrice implements Serializable {
 	public double getStockChange() {
 		return exportAmountBeforeLoss - transportLosses - importAmount;
 	}
+
+	public double getProduction() {
+		return production;
+	}
 	
 	public GlobalPrice createPriceAdjustedByFactor(double factor) {
 		return new GlobalPrice(timestep, exportPrice * factor, stockLevel, importAmount, exportAmountBeforeLoss,
-				transportLosses, referencePrice, stockUseRatio, prevStockChange);
+				transportLosses, referencePrice, stockUseRatio, prevStockChange, production);
 	}
 
 	public void resetStock(Double stock) {
 		this.stockLevel = stock;
 	}
 
-	public boolean isStockCritical() { // is stock projected to go below 0?
-		double stockChange = exportAmountBeforeLoss - transportLosses - importAmount;
-		return stockLevel + stockChange + (stockChange - prevStockChange) <= 0;
+	public double getTradeChangeUp(double maxOfProdOrSupply) {
+		// plus 5% of global production divided by number of countries in case country has 0 supply
+		maxOfProdOrSupply += 0.05 * production / CountryManager.getInstance().getCountryCount();
+		double supplyFraction = maxOfProdOrSupply / (production * 1.05);
+		// ensures stocks don't go below zero even if no exports
+		return (0.95 * stockLevel - (importAmount - getExportsAfterTransportLosses())) * supplyFraction;
 	}
 
-	public double getExportAdjustment(double targetStockUseRatio) {
-		double targetStockLevel = (stockLevel / stockUseRatio) * targetStockUseRatio;
-		double surplus = stockLevel - targetStockLevel;
-		double targetExport = exportAmountBeforeLoss
-				- ModelConfig.EXPORT_ADJUSTMENT_RATE * (exportAmountBeforeLoss - importAmount)
-				- ModelConfig.STOCK_ADJUSTMENT_RATE * surplus;
-		targetExport = Math.max(targetExport, 0); // can't be negative
-		double ratio = exportAmountBeforeLoss > 0 ? targetExport / exportAmountBeforeLoss : 1;
-		return ratio;
+	public double getTradeChangeDown(double maxOfProdOrSupply, double targetStockUseRatio) {
+		// plus 5% of global production divided by number of countries in case country has 0 supply
+		maxOfProdOrSupply += 0.05 * production / CountryManager.getInstance().getCountryCount();
+		double supplyFraction = maxOfProdOrSupply / (production * 1.05);
+		// Allow exports to increase until stocks 50% above target stock level
+		return Math.max(targetStockUseRatio * production * 1.5 - stockLevel, 0) * supplyFraction;
 	}
 
-
 }
\ No newline at end of file
diff --git a/src/ac/ed/lurg/types/CropType.java b/src/ac/ed/lurg/types/CropType.java
index 47a13ab1..2019fb77 100644
--- a/src/ac/ed/lurg/types/CropType.java
+++ b/src/ac/ed/lurg/types/CropType.java
@@ -11,18 +11,18 @@ import ac.ed.lurg.utils.LogWriter;
 
 public enum CropType {
 
-	WHEAT("WheatBarleyOats", "wheat", true, false, 9.5, ModelConfig.INITAL_PRICE_WHEAT, 0.45),
-	MAIZE("MaizeMilletSorghum", "maize", true, false, 5.1, ModelConfig.INITAL_PRICE_MAIZE, 0.35),
-	RICE("Rice (Paddy Equivalent)", "rice", true, false, 8.3, ModelConfig.INITAL_PRICE_RICE, 0.50),
-	OILCROPS("Oilcrops", "oilcrops", true, false, 4.4, ModelConfig.INITAL_PRICE_OILCROPS, 0.25),
-	PULSES("Pulses", "pulses", true, false, 10.8, ModelConfig.INITAL_PRICE_PULSES, 0.45),
-	STARCHY_ROOTS("Starchy Roots", "starchyRoots", true, false, 14.3, ModelConfig.INITAL_PRICE_STARCHYROOTS, 0.10),
-	ENERGY_CROPS("Energy crops", "energycrops", true, false, 5, ModelConfig.INITAL_PRICE_ENERGYCROPS, 0.2),
+	WHEAT("WheatBarleyOats", "wheat", true, false, 9.5, ModelConfig.INITAL_PRICE_WHEAT, 0.3),
+	MAIZE("MaizeMilletSorghum", "maize", true, false, 5.1, ModelConfig.INITAL_PRICE_MAIZE, 0.3),
+	RICE("Rice (Paddy Equivalent)", "rice", true, false, 8.3, ModelConfig.INITAL_PRICE_RICE, 0.3),
+	OILCROPS("Oilcrops", "oilcrops", true, false, 4.4, ModelConfig.INITAL_PRICE_OILCROPS, 0.3),
+	PULSES("Pulses", "pulses", true, false, 10.8, ModelConfig.INITAL_PRICE_PULSES, 0.3),
+	STARCHY_ROOTS("Starchy Roots", "starchyRoots", true, false, 14.3, ModelConfig.INITAL_PRICE_STARCHYROOTS, 0.3),
+	ENERGY_CROPS("Energy crops", "energycrops", true, false, 5, ModelConfig.INITAL_PRICE_ENERGYCROPS, 0.3),
 	SETASIDE("setaside", "setaside", false, false, 0, 0, 0),
-	MONOGASTRICS("Monogastrics", "monogastrics", true, true, 3.1, ModelConfig.INITAL_PRICE_MONOGASTRICS, 0.35),
-	RUMINANTS("Ruminants", "ruminants", true, true, 2.2, ModelConfig.INITAL_PRICE_RUMINANTS, 0.10),
-	FRUITVEG("FruitVeg", "fruitveg", true, false, 8.9, ModelConfig.INITAL_PRICE_FRUITVEG, 0.10),
-	SUGAR("Sugar", "sugar", true, false, 0.1, ModelConfig.INITAL_PRICE_SUGAR, 0.40),
+	MONOGASTRICS("Monogastrics", "monogastrics", true, true, 3.1, ModelConfig.INITAL_PRICE_MONOGASTRICS, 0.3),
+	RUMINANTS("Ruminants", "ruminants", true, true, 2.2, ModelConfig.INITAL_PRICE_RUMINANTS, 0.3),
+	FRUITVEG("FruitVeg", "fruitveg", true, false, 8.9, ModelConfig.INITAL_PRICE_FRUITVEG, 0.3),
+	SUGAR("Sugar", "sugar", true, false, 0.1, ModelConfig.INITAL_PRICE_SUGAR, 0.3),
 	PASTURE("pasture", "pasture", false, false, 0, 0, 0);  // confusing having a land cover and a crop type.  Needed here for yields (but not used for cropland area fractions).
 
 	private String faoName;
-- 
GitLab