{"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): Ethereum 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":"markdown","source":"# 1- Import Libs","metadata":{}},{"cell_type":"code","source":"!pip install gluonts","metadata":{"execution":{"iopub.status.busy":"2021-11-10T13:07:46.368803Z","iopub.execute_input":"2021-11-10T13:07:46.369524Z","iopub.status.idle":"2021-11-10T13:07:55.754246Z","shell.execute_reply.started":"2021-11-10T13:07:46.369421Z","shell.execute_reply":"2021-11-10T13:07:55.753177Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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\nfrom gluonts.dataset.common import ListDataset\n\nimport math\nfrom math import ceil\n# Error libs\nfrom sklearn.metrics import r2_score\nfrom sklearn.metrics import explained_variance_score\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.metrics import mean_squared_log_error\nfrom sklearn.metrics import max_error\nfrom sklearn.metrics import median_absolute_error\n#from sklearn.metrics import mean_poisson_deviance\nfrom sklearn.metrics import mean_gamma_deviance\nfrom sklearn.metrics import mean_tweedie_deviance\nfrom math import sqrt","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-10T13:07:55.757844Z","iopub.execute_input":"2021-11-10T13:07:55.758323Z","iopub.status.idle":"2021-11-10T13:07:56.151089Z","shell.execute_reply.started":"2021-11-10T13:07:55.758291Z","shell.execute_reply":"2021-11-10T13:07:56.150087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2- Import Dataset","metadata":{}},{"cell_type":"code","source":"train = pd.read_csv(\"../input/g-research-crypto-forecasting/train.csv\")\n# extract the data corresponding to ethereum (Asset_ID = 6)\nethereum = train.query(\"Asset_ID == 6\").reset_index(drop = True)\nethereum['timestamp'] = pd.to_datetime(ethereum['timestamp'], unit='s')\nethereum = ethereum.set_index('timestamp')\n# extract the \"High\" value at 9:00 a.m. daily\nindexer_9am = ethereum.index.indexer_at_time('9:00:00')\nvalues_at_9am = ethereum.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-10T13:08:19.719164Z","iopub.execute_input":"2021-11-10T13:08:19.720152Z","iopub.status.idle":"2021-11-10T13:08:53.963970Z","shell.execute_reply.started":"2021-11-10T13:08:19.720100Z","shell.execute_reply":"2021-11-10T13:08:53.963070Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3- Create Model","metadata":{}},{"cell_type":"markdown","source":"We shall forecast three months (30 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 =  30 # predict these many days, these are removed from the end of the training data","metadata":{"execution":{"iopub.status.busy":"2021-11-10T13:08:58.386293Z","iopub.execute_input":"2021-11-10T13:08:58.386808Z","iopub.status.idle":"2021-11-10T13:08:58.392270Z","shell.execute_reply.started":"2021-11-10T13:08:58.386773Z","shell.execute_reply":"2021-11-10T13:08:58.390849Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Convert our dataframe into a GluonTS dataset","metadata":{}},{"cell_type":"code","source":"data_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=\"1D\")","metadata":{"execution":{"iopub.status.busy":"2021-11-10T13:09:00.157491Z","iopub.execute_input":"2021-11-10T13:09:00.158065Z","iopub.status.idle":"2021-11-10T13:09:00.164433Z","shell.execute_reply.started":"2021-11-10T13:09:00.158013Z","shell.execute_reply":"2021-11-10T13:09:00.163265Z"},"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=[60],\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-10T13:09:08.764648Z","iopub.execute_input":"2021-11-10T13:09:08.764970Z","iopub.status.idle":"2021-11-10T13:10:06.610203Z","shell.execute_reply.started":"2021-11-10T13:09:08.764898Z","shell.execute_reply":"2021-11-10T13:10:06.609064Z"},"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-10T13:10:06.612650Z","iopub.execute_input":"2021-11-10T13:10:06.613228Z","iopub.status.idle":"2021-11-10T13:10:06.774107Z","shell.execute_reply.started":"2021-11-10T13:10:06.613181Z","shell.execute_reply":"2021-11-10T13:10:06.773037Z"},"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, 95.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-10T13:10:06.775798Z","iopub.execute_input":"2021-11-10T13:10:06.776145Z","iopub.status.idle":"2021-11-10T13:10:07.307666Z","shell.execute_reply.started":"2021-11-10T13:10:06.776101Z","shell.execute_reply":"2021-11-10T13:10:07.306659Z"},"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-10T13:10:07.309542Z","iopub.execute_input":"2021-11-10T13:10:07.309876Z","iopub.status.idle":"2021-11-10T13:10:07.319147Z","shell.execute_reply.started":"2021-11-10T13:10:07.309842Z","shell.execute_reply":"2021-11-10T13:10:07.317941Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ControlTrain=high_values_at_9am.iloc[-30:,:]\nControlTrain[\"Prediction\"]=forecast_entry.mean\nControlTrain","metadata":{"execution":{"iopub.status.busy":"2021-11-10T13:10:07.321007Z","iopub.execute_input":"2021-11-10T13:10:07.321575Z","iopub.status.idle":"2021-11-10T13:10:07.353518Z","shell.execute_reply.started":"2021-11-10T13:10:07.321527Z","shell.execute_reply":"2021-11-10T13:10:07.352583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Evaluate the results","metadata":{}},{"cell_type":"code","source":"class errors:\n      \n    def mean_absolute_percentage(y, y_hat): \n        \n        #y_true=y_true.replace(0,0.001)\n        #y_pred=y_pred.replace(0,0.001)\n        y, y_hat = np.array(y), np.array(y_hat)\n        return np.mean(np.abs((y - y_hat) / y)) * 100    \n    \n    def mean_absolute_percentage_notabs(y, y_hat): \n        \n        #y_true=y_true.replace(0,0.001)\n        #y_pred=y_pred.replace(0,0.001)\n        #y, y_hat = np.array(y), np.array(y_hat)\n        return np.mean(((y - y_hat) / y)) * 100    \n    \n    def bias(y,yhat):    \n        y, y_hat = np.array(y), np.array(y_hat)\n        return np.sum(y - y_hat)\n        \n    def smape(y_true, y_pred):\n        denominator = (np.abs(y_true) + np.abs(y_pred))\n        diff = np.abs(y_true - y_pred) / denominator\n        diff[denominator == 0] = 0.0\n        return 200 * np.mean(diff)\n    def performance_metrics (y,yhat):\n        #print(\"y shape : \",y.shape , \" , y_hat shape : \", yhat.shape)\n        try:\n            '''\n            R² is the ratio between how good our model is vs how good is the naive mean model.\n            '''\n            R_2 = round(r2_score(y, yhat),3)                \n        except:\n            R_2 = \"EXC\" \n            \n        '''\n        Adjusted R2 required Regressor Number\n        '''\n\n        try:\n            '''\n            Explained variance regression score function\n            '''\n            EVS = round(explained_variance_score(y, yhat),3)\n        except:\n            EVS = \"EXC\"  \n            \n        try:\n            '''\n            Mean absolute error regression loss\n            it’s not that sensitive to outliers as mean square error. \n            MAE is more robust (less sensitive to outliers) than MSE\n            The measure gives an idea of the magnitude of the error, but no idea of the direction (e.g. over or under predicting).\n            '''\n            MAE = int(mean_absolute_error(y, yhat))                   \n        except:\n            MAE = \"EXC\" \n            \n        try:\n            '''\n            Mean squared error regression loss\n            MSE basically measures average squared error of our predictions. For each point, it calculates square difference between the predictions and the target and then average those values.\n            The higher this value, the worse the model is. It is never negative, since we’re squaring the individual prediction-wise errors before summing them, but would be zero for a perfect model .\n            Advantage: Useful if we have unexpected values that we should care about. Vey high or low value that we should pay attention.\n            Disadvantage: If we make a single very bad prediction, the squaring will make the error even worse and it may skew the metric towards overestimating the model’s badness. That is a particularly problematic behaviour if we have noisy data (that is, data that\n            for whatever reason is not entirely reliable) — even a “perfect” model may have a high MSE in that situation, so it becomes hard to judge how well the model is performing. On the other hand, if all the errors are small, or rather, smaller than 1, than the\n            opposite effect is felt: we may underestimate the model’s badness.\n            Note that if we want to have a constant prediction the best one will be the mean value of the target values. It can be found by setting the derivative of our total error with respect to that constant to zero, and find it from this equation.\n            '''\n            MSE = int(mean_squared_error(y, yhat))                   \n        except:\n            MSE = \"EXC\" \n            \n        try:    \n            '''\n            RMSE\n            '''\n            RMSE = round(mean_squared_error(y, yhat, squared=False),2)\n        except: \n            RMSE = \"EXC\" \n            \n        try:\n            '''\n            The bias is defined as the average error.\n            '''\n            BIAS = int (bias(y,yhat))\n        except:\n            BIAS = \"EXC\" \n            \n        try:\n            \n            MAPE= round(errors.mean_absolute_percentage(y, yhat),2)         \n        except:\n            MAPE= \"EXC\" \n            \n        try:\n            '''  \t\n            max_error metric calculates the maximum residual error.\n            '''            \n            MAXE = round (max_error(y, yhat),2)\n        except:            \n            MAXE =\"EXC\"\n            \n        try:            \n            '''\n            Mean squared logarithmic error regression loss\n            '''\n            MSLE = round(mean_squared_log_error(y, yhat),2)\n        except:            \n            MSLE =\"EXC\"\n        try:            \n            '''\n            Median absolute error regression loss\n            '''\n            MDAE = round(median_absolute_error(y, yhat),2)\n        except:            \n            MDAE =\"EXC\"\n        try:            \n            '''\n            Mean Poisson Deviance\n            '''\n            MPD = 0 # round(mean_poisson_deviance(y, yhat),2)\n        except:            \n            MPD =\"EXC\"\n            \n        try:            \n            '''\n            Mean Gamma Deviance\n            '''\n            MGD = round(mean_gamma_deviance(y, yhat),2)\n        except:            \n            MGD =\"EXC\"\n            \n        try:            \n            '''\n            Mean Tweedie Deviance\n            '''\n            MTD = round(mean_tweedie_deviance(y, yhat),2)\n        except:            \n            MTD =\"EXC\"\n        try:            \n            '''\n            SMAPE\n            '''\n            SMAPE = round(smape(y, yhat),2)\n        except:            \n            SMAPE =\"EXC\"     \n        \n        return R_2, EVS, MAE, MSE, RMSE, BIAS, MAPE, MTD, MGD, MPD, MSLE, MDAE, MAXE, SMAPE\n    \n    \n    \n    def results_dict (y,yhat):\n        \n        R_2, EVS, MAE, MSE, RMSE, BIAS, MAPE, MTD, MGD, MPD, MSLE, MDAE, MAXE,SMAPE= errors.performance_metrics (y,yhat)\n        \n        return {'r2' : R_2 , 'evs' : EVS , 'mae' : MAE , 'mse' : MSE , 'mape' : MAPE, 'rmse' : RMSE, 'bias' : BIAS, 'mtd' : MTD, 'mgd' : MGD, 'mpd' : MPD, 'msle' : MSLE, 'mdae' : MDAE, 'maxe' : MAXE, \"smape\": SMAPE}\n    \n    \n    \n    def print_results(y,y_hat):\n        \n        R_2, EVS, MAE, MSE, RMSE, BIAS, MAPE, MTD, MGD, MPD, MSLE, MDAE, MAXE,SMAPE = errors.performance_metrics (y,y_hat)\n        print ('----------------------------------------------------------------')\n        print ('                        MODEL ERROR METRICS                     ')\n        print ('----------------------------------------------------------------')\n        print (\"- R2 Score                       (R^2) = \" + str (R_2) )\n        print (\"- Explained Variance Score       (EVS) = \" + str (EVS) )\n        print (\"- Mean Absolute Error            (MAE) = \" + str (MAE) )\n        print (\"- Mean Squared Error             (MSE) = \" + str (MSE) )\n        print (\"- Mean Absolute Percentage Error (MAPE)= %\" + str (MAPE))\n        print (\"- Root Mean Squared Error        (RMSE)= \" + str (RMSE))\n        print (\"- BIAS...SUM(y-yhat)             (BIAS)= \" + str (BIAS))\n        print (\"- Mean Tweedie Deviance          (MTD) = \" + str (MTD) )\n        print (\"- Mean Gamma Deviance            (MGD) = \" + str (MGD) )\n        print (\"- Mean Poisson Deviance          (MPD) = \" + str (MPD) )\n        print (\"- Mean squared logarithmic Error (MSLE)= \" + str (MSLE))     \n        print (\"- Median absolute Error          (MDAE)= \" + str (MDAE))        \n        print (\"- Maximum Residual Error         (MAXE)= \" + str (MAXE))\n        print (\"-  S Mean Absolute Percentage Error (SMAPE)= %\" + str (SMAPE))\n        \n        print ('----------------------------------------------------------------')\n        print ('\\n')","metadata":{"execution":{"iopub.status.busy":"2021-11-10T13:10:16.966802Z","iopub.execute_input":"2021-11-10T13:10:16.967240Z","iopub.status.idle":"2021-11-10T13:10:16.999124Z","shell.execute_reply.started":"2021-11-10T13:10:16.967209Z","shell.execute_reply":"2021-11-10T13:10:16.997646Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"errors.print_results(ControlTrain.High,ControlTrain.Prediction)","metadata":{"execution":{"iopub.status.busy":"2021-11-10T13:10:20.132547Z","iopub.execute_input":"2021-11-10T13:10:20.132877Z","iopub.status.idle":"2021-11-10T13:10:20.151881Z","shell.execute_reply.started":"2021-11-10T13:10:20.132842Z","shell.execute_reply":"2021-11-10T13:10:20.150764Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}