{"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":"# ARIMA Model with Exponential Smoothing","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n# https://www.kaggle.com/sohier/detailed-api-introduction\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\n# import os\n# for dirname, _, filenames in os.walk('/kaggle/input'):\n#     for filename in filenames:\n#         print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-12-30T00:29:45.001155Z","iopub.execute_input":"2021-12-30T00:29:45.00143Z","iopub.status.idle":"2021-12-30T00:29:45.005369Z","shell.execute_reply.started":"2021-12-30T00:29:45.00139Z","shell.execute_reply":"2021-12-30T00:29:45.00474Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load dependencies","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport warnings    # `do not disturb`\nwarnings.filterwarnings('ignore')\nimport pickle\nimport time\n\n# plots\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\n# statistics and econometrics\nimport statsmodels.formula.api as smf\nimport statsmodels.tsa.api as smt\nimport statsmodels.tsa.arima.model\nimport statsmodels.api as sm\nimport scipy.stats as scs\nfrom scipy.stats import pearsonr\nfrom scipy.stats import spearmanr\n\n# model performance\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.preprocessing import StandardScaler\nfrom statsmodels.tsa.arima_model import ARIMAResults","metadata":{"execution":{"iopub.status.busy":"2021-12-30T17:37:59.92928Z","iopub.execute_input":"2021-12-30T17:37:59.929663Z","iopub.status.idle":"2021-12-30T17:38:01.616885Z","shell.execute_reply.started":"2021-12-30T17:37:59.929564Z","shell.execute_reply":"2021-12-30T17:38:01.615926Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load data","metadata":{}},{"cell_type":"code","source":"# # Bring in training data and asset information\ndata_folder = \"/kaggle/input/g-research-crypto-forecasting/\"\nasset_details = pd.read_csv(data_folder + 'asset_details.csv')\ntrain = pd.read_csv(data_folder + \"train.csv\")","metadata":{"execution":{"iopub.status.busy":"2021-12-30T17:38:01.619343Z","iopub.execute_input":"2021-12-30T17:38:01.620007Z","iopub.status.idle":"2021-12-30T17:38:56.703285Z","shell.execute_reply.started":"2021-12-30T17:38:01.619958Z","shell.execute_reply":"2021-12-30T17:38:56.702257Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Define custom functions","metadata":{}},{"cell_type":"code","source":"def mean_absolute_percentage_error(y_true, y_pred): \n    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100\n\ndef moving_average(series, n=15):\n    \"\"\"\n    Calculate average of last n observations\n    \"\"\"\n    return np.average(series[-n:])\n\ndef weighted_average(series, weights):\n    \"\"\"\n    Calculate weighted average on the series.\n    Assuming weights are sorted in descending order\n    (larger weights are assigned to more recent observations).\n    \"\"\"\n    result = 0.0\n    for n in range(len(weights)):\n        result += series.iloc[-n-1] * weights[n]\n    return float(result)\n\ndef exponential_smoothing(series, alpha):\n    \"\"\"\n    series - dataset with timestamps\n    alpha - float [0.0, 1.0], smoothing parameter\n    \"\"\"\n    values = series.values\n    result = [values[0]] # first value is same as series\n    for n in range(1, len(series)):\n        result.append(alpha * values[n] + (1 - alpha) * result[n-1])\n    return result\n\ndef plotExponentialSmoothing(series, alphas, plot_intervals = False, scale = 1.96, plot_anomalies=False):\n    \"\"\"\n    Plots exponential smoothing with different alphas\n\n    series - dataset with timestamps\n    alphas - list of floats, smoothing parameters\n        \n    \"\"\"  \n    plt.figure(figsize=(15, 5))\n    for alpha in alphas:\n        plt.plot(exponential_smoothing(series, alpha), c=\"orange\", label=\"Alpha {}\".format(alpha))\n    plt.plot(series.values, \"b\", label = \"Actual\")\n    plt.title(\"Moving Exponential Average\\n Alpha = {}\".format(alphas))\n    plt.legend(loc=\"best\")\n    plt.axis('tight')\n    \n    if plot_intervals:\n        mae = mean_absolute_error(\n            series, exponential_smoothing(series, alpha)\n        )\n        \n        deviation = np.std(series - exponential_smoothing(series, alpha))\n        lower_bound = exponential_smoothing(series, alpha) - (mae + scale * deviation)\n        upper_bound = exponential_smoothing(series, alpha) + (mae + scale * deviation)\n        plt.plot(upper_bound, \"r--\", label=\"Upper Bound / Lower Bound\")\n        plt.plot(lower_bound, \"r--\")\n        \n        # Having the intervals, find abnormal values\n        if plot_anomalies:\n            anomalies = pd.Series(index=series.index, name=series.name)\n            anomalies[series<lower_bound] = series[series<lower_bound]\n            print(str(len(anomalies[series<lower_bound])) + str(\" lower-bound anomalies found\"))\n            anomalies[series>upper_bound] = series[series>upper_bound]\n            print(str(len(anomalies[series<lower_bound])) + str(\" upper-bound anomalies found\"))\n            plt.plot(anomalies, \"black\", markersize=10)\n        \n    plt.legend(loc=\"best\")\n    plt.grid(b=True)\n    \ndef tsplot(y, lags=None, figsize=(12, 7), style='bmh'):\n    \"\"\"\n    Plot time series, its ACF and PACF, calculate Dickey–Fuller test\n\n    y - timeseries\n    lags - how many lags to include in ACF, PACF calculation\n    \"\"\"\n    if not isinstance(y, pd.Series):\n        y = pd.Series(y)\n    \n    \n    with plt.style.context(style):   \n        fig = plt.figure(figsize=figsize)\n        layout = (2, 2)\n        ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)\n        acf_ax = plt.subplot2grid(layout, (1, 0))\n        pacf_ax = plt.subplot2grid(layout, (1, 1))\n        \n        y.plot(ax=ts_ax)\n        p_value = sm.tsa.stattools.adfuller(y)[1]\n        ts_ax.set_title(f'Time Series Analysis Plots\\n Dickey-Fuller: p={p_value:.5f}')\n        smt.graphics.plot_acf(y, lags=lags, ax=acf_ax)\n        smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax)\n        plt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2021-12-30T17:38:56.705329Z","iopub.execute_input":"2021-12-30T17:38:56.705655Z","iopub.status.idle":"2021-12-30T17:38:56.729719Z","shell.execute_reply.started":"2021-12-30T17:38:56.705612Z","shell.execute_reply":"2021-12-30T17:38:56.728637Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training Loop\n- Exponential smoothing (alpha=0.005)\n- ARIMA model with p=2, d=0, q=13","metadata":{}},{"cell_type":"code","source":"n_steps_to_train = 500\nn_steps_to_plot = 100\nn_steps_to_predict = 3\n\n# ARIMA Model params\np=2\nd=0\nq=13\nalpha=0.005\n\n# Initialize dict objects to store models, errors, and forecasts\nmodels = dict()\n\n# Create separate groups for each asset id/name\nfor i in range(len(asset_details)):\n    current_asset_id = asset_details.Asset_ID[i]\n    current_asset_name = asset_details.Asset_Name[i]\n    current_dataset = train[train[\"Asset_ID\"] == current_asset_id].set_index(\"timestamp\")\n    current_dataset.dropna(axis=0, inplace=True)\n    \n    current_dataset = current_dataset.reindex(\n        range(\n            current_dataset.index[0], current_dataset.index[-1] + 60, 60\n        ), method='pad'\n    )\n    \n    # Ensure consistency in distribution of newly created datasets\n    comparison = train[train[\"Asset_ID\"] == current_asset_id].set_index(\"timestamp\")\n    \n    n_steps_to_correlate = 50000\n    first_apple = current_dataset.Close.values[-n_steps_to_correlate:]\n    second_apple = comparison.Close.values[-n_steps_to_correlate:]\n    \n    assert len(first_apple) == len(second_apple), \"The input lengths do not match\"\n    \n    print(f\"Calculating correlations for: {current_asset_name}\")\n    \n    corr, _ = pearsonr(first_apple, second_apple)\n    corr, _ = spearmanr(first_apple, second_apple)\n    \n    print(\"Pearsons: %.3f\" % corr, \"Spearmans: %.3f \\n\" % corr)\n\n    # Calculate exponential smoothing\n    current_dataset_mini = current_dataset.copy()[-n_steps_to_train:]\n    current_dataset_mini[\"Smoothed\"] = exponential_smoothing(current_dataset_mini.Target, alpha=alpha)\n    data = current_dataset_mini[\"Smoothed\"]\n    \n    # Train model\n    model = statsmodels.tsa.arima.model.ARIMA(endog=data.values, exog=None, order=(p,d,q)).fit()\n    aic = model.aic\n    actual = np.array(data.values)\n    \n    # Save model\n    models[current_asset_id] = model","metadata":{"execution":{"iopub.status.busy":"2021-12-30T17:38:56.731628Z","iopub.execute_input":"2021-12-30T17:38:56.732293Z","iopub.status.idle":"2021-12-30T17:40:40.743185Z","shell.execute_reply.started":"2021-12-30T17:38:56.732244Z","shell.execute_reply":"2021-12-30T17:40:40.742204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prediction Loop","metadata":{}},{"cell_type":"code","source":"# for i, (test_df, sample_prediction_df) in enumerate(iter_test):\n#     sample_prediction_df[\"Target\"] = 0.0\n#     env.predict(sample_prediction_df)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gresearch_crypto\nenv = gresearch_crypto.make_env()\niter_test = env.iter_test()\n# from tqdm import tqdm\n\nprediction_step = 0\nfor (test_df, sample_prediction_df) in iter_test:\n    \n    # initialize predictions list\n    y_preds = []\n    \n    # always sort by time, then by Asset_ID\n    test_df = test_df.sort_values(by=['timestamp', 'Asset_ID'])\n    \n    # Map asset_id to row_id\n    asset_id_row_id_map = {asset_id: row_id for asset_id, row_id in test_df[[\"Asset_ID\", \"row_id\"]].values}\n    \n    # Loop through mapping\n    for key in asset_id_row_id_map.keys():\n        model = models[key]\n        y_preds.append(\n            float(\n                model.predict(\n                    start = n_steps_to_train + prediction_step, end = n_steps_to_train + prediction_step\n                )\n            )\n        )\n    \n    # Increment prediction step\n    prediction_step += 1\n    \n    # Update targets and submit\n    sample_prediction_df.Target = y_preds\n    env.predict(sample_prediction_df)","metadata":{"execution":{"iopub.status.busy":"2021-12-30T17:41:56.750535Z","iopub.execute_input":"2021-12-30T17:41:56.750973Z","iopub.status.idle":"2021-12-30T17:41:57.057823Z","shell.execute_reply.started":"2021-12-30T17:41:56.750936Z","shell.execute_reply":"2021-12-30T17:41:57.056918Z"},"trusted":true},"execution_count":null,"outputs":[]}]}