def elasticityOfSubstitution (ρ : ℝ) : ℝ := 1 / (1 - ρ)

theorem sigma_rho_inverse {ρ : ℝ} (hρ : ρ ≠ 1) :
    ρ = 1 - 1 / elasticityOfSubstitution ρ := by
  simp only [elasticityOfSubstitution]
  have h : (1 : ℝ) - ρ ≠ 0 := by intro h; apply hρ; linarith
  field_simp
  linarith

theorem factorShare_sum_eq_one {a : Fin J → ℝ} {ρ : ℝ} {x : Fin J → ℝ}
    (hden : 0 < ∑ i : Fin J, a i * (x i) ^ ρ) :
    ∑ j : Fin J, factorShare J a ρ x j = 1 := by
  simp only [factorShare]
  simp_rw [div_eq_mul_inv]
  rw [← Finset.sum_mul]
  exact mul_inv_cancel₀ (ne_of_gt hden)

theorem shephards_lemma {J : ℕ} (hJ : 0 < J) {ρ : ℝ} (hρ : ρ ≠ 0) (hρ1 : ρ ≠ 1)
    {p : Fin J → ℝ} (hp : ∀ j, 0 < p j) (k : Fin J) :
    HasDerivAt (fun t => cesUnitCost J ρ (Function.update p k t))
               (cesDemand J ρ p k 1) (p k) := by
  simp only [cesUnitCost, cesPriceIndex, cesDemand, unnormCES, dualExponent]
  set r := ρ / (ρ - 1)
  have hr : r ≠ 0 := div_ne_zero hρ (sub_ne_zero.mpr hρ1)
  set S := ∑ j : Fin J, (p j) ^ r
  set C := ∑ j ∈ Finset.univ.erase k, (p j) ^ r
  have hSpk : (p k) ^ r + C = S := by
    simp only [S, C]
    have h := Finset.sum_erase_add Finset.univ (fun j => (p j) ^ r) (Finset.mem_univ k)
    linarith
  have hSpos : 0 < S := by
    simp only [S]
    exact Finset.sum_pos (fun j _ => rpow_pos_of_pos (hp j) r)
      ⟨⟨0, hJ⟩, Finset.mem_univ _⟩
  -- Rewrite the function: updating p k to t gives inner sum = t^r + C
  have hfun_eq : (fun t : ℝ =>
      (∑ j : Fin J, (Function.update p k t j) ^ r) ^ (1 / r)) =
      (fun t => (t ^ r + C) ^ (1 / r)) := funext (fun t => by
    congr 1
    conv_lhs => rw [← Finset.insert_erase (Finset.mem_univ k)]
    rw [Finset.sum_insert (Finset.notMem_erase k _)]
    simp only [Function.update_self]; congr 1
    apply Finset.sum_congr rfl; intro j hj
    simp [Function.update_of_ne ((Finset.mem_erase.mp hj).1)])
  rw [hfun_eq]
  -- Inner derivative: d/dt [t^r + C] = r * t^{r-1}
  have hAd : HasDerivAt (fun t : ℝ => t ^ r + C)
      (r * (p k) ^ (r - 1)) (p k) := by
    have h := (hasDerivAt_id (p k)).rpow_const (p := r) (Or.inl (hp k).ne')
    simp only [id, one_mul] at h
    exact h.add_const C
  -- Chain rule: d/dt [(t^r + C)^{1/r}]
  have h_rpow := hAd.rpow_const (p := 1 / r) (Or.inl (hSpk ▸ hSpos).ne')
  have hOuter := hSpk ▸ h_rpow
  -- Algebraic identity: chain rule derivative = cesDemand form
  have hkey : r * (p k) ^ (r - 1) * (1 / r) * S ^ (1 / r - 1) =
      1 * (p k) ^ (r - 1) * (S ^ (1 / r)) ^ (1 - r) := by
    rw [← Real.rpow_mul (le_of_lt hSpos)]
    have h1 : 1 / r * (1 - r) = 1 / r - 1 := by field_simp [hr]
    rw [h1]
    have h2 : r * (p k) ^ (r - 1) * (1 / r) = 1 * (p k) ^ (r - 1) := by
      field_simp [hr]
    rw [h2]
  exact hkey ▸ hOuter

theorem ces_limit_cobbDouglas {J : ℕ} (hJ : 0 < J) {x : Fin J → ℝ} (hx : ∀ j, 0 < x j) :
    Filter.Tendsto
      (fun ρ : ℝ => ((1 / (↑J : ℝ)) * ∑ j : Fin J, (x j) ^ ρ) ^ (1 / ρ))
      (nhdsWithin 0 (Set.Ioi 0))
      (nhds (geometricMean J x)) := by
  set A : ℝ → ℝ := fun ρ => (1 / (↑J : ℝ)) * ∑ j : Fin J, (x j) ^ ρ with hA_def
  set D := (1 / (↑J : ℝ)) * ∑ j : Fin J, Real.log (x j)
  have hJpos : (0 : ℝ) < ↑J := Nat.cast_pos.mpr hJ
  have hAd := ces_inner_hasDerivAt hJ hx
  have hA0 : A 0 = 1 := by simp [A, rpow_zero, hJ.ne']
  have hApos : ∀ ρ, 0 < A ρ := fun ρ =>
    mul_pos (div_pos one_pos hJpos)
      (Finset.sum_pos (fun j _ => rpow_pos_of_pos (hx j) ρ) ⟨⟨0, hJ⟩, Finset.mem_univ _⟩)
  -- log A(ρ) has derivative D at 0 (chain rule: (log A)'(0) = A'(0)/A(0) = D/1 = D)
  have hlogAd : HasDerivAt (fun ρ => Real.log (A ρ)) D 0 := by
    have h := (hasDerivAt_log (hApos 0).ne').comp 0 hAd
    simp only [hA0, inv_one, one_mul] at h
    exact h
  have hlogA0 : Real.log (A 0) = 0 := by rw [hA0]; exact log_one
  -- log A(ρ)/ρ → D as ρ → 0+ (this is the "0/0" limit resolved by the derivative)
  have hslope : Filter.Tendsto (fun ρ => Real.log (A ρ) / ρ)
      (nhdsWithin 0 (Set.Ioi 0)) (nhds D) := by
    have hts := hasDerivAt_iff_tendsto_slope.mp hlogAd
    have hmono : nhdsWithin (0 : ℝ) (Set.Ioi 0) ≤ nhdsWithin 0 {0}ᶜ :=
      nhdsWithin_mono 0 (fun _ hρ => hρ.ne')
    exact (hts.mono_left hmono).congr'
      (eventually_nhdsWithin_of_forall (fun ρ _ => by
        simp only [slope_def_field, hlogA0, sub_zero]))
  -- exp(D) = geometricMean J x
  have hexpD : Real.exp D = geometricMean J x := by
    simp only [geometricMean]
    rw [show D = (1 / (↑J : ℝ)) * ∑ j : Fin J, Real.log (x j) from rfl,
        ← Real.log_prod (fun j _ => (hx j).ne'),
        Real.rpow_def_of_pos (Finset.prod_pos (fun j _ => hx j))]
    congr 1; ring
  -- Conclude: A(ρ)^{1/ρ} = exp(log A(ρ)/ρ) → exp(D) = geometricMean
  rw [← hexpD]
  exact ((Real.continuous_exp.continuousAt.tendsto).comp hslope).congr'
    (eventually_nhdsWithin_of_forall (fun ρ _ => by
      simp only [Function.comp, A, rpow_def_of_pos (hApos _)]
      congr 1; ring))

theorem logit_is_ces_at_q_one (T : ℝ) (_hT : T ≠ 0) (ε : Fin J → ℝ) (j : Fin J) :
    logitProbability J T ε j =
    escortDistribution J 1 (fun k => Real.exp (ε k / T)) j := by
  simp only [logitProbability, escortDistribution, rpow_one]

theorem fisherInfoRho_zero_at_symmetry [NeZero J]
    {c : ℝ} (hc : 0 < c) (ρ : ℝ) :
    fisherInfoRho (fun _ : Fin J => c) ρ = 0 :=
  escort_fisher_zero_at_symmetry hc ρ

theorem cramerRao_diverges_at_symmetry [NeZero J]
    {c : ℝ} (hc : 0 < c) (ρ : ℝ) (N : ℕ) (_hN : 0 < N) :
    cramerRaoBound N (fun _ : Fin J => c) ρ =
    1 / (↑N * 0) := by
  simp only [cramerRaoBound, fisherInfoRho,
             escort_fisher_zero_at_symmetry hc]