def _make_multiple_historical_forecasts(
        self,
        train_n_points: int,
        series: Union[TimeSeries, Sequence[TimeSeries]],
        direct_predictions: Union[TimeSeries, Sequence[TimeSeries]],
        past_covariates: Optional[Union[TimeSeries, Sequence[TimeSeries]]] = None,
        future_covariates: Optional[Union[TimeSeries, Sequence[TimeSeries]]] = None,
        num_samples: int = 1,
    ) -> Union[TimeSeries, Sequence[TimeSeries]]:
        """
        For GlobalForecastingModel, when predicting n > output_chunk_length, `historical_forecasts()` generally
        produce better forecasts than `predict()`.

        To get as close as possible to the predictions generated by the forecasting models during inference,
        `historical_forecasts` forecast horizon is equal to each model output_chunk_length.

        train_n_points are generated, starting from the end of the series.
        """
        is_single_series = isinstance(series, TimeSeries)
        series = series2seq(series)
        direct_predictions = series2seq(direct_predictions)
        past_covariates = series2seq(past_covariates)
        future_covariates = series2seq(future_covariates)

        n_components = series[0].n_components
        model_predict_cols = direct_predictions[0].columns.tolist()

        predictions = []
        for m_idx, model in enumerate(self.forecasting_models):
            # we start historical fc at multiple of the output length before the end.
            n_ocl_back = train_n_points // model.output_chunk_length
            start_hist_forecasts = n_ocl_back * model.output_chunk_length

            # we use the precomputed `direct_prediction` to fill any missing prediction
            # timesteps at the beginning (if train_n_points is not perfectly divisible by output length)
            missing_steps = train_n_points % model.output_chunk_length

            tmp_pred = model.historical_forecasts(
                series=series,
                past_covariates=(
                    past_covariates if model.supports_past_covariates else None
                ),
                future_covariates=(
                    future_covariates if model.supports_future_covariates else None
                ),
                forecast_horizon=model.output_chunk_length,
                stride=model.output_chunk_length,
                num_samples=(
                    num_samples if model.supports_probabilistic_prediction else 1
                ),
                start=-start_hist_forecasts,
                start_format="position",
                retrain=False,
                overlap_end=False,
                last_points_only=False,
                show_warnings=self.show_warnings,
                predict_likelihood_parameters=False,
            )
            # concatenate the strided predictions of output_chunk_length values each
            tmp_pred = [concatenate(sub_pred, axis=0) for sub_pred in tmp_pred]

            # add the missing steps at beginning by taking the first values of precomputed predictions
            if missing_steps:
                # add the missing steps at beginning by taking the first values of precomputed predictions
                # get the model's direct (uni/multivariate) predictions
                pred_cols = model_predict_cols[
                    m_idx * n_components : (m_idx + 1) * n_components
                ]
                hfc_cols = tmp_pred[0].columns.tolist()
                tmp_pred = [
                    concatenate(
                        [
                            preds_dir[:missing_steps][pred_cols].with_columns_renamed(
                                pred_cols, hfc_cols
                            ),
                            preds_hfc,
                        ],
                        axis=0,
                    )
                    for preds_dir, preds_hfc in zip(direct_predictions, tmp_pred)
                ]
            predictions.append(tmp_pred)

        tmp_predictions = []
        # slice the forecasts, training series-wise, to align them
        for prediction in predictions:
            tmp_predictions.append([ts for idx, ts in enumerate(prediction)])
        predictions = [seq2series(prediction) for prediction in tmp_predictions]

        # reduce the probabilistics series
        if self.train_samples_reduction is not None and self.train_num_samples > 1:
            predictions = [
                self._predictions_reduction(prediction) for prediction in predictions
            ]

        return (
            self._stack_ts_seq(predictions)
            if is_single_series
            else self._stack_ts_multiseq(predictions)
        )

        For GlobalForecastingModel, when predicting n > output_chunk_length, `historical_forecasts()` generally
        produce better forecasts than `predict()`.

        To get as close as possible to the predictions generated by the forecasting models during inference,
        `historical_forecasts` forecast horizon is equal to each model output_chunk_length.

        train_n_points are generated, starting from the end of the series.
        
darts/models/forecasting/regression_ensemble_model.py