@@ -457,18 +457,18 @@ def from_model(
457457 # Only constant breaks (mean-shift model)
458458 # exog contains all regressors (non-breaking)
459459 # exog_break defaults to constant
460-
460+
461461 # We need to separate the constant from other regressors
462462 # Check if constant exists in exog_all
463463 # Find ALL constant columns (there may be duplicates)
464464 const_indices : list [int ] = []
465-
465+
466466 for i in range (exog_all .shape [1 ]):
467467 # Check for constant column (all ones)
468468 # We use a loose tolerance because of potential floating point issues
469469 if np .allclose (exog_all [:, i ], 1.0 , atol = 1e-5 ):
470470 const_indices .append (i )
471-
471+
472472 if const_indices :
473473 # Remove all constant columns from exog (non-breaking)
474474 non_breaking_indices = [
@@ -478,10 +478,10 @@ def from_model(
478478 exog_non_break = exog_all [:, non_breaking_indices ]
479479 else :
480480 exog_non_break = None
481-
481+
482482 # exog_break will be a single constant column
483483 exog_break = np .ones ((len (endog ), 1 ))
484-
484+
485485 return cls (endog , exog = exog_non_break , exog_break = exog_break )
486486 else :
487487 # No constant found, so we just use exog_all as non-breaking
@@ -542,7 +542,7 @@ def _compute_ssr_segment(
542542 X_list .append (self .exog_break [start :end ])
543543 if self .exog is not None :
544544 X_list .append (self .exog [start :end ])
545-
545+
546546 if not X_list :
547547 X_seg = None
548548 else :
@@ -855,101 +855,105 @@ def fit(
855855 #
856856 # When p = 0 (pure structural change), no iteration is needed and the
857857 # standard Bai-Perron dynamic programming algorithm is used directly.
858-
858+
859859 if self .p > 0 :
860860 # Partial structural change model
861861 # Iterative procedure per Bai & Perron (2003), Section 3.
862-
862+
863863 # Initial estimate of fixed coefficients (assuming no breaks)
864864 X_full = []
865865 if self .exog_break is not None :
866866 X_full .append (self .exog_break )
867867 if self .exog is not None :
868868 X_full .append (self .exog )
869-
869+
870870 X_mat = np .column_stack (X_full )
871871 # Use OLS to get initial beta
872872 beta_full = np .linalg .lstsq (X_mat , self .endog , rcond = None )[0 ]
873-
873+
874874 # Extract fixed coefficients (last p)
875- beta_fixed = beta_full [- self .p :]
876-
875+ beta_fixed = beta_full [- self .p :]
876+
877877 ssr_vals = {}
878878 breaks_by_m = {}
879-
879+
880880 # m=0 case: No breaks
881881 # SSR is just the SSR from the full model with no breaks
882882 resid_0 = self .endog - X_mat @ beta_full
883883 ssr_vals [0 ] = float (np .sum (resid_0 ** 2 ))
884884 breaks_by_m [0 ] = []
885-
885+
886886 for m in range (1 , max_breaks + 1 ):
887887 # Initialize beta_fixed from m=0 (or could use previous m's result)
888888 # We restart from the no-break estimate to avoid getting stuck in local optima
889889 curr_beta_fixed = beta_fixed .copy ()
890-
890+
891891 # Iterative procedure for specific m
892- # We need to define Full_Design and full_beta outside the loop
892+ # We need to define Full_Design and full_beta outside the loop
893893 # in case the loop doesn't run (though it always runs at least once)
894894 Full_Design = None
895895 full_beta = None
896896 max_iter = 50
897897 converged = False
898-
898+
899899 for _iter in range (max_iter ):
900900 # 1. Partial out fixed regressors
901901 if self .exog is not None :
902902 y_star = self .endog - self .exog @ curr_beta_fixed
903903 else :
904904 y_star = self .endog
905-
905+
906906 # 2. Find optimal breaks for y_star on exog_break
907907 # We pass x=self.exog_break to override default behavior
908- ssr_matrix_m = self ._build_ssr_matrix (h , y = y_star , x = self .exog_break )
908+ ssr_matrix_m = self ._build_ssr_matrix (
909+ h , y = y_star , x = self .exog_break
910+ )
909911 ssr_m , breaks_m = self ._dynamic_programming (m , h , ssr_matrix_m )
910-
912+
911913 # 3. Re-estimate beta_fixed given these breaks
912914 # Construct full design matrix with breaks
913-
915+
914916 # Z_bar columns (breaking regressors)
915917 boundaries = [0 ] + breaks_m + [T ]
916918 Z_cols = []
917919 for i in range (len (boundaries ) - 1 ):
918- start , end = boundaries [i ], boundaries [i + 1 ]
920+ start , end = boundaries [i ], boundaries [i + 1 ]
919921 for col in range (self .q ):
920922 z_col = np .zeros (T )
921923 z_col [start :end ] = self .exog_break [start :end , col ]
922924 Z_cols .append (z_col )
923-
925+
924926 if not Z_cols :
925- # Should not happen if q > 0
926- Z_bar = np .zeros ((T , 0 ))
927+ # Should not happen if q > 0
928+ Z_bar = np .zeros ((T , 0 ))
927929 else :
928930 Z_bar = np .column_stack (Z_cols )
929-
931+
930932 Full_Design = np .column_stack ([Z_bar , self .exog ])
931-
933+
932934 # Estimate
933935 full_beta = np .linalg .lstsq (Full_Design , self .endog , rcond = None )[0 ]
934-
936+
935937 # Update beta_fixed (last p coeffs)
936- new_beta_fixed = full_beta [- self .p :]
937-
938+ new_beta_fixed = full_beta [- self .p :]
939+
938940 # Check convergence
939- if np .allclose (curr_beta_fixed , new_beta_fixed , rtol = 1e-4 , atol = 1e-6 ):
941+ if np .allclose (
942+ curr_beta_fixed , new_beta_fixed , rtol = 1e-4 , atol = 1e-6
943+ ):
940944 curr_beta_fixed = new_beta_fixed
941945 converged = True
942946 break
943947 curr_beta_fixed = new_beta_fixed
944-
948+
945949 if not converged :
946950 warnings .warn (
947951 f"Partial structural change iterative procedure did not "
948952 f"converge for m={ m } after { max_iter } iterations. "
949953 f"Results may be unreliable." ,
950954 stacklevel = 2 ,
951955 )
952-
956+
953957 # Store results for this m
954958 # Re-calculate SSR with final beta
955959 if Full_Design is not None and full_beta is not None :
@@ -961,19 +965,18 @@ def fit(
961965 ssr_vals [m ] = np .inf
962966 breaks_by_m [m ] = []
963967
964-
965968 else :
966969 # Pure structural change (original code)
967970 ssr_matrix = self ._build_ssr_matrix (h )
968-
971+
969972 ssr_vals = {}
970973 breaks_by_m = {}
971-
974+
972975 # No breaks case
973976 ssr_0 = float (ssr_matrix [0 , T - 1 ])
974977 ssr_vals [0 ] = ssr_0
975978 breaks_by_m [0 ] = []
976-
979+
977980 for m in range (1 , max_breaks + 1 ):
978981 ssr_m , breaks_m = self ._dynamic_programming (m , h , ssr_matrix )
979982 ssr_vals [m ] = ssr_m
@@ -987,27 +990,29 @@ def fit(
987990 seqf_critical : dict [int , float ] = {}
988991 bic_vals : dict [int , float ] = {}
989992 lwz_vals : dict [int , float ] = {}
990-
993+
991994 # Compute Information Criteria and SupF stats
992995 for m in range (max_breaks + 1 ):
993996 if m == 0 :
994- bic_vals [0 ], lwz_vals [0 ] = self ._compute_information_criteria (ssr_vals [0 ], 0 )
997+ bic_vals [0 ], lwz_vals [0 ] = self ._compute_information_criteria (
998+ ssr_vals [0 ], 0
999+ )
9951000 continue
996-
1001+
9971002 # Sup-F test
9981003 # F = ((SSR_0 - SSR_m) / (m * q)) / (SSR_m / (T - (m + 1) * q - p))
9991004 df1 = m * q
10001005 df2 = T - (m + 1 ) * q - self .p
1001-
1006+
10021007 ssr_m = ssr_vals [m ]
10031008 ssr_0 = ssr_vals [0 ]
1004-
1009+
10051010 if df2 <= 0 or ssr_m <= 0 :
10061011 supf = 0.0
10071012 else :
10081013 supf = ((ssr_0 - ssr_m ) / df1 ) / (ssr_m / df2 )
10091014 supf = max (0.0 , supf )
1010-
1015+
10111016 supf_stats [m ] = supf
10121017
10131018 # Critical value (use approximation)
0 commit comments