bootstrap_algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bootstrapped quantiles and standard errors for 90% confidence intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from sklearn.utils import resample\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import rpy2.robjects.packages as rpackages\n",
    "from rpy2.robjects.packages import importr\n",
    "import rpy2.rinterface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "KernSmooth = rpackages.importr('KernSmooth')\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import and prepare data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def process_data(df, output_file):\n",
    "    \"\"\"This function adds squared and interaction terms to the Carneiro data set.\"\"\"\n",
    "\n",
    "    # Delete redundant columns\\n\",\n",
    "    for key_ in [\"newid\", \"caseid\"]:\n",
    "        del df[key_]\n",
    "\n",
    "    # Add squared terms\n",
    "    for key_ in [\"mhgc\", \"cafqt\", \"avurate\", \"lurate_17\", \"numsibs\", \"lavlocwage17\"]:\n",
    "        str_ = key_ + \"sq\"\n",
    "        df[str_] = df[key_] ** 2\n",
    "\n",
    "    # Add interaction terms\n",
    "    for j in [\"pub4\", \"lwage5_17\", \"lurate_17\", \"tuit4c\"]:\n",
    "        for i in [\"cafqt\", \"mhgc\", \"numsibs\"]:\n",
    "            df[j + i] = df[j] * df[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "basic = pd.read_stata('data/basicvariables.dta')\n",
    "local = pd.read_stata('data/localvariables.dta') \n",
    "df = pd.concat([basic, local], axis = 1)\n",
    "process_data(df,'data/aer-replication-mock')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Z_list = ['const', \n",
    " 'cafqt',\n",
    " 'cafqtsq',\n",
    " 'mhgc',\n",
    " 'mhgcsq',\n",
    " 'numsibs',\n",
    " 'numsibssq',\n",
    " 'urban14',\n",
    " 'lavlocwage17',\n",
    " 'lavlocwage17sq',\n",
    " 'avurate',\n",
    " 'avuratesq',\n",
    " 'd57',\n",
    " 'd58',\n",
    " 'd59',\n",
    " 'd60',\n",
    " 'd61',\n",
    " 'd62',\n",
    " 'd63',\n",
    " 'lwage5_17numsibs',\n",
    " 'lwage5_17mhgc',\n",
    " 'lwage5_17cafqt',\n",
    " 'lwage5_17',\n",
    " 'lurate_17',\n",
    " 'lurate_17numsibs',\n",
    " 'lurate_17mhgc',\n",
    " 'lurate_17cafqt',\n",
    " 'tuit4c',\n",
    " 'tuit4cnumsibs',\n",
    " 'tuit4cmhgc',\n",
    " 'tuit4ccafqt',\n",
    " 'pub4',\n",
    " 'pub4numsibs',\n",
    " 'pub4mhgc',\n",
    " 'pub4cafqt']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "X_list = ['const',\n",
    "'exp', \n",
    "'expsq', \n",
    "'lwage5', \n",
    "'lurate', \n",
    "'cafqt', \n",
    "'cafqtsq', \n",
    "'mhgc', \n",
    "'mhgcsq', \n",
    "'numsibs', \n",
    "'numsibssq', \n",
    "'urban14', \n",
    "'lavlocwage17', \n",
    "'lavlocwage17sq', \n",
    "'avurate', \n",
    "'avuratesq', \n",
    "'d57', \n",
    "'d58', \n",
    "'d59', \n",
    "'d60', \n",
    "'d61', \n",
    "'d62', \n",
    "'d63']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bootstrap algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, draw $B = 250$ independent bootstrap samples and compute the MTE* $ \\ \\hat{\\theta}^{MTE}(b)$ for each bootstrap replication $b = 1,..,B$, where $B$ is the number of bootstrap samples used. Second, estimate the standard error of $\\hat{\\theta}^{MTE}$, denoted by $\\widehat{\\mathbf{se}}_B$, separately for each gridpoint $u_D$:\n",
    "\n",
    "$$ \\widehat{\\mathbf{se}}_B = \\frac{1}{B-1} \\sum_{b = 1}^{B}{\\{ \\hat{\\theta}}^{MTE}(b) \\ - \\ \\hat{\\theta}^{MTE}(\\cdot)    \\}^2 $$\n",
    "\n",
    "where $\\hat{\\theta}^{MTE}(\\cdot) = B^{-1} \\sum_{b = 1}^{B}{ \\hat{\\theta}}^{MTE}(b)$.\n",
    "\n",
    "With every iteration of the bootstrap, $\\hat{P}(z)$ is reestimated, so the evaluation points $u_D$ over which to compute the *MTE* vary with each iteration. I compute the mean for each row of the 500 evaluation points that are reproduced with each bootstrap. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Number of iterations\n",
    "B = 250"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Prepare empty arrays to store output values\n",
    "mte_4_R = np.zeros([500, B])\n",
    "quantiles_u_R = np.zeros([500, B])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for i in range(B): \n",
    "    boot = resample(df, replace=True, n_samples=len(df), random_state=None)\n",
    "\n",
    "    # define Z and estimate propensity score via logit\n",
    "    Z_boot = boot[Z_list]\n",
    "    logitRslt_boot = sm.Logit(boot['state'], Z_boot).fit(disp=0)\n",
    "    \n",
    "    ps = logitRslt_boot.predict(Z_boot)\n",
    "    \n",
    "    # got an error before, had to change \"boot['ps'] = ps\" to the code below\n",
    "    boot.insert(len(df.columns), 'ps', ps)\n",
    "    \n",
    "    treated = boot[['state', 'ps']][boot['state'] == 1]\n",
    "    untreated = boot[['state', 'ps']][boot['state'] == 0]\n",
    "    \n",
    "    boot_trim = boot[(boot.ps >= np.min(treated['ps'])) & (boot.ps <= np.max(untreated['ps']))]\n",
    "    ps_trim = ps[(ps >= np.min(treated['ps'])) & (ps <= np.max(untreated['ps']))]\n",
    "    \n",
    "    # sort by ps\n",
    "    boot_trim = boot_trim.sort_values(by='ps', ascending=True)\n",
    "    ps_trim = np.sort(ps_trim)\n",
    "    ps_trim = pd.Series(ps_trim)\n",
    "    \n",
    "    X_trim = boot_trim[X_list]\n",
    "    X_trim.reset_index(drop=True, inplace=True)\n",
    "    \n",
    "    # construct X_k * P(z)\n",
    "    N = len(X_trim)\n",
    "    P_z = pd.concat([ps_trim]*len(X_trim.columns), axis=1, ignore_index=True)\n",
    "    Xp_trim= pd.DataFrame(X_trim.values*P_z.values, columns=[key_ + \"_ps\" for key_ in list(X_trim)], index=X_trim.index)\n",
    "    \n",
    "    \n",
    "    Y_trim = boot_trim[['wage']]\n",
    "\n",
    "    %reload_ext rpy2.ipython\n",
    "        \n",
    "    %R -i X_trim\n",
    "    %R -i Xp_trim\n",
    "    %R -i Y_trim\n",
    "    %R -i ps_trim\n",
    "\n",
    "    %R X <- as.matrix(X_trim)\n",
    "    %R Xp <- as.matrix(Xp_trim)\n",
    "    %R Y <- as.matrix(Y_trim)\n",
    "    %R ps <- as.matrix(ps_trim)\n",
    "\n",
    "    # Compute the lowess residuals of X, Xp and Y    \n",
    "    %R res_Y_R <- loess(Y ~ ps, span = 0.05, degree = 1)$res\n",
    "    %R res_X_R <- apply(X, 2, function(x) loess(x ~ ps, span = 0.05, degree = 1)$res)\n",
    "    %R res_Xp_R <- apply(Xp, 2, function(x) loess(x ~ ps, span = 0.05, degree = 1)$res)\n",
    "\n",
    "    %R res_Y_R <- data.frame(res_Y_R)\n",
    "    %R res_X_R <- data.frame(res_X_R)\n",
    "    %R res_Xp_R <- data.frame(res_Xp_R)\n",
    "\n",
    "    %R -o res_Y_R\n",
    "    %R -o res_X_R\n",
    "    %R -o res_Xp_R\n",
    "\n",
    "    # Append res_X and res_Xp. \n",
    "    res_X_Xp_R = np.append(res_X_R, res_Xp_R, axis=1)\n",
    "\n",
    "    model_R = sm.OLS(res_Y_R, res_X_Xp_R)\n",
    "    results_R = model_R.fit()\n",
    "    b1_R = results_R.params[:len(list(X_trim))]\n",
    "    b1_b0_R = results_R.params[len((list(X_trim))):]\n",
    "\n",
    "    # prepare arrays\n",
    "    ps_arr_R = np.array(ps_trim)\n",
    "    X_arr_R = np.array(X_trim)\n",
    "    Xp_arr_R = np.array(Xp_trim)\n",
    "    Y_arr_R = boot_trim['wage']\n",
    "\n",
    "    # compute the unobserved part of Y\n",
    "    Y_tilde_R = Y_arr_R - np.dot(X_arr_R, b1_R) - np.dot(Xp_arr_R, b1_b0_R)\n",
    "        \n",
    "    %reload_ext rpy2.ipython\n",
    "    \n",
    "    %R -i Y_tilde_R\n",
    "    %R -i ps_arr_R\n",
    "    %R Y_tilde_R <- as.matrix(Y_tilde_R)\n",
    "    %R ps_arr_R <- as.matrix(ps_arr_R)\n",
    "\n",
    "    %R mte_u_poly_R <- locpoly(ps_arr_R, Y_tilde_R, drv = 1L, degree=2, bandwidth = 0.322, gridsize = 500L, range.x = c(0.005, 0.995))\n",
    "\n",
    "    # Return the mte_u component and the quantiles of u_D as data frames and export them back into Python.\n",
    "    %R mte_u_R <- mte_u_poly_R$y\n",
    "    %R mte_u_R <- data.frame(mte_u_R)\n",
    "    %R quantiles_R <- mte_u_poly_R$x\n",
    "    %R quantiles_R <- data.frame(quantiles_R)\n",
    "\n",
    "    %R -o mte_u_R\n",
    "    %R -o quantiles_R\n",
    "    \n",
    "    # Turn data frames into arrays\n",
    "    mte_u_R = np.array(mte_u_R)\n",
    "    mte_u_R = np.concatenate(mte_u_R, axis=0)\n",
    "    \n",
    "    quantiles_R = np.array(quantiles_R)\n",
    "    quantiles_R = np.concatenate(quantiles_R, axis=0)\n",
    "    \n",
    "    mte_x_R = np.dot(X_trim, b1_b0_R).mean(axis=0)\n",
    "    mte_R = mte_x_R + mte_u_R\n",
    "    \n",
    "    mte_4_R[:,i] = mte_R/4\n",
    "        \n",
    "    quantiles_u_R[:,i] = quantiles_R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate mean of quantiles $u_D$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "quantiles_R_mean = quantiles_u_R.mean(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Turn into data frame to export as csv\n",
    "quantiles_R_mean = pd.DataFrame(quantiles_R_mean)\n",
    "quantiles_R_mean.to_csv(r'/Users/sebastiangsell/Documents/MSc_Economics_Bonn/Thesis/replication-semipar-carneiro2011/quantiles_R_mean.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute separate standard error for each gridpoint $u_D$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mte_R_mean = mte_4_R.mean(axis=1)\n",
    "\n",
    "no_sum = np.zeros([len(mte_4_R), B])\n",
    "\n",
    "for row in range(len(mte_4_R)):\n",
    "    for column in range(B):\n",
    "        no_sum[row,column] = ((mte_4_R[row,column] - mte_R_mean[row])**2)\n",
    "\n",
    "se_R = (1/(B-1)) * np.sum(no_sum, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Turn into data frame to export as csv\n",
    "se_R = pd.DataFrame(se_R)\n",
    "se_R.to_csv(r'/Users/sebastiangsell/Documents/MSc_Economics_Bonn/Thesis/replication-semipar-carneiro2011/se_R.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}