MQO Ceteris Paribus.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MQO Ceteris Paribus\n",
    "\n",
    "### Lendo os Dados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<table border=\"1\" class=\"dataframe\">\n",
      "  <thead>\n",
      "    <tr style=\"text-align: right;\">\n",
      "      <th></th>\n",
      "      <th>age</th>\n",
      "      <th>female</th>\n",
      "      <th>education</th>\n",
      "      <th>earnings</th>\n",
      "      <th>hours</th>\n",
      "      <th>week</th>\n",
      "      <th>tenure</th>\n",
      "    </tr>\n",
      "  </thead>\n",
      "  <tbody>\n",
      "    <tr>\n",
      "      <th>0</th>\n",
      "      <td>52</td>\n",
      "      <td>0</td>\n",
      "      <td>12</td>\n",
      "      <td>62.393162</td>\n",
      "      <td>45</td>\n",
      "      <td>52</td>\n",
      "      <td>34</td>\n",
      "    </tr>\n",
      "    <tr>\n",
      "      <th>1</th>\n",
      "      <td>38</td>\n",
      "      <td>0</td>\n",
      "      <td>18</td>\n",
      "      <td>21.367521</td>\n",
      "      <td>45</td>\n",
      "      <td>52</td>\n",
      "      <td>14</td>\n",
      "    </tr>\n",
      "    <tr>\n",
      "      <th>2</th>\n",
      "      <td>38</td>\n",
      "      <td>0</td>\n",
      "      <td>14</td>\n",
      "      <td>15.686275</td>\n",
      "      <td>40</td>\n",
      "      <td>51</td>\n",
      "      <td>18</td>\n",
      "    </tr>\n",
      "    <tr>\n",
      "      <th>3</th>\n",
      "      <td>41</td>\n",
      "      <td>1</td>\n",
      "      <td>13</td>\n",
      "      <td>22.596154</td>\n",
      "      <td>40</td>\n",
      "      <td>52</td>\n",
      "      <td>22</td>\n",
      "    </tr>\n",
      "    <tr>\n",
      "      <th>4</th>\n",
      "      <td>42</td>\n",
      "      <td>0</td>\n",
      "      <td>13</td>\n",
      "      <td>62.125000</td>\n",
      "      <td>50</td>\n",
      "      <td>52</td>\n",
      "      <td>23</td>\n",
      "    </tr>\n",
      "  </tbody>\n",
      "</table>\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from numpy import dot\n",
    "from numpy.linalg import inv\n",
    "\n",
    "variables = {0: \"age\",\n",
    "             1: \"female\",\n",
    "             3: \"education\",\n",
    "             4: \"earnings\",\n",
    "             5: \"hours\",\n",
    "             6: \"week\"}\n",
    "\n",
    "\n",
    "dataset = (pd.read_csv(\"./data/cps09mar.txt\", sep = \"\\t\", header=None)\n",
    "            .rename(index=str, columns=variables)\n",
    "           .assign(earnings=lambda df: df[\"earnings\"] / (df[\"hours\"] * df[\"week\"]),\n",
    "                   tenure=lambda df: df[\"age\"] - df[\"education\"] - 6)\n",
    "           .query(\"earnings>0\")\n",
    "           .loc[:, list(variables.values()) + [\"tenure\"]]\n",
    "          )\n",
    "\n",
    "print(dataset.head().to_html())\n",
    "# dataset.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OLS Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>    1.1501</td> <td>    0.016</td> <td>   71.436</td> <td> 0.000</td> <td>    1.119</td> <td>    1.182</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>education</th> <td>    0.1082</td> <td>    0.001</td> <td>  114.390</td> <td> 0.000</td> <td>    0.106</td> <td>    0.110</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age</th>       <td>    0.0095</td> <td>    0.000</td> <td>   42.276</td> <td> 0.000</td> <td>    0.009</td> <td>    0.010</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>female</th>    <td>   -0.2629</td> <td>    0.005</td> <td>  -50.165</td> <td> 0.000</td> <td>   -0.273</td> <td>   -0.253</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import statsmodels.formula.api as smf\n",
    "\n",
    "results = smf.ols('np.log(earnings) ~ education + age + female', data=dataset).fit()\n",
    "\n",
    "results.summary().tables[1]\n",
    "# print(results.summary().tables[1].as_html())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = (dataset\n",
    "     .loc[:, [\"education\", \"age\", \"female\"]]\n",
    "     .assign(intercept=1))\n",
    "\n",
    "y = (dataset\n",
    "     .assign(log_earnings= lambda df: np.log(df[\"earnings\"]))\n",
    "     .loc[:, [\"log_earnings\"]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>beta_hat</th>\n",
       "      <th>variables</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.108156</td>\n",
       "      <td>education</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.009540</td>\n",
       "      <td>age</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.262894</td>\n",
       "      <td>female</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.150116</td>\n",
       "      <td>intercept</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   beta_hat  variables\n",
       "0  0.108156  education\n",
       "1  0.009540        age\n",
       "2 -0.262894     female\n",
       "3  1.150116  intercept"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def OLS_estimator(X, y):\n",
    "    variables = X.columns\n",
    "    beta_hat = inv(X.T.dot(X)).dot(X.T).dot(y).reshape(-1,)\n",
    "    return pd.DataFrame({\"variables\": variables, \"beta_hat\":beta_hat})\n",
    "\n",
    "def OLS_predictor(X, y):\n",
    "    variables = X.columns\n",
    "    beta_hat = inv(X.T.dot(X)).dot(X.T).dot(y).reshape(-1,)\n",
    "    return X.values.dot(beta_hat)\n",
    "\n",
    "OLS_estimator(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Least Squares Residuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = (dataset\n",
    "     .loc[:, [\"education\", \"age\", \"female\"]]\n",
    "     .assign(intercept=1))\n",
    "\n",
    "y = (dataset\n",
    "     .assign(log_earnings= lambda df: np.log(df[\"earnings\"]))\n",
    "     .loc[:, [\"log_earnings\"]])\n",
    "\n",
    "beta_hat = OLS_estimator(X, y).beta_hat\n",
    "y_hat = X.values.dot(beta_hat)\n",
    "e_hat = y.values.reshape(-1,) - y_hat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.round(e_hat.mean(), 6) # mean zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0., -0., -0., -0.])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.round(X.values.T.dot(e_hat), 6) # no correlation with X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Projection Matirx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample = dataset.sample(1000)\n",
    "\n",
    "X = (sample\n",
    "     .loc[:, [\"education\", \"age\", \"female\"]]\n",
    "     .assign(intercept=1)\n",
    "     .values)\n",
    "\n",
    "y = (sample\n",
    "     .assign(log_earnings= lambda df: np.log(df[\"earnings\"]))\n",
    "     .loc[:, [\"log_earnings\"]]\n",
    "     .values)\n",
    "\n",
    "beta_hat = inv(X.T.dot(X)).dot(X.T).dot(y)\n",
    "\n",
    "\n",
    "P = X.dot(inv(X.T.dot(X))).dot(X.T) # this is huge. Watch out!!! (NxN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_hat1 = P.dot(y) \n",
    "y_hat2 = X.dot(beta_hat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.all(np.isclose(y_hat2, y_hat1, 1e-6)) # both are the same"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ortogonal Projection Matirx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = np.eye(P.shape[0]) - P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.],\n",
       "       ...,\n",
       "       [ 0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.],\n",
       "       [-0., -0., -0., -0.]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.round(M.dot(X), 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.361755153084232e-14, 1.3500311979441903e-14)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e_hat1 = M.dot(y)\n",
    "e_hat2 = y - y_hat1\n",
    "e_hat1.mean(), e_hat2.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.all(np.isclose(e_hat1, e_hat2, 1e-6)) # both are the same"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression on Components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def regress(y, X):\n",
    "    return inv(X.T.dot(X)).dot(X.T).dot(y)\n",
    "\n",
    "X = (dataset\n",
    "     .loc[:, [\"education\", \"age\", \"female\"]]\n",
    "     .assign(intercept=1))\n",
    "\n",
    "y = (dataset\n",
    "     .assign(log_earnings= lambda df: np.log(df[\"earnings\"]))\n",
    "     .loc[:, [\"log_earnings\"]]\n",
    "     .values)\n",
    "\n",
    "X1 = X.loc[:, [\"intercept\", \"age\", \"female\"]].values\n",
    "X2 = X.loc[:, [\"education\"]].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.10815617]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta1_tilde = regress(y, X1)\n",
    "y1_hat_tilde = X1.dot(beta1_tilde)\n",
    "e1_tilde = y - y1_hat_tilde\n",
    "\n",
    "beta2_tilde = regress(X2, X1)\n",
    "X2_hat_tilde = X1.dot(beta2_tilde)\n",
    "X2e_tilde = X2 - X2_hat_tilde\n",
    "\n",
    "beta2_hat = regress(e1_tilde, X2e_tilde)\n",
    "\n",
    "beta2_hat # coef for education"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.10815617],\n",
       "       [ 0.00954044],\n",
       "       [-0.26289372],\n",
       "       [ 1.15011621]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regress(y, X.values)"
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}