{"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":"# M5 using mlforecast\n\n[mlforecast](https://nixtla.github.io/mlforecast/) is a framework to perform time series forecasting using machine learning models. It abstracts away most of the details and tries to mimic the scikit-learn API.\n\nThis notebook is inspired by https://www.kaggle.com/kneroma/m5-first-public-notebook-under-0-50.","metadata":{}},{"cell_type":"code","source":"%pip install mlforecast","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:37:32.223783Z","iopub.execute_input":"2022-06-30T02:37:32.224454Z","iopub.status.idle":"2022-06-30T02:37:49.208972Z","shell.execute_reply.started":"2022-06-30T02:37:32.224324Z","shell.execute_reply":"2022-06-30T02:37:49.207832Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Libraries","metadata":{}},{"cell_type":"code","source":"from pathlib import Path\n\nimport lightgbm as lgb\nimport numpy as np\nimport pandas as pd\nfrom mlforecast.core import TimeSeries\nfrom mlforecast.forecast import Forecast\nfrom window_ops.rolling import rolling_mean","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:37:49.211539Z","iopub.execute_input":"2022-06-30T02:37:49.211995Z","iopub.status.idle":"2022-06-30T02:37:52.564421Z","shell.execute_reply.started":"2022-06-30T02:37:49.211942Z","shell.execute_reply":"2022-06-30T02:37:52.563124Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data loading","metadata":{}},{"cell_type":"code","source":"input_path = Path('../input/m5-preprocess/processed/')\n\ndata = pd.read_parquet(input_path / 'sales.parquet')\ndata","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:37:52.566334Z","iopub.execute_input":"2022-06-30T02:37:52.566654Z","iopub.status.idle":"2022-06-30T02:38:02.525102Z","shell.execute_reply.started":"2022-06-30T02:37:52.566623Z","shell.execute_reply":"2022-06-30T02:38:02.523912Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"These are all sales in the dataset, however due to memory limitations we'll take from the 350th day onwards.","metadata":{}},{"cell_type":"code","source":"dates = sorted(data['date'].unique())\ndata = data[data['date'] >= dates[349]]\ndata.shape","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:38:02.527016Z","iopub.execute_input":"2022-06-30T02:38:02.527450Z","iopub.status.idle":"2022-06-30T02:38:06.667054Z","shell.execute_reply.started":"2022-06-30T02:38:02.527409Z","shell.execute_reply":"2022-06-30T02:38:06.665827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"mlforecast requires a dataframe with an index named **unique_id** which identifies each time serie, a column **ds** containing the datestamps and a column **y** with the series values.","metadata":{}},{"cell_type":"code","source":"data = data.rename(columns={'id': 'unique_id', 'date': 'ds'})\ndata = data.set_index('unique_id')\ndata","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:38:06.668495Z","iopub.execute_input":"2022-06-30T02:38:06.668866Z","iopub.status.idle":"2022-06-30T02:38:07.891644Z","shell.execute_reply.started":"2022-06-30T02:38:06.668832Z","shell.execute_reply":"2022-06-30T02:38:07.890633Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Metadata for predictions","metadata":{}},{"cell_type":"code","source":"prices = pd.read_parquet(input_path / 'prices.parquet')\nprices","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:38:07.892850Z","iopub.execute_input":"2022-06-30T02:38:07.893125Z","iopub.status.idle":"2022-06-30T02:38:07.966217Z","shell.execute_reply.started":"2022-06-30T02:38:07.893098Z","shell.execute_reply":"2022-06-30T02:38:07.965081Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cal = pd.read_parquet(input_path / 'calendar.parquet')\ncal = cal.rename(columns={'date': 'ds'})\ncal.head()","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:38:07.967871Z","iopub.execute_input":"2022-06-30T02:38:07.968282Z","iopub.status.idle":"2022-06-30T02:38:07.999216Z","shell.execute_reply.started":"2022-06-30T02:38:07.968240Z","shell.execute_reply":"2022-06-30T02:38:07.998139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Forecast setup","metadata":{}},{"cell_type":"markdown","source":"There are two inputs needed: a regressor that follows the scikit-learn API and a time series object which defines the features to be computed.","metadata":{}},{"cell_type":"markdown","source":"### Model","metadata":{}},{"cell_type":"code","source":"lgb_params = {\n    'objective': 'poisson',\n    'metric': 'rmse',\n    'learning_rate': 0.075,\n    'bagging_freq': 1,\n    'bagging_fraction': 0.75,\n    'lambda_l2': 0.1,\n    'n_estimators': 1200,\n    'num_leaves': 128,\n    'min_data_in_leaf': 100,\n}\n\nmodel = lgb.LGBMRegressor(**lgb_params)\nmodel","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:38:08.002397Z","iopub.execute_input":"2022-06-30T02:38:08.003080Z","iopub.status.idle":"2022-06-30T02:38:08.014553Z","shell.execute_reply.started":"2022-06-30T02:38:08.003023Z","shell.execute_reply":"2022-06-30T02:38:08.013804Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### TimeSeries\nThis is where we define the features. A brief description of each argument:\n\n* **freq**: frequency of our time series. This is a pandas abbreviation and is used to get the next dates when computing the predictions.\n* **lags**: lags that we want to use as features.\n* **lag_transforms**: dictionary where the keys are the lags that we want to use and the values are a list of transformations to apply to them. The transformations are defined as `numba` jitted functions. If the function takes more arguments than the input array, these are passed as a tuple `(func, arg1, arg2, ...)`.\n* **date_features**: date attributes to use for training. These are computed from the `ds` column and are updated in each timestep.\n* **num_threads**: number of threads to use in preprocessing and updates, defaults to all cpus. Since the transformations are `numba` jitted functions, we can use multithreading to compute our features.","metadata":{}},{"cell_type":"code","source":"ts = TimeSeries(\n    freq='D',\n    lags=[7, 28],\n    lag_transforms = {\n        7:  [(rolling_mean, 7), (rolling_mean, 28)],\n        28: [(rolling_mean, 7), (rolling_mean, 28)],\n    },\n    date_features=['year', 'month', 'day', 'dayofweek', 'quarter', 'week'],\n    num_threads=4,\n)\nts","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:38:08.016087Z","iopub.execute_input":"2022-06-30T02:38:08.016370Z","iopub.status.idle":"2022-06-30T02:38:08.026240Z","shell.execute_reply.started":"2022-06-30T02:38:08.016343Z","shell.execute_reply":"2022-06-30T02:38:08.025290Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Define forecaster\nOnce we have our model and flow configuration setup, we instantiate a `Forecast` object with them.","metadata":{}},{"cell_type":"code","source":"fcst = Forecast(model, ts)","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:38:08.027756Z","iopub.execute_input":"2022-06-30T02:38:08.028425Z","iopub.status.idle":"2022-06-30T02:38:08.036113Z","shell.execute_reply.started":"2022-06-30T02:38:08.028379Z","shell.execute_reply":"2022-06-30T02:38:08.035075Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training","metadata":{}},{"cell_type":"markdown","source":"If we only want to preprocess our data and train on all of it we can just call `Forecast.fit`. In this case we're going to make a train-valid split to get some information on the training, so we instead call `Forecast.preprocess` to get the dataframe with all our features and (internally) save the information for the forecasting step. `Forecast.preprocess` takes the following additional arguments:\n\n* **dropna**: whether or not to drop rows with null values after building all the features. Using lags and transformations on the lags generates many rows with `np.nan`s, this is a flag to indicate whether we want to drop them when we're done.\n* **keep_last_n**: keep only last `n` samples from each time serie after computing the features. The updates are performed by applying the transformations on the series again and taking only the last value. This can save memory if you have very long series and your transformations only use a small window, like in this case where we have series with thousands of data points and our transformations require only 28 (lag) + 27 (window) samples.\n* **static_features**: define which features are static. By default all extra columns (other than **ds** and **y**) are considered static and are replicated when building the features for the next timestep, setting this overrides that and repeats only the ones defined here.","metadata":{}},{"cell_type":"code","source":"%%time\nfeatures_df = fcst.preprocess(\n    data,\n    dropna=True,\n    keep_last_n=28+27,\n    static_features=['item_id', 'dept_id', 'cat_id', 'store_id', 'state_id'],\n)\ndel data","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:38:08.037560Z","iopub.execute_input":"2022-06-30T02:38:08.037888Z","iopub.status.idle":"2022-06-30T02:39:13.937633Z","shell.execute_reply.started":"2022-06-30T02:38:08.037859Z","shell.execute_reply":"2022-06-30T02:39:13.936255Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"features_df.columns","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:39:13.939079Z","iopub.execute_input":"2022-06-30T02:39:13.939471Z","iopub.status.idle":"2022-06-30T02:39:13.947306Z","shell.execute_reply.started":"2022-06-30T02:39:13.939437Z","shell.execute_reply":"2022-06-30T02:39:13.946056Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The names of the transformations are built using the name of the function and its arguments.\n\nThe order is always:\n1. Static features\n2. Lags\n3. Lag transforms\n4. Date features","metadata":{}},{"cell_type":"markdown","source":"Perform a train-valid split with 95% on train and 5% on valid.","metadata":{}},{"cell_type":"code","source":"np.random.seed(11)\ntrain_mask = np.random.rand(features_df.shape[0]) < 0.95\ntrain, valid = features_df[train_mask], features_df[~train_mask]\nX_train, y_train = train.drop(columns=['ds', 'y']), train.y\nX_valid, y_valid = valid.drop(columns=['ds', 'y']), valid.y\ndel features_df, train, valid","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:39:13.948811Z","iopub.execute_input":"2022-06-30T02:39:13.949221Z","iopub.status.idle":"2022-06-30T02:39:22.881952Z","shell.execute_reply.started":"2022-06-30T02:39:13.949177Z","shell.execute_reply":"2022-06-30T02:39:22.880896Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Calling `Forecast.fit` performs the preprocessing step as well. If we've already done that we just call `Forecast.model.fit` instead.","metadata":{}},{"cell_type":"code","source":"%time fcst.model.fit(X_train, y_train, eval_set=[(X_train, y_train), (X_valid, y_valid)], verbose=20)","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:39:22.883281Z","iopub.execute_input":"2022-06-30T02:39:22.883583Z","iopub.status.idle":"2022-06-30T02:40:41.225590Z","shell.execute_reply.started":"2022-06-30T02:39:22.883554Z","shell.execute_reply":"2022-06-30T02:40:41.224528Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Predictions","metadata":{}},{"cell_type":"markdown","source":"By default the predictions are computed repeating the static features and updating the transformations and the date features. If you want to do something different you can define your own predict function as explained [here](https://nixtla.github.io/mlforecast/forecast.html#Custom-predictions).","metadata":{}},{"cell_type":"code","source":"def my_predict_fn(\n    model,\n    new_x,\n    dynamic_dfs,\n    features_order,\n    alpha,\n) -> np.ndarray:\n    for df in dynamic_dfs:\n        new_x = new_x.merge(df, how='left')\n    predictions = model.predict(new_x[features_order])\n    return alpha * predictions","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:40:41.226887Z","iopub.execute_input":"2022-06-30T02:40:41.227178Z","iopub.status.idle":"2022-06-30T02:40:41.233285Z","shell.execute_reply.started":"2022-06-30T02:40:41.227151Z","shell.execute_reply":"2022-06-30T02:40:41.232148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Calling `Forecast.predict(horizon)` computes the predictions for the next `horizon` steps. We can also provide a custom `predict_fn` like we do in this case, using `my_predict_fn` defined above. This step uses multithreading if `num_threads` was set to a value greater than 1 or was left empty and you have more than 1 cpu (here we have 4).","metadata":{}},{"cell_type":"code","source":"fcst.ts.num_threads","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:40:41.234481Z","iopub.execute_input":"2022-06-30T02:40:41.234855Z","iopub.status.idle":"2022-06-30T02:40:41.248947Z","shell.execute_reply.started":"2022-06-30T02:40:41.234821Z","shell.execute_reply":"2022-06-30T02:40:41.247639Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nalphas = [1.028, 1.023, 1.018]\npreds = None\nfor alpha in alphas:\n    alpha_preds = fcst.predict(28, dynamic_dfs=[cal, prices], predict_fn=my_predict_fn, alpha=alpha)\n    alpha_preds = alpha_preds.set_index('ds', append=True)\n    if preds is None:\n        preds = 1 / 3 * alpha_preds\n    else:\n        preds += 1 / 3 * alpha_preds\npreds","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:40:41.251769Z","iopub.execute_input":"2022-06-30T02:40:41.252347Z","iopub.status.idle":"2022-06-30T02:41:05.435961Z","shell.execute_reply.started":"2022-06-30T02:40:41.252309Z","shell.execute_reply":"2022-06-30T02:41:05.434805Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Submission","metadata":{}},{"cell_type":"code","source":"wide = preds.reset_index().pivot_table(index='unique_id', columns='ds')\nwide.columns = [f'F{i+1}' for i in range(28)]\nwide.columns.name = None\nwide.index.name = 'id'\nwide","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:41:05.437387Z","iopub.execute_input":"2022-06-30T02:41:05.437670Z","iopub.status.idle":"2022-06-30T02:41:06.439141Z","shell.execute_reply.started":"2022-06-30T02:41:05.437642Z","shell.execute_reply":"2022-06-30T02:41:06.437835Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_sub = pd.read_csv(\n    '../input/m5-forecasting-accuracy/sample_submission.csv', index_col='id'\n)\nsample_sub.update(wide)\nnp.testing.assert_allclose(sample_sub.sum().sum(), preds['y_pred'].sum())\nsample_sub.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2022-06-30T02:41:06.440889Z","iopub.execute_input":"2022-06-30T02:41:06.441316Z","iopub.status.idle":"2022-06-30T02:41:09.657470Z","shell.execute_reply.started":"2022-06-30T02:41:06.441274Z","shell.execute_reply":"2022-06-30T02:41:09.656305Z"},"trusted":true},"execution_count":null,"outputs":[]}]}