1"""Balanced-growth path of the economy with frozen AI (alpha fixed): the reference state.23With Harrod-neutral efficiency growth g on labor, land and the energy resource, all stocks, output and4consumption grow at g, prices per efficiency unit are constant, and the interest rate is5r* = (1+g)^gamma / beta - 1 (Euler with a constant price index). Each capital stock satisfies the6steady-state version of its Euler equation (investment rate x = delta + g):7 R_k = q_k (r* + delta_k) - p_k phi_k g (delta_k + g/2), q_k = p_k (1 + phi_k g).8"""9from __future__ import annotations1011import numpy as np12from scipy.optimize import fsolve1314from .production import housing_rent, static_period15from .tasks import technology_from_params16from .transition import KINDS171819def bgp_interest_rate(p):20 return (1.0 + p.g_Z) ** p.gamma_u / p.beta - 1.0212223def target_rentals(p, r, pQ=None):24 g = p.g_Z25 price = {"N": 1.0, "Q": p.pQ0 if pQ is None else pQ, "E": 1.0, "H": 1.0}26 R = {}27 for k in KINDS:28 d, ph, pk = p.delta[k], p.phi[k], price[k]29 q = pk * (1 + ph * g)30 R[k] = q * (r + d) - pk * ph * g * (d + g / 2)31 return R323334def solve_bgp(p, alpha=None, eps=None, K0=None, KQ=None, verbose=False):35 """Stationary (detrended) stocks K_N, K_E, K_H such that every capital Euler holds at r*, holding the36 compute stock at KQ (default p.KQ0) and solving the compute Euler for the price p^Q_0 in closed form:37 p^Q_0 = R_Q / [ (1+phi g)(r*+delta_Q) - phi g (delta_Q + g/2) ].38 Returns the stocks, r*, the implied pQ0 and the period dictionary (per efficiency unit, gf = 1).39 """40 tech = technology_from_params(p)41 alpha = p.alpha_path()[0] if alpha is None else alpha42 eps = p.eps0 if eps is None else eps43 KQ = p.KQ0 if KQ is None else KQ44 r = bgp_interest_rate(p)45 Rstar = target_rentals(p, r, 1.0)46 g = p.g_Z47 x = {k: p.delta[k] + g for k in KINDS}48 kk = ("N", "E", "H")4950 def period(K, pQ):51 s = static_period(p, tech, K["N"], K["Q"], K["E"], K["H"], alpha, p.Z0, eps, gf=1.0)52 inv = sum((pQ if k == "Q" else 1.0) * (x[k] * K[k] + 0.5 * p.phi[k] * g ** 2 * K[k]) for k in KINDS)53 C = s["Y"] - inv54 housing_rent(p, s, C)55 return s, C5657 def pQ_of(s):58 d, ph = p.delta["Q"], p.phi["Q"]59 return s["RQ"] / ((1 + ph * g) * (r + d) - ph * g * (d + g / 2))6061 def resid(logK):62 K = dict(zip(kk, np.exp(logK)))63 K["Q"] = KQ64 s, C = period(K, 1.0)65 pQ = max(pQ_of(s), 1e-9)66 s, C = period(K, pQ)67 Rk = {"N": s["RN"], "E": s["RE"], "H": s["RHK"]}68 return np.array([np.log(max(Rk[k], 1e-300) / Rstar[k]) for k in kk])6970 K0 = {"N": p.KN0, "E": p.KE0, "H": p.KH0} if K0 is None else K071 lk0 = np.log([K0[k] for k in kk])72 sol, info, ier, msg = fsolve(resid, lk0, full_output=True, xtol=1e-13)73 K = dict(zip(kk, np.exp(sol)))74 K["Q"] = KQ75 s, C = period(K, 1.0)76 pQ = pQ_of(s)77 s, C = period(K, pQ)78 res = resid(sol)79 if verbose:80 print("bgp fsolve:", msg.strip(), "max|res|", np.abs(res).max(), "pQ0", pQ)81 return dict(K=K, r=r, pQ0=pQ, period=s, C=C, resid=res, ok=(ier == 1 and np.abs(res).max() < 1e-8))828384def apply_bgp(p, b):85 """Return a copy of p whose initial stocks and compute price are the BGP values."""86 return p.replace(KN0=b["K"]["N"], KQ0=b["K"]["Q"], KE0=b["K"]["E"], KH0=b["K"]["H"], pQ0=b["pQ0"])87