Какими методами можно улучшить предсказания временного ряда с множеством аномальных значений

Есть датасет с некоторым количеством ежедневных событий. Ссылка https://disk.yandex.ru/d/BDsxHn15sU7hMA

Нужно спрогнозировать события на следующие n дней. Данные имеют следующий вид: )

Декомпозиция данных на тренд, сезонность и остатки:

В данных есть недельная сезонность.

Делал прогноз при помощи МL моделей (Линейные с регуляризацией, LightGBM, Catboost), но результат не очень хороший, вероятно из-за множества аномалий. В качестве метрик использовал RMSE, MAE. Для моделей использовал генерацию новых признаков.

def make_features(data, max_lag, rolling_mean_size):
    data['year'] = data.index.month
    data['month'] = data.index.month
    data['day'] = data.index.day
    data['dayofweek'] = data.index.dayofweek
    for lag in range(1, max_lag+1):
        data['lag_{}'.format(lag)]= data['events'].shift(lag)
    data['rolling_mean'] = data['events'].shift().rolling(rolling_mean_size).mean()
    return data
 

def model_linear(model):
    
    best_lag = 0
    best_roll_size = 0
    best_alpha = 0
    best_rmse = 10
    
    for lag in range(1, 100, 5):
        for roll_size in range(1, 100, 5):
            for alpha in np.logspace(-3, 1, 5):

                df = make_features(data_resample, lag, roll_size)
                train, test = train_test_split(df, shuffle=False, test_size=0.1, random_state=2021)
                train = train.dropna()
                
                X_train = train.drop('events', axis = 1)
                y_train = train['events']
                X_test = test.drop('events', axis = 1)
                y_test = test['events']
                
                if model == 'Lasso' :
                    lm = Lasso(alpha=alpha)
                elif model == 'Ridge' :
                    lm = Ridge(alpha=alpha)
                lm.fit(X_train, y_train)

                preds_test = lm.predict(X_test)
                rmse = mean_squared_error(y_test, preds_test) ** 0.5
                if rmse < best_rmse :
                    best_rmse = rmse
                    best_lag = lag
                    best_roll_size = roll_size
                    best_alpha = alpha
    return best_rmse, best_lag, best_roll_size, best_alpha, preds_test

model_lgb = lgb.LGBMRegressor(random_state=12345)
parameters_lgb = {'max_depth': range(5, 11),
                 'learning_rate': [0.001, 0.01, 0.05, 0.1],
                 'n_estimators': range(100, 300, 50)}
search_lgb = RandomizedSearchCV(model_lgb, parameters_lgb, cv = ts_cv, n_jobs = -1, random_state = 12345)
search_lgb.fit(features_train, target_train)
best_lgb = search_lgb.best_estimator_
predict_lgb_valid = best_lgb.predict(features_valid)
rmse_lgb_valid = mean_squared_error(np.array(target_valid), predict_lgb_valid)**0.5
print('Качество модели LighGBM на валидационной выборке:', rmse_lgb_valid)

Метрика RMSE на тестовом датасете получалась где-то 4-5, при том что в 50% значений имеют среднее значение 5.

На картинке видно, что модель вообще не учла разброс. Пробовал также модель SARIMA, но результат примерно такой же.

Какие есть еще способы улучшить модель?