{"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":"**This notebook is an exercise in the [Time Series](https://www.kaggle.com/learn/time-series) course.  You can reference the tutorial at [this link](https://www.kaggle.com/ryanholbrook/trend).**\n\n---\n","metadata":{}},{"cell_type":"markdown","source":"\n# Introduction #","metadata":{}},{"cell_type":"markdown","source":"Run this cell to set everything up!","metadata":{}},{"cell_type":"code","source":"# Setup feedback system\nfrom learntools.core import binder\nbinder.bind(globals())\nfrom learntools.time_series.ex2 import *\n\n# Setup notebook\nfrom pathlib import Path\nfrom learntools.time_series.style import *  # plot style settings\n\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport seaborn as sns\nfrom sklearn.linear_model import LinearRegression\n\ndata_dir = Path('../input/ts-course-data/')\ncomp_dir = Path('../input/store-sales-time-series-forecasting')\n\nretail_sales = pd.read_csv(\n    data_dir / \"us-retail-sales.csv\",\n    parse_dates=['Month'],\n    index_col='Month',\n).to_period('D')\nfood_sales = retail_sales.loc[:, 'FoodAndBeverage']\nauto_sales = retail_sales.loc[:, 'Automobiles']\n\ndtype = {\n    'store_nbr': 'category',\n    'family': 'category',\n    'sales': 'float32',\n    'onpromotion': 'uint64',\n}\nstore_sales = pd.read_csv(\n    comp_dir / 'train.csv',\n    dtype=dtype,\n    parse_dates=['date'],\n    infer_datetime_format=True,\n)\nstore_sales = store_sales.set_index('date').to_period('D')\nstore_sales = store_sales.set_index(['store_nbr', 'family'], append=True)\naverage_sales = store_sales.groupby('date').mean()['sales']","metadata":{"lines_to_next_cell":0,"execution":{"iopub.status.busy":"2022-07-13T19:26:03.474278Z","iopub.execute_input":"2022-07-13T19:26:03.475262Z","iopub.status.idle":"2022-07-13T19:26:12.800399Z","shell.execute_reply.started":"2022-07-13T19:26:03.475138Z","shell.execute_reply":"2022-07-13T19:26:12.799332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"-------------------------------------------------------------------------------","metadata":{}},{"cell_type":"markdown","source":"# 1) Determine trend with a moving average plot\n\nThe *US Retail Sales* dataset contains monthly sales data for a number of retail industries in the United States. Run the next cell to see a plot of the *Food and Beverage* series.","metadata":{}},{"cell_type":"code","source":"ax = food_sales.plot(**plot_params)\nax.set(title=\"US Food and Beverage Sales\", ylabel=\"Millions of Dollars\");","metadata":{"execution":{"iopub.status.busy":"2022-07-13T19:26:32.985500Z","iopub.execute_input":"2022-07-13T19:26:32.985864Z","iopub.status.idle":"2022-07-13T19:26:33.326466Z","shell.execute_reply.started":"2022-07-13T19:26:32.985834Z","shell.execute_reply":"2022-07-13T19:26:33.325329Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"food_sales.head()","metadata":{"execution":{"iopub.status.busy":"2022-07-13T19:29:57.947726Z","iopub.execute_input":"2022-07-13T19:29:57.948118Z","iopub.status.idle":"2022-07-13T19:29:57.958052Z","shell.execute_reply.started":"2022-07-13T19:29:57.948085Z","shell.execute_reply":"2022-07-13T19:29:57.956696Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now make a moving average plot to estimate the trend for this series.","metadata":{}},{"cell_type":"code","source":"average_sales","metadata":{"execution":{"iopub.status.busy":"2022-07-13T19:38:33.707083Z","iopub.execute_input":"2022-07-13T19:38:33.707650Z","iopub.status.idle":"2022-07-13T19:38:33.718365Z","shell.execute_reply.started":"2022-07-13T19:38:33.707607Z","shell.execute_reply":"2022-07-13T19:38:33.717292Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# YOUR CODE HERE: Add methods to `food_sales` to compute a moving\n# average with appropriate parameters for trend estimation.\ntrend = food_sales.rolling(\n    window=12,\n    center=True,\n    min_periods=6    \n).mean()\n\n# Check your answer\nq_1.check()\n\n# Make a plot\nax = food_sales.plot(**plot_params, alpha=0.5)\nax = trend.plot(ax=ax, linewidth=3)","metadata":{"lines_to_next_cell":0,"execution":{"iopub.status.busy":"2022-07-13T19:31:36.809837Z","iopub.execute_input":"2022-07-13T19:31:36.810263Z","iopub.status.idle":"2022-07-13T19:31:37.137508Z","shell.execute_reply.started":"2022-07-13T19:31:36.810213Z","shell.execute_reply":"2022-07-13T19:31:37.136342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Uncomment to get a hint or solution\n#q_1.hint()\n#q_1.solution()","metadata":{"lines_to_next_cell":0},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"-------------------------------------------------------------------------------\n\n# 2) Identify trend\n\nWhat order polynomial trend might be appropriate for the *Food and Beverage Sales* series? Can you think of a non-polynomial curve that might work even better?\n\nOnce you've thought about it, run this cell for some discussion.","metadata":{}},{"cell_type":"markdown","source":"First order? ","metadata":{}},{"cell_type":"code","source":"# View the solution (Run this cell to receive credit!)\nq_2.check()","metadata":{"execution":{"iopub.status.busy":"2022-07-13T19:33:18.075909Z","iopub.execute_input":"2022-07-13T19:33:18.076312Z","iopub.status.idle":"2022-07-13T19:33:18.086477Z","shell.execute_reply.started":"2022-07-13T19:33:18.076280Z","shell.execute_reply":"2022-07-13T19:33:18.085088Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"-------------------------------------------------------------------------------\n\nWe'll continue using the time series of average sales in this lesson. Run this cell to see a moving average plot of `average_sales` estimating the trend.","metadata":{}},{"cell_type":"code","source":"trend = average_sales.rolling(\n    window=365,\n    center=True,\n    min_periods=183,\n).mean()\n\nax = average_sales.plot(**plot_params, alpha=0.5)\nax = trend.plot(ax=ax, linewidth=3)","metadata":{"execution":{"iopub.status.busy":"2022-07-13T19:34:08.731175Z","iopub.execute_input":"2022-07-13T19:34:08.731584Z","iopub.status.idle":"2022-07-13T19:34:09.040001Z","shell.execute_reply.started":"2022-07-13T19:34:08.731553Z","shell.execute_reply":"2022-07-13T19:34:09.038734Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3) Create a Trend Feature\n\nUse `DeterministicProcess` to create a feature set for a cubic trend model. Also create features for a 90-day forecast.","metadata":{}},{"cell_type":"code","source":"from statsmodels.tsa.deterministic import DeterministicProcess\n\ny = average_sales.copy()  # the target\n\n# YOUR CODE HERE: Instantiate `DeterministicProcess` with arguments\n# appropriate for a cubic trend model\ndp = DeterministicProcess(\n    index=average_sales.index,\n    constant=False,\n    order=3,\n    drop=True\n)\n\n# YOUR CODE HERE: Create the feature set for the dates given in y.index\nX = dp.in_sample()\n\n# YOUR CODE HERE: Create features for a 90-day forecast.\nX_fore = dp.out_of_sample(steps=90)\n\n\n# Check your answer\nq_3.check()","metadata":{"execution":{"iopub.status.busy":"2022-07-13T19:42:24.339106Z","iopub.execute_input":"2022-07-13T19:42:24.339465Z","iopub.status.idle":"2022-07-13T19:42:24.357944Z","shell.execute_reply.started":"2022-07-13T19:42:24.339436Z","shell.execute_reply":"2022-07-13T19:42:24.356721Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Lines below will give you a hint or solution code\n#q_3.hint()\n#q_3.solution()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"You can see the a plot of the result by running the next cell.","metadata":{}},{"cell_type":"code","source":"model = LinearRegression()\nmodel.fit(X, y)\n\ny_pred = pd.Series(model.predict(X), index=X.index)\ny_fore = pd.Series(model.predict(X_fore), index=X_fore.index)\n\nax = y.plot(**plot_params, alpha=0.5, title=\"Average Sales\", ylabel=\"items sold\")\nax = y_pred.plot(ax=ax, linewidth=3, label=\"Trend\", color='C0')\nax = y_fore.plot(ax=ax, linewidth=3, label=\"Trend Forecast\", color='C3')\nax.legend();","metadata":{"execution":{"iopub.status.busy":"2022-07-13T19:42:58.749931Z","iopub.execute_input":"2022-07-13T19:42:58.750314Z","iopub.status.idle":"2022-07-13T19:42:59.120315Z","shell.execute_reply.started":"2022-07-13T19:42:58.750280Z","shell.execute_reply":"2022-07-13T19:42:59.119176Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"--------------------------------------------------------------------------------\n\nOne way to fit more complicated trends is to increase the order of the polynomial you use. To get a better fit to the somewhat complicated trend in *Store Sales*, we could try using an order 11 polynomial.","metadata":{}},{"cell_type":"code","source":"from statsmodels.tsa.deterministic import DeterministicProcess\n\ndp = DeterministicProcess(index=y.index, order=11)\nX = dp.in_sample()\n\nmodel = LinearRegression()\nmodel.fit(X, y)\n\ny_pred = pd.Series(model.predict(X), index=X.index)\n\nax = y.plot(**plot_params, alpha=0.5, title=\"Average Sales\", ylabel=\"items sold\")\nax = y_pred.plot(ax=ax, linewidth=3, label=\"Trend\", color='C0')\nax.legend();","metadata":{"execution":{"iopub.status.busy":"2022-07-13T19:43:32.865714Z","iopub.execute_input":"2022-07-13T19:43:32.866206Z","iopub.status.idle":"2022-07-13T19:43:33.329231Z","shell.execute_reply.started":"2022-07-13T19:43:32.866145Z","shell.execute_reply":"2022-07-13T19:43:33.327789Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4) Understand risks of forecasting with high-order polynomials\n\nHigh-order polynomials are generally not well-suited to forecasting, however. Can you guess why?","metadata":{}},{"cell_type":"code","source":"# View the solution (Run this cell to receive credit!)\nq_4.check()","metadata":{"execution":{"iopub.status.busy":"2022-07-13T19:43:52.751071Z","iopub.execute_input":"2022-07-13T19:43:52.752063Z","iopub.status.idle":"2022-07-13T19:43:52.761451Z","shell.execute_reply.started":"2022-07-13T19:43:52.752018Z","shell.execute_reply":"2022-07-13T19:43:52.760179Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Uncomment the next line for a hint\n#q_4.hint()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Run this cell to see the same 90-day forecast using an order 11 polynomial. Does it confirm your intuition?","metadata":{}},{"cell_type":"code","source":"X_fore = dp.out_of_sample(steps=90)\ny_fore = pd.Series(model.predict(X_fore), index=X_fore.index)\n\nax = y.plot(**plot_params, alpha=0.5, title=\"Average Sales\", ylabel=\"items sold\")\nax = y_pred.plot(ax=ax, linewidth=3, label=\"Trend\", color='C0')\nax = y_fore.plot(ax=ax, linewidth=3, label=\"Trend Forecast\", color='C3')\nax.legend();","metadata":{"execution":{"iopub.status.busy":"2022-07-13T19:44:09.207800Z","iopub.execute_input":"2022-07-13T19:44:09.208212Z","iopub.status.idle":"2022-07-13T19:44:09.909827Z","shell.execute_reply.started":"2022-07-13T19:44:09.208178Z","shell.execute_reply":"2022-07-13T19:44:09.908373Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"--------------------------------------------------------------------------------\n\n# (Optional) Fit trend with splines\n\n*Splines* are a nice alternative to polynomials when you want to fit a trend. The *Multivariate Adaptive Regression Splines* (MARS) algorithm in the `pyearth` library is powerful and easy to use. There are a lot of hyperparameters you may want to investigate.","metadata":{}},{"cell_type":"code","source":"from pyearth import Earth\n\n# Target and features are the same as before\ny = average_sales.copy()\ndp = DeterministicProcess(index=y.index, order=1)\nX = dp.in_sample()\n\n# Fit a MARS model with `Earth`\nmodel = Earth()\nmodel.fit(X, y)\n\ny_pred = pd.Series(model.predict(X), index=X.index)\n\nax = y.plot(**plot_params, title=\"Average Sales\", ylabel=\"items sold\")\nax = y_pred.plot(ax=ax, linewidth=3, label=\"Trend\")","metadata":{"execution":{"iopub.status.busy":"2022-07-13T19:44:44.625919Z","iopub.execute_input":"2022-07-13T19:44:44.626346Z","iopub.status.idle":"2022-07-13T19:44:51.648783Z","shell.execute_reply.started":"2022-07-13T19:44:44.626312Z","shell.execute_reply":"2022-07-13T19:44:51.647301Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Forecasting complicated trends like this will typically be difficult (if not impossible). With historical data, however, you can use splines to isolate other patterns in a time series by *detrending*.","metadata":{}},{"cell_type":"code","source":"y_detrended = y - y_pred   # remove the trend from store_sales\n\ny_detrended.plot(**plot_params, title=\"Detrended Average Sales\");","metadata":{"execution":{"iopub.status.busy":"2022-07-13T19:45:06.199666Z","iopub.execute_input":"2022-07-13T19:45:06.200077Z","iopub.status.idle":"2022-07-13T19:45:06.495007Z","shell.execute_reply.started":"2022-07-13T19:45:06.200042Z","shell.execute_reply":"2022-07-13T19:45:06.493765Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Keep Going #\n\n[**Model seasonality**](https://www.kaggle.com/ryanholbrook/seasonality), another common type of time dependence, with indicators and Fourier features.","metadata":{}},{"cell_type":"markdown","source":"---\n\n\n\n\n*Have questions or comments? Visit the [course discussion forum](https://www.kaggle.com/learn/time-series/discussion) to chat with other learners.*","metadata":{}}]}