################################################################################
#                        skforecast.model_selection                            #
#                                                                              #
# This work by Joaquín Amat Rodrigo is licensed under a Creative Commons       #
# Attribution 4.0 International License.                                       #
################################################################################
# coding=utf-8


import typing
from typing import Union, Dict, List, Tuple
import numpy as np
import pandas as pd
import logging
import tqdm
from sklearn.metrics import mean_squared_error 
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_absolute_percentage_error
from sklearn.model_selection import ParameterGrid

from .ForecasterAutoreg import ForecasterAutoreg
from .ForecasterCustom import ForecasterCustom
from .ForecasterAutoregMultiOutput import ForecasterAutoregMultiOutput

logging.basicConfig(
    format = '%(asctime)-5s %(name)-10s %(levelname)-5s %(message)s', 
    level  = logging.INFO,
)


def time_series_spliter(y: Union[np.ndarray, pd.Series],
                        initial_train_size: int, steps: int,
                        allow_incomplete_fold: bool=True,
                        verbose: bool=True):
    '''
    
    Split indices of a time series into multiple train-test pairs. The order of
    is maintained and the training set increases in each iteration.
    
    Parameters
    ----------        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    verbose : bool, default `True`
        Print number of splits created.

    Yields
    ------
    train : 1D np.ndarray
        Training indices.
        
    test : 1D np.ndarray
        Test indices.
        
    '''
    
    if not isinstance(y, (np.ndarray, pd.Series)):

        raise Exception('`y` must be `1D np.ndarray` o `pd.Series`.')

    elif isinstance(y, np.ndarray) and y.ndim != 1:

        raise Exception(
            f"`y` must be `1D np.ndarray` o `pd.Series`, "
            f"got `np.ndarray` with {y.ndim} dimensions."
        )
        
    if isinstance(y, pd.Series):
        y = y.to_numpy().copy()
    
  
    folds = (len(y) - initial_train_size) // steps  + 1
    # +1 fold is needed to allow including the remainder in the last iteration.
    remainder = (len(y) - initial_train_size) % steps   
    
    if verbose:
        if folds == 1:
            print(f"Number of folds: {folds - 1}")
            print("Not enought observations in `y` to create to create even a complete fold.")

        elif remainder == 0:
            print(f"Number of folds: {folds - 1}")

        elif remainder != 0 and allow_incomplete_fold:
            print(f"Number of folds: {folds}")
            print(
                f"Since `allow_incomplete_fold=True`, "
                f"last fold only includes {remainder} observations instead of {steps}."
            )
            print(
                'Incomplete folds with few observations could overestimate or ',
                'underestimate validation metrics.'
            )
        elif remainder != 0 and not allow_incomplete_fold:
            print(f"Number of folds: {folds - 1}")
            print(
                f"Since `allow_incomplete_fold=False`, "
                f"last {remainder} observations are descarted."
            )

    if folds == 1:
        # There are no observations to create even a complete fold
        return []
    
    for i in range(folds):
          
        if i < folds - 1:
            train_end     = initial_train_size + i * steps    
            train_indices = range(train_end)
            test_indices  = range(train_end, train_end + steps)
            
        else:
            if remainder != 0 and allow_incomplete_fold:
                train_end     = initial_train_size + i * steps  
                train_indices = range(train_end)
                test_indices  = range(train_end, len(y))
            else:
                break
        
        yield train_indices, test_indices
        
        

def cv_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
                  initial_train_size: int, steps: int,
                  metric: str, exog: Union[np.ndarray, pd.Series]=None,
                  allow_incomplete_fold: bool=True, verbose: bool=True) -> np.array:
    '''
    Cross-validation of `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput`
    object. The order of is maintained and the training set increases in each iteration.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
            
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    verbose : bool, default `True`
        Print number of folds used for cross validation.

    Returns 
    -------
    cv_results: 1D np.ndarray
        Value of the metric for each partition.

    '''
    
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)

    if metric not in ['mean_squared_error', 'mean_absolute_error',
                      'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and "
            f"'mean_absolute_percentage_error'. Got {metric}."
        )
    
    cv_results = []
    
    metrics = {
        'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error
    }
    
    metric = metrics[metric]
    
    splits = time_series_spliter(
                y                     = y,
                initial_train_size    = initial_train_size,
                steps                 = steps,
                allow_incomplete_fold = allow_incomplete_fold,
                verbose               = verbose
             )
    
    for train_index, test_index in splits:
        
        if exog is None:
            forecaster.fit(y=y[train_index])      
            pred = forecaster.predict(steps=len(test_index))
            
        else:
            
            forecaster.fit(y=y[train_index], exog=exog[train_index])      
            pred = forecaster.predict(steps=len(test_index), exog=exog[test_index])
               
        metric_value = metric(
                            y_true = y[test_index],
                            y_pred = pred
                       )
        
        cv_results.append(metric_value)
                          
    return np.array(cv_results)


def backtesting_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
                           initial_train_size: int, steps: int,
                           metric: str, exog: Union[np.ndarray, pd.Series]=None,
                           verbose: bool=False) -> Tuple[np.array, np.array]:
    '''
    Backtesting (validation) of `ForecasterAutoreg`, `ForecasterCustom` or
    `ForecasterAutoregMultiOutput` object.
    The model is trained only once using the `initial_train_size` first observations.
    In each iteration, a number of `steps` predictions are evaluated.
    
    This evaluation is much faster than `cv_forecaster()` since the model is
    trained only once.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
            
    verbose : bool, default `False`
        Print number of folds used for backtesting.

    Returns 
    -------
    backtest_predictions: 1D np.ndarray
        Value of predictions.
        
    metric_value: np.ndarray shape (1,)
        Value of the metric.

    '''
    
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)

    if metric not in ['mean_squared_error', 'mean_absolute_error',
                      'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and "
            f"'mean_absolute_percentage_error'. Got {metric}."
        )
    
    backtest_predictions = []
    
    metrics = {
        'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error
    }
    
    metric = metrics[metric]
        
    if exog is None:
        forecaster.fit(y=y[:initial_train_size])      
    else:
        forecaster.fit(y=y[:initial_train_size], exog=exog[:initial_train_size])     
    
    
    folds     = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    window_size = len(forecaster.last_window)

    if isinstance(forecaster, ForecasterAutoregMultiOutput) and remainder != 0:
        # In ForecasterAutoregMultiOutput predictions are not iterative,
        # therefore no remainder is allowed.
        logging.warning(
                f"Backtesting `ForecasterAutoregMultiOutput` only allow completed "
                f"folds. Last {remainder} observations are excluded."
        )
        remainder=0
    
    if verbose:
        print(f"Number of observations used for training: {initial_train_size}")
        print(f"Number of observations used for testing: {len(y) - initial_train_size}")
        print(f"    Number of folds: {folds - 1 * (remainder == 0)}")
        print(f"    Number of steps per fold: {steps}")
        if remainder != 0:
            print(f"    Last fold only includes {remainder} observations")
      
    for i in range(folds):
        last_window_end   = initial_train_size + i * steps
        last_window_start = (initial_train_size + i * steps) - window_size 
        last_window       = y[last_window_start:last_window_end]
                
        if i < folds - 1:
            if exog is None:
                pred = forecaster.predict(
                            steps       = steps,
                            last_window = last_window
                        )
            else:
                pred = forecaster.predict(
                            steps       = steps,
                            last_window = last_window,
                            exog        = exog[last_window_end:last_window_end + steps]
                        )
        elif remainder != 0:
            steps = remainder 
            if exog is None:
                pred = forecaster.predict(
                            steps       = steps,
                            last_window = last_window
                        )
            else:
                pred = forecaster.predict(
                            steps       = steps,
                            last_window = last_window,
                            exog        = exog[last_window_end:last_window_end + steps]
                        )
        else:
            continue
        
        backtest_predictions.append(pred)
    
    backtest_predictions = np.concatenate(backtest_predictions)
    metric_value = metric(
                        y_true = y[initial_train_size: initial_train_size + len(backtest_predictions)],
                        y_pred = backtest_predictions
                   )

    return np.array([metric_value]), backtest_predictions



def grid_search_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
                           param_grid: dict, initial_train_size: int, steps: int,
                           metric: str, exog: Union[np.ndarray, pd.Series]=None,
                           lags_grid: list=None, method: str='cv',
                           allow_incomplete_fold: bool=True, return_best: bool=True,
                           verbose: bool=True) -> pd.DataFrame:
    '''
    Exhaustive search over specified parameter values for a Forecaster object.
    Validation is done using time series cross-validation or backtesting.
    
    Parameters
    ----------
    
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        ForecasterAutoreg, ForecasterCustom or ForecasterAutoregMultiOutput object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
        
    param_grid : dict
        Dictionary with parameters names (`str`) as keys and lists of parameter
        settings to try as values.
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
           
    lags_grid : list of int, lists, np.narray or range. 
        Lists of `lags` to try. Only used if forecaster is an instance of 
        `ForecasterCustom`.
        
    method : {'cv', 'backtesting'}
        Method used to estimate the metric for each parameter combination.
        'cv' for time series crosvalidation and 'backtesting' for simple
        backtesting. 'backtesting' is much faster since the model is fitted only
        once.
        
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    return_best : bool
        Refit the `forecaster` using the best found parameters on the whole data.
        
    verbose : bool, default `True`
        Print number of folds used for cv or backtesting.

    Returns 
    -------
    results: pandas.DataFrame
        Metric value estimated for each combination of parameters.

    '''
    
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    
    if isinstance(forecaster, ForecasterCustom):
        if lags_grid is not None:
            logging.warning(
                '`lags_grid` ignored if forecaster is an instance of `ForecasterCustom`.'
            )
        lags_grid = ['custom predictors']
        
    elif lags_grid is None:
        lags_grid = [forecaster.lags]
        
    if isinstance(forecaster, ForecasterAutoregMultiOutput):
        param_grid_renamed = {'estimator__'+key: value for key, value in param_grid.items()}
        param_grid = param_grid_renamed
        
    lags_list = []
    params_list = []
    metric_list = []
    
    param_grid =  list(ParameterGrid(param_grid))

    logging.info(
        f"Number of models to fit: {len(param_grid)*len(lags_grid)}"
    )
    
    for lags in tqdm.tqdm(lags_grid, desc='loop lags_grid', position=0):
        
        if isinstance(forecaster, (ForecasterAutoreg, ForecasterAutoregMultiOutput)):
            forecaster.set_lags(lags)
            lags = forecaster.lags.copy()
        
        for params in tqdm.tqdm(param_grid, desc='loop param_grid', position=1, leave=False):

            forecaster.set_params(**params)
            
            if method == 'cv':
                metrics = cv_forecaster(
                                forecaster = forecaster,
                                y          = y,
                                exog       = exog,
                                initial_train_size = initial_train_size,
                                steps  = steps,
                                metric = metric,
                                allow_incomplete_fold = allow_incomplete_fold,
                                verbose = verbose
                             )
            else:
                metrics = backtesting_forecaster(
                                forecaster = forecaster,
                                y          = y,
                                exog       = exog,
                                initial_train_size = initial_train_size,
                                steps  = steps,
                                metric = metric,
                                verbose = verbose
                             )[0]

            lags_list.append(lags)
            params_list.append(params)
            metric_list.append(metrics.mean())
            
    results = pd.DataFrame({
                'lags'  : lags_list,
                'params': params_list,
                'metric': metric_list})
    
    results = results.sort_values(by='metric', ascending=True)
    
    if return_best:
        
        best_lags = results['lags'].iloc[0]
        best_params = results['params'].iloc[0]
        logging.info(
            f"Refitting `forecaster` using the best found parameters and the whole data set: \n"
            f"lags: {best_lags} \n"
            f"params: {best_params}\n"
        )
        
        if isinstance(forecaster, (ForecasterAutoreg, ForecasterAutoregMultiOutput)):
            forecaster.set_lags(best_lags)
            
        forecaster.set_params(**best_params)
        forecaster.fit(y=y, exog=exog)
            
    return results


def backtesting_forecaster_intervals(forecaster, y: Union[np.ndarray, pd.Series],
                           initial_train_size: int, steps: int,
                           metric: str, exog: Union[np.ndarray, pd.Series]=None,
                           interval: list=[5, 95], n_boot: int=500,
                           in_sample_residuals: bool=True,
                           verbose: bool=False) -> Tuple[np.array, np.array]:
    '''
    Backtesting (validation) of `ForecasterAutoreg` or `ForecasterCustom` object.
    The model is trained only once using the `initial_train_size` first observations.
    In each iteration, a number of `steps` predictions are evaluated. Both, 
    predictions and intervals, are calculated.
    
    This evaluation is much faster than `cv_forecaster()` since the model is
    trained only once.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom
        `ForecasterAutoreg` or `ForecasterCustom` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
        
    interval: list, default `[5, 100]`
            Confidence of the prediction interval estimated. Sequence of percentiles
            to compute, which must be between 0 and 100 inclusive.
            
    n_boot: int, default `500`
        Number of bootstrapping iterations used to estimate prediction
        intervals.

    in_sample_residuals: bool, default `True`
        If `True`, residuals from the training data are used as proxy of
        prediction error to create prediction intervals.
        
    verbose : bool, default `True`
        Print number of folds used for backtesting.

    Returns 
    -------
    backtest_predictions: np.ndarray
        If `include_intervals=True`, 2D np.ndarray shape(steps, 3) with predicted
        value and their estimated interval.
            Column 0 = predictions
            Column 1 = lower bound interval
            Column 2 = upper bound interval
        
    metric_value: np.ndarray shape (1,)
        Value of the metric.

    Notes
    -----
    More information about prediction intervals in forecasting:
    https://otexts.com/fpp2/prediction-intervals.html
    Forecasting: Principles and Practice (2nd ed) Rob J Hyndman and
    George Athanasopoulos.

    '''
    
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)

    if metric not in ['mean_squared_error', 'mean_absolute_error',
                      'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and "
            f"'mean_absolute_percentage_error'. Got {metric}."
        )
    
    backtest_predictions = []
    
    metrics = {
        'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error
    }
    
    metric = metrics[metric]
        
    if exog is None:
        forecaster.fit(y=y[:initial_train_size])   
    else:
        forecaster.fit(y=y[:initial_train_size], exog=exog[:initial_train_size])     
    
    
    folds     = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    window_size = len(forecaster.last_window)
    
    if verbose:
        print(f"Number of observations used for training: {initial_train_size}")
        print(f"Number of observations used for testing: {len(y) - initial_train_size}")
        print(f"    Number of folds: {folds - 1 * (remainder == 0)}")
        print(f"    Number of steps per fold: {steps}")
        if remainder != 0:
            print(f"    Last fold only includes {remainder} observations")
    
    for i in range(folds):
        
        last_window_end   = initial_train_size + i * steps
        last_window_start = (initial_train_size + i * steps) - window_size 
        last_window       = y[last_window_start:last_window_end]
                
        if i < folds - 1:
            if exog is None:
                pred = forecaster.predict_interval(
                            steps       = steps,
                            last_window = last_window,
                            interval    = interval,
                            n_boot      = n_boot,
                            in_sample_residuals = in_sample_residuals
                        )
            else:
                pred = forecaster.predict_interval(
                            steps       = steps,
                            last_window = last_window,
                            exog        = exog[last_window_end:last_window_end + steps],
                            interval    = interval,
                            n_boot      = n_boot,
                            in_sample_residuals = in_sample_residuals
                        )
        elif remainder != 0:
            steps = remainder 
            if exog is None:
                pred = forecaster.predict_interval(
                            steps       = steps,
                            last_window = last_window,
                            interval    = interval,
                            n_boot      = n_boot,
                            in_sample_residuals = in_sample_residuals
                        )
            else:
                pred = forecaster.predict_interval(
                            steps       = steps,
                            last_window = last_window,
                            exog        = exog[last_window_end:last_window_end + steps],
                            interval    = interval,
                            n_boot      = n_boot,
                            in_sample_residuals = in_sample_residuals
                        )
        else:
            continue
        
        backtest_predictions.append(pred)
    
    backtest_predictions = np.concatenate(backtest_predictions)
    metric_value = metric(
                        y_true = y[initial_train_size:],
                        y_pred = backtest_predictions[:, 0]
                   )

    return np.array([metric_value]), backtest_predictions
f57437e96df67467652493c5955df54619ee848f
<mask token>


def time_series_spliter(y: Union[np.ndarray, pd.Series], initial_train_size:
    int, steps: int, allow_incomplete_fold: bool=True, verbose: bool=True):
    """
    
    Split indices of a time series into multiple train-test pairs. The order of
    is maintained and the training set increases in each iteration.
    
    Parameters
    ----------        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    verbose : bool, default `True`
        Print number of splits created.

    Yields
    ------
    train : 1D np.ndarray
        Training indices.
        
    test : 1D np.ndarray
        Test indices.
        
    """
    if not isinstance(y, (np.ndarray, pd.Series)):
        raise Exception('`y` must be `1D np.ndarray` o `pd.Series`.')
    elif isinstance(y, np.ndarray) and y.ndim != 1:
        raise Exception(
            f'`y` must be `1D np.ndarray` o `pd.Series`, got `np.ndarray` with {y.ndim} dimensions.'
            )
    if isinstance(y, pd.Series):
        y = y.to_numpy().copy()
    folds = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    if verbose:
        if folds == 1:
            print(f'Number of folds: {folds - 1}')
            print(
                'Not enought observations in `y` to create to create even a complete fold.'
                )
        elif remainder == 0:
            print(f'Number of folds: {folds - 1}')
        elif remainder != 0 and allow_incomplete_fold:
            print(f'Number of folds: {folds}')
            print(
                f'Since `allow_incomplete_fold=True`, last fold only includes {remainder} observations instead of {steps}.'
                )
            print(
                'Incomplete folds with few observations could overestimate or '
                , 'underestimate validation metrics.')
        elif remainder != 0 and not allow_incomplete_fold:
            print(f'Number of folds: {folds - 1}')
            print(
                f'Since `allow_incomplete_fold=False`, last {remainder} observations are descarted.'
                )
    if folds == 1:
        return []
    for i in range(folds):
        if i < folds - 1:
            train_end = initial_train_size + i * steps
            train_indices = range(train_end)
            test_indices = range(train_end, train_end + steps)
        elif remainder != 0 and allow_incomplete_fold:
            train_end = initial_train_size + i * steps
            train_indices = range(train_end)
            test_indices = range(train_end, len(y))
        else:
            break
        yield train_indices, test_indices


def cv_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
    initial_train_size: int, steps: int, metric: str, exog: Union[np.
    ndarray, pd.Series]=None, allow_incomplete_fold: bool=True, verbose:
    bool=True) ->np.array:
    """
    Cross-validation of `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput`
    object. The order of is maintained and the training set increases in each iteration.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
            
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    verbose : bool, default `True`
        Print number of folds used for cross validation.

    Returns 
    -------
    cv_results: 1D np.ndarray
        Value of the metric for each partition.

    """
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    if metric not in ['mean_squared_error', 'mean_absolute_error',
        'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and 'mean_absolute_percentage_error'. Got {metric}."
            )
    cv_results = []
    metrics = {'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error}
    metric = metrics[metric]
    splits = time_series_spliter(y=y, initial_train_size=initial_train_size,
        steps=steps, allow_incomplete_fold=allow_incomplete_fold, verbose=
        verbose)
    for train_index, test_index in splits:
        if exog is None:
            forecaster.fit(y=y[train_index])
            pred = forecaster.predict(steps=len(test_index))
        else:
            forecaster.fit(y=y[train_index], exog=exog[train_index])
            pred = forecaster.predict(steps=len(test_index), exog=exog[
                test_index])
        metric_value = metric(y_true=y[test_index], y_pred=pred)
        cv_results.append(metric_value)
    return np.array(cv_results)


def backtesting_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
    initial_train_size: int, steps: int, metric: str, exog: Union[np.
    ndarray, pd.Series]=None, verbose: bool=False) ->Tuple[np.array, np.array]:
    """
    Backtesting (validation) of `ForecasterAutoreg`, `ForecasterCustom` or
    `ForecasterAutoregMultiOutput` object.
    The model is trained only once using the `initial_train_size` first observations.
    In each iteration, a number of `steps` predictions are evaluated.
    
    This evaluation is much faster than `cv_forecaster()` since the model is
    trained only once.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
            
    verbose : bool, default `False`
        Print number of folds used for backtesting.

    Returns 
    -------
    backtest_predictions: 1D np.ndarray
        Value of predictions.
        
    metric_value: np.ndarray shape (1,)
        Value of the metric.

    """
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    if metric not in ['mean_squared_error', 'mean_absolute_error',
        'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and 'mean_absolute_percentage_error'. Got {metric}."
            )
    backtest_predictions = []
    metrics = {'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error}
    metric = metrics[metric]
    if exog is None:
        forecaster.fit(y=y[:initial_train_size])
    else:
        forecaster.fit(y=y[:initial_train_size], exog=exog[:initial_train_size]
            )
    folds = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    window_size = len(forecaster.last_window)
    if isinstance(forecaster, ForecasterAutoregMultiOutput) and remainder != 0:
        logging.warning(
            f'Backtesting `ForecasterAutoregMultiOutput` only allow completed folds. Last {remainder} observations are excluded.'
            )
        remainder = 0
    if verbose:
        print(f'Number of observations used for training: {initial_train_size}'
            )
        print(
            f'Number of observations used for testing: {len(y) - initial_train_size}'
            )
        print(f'    Number of folds: {folds - 1 * (remainder == 0)}')
        print(f'    Number of steps per fold: {steps}')
        if remainder != 0:
            print(f'    Last fold only includes {remainder} observations')
    for i in range(folds):
        last_window_end = initial_train_size + i * steps
        last_window_start = initial_train_size + i * steps - window_size
        last_window = y[last_window_start:last_window_end]
        if i < folds - 1:
            if exog is None:
                pred = forecaster.predict(steps=steps, last_window=last_window)
            else:
                pred = forecaster.predict(steps=steps, last_window=
                    last_window, exog=exog[last_window_end:last_window_end +
                    steps])
        elif remainder != 0:
            steps = remainder
            if exog is None:
                pred = forecaster.predict(steps=steps, last_window=last_window)
            else:
                pred = forecaster.predict(steps=steps, last_window=
                    last_window, exog=exog[last_window_end:last_window_end +
                    steps])
        else:
            continue
        backtest_predictions.append(pred)
    backtest_predictions = np.concatenate(backtest_predictions)
    metric_value = metric(y_true=y[initial_train_size:initial_train_size +
        len(backtest_predictions)], y_pred=backtest_predictions)
    return np.array([metric_value]), backtest_predictions


<mask token>


def backtesting_forecaster_intervals(forecaster, y: Union[np.ndarray, pd.
    Series], initial_train_size: int, steps: int, metric: str, exog: Union[
    np.ndarray, pd.Series]=None, interval: list=[5, 95], n_boot: int=500,
    in_sample_residuals: bool=True, verbose: bool=False) ->Tuple[np.array,
    np.array]:
    """
    Backtesting (validation) of `ForecasterAutoreg` or `ForecasterCustom` object.
    The model is trained only once using the `initial_train_size` first observations.
    In each iteration, a number of `steps` predictions are evaluated. Both, 
    predictions and intervals, are calculated.
    
    This evaluation is much faster than `cv_forecaster()` since the model is
    trained only once.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom
        `ForecasterAutoreg` or `ForecasterCustom` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
        
    interval: list, default `[5, 100]`
            Confidence of the prediction interval estimated. Sequence of percentiles
            to compute, which must be between 0 and 100 inclusive.
            
    n_boot: int, default `500`
        Number of bootstrapping iterations used to estimate prediction
        intervals.

    in_sample_residuals: bool, default `True`
        If `True`, residuals from the training data are used as proxy of
        prediction error to create prediction intervals.
        
    verbose : bool, default `True`
        Print number of folds used for backtesting.

    Returns 
    -------
    backtest_predictions: np.ndarray
        If `include_intervals=True`, 2D np.ndarray shape(steps, 3) with predicted
        value and their estimated interval.
            Column 0 = predictions
            Column 1 = lower bound interval
            Column 2 = upper bound interval
        
    metric_value: np.ndarray shape (1,)
        Value of the metric.

    Notes
    -----
    More information about prediction intervals in forecasting:
    https://otexts.com/fpp2/prediction-intervals.html
    Forecasting: Principles and Practice (2nd ed) Rob J Hyndman and
    George Athanasopoulos.

    """
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    if metric not in ['mean_squared_error', 'mean_absolute_error',
        'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and 'mean_absolute_percentage_error'. Got {metric}."
            )
    backtest_predictions = []
    metrics = {'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error}
    metric = metrics[metric]
    if exog is None:
        forecaster.fit(y=y[:initial_train_size])
    else:
        forecaster.fit(y=y[:initial_train_size], exog=exog[:initial_train_size]
            )
    folds = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    window_size = len(forecaster.last_window)
    if verbose:
        print(f'Number of observations used for training: {initial_train_size}'
            )
        print(
            f'Number of observations used for testing: {len(y) - initial_train_size}'
            )
        print(f'    Number of folds: {folds - 1 * (remainder == 0)}')
        print(f'    Number of steps per fold: {steps}')
        if remainder != 0:
            print(f'    Last fold only includes {remainder} observations')
    for i in range(folds):
        last_window_end = initial_train_size + i * steps
        last_window_start = initial_train_size + i * steps - window_size
        last_window = y[last_window_start:last_window_end]
        if i < folds - 1:
            if exog is None:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
            else:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, exog=exog[last_window_end:last_window_end +
                    steps], interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
        elif remainder != 0:
            steps = remainder
            if exog is None:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
            else:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, exog=exog[last_window_end:last_window_end +
                    steps], interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
        else:
            continue
        backtest_predictions.append(pred)
    backtest_predictions = np.concatenate(backtest_predictions)
    metric_value = metric(y_true=y[initial_train_size:], y_pred=
        backtest_predictions[:, 0])
    return np.array([metric_value]), backtest_predictions

<mask token>


def time_series_spliter(y: Union[np.ndarray, pd.Series], initial_train_size:
    int, steps: int, allow_incomplete_fold: bool=True, verbose: bool=True):
    """
    
    Split indices of a time series into multiple train-test pairs. The order of
    is maintained and the training set increases in each iteration.
    
    Parameters
    ----------        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    verbose : bool, default `True`
        Print number of splits created.

    Yields
    ------
    train : 1D np.ndarray
        Training indices.
        
    test : 1D np.ndarray
        Test indices.
        
    """
    if not isinstance(y, (np.ndarray, pd.Series)):
        raise Exception('`y` must be `1D np.ndarray` o `pd.Series`.')
    elif isinstance(y, np.ndarray) and y.ndim != 1:
        raise Exception(
            f'`y` must be `1D np.ndarray` o `pd.Series`, got `np.ndarray` with {y.ndim} dimensions.'
            )
    if isinstance(y, pd.Series):
        y = y.to_numpy().copy()
    folds = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    if verbose:
        if folds == 1:
            print(f'Number of folds: {folds - 1}')
            print(
                'Not enought observations in `y` to create to create even a complete fold.'
                )
        elif remainder == 0:
            print(f'Number of folds: {folds - 1}')
        elif remainder != 0 and allow_incomplete_fold:
            print(f'Number of folds: {folds}')
            print(
                f'Since `allow_incomplete_fold=True`, last fold only includes {remainder} observations instead of {steps}.'
                )
            print(
                'Incomplete folds with few observations could overestimate or '
                , 'underestimate validation metrics.')
        elif remainder != 0 and not allow_incomplete_fold:
            print(f'Number of folds: {folds - 1}')
            print(
                f'Since `allow_incomplete_fold=False`, last {remainder} observations are descarted.'
                )
    if folds == 1:
        return []
    for i in range(folds):
        if i < folds - 1:
            train_end = initial_train_size + i * steps
            train_indices = range(train_end)
            test_indices = range(train_end, train_end + steps)
        elif remainder != 0 and allow_incomplete_fold:
            train_end = initial_train_size + i * steps
            train_indices = range(train_end)
            test_indices = range(train_end, len(y))
        else:
            break
        yield train_indices, test_indices


def cv_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
    initial_train_size: int, steps: int, metric: str, exog: Union[np.
    ndarray, pd.Series]=None, allow_incomplete_fold: bool=True, verbose:
    bool=True) ->np.array:
    """
    Cross-validation of `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput`
    object. The order of is maintained and the training set increases in each iteration.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
            
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    verbose : bool, default `True`
        Print number of folds used for cross validation.

    Returns 
    -------
    cv_results: 1D np.ndarray
        Value of the metric for each partition.

    """
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    if metric not in ['mean_squared_error', 'mean_absolute_error',
        'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and 'mean_absolute_percentage_error'. Got {metric}."
            )
    cv_results = []
    metrics = {'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error}
    metric = metrics[metric]
    splits = time_series_spliter(y=y, initial_train_size=initial_train_size,
        steps=steps, allow_incomplete_fold=allow_incomplete_fold, verbose=
        verbose)
    for train_index, test_index in splits:
        if exog is None:
            forecaster.fit(y=y[train_index])
            pred = forecaster.predict(steps=len(test_index))
        else:
            forecaster.fit(y=y[train_index], exog=exog[train_index])
            pred = forecaster.predict(steps=len(test_index), exog=exog[
                test_index])
        metric_value = metric(y_true=y[test_index], y_pred=pred)
        cv_results.append(metric_value)
    return np.array(cv_results)


def backtesting_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
    initial_train_size: int, steps: int, metric: str, exog: Union[np.
    ndarray, pd.Series]=None, verbose: bool=False) ->Tuple[np.array, np.array]:
    """
    Backtesting (validation) of `ForecasterAutoreg`, `ForecasterCustom` or
    `ForecasterAutoregMultiOutput` object.
    The model is trained only once using the `initial_train_size` first observations.
    In each iteration, a number of `steps` predictions are evaluated.
    
    This evaluation is much faster than `cv_forecaster()` since the model is
    trained only once.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
            
    verbose : bool, default `False`
        Print number of folds used for backtesting.

    Returns 
    -------
    backtest_predictions: 1D np.ndarray
        Value of predictions.
        
    metric_value: np.ndarray shape (1,)
        Value of the metric.

    """
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    if metric not in ['mean_squared_error', 'mean_absolute_error',
        'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and 'mean_absolute_percentage_error'. Got {metric}."
            )
    backtest_predictions = []
    metrics = {'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error}
    metric = metrics[metric]
    if exog is None:
        forecaster.fit(y=y[:initial_train_size])
    else:
        forecaster.fit(y=y[:initial_train_size], exog=exog[:initial_train_size]
            )
    folds = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    window_size = len(forecaster.last_window)
    if isinstance(forecaster, ForecasterAutoregMultiOutput) and remainder != 0:
        logging.warning(
            f'Backtesting `ForecasterAutoregMultiOutput` only allow completed folds. Last {remainder} observations are excluded.'
            )
        remainder = 0
    if verbose:
        print(f'Number of observations used for training: {initial_train_size}'
            )
        print(
            f'Number of observations used for testing: {len(y) - initial_train_size}'
            )
        print(f'    Number of folds: {folds - 1 * (remainder == 0)}')
        print(f'    Number of steps per fold: {steps}')
        if remainder != 0:
            print(f'    Last fold only includes {remainder} observations')
    for i in range(folds):
        last_window_end = initial_train_size + i * steps
        last_window_start = initial_train_size + i * steps - window_size
        last_window = y[last_window_start:last_window_end]
        if i < folds - 1:
            if exog is None:
                pred = forecaster.predict(steps=steps, last_window=last_window)
            else:
                pred = forecaster.predict(steps=steps, last_window=
                    last_window, exog=exog[last_window_end:last_window_end +
                    steps])
        elif remainder != 0:
            steps = remainder
            if exog is None:
                pred = forecaster.predict(steps=steps, last_window=last_window)
            else:
                pred = forecaster.predict(steps=steps, last_window=
                    last_window, exog=exog[last_window_end:last_window_end +
                    steps])
        else:
            continue
        backtest_predictions.append(pred)
    backtest_predictions = np.concatenate(backtest_predictions)
    metric_value = metric(y_true=y[initial_train_size:initial_train_size +
        len(backtest_predictions)], y_pred=backtest_predictions)
    return np.array([metric_value]), backtest_predictions


def grid_search_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
    param_grid: dict, initial_train_size: int, steps: int, metric: str,
    exog: Union[np.ndarray, pd.Series]=None, lags_grid: list=None, method:
    str='cv', allow_incomplete_fold: bool=True, return_best: bool=True,
    verbose: bool=True) ->pd.DataFrame:
    """
    Exhaustive search over specified parameter values for a Forecaster object.
    Validation is done using time series cross-validation or backtesting.
    
    Parameters
    ----------
    
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        ForecasterAutoreg, ForecasterCustom or ForecasterAutoregMultiOutput object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
        
    param_grid : dict
        Dictionary with parameters names (`str`) as keys and lists of parameter
        settings to try as values.
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
           
    lags_grid : list of int, lists, np.narray or range. 
        Lists of `lags` to try. Only used if forecaster is an instance of 
        `ForecasterCustom`.
        
    method : {'cv', 'backtesting'}
        Method used to estimate the metric for each parameter combination.
        'cv' for time series crosvalidation and 'backtesting' for simple
        backtesting. 'backtesting' is much faster since the model is fitted only
        once.
        
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    return_best : bool
        Refit the `forecaster` using the best found parameters on the whole data.
        
    verbose : bool, default `True`
        Print number of folds used for cv or backtesting.

    Returns 
    -------
    results: pandas.DataFrame
        Metric value estimated for each combination of parameters.

    """
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    if isinstance(forecaster, ForecasterCustom):
        if lags_grid is not None:
            logging.warning(
                '`lags_grid` ignored if forecaster is an instance of `ForecasterCustom`.'
                )
        lags_grid = ['custom predictors']
    elif lags_grid is None:
        lags_grid = [forecaster.lags]
    if isinstance(forecaster, ForecasterAutoregMultiOutput):
        param_grid_renamed = {('estimator__' + key): value for key, value in
            param_grid.items()}
        param_grid = param_grid_renamed
    lags_list = []
    params_list = []
    metric_list = []
    param_grid = list(ParameterGrid(param_grid))
    logging.info(f'Number of models to fit: {len(param_grid) * len(lags_grid)}'
        )
    for lags in tqdm.tqdm(lags_grid, desc='loop lags_grid', position=0):
        if isinstance(forecaster, (ForecasterAutoreg,
            ForecasterAutoregMultiOutput)):
            forecaster.set_lags(lags)
            lags = forecaster.lags.copy()
        for params in tqdm.tqdm(param_grid, desc='loop param_grid',
            position=1, leave=False):
            forecaster.set_params(**params)
            if method == 'cv':
                metrics = cv_forecaster(forecaster=forecaster, y=y, exog=
                    exog, initial_train_size=initial_train_size, steps=
                    steps, metric=metric, allow_incomplete_fold=
                    allow_incomplete_fold, verbose=verbose)
            else:
                metrics = backtesting_forecaster(forecaster=forecaster, y=y,
                    exog=exog, initial_train_size=initial_train_size, steps
                    =steps, metric=metric, verbose=verbose)[0]
            lags_list.append(lags)
            params_list.append(params)
            metric_list.append(metrics.mean())
    results = pd.DataFrame({'lags': lags_list, 'params': params_list,
        'metric': metric_list})
    results = results.sort_values(by='metric', ascending=True)
    if return_best:
        best_lags = results['lags'].iloc[0]
        best_params = results['params'].iloc[0]
        logging.info(
            f"""Refitting `forecaster` using the best found parameters and the whole data set: 
lags: {best_lags} 
params: {best_params}
"""
            )
        if isinstance(forecaster, (ForecasterAutoreg,
            ForecasterAutoregMultiOutput)):
            forecaster.set_lags(best_lags)
        forecaster.set_params(**best_params)
        forecaster.fit(y=y, exog=exog)
    return results


def backtesting_forecaster_intervals(forecaster, y: Union[np.ndarray, pd.
    Series], initial_train_size: int, steps: int, metric: str, exog: Union[
    np.ndarray, pd.Series]=None, interval: list=[5, 95], n_boot: int=500,
    in_sample_residuals: bool=True, verbose: bool=False) ->Tuple[np.array,
    np.array]:
    """
    Backtesting (validation) of `ForecasterAutoreg` or `ForecasterCustom` object.
    The model is trained only once using the `initial_train_size` first observations.
    In each iteration, a number of `steps` predictions are evaluated. Both, 
    predictions and intervals, are calculated.
    
    This evaluation is much faster than `cv_forecaster()` since the model is
    trained only once.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom
        `ForecasterAutoreg` or `ForecasterCustom` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
        
    interval: list, default `[5, 100]`
            Confidence of the prediction interval estimated. Sequence of percentiles
            to compute, which must be between 0 and 100 inclusive.
            
    n_boot: int, default `500`
        Number of bootstrapping iterations used to estimate prediction
        intervals.

    in_sample_residuals: bool, default `True`
        If `True`, residuals from the training data are used as proxy of
        prediction error to create prediction intervals.
        
    verbose : bool, default `True`
        Print number of folds used for backtesting.

    Returns 
    -------
    backtest_predictions: np.ndarray
        If `include_intervals=True`, 2D np.ndarray shape(steps, 3) with predicted
        value and their estimated interval.
            Column 0 = predictions
            Column 1 = lower bound interval
            Column 2 = upper bound interval
        
    metric_value: np.ndarray shape (1,)
        Value of the metric.

    Notes
    -----
    More information about prediction intervals in forecasting:
    https://otexts.com/fpp2/prediction-intervals.html
    Forecasting: Principles and Practice (2nd ed) Rob J Hyndman and
    George Athanasopoulos.

    """
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    if metric not in ['mean_squared_error', 'mean_absolute_error',
        'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and 'mean_absolute_percentage_error'. Got {metric}."
            )
    backtest_predictions = []
    metrics = {'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error}
    metric = metrics[metric]
    if exog is None:
        forecaster.fit(y=y[:initial_train_size])
    else:
        forecaster.fit(y=y[:initial_train_size], exog=exog[:initial_train_size]
            )
    folds = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    window_size = len(forecaster.last_window)
    if verbose:
        print(f'Number of observations used for training: {initial_train_size}'
            )
        print(
            f'Number of observations used for testing: {len(y) - initial_train_size}'
            )
        print(f'    Number of folds: {folds - 1 * (remainder == 0)}')
        print(f'    Number of steps per fold: {steps}')
        if remainder != 0:
            print(f'    Last fold only includes {remainder} observations')
    for i in range(folds):
        last_window_end = initial_train_size + i * steps
        last_window_start = initial_train_size + i * steps - window_size
        last_window = y[last_window_start:last_window_end]
        if i < folds - 1:
            if exog is None:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
            else:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, exog=exog[last_window_end:last_window_end +
                    steps], interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
        elif remainder != 0:
            steps = remainder
            if exog is None:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
            else:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, exog=exog[last_window_end:last_window_end +
                    steps], interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
        else:
            continue
        backtest_predictions.append(pred)
    backtest_predictions = np.concatenate(backtest_predictions)
    metric_value = metric(y_true=y[initial_train_size:], y_pred=
        backtest_predictions[:, 0])
    return np.array([metric_value]), backtest_predictions

<mask token>
logging.basicConfig(format=
    '%(asctime)-5s %(name)-10s %(levelname)-5s %(message)s', level=logging.INFO
    )


def time_series_spliter(y: Union[np.ndarray, pd.Series], initial_train_size:
    int, steps: int, allow_incomplete_fold: bool=True, verbose: bool=True):
    """
    
    Split indices of a time series into multiple train-test pairs. The order of
    is maintained and the training set increases in each iteration.
    
    Parameters
    ----------        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    verbose : bool, default `True`
        Print number of splits created.

    Yields
    ------
    train : 1D np.ndarray
        Training indices.
        
    test : 1D np.ndarray
        Test indices.
        
    """
    if not isinstance(y, (np.ndarray, pd.Series)):
        raise Exception('`y` must be `1D np.ndarray` o `pd.Series`.')
    elif isinstance(y, np.ndarray) and y.ndim != 1:
        raise Exception(
            f'`y` must be `1D np.ndarray` o `pd.Series`, got `np.ndarray` with {y.ndim} dimensions.'
            )
    if isinstance(y, pd.Series):
        y = y.to_numpy().copy()
    folds = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    if verbose:
        if folds == 1:
            print(f'Number of folds: {folds - 1}')
            print(
                'Not enought observations in `y` to create to create even a complete fold.'
                )
        elif remainder == 0:
            print(f'Number of folds: {folds - 1}')
        elif remainder != 0 and allow_incomplete_fold:
            print(f'Number of folds: {folds}')
            print(
                f'Since `allow_incomplete_fold=True`, last fold only includes {remainder} observations instead of {steps}.'
                )
            print(
                'Incomplete folds with few observations could overestimate or '
                , 'underestimate validation metrics.')
        elif remainder != 0 and not allow_incomplete_fold:
            print(f'Number of folds: {folds - 1}')
            print(
                f'Since `allow_incomplete_fold=False`, last {remainder} observations are descarted.'
                )
    if folds == 1:
        return []
    for i in range(folds):
        if i < folds - 1:
            train_end = initial_train_size + i * steps
            train_indices = range(train_end)
            test_indices = range(train_end, train_end + steps)
        elif remainder != 0 and allow_incomplete_fold:
            train_end = initial_train_size + i * steps
            train_indices = range(train_end)
            test_indices = range(train_end, len(y))
        else:
            break
        yield train_indices, test_indices


def cv_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
    initial_train_size: int, steps: int, metric: str, exog: Union[np.
    ndarray, pd.Series]=None, allow_incomplete_fold: bool=True, verbose:
    bool=True) ->np.array:
    """
    Cross-validation of `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput`
    object. The order of is maintained and the training set increases in each iteration.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
            
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    verbose : bool, default `True`
        Print number of folds used for cross validation.

    Returns 
    -------
    cv_results: 1D np.ndarray
        Value of the metric for each partition.

    """
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    if metric not in ['mean_squared_error', 'mean_absolute_error',
        'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and 'mean_absolute_percentage_error'. Got {metric}."
            )
    cv_results = []
    metrics = {'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error}
    metric = metrics[metric]
    splits = time_series_spliter(y=y, initial_train_size=initial_train_size,
        steps=steps, allow_incomplete_fold=allow_incomplete_fold, verbose=
        verbose)
    for train_index, test_index in splits:
        if exog is None:
            forecaster.fit(y=y[train_index])
            pred = forecaster.predict(steps=len(test_index))
        else:
            forecaster.fit(y=y[train_index], exog=exog[train_index])
            pred = forecaster.predict(steps=len(test_index), exog=exog[
                test_index])
        metric_value = metric(y_true=y[test_index], y_pred=pred)
        cv_results.append(metric_value)
    return np.array(cv_results)


def backtesting_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
    initial_train_size: int, steps: int, metric: str, exog: Union[np.
    ndarray, pd.Series]=None, verbose: bool=False) ->Tuple[np.array, np.array]:
    """
    Backtesting (validation) of `ForecasterAutoreg`, `ForecasterCustom` or
    `ForecasterAutoregMultiOutput` object.
    The model is trained only once using the `initial_train_size` first observations.
    In each iteration, a number of `steps` predictions are evaluated.
    
    This evaluation is much faster than `cv_forecaster()` since the model is
    trained only once.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
            
    verbose : bool, default `False`
        Print number of folds used for backtesting.

    Returns 
    -------
    backtest_predictions: 1D np.ndarray
        Value of predictions.
        
    metric_value: np.ndarray shape (1,)
        Value of the metric.

    """
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    if metric not in ['mean_squared_error', 'mean_absolute_error',
        'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and 'mean_absolute_percentage_error'. Got {metric}."
            )
    backtest_predictions = []
    metrics = {'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error}
    metric = metrics[metric]
    if exog is None:
        forecaster.fit(y=y[:initial_train_size])
    else:
        forecaster.fit(y=y[:initial_train_size], exog=exog[:initial_train_size]
            )
    folds = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    window_size = len(forecaster.last_window)
    if isinstance(forecaster, ForecasterAutoregMultiOutput) and remainder != 0:
        logging.warning(
            f'Backtesting `ForecasterAutoregMultiOutput` only allow completed folds. Last {remainder} observations are excluded.'
            )
        remainder = 0
    if verbose:
        print(f'Number of observations used for training: {initial_train_size}'
            )
        print(
            f'Number of observations used for testing: {len(y) - initial_train_size}'
            )
        print(f'    Number of folds: {folds - 1 * (remainder == 0)}')
        print(f'    Number of steps per fold: {steps}')
        if remainder != 0:
            print(f'    Last fold only includes {remainder} observations')
    for i in range(folds):
        last_window_end = initial_train_size + i * steps
        last_window_start = initial_train_size + i * steps - window_size
        last_window = y[last_window_start:last_window_end]
        if i < folds - 1:
            if exog is None:
                pred = forecaster.predict(steps=steps, last_window=last_window)
            else:
                pred = forecaster.predict(steps=steps, last_window=
                    last_window, exog=exog[last_window_end:last_window_end +
                    steps])
        elif remainder != 0:
            steps = remainder
            if exog is None:
                pred = forecaster.predict(steps=steps, last_window=last_window)
            else:
                pred = forecaster.predict(steps=steps, last_window=
                    last_window, exog=exog[last_window_end:last_window_end +
                    steps])
        else:
            continue
        backtest_predictions.append(pred)
    backtest_predictions = np.concatenate(backtest_predictions)
    metric_value = metric(y_true=y[initial_train_size:initial_train_size +
        len(backtest_predictions)], y_pred=backtest_predictions)
    return np.array([metric_value]), backtest_predictions


def grid_search_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
    param_grid: dict, initial_train_size: int, steps: int, metric: str,
    exog: Union[np.ndarray, pd.Series]=None, lags_grid: list=None, method:
    str='cv', allow_incomplete_fold: bool=True, return_best: bool=True,
    verbose: bool=True) ->pd.DataFrame:
    """
    Exhaustive search over specified parameter values for a Forecaster object.
    Validation is done using time series cross-validation or backtesting.
    
    Parameters
    ----------
    
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        ForecasterAutoreg, ForecasterCustom or ForecasterAutoregMultiOutput object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
        
    param_grid : dict
        Dictionary with parameters names (`str`) as keys and lists of parameter
        settings to try as values.
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
           
    lags_grid : list of int, lists, np.narray or range. 
        Lists of `lags` to try. Only used if forecaster is an instance of 
        `ForecasterCustom`.
        
    method : {'cv', 'backtesting'}
        Method used to estimate the metric for each parameter combination.
        'cv' for time series crosvalidation and 'backtesting' for simple
        backtesting. 'backtesting' is much faster since the model is fitted only
        once.
        
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    return_best : bool
        Refit the `forecaster` using the best found parameters on the whole data.
        
    verbose : bool, default `True`
        Print number of folds used for cv or backtesting.

    Returns 
    -------
    results: pandas.DataFrame
        Metric value estimated for each combination of parameters.

    """
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    if isinstance(forecaster, ForecasterCustom):
        if lags_grid is not None:
            logging.warning(
                '`lags_grid` ignored if forecaster is an instance of `ForecasterCustom`.'
                )
        lags_grid = ['custom predictors']
    elif lags_grid is None:
        lags_grid = [forecaster.lags]
    if isinstance(forecaster, ForecasterAutoregMultiOutput):
        param_grid_renamed = {('estimator__' + key): value for key, value in
            param_grid.items()}
        param_grid = param_grid_renamed
    lags_list = []
    params_list = []
    metric_list = []
    param_grid = list(ParameterGrid(param_grid))
    logging.info(f'Number of models to fit: {len(param_grid) * len(lags_grid)}'
        )
    for lags in tqdm.tqdm(lags_grid, desc='loop lags_grid', position=0):
        if isinstance(forecaster, (ForecasterAutoreg,
            ForecasterAutoregMultiOutput)):
            forecaster.set_lags(lags)
            lags = forecaster.lags.copy()
        for params in tqdm.tqdm(param_grid, desc='loop param_grid',
            position=1, leave=False):
            forecaster.set_params(**params)
            if method == 'cv':
                metrics = cv_forecaster(forecaster=forecaster, y=y, exog=
                    exog, initial_train_size=initial_train_size, steps=
                    steps, metric=metric, allow_incomplete_fold=
                    allow_incomplete_fold, verbose=verbose)
            else:
                metrics = backtesting_forecaster(forecaster=forecaster, y=y,
                    exog=exog, initial_train_size=initial_train_size, steps
                    =steps, metric=metric, verbose=verbose)[0]
            lags_list.append(lags)
            params_list.append(params)
            metric_list.append(metrics.mean())
    results = pd.DataFrame({'lags': lags_list, 'params': params_list,
        'metric': metric_list})
    results = results.sort_values(by='metric', ascending=True)
    if return_best:
        best_lags = results['lags'].iloc[0]
        best_params = results['params'].iloc[0]
        logging.info(
            f"""Refitting `forecaster` using the best found parameters and the whole data set: 
lags: {best_lags} 
params: {best_params}
"""
            )
        if isinstance(forecaster, (ForecasterAutoreg,
            ForecasterAutoregMultiOutput)):
            forecaster.set_lags(best_lags)
        forecaster.set_params(**best_params)
        forecaster.fit(y=y, exog=exog)
    return results


def backtesting_forecaster_intervals(forecaster, y: Union[np.ndarray, pd.
    Series], initial_train_size: int, steps: int, metric: str, exog: Union[
    np.ndarray, pd.Series]=None, interval: list=[5, 95], n_boot: int=500,
    in_sample_residuals: bool=True, verbose: bool=False) ->Tuple[np.array,
    np.array]:
    """
    Backtesting (validation) of `ForecasterAutoreg` or `ForecasterCustom` object.
    The model is trained only once using the `initial_train_size` first observations.
    In each iteration, a number of `steps` predictions are evaluated. Both, 
    predictions and intervals, are calculated.
    
    This evaluation is much faster than `cv_forecaster()` since the model is
    trained only once.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom
        `ForecasterAutoreg` or `ForecasterCustom` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
        
    interval: list, default `[5, 100]`
            Confidence of the prediction interval estimated. Sequence of percentiles
            to compute, which must be between 0 and 100 inclusive.
            
    n_boot: int, default `500`
        Number of bootstrapping iterations used to estimate prediction
        intervals.

    in_sample_residuals: bool, default `True`
        If `True`, residuals from the training data are used as proxy of
        prediction error to create prediction intervals.
        
    verbose : bool, default `True`
        Print number of folds used for backtesting.

    Returns 
    -------
    backtest_predictions: np.ndarray
        If `include_intervals=True`, 2D np.ndarray shape(steps, 3) with predicted
        value and their estimated interval.
            Column 0 = predictions
            Column 1 = lower bound interval
            Column 2 = upper bound interval
        
    metric_value: np.ndarray shape (1,)
        Value of the metric.

    Notes
    -----
    More information about prediction intervals in forecasting:
    https://otexts.com/fpp2/prediction-intervals.html
    Forecasting: Principles and Practice (2nd ed) Rob J Hyndman and
    George Athanasopoulos.

    """
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    if metric not in ['mean_squared_error', 'mean_absolute_error',
        'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and 'mean_absolute_percentage_error'. Got {metric}."
            )
    backtest_predictions = []
    metrics = {'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error}
    metric = metrics[metric]
    if exog is None:
        forecaster.fit(y=y[:initial_train_size])
    else:
        forecaster.fit(y=y[:initial_train_size], exog=exog[:initial_train_size]
            )
    folds = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    window_size = len(forecaster.last_window)
    if verbose:
        print(f'Number of observations used for training: {initial_train_size}'
            )
        print(
            f'Number of observations used for testing: {len(y) - initial_train_size}'
            )
        print(f'    Number of folds: {folds - 1 * (remainder == 0)}')
        print(f'    Number of steps per fold: {steps}')
        if remainder != 0:
            print(f'    Last fold only includes {remainder} observations')
    for i in range(folds):
        last_window_end = initial_train_size + i * steps
        last_window_start = initial_train_size + i * steps - window_size
        last_window = y[last_window_start:last_window_end]
        if i < folds - 1:
            if exog is None:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
            else:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, exog=exog[last_window_end:last_window_end +
                    steps], interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
        elif remainder != 0:
            steps = remainder
            if exog is None:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
            else:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, exog=exog[last_window_end:last_window_end +
                    steps], interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
        else:
            continue
        backtest_predictions.append(pred)
    backtest_predictions = np.concatenate(backtest_predictions)
    metric_value = metric(y_true=y[initial_train_size:], y_pred=
        backtest_predictions[:, 0])
    return np.array([metric_value]), backtest_predictions

import typing
from typing import Union, Dict, List, Tuple
import numpy as np
import pandas as pd
import logging
import tqdm
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_absolute_percentage_error
from sklearn.model_selection import ParameterGrid
from .ForecasterAutoreg import ForecasterAutoreg
from .ForecasterCustom import ForecasterCustom
from .ForecasterAutoregMultiOutput import ForecasterAutoregMultiOutput
logging.basicConfig(format=
    '%(asctime)-5s %(name)-10s %(levelname)-5s %(message)s', level=logging.INFO
    )


def time_series_spliter(y: Union[np.ndarray, pd.Series], initial_train_size:
    int, steps: int, allow_incomplete_fold: bool=True, verbose: bool=True):
    """
    
    Split indices of a time series into multiple train-test pairs. The order of
    is maintained and the training set increases in each iteration.
    
    Parameters
    ----------        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    verbose : bool, default `True`
        Print number of splits created.

    Yields
    ------
    train : 1D np.ndarray
        Training indices.
        
    test : 1D np.ndarray
        Test indices.
        
    """
    if not isinstance(y, (np.ndarray, pd.Series)):
        raise Exception('`y` must be `1D np.ndarray` o `pd.Series`.')
    elif isinstance(y, np.ndarray) and y.ndim != 1:
        raise Exception(
            f'`y` must be `1D np.ndarray` o `pd.Series`, got `np.ndarray` with {y.ndim} dimensions.'
            )
    if isinstance(y, pd.Series):
        y = y.to_numpy().copy()
    folds = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    if verbose:
        if folds == 1:
            print(f'Number of folds: {folds - 1}')
            print(
                'Not enought observations in `y` to create to create even a complete fold.'
                )
        elif remainder == 0:
            print(f'Number of folds: {folds - 1}')
        elif remainder != 0 and allow_incomplete_fold:
            print(f'Number of folds: {folds}')
            print(
                f'Since `allow_incomplete_fold=True`, last fold only includes {remainder} observations instead of {steps}.'
                )
            print(
                'Incomplete folds with few observations could overestimate or '
                , 'underestimate validation metrics.')
        elif remainder != 0 and not allow_incomplete_fold:
            print(f'Number of folds: {folds - 1}')
            print(
                f'Since `allow_incomplete_fold=False`, last {remainder} observations are descarted.'
                )
    if folds == 1:
        return []
    for i in range(folds):
        if i < folds - 1:
            train_end = initial_train_size + i * steps
            train_indices = range(train_end)
            test_indices = range(train_end, train_end + steps)
        elif remainder != 0 and allow_incomplete_fold:
            train_end = initial_train_size + i * steps
            train_indices = range(train_end)
            test_indices = range(train_end, len(y))
        else:
            break
        yield train_indices, test_indices


def cv_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
    initial_train_size: int, steps: int, metric: str, exog: Union[np.
    ndarray, pd.Series]=None, allow_incomplete_fold: bool=True, verbose:
    bool=True) ->np.array:
    """
    Cross-validation of `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput`
    object. The order of is maintained and the training set increases in each iteration.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
            
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    verbose : bool, default `True`
        Print number of folds used for cross validation.

    Returns 
    -------
    cv_results: 1D np.ndarray
        Value of the metric for each partition.

    """
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    if metric not in ['mean_squared_error', 'mean_absolute_error',
        'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and 'mean_absolute_percentage_error'. Got {metric}."
            )
    cv_results = []
    metrics = {'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error}
    metric = metrics[metric]
    splits = time_series_spliter(y=y, initial_train_size=initial_train_size,
        steps=steps, allow_incomplete_fold=allow_incomplete_fold, verbose=
        verbose)
    for train_index, test_index in splits:
        if exog is None:
            forecaster.fit(y=y[train_index])
            pred = forecaster.predict(steps=len(test_index))
        else:
            forecaster.fit(y=y[train_index], exog=exog[train_index])
            pred = forecaster.predict(steps=len(test_index), exog=exog[
                test_index])
        metric_value = metric(y_true=y[test_index], y_pred=pred)
        cv_results.append(metric_value)
    return np.array(cv_results)


def backtesting_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
    initial_train_size: int, steps: int, metric: str, exog: Union[np.
    ndarray, pd.Series]=None, verbose: bool=False) ->Tuple[np.array, np.array]:
    """
    Backtesting (validation) of `ForecasterAutoreg`, `ForecasterCustom` or
    `ForecasterAutoregMultiOutput` object.
    The model is trained only once using the `initial_train_size` first observations.
    In each iteration, a number of `steps` predictions are evaluated.
    
    This evaluation is much faster than `cv_forecaster()` since the model is
    trained only once.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
            
    verbose : bool, default `False`
        Print number of folds used for backtesting.

    Returns 
    -------
    backtest_predictions: 1D np.ndarray
        Value of predictions.
        
    metric_value: np.ndarray shape (1,)
        Value of the metric.

    """
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    if metric not in ['mean_squared_error', 'mean_absolute_error',
        'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and 'mean_absolute_percentage_error'. Got {metric}."
            )
    backtest_predictions = []
    metrics = {'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error}
    metric = metrics[metric]
    if exog is None:
        forecaster.fit(y=y[:initial_train_size])
    else:
        forecaster.fit(y=y[:initial_train_size], exog=exog[:initial_train_size]
            )
    folds = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    window_size = len(forecaster.last_window)
    if isinstance(forecaster, ForecasterAutoregMultiOutput) and remainder != 0:
        logging.warning(
            f'Backtesting `ForecasterAutoregMultiOutput` only allow completed folds. Last {remainder} observations are excluded.'
            )
        remainder = 0
    if verbose:
        print(f'Number of observations used for training: {initial_train_size}'
            )
        print(
            f'Number of observations used for testing: {len(y) - initial_train_size}'
            )
        print(f'    Number of folds: {folds - 1 * (remainder == 0)}')
        print(f'    Number of steps per fold: {steps}')
        if remainder != 0:
            print(f'    Last fold only includes {remainder} observations')
    for i in range(folds):
        last_window_end = initial_train_size + i * steps
        last_window_start = initial_train_size + i * steps - window_size
        last_window = y[last_window_start:last_window_end]
        if i < folds - 1:
            if exog is None:
                pred = forecaster.predict(steps=steps, last_window=last_window)
            else:
                pred = forecaster.predict(steps=steps, last_window=
                    last_window, exog=exog[last_window_end:last_window_end +
                    steps])
        elif remainder != 0:
            steps = remainder
            if exog is None:
                pred = forecaster.predict(steps=steps, last_window=last_window)
            else:
                pred = forecaster.predict(steps=steps, last_window=
                    last_window, exog=exog[last_window_end:last_window_end +
                    steps])
        else:
            continue
        backtest_predictions.append(pred)
    backtest_predictions = np.concatenate(backtest_predictions)
    metric_value = metric(y_true=y[initial_train_size:initial_train_size +
        len(backtest_predictions)], y_pred=backtest_predictions)
    return np.array([metric_value]), backtest_predictions


def grid_search_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
    param_grid: dict, initial_train_size: int, steps: int, metric: str,
    exog: Union[np.ndarray, pd.Series]=None, lags_grid: list=None, method:
    str='cv', allow_incomplete_fold: bool=True, return_best: bool=True,
    verbose: bool=True) ->pd.DataFrame:
    """
    Exhaustive search over specified parameter values for a Forecaster object.
    Validation is done using time series cross-validation or backtesting.
    
    Parameters
    ----------
    
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        ForecasterAutoreg, ForecasterCustom or ForecasterAutoregMultiOutput object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
        
    param_grid : dict
        Dictionary with parameters names (`str`) as keys and lists of parameter
        settings to try as values.
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
           
    lags_grid : list of int, lists, np.narray or range. 
        Lists of `lags` to try. Only used if forecaster is an instance of 
        `ForecasterCustom`.
        
    method : {'cv', 'backtesting'}
        Method used to estimate the metric for each parameter combination.
        'cv' for time series crosvalidation and 'backtesting' for simple
        backtesting. 'backtesting' is much faster since the model is fitted only
        once.
        
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    return_best : bool
        Refit the `forecaster` using the best found parameters on the whole data.
        
    verbose : bool, default `True`
        Print number of folds used for cv or backtesting.

    Returns 
    -------
    results: pandas.DataFrame
        Metric value estimated for each combination of parameters.

    """
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    if isinstance(forecaster, ForecasterCustom):
        if lags_grid is not None:
            logging.warning(
                '`lags_grid` ignored if forecaster is an instance of `ForecasterCustom`.'
                )
        lags_grid = ['custom predictors']
    elif lags_grid is None:
        lags_grid = [forecaster.lags]
    if isinstance(forecaster, ForecasterAutoregMultiOutput):
        param_grid_renamed = {('estimator__' + key): value for key, value in
            param_grid.items()}
        param_grid = param_grid_renamed
    lags_list = []
    params_list = []
    metric_list = []
    param_grid = list(ParameterGrid(param_grid))
    logging.info(f'Number of models to fit: {len(param_grid) * len(lags_grid)}'
        )
    for lags in tqdm.tqdm(lags_grid, desc='loop lags_grid', position=0):
        if isinstance(forecaster, (ForecasterAutoreg,
            ForecasterAutoregMultiOutput)):
            forecaster.set_lags(lags)
            lags = forecaster.lags.copy()
        for params in tqdm.tqdm(param_grid, desc='loop param_grid',
            position=1, leave=False):
            forecaster.set_params(**params)
            if method == 'cv':
                metrics = cv_forecaster(forecaster=forecaster, y=y, exog=
                    exog, initial_train_size=initial_train_size, steps=
                    steps, metric=metric, allow_incomplete_fold=
                    allow_incomplete_fold, verbose=verbose)
            else:
                metrics = backtesting_forecaster(forecaster=forecaster, y=y,
                    exog=exog, initial_train_size=initial_train_size, steps
                    =steps, metric=metric, verbose=verbose)[0]
            lags_list.append(lags)
            params_list.append(params)
            metric_list.append(metrics.mean())
    results = pd.DataFrame({'lags': lags_list, 'params': params_list,
        'metric': metric_list})
    results = results.sort_values(by='metric', ascending=True)
    if return_best:
        best_lags = results['lags'].iloc[0]
        best_params = results['params'].iloc[0]
        logging.info(
            f"""Refitting `forecaster` using the best found parameters and the whole data set: 
lags: {best_lags} 
params: {best_params}
"""
            )
        if isinstance(forecaster, (ForecasterAutoreg,
            ForecasterAutoregMultiOutput)):
            forecaster.set_lags(best_lags)
        forecaster.set_params(**best_params)
        forecaster.fit(y=y, exog=exog)
    return results


def backtesting_forecaster_intervals(forecaster, y: Union[np.ndarray, pd.
    Series], initial_train_size: int, steps: int, metric: str, exog: Union[
    np.ndarray, pd.Series]=None, interval: list=[5, 95], n_boot: int=500,
    in_sample_residuals: bool=True, verbose: bool=False) ->Tuple[np.array,
    np.array]:
    """
    Backtesting (validation) of `ForecasterAutoreg` or `ForecasterCustom` object.
    The model is trained only once using the `initial_train_size` first observations.
    In each iteration, a number of `steps` predictions are evaluated. Both, 
    predictions and intervals, are calculated.
    
    This evaluation is much faster than `cv_forecaster()` since the model is
    trained only once.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom
        `ForecasterAutoreg` or `ForecasterCustom` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
        
    interval: list, default `[5, 100]`
            Confidence of the prediction interval estimated. Sequence of percentiles
            to compute, which must be between 0 and 100 inclusive.
            
    n_boot: int, default `500`
        Number of bootstrapping iterations used to estimate prediction
        intervals.

    in_sample_residuals: bool, default `True`
        If `True`, residuals from the training data are used as proxy of
        prediction error to create prediction intervals.
        
    verbose : bool, default `True`
        Print number of folds used for backtesting.

    Returns 
    -------
    backtest_predictions: np.ndarray
        If `include_intervals=True`, 2D np.ndarray shape(steps, 3) with predicted
        value and their estimated interval.
            Column 0 = predictions
            Column 1 = lower bound interval
            Column 2 = upper bound interval
        
    metric_value: np.ndarray shape (1,)
        Value of the metric.

    Notes
    -----
    More information about prediction intervals in forecasting:
    https://otexts.com/fpp2/prediction-intervals.html
    Forecasting: Principles and Practice (2nd ed) Rob J Hyndman and
    George Athanasopoulos.

    """
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    if metric not in ['mean_squared_error', 'mean_absolute_error',
        'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and 'mean_absolute_percentage_error'. Got {metric}."
            )
    backtest_predictions = []
    metrics = {'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error}
    metric = metrics[metric]
    if exog is None:
        forecaster.fit(y=y[:initial_train_size])
    else:
        forecaster.fit(y=y[:initial_train_size], exog=exog[:initial_train_size]
            )
    folds = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    window_size = len(forecaster.last_window)
    if verbose:
        print(f'Number of observations used for training: {initial_train_size}'
            )
        print(
            f'Number of observations used for testing: {len(y) - initial_train_size}'
            )
        print(f'    Number of folds: {folds - 1 * (remainder == 0)}')
        print(f'    Number of steps per fold: {steps}')
        if remainder != 0:
            print(f'    Last fold only includes {remainder} observations')
    for i in range(folds):
        last_window_end = initial_train_size + i * steps
        last_window_start = initial_train_size + i * steps - window_size
        last_window = y[last_window_start:last_window_end]
        if i < folds - 1:
            if exog is None:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
            else:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, exog=exog[last_window_end:last_window_end +
                    steps], interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
        elif remainder != 0:
            steps = remainder
            if exog is None:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
            else:
                pred = forecaster.predict_interval(steps=steps, last_window
                    =last_window, exog=exog[last_window_end:last_window_end +
                    steps], interval=interval, n_boot=n_boot,
                    in_sample_residuals=in_sample_residuals)
        else:
            continue
        backtest_predictions.append(pred)
    backtest_predictions = np.concatenate(backtest_predictions)
    metric_value = metric(y_true=y[initial_train_size:], y_pred=
        backtest_predictions[:, 0])
    return np.array([metric_value]), backtest_predictions

################################################################################
#                        skforecast.model_selection                            #
#                                                                              #
# This work by Joaquín Amat Rodrigo is licensed under a Creative Commons       #
# Attribution 4.0 International License.                                       #
################################################################################
# coding=utf-8


import typing
from typing import Union, Dict, List, Tuple
import numpy as np
import pandas as pd
import logging
import tqdm
from sklearn.metrics import mean_squared_error 
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_absolute_percentage_error
from sklearn.model_selection import ParameterGrid

from .ForecasterAutoreg import ForecasterAutoreg
from .ForecasterCustom import ForecasterCustom
from .ForecasterAutoregMultiOutput import ForecasterAutoregMultiOutput

logging.basicConfig(
    format = '%(asctime)-5s %(name)-10s %(levelname)-5s %(message)s', 
    level  = logging.INFO,
)


def time_series_spliter(y: Union[np.ndarray, pd.Series],
                        initial_train_size: int, steps: int,
                        allow_incomplete_fold: bool=True,
                        verbose: bool=True):
    '''
    
    Split indices of a time series into multiple train-test pairs. The order of
    is maintained and the training set increases in each iteration.
    
    Parameters
    ----------        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    verbose : bool, default `True`
        Print number of splits created.

    Yields
    ------
    train : 1D np.ndarray
        Training indices.
        
    test : 1D np.ndarray
        Test indices.
        
    '''
    
    if not isinstance(y, (np.ndarray, pd.Series)):

        raise Exception('`y` must be `1D np.ndarray` o `pd.Series`.')

    elif isinstance(y, np.ndarray) and y.ndim != 1:

        raise Exception(
            f"`y` must be `1D np.ndarray` o `pd.Series`, "
            f"got `np.ndarray` with {y.ndim} dimensions."
        )
        
    if isinstance(y, pd.Series):
        y = y.to_numpy().copy()
    
  
    folds = (len(y) - initial_train_size) // steps  + 1
    # +1 fold is needed to allow including the remainder in the last iteration.
    remainder = (len(y) - initial_train_size) % steps   
    
    if verbose:
        if folds == 1:
            print(f"Number of folds: {folds - 1}")
            print("Not enought observations in `y` to create to create even a complete fold.")

        elif remainder == 0:
            print(f"Number of folds: {folds - 1}")

        elif remainder != 0 and allow_incomplete_fold:
            print(f"Number of folds: {folds}")
            print(
                f"Since `allow_incomplete_fold=True`, "
                f"last fold only includes {remainder} observations instead of {steps}."
            )
            print(
                'Incomplete folds with few observations could overestimate or ',
                'underestimate validation metrics.'
            )
        elif remainder != 0 and not allow_incomplete_fold:
            print(f"Number of folds: {folds - 1}")
            print(
                f"Since `allow_incomplete_fold=False`, "
                f"last {remainder} observations are descarted."
            )

    if folds == 1:
        # There are no observations to create even a complete fold
        return []
    
    for i in range(folds):
          
        if i < folds - 1:
            train_end     = initial_train_size + i * steps    
            train_indices = range(train_end)
            test_indices  = range(train_end, train_end + steps)
            
        else:
            if remainder != 0 and allow_incomplete_fold:
                train_end     = initial_train_size + i * steps  
                train_indices = range(train_end)
                test_indices  = range(train_end, len(y))
            else:
                break
        
        yield train_indices, test_indices
        
        

def cv_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
                  initial_train_size: int, steps: int,
                  metric: str, exog: Union[np.ndarray, pd.Series]=None,
                  allow_incomplete_fold: bool=True, verbose: bool=True) -> np.array:
    '''
    Cross-validation of `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput`
    object. The order of is maintained and the training set increases in each iteration.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
            
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    verbose : bool, default `True`
        Print number of folds used for cross validation.

    Returns 
    -------
    cv_results: 1D np.ndarray
        Value of the metric for each partition.

    '''
    
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)

    if metric not in ['mean_squared_error', 'mean_absolute_error',
                      'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and "
            f"'mean_absolute_percentage_error'. Got {metric}."
        )
    
    cv_results = []
    
    metrics = {
        'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error
    }
    
    metric = metrics[metric]
    
    splits = time_series_spliter(
                y                     = y,
                initial_train_size    = initial_train_size,
                steps                 = steps,
                allow_incomplete_fold = allow_incomplete_fold,
                verbose               = verbose
             )
    
    for train_index, test_index in splits:
        
        if exog is None:
            forecaster.fit(y=y[train_index])      
            pred = forecaster.predict(steps=len(test_index))
            
        else:
            
            forecaster.fit(y=y[train_index], exog=exog[train_index])      
            pred = forecaster.predict(steps=len(test_index), exog=exog[test_index])
               
        metric_value = metric(
                            y_true = y[test_index],
                            y_pred = pred
                       )
        
        cv_results.append(metric_value)
                          
    return np.array(cv_results)


def backtesting_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
                           initial_train_size: int, steps: int,
                           metric: str, exog: Union[np.ndarray, pd.Series]=None,
                           verbose: bool=False) -> Tuple[np.array, np.array]:
    '''
    Backtesting (validation) of `ForecasterAutoreg`, `ForecasterCustom` or
    `ForecasterAutoregMultiOutput` object.
    The model is trained only once using the `initial_train_size` first observations.
    In each iteration, a number of `steps` predictions are evaluated.
    
    This evaluation is much faster than `cv_forecaster()` since the model is
    trained only once.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        `ForecasterAutoreg`, `ForecasterCustom` or `ForecasterAutoregMultiOutput` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
            
    verbose : bool, default `False`
        Print number of folds used for backtesting.

    Returns 
    -------
    backtest_predictions: 1D np.ndarray
        Value of predictions.
        
    metric_value: np.ndarray shape (1,)
        Value of the metric.

    '''
    
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)

    if metric not in ['mean_squared_error', 'mean_absolute_error',
                      'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and "
            f"'mean_absolute_percentage_error'. Got {metric}."
        )
    
    backtest_predictions = []
    
    metrics = {
        'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error
    }
    
    metric = metrics[metric]
        
    if exog is None:
        forecaster.fit(y=y[:initial_train_size])      
    else:
        forecaster.fit(y=y[:initial_train_size], exog=exog[:initial_train_size])     
    
    
    folds     = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    window_size = len(forecaster.last_window)

    if isinstance(forecaster, ForecasterAutoregMultiOutput) and remainder != 0:
        # In ForecasterAutoregMultiOutput predictions are not iterative,
        # therefore no remainder is allowed.
        logging.warning(
                f"Backtesting `ForecasterAutoregMultiOutput` only allow completed "
                f"folds. Last {remainder} observations are excluded."
        )
        remainder=0
    
    if verbose:
        print(f"Number of observations used for training: {initial_train_size}")
        print(f"Number of observations used for testing: {len(y) - initial_train_size}")
        print(f"    Number of folds: {folds - 1 * (remainder == 0)}")
        print(f"    Number of steps per fold: {steps}")
        if remainder != 0:
            print(f"    Last fold only includes {remainder} observations")
      
    for i in range(folds):
        last_window_end   = initial_train_size + i * steps
        last_window_start = (initial_train_size + i * steps) - window_size 
        last_window       = y[last_window_start:last_window_end]
                
        if i < folds - 1:
            if exog is None:
                pred = forecaster.predict(
                            steps       = steps,
                            last_window = last_window
                        )
            else:
                pred = forecaster.predict(
                            steps       = steps,
                            last_window = last_window,
                            exog        = exog[last_window_end:last_window_end + steps]
                        )
        elif remainder != 0:
            steps = remainder 
            if exog is None:
                pred = forecaster.predict(
                            steps       = steps,
                            last_window = last_window
                        )
            else:
                pred = forecaster.predict(
                            steps       = steps,
                            last_window = last_window,
                            exog        = exog[last_window_end:last_window_end + steps]
                        )
        else:
            continue
        
        backtest_predictions.append(pred)
    
    backtest_predictions = np.concatenate(backtest_predictions)
    metric_value = metric(
                        y_true = y[initial_train_size: initial_train_size + len(backtest_predictions)],
                        y_pred = backtest_predictions
                   )

    return np.array([metric_value]), backtest_predictions



def grid_search_forecaster(forecaster, y: Union[np.ndarray, pd.Series],
                           param_grid: dict, initial_train_size: int, steps: int,
                           metric: str, exog: Union[np.ndarray, pd.Series]=None,
                           lags_grid: list=None, method: str='cv',
                           allow_incomplete_fold: bool=True, return_best: bool=True,
                           verbose: bool=True) -> pd.DataFrame:
    '''
    Exhaustive search over specified parameter values for a Forecaster object.
    Validation is done using time series cross-validation or backtesting.
    
    Parameters
    ----------
    
    forecaster : ForecasterAutoreg, ForecasterCustom, ForecasterAutoregMultiOutput
        ForecasterAutoreg, ForecasterCustom or ForecasterAutoregMultiOutput object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
        
    param_grid : dict
        Dictionary with parameters names (`str`) as keys and lists of parameter
        settings to try as values.
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
           
    lags_grid : list of int, lists, np.narray or range. 
        Lists of `lags` to try. Only used if forecaster is an instance of 
        `ForecasterCustom`.
        
    method : {'cv', 'backtesting'}
        Method used to estimate the metric for each parameter combination.
        'cv' for time series crosvalidation and 'backtesting' for simple
        backtesting. 'backtesting' is much faster since the model is fitted only
        once.
        
    allow_incomplete_fold : bool, default `True`
        The last test set is allowed to be incomplete if it does not reach `steps`
        observations. Otherwise, the latest observations are discarded.
        
    return_best : bool
        Refit the `forecaster` using the best found parameters on the whole data.
        
    verbose : bool, default `True`
        Print number of folds used for cv or backtesting.

    Returns 
    -------
    results: pandas.DataFrame
        Metric value estimated for each combination of parameters.

    '''
    
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)
    
    if isinstance(forecaster, ForecasterCustom):
        if lags_grid is not None:
            logging.warning(
                '`lags_grid` ignored if forecaster is an instance of `ForecasterCustom`.'
            )
        lags_grid = ['custom predictors']
        
    elif lags_grid is None:
        lags_grid = [forecaster.lags]
        
    if isinstance(forecaster, ForecasterAutoregMultiOutput):
        param_grid_renamed = {'estimator__'+key: value for key, value in param_grid.items()}
        param_grid = param_grid_renamed
        
    lags_list = []
    params_list = []
    metric_list = []
    
    param_grid =  list(ParameterGrid(param_grid))

    logging.info(
        f"Number of models to fit: {len(param_grid)*len(lags_grid)}"
    )
    
    for lags in tqdm.tqdm(lags_grid, desc='loop lags_grid', position=0):
        
        if isinstance(forecaster, (ForecasterAutoreg, ForecasterAutoregMultiOutput)):
            forecaster.set_lags(lags)
            lags = forecaster.lags.copy()
        
        for params in tqdm.tqdm(param_grid, desc='loop param_grid', position=1, leave=False):

            forecaster.set_params(**params)
            
            if method == 'cv':
                metrics = cv_forecaster(
                                forecaster = forecaster,
                                y          = y,
                                exog       = exog,
                                initial_train_size = initial_train_size,
                                steps  = steps,
                                metric = metric,
                                allow_incomplete_fold = allow_incomplete_fold,
                                verbose = verbose
                             )
            else:
                metrics = backtesting_forecaster(
                                forecaster = forecaster,
                                y          = y,
                                exog       = exog,
                                initial_train_size = initial_train_size,
                                steps  = steps,
                                metric = metric,
                                verbose = verbose
                             )[0]

            lags_list.append(lags)
            params_list.append(params)
            metric_list.append(metrics.mean())
            
    results = pd.DataFrame({
                'lags'  : lags_list,
                'params': params_list,
                'metric': metric_list})
    
    results = results.sort_values(by='metric', ascending=True)
    
    if return_best:
        
        best_lags = results['lags'].iloc[0]
        best_params = results['params'].iloc[0]
        logging.info(
            f"Refitting `forecaster` using the best found parameters and the whole data set: \n"
            f"lags: {best_lags} \n"
            f"params: {best_params}\n"
        )
        
        if isinstance(forecaster, (ForecasterAutoreg, ForecasterAutoregMultiOutput)):
            forecaster.set_lags(best_lags)
            
        forecaster.set_params(**best_params)
        forecaster.fit(y=y, exog=exog)
            
    return results


def backtesting_forecaster_intervals(forecaster, y: Union[np.ndarray, pd.Series],
                           initial_train_size: int, steps: int,
                           metric: str, exog: Union[np.ndarray, pd.Series]=None,
                           interval: list=[5, 95], n_boot: int=500,
                           in_sample_residuals: bool=True,
                           verbose: bool=False) -> Tuple[np.array, np.array]:
    '''
    Backtesting (validation) of `ForecasterAutoreg` or `ForecasterCustom` object.
    The model is trained only once using the `initial_train_size` first observations.
    In each iteration, a number of `steps` predictions are evaluated. Both, 
    predictions and intervals, are calculated.
    
    This evaluation is much faster than `cv_forecaster()` since the model is
    trained only once.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg, ForecasterCustom
        `ForecasterAutoreg` or `ForecasterCustom` object.
        
    y : 1D np.ndarray, pd.Series
        Training time series values. 
    
    initial_train_size: int 
        Number of samples in the initial train split.
        
    steps : int
        Number of steps to predict.
        
    metric : {'mean_squared_error', 'mean_absolute_error', 'mean_absolute_percentage_error'}
        Metric used to quantify the goodness of fit of the model.
        
    exog : np.ndarray, pd.Series, default `None`
            Exogenous variable/s included as predictor/s. Must have the same
            number of observations as `y` and should be aligned so that y[i] is
            regressed on exog[i].
        
    interval: list, default `[5, 100]`
            Confidence of the prediction interval estimated. Sequence of percentiles
            to compute, which must be between 0 and 100 inclusive.
            
    n_boot: int, default `500`
        Number of bootstrapping iterations used to estimate prediction
        intervals.

    in_sample_residuals: bool, default `True`
        If `True`, residuals from the training data are used as proxy of
        prediction error to create prediction intervals.
        
    verbose : bool, default `True`
        Print number of folds used for backtesting.

    Returns 
    -------
    backtest_predictions: np.ndarray
        If `include_intervals=True`, 2D np.ndarray shape(steps, 3) with predicted
        value and their estimated interval.
            Column 0 = predictions
            Column 1 = lower bound interval
            Column 2 = upper bound interval
        
    metric_value: np.ndarray shape (1,)
        Value of the metric.

    Notes
    -----
    More information about prediction intervals in forecasting:
    https://otexts.com/fpp2/prediction-intervals.html
    Forecasting: Principles and Practice (2nd ed) Rob J Hyndman and
    George Athanasopoulos.

    '''
    
    forecaster._check_y(y=y)
    y = forecaster._preproces_y(y=y)
    
    if exog is not None:
        forecaster._check_exog(exog=exog)
        exog = forecaster._preproces_exog(exog=exog)

    if metric not in ['mean_squared_error', 'mean_absolute_error',
                      'mean_absolute_percentage_error']:
        raise Exception(
            f"Allowed metrics are: 'mean_squared_error', 'mean_absolute_error' and "
            f"'mean_absolute_percentage_error'. Got {metric}."
        )
    
    backtest_predictions = []
    
    metrics = {
        'mean_squared_error': mean_squared_error,
        'mean_absolute_error': mean_absolute_error,
        'mean_absolute_percentage_error': mean_absolute_percentage_error
    }
    
    metric = metrics[metric]
        
    if exog is None:
        forecaster.fit(y=y[:initial_train_size])   
    else:
        forecaster.fit(y=y[:initial_train_size], exog=exog[:initial_train_size])     
    
    
    folds     = (len(y) - initial_train_size) // steps + 1
    remainder = (len(y) - initial_train_size) % steps
    window_size = len(forecaster.last_window)
    
    if verbose:
        print(f"Number of observations used for training: {initial_train_size}")
        print(f"Number of observations used for testing: {len(y) - initial_train_size}")
        print(f"    Number of folds: {folds - 1 * (remainder == 0)}")
        print(f"    Number of steps per fold: {steps}")
        if remainder != 0:
            print(f"    Last fold only includes {remainder} observations")
    
    for i in range(folds):
        
        last_window_end   = initial_train_size + i * steps
        last_window_start = (initial_train_size + i * steps) - window_size 
        last_window       = y[last_window_start:last_window_end]
                
        if i < folds - 1:
            if exog is None:
                pred = forecaster.predict_interval(
                            steps       = steps,
                            last_window = last_window,
                            interval    = interval,
                            n_boot      = n_boot,
                            in_sample_residuals = in_sample_residuals
                        )
            else:
                pred = forecaster.predict_interval(
                            steps       = steps,
                            last_window = last_window,
                            exog        = exog[last_window_end:last_window_end + steps],
                            interval    = interval,
                            n_boot      = n_boot,
                            in_sample_residuals = in_sample_residuals
                        )
        elif remainder != 0:
            steps = remainder 
            if exog is None:
                pred = forecaster.predict_interval(
                            steps       = steps,
                            last_window = last_window,
                            interval    = interval,
                            n_boot      = n_boot,
                            in_sample_residuals = in_sample_residuals
                        )
            else:
                pred = forecaster.predict_interval(
                            steps       = steps,
                            last_window = last_window,
                            exog        = exog[last_window_end:last_window_end + steps],
                            interval    = interval,
                            n_boot      = n_boot,
                            in_sample_residuals = in_sample_residuals
                        )
        else:
            continue
        
        backtest_predictions.append(pred)
    
    backtest_predictions = np.concatenate(backtest_predictions)
    metric_value = metric(
                        y_true = y[initial_train_size:],
                        y_pred = backtest_predictions[:, 0]
                   )

    return np.array([metric_value]), backtest_predictions