{"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":"#library\nimport os\nfrom datetime import datetime\nimport time\nfrom sklearn.preprocessing import StandardScaler\nimport plotly.express as px\nimport plotly.graph_objs as go\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport math\n\nimport statsmodels.api as sm\nfrom statsmodels.tsa.ar_model import AutoReg\nfrom statsmodels.tsa.arima_model import ARMA\nfrom statsmodels.tsa.stattools import adfuller\nfrom statsmodels.graphics.tsaplots import plot_acf\nfrom statsmodels.graphics.tsaplots import plot_pacf\nfrom statsmodels.tsa.arima_model import ARIMA\nfrom statsmodels.tsa.statespace.sarimax import SARIMAX\nfrom statsmodels.tsa.seasonal import seasonal_decompose\nfrom scipy import stats\nfrom itertools import product\n\nfrom plotly.subplots import make_subplots\n\nfrom statsmodels.graphics.tsaplots import plot_acf, plot_pacf\n\n\n\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n        \n        \n        \ninfo = pd.read_csv(\"/kaggle/input/g-research-crypto-forecasting/asset_details.csv\")\nctrain = pd.read_csv(\"/kaggle/input/g-research-crypto-forecasting/train.csv\")\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-01-04T08:03:33.022211Z","iopub.execute_input":"2022-01-04T08:03:33.022883Z","iopub.status.idle":"2022-01-04T08:04:34.645388Z","shell.execute_reply.started":"2022-01-04T08:03:33.022846Z","shell.execute_reply":"2022-01-04T08:04:34.644612Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ctrain.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:34.647067Z","iopub.execute_input":"2022-01-04T08:04:34.647327Z","iopub.status.idle":"2022-01-04T08:04:34.671108Z","shell.execute_reply.started":"2022-01-04T08:04:34.647294Z","shell.execute_reply":"2022-01-04T08:04:34.670367Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"info.Asset_Name.unique()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:34.672332Z","iopub.execute_input":"2022-01-04T08:04:34.672688Z","iopub.status.idle":"2022-01-04T08:04:34.684847Z","shell.execute_reply.started":"2022-01-04T08:04:34.672652Z","shell.execute_reply":"2022-01-04T08:04:34.68408Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ctrain.columns","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:34.687344Z","iopub.execute_input":"2022-01-04T08:04:34.687604Z","iopub.status.idle":"2022-01-04T08:04:34.693104Z","shell.execute_reply.started":"2022-01-04T08:04:34.687571Z","shell.execute_reply":"2022-01-04T08:04:34.692377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#time frame that we have\npd.to_datetime(ctrain['timestamp'], unit='s')","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:34.694264Z","iopub.execute_input":"2022-01-04T08:04:34.6949Z","iopub.status.idle":"2022-01-04T08:04:36.051708Z","shell.execute_reply.started":"2022-01-04T08:04:34.694863Z","shell.execute_reply":"2022-01-04T08:04:36.051052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Candlestick charts","metadata":{}},{"cell_type":"code","source":"#Dictionary for assets\nrename_dict = {}\nasset_details = info \nfor a in asset_details['Asset_ID']:\n    rename_dict[a] = asset_details[asset_details.Asset_ID == a].Asset_Name.values[0]\n\ndisplay(rename_dict)","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:36.053032Z","iopub.execute_input":"2022-01-04T08:04:36.053299Z","iopub.status.idle":"2022-01-04T08:04:36.070815Z","shell.execute_reply.started":"2022-01-04T08:04:36.053266Z","shell.execute_reply":"2022-01-04T08:04:36.07006Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Convert timestamp\ntrain  = ctrain.copy()\ntrain['timestamp'] = train['timestamp'].astype('datetime64[s]')\n#supplemental_train['timestamp'] = supplemental_train['timestamp'].astype('datetime64[s]')\n#example_test['timestamp'] = example_test['timestamp'].astype('datetime64[s]')\n\n#train['date'] = read_DatetimeIndex(train['timestamp']).date\n#supplemental_train['date'] = pd.DatetimeIndex(supplemental_train['timestamp']).date\n#example_test['date'] = pd.DatetimeIndex(example_test['timestamp']).date\n\n# Resample\ntrain_daily = pd.DataFrame()\n\nfor asset_id in asset_details.Asset_ID:\n    train_single = train[train.Asset_ID == asset_id].copy()\n\n    train_single_new = train_single[['timestamp','Count']].resample('D', on='timestamp').sum()\n    train_single_new['Open'] = train_single[['timestamp','Open']].resample('D', on='timestamp').first()['Open']\n    train_single_new['High'] = train_single[['timestamp','High']].resample('D', on='timestamp').max()['High']\n    train_single_new['Low'] = train_single[['timestamp','Low']].resample('D', on='timestamp').min()['Low']\n    train_single_new['Close'] = train_single[['timestamp','Close']].resample('D', on='timestamp').last()['Close']\n    train_single_new['Volume'] = train_single[['timestamp','Volume']].resample('D', on='timestamp').sum()['Volume']\n    # train_single_new['VWAP']\n    #train_single_new['Target'] = train_single[['timestamp','Target']].resample('D', on='timestamp').mean()['Target']\n    train_single_new['Asset_ID'] = asset_id\n\n    train_daily = train_daily.append(train_single_new.reset_index(drop=False))\ntrain_daily = train_daily.sort_values(by = ['timestamp', 'Asset_ID']).reset_index(drop=True)\n\ntrain_daily = train_daily.pivot(index='timestamp', columns='Asset_ID')[['Count', 'Open', 'High', 'Low', 'Close', 'Volume']]\ntrain_daily = train_daily.reset_index(drop=False)\n\ndisplay(train_daily.head(10))","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:36.072344Z","iopub.execute_input":"2022-01-04T08:04:36.07261Z","iopub.status.idle":"2022-01-04T08:04:47.505091Z","shell.execute_reply.started":"2022-01-04T08:04:36.072576Z","shell.execute_reply":"2022-01-04T08:04:47.50435Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#visualize Bitcoin for recent data rows - last 200 rows\ncrypto_df = ctrain\n\ncrypto_df.index = pd.to_datetime(crypto_df.timestamp, unit='s')\nbtc = crypto_df[crypto_df[\"Asset_ID\"]==1] # Asset_ID = 1 for Bitcoin\nbtc_mini = btc.iloc[-200:] # Select recent data rows\n\n\nimport plotly.graph_objects as go\n\nfig = go.Figure(data=[go.Candlestick(x=btc_mini.index, open=btc_mini['Open'], high=btc_mini['High'], low=btc_mini['Low'], close=btc_mini['Close'])])\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:47.506301Z","iopub.execute_input":"2022-01-04T08:04:47.506668Z","iopub.status.idle":"2022-01-04T08:04:49.091861Z","shell.execute_reply.started":"2022-01-04T08:04:47.506625Z","shell.execute_reply":"2022-01-04T08:04:49.091173Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training Data Distribution among differnet Assets ","metadata":{}},{"cell_type":"code","source":"data = ctrain\n\nasset_count= []\nfor i in range(14):\n    count = (data[\"Asset_ID\"]==i).sum()\n    asset_count.append(count)\nfig = px.bar(x = asset_details.sort_values(\"Asset_ID\")[\"Asset_Name\"],\n             y = asset_count , \n             color = asset_count ,\n             color_continuous_scale=\"Emrld\") \nfig.update_xaxes(title=\"Assets\")\nfig.update_yaxes(title = \"Number of Rows\")\nfig.update_layout(showlegend = True,\n    title = {\n        'text': 'Data Distribution ',\n        'y':0.95,\n        'x':0.5,\n        'xanchor': 'center',\n        'yanchor': 'top'})\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:49.09307Z","iopub.execute_input":"2022-01-04T08:04:49.093432Z","iopub.status.idle":"2022-01-04T08:04:50.438403Z","shell.execute_reply.started":"2022-01-04T08:04:49.093395Z","shell.execute_reply":"2022-01-04T08:04:50.437602Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"1. Data visualtion for 14 popular cryptocurrency\n2. Price History for selected individual cryptocurrency\n3. Basic Arima Model for price prediction","metadata":{}},{"cell_type":"markdown","source":"## Time history of 3 coins and its return\n\nref: https://www.kaggle.com/fangya/cryptocurrency-data-visualization-arima","metadata":{}},{"cell_type":"code","source":"# Impute Missing Time Value\n\ndef c_time_sub(asset_id, data=ctrain):\n    df=data[ctrain[\"Asset_ID\"]==asset_id].set_index(\"timestamp\")\n    df=df.reindex(range(df.index[0],df.index[-1]+60,60), method=\"pad\")\n    df.index = pd.to_datetime(df.index, unit='s')\n    return df\n\n\n#  Bitcoin\nbtc=c_time_sub(asset_id=1)\n\n#  Ethereum\neth=c_time_sub(asset_id=6)\n\n#  Cardano\nada=c_time_sub(asset_id=3)\n","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:50.443496Z","iopub.execute_input":"2022-01-04T08:04:50.444004Z","iopub.status.idle":"2022-01-04T08:04:52.319107Z","shell.execute_reply.started":"2022-01-04T08:04:50.443952Z","shell.execute_reply":"2022-01-04T08:04:52.318386Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"btc.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:52.320388Z","iopub.execute_input":"2022-01-04T08:04:52.320639Z","iopub.status.idle":"2022-01-04T08:04:52.336758Z","shell.execute_reply.started":"2022-01-04T08:04:52.320607Z","shell.execute_reply":"2022-01-04T08:04:52.335813Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"raw","source":"","metadata":{}},{"cell_type":"markdown","source":"## Cryptocurrency Log Return Correlation Plot for 2021","metadata":{}},{"cell_type":"code","source":"import time\n\n# define function to compute log returns\ndef log_return(series, periods=1):\n    return np.log(series).diff(periods=periods)\ndef plot_corr_2021(ctrain, info):\n    # auxiliary function, from datetime to timestamp\n    totimestamp = lambda s: np.int32(time.mktime(datetime.strptime(s, \"%d/%m/%Y\").timetuple()))\n    # create dataframe with returns for all assets\n    \n    crypto_df = ctrain\n    all_assets_2021 = pd.DataFrame([])\n    \n    for asset_id, asset_name in rename_dict.items():\n        print(asset_name)\n        asset = crypto_df[crypto_df[\"Asset_ID\"]==asset_id].set_index(\"timestamp\")\n        asset = asset.loc[totimestamp('01/01/2021'):]\n        asset = asset.reindex(range(asset.index[0],asset.index[-1]+60,60),method='pad')\n        lret = log_return(asset.Close.fillna(0))[1:]\n        \n        all_assets_2021 = all_assets_2021.join(lret, rsuffix=asset_name, how=\"outer\")\n    all_assets_2021.rename(columns = {'Close':'CloseBitcoin Cash'}, inplace=True)\n    sns.heatmap(all_assets_2021.corr());\n    \"\"\"plt.imshow(all_assets_2021.corr());\n    plt.yticks(info.Asset_ID.values, info.Asset_Name.values);\n    plt.xticks(info.Asset_ID.values, info.Asset_Name.values, rotation='vertical');\n    plt.colorbar();\n    plt.savefig('corr.jpg')\"\"\"\n    all_assets_2021_corr =  all_assets_2021.corr()\n    \"\"\"all_assets_2021_corr.columns = info.Asset_Name.values\n    all_assets_2021_corr.index = info.Asset_Name.values\"\"\"\n    return all_assets_2021, all_assets_2021_corr\nall_assets_2021, all_assets_2021_corr = plot_corr_2021(ctrain, info)","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:52.338557Z","iopub.execute_input":"2022-01-04T08:04:52.33978Z","iopub.status.idle":"2022-01-04T08:04:57.69369Z","shell.execute_reply.started":"2022-01-04T08:04:52.339738Z","shell.execute_reply":"2022-01-04T08:04:57.693042Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_assets_2021_corr.columns","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:57.694991Z","iopub.execute_input":"2022-01-04T08:04:57.69545Z","iopub.status.idle":"2022-01-04T08:04:57.701319Z","shell.execute_reply.started":"2022-01-04T08:04:57.695412Z","shell.execute_reply":"2022-01-04T08:04:57.700618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#the top correlation\ndef get_redundant_pairs(df):\n    '''Get diagonal and lower triangular pairs of correlation matrix'''\n    pairs_to_drop = set()\n    cols = df.columns\n    for i in range(0, df.shape[1]):\n        for j in range(0, i+1):\n            pairs_to_drop.add((cols[i], cols[j]))\n    return pairs_to_drop\n\ndef get_top_abs_correlations(df, n=5, ascending=False):\n    au_corr = df.corr().abs().unstack()\n    labels_to_drop = get_redundant_pairs(df)\n    au_corr = au_corr.drop(labels=labels_to_drop).sort_values(ascending=ascending)\n    return au_corr[0:n]\npd.DataFrame(get_top_abs_correlations(all_assets_2021_corr, n=10, ascending=False), columns=['Correlation'])\n","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:57.702638Z","iopub.execute_input":"2022-01-04T08:04:57.703085Z","iopub.status.idle":"2022-01-04T08:04:57.726433Z","shell.execute_reply.started":"2022-01-04T08:04:57.703051Z","shell.execute_reply":"2022-01-04T08:04:57.725444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#least correlated\npd.DataFrame(get_top_abs_correlations(all_assets_2021_corr, n=10, ascending=True), columns=['Correlation'])\n","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:57.727983Z","iopub.execute_input":"2022-01-04T08:04:57.728313Z","iopub.status.idle":"2022-01-04T08:04:57.748958Z","shell.execute_reply.started":"2022-01-04T08:04:57.728248Z","shell.execute_reply":"2022-01-04T08:04:57.748052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rename_dict[2]\n","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:57.750357Z","iopub.execute_input":"2022-01-04T08:04:57.750642Z","iopub.status.idle":"2022-01-04T08:04:57.760043Z","shell.execute_reply.started":"2022-01-04T08:04:57.750605Z","shell.execute_reply":"2022-01-04T08:04:57.759036Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#correlations over time\n# join two asset in single DataFrame\ndef corr_over_time(crypto_df, asset_id_1 =1, asset_id_2 = 6):\n    btc = crypto_df[crypto_df[\"Asset_ID\"]==asset_id_1].set_index(\"timestamp\")\n    eth= crypto_df[crypto_df[\"Asset_ID\"]==asset_id_2].set_index(\"timestamp\")\n\n    lret_btc_long = log_return(btc.Close)[1:]\n    lret_eth_long = log_return(eth.Close)[1:]\n    lret_btc_long.rename('lret_btc', inplace=True)\n    lret_eth_long.rename('lret_eth', inplace=True)\n    two_assets = pd.concat([lret_btc_long, lret_eth_long], axis=1)\n\n    # group consecutive rows and use .corr() for correlation between columns\n    corr_time = two_assets.groupby(two_assets.index//(10000*60)).corr().loc[:,\"lret_btc\"].loc[:,\"lret_eth\"]\n    return corr_time\n\n    #corr_time.plot();\n    \"\"\"plt.xticks([])\n    plt.ylabel(\"Correlation\")\n    plt.title(\"Correlation between BTC and ETH over time\");\"\"\"\n\nfig, ax = plt.subplots(nrows=14, ncols=14, figsize=(30, 25))\nfor i in range(0,13):\n    for j in range(i+1,14):\n        g = sns.lineplot(data = corr_over_time(ctrain,asset_id_1 =i, asset_id_2 = j),ax=ax[i,j])\n        g.set(xlabel=rename_dict[i], ylabel=rename_dict[j])\n","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:04:57.761152Z","iopub.execute_input":"2022-01-04T08:04:57.761854Z","iopub.status.idle":"2022-01-04T08:07:19.418649Z","shell.execute_reply.started":"2022-01-04T08:04:57.761814Z","shell.execute_reply":"2022-01-04T08:07:19.418004Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The illustration of BTC and ETH prices groupped by year [9] https://www.kaggle.com/iamleonie/to-the-moon-g-research-crypto-forecasting-eda#Data-Overview","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"The corr changes from time ","metadata":{}},{"cell_type":"markdown","source":"### Closing Price for BTC, ETH, ADA","metadata":{"execution":{"iopub.status.busy":"2021-12-25T06:57:07.881708Z","iopub.execute_input":"2021-12-25T06:57:07.882004Z","iopub.status.idle":"2021-12-25T06:57:07.887355Z","shell.execute_reply.started":"2021-12-25T06:57:07.881976Z","shell.execute_reply":"2021-12-25T06:57:07.886118Z"}}},{"cell_type":"markdown","source":"","metadata":{"execution":{"iopub.status.busy":"2021-12-24T06:58:20.162696Z","iopub.execute_input":"2021-12-24T06:58:20.162993Z","iopub.status.idle":"2021-12-24T06:58:20.170021Z","shell.execute_reply.started":"2021-12-24T06:58:20.162959Z","shell.execute_reply":"2021-12-24T06:58:20.168967Z"}}},{"cell_type":"code","source":"# Plot the Closing Price for BTC, ETH, ADA (Cardano)\nf= plt.figure(figsize=(10,12))  \n\ndef gplot(no , data, price, label, ylabel, color):\n    ax=f.add_subplot(no)\n    plt.plot(data[price], label=label, color=color)\n    plt.legend()\n    plt.xlabel(\"Time\")\n    plt.ylabel(ylabel)\n    return plt\n\ngplot(no=311, data=btc, price=\"Close\" , label=\"BTC 2021 Overall Performance\", ylabel=\"BTC Closing Price\", color=\"Lightskyblue\")\ngplot(no=312, data=eth, price=\"Close\" ,label=\"ETH 2021 Overall Performance\", ylabel=\"ETH Closing Price\", color=\"Coral\")\ngplot(no=313, data=ada, price=\"Close\" ,label=\"Cardano 2021 Overall Performance\", ylabel=\"ADA Closing Price\", color=\"khaki\")\n\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:07:19.419609Z","iopub.execute_input":"2022-01-04T08:07:19.419824Z","iopub.status.idle":"2022-01-04T08:07:24.875131Z","shell.execute_reply.started":"2022-01-04T08:07:19.419797Z","shell.execute_reply":"2022-01-04T08:07:24.874297Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Target : 15 minute resudualized returns\n#Residual Return: An asset's residual return equals its excess return minus beta times the benchmark excess return.\n\nf= plt.figure(figsize=(10,12))  \ngplot(no=311, data=btc, price=\"Target\" , label=\"BTC 2021 15min Return Residue\", ylabel=\"BTC residual return\", color=\"Aqua\")\ngplot(no=312, data=eth, price=\"Target\" ,label=\"ETH 2021 15min Return Residue\", ylabel=\"ETH residual return\", color=\"Pink\")\ngplot(no=313, data=ada, price=\"Target\" ,label=\"ADA 2021 15min Return Residue\", ylabel=\"ADA residual return\", color=\"gold\")\n\n\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:07:24.876486Z","iopub.execute_input":"2022-01-04T08:07:24.876891Z","iopub.status.idle":"2022-01-04T08:07:29.999181Z","shell.execute_reply.started":"2022-01-04T08:07:24.876854Z","shell.execute_reply":"2022-01-04T08:07:29.998526Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As we can see ETH, and BTC residual return are relatively stable compared to ADA. This might be a good implication that if the investor would take short time trading opportunities, ADA is a better choice.\n\nIf the investor is risk averse, BTC or ETH will be a better fit","metadata":{}},{"cell_type":"code","source":"#candlestick\n\ndef c_chart(data,label):\n    candlestick = go.Figure(data = [go.Candlestick(x =data.index, \n                                               open = data[('Open')], \n                                               high = data[('High')], \n                                               low = data[('Low')], \n                                               close = data[('Close')])])\n    candlestick.update_xaxes(title_text = 'Time',\n                             rangeslider_visible = True)\n\n    candlestick.update_layout(\n    title = {\n        'text': '{:} Candelstick Chart'.format(label),\n        \"y\":0.8,\n        \"x\":0.5,\n        'xanchor': 'center',\n        'yanchor': 'top'})\n\n    candlestick.update_yaxes(title_text = 'Price in USD', ticksuffix = '$')\n    return candlestick\n\n%matplotlib inline\nbtc_candle=c_chart(btc[-90:], label=\"BTC Price\")\nbtc_candle.show()\n\neth_candle=c_chart(eth[-90:], label=\"ETH Price \")\neth_candle.show()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:07:30.000489Z","iopub.execute_input":"2022-01-04T08:07:30.000869Z","iopub.status.idle":"2022-01-04T08:07:30.049787Z","shell.execute_reply.started":"2022-01-04T08:07:30.000834Z","shell.execute_reply":"2022-01-04T08:07:30.049091Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"## Technical Analysis charts\nAcross the industry, there are hundreds of patterns and signals that have been developed by researchers to support technical analysis trading, it would be a bit tedious to code them. \nFortunately there is a technical analysis library in python available(created by <a href=\"https://github.com/bukosabino\" >bukosabino</a>), the library has implemented 42 indicators:\n\n\n<b>Volume:</b>\n<ul>\n<li>Money Flow Index (MFI)</li>\n<li>Accumulation/Distribution Index (ADI)</li>\n<li>On-Balance Volume (OBV)</li>\n<li>Chaikin Money Flow (CMF)</li>\n<li>Force Index (FI)</li>\n<li>Ease of Movement (EoM, EMV)</li>\n<li>Volume-price Trend (VPT)</li>\n<li>Negative Volume Index (NVI)</li>\n<li>Volume Weighted Average Price (VWAP)</li>\n</ul>\n<b>Volatility:</b>\n<ul>\n<li>Average True Range (ATR)</li>\n<li>Bollinger Bands (BB)</li>\n<li>Keltner Channel (KC)</li>\n<li>Donchian Channel (DC)</li>\n<li>Ulcer Index (UI)</li>\n</ul>\n\n<b>Trend:</b>\n<ul>\n<li>Simple Moving Average (SMA)</li>\n<li>Exponential Moving Average (EMA)</li>\n<li>Weighted Moving Average (WMA)</li>\n<li>Moving Average Convergence Divergence (MACD)</li>\n<li>Average Directional Movement Index (ADX)</li>\n<li>Vortex Indicator (VI)</li>\n<li>Trix (TRIX)</li>\n<li>Mass Index (MI)</li>\n<li>Commodity Channel Index (CCI)</li>\n<li>Detrended Price Oscillator (DPO)</li>\n<li>KST Oscillator (KST)</li>\n<li>Ichimoku Kinkō Hyō (Ichimoku)</li>\n<li>Parabolic Stop And Reverse (Parabolic SAR)</li>\n<li>Schaff Trend Cycle (STC)</li>\n</ul>\n\n<b>Momentum:</b>\n<ul>\n<li>Relative Strength Index (RSI)</li>\n<li>Stochastic RSI (SRSI)</li>\n<li>True strength index (TSI)</li>\n<li>Ultimate Oscillator (UO)</li>\n<li>Stochastic Oscillator (SR)</li>\n<li>Williams %R (WR)</li>\n<li>Awesome Oscillator (AO)</li>\n<li>Kaufman's Adaptive Moving Average (KAMA)</li>\n<li>Rate of Change (ROC)</li>\n<li>Percentage Price Oscillator (PPO)</li>\n<li>Percentage Volume Oscillator (PVO)</li>\n</ul>\n<b>Others:</b>\n<ul>\n<li>Daily Return (DR)</li>\n<li>Daily Log Return (DLR)</li>\n<li>Cumulative Return (CR)</li>\n</ul>\n\n\n\n<b>Documentation:</b>\n<ul>\n    <li><a href = \"https://technical-analysis-library-in-python.readthedocs.io/en/latest/\">Technical Analysis library documentation</a></li>\n</ul>\n</div>\n\n","metadata":{}},{"cell_type":"markdown","source":"## Moving average","metadata":{"execution":{"iopub.status.busy":"2021-12-25T07:11:28.683455Z","iopub.execute_input":"2021-12-25T07:11:28.684183Z","iopub.status.idle":"2021-12-25T07:11:28.68863Z","shell.execute_reply.started":"2021-12-25T07:11:28.684141Z","shell.execute_reply":"2021-12-25T07:11:28.687804Z"}}},{"cell_type":"code","source":"ctrain.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:07:30.050863Z","iopub.execute_input":"2022-01-04T08:07:30.051125Z","iopub.status.idle":"2022-01-04T08:07:30.068633Z","shell.execute_reply.started":"2022-01-04T08:07:30.05109Z","shell.execute_reply":"2022-01-04T08:07:30.06793Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"btc.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:07:30.069734Z","iopub.execute_input":"2022-01-04T08:07:30.070782Z","iopub.status.idle":"2022-01-04T08:07:30.085433Z","shell.execute_reply.started":"2022-01-04T08:07:30.070745Z","shell.execute_reply":"2022-01-04T08:07:30.084701Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rename_dict","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:07:30.08655Z","iopub.execute_input":"2022-01-04T08:07:30.086921Z","iopub.status.idle":"2022-01-04T08:07:30.096621Z","shell.execute_reply.started":"2022-01-04T08:07:30.086886Z","shell.execute_reply":"2022-01-04T08:07:30.095855Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ndef mov_avg(df, a = 'W',Asset_ID =1, N = 10):\n    \n    df = df[df.Asset_ID == Asset_ID]\n    \n    df.set_index(df['timestamp'], inplace=True)\n    df.drop('timestamp', axis=1, inplace=True)\n    df.index = pd.to_datetime(df.index, unit='s')\n\n    #weekly\n    convertion={\n        'Open':'first',\n        'High':'max',\n        'Low':'min',\n        'Close':'mean',\n        'Volume':'sum',    \n    }\n    ds_df = df.resample(a).apply(convertion)\n    #avg # of mins in a week\n\n    \n    \n    #Simple ma\n    ds_df['rolling_mean' + str(a) + '_' + str(5)] = ds_df.Close.rolling(window=5).mean()\n    ds_df['rolling_mean' + str(a) + '_' + str(10)] = ds_df.Close.rolling(window=10).mean()\n    fig = go.Figure(go.Candlestick(x=ds_df.index,open=ds_df['Open'],high=ds_df['High'],low=ds_df['Low'],close=ds_df['Close']))\n    fig.update_layout(title='SMA '+str(rename_dict[Asset_ID])+' Close', yaxis_title=str(rename_dict[Asset_ID]))\n    fig.update_yaxes(type=\"log\")\n    fig.add_trace(go.Scatter(x=ds_df.index, y=ds_df['Close'],mode='lines',name='Close'))\n    fig.add_trace(go.Scatter(x=ds_df.index, y=ds_df['rolling_mean' + str(a) + '_' + str(5)], mode='lines', name='SMA MEAN of 5 ' + str(a),line=dict(color='royalblue', width=2)))\n    fig.add_trace(go.Scatter(x=ds_df.index, y=ds_df['rolling_mean' + str(a) + '_' + str(10)], mode='lines', name='SMA MEAN of 10 ' + str(a), line=dict(color='#555555', width=2)))\n    fig.show()\n    \n    \n    #Exp ma\n    ewma = pd.Series.ewm\n    ds_df['rolling_ema_'+ str(N)]  = ds_df.Close.ewm(min_periods=N, span=N).mean()\n    fig = go.Figure(go.Candlestick(x=ds_df.index,open=ds_df['Open'],high=ds_df['High'],low=ds_df['Low'],close=ds_df['Close']))\n    fig.update_layout(title='EMA ' +str(rename_dict[Asset_ID])+' Close', yaxis_title=str(rename_dict[Asset_ID]))\n    fig.update_yaxes(type=\"log\")\n    fig.add_trace(go.Scatter(x=ds_df.index, y=ds_df['Close'],mode='lines',name='Close'))\n    fig.add_trace(go.Scatter(x=ds_df.index, y=ds_df['rolling_ema_' + str(N)], mode='lines', name='EMA MEAN span of '+str(N)+' rolled by '+ str(a),line=dict(color='royalblue', width=2)))\n    fig.show()\n    \n    \n \n\n    \n    \n    \n    \n    \n    \nmov_avg(ctrain, Asset_ID =1 )\n","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:07:30.098078Z","iopub.execute_input":"2022-01-04T08:07:30.098429Z","iopub.status.idle":"2022-01-04T08:07:30.560691Z","shell.execute_reply.started":"2022-01-04T08:07:30.098393Z","shell.execute_reply":"2022-01-04T08:07:30.559938Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#ETH\nmov_avg(ctrain, Asset_ID =6 )","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:07:30.562106Z","iopub.execute_input":"2022-01-04T08:07:30.562364Z","iopub.status.idle":"2022-01-04T08:07:30.989266Z","shell.execute_reply.started":"2022-01-04T08:07:30.562331Z","shell.execute_reply":"2022-01-04T08:07:30.988391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Other technical analysis could be found here https://www.kaggle.com/shakshyathedetector/crypto-prediction-technical-analysis-features#Moving-average\nThere are:\n- Moving average convergence divergence (MACD) is a trend-following momentum indicator that shows the relationship between two moving averages of prices. The MACD is calculated by subtracting the 26-day exponential moving average (EMA) from the 12-day EMA\n- Bollinger Bands are a type of statistical chart characterizing the prices and volatility over time of a financial instrument or commodity, using a formulaic method propounded by John Bollinger in the 1980s. Financial traders employ these charts as a methodical tool to inform trading decisions, control automated trading systems, or as a component of technical analysis. Bollinger Bands display a graphical band (the envelope maximum and minimum of moving averages, similar to Keltner or Donchian channels) and volatility (expressed by the width of the envelope) in one two-dimensional chart.\n- A Volume Moving Average is the simplest volume-based technical indicator. Similar to a price moving average, a VMA is an average volume of a security (stock), commodity, index or exchange over a selected period of time. Volume Moving Averages are used in charts and in technical analysis to smooth and describe a volume trend by filtering short term spikes and gaps.","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Statistical analysis","metadata":{}},{"cell_type":"code","source":"info = pd.read_csv(\"/kaggle/input/g-research-crypto-forecasting/asset_details.csv\")\nctrain = pd.read_csv(\"/kaggle/input/g-research-crypto-forecasting/train.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:07:30.990774Z","iopub.execute_input":"2022-01-04T08:07:30.991265Z","iopub.status.idle":"2022-01-04T08:08:04.530912Z","shell.execute_reply.started":"2022-01-04T08:07:30.991214Z","shell.execute_reply":"2022-01-04T08:08:04.526647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ctrain.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:08:04.534632Z","iopub.execute_input":"2022-01-04T08:08:04.534891Z","iopub.status.idle":"2022-01-04T08:08:04.555773Z","shell.execute_reply.started":"2022-01-04T08:08:04.534857Z","shell.execute_reply":"2022-01-04T08:08:04.555079Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import time\n\n# define function to compute log returns\ndef log_return(series, periods=1):\n    return np.log(series).diff(periods=periods)\ndef logret_all_assets(ctrain, info):\n    # auxiliary function, from datetime to timestamp\n    totimestamp = lambda s: np.int32(time.mktime(datetime.strptime(s, \"%d/%m/%Y\").timetuple()))\n    # create dataframe with returns for all assets\n    \n    crypto_df = ctrain\n    all_assets_2021 = pd.DataFrame([])\n    \n    for asset_id, asset_name in rename_dict.items():\n        #print(asset_name)\n        asset = crypto_df[crypto_df[\"Asset_ID\"]==asset_id].set_index(\"timestamp\")\n        #asset = asset.loc[totimestamp('01/01/2021'):]\n        asset = asset.reindex(range(asset.index[0],asset.index[-1]+60,60),method='pad')\n        lret = log_return(asset.Close.fillna(0))[1:]\n        \n        all_assets_2021 = all_assets_2021.join(lret, rsuffix=asset_name, how=\"outer\")\n    all_assets_2021.rename(columns = {'Close':'Close Bitcoin Cash'}, inplace=True)\n    all_assets_2021.index  = pd.to_datetime(all_assets_2021.index, unit='s')\n    return all_assets_2021\nall_assets = logret_all_assets(ctrain, info)","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:08:04.560909Z","iopub.execute_input":"2022-01-04T08:08:04.561196Z","iopub.status.idle":"2022-01-04T08:08:13.646659Z","shell.execute_reply.started":"2022-01-04T08:08:04.561164Z","shell.execute_reply":"2022-01-04T08:08:13.64587Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#logret of all assets\nall_assets.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:08:13.647941Z","iopub.execute_input":"2022-01-04T08:08:13.6482Z","iopub.status.idle":"2022-01-04T08:08:13.669547Z","shell.execute_reply.started":"2022-01-04T08:08:13.648168Z","shell.execute_reply":"2022-01-04T08:08:13.668943Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_assets.to_csv('logret_all_years_assets.csv')","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:08:13.670688Z","iopub.execute_input":"2022-01-04T08:08:13.671048Z","iopub.status.idle":"2022-01-04T08:09:07.720885Z","shell.execute_reply.started":"2022-01-04T08:08:13.670996Z","shell.execute_reply":"2022-01-04T08:09:07.720092Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.boxplot(data = all_assets)","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:09:07.722225Z","iopub.execute_input":"2022-01-04T08:09:07.722508Z","iopub.status.idle":"2022-01-04T08:09:20.902655Z","shell.execute_reply.started":"2022-01-04T08:09:07.722473Z","shell.execute_reply":"2022-01-04T08:09:20.901975Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_assets.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:09:20.904016Z","iopub.execute_input":"2022-01-04T08:09:20.904282Z","iopub.status.idle":"2022-01-04T08:09:20.923269Z","shell.execute_reply.started":"2022-01-04T08:09:20.904247Z","shell.execute_reply":"2022-01-04T08:09:20.922343Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"stat = all_assets.describe()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:09:20.924974Z","iopub.execute_input":"2022-01-04T08:09:20.925307Z","iopub.status.idle":"2022-01-04T08:09:21.866454Z","shell.execute_reply.started":"2022-01-04T08:09:20.925263Z","shell.execute_reply":"2022-01-04T08:09:21.865707Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"stat","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:09:21.86785Z","iopub.execute_input":"2022-01-04T08:09:21.86811Z","iopub.status.idle":"2022-01-04T08:09:21.894875Z","shell.execute_reply.started":"2022-01-04T08:09:21.868074Z","shell.execute_reply":"2022-01-04T08:09:21.89427Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"stat.loc['mean',:].sort_values()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:09:21.895993Z","iopub.execute_input":"2022-01-04T08:09:21.896264Z","iopub.status.idle":"2022-01-04T08:09:21.905841Z","shell.execute_reply.started":"2022-01-04T08:09:21.896216Z","shell.execute_reply":"2022-01-04T08:09:21.905091Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"https://www.investopedia.com/articles/investing/102715/computing-historical-volatility-excel.asp\n r (t) = ln (P (t) / P (t-1)) ->  \n \nThe total return R at time t is:  \nR = r1 + r2 + r3 + 2 + ... +rt-1+ rt,\n\n","metadata":{}},{"cell_type":"code","source":"\n\nstat.loc['std',:].sort_values()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:09:21.907105Z","iopub.execute_input":"2022-01-04T08:09:21.907814Z","iopub.status.idle":"2022-01-04T08:09:21.917071Z","shell.execute_reply.started":"2022-01-04T08:09:21.907774Z","shell.execute_reply":"2022-01-04T08:09:21.916342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For all time minutes logreturn  ","metadata":{}},{"cell_type":"code","source":"all_assets.describe()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:09:21.918396Z","iopub.execute_input":"2022-01-04T08:09:21.918859Z","iopub.status.idle":"2022-01-04T08:09:22.818781Z","shell.execute_reply.started":"2022-01-04T08:09:21.918821Z","shell.execute_reply":"2022-01-04T08:09:22.818054Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_assets.columns","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:09:22.820187Z","iopub.execute_input":"2022-01-04T08:09:22.82068Z","iopub.status.idle":"2022-01-04T08:09:22.8263Z","shell.execute_reply.started":"2022-01-04T08:09:22.820642Z","shell.execute_reply":"2022-01-04T08:09:22.825565Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CREATE A FUNCTION THAT CALCULATE REALIZED VOLATILITY for btc\n# FROM minute-by-minute LOG RETURNS\ndef realized_volatility_daily(series_log_return):\n    \"\"\"\n    Get the daily realized volatility which is calculated as the square root\n    of sum of squares of log returns within a specific window interval \n    \"\"\"\n    n = len(series_log_return)\n    series_log_return = series_log_return.fillna(0)\n    return np.sqrt(np.sum(series_log_return**2)/(n - 1))\n\nk =24\nintervals = [ 7*k]\n#, 60 *k, 180 *k, 365*k\nvols_df = {}\n\n# ITERATE OVER intervals LIST\nfor i in intervals:\n    # GET DAILY LOG RETURNS USING THAT INTERVAL\n    vols = all_assets.CloseBitcoin.rolling(window=i)\\\n                         .apply(realized_volatility_daily).values\n\n    vols_df[i] = vols\n\n# CONVERT vols_df FROM DICTIONARY TO PANDAS DATAFRAME\nvols_df = pd.DataFrame(vols_df, columns=intervals, index=df.index)","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:09:22.82772Z","iopub.execute_input":"2022-01-04T08:09:22.828197Z","iopub.status.idle":"2022-01-04T08:13:22.46882Z","shell.execute_reply.started":"2022-01-04T08:09:22.82816Z","shell.execute_reply":"2022-01-04T08:13:22.467506Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"vols_df","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:13:22.469727Z","iopub.status.idle":"2022-01-04T08:13:22.47013Z","shell.execute_reply.started":"2022-01-04T08:13:22.469902Z","shell.execute_reply":"2022-01-04T08:13:22.469926Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# CHANGING MATPLOTLIB STYLE\nplt.style.use(['fivethirtyeight'])\n\nfig, ax = plt.subplots(figsize=(18,7))\n\nfor i in intervals:\n    if i == 7:\n        alpha = 0.5\n        lw = 1\n    else:\n        alpha = 1.0\n        lw = 2\n    ax.plot(vols_df[i], label=f'{i}-Day Interval Realized Volatility', \n            alpha=alpha, lw=lw)\n\nax.set_title('Realized Volatility Using Different Interval Windows', fontsize=21)\n\nplt.legend(loc='best', prop={'size': 14})\nplt.savefig(os.path.join(directory_to_img, 'diff_intervals.png'), \n            dpi=300, bbox_inches='tight')\nplt.show();","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:13:41.155661Z","iopub.execute_input":"2022-01-04T08:13:41.155967Z","iopub.status.idle":"2022-01-04T08:13:41.371856Z","shell.execute_reply.started":"2022-01-04T08:13:41.155935Z","shell.execute_reply":"2022-01-04T08:13:41.371069Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.boxplot(data = all_assets_2021)","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:13:43.81685Z","iopub.execute_input":"2022-01-04T08:13:43.817471Z","iopub.status.idle":"2022-01-04T08:13:45.789754Z","shell.execute_reply.started":"2022-01-04T08:13:43.817419Z","shell.execute_reply":"2022-01-04T08:13:45.789092Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"stat_2021 = all_assets_2021.describe()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:13:45.992942Z","iopub.execute_input":"2022-01-04T08:13:45.993293Z","iopub.status.idle":"2022-01-04T08:13:46.183064Z","shell.execute_reply.started":"2022-01-04T08:13:45.993258Z","shell.execute_reply":"2022-01-04T08:13:46.182349Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"stat_2021.loc['mean',:].sort_values()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:13:47.162803Z","iopub.execute_input":"2022-01-04T08:13:47.163058Z","iopub.status.idle":"2022-01-04T08:13:47.169869Z","shell.execute_reply.started":"2022-01-04T08:13:47.163021Z","shell.execute_reply":"2022-01-04T08:13:47.169268Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"stat_2021.loc['mean',:].sort_values()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:13:47.566628Z","iopub.execute_input":"2022-01-04T08:13:47.567204Z","iopub.status.idle":"2022-01-04T08:13:47.575589Z","shell.execute_reply.started":"2022-01-04T08:13:47.56717Z","shell.execute_reply":"2022-01-04T08:13:47.574534Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Dodgecoin has ","metadata":{}},{"cell_type":"code","source":"# sns.pairplot(all_assets_2021)","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:13:49.185218Z","iopub.execute_input":"2022-01-04T08:13:49.185938Z","iopub.status.idle":"2022-01-04T08:13:49.189628Z","shell.execute_reply.started":"2022-01-04T08:13:49.185903Z","shell.execute_reply":"2022-01-04T08:13:49.188584Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Imports\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom sklearn.cluster import KMeans\nfrom sklearn.decomposition import PCA\nfrom sklearn.preprocessing import StandardScaler\n%config InlineBackend.figure_format='retina'\nwarnings.filterwarnings( \"ignore\", module = \"matplotlib\\..*\" )\n\n\ndef vis_pca(all_assets):\n    X_std = StandardScaler().fit_transform(all_assets.fillna(0))\n    pca = PCA()\n    X_std_pca = pca.fit_transform(X_std)\n    plot = plt.scatter(X_std_pca[:,0], X_std_pca[:,1])\n    #plt.legend(handles=plot.legend_elements()[0])\n    plt.show()\n\n    # Plot the explained variances\n    features = range(pca.n_components_)\n    plt.bar(features, pca.explained_variance_ratio_, color='black')\n    plt.xlabel('PCA features')\n    plt.ylabel('variance %')\n    plt.xticks(features)\n    plt.show()\nvis_pca(all_assets)","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:13:51.405379Z","iopub.execute_input":"2022-01-04T08:13:51.405659Z","iopub.status.idle":"2022-01-04T08:13:59.787414Z","shell.execute_reply.started":"2022-01-04T08:13:51.405628Z","shell.execute_reply":"2022-01-04T08:13:59.786108Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Fig shows that the first one component explain the majority of the variance in our data. \n","metadata":{}},{"cell_type":"code","source":"\n\nvis_pca(all_assets.resample('M').mean())","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:13:22.486029Z","iopub.status.idle":"2022-01-04T08:13:22.486881Z","shell.execute_reply.started":"2022-01-04T08:13:22.486653Z","shell.execute_reply":"2022-01-04T08:13:22.486677Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"vis_pca(all_assets.resample('W').mean())","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:13:22.487933Z","iopub.status.idle":"2022-01-04T08:13:22.488791Z","shell.execute_reply.started":"2022-01-04T08:13:22.488562Z","shell.execute_reply":"2022-01-04T08:13:22.488586Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n\n\n\nX_std = StandardScaler().fit_transform(all_assets.fillna(0))\npca = PCA()\nX_std_pca = pca.fit_transform(X_std)\nplot = plt.scatter(X_std_pca[:,0], X_std_pca[:,1])\n#plt.legend(handles=plot.legend_elements()[0])\nplt.show()\n\n# Plot the explained variances\nfeatures = range(pca.n_components_)\nplt.bar(features, pca.explained_variance_ratio_, color='black')\nplt.xlabel('PCA features')\nplt.ylabel('variance %')\nplt.xticks(features)\n\n","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:13:22.489905Z","iopub.status.idle":"2022-01-04T08:13:22.490725Z","shell.execute_reply.started":"2022-01-04T08:13:22.490502Z","shell.execute_reply":"2022-01-04T08:13:22.490526Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"## ARIMA Model \n","metadata":{}},{"cell_type":"markdown","source":"Ref: https://www.kaggle.com/fangya/cryptocurrency-data-visualization-arima\n\nhttps://www.kaggle.com/iamleonie/to-the-moon-g-research-crypto-forecasting-eda#Data-Overview","metadata":{}},{"cell_type":"code","source":"btc.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### BTC Seasonal Decomposition\n","metadata":{}},{"cell_type":"code","source":"#freq https://stackoverflow.com/questions/17001389/pandas-resample-documentation\n#explanation https://towardsdatascience.com/time-series-decomposition-and-statsmodels-parameters-69e54d035453\nplt.rcParams[\"figure.figsize\"]=(15,7)\n\ndef season_df(data, label, a):\n    df=data.resample(a).mean()\n    seasonal_decompose(df.Close).plot()\n    print(label)\n    return plt.show()\nseason_df(data=btc, label=\"BTC Seasonal Decomposition\",a =\"M\")\n","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:07.295354Z","iopub.execute_input":"2022-01-04T08:14:07.295887Z","iopub.status.idle":"2022-01-04T08:14:08.294718Z","shell.execute_reply.started":"2022-01-04T08:14:07.29585Z","shell.execute_reply":"2022-01-04T08:14:08.294051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"season_df(data=btc, label=\"BTC Seasonal Decomposition\",a =\"D\")","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:10.835375Z","iopub.execute_input":"2022-01-04T08:14:10.835833Z","iopub.status.idle":"2022-01-04T08:14:12.112635Z","shell.execute_reply.started":"2022-01-04T08:14:10.835798Z","shell.execute_reply":"2022-01-04T08:14:12.110915Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"season_df(data=btc, label=\"BTC Seasonal Decomposition\",a =\"H\")","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:12.114318Z","iopub.execute_input":"2022-01-04T08:14:12.114749Z","iopub.status.idle":"2022-01-04T08:14:14.141857Z","shell.execute_reply.started":"2022-01-04T08:14:12.114715Z","shell.execute_reply":"2022-01-04T08:14:14.140902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n# CALCULATE PRICE RETURNS AS DAILY PERCENTAGE CHANGE USING pct_change()\nbtc['returns'] = 100 * btc.Close.pct_change().dropna()\n\n# CALCULATE LOG RETURNS BASED ON ABOVE FORMULA\nbtc['log_returns'] = np.log(btc.Close/btc.Close.shift(1))\nbtc.dropna(inplace=True)\n\n# PLOT DISTRIBUTION PLOTS OF RETURNS & LOG RETURNS\n# AND VISUALLY COMPARE THEM WITH THE STANDARD NORMAL DISTRIBUTION\ndf = btc.copy()\nwith sns.axes_style(\"darkgrid\"):\n    fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(18,12))\n\n    axes[0][0].plot(df.returns, color='blue')\n    axes[0][0].set_title('Returns')\n\n    sns.distplot(df.returns, norm_hist=True, fit=stats.norm, color='blue',\n                bins=50, ax=axes[0][1])\n    axes[0][1].set_title('Returns')\n\n    axes[1][0].plot(df.log_returns, color='green')\n    axes[1][0].set_title('Log Returns')\n\n    sns.distplot(df.log_returns, norm_hist=True, fit=stats.norm, color='green',\n                bins=50, ax=axes[1][1])\n    axes[1][1].set_title('Log Returns')\n    plt.tight_layout()\n\n    fig.show();","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy.stats import kurtosis, skew\n\n#returns\nprint( 'returns: excess kurtosis of normal distribution (should be 0): {}'.format( kurtosis(df.returns) ))\nprint( 'returns: skewness of normal distribution (should be 0): {}'.format( skew(df.returns) ))\n\n#log_returns\nprint( 'log_returns: excess kurtosis of normal distribution (should be 0): {}'.format( kurtosis(df.log_returns) ))\nprint( 'log_returns: skewness of normal distribution (should be 0): {}'.format( skew(df.log_returns) ))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- skewness:\nreturns - positive, log_returns - negative\n\n- positive kurtosis (leptokurtic) - Both Returns & Log Returns show higher peak with thicker tails than the standard normal distribution.\n","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Box-Cox Transformation for Arima\nSometimes after applying Box-Cox with a particular value of lambda the process may look stationary. It is sometimes possible that even if after applying the Box-Cox transformation the series does not appear to be stationary, diagnostics from ARIMA modeling can then be used to decide if differencing or seasonal differencing might be useful to to remove polynomial trends or seasonal trends respectively. After that the result might be an ARMA model that is stationary. If diagnostics confirm the orders p an q for the ARMA model, the AR and MA parameters can then be estimated.\n\nhttps://www.kaggle.com/taniaj/cryptocurrency-predictions-with-arima\n","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"btc.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:15.664458Z","iopub.execute_input":"2022-01-04T08:14:15.664721Z","iopub.status.idle":"2022-01-04T08:14:15.67885Z","shell.execute_reply.started":"2022-01-04T08:14:15.664693Z","shell.execute_reply":"2022-01-04T08:14:15.678182Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Box-Cox Transformation\nbtc_month=btc.resample(\"M\").mean()\nbtc_month[\"close_box\"], lmbda=stats.boxcox(btc_month.Close)\nprint(\"Dickey–Fuller test: p=%f\" % adfuller(btc_month.close_box)[1])\n","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:16.064465Z","iopub.execute_input":"2022-01-04T08:14:16.064676Z","iopub.status.idle":"2022-01-04T08:14:16.151364Z","shell.execute_reply.started":"2022-01-04T08:14:16.064651Z","shell.execute_reply":"2022-01-04T08:14:16.150609Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The p-value indicates that series is still not stationary.\n\n","metadata":{}},{"cell_type":"markdown","source":"# Seasonal differentiation \n\nOne method of differencing data is seasonal differencing, which involves computing the difference between an observation and the corresponding observation in the previous year.","metadata":{}},{"cell_type":"code","source":"# Seasonal differentiation (12 months)\nbtc_month['box_diff_seasonal_12'] = btc_month.close_box - btc_month.close_box.shift(12)\nprint(\"Dickey–Fuller test: p=%f\" % adfuller(btc_month.box_diff_seasonal_12[12:])[1])","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:20.328375Z","iopub.execute_input":"2022-01-04T08:14:20.328636Z","iopub.status.idle":"2022-01-04T08:14:20.344181Z","shell.execute_reply.started":"2022-01-04T08:14:20.328607Z","shell.execute_reply":"2022-01-04T08:14:20.343302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The p-value indicates that series is still not stationary.\n\n","metadata":{}},{"cell_type":"code","source":"# Seasonal differentiation (3 months)\nbtc_month['box_diff_seasonal_3'] = btc_month.close_box - btc_month.close_box.shift(3)\nprint(\"Dickey–Fuller test: p=%f\" % adfuller(btc_month.box_diff_seasonal_3[3:])[1])","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:21.572684Z","iopub.execute_input":"2022-01-04T08:14:21.573074Z","iopub.status.idle":"2022-01-04T08:14:21.588539Z","shell.execute_reply.started":"2022-01-04T08:14:21.573032Z","shell.execute_reply":"2022-01-04T08:14:21.587509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The p-value indicates that series is stationary as the computed p-value is lower than the significance level alpha = 0.05.","metadata":{}},{"cell_type":"markdown","source":"### Regular differentiation \nSometimes it may be necessary to difference the data a second time to obtain a stationary time series, which is referred to as second order differencing","metadata":{}},{"cell_type":"code","source":"btc_month.box_diff_seasonal_12.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:24.964435Z","iopub.execute_input":"2022-01-04T08:14:24.964697Z","iopub.status.idle":"2022-01-04T08:14:24.970883Z","shell.execute_reply.started":"2022-01-04T08:14:24.96467Z","shell.execute_reply":"2022-01-04T08:14:24.970144Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Regular differentiation (second diff for 12 months)\nbtc_month['box_diff2'] = btc_month.box_diff_seasonal_12 - btc_month.box_diff_seasonal_12.shift(1)\n\n# STL-decomposition\nseasonal_decompose(btc_month.box_diff2[13:]).plot()   \nprint(\"Dickey–Fuller test: p=%f\" % adfuller(btc_month.box_diff2[13:])[1])\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:25.808224Z","iopub.execute_input":"2022-01-04T08:14:25.808782Z","iopub.status.idle":"2022-01-04T08:14:26.670617Z","shell.execute_reply.started":"2022-01-04T08:14:25.808745Z","shell.execute_reply":"2022-01-04T08:14:26.668895Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The p-value indicates that series is stationary as the computed p-value is lower than the significance level alpha = 0.05.","metadata":{}},{"cell_type":"markdown","source":"### Autocorrelation \nAutocorrelation is the correlation of a time series with the same time series lagged. It summarizes the strength of a relationship with an observation in a time series with observations at prior time steps.\n\nWe create autocorrelation factor (ACF) and partial autocorrelation factor (PACF) plots to identify patterns in the above data which is stationary on both mean and variance. The idea is to identify presence of AR and MA components in the residuals.","metadata":{}},{"cell_type":"code","source":"#autocorrelation_plot(btc_month.close)\nplot_acf(btc_month.Close[13:].values.squeeze(), lags=12)\n\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:29.734759Z","iopub.execute_input":"2022-01-04T08:14:29.735352Z","iopub.status.idle":"2022-01-04T08:14:30.092251Z","shell.execute_reply.started":"2022-01-04T08:14:29.735315Z","shell.execute_reply":"2022-01-04T08:14:30.091549Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There is a positive correlation with the first 9 lags that is perhaps significant for the first 3 lags.\n\nA good starting point for the AR parameter of the model may be 3.\n\nLets try out autocorrelation on the differences","metadata":{}},{"cell_type":"code","source":"# Initial approximation of parameters using Autocorrelation and Partial Autocorrelation Plots\nax = plt.subplot(211)\n# Plot the autocorrelation function\n#sm.graphics.tsa.plot_acf(btc_month.box_diff2[13:].values.squeeze(), lags=48, ax=ax)\nplot_acf(btc_month.box_diff2[13:].values.squeeze(), lags=12, ax=ax)\nax = plt.subplot(212)\n#sm.graphics.tsa.plot_pacf(btc_month.box_diff2[13:].values.squeeze(), lags=48, ax=ax)\nplot_pacf(btc_month.box_diff2[13:].values.squeeze(), lags=12, ax=ax)\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:31.799047Z","iopub.execute_input":"2022-01-04T08:14:31.799695Z","iopub.status.idle":"2022-01-04T08:14:32.257732Z","shell.execute_reply.started":"2022-01-04T08:14:31.799657Z","shell.execute_reply":"2022-01-04T08:14:32.257034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There are not many spikes in the plots outside the insignificant zone (shaded) so there may not be enough information available in the residuals to be extracted by AR and MA models.","metadata":{}},{"cell_type":"markdown","source":"### Arima","metadata":{}},{"cell_type":"markdown","source":"AutoRegressive Integrated Moving Average\n\nARIMA models are denoted with the notation ARIMA(p, d, q). These parameters account for seasonality, trend, and noise in datasets:","metadata":{}},{"cell_type":"markdown","source":"p - the number of lag observations to include in the model, or lag order. (AR)\nd - the number of times that the raw observations are differenced, or the degree of differencing. (I)\nq - the size of the moving average window, also called the order of moving average.(MA)","metadata":{}},{"cell_type":"markdown","source":"A linear regression model is constructed including the specified number and type of terms, and the data is prepared by a degree of differencing in order to make it stationary, i.e. to remove trend and seasonal structures that negatively affect the regression model. A value of 0 for a parameter indicates to not use that element of the model.","metadata":{}},{"cell_type":"markdown","source":"### Parameter Selection\nWe will iteratively explore different combinations of parameters. For each combination we fit a new ARIMA model with SARIMAX() and assess its overall quality.\n\nWe will use the AIC (Akaike Information Criterion) value, returned with ARIMA models fitted using statsmodels. The AIC measures how well a model fits the data while taking into account the overall complexity of the model. A model that fits the data very well while using lots of features will be assigned a larger AIC score than a model that uses fewer features to achieve the same goodness-of-fit. Therefore, we are interested in finding the model that yields the lowest AIC value.","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#start=datetime(2018, 1, 31), end=datetime(2021, 5, 31)\n\nbtc_month.index","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:36.917038Z","iopub.execute_input":"2022-01-04T08:14:36.917324Z","iopub.status.idle":"2022-01-04T08:14:36.928012Z","shell.execute_reply.started":"2022-01-04T08:14:36.917294Z","shell.execute_reply":"2022-01-04T08:14:36.927188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"qs = range(0, 3)\nps = range(0, 3)\nd=1\nparameters = product(ps, qs)\nparameters_list = list(parameters)\nlen(parameters_list)\n\n# Model Selection\nresults = []\nbest_aic = float(\"inf\")\nwarnings.filterwarnings('ignore')\nfor param in parameters_list:\n    try:\n        model = SARIMAX(btc_month['2018-01-31':'2021-05-31'].close_box, order=(param[0], d, param[1])).fit(disp=-1)\n    except ValueError:\n        print('bad parameter combination:', param)\n        continue\n    aic = model.aic\n    if aic < best_aic:\n        best_model = model\n        best_aic = aic\n        best_param = param\n    results.append([param, model.aic])","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:38.691175Z","iopub.execute_input":"2022-01-04T08:14:38.691641Z","iopub.status.idle":"2022-01-04T08:14:40.196031Z","shell.execute_reply.started":"2022-01-04T08:14:38.691606Z","shell.execute_reply":"2022-01-04T08:14:40.195212Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"result_table = pd.DataFrame(results)\nresult_table.columns = ['parameters', 'aic']\nprint(result_table.sort_values(by = 'aic', ascending=True).head())\nprint(best_model.summary())","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:42.189928Z","iopub.execute_input":"2022-01-04T08:14:42.19019Z","iopub.status.idle":"2022-01-04T08:14:42.20681Z","shell.execute_reply.started":"2022-01-04T08:14:42.190162Z","shell.execute_reply":"2022-01-04T08:14:42.206102Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Dickey–Fuller test:: p=%f\" % adfuller(best_model.resid[13:])[1])\n","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:44.112931Z","iopub.execute_input":"2022-01-04T08:14:44.113492Z","iopub.status.idle":"2022-01-04T08:14:44.129064Z","shell.execute_reply.started":"2022-01-04T08:14:44.113452Z","shell.execute_reply":"2022-01-04T08:14:44.128391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Analysis of Results \nThe coef column shows the weight (i.e. importance) of each feature and how each one impacts the time series. The P>|z| column informs us of the significance of each feature weight. Here, each weight has a p-value lower or close to 0.05, so it is reasonable to retain all of them in our model.\n\nWhen fitting seasonal ARIMA models (and any other models for that matter), it is important to run model diagnostics to ensure that none of the assumptions made by the model have been violated. The plot_diagnostics object allows us to quickly generate model diagnostics and investigate for any unusual behavior.","metadata":{}},{"cell_type":"markdown","source":"Our primary concern is to ensure that the residuals of our model are uncorrelated and normally distributed with zero-mean. If the seasonal ARIMA model does not satisfy these properties, it is a good indication that it can be further improved.\nIn the histogram (top right), the KDE line should follow the N(0,1) line (normal distribution with mean 0, standard deviation 1) closely. This is an indication whether the residuals are normally distributed or not.\nIn the Q-Q-plot the ordered distribution of residuals (blue dots) should follow the linear trend of the samples taken from a standard normal distribution with N(0, 1). Again, this is an indication whether the residuals are normally distributed.\n\nThe standardized residual(top left) plot doesn't display any obvious seasonality. \nThis is confirmed by the Correlogram (bottom right), which shows that the time series residuals have low correlation with lagged versions of itself. (?)","metadata":{}},{"cell_type":"code","source":"best_model.plot_diagnostics(figsize=(15, 12))\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:46.145363Z","iopub.execute_input":"2022-01-04T08:14:46.145882Z","iopub.status.idle":"2022-01-04T08:14:47.030732Z","shell.execute_reply.started":"2022-01-04T08:14:46.145844Z","shell.execute_reply":"2022-01-04T08:14:47.030056Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"btc_month.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:48.952891Z","iopub.execute_input":"2022-01-04T08:14:48.953583Z","iopub.status.idle":"2022-01-04T08:14:48.970038Z","shell.execute_reply.started":"2022-01-04T08:14:48.953542Z","shell.execute_reply":"2022-01-04T08:14:48.969098Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Conclusion: We may consider trying to standardise the distribution further. But lets go ahead and do a prediction anyway...","metadata":{}},{"cell_type":"code","source":"#one-step ahead predictions\ndef invboxcox(y,lmbda):\n    if lmbda == 0:\n        return(np.exp(y))\n    else:\n        return(np.exp(np.log(lmbda*y+1)/lmbda))\n# Prediction\nbtc_month_pred = btc_month[['Close']]\n\nbtc_month_pred['forecast'] = invboxcox(best_model.predict(start=datetime(2018, 1, 31), end=datetime(2021, 9, 21)), lmbda)\nplt.figure(figsize=(18,10))\nbtc_month_pred.Close.plot()\nbtc_month_pred.forecast.plot(color='r', ls='--', label='Predicted Close')\nplt.legend()\nplt.title('Bitcoin monthly forecast')\nplt.ylabel('USD')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-01-04T08:14:51.066406Z","iopub.execute_input":"2022-01-04T08:14:51.066851Z","iopub.status.idle":"2022-01-04T08:14:51.57674Z","shell.execute_reply.started":"2022-01-04T08:14:51.066816Z","shell.execute_reply":"2022-01-04T08:14:51.576109Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## SARIMAX","metadata":{}},{"cell_type":"code","source":"##Addd SARIMAX (week 6)\n##Add rmse\n##Add other metric that is used in the competition","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Initial approximation of parameters\nfrom itertools import product\nimport warnings\nwarnings.filterwarnings('ignore')\nfrom statsmodels.tsa.statespace.sarimax import SARIMAX\n\n\nQs = range(0, 2)\nqs = range(0, 3)\nPs = range(0, 3)\nps = range(0, 3)\nD=1\nd=1\nparameters = product(ps, qs, Ps, Qs)\nparameters_list = list(parameters)\nlen(parameters_list)\n\n# Model Selection\nresults = []\nbest_aic = float(\"inf\")\nwarnings.filterwarnings('ignore')\nfor param in parameters_list:\n    try:\n#        model = SARIMAX(btc_month.close_box, order=(param[0], d, param[1]), seasonal_order=(param[2], D, param[3], 12)).fit(disp=-1)\n        model = SARIMAX(btc_month.close_box, order=(param[0], d, param[1]), seasonal_order=(param[2], D, param[3], 4)).fit(disp=-1)\n    except ValueError:\n        print('bad parameter combination:', param)\n        continue\n    aic = model.aic\n    if aic < best_aic:\n        best_model = model\n        best_aic = aic\n        best_param = param\n    results.append([param, model.aic])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Best Models\nresult_table = pd.DataFrame(results)\nresult_table.columns = ['parameters', 'aic']\nprint(result_table.sort_values(by = 'aic', ascending=True).head())\nprint(best_model.summary())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Dickey–Fuller test:: p=%f\" % adfuller(best_model.resid[13:])[1])\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"best_model.plot_diagnostics(figsize=(15, 12))\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# STL-decomposition\nplt.figure(figsize=(15,7))\nplt.subplot(211)\nbest_model.resid[13:].plot()\nplt.ylabel(u'Residuals')\nax = plt.subplot(212)\n#sm.graphics.tsa.plot_acf(best_model.resid[13:].values.squeeze(), lags=48, ax=ax)\nplot_acf(best_model.resid[13:].values.squeeze(), lags=12, ax=ax)\n\nprint(\"Dickey–Fuller test:: p=%f\" % adfuller(best_model.resid[13:])[1])\n\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"btc_month.columns","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"btc_month.index","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Prediction\nbtc_month2 = btc_month[['Close']]\ndate_list = [datetime(2018, 1, 31), datetime(2018, 4, 30), datetime(2018, 5, 31), datetime(2018, 6, 30)]\nfuture = pd.DataFrame(index=date_list, columns= btc_month.columns)\nbtc_month2 = pd.concat([btc_month2, future])\n\nbtc_month2['forecast'] = invboxcox(best_model.predict(start=0, end=75), lmbda)\n\nplt.figure(figsize=(15,7))\nbtc_month2.Close.plot()\nbtc_month2.forecast.plot(color='r', ls='--', label='forecast')\nplt.legend()\nplt.title('Bitcoin Monthly Close Forecast')\nplt.ylabel('USD')\nplt.savefig('bitcoin_monthly_forecast.png')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_forecasted = btc_month2.forecast\ny_truth = btc_month2['2015-01-01':'2017-01-01'].Close\n\n# Compute the root mean square error\nrmse = np.sqrt(((y_forecasted - y_truth) ** 2).mean())\nprint('Mean Squared Error: {}'.format(round(rmse, 2)))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### LSTM\n","metadata":{}},{"cell_type":"code","source":"def get_features(df):\n    df_feat = df[['Count', 'Open', 'High', 'Low', 'Close', 'Volume', 'VWAP', 'Target']].copy()\n    return df_feat","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:08:21.401741Z","iopub.execute_input":"2022-01-03T14:08:21.402283Z","iopub.status.idle":"2022-01-03T14:08:21.428059Z","shell.execute_reply.started":"2022-01-03T14:08:21.402185Z","shell.execute_reply":"2022-01-03T14:08:21.427323Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport gc\nimport traceback\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport gresearch_crypto\nimport tensorflow as tf\nfrom tensorflow import keras\nimport matplotlib.pyplot as plt\nimport pandas as pd, numpy as np\nfrom tensorflow.keras import layers\nimport tensorflow_probability as tfp\nimport tensorflow.keras.backend as K\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.preprocessing import RobustScaler\nfrom tensorflow.keras.preprocessing.sequence import TimeseriesGenerator\npd.set_option('display.max_columns', None)","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:08:21.429954Z","iopub.execute_input":"2022-01-03T14:08:21.430219Z","iopub.status.idle":"2022-01-03T14:08:28.548909Z","shell.execute_reply.started":"2022-01-03T14:08:21.430182Z","shell.execute_reply":"2022-01-03T14:08:28.548147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df = pd.read_csv('/kaggle/input/g-research-crypto-forecasting/train.csv')","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:08:28.550205Z","iopub.execute_input":"2022-01-03T14:08:28.550441Z","iopub.status.idle":"2022-01-03T14:09:21.837678Z","shell.execute_reply.started":"2022-01-03T14:08:28.550408Z","shell.execute_reply":"2022-01-03T14:09:21.836766Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Impute Missing Time Value\n\"\"\"def c_time_sub(asset_id, data=ctrain):\n    df=data[ctrain[\"Asset_ID\"]==asset_id].set_index(\"timestamp\")\n    df=df.reindex(range(df.index[0],df.index[-1]+60,60), method=\"pad\")\n    return df\"\"\"\n\ntrain_df.index=train_df['timestamp']\ntrain = train_df.sort_index()\nind = train.index.unique()\ndef reindex(df):\n    df = df.reindex(range(ind[0],ind[-1]+60,60),method='nearest')\n    df = df.fillna(method=\"ffill\").fillna(method=\"bfill\")\n    return df\ntrain = train.groupby('Asset_ID').apply(reindex).reset_index(0, drop=True).sort_index()\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:09:21.839941Z","iopub.execute_input":"2022-01-03T14:09:21.840196Z","iopub.status.idle":"2022-01-03T14:09:43.677888Z","shell.execute_reply.started":"2022-01-03T14:09:21.840162Z","shell.execute_reply":"2022-01-03T14:09:43.677133Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = train.reset_index(drop = True)\ntrain[['Count', 'Open', 'High', 'Low', 'Close', 'Volume', 'VWAP', 'Target']] = train[['Count', 'Open', 'High', 'Low', 'Close', 'Volume', 'VWAP', 'Target']].astype(np.float32)\ntrain['Target'] = train['Target'].fillna(0)\n\nVWAP_max = np.max(train[np.isfinite(train.VWAP)].VWAP)\nVWAP_min = np.min(train[np.isfinite(train.VWAP)].VWAP)\ntrain['VWAP'] = np.nan_to_num(train.VWAP, posinf=VWAP_max, neginf=VWAP_min)\ndf = train[['Asset_ID', 'Target']].copy()\n\ntimes = dict((t,i) for i,t in enumerate(df.index.unique()))\n\ndf['id'] = df.index.map(times)\ndf['id'] = df['id'].astype(str) + '_' + df['Asset_ID'].astype(str)\nids = df.id.copy()\ndel df","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:09:43.679426Z","iopub.execute_input":"2022-01-03T14:09:43.679957Z","iopub.status.idle":"2022-01-03T14:11:38.149329Z","shell.execute_reply.started":"2022-01-03T14:09:43.679917Z","shell.execute_reply":"2022-01-03T14:11:38.144275Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = get_features(train)","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:11:38.152626Z","iopub.execute_input":"2022-01-03T14:11:38.154122Z","iopub.status.idle":"2022-01-03T14:11:38.847372Z","shell.execute_reply.started":"2022-01-03T14:11:38.154078Z","shell.execute_reply":"2022-01-03T14:11:38.846587Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:11:38.848755Z","iopub.execute_input":"2022-01-03T14:11:38.849024Z","iopub.status.idle":"2022-01-03T14:11:38.874662Z","shell.execute_reply.started":"2022-01-03T14:11:38.84899Z","shell.execute_reply":"2022-01-03T14:11:38.873803Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_features = [i for i in train.columns if i not in ['Target', 'date', 'timestamp', 'Asset_ID', 'groups']]\ntrain.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:11:38.876293Z","iopub.execute_input":"2022-01-03T14:11:38.876583Z","iopub.status.idle":"2022-01-03T14:11:38.891237Z","shell.execute_reply.started":"2022-01-03T14:11:38.876545Z","shell.execute_reply":"2022-01-03T14:11:38.890512Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_features","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:11:38.892468Z","iopub.execute_input":"2022-01-03T14:11:38.892926Z","iopub.status.idle":"2022-01-03T14:11:38.903385Z","shell.execute_reply.started":"2022-01-03T14:11:38.89288Z","shell.execute_reply":"2022-01-03T14:11:38.902496Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"train['Target'] = train['Target'].fillna(0)\n","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:11:38.906396Z","iopub.execute_input":"2022-01-03T14:11:38.906917Z","iopub.status.idle":"2022-01-03T14:11:39.047065Z","shell.execute_reply.started":"2022-01-03T14:11:38.906875Z","shell.execute_reply":"2022-01-03T14:11:39.04628Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train['Target'].head()","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:03:13.740617Z","iopub.status.idle":"2022-01-03T14:03:13.741658Z","shell.execute_reply.started":"2022-01-03T14:03:13.741401Z","shell.execute_reply":"2022-01-03T14:03:13.741427Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DEVICE = \"TPU\" #or \"GPU\"\n\nSEED = 42\n\nEPOCHS = 10\nDEBUG = True\nN_ASSETS = 1\nWINDOW_SIZE = 15\nBATCH_SIZE = 1024\nPCT_VALIDATION = 10 # last 10% of the data are used as validation set","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:12:00.262962Z","iopub.execute_input":"2022-01-03T14:12:00.263449Z","iopub.status.idle":"2022-01-03T14:12:00.267968Z","shell.execute_reply.started":"2022-01-03T14:12:00.263411Z","shell.execute_reply":"2022-01-03T14:12:00.26706Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = train_df.copy()\n\ntargets = train['Target'].to_numpy().reshape(-1, N_ASSETS)\nfeatures = train.columns.drop(['Asset_ID', 'Target','timestamp'])\ntrain = train[features]\ntrain = train.values\ntrain = train.reshape(-1, N_ASSETS, train.shape[-1])","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:12:00.787531Z","iopub.execute_input":"2022-01-03T14:12:00.788112Z","iopub.status.idle":"2022-01-03T14:12:01.944651Z","shell.execute_reply.started":"2022-01-03T14:12:00.788071Z","shell.execute_reply":"2022-01-03T14:12:01.943908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Dataset creation\nSamples with a duration of WINDOW_SIZE records (minutes) will be formed from the train array. Each sample has a target vector corresponding to the final index if WINDOW_SIZE record.","metadata":{}},{"cell_type":"code","source":"class sample_generator(keras.utils.Sequence):\n    def __init__(self, x_set, y_set, batch_size, length):\n        self.x, self.y = x_set, y_set\n        self.batch_size = batch_size\n        self.length = length\n        self.size = len(x_set)\n    def __len__(self): return int(np.ceil(len(self.x) / float(self.batch_size)))\n    def __getitem__(self, idx):\n        batch_x=[]\n        batch_y=[]\n        for i in range(self.batch_size):\n            start_ind = self.batch_size*idx + i\n            end_ind = start_ind + self.length \n            if end_ind <= self.size:\n                batch_x.append(self.x[start_ind : end_ind])\n                batch_y.append(self.y[end_ind -1])\n        return np.array(batch_x), np.array(batch_y)","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:12:06.216332Z","iopub.execute_input":"2022-01-03T14:12:06.216834Z","iopub.status.idle":"2022-01-03T14:12:06.23936Z","shell.execute_reply.started":"2022-01-03T14:12:06.21677Z","shell.execute_reply":"2022-01-03T14:12:06.237091Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train, X_test = train[:-len(train)//PCT_VALIDATION], train[-len(train)//PCT_VALIDATION:]\ny_train, y_test = targets[:-len(train)//PCT_VALIDATION], targets[-len(train)//PCT_VALIDATION:]","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:12:08.607638Z","iopub.execute_input":"2022-01-03T14:12:08.607915Z","iopub.status.idle":"2022-01-03T14:12:08.612293Z","shell.execute_reply.started":"2022-01-03T14:12:08.607882Z","shell.execute_reply":"2022-01-03T14:12:08.611491Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_generator = sample_generator(X_train, y_train, length = WINDOW_SIZE, batch_size = BATCH_SIZE)\nval_generator = sample_generator(X_test, y_test, length = WINDOW_SIZE, batch_size = BATCH_SIZE)\nprint(f'Sample shape: {train_generator[0][0].shape}')\nprint(f'Target shape: {train_generator[0][1].shape}')","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:12:09.340211Z","iopub.execute_input":"2022-01-03T14:12:09.340704Z","iopub.status.idle":"2022-01-03T14:12:09.35733Z","shell.execute_reply.started":"2022-01-03T14:12:09.340666Z","shell.execute_reply":"2022-01-03T14:12:09.356538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training\nOur model will be trained for the number of FOLDS and EPOCHS you chose in the configuration above. Each fold the model with lowest validation loss will be saved and used to predict OOF and test. Adjust the variable VERBOSE. The variable VERBOSE=1 or 2 will display the training and validation loss for each epoch as text.","metadata":{}},{"cell_type":"code","source":"#Correlations for predicted and real\ndef MaxCorrelation(y_true,y_pred): \n    return -tf.math.abs(tfp.stats.correlation(y_pred,y_true, sample_axis=None, event_axis=None))\ndef Correlation(y_true,y_pred): \n    return tf.math.abs(tfp.stats.correlation(y_pred,y_true, sample_axis=None, event_axis=None))\n\n#Masked losses\ndef masked_mse(y_true, y_pred):\n    mask = tf.math.not_equal(y_true, 0.)\n    y_true_masked = tf.boolean_mask(y_true, mask)\n    y_pred_masked = tf.boolean_mask(y_pred, mask)\n    return tf.keras.losses.mean_squared_error(y_true = y_true_masked, y_pred = y_pred_masked)\n\n\ndef masked_mae(y_true, y_pred):\n    mask = tf.math.not_equal(y_true, 0.)\n    y_true_masked = tf.boolean_mask(y_true, mask)\n    y_pred_masked = tf.boolean_mask(y_pred, mask)\n    return tf.keras.losses.mean_absolute_error(y_true = y_true_masked, y_pred = y_pred_masked)\n\ndef masked_cosine(y_true, y_pred):\n    mask = tf.math.not_equal(y_true, 0.)\n    y_true_masked = tf.boolean_mask(y_true, mask)\n    y_pred_masked = tf.boolean_mask(y_pred, mask)\n    return tf.keras.losses.cosine_similarity(y_true_masked, y_pred_masked)","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:12:10.838116Z","iopub.execute_input":"2022-01-03T14:12:10.838924Z","iopub.status.idle":"2022-01-03T14:12:10.848071Z","shell.execute_reply.started":"2022-01-03T14:12:10.838851Z","shell.execute_reply":"2022-01-03T14:12:10.847281Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Long Short-Term Memory layer\ndef get_squence_model(x):\n    x = layers.LSTM(units=32, return_sequences=True)(x)\n    return x\n\n#Model\ndef get_model(n_assets = N_ASSETS):\n    x_input = keras.Input(shape=(train_generator[0][0].shape[1], n_assets, train_generator[0][0].shape[-1]))\n    branch_outputs = []\n    for i in range(n_assets):\n        a = layers.Lambda(lambda x: x[:,:, i])(x_input) # Slicing the ith asset:\n        a = layers.Masking(mask_value = 0., )(a)\n        a = get_squence_model(a)\n        a = layers.GlobalAvgPool1D()(a)\n        branch_outputs.append(a)\n    x = layers.Concatenate()(branch_outputs)\n    x = layers.Dense(units = 128)(x)\n    out = layers.Dense(units = n_assets)(x)\n    model = keras.Model(inputs=x_input, outputs=out)\n    model.compile(optimizer = tf.keras.optimizers.Adam(learning_rate=1e-3), loss = masked_cosine, metrics=[Correlation])\n    return model\n\nmodel = get_model()\nmodel.summary()","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:12:12.720709Z","iopub.execute_input":"2022-01-03T14:12:12.721253Z","iopub.status.idle":"2022-01-03T14:12:16.148843Z","shell.execute_reply.started":"2022-01-03T14:12:12.721212Z","shell.execute_reply":"2022-01-03T14:12:16.14807Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tf.keras.utils.plot_model(get_model(n_assets=1), show_shapes=True)","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:12:16.151057Z","iopub.execute_input":"2022-01-03T14:12:16.151473Z","iopub.status.idle":"2022-01-03T14:12:18.39268Z","shell.execute_reply.started":"2022-01-03T14:12:16.151435Z","shell.execute_reply":"2022-01-03T14:12:18.391741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(features)","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:12:18.394174Z","iopub.execute_input":"2022-01-03T14:12:18.394423Z","iopub.status.idle":"2022-01-03T14:12:18.400476Z","shell.execute_reply.started":"2022-01-03T14:12:18.394393Z","shell.execute_reply":"2022-01-03T14:12:18.399267Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tf.random.set_seed(0)\nestop = keras.callbacks.EarlyStopping(monitor = 'val_loss', patience = 7, verbose = 0, mode = 'min',restore_best_weights = True)\nscheduler = keras.optimizers.schedules.ExponentialDecay(1e-3, (0.5 * len(X_train) / BATCH_SIZE), 1e-3)\nlr = keras.callbacks.LearningRateScheduler(scheduler, verbose = 1)\n\nhistory = model.fit(train_generator, validation_data = (val_generator), epochs = EPOCHS, callbacks = [lr, estop])\n","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:13:36.382812Z","iopub.execute_input":"2022-01-03T14:13:36.383557Z","iopub.status.idle":"2022-01-03T14:45:32.290069Z","shell.execute_reply.started":"2022-01-03T14:13:36.383518Z","shell.execute_reply":"2022-01-03T14:45:32.289314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(1, 2, figsize=(16, 8))\nhistories = pd.DataFrame(history.history)\nepochs = list(range(1,len(histories)+1))\nloss = histories['loss']\nval_loss = histories['val_loss']\nCorrelation = histories['Correlation']\nval_Correlation = histories['val_Correlation']\nax[0].plot(epochs, loss, label = 'Train Loss')\nax[0].plot(epochs, val_loss, label = 'Val Loss')\nax[0].set_title('Losses')\nax[0].set_xlabel('Epoch')\nax[0].legend(loc='upper right')\nax[1].plot(epochs, Correlation, label = 'Train Correlation')\nax[1].plot(epochs, val_Correlation, label = 'Val Correlation')\nax[1].set_title('Correlations')\nax[1].set_xlabel('Epoch')\nax[1].legend(loc='upper right')\nfig.show()\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:46:50.119045Z","iopub.execute_input":"2022-01-03T14:46:50.119799Z","iopub.status.idle":"2022-01-03T14:46:50.722427Z","shell.execute_reply.started":"2022-01-03T14:46:50.119759Z","shell.execute_reply":"2022-01-03T14:46:50.721731Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Submit to Kaggle","metadata":{}},{"cell_type":"code","source":"def get_features(df):\n    df_feat = df[['Count', 'Open', 'High', 'Low', 'Close', 'Volume', 'VWAP']].copy()\n    return df_feat","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:47:46.994671Z","iopub.execute_input":"2022-01-03T14:47:46.994974Z","iopub.status.idle":"2022-01-03T14:47:46.999375Z","shell.execute_reply.started":"2022-01-03T14:47:46.99494Z","shell.execute_reply":"2022-01-03T14:47:46.998379Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"env = gresearch_crypto.make_env()\niter_test = env.iter_test()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for (test_df, sample_prediction_df) in iter_test:\n    test_df = get_features(test_df)\n    test_data = test_df.merge(example, how='outer', on='Asset_ID').sort_values('asset_order')\n    test = np.array(test_data[features].fillna(0))\n    test = test.reshape(-1, 1, N_ASSETS, test.shape[-1])\n    test_sample = np.hstack([test_sample, test])[:,-1 * WINDOW_SIZE:]\n    y_pred = model.predict(test_sample).squeeze().reshape(-1, 1).squeeze()\n    test_data['Target'] = y_pred\n    for _, row in test_df.iterrows():\n        try: sample_prediction_df.loc[sample_prediction_df['row_id'] == row['row_id'], 'Target'] = test_data.loc[test_data['row_id'] == row['row_id'], 'Target'].item()\n        except: sample_prediction_df.loc[sample_prediction_df['row_id'] == row['row_id'], 'Target'] = 0\n    env.predict(sample_prediction_df)","metadata":{"execution":{"iopub.status.busy":"2022-01-03T14:48:04.616977Z","iopub.execute_input":"2022-01-03T14:48:04.617564Z","iopub.status.idle":"2022-01-03T14:48:04.657987Z","shell.execute_reply.started":"2022-01-03T14:48:04.617511Z","shell.execute_reply":"2022-01-03T14:48:04.656885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#library\nimport os\nfrom datetime import datetime\nimport time\nfrom sklearn.preprocessing import StandardScaler\nimport plotly.express as px\nimport plotly.graph_objs as go\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport math\n\nimport statsmodels.api as sm\nfrom statsmodels.tsa.ar_model import AutoReg\nfrom statsmodels.tsa.arima_model import ARMA\nfrom statsmodels.tsa.stattools import adfuller\nfrom statsmodels.graphics.tsaplots import plot_acf\nfrom statsmodels.graphics.tsaplots import plot_pacf\nfrom statsmodels.tsa.arima_model import ARIMA\nfrom statsmodels.tsa.statespace.sarimax import SARIMAX\nfrom statsmodels.tsa.seasonal import seasonal_decompose\nfrom scipy import stats\nfrom itertools import product\n\n\nfrom statsmodels.graphics.tsaplots import plot_acf, plot_pacf\n\n\n\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n        \n        \n        \ninfo = pd.read_csv(\"/kaggle/input/g-research-crypto-forecasting/asset_details.csv\")\nctrain = pd.read_csv(\"/kaggle/input/g-research-crypto-forecasting/train.csv\")\n","metadata":{"execution":{"iopub.status.busy":"2022-01-03T13:24:53.547735Z","iopub.execute_input":"2022-01-03T13:24:53.548112Z","iopub.status.idle":"2022-01-03T13:25:52.174912Z","shell.execute_reply.started":"2022-01-03T13:24:53.548002Z","shell.execute_reply":"2022-01-03T13:25:52.17376Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"info = pd.read_csv(\"/kaggle/input/g-research-crypto-forecasting/asset_details.csv\")\nctrain = pd.read_csv(\"/kaggle/input/g-research-crypto-forecasting/train.csv\")\n","metadata":{"execution":{"iopub.status.busy":"2022-01-03T13:25:52.177091Z","iopub.execute_input":"2022-01-03T13:25:52.177454Z","iopub.status.idle":"2022-01-03T13:26:23.7793Z","shell.execute_reply.started":"2022-01-03T13:25:52.177407Z","shell.execute_reply":"2022-01-03T13:26:23.77819Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def preprocess_data(df,asset_id =1):    \n    df=df[df[\"Asset_ID\"]==asset_id].set_index(\"timestamp\")\n    df=df.reindex(range(df.index[0],df.index[-1]+60,60), method=\"pad\")\n    df.index = pd.to_datetime(df.index, unit='s')\n\n    # CALCULATE PRICE RETURNS AS DAILY PERCENTAGE CHANGE USING pct_change()\n    df['returns'] = 100 * df.Close.pct_change().dropna()\n\n    # CALCULATE LOG RETURNS BASED ON ABOVE FORMULA\n    df['log_returns'] = np.log(df.Close/df.Close.shift(1))\n\n    # DROPPING THE 1ST ROW OF DATA \n    # BECAUSE I SHIFTED IT FORWARD TO CALCULATE RETURNS/LOG RETURNS\n    df.dropna(inplace=True)\n    \n    #The Open/High/Low/Close prices are usually very similar and highly correlated to each other. Therefore, instead of keeping all of them in the dataset, I would add 2 more features:\n    df['HL_sprd'] = np.log((df.High - df.Low) / df.Close)\n    df['CO_sprd'] = (df.Close - df.Open) / df.Open\n    df['Volume'] = np.log(df.Volume)\n    features = ['HL_sprd', 'CO_sprd', 'Volume','log_returns']\n    x_train_processed = df[features]     \n    y_train_processed = df['Target']\n    print()\n    return df,x_train_processed,y_train_processed","metadata":{"execution":{"iopub.status.busy":"2022-01-03T13:26:50.239614Z","iopub.execute_input":"2022-01-03T13:26:50.240079Z","iopub.status.idle":"2022-01-03T13:26:50.252357Z","shell.execute_reply.started":"2022-01-03T13:26:50.240035Z","shell.execute_reply":"2022-01-03T13:26:50.251254Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df, X_btc, y_btc  = preprocess_data(ctrain,1)","metadata":{"execution":{"iopub.status.busy":"2022-01-03T13:26:50.958242Z","iopub.execute_input":"2022-01-03T13:26:50.959082Z","iopub.status.idle":"2022-01-03T13:26:52.788039Z","shell.execute_reply.started":"2022-01-03T13:26:50.959045Z","shell.execute_reply":"2022-01-03T13:26:52.786911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_btc.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-03T13:26:55.018134Z","iopub.execute_input":"2022-01-03T13:26:55.018882Z","iopub.status.idle":"2022-01-03T13:26:55.038228Z","shell.execute_reply.started":"2022-01-03T13:26:55.01883Z","shell.execute_reply":"2022-01-03T13:26:55.037281Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_btc.head()","metadata":{"execution":{"iopub.status.busy":"2022-01-03T13:26:57.291519Z","iopub.execute_input":"2022-01-03T13:26:57.291842Z","iopub.status.idle":"2022-01-03T13:26:57.301301Z","shell.execute_reply.started":"2022-01-03T13:26:57.291807Z","shell.execute_reply":"2022-01-03T13:26:57.300175Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In statistics and econometrics, an augmented Dickey–Fuller test (ADF) tests the null hypothesis that a unit root is present in a time series sample.","metadata":{}},{"cell_type":"code","source":"btc = ctrain[ctrain[\"Asset_ID\"]==1].set_index(\"timestamp\")","metadata":{"execution":{"iopub.status.busy":"2022-01-03T13:31:11.180346Z","iopub.execute_input":"2022-01-03T13:31:11.181277Z","iopub.status.idle":"2022-01-03T13:31:11.48436Z","shell.execute_reply.started":"2022-01-03T13:31:11.181209Z","shell.execute_reply":"2022-01-03T13:31:11.483306Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"    # CALCULATE PRICE RETURNS AS DAILY PERCENTAGE CHANGE USING pct_change()\nbtc['returns'] = 100 * btc.Close.pct_change().dropna()\n\n    # CALCULATE LOG RETURNS BASED ON ABOVE FORMULA\nbtc['log_returns'] = np.log(btc.Close/btc.Close.shift(1))","metadata":{"execution":{"iopub.status.busy":"2022-01-03T13:32:13.812359Z","iopub.execute_input":"2022-01-03T13:32:13.812679Z","iopub.status.idle":"2022-01-03T13:32:14.429134Z","shell.execute_reply.started":"2022-01-03T13:32:13.812644Z","shell.execute_reply":"2022-01-03T13:32:14.428152Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"    # CALCULATE PRICE RETURNS AS DAILY PERCENTAGE CHANGE USING pct_change()\ndf['returns'] = 100 * df.Close.pct_change().dropna()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from statsmodels.tsa.stattools import adfuller\n# LOG RETURNS\nadfuller_results = adfuller(btc.log_returns.dropna())\n\nprint(f'ADF Statistic: {adfuller_results[0]}')\nprint(f'p-value: {adfuller_results[1]}')\nprint('Critical Values:')\nfor key, value in adfuller_results[4].items():\n    print(f'{key}: {value:.4f}')","metadata":{"execution":{"iopub.status.busy":"2022-01-03T13:32:15.248833Z","iopub.execute_input":"2022-01-03T13:32:15.249482Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# RETURNS\nadfuller_results = adfuller(df.returns.dropna())\n\nprint(f'ADF Statistic: {adfuller_results[0]}')\nprint(f'p-value: {adfuller_results[1]}')\nprint('Critical Values:')\nfor key, value in adfuller_results[4].items():\n    print(f'{key}: {value:.4f}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Split","metadata":{}},{"cell_type":"code","source":"#SPLIT\ndef custom_split(X_btc,y_btc):   \n    PCT_VALIDATION = 10 # last 10% of the data are used as validation set\n\n\n    X_train, X_test = X_btc[:-len(X_btc)//PCT_VALIDATION], X_btc[-len(X_btc)//PCT_VALIDATION:]\n    y_train, y_test = y_btc[:-len(y_btc)//PCT_VALIDATION], y_btc[-len(y_btc)//PCT_VALIDATION:]\n    return X_train, X_test, y_train, y_test\n\nX_train, X_test, y_train, y_test = custom_split(X_btc,y_btc)   ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Scaler","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"##think about it later\n\n\nfrom sklearn.preprocessing import MinMaxScaler\n\n\n\ndef transform_volatility_to_scaler(scaler, tf_series):\n    '''\n    Transform a series to a fitted scaler\n    '''\n    idx = tf_series.index\n    output = pd.Series(scaler.transform(tf_series.values.reshape(-1,1))[:,0],\n                       index=idx)\n    return output\n\n# INITIATE SCALER\nscaler = MinMaxScaler()\n\n# FIT SCALER TO CURRENT VOLATILITY IN TRAINING SET\nscaler = scaler.fit(X_train.values.reshape(-1,1))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\"# TRANSFORM TRAINING CURRENT & FUTURE VOLATILITIES \nx_train_scaled = transform_volatility_to_scaler(scaler_vol, x_train)\ny_train_scaled = transform_volatility_to_scaler(scaler_vol, y_train)\n\n# TRANSFORMING VALIDATION CURRENT & FUTURE VOLATILITIES\nx_val_scaled = transform_volatility_to_scaler(scaler_vol, x_val)\ny_val_scaled = transform_volatility_to_scaler(scaler_vol, y_val)\n\n# TRANSFORMING TEST CURRENT & FUTURE VOLATILITIES\nx_test_scaled = transform_volatility_to_scaler(scaler_vol, x_test)\ny_test_scaled = transform_volatility_to_scaler(scaler_vol, y_test)\"\"\"\"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# DEFINE ROOT MEAN SQUARED PERCENTAGE ERROR FUNCTION\ndef RMSPE(y_true, y_pred):\n    \"\"\"\n    Compute Root Mean Squared Percentage Error between 2 arrays\n    \"\"\"\n    output = np.sqrt(np.mean(np.square((y_true - y_pred) / y_true)))\n    return output\n\n\n\n# DEFINE ROOT MEAN SQUARED ERROR FUNCTION\ndef RMSE(y_true, y_pred):\n    \"\"\"\n    Compute Root Mean Squared Error between 2 arrays\n    \"\"\"\n    output = np.sqrt(mse(y_true, y_pred))\n    return output\n\n# INITIATE A DATAFRAME FOR MODEL PERFORMANCE TRACKING & COMPARISON \nperf_df = pd.DataFrame(columns=['Model', 'Validation RMSPE', 'Validation RMSE'])\n\n# A FUNCTION THAT LOGS MODEL NAME, RMSE AND RMPSE INTO perf_df\n# FOR EASY COMPARISON LATER\ndef log_perf(y_true, y_pred, model_name):\n    perf_df.loc[len(perf_df.index)] = [model_name, \n                                       RMSPE(y_true, y_pred), \n                                       RMSE(y_true, y_pred)]\n    return perf_df\n\n\n# PLOTTING MODEL PREDICTIONS VS. TARGET VALUES\ndef viz_model(y_true, y_pred, model_name):\n    sns.set_context(\"paper\", font_scale=1.7)\n    plt.rcParams[\"axes.grid\"] = False\n\n    with sns.axes_style(\"whitegrid\"):\n        plt.figure(figsize=(18,7))\n        plt.plot(x_val_scaled, color='gray',  ls=':',\n                label=f\"Scaled Current Daily Volatility\")\n        \n        plt.plot(y_true, color='blue', lw=2, \n                label=f\"Target Volatility\")\n        plt.plot(y_pred, color='orange', lw=2.5,\n                label=f'Forecasted Volatility')\n        \n        # plt.plot(lr_val, color='gray', alpha=0.4,\n        #         label='Daily Log Returns')\n\n        plt.title(f'{model_name} \\non Validation Data')\n        plt.legend(loc='best', frameon=True)\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 4.7. Multivariate 2-Layered Bidirectional LSTM with Dropout\n\nMy first Multivariate LSTM model would be relatively simple with 2 hidden layers of Bidirectional LSTM (quite similar to the Univariate structure above). However, I will also include Dropout layers in between. Having more features means the model would be more prone to **overfitting**, and Dropout layers can help reduce that.","metadata":{}},{"cell_type":"code","source":"import tensorflow as tf\nfrom tensorflow.keras.callbacks import ModelCheckpoint, EarlyStopping\nfrom tensorflow.keras import backend as K\nfrom decimal import Decimal\n# CONVERT A SERIES TO MATRIX DATASET\n\nseed = 20\n\n\ndef windowed_dataset(x_series, y_series, lookback_window):\n    dataX, dataY = [], []\n    for i in range((lookback_window-1), len(x_series)):\n        start_idx = x_series.index[i-lookback_window+1]\n        end_idx = x_series.index[i]\n        a = x_series[start_idx:end_idx].values\n        dataX.append(a)\n        dataY.append(y_series[end_idx])\n\n    return np.array(dataX), np.array(dataY)\n\n\ntf.keras.backend.clear_session()\n\n# SET SEED FOR REPRODUCIBILITY\nnp.random.seed(seed)\n\nn_past = 30\nbatch_size = 64\nn_dims = X_btc.shape[1]\n\nmat_X_train, mat_y_train = windowed_dataset(X_train, y_train, n_past)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### LSTM\nversion on git https://github.com/ayman-mezghani/ml4f-project/blob/master/g-research.ipynb","metadata":{}},{"cell_type":"code","source":"import os\nimport gc\nimport traceback\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport gresearch_crypto\nimport tensorflow as tf\nfrom tensorflow import keras\nimport matplotlib.pyplot as plt\nimport pandas as pd, numpy as np\nfrom tensorflow.keras import layers\nimport tensorflow_probability as tfp\nimport tensorflow.keras.backend as K\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.preprocessing import RobustScaler\nfrom tensorflow.keras.preprocessing.sequence import TimeseriesGenerator\npd.set_option('display.max_columns', None)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"        \nasset_details = pd.read_csv(\"/kaggle/input/g-research-crypto-forecasting/asset_details.csv\")\ntrain = pd.read_csv(\"/kaggle/input/g-research-crypto-forecasting/train.csv\")\ntest = pd.read_csv(\"/kaggle/input/g-research-crypto-forecasting/example_test.csv\")\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n","metadata":{"execution":{"iopub.status.busy":"2021-12-28T07:15:36.872515Z","iopub.execute_input":"2021-12-28T07:15:36.872718Z","iopub.status.idle":"2021-12-28T07:15:38.490051Z","shell.execute_reply.started":"2021-12-28T07:15:36.872692Z","shell.execute_reply":"2021-12-28T07:15:38.489197Z"},"jupyter":{"outputs_hidden":true}}},{"cell_type":"markdown","source":"Feature eng","metadata":{}},{"cell_type":"code","source":"#Params\nDEBUG = True\nEPOCHS = 200\nN_ASSETS = 14\nWINDOW_SIZE = 15\nBATCH_SIZE = 1024\nPCT_VALIDATION = 10 # last 10% of the data are used as validation set","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gc \n\ndef time_change(train):\n    # Convert timestamp\n    train['timestamp'] = pd.to_datetime(train['timestamp'], unit='s')\n\n    train.index = train.timestamp\n    train.sort_index(inplace=True)\n    train['is_real'] = True\n    ind = train.index.unique()\n\n    def reindex(df):\n        res = df.reindex(pd.date_range(ind.min(), ind.max(), freq='min'))\n        res['is_real'].fillna(False, inplace=True)\n        res['timestamp'] = res.index\n        res = res.fillna(method=\"ffill\").fillna(method=\"bfill\")\n        return res\n    train = train.groupby('Asset_ID').apply(reindex).reset_index(0, drop=True).sort_values(by=['timestamp', 'Asset_ID'])\n    gc.collect()\n    return train\n\ndef vmap_change(train):\n    #VWAP  - The average price of the asset over the time interval, weighted by volume.\n    VWAP_max = np.max(train[np.isfinite(train.VWAP)].VWAP)\n    VWAP_min = np.min(train[np.isfinite(train.VWAP)].VWAP)\n    train['VWAP'] = np.nan_to_num(train.VWAP, posinf=VWAP_max, neginf=VWAP_min)\n    \n    if DEBUG:\n        train = train[-500010:]\n    return train\n\ndef target_train_shape_change(train):\n    targets = train['Target'].to_numpy().reshape(-1, N_ASSETS)\n    feature_cols = train.columns.drop(['Asset_ID', 'Target', 'timestamp', 'is_real'])\n\n    #non-real data features are set to 0 which are then masked by the model\n    for col in feature_cols:\n        train[col] = train[col] * train.is_real\n        \n    train_data = train[feature_cols].values\n    train_data = train_data.reshape(-1, N_ASSETS, train_data.shape[-1])\n    return targets, train_data\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.columns","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_features(df):\n    df_feat = df[['Count', 'Open', 'High', 'Low', 'Close', 'Volume', 'VWAP']].copy()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"targets, train_data = target_train_shape_change(vmap_change(time_change(train)))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Samples with a duration of WINDOW_SIZE records (minutes) will be formed from the train array. \n#Each sample has a target vector corresponding to the final index if WINDOW_SIZE record.\n\nclass sample_generator(keras.utils.Sequence):\n    def __init__(self, x_set, y_set, batch_size, length):\n        self.x, self.y = x_set, y_set\n        self.batch_size = batch_size\n        self.length = length\n        self.size = len(x_set)\n    def __len__(self): return int(np.ceil(len(self.x) / float(self.batch_size)))\n    def __getitem__(self, idx):\n        batch_x = []\n        batch_y = []\n        for i in range(self.batch_size):\n            start_ind = self.batch_size * idx + i\n            end_ind = start_ind + self.length \n            if end_ind <= self.size:\n                batch_x.append(self.x[start_ind : end_ind])\n                batch_y.append(self.y[end_ind -1])\n        return np.array(batch_x), np.array(batch_y)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train, X_test = train_data[:-len(train_data)//PCT_VALIDATION], train_data[-len(train_data)//PCT_VALIDATION:]\ny_train, y_test = targets[:-len(train_data)//PCT_VALIDATION], targets[-len(train_data)//PCT_VALIDATION:]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_generator = sample_generator(X_train, y_train, length=WINDOW_SIZE, batch_size=BATCH_SIZE)\nval_generator = sample_generator(X_test, y_test, length=WINDOW_SIZE, batch_size=BATCH_SIZE)\nprint(f'Sample shape: {train_generator[0][0].shape}')\nprint(f'Target shape: {train_generator[0][1].shape}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Metrics and Loss functions","metadata":{}},{"cell_type":"code","source":"#Correlations for predicted and real\ndef MaxCorrelation(y_true,y_pred): \n    return -tf.math.abs(tfp.stats.correlation(y_pred, y_true, sample_axis=None, event_axis=None))\n\ndef Correlation(y_true,y_pred): \n    return tf.math.abs(tfp.stats.correlation(y_pred, y_true, sample_axis=None, event_axis=None))\n\n#Masked losses\ndef masked_mse(y_true, y_pred):\n    mask = tf.math.not_equal(y_true, 0.)\n    y_true_masked = tf.boolean_mask(y_true, mask)\n    y_pred_masked = tf.boolean_mask(y_pred, mask)\n    return tf.keras.losses.mean_squared_error(y_true=y_true_masked, y_pred=y_pred_masked)\n\ndef masked_mae(y_true, y_pred):\n    mask = tf.math.not_equal(y_true, 0.)\n    y_true_masked = tf.boolean_mask(y_true, mask)\n    y_pred_masked = tf.boolean_mask(y_pred, mask)\n    return tf.keras.losses.mean_absolute_error(y_true=y_true_masked, y_pred=y_pred_masked)\n\ndef masked_cosine(y_true, y_pred):\n    mask = tf.math.not_equal(y_true, 0.)\n    y_true_masked = tf.boolean_mask(y_true, mask)\n    y_pred_masked = tf.boolean_mask(y_pred, mask)\n    return tf.keras.losses.cosine_similarity(y_true_masked, y_pred_masked)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Multivariate 1-Layered LSTM with Gloobal Average Pool","metadata":{}},{"cell_type":"code","source":"train_generator[0][0].shape[1]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Model\ndef get_model_1_layer(n_assets=N_ASSETS):\n    x_input = keras.Input(shape=(train_generator[0][0].shape[1], n_assets, train_generator[0][0].shape[-1]))\n    branch_outputs = []\n    \n    for i in range(n_assets):\n        a = layers.Lambda(lambda x: x[:,:, i])(x_input) # Slicing the ith asset:\n        a = layers.Masking(mask_value = 0., )(a)\n        a = layers.LSTM(units=32, return_sequences=True)(a)\n        a = layers.GlobalAvgPool1D()(a)\n        branch_outputs.append(a)\n        \n    x = layers.Concatenate()(branch_outputs)\n    x = layers.Dense(units = 128)(x)\n    out = layers.Dense(units = n_assets)(x)\n    model = keras.Model(inputs=x_input, outputs=out)\n    model.compile(optimizer = tf.keras.optimizers.Adam(learning_rate=1e-3), loss=masked_mse, metrics=[Correlation])\n    return model\n\nmodel_1_layer = get_model_1_layer()\nmodel_1_layer.summary()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tf.keras.utils.plot_model(get_model_1_layer(n_assets=1), show_shapes=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tf.random.set_seed(0)\n\n#Stop training when a monitored metric has stopped improving\nestop = keras.callbacks.EarlyStopping(monitor='val_loss', patience=7, verbose=0, mode='min', restore_best_weights=True)\n\nscheduler = keras.optimizers.schedules.ExponentialDecay(1e-3, (0.5 * len(X_train) / BATCH_SIZE), 1e-3)\nlr = keras.callbacks.LearningRateScheduler(scheduler, verbose=1)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"history_1_layer = model_1_layer.fit(train_generator, validation_data=(val_generator), epochs=EPOCHS, callbacks=[lr, estop])\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_training_history(history):\n    fig, ax = plt.subplots(1, 2, figsize=(16, 8))\n    histories = pd.DataFrame(history.history)\n    epochs = list(range(1,len(histories)+1))\n    loss = histories['loss']\n    val_loss = histories['val_loss']\n    correlation = histories['Correlation']\n    val_correlation = histories['val_Correlation']\n    ax[0].plot(epochs, loss, label='Train Loss')\n    ax[0].plot(epochs, val_loss, label='Val Loss')\n    ax[0].set_title('Losses')\n    ax[0].set_xlabel('Epoch')\n    ax[0].legend(loc='upper right')\n    ax[1].plot(epochs, correlation, label='Train Correlation')\n    ax[1].plot(epochs, val_correlation, label='Val Correlation')\n    ax[1].set_title('Correlations')\n    ax[1].set_xlabel('Epoch')\n    ax[1].legend(loc='upper right')\n    fig.show()\ngc.collect()\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_training_history(history_1_layer)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions_1_layer = model_1_layer.predict(val_generator)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"perf_df = pd.DataFrame(columns=['Model', 'asset', 'corr','weights'])# A FUNCTION THAT LOGS MODEL NAME, RMSE AND RMPSE INTO perf_df\n\nprint('Asset:    Corr. coef.')\nprint('---------------------')\nfor i in range(N_ASSETS):\n    # drop first 14 values in the y_test, since they are absent in val_generator labels\n    y_true = np.squeeze(y_test[WINDOW_SIZE - 1:, i])\n    y_pred = np.squeeze(predictions_1_layer[:, i])\n    real_target_ind = np.argwhere(y_true != 0)\n    asset_name = asset_details[asset_details.Asset_ID == i]['Asset_Name'].item()\n    weights = asset_details[asset_details.Asset_ID == i]['Weight'].item()\n    correl = np.corrcoef(y_pred[real_target_ind].flatten(), y_true[real_target_ind].flatten())[0,1]\n    perf_df = perf_df.append([pd.Series(['1 layer LSTM',asset_name,correl,weights ],index = perf_df.columns)], ignore_index=True)\n    \n    print(f\"{asset_name}: {correl:.4f}\")\n    plt.plot(y_true, label='Target')\n    plt.plot(y_pred, label='Prediction')\n    plt.xlabel('Time')\n    plt.ylabel('Target')\n    plt.legend()\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"perf_df","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"perf_df = perf_df.append([pd.Series(['1 layer LSTM', 'WEIGTED', (perf_df['corr']*perf_df['weights']).sum(), 1 ],index = perf_df.columns)], ignore_index=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"perf_df.to_csv('Performance.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"perf_df","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import datetime\ntimestamp = datetime.datetime.now()\nmodel_1_layer.save(f\"model_1_layer{timestamp.strftime('%s')}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model_1_layer.save('model_1_layer.h5')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# load data that is just before the test set\nsup = pd.read_csv('/kaggle/input/g-research-crypto-forecasting/supplemental_train.csv').sort_values(by=['timestamp', 'Asset_ID'])\n\n# only the last 15 min are useful\nsup = sup[-WINDOW_SIZE * (N_ASSETS):]\n\n# prepare the sup data. Test data will be appended\ntest_sample = np.array(get_features(sup))\n#test_sample = np.array(sup[feature_cols])\ntest_sample = test_sample.reshape(-1, (N_ASSETS), test_sample.shape[-1])\ntest_sample = np.expand_dims(test_sample, axis=0)\n\n\nenv = gresearch_crypto.make_env()\niter_test = env.iter_test()\n\n\"\"\"\nenv = gresearch_crypto.make_env()\niter_test = env.iter_test()\n\"\"\"\nfor (test_df, sample_prediction_df) in iter_test:\n    # get the features\n    test_df = get_features(test_df)\n    \n    # always sort by time, then by Asset_ID\n    test_df = test_df.sort_values(by=['timestamp', 'Asset_ID'])\n    \n    # to map to the corresponding rows    \n    asset_id_row_id_map = {a_id: r_id for a_id, r_id in test_df[['Asset_ID', 'row_id']].values}\n    \n    # matrix of features\n    test = np.array(test_df[feature_cols].fillna(0))\n    \n    # reshaping (similar to the train)\n    test = test.reshape(-1, 1, N_ASSETS, test.shape[-1])\n    \n    # Stack the test data to previous data, keep only last 15 min \n    test_sample = np.hstack([test_sample, test])[:,-WINDOW_SIZE:]\n    \n    y_pred = model.predict(test_sample).squeeze().reshape(-1, 1).squeeze()\n    \n    for i, p in enumerate(y_pred):\n        sample_prediction_df.loc[sample_prediction_df['row_id'] == asset_id_row_id_map[i], 'Target'] = p\n    env.predict(sample_prediction_df)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Multivariate 2-Layered Bidirectional LSTM Model","metadata":{}},{"cell_type":"code","source":"#Model\ndef get_model(n_assets=N_ASSETS):\n    x_input = keras.Input(shape=(train_generator[0][0].shape[1], n_assets, train_generator[0][0].shape[-1]))\n    branch_outputs = []\n\n    for i in range(n_assets):\n        a = layers.Lambda(lambda x: x[:,:, i])(x_input) # Slicing the ith asset:\n        a = layers.Masking(mask_value = 0., )(a)\n        a = layers.Bidirectional(layers.LSTM(32, return_sequences=True))(a)\n        a = layers.Bidirectional(layers.LSTM(16))(a)\n        #a = layers.GlobalAvgPool1D()(a)\n        branch_outputs.append(a)\n    \n    x = layers.Concatenate()(branch_outputs)\n    x = layers.Dense(units = 128)(x)\n    out = layers.Dense(units = n_assets)(x)\n    model = keras.Model(inputs=x_input, outputs=out)\n    model.compile(optimizer = tf.keras.optimizers.Adam(learning_rate=1e-3), loss=masked_mse, metrics=[Correlation])\n    return model\n\nmodel = get_model()\nmodel.summary()\n\ntf.keras.utils.plot_model(get_model(n_assets=1), show_shapes=True)\n\n#print(features)\n\ntf.random.set_seed(0)\nestop = keras.callbacks.EarlyStopping(monitor='val_loss', patience=7, verbose=0, mode='min', restore_best_weights=True)\nscheduler = keras.optimizers.schedules.ExponentialDecay(1e-3, (0.5 * len(X_train) / BATCH_SIZE), 1e-3)\nlr = keras.callbacks.LearningRateScheduler(scheduler, verbose=1)\n\n\n#train and test\n\n\nhistory = model.fit(train_generator, validation_data=(val_generator), epochs=EPOCHS, callbacks=[lr, estop])\n\nplot_training_history(history)\n\npredictions = model.predict(val_generator)\n\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_training_history(history)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# DEFINE ROOT MEAN SQUARED ERROR FUNCTION\ndef RMSE(y_true, y_pred):\n    \"\"\"\n    Compute Root Mean Squared Error between 2 arrays\n    \"\"\"\n    output = np.sqrt(mse(y_true, y_pred))\n    return output\n# INITIATE A DATAFRAME FOR MODEL PERFORMANCE TRACKING & COMPARISON \nperf_df = pd.DataFrame(columns=['Model', 'Validation RMSPE', 'Validation RMSE'])\n\n# A FUNCTION THAT LOGS MODEL NAME, RMSE AND RMPSE INTO perf_df\n# FOR EASY COMPARISON LATER\ndef log_perf(y_true, y_pred, model_name):\n    perf_df.loc[len(perf_df.index)] = [model_name, \n                                       RMSPE(y_true, y_pred), \n                                       RMSE(y_true, y_pred)]\n    return perf_df\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_test","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions[:,0]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Asset:    Corr. coef.')\nprint('---------------------')\nfor i in range(N_ASSETS):\n    # drop first 14 values in the y_test, since they are absent in val_generator labels\n    y_true = np.squeeze(y_test[WINDOW_SIZE - 1:, i])\n    y_pred = np.squeeze(predictions[:, i])\n    real_target_ind = np.argwhere(y_true != 0)\n    asset_name = asset_details[asset_details.Asset_ID == i]['Asset_Name'].item()\n    print(f\"{asset_name}: {np.corrcoef(y_pred[real_target_ind].flatten(), y_true[real_target_ind].flatten())[0,1]:.4f}\")\n    plt.plot(y_true, label='Target')\n    plt.plot(y_pred, label='Prediction')\n    plt.xlabel('Time')\n    plt.ylabel('Target')\n    plt.legend()\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# load data that is just before the test set\nsup = pd.read_csv('/kaggle/input/g-research-crypto-forecasting/supplemental_train.csv').sort_values(by=['timestamp', 'Asset_ID'])\n\n# only the last 15 min are useful\nsup = sup[-WINDOW_SIZE * (N_ASSETS):]\n\n# prepare the sup data. Test data will be appended\ntest_sample = np.array(get_features(sup))\n#test_sample = np.array(sup[feature_cols])\ntest_sample = test_sample.reshape(-1, (N_ASSETS), test_sample.shape[-1])\ntest_sample = np.expand_dims(test_sample, axis=0)\n\n\nenv = gresearch_crypto.make_env()\niter_test = env.iter_test()\n\n\"\"\"\nenv = gresearch_crypto.make_env()\niter_test = env.iter_test()\n\"\"\"\nfor (test_df, sample_prediction_df) in iter_test:\n    # get the features\n    test_df = get_features(test_df)\n    \n    # always sort by time, then by Asset_ID\n    test_df = test_df.sort_values(by=['timestamp', 'Asset_ID'])\n    \n    # to map to the corresponding rows    \n    asset_id_row_id_map = {a_id: r_id for a_id, r_id in test_df[['Asset_ID', 'row_id']].values}\n    \n    # matrix of features\n    test = np.array(test_df[feature_cols].fillna(0))\n    \n    # reshaping (similar to the train)\n    test = test.reshape(-1, 1, N_ASSETS, test.shape[-1])\n    \n    # Stack the test data to previous data, keep only last 15 min \n    test_sample = np.hstack([test_sample, test])[:,-WINDOW_SIZE:]\n    \n    y_pred = model.predict(test_sample).squeeze().reshape(-1, 1).squeeze()\n    \n    for i, p in enumerate(y_pred):\n        sample_prediction_df.loc[sample_prediction_df['row_id'] == asset_id_row_id_map[i], 'Target'] = p\n    env.predict(sample_prediction_df)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_df.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"iter_test = env.iter_test()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}