def analyze_residuals(y_true, y_pred):
    """
    Усовершенствованный анализ остатков для временных рядов:
    - Проверка нормальности (тест Жарка-Бера)
    - Модифицированный тест на гетероскедастичность
    - Визуализация распределения остатков"""
 
    residuals = y_true - y_pred
    flat_residuals = residuals.flatten()
    
    # 1. Тест на нормальность (Жарке-Бера)
    jb_test = jarque_bera(flat_residuals)
    print("\nТест Жарка-Бера на нормальность остатков:")
    print(f"  JB статистика: {jb_test[0]:.4f}")
    print(f"  p-value: {jb_test[1]:.4f}")
    print("  Вывод: Остатки" + (" не" if jb_test[1] < 0.05 else "") + " нормально распределены")
    
    # 2. Модифицированный тест на гетероскедастичность
    try:
        # Для временных рядов используем лагированные остатки в качестве exog
        lagged_residuals = np.column_stack([
            np.ones(len(flat_residuals)),
            np.roll(flat_residuals, 1),  # Лаг 1
            np.roll(flat_residuals, 2)   # Лаг 2
        ])[2:]  # Отбрасываем первые два наблюдения
        
        current_residuals = flat_residuals[2:]
        
        if len(current_residuals) > 10:  # Минимальное количество наблюдений
            bp_test = het_breuschpagan(current_residuals, lagged_residuals)
            print("\nМодифицированный тест на гетероскедастичность:")
            print(f"  LM статистика: {bp_test[0]:.4f}")
            print(f"  p-value: {bp_test[1]:.4f}")
            print("  Вывод: Остатки" + (" имеют" if bp_test[1] < 0.05 else " не имеют") + " гетероскедастичность")
        else:
            print("\nНедостаточно данных для теста на гетероскедастичность")
    except Exception as e:
        print(f"\nОшибка при выполнении теста на гетероскедастичность: {str(e)}")
    
    # 3. Визуализация распределения остатков
    plt.figure(figsize=(12, 5))
    
    # Гистограмма с нормальным распределением
    plt.subplot(1, 2, 1)
    plt.hist(flat_residuals, bins=30, density=True, alpha=0.6, color='g')
    xmin, xmax = plt.xlim()
    x = np.linspace(xmin, xmax, 100)
    p = stats.norm.pdf(x, np.mean(flat_residuals), np.std(flat_residuals))
    plt.plot(x, p, 'k', linewidth=2)
    plt.title('Распределение остатков')
    
    # График остатков во времени
    plt.subplot(1, 2, 2)
    plt.plot(flat_residuals, alpha=0.7)
    plt.axhline(y=0, color='r', linestyle='--')
    plt.xlabel('Временной период')
    plt.ylabel('Остатки')
    plt.title('Динамика остатков')
    
    plt.tight_layout()
    plt.show()
    
    # Дополнительный график: остатки vs предсказания
    plt.figure(figsize=(6, 4))
    plt.scatter(y_pred.flatten(), flat_residuals, alpha=0.5)
    plt.axhline(y=0, color='r', linestyle='--')
    plt.xlabel('Предсказанные значения')
    plt.ylabel('Остатки')
    plt.title('Остатки vs Предсказания')
    plt.tight_layout()
    plt.show()

    


def calculate_pfi_sarimax(model, y_test, X_test, metric_func, n_repeats=10):
    """
    Расчет Permutation Feature Importance (PFI) для SARIMAX модели.
    Возвращает:
    - pfi_abs: абсолютные значения важности (разница метрики)
    - pfi_norm: нормализованные значения (в % от общей важности)
    """
    if X_test is None or len(X_test) == 0:
        return {}, {}

    baseline_predictions = model.predict(n_periods=len(y_test), X=X_test)
    baseline_metric = metric_func(y_test, baseline_predictions)
    
    feature_importance = {i: [] for i in range(X_test.shape[1])}
    
    for _ in range(n_repeats):
        for feature_idx in range(X_test.shape[1]):
            X_test_permuted = X_test.copy()
            np.random.shuffle(X_test_permuted[:, feature_idx])
            permuted_predictions = model.predict(n_periods=len(y_test), X=X_test_permuted)
            permuted_metric = metric_func(y_test, permuted_predictions)
            feature_importance[feature_idx].append(permuted_metric - baseline_metric)
    
    pfi_abs = {feature_idx: np.mean(scores) for feature_idx, scores in feature_importance.items()}
    
    # Нормализация (сумма модулей = 100%)
    total_importance = sum(abs(val) for val in pfi_abs.values())
    pfi_norm = {}
    if total_importance > 0:
        pfi_norm = {feature_idx: (abs(val) / total_importance) * 100 
                   for feature_idx, val in pfi_abs.items()}
    
    return pfi_abs, pfi_norm

def process_targets_sarimax(data, data_for_test, target_cols_list, base_covariates=None, horizon=8):
    all_results = []
    aggregated_feature_importance_abs = {}
    aggregated_feature_importance_norm = {} 
    price_importance_values_abs = []
    price_importance_values_norm = []
    
    if base_covariates is None:
        base_covariates = []
    
    for target in target_cols_list:
        print(f"\n{'='*50}\nОбработка таргета: {target}\n{'='*50}")
        
        try:
            # 1. Формирование ковариат
            current_target_col = target
            product_id = target.split('_')[-1]
            covariates = base_covariates.copy()
            
            possible_specific = [f'Цен_{product_id}', f'Марж руб_{product_id}', f'Марж %_{product_id}']
            covariates.extend([cov for cov in possible_specific if cov in data.columns])
            covariates = list(dict.fromkeys(covariates))
            
            print(f"Используемые ковариаты: {covariates}")
            
            # 2. Подготовка данных
            y_train, X_train, y_train_targets = prepare_sarimax_data(data, target, covariates, horizon)
            y_test, X_test, y_test_targets = prepare_sarimax_data(data_for_test, target, covariates, horizon)
            
            # 3. Масштабирование
            scaler = None
            if covariates and X_train is not None:
                scaler = MinMaxScaler() 
                X_train = scaler.fit_transform(X_train)
                X_test = scaler.transform(X_test) if X_test is not None else None
            
            # 4. Подбор модели
            print("Подбор SARIMAX...")
            model = auto_arima(
                y=y_train,
                X=X_train,
                seasonal=True,
                m=2,
                max_p=3,
                max_q=3,
                max_P=2,
                max_Q=2,
                d=None,
                D=None,
                stepwise=False,
                with_intercept=True,
                information_criterion='aic',
                test='adf',
                seasonal_test='ocsb',
                maxiter=10,
                suppress_warnings=True,
                trace=True
)
            
            print(f"Лучшая модель: {model.order}, {model.seasonal_order}")
            
            # 5. Прогнозирование
            predictions = []
            for i in range(len(y_test_targets)):
                try:
                    exog = None
                    if covariates and X_test is not None:
                        exog = X_test[i:i+horizon] if i + horizon <= len(X_test) else X_test[-horizon:]
        
                    forecast = model.predict(n_periods=horizon, X=exog)
        
                    if np.isscalar(forecast):
                        forecast = np.array([forecast] * horizon)
                    elif len(forecast) < horizon:
                        forecast = np.pad(forecast, (0, horizon - len(forecast)), 'edge')
                    elif len(forecast) > horizon:
                        forecast = forecast[:horizon]
            
                    predictions.append(forecast)
        
                    if i < len(y_test):
                        update_X = X_test[i].reshape(1, -1) if (covariates and X_test is not None) else None
                        model.update([y_test[i]], X=update_X)
    
                except Exception as e:
                    print(f"Ошибка прогнозирования на шаге {i}: {str(e)}")
                    predictions.append(np.full(horizon, y_test[i]))

            predictions = np.array(predictions)
            
            # 6. Метрики и анализ остатков
            current_rmse = rmse(y_test_targets, predictions)
            current_mase = mase(y_test_targets, predictions, y_train)
            
            diff = predictions - y_test_targets
            under_forecast_count = np.sum(diff < 0)
            total_forecasts = len(y_test_targets) * horizon
            under_forecast_ratio = under_forecast_count / total_forecasts * 100

            # Анализ остатков
            analyze_residuals(y_test_targets, predictions)

            # 7. Расчёт PFI
            pfi_abs, pfi_norm = {}, {}
            if covariates and X_test is not None:
                pfi_abs, pfi_norm = calculate_pfi_sarimax(
                    model=model,
                    y_test=y_test,
                    X_test=X_test,
                    metric_func=rmse,
                    n_repeats=20
                )
                
                print("\nВажность ковариат (абсолютная):")
                for idx, importance in sorted(pfi_abs.items(), key=lambda x: x[1], reverse=True):
                    print(f"  {covariates[idx]}: {importance:.4f}")
                
                print("\nВажность ковариат (нормализованная %):")
                for idx, importance in sorted(pfi_norm.items(), key=lambda x: x[1], reverse=True):
                    print(f"  {covariates[idx]}: {importance:.2f}%")
                
                # Агрегация важности
                for idx, importance in pfi_abs.items():
                    feature_name = covariates[idx]
                    if feature_name.startswith('Маржа руб_'):
                        price_importance_values_abs.append(importance)
                    else:
                        if feature_name not in aggregated_feature_importance_abs:
                            aggregated_feature_importance_abs[feature_name] = []
                        aggregated_feature_importance_abs[feature_name].append(importance)
                
                for idx, importance in pfi_norm.items():
                    feature_name = covariates[idx]
                    if feature_name.startswith('Маржа руб_'):
                        price_importance_values_norm.append(importance)
                    else:
                        if feature_name not in aggregated_feature_importance_norm:
                            aggregated_feature_importance_norm[feature_name] = []
                        aggregated_feature_importance_norm[feature_name].append(importance)

            all_results.append({
                'target': target,
                'RMSE': current_rmse,
                'MASE': current_mase,
                'under_forecast_ratio': under_forecast_ratio,
                'model': model,
                'order': model.order,
                'seasonal_order': model.seasonal_order,
                'covariates': covariates,
                'pfi_abs': pfi_abs,
                'pfi_norm': pfi_norm
            })

            print(f"\nRMSE: {current_rmse:.4f}, MASE: {current_mase:.4f}, Under-forecast: {under_forecast_ratio:.2f}%")

        except Exception as e:
            print(f"Ошибка: {str(e)}")
            all_results.append({'target': target, 'error': str(e)})

    # Сводные результаты
    successful_results = [r for r in all_results if 'RMSE' in r and 'MASE' in r]

    if successful_results:
        print("\n" + "="*50)
        print("Сводные результаты по всем моделям:")
        print("="*50)
        
        avg_rmse = np.mean([r['RMSE'] for r in successful_results])
        avg_mase = np.mean([r['MASE'] for r in successful_results])
        avg_under_forecast = np.mean([r['under_forecast_ratio'] for r in successful_results])
        
        print(f"\nКоличество успешных моделей: {len(successful_results)}")
        print(f"\nСредние метрики:")
        print(f"  RMSE: {avg_rmse:.4f}")
        print(f"  MASE: {avg_mase:.4f}")
        print(f"  Under-forecast ratio: {avg_under_forecast:.2f}%")
        
        if aggregated_feature_importance_abs:
            mean_feature_importance_abs = {feature: np.mean(values) for feature, values in aggregated_feature_importance_abs.items()}
            mean_price_importance_abs = np.mean(price_importance_values_abs) if price_importance_values_abs else 0.0
            
            print("\nСредняя важность ковариат (абсолютная):")
            for feature, importance in sorted(mean_feature_importance_abs.items(), key=lambda x: x[1], reverse=True):
                print(f"  {feature}: {importance:.4f}")
            
            print(f"\nСредняя важность цен (все признаки Цена_): {mean_price_importance_abs:.4f}")
            
            # Визуализация
            if mean_feature_importance_abs:
                features, importances = zip(*sorted(mean_feature_importance_abs.items(), key=lambda x: x[1], reverse=True))
                if price_importance_values_abs:
                    features = list(features) + ['Маржа руб']
                    importances = list(importances) + [mean_price_importance_abs]
                
                plt.figure(figsize=(12, 6))
                bars = plt.bar(features, importances, color=['blue' if 'Маржа руб' not in str(f) else 'salmon' for f in features])
                plt.xlabel('Факторы')
                plt.ylabel('Важность (абсолютная, ΔRMSE)')
                plt.title('Средняя важность факторов (абсолютный PFI)')
                plt.xticks(rotation=45, ha='right')
                
                for bar in bars:
                    height = bar.get_height()
                    plt.text(bar.get_x() + bar.get_width()/2., height,
                            f'{height:.3f}',
                            ha='center', va='bottom')
                
                plt.tight_layout()
                plt.show()

        if aggregated_feature_importance_norm:
            mean_feature_importance_norm = {feature: np.mean(values) for feature, values in aggregated_feature_importance_norm.items()}
            mean_price_importance_norm = np.mean(price_importance_values_norm) if price_importance_values_norm else 0.0
            
            print("\nСредняя важность ковариат (нормализованная %):")
            for feature, importance in sorted(mean_feature_importance_norm.items(), key=lambda x: x[1], reverse=True):
                print(f"  {feature}: {importance:.2f}%")
            
            print(f"\nСредняя важность цен (все признаки Маржа руб_): {mean_price_importance_norm:.2f}%")

    return all_results correct mistakes