{"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":"code","source":"import pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport matplotlib.pyplot as plt","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-11-07T13:29:49.952842Z","iopub.execute_input":"2021-11-07T13:29:49.953552Z","iopub.status.idle":"2021-11-07T13:29:49.959246Z","shell.execute_reply.started":"2021-11-07T13:29:49.953514Z","shell.execute_reply":"2021-11-07T13:29:49.956828Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data","metadata":{}},{"cell_type":"code","source":"df=pd.read_csv('../input/g-research-crypto-forecasting/train.csv')\ndf.replace([np.inf, -np.inf], np.nan, inplace=True)","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:29:50.301822Z","iopub.execute_input":"2021-11-07T13:29:50.302320Z","iopub.status.idle":"2021-11-07T13:30:27.190455Z","shell.execute_reply.started":"2021-11-07T13:29:50.302284Z","shell.execute_reply":"2021-11-07T13:30:27.185892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df=df.dropna(how='any')","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:27.196492Z","iopub.execute_input":"2021-11-07T13:30:27.198237Z","iopub.status.idle":"2021-11-07T13:30:29.119239Z","shell.execute_reply.started":"2021-11-07T13:30:27.198193Z","shell.execute_reply":"2021-11-07T13:30:29.118474Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.isna().value_counts()","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:29.120487Z","iopub.execute_input":"2021-11-07T13:30:29.120743Z","iopub.status.idle":"2021-11-07T13:30:33.735169Z","shell.execute_reply.started":"2021-11-07T13:30:29.120709Z","shell.execute_reply":"2021-11-07T13:30:33.734463Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* Asset","metadata":{}},{"cell_type":"code","source":"asset=pd.read_csv('../input/g-research-crypto-forecasting/asset_details.csv')","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:33.737162Z","iopub.execute_input":"2021-11-07T13:30:33.737492Z","iopub.status.idle":"2021-11-07T13:30:33.748534Z","shell.execute_reply.started":"2021-11-07T13:30:33.737454Z","shell.execute_reply":"2021-11-07T13:30:33.747789Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"asset_map=lambda x : asset[asset.Asset_ID==x].Asset_Name.tolist()[0]","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:33.749814Z","iopub.execute_input":"2021-11-07T13:30:33.751145Z","iopub.status.idle":"2021-11-07T13:30:33.755367Z","shell.execute_reply.started":"2021-11-07T13:30:33.751108Z","shell.execute_reply":"2021-11-07T13:30:33.754694Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* Group\n\n","metadata":{}},{"cell_type":"code","source":"dfs={}\nfor i in np.unique(df.Asset_ID):\n    dfs[i]=df[df.Asset_ID==i].reset_index(drop=True)\n    dfs[i]=dfs[i].drop(['Asset_ID'],axis=1)\n    dfs[i]['date']=pd.to_datetime(dfs[i].timestamp,unit='s')\n    dfs[i]=dfs[i].set_index('date')\n","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:33.756553Z","iopub.execute_input":"2021-11-07T13:30:33.756916Z","iopub.status.idle":"2021-11-07T13:30:39.539620Z","shell.execute_reply.started":"2021-11-07T13:30:33.756883Z","shell.execute_reply":"2021-11-07T13:30:39.538892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dfs[0].head(3)","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:39.540855Z","iopub.execute_input":"2021-11-07T13:30:39.541298Z","iopub.status.idle":"2021-11-07T13:30:39.559659Z","shell.execute_reply.started":"2021-11-07T13:30:39.541260Z","shell.execute_reply":"2021-11-07T13:30:39.558846Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dfs[0].tail(3)","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:39.561136Z","iopub.execute_input":"2021-11-07T13:30:39.561489Z","iopub.status.idle":"2021-11-07T13:30:39.578640Z","shell.execute_reply.started":"2021-11-07T13:30:39.561449Z","shell.execute_reply":"2021-11-07T13:30:39.577845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* Plot trend (There's some data missing)","metadata":{}},{"cell_type":"code","source":"num_forecast_steps=60*24*3  #3 days\n\nyear='2020'","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:39.581130Z","iopub.execute_input":"2021-11-07T13:30:39.582082Z","iopub.status.idle":"2021-11-07T13:30:39.585527Z","shell.execute_reply.started":"2021-11-07T13:30:39.582047Z","shell.execute_reply":"2021-11-07T13:30:39.584678Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(12, 6))\nax = fig.add_subplot(1, 1, 1)\nfor k,v in dfs.items():\n    sns.lineplot(x=dfs[k].loc[year].index[:num_forecast_steps],y='Target',\n                 data=dfs[k].loc[year][:num_forecast_steps],ax=ax,label=asset_map(k))\nax.set_ylabel(\"Return\")\nax.set_xlabel(\"time\")\nfig.suptitle('A day of returns',fontsize=15)\n\nfig.autofmt_xdate()\nplt.legend(loc='upper right')","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:39.589902Z","iopub.execute_input":"2021-11-07T13:30:39.590381Z","iopub.status.idle":"2021-11-07T13:30:43.472372Z","shell.execute_reply.started":"2021-11-07T13:30:39.590337Z","shell.execute_reply":"2021-11-07T13:30:43.471675Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Modeling","metadata":{}},{"cell_type":"code","source":"import tensorflow as tf\nimport tensorflow_probability as tfp\n\nfrom tensorflow_probability import distributions as tfd  #distribution\nfrom tensorflow_probability import sts  #Bayesian structural time series models.","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:43.473879Z","iopub.execute_input":"2021-11-07T13:30:43.474309Z","iopub.status.idle":"2021-11-07T13:30:43.479626Z","shell.execute_reply.started":"2021-11-07T13:30:43.474253Z","shell.execute_reply":"2021-11-07T13:30:43.478908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Functions","metadata":{}},{"cell_type":"markdown","source":"* Plot a forecast distribution against the 'true' time series","metadata":{}},{"cell_type":"code","source":"def plot_forecast(x, y,\n                  forecast_mean, forecast_scale, forecast_samples,\n                  title, x_locator=None, x_formatter=None):\n    \n    #x:dates\n    \n    #y:target\n    colors = sns.color_palette()\n    c1, c2 = colors[0], colors[1]\n    fig = plt.figure(figsize=(12, 6))\n    ax = fig.add_subplot(1, 1, 1)\n    num_steps = len(y)\n    num_steps_forecast = forecast_mean.shape[-1]\n    \n    num_steps_train = num_steps - num_steps_forecast\n    \n    ax.plot(x, y, lw=2, color=c1, label='ground truth')\n    \n    forecast_steps=x[num_steps_train:]\n    \n    ax.plot(forecast_steps, forecast_samples.T, lw=1, color=c2, alpha=0.1)\n\n    ax.plot(forecast_steps, forecast_mean, lw=2, ls='--', color=c2,\n           label='forecast')\n    ax.fill_between(forecast_steps,\n                   forecast_mean-2*forecast_scale,\n                   forecast_mean+2*forecast_scale, color=c2, alpha=0.2)\n\n    ymin, ymax = min(np.min(forecast_samples), np.min(y)), max(np.max(forecast_samples), np.max(y))\n    yrange = ymax-ymin\n    ax.set_ylim([ymin - yrange*0.1, ymax + yrange*0.1])\n    ax.set_title(\"{}\".format(title))\n    ax.legend()\n\n    if x_locator is not None:\n        ax.xaxis.set_major_locator(x_locator)\n        ax.xaxis.set_major_formatter(x_formatter)\n        fig.autofmt_xdate()\n\n    return fig, ax","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:43.480867Z","iopub.execute_input":"2021-11-07T13:30:43.481637Z","iopub.status.idle":"2021-11-07T13:30:43.493518Z","shell.execute_reply.started":"2021-11-07T13:30:43.481602Z","shell.execute_reply":"2021-11-07T13:30:43.492877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* Plot the contributions of posterior components in a single figure.","metadata":{}},{"cell_type":"code","source":"def plot_components(dates,\n                    component_means_dict,\n                    component_stddevs_dict,\n                    x_locator=None,\n                    x_formatter=None):\n    colors = sns.color_palette()\n    c1, c2 = colors[0], colors[1]\n\n    axes_dict = collections.OrderedDict()\n    num_components = len(component_means_dict)\n    fig = plt.figure(figsize=(12, 2.5 * num_components))\n    for i, component_name in enumerate(component_means_dict.keys()):\n        component_mean = component_means_dict[component_name]\n        component_stddev = component_stddevs_dict[component_name]\n\n        ax = fig.add_subplot(num_components,1,1+i)\n        ax.plot(dates, component_mean, lw=2)\n        ax.fill_between(dates,\n                     component_mean-2*component_stddev,\n                     component_mean+2*component_stddev,\n                     color=c2, alpha=0.5)\n        ax.set_title(component_name)\n        if x_locator is not None:\n            ax.xaxis.set_major_locator(x_locator)\n            ax.xaxis.set_major_formatter(x_formatter)\n        axes_dict[component_name] = ax\n    fig.autofmt_xdate()\n    fig.tight_layout()\n    return fig, axes_dict","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:43.494936Z","iopub.execute_input":"2021-11-07T13:30:43.495286Z","iopub.status.idle":"2021-11-07T13:30:43.505620Z","shell.execute_reply.started":"2021-11-07T13:30:43.495252Z","shell.execute_reply":"2021-11-07T13:30:43.504914Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* Plot a time series against a model's one-step predictions","metadata":{}},{"cell_type":"code","source":"def plot_one_step_predictive(dates, observed_time_series,\n                             one_step_mean, one_step_scale,\n                             x_locator=None, x_formatter=None):\n    colors = sns.color_palette()\n    c1, c2 = colors[0], colors[1]\n\n    fig=plt.figure(figsize=(12, 6))\n    ax = fig.add_subplot(1,1,1)\n    num_timesteps = one_step_mean.shape[-1]\n    ax.plot(dates, observed_time_series, label=\"observed time series\", color=c1)\n    ax.plot(dates, one_step_mean, label=\"one-step prediction\", color=c2)\n    ax.fill_between(dates,\n                  one_step_mean - one_step_scale,\n                  one_step_mean + one_step_scale,\n                  alpha=0.1, color=c2)\n    ax.legend()\n\n    if x_locator is not None:\n        ax.xaxis.set_major_locator(x_locator)\n        ax.xaxis.set_major_formatter(x_formatter)\n        fig.autofmt_xdate()\n    fig.tight_layout()\n    return fig, ax","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:43.506699Z","iopub.execute_input":"2021-11-07T13:30:43.506956Z","iopub.status.idle":"2021-11-07T13:30:43.517587Z","shell.execute_reply.started":"2021-11-07T13:30:43.506922Z","shell.execute_reply":"2021-11-07T13:30:43.516911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* Build model\n\n    Local-Linear Trend model + minute-of-hour seasonal effect","metadata":{}},{"cell_type":"code","source":"def build_model(observed_time_series,num_seasons):  \n    trend = sts.LocalLinearTrend(observed_time_series=observed_time_series)\n    \n    seasonal = tfp.sts.Seasonal(num_seasons=num_seasons, observed_time_series=observed_time_series)\n    \n    model = sts.Sum([trend, seasonal], observed_time_series=observed_time_series)\n    \n    return model","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:43.520305Z","iopub.execute_input":"2021-11-07T13:30:43.520487Z","iopub.status.idle":"2021-11-07T13:30:43.530144Z","shell.execute_reply.started":"2021-11-07T13:30:43.520460Z","shell.execute_reply":"2021-11-07T13:30:43.529437Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Example :  LocalLinearTrend Model\n\n   * Forecast only by single variable \n\n   * Given Target Observation (2020 , first 600 minutes data)\n   \n   * Forecast next unseen Target (2020 ,next 300 minutes data)","metadata":{}},{"cell_type":"code","source":"#since data  is minutes wise, season is an hour\n\nasset=1 #bitcoin\nyear='2020'\nnum_seasons=60 #60 minutes--> 1 hours\nnum_obs=600 #600 minutes observation --> 10 hours\n\nmodel=build_model(dfs[asset].loc[year][:num_obs].Target,num_seasons)\n\n\n# Build the variational surrogate posteriors `qs`.\nvariational_posteriors = sts.build_factored_surrogate_posterior(model=model)","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:43.531468Z","iopub.execute_input":"2021-11-07T13:30:43.531942Z","iopub.status.idle":"2021-11-07T13:30:43.773150Z","shell.execute_reply.started":"2021-11-07T13:30:43.531905Z","shell.execute_reply":"2021-11-07T13:30:43.772423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"num_variational_steps = 100  #update steps\n\n#optimize observation \n\nloss_curve = tfp.vi.fit_surrogate_posterior(\n    \n    target_log_prob_fn=model.joint_log_prob(observed_time_series=dfs[asset].loc[year].Target[:num_obs]),\n    \n    surrogate_posterior=variational_posteriors,\n    \n    optimizer=tf.optimizers.Adam(learning_rate=0.1),\n    \n    num_steps=num_variational_steps)\n\nplt.plot(loss_curve)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:30:43.774326Z","iopub.execute_input":"2021-11-07T13:30:43.774594Z","iopub.status.idle":"2021-11-07T13:33:14.924682Z","shell.execute_reply.started":"2021-11-07T13:30:43.774542Z","shell.execute_reply":"2021-11-07T13:33:14.923899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Draw samples from the variational posterior.\nq_samples = variational_posteriors.sample(1000)","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:33:14.925790Z","iopub.execute_input":"2021-11-07T13:33:14.926058Z","iopub.status.idle":"2021-11-07T13:33:14.959322Z","shell.execute_reply.started":"2021-11-07T13:33:14.926021Z","shell.execute_reply":"2021-11-07T13:33:14.958652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Inferred parameters:\")\nfor param in model.parameters:\n    print(\"{}: {} +- {}\".format(param.name,\n                              np.mean(q_samples[param.name], axis=0),\n                              np.std(q_samples[param.name], axis=0)))","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:33:14.962160Z","iopub.execute_input":"2021-11-07T13:33:14.962351Z","iopub.status.idle":"2021-11-07T13:33:14.972330Z","shell.execute_reply.started":"2021-11-07T13:33:14.962327Z","shell.execute_reply":"2021-11-07T13:33:14.971588Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"num_forecast_steps=60*5 #5 hour forecast\n\nforecast_dist = sts.forecast(\n    model,\n    observed_time_series=dfs[asset].loc[year].Target[:num_obs],\n    parameter_samples=q_samples,\n    num_steps_forecast=num_forecast_steps)","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:33:14.973853Z","iopub.execute_input":"2021-11-07T13:33:14.974268Z","iopub.status.idle":"2021-11-07T13:33:18.121699Z","shell.execute_reply.started":"2021-11-07T13:33:14.974232Z","shell.execute_reply":"2021-11-07T13:33:18.120732Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"num_samples=100  #generate n forecast samples\n\nforecast_mean, forecast_scale, forecast_samples = (\n    forecast_dist.mean().numpy()[..., 0],\n    forecast_dist.stddev().numpy()[..., 0],\n    forecast_dist.sample(num_samples).numpy()[..., 0])","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:40:22.153731Z","iopub.execute_input":"2021-11-07T13:40:22.154312Z","iopub.status.idle":"2021-11-07T13:40:29.462721Z","shell.execute_reply.started":"2021-11-07T13:40:22.154273Z","shell.execute_reply":"2021-11-07T13:40:29.462017Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* All 100 series and mean and std of each time","metadata":{}},{"cell_type":"code","source":"forecast_mean.shape,forecast_scale.shape","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:40:33.082792Z","iopub.execute_input":"2021-11-07T13:40:33.083061Z","iopub.status.idle":"2021-11-07T13:40:33.087791Z","shell.execute_reply.started":"2021-11-07T13:40:33.083029Z","shell.execute_reply":"2021-11-07T13:40:33.087194Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"forecast_samples.shape","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:40:33.652654Z","iopub.execute_input":"2021-11-07T13:40:33.653130Z","iopub.status.idle":"2021-11-07T13:40:33.658406Z","shell.execute_reply.started":"2021-11-07T13:40:33.653095Z","shell.execute_reply":"2021-11-07T13:40:33.657599Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plot_forecast(\n    dfs[asset].loc[year].index[:num_obs+num_forecast_steps],\n    dfs[asset].loc[year].Target[:num_obs+num_forecast_steps],\n    forecast_mean, forecast_scale, forecast_samples,\n    x_locator=None,\n    x_formatter=None,\n    title=f\"Returns in {year}\")\n\n\nax.legend(loc=\"upper left\")\nax.set_ylabel(f\"Returns\")\nax.set_xlabel(\"minutes\")\nfig.autofmt_xdate()","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:40:36.183120Z","iopub.execute_input":"2021-11-07T13:40:36.183504Z","iopub.status.idle":"2021-11-07T13:40:36.598482Z","shell.execute_reply.started":"2021-11-07T13:40:36.183469Z","shell.execute_reply":"2021-11-07T13:40:36.597817Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_true= dfs[asset].loc[year].Target[num_obs:num_obs+num_forecast_steps]\ny_pred=forecast_mean\n\ntf.keras.metrics.mean_squared_error(y_true,y_pred)","metadata":{"execution":{"iopub.status.busy":"2021-11-07T13:40:38.321803Z","iopub.execute_input":"2021-11-07T13:40:38.322526Z","iopub.status.idle":"2021-11-07T13:40:38.332995Z","shell.execute_reply.started":"2021-11-07T13:40:38.322487Z","shell.execute_reply":"2021-11-07T13:40:38.332099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* I will learn more about tfp library and try more statistical model then update the version","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}