To do so we look at (i) how one may check for stationarity graphically and via (ADF, KPSS & PP) test statistics, (ii) differencing, (iii) ex-ante (before modelling) parameter identification via autocorrelation & partial autocorrelation functions, (iv) the difference between Autoregressive Moving Average (ARMA) & Seasonal Integrated ARMA (SARIMA) models, (v) ex-post (after modelling) parameter identification via information criteria and mean of absolute & squared errors (as well as the use of Python's 'pickle' library), (vi) how one may choose an optimal model specification - reducing errors, (vii) recursive one-step-ahead out-of-sample forecasts, and finally (iix) corresponding month models before concluding.

As a precursor to ARMA modelling, the reader is advised to read 'Information Demand and Stock Return Predictability (Coded in R): Part 2: Econometric Modelling of Non-Linear Variances' - an article without coding that looks into our methods of modelling.

Webinar

Contents

• What are SARIMA Models?
• Get to Coding
  • Development Tools & Resources
  • Import Libraries
• Data Sets
  • Collecting Datastream Data
  • Collecting Eikon data
  • Joining the Datastream and Eikon data-frames together
• Stationarity
  • Checking For Stationarity Via Graphical Analysis
  • Checking For Stationarity Via Test Statistics
    • Interpreting the stationarity table
  • Differencing
• Modelling
  • Ex-Ante Parameter Identification: Autocorrelation and Partial Autocorrelation Functions
    • ARMA(5,2)
    • SARIMA(5,0,2)(1, 0, 0)12
  • Ex-Post Parameter Identification: Information Criteria and Mean of Absolute & squared Errors
    • Graphical Analysis
    • Information Criteria
    • Mean Squared Errors and Sum of squared Errors
      • Pickle
      • pmdarima
  • Optimal Model
    • Fitted values of our 'Optimal Model'
    • Recursive One-Step-Ahead Out-of-Sample Forecasts
  • corresponding Month Model
    • Check the best SAMIMA model with ARMA(0,0) imposed
    • Fit the sarima0,0,0,3,0,0,12 model
    • Recursive One-Step-Ahead Out-of-Sample Forecasts
  • Comparing Models
• Conclusion
• Webinar Recording
• References

What are SARIMA Models? 

Seasonal AutoRegressive Integrated Moving Average (SARIMA) models are best explained in a separate article. "Time Series Forecasting with SARIMA in Python" is great at showing both the mathematics and fundamental Python basics at play below. I would suggest reading this article 1st that goes through the slightly simpler ARMA model. The only difference is the inclusion of the seasonal factor in:

Get to Coding

Development Tools & Resources

The example code demonstrating the use case is based on the following development tools and resources:

• Refinitiv's DataStream Web Services (DSWS): Access to DataStream data. A DataStream or Refinitiv Workspace IDentification (ID) will be needed to run the code below.
• Refinitiv's Eikon API: To access Eikon data. An Eikon API key will be needed to use this API.

Import Libraries

We need to gather our data. Since Refinitiv's DataStream Web Services (DSWS) allows for access to ESG data covering nearly 70% of global market cap and over 400 metrics, naturally it is more than appropriate. We can access DSWS via the Python library "DatastreamDSWS" that can be installed simply by using pip install.

We use Refinitiv's Eikon Python Application Programming Interface (API) to access financial data such as Clc1. We can access it via the Python library "eikon" that can be installed simply by using pip install.

The library 'eikon' in this notebook is version 1.1.8
The library 'numpy' in this notebook is version 1.18.2
The library 'pandas' in this notebook is version 1.2.4
The library 'plotly' in this notebook is version 4.14.3
The library 'cufflinks' in this notebook is version 0.17.3
The library 'seaborn' in this notebook is version 0.11.1

Data-sets 

• Dependent variables (Y):

  • RMCPANNL: Romania, Consumer Prices, by Commodity, All Items, Total, Change y/y

    Annual inflation rate – the increase (in percentage term) in consumer prices in one month of the current year compared to the same month of the previous year. This rate is calculated as a ratio, expressed as a percentage, between the price index of one month of the current year and that of the same month of the previous year, calculated against the same base, from which 100 is subtracted.

    This data-set will be differenced once to be stationary.

  • RMCONPRCF

• Independent / Explanatory variables (X):

  Each of the following datasetes can be used in 'raw' (id est (i.e.): without changing the data retrieved from Datastream) or adjusted like RMCPANNL (i.e.: looking at that metric's increase in one month compared to the same month of the previous year and then differenced once to be stationary) each defined as 'raw' or 'yoy' variables respectively.

  • RMUNPTOTP: Romania, Unemployed, Overall, Registered.
  • CLc1: 'CLc1 NYMEX Light Sweet Crude Oil (WTI) Electronic Energy Future Continuation 1' on eikon (which we will use) and 'NYMEX - Light Crude Oil TRC1' on Datastream: Crude oil prices from the US market's NYMEX (New York Mercantile Exchange).
  • RMXRUSD.: Romania, New Romanian Leu Per US Dollar, RON: Foreign Exchange (FX) rate of 1USD to 1Leu.
  • RMXREUR.: Romania, New Romanian Leu Per Euro, RON: Foreign Exchange (FX) rate of 1EUR to 1Leu.
  • RMGOVBALA: Romania, Deficit / Surplus, GFS2001, Total (New Methodology Since January 2006), Cumulative, RON: Romanian Government's National Debt level. If this value is positive, it indicates a surplus; if it is negative, it indicates a deficit.

Due to popular demand, we decided not to use Consumer Prices (Romania, Consumer Prices, by Commodity, All Items, Total, Index, Corresponding Period of the Previous Year = 100) as a dependent variable, but RMCPANNL instead.

Collecting Datastream Data 

We can go ahead and collect our data from Datastream:

Note that RMCPANNL's values are percentages, we thus ought to divide it by 100 to get proper fractions:

Let's have a look at our Pandas data-frame ' df '

256 rows × 5 columns

Note on timing: I aim to make this article as realistic as possible, putting myself in the shoes of a professional using this workflow to predict inflation. I do not expect our SARIMA models to perform particularly well, mainly because we did not do extensive research on the best exogenous variables to use - I leave that to you to decide. But while I do not expect this model to perform well, the workflow ought to be at its best. To use it optimally, please try to use exogenous variables (such as RMXRUSD.) as soon as their monthly values are released. It is possible for the values as of the - 15th of the month - of one variable to be published at another time to another (exempli gratia (e.g.): values for RMXRUSD. and RMGOVBALA could be published on the 20th and 25th of the month respectively, even though they are both representing values for the 15th of the month).

The 'X' is the default value for each of the items we requested. Let's change those with relevant value names:

Collecting Eikon data 

Now let's add our CLc1 data from Eikon. Note that we want monthly data, so let's collect monthly (30 day) moving average data (MAVG(TR.CLOSEPRICE,-29)). 1st let's collect that data:

5597 rows × 4 columns

But note that the last few moving average data points are NaN. Just in case these last few days are useful, let's make our own moving average:

5597 rows × 4 columns

Joining the Datastream and Eikon data-frames together 

Now, if we want to combine the data from df and CLc1, we may want to do so joining the two data-frames on their dates; but one may notice that CLc1's Timestamp and df's index (Dates) portray time in different ways:

'1999-01-01T00:00:00Z'

'2000-01-15'

To harmonize the two, I will use df's way:

'2000-01-31'

5597 rows × 4 columns

Now we are encountering the issue whereby the 15th of the month every month is not necessarily a trading day. But we want to normalise the two data-frames to have monthly data on the 15th of every month (since that's how Datastream returns data). What we can do is use CLc1 30D Moving Average on the 15th of every month and (if it's not available) we'll use the 14th's; if that's unavailable we'll use the 13th; and so on until a set limit - say the 12th of every month (repeating this process 3 times). This is what the function below does:

As per the above, we can see the difference between the index and Timestamp in row 259, indicating that our function above worked.

265 rows × 2 columns

Now we can join our two data-frames:

256 rows × 7 columns

Now, we want comparable data, and RMCPANNL is the only series depicting the increase in one month compared to the same month of the previous year. So let's transform our other series to depict the data similarly:

242 rows × 12 columns

Stationarity 

Checking For Stationarity Via Graphical Analysis 

We must admit that our RMCPANNL dataset does not look stationary

RMCPANNL's 1st difference (1d) does look rather stationary:

Upon quick graphical visualisation, it seems as though we ought to choose RMCPANNL's 1st difference

Checking For Stationarity Via Test Statistics 

We can create a Python function for the three aforementioned test-statistics for stationarity: Augmented Dickey-Fuller (ADF), Phillips-Perron (PP) and the Kwiatkowski–Phillips–Schmidt–Shin (KPSS) ones.

For the ADF test, one may use the adfuller function from the statsmodels library (imported via the line from statsmodels.tsa.stattools import adfuller), but that library doesn't have functions for all the tests in question. the arch library does - however. So I will define a Python function Stationarity_Table returning all relevant arch stationarity results of interest:

4.19

Interpreting the stationarity table: As per the above, the null hypothesis of the ADF is that there is a unit root, with the alternative that there is no unit root. If the P-value is above a critical size, then we cannot reject that there is a unit root. If the process has a unit root, then it is a non-stationary time series. If the p-value is close to significant, then the critical values should be used to judge whether to reject the null.

Nonetheless, analysing the table's figures for RMUNPTOTP, RMXRUSD, RMXREUR, and RMGOVBALA shows us that they're non-stationary (for the most part - especially with the ADF test). We ought to difference them once and see if they are stationary then.

Differencing 

241 rows × 18 columns

Now they're all stationary. We ought to use datasets treated similarly for consistency (e.g.: compare 'RMCPANNL yoy1d' with 'RMUNPTOTP yoy1d' rather than with 'RMUNPTOTP yoy'), so we will use the ones shown in the last table (i.e.: 'yoy1d' datasets).

N.B.: The RMGOVBALA yoy1d dataset is causing issues with the KPSS test, but that one is only for weak stationarity, and is not of as much interest to us in our study - I just wanted to include it to portray the extent of the Stationarity_Table function. The RMGOVBALA yoy1d dataset is causing issues with the KPSS test because of its large variance, as can be seen in this graph:

Modelling 

Now we will look at which ARMA-family model to implement. We'll look at the best model orders to choose before (ex-ante) and after (ex-post) going through modelling procedures.

Ex-Ante Parameter Identification: Autocorrelation and Partial Autocorrelation Functions 

To figure out our parameter identification (i.e.: what lags to use in our model), we may want to plot ACF and PACF diagrams:

N.B.: the semicolon at the end of the code for our plots below is there to avoid duplicates.

From these graphs (also known as (a.k.a.): correlograms), we could identify an AR(1,4,5) (non-consecutive) model; but the MA part is harder to decipher. Let's have a look at how those models would fair:

• 1st: a simple ARMA(5,2) model (i.e.: a SARIMA(5,0,2) model), then
• 2nd: a SARIMA(5,0,2)(1, 0, 0)12 model:

ARMA(5,2)

Warnings:

[1] Covariance matrix calculated using the outer product of gradients (complex-step).

SARIMA(5,0,2)(1, 0, 0)12 

Warnings:

[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Ex-Post Parameter Identification: Information Criteria and Mean of Absolute & squared Errors 

Graphical Analysis 

After having fit a model, one may try to look at the residuals to compare models:

Looking at such graphical information is not simple. (i) Not only is it arbitrary to chose a model over another from graphical information (ii) it quickly becomes difficult when comparing many models together. That's why we may compare them instead with performance measurements such as Information Criteria:

Information Criteria 

We can similarly go through each permutation within limits of our SARIMA model's parameters, and go beyond just recording IC information - including mean squared error (MSE) and sum of squared error (SSE, a.k.a.: RSS) information too.

Mean Squared Errors and Sum of squared Errors 

We can define the below function to return a table 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, and then choose the one providing the lowest (of any one chosen error measure above) as the optimal model of choice:

N.B.: Note that we need to lag our exogenous variables since we are putting ourselves in the shoes of a programmer who is using past values to predict tomorrow's.

Now let's put this Python function into practice:

Let's have a look at that table returned:

3875 rows × 6 columns

Let's have a look at that table lowest returned:

Pickle 

It is quite time consuming to request DSWS for all these codes. Furthermore, the table test1_df is rather large, and we only accounted for p, q, P, Q values up to (but not including) 3. What we could do is save our table in a file outside our code, then delete our table in our code, continue our analysis with test2_df and compare the two afterwards. Let's save our progress this far using Pickle. Since this is a step that can overwrite existing files, I commented them out. Feel free to uncomment them back in if you would like to save your data:

As aforementioned, let's delete our tables in our code for now - after having saved them outside our code:

Now let's continue with our examination:

pmdarima 

If we knew which exogenous variables to use, we could have used pmdarima to find the best model specification too:

ARIMA(maxiter=40, order=(0, 1, 1), scoring_args={},
seasonal_order=(0, 0, 1, 12), suppress_warnings=True)

Unfortunately, we don't know what permutation of exogenous variables is best, and Information-Criteria-optimisation methods are thus not sufficient.

Optimal model 

Finally, let's compare our two tables to check the lowest MSE and SSE:

Thus, the best model out of our choices is the 'sarima5,0,4,4,0,0,12,c,exo=('RMUNPTOTP', 'RMXRUSD.')' model (best as in SSE and MSE reducing):

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Looking at P-Values, we can see that the only coefficients with statistical significance are ar.L5 and ar.S.L12 to ar.S.L36.

Fitted values of our 'Optimal Model' 

Recursive One-Step-Ahead Out-of-Sample Forecasts 

Now we can simply apply our model but train it on a smaller dataset (say 80% of our data, i.e.: approx. 200 first observations) to simulate an out-of-sample test:

200
210
220
230

Let's save our data in s4arima504c_osaoosf:

242 rows × 19 columns

As you can see, we are missing the last row of data in df2. We can remedy this by appending it with the last row of s4arima504c_osaoosf_plot_df:

243 rows × 19 columns

Corresponding Month Model 

Many in the industry only focus on the same month of previous years in their predicting models (e.g.: only predict this year's March inflation value with previous March values, similarly so with each month of the year). There is some merit to the thought, as other datapoints may sway our data with irrelevant variables. The reason is that demand/supply profiles might be different each month of the day and it would be incorrect to compare months in a continuous series (i.e.: January to February et cetera (etc.)). To do so, we can exploit the SARIMA model and only use Seasonal AutoRegressive variables (in addition to our exogenous variables):

1st: Let's check the best SAMIMA model with ARMA(0,0) imposed:

775 rows × 6 columns

Index(['sarima0,0,0,3,0,0,12,c,exo=('RMUNPTOTP', 'RMXREUR.', 'RMXRUSD.')'], dtype='object')

Note that our Reg_Perm_Table function's default Seasonal parameters (P and Q) range from 0 to 5 and that the optimum (i.e.: lowest SSE producing one) is lower than 5; this would indicate that we may very well have our optimal SAMIMA model with ARMA(0,0) imposed

2nd: Fit the sarima0,0,0,3,0,0,12 model:

3rd: Recursive One-Step-Ahead Out-of-Sample Forecasts

200
210
220
230

Let's save our data in s3arima000c_osaoosf:

Comparing Models 

0.0013175028279920704 0.0012139805753142835

From the two SSE values above, we can determine that our s4arima504c model is the better of the two.

Conclusion 

The fact that a model with seasonal components (SARIMA(5,0,4)(4,0,0)1212) was the optimal one (even before we investigate Corresponding Month Models) suggests that there was a lot of information in corresponding months in previous years. As a matter of fact, even the directly previous months' coefficients seemed of little relevance. Indeed, s4arima504c_fit P-values of these cofactors from lag 1 to 4 (ar.L1 to ar.L4) never went below 0.1; i.e.: they were never statistically significant. The 5th lag however (ar.L5, with value 0.2342) had a P-value lower than 0.05, suggesting it to be statistically significant at the 5% Confidence Level (CL). Moving Average elements (ma.L1 to ma.L4) never showed statistical significance (as per their P-values). (Statistical significance is only really considered at the 10% CL, when the P-value is lower than 0.1.)

Seasonal coefficients 1 to 3 (ar.S.L12 to ar.S.L36) proved significant at the 5% CL. The 4th (ar.S.L48) did not.

It would be interesting to look at a non-consecutive lags model with only ar.L5 & ar.S.L12 to ar.S.L36 and investigate it's out-of-sample one-step-ahead forecasts as a back-test.

Additionally, what make this model interesting is that one can use it for more than one-step-ahead forecasts depending on the frequency of exogenous data provision. Say that the exogenous RMXRUSD variable is only published on the 30th of the month (showing data for the 15th of that month); than one can use data published a month prior to forecast two-steps-ahead, alleviating this issue - although not solving it per se (by itself).

Romanian inflation is a particular case study. High inflation levels a the start of our data-set might have skewed results significantly. A look at models such as the ones above in predicting other countries' inflation levels could be more fruitful. Alternatively, looking into more differenced data could be equally interesting. (Indeed, using Stationarity_Table with more differenced data seem to yield more stationary data.)

References

Model Specifications

• Time Series Forecasting with SARIMA in Python
• Problem defining ARIMA order
• Model selection by The Akaike’s Information Criterion (AIC) what is common practice?
• pmdarima Python Library
  • pmdarima.arima
    • Tips to using auto_arima
  • pmdarima.arima.auto_arima
• In general, choosing the order of the model (e.g. p, q and P, Q) via information criteria is a good idea, but it is not a good idea to choose the differencing order (d or D) that way.
• Time Series Talk : Seasonal ARIMA Model

Model Performance

• How to find relative influence of Exogenous variables from a ARIMA model?
• Out of sample risk estimates are the gold standard for performance evaluation

Academic Refferences

• FORECASTING INFLATION THROUGH ECONOMETRIC MODELS: AN EMPIRICAL STUDY ON PAKISTANI DATA
• Quarterly inflation rate target and forecasts in Romania
• THE EVOLUTION OF CONSUMER PRICES IN JANUARY 202
• Concepte de bază și definiții

Multi-Processing

• Python Multiprocessing Tutorial: Run Code in Parallel Using the Multiprocessing Module
• Parallel Processing in Python – A Practical Guide with Examples
• Multiprocessing in Python on Windows and Jupyter/Ipython — Making it work

• SEM - Comparing Models with Different Variables

Misalanious

• Display number with leading zeros
• Ordered_unique_values_in_list
• Select by partial string from a pandas DataFrame
• Jupyter notebook not running code. Stuck on In [*]
• Modelling