{"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":"### Probabilistic forecasting using [GluonTS](https://ts.gluon.ai/index.html): Bitcoin example\n[Probabilistic forecasting](https://en.wikipedia.org/wiki/Probabilistic_forecasting), rather than providing a single point prediction, provides a probability distribution as the outcome.\nTo do this we shall be using the [GluonTS - Probabilistic Time Series Modeling](https://ts.gluon.ai/index.html) package. This notebook is heavily based on the [Quick Start Tutorial](https://ts.gluon.ai/tutorials/forecasting/quick_start_tutorial.html) and the [Extended Forecasting Tutorial](https://ts.gluon.ai/tutorials/forecasting/extended_tutorial.html) that are both provided with the package.","metadata":{}},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd\nimport matplotlib.pyplot as plt\n#plt.style.use('fivethirtyeight')\nplt.rcParams.update({'font.size': 12})\n#plt.rcParams[\"figure.figsize\"] = (12, 4)\nfrom datetime import date","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-11T20:29:02.079638Z","iopub.execute_input":"2021-11-11T20:29:02.080294Z","iopub.status.idle":"2021-11-11T20:29:02.103455Z","shell.execute_reply.started":"2021-11-11T20:29:02.080167Z","shell.execute_reply":"2021-11-11T20:29:02.102784Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Install GluonTS","metadata":{}},{"cell_type":"code","source":"!pip install -q gluonts","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-11-11T20:29:02.104707Z","iopub.execute_input":"2021-11-11T20:29:02.106464Z","iopub.status.idle":"2021-11-11T20:29:12.787013Z","shell.execute_reply.started":"2021-11-11T20:29:02.106420Z","shell.execute_reply":"2021-11-11T20:29:12.786068Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For our data we shall be using the dataset provided by the kaggle [G-Research Crypto Forecasting](https://www.kaggle.com/c/g-research-crypto-forecasting) competition. We shall extract the data corresponding to the Bitcoin cryptocurrency. From this Bitcoin data we shall use only one datapoint each day, arbitrarily chosen to be the **High** value at 9:00 a.m.","metadata":{}},{"cell_type":"code","source":"train = pd.read_csv(\"../input/g-research-crypto-forecasting/train.csv\")\n# extract the data corresponding to Bitcoin (Asset_ID = 1)\nbitcoin = train.query(\"Asset_ID == 1\").reset_index(drop = True)\nbitcoin['timestamp'] = pd.to_datetime(bitcoin['timestamp'], unit='s')\nbitcoin = bitcoin.set_index('timestamp')\n# extract the \"High\" value at 9:00 a.m. daily\nindexer_9am = bitcoin.index.indexer_at_time('9:00:00')\nvalues_at_9am = bitcoin.iloc[indexer_9am]\nhigh_values_at_9am = values_at_9am[[\"High\"]]\n# take a look\nhigh_values_at_9am","metadata":{"execution":{"iopub.status.busy":"2021-11-11T20:29:12.788668Z","iopub.execute_input":"2021-11-11T20:29:12.788891Z","iopub.status.idle":"2021-11-11T20:30:19.999748Z","shell.execute_reply.started":"2021-11-11T20:29:12.788864Z","shell.execute_reply":"2021-11-11T20:30:19.998963Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We shall forecast three months (90 days) worth of data, based on the 180 days prior to the start of the forecasting period","metadata":{}},{"cell_type":"code","source":"freq = \"1D\"             # the frequency of our data, here daily\ncontext_length    = 180 # train on this number of days\nprediction_length =  90 # predict these many days, these are removed from the end of the training data","metadata":{"execution":{"iopub.status.busy":"2021-11-11T20:30:20.002112Z","iopub.execute_input":"2021-11-11T20:30:20.002624Z","iopub.status.idle":"2021-11-11T20:30:20.006749Z","shell.execute_reply.started":"2021-11-11T20:30:20.002580Z","shell.execute_reply":"2021-11-11T20:30:20.006059Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Convert our dataframe into a GluonTS dataset","metadata":{}},{"cell_type":"code","source":"from gluonts.dataset.common import ListDataset\n\ndata_list = [{\"start\": \"2018-01-01 09:00:00\", \"target\": high_values_at_9am[c].values} for c in high_values_at_9am.columns]\ntrain_ds  = ListDataset(data_iter=data_list,freq=freq)","metadata":{"execution":{"iopub.status.busy":"2021-11-11T20:30:20.007822Z","iopub.execute_input":"2021-11-11T20:30:20.008499Z","iopub.status.idle":"2021-11-11T20:30:20.189014Z","shell.execute_reply.started":"2021-11-11T20:30:20.008471Z","shell.execute_reply":"2021-11-11T20:30:20.188299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Training\nHere we use the GulonTS’s pre-built [feedforward neural network estimator](https://ts.gluon.ai/api/gluonts/gluonts.model.simple_feedforward.html) `SimpleFeedForwardEstimator` in conjunction with a [trainer](https://ts.gluon.ai/api/gluonts/gluonts.mx.trainer.html). This feedforward network can be [substituted for a recurrent neural network (RNN)](https://ts.gluon.ai/tutorials/forecasting/extended_tutorial.html#From-feedforward-to-RNN), such as the  [DeepAREstimator](https://ts.gluon.ai/api/gluonts/gluonts.model.deepar.html).","metadata":{}},{"cell_type":"code","source":"from gluonts.model.simple_feedforward  import SimpleFeedForwardEstimator\n# to use the DeepAREstimator \n# from gluonts.model.deepar import DeepAREstimator\nfrom gluonts.mx.distribution.student_t import StudentTOutput\nfrom gluonts.mx import Trainer\n\nestimator = SimpleFeedForwardEstimator(num_hidden_dimensions=[50],\n                                       freq=freq,\n                                       context_length=context_length,\n                                       prediction_length=prediction_length,\n                                       distr_output=StudentTOutput(),\n                                       trainer=Trainer(epochs=50,\n                                                       learning_rate=1e-3,\n                                                       num_batches_per_epoch=100,\n                                                       patience=10))\n\npredictor = estimator.train(train_ds)\nprint(\"Done\")","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-11-11T20:30:20.190113Z","iopub.execute_input":"2021-11-11T20:30:20.190339Z","iopub.status.idle":"2021-11-11T20:31:19.754304Z","shell.execute_reply.started":"2021-11-11T20:30:20.190313Z","shell.execute_reply":"2021-11-11T20:31:19.753688Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Prediction\nWe shall now make (here 2000) predictions using our model","metadata":{}},{"cell_type":"code","source":"from gluonts.evaluation import make_evaluation_predictions\n\nforecast_it, ts_it = make_evaluation_predictions(\n    dataset=train_ds,  # dataset\n    predictor=predictor,  # predictor\n    num_samples=2000,  # number of sample paths we want for evaluation\n)\n\nforecasts = list(forecast_it)\ntss = list(ts_it)","metadata":{"execution":{"iopub.status.busy":"2021-11-11T20:31:19.755480Z","iopub.execute_input":"2021-11-11T20:31:19.755846Z","iopub.status.idle":"2021-11-11T20:31:19.870875Z","shell.execute_reply.started":"2021-11-11T20:31:19.755816Z","shell.execute_reply":"2021-11-11T20:31:19.870185Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Visualization\nAnd finally we shall now plot our probabilistic forecast (here with the 50% and 90% prediction intervals shown) along with the ground truth values","metadata":{}},{"cell_type":"code","source":"def plot_prob_forecasts(ts_entry, forecast_entry):\n    plot_length = context_length + prediction_length\n    prediction_intervals = (50.0, 90.0)\n    legend = [\"ground truth\", \"median prediction\"] + [f\"{k}% prediction interval\" for k in prediction_intervals][::-1]\n\n    fig, ax = plt.subplots(1, 1, figsize=(18, 7))\n    ts_entry[-plot_length:].plot(ax=ax)  # plot the time series\n    forecast_entry.plot(prediction_intervals=prediction_intervals, color='cadetblue')\n    plt.axvline(forecast_entry.start_date, color='g', lw=1) # end of train dataset\n    plt.grid(which=\"major\")\n    plt.legend(legend, loc=\"upper left\")\n    plt.show();\n    \nplot_prob_forecasts(tss[0], forecasts[0])","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-11T20:31:19.871849Z","iopub.execute_input":"2021-11-11T20:31:19.872035Z","iopub.status.idle":"2021-11-11T20:31:20.395990Z","shell.execute_reply.started":"2021-11-11T20:31:19.872012Z","shell.execute_reply":"2021-11-11T20:31:20.395285Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can extract the mean and median values of our forecast as follows:","metadata":{}},{"cell_type":"code","source":"# first entry of the forecast list\nforecast_entry = forecasts[0]\nprint(f\"Mean of the future window:\\n {forecast_entry.mean}\")\nprint(f\"0.5-quantile (median) of the future window:\\n {forecast_entry.quantile(0.5)}\")","metadata":{"execution":{"iopub.status.busy":"2021-11-11T20:31:20.397533Z","iopub.execute_input":"2021-11-11T20:31:20.398035Z","iopub.status.idle":"2021-11-11T20:31:20.407274Z","shell.execute_reply.started":"2021-11-11T20:31:20.397995Z","shell.execute_reply":"2021-11-11T20:31:20.406315Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It is worth mentioning that if desired [point forecasts](https://ts.gluon.ai/tutorials/forecasting/extended_tutorial.html#Point-forecasts-with-a-simple-feedforward-network) can also be made using GluonTS.\n\n### What about the minutes?\nIn the above example we have taken one sample per day, however, we have much more fine grained data at our disposal. Would we do better if we looked at the data on the scale of minutes? It was [Benoit Mandelbrot](https://en.wikipedia.org/wiki/Benoit_Mandelbrot) that was the first to suggest that stock data may be [fractal](https://en.wikipedia.org/wiki/Fractal) in nature; the data looking similar no matter what the scale. To quantify this somewhat here we shall calculate the [Higuchi fractal dimension](https://en.wikipedia.org/wiki/Higuchi_dimension) for the daily sample above, and for a similar sample size on the minute scale. To do this we shall use the [HFDA](https://github.com/hiroki-kojima/HFDA) package written by Hiroki Kojima. ","metadata":{}},{"cell_type":"code","source":"!pip install -q hfda\nimport hfda","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-11-11T20:31:20.410280Z","iopub.execute_input":"2021-11-11T20:31:20.410886Z","iopub.status.idle":"2021-11-11T20:31:30.443060Z","shell.execute_reply.started":"2021-11-11T20:31:20.410845Z","shell.execute_reply":"2021-11-11T20:31:30.442270Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"k_max = 15\nD = hfda.measure(np.concatenate(high_values_at_9am.values), k_max)\nprint(\"Higuchi fractal dimension = %.2f\" % D)","metadata":{"execution":{"iopub.status.busy":"2021-11-11T20:31:30.444338Z","iopub.execute_input":"2021-11-11T20:31:30.444554Z","iopub.status.idle":"2021-11-11T20:31:30.467391Z","shell.execute_reply.started":"2021-11-11T20:31:30.444528Z","shell.execute_reply":"2021-11-11T20:31:30.466749Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we shall create the minutes dataset","metadata":{}},{"cell_type":"code","source":"bitcoin_minutes = bitcoin[[\"High\"]].tail(1358)\n\nD = hfda.measure(np.concatenate(bitcoin_minutes.values), k_max)\nprint(\"Higuchi fractal dimension = %.2f\" %D)","metadata":{"execution":{"iopub.status.busy":"2021-11-11T20:31:30.468528Z","iopub.execute_input":"2021-11-11T20:31:30.468730Z","iopub.status.idle":"2021-11-11T20:31:30.489131Z","shell.execute_reply.started":"2021-11-11T20:31:30.468703Z","shell.execute_reply":"2021-11-11T20:31:30.488118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Indeed the values seem similar, with only a slightly lower value for the minute data . Let us repeat a probabilistic forecast, now for the minutes data","metadata":{}},{"cell_type":"code","source":"freq = \"1min\"           # the frequency of our data, now in minutes\ncontext_length    = 180 # train on this number of days\nprediction_length =  90 # predict these many days, these are removed from the end of the training data\n\ndata_list = [{\"start\": \"2021-09-20 01:23:00\", \"target\": bitcoin_minutes[c].values} for c in bitcoin_minutes.columns]\ntrain_ds  = ListDataset(data_iter=data_list,freq=freq)\n\nestimator = SimpleFeedForwardEstimator(num_hidden_dimensions=[50],\n                                       freq=freq,\n                                       context_length=context_length,\n                                       prediction_length=prediction_length,\n                                       distr_output=StudentTOutput(),\n                                       trainer=Trainer(epochs=50,\n                                                       learning_rate=1e-3,\n                                                       num_batches_per_epoch=100,\n                                                       patience=10))\n\npredictor = estimator.train(train_ds)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-11-11T20:31:30.490161Z","iopub.execute_input":"2021-11-11T20:31:30.490363Z","iopub.status.idle":"2021-11-11T20:32:28.316954Z","shell.execute_reply.started":"2021-11-11T20:31:30.490337Z","shell.execute_reply":"2021-11-11T20:32:28.315286Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"forecast_it, ts_it = make_evaluation_predictions(\n    dataset=train_ds,  # dataset\n    predictor=predictor,  # predictor\n    num_samples=2000,  # number of sample paths we want for evaluation\n)\n\nforecasts = list(forecast_it)\ntss = list(ts_it)\n\nplot_prob_forecasts(tss[0], forecasts[0])","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-11T20:32:28.318038Z","iopub.execute_input":"2021-11-11T20:32:28.318736Z","iopub.status.idle":"2021-11-11T20:32:28.909790Z","shell.execute_reply.started":"2021-11-11T20:32:28.318695Z","shell.execute_reply":"2021-11-11T20:32:28.909028Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Indeed it does seem that forecasting on the minute scale is just as challenging as forecasting on the daily scale.\n\n### Related packages\n* [NGBoost: Natural Gradient Boosting for Probabilistic Prediction](https://stanfordmlgroup.github.io/projects/ngboost/)\n\n### Related kaggle notebooks\n* [\"M5 - Sales Uncertainty Prediction\"](https://www.kaggle.com/allunia/m5-sales-uncertainty-prediction) written by [Laura Fink](https://www.kaggle.com/allunia)\n* [\"M5 Forecasting Competition GluonTS Template\"](https://www.kaggle.com/steverab/m5-forecasting-competition-gluonts-template) written by [Stephan Rabanser](https://www.kaggle.com/steverab)\n\n### Related reading\n* [Rob J. Hyndman and George Athanasopoulos \"*Forecasting: Principles and Practice*\", (3rd Edition)](https://otexts.com/fpp3/)\n* [Fotios Petropoulos, *et al. \"Forecasting: Theory and practice*\", arXiv:2012.03854 (2020)](https://arxiv.org/pdf/2012.03854.pdf)\n* [DeepAR Forecasting Algorithm](https://docs.aws.amazon.com/sagemaker/latest/dg/deepar.html)\n* [Benoit B. Mandelbrot \"*The Variation of Certain Speculative Prices*\", The Journal of Business **Vol. 36** pp. 394-419 (1963)](https://www.jstor.org/stable/2350970)\n* [Benoit B. Mandelbrot \"*The Variation of Some Other Speculative Prices*\", The Journal of Business **Vol. 40** pp. 393-413 (1967)](https://www.jstor.org/stable/2351623)\n* [Benoit B. Mandelbrot, Adlai J. Fisher and Laurent E. Calvet \"*A Multifractal Model of Asset Returns*\", Cowles Foundation Discussion Paper No. 1164 (1997)](https://ssrn.com/abstract=78588)\n* [T. Higuchi \"*Approach to an irregular time series on the basis of the fractal theory*\", Physica D: Nonlinear Phenomena **Vol. 31** pp. 277-283 (1988)](https://doi.org/10.1016/0167-2789(88)90081-4)","metadata":{}}]}