Merge pull request #221 from centreformicrosimulation/experimental/lifetime_draws

Experimental/lifetime draws

Author: andrewbaxter439

input/reg_labourSupplyUtility.xlsx

@@ -18,10 +18,7 @@
 import org.apache.commons.math3.distribution.MultivariateNormalDistribution;
 import org.apache.commons.math3.distribution.NormalDistribution;
 import org.apache.commons.math3.util.Pair;
-import simpaths.data.startingpop.DataParser;
-import simpaths.model.AnnuityRates;
 import simpaths.model.decisions.Grids;
-import simpaths.model.enums.*;
 import simpaths.model.taxes.DonorTaxUnit;
 import simpaths.model.taxes.MatchFeature;
 import simpaths.model.taxes.database.TaxDonorDataParser;
@@ -241,6 +238,10 @@ else if(numberOfChildren <= 5) {
     public static final double WEEKS_PER_MONTH = 365.25/(7.*12.);	// = 4.348214286
     public static final double WEEKS_PER_YEAR = 365.25 / 7.;
 
+    // Determine probability of yearly labour supply matches persisting from previous year
+    public static double labour_innovation_employment_persistence_probability = 0.9;
+    public static double labour_innovation_notinemployment_persistence_probability = 0.1;
+
     public static final int HOURS_IN_WEEK = 24 * 7; //This is used to calculate leisure in labour supply
     //Is it possible for people to start going to the labour module (e.g. age 17) while they are living with parents (until age 18)?
     //Cannot see how its possible if it is the household that decides how much labour to supply.  If someone finishes school at 17, they need to leave home before they can enter the labour market.  So set age for finishing school and leaving home to 18.
@@ -277,8 +278,8 @@ else if(numberOfChildren <= 5) {
 
     //public static int MAX_AGE_IN_EDUCATION;// = MAX_AGE;//30;			// Max age a person can stay in education	//Cannot set here, as MAX_AGE is not known yet.  Now set to MAX_AGE in buildObjects in Model class.
     //public static int MAX_AGE_MARRIAGE;// = MAX_AGE;//75;  			// Max age a person can marry		//Cannot set here, as MAX_AGE is not known yet.  Now set to MAX_AGE in buildObjects in Model class.
-    private static int MIN_START_YEAR = 2011; //Minimum allowed starting point. Should correspond to the oldest initial population.
-    private static int MAX_START_YEAR = 2021; //Maximum allowed starting point. Should correspond to the most recent initial population.
+    private static int MIN_START_YEAR = 2015; //Minimum allowed starting point. Should correspond to the oldest initial population.
+    private static int MAX_START_YEAR = 2019; //Maximum allowed starting point. Should correspond to the most recent initial population.
     public static int startYear;
     public static int endYear;
     private static final int MIN_START_YEAR_TESTING = 2019;
@@ -1139,10 +1140,10 @@ else if(country.equals(Country.UK)) {
             columnsEmploymentSelectionFemalesE = 29;
             columnsLabourSupplyUtilityMales = 19;
             columnsLabourSupplyUtilityFemales = 12;
-            columnsLabourSupplyUtilityMalesWithDependent = 23;
-            columnsLabourSupplyUtilityFemalesWithDependent = 23;
-            columnsLabourSupplyUtilityACMales = 17;
-            columnsLabourSupplyUtilityACFemales = 17;
+            columnsLabourSupplyUtilityMalesWithDependent = 15;
+            columnsLabourSupplyUtilityFemalesWithDependent = 15;
+            columnsLabourSupplyUtilityACMales = 31;
+            columnsLabourSupplyUtilityACFemales = 31;
             columnsLabourSupplyUtilityCouples = 64;
             columnsLabourCovid19_SE = 1;
             columnsLabourCovid19_2a_processes = 1;

src/main/java/simpaths/data/Parameters.java

@@ -95,6 +95,8 @@ public class BenefitUnit implements EventListener, IDoubleSource, Weight, Compar
     @Transient ArrayList<Triple<Les_c7_covid, Double, Integer>> covid19MonthlyStateAndGrossIncomeAndWorkHoursTripleFemale = new ArrayList<>(); // This ArrayList stores monthly values of labour market states and gross incomes, to be sampled from by the LabourMarket class, for the female member of the benefit unit
 
     @Transient Innovations innovations;
+    @Transient private Double lastLabourInnov;
+    @Transient private Integer lastYear;
 
     @Transient private Integer yearLocal;
     @Transient private Occupancy occupancyLocal;
@@ -153,7 +155,7 @@ public BenefitUnit(Long id, long seed) {
         key  = new PanelEntityKey(id);        //Sets up key
 
         this.seed = seed;
-        innovations = new Innovations(9, seed);
+        innovations = new Innovations(9, 1, seed);
 
         this.numberChildrenAll_lag1 = 0;
         this.numberChildren02_lag1 = 0;
@@ -1108,7 +1110,28 @@ protected void updateLabourSupplyAndIncome() {
             }
 
             //Sample labour supply from possible labour (pairs of) values
-            double labourInnov = innovations.getDoubleDraw(5);
+            double labourInnov = 0;
+
+            // Check if:
+            // - single occupant is employed
+            // - or male of couple is employed
+            // - or male is not at risk of work and female is employed
+            // -> then persist employment at prob_{persist_employment}
+            // Otherwise
+            // -> persist unemployment at prob_{persist_unemployment}
+
+            if (occupancy.Single_Male.equals(occupancy) && male.atRiskOfWork() && male.getEmployed() == 1) {
+                    labourInnov = getLabourInnovation(Parameters.labour_innovation_employment_persistence_probability);
+            } else if (occupancy.Single_Female.equals(occupancy) && female.atRiskOfWork() && female.getEmployed() == 1) {
+                    labourInnov = getLabourInnovation(Parameters.labour_innovation_employment_persistence_probability);
+            } else if (occupancy.equals(Occupancy.Couple) &&
+                    ((male.atRiskOfWork() && male.getEmployed() == 1) ||
+                            (!male.atRiskOfWork() && female.atRiskOfWork() && female.getEmployed() == 1))) {
+                labourInnov = getLabourInnovation(Parameters.labour_innovation_employment_persistence_probability);
+            } else {
+                labourInnov = getLabourInnovation(Parameters.labour_innovation_notinemployment_persistence_probability);
+            }
+
             try {
                 MultiKeyMap<Labour, Double> labourSupplyUtilityRegressionProbabilitiesByLabourPairs = convertRegressionScoresToProbabilities(labourSupplyUtilityRegressionScoresByLabourPairs);
                 labourSupplyChoice = ManagerRegressions.multiEvent(labourSupplyUtilityRegressionProbabilitiesByLabourPairs, labourInnov);
@@ -1529,6 +1552,30 @@ public enum Regressors {
         MaleEduH_30,
         MaleEduM_40,
         MaleEduH_40,
+        MaleEduM_1,
+        MaleEduH_1,
+        MaleEduM_2,
+        MaleEduH_2,
+        MaleEduM_3,
+        MaleEduH_3,
+        MaleEduM_4,
+        MaleEduH_4,
+        FemaleEduM_10,
+        FemaleEduH_10,
+        FemaleEduM_20,
+        FemaleEduH_20,
+        FemaleEduM_30,
+        FemaleEduH_30,
+        FemaleEduM_40,
+        FemaleEduH_40,
+        FemaleEduM_1,
+        FemaleEduH_1,
+        FemaleEduM_2,
+        FemaleEduH_2,
+        FemaleEduM_3,
+        FemaleEduH_3,
+        FemaleEduM_4,
+        FemaleEduH_4,
         MaleLeisure_dnc,
         FemaleLeisure_dnc,
         MaleLeisure_dnc02,
@@ -2416,6 +2463,78 @@ public double getDoubleValue(Enum<?> variableID) {
             case MaleEduH_40 -> {
                 return (getMale() != null && getMale().getLabourSupplyWeekly().equals(Labour.FORTY) && getMale().getDeh_c3().equals(Education.High)) ? 1. : 0.;
             }
+            case MaleEduM_1 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.TEN) && getMale().getDeh_c3().equals(Education.Medium)) ? 1. : 0.;
+            }
+            case MaleEduH_1 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.TEN) && getMale().getDeh_c3().equals(Education.High)) ? 1. : 0.;
+            }
+            case MaleEduM_2 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.TWENTY) && getMale().getDeh_c3().equals(Education.Medium)) ? 1. : 0.;
+            }
+            case MaleEduH_2 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.TWENTY) && getMale().getDeh_c3().equals(Education.High)) ? 1. : 0.;
+            }
+            case MaleEduM_3 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.THIRTY) && getMale().getDeh_c3().equals(Education.Medium)) ? 1. : 0.;
+            }
+            case MaleEduH_3 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.THIRTY) && getMale().getDeh_c3().equals(Education.High)) ? 1. : 0.;
+            }
+            case MaleEduM_4 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.FORTY) && getMale().getDeh_c3().equals(Education.Medium)) ? 1. : 0.;
+            }
+            case MaleEduH_4 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.FORTY) && getMale().getDeh_c3().equals(Education.High)) ? 1. : 0.;
+            }
+            case FemaleEduM_10 -> {
+                return (getMale() != null && getFemale() != null && getMale().getLabourSupplyWeekly().equals(Labour.TEN) && getFemale().getDeh_c3().equals(Education.Medium)) ? 1. : 0.;
+            }
+            case FemaleEduH_10 -> {
+                return (getMale() != null && getFemale() != null && getMale().getLabourSupplyWeekly().equals(Labour.TEN) && getFemale().getDeh_c3().equals(Education.High)) ? 1. : 0.;
+            }
+            case FemaleEduM_20 -> {
+                return (getMale() != null && getFemale() != null && getMale().getLabourSupplyWeekly().equals(Labour.TWENTY) && getFemale().getDeh_c3().equals(Education.Medium)) ? 1. : 0.;
+            }
+            case FemaleEduH_20 -> {
+                return (getMale() != null && getFemale() != null && getMale().getLabourSupplyWeekly().equals(Labour.TWENTY) && getFemale().getDeh_c3().equals(Education.High)) ? 1. : 0.;
+            }
+            case FemaleEduM_30 -> {
+                return (getMale() != null && getFemale() != null && getMale().getLabourSupplyWeekly().equals(Labour.THIRTY) && getFemale().getDeh_c3().equals(Education.Medium)) ? 1. : 0.;
+            }
+            case FemaleEduH_30 -> {
+                return (getMale() != null && getFemale() != null && getMale().getLabourSupplyWeekly().equals(Labour.THIRTY) && getFemale().getDeh_c3().equals(Education.High)) ? 1. : 0.;
+            }
+            case FemaleEduM_40 -> {
+                return (getMale() != null && getFemale() != null && getMale().getLabourSupplyWeekly().equals(Labour.FORTY) && getFemale().getDeh_c3().equals(Education.Medium)) ? 1. : 0.;
+            }
+            case FemaleEduH_40 -> {
+                return (getMale() != null && getFemale() != null && getMale().getLabourSupplyWeekly().equals(Labour.FORTY) && getFemale().getDeh_c3().equals(Education.High)) ? 1. : 0.;
+            }
+            case FemaleEduM_1 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.TEN) && getFemale().getDeh_c3().equals(Education.Medium)) ? 1. : 0.;
+            }
+            case FemaleEduH_1 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.TEN) && getFemale().getDeh_c3().equals(Education.High)) ? 1. : 0.;
+            }
+            case FemaleEduM_2 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.TWENTY) && getFemale().getDeh_c3().equals(Education.Medium)) ? 1. : 0.;
+            }
+            case FemaleEduH_2 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.TWENTY) && getFemale().getDeh_c3().equals(Education.High)) ? 1. : 0.;
+            }
+            case FemaleEduM_3 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.THIRTY) && getFemale().getDeh_c3().equals(Education.Medium)) ? 1. : 0.;
+            }
+            case FemaleEduH_3 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.THIRTY) && getFemale().getDeh_c3().equals(Education.High)) ? 1. : 0.;
+            }
+            case FemaleEduM_4 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.FORTY) && getFemale().getDeh_c3().equals(Education.Medium)) ? 1. : 0.;
+            }
+            case FemaleEduH_4 -> {
+                return (getMale() != null  && getFemale() != null && getFemale().getLabourSupplyWeekly().equals(Labour.FORTY) && getFemale().getDeh_c3().equals(Education.High)) ? 1. : 0.;
+            }
             case MaleLeisure_dnc02 -> {
                 return (Parameters.HOURS_IN_WEEK - getMale().getLabourSupplyHoursWeekly()) * getIndicatorChildren(0,1).ordinal();
             }
@@ -3646,4 +3765,55 @@ public void removeMember(Person member) {
     }
 
     public static void setBenefitUnitIdCounter(long id) {benefitUnitIdCounter = id;}
+
+        /**
+         * Calculates the labour innovation value for the current period with persistence.
+         *
+         * This method implements a persistence mechanism for labour innovations where there's
+         * a probability that the previous period's innovation will persist into the current period.
+         * A new innovation is drawn each period, but it's only used with probability (1 - persistenceProbability).
+         *
+         * The persistence mechanism works as follows:
+         * 1. Draw a new potential innovation value
+         * 2. Draw a random number to determine if we should use the new value or persist with the old one
+         * 3. If it's the first period or there was a gap in years, use the new innovation
+         * 4. Otherwise, use the persistence probability to decide between new and old values
+         *
+         * @param persistenceProbability The probability (between 0 and 1) that the previous period's
+         *                              innovation will persist. A value of 0 means always use new innovations,
+         *                              while 1 means always keep the old innovation.
+         * @return The labour innovation value for this period
+         */
+    private double getLabourInnovation(double persistenceProbability) {
+
+        double newLabourInnovation = innovations.getDoubleDraw(5);
+
+        // persistenceRandomDraw - random value to determine if we keep old innovation
+        double persistenceRandomDraw = innovations.getDoubleDraw(6);
+
+        int currentYear = model.getYear();
+
+        // Initialize on first call
+        if (lastLabourInnov == null || lastYear == null) {
+            lastLabourInnov = newLabourInnovation;
+            lastYear = currentYear;
+            return newLabourInnovation;
+        }
+
+        double labourInnov;
+        if (persistenceRandomDraw < persistenceProbability) {
+            labourInnov = lastLabourInnov;
+        } else {
+            labourInnov = newLabourInnovation;
+        }
+
+        // Store values for next period
+        lastLabourInnov = labourInnov;
+        lastYear = currentYear;
+
+        return labourInnov;
+
+
+
+        }
 }

src/main/java/simpaths/model/BenefitUnit.java

@@ -28,6 +28,14 @@ public Innovations(int nDoubleInnovs, int nSingleDrawDoubleInnovs, int nSingleDr
         }
     }
 
+    public Innovations(int nDoubleInnovs, int nSingleDrawDoubleInnovs, long seed) {
+        this(nDoubleInnovs, seed);
+        singleDrawDoubleInnovs = new double[nSingleDrawDoubleInnovs];
+        for (int ii = 0; ii < nSingleDrawDoubleInnovs; ii++) {
+            singleDrawDoubleInnovs[ii] = generator.nextDouble();
+        }
+    }
+
     public void getNewDoubleDraws() {
         for (int ii = 0; ii < doubleInnovs.length; ii++) {
             doubleInnovs[ii] = generator.nextDouble();