<jupyter_start><jupyter_text>ARIMA ARIMA, acronym of ‘Auto Regressive Integrated Moving Average’ is actually a class of models that ‘explains’ a given time series based on its own past values, that is, its own lags and the lagged forecast errors. ARIMA models aim to describe the autocorrelations in the data.  Data transformationIf the data show variation that increases or decreases with the level of the series, then a transformation can be useful. A useful family of transformations, that includes both logarithms and power transformations, is the family of Box-Cox transformations, which depend on the parameter $\lambda$ and are defined as follows:* $\lambda=0, w_t=log( y_t)$* $\lambda\neq0, w_t=(y^\lambda_t-1)/\lambda$A good value of $\lambda$ is one which makes the size of the seasonal variation about the same across the whole series, as that makes the forecasting model simpler.We will find the best $\lambda$ to transform Australian monthly electricity production:<jupyter_code>!pip install tscv
!pip install pmdarima
from datetime import datetime
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import tscv
from rpy2.robjects.packages import importr
import rpy2.robjects as ro
from statsmodels.tsa.holtwinters import SimpleExpSmoothing, Holt, ExponentialSmoothing
from statsmodels.tsa.stattools import kpss
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.graphics.tsaplots import plot_acf
from pmdarima.arima import auto_arima
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.arima_model import ARIMA
from tscv import GapWalkForward
from rpy2.robjects import pandas2ri, r
from scipy.stats import boxcox
pandas2ri.activate()
# install fpp2 package
utils = importr('utils')
utils.install_packages('fpp2')
# get the data and create dataframe
ro.r('library(fpp2)')
elec_np = r('elec')
elec_df = pd.DataFrame(elec_np, index=pd.date_range(start='1/1/1956',  periods=476, freq='M'), columns=['production'])
# plot data
elec_df.plot(figsize=(10,8))
plt.xlabel('year')
plt.ylabel('Australian monthly electricity production')
plt.show()<jupyter_output><empty_output><jupyter_text>The data show variation that increases with the years. We use the function **boxcox** to transform the data:<jupyter_code># estimate lambda and transform the data
elec_transf = boxcox(elec_df.production)
elec_df['boxcox'] = elec_transf[0]

# plot data
elec_df.boxcox.plot(figsize=(10,8))
plt.xlabel('year')
plt.ylabel('boxcox(elect)')
plt.show()

print('Lambda value is: ', elec_transf[1])<jupyter_output><empty_output><jupyter_text>Non-Seasonal ARIMA models ARIMA models combine differencing, autoregresion and moving average model. The full model can be written as: $y^{'}_{t}=c+\phi_{1} y^{'}_{t-1}+...+\phi_{p} y^{'}_{t-p}+\theta_{1} \epsilon_{t-1} +...+ \theta_{q} \epsilon_{t-q}+\epsilon_{t}$ Where $y^{'}_{t}$ is the differenced series. This model is called ARIMA($p,d,q$) model where: * $p$ = order of the autoregressive part * $d$ = degree of first differencing involved * $q$ = order of the moving average part We will use the quarterly percentage changes in US consumption expenditure to fit a non-seasonal ARIMA.<jupyter_code># get the data and create dataframe

us_np = r('uschange[,"Consumption"]')
index= pd.date_range(start='1970', end='2017', freq='Q')[0:-1]
us_df = pd.DataFrame(us_np, index)
us_df.columns = ['consumption']
# plot data
us_df.plot(figsize=(10,8))
plt.xlabel('quarter')
plt.ylabel('Quarterly percentange change')
plt.show()<jupyter_output><empty_output><jupyter_text>Pmdarima package brings R's auto.arima functionality to Python by wrapping statsmodel ARIMA and SARIMAX models. We will search the best model by minimizing the Akaike Information Criterion. We define the search range of the parameters $p$ and $q$.<jupyter_code># fitting a stepwise model:
stepwise_fit = auto_arima(us_df.consumption, start_p=1, start_q=1, max_p=3, max_q=3, m=4,
                             seasonal=False, trace=True,
                             error_action='ignore',  # don't want to know if an order does not work
                             suppress_warnings=True,  # don't want convergence warnings
                             information_criterion='aicc', # optimize aicc
                             stepwise=True)  # set to stepwise

stepwise_fit.summary()<jupyter_output>Performing stepwise search to minimize aicc
Fit ARIMA(1,0,1)x(0,0,0,0) [intercept=True]; AIC=343.815, BIC=356.740, Time=0.127 seconds
Fit ARIMA(0,0,0)x(0,0,0,0) [intercept=True]; AIC=375.632, BIC=382.095, Time=0.042 seconds
Fit ARIMA(1,0,0)x(0,0,0,0) [intercept=True]; AIC=353.332, BIC=363.025, Time=0.052 seconds
Fit ARIMA(0,0,1)x(0,0,0,0) [intercept=True]; AIC=360.371, BIC=370.064, Time=0.126 seconds
Fit ARIMA(0,0,0)x(0,0,0,0) [intercept=False]; AIC=529.744, BIC=532.975, Time=0.023 seconds
Fit ARIMA(2,0,1)x(0,0,0,0) [intercept=True]; AIC=344.493, BIC=360.649, Time=0.163 seconds
Fit ARIMA(1,0,2)x(0,0,0,0) [intercept=True]; AIC=343.330, BIC=359.485, Time=0.139 seconds
Fit ARIMA(0,0,2)x(0,0,0,0) [intercept=True]; AIC=356.797, BIC=369.721, Time=0.080 seconds
Fit ARIMA(2,0,2)x(0,0,0,0) [intercept=True]; AIC=342.284, BIC=361.671, Time=0.255 seconds
Fit ARIMA(3,0,2)x(0,0,0,0) [intercept=True]; AIC=342.536, BIC=365.154, Time=0.299 seconds
Fit ARIMA(2,0,3)x(0,0,0,0) [intercept=True]; AIC=343.17[...]<jupyter_text>The best model found is ARIMA(1,0,3). Attending the parameters provided (**intercept**, **ar.L1**, **ma.L1**, **ma.L2**, **ma.L3**), the model can be expressed as:$y^{'}_{t}= 0.307 + 0.589 y^{'}_{t-1}-0.353 \epsilon_{t-1} +0.085\epsilon_{t-2} + 0.174 \epsilon_{t-3}+\epsilon_{t}$If we know in advance the parameters, we can call directly the function _ARIMA_ from the package **statsmodels**.<jupyter_code># fit ARIMA model

arima_model = ARIMA(us_df.consumption, order=(1,0,3))
arima_fit = arima_model.fit()
arima_fit.summary()<jupyter_output><empty_output><jupyter_text>We predict future values as usual.<jupyter_code># predicting the next 10 values

predict_arima = arima_fit.forecast(10)[0]
predict_arima

fig, ax = plt.subplots(figsize=(16, 8))
ax.plot(us_df.index, us_df.consumption.values, label='consumption')
ax.plot(pd.date_range(start='31/12/2016',  periods=10, freq='Q'), predict_arima, color="green", label='10 quarters ahead')
plt.ylabel('Quarterly percentange change') 
plt.xlabel('year') 
plt.legend();<jupyter_output><empty_output><jupyter_text>The function **auto_arima** does not consider all possible models in its search. We can make it work harder by using the argument _stepwise_=FALSE. However, the stepwise algorithm can be significantly faster than fitting all hyper-parameter combinations and is less likely to over-fit the model so it would be our first choice.<jupyter_code># fitting a grid search model:
grid_fit = auto_arima(us_df.consumption, start_p=1, start_q=1, max_p=3, max_q=3, m=4,
                             seasonal=False, trace=True,
                             error_action='ignore',  # don't want to know if an order does not work
                             suppress_warnings=True,  # don't want convergence warnings
                             information_criterion='aicc', # optimize aicc
                             stepwise=False)  # set to stepwise

grid_fit.summary()<jupyter_output>Fit ARIMA(0,0,0)x(0,0,0,0) [intercept=True]; AIC=375.632, BIC=382.095, Time=0.038 seconds
Fit ARIMA(0,0,1)x(0,0,0,0) [intercept=True]; AIC=360.371, BIC=370.064, Time=0.126 seconds
Fit ARIMA(0,0,2)x(0,0,0,0) [intercept=True]; AIC=356.797, BIC=369.721, Time=0.078 seconds
Fit ARIMA(0,0,3)x(0,0,0,0) [intercept=True]; AIC=342.758, BIC=358.914, Time=0.099 seconds
Fit ARIMA(1,0,0)x(0,0,0,0) [intercept=True]; AIC=353.332, BIC=363.025, Time=0.050 seconds
Fit ARIMA(1,0,1)x(0,0,0,0) [intercept=True]; AIC=343.815, BIC=356.740, Time=0.116 seconds
Fit ARIMA(1,0,2)x(0,0,0,0) [intercept=True]; AIC=343.330, BIC=359.485, Time=0.135 seconds
Fit ARIMA(1,0,3)x(0,0,0,0) [intercept=True]; AIC=341.613, BIC=361.000, Time=0.159 seconds
Fit ARIMA(2,0,0)x(0,0,0,0) [intercept=True]; AIC=346.268, BIC=359.192, Time=0.074 seconds
Fit ARIMA(2,0,1)x(0,0,0,0) [intercept=True]; AIC=344.493, BIC=360.649, Time=0.164 seconds
Fit ARIMA(2,0,2)x(0,0,0,0) [intercept=True]; AIC=342.284, BIC=361.671, Time=0.267 seconds
Fit ARIMA([...]<jupyter_text>The model found now ARIMA(3,0,0) is slightly better than the previous one, attending to the AIC.  Check the residualsIt is recommended to check if the residuals behave like white noise. If not the model is not considered good enough and we need to try a better adjustment.We will check the residuals from the model fitted in the last section.<jupyter_code># ACF plot

plot_acf(stepwise_fit.resid(), zero=False)
plt.show()<jupyter_output><empty_output><jupyter_text>The ACF plot shows that almost all autocorrelations are within the threshold limits except at 8 and 22 lags.Ljung-Box test helps to determine if the residuals can be considered as white noise.<jupyter_code>acorr_ljungbox(stepwise_fit.resid(), lags=22)[1]<jupyter_output>/usr/local/lib/python3.6/dist-packages/statsmodels/stats/diagnostic.py:524: FutureWarning: The value returned will change to a single DataFrame after 0.12 is released.  Set return_df to True to use to return a DataFrame now.  Set return_df to False to silence this warning.
  warnings.warn(msg, FutureWarning)<jupyter_text>Ljung-Box test returns large p-values for all lags tested, also suggesting that the residuals are white noise.  Seasonal ARIMA modelsARIMA models are also capable of modelling a wide range of seasonal data. A seasonal ARIMA model is formed by including additional seasonal terms in the ARIMA:ARIMA(p,d,q)(P,D,Q)mwhere m = number of observations per year. The seasonal part of the model consists of terms that are similar to the non-seasonal components of the model, but involve backshifts of the seasonal period.We will fit a seasonal ARIMA model using quarterly European retail trade data from 1996 to 2011.<jupyter_code># get the data and create dataframe

eu_np = r('euretail')
index= pd.date_range(start='1996', end='2012', freq='Q')
eu_df = pd.DataFrame(eu_np, index)
eu_df.columns = ['retail']
# plot data

eu_df.plot(figsize=(10,8))
plt.xlabel('year')
plt.ylabel('Quarterly retail index')
plt.show()<jupyter_output><empty_output><jupyter_text>Now, we will find the best ARIMA model.<jupyter_code># fitting a grid search model:
eu_fit_p = auto_arima(eu_df.retail, start_p=0, start_q=0, max_p=3, max_q=3, m=4,
                          seasonal=True, trace=True, start_P=0, max_P=3,
                          start_Q=0, max_Q=3,
                          information_criterion='aicc', # optimize by aicc
                          error_action='ignore',  # don't want to know if an order does not work
                          suppress_warnings=True,  # don't want convergence warnings
                          stepwise=False)  # set to stepwise

eu_fit_p.summary()<jupyter_output>Fit ARIMA(0,2,0)x(0,0,0,4) [intercept=True]; AIC=127.888, BIC=132.142, Time=0.040 seconds
Fit ARIMA(0,2,0)x(0,0,1,4) [intercept=True]; AIC=122.289, BIC=128.670, Time=0.044 seconds
Fit ARIMA(0,2,0)x(0,0,2,4) [intercept=True]; AIC=108.684, BIC=117.192, Time=0.082 seconds
Fit ARIMA(0,2,0)x(0,0,3,4) [intercept=True]; AIC=110.010, BIC=120.646, Time=0.239 seconds
Fit ARIMA(0,2,0)x(1,0,0,4) [intercept=True]; AIC=113.800, BIC=120.181, Time=0.051 seconds
Fit ARIMA(0,2,0)x(1,0,1,4) [intercept=True]; AIC=96.691, BIC=105.200, Time=0.288 seconds
Fit ARIMA(0,2,0)x(1,0,2,4) [intercept=True]; AIC=98.605, BIC=109.241, Time=0.488 seconds
Fit ARIMA(0,2,0)x(1,0,3,4) [intercept=True]; AIC=98.599, BIC=111.362, Time=0.746 seconds
Fit ARIMA(0,2,0)x(2,0,0,4) [intercept=True]; AIC=101.061, BIC=109.569, Time=0.108 seconds
Fit ARIMA(0,2,0)x(2,0,1,4) [intercept=True]; AIC=98.576, BIC=109.211, Time=0.621 seconds
Fit ARIMA(0,2,0)x(2,0,2,4) [intercept=True]; AIC=97.829, BIC=110.592, Time=0.617 seconds
Fit ARIMA(0,2,0[...]<jupyter_text>The model found ARIMA(1,2,1)(1,0,1)4 differs from the model reported in https://otexts.com/fpp2/seasonal-arima.html (ARIMA(0,1,3)(0,1,1)4). The parameters of differencing used in Python, during the search, are: d=2 and D=0. The best model uses d=0 and D=1 instead. We call the R code to get the best combination reported for this problem.<jupyter_code>ro.r('library(forecast)')
ro.r('print(auto.arima(euretail))')<jupyter_output>Series: euretail 
ARIMA(0,1,3)(0,1,1)[4] 

Coefficients:
         ma1     ma2     ma3     sma1
      0.2630  0.3694  0.4200  -0.6636
s.e.  0.1237  0.1255  0.1294   0.1545

sigma^2 estimated as 0.156:  log likelihood=-28.63
AIC=67.26   AICc=68.39   BIC=77.65<jupyter_text>This new model shows a best value of AIC=67.26 in compare with the AIC=78.30 obtained previously. Now we fit this model using the Python SARIMAX function from the package **statsmodels**:<jupyter_code># create the seasonal arima model

arima_model = SARIMAX(eu_df.retail, order=(0,1,3), seasonal_order=(0,1,1,4))
eu_fit_r = arima_model.fit()
eu_fit_r.summary()<jupyter_output><empty_output><jupyter_text>We will predict the next 12 values.<jupyter_code># predicting the next 12 values

predict_arima = eu_fit_r.forecast(steps=12)
predict_arima

fig, ax = plt.subplots(figsize=(16, 8))
ax.plot(eu_df.index, eu_df.retail.values, label='retail')
ax.plot(pd.date_range(start='31/1/2012',  periods=12, freq='Q'), predict_arima, color="green", label='12 quarters ahead')
plt.ylabel('Quarterly retail index') 
plt.xlabel('year') 
plt.legend();<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13309750165819145, 'fraction_numerical': 0.07561353084236126, 'mean_word_length': 5.1043634727845255, 'pattern_counts': {'":': 0, '<': 52, '<?xml version=': 0, '>': 52, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '579845', 'license': 'Apache-2.0', 'n_tokens_mistral': 5547, 'n_tokens_neox': 5012, 'n_words': 1310, 'path': 'notebooks/arima.ipynb', 'repo_name': 'mabelc/forcasting', 'size': '13569'}