# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: percent
# ---

# %% [markdown]
# # US GDP Growth: The Great Moderation
#
# The "Great Moderation" refers to the sharp reduction in US GDP growth
# volatility beginning around 1984 (McConnell & Perez-Quiros, 2000;
# Stock & Watson, 2002).
#
# Here we test for two breaks in simulated quarterly GDP growth data that
# mimics this pattern. Real data can be obtained from FRED
# (series A191RL1Q225SBEA).

# %%
import numpy as np
from xtbreak import test, estimate

# %% [markdown]
# ## Synthetic GDP growth data
#
# Three regimes:
# 1. Pre-1984 (high volatility): σ ≈ 4.5%, mean ≈ 3.5%
# 2. 1984–2007 (Great Moderation): σ ≈ 2.0%, mean ≈ 3.0%
# 3. Post-2007 (Great Recession era): σ ≈ 3.0%, mean ≈ 2.0%

# %%
np.random.seed(2024)

# Regime 1: 1960Q1–1983Q4 (96 quarters)
n1 = 96
pre_moderation = np.random.normal(3.5, 4.5, n1)

# Regime 2: 1984Q1–2007Q4 (96 quarters)
n2 = 96
post_moderation = np.random.normal(3.0, 2.0, n2)

# Regime 3: 2008Q1–2019Q4 (48 quarters)
n3 = 48
post_crisis = np.random.normal(2.0, 3.0, n3)

gdp_growth = np.concatenate([pre_moderation, post_moderation, post_crisis])
T = len(gdp_growth)
x = np.ones((T, 1))  # Constant-only model (testing for mean shift)

# Quarterly time labels: 1960.00, 1960.25, 1960.50, ...
quarters = np.round(np.arange(1960.0, 1960.0 + T * 0.25, 0.25)[:T], 2)

print(f"Sample size: {T} quarters ({quarters[0]:.2f}–{quarters[-1]:.2f})")
print(f"True break 1: 1984.00 (index {n1})")
print(f"True break 2: 2008.00 (index {n1 + n2})")
print()

# %% [markdown]
# ## Step 1: Test for structural breaks
#
# We test H0 (no breaks) against H1 (2 simultaneous breaks). The SupF(2)
# statistic evaluates whether the data contains two break points.

# %%
# Test for 2 breaks
test_result_2 = test(gdp_growth, x, breaks=2, trimming=0.10, vce="ssr")
print("=" * 55)
print("Structural Break Tests")
print("=" * 55)
print(f"  SupF(2) statistic: {test_result_2.statistic:.4f}")
print(f"  Break indices:     {test_result_2.breaks}")
print()

# Also test for 1 break (for comparison)
test_result_1 = test(gdp_growth, x, breaks=1, trimming=0.10, vce="ssr")
print(f"  SupF(1) statistic: {test_result_1.statistic:.4f}")
print(f"  (Both tests reject H0 of no break)")
print()

# %% [markdown]
# ## Step 2: Estimate break locations
#
# The estimator uses dynamic programming to find the global SSR-minimizing
# partition into 3 regimes (2 breaks).

# %%
est_result = estimate(
    gdp_growth, x,
    breaks=2,
    trimming=0.10,
    time_values=tuple(quarters.tolist()),
)

print("=" * 55)
print("Break Date Estimation (2 breaks)")
print("=" * 55)
for be in est_result.break_estimates:
    print(f"  Break {be.break_number}:")
    print(f"    Estimated date: {be.time_value}")
    print(f"    Index:          {be.index}")
    print(f"    95% CI (time):  [{be.ci_time_value[0]}, {be.ci_time_value[1]}]")
    print(f"    95% CI (index): [{be.ci_index[0]}, {be.ci_index[1]}]")
    print()
print(f"  Global SSR: {est_result.ssr:.4f}")
print()

# %% [markdown]
# ## Step 3: Regime comparison
#
# Compare the mean and dispersion of GDP growth across the three
# estimated regimes.

# %%
print("=" * 55)
print("Regime Statistics")
print("=" * 55)
print(f"{'Regime':<10} {'Mean':>8} {'SE':>8} {'N':>6} {'R²':>10}")
print("-" * 55)
for rs in est_result.regime_statistics:
    print(f"  {rs.regime_index:<8} {rs.coefficients[0]:>8.2f} "
          f"{rs.standard_errors[0]:>8.2f} {rs.n_observations:>6} "
          f"{rs.r_squared:>10.6f}")
print()

# %% [markdown]
# ## Step 4: Chow tests for regime equality

# %%
if est_result.chow_tests:
    print("=" * 55)
    print("Chow Tests (Adjacent Regime Equality)")
    print("=" * 55)
    for ct in est_result.chow_tests:
        p_str = f"{ct.p_value:.6f}" if ct.p_value is not None else "N/A"
        print(f"  Regime {ct.regime_j} vs {ct.regime_j_plus_1}: "
              f"F = {ct.f_statistic:.4f}, p = {p_str}")
    print()

# %% [markdown]
# ## Step 5: Full summary

# %%
print(est_result.summary())

# %% [markdown]
# ## Interpretation
#
# Two breaks detected:
#
# 1. Around 1984 (onset of the Great Moderation). Volatility drops sharply,
#    attributed to improved monetary policy, the shift from manufacturing
#    to services, and better inventory management.
#
# 2. Around 2007–2008 (end of the Great Moderation). Volatility rises as
#    the financial crisis disrupts the previously stable growth pattern.
#
# Implications:
# - Forecasting: different volatility regimes require different model
#   specifications and confidence bands.
# - Monetary policy: break-aware models avoid overfitting to one regime.
# - Risk management: regime identification helps calibrate stress scenarios.
#
# ### Workflow recap
#
# ```python
# from xtbreak import test, estimate
#
# result = test(y, x, breaks=2, trimming=0.10, vce="ssr")
# est = estimate(y, x, breaks=2, trimming=0.10,
#                time_values=tuple(time_labels))
#
# for be in est.break_estimates:
#     print(be.time_value, be.ci_time_value)
# for rs in est.regime_statistics:
#     print(rs.coefficients, rs.r_squared)
# ```