tcn_signals_research.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f4f2af90",
   "metadata": {
    "papermill": {
     "duration": 0.00682,
     "end_time": "2023-11-05T17:59:38.297125",
     "exception": false,
     "start_time": "2023-11-05T17:59:38.290305",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Temporal Convolution Neural Network with Conditioning for Broad Market Signals\n",
    "\n",
    "\n",
    "<a href=\"https://www.kaggle.com/addarm/temporal-convolution-neural-network-for-markets\" target=\"_blank\"><img align=\"left\" alt=\"Kaggle\" title=\"Open in Kaggle\" src=\"https://kaggle.com/static/images/open-in-kaggle.svg\"></a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7699b78c",
   "metadata": {},
   "source": [
    "![image](./images/tcn_banner.PNG)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c10ced8",
   "metadata": {},
   "source": [
    "There is a naive rule-of-thumb strategy amongst the day-trading, that if you buy a levered ETF at open and sell at close, and applying the law of large numbers, you'd be wealthy in no time. This is an interesting problem to solve with deep learning networks. In this article, we will build one called a Temporal Convolusion Network (TCN) for binary classification, which will learn and predict the market direction day by day, if its bullish (and profitable) or bearish.\n",
    "\n",
    "Convolutional neural networks (CNNs) are commonly used for image classification, as the spectral convolutions can extract features and patterns. The same is applicable to timeseries, as proven by the seminal paper from Oord et al. describing Google Deepmind's WaveNet architecture (paper and code in the reference section). We will borrow concepts from WaveNet to teach it to identify deep structures in the market, and try to emmulate a profitable day trader."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a827e09",
   "metadata": {
    "_cell_guid": "b1076dfc-b9ad-4769-8c92-a6c4dae69d19",
    "_uuid": "8f2839f25d086af736a60e9eeb907d3b93b6e0e5",
    "execution": {
     "iopub.execute_input": "2023-11-05T17:59:38.335093Z",
     "iopub.status.busy": "2023-11-05T17:59:38.334737Z",
     "iopub.status.idle": "2023-11-05T18:01:02.444324Z",
     "shell.execute_reply": "2023-11-05T18:01:02.443308Z"
    },
    "papermill": {
     "duration": 84.119909,
     "end_time": "2023-11-05T18:01:02.447010",
     "exception": false,
     "start_time": "2023-11-05T17:59:38.327101",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "from tqdm import tqdm\n",
    "\n",
    "import dotenv\n",
    "%load_ext dotenv\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "IS_KAGGLE = os.getenv('IS_KAGGLE', 'True') == 'True'\n",
    "\n",
    "if IS_KAGGLE:\n",
    "    # Kaggle confgs\n",
    "    print('Running in Kaggle...')\n",
    "    %pip install yfinance\n",
    "    %pip install statsmodels\n",
    "    %pip install seaborn\n",
    "    %pip install itertools\n",
    "    %pip install scikit-learn\n",
    "    %pip install tensorflow\n",
    "    %pip install shutils\n",
    "    %pip install seaborn\n",
    "    %pip install tensorboard\n",
    "    %pip install tqdm\n",
    "    %pip install matplotlib\n",
    "    %pip install python-dotenv\n",
    "    for dirname, _, filenames in os.walk('/kaggle/input'):\n",
    "        for filename in filenames:\n",
    "            print(os.path.join(dirname, filename))\n",
    "else:\n",
    "    print('Running Local...')\n",
    "\n",
    "import yfinance as yf\n",
    "\n",
    "\n",
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f226f2c",
   "metadata": {},
   "source": [
    "# Financial Data\n",
    "\n",
    "Using `yfinance`, we will pull 20 years of market data, specifically: \n",
    "- **SPY** ETF - which will be our target.\n",
    "- CBoE's **VIX** as a proxy for market senitment and volatility.\n",
    "- The **Russell500** small caps index, as a proxy for speculation.\n",
    "- The **Gold index** as proxy for investor uncertainty.\n",
    "- The **10 year US treasury note** as proxy for inflation and rates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38727fc7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-05T18:01:02.470424Z",
     "iopub.status.busy": "2023-11-05T18:01:02.469635Z",
     "iopub.status.idle": "2023-11-05T18:01:03.256388Z",
     "shell.execute_reply": "2023-11-05T18:01:03.255137Z"
    },
    "papermill": {
     "duration": 0.801372,
     "end_time": "2023-11-05T18:01:03.258937",
     "exception": false,
     "start_time": "2023-11-05T18:01:02.457565",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "from scipy.stats import skew, kurtosis\n",
    "from pandas.tseries.offsets import BDay\n",
    "\n",
    "START_DATE = \"1999-01-01\"\n",
    "END_DATE = pd.Timestamp(datetime.now() - BDay(1)).strftime('%Y-%m-%d')\n",
    "DATA_DIR = \"data\"\n",
    "INDEX = \"Date\"\n",
    "TARGET_ETF = \"SPY\"  # S&P 500\n",
    "RATES_INDEX = \"^TNX\"  # 10 Year Treasury Note Yield\n",
    "VOLATILITY_INDEX = \"^VIX\"  # CBOE Volatility Index\n",
    "SMALLCAP_INDEX = \"^RUT\" # Russell 2000 Index\n",
    "GOLD_INDEX = \"GC=F\" # Gold futures\n",
    "tickers_symbols = [\n",
    "    TARGET_ETF,\n",
    "    VOLATILITY_INDEX,\n",
    "    RATES_INDEX,\n",
    "    SMALLCAP_INDEX,\n",
    "    GOLD_INDEX,\n",
    "]\n",
    "INTERVAL = \"1d\"\n",
    "\n",
    "def get_tickerdata(tickers_symbols, start=START_DATE, end=END_DATE, interval=INTERVAL, datadir=DATA_DIR):\n",
    "    tickers = {}\n",
    "    earliest_end= datetime.strptime(end,'%Y-%m-%d')\n",
    "    latest_start = datetime.strptime(start,'%Y-%m-%d')\n",
    "    os.makedirs(DATA_DIR, exist_ok=True)\n",
    "    for symbol in tickers_symbols:\n",
    "        cached_file_path = f\"{datadir}/{symbol}-{start}-{end}-{interval}.csv\"\n",
    "\n",
    "        try:\n",
    "            if os.path.exists(cached_file_path):\n",
    "                df = pd.read_csv(cached_file_path, index_col=INDEX)\n",
    "                df.index = pd.to_datetime(df.index)\n",
    "                assert len(df) > 0\n",
    "            else:\n",
    "                df = yf.download(\n",
    "                    symbol,\n",
    "                    start=START_DATE,\n",
    "                    end=END_DATE,\n",
    "                    progress=False,\n",
    "                    interval=INTERVAL,\n",
    "                )\n",
    "                assert len(df) > 0\n",
    "                df.to_csv(cached_file_path)\n",
    "            min_date = df.index.min()\n",
    "            max_date = df.index.max()\n",
    "            nan_count = df[\"Close\"].isnull().sum()\n",
    "            skewness = round(skew(df[\"Close\"].dropna()), 2)\n",
    "            kurt = round(kurtosis(df[\"Close\"].dropna()), 2)\n",
    "            outliers_count = (df[\"Close\"] > df[\"Close\"].mean() + (3 * df[\"Close\"].std())).sum()\n",
    "            print(\n",
    "                f\"{symbol} => min_date: {min_date}, max_date: {max_date}, kurt:{kurt}, skewness:{skewness}, outliers_count:{outliers_count},  nan_count: {nan_count}\"\n",
    "            )\n",
    "            tickers[symbol] = df\n",
    "\n",
    "            if min_date > latest_start:\n",
    "                latest_start = min_date\n",
    "            if max_date < earliest_end:\n",
    "                earliest_end = max_date\n",
    "        except Exception as e:\n",
    "            print(f\"Error with {symbol}: {e}\")\n",
    "\n",
    "    return tickers, latest_start, earliest_end\n",
    "\n",
    "tickers, latest_start, earliest_end = get_tickerdata(tickers_symbols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "946b0d99",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Date ranges {latest_start} - {earliest_end}\")\n",
    "for ticker in tickers_symbols:\n",
    "    df = tickers.get(ticker)\n",
    "    df = df[(df.index >= latest_start) & (df.index <= earliest_end)]\n",
    "\n",
    "    assert len(df) > 0 and not df.isna().any().any()\n",
    "\n",
    "    tickers[ticker] = df\n",
    "\n",
    "tickers.get(TARGET_ETF).head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "122af07e",
   "metadata": {},
   "outputs": [],
   "source": [
    "percentage_returns = {}\n",
    "\n",
    "for ticker in tickers_symbols:\n",
    "    df = tickers.get(ticker)\n",
    "    percentage_returns[ticker] = ((1 + df[\"Close\"].pct_change()).cumprod() - 1) * 100\n",
    "plt.figure(figsize=(16, 6))\n",
    "for ticker, pr in percentage_returns.items():\n",
    "    plt.plot(pr, label=ticker, alpha=0.75)\n",
    "\n",
    "plt.legend(loc=\"upper left\")\n",
    "plt.grid(True)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83e06e38",
   "metadata": {},
   "source": [
    "# Naive Baseline Strategy\n",
    "\n",
    "Buy at open, sell at close - easy right? Or that is what the financial gurus on youtube want to persuade us to believe and buy their courses.\n",
    "\n",
    "Let's do a quich and dirty test for this null hypothesis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34c42ed0",
   "metadata": {},
   "outputs": [],
   "source": [
    "TARGET_FACTOR = 55/100/100  # assume 50BP spread on average\n",
    "print(f\"target factor used: {TARGET_FACTOR} for window: {WINDOW_SIZE} predicting: {PREDICTION_HORIZON} step\")\n",
    "\n",
    "target_etf = tickers.get(TARGET_ETF)\n",
    "\n",
    "next_close = (target_etf[\"Close\"].shift(-1) * (1- TARGET_FACTOR))\n",
    "days_rising = target_etf[next_close > target_etf[\"Close\"]]\n",
    "days_not_rising = target_etf[next_close <= target_etf[\"Close\"]]\n",
    "days_rising_count = days_rising.groupby(days_rising.index.year).size()\n",
    "days_not_rising_count = days_not_rising.groupby(days_not_rising.index.year).size()\n",
    "\n",
    "total_days = target_etf.groupby(target_etf.index.year).size()\n",
    "percentage_rising = (days_rising_count / total_days) * 100\n",
    "\n",
    "yearly_counts = pd.DataFrame(\n",
    "    {\"Rising\": days_rising_count, \"Not Rising\": days_not_rising_count}\n",
    ")\n",
    "yearly_counts.plot(kind=\"bar\", color=[\"green\", \"red\"], figsize=(16, 4))\n",
    "plt.title(\"Days Rising vs Not Rising per Year\")\n",
    "plt.xlabel(\"Year\")\n",
    "plt.ylabel(\"Count\")\n",
    "plt.tick_params(axis=\"x\", rotation=45)\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(f\"Baseline Accuracy: {percentage_rising.mean():0.2f}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfdd0638",
   "metadata": {},
   "source": [
    "Looking at these numbers, we have 29% of making a profit with this strategy, such a low ROI is true because we are factoring in a 50 basis points spread, and various fees and commisions the broker would charge us. Therefore the hypothesis is rejected, and the alternative is true, you'll be unprofitabled - unless you leverage some state-of-the-art neural network (NN)\n",
    "\n",
    "We know that this NN precision score should outperform the naive baseline's 29% to be succesful.\n",
    "\n",
    "# Feature Selection and Engineering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e5751b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import FunctionTransformer, MinMaxScaler, normalize\n",
    "\n",
    "MONTH_SINE = \"month_sin\"\n",
    "MONTH_COS = \"month_cos\"\n",
    "MONTH_RBF = \"month_rbf\"\n",
    "DAY_SINE = \"day_sin\"\n",
    "DAY_COS = \"day_cos\"\n",
    "DAY_RBF = \"day_rbf\"\n",
    "Q_SINE = \"quart_sin\"\n",
    "Q_COS = \"quart_cos\"\n",
    "Q_RBF = \"quart_rbf\"\n",
    "BIZ_SINE = \"biz_sin\"\n",
    "BIZ_COS = \"biz_cos\"\n",
    "BIZ_RBF = \"biz_rbf\"\n",
    "TARGET_LABEL = \"Close_target\"\n",
    "TARGET_TS = \"Close\"  # The TS which we will condition.\n",
    "EXOG_TS = tickers_symbols.copy() # The TS to use for conditioning\n",
    "EXOG_TS.remove(TARGET_ETF)\n",
    "PRICE_FEATURES = [\"Open\", \"High\", \"Low\", \"Close\"]\n",
    "TIME_FEATURES = [\n",
    "    DAY_RBF,\n",
    "    MONTH_RBF,\n",
    "    Q_RBF,\n",
    "    BIZ_RBF,\n",
    "]\n",
    "TS_FEATURES = [\"Volume\"]\n",
    "TARGET_TS_FEATURES_NOTIME = EXOG_TS + TS_FEATURES + PRICE_FEATURES\n",
    "TARGET_TS_FEATURES = TARGET_TS_FEATURES_NOTIME + TIME_FEATURES\n",
    "WINDOW_SIZE = 252 // 4  # 1 years trading, sampled across days\n",
    "PREDICTION_HORIZON = 1  # next 1 trading day"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a378eb58",
   "metadata": {},
   "source": [
    "Before even consuming the financial timeserieses, we need to do some wrangling and feature engineering. Neural Networks (NN) see numbers differently than us, they are interested in magnitudes, cycles and distances, this means we have to encode our data in a way that makes a NN learn better and be more efficient, and accurate.\n",
    "\n",
    "Throughout the experiments we used log returns (capturing percentage moves): $r_t = \\log\\left(\\frac{P_t}{P_{t-1}}\\right)$, first order of integration (differencing the values of today with yesterday): $Y_t = \\sum_{i=0}^{t} \\Delta X_i$, and finally only normalizing the raw data: $X_{\\text{norm}} = \\frac{X - \\mu}{\\sigma}$. \n",
    "\n",
    "We chose the first order differencing - normalized - we'll revisit this in the architecture below.\n",
    "\n",
    "For temporal data, we took a different approach. Since time is cyclic (even years, you can think of these as business cycles) it can be encoded into sin/cosine waves, for example days can be encoded in the following equation:\n",
    "\n",
    "$$\n",
    "\\text{sin}_{d} = \\sin\\left(\\frac{2\\pi d}{365}\\right)\n",
    "\\newline\n",
    "\\text{cos}_{d} = \\cos\\left(\\frac{2\\pi d}{365}\\right)\n",
    "$$\n",
    "\n",
    "Though a more accurate representation of time is available through Radial Basis Functions (RBF) $K(d, d') = \\exp\\left(-\\frac{(d - d')^2}{2\\sigma^2}\\right)$. This is more accurate as distances are measured from a focal point, therefore Sunday is larger than Monday although in the days' sequence they are side-by-side. In a cyclic function, the distance of Sunday is not well represented, this while the RBF one gives the right distance representation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1f1f024",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def create_time_features(data_df):\n",
    "    \"\"\"\n",
    "    Encodes time cyclic features for a dataset with monthly sampling.\n",
    "    Including cyclic encoding for day and year.\n",
    "    :param data_df: The timeseries with a date in the format YYYY-MM-DD as index.\n",
    "    :return: data_df with added wave features for month, day, and year.\n",
    "    \"\"\"\n",
    "    if not isinstance(data_df.index, pd.DatetimeIndex):\n",
    "        raise ValueError(\"The DataFrame index must be a DateTimeIndex.\")\n",
    "\n",
    "    def _sin_transformer(period):\n",
    "        return FunctionTransformer(lambda x: np.sin(x / period * 2 * np.pi))\n",
    "\n",
    "    def _cos_transformer(period):\n",
    "        return FunctionTransformer(lambda x: np.cos(x / period * 2 * np.pi))\n",
    "\n",
    "    def rbf_transform(x, period, input_range):\n",
    "        x_normalized = (x - input_range[0]) / (input_range[1] - input_range[0]) * period\n",
    "        return np.exp(-0.5 * ((x_normalized - period / 2) / 1.0) ** 2)\n",
    "\n",
    "    data_df[DAY_RBF] = rbf_transform(data_df.index.day, 31, (1, 31))\n",
    "    data_df[MONTH_RBF] = rbf_transform(data_df.index.month, 12, (1, 12))\n",
    "    data_df[Q_RBF] = rbf_transform((data_df.index.month - 1) // 3 + 1, 4, (1, 4))\n",
    "    min_year, max_year = data_df.index.year.min(), data_df.index.year.max()\n",
    "    data_df[BIZ_RBF] = rbf_transform(\n",
    "        data_df.index.year, max_year - min_year, (min_year, max_year)\n",
    "    )\n",
    "\n",
    "    return data_df\n",
    "\n",
    "def create_features_df(tickers, data_df, target_label=TARGET_ETF):\n",
    "    \"\"\"\n",
    "    Create all exogenous features that lead to our target etf.\n",
    "        - if the trading day close is higher than the open.\n",
    "        - price log returns or 1D integrations.\n",
    "    :param tickers: All the timeseries with a date in the format YYYY-MM-DD as index.\n",
    "    :param data_df: Any pre-engineered features.\n",
    "    :return: data_df With TARGET_LABEL timeseries at column 0, and the rest are for conditioning.\n",
    "    \"\"\"\n",
    "\n",
    "    def _get_first_difference(data_df):\n",
    "        return data_df.pct_change().fillna(0)\n",
    "\n",
    "    def _get_log_returns(data_df):\n",
    "        return np.log(data_df / data_df.shift(1)).fillna(0)\n",
    "\n",
    "    IDX_COL = \"Open\"\n",
    "    price_transform = FunctionTransformer(_get_first_difference)\n",
    "\n",
    "    data_df[PRICE_FEATURES] = price_transform.fit_transform(data_df[PRICE_FEATURES])\n",
    "    rates_df = tickers.get(RATES_INDEX)\n",
    "\n",
    "    for index in EXOG_TS:\n",
    "        if index == target_label:\n",
    "            continue\n",
    "        index_df = tickers.get(index)\n",
    "        transformed_data = price_transform.fit_transform(index_df[[IDX_COL]])\n",
    "        data_df[index] = transformed_data\n",
    "\n",
    "    data_df = data_df.fillna(0)\n",
    "\n",
    "    return data_df\n",
    "\n",
    "IS_CLASSIFICATION = True\n",
    "\n",
    "def prepare_data(tickers, target_label=TARGET_ETF, classify=IS_CLASSIFICATION, to_normalize=True):\n",
    "    \"\"\"\n",
    "    Utility function to prepare the data.\n",
    "    :data_df dataframe: dataframe with `window_size` months of data to predict the `window_size`+`horizon`.\n",
    "    :param window_size: int, length of the input sequence\n",
    "    :param horizon: int, forecasting horizon, defaults to 1\n",
    "    :return: Array in the shape of (n_samples, n_steps, n_features)\n",
    "    \"\"\"\n",
    "    y_scaler = None, None\n",
    "    data_df = tickers.get(target_label).copy()\n",
    "    close_tomorrow_df = (data_df['Close'].shift(-1) * (1- TARGET_FACTOR))\n",
    "    if classify:\n",
    "        today_df = data_df['Close'] # data_df['Open'].shift(-1)\n",
    "        # Calculate the target label: 1 if next day's close is higher than its open.\n",
    "        data_df[TARGET_LABEL] = (close_tomorrow_df > today_df).astype(int)\n",
    "    else:\n",
    "        # A target variable with a large spread of values, in turn, may result in large error gradient values\n",
    "        # causing weight values to change dramatically. We predict tomorrows close.\n",
    "        y_scaler = MinMaxScaler(feature_range=(-1, 1))\n",
    "        y_scaled = y_scaler.fit_transform(close_tomorrow_df)\n",
    "        data_df[TARGET_LABEL] = y_scaled.flatten()\n",
    "\n",
    "    data_df = create_features_df(tickers, data_df, target_label=target_label)\n",
    "    data_df_normalized = None\n",
    "    if to_normalize:\n",
    "        data_df_normalized = normalize(data_df[TARGET_TS_FEATURES_NOTIME], norm=\"l2\")\n",
    "    else:\n",
    "        data_df_normalized = data_df[TARGET_TS_FEATURES_NOTIME]\n",
    "    data_df_normalized = pd.DataFrame(\n",
    "        data_df_normalized, columns=TARGET_TS_FEATURES_NOTIME\n",
    "    )\n",
    "    data_df_normalized = pd.concat(\n",
    "        [\n",
    "            data_df[TARGET_LABEL].reset_index(drop=True),\n",
    "            data_df_normalized.reset_index(drop=True),\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "    if any(feature in TIME_FEATURES for feature in TARGET_TS_FEATURES):\n",
    "        # Don't normalize these.\n",
    "        data_df = create_time_features(data_df)\n",
    "        data_df_normalized = pd.concat(\n",
    "            [\n",
    "                data_df[TIME_FEATURES].reset_index(drop=True),\n",
    "                data_df_normalized.reset_index(drop=True),\n",
    "            ],\n",
    "            axis=1,\n",
    "        )\n",
    "    # we might drop the first and last row of data (shifts and difference => NAs)\n",
    "    return data_df_normalized.dropna(axis=0), y_scaler\n",
    "\n",
    "data_df, y_scaler = prepare_data(tickers, to_normalize=True)\n",
    "data_df.head(15)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1c200b8",
   "metadata": {},
   "source": [
    "## Data Analyis & Feature Selection\n",
    "\n",
    "let's start with a simple check, if we are to predict the number of days the target instrument (SPY in our case) closes higher than the start (we know already its 29% of the time) and factor in spreads and fees, is our data set balanced?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bce56bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "label_counts = data_df[TARGET_LABEL].value_counts()\n",
    "\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.bar(label_counts.index.astype(str), label_counts.values, color=['blue', 'red'])\n",
    "plt.xlabel('Target Label')\n",
    "plt.ylabel('Count')\n",
    "plt.title('Distribution of Target Labels (Close>Open)')\n",
    "plt.xticks([0, 1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abad8572",
   "metadata": {},
   "source": [
    "No its not - this is can be a problem for the NN as it can get biased on the 0 label (not higher than the start). In this case we have to use a specific loss function called Focal Loss as an alternative to the standard Cross-Enthropy, more on this in the model training section.\n",
    "\n",
    "### Principla Component Analysis (PCA)\n",
    "\n",
    "We carry out a PCA to assert which features have most info to our target, and attempt to identify anomalies in our timeseries. Let's start with a timeseries definition: An interdependant multiplicative timeseries (as most financial timeseries are) is comprised of 4 components:\n",
    "\n",
    "$Y(t)=T(t)×S(t)×C(t)×\\epsilon(t)$\n",
    "\n",
    "Where $T$ is the trend, $S$ is the seasonal trend, $C$ is the cyclical trend and $\\epsilon$ is noise in any point in time $t$. We know that all our data is:\n",
    "1. Autoregressive and probably multicollinear\n",
    "2. Not stationary, the mean moves.\n",
    "3. Has trend, seasonality and all manner of noisome structures.\n",
    "\n",
    "With this in mind, PCA will decompose these into components of meaningful variance. The cummulative variance shows which features are of importance in our data and which will contribute most information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0262feac",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "MAX_VARIANCE = 0.95\n",
    "\n",
    "\n",
    "data_df_pca = data_df.drop(columns=[TARGET_LABEL])\n",
    "pca = PCA()\n",
    "xdata = pca.fit_transform(data_df_pca)\n",
    "\n",
    "cum_var_exp = np.cumsum(pca.explained_variance_ratio_)\n",
    "num_components = np.argmax(cum_var_exp >= MAX_VARIANCE) + 1\n",
    "print(f\"Max components for {MAX_VARIANCE*100}% variance: {num_components} out of {data_df.shape[1]}\")\n",
    "\n",
    "eigenvectors = pca.components_\n",
    "loadings_df = pd.DataFrame(eigenvectors, columns=data_df_pca.columns).T\n",
    "summed_loadings = np.sum(np.abs(eigenvectors), axis=0)\n",
    "summed_loadings_df = pd.DataFrame(summed_loadings, index=data_df_pca.columns, columns=[\"Sum\"]).sort_values(by=\"Sum\", ascending=False)\n",
    "\n",
    "fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(18, 10), gridspec_kw={'height_ratios': [1, 2, 4]})\n",
    "\n",
    "summed_loadings_df.plot(kind=\"bar\", legend=False, ax=axes[0])\n",
    "axes[0].set_title(\"Summed Loadings Across All Principal Components\")\n",
    "axes[0].set_ylabel(\"Summed Loadings\")\n",
    "axes[0].set_xlabel(\"Features\")\n",
    "axes[0].tick_params(axis='x', labelrotation=45)\n",
    "\n",
    "\n",
    "loadings_df.plot(kind=\"bar\", legend=False, ax=axes[1])\n",
    "axes[1].set_title(f\"Loadings for Principal Components\")\n",
    "axes[1].set_ylabel(\"Loadings\")\n",
    "axes[1].set_xlabel(\"Features\")\n",
    "axes[1].tick_params(axis='x', labelrotation=45)\n",
    "axes[1].legend(loc='upper left', bbox_to_anchor=(1, 1))\n",
    "\n",
    "pca_reduced = PCA(n_components=num_components)\n",
    "xdata_reduced = pca_reduced.fit_transform(data_df_pca)\n",
    "xdata_projected = pca_reduced.inverse_transform(xdata_reduced)\n",
    "error = np.sum((data_df_pca - xdata_projected) ** 2, axis=1)\n",
    "anomaly_threshold = np.percentile(error, MAX_VARIANCE*100)\n",
    "anomalies = error > anomaly_threshold\n",
    "\n",
    "data_df_pca = data_df.copy()\n",
    "anomaly_dates = data_df.index[np.where(anomalies)]\n",
    "for i, column in enumerate(data_df_pca.columns):\n",
    "    if column == TARGET_LABEL or column in TIME_FEATURES or column == \"Volume\":\n",
    "        continue\n",
    "    axes[2].scatter(anomaly_dates, data_df_pca.loc[anomaly_dates, column], label='Anomaly' if i == 0 else \"__nolegend__\", color='k')\n",
    "    axes[2].plot(data_df_pca.index, data_df_pca[column], label=f'{column}', alpha=0.8)\n",
    "axes[2].set_title('Time Series with Anomalies Highlighted')\n",
    "axes[2].legend(loc='upper left', bbox_to_anchor=(1, 1))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9b3d01c",
   "metadata": {},
   "source": [
    "By using PCA we can deconstruct and reconstruct with reduced components to indentify variance, and therefore anomalies. These are plotted in the chart above, with the black circles signalling outliers above the 95th percentile\n",
    "\n",
    "The chart is interesting, the market movements during the 2008 global recession were relatively uniform across all sectors (broad market downturn), returns might move together in a way that **PCA considers normal**. In contrast, the periods of 2000-2005 (DotCom bubble) and post-2019 (Covid) have exhibited more varied behaviors across these indices, leading to higher reconstruction errors for those periods.\n",
    "\n",
    "We also use the Information Coefficient (IC), calculated here using **Spearman's rank correlation** to understand the features' importance. It measures the strength and direction of a monotonic relationship between two variables. IC is useful because it can capture nonlinear relationships between variables. A high absolute value of the IC indicates a strong relationship, which could be either positive or negative. Spearman's rank is represented by the following formula:\n",
    "\n",
    "$$\n",
    "\\rho = 1 - \\frac{6 \\sum d_i^2}{n(n^2 - 1)}\n",
    "$$\n",
    "\n",
    "- Where the sum of d represents the sum of the squared differences in ranks between each pair of observations.\n",
    "- n is the number of observations.\n",
    "- The constant 6 is part of the normalization factor that scales the sum of squared rank differences.\n",
    "\n",
    "The **Mutual Information (MI)** also quantifies the amount of information obtained about one random variable through observing the other random variable. It's measuring how much information each feature in our dataset provides about the target variable. Unlike correlation, MI can capture any kind of relationship between variables, not just linear or monotonic. MI between X and Y can be represented by the formula:\n",
    "\n",
    "$$\n",
    "I(X; Y) = \\sum_{y \\in Y} \\sum_{x \\in X} p(x, y) \\log \\left( \\frac{p(x, y)}{p(x)p(y)} \\right)\n",
    "$$\n",
    "\n",
    "- p(x,y) is the joint probability distribution function of X and Y, indicating the probability that X takes on value x and Y takes on value y simultaneously.\n",
    "- p(x) and p(y) are the marginal probability distribution functions of X and Y, respectively, indicating the probabilities of X taking on value x and Y taking on value y independently.\n",
    "\n",
    "Many good things come from feature selection: less chance of over-fitting, faster training times, more accurate models - though sometimes at the cost of interpretability. In our case all features are deemed important within an economic context, though we ought to finetune a bit in the next section. The graph below represents the IC and MI of our current features:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26f5cc4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import spearmanr\n",
    "from sklearn.feature_selection import mutual_info_regression\n",
    "\n",
    "ic = {}\n",
    "for column in data_df.columns:\n",
    "    if column != TARGET_LABEL:\n",
    "        corr, p_val = spearmanr(data_df[TARGET_LABEL], data_df[column])\n",
    "        ic[column] = [corr, p_val]\n",
    "\n",
    "ic_df = pd.DataFrame(ic, index=[\"IC\", \"p-value\"]).T\n",
    "\n",
    "mi = mutual_info_regression(X=data_df.drop(columns=[TARGET_LABEL]), y=data_df[TARGET_LABEL])\n",
    "mi_series = pd.Series(mi, index=data_df.drop(columns=[TARGET_LABEL]).columns)\n",
    "metrics = pd.concat(\n",
    "    [\n",
    "        mi_series.to_frame(\"Mutual Information\"),\n",
    "        ic_df[\"IC\"].to_frame(\"Information Coefficient\"),\n",
    "    ],\n",
    "    axis=1,\n",
    ")\n",
    "metrics = metrics.sort_values(by=\"Mutual Information\", ascending=False)\n",
    "ax = metrics.plot.bar(figsize=(18, 4), rot=45)\n",
    "ax.set_xlabel(\"Features\")\n",
    "ax.set_ylabel(\"Scores\")\n",
    "ax.set_title(\"Feature Evaluation: MI & IC\")\n",
    "sns.despine()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fee3767",
   "metadata": {},
   "source": [
    "## Feature Selection\n",
    "\n",
    "Let's do a union of the important features from both PCA and IC/MI methods, with the assumption that only these are what the model needs.\n",
    "\n",
    "There are other checks we can do, such as correlation and cointegration, plus checking their T-scores. Sadly in earlier experiments, this feature set was too simple for the model, causing it to have high bias. We will train it with the full load."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2fcad60",
   "metadata": {},
   "outputs": [],
   "source": [
    "MAX_FEATURES_COUNT = 10\n",
    "\n",
    "features_pca = summed_loadings_df.head(MAX_FEATURES_COUNT).index.tolist()\n",
    "features_miic = (metrics.head(MAX_FEATURES_COUNT).index.tolist())\n",
    "print(F\"Top {MAX_FEATURES_COUNT} PCA Loadings: {features_pca}\")\n",
    "print(F\"Top {MAX_FEATURES_COUNT} MI/IC: {features_miic}\")\n",
    "SELECTED_FEATURES = list(set(features_pca) & set(features_miic))\n",
    "\n",
    "print(f\"Selected {len(SELECTED_FEATURES)} features: {SELECTED_FEATURES}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dc9a383",
   "metadata": {},
   "outputs": [],
   "source": [
    "SELECTED_EXOG =  EXOG_TS # [feature for feature in SELECTED_FEATURES if feature in EXOG_TS] #\n",
    "SELECTED_FEATURES = data_df.drop(columns=EXOG_TS+[\"Open\", \"Volume\"], axis=1).columns.to_list() # [feature for feature in SELECTED_FEATURES if feature not in EXOG_TS +[\"Open\", \"Volume\"])] #\n",
    "\n",
    "SELECTED_EXOG, SELECTED_FEATURES"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af3cb19b",
   "metadata": {},
   "source": [
    "Finally, to make sure we have only features that contribute new information, we test for multicolinearity with **Variance Inflation Factor** (VIF): $VIF_i = \\frac{1}{1 - R_i^2}$ where Ri coefficient of determination of a regression. For feature not to be colinear, they have to score between 1 to 5 (ignoring the constant term which we placed there for an intercept to help in the variance calculation):\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "489e9928",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.stats.outliers_influence import variance_inflation_factor\n",
    "from statsmodels.tools.tools import add_constant\n",
    "\n",
    "data_with_constant = add_constant(data_df[SELECTED_FEATURES + SELECTED_EXOG])\n",
    "vif_data = pd.DataFrame()\n",
    "vif_data[\"feature\"] = data_with_constant.columns\n",
    "vif_data[\"VIF\"] = [\n",
    "    variance_inflation_factor(data_with_constant.values, i)\n",
    "    for i in range(data_with_constant.shape[1])\n",
    "]\n",
    "\n",
    "print(vif_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05e574d6",
   "metadata": {},
   "source": [
    "# Temporal Convolution Neural Network (TCN)\n",
    "\n",
    "CNNs, RNNs and other stateful deeplearning models are used for timeseries (including transformers, as language is also a timeseries). The TCN is made specifically for timeseries problems, following Van den Oord et al. architecture, we will build one that is conditioned for the capital markets:\n",
    "1. 2 Input 3D Tensors of shape `(batch_size, window_size, n_features)` for the target timeseries, and `(batch_size, window_size, n_timeseries)` for the conditioning timeseries sets.\n",
    "2. Output 2D tensor of shape `(batch_size, horizon)`\n",
    "3. 2 conditioning layers\n",
    "   1. x1 Dilated 1D convolutions for the exogenous timeseries to condition the main with.\n",
    "   2. x1 Dilated 1D convolutions for the main timeseries.\n",
    "   3. x2 residual layers to add back historic data to the x2 convolution outputs.\n",
    "   4. x2 Add layers to reinstate residuals.\n",
    "   5. x1 contenate layer to combine and output the conditioned main timeseries and residuals.\n",
    "4. 6 hidenn layers (from the paper's code) comprised of:\n",
    "   1. x2 1D [convulations](https://www.tensorflow.org/tutorials/structured_data/time_series#cnn) with relu activation and a spatial dropouts.\n",
    "   2. 1D Dilated Convolution to capture residuals with linear activation.\n",
    "   3. An addition layer to add back the [residuals](https://www.tensorflow.org/tutorials/structured_data/time_series#advanced_residual_connections) into the next layers input\n",
    "5. a single dense layer to output the next timestep according to the given horizon.\n",
    "6. Training is done with ADAM learning optimizer, early stop function, and BinaryFocalCrossentropy to minimize.\n",
    "\n",
    "The architecture is visualized in the graph below, taken from tenserboard's graph plot (the arrows that point to the center, are going out of the block to the next one):\n",
    "\n",
    "\n",
    "![NNGraph](./images/fullgraph_tcn.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "421953b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "from tensorflow.keras.layers import (\n",
    "    SpatialDropout1D,\n",
    "    Dense,\n",
    "    Conv1D,\n",
    "    Layer,\n",
    "    Add,\n",
    "    Input,\n",
    "    Concatenate,\n",
    "    Flatten,\n",
    ")\n",
    "from tensorflow.keras import Model\n",
    "\n",
    "assert tf.__version__ >= \"2.12\"\n",
    "\n",
    "NONLINEAR_ACTIVATION = \"tanh\" # LeakyReLU(alpha=0.01)\n",
    "\n",
    "\n",
    "class TCNBlock(Layer):\n",
    "    \"\"\"\n",
    "    TCN Residual Block that uses zero-padding to maintain `steps` value of the ouput equal to the one in the input.\n",
    "    Residual Block is obtained by stacking togeather (2x) the following:\n",
    "        - 1D Dilated Convolution\n",
    "        - ReLu\n",
    "        - Spatial Dropout\n",
    "    And adding the input after trasnforming it with a 1x1 Conv\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(\n",
    "        self,\n",
    "        filters=1,\n",
    "        kernel_size=2,\n",
    "        dilation_rate=1,\n",
    "        kernel_initializer=\"glorot_normal\",\n",
    "        bias_initializer=\"glorot_normal\",\n",
    "        kernel_regularizer=None,\n",
    "        bias_regularizer=None,\n",
    "        use_bias=False,\n",
    "        dropout_rate=0.0,\n",
    "        layer_id=None,\n",
    "        **kwargs,\n",
    "    ):\n",
    "        \"\"\" \"\n",
    "        Arguments\n",
    "            filters: Integer, the dimensionality of the output space\n",
    "                (i.e. the number of output filters in the convolution).\n",
    "            kernel_size: An integer or tuple/list of a single integer,\n",
    "                specifying the length of the 1D convolution window.\n",
    "            dilation_rate: an integer or tuple/list of a single integer, specifying\n",
    "                the dilation rate to use for dilated convolution.\n",
    "                Usually dilation rate increases exponentially with the depth of the network.\n",
    "            activation: Activation function to use\n",
    "                If you don't specify anything, no activation is applied\n",
    "                (ie. \"linear\" activation: `a(x) = x`).\n",
    "            use_bias: Boolean, whether the layer uses a bias vector.\n",
    "            kernel_initializer: Initializer for the `kernel` weights matrix\n",
    "            bias_initializer: Initializer for the bias vector\n",
    "            kernel_regularizer: Regularizer function applied to the `kernel` weights matrix\n",
    "            bias_regularizer: Regularizer function applied to the bias vector\n",
    "                (see [regularizer](../regularizers.md)).\n",
    "        # Input shape\n",
    "            3D tensor with shape: `(batch, steps, n_features)`\n",
    "        # Output shape\n",
    "            3D tensor with shape: `(batch, steps, filters)`\n",
    "        \"\"\"\n",
    "        super(TCNBlock, self).__init__(**kwargs)\n",
    "        assert dilation_rate is not None and dilation_rate > 0 and filters > 0 and kernel_size > 0\n",
    "\n",
    "        self.filters = filters\n",
    "        self.kernel_size = kernel_size\n",
    "        self.dilation_rate = dilation_rate\n",
    "        self.kernel_initializer = kernel_initializer\n",
    "        self.bias_initializer = bias_initializer\n",
    "        self.kernel_regularizer = kernel_regularizer\n",
    "        self.bias_regularizer = bias_regularizer\n",
    "        self.use_bias = use_bias\n",
    "        self.dropout_rate = dropout_rate\n",
    "        self.layer_id = str(layer_id)\n",
    "\n",
    "    def get_config(self):\n",
    "        config = super(TCNBlock, self).get_config()\n",
    "        config.update({\n",
    "            'filters': self.filters,\n",
    "            'kernel_size': self.kernel_size,\n",
    "            'dilation_rate': self.dilation_rate,\n",
    "            'kernel_initializer': self.kernel_initializer,\n",
    "            'bias_initializer': self.bias_initializer,\n",
    "            'kernel_regularizer': self.kernel_regularizer,\n",
    "            'bias_regularizer': self.bias_regularizer,\n",
    "            'use_bias': self.use_bias,\n",
    "            'dropout_rate': self.dropout_rate,\n",
    "        })\n",
    "        return config\n",
    "\n",
    "    def build(self, inputs):\n",
    "        # Capture feature set from the input\n",
    "        self.conv1 = Conv1D(\n",
    "            filters=self.filters,\n",
    "            kernel_size=self.kernel_size,\n",
    "            use_bias=self.use_bias,\n",
    "            bias_initializer=self.bias_initializer,\n",
    "            bias_regularizer=self.bias_regularizer,\n",
    "            kernel_initializer=self.kernel_initializer,\n",
    "            kernel_regularizer=self.kernel_regularizer,\n",
    "            padding=\"causal\",\n",
    "            dilation_rate=self.dilation_rate,\n",
    "            activation=NONLINEAR_ACTIVATION,\n",
    "            name=f\"Conv1D_1_{self.layer_id}\"\n",
    "        )\n",
    "        # Spatial dropout is specific to convolutions by dropping an entire timewindow,\n",
    "        # not to rely too heavily on specific features detected by the kernels.\n",
    "        self.dropout1 = SpatialDropout1D(\n",
    "            self.dropout_rate, trainable=True, name=f\"SpatialDropout1D_1_{self.layer_id}\"\n",
    "        )\n",
    "        # Capture a higher order feature set from the previous convolution\n",
    "        self.conv2 = Conv1D(\n",
    "            filters=self.filters,\n",
    "            kernel_size=self.kernel_size,\n",
    "            use_bias=self.use_bias,\n",
    "            bias_initializer=self.bias_initializer,\n",
    "            bias_regularizer=self.bias_regularizer,\n",
    "            kernel_initializer=self.kernel_initializer,\n",
    "            kernel_regularizer=self.kernel_regularizer,\n",
    "            padding=\"causal\",\n",
    "            dilation_rate=self.dilation_rate,\n",
    "            activation=NONLINEAR_ACTIVATION,\n",
    "            name=f\"Conv1D_2_{self.layer_id}\"\n",
    "        )\n",
    "        self.dropout2 = SpatialDropout1D(\n",
    "            self.dropout_rate, trainable=True, name=f\"SpatialDropout1D_2_{self.layer_id}\"\n",
    "        )\n",
    "        # The skip connection is an addition of the input to the block with the output of the second dropout layer.\n",
    "        # Solves vanishing gradient, carries info from earlier layers to later layers, allowing gradients to flow across this alternative path.\n",
    "        # Does not learn direct mappings, but differences (residuals) while keeping temporal context.\n",
    "        # Note how it keeps dims intact with kernel 1.\n",
    "        self.skip_out = Conv1D(\n",
    "            filters=self.filters,\n",
    "            kernel_size=1,\n",
    "            activation=\"relu\",\n",
    "            name=f\"Conv1D_skipconnection_{self.layer_id}\",\n",
    "        )\n",
    "        # This is the elementwise add for the residual connection and Conv1d 2's output\n",
    "        self.residual_out = Add(name=f\"residual_Add_{self.layer_id}\")\n",
    "\n",
    "    def call(self, inputs):\n",
    "        x = self.conv1(inputs)\n",
    "        x = self.dropout1(x)\n",
    "        x = self.conv2(x)\n",
    "        x = self.dropout2(x)\n",
    "\n",
    "        # Residual output by adding the inputs back\n",
    "        skip_out_x = self.skip_out(inputs)\n",
    "        x = self.residual_out([x, skip_out_x])\n",
    "        return x, skip_out_x\n",
    "\n",
    "class ConditionalBlock(Layer):\n",
    "    \"\"\"\n",
    "    TCN condtioning Block that conditions a target timeseries to exogenous timeserieses.\n",
    "    The Block is obtained by stacking togeather the following:\n",
    "        - 1D Dilated Convolution for the main TS.\n",
    "        - 1D Dilated Convolution for the exog TSs.\n",
    "        - 1D Dilated skip layer for both to retain history.\n",
    "        - ReLu\n",
    "        - Spatial Dropout\n",
    "    And adding the input after trasnforming it with a 1x1 Conv\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(\n",
    "        self,\n",
    "        filters=1,\n",
    "        kernel_size=2,\n",
    "        kernel_initializer=\"glorot_normal\",\n",
    "        bias_initializer=\"glorot_normal\",\n",
    "        kernel_regularizer=None,\n",
    "        bias_regularizer=None,\n",
    "        use_bias=False,\n",
    "        dropout_rate=0.01,\n",
    "        layer_id=None,\n",
    "        **kwargs,\n",
    "    ):\n",
    "        super(ConditionalBlock, self).__init__(**kwargs)\n",
    "\n",
    "        assert filters > 0 and kernel_size > 0\n",
    "\n",
    "        self.filters = filters\n",
    "        self.kernel_size = kernel_size\n",
    "        self.kernel_initializer = kernel_initializer\n",
    "        self.bias_initializer = bias_initializer\n",
    "        self.kernel_regularizer = kernel_regularizer\n",
    "        self.bias_regularizer = bias_regularizer\n",
    "        self.use_bias = use_bias\n",
    "        self.dropout_rate = dropout_rate\n",
    "        self.layer_id = str(layer_id)\n",
    "\n",
    "    def get_config(self):\n",
    "        config = super(ConditionalBlock, self).get_config()\n",
    "        config.update({\n",
    "            'filters': self.filters,\n",
    "            'kernel_size': self.kernel_size,\n",
    "            'kernel_initializer': self.kernel_initializer,\n",
    "            'bias_initializer': self.bias_initializer,\n",
    "            'kernel_regularizer': self.kernel_regularizer,\n",
    "            'bias_regularizer': self.bias_regularizer,\n",
    "            'use_bias': self.use_bias,\n",
    "            'dropout_rate': self.dropout_rate,\n",
    "            'id': self.layer_id\n",
    "        })\n",
    "        return config\n",
    "\n",
    "    def build(self, inputs):\n",
    "        self.main_conv = Conv1D(\n",
    "            filters=self.filters,\n",
    "            kernel_size=self.kernel_size,\n",
    "            use_bias=self.use_bias,\n",
    "            bias_initializer=self.bias_initializer,\n",
    "            bias_regularizer=self.bias_regularizer,\n",
    "            kernel_initializer=self.kernel_initializer,\n",
    "            kernel_regularizer=self.kernel_regularizer,\n",
    "            padding=\"causal\",\n",
    "            activation=NONLINEAR_ACTIVATION,\n",
    "            name=f\"Conv1D_Conditional_1\",\n",
    "        )\n",
    "        self.dropout1 = SpatialDropout1D(\n",
    "            self.dropout_rate, trainable=True, name=f\"SpatialDropout1D_1_{self.layer_id}\"\n",
    "        )\n",
    "        self.main_skip_conn = Conv1D(\n",
    "            filters=self.filters,\n",
    "            kernel_size=1,\n",
    "            activation=\"relu\",\n",
    "            name=f\"Skip_Conditional_1\",\n",
    "        )\n",
    "        self.cond_conv = Conv1D(\n",
    "            filters=self.filters,\n",
    "            kernel_size=self.kernel_size,\n",
    "            use_bias=self.use_bias,\n",
    "            bias_initializer=self.bias_initializer,\n",
    "            bias_regularizer=self.bias_regularizer,\n",
    "            kernel_initializer=self.kernel_initializer,\n",
    "            kernel_regularizer=self.kernel_regularizer,\n",
    "            padding=\"causal\",\n",
    "            activation=NONLINEAR_ACTIVATION,\n",
    "            name=f\"Conv1D_Conditional_2\",\n",
    "        )\n",
    "        self.cond_skip_conn = Conv1D(\n",
    "            filters=self.filters,\n",
    "            kernel_size=1,\n",
    "            activation=\"relu\",\n",
    "            name=f\"Skip_Conditional_2\",\n",
    "        )\n",
    "        self.dropout2 = SpatialDropout1D(\n",
    "            self.dropout_rate, trainable=True, name=f\"SpatialDropout1D_2_{self.layer_id}\"\n",
    "        )\n",
    "\n",
    "    def call(self, inputs):\n",
    "        \"\"\"\n",
    "        Will apply causal convolutions to every TS and concatenate the results.\n",
    "        :param inputs: Array\n",
    "            A list where inputs[0] is the main input and inputs[1] are the conditional inputs\n",
    "        :return: Array\n",
    "            Tensor of concatenated results.\n",
    "        \"\"\"\n",
    "        main_input, cond_input = inputs[0], inputs[1] if len(inputs) > 1 else None\n",
    "\n",
    "        x = self.main_conv(main_input)\n",
    "        x = self.dropout1(x)\n",
    "        skip_out_x = self.main_skip_conn(main_input)\n",
    "        x = Add()([x, skip_out_x])\n",
    "        if cond_input is not None:\n",
    "            cond_x = self.cond_conv(cond_input)\n",
    "            cond_x = self.dropout2(cond_x)\n",
    "            cond_skip_out_x = self.cond_skip_conn(cond_input)\n",
    "            cond_x = Add()([cond_x, cond_skip_out_x])\n",
    "\n",
    "            x = Concatenate(axis=-1)([x, cond_x])\n",
    "        return x\n",
    "\n",
    "\n",
    "def TCN(\n",
    "    input_shape,\n",
    "    dense_units=None,\n",
    "    conditioning_shapes=None,\n",
    "    output_horizon=1,\n",
    "    filters=[32],\n",
    "    kernel_size=2,\n",
    "    dilation_rate=2,\n",
    "    kernel_initializer=\"glorot_normal\",\n",
    "    bias_initializer=\"glorot_normal\",\n",
    "    kernel_regularizer=None,\n",
    "    bias_regularizer=None,\n",
    "    use_bias=False,\n",
    "    dropout_rate=0.01,\n",
    "):\n",
    "    \"\"\"\n",
    "    Tensorflow TCN Model builder.\n",
    "    see: https://www.tensorflow.org/api_docs/python/tf/keras/Model\n",
    "    see: https://www.tensorflow.org/guide/keras/making_new_layers_and_models_via_subclassing#the_model_class\n",
    "    see: https://www.tensorflow.org/api_docs/python/tf/keras/regularizers/L2\n",
    "\n",
    "    :param layers: int\n",
    "        Number of layers for the network. Defaults to 1 layer.\n",
    "    :param filters: int\n",
    "        the number of output filters in the convolution. Defaults to 32.\n",
    "    :param kernel_size: int or tuple\n",
    "        the length of the 1D convolution window\n",
    "    :param dilation_rate: int\n",
    "        the dilation rate to use for dilated convolution. Defaults to 1.\n",
    "    :param output_horizon: int\n",
    "        the output horizon.\n",
    "    \"\"\"\n",
    "    main_input = Input(shape=input_shape, name=\"main_input\")\n",
    "    cond_input = (\n",
    "        Input(shape=conditioning_shapes, name=\"exog_input\")\n",
    "        if conditioning_shapes is not None and len(conditioning_shapes) > 0\n",
    "        else None\n",
    "    )\n",
    "    x = main_input\n",
    "    if cond_input is not None:\n",
    "        x = ConditionalBlock(\n",
    "            filters=filters[0],\n",
    "            kernel_size=kernel_size,\n",
    "            kernel_initializer=kernel_initializer,\n",
    "            bias_initializer=bias_initializer,\n",
    "            kernel_regularizer=kernel_regularizer,\n",
    "            bias_regularizer=bias_regularizer,\n",
    "            use_bias=use_bias,\n",
    "            dropout_rate=dropout_rate,\n",
    "        )([main_input] + [cond_input])\n",
    "    skip_connections = []\n",
    "    for i, filter in enumerate(filters):\n",
    "        x, x_skip = TCNBlock(\n",
    "            filters=filter,\n",
    "            kernel_size=kernel_size,\n",
    "            dilation_rate=dilation_rate ** (i + 1),\n",
    "            kernel_initializer=kernel_initializer,\n",
    "            bias_initializer=bias_initializer,\n",
    "            kernel_regularizer=kernel_regularizer,\n",
    "            bias_regularizer=bias_regularizer,\n",
    "            use_bias=use_bias,\n",
    "            dropout_rate=dropout_rate,\n",
    "            layer_id=i,\n",
    "        )(x)\n",
    "        skip_connections.append(x_skip)\n",
    "    if skip_connections:\n",
    "        skip_connections.append(x)\n",
    "        aggregated = Concatenate(axis=-1, name=f\"Final_Residuals\")(skip_connections)\n",
    "        aggregated = Conv1D(filters[-1], kernel_size=1, activation='relu', padding='same')(aggregated)\n",
    "\n",
    "\n",
    "    if dense_units:\n",
    "        # Dense networks for deep learning ifrequired.\n",
    "        x = Flatten()(x)\n",
    "        # First layer\n",
    "        x = Dense(dense_units[0], input_shape=input_shape, activation=\"LeakyReLU\", name=f\"Dense_0\")(x)\n",
    "        for i, units  in enumerate(dense_units, start=1):\n",
    "            x = Dense(units , activation=\"LeakyReLU\", name=f\"Dense__{i}\")(x)\n",
    "        # Last layer\n",
    "        x = Dense(input_shape[0], activation=\"sigmoid\", name=f\"Dense_Classifier\")(x)\n",
    "    else:\n",
    "        x = Conv1D(filters=output_horizon, kernel_size=1, padding=\"causal\", activation=\"sigmoid\",name=f\"Conv_Classifier\")(x)\n",
    "    model = Model(\n",
    "        inputs=[main_input, cond_input] if cond_input is not None else [main_input],\n",
    "        outputs=x,\n",
    "        name=\"TCN_Conditional_Model\",\n",
    "    )\n",
    "\n",
    "    return model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30e8329e",
   "metadata": {},
   "source": [
    "With tensorflow2 we built the Temporal Convolution Network (TCN), a convolution network’s operations over consecutive layers can be represented with the following :\n",
    "$$\n",
    "x_{l}^{t} = g\\left(\\sum_{k=0}^{K-1} w_{lk} \\cdot x_{l-1}^{t-(k \\times d)} + b_l\\right)\n",
    "$$\n",
    "Where:\n",
    "- xt_l is the output of the neuron at position t in the l-th layer.\n",
    "- K is the kernel’s with, determining the number of past time windows considered.\n",
    "- wlk is the weight for the k-th position in the kernel used to give importance to past data.\n",
    "- d is the dilation factor, or the space between inputs allowing the network to integrate historical information.\n",
    "- bl is the bias term.\n",
    "- g is ReLU defined as $ g(x) = \\max(0, x) $.\n",
    "\n",
    "The NN has some special attributes worth talking about, these being dilations and causal layers, which can be visualized by this image taken from the WaveNet paper:\n",
    "\n",
    "![image](./images/wavenet.png)\n",
    "\n",
    "and compare it to what a typical convolution NN would look like:\n",
    "\n",
    "![image](./images/wavenet_normal.png)\n",
    "\n",
    "## Dilated Causal Convolutions\n",
    "\n",
    "From the images above, note how at each layer of neurons, the connections skip a neuron or two - this a **dilation factor**. Dilations in the TCN increase exponentially its spectral field within the depth of the network, allowing the neurons to capture history, without leakage.\n",
    "\n",
    "Additionally, we observe how the neurons are not fully-connected, or they are connected to the previous neuron but not to the next - this is what is know as a **causal padding**, meaning no future information is leaked in to the current neuron that is the 'now'.\n",
    "\n",
    "## Residuals and Skip Connections\n",
    "\n",
    "Skip connections are a technique from Residual Networks (ResNets), that were used to mitigate the vanishing gradient problem and improve model performance. A skip connection allows the gradient to be directly backpropagated to earlier layers. In our TCN architecture, skip connections are used to connect the output of each convolutional block directly to a final layer in the network (before the output layer) allowing the network to make use of features learned at various depths directly when making predictions. This is particularly useful for time series data where both recent, and more distant past information might be relevant for prediction.\n",
    "\n",
    "The image from WaveNet illustrates this concept:\n",
    "\n",
    "\n",
    "![image](./images/skip_layers.png)\n",
    "\n",
    "While their function are similar, these connections have slightly different roles:\n",
    "- Residual Connections is used within blocks to add the block's input to its output for gradient flow.\n",
    "- Skip Connections skip multiple layers or blocks and are used to directly connect early layers to later ones, allowing the network to leverage features learned at various stages of the network directly in the output.\n",
    "\n",
    "\n",
    "# Train Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ba82142",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import shutil\n",
    "from tensorflow.config.experimental import list_physical_devices, set_memory_growth\n",
    "\n",
    "tf.keras.backend.clear_session()\n",
    "\n",
    "# https://learn.microsoft.com/en-us/windows/ai/directml/gpu-tensorflow-plugin\n",
    "gpus = list_physical_devices(\"GPU\")\n",
    "print(\"Num GPUs Available: \", len(gpus))\n",
    "if gpus:\n",
    "    try:\n",
    "        for gpu in gpus:\n",
    "            set_memory_growth(gpu, True)\n",
    "    except RuntimeError as e:\n",
    "        print(e)\n",
    "\n",
    "MODEL_DIR = f\"./models/{datetime.now().strftime('%Y%m%d')}\"\n",
    "DELETE_OLD_LOGS = True\n",
    "\n",
    "LOG_BASEPATH = \"./logs\"\n",
    "if DELETE_OLD_LOGS and os.path.exists(LOG_BASEPATH):\n",
    "    assert os.path.isdir(LOG_BASEPATH)\n",
    "    shutil.rmtree(LOG_BASEPATH)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea9a3ae9",
   "metadata": {},
   "source": [
    "## Tensorboard for profiling\n",
    "\n",
    "Time to train and evaluate our model. We'll use tensorbards, and run it with the command: `tensorboard --logdir ./logs`, or run the following code in your notebook `%tensorboard --logdir logs` for in-notebook magic (which I don't recommend outside of a Collab environment). Addtionally, we have to provide calls backs for tensorflow to feed logs to tensorboards."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06a2ca2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tensorboard import program\n",
    "\n",
    "tb = program.TensorBoard()\n",
    "tb.configure(argv=[None, '--logdir', LOG_BASEPATH, '--bind_all'])\n",
    "url = tb.launch()\n",
    "print(f\"TensorBoard started at {url}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f8ff228",
   "metadata": {},
   "source": [
    "## Hyper Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22d53654",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tensorboard.plugins.hparams import api as hp\n",
    "from tensorflow.keras.losses import BinaryCrossentropy, BinaryFocalCrossentropy\n",
    "from tensorflow.keras.metrics import AUC, BinaryAccuracy\n",
    "from datetime import datetime\n",
    "from sklearn.model_selection import ParameterGrid, TimeSeriesSplit\n",
    "import json\n",
    "import pickle\n",
    "\n",
    "OOS_SPLIT = 0.1\n",
    "VAL_SPLIT = 0.20\n",
    "EPOCHS = 300\n",
    "PATIENCE_EPOCHS = 15\n",
    "BATCH_SIZE = 124\n",
    "MAX_FILTER = 512\n",
    "MIN_FILTER = 32\n",
    "FILTERS = [MIN_FILTER, MIN_FILTER*4] # [WINDOW_SIZE//2, WINDOW_SIZE, WINDOW_SIZE*2]\n",
    "HIDDEN_DENSE = [WINDOW_SIZE] # [WINDOW_SIZE//2, PREDICTION_HORIZON]\n",
    "BIAS = True\n",
    "DROPRATE = 0.05\n",
    "POOL_SIZE = 8\n",
    "KERNEL_SIZE = 2\n",
    "DILATION_RATE = 1\n",
    "MAX_LAYERS = 2\n",
    "REG_WEIGHTS = 0.0001\n",
    "LEARN_RATE = 0.0025\n",
    "MODEL_LOG_DIR = f'{LOG_BASEPATH}/{datetime.now().strftime(\"%d%H%M%S\")}'\n",
    "TARGET_METRIC = \"auc\"\n",
    "LOSS = BinaryFocalCrossentropy(apply_class_balancing=True, from_logits=True)\n",
    "METRICS = [AUC(name=TARGET_METRIC), BinaryCrossentropy(from_logits=True), BinaryAccuracy(), 'accuracy']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82617c46",
   "metadata": {},
   "source": [
    "## Model Training\n",
    "\n",
    "Some noteworthy configrations in the training:\n",
    "- `EarlyStopping`: We configured a long 300-epoch horizon for training, this will not be used in its full size for our training. If EarlyStopping detects a flattening gradient, it will stop the training and re-esteblish the better weights. In most of our runs it stopped between 15 to 20 epochs.\n",
    "- `ReduceLROnPlateau`: In addition to early stopping, we also want to decrease the learning alpha as the gradient start to plateauing to help it converge. We start with a large learning rate, and quarter it everytime progress stops, until the early stopping kicks in.\n",
    "- `TensorBoard`: The callback to give logs and data to tensorboard.\n",
    "- And a final `LambdaCallback` to print a confusion matrix on every epoch, to allow us to gauge the learning and classification errors.\n",
    "\n",
    "Let's talk about the error function we chose to optimize: `BinaryFocalCrossentropy`.\n",
    "\n",
    "The idea behind **Focal Loss, and Binary Crossentropy**, is to dynamically scale the loss based on the correctness of the classification. It applies a modulating factor to the cross-entropy loss, reducing the contributions from easy examples and focusing the model on hard examples. For a binary classification problem, the Binary Focal Loss is defined as:\n",
    "\n",
    "$$\n",
    "FL(p_t) = -\\alpha_t (1 - p_t)^\\gamma \\log(p_t)\n",
    "$$\n",
    "Where:\n",
    "- \\(p_t\\) is the model's predicted probability for the class with label \\(y=1\\). Specifically, \\(p_t = p\\) if \\(y=1\\), and \\(p_t = 1 - p\\) otherwise, where \\(p\\) is the predicted probability for the class with label \\(1\\).\n",
    "- \\(\\alpha_t\\) is a weighting factor for the class with label \\(1\\), which helps to manage class imbalance by assigning more weight to the less frequent class.\n",
    "- \\(\\gamma\\) is the focusing parameter that adjusts the rate at which easy examples are down-weighted. When \\(\\gamma = 0\\), Focal Loss reduces to binary cross-entropy loss.\n",
    "- The term \\((1 - p_t)^\\gamma\\) effectively reduces the loss contribution from well-classified examples (\\(p_t \\approx 1\\) for true positives and \\(p_t \\approx 0\\) for true negatives), putting more focus on examples that are misclassified or classified with low confidence.\n",
    "- The \\(\\log(p_t)\\) term is the log-likelihood of the true class, similar to what is used in cross-entropy loss.\n",
    "- \n",
    "The loss function focuses on Hard Examples only, by reducing the loss contribution from easy-to-classify examples, these being Bearish signals (which are a majority class in our imbalanced dataset), and focus on the hard, often misclassified examples (typically the minority class) which are our profitable bullish signals. This helps in learning discriminative features for the minority class without being overwhelmed by the majority class.\n",
    "\n",
    "Another important metric that we'll include is the **Receiver Operating Characteristic (ROC) Curve**. The models' training process will provide the area under the ROC curve, which is the True Positive Rate (TPR, also known as sensitivity or recall) against the False Positive Rate (FPR, 1 - specificity) at various threshold settings. The image from wikipedia can help in disambiguating these concepts:\n",
    "\n",
    "![image](./images/Precisionrecall.svg.png)\n",
    "\n",
    "The integration of the following make up the ROC AUC metric:\n",
    "$$\n",
    "\\text{True Positive Rate (TPR)} = \\text{Recall} = \\frac{TP}{TP + FN}\n",
    "\\newline\n",
    "\\text{False Positive Rate (FPR)} = \\frac{FP}{FP + TN}\n",
    "$$\n",
    "\n",
    "The AUC measures the entire two-dimensional area underneath the curve from (0,0) to (1,1):\n",
    "- AUC = 1: The classifier is perfect, able to distinguish all positive and negative instances correctly.\n",
    "- 0.5 < AUC < 1: The classifier is better than random guessing, with higher values indicating better classification performance.\n",
    "- AUC = 0.5: The classifier's performance is no better than random guessing.\n",
    "- AUC < 0.5: The classifier is performing worse than random guessing. In practice, this situation is uncommon since it indicates that one could improve the classifier by simply inverting its decisions.\n",
    "\n",
    "It's important that True Positive (TP) accuracy is high enough, since we cannot offord monetary losses on misclassification. This metric is common in financial and medical machine learning problems, where mistakes are costly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f481f001",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tensorflow.keras.regularizers import L2, L1L2\n",
    "from tensorflow.keras.callbacks import EarlyStopping, TensorBoard,ModelCheckpoint,ReduceLROnPlateau,LambdaCallback\n",
    "from tensorflow.keras.optimizers import Adam\n",
    "from tensorflow.summary import create_file_writer\n",
    "from tensorflow.debugging.experimental import enable_dump_debug_info\n",
    "from tensorflow.math import confusion_matrix\n",
    "import io\n",
    "\n",
    "def plot_confusion_matrix(cm, labels, cm2=None, labels2=None):\n",
    "        plt.figure(figsize=(8 if cm2 is not None else 4, 4))\n",
    "        if cm2 is not None:\n",
    "            plt.subplot(1, 2, 1)\n",
    "        plt.imshow(cm, interpolation='nearest', cmap=plt.cm.Accent)\n",
    "\n",
    "        df_cm = pd.DataFrame((cm / np.sum(cm, axis=1)[:, None])*100, index=[i for i in labels], columns=[i for i in labels])\n",
    "        cm_plot1 = sns.heatmap(df_cm, annot=True,  fmt=\".2f\", cmap='Blues', xticklabels=labels, yticklabels=labels).get_figure()\n",
    "        plt.xlabel('Predicted Labels')\n",
    "        plt.ylabel('True Labels')\n",
    "        plt.title('Confusion Matrix 1')\n",
    "        tick_marks = np.arange(len(labels))\n",
    "        plt.xticks(tick_marks, labels, rotation=45)\n",
    "        plt.yticks(tick_marks, labels)\n",
    "\n",
    "        cm_plot2=None\n",
    "        if cm2 is not None:\n",
    "            plt.subplot(1, 2, 2)\n",
    "            df_cm = pd.DataFrame((cm2 / np.sum(cm2, axis=1)[:, None])*100, index=[i for i in labels2], columns=[i for i in labels2])\n",
    "            cm_plot12 = sns.heatmap(df_cm, annot=True,  fmt=\".2f\", cmap='Reds', xticklabels=labels, yticklabels=labels).get_figure()\n",
    "            plt.xlabel('Predicted Labels')\n",
    "            plt.title('Confusion Matrix 2')\n",
    "        plt.tight_layout()\n",
    "\n",
    "        return cm_plot1, cm_plot2\n",
    "\n",
    "# enable_dump_debug_info(LOG_BASEPATH, tensor_debug_mode=\"FULL_HEALTH\", circular_buffer_size=-1)\n",
    "\n",
    "def build_tcn(\n",
    "    input_shape,\n",
    "    X, y,\n",
    "    Xt=None, yt=None,\n",
    "    conditioning_shapes=None,\n",
    "    val_split=VAL_SPLIT,\n",
    "    output_horizon=PREDICTION_HORIZON,\n",
    "    filters=FILTERS,\n",
    "    kernel_size=KERNEL_SIZE,\n",
    "    dilation_rate=DILATION_RATE,\n",
    "    kernel_regularizer=L1L2(l1=REG_WEIGHTS, l2=REG_WEIGHTS//10),\n",
    "    bias_regularizer=L1L2(l1=REG_WEIGHTS, l2=REG_WEIGHTS//10),\n",
    "    dropout_rate=DROPRATE,\n",
    "    dense_units=HIDDEN_DENSE,\n",
    "    lr=LEARN_RATE,\n",
    "    patience=PATIENCE_EPOCHS,\n",
    "    epochs=EPOCHS,\n",
    "    batch_size=BATCH_SIZE,\n",
    "    use_bias=BIAS,\n",
    "    loss=LOSS,\n",
    "    tb=True,\n",
    "):\n",
    "    def log_confusion_matrix(epoch, logs):\n",
    "        def _plot_to_image(figure):\n",
    "            \"\"\"Converts the matplotlib plot specified by 'figure' to a PNG image and\n",
    "            returns it. The supplied figure is closed and inaccessible after this call.\"\"\"\n",
    "            buf = io.BytesIO()\n",
    "            plt.savefig(buf, format='png')\n",
    "            plt.close(figure)\n",
    "            buf.seek(0)\n",
    "            image = tf.image.decode_png(buf.getvalue(), channels=4)\n",
    "            image = tf.expand_dims(image, 0)\n",
    "            return image\n",
    "        # model is global as is XT and yt\n",
    "        assert Xt is not None and len(Xt) > 1 and len(yt) > 1\n",
    "        ypred = model.predict(Xt)\n",
    "        cm = confusion_matrix(yt.flatten(), ypred.flatten())\n",
    "        figure, _ = plot_confusion_matrix(cm, labels=[1,0])\n",
    "        cm_image = _plot_to_image(figure)\n",
    "\n",
    "        file_writer_cm = create_file_writer(LOG_BASEPATH)\n",
    "        with file_writer_cm.as_default():\n",
    "            tf.summary.image(\"Confusion Matrix\", cm_image, step=epoch)\n",
    "\n",
    "    assert len(X) > 1 and len(y) > 1 and input_shape is not None\n",
    "    global model, globalXt, globalyt\n",
    "    globalXt = Xt\n",
    "    globalyt = yt\n",
    "\n",
    "    model = TCN(\n",
    "        input_shape=input_shape,\n",
    "        conditioning_shapes=conditioning_shapes,\n",
    "        dense_units=dense_units,\n",
    "        output_horizon=output_horizon,\n",
    "        filters=filters,\n",
    "        kernel_size=kernel_size,\n",
    "        dilation_rate=dilation_rate,\n",
    "        kernel_regularizer=kernel_regularizer,\n",
    "        bias_regularizer=bias_regularizer,\n",
    "        use_bias=use_bias,\n",
    "        dropout_rate=dropout_rate,\n",
    "    )\n",
    "\n",
    "    model.compile(loss=loss, optimizer=Adam(lr), metrics=METRICS)\n",
    "    callbacks = [EarlyStopping(\n",
    "                    patience=patience,\n",
    "                    monitor=f\"val_{TARGET_METRIC}\",\n",
    "                    restore_best_weights=True,\n",
    "                ),\n",
    "                ReduceLROnPlateau(\n",
    "                    monitor=f\"val_{TARGET_METRIC}\",\n",
    "                    factor=0.3,\n",
    "                    patience=patience//2,\n",
    "                    verbose=1,\n",
    "                    min_delta=0.00001,\n",
    "                )]\n",
    "    if tb:\n",
    "        callbacks.append(TensorBoard(log_dir=MODEL_LOG_DIR,\n",
    "                                    histogram_freq=1,\n",
    "                                    write_graph=True,\n",
    "                                    write_images=True,\n",
    "                                    update_freq='epoch',\n",
    "                                    profile_batch=2,\n",
    "                                    embeddings_freq=1))\n",
    "    if tb and IS_CLASSIFICATION:\n",
    "        callbacks.append(LambdaCallback(on_epoch_end=log_confusion_matrix))\n",
    "    if Xt is not None:\n",
    "        assert len(Xt) > 1 and len(yt) > 1\n",
    "        history = model.fit(\n",
    "            X,\n",
    "            y,\n",
    "            validation_data=(Xt, yt),\n",
    "            epochs=epochs,\n",
    "            batch_size=batch_size,\n",
    "            callbacks=callbacks,\n",
    "            verbose=1,\n",
    "        )\n",
    "    else:\n",
    "        assert val_split > 0 and  Xt is None and yt is None\n",
    "        history = model.fit(\n",
    "            X,\n",
    "            y,\n",
    "            validation_split=val_split,\n",
    "            epochs=epochs,\n",
    "            batch_size=batch_size,\n",
    "            callbacks=callbacks,\n",
    "            verbose=0,\n",
    "        )\n",
    "    return model, history"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7485c67",
   "metadata": {},
   "source": [
    "### Constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d5eca9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "LOAD_MODEL_CACHE = False\n",
    "\n",
    "OOS_SIZE = round(len(data_df) * OOS_SPLIT)\n",
    "\n",
    "x_oos = data_df[-OOS_SIZE:]\n",
    "data_is_df = data_df[:-OOS_SIZE]\n",
    "\n",
    "VAL_SIZE = round(len(data_is_df) * VAL_SPLIT)\n",
    "\n",
    "data_x_df = data_is_df[:-VAL_SIZE]\n",
    "data_t_df = data_is_df[-VAL_SIZE:]\n",
    "\n",
    "print(f\"{x_oos.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01f5fc3c",
   "metadata": {},
   "source": [
    "## Encoding Time Windows\n",
    "\n",
    "Given all the special attributes we described for the TCN, you shouldn't be surpised that the data needs special threatment, it needs to be encoded into what is known as **time windows or time steps**. The image from Gasparin et al. gives a good visualization of these time steps:\n",
    "\n",
    "![image](./images/timewindows.jpg)\n",
    "\n",
    "The data was normalised and encoded into 3months-long windows, with a one-day lag from so we can predict the market direction of the next day:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54750ecd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def prepare_windows(\n",
    "    data_df,\n",
    "    target_df,\n",
    "    prime_ts=SELECTED_FEATURES,\n",
    "    exog_ts=SELECTED_EXOG,\n",
    "    window_size=WINDOW_SIZE,\n",
    "    horizon=PREDICTION_HORIZON,\n",
    "):\n",
    "    \"\"\"\n",
    "    Create input and target windows suitable for TCN model.\n",
    "    :param data: DataFrame with shape (n_samples, n_features)\n",
    "    :param features: List of strings, names of the feature columns\n",
    "    :param target_df: Optional the labels if this encoding is for training.\n",
    "    :param window_size: int, length of the input sequence.\n",
    "    :param horizon: int, forecasting horizon.\n",
    "    :return: 3 Arrays in the shape of (n_samples, n_steps, n_features) for the training data, the exogenous data, and the labels (this last is optional)\n",
    "    \"\"\"\n",
    "    X, Xexog, y = [], [], []\n",
    "    for i in tqdm(\n",
    "        range(len(data_df) - window_size - horizon + 1), desc=f\"Encoding Widows of {window_size}\"\n",
    "    ):\n",
    "        input_window = data_df[prime_ts].iloc[i : i + window_size].values\n",
    "        X.append(input_window)\n",
    "        input_window = data_df[exog_ts].iloc[i : i + window_size].values\n",
    "        Xexog.append(input_window)\n",
    "        if target_df is not None:\n",
    "            target_window = target_df.iloc[i  : i + window_size].values\n",
    "            y.append(target_window)\n",
    "    return np.array(X), np.array(Xexog), np.array(y)\n",
    "\n",
    "train_data, train_exog_data, ytrain_data = prepare_windows(data_x_df, data_x_df[TARGET_LABEL])\n",
    "test_data, test_exog_data, ytest_data = prepare_windows(data_t_df, data_t_df[TARGET_LABEL])\n",
    "\n",
    "assert not np.any(pd.isna(train_data)) and not np.any(pd.isna(train_exog_data))\n",
    "\n",
    "print(f\"Label shape encoded: {ytrain_data.shape}\")\n",
    "print(f\"Data shapes for prime TS: {train_data.shape}, exog TS: {train_exog_data.shape}\")\n",
    "print(f\"First window: {train_exog_data[:1][0]}\")\n",
    "print(f\"First window: {train_data[:1][0]}\")\n",
    "print(f\"First window targets: {ytrain_data[:1][0]}\")\n",
    "\n",
    "input_shape = (\n",
    "    WINDOW_SIZE,\n",
    "    1 if len(train_data.shape) < 3 else train_data.shape[2],\n",
    ")  # if we have no additonal features X.shape[1]\n",
    "conditioning_shapes = (WINDOW_SIZE, train_exog_data.shape[2])\n",
    "print(f\"Model logs for Tensorboard available here: {MODEL_LOG_DIR}\")\n",
    "print(f\"Input Shape: {input_shape} and Condtioning shapes: {conditioning_shapes}\")\n",
    "\n",
    "assert not np.any(np.isnan(train_data))\n",
    "assert not np.any(np.isnan(train_exog_data))\n",
    "assert not np.any(np.isnan(ytrain_data))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d18f078d",
   "metadata": {},
   "source": [
    "Before we train, we will set aside some out-of-sample data which the model will not see. This will be used to test the model's generalization later on:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91b14747",
   "metadata": {},
   "outputs": [],
   "source": [
    "xtrain_oos, Xexog_oos, y_oos = prepare_windows(x_oos, x_oos[TARGET_LABEL])\n",
    "y_oos"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a14dda81",
   "metadata": {},
   "source": [
    "### GridSearch and HyperParams\n",
    "\n",
    "The model is complex, and to manually fine-tune every parameter is impossible. We'll run a grid search, which after a day or so of searching, should identify good enough hyperparameters. In our tests, we got this setup:\n",
    "\n",
    "```json\n",
    "{\n",
    "    \"best_params\": {\n",
    "        \"batch_size\": 124,\n",
    "        \"bias\": true,\n",
    "        \"dense_units\": \"126\",\n",
    "        \"dilation_rate\": 1,\n",
    "        \"dropout_rate\": 0.05,\n",
    "        \"epochs\": 300,\n",
    "        \"filters\": \"32_128\",\n",
    "        \"kernel_size\": 2,\n",
    "        \"learning_rate\": 0.0025,\n",
    "        \"patience\": 15,\n",
    "        \"reg_weight\": 0.0002\n",
    "    },\n",
    "    \"best_loss\": 0.00413,\n",
    "    \"best_metric\": 0.8180,\n",
    "}\n",
    "```\n",
    "\n",
    "The x_x_x hparams are our way of encoding an array in tensorflow's hparam framework. These need to be decoded before being passed to the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d598dfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# See: https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.ParameterGrid.html\n",
    "# See paper: https://www.mdpi.com/2076-3417/10/7/2322\n",
    "\n",
    "HP_KERNEL_SIZE = hp.HParam(\"kernel_size\", hp.Discrete([KERNEL_SIZE * 2, KERNEL_SIZE]))\n",
    "HP_BATCH_SIZE = hp.HParam(\"batch_size\", hp.Discrete([BATCH_SIZE]))\n",
    "HP_EPOCHS = hp.HParam(\"epochs\", hp.Discrete([EPOCHS]))\n",
    "HP_DILATION_RATE = hp.HParam(\"dilation_rate\", hp.Discrete([DILATION_RATE]))\n",
    "HP_DROPOUT_RATE = hp.HParam(\"dropout_rate\", hp.Discrete([DROPRATE, DROPRATE*2]))\n",
    "HP_REG_WEIGHTS = hp.HParam(\"reg_weight\", hp.Discrete([REG_WEIGHTS, REG_WEIGHTS*2]))\n",
    "HP_LEARNING_RATE = hp.HParam(\"learning_rate\", hp.Discrete([LEARN_RATE]))\n",
    "HP_PATIENCE = hp.HParam(\"patience\", hp.Discrete([PATIENCE_EPOCHS]))\n",
    "HP_BIAS = hp.HParam(\"bias\", hp.Discrete([BIAS, False]))\n",
    "HP_HIDDEN_DENSE = hp.HParam(\"dense_units\", hp.Discrete([\n",
    "    f\"{WINDOW_SIZE}\",\n",
    "    f\"{WINDOW_SIZE*2}_{WINDOW_SIZE}\",\n",
    "]))\n",
    "HP_FILTERS = hp.HParam(\"filters\", hp.Discrete([\n",
    "    f\"{MIN_FILTER}\",\n",
    "    f\"{MIN_FILTER}_{MIN_FILTER*2}\",\n",
    "    f\"{MIN_FILTER*2}_{MIN_FILTER*4}\",\n",
    "]))\n",
    "\n",
    "HPARAMS = [\n",
    "    HP_FILTERS,\n",
    "    HP_KERNEL_SIZE,\n",
    "    HP_BATCH_SIZE,\n",
    "    HP_EPOCHS,\n",
    "    HP_DILATION_RATE,\n",
    "    HP_DROPOUT_RATE,\n",
    "    HP_REG_WEIGHTS,\n",
    "    HP_LEARNING_RATE,\n",
    "    HP_PATIENCE,\n",
    "    HP_BIAS,\n",
    "    HP_HIDDEN_DENSE\n",
    "]\n",
    "\n",
    "GRID_SEARCH_TRAIN = False\n",
    "\n",
    "def grid_search_build_tcn(\n",
    "    input_shape, X, y, Xt=None, yt=None, hparams=HPARAMS, conditioning_shapes=None, file_name=\"best_params.json\"\n",
    "):\n",
    "    def _decode_arrays(config_str):\n",
    "        return [int(unit) for unit in config_str.split('_')]\n",
    "\n",
    "    def _save_best_params(best_params, best_loss, best_metric, file_name=\"best_params.json\"):\n",
    "        os.makedirs(MODEL_DIR, exist_ok=True)\n",
    "        with open(f\"{MODEL_DIR}/{file_name}\", \"w\") as file:\n",
    "            json.dump({\"best_params\": best_params, \"best_loss\": best_loss, \"best_metric\": best_metric}, file)\n",
    "\n",
    "    with create_file_writer(f\"{MODEL_LOG_DIR}/hparam_tuning\").as_default():\n",
    "        hp.hparams_config(\n",
    "            hparams=hparams,\n",
    "            metrics=[hp.Metric(TARGET_METRIC, display_name=TARGET_METRIC)],\n",
    "        )\n",
    "    grid = list(ParameterGrid({h.name: h.domain.values for h in hparams}))\n",
    "    best_model = None\n",
    "    best_loss = np.inf\n",
    "    best_metric = -np.inf\n",
    "    best_params = None\n",
    "    best_history = None\n",
    "    for hp_values in tqdm(grid, desc=\"Grid Search..\"):\n",
    "        dense_units = _decode_arrays(hp_values[\"dense_units\"])\n",
    "        filters = _decode_arrays(hp_values[\"filters\"])\n",
    "        k = hp_values[\"kernel_size\"]\n",
    "        d = hp_values[\"dilation_rate\"]\n",
    "        rw = hp_values[\"reg_weight\"]\n",
    "        drop = hp_values[\"dropout_rate\"]\n",
    "        print(f\"Shapes{input_shape}{conditioning_shapes}: x{X[0].shape}xg{X[1].shape}y{y.shape}, filters {filters}, dense {dense_units}, k: {k}, d: {d}, rw: {rw}, drop: {drop}\")\n",
    "        print(f\"\")\n",
    "\n",
    "        model, history = build_tcn(input_shape, X, y,\n",
    "                                    conditioning_shapes=conditioning_shapes,\n",
    "                                    output_horizon=PREDICTION_HORIZON,\n",
    "                                    Xt=Xt, yt=yt,\n",
    "                                    filters=filters,\n",
    "                                    kernel_size=k,\n",
    "                                    dilation_rate=d,\n",
    "                                    kernel_regularizer=L1L2(l1=rw, l2=rw),\n",
    "                                    bias_regularizer=L1L2(l1=rw, l2=rw),\n",
    "                                    dropout_rate=drop,\n",
    "                                    dense_units=dense_units)\n",
    "        loss = np.min(history.history[f\"val_loss\"])\n",
    "        metric = np.min(history.history[f\"val_{TARGET_METRIC}\"])\n",
    "        if (loss < best_loss) and (best_metric > metric):\n",
    "            print(f\"best metric: {metric}\")\n",
    "            print(f\"best loss: {loss}\")\n",
    "            print(f\"best params: {hp_values}\")\n",
    "            best_history = history\n",
    "            best_loss = loss\n",
    "            best_metric = metric\n",
    "            best_model = model\n",
    "            best_params = hp_values\n",
    "            _save_best_params(best_params, best_loss, best_metric, file_name)\n",
    "    return best_model, best_history\n",
    "\n",
    "if GRID_SEARCH_TRAIN:\n",
    "    model, history = grid_search_build_tcn(input_shape, [train_data, train_exog_data], ytrain_data, Xt=[test_data, test_exog_data], yt=ytest_data, conditioning_shapes=conditioning_shapes)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8007549",
   "metadata": {},
   "source": [
    "## Load or Train\n",
    "\n",
    "With the hyper parameter search complete, we inject our optimal parameters and train the model in one run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92eb87e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tensorflow.keras.models import load_model\n",
    "\n",
    "os.makedirs(MODEL_DIR, exist_ok=True)\n",
    "assert os.path.exists(MODEL_DIR)\n",
    "\n",
    "if LOAD_MODEL_CACHE:\n",
    "    model = load_model(f\"{MODEL_DIR}/tcn.h5\")\n",
    "    with open(f\"{MODEL_DIR}/history.pkl\", 'rb') as file:\n",
    "        loaded_history = pickle.load(file)\n",
    "elif not GRID_SEARCH_TRAIN:\n",
    "     model, history = build_tcn( # Xexog\n",
    "            input_shape, X=[train_data, train_exog_data], y=ytrain_data, Xt=[test_data, test_exog_data], yt=ytest_data, conditioning_shapes=conditioning_shapes\n",
    "        )\n",
    "model.save(f\"{MODEL_DIR}/tcn.h5\")\n",
    "tf.keras.utils.plot_model(\n",
    "    model,\n",
    "    to_file='./images/tcn_model.png',\n",
    "    show_shapes=True,\n",
    "    show_dtype=False,\n",
    "    show_layer_names=True,\n",
    "    rankdir='TB',\n",
    "    expand_nested=True,\n",
    "    dpi=102,\n",
    "    layer_range=None,\n",
    "    show_layer_activations=True,\n",
    "    show_trainable=True,\n",
    ")\n",
    "\n",
    "# model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d34348cd",
   "metadata": {},
   "source": [
    "# Evaluation\n",
    "\n",
    "Initially, we discussed how precision is our most important metric. Another good method of visualizing this, apart from the RoC AUC, is the confussion matrix. In our binary classificiation context, the confusion matrix consists of 2 rows and 2 columns, each representing the instances in a predicted class versus instances in an actual class. Envision a square matrix of these:\n",
    "\n",
    "- True Positives (TP): The number of positive instances correctly classified as positive.\n",
    "- False Positives (FP): The number of negative instances incorrectly classified as positive.\n",
    "- True Negatives (TN): The number of negative instances correctly classified as negative.\n",
    "- False Negatives (FN): The number of positive instances incorrectly classified as negative.\n",
    "\n",
    "And the results of our model below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b862085f",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert not np.array_equal(train_data, test_data) and not np.array_equal(train_exog_data, test_exog_data), \"Training and test should not be identical.\"\n",
    "\n",
    "y_pred = model.predict([train_data, train_exog_data])\n",
    "yt_pred = model.predict([test_data, test_exog_data])\n",
    "\n",
    "y_pred_orig = y_pred.copy()\n",
    "yt_pred_orig = yt_pred.copy()\n",
    "y_pred = (y_pred > 0.5).astype(int)\n",
    "yt_pred = (yt_pred > 0.5).astype(int)\n",
    "\n",
    "print(f\"Prediction shape: {yt_pred.shape} vs test data shape: {ytest_data.shape}\")\n",
    "print(f\"Test data 1 horizon sample: {ytest_data[0]}\")\n",
    "print(f\"Prediction data 1 horizon sample: {yt_pred[0].T}\")\n",
    "print(f\"Prediction 1 horizon sample: {yt_pred.flatten()[0]} VS {ytest_data.flatten()[0]}\")\n",
    "\n",
    "cm_train = confusion_matrix(ytrain_data.flatten(), y_pred.flatten())\n",
    "cm_test = confusion_matrix(ytest_data.flatten(), yt_pred.flatten())\n",
    "\n",
    "cm_train_np = cm_train.numpy()\n",
    "cm_test_np = cm_test.numpy()\n",
    "\n",
    "plot_confusion_matrix(cm_train_np, [0,1], cm_test_np, [0,1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "130e7627",
   "metadata": {},
   "source": [
    "We also look at other metrics, most already mentioned (AUC, TP, TN, etc.). Another metric we'll look at is the F1 Score, which is a statistical measure used to evaluate the performance of binary classification models in situations where the class distribution is imbalanced. \n",
    "It is the harmonic mean of precision and recall which signal the robustness of our classifier. The F1 Score reaches its best value at 1 (perfect precision and recall) and worst at 0, and is calculated as follows:\n",
    "$$\n",
    "F1 = 2 \\cdot \\frac{\\text{precision} \\cdot \\text{recall}}{\\text{precision} + \\text{recall}}\n",
    "$$\n",
    "\n",
    "An out of the box F1 score assumes both precision and recall are important, another variant called the F-Beta Score, defined as follows:\n",
    "\n",
    "$$\n",
    "F_{\\beta} = (1 + \\beta^2) \\cdot \\frac{\\text{precision} \\cdot \\text{recall}}{(\\beta^2 \\cdot \\text{precision}) + \\text{recall}}\n",
    "$$\n",
    "\n",
    "Where beta can be altered to give more importance to either of the metrics, weighting precision (<1) or recall (β>1).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cf9703b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import (\n",
    "    accuracy_score, precision_score, recall_score,\n",
    "    fbeta_score, roc_auc_score, roc_curve, auc\n",
    ")\n",
    "\n",
    "def directional_accuracy(y, y_pred):\n",
    "    a = np.array(y).flatten()\n",
    "    p = np.array(y_pred).flatten()\n",
    "\n",
    "    a_dir = np.sign(np.diff(a))\n",
    "    p_dir = np.sign(np.diff(p))\n",
    "    correct_dirs = np.sum(a_dir == p_dir)\n",
    "    acc = correct_dirs / len(a_dir)\n",
    "\n",
    "    return acc\n",
    "\n",
    "\n",
    "print(f\"shapes y_pred: {y_pred.shape} and yt_pred: {yt_pred.shape}\")\n",
    "\n",
    "accuracy_test = accuracy_score(ytest_data.flatten(), yt_pred.flatten())\n",
    "precision_test = precision_score(ytest_data.flatten(), yt_pred.flatten())\n",
    "recall_test = recall_score(ytest_data.flatten(), yt_pred.flatten())\n",
    "f1_test = fbeta_score(ytest_data.flatten(), yt_pred.flatten(), average=\"weighted\", beta=0.5)\n",
    "roc_auc_test = roc_auc_score(ytest_data.flatten(), yt_pred_orig.flatten(), average='weighted')\n",
    "metrics_df = pd.DataFrame({\n",
    "    \"Accuracy\": [accuracy_test],\n",
    "    \"Precision\": [precision_test],\n",
    "    \"Recall\": [recall_test],\n",
    "    \"F1b Score\": [f1_test],\n",
    "    \"ROC AUC\": [roc_auc_test],\n",
    "}, index=[\"Test\"])\n",
    "\n",
    "fig, axs = plt.subplots(3 if IS_CLASSIFICATION else 2, 1, figsize=(12, 10))\n",
    "axs[0].plot(history.history[\"loss\"], label=\"Train loss\")\n",
    "axs[0].plot(history.history[\"val_loss\"], label=\"Validation loss\")\n",
    "axs[0].set_xlabel(\"Epochs\")\n",
    "axs[0].set_ylabel(\"Loss\")\n",
    "axs[0].legend()\n",
    "axs[1].plot(\n",
    "    history.history[TARGET_METRIC],\n",
    "    label=f\"Train {TARGET_METRIC}\",\n",
    ")\n",
    "axs[1].plot(\n",
    "    history.history[f\"val_{TARGET_METRIC}\"],\n",
    "    label=f\"Test {TARGET_METRIC}\",\n",
    "    color=\"k\",\n",
    ")\n",
    "axs[1].axhline(accuracy_test, color=\"r\", linestyle=\"-.\", label=f\"Acc: {accuracy_test:.2f}\")\n",
    "axs[1].axhline(precision_test, color=\"g\", linestyle=\"--\", label=f\"Prec: {precision_test:.2f}\")\n",
    "axs[1].axhline(recall_test, color=\"g\", linestyle=\"-.\", label=f\"Rec: {recall_test:.2f}\")\n",
    "axs[1].axhline(f1_test, color=\"c\", linestyle=\"--\", label=f\"Fbeta: {f1_test:.2f}\")\n",
    "axs[1].set_xlabel(\"Epochs\")\n",
    "axs[1].set_ylabel(\"Score\")\n",
    "axs[1].legend()\n",
    "\n",
    "\n",
    "fpr, tpr, _ = roc_curve(ytest_data.flatten(), yt_pred_orig.flatten())\n",
    "roc_auc_test = auc(fpr, tpr)\n",
    "axs[2].plot(fpr, tpr, color='darkorange', lw=2, label=f'ROC curve (area = {roc_auc_test:.2f})')\n",
    "axs[2].plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--', label='ROC Line')\n",
    "axs[2].set_xlabel('False Positive Rate')\n",
    "axs[2].set_ylabel('True Positive Rate')\n",
    "axs[2].legend()\n",
    "\n",
    "plt.show()\n",
    "\n",
    "metrics_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b4b4c45",
   "metadata": {},
   "source": [
    "# Walk Forwards\n",
    "\n",
    "Very good metrics, with an AUC of almost 1 - which is a must in this domain of prediction if we wish to avoid financial losses. let's test the model on data it never saw with **walkforwards on out-of-sample data** which we have withheld before the whole training and validation cycle: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bf8ad22",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_oos_pred_raw = model.predict([xtrain_oos, Xexog_oos])\n",
    "y_oos_pred = (y_oos_pred_raw > 0.5).astype(int)\n",
    "\n",
    "print(f\"Prediction shape: {y_oos_pred.shape} vs test data shape: {y_oos.shape}\")\n",
    "print(f\"Test data 1 horizon sample: {y_oos[0]}\")\n",
    "print(f\"Predicted data 1 horizon sample: {y_oos_pred[0].T}\")\n",
    "print(f\"Prediction 1 horizon sample: {y_oos_pred.flatten()[0]} VS {y_oos.flatten()[0]}\")\n",
    "\n",
    "metrics_oos_df = None\n",
    "if IS_CLASSIFICATION:\n",
    "    accuracy_test = accuracy_score(y_oos.flatten(), y_oos_pred.flatten())\n",
    "    precision_test = precision_score(y_oos.flatten(), y_oos_pred.flatten())\n",
    "    recall_test = recall_score(y_oos.flatten(), y_oos_pred.flatten())\n",
    "    f1_test = fbeta_score(y_oos.flatten(), y_oos_pred.flatten(), average=\"weighted\", beta=0.5)\n",
    "    roc_auc_test = roc_auc_score(y_oos.flatten(), y_oos_pred_raw.flatten())\n",
    "    metrics_oos_df = pd.DataFrame({\n",
    "        \"Accuracy\": [accuracy_test],\n",
    "        \"Precision\": [precision_test],\n",
    "        \"Recall\": [recall_test],\n",
    "        \"F1b Score\": [f1_test],\n",
    "        \"ROC AUC\": [roc_auc_test],\n",
    "    }, index=[\"Test\"])\n",
    "\n",
    "metrics_oos_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8043ed7",
   "metadata": {},
   "source": [
    "## Cross Validation\n",
    "\n",
    "Good training/valdation metrics and good out-of-sample metrics, the model is too good to be trusted!\n",
    "\n",
    "We'll need more rigourous testing, starting with a **3 fold cross validation using a TimeSeries splitte**r from Scikit-Learn To increase trust in the output, we ill **perturb our data with faint gaussian noise**, this should destabilize a model without a robust classifier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e08d321",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import TimeSeriesSplit\n",
    "import json\n",
    "import os\n",
    "\n",
    "CV_MODEL = True\n",
    "CV_SPLITS = 3\n",
    "\n",
    "def train_cv_model(X, y, input_shape, conditioning_shapes=None, n_splits=5, perturb=True):\n",
    "    def _save_cv(results_df, file_name=\"cv_results.json\"):\n",
    "        os.makedirs(MODEL_DIR, exist_ok=True)\n",
    "        file_path = os.path.join(MODEL_DIR, file_name)\n",
    "        with open(file_path, \"w\") as file:\n",
    "            json.dump({\"CV results\": results_df.to_dict(orient=\"records\")}, file)\n",
    "\n",
    "    def _perturb_gaussiannoise(X, noise_level=0.1):\n",
    "        sigma = noise_level * np.std(X)\n",
    "        noise = np.random.normal(0, sigma, X.shape)\n",
    "        return X + noise\n",
    "\n",
    "    if perturb:\n",
    "        X = _perturb_gaussiannoise(X)\n",
    "\n",
    "    results = []\n",
    "    tscv = TimeSeriesSplit(n_splits=n_splits)\n",
    "\n",
    "    for train_index, test_index in tqdm(tscv.split(X), desc=f\"CV Testing for n_splits: {n_splits}\"):\n",
    "        X_train = X.iloc[train_index]\n",
    "        y_train = y.iloc[train_index]\n",
    "        X_test = X.iloc[test_index]\n",
    "        y_test = y.iloc[test_index]\n",
    "\n",
    "        X_train_windows, X_Exog_train, y_train_windows = prepare_windows(X_train, y_train)\n",
    "        X_test_windows, X_Exog_test, y_test_windows = prepare_windows(X_test, y_test)\n",
    "\n",
    "        try:\n",
    "            cv_model, _ = build_tcn(input_shape, [X_train_windows, X_Exog_train], y_train_windows, [X_test_windows, X_Exog_test], y_test_windows, conditioning_shapes=conditioning_shapes, tb=False)\n",
    "            result = cv_model.evaluate([X_test_windows, X_Exog_test], y_test_windows, verbose=0)\n",
    "            results.append(result)\n",
    "        except Exception as e:\n",
    "            print(f\"CV error on fold with exception: {e}\")\n",
    "\n",
    "    metrics_names = [metric.name for metric in cv_model.metrics]\n",
    "    results_df = pd.DataFrame(results, columns=metrics_names)\n",
    "    _save_cv(results_df)\n",
    "\n",
    "    return results_df\n",
    "\n",
    "\n",
    "if CV_MODEL:\n",
    "    results_df = train_cv_model(data_df, data_df[TARGET_LABEL], input_shape, conditioning_shapes=conditioning_shapes)\n",
    "    print(results_df)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77790620",
   "metadata": {},
   "source": [
    "# Predict Today"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66ac43c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "from pandas.tseries.offsets import BDay\n",
    "\n",
    "assert model is not None\n",
    "\n",
    "LAST_TRADING_DAY = pd.Timestamp(datetime.now() - BDay(1))\n",
    "FIRST_WINDOW_DAY = pd.Timestamp(LAST_TRADING_DAY - BDay(WINDOW_SIZE + 5))\n",
    "assert len( pd.bdate_range(start=FIRST_WINDOW_DAY, end=LAST_TRADING_DAY)) >= WINDOW_SIZE, f\"Expected larger than {WINDOW_SIZE} business days\"\n",
    "print(f\"Date ranges {FIRST_WINDOW_DAY} - {LAST_TRADING_DAY}\")\n",
    "\n",
    "pred_tickers = tickers.copy()\n",
    "for ticker in tickers_symbols:\n",
    "    df = tickers.get(ticker)\n",
    "    df = df[(df.index >= latest_start) & (df.index <= earliest_end)]\n",
    "    assert len(df) > 0 and not df.isna().any().any()\n",
    "    pred_tickers[ticker] = df\n",
    "assert len(pred_tickers) > 0 and len(pred_tickers[TARGET_ETF]) > 0\n",
    "\n",
    "preddata_df, _ = prepare_data(pred_tickers, to_normalize=True)\n",
    "print(f\"Pred Data Shape {preddata_df.shape}\")\n",
    "\n",
    "assert not np.any(pd.isna(preddata_df)) and not np.any(pd.isna(preddata_df))\n",
    "pred_X, pred_Xexog, _ = prepare_windows(preddata_df, None)\n",
    "\n",
    "today_pred = model.predict([pred_X, pred_Xexog])\n",
    "print(f\"today_pred Data Shape {today_pred.shape}\")\n",
    "# Remember shape: (Window, Days, features) - in label's case only 1 feature.\n",
    "print(f\"Predicting the {WINDOW_SIZE}th which was yesterday: {today_pred[-1][-2]}\")\n",
    "print(f\"Predicting the {WINDOW_SIZE+1}th which is today: {today_pred[-1][-1]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8950717",
   "metadata": {},
   "source": [
    "# Unseen Securities Performance\n",
    "\n",
    "The CV metrics are very promising, it does miss a lot of predictions, but those that it gets, is almost 100% correct! Still not to be trusted until we test the model in the real world now.\n",
    "\n",
    "Testing on the target index constituents, we assume the model has generalized well enough to have good performance, and we will not be suprised if it has excellent performance because all financial data is serially corellated and multtcollinear."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5f44935",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert model is not None\n",
    "\n",
    "OHTER_ETFS = [\"SSO\", \"SPXL\", \"BAC\", \"TQQQ\", \"AAPL\", \"HSBC\", \"TM\"]  # These are leveraged & constituents\n",
    "LAST_TRADING_DAY_STR = LAST_TRADING_DAY.strftime('%Y-%m-%d')\n",
    "FIRST_WINDOW_DAY_STR = FIRST_WINDOW_DAY.strftime('%Y-%m-%d')\n",
    "\n",
    "def generalization_test(unseen_tickers_list, start, end, gen_model=None, input_shape=input_shape, conditioning_shapes=conditioning_shapes):\n",
    "    metrics_unseen = []\n",
    "    for ticker in tqdm(unseen_tickers_list, desc=\"Testing on Unseen TS\"):\n",
    "        current_tickers = tickers.copy()\n",
    "        other_etf, latest_start, earliest_end = get_tickerdata([ticker], start=start, end=end)\n",
    "        current_tickers[TARGET_ETF] = other_etf[ticker]\n",
    "        print(f\"Date ranges for {ticker} are {latest_start.strftime('%Y-%m-%d')} - {earliest_end.strftime('%Y-%m-%d')}\")\n",
    "\n",
    "        for ticker_symbol in tickers_symbols:\n",
    "            df = current_tickers.get(ticker_symbol)\n",
    "            df = df[(df.index >= latest_start) & (df.index <= earliest_end)]\n",
    "            assert not df.empty and not df.isna().any().any() and len(df) > WINDOW_SIZE//2, f\"Data validation failed for {ticker_symbol}\"\n",
    "            current_tickers[ticker_symbol] = df\n",
    "        unseendata_df, _ = prepare_data(current_tickers, to_normalize=True)\n",
    "        assert not np.any(pd.isna(unseendata_df)) and len(unseendata_df) > WINDOW_SIZE//2, \"Unseen data preparation failed\"\n",
    "        unseen_X, unseen_Xexog, unseen_y = prepare_windows(unseendata_df, unseendata_df[TARGET_LABEL])\n",
    "\n",
    "        unseen_ypred_raw = gen_model.predict([unseen_X, unseen_Xexog])\n",
    "        unseen_pred = (unseen_ypred_raw > 0.5).astype(int)\n",
    "        metrics = {}\n",
    "        if IS_CLASSIFICATION:\n",
    "            metrics = {\n",
    "                \"Accuracy\": accuracy_score(unseen_y.flatten(), unseen_pred.flatten()),\n",
    "                \"Precision\": precision_score(unseen_y.flatten(), unseen_pred.flatten()),\n",
    "                \"Recall\": recall_score(unseen_y.flatten(), unseen_pred.flatten()),\n",
    "                \"F1b Score\": fbeta_score(y_oos.flatten(), y_oos_pred.flatten(), average=\"weighted\", beta=0.5),\n",
    "                \"ROC AUC\": roc_auc_score(unseen_y.flatten(), unseen_pred.flatten(), average='weighted')\n",
    "            }\n",
    "        metrics_unseen.append({ticker: metrics})\n",
    "    metrics_unseen_df = pd.DataFrame.from_dict({list(d.keys())[0]: list(d.values())[0] for d in metrics_unseen}, orient='index')\n",
    "    return metrics_unseen_df\n",
    "\n",
    "metrics_unseen_df = generalization_test(OHTER_ETFS, FIRST_WINDOW_DAY_STR, LAST_TRADING_DAY_STR, gen_model=model)\n",
    "metrics_unseen_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "940a6579",
   "metadata": {},
   "source": [
    "# Black Swan Performance\n",
    "\n",
    "Let's also do some exotic tests, we will train and apply our architecture to a black swan event, as we assume that an economy that is always on the rise is easy to predict, but in market turmoil, there should be some failure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f61177c",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert model is not None\n",
    "\n",
    "COVID_REGIME_START = \"2019-10-01\"\n",
    "COVID_REGIME_END = \"2021-01-01\"\n",
    "\n",
    "current_tickers = tickers.copy()\n",
    "other_etf, latest_start, earliest_end = get_tickerdata([ticker], start=START_DATE, end=COVID_REGIME_END)\n",
    "current_tickers[TARGET_ETF] = other_etf[ticker]\n",
    "print(f\"Date ranges for {ticker} are {latest_start.strftime('%Y-%m-%d')} - {earliest_end.strftime('%Y-%m-%d')}\")\n",
    "\n",
    "for ticker_symbol in tickers_symbols:\n",
    "    df = current_tickers.get(ticker_symbol)\n",
    "    df = df[(df.index >= latest_start) & (df.index <= earliest_end)]\n",
    "    assert not df.empty and not df.isna().any().any() and len(df) > WINDOW_SIZE//2, f\"Data validation failed for {ticker_symbol}\"\n",
    "    current_tickers[ticker_symbol] = df\n",
    "\n",
    "unseendata_df, _ = prepare_data(current_tickers, to_normalize=True)\n",
    "unseen_X, unseen_Xexog, unseen_y = prepare_windows(unseendata_df, unseendata_df[TARGET_LABEL])\n",
    "gen_model, gen_hist = build_tcn(input_shape, [unseen_X, unseen_Xexog], unseen_y, conditioning_shapes=conditioning_shapes, tb=False)\n",
    "print(f'{TARGET_ETF} val AUC: {gen_hist.history[f\"val_{TARGET_METRIC}\"][-1]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2565329b",
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics_unseen_df = generalization_test(OHTER_ETFS, START_DATE, COVID_REGIME_END, gen_model=gen_model)\n",
    "metrics_unseen_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7dde300",
   "metadata": {},
   "source": [
    "## Low Correlation Datasets.\n",
    "\n",
    "Performance is unaffected in a market downturn. High accuracy, high precision - what more do we want!\n",
    "\n",
    "Let's truly stress the model with predictions in the same chaotic conditions, on market-uncorrelated timeseries - specifically the meme-stocks that came to existance in the covid regime:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99943f7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "OHTER_ETFS = [\"ZM\", \"SHOP\", \"VIG\", \"VOO\", \"TSLA\", \"WDI.HM\"]  # Including meme stocks and frauds\n",
    "\n",
    "metrics_unseen_df = generalization_test(OHTER_ETFS, COVID_REGIME_START, COVID_REGIME_END, gen_model=gen_model)\n",
    "metrics_unseen_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78e9a9c5",
   "metadata": {
    "papermill": {
     "duration": 0.037345,
     "end_time": "2023-11-05T18:01:12.947046",
     "exception": false,
     "start_time": "2023-11-05T18:01:12.909701",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The model really seems to be unchallenged! Is the problem at hand that easy to solve for a moderately complex NN? Or are we not seeing some deep issue with our data and architecture?\n",
    "\n",
    "# Conclusion\n",
    "\n",
    "In this article, we tested a day-trading assumption that buying at start and selling at close is profitable - which we asserted that it has 29% chance of success. We outperformed this naive strategy using deep learning, which had a true positive rate of 99%.\n",
    "\n",
    "The deap learning model is one inspired by Google's WaveNet from Van den Oord et al., which uses a temporal convolution network, with causal dilating layers, skip layers to preserve historic residuals & expedite gradient intercept, and finally conditioning the input timeseries with multiple timeseries before the TCN's main convolutions are done.\n",
    "\n",
    "This is model, dataset and arcticle are for research purposes only, but here is GPT4's depiction of it beating the finfluencing youtuber, have fun!\n",
    "\n",
    "![image](./images/nn_vs_youtuber.PNG)\n",
    "\n",
    "## References\n",
    "\n",
    "- [YFinance Github](https://github.com/ranaroussi/yfinance)\n",
    "- [WaveNet: A Generative Model for Raw Audio](https://paperswithcode.com/paper/wavenet-a-generative-model-for-raw-audio)\n",
    "- [Time series forecasting Tensorflow](https://www.tensorflow.org/tutorials/structured_data/time_series)\n",
    "- [TensorBoard: TensorFlow's visualization toolkit](https://www.tensorflow.org/tensorboard)\n",
    "- [Principal component analysis](https://en.wikipedia.org/wiki/Principal_component_analysis)\n",
    "- [Spearman's rank correlation coefficient](https://en.wikipedia.org/wiki/Spearman%27s_rank_correlation_coefficient)\n",
    "- [Deep Learning for Time Series Forecasting: The Electric Load Case](https://arxiv.org/abs/1907.09207)\n",
    "- [Hyperparameter Tuning with the HParams Dashboard ](https://www.tensorflow.org/tensorboard/hyperparameter_tuning_with_hparams#1_experiment_setup_and_the_hparams_experiment_summary)\n",
    "- [fbeta_score](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.fbeta_score.html)\n",
    "- [Advances in Financial Machine Learning](https://www.amazon.com/Advances-Financial-Machine-Learning-Marcos/dp/1119482089)\n",
    "\n",
    "\n",
    "## Github\n",
    "\n",
    "Article here is also available on [Github](https://github.com/adamd1985/quant_research/blob/main/tcn_signals_research.ipynb)\n",
    "\n",
    "Kaggle notebook available [here](https://www.kaggle.com/addarm/temporal-convolution-neural-network-for-markets)\n",
    "\n",
    "\n",
    "## Media\n",
    "\n",
    "All media used (in the form of code or images) are either solely owned by me, acquired through licensing, or part of the Public Domain and granted use through Creative Commons License.\n",
    "\n",
    "## CC Licensing and Use\n",
    "\n",
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc/4.0/\"><img alt=\"Creative Commons License\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by-nc/4.0/88x31.png\" /></a><br />This work is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc/4.0/\">Creative Commons Attribution-NonCommercial 4.0 International License</a>."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  },
  "papermill": {
   "default_parameters": {},
   "duration": 98.905855,
   "end_time": "2023-11-05T18:01:13.607303",
   "environment_variables": {},
   "exception": null,
   "input_path": "__notebook__.ipynb",
   "output_path": "__notebook__.ipynb",
   "parameters": {},
   "start_time": "2023-11-05T17:59:34.701448",
   "version": "2.4.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}