ps7.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "16c41db4",
   "metadata": {},
   "source": [
    "# Problem Set 7\n",
    "## Daryl Larsen\n",
    "### 3. (Monte Carlo simulations for Probit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a6136fa5",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra, Random, Distributions, PrettyTables, Optim, Plots"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57754248",
   "metadata": {},
   "source": [
    "#### (a) Generate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "847c8ae2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "data (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Set some seed (π sounds good)\n",
    "Random.seed!(314159265358)\n",
    "\n",
    "function data(n; β=[-0.5 , 1], ρ=0.9, δ=[0 0])\n",
    "    # generate standard normal X and ones\n",
    "    X1 = ones(n, 1)\n",
    "    X2 = rand(Normal(0,1),n)\n",
    "    X = [X1 X2]\n",
    "    \n",
    "    # generate standard normal u and Y according to formula\n",
    "    U = rand(Normal(0,1),n)\n",
    "    Y = 1 .* (X*β .+ U .>= 0)\n",
    "    \n",
    "    return ( Y = Y, X = X, U = U)  \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "722614a4",
   "metadata": {},
   "source": [
    "#### (b) Construct sample criterion function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "aae5f12e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sampleCriterion (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function sampleCriterion(n,Y,b,X)\n",
    "    xb = X * b\n",
    "    Φ = zeros(n,1)\n",
    "    Φ1 = zeros(n,1)\n",
    "    for j = 1:n\n",
    "        Φ[j,1] = log(cdf(Normal(0,1),xb[j,1]))\n",
    "        Φ1[j,1] = log(1 - cdf(Normal(0,1),xb[j,1]))\n",
    "    end\n",
    "    Q = - (Y' * Φ) ./n - ((1 .- Y)' * Φ1) ./n\n",
    "    \n",
    "    return Q\n",
    "end "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a60cdba",
   "metadata": {},
   "source": [
    "#### (c) Compute MLE of β0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3604f02c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " -0.5881186381151735\n",
       "  1.0361748198847573"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# generate data\n",
    "n = 100\n",
    "(Y, X, U) = data(n)\n",
    "\n",
    "result = optimize(b->sampleCriterion(n,Y,b,X)[1,1],[0.0 ; 0.0])\n",
    "βprobit = Optim.minimizer(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "304ec14e",
   "metadata": {},
   "source": [
    "#### (d) Logit criterion function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4a1c6951",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "logitCriterion (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function logitCriterion(n,Y,b,X)\n",
    "    xb = X * b\n",
    "    Λ = zeros(n,1)\n",
    "    Λ1 = zeros(n,1)\n",
    "    for j = 1:n\n",
    "        Λ[j,1] = log( exp(xb[j,1]) ./ (1 .+ exp(xb[j,1])))\n",
    "        Λ1[j,1] = log(1 - exp(xb[j,1]) ./ (1 .+ exp(xb[j,1])))\n",
    "    end\n",
    "    Q = - (Y' * Λ) ./n - ((1 .- Y)' * Λ1) ./n\n",
    "    \n",
    "    return Q\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "08c34998",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " -1.0159651032899697\n",
       "  1.8127426924132513"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = optimize(b->logitCriterion(n,Y,b,X)[1,1],[0.0 ; 0.0])\n",
    "βlogit = Optim.minimizer(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62ff6958",
   "metadata": {},
   "source": [
    "#### (e) OLS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "65d764b3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " 0.3456312177026701\n",
       " 0.25026850544971313"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "βOLS = (X'*X)\\X'*Y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b416f2d3",
   "metadata": {},
   "source": [
    "#### (f) Monte Carlo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "63ad7d96",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "┌──────┬───────────┬──────────┬───────────┐\n",
      "│\u001b[1m      \u001b[0m│\u001b[1m    Probit \u001b[0m│\u001b[1m    Logit \u001b[0m│\u001b[1m       OLS \u001b[0m│\n",
      "├──────┼───────────┼──────────┼───────────┤\n",
      "│ Bias │ 0.0480279 │ 0.795611 │ -0.734105 │\n",
      "│ RMSE │  0.228026 │ 0.891947 │   0.73481 │\n",
      "└──────┴───────────┴──────────┴───────────┘\n"
     ]
    }
   ],
   "source": [
    "rep = 10000\n",
    "biasProbit=0\n",
    "biasLogit=0\n",
    "biasOLS=0\n",
    "rmseProbit=0\n",
    "rmseLogit=0\n",
    "rmseOLS=0\n",
    "for i in 1:rep\n",
    "    Random.seed!(i);\n",
    "    (Y, X, U) = data(n)\n",
    "    result = optimize(b->sampleCriterion(n,Y,b,X)[1,1],[0.0 ; 0.0])\n",
    "    βProbit = Optim.minimizer(result)\n",
    "    \n",
    "    result = optimize(b->logitCriterion(n,Y,b,X)[1,1],[0.0 ; 0.0])\n",
    "    βLogit = Optim.minimizer(result)\n",
    "    \n",
    "    βOLS = (X'*X)\\X'*Y\n",
    "    \n",
    "    biasProbit += βProbit[2] -1;\n",
    "    rmseProbit += (βProbit[2] -1)^2;\n",
    "    biasLogit += βLogit[2] -1;\n",
    "    rmseLogit += (βLogit[2] -1)^2\n",
    "    biasOLS += βOLS[2] -1;\n",
    "    rmseOLS += (βOLS[2] -1)^2;\n",
    "end\n",
    "\n",
    "biasProbit = biasProbit/rep\n",
    "biasLogit = biasLogit/rep\n",
    "biasOLS = biasOLS/rep\n",
    "rmseProbit = sqrt(rmseProbit/rep)\n",
    "rmseLogit = sqrt(rmseLogit/rep)\n",
    "rmseOLS = sqrt(rmseOLS/rep)\n",
    "\n",
    "table = [ \"Bias\" biasProbit biasLogit biasOLS;\n",
    "\"RMSE\" rmseProbit rmseLogit rmseOLS]\n",
    "pretty_table(table, [\"\" \"Probit\" \"Logit\" \"OLS\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca25284d",
   "metadata": {},
   "source": [
    "#### (g) Which is best?\n",
    "Probit has a much smaller bias and RMSE than both Logit and OLS. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f514a905",
   "metadata": {},
   "source": [
    "#### (h) Do it again but bigger"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4afc6f01",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "┌──────┬────────────┬──────────┬───────────┐\n",
      "│\u001b[1m      \u001b[0m│\u001b[1m     Probit \u001b[0m│\u001b[1m    Logit \u001b[0m│\u001b[1m       OLS \u001b[0m│\n",
      "├──────┼────────────┼──────────┼───────────┤\n",
      "│ Bias │ 0.00347562 │ 0.718006 │ -0.734958 │\n",
      "│ RMSE │  0.0629207 │ 0.727009 │  0.735026 │\n",
      "└──────┴────────────┴──────────┴───────────┘\n"
     ]
    }
   ],
   "source": [
    "n=1000\n",
    "rep = 10000\n",
    "biasProbit=0\n",
    "biasLogit=0\n",
    "biasOLS=0\n",
    "rmseProbit=0\n",
    "rmseLogit=0\n",
    "rmseOLS=0\n",
    "for i in 1:rep\n",
    "    Random.seed!(i);\n",
    "    (Y, X, U) = data(n)\n",
    "    result = optimize(b->sampleCriterion(n,Y,b,X)[1,1],[0.0 ; 0.0])\n",
    "    βProbit = Optim.minimizer(result)\n",
    "    \n",
    "    result = optimize(b->logitCriterion(n,Y,b,X)[1,1],[0.0 ; 0.0])\n",
    "    βLogit = Optim.minimizer(result)\n",
    "    \n",
    "    βOLS = (X'*X)\\X'*Y\n",
    "    \n",
    "    biasProbit += βProbit[2] -1;\n",
    "    rmseProbit += (βProbit[2] -1)^2;\n",
    "    biasLogit += βLogit[2] -1;\n",
    "    rmseLogit += (βLogit[2] -1)^2\n",
    "    biasOLS += βOLS[2] -1;\n",
    "    rmseOLS += (βOLS[2] -1)^2;\n",
    "end\n",
    "\n",
    "biasProbit = biasProbit/rep\n",
    "biasLogit = biasLogit/rep\n",
    "biasOLS = biasOLS/rep\n",
    "rmseProbit = sqrt(rmseProbit/rep)\n",
    "rmseLogit = sqrt(rmseLogit/rep)\n",
    "rmseOLS = sqrt(rmseOLS/rep)\n",
    "\n",
    "table = [ \"Bias\" biasProbit biasLogit biasOLS;\n",
    "\"RMSE\" rmseProbit rmseLogit rmseOLS]\n",
    "pretty_table(table, [\"\" \"Probit\" \"Logit\" \"OLS\"])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.0",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}