import numpy as np
import pandas as pd
from statsmodels.api import OLS, add_constant
from scipy import stats
def chow_test_tax_series(
df_annual: pd.DataFrame,
y_col: str,
break_year: int = 1971,
min_obs_each_side: int = 8,
save_csv_path: str = None
):
"""
Chow test for a structural break at a known year on an annual series.
Model: y ~ const + t
where t is a normalized year index for numerical stability.
H0: No structural break (β_pre = β_post)
H1: Structural break exists
"""
print("\n" + "=" * 80)
print("CHOW TEST FOR STRUCTURAL BREAK")
print("=" * 80)
# Clean and keep what we need
df_clean = df_annual[["Year", y_col]].dropna().copy()
df_clean = df_clean.sort_values("Year").reset_index(drop=True)
# Split (match your original: pre < break_year, post >= break_year)
df_pre = df_clean[df_clean["Year"] < break_year].copy()
df_post = df_clean[df_clean["Year"] >= break_year].copy()
n_pre, n_post = len(df_pre), len(df_post)
n_total = n_pre + n_post
print(f"\nBreak year: {break_year}")
print(f"Observations pre-break: {n_pre}")
print(f"Observations post-break: {n_post}")
print(f"Total observations: {n_total}")
if n_pre < min_obs_each_side or n_post < min_obs_each_side:
raise ValueError(
f"Not enough annual observations on each side of {break_year}. "
f"Need at least {min_obs_each_side} each side. Got pre={n_pre}, post={n_post}."
)
# Normalize time index for stability (use the same reference for all)
base_year = df_clean["Year"].min()
df_clean["t"] = df_clean["Year"] - base_year
df_pre["t"] = df_pre["Year"] - base_year
df_post["t"] = df_post["Year"] - base_year
# Pooled regression
X_pooled = add_constant(df_clean["t"].astype(float))
y_pooled = df_clean[y_col].astype(float)
model_pooled = OLS(y_pooled, X_pooled).fit()
RSS_pooled = float(model_pooled.ssr)
# Separate regressions
X_pre = add_constant(df_pre["t"].astype(float))
y_pre = df_pre[y_col].astype(float)
model_pre = OLS(y_pre, X_pre).fit()
RSS_pre = float(model_pre.ssr)
X_post = add_constant(df_post["t"].astype(float))
y_post = df_post[y_col].astype(float)
model_post = OLS(y_post, X_post).fit()
RSS_post = float(model_post.ssr)
RSS_separate = RSS_pre + RSS_post
# Chow statistic
k = 2 # intercept + slope
df_num = k
df_den = n_total - 2 * k
F_stat = ((RSS_pooled - RSS_separate) / k) / (RSS_separate / df_den)
p_value = 1 - stats.f.cdf(F_stat, df_num, df_den)
critical_05 = stats.f.ppf(0.95, df_num, df_den)
critical_01 = stats.f.ppf(0.99, df_num, df_den)
print("\n" + "-" * 80)
print("RESULTS:")
print("-" * 80)
print(f"\nRSS Pooled (restricted): {RSS_pooled:.4f}")
print(f"RSS Separate (unrestricted): {RSS_separate:.4f}")
print(f" RSS Pre-break: {RSS_pre:.4f}")
print(f" RSS Post-break: {RSS_post:.4f}")
print(f"\nChow F-statistic: {F_stat:.4f}")
print(f"Degrees of freedom: ({df_num}, {df_den})")
print(f"P-value: {p_value:.6f}")
print(f"\nCritical values:")
print(f" 5% level: {critical_05:.4f}")
print(f" 1% level: {critical_01:.4f}")
print("\nConclusion:")
if p_value < 0.01:
result = "Reject H0 (1%)"
print(f" *** REJECT NULL HYPOTHESIS at 1% level ***")
print(f" Strong evidence of structural break at {break_year}")
elif p_value < 0.05:
result = "Reject H0 (5%)"
print(f" ** REJECT NULL HYPOTHESIS at 5% level **")
print(f" Significant evidence of structural break at {break_year}")
elif p_value < 0.10:
result = "Reject H0 (10%)"
print(f" * REJECT NULL HYPOTHESIS at 10% level *")
print(f" Moderate evidence of structural break at {break_year}")
else:
result = "Fail to Reject H0"
print(f" FAIL TO REJECT NULL HYPOTHESIS")
print(f" Insufficient evidence of structural break at {break_year}")
# Coefficients
print("\n" + "-" * 80)
print("REGRESSION COEFFICIENTS:")
print("-" * 80)
print(f"\nPre-{break_year} period:")
print(f" Intercept: {model_pre.params.iloc[0]:.4f} (se: {model_pre.bse.iloc[0]:.4f})")
print(f" Slope: {model_pre.params.iloc[1]:.4f} (se: {model_pre.bse.iloc[1]:.4f})")
print(f" R-squared: {model_pre.rsquared:.4f}")
print(f"\nPost-{break_year} period:")
print(f" Intercept: {model_post.params.iloc[0]:.4f} (se: {model_post.bse.iloc[0]:.4f})")
print(f" Slope: {model_post.params.iloc[1]:.4f} (se: {model_post.bse.iloc[1]:.4f})")
print(f" R-squared: {model_post.rsquared:.4f}")
print(f"\nPooled (no break):")
print(f" Intercept: {model_pooled.params.iloc[0]:.4f} (se: {model_pooled.bse.iloc[0]:.4f})")
print(f" Slope: {model_pooled.params.iloc[1]:.4f} (se: {model_pooled.bse.iloc[1]:.4f})")
print(f" R-squared: {model_pooled.rsquared:.4f}")
results_dict = {
"Test": "Chow Test",
"Series": y_col,
"Break Year": break_year,
"F-statistic": float(F_stat),
"P-value": float(p_value),
"Critical Value (5%)": float(critical_05),
"Critical Value (1%)": float(critical_01),
"Result": result,
"Slope Pre": float(model_pre.params.iloc[1]),
"Slope Post": float(model_post.params.iloc[1]),
"Slope Change": float(model_post.params.iloc[1] - model_pre.params.iloc[1]),
"n_pre": int(n_pre),
"n_post": int(n_post),
}
if save_csv_path is not None:
pd.DataFrame([results_dict]).to_csv(save_csv_path, index=False)
return results_dict
