{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## メモ\n\nラグやマーケットリターンなどを組み合わせて適当に特徴量作り、Ridge回帰でTargetを予測したnotebook。\n提出がタイムアウトになるので、提出部分の高速化が必要。\n\nコードは公開していないが、俺が運用している仮想通貨ボットで使っている特徴量 + ridge回帰(alphaは適当)のスコア(相関)は0.013程度だった。\nこのnotebookだと相関が0.04くらい。\n\n## Notes\n\nA notebook that predicts the target by ridge regression by creating some features by combining lag and market returns.\nSince the submission will time out, it is necessary to speed up the submission part.\n\nAlthough the code is not disclosed, the score (correlation) when using the features used in my crypto trading bot was about 0.013.\nWith this notebook, the correlation is about 0.04.\n\n## TODO\n\n- ffill missing bars\n- evaluation metrics (use weights. other metrics like sharpe, double sharpe)\n- model and feature improvement\n- robust cv (nested cv or bbc-cv) + hyper parameter tuning","metadata":{}},{"cell_type":"code","source":"import os\nimport lzma\n\nimport cloudpickle\nimport datatable as dt\nimport gresearch_crypto\nimport joblib\nimport lightgbm as lgb\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nfrom scipy.stats import pearsonr, spearmanr\nfrom sklearn.base import BaseEstimator, clone\nfrom sklearn.ensemble import VotingRegressor, BaggingRegressor\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import r2_score\nfrom sklearn.model_selection import cross_val_predict, cross_val_score\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.preprocessing import QuantileTransformer, StandardScaler","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-11-08T02:46:27.932720Z","iopub.execute_input":"2021-11-08T02:46:27.933210Z","iopub.status.idle":"2021-11-08T02:46:27.942525Z","shell.execute_reply.started":"2021-11-08T02:46:27.933155Z","shell.execute_reply":"2021-11-08T02:46:27.941149Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_feature_columns(df):\n    features = df.columns[df.columns.str.startswith('feature')]\n    return sorted(list(features))\n\ndef save_model(model, path):\n    data = cloudpickle.dumps(model)\n    data = lzma.compress(data)\n    with open(path, 'wb') as f:\n        f.write(data)\n        \ndef process_data(df, df_asset):\n    df = df.copy()\n    df['timestamp'] = pd.to_datetime(df['timestamp'], unit='s', utc=True)\n    df = df.rename(columns={\n        'Asset_ID': 'market',\n        'Open': 'op',\n        'High': 'hi',\n        'Low': 'lo',\n        'Close': 'cl',\n        'Volume': 'volume',\n        'VWAP': 'vwap',\n        'Count': 'trade_count',\n        'Target': 'target',\n    })\n    df = df.join(df_asset[['weight']], on='market', how='left')\n    df = df.set_index(['timestamp', 'market'])\n    return df\n\ndef sort_and_remove_duplicates(df):\n    df = df.sort_index(kind='mergesort')\n    # https://stackoverflow.com/questions/13035764/remove-rows-with-duplicate-indices-pandas-dataframe-and-timeseries\n    df = df.loc[~df.index.duplicated(keep='last')]\n    return df\n\ndef my_purge_kfold(n, n_splits=5, purge=3750 * 14):\n    idx = np.arange(n)\n    cv = []\n    for i in range(n_splits):\n        val_start = i * n // n_splits\n        val_end = (i + 1) * n // n_splits\n        val_idx = idx[val_start:val_end]\n        train_idx = idx[(idx < val_start - purge) | (val_end + purge <= idx)]\n        cv.append((\n            train_idx,\n            val_idx,\n        ))\n    return cv","metadata":{"execution":{"iopub.status.busy":"2021-11-08T02:46:27.945707Z","iopub.execute_input":"2021-11-08T02:46:27.946115Z","iopub.status.idle":"2021-11-08T02:46:27.963503Z","shell.execute_reply.started":"2021-11-08T02:46:27.946064Z","shell.execute_reply":"2021-11-08T02:46:27.962727Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# preprocess asset data\n\ndf = dt.fread('../input/g-research-crypto-forecasting/asset_details.csv').to_pandas()\ndf = df.rename(columns={\n    'Asset_ID': 'market',\n    'Weight': 'weight',\n    'Asset_Name': 'name',\n})\ndf = df.set_index('market')\ndf = df.sort_values('market')\ndf.to_pickle('/tmp/df_asset.pkl')\ndisplay(df)","metadata":{"execution":{"iopub.status.busy":"2021-11-08T02:46:27.964865Z","iopub.execute_input":"2021-11-08T02:46:27.965190Z","iopub.status.idle":"2021-11-08T02:46:28.007483Z","shell.execute_reply.started":"2021-11-08T02:46:27.965146Z","shell.execute_reply":"2021-11-08T02:46:28.006459Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# preprocess train data\n\n# supplemental_trainは提出後に増えるデータらしい\n# trainとsupplemental_trainをくっつけてtrainにすれば良さそう\ndf = pd.concat([\n    dt.fread('../input/g-research-crypto-forecasting/train.csv').to_pandas(),\n    dt.fread('../input/g-research-crypto-forecasting/supplemental_train.csv').to_pandas(),\n])\ndf_asset = pd.read_pickle('/tmp/df_asset.pkl')\ndf = process_data(df, df_asset)\ndf = sort_and_remove_duplicates(df)\ndf.to_pickle('/tmp/df.pkl')\ndisplay(df)","metadata":{"execution":{"iopub.status.busy":"2021-11-08T02:46:28.012294Z","iopub.execute_input":"2021-11-08T02:46:28.012599Z","iopub.status.idle":"2021-11-08T02:46:59.414982Z","shell.execute_reply.started":"2021-11-08T02:46:28.012567Z","shell.execute_reply":"2021-11-08T02:46:59.414328Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# check interval is 1 min\n\ndf = pd.read_pickle('/tmp/df.pkl')\ndf = df.reset_index()\ndf['timestamp'] = df['timestamp'].view(int) / 10 ** 9\ndf['interval'] = df['timestamp'] - df.groupby('market')['timestamp'].shift(1)\ndf = df.dropna()\nfor interval in range(60, 361, 60):\n    print('{} {}'.format(interval, np.mean(df['interval'] == interval)))","metadata":{"execution":{"iopub.status.busy":"2021-11-08T02:46:59.416136Z","iopub.execute_input":"2021-11-08T02:46:59.416510Z","iopub.status.idle":"2021-11-08T02:47:09.531577Z","shell.execute_reply.started":"2021-11-08T02:46:59.416478Z","shell.execute_reply":"2021-11-08T02:47:09.530421Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# calc features\n\ndef calc_features(df):\n    df = df.copy()\n    \n    df['ln_cl'] = np.log(df['cl'])\n    \n#     df['feature_upper_shadow'] = df['hi'] - np.maximum(df['op'], df['cl'])\n#     df['feature_lower_shadow'] = np.minimum(df['cl'], df['op']) - df['lo']\n\n    # shift is faster than diff\n    df['feature_cl_diff1'] = df['ln_cl'] - df.groupby('market')['ln_cl'].shift(15)\n    df['raw_return_causal'] = df['ln_cl'] - df.groupby('market')['ln_cl'].shift(15)\n    \n    inv_weight_sum = 1.0 / df.groupby('timestamp')['weight'].transform('sum')\n    \n    df['market_return_causal'] = (df['raw_return_causal'] * df['weight']).groupby('timestamp').transform('sum') * inv_weight_sum\n    \n    df['beta_causal'] = (\n        (df['raw_return_causal'] * df['market_return_causal']).groupby('market').transform(lambda x: x.rolling(3750, 1).mean())\n        / (df['market_return_causal'] ** 2).groupby('market').transform(lambda x: x.rolling(3750, 1).mean())\n    )\n    \n    df['feature_cl_diff1_mean_simple'] = df['feature_cl_diff1'].groupby('timestamp').transform('mean')\n    df['feature_cl_diff1_mean_weight'] = (df['feature_cl_diff1'] * df['weight']).groupby('timestamp').transform('sum') * inv_weight_sum\n    df['feature_cl_diff1_resid'] = df['feature_cl_diff1'] - df['beta_causal'] * df['feature_cl_diff1_mean_weight']\n    \n    df['feature_cl_diff1_rank'] = df.groupby('timestamp')['feature_cl_diff1'].transform('rank')\n    \n    df = df.rename(columns={\n        'beta_causal': 'feature_beta_causal',\n    })\n    \n    return df\n\ndf = pd.read_pickle('/tmp/df.pkl')\ndf = calc_features(df)\nfeatures = get_feature_columns(df)\ndf = df[features + ['target', 'weight']]\ndf.to_pickle('/tmp/df_features.pkl')","metadata":{"execution":{"iopub.status.busy":"2021-11-08T02:47:09.535296Z","iopub.execute_input":"2021-11-08T02:47:09.537149Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.read_pickle('/tmp/df_features.pkl')\ndf = df.dropna()\nfeatures = get_feature_columns(df)\nfor feature in features:\n    print('{} {}'.format(feature, pearsonr(df[feature], df['target'])))\n    \nfor market, df_market in df.groupby('market'):\n    for feature in features:\n        print('{} {} {}'.format(market, feature, pearsonr(df_market[feature], df_market['target'])))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# cv\ndf = pd.read_pickle('/tmp/df_features.pkl')\nfeatures = get_feature_columns(df)\ndf = df.dropna()\n# df = df.loc[df.index.get_level_values(0) < pd.to_datetime('2019-01-01 00:00:00Z')]\ndf = df.loc[df.index.get_level_values(0) < pd.to_datetime('2021-01-01 00:00:00Z')]\n\nmodel = Ridge()\n# model = lgb.LGBMRegressor(n_jobs=-1, random_state=1)\n\nmodel = Pipeline([\n    ('scaler', StandardScaler()),\n    ('model', model)\n])\n\n# model = BaggingRegressor(\n#     model,\n#     n_estimators=1,\n#     random_state=1,\n# )\n\ncv = my_purge_kfold(df.shape[0])\ndf['y_pred'] = cross_val_predict(\n    model,\n    X=df[features],\n    y=df['target'], \n    cv=cv, \n#     n_jobs=-1,\n)\n\nprint(r2_score(df['target'], df['y_pred']))\nprint(pearsonr(df['target'], df['y_pred']))\nprint(df['target'].std())\n\nprint('pearsonr by market')\ndisplay(df.groupby('market').apply(lambda x: pearsonr(x['target'], x['y_pred'])[0]))\n\ndf2 = df.reset_index().set_index('timestamp')\nmarket_count = df2['market'].unique().size\ndf2['target'].rolling(3 * 30 * 24 * 60 * market_count).corr(df2['y_pred']).iloc[::24 * 60 * market_count].plot()\nplt.title('3 month rolling pearsonr')\nplt.show()\n\nif False:\n    for market, df_market in df.groupby('market'):\n        df2 = df_market.reset_index().set_index('timestamp')\n        df2['target'].rolling(3 * 30 * 24 * 60).corr(df2['y_pred']).iloc[::24 * 60].plot()\n        plt.title('3 month rolling pearsonr {}'.format(market))\n        plt.show()\n\n# ボットで使っている特徴量(重いのは削除) + ridgeのスコア\n# r2 0.00012130795512599324\n# pearsonr (0.013457722434615444, 0.0)\n# target std 0.005677241985371386","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# refit\nmodel.fit(df[get_feature_columns(df)], df['target'])\nsave_model(model, 'model.xz')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# submit\n\nrecent_sec = 4000 * (5 * 60)\ndf = pd.read_pickle('/tmp/df.pkl')\ndf_asset = pd.read_pickle('/tmp/df_asset.pkl')\nmodel = joblib.load('model.xz')\n\nenv = gresearch_crypto.make_env() \niter_test = env.iter_test()\nfor (test_df, sample_prediction_df) in iter_test:\n    # 最新データを追加\n    test_df2 = process_data(test_df, df_asset)\n    if 'row_id' in df.columns:\n        df = df.drop(columns=['row_id']) # 念の為row_id重複に対応\n    df = df.append(test_df2)\n    \n    # 最近のデータだけにする(for performance)\n    test_min_timestamp = test_df2.index.get_level_values(0).min()\n    df = df.loc[test_min_timestamp - pd.to_timedelta(recent_sec, unit='s') <= df.index.get_level_values(0)]\n    df = sort_and_remove_duplicates(df)\n    \n    # 特徴量計算\n    df_features = calc_features(df)\n    \n    # 予測\n    df_features = df_features.loc[~df_features['row_id'].isna()]\n    df_features['Target'] = model.predict(df_features[get_feature_columns(df_features)].values)\n    sample_prediction_df = sample_prediction_df.merge(df_features[['row_id', 'Target']], how='left', on='row_id')\n    \n    if False:\n        display(test_df)\n        display(sample_prediction_df)\n    \n    env.predict(sample_prediction_df)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}