# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + id="Jwcln1tUyyHP"
import warnings
import itertools
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller
from sklearn.metrics import mean_squared_error
from pmdarima import auto_arima


# + [markdown] id="iHRDPm_4zIGz"
# # Auto-regressive integrated moving average models

# + [markdown] id="VEE0mD-9wIne"
# ## Introduction
# ARIMA stands for Auto-Regressive Integrated Moving Average. Let's go through each one of these concepts one by one. 
#
# ### Auto-regressive
# $$y_t = c + \rho_1  y_{t-1} + \rho_2  y_{t-2} + e_t$$
#
# The autoregressive part of the model assumes that the best way to explain the current value in a time series is using a set number of the past observations. In this case, we are looking at a AR(2) model, and $\rho_1$ and $\rho_2$ are the parameters of interest for the model. 
#
# ### Moving average
# The moving average component assumes that the best explanation for a time series is the error term of the previous values in the series:
# $$y_t = c + \mu_1  e_{t-1} + \mu_1  e_{t-2} + e_t$$
#
# This would be a MA(2) model, and $\mu_1$ and $\mu_2$ are the parameters of interest of the model. 
#
# ### Integrated
# This value determines the amount of times we differenciate the series ($y_t-y_{t-1}$) for the AR and MA models above. 
#
# The final equation depends on the integrated term. 
# For more information on the functional form of the seasonal ARIMA, please check the `statsmodel` documentation. 
#
#
#
#
#
#

# + [markdown] id="vR-JI2-uIEER"
# The parameters we are going to use below are $p$ for the AR component, $d$ for the integrated component, and $q$ for the MA component. 
#
#

# + [markdown] id="9uTlwgCwzPZi"
# # Data processing

# + colab={"base_uri": "https://localhost:8080/", "height": 833} id="ZEHZ3sfKy2Ck" outputId="98b19c58-37d4-4603-b899-cb6a5d6f3530"
# For this tutorial we are going to use the CO2 data available in the 
# Statsmodel package
data = sm.datasets.co2.load_pandas()
y = data.data
y.head(25)

# + colab={"base_uri": "https://localhost:8080/"} id="caQOBGO0y5XW" outputId="3f1fd648-255e-4ee8-81e8-2ae0a2315caf"
# We would like to reorganise the dataset in monthly data.
# The .resample() function indexes the data 
# M for monthly data, S for start of the month
# Since there is more than one observation for every month, we consider the mean
y = y['co2'].resample('MS').mean()
y.head(25)

# + colab={"base_uri": "https://localhost:8080/"} id="h1-HFwLEtbzi" outputId="58038586-559e-44be-9821-f51038db4bfb"
# There are still some missing values. We are going to fill them 
# using the last available value. There are other possibilities, but that 
# is not the focus of this session
y = y.fillna(y.bfill())

y.head()

# + colab={"base_uri": "https://localhost:8080/", "height": 592} id="XcTBxQm6y9f0" outputId="8e730bcf-79cc-4f23-93bb-74d21b0aff8e"
y.plot(figsize=(20, 10))
plt.show()

# + colab={"base_uri": "https://localhost:8080/"} id="25Dw5F2a-d2L" outputId="36c0cdc6-5979-4a2c-cae1-e29d4611e616"
print("p-value:", adfuller(y.dropna())[1])

# + colab={"base_uri": "https://localhost:8080/"} id="D0gdgFY92fU0" outputId="c904891f-5c1d-4acb-ccce-2797470cea9b"
print("p-value:", adfuller(y.diff().dropna())[1])

# + colab={"base_uri": "https://localhost:8080/", "height": 592} id="6Vq5_rxa2sjq" outputId="88784092-dcd4-4028-a2e6-5adaf423c772"
y.diff().plot(figsize=(20, 10))
plt.show()

# + [markdown] id="6A5RjeckznKJ"
# ## Grid search

# + [markdown] id="tZHwdzyvzyG-"
# Since we have little information about the data generating process that generates the data we are going to study, we can use a grid search to determine the set of parameters that best explains the data using a simple metric: the Akaike Information Criteria. 
#
# The AIC is a tool we use to measure the performance of a model. By itself it does not really provide much information, but it shines when it is used to compare two or more models. 

# + colab={"base_uri": "https://localhost:8080/", "height": 409} id="keyC6Zei09qC" outputId="06ae9d03-ac65-4940-a133-a17f708f11b8"
# Trying out different models
m = ARIMA(y, order = (0, 1, 1))
m.fit().summary()


# + colab={"base_uri": "https://localhost:8080/", "height": 430} id="CJBRbiY03yPd" outputId="f83311c1-8548-4ff0-a724-6347f7523d7a"
m = ARIMA(y.diff()[1:], order = (0, 0, 1))
m.fit().summary()

# + id="0eeVkBj5AC8L"
p = d = q = range(0, 3)
pdq = list(itertools.product(p, d, q))


# + colab={"base_uri": "https://localhost:8080/"} id="GbDxESJ1AFNG" outputId="539ca014-5e7e-4101-e019-0490d963fdd7"
pdq[1]

# + colab={"base_uri": "https://localhost:8080/"} id="MhnNalR004ir" outputId="c15747e3-f34d-47ea-9abe-5c177ea857be"
for param in pdq:
        try:
            mod = ARIMA(y, order = param)

            results = mod.fit()

            print('ARIMA{}'.format(param), results.aic)
        except:
            continue

# + [markdown] id="k2Ojk0gqEKu_"
# According to this calculation, ARIMA(2,1,1) is the best performing model in terms of AIC out of the set of possible parameters analysed. 

# + colab={"base_uri": "https://localhost:8080/"} id="DK40QhiiEJXK" outputId="d3f9c665-255e-4783-d3ca-47866ce377d1"

chosenarimamodel = ARIMA(y, order = (2, 1, 1))
print(chosenarimamodel.fit().summary().tables[1])

# + id="XIvw2Zw5zdJa"
trainsize = int(len(y) * 0.90)
train, test = y[0:trainsize], y[trainsize:]
history = [x for x in train]
predictions = list()

# + colab={"base_uri": "https://localhost:8080/"} id="z_plOc7N_F_c" outputId="cd41b576-f137-420a-e342-1b9f0cda601c"
for t in range(len(test)):
	model = ARIMA(history, order = (2,1,1))
	model_fit = model.fit()
	yhat = model_fit.forecast()[0]
	predictions.append(yhat)
	history.append(test[t])
	# calculate out of sample error
error = mean_squared_error(test, predictions)
print(error)

# + colab={"base_uri": "https://localhost:8080/", "height": 283} id="d5CTCG5l_hEd" outputId="262cf8fb-1f22-4f01-8d93-083c018bbf2b"
train = train[1:]
start=len(train)
end=len(train)+len(test)-1

testarima = ARIMA(train, order = (2, 1, 1))

pred=testarima.fit().predict(start=start,end=end).rename('ARIMA Predictions')
pred.plot(legend=True)
test.plot(legend=True)


# + [markdown] id="HwYwAu58NWIJ"
# The forecast is not great, primarily because our model does not take into consideration the clearly seasonal model of the data. 

# + [markdown] id="zimQKkOwDt2X"
# We can expand this experiment to using the `SARIMAX` function, which allows for a seasonal component to be added to the modelling. With long and clearly seasonal datasets like this one, it might generate better results. 

# + [markdown] id="Dfmu0QJqGDDb"
# ## Seasonal ARIMA 

# + [markdown] id="-pMX4dTuGG_F"
# We can include a seasonal component to the estimation of the time series to incorporate the typical monthly behaviour not included in traditional ARIMA analysis. 
#
# SARIMAX stands for _Seasonal ARIMA with Exogenous variables_. In this case we did not add anu exogenous variables to the calculations, but it's relatively straightforward to do (see Example 4 of the documentation).  

# + id="YrLcUhWVzcFc"
p = d = q = range(0, 2)

pdq = list(itertools.product(p, d, q))
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]


# + colab={"base_uri": "https://localhost:8080/"} id="iWpIOFS9FIV9" outputId="3da76cf7-2d75-4cc1-ecc8-49b5b6e1a940"
warnings.filterwarnings("ignore") # specify to ignore warning messages

for param in pdq:
    for param_seasonal in seasonal_pdq:
        try:
            mod = sm.tsa.statespace.SARIMAX(y,
                                            order=param,
                                            seasonal_order=param_seasonal,
                                            enforce_stationarity=False,
                                            enforce_invertibility=False)

            results = mod.fit()

            print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
        except:
            continue

# + colab={"base_uri": "https://localhost:8080/"} id="mplMS-4EIEiU" outputId="b55ef64e-c57d-4fbf-f68b-9aff6b172128"
mod = sm.tsa.statespace.SARIMAX(y,
                                order=(1, 1, 1),
                                seasonal_order=(1, 1, 1, 12),
                                enforce_stationarity=False,
                                enforce_invertibility=False)

results = mod.fit()

print(results.summary().tables[1])

# + colab={"base_uri": "https://localhost:8080/", "height": 610} id="2p0ak5GrIzQr" outputId="225fd6cf-c23e-4761-82a1-25912580e8d7"
train = train[1:]
start=len(train)
end=len(train)+len(test)-1

testarima = sm.tsa.statespace.SARIMAX(y,
                                order=(1, 1, 1),
                                seasonal_order=(1, 1, 1, 12),
                                enforce_stationarity=False,
                                enforce_invertibility=False)

pred=testarima.fit().predict(start=start,end=end).rename('ARIMA Predictions')
pred.plot(legend=True, figsize=(20, 10))
test.plot(legend=True)

# + [markdown] id="3VhUVTiUGlXE"
# ## Some notes
#
# - The `auto_arima` function allows for automatic grid search. Be careful of the parameters used for the estimation, since controls need to be fined tuned to replicate the results we have above. 
# - We determined the best model using the AIC. There are other target functions you can use to choose your $p$, $d$, and $q$ parameters.  One possibility is the _mean square error_ between the predicted values given a set of parameters, and the test set. <NAME> has excellent tutorials on how to do this. 
# - An alternative way to check whether we have the right set of parameters is to check the ACF of the errors.
# - We must not forget to check whether the chosen parameters to explain our model fit our understanding of the underlying data. 

# + id="01ebmX4FGlJx"


# + [markdown] id="jlmyA36BMpoP"
# ## Sources

# + id="TuE7DYfs1nVM"
# https://www.digitalocean.com/community/tutorials/a-guide-to-time-series-forecasting-with-arima-in-python-3
# https://machinelearningmastery.com/make-sample-forecasts-arima-python/
# https://machinelearningmastery.com/grid-search-arima-hyperparameters-with-python/
# https://www.pluralsight.com/guides/advanced-time-series-modeling-(arima)-models-in-python
# https://people.duke.edu/~rnau/411arim.htm
# https://medium.com/swlh/temperature-forecasting-with-arima-model-in-python-427b2d3bcb53
# https://www.statsmodels.org/stable/index.html

StarcoderdataJupyterScriptsDedupFiltered