alt_data_factor_research.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "\"\"\"\n",
    "Further suppose that the dataset is ticker-tagged for 200 equities which respond to the dataset in roughly similar ways. Maybe the economic basis of the data is shared by all the equities, or maybe investor trading behavior on the data unites them.\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# simulate the returns and signal\n",
    "r = 0 # avg return\n",
    "r_sig = 1 # avg return conditioned on signal\n",
    "sigma_alpha = 0.5 # alpha std dev\n",
    "sigma_beta = 0.3 # beta std dev\n",
    "rho = 0.4 # correlation between alpha and beta\n",
    "\n",
    "# Means\n",
    "Mu = np.array([r, r_sig])\n",
    "\n",
    "# covariance\n",
    "Cov = np.diag([r, r_sig]) @ np.array([1,rho,rho,1]).reshape(2,2) @ np.diag([r, r_sig])\n",
    "\n",
    "print(Mu)\n",
    "print(Cov)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "\"\"\"\n",
    "More formally we can say the alt dataset comprises an explanatory factor portfolio, and the equity returns have betas to that factor. That implies a covariance structure between the equities we can exploit for superior signal construction. So we continue simulating this data...\n",
    "\"\"\"\n",
    "\n",
    "# generate population equity data\n",
    "m_symbols = 200\n",
    "\n",
    "np.random.seed(60)\n",
    "mvr = np.random.multivariate_normal(Mu, Cov, m_symbols)\n",
    "\n",
    "alpha_ret = mvr[:,0]\n",
    "beta_ret = mvr[:,1]\n",
    "\n",
    "print(mvr)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "# simulate 50 return samples for each symbol - train set\n",
    "\n",
    "np.random.seed(64)\n",
    "n_samples = 50\n",
    "\n",
    "event = np.random.uniform(0,1, int(n_samples * m_symbols))\n",
    "sym_id = np.repeat(list(range(1, m_symbols + 1)), n_samples)\n",
    "mu = alpha_ret[sym_id - 1] + beta_ret[sym_id - 1] * event\n",
    "sigma = 0.5\n",
    "\n",
    "ret = np.random.normal(mu, sigma, m_symbols * n_samples)\n",
    "train_df = pd.DataFrame({'symbol': sym_id, 'event':event, 'ret': ret})\n",
    "train_df"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "\"\"\"\n",
    "Now we have a train set and a test set. First let's try fitting an OLS model to each equity based on only its sample of the data, individually. We'll estimate the model on the train set, then evaluate it via MSE on the test set (MSE chosen over MAE to penalize very high errors).\n",
    "\"\"\"\n",
    "\n",
    "# simulate 20 return samples for each symbol - test set\n",
    "\n",
    "np.random.seed(32)\n",
    "n_samples = 20\n",
    "\n",
    "event = np.random.uniform(0,1, int(n_samples * m_symbols))\n",
    "sym_id = np.repeat(list(range(1, m_symbols + 1)), n_samples)\n",
    "mu = alpha_ret[sym_id - 1] + beta_ret[sym_id - 1] * event\n",
    "sigma = 0.5\n",
    "\n",
    "ret = np.random.normal(mu, sigma, m_symbols * n_samples)\n",
    "test_df = pd.DataFrame({'symbol': sym_id, 'event':event, 'ret': ret})\n",
    "test_df"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "\"\"\"\n",
    "We evaluate the per-symbol models (\"no pooling\") and average their MSEs, then evaluate out of sample on the entire dataset (\"complete pooling\"). Per-symbol the average MSE is ~0.604, total is 1.21. Can we do better than this? Let's try estimating the model with complete pooling.\n",
    "\"\"\"\n",
    "\n",
    "import statsmodels.api as sm\n",
    "\n",
    "symbol_mse = []\n",
    "total_y = []\n",
    "total_y_hat = []\n",
    "rsq = []\n",
    "pval = []\n",
    "\n",
    "for symbol in train_df['symbol'].drop_duplicates().to_list():\n",
    "    train_data = train_df.loc[train_df['symbol'] == symbol][['event', 'ret']].to_numpy()\n",
    "    y = train_data.T[1].reshape(-1,1)\n",
    "    X = sm.add_constant(train_data.T[0])\n",
    "    model = sm.OLS(y, X).fit()\n",
    "    rsq.append(model.rsquared)\n",
    "    pval.append(model.pvalues)\n",
    "\n",
    "    # evaluate the model via MSE\n",
    "    test_data = test_df.loc[test_df['symbol'] == symbol][['event', 'ret']].to_numpy()\n",
    "    y = test_data.T[1].reshape(-1,1)\n",
    "    X = sm.add_constant(test_data.T[0])\n",
    "\n",
    "    y_hat = model.predict(exog=X)\n",
    "\n",
    "    mse = np.square(y - y_hat).mean()\n",
    "    total_y += [k for k in y]\n",
    "    total_y_hat += [k for k in y_hat]\n",
    "    symbol_mse.append(mse)\n",
    "\n",
    "# measure avg of MSE across per-symbol models\n",
    "np.mean(symbol_mse)\n"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "# also measure MSE across entire dataset - no pooling performs much worse on entire dataset than per-symbol\n",
    "\n",
    "np.square(np.array(total_y) - np.array(total_y_hat)).mean()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "\"\"\"\n",
    "The model estimates the alt dataset's signal across the entire sample of equities without any conditioning on equities. MSEs on a per-symbol basis are worse than before, but the total dataset MSE has improved. Can we do better than this?\n",
    "\"\"\"\n",
    "\n",
    "pooled_mse = []\n",
    "\n",
    "train_data = train_df[['event', 'ret']].to_numpy()\n",
    "y = train_data.T[1].reshape(-1,1)\n",
    "X = sm.add_constant(train_data.T[0])\n",
    "model = sm.OLS(y, X).fit()\n",
    "print(model.summary())"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "for symbol in train_df['symbol'].drop_duplicates().to_list():\n",
    "    # evaluate the model via MSE\n",
    "\n",
    "    test_data = test_df.loc[test_df['symbol'] == symbol][['event', 'ret']].to_numpy()\n",
    "    y = test_data.T[1].reshape(-1,1)\n",
    "    X = sm.add_constant(test_data.T[0])\n",
    "\n",
    "    #y_hat = X @ model.params.T\n",
    "    y_hat = model.predict(exog=X)\n",
    "\n",
    "    mse = np.square(y - y_hat).mean()\n",
    "    pooled_mse.append(mse)\n",
    "\n",
    "np.mean(pooled_mse)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "# also measure MSE across entire dataset - total pooling performs much better on entire dataset than per-symbol\n",
    "\n",
    "X = sm.add_constant(test_df['event'].to_numpy())\n",
    "y_hat = model.predict(X)\n",
    "np.square(y-y_hat).mean()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "\"\"\"\n",
    "In fact, we can fit a multilevel regression model instead. The \"fixed effects\" have an interpretation as the completely pooled regression coefficient, and the \"random effects\" are per-symbol deviations. To make a per-symbol prediction, we sum the effects. The MSEs improve.\n",
    "\"\"\"\n",
    "\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "# build mixed effects model - partial pooling on symbols\n",
    "\n",
    "model = smf.mixedlm('ret ~ event', train_df, groups=train_df['symbol'])\n",
    "md = model.fit()\n",
    "re = {symbol: md.random_effects[symbol][0] for symbol in md.random_effects}\n",
    "test_df['fe'] = md.predict(exog=test_df)\n",
    "test_df['re'] = test_df['symbol'].map(re)\n",
    "test_df['ret_hat'] = test_df['fe'] + test_df['re']\n",
    "\n",
    "# evaluate the model per-symbol - it outperforms both the unpooled and completely pooled versions\n",
    "\n",
    "mse = []\n",
    "\n",
    "for symbol in test_df['symbol'].drop_duplicates().to_list():\n",
    "    test_data = test_df.loc[test_df['symbol'] == symbol]\n",
    "    sym_mse = np.square(test_data['ret'] - test_data['ret_hat']).mean()\n",
    "    mse.append(sym_mse)\n",
    "\n",
    "np.mean(mse)"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "# evaluate the MSE across the entire dataset - it outperforms both the unpooled and completely pooled versions\n",
    "\n",
    "me_mse = np.square(test_df['ret'] - test_df['ret_hat']).mean()\n",
    "me_mse"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "\"\"\"\n",
    "Aside from improved out-of-sample accuracy, this comes with an elegant basis for portfolio construction. The random effects have a natural interpretation as the relative out (under) performance of each symbol, neutral to the group/cluster on which the model was estimated.\n",
    "\n",
    "So for each alt dataset with which this modeling approach works, you can simple sort by random effects in each group, long the top decile of random effects, short the bottom decile of random effects, and hedge the remainder to trade a neutral basket.\n",
    "\"\"\"\n",
    "\n",
    "import quantstats as qs\n",
    "\n",
    "# Sort within each symbol group by random effects\n",
    "train_df = train_df.sort_values(by=['event'])\n",
    "\n",
    "# Define function to classify deciles\n",
    "def classify_decile(group):\n",
    "    group['decile'] = pd.qcut(group['event'], 10, labels=False)\n",
    "    return group\n",
    "\n",
    "# Apply the function to the DataFrame\n",
    "train_df = train_df.groupby('symbol', group_keys=False).apply(classify_decile)\n",
    "\n",
    "# Create long and short positions\n",
    "train_df['position'] = 0\n",
    "train_df.loc[train_df['decile'] == 9, 'position'] = 1  # Long top decile\n",
    "train_df.loc[train_df['decile'] == 0, 'position'] = -1  # Short bottom decile\n",
    "\n",
    "# Hedge the remainder to trade a neutral basket\n",
    "train_df['hedged_position'] = train_df.groupby('symbol')['position'].transform(lambda x: x - x.mean())\n",
    "\n",
    "# Calculate the strategy returns\n",
    "train_df['strategy_ret'] = train_df['hedged_position'] * train_df['ret']\n",
    "\n",
    "# Aggregate strategy returns\n",
    "strategy_returns = train_df.groupby('symbol')['strategy_ret'].sum()\n",
    "\n",
    "# Convert to cumulative returns for QuantStats analysis\n",
    "strategy_returns"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [
    "strategy_returns.plot()"
   ],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "source": [],
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "execution_count": null
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}