Bug in computing the pre-time-path population

[jdebacker]

The ogcore.demographics module computes the population over the time path, but also needs the population distribution in the year before the time path begins. This is inferred and there seems to be an error constructing this when using certain population distributions (e.g., that from Ethiopia).

@claude is suggesting the following issues and fixes:
I've found the root cause of the infeasible growth rates. Here's
the complete analysis:

Root Cause: Backward Population Solver Produces Negative Values

The Problem:
When ogcore.demographics calculates the pre-period population (2023) by
solving backwards from 2024, it encounters a mathematical impossibility
with Ethiopia's demographic structure, particularly at very old ages
(93-99).

Why This Happens

The demographic equation is:
pop_current[age] = pop_previous[age-1] * survival_rate + pop_previous[age]

• imm_rate

When solving backwards for pop_previous[age]:
pop_previous[age] = (pop_current[age] - survived_from_previous_age) /
imm_rate

The bug occurs when:

• More people survive from age N-1 than exist at age N in the current
  period
• Immigration rate is positive
• Result: pop_previous[age] becomes hugely negative

Specific Examples from Ethiopia

At age 99:

• Current population (2024): 189 people
• People surviving from age 98: 365,692,350 (366 million!)
• Immigration rate: +19.9% (positive)
• Solved pre-population: -1.8 BILLION ❌

This creates a mathematical impossibility. The equation needs a negative
population to make the immigration math work, but negative populations
don't exist.

The Impact

This causes:

1. 76 out of 80 working-age groups to have negative omega_S_preTP values
2. omega_S_preTP ranges from -0.59 to +1.49 (should be all positive,
   summing to 1.0)
3. Growth rate g_n[0] = -105.3% (< -100%, which is impossible)

Why Ethiopia but not South Africa?

Ethiopia has:

• Very young, rapidly growing population (high fertility)
• Steep population decline at old ages
• The backward solver hits this edge case

South Africa has:

• More balanced age distribution
• Slower population changes
• The backward solver works fine

Recommendations

demographics.py:1006-1033 - The pre_pop_solve() function needs to be
fixed. Options:

1. Use actual historical data instead of solving backwards (best solution)
2. Add constraints to prevent negative population values in the solver
3. Use different immigration assumptions for backward vs forward solving
4. Skip the backward solve and use forward data from get_pop() which
   already retrieves 2023 data

The good news: The UN data itself is perfectly fine. This is purely an
ogcore.demographics implementation issue in how it initializes the
pre-transition period population.

cc @rickecon

[jdebacker]

Proposed patch to ogcore.demographics -- don't solve for pre-TP population, just pull it from the UN data:

--- a/ogcore/demographics.py
+++ b/ogcore/demographics.py
@@ -990,56 +990,35 @@ def get_pop_objs(
     # Check that counterfactual pop dist is close to pop dist given
     # assert np.allclose(pop_counter_2D, pop_dist)
     assert np.allclose(pop_counter_2D, pop_2D)
+
+    # Get percentage distribution for S periods for pre-TP period
+    # Use the actual UN historical data instead of solving backwards
+    pre_pop_EpS = pop_rebin(pre_pop, E + S)

-    """"
-    CHANGE - in OG-Core, we are implicitLy assuming pre-TP rates of mortality,
-    fertility, and immigration are the same as the period 0 rates.
-
-    So let's just infer the pre-pop_dist from those.
-    """
-    pop1 = pop_2D[0, :]
-    fert0 = fert_rates[0, :]
-    mort0 = mort_rates[0, :]
-    infmort0 = infmort_rates[0]
-    imm0 = imm_rates_orig[0, :]
-    pre_pop_guess = pop1.copy()
-
-    # I can't solve this analytically, so set up a system of equation
-    # to solve
-    def pre_pop_solve(pre_pop_guess, pop1, fert0, mort0, infmort0, imm0):
-        pre_pop = pre_pop_guess
-        errors = np.zeros(E + S)
-        errors[0] = pop1[0] - (
-            (1 - infmort0) * (fert0 * pre_pop).sum() + imm0[0] * pre_pop[0]
-        )
-        errors[1:] = pop1[1:] - (
-            pre_pop[:-1] * (1 - mort0[:-1]) + pre_pop[1:] * imm0[1:]
-        )
-        # print("Max error = ", np.abs(errors).max())
-        return errors
-
-    opt_res = opt.root(
-        pre_pop_solve,
-        pre_pop_guess,
-        args=(pop1, fert0, mort0, infmort0, imm0),
-        method="lm",
-    )
-    pre_pop = opt_res.x
-    print(
-        "Success? ",
-        opt_res.success,
-        ", Max diff = ",
-        np.abs(opt_res.fun).max(),
-    )
-    pre_pop_EpS = pop_rebin(pre_pop, E + S)
-
-    # Check result
+    # Verify that the UN pre-period data is reasonably consistent
+    # with the period 0 population using the demographic transition equations
     initial_pop_counter = np.zeros(E + S)
-    newborns = (fert_rates[0, :] * pre_pop[:]).sum()
+    newborns = (fert_rates[0, :] * pre_pop_EpS[:]).sum()
     initial_pop_counter[0] = (
         1 - infmort_rates[0]
-    ) * newborns + imm_rates_orig[0, 0] * pre_pop[0]
+    ) * newborns + imm_rates_orig[0, 0] * pre_pop_EpS[0]
     initial_pop_counter[1:] = (
-        pre_pop[:-1] * (1 - mort_rates[0, :-1])
-        + pre_pop[1:] * imm_rates_orig[0, 1:]
+        pre_pop_EpS[:-1] * (1 - mort_rates[0, :-1])
+        + pre_pop_EpS[1:] * imm_rates_orig[0, 1:]
     )
-    # Test that using pre pop get to pop in period 1
-    print("Max diff = ", np.abs(pop_2D[0, :] - initial_pop_counter).max())
-    # assert np.allclose(initial_pop_counter, pop_2D[0, :])
+
+    max_diff = np.abs(pop_2D[0, :] - initial_pop_counter).max()
+    print(
+        "Pre-period population verification: Max diff = ",
+        max_diff
+    )
+
+    if max_diff > 100000:
+        print(
+            "WARNING: Large difference between UN pre-period population "
+            + "and period 0 population ({:.2f}). ".format(max_diff)
+            + "This may indicate inconsistencies in the data or "
+            + "immigration rate calculations, but using UN historical "
+            + "data as it is more reliable than backward-solved estimates."
+        )

     # Create the transition matrix for the population distribution
     # from T0 going forward (i.e., past when we have data on forecasts)