NYU Quant Macroeconomics Note.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is created for self-studying for Quantitative Macroeconomics course at NYU, taught by Gianluca Violante. Course website: https://sites.google.com/a/nyu.edu/glviolante/teaching/quantmacro"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "addpath('C:\\Users\\barry\\Documents\\Academics\\Quantecon-Notes\\functions')\n",
    "cd 'C:\\Users\\barry\\Documents\\Academics\\Quantecon-Notes'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 2\n",
    "\n",
    "## Differentiation\n",
    "Finite difference: $f'(x_0) \\simeq \\frac{f(x_0+h)-f(x_0)}{h}$\n",
    "\n",
    "Two-sided derivative: $f'(x_0) \\simeq \\frac{f(x_0+h)-f(x_0-h)}{2h}$. It is better as long as the fourth order residual on both sides of $x_0$ are of similar magnitude (prove using Taylor expansion)\n",
    "\n",
    "In practice, compute TS as $f'(x_0) \\simeq \\frac{f(x_0+h)-f(x_0-h)}{x_0+h - (x_0-h)}$ so both have same argument. \n",
    "\n",
    "Miranda-Fackler: set $h = max(|x_0|,1)\\sqrt[3]{\\epsilon} \\simeq \\max(|x_0,1|)\\times 6 \\times 10^{-6}$, where $\\epsilon$ is margion epsilon (for 64-bit machine, $\\epsilon = 2^{-52}$).\n",
    "\n",
    "**Example: Two-sided finite difference for single variable function**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created file 'C:\\Users\\barry\\Documents\\Academics\\Quantecon-Notes\\TS_diff.m'.\n"
     ]
    }
   ],
   "source": [
    "%%file TS_diff.m\n",
    "function TS = TS_diff(f,x0)\n",
    "    \n",
    "    h = max(abs(x0),1)*6.0*1e-6;\n",
    "    \n",
    "    TS = (f(x0+h)-f(x0-h))/(2*h);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ans =\n",
      "\n",
      "    7.3891\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "TS_diff(@exp, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example: Two-sided finite difference for Jacobian of an arbitrary function**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created file 'C:\\Users\\barry\\Documents\\Academics\\Quantecon-Notes\\functions\\TS_jac.m'.\n"
     ]
    }
   ],
   "source": [
    "%%file functions/TS_jac.m\n",
    "function fjac = TS_jac(f,x)\n",
    "    \n",
    "    h = max(abs(x),1)*6.0*1e-6;\n",
    "    xh1 = x+h;\n",
    "    xh0 = x-h;\n",
    "    h = xh1-xh0;\n",
    "    \n",
    "    for j =1:length(x)\n",
    "        xx = x; xx(j) = xh1(j); f1 = feval(f,xx);\n",
    "        xx(j) = xh0(j); f0 = feval(f,xx);\n",
    "        fjac(:,j) = (f1-f0)/h(j);\n",
    "    end\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ans =\n",
      "\n",
      "    0.5403         0         0         0         0\n",
      "         0   -0.4161         0         0         0\n",
      "         0         0   -0.9900         0         0\n",
      "         0         0         0    0.2837         0\n",
      "         0         0         0         0    0.5403\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "x = [1,2,3,5,1];\n",
    "TS_jac(@sin, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Root Finding\n",
    "\n",
    "(See other notebook)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The matching model of the labor market\n",
    "\n",
    "Model setup: risk-neutral households that discount future at rate $\\beta = 1/R$. Linear production function, which we can therefore think 1 firm only hires 1 worker (homogeneity). Aggregate productivity $z$ draws from $\\Gamma(z',z)$ with $N$ values ($z$ is current state and $z'$ is any possible future state. $\\Gamma$ is transition matrix). Employed workers receive a wage $w(z)$ and separate at rate $\\delta$ (rate for leaving the firm). \n",
    "\n",
    "Firm's profit: $z-w(z)$\n",
    "\n",
    "Unemployed workers receive benefit $b$ and search for vacancies.\n",
    "\n",
    "In each period there are fraction $u$ of unemployed workers, and $v$ vacant jobs. Number of matches taking place per unit of time is given by aggregate CRS (constant return to scale) matching function: $m = m(u,v)$.\n",
    "\n",
    "Unemployed worker moving to employment as a Poisson process with rate $p(\\theta) = m/u$, and job vacancy filled with rate $q(\\theta) = m/v$. $\\theta = v/u$ is called labor market tightness. Flow cost for firm of posting a vacancy is $\\kappa$\n",
    "\n",
    "Value of unemployed and employed workers at state $z$: \n",
    "\n",
    "$$U(z) = b + \\beta \\mathbb{E}_z[p(\\theta)W(z') + (1-p(\\theta))U(z')]$$\n",
    "\n",
    "$$W(z) = w(z) + \\beta \\mathbb{E}_z[\\delta U(z') + (1-\\delta)W(z')]$$\n",
    "\n",
    "Value of vacant and matched firms at state $z$:\n",
    "\n",
    "$$V(z) = -\\kappa + \\beta \\mathbb{E}_z[q(\\theta)J(v') + (1-q(\\theta))V(z')]$$\n",
    "\n",
    "$$J(z) = z - w(z) + \\beta \\mathbb{E}_z[\\delta V(z') + (1-\\delta) J(z')]$$\n",
    "\n",
    "Free entry of the firm: $V(z) = 0 \\forall z$\n",
    "\n",
    "**Solving for $\\theta$:**\n",
    "\n",
    "Wage $w(z)$ is determined by Nash bargaining between worker and the firm. Wage solves $\\eta$ is bargaining power for worker\n",
    "\n",
    "$$\\max_{w(z)} [W(z)-U(z)]^\\eta[J(z)-V(z)]^{1-\\eta}$$\n",
    "\n",
    "Solution is $w(z) = \\eta z + (1-\\eta)b + \\eta \\kappa \\theta(z)$ (use chain rule to take FOC). \n",
    "\n",
    "After some algebra (use V(z) = 0):\n",
    "\n",
    "$$\\frac{\\kappa}{\\beta q(\\theta(z_i))} = \\sum_{j=1}^N [(1-\\eta)(z_i-b)-\\eta\\kappa\\theta(z_j)+\\frac{(1-\\delta)\\kappa}{q(\\theta(z_j))}]\\Gamma(z_j,z_i)$$\n",
    "\n",
    "which is a system of $N$ nonlinear equations in $N$ unknowns $\\theta(z_i)$. Given $u_t$ and solve for $\\theta_t$, we can simulate $v_t = u_t \\theta_t$, $w_t$, and recall the law of motion for the endogenous state $u_{t+1} = u_t + \\delta(1-u_t) - p(\\theta_t)u(t)$   \n",
    "\n",
    "**Full solution in homework (Todo)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example: Cournot duopoly**\n",
    "\n",
    "The inverse demand function is $p=q^{-1/\\eta}$ and the two firms faces cost function $C_i(q_i)=\\frac{1}{2}c_iq_i^2$ for $i=1,2$. The profit of the firm is $\\pi_i(q_1, q_2) = P(q_1+q_2)q_i - C(q_i)$ and it solves the FOC: $\\frac{\\partial \\pi_i}{\\partial q_i} = 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created file 'C:\\Users\\barry\\Documents\\Academics\\Quantecon-Notes\\functions\\cournot.m'.\n"
     ]
    }
   ],
   "source": [
    "%%file functions/cournot.m\n",
    "\n",
    "function fval = cournot(q)\n",
    "    c = [0.6;0.8]; eta = 1.6;\n",
    "    e = -1/eta;\n",
    "    fval = sum(q)^e + e*sum(q)^(e-1)*q - diag(c)*q;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Equation solved.\n",
      "\n",
      "fsolve completed because the vector of function values is near zero\n",
      "as measured by the value of the function tolerance, and\n",
      "the problem appears regular as measured by the gradient.\n",
      "\n",
      "\n",
      "ans =\n",
      "\n",
      "    0.8396\n",
      "    0.6888\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "fsolve(@cournot,[0.2;0.2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integration\n",
    "\n",
    "Overview: three methods\n",
    "1. Newton-Cotes methods: employ piecewise polynomial approximations to the integrand with evenly spaced nodes\n",
    "2. Gaussian Quadrature: choose modes and weights that satisfy some moment-matching conditions\n",
    "3. Monte-Carlo methods: use equally weighted random nodes\n",
    "\n",
    "### Newton-Cotes: \n",
    "Basic idea: approximate $f$ with a piecewise linear polynomial $\\tilde{f}$ whose integral is easy to compute. \n",
    "\n",
    "Trapezoid rule: $f(x) \\simeq f(x_i) + \\frac{f(x_{i+1})-f(x_i)}{x_{i+1} - x_i}(x-x_i)$ and $\\int_a^b f(x)dx = \\sum_{i=0}^N w_i \\frac{f(x_{i+1})+f(x_i)}{2}$ where $w_i = x_{i+1} - x_i$. At the end use $w_i = 1/2$\n",
    "\n",
    "Simpson rule: based on quadratic approximation. For smooth function, Simpson rule is better than trapezoid rule. But for nondifferentiable function, use trapezoid.\n",
    "\n",
    "### Gaussian quadrature\n",
    "\n",
    "A function can be approximated by a polynomial, and weights and quadrature points can be chosen efficiently. \n",
    "\n",
    "Let $\\Phi = {\\phi_k(x)}^\\infty_{k=0}$ be a family of polynomials defined on $[a,b]$, $\\phi_k(x) = \\sum_{i=0}^k \\alpha_{k,i}x^i$.\n",
    "\n",
    "Definition: $\\Phi$ is orghogonal wrt the weight function $w$ if $<\\phi_k,\\phi_j> = \\int_a^b \\phi_k(x)\\phi_j(x)w(x)dx = 0$ for all $k \\neq j$. If $<\\phi_k,\\phi_k> = 1$, then it is called orthonormal. \n",
    "\n",
    "Theorem: Suppose $\\Phi$ is orghonormal wrt $w(x)$ and ${x_i}_1^n$ be the zeros of $\\phi_n(x)$. Then $$\\int_a^b f(x)w(x)dx = \\sum_{i=1}^n w_if(x_i) + \\frac{f(\\xi)^{2n}}{2n!}$$\n",
    "for some $\\xi \\in [a,b]$ and the weights are given by $$w_i = -\\frac{1}{\\phi_n'(x_i)\\phi_{n+1}(x_i)} > 0$$\n",
    "\n",
    "The thoerem gives zeros of orthonormal polynomials as quadrature nodes weights. In practice, use well-known polynomial families. \n",
    "\n",
    "<img src=\"misc/polynomial quadrature.PNG\" style=\"width: 500px;\">\n",
    "\n",
    "### Monte-Carlo\n",
    "\n",
    "Useful for high dimensional integration. \n",
    "\n",
    "Basic idea: from LLN, $\\lim_{n\\rightarrow\\infty}\\frac{1}{n}\\sum_{i=1}^n f(x_i) = \\int_0^1 f(x)dx$. \n",
    "\n",
    "Method: draw $n$ realizations of $x$ and compute the sum $I_n = \\frac{1}{n}\\sum_{i=1}^n f(x_i)$\n",
    "\n",
    "Variance-reduction techniques:\n",
    "1. Stratified sampling: divide $[0,1]$ into two subintervals $[0,\\alpha]$ and $[\\alpha, 1]$ in order to avoid concentration of draws in one subinterval from bad luck and use $I_n = \\frac{\\alpha}{n} \\sum_{i=1}^n f(x_{1i}) + \\frac{1-\\alpha}{n}\\sum_{i=1}^n f(x_{2i})$\n",
    "2. Antithetic variance: if f is weakly increasing then $f(x)$ and $f(1-x)$ are negatively correlated. Therefore, compute $I_n = \\frac{1}{2n}\\sum_{i=1}^n [f(x_i) + f(1-x_i)]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Matlab",
   "language": "matlab",
   "name": "matlab"
  },
  "language_info": {
   "codemirror_mode": "octave",
   "file_extension": ".m",
   "help_links": [
    {
     "text": "MetaKernel Magics",
     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "matlab",
   "version": "0.16.7"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}