Forecasting Inflation: Romanian Case Study Using SARIMA Models

    	
            
import DatastreamDSWS as dsws
 
# We can use our Refinitiv's Datastream Web Socket (DSWS) API keys that allows us to be identified by Refinitiv's back-end services and enables us to request (and fetch) data: Credentials are placed in a text file so that it may be used in this code without showing it itself.
(dsws_username, dsws_password) = (open("Datastream_username.txt","r"),
                                  open("Datastream_password.txt","r"))
 
ds = dsws.Datastream(username = str(dsws_username.read()),
                     password = str(dsws_password.read()))
 
# It is best to close the files we opened in order to make sure that we don't stop any other services/programs from accessing them if they need to.
dsws_username.close()
dsws_password.close()
 
 
# # Alternatively one can use the following:
# import getpass
# dsusername = input()
# dspassword = getpass.getpass()
# ds = dsws.Datastream(username = dsusername, password = dspassword)

        
        Copy
    

    	
            
import eikon as ek
 
# The key is placed in a text file so that it may be used in this code without showing it itself:
eikon_key = open("eikon.txt","r")
ek.set_app_key(str(eikon_key.read()))
# It is best to close the files we opened in order to make sure that we don't stop any other services/programs from accessing them if they need to:
eikon_key.close()

        
        Copy
    

    	
            
import numpy as np
import pandas as pd
import pickle # need to ' pip install pickle-mixin '. This library is native to Python and therefore doesn't have a version of its own.
from datetime import datetime, timedelta # ' datetime ' is native to Python and therefore doesn't have a version of its own.

        
        Copy
    

    	
            
# Import the relevant plotting libraries:
 
import plotly # Needed to show plots in line in Notebook.
from plotly.offline import init_notebook_mode, iplot
init_notebook_mode(connected = False) # Initialize plotly.js in the browser if it hasn't been loaded into the Document Object Model (DOM) yet.
 
import cufflinks, seaborn # Needed to show plots in line in Notebook.
cufflinks.go_offline() # plotly & cufflinks are online platforms but we want to stay offline.

        
        Copy
    

    	
            
for i,j in zip([ek, np, pd, plotly, cufflinks, seaborn],
               ["eikon", "numpy", "pandas", "plotly", "cufflinks", "seaborn"]):
    print("The library '" + j + "' in this notebook is version " + i.__version__)

        
        Copy
    

    	
            
df = ds.get_data(tickers = "RMCPANNL,RMUNPTOTP,RMXRUSD.,RMXREUR.,RMGOVBALA", # be carefulnot to put spaces in between elements here, or else these spaces will be included in column names.
                 fields = "X",
                 start = '2000-01-01',
                 end = '2021-04-06',
                 freq = 'M')

        
        Copy
    

    	
            df["RMCPANNL", "X"] = df["RMCPANNL"]["X"]/100
        
        Copy
    

    	
            df
        
        Copy
    

    	
            
df.columns = pd.MultiIndex.from_tuples(
    [('RMCPANNL', 'yoy'), # yoy will be the increase in one month compared to the same month of the previous year
     ('RMUNPTOTP', 'raw'),
     ('RMXRUSD.', 'raw'),
     ('RMXREUR.', 'raw'),
     ('RMGOVBALA', 'raw')],
    names=['Instrument', 'Field'])

        
        Copy
    

    	
            
CLc1, err = ek.get_data(
    instruments = "CLc1",
    fields = ["TR.CLOSEPRICE.timestamp",
              "TR.CLOSEPRICE",
              "MAVG(TR.CLOSEPRICE,-29)"],
    parameters = {'SDate': '1999-01-01', # We want to start 1 year earlier because we're about to take moving averages which makes us loose degrees of freedom.
                  'EDate': str(datetime.today())[:10], # ' df.index[-1] ' picks the last date in our previously defined data-frame ' df '. You may want to use 'str(datetime.today())[:10]'.
                  'FRQ': 'D'})
CLc1

        
        Copy
    

    	
            
# This takes the moving average of the last 30 datapoints, not the last 30 days (since not all of the last 30 days might have been trading days with price data).
CLc1["CLc1 30D Moving Average"] = CLc1["Close Price"].dropna().rolling(30).mean()
# Let's remove the ' MAVG(TR.CLOSEPRICE,-29) ' column, it was only there to display the issue outlined.
CLc1.drop('MAVG(TR.CLOSEPRICE,-29)', axis = 1, inplace = True)
CLc1 # Let's have a look at our data-frame now.

        
        Copy
    

    	
            CLc1.Timestamp[0] # Note that ' CLc1.Timestamp ' is the same as CLc1["Timestamp"]
        
        Copy
    

    	
            df.index[0]
        
        Copy
    

    	
            '2000-01-31T00:00:00Z'[0:10]
        
        Copy
    

    	
            
CLc1.index = [i[0:10] for i in CLc1.Timestamp]
CLc1

        
        Copy
    

    	
            
def Match_previous_days(df1, df2, common_day = 15, back = 3):
    """Match_previous_days Version 1.0:
    This function compares two Pandas data-frames (df1 and df2, with monthly and daily
    data respectively and dates in its index), and checks if any index from df1 is
    missing in df2. If one for a day set in 'common_day' is missing, it changes it to
    the previous day. It does so as many times as is set in 'back'.
    It is created to ease the joining of Datastream API (DSWS) and Eikon API
    retrieved Pandas data-frames (df1 and df2 respectively). E.g.:
    Setting ' common_day ' to ' 15 ' and back to ' 3 ', our function will scan for
    any row-name (called an index) in df1 for the 15th of the month that does not
    appear in df2 and change any such index name to the 14th of the month (essentially
    showing df2 data for the 14th of the month for those rows specifically as 15th of
    the month). It can do it again for the 13th of the month, and then 12th; 3 times 
    in total (as specified in ' back ').
    
    
    Parameters:
    ----------
     
    df1: Pandas data-frame
        A monthly data-frame from DSWS - or formatted similarly
        (i.e.: with Multi-Index'ed columns, etc.).
     
    df2: Pandas data-frame
        A daily data-frame from eikon's API - or formatted
        similarly - with one exception that its index is made of strings for dates in the
        "yyyy-mm-dd" format.
    
    common_day: int
        The day of the month to start from. For an example, see 1st
        paragraph above.
        Default: common_day = 15
    
    common_day: int
        The number of times we are ready to run the procedure for. For an
        example, see 1st paragraph above.
        Default: back = 3
    
    
    Dependencies:
    ----------
    
    pandas 1.2.3
    numpy 1.19.5
    
    
    Examples:
    --------
    
    >>> import DatastreamDSWS as dsws
    >>> ds = dsws.Datastream(username = "insert dsws username here", password = "insert dsws password here")
    >>> 
    >>> import eikon as ek
    >>> ek.set_app_key("insert eikon key here")
    >>> 
    >>> df_1 = ds.get_data(tickers = "RMCPANNL,RMUNPTOTP,RMXRUSD.,RMXREUR.,RMGOVBALA", # be carefulnot to put spaces in between elements here, or else these spaces will be included in column names.
    >>>                    fields = "X", start = '2000-01-01', freq = 'M')
    >>> 
    >>> df_2, err = ek.get_data(instruments = "CLc1", fields = ["TR.CLOSEPRICE.timestamp", "TR.CLOSEPRICE"],
    >>>                         parameters = {'SDate': '1999-01-01', # We want to start 1 year earlier because we're about to take moving averages which makes us loose degrees of freedom.
    >>>                                       'EDate': str(datetime.today())[:10], # ' df_1.index[-1] ' picks the last date in our previously defined data-frame ' df_1 '. You may want to use 'str(datetime.today())[:10]'.
    >>>                                       'FRQ': 'D'})
    >>> df_2["CLc1 30D Moving Average"] = df_2["Close Price"].rolling(30).mean() # This takes the moving average of the last 30 datapoints, not the last 30 days (since not all of the last 30 days might have been trading days with price data).
    >>> df_2.index = [i[0:10] for i in CLc1.Timestamp]
    >>> 
    >>> Match_previous_days(df1 = df_1, df2 = df_2, common_day = 15, back = 3) # Note that this is the same as ' Match_previous_days(df1 = df_1, df2 = df_2) '.
    """
    
    _common_day = "-" + "{:02d}".format(common_day)
    df2_list = []
    df2_df1_diff_index1 = []
    df2_df1_diff_index2 = []
    index_dictionary = []
    
    df2_list.append(df2.copy())
    
    for k in range(back):
        
        if k > 0:
            df2_list.append(df2_list[-1])
        
        df2_df1_diff_index1.append(np.setdiff1d(
            df1.index.tolist(),
            df2_list[k][df2_list[k].index.str.contains(_common_day)].index.tolist()))
        
        df2_df1_diff_index2.append([
            i.replace("-15",  "-" + "{:02d}".format(common_day - 1 - k)) # ' -1 ' because k starts at 0.
            for i in df2_df1_diff_index1[k]])
        
        index_dictionary.append({})
#         display(index_dictionary[k])
        for i,j in zip(df2_df1_diff_index2[k], df2_df1_diff_index1[k]):
            index_dictionary[k][i] = j
#         display(index_dictionary[k])
        
        df2_list[k].rename(index = index_dictionary[k], inplace = True)
        
    return df2_list[-1]

        
        Copy
    

    	
            CLc1_of_interest0 = Match_previous_days(df1 = df, df2 = CLc1, common_day = 15, back = 4) # Applying the previously defined function.
        
        Copy
    

    	
            CLc1_of_interest0.iloc[258:261]
        
        Copy
    

    	
            
CLc1_of_interest1 = CLc1_of_interest0[CLc1_of_interest0.index.str.contains("-15")] # We're only interested in the 15th of each month.
 
CLc1_of_interest = pd.DataFrame(
    data = CLc1_of_interest1[
        ["Close Price", "CLc1 30D Moving Average"] # We're only interested in those two columns.
    ].values,
    index = CLc1_of_interest1.index,
    columns = pd.MultiIndex.from_tuples(
        [('CLc1', 'Close Price'),
         ('CLc1', 'Close Price 30D Moving Average')]))
 
CLc1_of_interest

        
        Copy
    

    	
            
df = pd.merge(df, CLc1_of_interest,
              left_index = True,
              right_index = True)
df

        
        Copy
    

    	
            
df["RMUNPTOTP", "yoy"] = df["RMUNPTOTP"]["raw"].pct_change(periods = 12)
df["RMXRUSD.", "yoy"] = df["RMXRUSD."]["raw"].pct_change(periods = 12)
df["RMXRUSD.", "yoy"] = df["RMXRUSD."]["raw"].pct_change(periods = 12)
df["RMXREUR.", "yoy"] = df["RMXREUR."]["raw"].pct_change(periods = 12)
df["RMGOVBALA", "yoy"] = df["RMGOVBALA"]["raw"].pct_change(periods = 12)
df["CLc1", "yoy"] = df["CLc1"][
    "Close Price 30D Moving Average"].pct_change(periods = 12)
df = df.sort_index(axis=1)
df = df.dropna()
df

        
        Copy
    

    	
            
df["RMCPANNL"]["yoy"].iplot(theme = "solar", title = "RMCPANNL", xTitle = "date",
                            yTitle = "% increase from same month in previous year (yoy)")

        
        Copy
    

    	
            
df["RMCPANNL"]["yoy"].diff(1).iplot(theme = "solar",
                                    title = "RMCPANNL's 1st difference",
                                    xTitle = "date",
                                    yTitle = "1st difference in yoy (yoy1d)")

        
        Copy
    

    	
            
import arch
from arch.unitroot import ADF, PhillipsPerron, KPSS
print(arch.__version__)

        
        Copy
    

    	
            
def Stationarity_Table(data = [df["RMCPANNL"]["yoy"].dropna(),
                               df["RMUNPTOTP"]["raw"].dropna()],
                       dataset_names = ["RMCPANNL", "RMUNPTOTP"],
                       tests = ["ADF", "PP", "KPSS"],
                       enumerate_data = False):
    """Stationarity_Table Version 1.0:
    This function returns a Pandas data-frame with the results of the
    Augmented Dickey-Fuller (ADF), Phillips-Perron (PP) and the
    Kwiatkowski–Phillips–Schmidt–Shin (KPSS) stationarity tests from
    the arch library.
    
    
    Parameters:
    ----------
    
    data: list
        list of nan-less Pandas data-frames of number-series on which the test will be performed
        Default: data = [df["RMCPANNL"]["yoy"].dropna(), df["RMUNPTOTP"]["raw"].dropna()]
    
    dataset_names: list
        List of strings of the names to show as columns in our resulted table
        Default: dataset_names = ["RMCPANNL", "RMUNPTOTP"]
    
    tests: list
        List of strings of the stationarity tests wished.
        You may choose between "ADF", "PP", and/or "KPSS".
        Default: tests = ["ADF", "PP", "KPSS"]
    
    
    Dependencies:
    ----------
    
    pandas 1.2.3
    numpy 1.19.5
    arch 4.15 as ' from arch.unitroot import ADF, PhillipsPerron, KPSS '
    
    
    Examples:
    --------
    
    >>> import DatastreamDSWS as dsws
    >>> ds = dsws.Datastream(username = "insert dsws username here", password = "insert dsws password here")
    >>> 
    >>> df_1 = ds.get_data(tickers = "RMCPANNL,RMUNPTOTP", # be carefulnot to put spaces in between elements here, or else these spaces will be included in column names.
    >>>                    fields = "X", start = '2000-01-01', freq = 'M')
    >>> 
    >>> Stationarity_Table(data = [df_1["RMCPANNL"]["X"].dropna(), df_1["RMUNPTOTP"]["X"].dropna()], dataset_names = ["RMCPANNL", "RMUNPTOTP"], tests = ["ADF", "PP", "KPSS"], enumerate_data = False)
    """
    
    _columns, stationarity_test = [], []
    for i in tests:
        _columns.append((i, "Test-statistic", "")) # The ', ""' is there for the Pandas MultiIndex
        _columns.append((i, "p-value", ""))
        _columns.append((i, "Number of lags used", ""))
        _columns.append(
            (i,
             "Number of observations used for the regression and critical values'calculation",
             ""))
        _columns.append((i, "Critical values", '1%'))
        _columns.append((i, "Critical values", '5%'))
        _columns.append((i, "Critical values", '10%'))
        _columns.append((i, "Null hypothesis", ""))
    
    for d, _d in zip([D.astype(np.float64) for D in data], enumerate(data)):
        
        if enumerate_data == True:
            print(_d[0])
        
        _tests, _data = [], []
        if "ADF" in tests: _tests.append(ADF(d))
        if "PP" in tests: _tests.append(PhillipsPerron(d))
        if "KPSS" in tests: _tests.append(KPSS(d))
        
        for i in _tests:
            _data.append(i.stat)
            _data.append(i.pvalue)
            _data.append(i.lags)
            _data.append(i.nobs)
            for k in [i.critical_values[j] for j in i.critical_values]:
                _data.append(k)
            _data.append(i.null_hypothesis)
        
        stationarity_test.append(
            pd.DataFrame(
                data = np.array(_data)[np.newaxis],
                columns = pd.MultiIndex.from_tuples(_columns)).T)
    
    stationarity_test_table = pd.concat(
        stationarity_test, axis = 1, join = "inner")
    
    stationarity_test_table = pd.DataFrame(
        data = stationarity_test_table.values,
        index = stationarity_test_table.index,
        columns = dataset_names)
    
    return stationarity_test_table

        
        Copy
    

    	
            
Stationarity_Table(
    data =  [df["RMCPANNL"]["yoy"].dropna(),
             df["RMUNPTOTP"]["raw"].dropna(),
             df["RMUNPTOTP"]["yoy"].dropna(),
             df["RMXRUSD."]["raw"].dropna(),
             df["RMXRUSD."]["yoy"].dropna(),
             df["RMXREUR."]["raw"].dropna(),
             df["RMXREUR."]["yoy"].dropna(),
             df["RMGOVBALA"]["raw"].dropna(),
             df["RMGOVBALA"]["yoy"].dropna(),
             df["CLc1"]["Close Price 30D Moving Average"].dropna(),
             df["CLc1"]["yoy"].dropna()],
    dataset_names = ["RMCPANNL yoy",
                     "RMUNPTOTP raw","RMUNPTOTP yoy",
                     "RMXRUSD. raw", "RMXRUSD. yoy",
                     "RMXREUR. raw", "RMXREUR. yoy",
                     "RMGOVBALA raw", "RMGOVBALA yoy",
                     "CLc1 Close Price 30D Moving Average",
                     "CLc1 Close Price yoy"],
    tests = ["ADF", "PP", "KPSS"],
    enumerate_data = False)
# # for ' data ', one could use:
# [df[i][j].dropna() for i,j in zip(
#     ["RMCPANNL", "RMUNPTOTP", "RMUNPTOTP",
#      "RMXRUSD.", "RMXRUSD.", "RMXREUR.", "RMXREUR.",
#      "RMGOVBALA", "RMGOVBALA", "CLc1", "CLc1"],
#     ["yoy", "raw", "yoy", "raw", "yoy", "raw", "yoy",
#      "raw", "yoy", "Close Price 30D Moving Average", "yoy"])]

        
        Copy
    

    	
            
df["RMCPANNL", "yoy1d"] = df["RMCPANNL"]["yoy"].diff(1)
df["RMUNPTOTP", "yoy1d"] = df["RMUNPTOTP"]["yoy"].diff(1)
df["RMXRUSD.", "yoy1d"] = df["RMXRUSD."]["yoy"].diff(1)
df["RMXRUSD.", "yoy1d"] = df["RMXRUSD."]["yoy"].diff(1)
df["RMXREUR.", "yoy1d"] = df["RMXREUR."]["yoy"].diff(1)
df["RMGOVBALA", "yoy1d"] = df["RMGOVBALA"]["yoy"].diff(1)
df["CLc1", "yoy1d"] = df["CLc1"]["yoy"].diff(1)
df = df.sort_index(axis = 1)
 
df.dropna()

        
        Copy
    

    	
            
Stationarity_Table(
    data = [
        df[i]["yoy1d"].dropna()
        for i in [
            "RMCPANNL", "RMUNPTOTP", "RMXRUSD.",
            "RMXREUR.", "RMGOVBALA", "CLc1"]],
    dataset_names = [
        "RMCPANNL yoy1d", "RMUNPTOTP yoy1d",
        "RMXRUSD. yoy1d", "RMXREUR. yoy1d",
        "RMGOVBALA yoy1d", "CLc1 Close Price yoy1d"],
    tests = ["ADF", "PP"])

        
        Copy
    

    	
            
df["RMGOVBALA"]["yoy1d"].iplot(
    theme = "solar", xTitle = "date", yTitle = "yoy1d",
    title = "RMGOVBALA's year on year % increase's 1st difference")

        
        Copy
    

    	
            
import matplotlib.pyplot as plt
import statsmodels.api as sm

        
        Copy
    

    	
            
sm.graphics.tsa.plot_acf(x = df["RMCPANNL"]["yoy1d"].dropna(), lags = 25);
sm.graphics.tsa.plot_pacf(x = df["RMCPANNL"]["yoy1d"].dropna(), lags = 25);

        
        Copy
    

    	
            from statsmodels.tsa.arima.model import ARIMA
        
        Copy
    

    	
            
arima502c = ARIMA(endog = np.array(df["RMCPANNL"]["yoy1d"].dropna()),
                  order = (5, 0, 2), # We don't need to difference our data again, so 'd' is set to 0.
                  trend = "c")
arima502c_fit = arima502c.fit()

        
        Copy
    

    	
            arima502c_fit.summary()
        
        Copy
    

    	
            
s1arima502c = ARIMA(endog = np.array(df["RMCPANNL"]["yoy1d"].dropna()),
                    order = (5, 0, 2),
                    trend = "c",
                    seasonal_order = (1, 0, 0, 12))
s1arima502c_fit = s1arima502c.fit()

        
        Copy
    

    	
            s1arima502c_fit.summary()
        
        Copy
    

    	
            
models_resid = pd.DataFrame(
    data = {'arima502c residuals' : arima502c_fit.resid,
            's1arima502c residuals' : s1arima502c_fit.resid})
models_resid.iplot(
    theme = "solar", xTitle = "date",
    yTitle = "residuals (%)",
    title = "arima502c and s1arima502c models' residuals")

        
        Copy
    

    	
            
def IC_Table(models_fit, models_name):
    return pd.DataFrame(
        data = [[i.aic, i.aicc, i.bic, i.hqic]
                for i in models_fit],
        columns = ['AIC', 'AICC', 'BIC', 'HQIC'],
        index = [str(i) for i in models_name])
 
# ' ic ' standing for 'information criteria':
ic0 = IC_Table(models_fit = [arima502c_fit, s1arima502c_fit],
                 models_name = ['arima502c', 's1arima502c'])
 
ic0

        
        Copy
    

    	
            import itertools # We need this library to iteratethrough each permutation of exogenous variable possible
        
        Copy
    

    	
            
def Reg_Perm_Table(endog = np.array(df.dropna()["RMCPANNL"]["yoy1d"][1:].values.tolist()), # the ' [1:] ' is here to get same number of figures as in our exogenous variables.
                   exogenous = [df.dropna()[i]["yoy1d"].shift(1)[1:].values.tolist() # ' .shift(1) ' is here to get lagged values; the ' [1:] ' is there to remove the 0th value which is nan due to the lagging.
                                for i in ["CLc1", "RMGOVBALA", "RMUNPTOTP", "RMXREUR.", "RMXRUSD."]],
                   exogenous_names = ["CLc1", "RMGOVBALA", "RMUNPTOTP", "RMXREUR.", "RMXRUSD."],
                   d = [0], D = [0], s = [12],
                   p = range(5), q = range(5),
                   P = range(5), Q = range(5),
                   trend = "c",
                   print_progress = False,
                   try_loop = False,
                   save_model_fit = False,
                   lowest_no = 3):
    """ Reg_Perm_Table Version 1.0:
    This Python function ' Reg_Perm_Table ' goes through each permutation of possible ARIMA orders
    (chosen in parameters d, D, s, p, P, q, and Q) and exogenous variables and returns two Pandas
    data-frames (table and lowest) to help decide (i) the optimal exogenous variables to use,
    (ii) the optimal order of our ARIMA model to use.
    ' table ' consists of the Akaike information criterion (AIC),
    the AICc (AIC with a correction for small sample sizes),
    the Bayesian information criterion (BIC) (aka Schwarz criterion (also SBC, SBIC)) and
    the Hannan-Quinn information criterion (HQC)
    as well as the Mean Squared Errors (MSE) and
    the Sum of Squared Errors (SSE) (aka: residual sum of squares (RSS)) of any specified ARIMAX model.
    ' lowest ' consists of the name of all models with lowest values of all the metrics in ' table '.
    
    For more information on ARIMA orders (both (i) orders p, d, and q and (ii) seasonal orders P, D,
    Q and s), look at the statsmodels website:
    https://www.statsmodels.org/stable/generated/statsmodels.tsa.arima.model.ARIMA.html
    
    
    Parameters:
    ----------
    
    endog: Numpy array
        Numpy array of a list of values used as endogenous variables in our ARIMAX model.
        Defaulted to: np.array(df.dropna()["RMCPANNL"]["yoy1d"][1:].values.tolist()) # The ' [1:] ' is here to get same number of figures as in our exogenous variables.
    
    exogenous: list
        List of lists of each of the exogenous variables wished. can be set to ' None ' BUT
        IT CANNOT CONSIST OF ONLY ONE ELEMENT.
        Defaulted to: exogenous = [df.dropna()[i]["yoy1d"].shift(1)[1:].values.tolist() for i in ["CLc1", "RMGOVBALA", "RMUNPTOTP", "RMXREUR.", "RMXRUSD."]] # ' .shift(1) ' is here to get lagged values; the ' [1:] ' is there to remove the 0th value which is nan due to the lagging.
    
    exogenous_names: list
        list of strings of the names of the exogenous variables wished to be seen in
        the outputed data-frames. They have to be in the same order as in the ' exogenous ' variables.
        Defaulted to: exogenous_names = ["CLc1", "RMGOVBALA", "RMUNPTOTP", "RMXREUR.", "RMXRUSD."]
    
    p: list
        list of the ARIMA order 'p' that the user wishes to iteratethrough. This can be a range.
        Defaulted to: p = [0]
    
    d: list
        list of the ARIMA order 'd' that the user wishes to iteratethrough. This can be a range.
        Defaulted to: d = [0]
    
    q: list
        list of the ARIMA order 'q' that the user wishes to iteratethrough. This can be a range.
        Defaulted to: q = [0]
    
    P: list
        list of the ARIMA seasonal order 'P' that the user wishes to iteratethrough. This can be a range.
        Defaulted to: P = [0]
    
    D: list
        list of the ARIMA seasonal order 'D' that the user wishes to iteratethrough. This can be a range.
        Defaulted to: D = [0]
    
    Q: list
        list of the ARIMA seasonal order 'Q' that the user wishes to iteratethrough. This can be a range.
        Defaulted to: Q = [0]
    
    s: list
        list of the ARIMA seasonal order 's' that the user wishes to iteratethrough. This can be a range.
        Defaulted to: s = [0]
    
    trend: str
        str{‘n’,’c’,’t’,’ct’} or 'iterate-able', optional as per statsmodels.tsa.arima.model.ARIMA.
        Defaulted to: trend = "c"
    
    print_progress: Boolean
        If True, then a summary of each fitted model will appear in order of execution.
        Defaulted to: print_progress = False
    
    try_loop: Boolean
        If true, a try and except (then pass) loop is used over the model specification and
        fitting. If that part of the code fails, it will just ignore the model that caused an error from the output.
        Defaulted to: try_loop = False
    
    save_model_fit: Boolean
        ARIMA model fit's Python objects take up a lot of memory, so they are not saved. If
        this argument is st to True, our function then returns a new list full of each fit objects in the order they
        were executed called ' model_fit_list '; the function returns: table, model_fit_list, lowest
        Defaulted to: save_model_fit = False
    
    lowest_no: int
        the returned table ' lowest ' provides the name of each model with lowest AIC, AICC, BIC,
        HQIC, MSE and SSE up to the 3rd lowest value by default, but you can change that here.
        Defaulted to: lowest_no = 3
    
    
    Dependencies:
    ----------
    
    itertools
    numpy as np: version 1.19.5
    pandas as pd: version 1.2.3
    from statsmodels.tsa.arima.model import ARIMA: version 0.11.1
    
    
    Examples:
    --------
    
    >>> import DatastreamDSWS as dsws # Using LSEG Refinitiv's Datastream Python API named DataStreem Web Service (DSWS).
    >>> ds = dsws.Datastream(username = "insert dsws username here", password = "insert dsws password here")
    >>> 
    >>> df_1 = ds.get_data(tickers = "RMCPANNL,RMGOVBALA,RMXRUSD.", # be carefulnot to put spaces in between elements here, or else these spaces will be included in column names.
    >>>                    fields = "X", start = '2000-01-01', freq = 'M')
    >>> 
    >>> Reg_Perm_Table(endog = np.array(df_1.dropna()["RMCPANNL"]["X"][1:].values.tolist()), # the ' [1:] ' is here to get same number of figures as in our exogenous variables.
    >>>                exogenous = [df_1.dropna()[i]["X"].shift(1)[1:].values.tolist() # ' .shift(1) ' is here to get lagged values; the ' [1:] ' is there to remove the 0th value which is nan due to the lagging.
    >>>                             for i in ["RMXRUSD.", "RMGOVBALA"]],
    >>>                exogenous_names = ["RMXRUSD.", "RMGOVBALA"])
    """
    
    model_table_data, model_fit_list, models_name_list = [], [], [] # Lists to be appended.
    
    for _p in p: # Going through each permutation of possible arguments.
        for _d in d:
            for _q in q:
                for _P in P:
                    for _D in D:
                        for _Q in Q:
                            for _s in s:
                                
                                def Modelling(exog1, endog1 = endog,
                                              order1 = (_p, _d, _q), trend1 = trend,
                                              seasonal_order1 = (_P, _D, _Q, _s)):
                                    # global?
                                    return ARIMA(endog = endog1, exog = exog1,
                                                 order = order1, trend = trend1,
                                                 seasonal_order = (_P, _D, _Q, _s)).fit()
                                if exogenous != None:
                                    for j in range(1,len(exogenous_names)+1):
                                        for i,i_name in zip(itertools.combinations(exogenous, r = j),
                                                            itertools.combinations(exogenous_names, r = j)):
                                            
                                            if try_loop == True:
                                                try: # statsmodels' ARIMA function does not work 100% of the time.
                                                    model = Modelling(exog1 = np.array(i).T)
                                                    model_table_data.append(
                                                        [model.aic, model.aicc, model.bic,
                                                         model.hqic, model.mse, model.sse])
                                                    if save_model_fit == True: model_fit_list.append(model)
                                                    models_name_list.append(
                                                        f"sarima{_p},{_d},{_q},{_P},{_D},{_Q},{_s},{trend},exo={i_name}")
                                                    if print_progress == True: # For debugging.
                                                        print(models_name_list[-1])
                                                        display(model_fit_list[-1].summary())
                                                except:
                                                    pass
                                            else: # i.e.: if try_loop != True
                                                model = Modelling(exog1 = np.array(i).T)
                                                model_table_data.append(
                                                    [model.aic, model.aicc, model.bic,
                                                     model.hqic, model.mse, model.sse])
                                                if save_model_fit == True: model_fit_list.append(model)
                                                models_name_list.append(
                                                    f"sarima{_p},{_d},{_q},{_P},{_D},{_Q},{_s},{trend},exo={i_name}")
                                                if print_progress == True: # For debugging.
                                                    print(models_name_list[-1])
                                                    display(model_fit_list[-1].summary())
                                
                                elif exogenous == None:
                                    if try_loop == True:
                                        try: # statsmodels' ARIMA function does not work 100% of the time.
                                            model = Modelling(exog1 = None)
                                            model_table_data.append(
                                                [model.aic, model.aicc, model.bic,
                                                 model.hqic, model.mse, model.sse])
                                            if save_model_fit == True: model_fit_list.append(model)
                                            models_name_list.append(
                                                f"sarima{_p},{_d},{_q},{_P},{_D},{_Q},{_s},{trend},exo={i_name}")
                                            if print_progress == True: # For debugging.
                                                print(models_name_list[-1])
                                                display(model_fit_list[-1].summary())
                                        except:
                                            pass
                                        
                                        else: # i.e.: if try_loop != True
                                            model = Modelling(exog1 = None)
                                            model_table_data.append(
                                                [model.aic, model.aicc, model.bic,
                                                 model.hqic, model.mse, model.sse])
                                            if save_model_fit == True: model_fit_list.append(model)
                                            models_name_list.append(
                                                f"sarima{_p},{_d},{_q},{_P},{_D},{_Q},{_s},{trend},exo={i_name}")
                                            if print_progress == True: # For debugging.
                                                print(models_name_list[-1])
                                                display(model_fit_list[-1].summary())
    
    table = pd.DataFrame(
        data = model_table_data,
        columns = ['AIC', 'AICC', 'BIC', 'HQIC', "MSE", "SSE"],
        index = [str(i) for i in models_name_list]) # Just like the previously defined 'IC_Table' function.
    
    if lowest_no == None or lowest_no == False:
        if save_model_fit == True:
            return table, model_fit_list
        else:
            return table
    else:
        if lowest_no == 1:
            lowest_index = ["1st lowest"]
        elif lowest_no == 2:
            lowest_index = [f"{i} lowest" for i in ['1st', '2nd']]
        elif lowest_no == 3:
            lowest_index = [f"{i} lowest" for i in ['1st', '2nd', '3rd']]
        elif lowest_no > 3:
            lowest_index = [f"{i} lowest"
                            for i in list(['1st', '2nd', '3rd'] +
                                          [f"{j}th" for j in range(4, lowest_no + 1)])]
        else:
            print(f"{lowest_no} should be a positive integer")
        lowest = pd.DataFrame({'AIC': table.nsmallest(lowest_no, 'AIC').index,
                               'AICC': table.nsmallest(lowest_no, 'AICC').index,
                               'BIC': table.nsmallest(lowest_no, 'BIC').index,
                               'HQIC': table.nsmallest(lowest_no, 'HQIC').index,
                               'MSE': table.nsmallest(lowest_no, 'MSE').index,
                               'SSE': table.nsmallest(lowest_no, 'SSE').index},
                              index = lowest_index)
        if save_model_fit == True:
            return table, model_fit_list, lowest
        else:
            return table, lowest

        
        Copy
    

    	
            
# We get a lot of 'Maximum Likelihood optimization failed to converge.
# Check mle_retvals' error messages, this cell will remove them.
# You don't need to remove them, it's just neater to.
import warnings
from statsmodels.tools.sm_exceptions import ConvergenceWarning
warnings.simplefilter('ignore', ConvergenceWarning)

        
        Copy
    

    	
            
# This step takes a long while
test1_df, test1_lowest = Reg_Perm_Table(p = [1])

        
        Copy
    

    	
            test1_lowest
        
        Copy
    

    	
            
# need to ' pip install pickle-mixin '
import pickle
 
# To save data out:
pickle_out = open("Regression_Permutation_Table1.pickle","wb") # This creates the '.pickle' file where our data of choice will be saved. ' wb ' stand for 'write bytes'.
pickle.dump((test1_df, test1_lowest), # ' (test1_df, test1_lowest) ' are the two data-frames we choose to save.
            pickle_out) # ' pickle_out ' specifies the '.pickle' file in which we want to write (actually: overwrite - everything previously in that file will be overwritten)
pickle_out.close() # We need to close this '.pickle' file; leaving it open could corrupt it.
 
# # To load data in:
# pickle_in = open("Regression_Permutation_Table2.pickle","rb") # ' rb ' stand for 'read bytes'.
# test1_df, test1_lowest = pickle.load(pickle_in)
# pickle_in.close() # We ought to close the file we opened to allow any other programs access if they need it.

        
        Copy
    

    	
            
del test1_df
del test1_lowest

        
        Copy
    

    	
            
test2_df, test2_lowest = Reg_Perm_Table(p = [2])
pickle_out = open("Regression_Permutation_Table2.pickle","wb")
pickle.dump((test2_df, test2_lowest), pickle_out)
pickle_out.close()
del test2_df
 
test3_df, test3_lowest = Reg_Perm_Table(p = [3])
pickle_out = open("Regression_Permutation_Table3.pickle","wb")
pickle.dump((test3_df, test3_lowest), pickle_out)
pickle_out.close()
del test3_df
 
test4_df, test4_lowest = Reg_Perm_Table(p = [4], try_loop = True) # We tripped on some problems with this one, so used the ' try_loop '.
pickle_out = open("Regression_Permutation_Table4.pickle","wb")
pickle.dump((test4_df, test4_lowest), pickle_out)
pickle_out.close()
del test4_df
 
test5_df, test5_lowest = Reg_Perm_Table(p = [5])
pickle_out = open("Regression_Permutation_Table5.pickle","wb")
pickle.dump((test5_df, test5_lowest), pickle_out)
pickle_out.close()
del test5_df

        
        Copy
    

    	
            import pmdarima # Install via ' pip install pmdarima '
        
        Copy
    

    	
            
pmdarima.arima.auto_arima(
    y = np.array(df.dropna()["RMCPANNL"]["yoy1d"][1:]),
    X = None,
    d = 1, # The order of first-differencing. If None (by default), the value will automatically be selected based on the results of the test (i.e., either the Kwiatkowski–Phillips–Schmidt–Shin, Augmented Dickey-Fuller or the Phillips–Perron test will be conducted to find the most probable value). Must be a positive integer or None. Note that if d is None, the runtime could be significantly longer.
    start_p = 0, max_p = 7,
    start_q = 0, max_q = 7,
    D = None, max_D = 1,
    start_P = 0, max_P = 2,
    start_Q = 0, max_Q = 2,
    m = 12, # The period for seasonal differencing, m refers to the number of periods in each season. For example, m is 4 for quarterly data, 12 for monthly data, or 1 for annual (non-seasonal) data. Default is 1. Note that if m == 1 (i.e., is non-seasonal), seasonal will be set to False. For more information on setting this parameter, see Setting m.
    seasonal = True,
    stationary = False,
    information_criterion = 'aic',
    test = 'kpss', # Type of unit root test to use in order to detect stationarity if stationary is False and d is None. Default is ‘kpss’ (Kwiatkowski–Phillips–Schmidt–Shin).
    seasonal_test = 'ocsb', # This determines which seasonal unit root test is used if seasonal is True and D is None. Default is ‘OCSB’.
    stepwise = True, # Whether to use the stepwise algorithm outlined in Hyndman and Khandakar (2008) to identify the optimal model parameters. The stepwise algorithm can be significantly faster than fitting all (or a random subset of) hyper-parameter combinations and is less likely to over-fit the model.
    method = 'lbfgs', # The method determines which solver from scipy.optimize is used
    maxiter = 40)

        
        Copy
    

    	
            
# To load data in:
for j in range(1,6):
    exec(f"pickle_in = open('Regression_Permutation_Table{j}.pickle','rb')") # ' rb ' stand for 'read bytes'.
    exec(f"test{j}_df, test{j}_lowest = pickle.load(pickle_in)")
    exec(f"pickle_in.close()") # We ought to close the file we opened to allow any other programs access if they need it.
 
for j in range(1,6):
    # ' list(set( ' is there to return unique elements in case the same model has both the lowest MSE and SSE
    #   if you think about it, they should always be the same.
    exec(f"display(test{j}_df.loc[list(set([i for i in test{j}_lowest[['MSE', 'SSE']].loc['1st lowest']]))])")

        
        Copy
    

    	
            
s4arima504c = ARIMA(endog = np.array(df.dropna()["RMCPANNL"]["yoy1d"][1:].values.tolist()),
                    exog = np.array([df.dropna()[i]["yoy1d"].shift(1)[1:].values.tolist()
                                     for i in ["RMUNPTOTP", "RMXRUSD."]]).T,
                    order = (5, 0, 4),
                    trend = "c",
                    seasonal_order = (4, 0, 0, 12))
s4arima504c_fit = s4arima504c.fit()
s4arima504c_fit.summary()

        
        Copy
    

    	
            
s4arima504c_fit_plot_dic = {
    'residuals' : s4arima504c_fit.resid,
    'fitted values' : s4arima504c_fit.fittedvalues,
    'endogenous variable (RMCPANNL)' : np.array(df.dropna()["RMCPANNL"]["yoy1d"][1:].values.tolist())}
 
exogenous_variables = ["RMUNPTOTP", "RMXRUSD."]
for j,k in zip([df.dropna()[i]["yoy1d"].shift(1)[1:].values.tolist()
                for i in exogenous_variables],
               exogenous_variables):
    s4arima504c_fit_plot_dic[k] = np.array(j).T
    
s4arima504c_fit_plot_df = pd.DataFrame(data = s4arima504c_fit_plot_dic,
                                       index = df.dropna()["RMCPANNL"]["yoy1d"][1:].index)
 
s4arima504c_fit_plot_df.iplot(theme = "solar",
                              title = "s4arima504c model",
                              xTitle = "date")

        
        Copy
    

    	
            
s4arima504c_osaoosf = []
 
for t in range(200, len(df.dropna()["RMCPANNL"]["yoy1d"][1:])):
    if t%10 == 0 : print(str(t)) # This line is to see the progression of our code.
    s4arima504c = ARIMA(
        endog = np.array(df.dropna()["RMCPANNL"]["yoy1d"][1:].values.tolist())[0:t-1],
        exog = np.array([df.dropna()[i]["yoy1d"].shift(1)[1:].values.tolist()[0:t-1]
                         for i in ["RMUNPTOTP", "RMXRUSD."]]).T,
        order = (5, 0, 4),
        trend = "c",
        seasonal_order = (4, 0, 0, 12))
    
    s4arima504c_isf = s4arima504c.fit() # 'isf' for in-sample-fit
    
    s4arima504c_osaoosf.append( # 'osaoosf' as in One-Step-Ahead Out-of-Sample Forecasts
        s4arima504c_isf.forecast(
            exog = np.array([df.dropna()[i]["yoy1d"].shift(1)[1:].values.tolist()[t]
                             for i in ["RMUNPTOTP", "RMXRUSD."]]).T)[0])

        
        Copy
    

    	
            
# saving data into a pickle file::
pickle_out = open("s4arima504c_osaoosf.pickle","wb")
pickle.dump(s4arima504c_osaoosf, pickle_out)
pickle_out.close()
 
# # loading data from the pickle file:
# pickle_in = open('s4arima504c_osaoosf.pickle','rb') # ' rb ' stand for 'read bytes'.
# s4arima504c_osaoosf = pickle.load(pickle_in)
# pickle_in.close() # We ought to close the file we opened to allow any other programs access if they need it.

        
        Copy
    

    	
            
RMCPANNL_array_wna = np.append( # 'wna' as in 'with na' value added at the end.
    np.array(df.dropna()["RMCPANNL"]["yoy1d"][1:].values.tolist())[200 + 1:],
    [np.nan])

        
        Copy
    

    	
            from datetime import datetime
        
        Copy
    

    	
            
# add a month to the index:
index_in_qu = df.dropna()["RMCPANNL"]["yoy1d"][1:][200:].index
index = [i for i in index_in_qu][1:] + [
    str(datetime.strptime(str(pd.to_datetime(
        [str(i) for i in  index_in_qu])[-1].replace(
        month = pd.to_datetime(str(df.index[-1])).month+1)),
                          '%Y-%m-%d %H:%M:%S'))[:10]]

        
        Copy
    

    	
            
s4arima504c_osaoosf_plot_dic = {
    'Recursive One-Step-Ahead Out-of-Sample Forecasts' : np.array(s4arima504c_osaoosf),
    'endogenous variable (RMCPANNL)' : RMCPANNL_array_wna}
 
exogenous_variables = ["RMUNPTOTP", "RMXRUSD."]
for j,k in zip([df.dropna()[i]["yoy1d"].shift(1)[1:].values.tolist()[200:]
                for i in exogenous_variables],
               exogenous_variables):
    s4arima504c_osaoosf_plot_dic[k] = np.array(j).T
 
s4arima504c_osaoosf_plot_df = pd.DataFrame(
    data = s4arima504c_osaoosf_plot_dic,
    index = index)

        
        Copy
    

    	
            
s4arima504c_osaoosf_plot_df.iplot(
    theme = "solar", xTitle = "date",
    title = "s4arima504c Recursive One-Step-Ahead Out-of-Sample Forecasts")

        
        Copy
    

    	
            df2 = df.copy()
        
        Copy
    

    	
            s4arima504c_osaoosf_plot_df.tail()
        
        Copy
    

    	
            
df2["RMCPANNL", "s4arima504c Recursive OSAOOSF"] = s4arima504c_osaoosf_plot_df["Recursive One-Step-Ahead Out-of-Sample Forecasts"]
 
df2 = df2.sort_index(axis=1) # This rearranges the columns and merges Fields under the same Instruments.
 
df2
 

        
        Copy
    

    	
            
df3 = df2.copy() # just to delimitate from before and after this step.
 
index_of_interest = s4arima504c_osaoosf_plot_df.tail(1).index[0]
df3 = df3.append(pd.Series(name = index_of_interest))
 
df3['RMCPANNL', 's4arima504c Recursive OSAOOSF'].loc[index_of_interest] = s4arima504c_osaoosf_plot_df["Recursive One-Step-Ahead Out-of-Sample Forecasts"][-1]
 
df3

        
        Copy
    

    	
            
test0_df, test0_lowest = Reg_Perm_Table(p = [0])
 
# To save data out:
pickle_out = open("Regression_Permutation_Table0.pickle","wb")
pickle.dump((test0_df, test0_lowest), pickle_out)
pickle_out.close()
 
# To load data in:
# pickle_in = open("Regression_Permutation_Table0.pickle","rb") # ' rb ' stand for 'read bytes'.
# test0_df, test0_lowest = pickle.load(pickle_in)
# pickle_in.close() # We ought to close the file we opened to allow any other programs access if they need it.
 
test0_df[test0_df.index.str.contains("sarima0,0,0")]

        
        Copy
    

    	
            test0_df[test0_df.index.str.contains("sarima0,0,0")].nsmallest(1, 'SSE').index
        
        Copy
    

    	
            
exogenous_variables = ['RMUNPTOTP', 'RMXREUR.', 'RMXRUSD.']
 
s3arima000c = ARIMA(
    endog = np.array(df.dropna()["RMCPANNL"]["yoy1d"][1:].values.tolist())[0:200],
    exog = np.array([df.dropna()[i]["yoy1d"].shift(1)[1:].values.tolist()[0:200]
                     for i in exogenous_variables]).T,
    order = (0, 0, 0),
    trend = "c",
    seasonal_order = (3, 0, 0, 12))
s3arima000c_isf = s3arima000c.fit() # 'isf' for in-sample-fit
 
s3arima000c_isf_plot_dic = {
    'residuals' : s3arima000c_isf.resid,
    'fitted values' : s3arima000c_isf.fittedvalues,
    'endogenous variable (RMCPANNL)' : np.array(
        df.dropna()["RMCPANNL"]["yoy1d"][1:].values.tolist())[0:200]}
 
for j,k in zip([df.dropna()[i]["yoy1d"].shift(1)[1:].values.tolist()[0:200]
                for i in exogenous_variables],
               exogenous_variables):
    s3arima000c_isf_plot_dic[k] = np.array(j).T
    
s3arima000c_isf_plot_df = pd.DataFrame(data = s3arima000c_isf_plot_dic,
                                       index = df.dropna()["RMCPANNL"]["yoy1d"][1:][0:200].index)
 
s3arima000c_isf_plot_df.iplot(theme = "solar",
                              title = "s3arima000c in-sample-forecasts model fit",
                              xTitle = "date")

        
        Copy
    

    	
            
s3arima000c_osaoosf = []
 
for t in range(200, len(df.dropna()["RMCPANNL"]["yoy1d"][1:])):
    if t%10 == 0 : print("Our code has created one-step-ahead out-of-sample forecasts from 200 to " str(t)) # This line is to see the progression of our code.
    s3arima000c = ARIMA(
        endog = np.array(df.dropna()["RMCPANNL"]["yoy1d"][1:].values.tolist())[0:t-1],
        exog = np.array([df.dropna()[i]["yoy1d"].shift(1)[1:].values.tolist()[0:t-1]
                         for i in exogenous_variables]).T,
        order = (0, 0, 0),
        trend = "c",
        seasonal_order = (3, 0, 0, 12))
    
    s3arima000c_isf = s3arima000c.fit() # 'isf' for in-sample-fit
    
    s3arima000c_osaoosf.append( # 'osaoosf' as in One-Step-Ahead Out-of-Sample Forecasts
        s3arima000c_isf.forecast(
            exog = np.array([df.dropna()[i]["yoy1d"].shift(1)[1:].values.tolist()[t]
                             for i in exogenous_variables]).T)[0])

        
        Copy
    

    	
            
# saving data into a pickle file::
pickle_out = open("s3arima000c_osaoosf.pickle","wb")
pickle.dump(s3arima000c_osaoosf, pickle_out)
pickle_out.close()
 
# # loading data from the pickle file:
# pickle_in = open('s3arima000c_osaoosf.pickle','rb') # ' rb ' stand for 'read bytes'.
# s3arima000c_osaoosf = pickle.load(pickle_in)
# pickle_in.close() # We ought to close the file we opened to allow any other programs access if they need it.

        
        Copy
    

    	
            
# add a month to the index:
index_in_qu = df.dropna()["RMCPANNL"]["yoy1d"][1:][200:].index
index = [i for i in index_in_qu][1:] + [
    str(datetime.strptime(str(pd.to_datetime(
        [str(i) for i in  index_in_qu])[-1].replace(
        month = pd.to_datetime(str(df.index[-1])).month+1)),
                          '%Y-%m-%d %H:%M:%S'))[:10]]

        
        Copy
    

    	
            
s3arima000c_osaoosf_plot_dic = {
    'Recursive One-Step-Ahead Out-of-Sample Forecasts' : np.array(s3arima000c_osaoosf),
    'endogenous variable (RMCPANNL)' : np.append(
        np.array(df.dropna()["RMCPANNL"]["yoy1d"][1:].values.tolist())[200 + 1:],
        [np.nan])}
 
for j,k in zip([df.dropna()[i]["yoy1d"].shift(1)[1:].values.tolist()[200:]
                for i in exogenous_variables],
               exogenous_variables):
    s3arima000c_osaoosf_plot_dic[k] = np.array(j).T
 
s3arima000c_osaoosf_plot_df = pd.DataFrame(
    data = s3arima000c_osaoosf_plot_dic,
    index = index)

        
        Copy
    

    	
            
s3arima000c_osaoosf_plot_df.iplot(
    theme = "solar", xTitle = "date",
    title = "s3arima000c Recursive One-Step-Ahead Out-of-Sample Forecasts")

        
        Copy
    

    	
            
for i,j in zip([s3arima000c_osaoosf_plot_df], [s4arima504c_osaoosf_plot_df]):
    s3arima000c_resid = i["endogenous variable (RMCPANNL)"] - i["Recursive One-Step-Ahead Out-of-Sample Forecasts"][:-1]
    s4arima504c_resid = j["endogenous variable (RMCPANNL)"] - j["Recursive One-Step-Ahead Out-of-Sample Forecasts"][:-1]
    s3arima000c_sse = (s3arima000c_resid**2).sum()
    s4arima504c_sse = (s4arima504c_resid**2).sum()

        
        Copy
    

    	
            print(s3arima000c_sse, s4arima504c_sse)
        
        Copy