EconometricsHW8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 (Remark: I used R instead of STATA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 403,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "library(mvtnorm)\n",
    "n =1000\n",
    "set.seed(10)\n",
    "mu = c(0, 0)\n",
    " A <- array(c(0.5, 0.25, 0.25, 4 ), dim=c(2,2))\n",
    "w = rmvnorm(n, mu, A)\n",
    "u = w[,1]\n",
    "v = w[,2]\n",
    "# plot(w)\n",
    "x = rnorm(n, mean = 0 ,sd = 1)\n",
    "z = rnorm(n, mean = 0 ,sd = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 404,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df    <- data.frame(u=u, v=v, x=x , z =z) \n",
    "# head(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 405,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "D = sign(1 + z + v)\n",
    "Ys = 1 + x + u\n",
    "df    <- data.frame(df ,dc =1 + z + v, D =D, Ys = Ys)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 406,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#replace the result of sign func.\n",
    "df$D[df$D == 0] = 1 \n",
    "df$D[df$D == -1] = 0 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 407,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>u</th><th scope=col>v</th><th scope=col>x</th><th scope=col>z</th><th scope=col>dc</th><th scope=col>D</th><th scope=col>Ys</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>1</th><td>-0.003926272</td><td>-0.3663731</td><td>-0.307865</td><td>0.4605151</td><td>1.094142</td><td>1</td><td>0.6882087</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>-1.016823</td><td>-1.324078</td><td>0.7580856</td><td>0.2350253</td><td>-0.08905249</td><td>0</td><td>0.7412631</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>0.2425872</td><td>0.8060361</td><td>-0.5738634</td><td>0.6432573</td><td>2.449293</td><td>1</td><td>0.6687238</td></tr>\n",
       "\t<tr><th scope=row>4</th><td>-0.8805651</td><td>-0.8384773</td><td>-0.9387445</td><td>0.9131981</td><td>1.074721</td><td>1</td><td>-0.8193096</td></tr>\n",
       "\t<tr><th scope=row>5</th><td>-1.164077</td><td>-0.6630874</td><td>-0.02769932</td><td>0.988286</td><td>1.325199</td><td>1</td><td>-0.1917762</td></tr>\n",
       "\t<tr><th scope=row>6</th><td>0.8423711</td><td>1.612</td><td>-1.066249</td><td>0.1127413</td><td>2.724742</td><td>1</td><td>0.7761224</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lllllll}\n",
       "  & u & v & x & z & dc & D & Ys\\\\\n",
       "\\hline\n",
       "\t1 & -0.003926272 & -0.3663731 & -0.307865 & 0.4605151 & 1.094142 & 1 & 0.6882087\\\\\n",
       "\t2 & -1.016823 & -1.324078 & 0.7580856 & 0.2350253 & -0.08905249 & 0 & 0.7412631\\\\\n",
       "\t3 & 0.2425872 & 0.8060361 & -0.5738634 & 0.6432573 & 2.449293 & 1 & 0.6687238\\\\\n",
       "\t4 & -0.8805651 & -0.8384773 & -0.9387445 & 0.9131981 & 1.074721 & 1 & -0.8193096\\\\\n",
       "\t5 & -1.164077 & -0.6630874 & -0.02769932 & 0.988286 & 1.325199 & 1 & -0.1917762\\\\\n",
       "\t6 & 0.8423711 & 1.612 & -1.066249 & 0.1127413 & 2.724742 & 1 & 0.7761224\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "             u          v           x         z          dc D         Ys\n",
       "1 -0.003926272 -0.3663731 -0.30786498 0.4605151  1.09414192 1  0.6882087\n",
       "2 -1.016822509 -1.3240778  0.75808558 0.2350253 -0.08905249 0  0.7412631\n",
       "3  0.242587220  0.8060361 -0.57386341 0.6432573  2.44929346 1  0.6687238\n",
       "4 -0.880565107 -0.8384773 -0.93874446 0.9131981  1.07472077 1 -0.8193096\n",
       "5 -1.164076881 -0.6630874 -0.02769932 0.9882860  1.32519861 1 -0.1917762\n",
       "6  0.842371051  1.6120003 -1.06624865 0.1127413  2.72474167 1  0.7761224"
      ]
     },
     "execution_count": 407,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "head(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (b) Heckit by hand  \n",
    " \n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Probit estimatin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 408,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "glm(formula = D ~ z, family = binomial(probit), data = df)\n",
       "\n",
       "Deviance Residuals: \n",
       "    Min       1Q   Median       3Q      Max  \n",
       "-2.6420  -1.1581   0.6316   0.8885   1.5408  \n",
       "\n",
       "Coefficients:\n",
       "            Estimate Std. Error z value Pr(>|z|)    \n",
       "(Intercept)  0.51812    0.04374   11.85   <2e-16 ***\n",
       "z            0.46218    0.04549   10.16   <2e-16 ***\n",
       "---\n",
       "Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n",
       "\n",
       "(Dispersion parameter for binomial family taken to be 1)\n",
       "\n",
       "    Null deviance: 1256.7  on 999  degrees of freedom\n",
       "Residual deviance: 1142.2  on 998  degrees of freedom\n",
       "AIC: 1146.2\n",
       "\n",
       "Number of Fisher Scoring iterations: 4\n"
      ]
     },
     "execution_count": 408,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/html": [
       "<strong>z:</strong> 0.462176132791821"
      ],
      "text/latex": [
       "\\textbf{z:} 0.462176132791821"
      ],
      "text/markdown": [
       "**z:** 0.462176132791821"
      ],
      "text/plain": [
       "        z \n",
       "0.4621761 "
      ]
     },
     "execution_count": 408,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res <- glm(D ~  z , data = df, family=binomial(probit))\n",
    "summary(res)\n",
    "    res$coefficients[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Inverse Mill's Ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 409,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>0.667508521484319</li>\n",
       "\t<li>0.730043593087919</li>\n",
       "\t<li>0.618751345207167</li>\n",
       "\t<li>0.550024750476546</li>\n",
       "\t<li>0.531629181306629</li>\n",
       "\t<li>0.765011290422085</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 0.667508521484319\n",
       "\\item 0.730043593087919\n",
       "\\item 0.618751345207167\n",
       "\\item 0.550024750476546\n",
       "\\item 0.531629181306629\n",
       "\\item 0.765011290422085\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 0.667508521484319\n",
       "2. 0.730043593087919\n",
       "3. 0.618751345207167\n",
       "4. 0.550024750476546\n",
       "5. 0.531629181306629\n",
       "6. 0.765011290422085\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 0.6675085 0.7300436 0.6187513 0.5500248 0.5316292 0.7650113"
      ]
     },
     "execution_count": 409,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lbd = dnorm(z*res$coefficients[2])/pnorm(z*res$coefficients[2])\n",
    "head(lbd)\n",
    "df    <- data.frame(df ,lbd= lbd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 410,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "lm(formula = Ys ~ x + lbd, data = df)\n",
       "\n",
       "Residuals:\n",
       "     Min       1Q   Median       3Q      Max \n",
       "-2.32577 -0.45767  0.01295  0.49208  2.51733 \n",
       "\n",
       "Coefficients:\n",
       "            Estimate Std. Error t value Pr(>|t|)    \n",
       "(Intercept)  1.06187    0.06723   15.79   <2e-16 ***\n",
       "x            1.02585    0.02249   45.62   <2e-16 ***\n",
       "lbd         -0.06732    0.07654   -0.88    0.379    \n",
       "---\n",
       "Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n",
       "\n",
       "Residual standard error: 0.7228 on 997 degrees of freedom\n",
       "Multiple R-squared:  0.6766,\tAdjusted R-squared:  0.6759 \n",
       "F-statistic:  1043 on 2 and 997 DF,  p-value: < 2.2e-16\n"
      ]
     },
     "execution_count": 410,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res <- lm(Ys ~  x + lbd, data = df)\n",
    "summary(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (c) 2-step estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 411,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "library(\"sampleSelection\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 412,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a <- heckit( D ~  z, Ys ~ x , data=df,method = \"2step\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 413,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>Estimate</th><th scope=col>Std. Error</th><th scope=col>t value</th><th scope=col>Pr(>|t|)</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>(Intercept)</th><td>5.181239e-01</td><td>4.371983e-02</td><td>1.185101e+01</td><td>2.123958e-30</td></tr>\n",
       "\t<tr><th scope=row>z</th><td>4.621762e-01</td><td>4.561822e-02</td><td>1.013140e+01</td><td>5.029723e-23</td></tr>\n",
       "\t<tr><th scope=row>(Intercept)</th><td>1.071536e+00</td><td>6.688249e-02</td><td>1.602117e+01</td><td>1.462510e-51</td></tr>\n",
       "\t<tr><th scope=row>x</th><td> 1.012247e+00</td><td> 2.712527e-02</td><td> 3.731750e+01</td><td>2.915670e-191</td></tr>\n",
       "\t<tr><th scope=row>invMillsRatio</th><td>-0.009687247</td><td> 0.127574552</td><td>-0.075934009</td><td> 0.939486879</td></tr>\n",
       "\t<tr><th scope=row>sigma</th><td>0.7163367</td><td>       NA</td><td>       NA</td><td>       NA</td></tr>\n",
       "\t<tr><th scope=row>rho</th><td>-0.01352332</td><td>         NA</td><td>         NA</td><td>         NA</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|llll}\n",
       "  & Estimate & Std. Error & t value & Pr(>|t|)\\\\\n",
       "\\hline\n",
       "\t(Intercept) & 5.181239e-01 & 4.371983e-02 & 1.185101e+01 & 2.123958e-30\\\\\n",
       "\tz & 4.621762e-01 & 4.561822e-02 & 1.013140e+01 & 5.029723e-23\\\\\n",
       "\t(Intercept) & 1.071536e+00 & 6.688249e-02 & 1.602117e+01 & 1.462510e-51\\\\\n",
       "\tx &  1.012247e+00 &  2.712527e-02 &  3.731750e+01 & 2.915670e-191\\\\\n",
       "\tinvMillsRatio & -0.009687247 &  0.127574552 & -0.075934009 &  0.939486879\\\\\n",
       "\tsigma & 0.7163367 &        NA &        NA &        NA\\\\\n",
       "\trho & -0.01352332 &          NA &          NA &          NA\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "1. 0.518123941194385\n",
       "2. 0.462176159516493\n",
       "3. 1.07153552462255\n",
       "4. 1.01224739002373\n",
       "5. -0.00968724720762969\n",
       "6. 0.716336679918967\n",
       "7. -0.0135233158920824\n",
       "8. 0.0437198273609888\n",
       "9. 0.0456182153400582\n",
       "10. 0.0668824927414191\n",
       "11. 0.0271252735521768\n",
       "12. 0.127574552070074\n",
       "13. NA\n",
       "14. NA\n",
       "15. 11.8510061102553\n",
       "16. 10.1313950155925\n",
       "17. 16.0211660884907\n",
       "18. 37.3174997876655\n",
       "19. -0.0759340091769141\n",
       "20. NA\n",
       "21. NA\n",
       "22. 2.12395783242378e-30\n",
       "23. 5.0297230434052e-23\n",
       "24. 1.46251047722009e-51\n",
       "25. 2.91567010804642e-191\n",
       "26. 0.939486879428983\n",
       "27. NA\n",
       "28. NA\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "                  Estimate Std. Error     t value      Pr(>|t|)\n",
       "(Intercept)    0.518123941 0.04371983 11.85100611  2.123958e-30\n",
       "z              0.462176160 0.04561822 10.13139502  5.029723e-23\n",
       "(Intercept)    1.071535525 0.06688249 16.02116609  1.462510e-51\n",
       "x              1.012247390 0.02712527 37.31749979 2.915670e-191\n",
       "invMillsRatio -0.009687247 0.12757455 -0.07593401  9.394869e-01\n",
       "sigma          0.716336680         NA          NA            NA\n",
       "rho           -0.013523316         NA          NA            NA"
      ]
     },
     "execution_count": 413,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coef(summary(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Std. error of (b): 0.02250  ・・・(1)  \n",
    "Std. error of (c): 0.02712  ・・・(2)**\n",
    "  \n",
    "Here we can see that two std. errors are different.  In this case (2) > (1) , however it is (1) > (2) in STATA as I saw in some of my friends result.  It is because  when we calculate standard error we must correct an asymptoticvariance.  In this case , heckit, we should replace 2-stage estimatour with GMM estimatour.  STATA does thi automatically, thus (1) > (2) happens in STATA.\n",
    "\n",
    " Remark: I have no idea why the opposite happened in R. Perhaps it is because glm() command in R uses more appropriate way to in its inside."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (a) Convergence in Probability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 421,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"N = 1\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>z:</strong> NA"
      ],
      "text/latex": [
       "\\textbf{z:} NA"
      ],
      "text/markdown": [
       "**z:** NA"
      ],
      "text/plain": [
       " z \n",
       "NA "
      ]
     },
     "execution_count": 421,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"N = 10\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>z:</strong> 1.04609604217411"
      ],
      "text/latex": [
       "\\textbf{z:} 1.04609604217411"
      ],
      "text/markdown": [
       "**z:** 1.04609604217411"
      ],
      "text/plain": [
       "       z \n",
       "1.046096 "
      ]
     },
     "execution_count": 421,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"N = 100\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>z:</strong> 0.372259396447066"
      ],
      "text/latex": [
       "\\textbf{z:} 0.372259396447066"
      ],
      "text/markdown": [
       "**z:** 0.372259396447066"
      ],
      "text/plain": [
       "        z \n",
       "0.3722594 "
      ]
     },
     "execution_count": 421,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"N = 10000\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>z:</strong> 0.462176132791821"
      ],
      "text/latex": [
       "\\textbf{z:} 0.462176132791821"
      ],
      "text/markdown": [
       "**z:** 0.462176132791821"
      ],
      "text/plain": [
       "        z \n",
       "0.4621761 "
      ]
     },
     "execution_count": 421,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"N = 100000\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>z:</strong> 0.50234199975966"
      ],
      "text/latex": [
       "\\textbf{z:} 0.50234199975966"
      ],
      "text/markdown": [
       "**z:** 0.50234199975966"
      ],
      "text/plain": [
       "       z \n",
       "0.502342 "
      ]
     },
     "execution_count": 421,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"N = 1000000\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>z:</strong> 0.499639769039458"
      ],
      "text/latex": [
       "\\textbf{z:} 0.499639769039458"
      ],
      "text/markdown": [
       "**z:** 0.499639769039458"
      ],
      "text/plain": [
       "        z \n",
       "0.4996398 "
      ]
     },
     "execution_count": 421,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"N = 10000000\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>z:</strong> 0.50233133991461"
      ],
      "text/latex": [
       "\\textbf{z:} 0.50233133991461"
      ],
      "text/markdown": [
       "**z:** 0.50233133991461"
      ],
      "text/plain": [
       "        z \n",
       "0.5023313 "
      ]
     },
     "execution_count": 421,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "test <- function(n) {\n",
    "set.seed(10)\n",
    "\n",
    "mu = c(0, 0)\n",
    "A <- array(c(0.5, 0.25, 0.25, 4 ), dim=c(2,2))\n",
    "w = rmvnorm(n, mu, A)\n",
    "u = w[,1]\n",
    "v = w[,2]\n",
    "x = rnorm(n, mean = 0 ,sd = 1)\n",
    "z = rnorm(n, mean = 0 ,sd = 1)\n",
    "df    <- data.frame(u=u, v=v, x=x , z =z) \n",
    "D = sign(1 + z + v)\n",
    "Ys = 1 + x + u\n",
    "df    <- data.frame(df ,dc =1 + z + v, D =D, Ys = Ys)\n",
    "df$D[df$D == 0] = 1 \n",
    "df$D[df$D == -1] = 0 \n",
    "res <- glm(D ~  z , data = df, family=binomial(probit))\n",
    "return(res)\n",
    "\n",
    "}\n",
    "\n",
    "print(\"N = 1\")\n",
    "test(1)$coefficients[2]\n",
    "print(\"N = 10\")\n",
    "test(10)$coefficients[2]\n",
    "print(\"N = 100\")\n",
    "test(100)$coefficients[2]\n",
    "print(\"N = 10000\")\n",
    "test(1000)$coefficients[2]\n",
    "print(\"N = 100000\")\n",
    "test(10000)$coefficients[2]\n",
    "print(\"N = 1000000\")\n",
    "test(100000)$coefficients[2]\n",
    "print(\"N = 10000000\")\n",
    "test(1000000)$coefficients[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "As shown in above, the coefficient converges to about 0.5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}