text: import warnings
  from math import sqrt
  from collections import OrderedDict
  from statistics import mean
  
  import numpy as np
  import statsmodels.api as sm
  from sklearn.metrics import mean_squared_error
  
  
  class ArimaModel:
      """Custom ArimaModel class"""
  
      def __init__(self, series_train_data, validation_criteria, data_frequency=0):
          """
  
          Parameters
          ----------
          series_train_data: list of values of time series
          validation_criteria: 'aic' for Akaike Information Criteria, 'cross' for cross validation with rmse
          data_frequency: interval of data (12 for monthly series, 4 for quaternal); if value bigger than 1 is provided,
                          then a seasonal model will be used; values lower than 2 are ignored
  
          """
          self.series_train_data = series_train_data
          self.validation_criteria = validation_criteria
          self.data_frequency = data_frequency
  
          self.series_size = len(self.series_train_data)
  
          # maximal values of standard arima model parameters
          self.max_p = 5
          self.max_d = 5
          self.max_q = 5
  
          # maximal values of seasonal arima model parameters
          self.max_P = 2
          self.max_D = 1
          self.max_Q = 2
  
          # significance level for KPSS test
          self.significance_level = 0.05
  
          if self.validation_criteria == 'aic':
              self.validation_function = self.__aic
          else:
              self.validation_function = self.__time_series_cross_validation
              self.max_nr_of_folds = 6
  
          # do not use seasonal component if series is of yearly frequency
          if data_frequency < 2:
              self.data_frequency = 0
  
          # parameters of arima model
          self.p, self.d, self.q, self.P, self.D, self.Q = self.__find_parameters()
  
          # fit a model with found parameters
          self.model = self.__fit_model()
  
      def __find_parameters(self, train_data=None):
          """Determine values of arima model parameters according to Hyndman-Khandakar procedure"""
          
          if train_data is None:
              train_data = self.series_train_data
  
          train_data_size = len(train_data)
  
          # seasonal models
          if self.data_frequency > 0:
              possible_P_values = range(self.max_P + 1)
              possible_D_values = range(self.max_D + 1)
              possible_Q_values = range(self.max_Q + 1)
              # initial sets of p, q, P, Q values in Hyndman-Khandakar procedure for seasonal models
              start_parameters_sets = [(0, 0, 0, 0), (1, 0, 1, 0), (0, 1, 0, 1), (2, 2, 1, 1)]
          # nonseasonal models
          else:
              possible_P_values = [0]
              possible_D_values = [0]
              possible_Q_values = [0]
              # initial sets of p, q, P, Q values in Hyndman-Khandakar procedure for nonseasonal models
              start_parameters_sets = [(0, 0, 0, 0), (1, 0, 0, 0), (0, 1, 0, 0), (2, 2, 0, 0)]
  
          possible_p_values = range(self.max_p + 1)
          possible_q_values = range(self.max_q + 1)
  
          best_parameters_list = []
  
          # since there is no extended Canova-Hansen test implemented in Python we will simply follow Hyndman-Khandakar
          # procedure for every possible D value and choose the best model according to chosen validation_criteria
          for D in possible_D_values:
  
              # set d parameter according to kpss test
              d = self.__find_d_parameter(D)
  
              validation_values_dict = {}
              # get only these initial sets of parameters for which there is sufficient amount of train data to fit model
              parameters_sets = [(p, q, P, Q) for p, q, P, Q in start_parameters_sets
                                 if self.__sufficient_train_data_size(train_data_size, p, d, q, P, D, Q)]
              for p, q, P, Q in parameters_sets:
                  if (p, q, P, Q) not in validation_values_dict:
                      validation_values_dict[(p, q, P, Q)] = \
                          self.validation_function(train_data, (p, d, q), (P, D, Q, self.data_frequency))
  
              best_p, best_q, best_P, best_Q = min(validation_values_dict, key=validation_values_dict.get)
              best_value = validation_values_dict[(best_p, best_q, best_P, best_Q)]
              parameters_sets = self.__next_sets_of_parameters(best_p, best_q, best_P, best_Q)
  
              # in case of nonseasonal models sets with P and Q other than 0 are dropped
              parameters_sets = [(p, q, P, Q) for p, q, P, Q in parameters_sets
                                 if p in possible_p_values and q in possible_q_values
                                 and P in possible_P_values and Q in possible_Q_values
                                 and self.__sufficient_train_data_size(train_data_size, p, d, q, P, D, Q)]
  
              # remove sets of parameters for which we have insufficient amount of train data to fit the model reasonably
  
              while parameters_sets:
                  for p, q, P, Q in parameters_sets:
                      if (p, q, P, Q) not in validation_values_dict:
                          try:
                              validation_values_dict[(p, q, P, Q)] = \
                                  self.validation_function(train_data, (p, d, q), (P, D, Q, self.data_frequency))
                          except Exception as exception:
                              print('Evaluating model with parameters: {} {} failed due to exception: {}. Skipping this '
                                    'set of parameters.'.format((p, d, q), (P, D, Q, self.data_frequency), exception))
                              # if evaluating a model with given parameters was unsuccessful then add them to dict
                              # with best_value so far -> they will not be used as smaller value is searched
                              validation_values_dict[(p, q, P, Q)] = best_value
  
                  min_value = min(validation_values_dict.values())
  
                  if min_value < best_value:
                      best_p, best_q, best_P, best_Q = min(validation_values_dict, key=validation_values_dict.get)
  
                      best_value = min_value
                      parameters_sets = self.__next_sets_of_parameters(best_p, best_q, best_P, best_Q)
  
                      parameters_sets = [(p, q, P, Q) for p, q, P, Q in parameters_sets
                                         if p in possible_p_values and q in possible_q_values
                                         and P in possible_P_values and Q in possible_Q_values
                                         and self.__sufficient_train_data_size(train_data_size, p, d, q, P, D, Q)]
  
                  else:
                      break
  
              best_parameters_list.append((best_value, best_p, d, best_q, best_P, D, best_Q))
  
          _, p, d, q, P, D, Q = min(best_parameters_list, key=lambda x: x[0])
  
          return p, d, q, P, D, Q
  
      def __find_d_parameter(self, D):
          """Find number of differences d according to KPSS statistical test"""
  
          # set the starting number of differences
          d = 0
          differentiated_data = self.series_train_data.copy()
  
          # if model has a seasonal component then apply seasonal differences first
          for _ in range(D):
              differentiated_data = [differentiated_data[i] - differentiated_data[i - self.data_frequency]
                                     for i in range(self.data_frequency, len(differentiated_data))]
  
          while d < self.max_d:
              # suppress warning about p-value outside the table of critical values
              # they mean that either p-value is greater than 0.1 or smaller than 0.01 which is irrelevant here
              with warnings.catch_warnings():
                  warnings.simplefilter("ignore")
                  _, kpss_p_value, _, _ = sm.tsa.stattools.kpss(differentiated_data)
              if kpss_p_value > self.significance_level:
                  break
              differentiated_data = [differentiated_data[i] - differentiated_data[i - 1]
                                     for i in range(1, len(differentiated_data))]
              d += 1
  
          return d
  
      def __next_sets_of_parameters(self, p, q, P, Q):
          """
          Return new variations of arima parameters according to Hyndman-Khandakar procedure:
              where one of p, q, P and Q is allowed to vary by ±1 from the current model or
              where p and q both vary by ±1 from the current model or
              where P and Q both vary by ±1 from the current model
          """
  
          return [(p + 1, q, P, Q),
                  (p - 1, q, P, Q),
                  (p, q + 1, P, Q),
                  (p, q - 1, P, Q),
                  (p, q, P + 1, Q),
                  (p, q, P - 1, Q),
                  (p, q, P, Q + 1),
                  (p, q, P, Q - 1),
                  (p + 1, q + 1, P, Q),
                  (p - 1, q - 1, P, Q),
                  (p, q, P + 1, Q + 1),
                  (p, q, P - 1, Q - 1)]
  
      def __sufficient_train_data_size(self, train_data_size, p, d, q, P, D, Q):
          """Check if train data size is sufficient to fit a model with given set of parameters"""
          # minimal size of series train data needed to fit a model with given parameters
          minimal_data_size = d + D * self.data_frequency + \
                              max(3 * q + 1, 3 * Q * self.data_frequency + 1, p, P * self.data_frequency) + 1
  
          if self.validation_criteria == 'aic':
              return minimal_data_size <= train_data_size
          # if cross validation is used for model validation check if we have at least twice as much data
          else:
              return minimal_data_size <= int(train_data_size / 2)
  
      def __aic(self, train_data, parameters, seasonal_parameters):
          """Fit a model with its train data and provided arguments and return Akaike Information Criteria"""
          fitted_model = self.__fit_model(train_data, parameters, seasonal_parameters)
          return fitted_model.aic
  
      def __time_series_cross_validation(self, train_data, parameters, seasonal_parameters):
          """Perform cross validation of time series and return mean of rmse of all models created for validation"""
          fold_nr = 0
          rmse_per_fold = []
  
          train_data_size = len(train_data)
  
          # set window_size so that max_nr_of_folds validations can be done
          # and model is fitted on at least half of train data
          window_size = int(train_data_size / 2 / self.max_nr_of_folds) or 1
          window_slide = window_size
  
          while fold_nr < self.max_nr_of_folds:
              index = train_data_size - window_size - window_slide * fold_nr
              # at least as much fit data as window_size
              if index < window_size:
                  break
              fit_data = train_data[:index]
              validation_data = train_data[index:index + window_size]
              fitted_model = self.__fit_model(fit_data, parameters, seasonal_parameters)
              validation_forecast = fitted_model.forecast(window_size).tolist()
              rmse_per_fold.append(sqrt(mean_squared_error(validation_data, validation_forecast)))
              fold_nr += 1
          return mean(rmse_per_fold)
  
      def __fit_model(self, fit_data=None, parameters=None, seasonal_parameters=None):
          """Fit arima model with given data and parameters"""
          
          if fit_data is None:
              fit_data=self.series_train_data
              
          if parameters is None:
              parameters=(self.p, self.d, self.q)
              
          if seasonal_parameters is None:
              seasonal_parameters=(self.P, self.D, self.Q, self.data_frequency)
              
          # suppress FutureWarning about wrong input data format
          with warnings.catch_warnings():
              warnings.filterwarnings("ignore", category=FutureWarning)
  
          # we use SARIMAX to avoid errors with non stationary series when fitting standard ARIMA model
              arima_model = sm.tsa.statespace.SARIMAX(fit_data,
                                                      order=parameters,
                                                      seasonal_order=seasonal_parameters,
                                                      enforce_stationarity=False,
                                                      enforce_invertibility=False)
              fitted_model = arima_model.fit(maxiter=200, method='powell', disp=False)
  
          return fitted_model
  
      def make_prediction(self, series_test_data, forecast_type):
          """Make prediction on train data and forecast of test data using either full or one-step method"""
  
          # prediction of train data
          prediction = self.model.predict()
          test_data_size = len(series_test_data)
  
          if forecast_type == 'full':
              # forecast of all test data
              forecast = self.model.forecast(test_data_size)
          elif forecast_type == 'one_step':
              refit_train_data = self.series_train_data.copy()
              # forecast of first value from test data
              forecast = self.model.forecast(1)
              # add expected value to train data
              refit_train_data.append(series_test_data[0])
  
              for i in range(1, test_data_size):
                  # reset arima parameters according to new train data
                  p, d, q, P, D, Q = self.__find_parameters(refit_train_data)
                  # fit new model on increased train_data
                  refitted_model = self.__fit_model(refit_train_data, (p, d, q), (P, D, Q, self.data_frequency))
                  # forecast single next value
                  forecast = np.append(forecast, refitted_model.forecast(1))
                  # append real value for next time step forecast
                  refit_train_data.append(series_test_data[i])
          return prediction.tolist(), forecast.tolist()
  
      def get_parameters(self):
          """Return list of model parameters"""
          return OrderedDict([('p', self.p), ('d', self.d), ('q', self.q), ('P', self.P), ('D', self.D), ('Q', self.Q)])
  