Lecture 6: Equities — Stock Valuation, the Dividend Discount Model, and Growth Opportunities

# ============================================================
# Setup
# ============================================================
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from scipy.optimize import brentq

plt.rcParams.update({
    'figure.figsize': (10, 5),
    'axes.spines.top': False,
    'axes.spines.right': False,
    'axes.grid': True,
    'grid.alpha': 0.3,
    'font.size': 12,
    'lines.linewidth': 2,
})

print("Libraries loaded successfully.")

# ============================================================
# Equity Market Scale — IPOs and Market Activity
# ============================================================
# Lo's slide data: IPO volume and total underwriting 1990-2004
years = list(range(1990, 2005))
ipo_volume = [5, 16, 24, 41, 28, 30, 50, 43, 37, 64, 76, 36, 26, 16, 48]  # $B
total_underwriting = [312, 587, 856, 1063, 716, 722, 979, 1317, 1868, 1960, 1851, 2535, 2581, 2890, 2859]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# IPOs
colors_ipo = ['#e74c3c' if v > 60 else '#3498db' for v in ipo_volume]
ax1.bar(years, ipo_volume, color=colors_ipo, alpha=0.8, edgecolor='white')
ax1.set_xlabel('Year', fontsize=12)
ax1.set_ylabel('IPO Volume ($B)', fontsize=12)
ax1.set_title('Initial Public Offerings\n(US Market)', fontsize=14, fontweight='bold')
ax1.tick_params(axis='x', rotation=45)

# Highlight dot-com boom
ax1.annotate('Dot-com\nboom', xy=(2000, 76), xytext=(1997, 80),
            fontsize=10, color='#e74c3c', fontweight='bold',
            arrowprops=dict(arrowstyle='->', color='#e74c3c'))

# Total underwriting
ax2.bar(years, [v/1000 for v in total_underwriting], color='#2c3e50', alpha=0.8, edgecolor='white')
ax2.set_xlabel('Year', fontsize=12)
ax2.set_ylabel('Total Underwriting ($T)', fontsize=12)
ax2.set_title('Total Securities Underwriting\n(US Market)', fontsize=14, fontweight='bold')
ax2.tick_params(axis='x', rotation=45)

plt.tight_layout()
plt.show()

# ============================================================
# One-Period vs. Multi-Period DDM — Equivalence Demonstration
# ============================================================

def ddm_general(dividends, terminal_price, r):
    """Value a stock from a finite list of dividends plus terminal price."""
    T = len(dividends)
    pv = sum(D / (1 + r)**t for t, D in enumerate(dividends, 1))
    pv += terminal_price / (1 + r)**T
    return pv


# Example: constant dividend stream, showing equivalence
r = 0.12
D = 2.00  # constant dividend forever

# Investor 1: holds for 1 year
# P1 = D/r (perpetuity from year 2 onward)
P1 = D / r
P0_1period = ddm_general([D], P1, r)

# Investor 2: holds for 5 years
P5 = D / r  # same perpetuity from year 6
P0_5period = ddm_general([D]*5, P5, r)

# Investor 3: holds forever (perpetuity formula)
P0_forever = D / r

print("=" * 60)
print("DDM EQUIVALENCE: Holding Period Doesn't Matter")
print(f"Constant dividend D = ${D:.2f}, discount rate r = {r:.0%}")
print("=" * 60)
print(f"1-year holding period:  P₀ = (D₁ + P₁)/(1+r) = ({D} + {P1:.2f})/1.12 = ${P0_1period:.4f}")
print(f"5-year holding period:  P₀ = Σ PV(D) + PV(P₅)                       = ${P0_5period:.4f}")
print(f"Infinite holding:       P₀ = D/r = {D}/{r}                            = ${P0_forever:.4f}")
print(f"\n→ All methods give P₀ = ${D/r:.2f} regardless of holding period ✓")
print(f"  The stock price today reflects ALL future dividends, not just")
print(f"  the dividends you personally will receive while holding it.")

# ============================================================
# DDM Convergence: How Much Value Comes From Near vs. Far Dividends?
# ============================================================
# ▶ MODIFY AND RE-RUN
D = 2.00
r = 0.10
g = 0.04  # growth rate
# ============================================================

# Gordon Growth Model value
P0 = D * (1 + g) / (r - g)

# Cumulative PV of dividends from year 1 to T
T_max = 100
cumulative_pv = np.zeros(T_max)
for t in range(1, T_max + 1):
    Dt = D * (1 + g)**t
    cumulative_pv[t-1] = cumulative_pv[t-2] + Dt / (1 + r)**t if t > 1 else Dt / (1 + r)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# Left: Individual PV contributions
T_show = 50
pvs = [D * (1 + g)**t / (1 + r)**t for t in range(1, T_show + 1)]
ax1.bar(range(1, T_show + 1), pvs, color='#3498db', alpha=0.7, width=0.8)
ax1.set_xlabel('Year', fontsize=12)
ax1.set_ylabel('PV of Dividend ($)', fontsize=12)
ax1.set_title(f'PV of Each Year\'s Dividend\n(D₀=${D}, g={g:.0%}, r={r:.0%})', fontsize=13, fontweight='bold')

# Right: Cumulative
years = np.arange(1, T_max + 1)
ax2.plot(years, cumulative_pv, color='#2c3e50', linewidth=3)
ax2.axhline(y=P0, color='#e74c3c', linewidth=2, linestyle='--', label=f'Gordon Model: P₀ = ${P0:.2f}')

# Milestones
for pct in [0.5, 0.75, 0.9]:
    idx = np.searchsorted(cumulative_pv, P0 * pct)
    if idx < T_max:
        ax2.plot(idx + 1, cumulative_pv[idx], 'o', color='#e67e22', markersize=8, zorder=5)
        ax2.annotate(f'{pct:.0%} in ~{idx+1} years', xy=(idx + 1, cumulative_pv[idx]),
                    xytext=(idx + 10, cumulative_pv[idx] - P0 * 0.05), fontsize=9,
                    arrowprops=dict(arrowstyle='->', color='#e67e22'))

ax2.set_xlabel('Years of Dividends Included', fontsize=12)
ax2.set_ylabel('Cumulative PV ($)', fontsize=12)
ax2.set_title('Convergence to Full DDM Value', fontsize=13, fontweight='bold')
ax2.legend(fontsize=11)

plt.tight_layout()
plt.show()

for pct in [0.5, 0.75, 0.9, 0.95]:
    idx = np.searchsorted(cumulative_pv, P0 * pct) + 1
    print(f"  {pct:.0%} of value captured in first {idx:>3d} years")

# ============================================================
# Gordon Growth Model — Core Calculations
# ============================================================

def gordon_growth_price(D0, r, g):
    """Gordon Growth Model: P0 = D0*(1+g)/(r-g). Requires r > g."""
    if r <= g:
        return float('inf')
    return D0 * (1 + g) / (r - g)


# Lo's example
D0 = 1.00
r = 0.20
g = 0.06
P0 = gordon_growth_price(D0, r, g)

print("=" * 60)
print("GORDON GROWTH MODEL — Lo's Example")
print("=" * 60)
print(f"Current dividend D₀ = ${D0:.2f}")
print(f"Growth rate g = {g:.0%}")
print(f"Required return r = {r:.0%}")
print(f"\nP₀ = D₀(1+g)/(r-g) = ${D0}×{1+g:.2f}/({r:.2f}-{g:.2f}) = ${P0:.2f}")

# Sensitivity table
print(f"\n{'':>8s}", end='')
for g_val in [0.02, 0.04, 0.06, 0.08, 0.10, 0.12]:
    print(f"{'g='+f'{g_val:.0%}':>9s}", end='')
print()
print("-" * 62)

for r_val in [0.12, 0.14, 0.16, 0.18, 0.20, 0.25]:
    print(f"r={r_val:.0%}  ", end='')
    for g_val in [0.02, 0.04, 0.06, 0.08, 0.10, 0.12]:
        if r_val > g_val:
            p = gordon_growth_price(1.0, r_val, g_val)
            print(f" ${p:>6.2f}", end='')
        else:
            print(f"     ∞ ", end='')
    print()

print("\nKey insight: Price is EXTREMELY sensitive to (r-g).")
print("A small change in growth or discount rate has a massive impact on valuation.")

# ============================================================
# Sensitivity Analysis: Price vs. Growth Rate and Discount Rate
# ============================================================
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# Left: Price vs. growth rate (fixed r)
g_range = np.linspace(0, 0.095, 200)
for r_val, color, ls in [(0.10, '#e74c3c', '-'), (0.12, '#e67e22', '--'),
                          (0.15, '#3498db', '-'), (0.20, '#2c3e50', '--')]:
    prices = [gordon_growth_price(1.0, r_val, g) for g in g_range]
    ax1.plot(g_range * 100, prices, color=color, linewidth=2, linestyle=ls,
            label=f'r = {r_val:.0%}')

ax1.set_xlabel('Growth Rate g (%)', fontsize=12)
ax1.set_ylabel('Stock Price ($)', fontsize=12)
ax1.set_title('Price vs. Growth Rate\n(D₀ = $1.00)', fontsize=14, fontweight='bold')
ax1.set_ylim(0, 60)
ax1.legend(fontsize=10)
ax1.annotate('Price explodes\nas g → r', xy=(8.5, 50), fontsize=10, color='#e74c3c',
            style='italic', bbox=dict(boxstyle='round', facecolor='#fdf2f2', alpha=0.9))

# Right: Price vs. discount rate (fixed g)
r_range = np.linspace(0.05, 0.25, 200)
for g_val, color, ls in [(0.02, '#2ecc71', '-'), (0.04, '#3498db', '--'),
                          (0.06, '#e67e22', '-'), (0.08, '#e74c3c', '--')]:
    prices = [gordon_growth_price(1.0, r, g_val) if r > g_val else np.nan for r in r_range]
    ax2.plot(r_range * 100, prices, color=color, linewidth=2, linestyle=ls,
            label=f'g = {g_val:.0%}')

ax2.set_xlabel('Required Return r (%)', fontsize=12)
ax2.set_ylabel('Stock Price ($)', fontsize=12)
ax2.set_title('Price vs. Discount Rate\n(D₀ = $1.00)', fontsize=14, fontweight='bold')
ax2.set_ylim(0, 60)
ax2.legend(fontsize=10)

plt.tight_layout()
plt.show()

print("The Gordon model has a singularity at g = r (price → ∞).")
print("In practice, g must be STRICTLY less than r for the model to be valid.")
print("For mature companies, g should not exceed nominal GDP growth (roughly 4-6%).")

# ============================================================
# Cost of Equity Estimation — Duke Power Example
# ============================================================
D0_P0 = 0.052  # current dividend yield

growth_estimates = {
    'Value Line': 0.050,
    'IBES Consensus': 0.055,
    'Historical': 0.041,
}

print("=" * 60)
print("COST OF EQUITY — Duke Power (September 1992)")
print("=" * 60)
print(f"Current dividend yield D₀/P₀ = {D0_P0:.1%}")
print()
print(f"{'Source':>18s} {'g':>8s} {'D₁/P₀':>10s} {'r = D₁/P₀ + g':>15s}")
print("-" * 55)

r_estimates = []
for source, g in growth_estimates.items():
    D1_P0 = D0_P0 * (1 + g)
    r = D1_P0 + g
    r_estimates.append(r)
    print(f"{source:>18s} {g:>8.1%} {D1_P0:>10.2%} {r:>15.2%}")

print(f"\nRange of cost of equity estimates: {min(r_estimates):.1%} to {max(r_estimates):.1%}")
print(f"Average: {np.mean(r_estimates):.1%}")
print("\nNote: Duke Power is a regulated utility — low risk, stable dividends.")
print("Its cost of equity is correspondingly modest (~10%) compared to")
print("high-growth tech firms (which might have r = 15-20%).")

# ============================================================
# Return Decomposition: Dividend Yield vs. Capital Gains
# ============================================================
# Show how return composition changes with growth rate

g_range = np.linspace(0, 0.09, 200)
r_fixed = 0.12
D0 = 1.0

div_yields = np.array([D0 * (1 + g) / gordon_growth_price(D0, r_fixed, g) if r_fixed > g else np.nan for g in g_range])
cap_gains = g_range

fig, ax = plt.subplots(figsize=(11, 6))

ax.fill_between(g_range * 100, 0, div_yields * 100, alpha=0.4, color='#3498db', label='Dividend Yield')
ax.fill_between(g_range * 100, div_yields * 100, (div_yields + cap_gains) * 100,
                alpha=0.4, color='#e67e22', label='Capital Gains (g)')

ax.axhline(y=r_fixed * 100, color='#2c3e50', linewidth=2, linestyle='--', label=f'Total Return = r = {r_fixed:.0%}')

ax.set_xlabel('Growth Rate g (%)', fontsize=13)
ax.set_ylabel('Return Component (%)', fontsize=13)
ax.set_title('Return Decomposition: r = Dividend Yield + Capital Gains',
             fontsize=14, fontweight='bold')
ax.legend(fontsize=11)
ax.set_ylim(0, 14)
ax.set_xlim(0, 9)

# Annotations
ax.annotate('Low-growth (utilities):\nhigh dividend yield', xy=(1.5, 10), fontsize=10,
            color='#3498db', fontweight='bold',
            bbox=dict(boxstyle='round', facecolor='#eaf2f8', alpha=0.9))
ax.annotate('High-growth (tech):\nlow dividend yield,\nhigh capital gains', xy=(6.5, 5), fontsize=10,
            color='#e67e22', fontweight='bold',
            bbox=dict(boxstyle='round', facecolor='#fdf6e3', alpha=0.9))

plt.tight_layout()
plt.show()

print("The total return is always r (by the Gordon model: r = D₁/P₀ + g).")
print("As growth increases, the mix shifts from dividend yield to capital gains.")
print("This is why growth stocks have low dividend yields — it's not a flaw!")

# ============================================================
# Two-Stage Growth DDM
# ============================================================

def two_stage_ddm(D0, r, g_high, T_high, g_low):
    """Two-stage DDM: high growth for T years, then stable growth forever."""
    # Phase 1: high-growth dividends
    phase1_pv = 0
    DT = D0
    for t in range(1, T_high + 1):
        DT = D0 * (1 + g_high)**t
        phase1_pv += DT / (1 + r)**t
    
    # Phase 2: terminal value using Gordon model
    D_terminal = DT * (1 + g_low)
    P_T = D_terminal / (r - g_low) if r > g_low else float('inf')
    phase2_pv = P_T / (1 + r)**T_high
    
    return phase1_pv, phase2_pv, phase1_pv + phase2_pv


# Lo's example
D0 = 1.00
r = 0.20
g_H = 0.06
T = 7
g_L = 0.00

pv1, pv2, P0 = two_stage_ddm(D0, r, g_H, T, g_L)

print("=" * 60)
print("TWO-STAGE DDM — Lo's Example")
print("=" * 60)
print(f"D₀ = ${D0:.2f}, r = {r:.0%}, g_H = {g_H:.0%} (years 1–{T}), g_L = {g_L:.0%} (year {T+1}+)")
print()

# Show dividend timeline
print("Dividend timeline:")
for t in range(1, T + 2):
    Dt = D0 * (1 + g_H)**min(t, T)
    if t > T:
        Dt = D0 * (1 + g_H)**T * (1 + g_L)
    phase = "← high growth" if t <= T else "← stable"
    print(f"  D_{t} = ${Dt:.4f}  {phase}")
    if t == T + 1:
        break

# Terminal price
D_T = D0 * (1 + g_H)**T
P_T = D_T * (1 + g_L) / (r - g_L)

print(f"\nTerminal value at T={T}: P_{T} = D_{T+1}/(r-g_L) = ${D_T*(1+g_L):.4f}/{r-g_L:.2f} = ${P_T:.2f}")
print(f"\nPV of high-growth dividends (years 1–{T}): ${pv1:.4f}")
print(f"PV of terminal value:                      ${pv2:.4f}")
print(f"Total stock price P₀:                      ${P0:.4f}")
print(f"\nTerminal value share: {pv2/(pv1+pv2)*100:.1f}% of total value")

# Compare with pure Gordon models
P0_allhigh = gordon_growth_price(D0, r, g_H)
P0_no_growth = D0 * (1 + 0) / (r - 0)  # g=0 forever

print(f"\nComparison:")
print(f"  If g = {g_H:.0%} forever:   P₀ = ${P0_allhigh:.2f}")
print(f"  Two-stage model:     P₀ = ${P0:.2f}")
print(f"  If g = {g_L:.0%} forever:   P₀ = ${P0_no_growth:.2f}")

# ============================================================
# Visualization: Two-Stage DDM Cash Flows and Valuation
# ============================================================
# ▶ MODIFY AND RE-RUN
D0 = 2.00
r = 0.12
g_high = 0.15     # high-growth phase
T_high = 5        # years of high growth
g_low = 0.04      # stable growth
# ============================================================

pv1, pv2, P0 = two_stage_ddm(D0, r, g_high, T_high, g_low)

# Build dividend series
T_plot = 25
divs = []
pvs = []
for t in range(1, T_plot + 1):
    if t <= T_high:
        Dt = D0 * (1 + g_high)**t
    else:
        Dt = D0 * (1 + g_high)**T_high * (1 + g_low)**(t - T_high)
    divs.append(Dt)
    pvs.append(Dt / (1 + r)**t)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# Left: Dividend timeline
years_plot = range(1, T_plot + 1)
colors_bar = ['#e74c3c' if t <= T_high else '#3498db' for t in years_plot]
ax1.bar(years_plot, divs, color=colors_bar, alpha=0.75, edgecolor='white')
ax1.axvline(x=T_high + 0.5, color='#2c3e50', linewidth=2, linestyle=':', alpha=0.7)
ax1.text(T_high / 2, max(divs) * 0.9, f'High Growth\n(g = {g_high:.0%})',
         ha='center', fontsize=11, color='#e74c3c', fontweight='bold')
ax1.text(T_high + (T_plot - T_high) / 2, max(divs) * 0.9, f'Stable Growth\n(g = {g_low:.0%})',
         ha='center', fontsize=11, color='#3498db', fontweight='bold')
ax1.set_xlabel('Year', fontsize=12)
ax1.set_ylabel('Dividend ($)', fontsize=12)
ax1.set_title('Expected Dividend Stream', fontsize=14, fontweight='bold')

# Right: PV contributions
ax2.bar(years_plot, pvs, color=colors_bar, alpha=0.75, edgecolor='white')
ax2.axvline(x=T_high + 0.5, color='#2c3e50', linewidth=2, linestyle=':', alpha=0.7)

# Show cumulative reaching total
cumulative = np.cumsum(pvs)
ax2_twin = ax2.twinx()
ax2_twin.plot(years_plot, cumulative, color='#2c3e50', linewidth=2.5, marker='o', markersize=3)
ax2_twin.axhline(y=P0, color='#e67e22', linewidth=2, linestyle='--')
ax2_twin.set_ylabel('Cumulative PV ($)', fontsize=12, color='#2c3e50')
ax2_twin.tick_params(axis='y', colors='#2c3e50')

ax2.set_xlabel('Year', fontsize=12)
ax2.set_ylabel('PV of Dividend ($)', fontsize=12)
ax2.set_title(f'PV of Each Dividend (P₀ = ${P0:.2f})', fontsize=14, fontweight='bold')

plt.tight_layout()
plt.show()

print(f"Two-stage DDM: P₀ = ${P0:.2f}")
print(f"  PV of high-growth dividends: ${pv1:.2f} ({pv1/P0*100:.1f}%)")
print(f"  PV of terminal value:        ${pv2:.2f} ({pv2/P0*100:.1f}%)")
print(f"\nThe terminal value typically dominates — most of a stock's value")
print(f"comes from the LONG-RUN growth phase, not the near-term high growth.")

# ============================================================
# Texas Western Example — EPS, Growth, and Valuation
# ============================================================

# Company parameters
EPS1 = 1.00
BV0 = 10.00
g = 0.08
r = 0.10
ROE = EPS1 / BV0  # = 10%
b = BV0 * g / EPS1  # plowback ratio
p = 1 - b             # payout ratio
D1 = EPS1 * p         # dividend per share

print("=" * 60)
print("TEXAS WESTERN — Lo's Example")
print("=" * 60)
print(f"EPS₁ = ${EPS1:.2f},  BV₀ = ${BV0:.2f}")
print(f"ROE = EPS/BV = {ROE:.0%}")
print(f"Growth rate g = {g:.0%}")
print(f"Plowback ratio b = {b:.2f}")
print(f"Payout ratio p = {p:.2f}")
print(f"DPS₁ = EPS × p = ${D1:.2f}")
print(f"Sustainable growth: b × ROE = {b:.2f} × {ROE:.2f} = {b * ROE:.0%} = g ✓")

# CASE 1: Expand at 8% forever
P0_case1 = D1 / (r - g)
print(f"\nCase 1: Expand at {g:.0%} forever")
print(f"  P₀ = D₁/(r-g) = {D1}/{r}-{g} = ${P0_case1:.2f}")
print(f"  P/E = {P0_case1/EPS1:.1f}×")

# CASE 2: Growth slows to 4% after year 5
g_high = 0.08
g_low = 0.04
T_trans = 5

# Build explicit year-by-year forecast
print(f"\nCase 2: {g_high:.0%} for {T_trans} years, then {g_low:.0%}")
print(f"\n{'Year':>6s} {'BV_start':>10s} {'EPS':>8s} {'Inv':>8s} {'DPS':>8s} {'b':>6s}")
print("-" * 55)

# Phase 1
BV_t = BV0
EPS_t = EPS1
pv_divs = 0
for t in range(1, T_trans + 1):
    EPS_t = BV_t * ROE
    inv_t = BV_t * g_high
    D_t = EPS_t - inv_t
    pv_divs += D_t / (1 + r)**t
    b_t = inv_t / EPS_t
    print(f"{t:>6d} {BV_t:>10.4f} {EPS_t:>8.4f} {inv_t:>8.4f} {D_t:>8.4f} {b_t:>6.2f}")
    BV_t = BV_t * (1 + g_high)

# Terminal value: from year T+1, ROE still 10% but growth drops to 4%
# New plowback: b_new = g_low / ROE
b_new = g_low / ROE
p_new = 1 - b_new
EPS_T1 = BV_t * ROE
D_T1 = EPS_T1 * p_new
P_T = D_T1 / (r - g_low)
pv_terminal = P_T / (1 + r)**T_trans

P0_case2 = pv_divs + pv_terminal

print(f"\nAt year {T_trans+1}: EPS = ${EPS_T1:.4f}, new payout = {p_new:.0%}, DPS = ${D_T1:.4f}")
print(f"Terminal price: P_{T_trans} = ${P_T:.2f}")
print(f"\nPV of dividends (years 1-{T_trans}): ${pv_divs:.2f}")
print(f"PV of terminal value:            ${pv_terminal:.2f}")
print(f"Stock price P₀:                  ${P0_case2:.2f}")
print(f"P/E = {P0_case2/EPS1:.1f}×")

print(f"\n→ BOTH cases give P₀ = ${P0_case1:.2f}! Why?")
print(f"  Because ROE = r = {ROE:.0%}. When the firm earns EXACTLY its cost of")
print(f"  capital on new investments, growth adds NO value.")
print(f"  Price = EPS₁/r = ${EPS1/r:.2f} regardless of growth strategy.")

# ============================================================
# Lo's ABC Software Example — PVGO Decomposition
# ============================================================
# From slide 18
EPS1_abc = 8.33
p_abc = 0.6        # payout ratio
ROE_abc = 0.25     # return on equity
r_abc = 0.15       # cost of capital

b_abc = 1 - p_abc  # plowback ratio = 0.4
g_abc = b_abc * ROE_abc  # sustainable growth = 10%
D1_abc = EPS1_abc * p_abc

print("=" * 65)
print("ABC SOFTWARE — Growth Opportunities (Lo's Example)")
print("=" * 65)
print(f"EPS₁ = ${EPS1_abc:.2f},  payout = {p_abc:.0%},  ROE = {ROE_abc:.0%},  r = {r_abc:.0%}")
print(f"Plowback b = {b_abc:.1f},  g = b × ROE = {b_abc:.1f} × {ROE_abc:.0%} = {g_abc:.0%}")
print(f"D₁ = EPS × p = ${D1_abc:.2f}")

# No-growth value
P0_nogrowth = EPS1_abc / r_abc
print(f"\nNo-growth value: EPS₁/r = ${EPS1_abc:.2f}/{r_abc:.2f} = ${P0_nogrowth:.2f}")

# Growth value
P0_growth = D1_abc / (r_abc - g_abc)
print(f"Growth value:    D₁/(r-g) = ${D1_abc:.2f}/({r_abc:.2f}-{g_abc:.2f}) = ${P0_growth:.2f}")

PVGO = P0_growth - P0_nogrowth
print(f"\nPVGO = P₀(growth) - P₀(no-growth) = ${P0_growth:.2f} - ${P0_nogrowth:.2f} = ${PVGO:.2f}")
print(f"PVGO share: {PVGO/P0_growth*100:.1f}% of stock price")

# Verify PVGO via NPV of growth investments
# At t=1: invest b*EPS = 0.4*8.33 = 3.33 at ROE=25% → perpetual CF of 0.25*3.33 = 0.83
inv_1 = b_abc * EPS1_abc
cf_from_inv1 = ROE_abc * inv_1
NPV_1 = -inv_1 + cf_from_inv1 / r_abc  # NPV at t=1

# PVGO = NPV_1 / (r - g) as a growing perpetuity
PVGO_check = NPV_1 / (r_abc - g_abc)

print(f"\nVerification via NPV of growth investments:")
print(f"  Year 1 investment: ${inv_1:.2f} at ROE={ROE_abc:.0%}")
print(f"  Perpetual CF from investment: ${cf_from_inv1:.2f}/year")
print(f"  NPV at t=1: -${inv_1:.2f} + ${cf_from_inv1:.2f}/{r_abc:.2f} = ${NPV_1:.2f}")
print(f"  PVGO = NPV₁/(r-g) = ${NPV_1:.2f}/({r_abc:.2f}-{g_abc:.2f}) = ${PVGO_check:.2f} ✓")

# ============================================================
# PVGO Decomposition — Visualization
# ============================================================

def pvgo_decomposition(EPS1, r, ROE, b):
    """Compute price, no-growth value, and PVGO."""
    g = b * ROE
    D1 = EPS1 * (1 - b)
    if r <= g:
        return float('inf'), EPS1/r, float('inf'), g
    P0 = D1 / (r - g)
    P_nogrowth = EPS1 / r
    PVGO = P0 - P_nogrowth
    return P0, P_nogrowth, PVGO, g

# Compare firms with different ROE
firms = [
    ('Value Destroyer',  8.33, 0.15, 0.10, 0.4),  # ROE < r
    ('Zero-Value Growth', 8.33, 0.15, 0.15, 0.4),  # ROE = r
    ('ABC Software',     8.33, 0.15, 0.25, 0.4),   # ROE > r
    ('Super Growth',     8.33, 0.15, 0.35, 0.4),   # ROE >> r
]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

names = []
no_growth_vals = []
pvgo_vals = []
pe_ratios = []

print("=" * 80)
print("PVGO DECOMPOSITION ACROSS FIRMS (all: EPS₁=$8.33, r=15%, b=0.4)")
print("=" * 80)
print(f"{'Firm':>20s} {'ROE':>6s} {'g':>6s} {'P₀':>8s} {'EPS/r':>8s} {'PVGO':>8s} {'P/E':>6s} {'PVGO%':>7s}")
print("-" * 80)

for name, eps, r, roe, b in firms:
    P0, P_ng, pvgo, g = pvgo_decomposition(eps, r, roe, b)
    pe = P0 / eps
    names.append(name.replace(' ', '\n'))
    no_growth_vals.append(P_ng)
    pvgo_vals.append(pvgo)
    pe_ratios.append(pe)
    print(f"{name:>20s} {roe:>6.0%} {g:>6.0%} ${P0:>7.2f} ${P_ng:>7.2f} ${pvgo:>+7.2f} {pe:>6.1f} {pvgo/P0*100:>6.1f}%")

# Left: Stacked bar chart
x = np.arange(len(names))
ax1.bar(x, no_growth_vals, color='#3498db', alpha=0.8, label='No-growth (EPS₁/r)', edgecolor='white')
for i, pv in enumerate(pvgo_vals):
    color = '#2ecc71' if pv >= 0 else '#e74c3c'
    ax1.bar(x[i], pv, bottom=no_growth_vals[i], color=color, alpha=0.8, edgecolor='white')

# Custom legend
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor='#3498db', alpha=0.8, label='No-growth value (EPS₁/r)'),
                   Patch(facecolor='#2ecc71', alpha=0.8, label='PVGO > 0 (value creating)'),
                   Patch(facecolor='#e74c3c', alpha=0.8, label='PVGO < 0 (value destroying)')]
ax1.legend(handles=legend_elements, fontsize=9, loc='upper left')

ax1.set_xticks(x)
ax1.set_xticklabels(names, fontsize=9)
ax1.set_ylabel('Stock Price ($)', fontsize=12)
ax1.set_title('Stock Price = No-Growth + PVGO', fontsize=14, fontweight='bold')
ax1.axhline(y=no_growth_vals[0], color='gray', linewidth=1, linestyle=':', alpha=0.5)
ax1.yaxis.set_major_formatter(mticker.FuncFormatter(lambda x, _: f'${x:.0f}'))

# Right: P/E ratio
colors_pe = ['#e74c3c', '#e67e22', '#2ecc71', '#27ae60']
ax2.bar(x, pe_ratios, color=colors_pe, alpha=0.8, edgecolor='white')
ax2.axhline(y=1/0.15, color='#2c3e50', linewidth=2, linestyle='--', 
            label=f'1/r = {1/0.15:.1f}× (no-growth P/E)')
ax2.set_xticks(x)
ax2.set_xticklabels(names, fontsize=9)
ax2.set_ylabel('P/E Ratio', fontsize=12)
ax2.set_title('P/E Ratio Reflects Growth Opportunities', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)

for i, pe in enumerate(pe_ratios):
    ax2.text(i, pe + 0.3, f'{pe:.1f}×', ha='center', fontsize=10, fontweight='bold')

plt.tight_layout()
plt.show()

# ============================================================
# When Does Growth Create/Destroy Value? — ROE vs. r
# ============================================================
EPS1 = 8.33
r = 0.15
b = 0.4

ROE_range = np.linspace(0.01, 0.40, 200)
prices = []
pvgos = []
pes = []

for roe in ROE_range:
    g = b * roe
    if r > g:
        D1 = EPS1 * (1 - b)
        P0 = D1 / (r - g)
        pvgo = P0 - EPS1/r
        prices.append(P0)
        pvgos.append(pvgo)
        pes.append(P0/EPS1)
    else:
        prices.append(np.nan)
        pvgos.append(np.nan)
        pes.append(np.nan)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# Left: PVGO vs ROE
ax1.plot(ROE_range * 100, pvgos, color='#2c3e50', linewidth=3)
ax1.axhline(y=0, color='#e74c3c', linewidth=2, linestyle='--')
ax1.axvline(x=r * 100, color='#e67e22', linewidth=2, linestyle=':', alpha=0.7)
ax1.fill_between(ROE_range * 100, 0, pvgos, where=np.array(pvgos) > 0,
                alpha=0.15, color='#2ecc71', label='PVGO > 0 (growth creates value)')
ax1.fill_between(ROE_range * 100, 0, pvgos, where=np.array(pvgos) < 0,
                alpha=0.15, color='#e74c3c', label='PVGO < 0 (growth destroys value)')

ax1.annotate(f'ROE = r = {r:.0%}\n(breakeven)', xy=(r*100, 0), xytext=(r*100 + 5, -15),
            fontsize=10, color='#e67e22', fontweight='bold',
            arrowprops=dict(arrowstyle='->', color='#e67e22', lw=1.5))

ax1.set_xlabel('Return on Equity (ROE) %', fontsize=12)
ax1.set_ylabel('PVGO ($)', fontsize=12)
ax1.set_title('PVGO vs. ROE\n(r = 15%, b = 0.4)', fontsize=14, fontweight='bold')
ax1.legend(fontsize=9)

# Right: P/E vs ROE
ax2.plot(ROE_range * 100, pes, color='#2c3e50', linewidth=3)
ax2.axhline(y=1/r, color='#e67e22', linewidth=2, linestyle='--', label=f'1/r = {1/r:.1f}× (no-growth P/E)')
ax2.axvline(x=r * 100, color='#e67e22', linewidth=1, linestyle=':', alpha=0.5)

ax2.set_xlabel('Return on Equity (ROE) %', fontsize=12)
ax2.set_ylabel('P/E Ratio', fontsize=12)
ax2.set_title('P/E Ratio vs. ROE\n(r = 15%, b = 0.4)', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.set_ylim(0, 30)

plt.tight_layout()
plt.show()

print("Critical insight:")
print("• ROE > r: Growth creates value → PVGO > 0 → P/E > 1/r (growth stock)")
print("• ROE = r: Growth is value-neutral → PVGO = 0 → P/E = 1/r")
print("• ROE < r: Growth DESTROYS value → PVGO < 0 → P/E < 1/r")
print("\nA 'growth stock' is defined by ROE > r, NOT by growing EPS!")

# Exercise 1 — Workspace
# D0 = 3.50; g = 0.03; r = 0.09
#
# (a) P0 = D0*(1+g)/(r-g) = 3.50*1.03/(0.09-0.03)
# P0 = ?
#
# (b) Trading at $65: r = D1/P0 + g = 3.50*1.03/65 + 0.03
# r = ?
#
# (c) New g=0.05, new r=0.11
# P0_new = 3.50*(1.05)/(0.11-0.05)
# Compare with original
#
# (d) 55 = 3.50*(1+g)/(0.09-g)  →  solve for g
# from scipy.optimize import brentq
# g_implied = brentq(lambda g: 3.50*(1+g)/(0.09-g) - 55, 0, 0.089)
# print(f"Implied g = {g_implied:.2%}")

# Exercise 2 — Workspace
# D0, r, g_H, T, g_L = 1.20, 0.14, 0.20, 5, 0.04
#
# (a)
# pv1, pv2, P0 = two_stage_ddm(D0, r, g_H, T, g_L)
# print(f"P₀ = ${P0:.2f}")
#
# (b)
# print(f"Terminal value share: {pv2/(pv1+pv2)*100:.1f}%")
#
# (c) Find g_L such that P0 = 35
# g_implied = brentq(lambda g: two_stage_ddm(D0, r, g_H, T, g)[2] - 35, 0, 0.139)
# print(f"Implied g_L = {g_implied:.2%}")

# Exercise 3 — Workspace
# EPS1, p, ROE, r = 5.00, 0.50, 0.18, 0.12
# b = 1 - p
# g = b * ROE
# D1 = EPS1 * p
# P0 = D1 / (r - g)
# P_nogrowth = EPS1 / r
# PVGO = P0 - P_nogrowth
# PE = P0 / EPS1
#
# print(f"g = {g:.2%}, D₁ = ${D1:.2f}")
# print(f"P₀ = ${P0:.2f}, P/E = {PE:.1f}×")
# print(f"No-growth: ${P_nogrowth:.2f}, PVGO: ${PVGO:.2f} ({PVGO/P0*100:.1f}%)")