{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"**If you think this notebook deserves an upvote, I'd love to have it. An upvote per view, its all I ask**\n(credit to [Dan Carlin](https://twitter.com/HardcoreHistory) for coining the phrase ;-) \n\n---------------------------------------\n\nThis is part of a series of notebooks about practical time series methods:\n\n* [Part 0: the basics](https://www.kaggle.com/konradb/ts-0-the-basics)\n* [Part 1a: smoothing methods](https://www.kaggle.com/konradb/ts-1a-smoothing-methods)\n* [Part 1b: Prophet](https://www.kaggle.com/konradb/ts-1b-prophet) \n* [Part 2: ARMA](https://www.kaggle.com/konradb/ts-2-arma-and-friends) - **this notebook**\n* [Part 3: Time series for finance](https://www.kaggle.com/konradb/ts-3-time-series-for-finance) \n* [Part 4: Sales and demand forecasting](https://www.kaggle.com/konradb/ts-4-sales-and-demand-forecasting)\n* [Part 5: Automatic for the people](https://www.kaggle.com/code/konradb/ts-5-automatic-for-the-people)\n* [Part 6: Deep learning for TS - sequences](https://www.kaggle.com/konradb/ts-6-deep-learning-for-ts-sequences) \n* [Part 7: Survival analysis](https://www.kaggle.com/konradb/ts-7-survival-analysis)\n* [Part 8: Hierarchical time series](https://www.kaggle.com/code/konradb/ts-8-hierarchical-time-series)\n\n\nThe series is accompanied by video presentations on the YouTube channel of [Abhishek](https://www.kaggle.com/abhishek):\n\n* [Talk 0](https://www.youtube.com/watch?v=cKzXOOtOXYY) \n* [Talk 1](https://www.youtube.com/watch?v=kAI67Sz92-s) - combining the content from parts 1a and 1b\n* [Talk 2](https://www.youtube.com/watch?v=LjV5DE3KR-U) - **based on this notebook**\n* [Talk 3](https://www.youtube.com/watch?v=74rDhJexmTg) \n* [Talk 4](https://www.youtube.com/watch?v=RdH8zd07u2E)  \n* [Talk 5](https://www.youtube.com/watch?v=wBP8Pc4Wxzs)\n* [Talk 6](https://www.youtube.com/watch?v=81AEI0tj0Kk)\n* [Talk 7](https://www.youtube.com/watch?v=m-8I_hkmz9o)\n* [Talk 8]()\n\n\n---------------------------------------\n\n\nThe notebook is split into three sections: we introduce the basic framework of linear processes, then present extensions and finally demonstrate how to solve a prediction problem from scratch. \n\n* [Basic linear processes](#section-one)\n* [Beyond ARMA](#section-two)\n* [Full pipeline](#section-three)\n* [A day at the races](#section-four)\n\nAs usual, we begin by importing the required packages","metadata":{}},{"cell_type":"code","source":"# It's a Surprise Tool That Will Help Us Later ;-)\n!pip install pmdarima","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:34.732382Z","iopub.execute_input":"2022-02-14T01:10:34.732943Z","iopub.status.idle":"2022-02-14T01:10:44.510691Z","shell.execute_reply.started":"2022-02-14T01:10:34.732831Z","shell.execute_reply":"2022-02-14T01:10:44.509219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import itertools\nimport pandas as pd\nimport numpy as np\nfrom random import gauss\n\nimport statsmodels.api as sm\nfrom statsmodels.graphics.tsaplots import plot_acf, plot_pacf\nfrom statsmodels.tsa.ar_model import AR\nfrom statsmodels.tsa.seasonal import seasonal_decompose\nfrom statsmodels.tsa.stattools import adfuller\nfrom statsmodels.tsa.arima_process import ArmaProcess\n\nfrom pandas.plotting import autocorrelation_plot\nfrom sklearn.metrics import mean_squared_error\n\nimport os\nfrom IPython.display import Image\n\nimport matplotlib.pyplot as plt\nimport warnings\nimport itertools\nimport statsmodels.api as sm\nfrom random import gauss\nfrom pandas.plotting import autocorrelation_plot\nfrom statsmodels.graphics.tsaplots import plot_acf, plot_pacf\n\nfrom pmdarima.arima import auto_arima\nfrom pmdarima import pipeline\nfrom pmdarima import model_selection\nfrom pmdarima import preprocessing as ppc\nfrom pmdarima import arima\n\nfrom fbprophet import Prophet \n\nimport matplotlib.pyplot as plt\nplt.style.use('fivethirtyeight') \n\nwarnings.simplefilter(action='ignore', category= FutureWarning)","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:44.513673Z","iopub.execute_input":"2022-02-14T01:10:44.514187Z","iopub.status.idle":"2022-02-14T01:10:46.832949Z","shell.execute_reply.started":"2022-02-14T01:10:44.514126Z","shell.execute_reply":"2022-02-14T01:10:46.831099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# general settings\nclass CFG:\n    data_folder = '../input/tsdata-1/'\n    img_dim1 = 20\n    img_dim2 = 10\n        \n# adjust the parameters for displayed figures    \nplt.rcParams.update({'figure.figsize': (CFG.img_dim1,CFG.img_dim2)})    ","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:46.835596Z","iopub.execute_input":"2022-02-14T01:10:46.835991Z","iopub.status.idle":"2022-02-14T01:10:46.843555Z","shell.execute_reply.started":"2022-02-14T01:10:46.835956Z","shell.execute_reply":"2022-02-14T01:10:46.842073Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"section-one\"></a>\n# Basic linear processes\n\nUntil 1980s, ARIMA processes were the workhorses of time series analysis, which for many people meant just finding right order of such model. This class of processes is defined through linear relationship between observation and noise factors.\n\nJust to indulge my inner nerd, I will go top-down in this one and start from a high level definition - because why not :-) A linear process is a time series $X_t$ defined by:\n\n\\begin{equation}\n    X_t = \\sum_{u = -\\infty}^{\\infty} \\psi_u \\epsilon_{t-u}\n\\end{equation}\n\nwhere $\\epsilon_t$ is a white noise series and \n\n\\begin{equation}\n\\sum_{u = -\\infty}^{\\infty} \\left|\\psi_u \\right|^2 < \\infty\n\\end{equation}\n\n\nTranslated to plain English:\n* (in our context) white noise a serially uncorrelated sequence of random variables with zero mean and finite variance. For a general definition (and the origin of the color convention) I can recommend the Wikipedia entry: https://en.wikipedia.org/wiki/White_noise\n\n* it is a linear combination - potentially infinite and depending on both past and present - of a white noise series\n\n* in the real world it has to be causal $\\implies$ depending on a finite number of past values\n\n* it is a unified framework for handling different types of data generating processes\n\n* for the mathematically inclided person reading this: huge parts of the theory behind ARMA processes pop up as special cases / corrolaries from results in functional analysis, specifically the [Hilbert projection theorem](https://en.wikipedia.org/wiki/Hilbert_projection_theorem) and [spectral theorem](https://en.wikipedia.org/wiki/Spectral_theorem); the requirement on sum of squared values being finite correspond to belonging to an $l_2$ space for our sequence: https://en.wikipedia.org/wiki/Lp_space\n\n* only needs second order statistics (variance) $\\rightarrow$ works best for elliptical distributions (e.g. Gaussian)\n","metadata":{}},{"cell_type":"markdown","source":"\n## AR processes","metadata":{}},{"cell_type":"markdown","source":"The simplest (non-trivial) example of a linear process is an autoregressive process of order $p$:\n\n\\begin{equation}\n    X_t = \\sum_{i= 1}^p \\phi_i X_{t-i} + \\epsilon_t\n\\end{equation}\n    \nUnpacking the formula:\n1. Natural extension of a multiple linear regression model: we forecast the variable of interest using a linear combination of past values of the variable (lagged values act as predictors). The term \\textit{auto}regression reflects the fact that we are regressing the variable on (version of) itself\n2. The single exponential smoothing method (described [here](https://www.kaggle.com/konradb/ts-1a-smoothing-methods) ) can be viewed as a special case of an autoregressive process of order 1\n3. Easy to estimate parameters and forecast \n\n\nNow that we have defined (and hopefully understood) what an AR process is, let's have a look at how to identify one. This is relevant if we need to decide whether an autoregressive process is a right kind of model for a particular dataset. For the sake of demonstration, we will use simulated data, where the data generating process is:\n\n\\begin{equation}\n    X_t = 0.9 X_{t-1} + \\epsilon_t\n\\end{equation}\n    \nObserve in the codeblock below that the coefficient changes sign - this is because of the notational convention in statsmodels, where we read the coefficients in order defined by the characteristic polynomial (more on that below):\n\n\\begin{equation}\n    X_t - 0.9 X_{t-1} = \\epsilon_t\n\\end{equation}\n","metadata":{}},{"cell_type":"markdown","source":"A fast technique for deciding if AR is the right kind of model for a given series is to examine the autocorrelation and partial autocorrelation functions (described [here](https://www.kaggle.com/konradb/practical-time-series-pt-1-the-basics#Dependence)). ","metadata":{}},{"cell_type":"code","source":"ar1 = np.array([1, -0.9])\nma1 = np.array([1])\nAR_object1 = ArmaProcess(ar1, ma1)\nsimulated_data_1 = AR_object1.generate_sample(nsample=1000)\nplt.plot(simulated_data_1)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:46.845425Z","iopub.execute_input":"2022-02-14T01:10:46.846009Z","iopub.status.idle":"2022-02-14T01:10:47.106894Z","shell.execute_reply.started":"2022-02-14T01:10:46.845966Z","shell.execute_reply":"2022-02-14T01:10:47.105829Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_acf(simulated_data_1, lags = 10)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:47.108369Z","iopub.execute_input":"2022-02-14T01:10:47.108831Z","iopub.status.idle":"2022-02-14T01:10:47.364973Z","shell.execute_reply.started":"2022-02-14T01:10:47.108784Z","shell.execute_reply":"2022-02-14T01:10:47.363839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_pacf(simulated_data_1, lags = 10)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:47.366602Z","iopub.execute_input":"2022-02-14T01:10:47.366981Z","iopub.status.idle":"2022-02-14T01:10:47.604766Z","shell.execute_reply.started":"2022-02-14T01:10:47.366948Z","shell.execute_reply":"2022-02-14T01:10:47.60379Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Intuition for identification:\n* process memory: direct and indirect dependence information\n* ACF for AR(p): strong until lag p, trailing off afterwards\n* PACF only describes the direct relationship between an observation and its lag. This would suggest that there would be no correlation for lag values beyond k.\n* analytical derivation is based on the Laurent series expansion and the backshift operator $B$\n","metadata":{}},{"cell_type":"markdown","source":"## MA processes","metadata":{}},{"cell_type":"markdown","source":"An autoregressive model expresses the forecast variable as a linear combination of past realizations of itself and the same idea can be applied to past forecast errors: a moving average process of order $q$ is defined by the relationship:\n    \n\\begin{equation}\nX_t = \\epsilon_t + \\sum_{i = 1}^q \\theta_i \\epsilon_{t-i}\n\\end{equation}\n\nwhere $\\epsilon_t$ is a white noise series.\n  \nLet's repeat the identification exercise from before:\n\n\\begin{equation}\n    X_t = \\epsilon_t + 0.9  \\epsilon_{t-1}\n\\end{equation}\n","metadata":{}},{"cell_type":"code","source":"ar1 = np.array([1])\nma1 = np.array([1, -0.9])\nMA_object1 = ArmaProcess(ar1, ma1)\nsimulated_data_1 = MA_object1.generate_sample(nsample=1000)\nplt.plot(simulated_data_1)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:47.606034Z","iopub.execute_input":"2022-02-14T01:10:47.606391Z","iopub.status.idle":"2022-02-14T01:10:47.828031Z","shell.execute_reply.started":"2022-02-14T01:10:47.606346Z","shell.execute_reply":"2022-02-14T01:10:47.826712Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_acf(simulated_data_1, lags = 10)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:47.832406Z","iopub.execute_input":"2022-02-14T01:10:47.832981Z","iopub.status.idle":"2022-02-14T01:10:48.086162Z","shell.execute_reply.started":"2022-02-14T01:10:47.832929Z","shell.execute_reply":"2022-02-14T01:10:48.085381Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_pacf(simulated_data_1, lags = 10)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:48.087841Z","iopub.execute_input":"2022-02-14T01:10:48.088322Z","iopub.status.idle":"2022-02-14T01:10:48.336275Z","shell.execute_reply.started":"2022-02-14T01:10:48.088273Z","shell.execute_reply":"2022-02-14T01:10:48.334882Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Intuition for identification:\n* MA model corrects future forecasts based on errors made on recent forecasts.\n* ACF for MA(k) series: to show a strong correlation with recent values up to the lag of k, then a sharp decline \n* PACF: strong relationship to the lag k and a trailing off of correlation from the lag onwards.","metadata":{}},{"cell_type":"markdown","source":"## ARMA","metadata":{}},{"cell_type":"markdown","source":"We have an autoregressive component and a moving average one, so it is quite natural to combine those two types of dynamics into a single model: ARMA(p,q) series satisfies the relationship:\n\n\\begin{equation}\nX_t = \\sum_{i=1}^p \\phi_i X_{t-i} + \\sum_{j=1}^q \\theta_j \\epsilon_{t-j} + \\epsilon_t\n\\end{equation}\n\nwhere our predictors on the right hand side include both lagged values of the series and lagged errors, $p$ is the order the autoregressive part and $q$ is the order of the moving average component.\n\nSome observations:\n1. if an ARMA(p,q) model is stationary (see [this one](https://www.kaggle.com/konradb/ts-0-the-basics#Stationarity) for a refresher), it can be represented as an infinite AR series\n\n\\begin{equation}    \n    X_t = \\sum_{u=1}^\\infty \\pi_u X_{t-u} + \\epsilon_t\n\\end{equation}\n\n2. this allows for estimation via maximum likelihood and a simple recursive forecast:\n\n\\begin{equation}\n    \\hat{X}_{T+1} = \\sum_{u=1}^\\infty \\hat{\\pi}_u X_{T+1 - u}\n\\end{equation}\n\nIn practice: Kalman filter.\n\n3. Stationarity of the ARMA(p,q) process is established by analyzing the characteristic polynomial. There is quite a substantial body of theory behind this idea, but we are focusing on practical applications here - so a crash-course argument would go like this:\n\n    * given ARMA(p,q) series defined by the equation above, we can write its characteristic polynomial \n    \n    \\begin{equation}\n    P(z) = 1 - \\phi_1 z - \\ldots - \\phi_p z^p\n    \\end{equation}\n    \n    * if we look for solutions in complex (as opposed to real numbers) domain, the equation P(z) = 0 has $p$ solutions $z_1, \\ldots, z_p$ (corollary from the Fundamental Theorem of Algebra: https://en.wikipedia.org/wiki/Fundamental_theorem_of_algebra)\n       \n    * if $|z_i| >1$ for all $i$ then the underlying model is stationary\n    * tests like Augmented Dickey-Fuller check for this behavior\n    ","metadata":{}},{"cell_type":"code","source":"\nImage(filename=\"../input/muh-images/unit_roots.png\", width= 600, height=200)\n ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-02-14T01:10:48.339046Z","iopub.execute_input":"2022-02-14T01:10:48.339702Z","iopub.status.idle":"2022-02-14T01:10:48.359136Z","shell.execute_reply.started":"2022-02-14T01:10:48.339643Z","shell.execute_reply":"2022-02-14T01:10:48.358138Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Forecasting with ARMA\n\nLet's put our newly introduced framework to the test and actually predict something. We look we will have a look at the quarterly change in the aggregate savings level in the United States - the data is sourced from FRED: https://fred.stlouisfed.org/series/A191RP1Q027SBEA.","metadata":{}},{"cell_type":"code","source":"xdat = pd.read_csv(CFG.data_folder + 'savings_change.csv')\n\nxdat.columns = ['date', 'value']\nxdat['date'] = pd.to_datetime(xdat['date'])\n\nxdat.set_index('date').plot()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:48.360616Z","iopub.execute_input":"2022-02-14T01:10:48.361209Z","iopub.status.idle":"2022-02-14T01:10:48.832814Z","shell.execute_reply.started":"2022-02-14T01:10:48.361148Z","shell.execute_reply":"2022-02-14T01:10:48.831275Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We start with a quick EDA: first, a seasonal decomposition","metadata":{}},{"cell_type":"code","source":"decomposition = sm.tsa.seasonal_decompose(xdat[\"value\"],period =12, model = 'additive') \nfigure = decomposition.plot()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:48.839824Z","iopub.execute_input":"2022-02-14T01:10:48.840458Z","iopub.status.idle":"2022-02-14T01:10:49.531947Z","shell.execute_reply.started":"2022-02-14T01:10:48.840405Z","shell.execute_reply":"2022-02-14T01:10:49.530597Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There is *a* seasonal component, but it is very small in magnitude; trend does not seem significant either - we can suspect the process is stationary and move to testing this formally (using our old acquaintance, the ADF test):","metadata":{}},{"cell_type":"code","source":"result = adfuller(xdat['value'])\nprint('ADF Statistic: %f' % result[0])\nprint('p-value: %f' % result[1])\nprint('Critical Values:')\nfor key, value in result[4].items():\n    print('\\t%s: %.3f' % (key, value))","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:49.533475Z","iopub.execute_input":"2022-02-14T01:10:49.533808Z","iopub.status.idle":"2022-02-14T01:10:49.570279Z","shell.execute_reply.started":"2022-02-14T01:10:49.533768Z","shell.execute_reply":"2022-02-14T01:10:49.568517Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"No red flags stopping us at this stage (reminder: this means we don't see an overall linear trend during the sample period) - we can proceed to model the series as ARMA(p,q). What about the values of p and q?","metadata":{}},{"cell_type":"code","source":"plot_acf(xdat['value'], lags = 25); print()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:49.572565Z","iopub.execute_input":"2022-02-14T01:10:49.57313Z","iopub.status.idle":"2022-02-14T01:10:49.87256Z","shell.execute_reply.started":"2022-02-14T01:10:49.573077Z","shell.execute_reply":"2022-02-14T01:10:49.871793Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_pacf(xdat['value'], lags = 25); print()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:49.873698Z","iopub.execute_input":"2022-02-14T01:10:49.874136Z","iopub.status.idle":"2022-02-14T01:10:50.159419Z","shell.execute_reply.started":"2022-02-14T01:10:49.874099Z","shell.execute_reply":"2022-02-14T01:10:50.158322Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"At this stage, we could make an educated guess at fitting an ARMA model - (1,0) seems like a reasonable first pass. We can do something better though and use the functionality provided by the `pmdarima` package. Its functionality modernizes the ARIMA functionality from `statsmodels`:\n* it allows working with pipelines combining preprocessing and fitting\n* enables cross-validation functionality akin to what we saw with`cross_validation` in [Prophet](https://www.kaggle.com/konradb/ts-1b-prophet#Performance-evaluation)\n* automates a (somewhat tedious) task of establishing the right order for an ARIMA model","metadata":{}},{"cell_type":"code","source":"# we split the data into training and validation, leaving out the last 3 years as holdout\nxtrain, xvalid = model_selection.train_test_split(xdat['value'], test_size = 12)","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:50.16088Z","iopub.execute_input":"2022-02-14T01:10:50.161202Z","iopub.status.idle":"2022-02-14T01:10:50.16737Z","shell.execute_reply.started":"2022-02-14T01:10:50.161169Z","shell.execute_reply":"2022-02-14T01:10:50.166296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"One of the useful features of `pmdarima` is a Fourier featurizer:\n- it allows us to capture seasonality in the model without worrying about the seasonal order (so use a pure ARIMA with `seasonality = False`) \n- its achieved by creating a set of covariates based on Fourier decomposition\n- we can include arbitrary length seasonal pattern\n- the short-term dynamics are easily handled with a simple ARMA error","metadata":{}},{"cell_type":"code","source":"# construct a pipeline, combining feature preprocessing and the model to fit\npipe = pipeline.Pipeline([\n    (\"fourier\", ppc.FourierFeaturizer(m=4)),\n    (\"arima\", arima.AutoARIMA(stepwise=True, trace=1, error_action=\"ignore\",\n                              seasonal=False,  # because we use Fourier\n                              suppress_warnings=True))\n])\n\npipe.fit(xtrain)","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:50.168929Z","iopub.execute_input":"2022-02-14T01:10:50.169269Z","iopub.status.idle":"2022-02-14T01:10:52.633139Z","shell.execute_reply.started":"2022-02-14T01:10:50.169237Z","shell.execute_reply":"2022-02-14T01:10:52.632123Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# we can inspect the results\npipe.summary()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:52.634467Z","iopub.execute_input":"2022-02-14T01:10:52.634806Z","iopub.status.idle":"2022-02-14T01:10:52.655908Z","shell.execute_reply.started":"2022-02-14T01:10:52.634765Z","shell.execute_reply":"2022-02-14T01:10:52.655151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Statsmodels provides us with a rich set of descriptive measures characterizing the fitted model:\n* the top table summarizes sample size, likelihood and the information criteria - this is helpful if we want to select models on parsimony. For those unfamiliar with the topic, information criteria (AIC, BIC, HQIC ) quantify the basic idea of balancing model complexity (number of parameters) and likelihood, so that we end up selecting a model which has the minimal number of parameters necessary to capture the details of our DGP, but not more than that:https://en.wikipedia.org/wiki/Akaike_information_criterion. An information criterion measures how well a model fits the data while taking into account the overall complexity - model that fits the data very well while using large number of parameters will be assigned a larger AIC score than a model that uses fewer features to achieve the same goodness-of-fit. \n\n* the middle table gives the estimated coefficients of the model: $ar.L_{i}$ corresponds to the coefficient of the AR part at lag $i$, which $sigma2$ is the variance of the noise component. The $Z$ statistic and the endpoints of a confidence interval allow for a quick assessment of statistical significance\n\n* finally, the box at the bottom summarizes some diagnostic test: Jarque-Bera is used for testing if the residuals of the model have Gaussian distribution, while Ljung-Box checks whether the results are serially independent. In our case the former shows no problem, while L-B indicates there is some serial dependence in the residuals of our model, so there mighte be some component of the dynamics we are not capturing. ","metadata":{}},{"cell_type":"code","source":"# we compute predictions from the pipeline object\npreds, conf_int = pipe.predict(n_periods= xvalid.shape[0], return_conf_int=True)\nprint(\"\\nForecasts:\")\nprint(preds)\n","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:52.657178Z","iopub.execute_input":"2022-02-14T01:10:52.657455Z","iopub.status.idle":"2022-02-14T01:10:52.676427Z","shell.execute_reply.started":"2022-02-14T01:10:52.657426Z","shell.execute_reply":"2022-02-14T01:10:52.675678Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"\\nConfidence intervals:\")\nprint(conf_int)\n","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:52.677908Z","iopub.execute_input":"2022-02-14T01:10:52.678187Z","iopub.status.idle":"2022-02-14T01:10:52.684056Z","shell.execute_reply.started":"2022-02-14T01:10:52.678158Z","shell.execute_reply":"2022-02-14T01:10:52.683256Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Unlike the uncertainty bands for Prophet (which are calculated non-parametrically by resampling the data with MCMC), ARMA confidence intervals are based on a parametric approach:\n* we assume the residual terms $\\epsilon_t$ have Gaussian distribution\n* Gaussian is closed under addition $\\implies$ $X_t$ is Gaussian\n* Hilbert projection theorem $\\implies$ expected value is our best forecast wrt MSE\n* confidence interval for the mean: $\\mu_t \\pm 1.96 \\; \\sigma_t$ with\n\\begin{equation}\n\\mu_t = E\\left[X_t | \\sigma(X_t) \\right]\n\\end{equation}\nand\n\\begin{equation}\n\\sigma^2_t = Var\\left[X_t | \\sigma(X_t) \\right]\n\\end{equation}","metadata":{}},{"cell_type":"code","source":"xvalid","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:52.685403Z","iopub.execute_input":"2022-02-14T01:10:52.685892Z","iopub.status.idle":"2022-02-14T01:10:52.699227Z","shell.execute_reply.started":"2022-02-14T01:10:52.68586Z","shell.execute_reply":"2022-02-14T01:10:52.698232Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# let's quickly visualize the forecast\nxvalid = pd.DataFrame(xvalid.values, columns = ['actual'])\nxvalid['predicted'] = preds\nxvalid.plot()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:52.700685Z","iopub.execute_input":"2022-02-14T01:10:52.701251Z","iopub.status.idle":"2022-02-14T01:10:52.953396Z","shell.execute_reply.started":"2022-02-14T01:10:52.701219Z","shell.execute_reply":"2022-02-14T01:10:52.952313Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In terms of forecast quality this model can be useful as a baseline, but overall it is clearly nothing to write home about:\n1. there is a predictable lag: one period after an increase in the original series the forecast goes up as well; this is due to autoregressive nature of the model\n2. the model captures the general dynamics of the model, but struggles with the range of values - this is a consequence of the constant variance assumption","metadata":{}},{"cell_type":"markdown","source":"<a id=\"section-two\"></a>\n# Beyond ARMA\n\nAs we can see from the examples above, there are issues when applying the baseline $ARMA(p,q)$ model to real life time series:\n* the models are only efficient if the series is stationary - peformance deteriorates rapidly otherwise \n* does not allow for trend\n* does not allow for seasonality\n\nThe last two issues can be handled by differencing the series at appropriate lags, and it is this idea that allows us to extend ARMA - first, we move towards ARIMA. \n\n### ARIMA\n\nA process $X_t$ is ARIMA(p,d,q) $\\iff$ $\\nabla^d X_t$ is a stationary ARMA(p,q)\n\nUnpacking the (slightly cryptic :-) formulation: \n* a series becomes a stationary ARMA(p,q) after differencing it $d$ times\n\n* differencing is a discrete version of differentiation: polynomial of order $d$ becomes a constant after taking $d$th derivative $\\implies$ a series stationary around a linear trend becomes stationary in first differences (quadratic trend - second differences etc)\n\n","metadata":{}},{"cell_type":"code","source":"# quick demo: Tesla daily closing price on NYSE \nxdat = pd.read_csv(CFG.data_folder + 'tesla_prices_5y.csv', usecols = ['Date', 'Close'])\nxdat.columns = ['date', 'value']\n\n# xdat['value'] = np.log1p(xdat['value'])\nxdat['date'] = pd.to_datetime(xdat['date'])\nxdat.set_index('date').plot()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:52.957135Z","iopub.execute_input":"2022-02-14T01:10:52.957462Z","iopub.status.idle":"2022-02-14T01:10:53.221384Z","shell.execute_reply.started":"2022-02-14T01:10:52.95743Z","shell.execute_reply":"2022-02-14T01:10:53.220219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Applying Dickey-Fuller test confirms that the series certainly is not stationary:","metadata":{}},{"cell_type":"code","source":"result = adfuller(xdat['value'])\nprint('ADF Statistic: %f' % result[0])\nprint('p-value: %f' % result[1])","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:53.223471Z","iopub.execute_input":"2022-02-14T01:10:53.223794Z","iopub.status.idle":"2022-02-14T01:10:53.289878Z","shell.execute_reply.started":"2022-02-14T01:10:53.223762Z","shell.execute_reply":"2022-02-14T01:10:53.288161Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"With apologies for the pun, applying `.diff()` operator makes all the difference:","metadata":{}},{"cell_type":"code","source":"ydat = xdat['value'].diff()\nydat.plot()\nprint()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:53.291895Z","iopub.execute_input":"2022-02-14T01:10:53.292612Z","iopub.status.idle":"2022-02-14T01:10:53.71278Z","shell.execute_reply.started":"2022-02-14T01:10:53.292553Z","shell.execute_reply":"2022-02-14T01:10:53.711665Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"result = adfuller(ydat.dropna())\nprint('ADF Statistic: %f' % result[0])\nprint('p-value: %f' % result[1])","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:53.714004Z","iopub.execute_input":"2022-02-14T01:10:53.714314Z","iopub.status.idle":"2022-02-14T01:10:53.786211Z","shell.execute_reply.started":"2022-02-14T01:10:53.714283Z","shell.execute_reply":"2022-02-14T01:10:53.784143Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As you can see, differencing the time series once removes a linear trend and produces a series for which the null is not rejected in Dickey-Fuller test - this is despite the fact the volatility clearly goes up in the latter part of the sample!","metadata":{}},{"cell_type":"markdown","source":"The basic idea is that if we are dealing with a series that can be represented as a stationary process with polynomial trend, we can model it jointly (instead of manually differencing $d$ times and then flipping back). \n\n\nSideline point worth mentioning is that a non-seasonal ARIMA encapsulates other models as special cases (to quote a line from several textbooks I used in university: verifying those assertions is left as an exercise to the reader ;-)\n\n- ARIMA(0,0,0) is white noise\n- ARIMA(0,1,0) is random walk\n- ARIMA(p,0,0) is AR(p)\n- ARIMA(0,0,q) is MA(q)\n- ARIMA(0, 1, 1) is equivalent to single (basic) exponential smoothing model\n- ARIMA(0, 2, 2) is equivalent to Holt's linear method with additive errors, i.e. [double exponential smoothing](https://www.kaggle.com/konradb/ts-1a-smoothing-methods#Popular-methods)\n","metadata":{}},{"cell_type":"markdown","source":"### SARIMA\n\nThe \"I\" in ARIMA corresponds to the \"integrated\" component, which is a formal way of saying we are incorporating the trend into our setup starting with a stationary $ARMA(p,q)$ process. Following a similar logic, we can incorporate the seasonal component and allow it to follow the same type of dynamic: a seasonal ARIMA (SARIMA) model can be denoted as ARIMA(p,d,q)(P,D,Q)m, where:\n- m refers to the number of periods in each season\n- (lowercase) p,d,q refer to the definition of the ARIMA part\n- (uppercase) P,D,Q refer to the autoregressive, differencing, and moving average terms for the seasonal part of the ARIMA model.\n\n","metadata":{}},{"cell_type":"markdown","source":"The last part is not entirely obvious at an intuitive model (at least it wasn't to me when I learned it). A good intuition to think about is that we want to identify components of the seasonal part the way we would for the series itself:\n* we isolate the seasonal component, e.g. through seasonal decomposition\n* to establish the value of $D$ we test the seasonal component - is it stationary? is the differenced version stationary? \n* the seasonal autoregressive value $P$ can be established by inspecting the PACF of the seasonal component\n* the seasonal moving average value $Q$ can be established by inspect the ACF of the seasonal component\n\nBelow, we demonstrate briefly how this analysis can be conducted in practice:","metadata":{}},{"cell_type":"code","source":"series = pd.read_csv(CFG.data_folder + 'passengers.csv')\nseries['date'] = pd.to_datetime(series['date'])\n\n\ndecomposition = seasonal_decompose(series['passengers'], period = 12, model = 'multiplicative')\nfigure = decomposition.plot()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:53.788101Z","iopub.execute_input":"2022-02-14T01:10:53.788837Z","iopub.status.idle":"2022-02-14T01:10:54.499706Z","shell.execute_reply.started":"2022-02-14T01:10:53.788782Z","shell.execute_reply":"2022-02-14T01:10:54.498992Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# check stationarity\nresult = adfuller(decomposition.seasonal)\nprint('ADF Statistic: %f' % result[0])\nprint('p-value: %f' % result[1])\n","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:54.500693Z","iopub.execute_input":"2022-02-14T01:10:54.501079Z","iopub.status.idle":"2022-02-14T01:10:54.517694Z","shell.execute_reply.started":"2022-02-14T01:10:54.501049Z","shell.execute_reply":"2022-02-14T01:10:54.516896Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Non-stationarity is rejected extremely strongly, do $D$ = 0 in our example. What about the other components?","metadata":{}},{"cell_type":"code","source":"plot_pacf(decomposition.seasonal, lags = 12); print()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:54.518732Z","iopub.execute_input":"2022-02-14T01:10:54.519192Z","iopub.status.idle":"2022-02-14T01:10:54.786806Z","shell.execute_reply.started":"2022-02-14T01:10:54.519159Z","shell.execute_reply":"2022-02-14T01:10:54.785827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Relevant PACF values up until lag 12 at least indicate $P$ = 1 is the right choice.","metadata":{"execution":{"iopub.status.busy":"2022-02-10T20:04:02.807986Z","iopub.execute_input":"2022-02-10T20:04:02.808593Z","iopub.status.idle":"2022-02-10T20:04:02.815801Z","shell.execute_reply.started":"2022-02-10T20:04:02.808546Z","shell.execute_reply":"2022-02-10T20:04:02.814432Z"}}},{"cell_type":"code","source":"plot_acf(decomposition.seasonal, lags = 12); print()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:54.788236Z","iopub.execute_input":"2022-02-14T01:10:54.788831Z","iopub.status.idle":"2022-02-14T01:10:55.066044Z","shell.execute_reply.started":"2022-02-14T01:10:54.788776Z","shell.execute_reply":"2022-02-14T01:10:55.065048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Based on the ACF graph above $Q=1$ seems like a reasonable choice.","metadata":{}},{"cell_type":"markdown","source":"\nFor the more mathematically members of the audience, the SARIMA model can be formulated quite compactly if we make use of the backshift (lag) operator (https://en.wikipedia.org/wiki/Lag_operator).  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. For example, an $ARIMA(1,1,1)(1,1,1)_4$ model for quarterly data (m = 4) series $y_t$ can be specified as \n\n\\begin{equation}       \n    (1 - \\phi_1 B)(1- \\Phi_1 B^4)(1- B)(1-B^4)y_t = (1 + \\theta_1 B)(1 + \\Theta_1 B)\\epsilon_t\n\\end{equation}\n\n### SARIMAX\n\nSometimes we do have data available for the forecast horizon (e.g. long term economic forecasts available in advance, or the lagged values of the series itself). This means we can use **exogenous** variables to improve the quality of the forecast - hence SARIMA**X**.","metadata":{}},{"cell_type":"markdown","source":"<a id=\"section-three\"></a>\n# Full pipeline\n\nWe are in a position to combine all the building blocks and solve a real-life problem in full generality. We will use the data from the Demand Forecasting competition: \nhttps://www.kaggle.com/c/demand-forecasting-kernels-only","metadata":{}},{"cell_type":"code","source":"# load the data\ntrain = pd.read_csv('../input/demand-forecasting-kernels-only/train.csv' ,parse_dates=['date'],index_col='date')\ntest = pd.read_csv('../input/demand-forecasting-kernels-only/test.csv', parse_dates=['date'],index_col='date')\ndf = pd.concat([train,test],sort=True)","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:55.067262Z","iopub.execute_input":"2022-02-14T01:10:55.067559Z","iopub.status.idle":"2022-02-14T01:10:55.974717Z","shell.execute_reply.started":"2022-02-14T01:10:55.067528Z","shell.execute_reply":"2022-02-14T01:10:55.97357Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# we subset to one item x store combination\nbuf = df[(df.item==1)&(df.store==1)].copy()\nbuf.head(10)","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:55.976402Z","iopub.execute_input":"2022-02-14T01:10:55.976851Z","iopub.status.idle":"2022-02-14T01:10:56.057735Z","shell.execute_reply.started":"2022-02-14T01:10:55.976803Z","shell.execute_reply":"2022-02-14T01:10:56.056703Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# what do the components look like? \ndecomposition = seasonal_decompose(buf.sales.dropna(),freq=365)\nfigure = decomposition.plot()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:56.059322Z","iopub.execute_input":"2022-02-14T01:10:56.059946Z","iopub.status.idle":"2022-02-14T01:10:57.243128Z","shell.execute_reply.started":"2022-02-14T01:10:56.059894Z","shell.execute_reply":"2022-02-14T01:10:57.241666Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"With four years worth of data, we can think about capturing seasonal patterns (always worth remembering: modeling seasonal effects makes sense if you have at least two complete data cycles in your dataset). We begin by splitting into training and test:","metadata":{}},{"cell_type":"code","source":"tr_start,tr_end = '2014-01-01','2017-09-30'\nte_start,te_end = '2017-10-01','2017-12-31'\nx0 = buf['sales'][tr_start:tr_end].dropna()\nx1 = buf['sales'][te_start:te_end].dropna()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:57.244532Z","iopub.execute_input":"2022-02-14T01:10:57.244863Z","iopub.status.idle":"2022-02-14T01:10:57.259168Z","shell.execute_reply.started":"2022-02-14T01:10:57.244828Z","shell.execute_reply":"2022-02-14T01:10:57.258024Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# examine autocorrelation\nplot_acf(x0, lags = 12); print()\nplot_pacf(x0, lags = 12); print()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:57.260821Z","iopub.execute_input":"2022-02-14T01:10:57.261169Z","iopub.status.idle":"2022-02-14T01:10:57.739941Z","shell.execute_reply.started":"2022-02-14T01:10:57.261135Z","shell.execute_reply":"2022-02-14T01:10:57.73887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Indeed it does. We will conduct a quick grid search across possible values for the model order coefficients. ","metadata":{}},{"cell_type":"code","source":"model_autoARIMA = auto_arima(x0, start_p=7, start_q=7 ,\n                      test='adf',       \n                      max_p= 7, max_q=7, \n                      m= 7,              \n                      d= 1,\n                      seasonal=True,   \n                      start_P=1, \n                      D=1, \n                      trace=True,\n                      error_action='ignore',  \n                      suppress_warnings=True, \n                      stepwise=True)\n\n\n","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:10:57.741419Z","iopub.execute_input":"2022-02-14T01:10:57.741782Z","iopub.status.idle":"2022-02-14T01:12:32.445628Z","shell.execute_reply.started":"2022-02-14T01:10:57.741724Z","shell.execute_reply":"2022-02-14T01:12:32.444211Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(model_autoARIMA.summary())\n","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:12:32.447901Z","iopub.execute_input":"2022-02-14T01:12:32.448362Z","iopub.status.idle":"2022-02-14T01:12:32.46581Z","shell.execute_reply.started":"2022-02-14T01:12:32.448315Z","shell.execute_reply":"2022-02-14T01:12:32.464398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model_autoARIMA.plot_diagnostics()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:12:32.467198Z","iopub.execute_input":"2022-02-14T01:12:32.467495Z","iopub.status.idle":"2022-02-14T01:12:33.178253Z","shell.execute_reply.started":"2022-02-14T01:12:32.467467Z","shell.execute_reply":"2022-02-14T01:12:33.176978Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pred = model_autoARIMA.predict(x1.shape[0])\n\npd.DataFrame({'test':x1,'pred':pred}).plot();plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:12:33.179692Z","iopub.execute_input":"2022-02-14T01:12:33.180012Z","iopub.status.idle":"2022-02-14T01:12:33.464641Z","shell.execute_reply.started":"2022-02-14T01:12:33.179981Z","shell.execute_reply":"2022-02-14T01:12:33.463863Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As you can see from the graph above (and the diagnostics summary that preceded it), automatically fitted SARIMA does a decent job at capturing the dynamics of the underlying process - but there are still issues with capturing the range of values, with the most likely reasons being time-varying variance or distribution asymmetry. \n\n\nLet's compare its performance against [Prophet](https://www.kaggle.com/konradb/ts-1b-prophet).\n\n<a id=\"section-four\"></a>\n# A day at the races\n\nWe will be using data from the Nifty-50 index, which is an aggregate for the National Stock Exchange in India (an updated dataset is maintained courtesy of [Rohan](https://www.kaggle.com/rohanrao)): https://www.kaggle.com/rohanrao/nifty50-stock-market-data. Both models will be fitted with default settings, without incorporating special days and the like.","metadata":{}},{"cell_type":"code","source":"stock_name = 'TATASTEEL'\n\n# load the data\ndf = pd.read_csv('../input/nifty50-stock-market-data/'+stock_name+'.csv')\ndf.set_index(\"Date\", drop=False, inplace=True)\n\ndf.head(5)","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:15:28.689788Z","iopub.execute_input":"2022-02-14T01:15:28.690421Z","iopub.status.idle":"2022-02-14T01:15:28.753194Z","shell.execute_reply.started":"2022-02-14T01:15:28.690382Z","shell.execute_reply":"2022-02-14T01:15:28.752367Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plot the target \ndf.VWAP.plot()","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:15:31.671849Z","iopub.execute_input":"2022-02-14T01:15:31.67267Z","iopub.status.idle":"2022-02-14T01:15:31.917986Z","shell.execute_reply.started":"2022-02-14T01:15:31.672623Z","shell.execute_reply":"2022-02-14T01:15:31.916806Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# we can construct some simple features to be used as explanatory variables: \n# rolling statistics over the last 3 days for the price and volume characteristics\n\nlag_features = [\"High\", \"Low\", \"Volume\"]\nwindow_size = 3\n\ndf_rolled = df[lag_features].rolling(window = window_size, min_periods=0)\ndf_mean = df_rolled.mean().shift(1).reset_index()\ndf_std = df_rolled.std().shift(1).reset_index()\n\nfor feature in lag_features:\n    df[feature + '_mean_lag' + str(window_size)] = df_mean[feature].values\n    df[feature + '_std_lag' + str(window_size)] = df_std[feature].values\n\n# Prophet can handle missing values, but ARIMA cannot\ndf.fillna(df.mean(), inplace=True)","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:15:35.108842Z","iopub.execute_input":"2022-02-14T01:15:35.109594Z","iopub.status.idle":"2022-02-14T01:15:35.43497Z","shell.execute_reply.started":"2022-02-14T01:15:35.10955Z","shell.execute_reply":"2022-02-14T01:15:35.433829Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# prepare the training / validation split\ndf_train = df[df.Date < \"2019\"]\ndf_valid = df[df.Date >= \"2019\"]\n\nexogenous_features = ['High_mean_lag3', 'High_std_lag3', 'Low_mean_lag3',\n       'Low_std_lag3', 'Volume_mean_lag3', 'Volume_std_lag3',]","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:15:37.822877Z","iopub.execute_input":"2022-02-14T01:15:37.823303Z","iopub.status.idle":"2022-02-14T01:15:37.834108Z","shell.execute_reply.started":"2022-02-14T01:15:37.823263Z","shell.execute_reply":"2022-02-14T01:15:37.83282Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# fit ARIMA in the same manner as before \nmodel_arima = auto_arima(df_train.VWAP, exogenous=df_train[exogenous_features], \n                   m = 7,  \n                   # max ranges for the p,q,P,Q parameters - can be extended\n                   max_p= 2, max_q=2,    \n                   max_P = 1, max_Q = 1,\n                   trace=True, error_action=\"ignore\", suppress_warnings=True)\nmodel_arima.fit(df_train.VWAP, exogenous=df_train[exogenous_features])\n\nforecast = model_arima.predict(n_periods=len(df_valid), exogenous=df_valid[exogenous_features])\ndf_valid[\"Forecast_ARIMAX\"] = forecast","metadata":{"execution":{"iopub.status.busy":"2022-02-14T00:11:34.856208Z","iopub.execute_input":"2022-02-14T00:11:34.856943Z","iopub.status.idle":"2022-02-14T00:21:59.528911Z","shell.execute_reply.started":"2022-02-14T00:11:34.856884Z","shell.execute_reply":"2022-02-14T00:21:59.527677Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# fit a Prophet model with default parameters\n\ndf = df_train[['Date', 'VWAP']].rename(columns={\"Date\": \"ds\", \"VWAP\": \"y\"})\n\nmodel_prophet = Prophet()\n\n# add the regressors to the dataframe holding the data\nfor f in exogenous_features:    \n    df[f] = df_train[f]\n    model_prophet.add_regressor(f)\n\n# the rest proceeds as before. \nmodel_prophet.fit(df)","metadata":{"execution":{"iopub.status.busy":"2022-02-14T01:15:42.100638Z","iopub.execute_input":"2022-02-14T01:15:42.101259Z","iopub.status.idle":"2022-02-14T01:15:50.152634Z","shell.execute_reply.started":"2022-02-14T01:15:42.101212Z","shell.execute_reply":"2022-02-14T01:15:50.15162Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"forecast = model_prophet.predict(df_valid[[\"Date\", \"VWAP\"] + exogenous_features].rename(columns={\"Date\": \"ds\"}))\ndf_valid[\"Forecast_Prophet\"] = forecast.yhat.values\n","metadata":{"execution":{"iopub.status.busy":"2022-02-14T00:26:16.985542Z","iopub.execute_input":"2022-02-14T00:26:16.985857Z","iopub.status.idle":"2022-02-14T00:26:21.599297Z","shell.execute_reply.started":"2022-02-14T00:26:16.985826Z","shell.execute_reply":"2022-02-14T00:26:21.598227Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_valid[[\"VWAP\", \"Forecast_ARIMAX\", \"Forecast_Prophet\"]].plot()\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def friendly_rmse(x,y):\n    return np.round(np.sqrt(mean_squared_error(x, y))  ,4)\n\nprint(\"RMSE of Auto ARIMAX:\", friendly_rmse(df_valid.VWAP, df_valid.Forecast_ARIMAX))\nprint(\"RMSE of Prophet:\", friendly_rmse(df_valid.VWAP, df_valid.Forecast_Prophet))","metadata":{"execution":{"iopub.status.busy":"2022-02-14T00:26:30.697509Z","iopub.execute_input":"2022-02-14T00:26:30.698023Z","iopub.status.idle":"2022-02-14T00:26:30.707756Z","shell.execute_reply.started":"2022-02-14T00:26:30.697991Z","shell.execute_reply":"2022-02-14T00:26:30.706444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see from this small comparison, a vintage method like SARIMAX achieves comparable performance with a modern solution like Prophet. While it takes a bit more effort to set up, a full probabilistic model allows for more detailed and informative inference.\n\n","metadata":{}}]}