Notebooks/SMM/SMMest.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulated Method of Moments Estimation\n",
    "### by [Richard W. Evans](https://sites.google.com/site/rickecon/), January 2018\n",
    "The code in this Jupyter notebook was written using Python 3.6. It also uses data file `data/Econ381totpts.txt` and the image file `images/MLEplots.png`. For the code to run properly, you should have the data file and the image fil in the same folder as the Jupyter notebook file. Otherwise, you will have to change the respective lines of the code that read in the data to reflect the location of that data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. The SMM estimator\n",
    "Simulated method of moments (SMM) is analogous to the generalized method of moments (GMM) estimator. SMM could really be thought of as a particular type of GMM estimator. The SMM estimator chooses model parameters $\\theta$ to make simulated model moments match data moments. Seminal papers developing SMM are McFadden (1989), Lee and Ingram (1991), and Duffie and Singleton (1993). Good textbook treatments of SMM are found in Adda and Cooper (2003, pp. 87-100) and Davidson and MacKinnon (2004, pp. 383-394).\n",
    "\n",
    "In ML estimation, we used data $x$ and model parameters $\\theta$ to maximize the likelihood of drawing that data $x$ from the model given parameters $\\theta$.\n",
    "\n",
    "$$ \\hat{\\theta}_{ML} = \\theta:\\quad \\max_{\\theta}\\ln\\mathcal{L} = \\sum_{i=1}^N\\ln\\Bigl(f(x_i|\\theta)\\Bigr) $$\n",
    "\n",
    "In GMM estimation, we used data $x$ and model parameters $\\theta$ to minimize the distance between model moments $m(x|\\theta)$ and data moments $m(x)$.\n",
    "\n",
    "$$ \\hat{\\theta}_{GMM} = \\theta:\\quad \\min_{\\theta}||m(x|\\theta) - m(x)|| $$\n",
    "\n",
    "The following difficulties can arise with GMM making it not possible or very difficult.\n",
    "\n",
    "* The model moment function $m(x|\\theta)$ is not known analytically.\n",
    "* The model moments $m(x|\\theta)$ are derived from *latent variables* that are not observed by the modeler. See Laroque and Salanie (1993).\n",
    "* The model moments $m(x|\\theta)$ are derived from *censored variables* that are only partially observed by the modeler.\n",
    "* The model moments $m(x|\\theta)$ are just difficult to derive analytically. Examples include moments that include multiple integrals over nonlinear functions as in McFadden (1989).\n",
    "\n",
    "SMM estimation is simply to simulate the model data $S$ times, and use the average values of the moments from the simulated data as the estimator for the model moments. Let $\\tilde{x}=\\{\\tilde{x}_1,\\tilde{x}_2,...\\tilde{x}_s,...\\tilde{x}_S\\}$ be the $S$ simulations of the model data.  \n",
    "\n",
    "$$ \\hat{m}\\left(\\tilde{x}|\\theta\\right) = \\frac{1}{S}\\sum_{s=1}^S m\\left(\\tilde{x}_s|\\theta\\right) $$\n",
    "\n",
    "Once we have an estimate of the model moments $\\hat{m}\\left(\\tilde{x}|\\theta\\right)$ from our $S$ simulations, SMM estimation is very similar to our presentation of GMM. The SMM approach of estimating the parameter vector $\\hat{\\theta}_{SMM}$ is to choose $\\theta$ to minimize some distance measure of the data moments $m(x)$ from the simulated model moments $\\hat{m}(\\tilde{x}|\\theta)$.\n",
    "\n",
    "$$ \\hat{\\theta}_{SMM}=\\theta:\\quad \\min_{\\theta}\\: ||\\hat{m}(\\tilde{x}|\\theta)-m(x)|| $$\n",
    "\n",
    "The distance measure $||\\hat{m}(\\tilde{x}|\\theta)-m(x)||$ can be any kind of norm. But it is important to recognize that your estimates $\\hat{\\theta}_{SMM}$ will be dependent on what distance measure (norm) you choose. The most widely studied and used distance metric in GMM and SMM estimation is the $L^2$ norm or the sum of squared errors in moments. Define the moment error function $e(\\tilde{x},x|\\theta)$ as the percent difference in the vector of simulated model moments from the data moments.\n",
    "\n",
    "$$ e(\\tilde{x},x|\\theta) \\equiv \\frac{\\hat{m}(\\tilde{x}|\\theta)-m(x)}{m(x)} $$\n",
    "\n",
    "It is important that the error function $e(\\tilde{x},x|\\theta)$ be a percent deviation of the moments (given that none of the data moments are 0). This puts all the moments in the same units, which helps make sure that no moments receive unintended weighting simply due to its units. This ensures that the problem is scaled properly and will suffer from as little as possible ill conditioning.\n",
    "\n",
    "In this case, the SMM estimator is the following,\n",
    "\n",
    "$$ \\hat{\\theta}_{SMM}=\\theta:\\quad \\min_{\\theta}\\:e(\\tilde{x},x|\\theta)^T \\, W \\, e(\\tilde{x},x|\\theta) $$\n",
    "\n",
    "where $W$ is a $R\\times R$ weighting matrix in the criterion function. For now, think of this weighting matrix as the identity matrix. But we will show in Section 2 a more optimal weighting matrix. We call the quadratic form expression $e(\\tilde{x},x|\\theta)^T \\, W \\, e(\\tilde{x},x|\\theta)$ the *criterion function* because it is a strictly positive scalar that is the object of the minimization in the SMM problem statement. The $R\\times R$ weighting matrix $W$ in the criterion function allows the econometrician to control how each moment is weighted in the minimization problem. For example, an $R\\times R$ identity matrix for $W$ would give each moment equal weighting, and the criterion function would be a simply sum of squared percent deviations (errors). Other weighting strategies can be dictated by the nature of the problem or model.\n",
    "\n",
    "One last item to emphasize with SMM, which we will highlight in the examples in this notebook, is that the errors that are drawn for the $S$ simulations of the model must be drawn only once so that the minimization problem for $\\hat{\\theta}_{SMM}$ does not have the underlying sampling changing for each guess of a value of $\\theta$. Put more simply, you want the random draws for all the simulations to be held constant so that the only thing changing in the minimization problem is the value of the vector of parameters $\\theta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. The Weighting Matrix (W)\n",
    "In the SMM criterion function in the problem statement above, some weighting matrices $W$ produce precise estimates while others produce poor estimates with large variances. We want to choose the optimal weighting matrix $W$ with the smallest possible asymptotic variance. This is an efficient or optimal SMM estimator. The optimal weighting matrix is the inverse variance covariance matrix of the moments at the optimal moments,\n",
    "\n",
    "$$ W^{opt} \\equiv \\Omega^{-1}(\\tilde{x},x|\\hat{\\theta}_{SMM}) $$\n",
    "\n",
    "where $\\Omega(\\tilde{x},x|\\theta)$ is the variance covariance matrix of the moment condition errors $e(\\tilde{x},x|\\theta)$. The intuition for using the inverse variance covariance matrix $\\Omega^{-1}$ as the optimal weighting matrix is the following. You want to downweight moments that have a high variance, and you want to weight more heavily the moments that are generated more precisely.\n",
    "\n",
    "Notice that this definition of the optimal weighting matrix is circular. $W^{opt}$ is a function of the SMM estimates $\\hat{\\theta}_{SMM}$, but the optimal weighting matrix is used in the estimation of $\\hat{\\theta}_{SMM}$. This means that one has to use some kind of iterative fixed point method to find the true optimal weighting matrix $W^{opt}$. Below are some examples of weighting matrices to use."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1. The identity matrix (W = I)\n",
    "Many times, you can get away with just using the identity matrix as your weighting matrix $W = I$. This changes the criterion function to a simple sum of squared error functions such that each moment has the same weight.\n",
    "\n",
    "$$ \\hat{\\theta}_{SMM}=\\theta:\\quad \\min_{\\theta}\\:e(\\tilde{x},x|\\theta)^T \\, e(\\tilde{x},x|\\theta) $$\n",
    "\n",
    "If the problem is well conditioned and well identified, then your SMM estimates $\\hat{\\theta}_{SMM}$ will not be greatly affected by this simplest of weighting matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2. Two-step variance covariance estimator of W\n",
    "The most common method of estimating the optimal weighting matrix for SMM estimates is the two-step variance covariance estimator. The name \"two-step\" refers to the two steps used to get the weighting matrix.\n",
    "\n",
    "The first step is to estimate the SMM parameter vector $\\hat{\\theta}_{1,SMM}$ using the simple identity matrix as the weighting matrix $W = I$.\n",
    "\n",
    "$$ \\hat{\\theta}_{1,SMM}=\\theta:\\quad \\min_{\\theta}\\:e(\\tilde{x},x|\\theta)^T \\, I \\, e(\\tilde{x},x|\\theta) $$\n",
    "\n",
    "You then use the vector of moment error functions from thet Step 1 SMM estimate $e(\\tilde{x},x|\\hat{\\theta}_{1,SMM})$ to get a new estimate of the variance covariance matrix of the moment error vector.\n",
    "\n",
    "$$ \\hat{\\Omega}_2 = \\frac{1}{N}e(\\tilde{x},x|\\hat{\\theta}_{1,SMM})\\,e(\\tilde{x},x|\\hat{\\theta}_{1,SMM})^T $$\n",
    "\n",
    "The optimal weighting matrix is the inverse of the two-step variance covariance matrix.\n",
    "\n",
    "$$ \\hat{W}^{two-step} \\equiv \\hat{\\Omega}_2^{-1} $$\n",
    "\n",
    "Lastly, re-estimate the SMM estimator using the optimal two-step weighting matrix.\n",
    "\n",
    "$$ \\hat{\\theta}_{2, SMM}=\\theta:\\quad \\min_{\\theta}\\:e(\\tilde{x},x|\\theta)^T \\, \\hat{W}^{two-step} \\, e(\\tilde{x},x|\\theta) $$\n",
    "\n",
    "$\\hat{\\theta}_{2, SMM}$ is called the two-step SMM estimator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3. Iterated variance covariance estimator of W\n",
    "The truly optimal weighting matrix $W^{opt}$ is the iterated variance-covariance estimator of $W$. This procedure is to just repeat the process described in the two-step SMM estimator until the estimated weighting matrix no longer changes between iterations. Let $i$ index the $i$th iterated SMM estimator,\n",
    "\n",
    "$$ \\hat{\\theta}_{i,SMM}=\\theta:\\quad \\min_{\\theta}\\:e(\\tilde{x},x|\\theta)^T \\, \\hat{W}_{i} \\, e(\\tilde{x},x|\\theta) $$\n",
    "\n",
    "and the $i+1$th estimate of the optimal weighting matrix is defined as the following.\n",
    "\n",
    "$$ \\hat{W}_{i+1} \\equiv \\hat{\\Omega}_{i+1} = \\frac{1}{N}e(\\tilde{x},x|\\hat{\\theta}_{i,SMM})\\,e(\\tilde{x},x|\\hat{\\theta}_{i,SMM})^T $$\n",
    "\n",
    "The iterated SMM estimator is the $\\hat{\\theta}_{i,SMM}$ such that $\\hat{W}_{i+1}$ is very close to $\\hat{W}_{i}$ for some distance metric (norm).\n",
    "\n",
    "$$ \\hat{\\theta}_{it,SMM} = \\hat{\\theta}_{i,SMM}: \\quad || \\hat{W}_{i+1} - \\hat{W}_{i} || < \\varepsilon $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4. Newey-West consistent estimator of $\\Omega$ and W\n",
    "The Newey-West estimator of the optimal weighting matrix and variance covariance matrix is consistent in the presence of heteroskedasticity and autocorrelation in the data (See Newey and West, 1987). Adda and Cooper (2003, p. 82, 89) have a nice exposition of how to compute the Newey-West weighting matrix $\\hat{W}_{nw}$ in the SMM case. See also Davidson and MacKinnon (2004, p. 364). The asymptotic representation of the optimal weighting matrix $\\hat{W}^{opt}$ is the following:\n",
    "\n",
    "$$ \\hat{W}^{opt} = \\lim_{N\\rightarrow\\infty}\\frac{1}{N}\\sum_{i=1}^N \\sum_{l=-\\infty}^\\infty e(x_i|\\theta)e(x_{i-l}|\\theta)^T $$\n",
    "\n",
    "The Newey-West consistend estimator of $\\hat{W}^{opt}$ is:\n",
    "\n",
    "$$ \\hat{W}_{nw} = \\Gamma_{0,N} + \\frac{1}{N}\\frac{1}{S}\\sum_{i=1}^N \\sum_{s=1}^S \\left[m\\left(\\tilde{x}^s_i|\\hat{\\theta}_{SMM}\\right) - \\frac{1}{L}\\sum_{l=1}^L m\\left(\\tilde{x}^l_i|\\hat{\\theta}_{SMM}\\right)\\right]\\cdot\\left[m\\left(\\tilde{x}^s_i|\\hat{\\theta}_{SMM}\\right) - \\frac{1}{L}\\sum_{l=1}^L m\\left(\\tilde{x}^l_i|\\hat{\\theta}_{SMM}\\right)\\right]^T $$\n",
    "\n",
    "where\n",
    "\n",
    "$$ \\Gamma_{v,N} = \\frac{1}{N}\\sum_{i=v+1}^N e(x_i|\\theta)e(x_{i-v}|\\theta)^T $$\n",
    "\n",
    "Of course, for autocorrelation, the subscript $i$ and the sample size $N$ can be changed to $t$ and $T$, respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Examples\n",
    "In this section, we will use SMM to estimate parameters of the models from the [maximum likelihood notebook](https://github.com/rickecon/StructEst_W17/blob/master/Notebooks/MLE/MLest.ipynb) and from the [GMM notebook](https://github.com/rickecon/StructEst_W17/blob/master/Notebooks/GMM/GMMest.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1. Fitting a truncated normal to intermediate macroeconomics test scores\n",
    "Let's revisit the problem from the ML and GMM notebooks of fitting a truncated normal distribution to intermediate macroeconomics test scores. The data are in the text file `Econ381totpts.txt`. Recall that these test scores are between 0 and 450. The figure below shows a histogram of the data, as well as three truncated normal PDF's. The black line is the ML estimate of $\\mu$ and $\\sigma$ of the truncated normal pdf. The red and the green lines are just the PDF's of two \"arbitrarily\" chosen combinations of the truncated normal parameters $\\mu$ and $\\sigma$.\n",
    "\n",
    "![title](images/MLEplots.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try estimating the parameters $\\mu$ and $\\sigma$ from the truncated normal distribution by SMM. What moments should we use? Let's try the mean and variance of the data. These two statistics of the data are defined by:\n",
    "\n",
    "$$ mean(scores_i) = \\frac{1}{N}\\sum_{i=1}^N scores_i $$\n",
    "\n",
    "$$ var(scores_i) = \\frac{1}{N-1}\\sum_{i=1}^{N} \\left(scores_i - mean(scores_i)\\right)^2 $$\n",
    "\n",
    "So the data moment vector $m(x)$ for SMM is the following.\n",
    "\n",
    "$$ m(scores_i) \\equiv \\begin{bmatrix} mean(scores_i) \\\\ var(scores_i) \\end{bmatrix} $$\n",
    "\n",
    "And the model moment vector $m(x|\\theta)$ for GMM is the following.\n",
    "\n",
    "$$ m(scores_i|\\mu,\\sigma) \\equiv \\begin{bmatrix} mean(scores_i|\\mu,\\sigma) \\\\ var(scores_i|\\mu,\\sigma) \\end{bmatrix} $$\n",
    "\n",
    "But let's assume that we need to simulate the data from the model (test scores) $S$ times in order to get the model moments. In this case, we don't need to simulate. But we will do so to show how SMM works."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What would one simulation (the $s$th simulation) of the test scores look like? There are 161 test score observations in the data file `Econ381totpts.txt`. So one simulation (the $s$th simulation) would be a draw of 161 test scores from a truncated normal distribution with parameters $\\mu$, $\\sigma$, and $cutoff=450$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import packages and load the data\n",
    "import numpy as np\n",
    "import numpy.random as rnd\n",
    "import numpy.linalg as lin\n",
    "import scipy.stats as sts\n",
    "import scipy.integrate as intgr\n",
    "import scipy.optimize as opt\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm\n",
    "cmap1 = matplotlib.cm.get_cmap('summer')\n",
    "# This next command is specifically for Jupyter Notebook\n",
    "%matplotlib notebook\n",
    "\n",
    "pts = np.loadtxt('data/Econ381totpts.txt')\n",
    "# pts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let random variable $y\\sim N(\\mu,\\sigma)$ be distributed normally with mean $\\mu$ and standard deviation $\\sigma$ with PDF given by $\\phi(y|\\mu,\\sigma)$ and CDF given by $\\Phi(y|\\mu,\\sigma)$. The truncated normal distribution of random variable $x\\in(a,b)$ based on $y$ but with cutoff values of $a\\geq -\\infty$ as a lower bound and $a < b\\leq\\infty$ as an upper bound has the following probability density function.\n",
    "\n",
    "$$ f(x|\\mu,\\sigma,a,b) = \\begin{cases} 0 \\quad\\text{if}\\quad x\\leq a \\\\ \\frac{\\phi(x|\\mu,\\sigma)}{\\Phi(b|\\mu,\\sigma) - \\Phi(a|\\mu,\\sigma)}\\quad\\text{if}\\quad a < x < b \\\\ 0 \\quad\\text{if}\\quad x\\geq b \\end{cases} $$\n",
    "\n",
    "The CDF of the truncated normal can be shown to be the following:\n",
    "\n",
    "$$ F(x|\\mu,\\sigma,a,b) = \\begin{cases} 0 \\quad\\text{if}\\quad x\\leq a \\\\ \\frac{\\Phi(x|\\mu,\\sigma) - \\Phi(a|\\mu,\\sigma)}{\\Phi(b|\\mu,\\sigma) - \\Phi(a|\\mu,\\sigma)}\\quad\\text{if}\\quad a < x < b \\\\ 0 \\quad\\text{if}\\quad x\\geq b \\end{cases} $$\n",
    "\n",
    "The inverse CDF of the truncated normal takes a value $p$ between 0 and 1 and solves for the value of $x$ for which $p=F(x|\\mu,\\sigma,a,b)$. The expression for the inverse CDF of the truncated normal is the following:\n",
    "\n",
    "$$ x = \\Phi^{-1}(z|\\mu,\\sigma) \\quad\\text{where}\\quad z = p\\Bigl[\\Phi(b|\\mu,\\sigma) - \\Phi(a|\\mu,\\sigma)\\Bigr] + \\Phi(a|\\mu,\\sigma) $$\n",
    "\n",
    "Note that $z$ is just a transformation of $p$ such that $z\\sim U\\Bigl(\\Phi^{-1}(a|\\mu,\\sigma), \\Phi^{-1}(b|\\mu,\\sigma)\\Bigr)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define function that gives PDF values from truncated normal distribution\n",
    "def trunc_norm_pdf(xvals, mu, sigma, cut_lb, cut_ub):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    Generate PDF values from a truncated normal distribution based on a\n",
    "    normal distribution with mean mu and standard deviation sigma and\n",
    "    cutoffs (cut_lb, cut_ub).\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    xvals  = (N, S) matrix, (N,) vector, or scalar in (cut_lb, cut_ub),\n",
    "             value(s) in the support of s~TN(mu, sig, cut_lb, cut_ub)\n",
    "    mu     = scalar, mean of the nontruncated normal distribution from\n",
    "             which the truncated normal is derived\n",
    "    sigma  = scalar > 0, standard deviation of the nontruncated normal\n",
    "             distribution from which the truncated normal is derived\n",
    "    cut_lb = scalar or string, ='None' if no lower bound cutoff is\n",
    "             given, otherwise is scalar lower bound value of\n",
    "             distribution. Values below this cutoff have zero\n",
    "             probability\n",
    "    cut_ub = scalar or string, ='None' if no upper bound cutoff is given\n",
    "             given, otherwise is scalar lower bound value of\n",
    "             distribution. Values below this cutoff have zero\n",
    "             probability\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        scipy.stats.norm()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    cut_ub_cdf = scalar in [0, 1], cdf of N(mu, sigma) at upper bound\n",
    "                 cutoff of truncated normal distribution\n",
    "    cut_lb_cdf = scalar in [0, 1], cdf of N(mu, sigma) at lower bound\n",
    "                 cutoff of truncated normal distribution\n",
    "    unif2_vals = (N, S) matrix, (N,) vector, or scalar in (0,1),\n",
    "                 rescaled uniform derived from original.\n",
    "    pdf_vals   = (N, S) matrix, (N,) vector, or scalar in (0,1), PDF\n",
    "                 values corresponding to xvals from truncated normal PDF\n",
    "                 with base normal normal distribution N(mu, sigma) and\n",
    "                 cutoffs (cut_lb, cut_ub)\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: pdf_vals\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    # No cutoffs: truncated normal = normal\n",
    "    if (cut_lb == None) & (cut_ub == None):\n",
    "        cut_ub_cdf = 1.0\n",
    "        cut_lb_cdf = 0.0\n",
    "    # Lower bound truncation, no upper bound truncation\n",
    "    elif (cut_lb != None) & (cut_ub == None):\n",
    "        cut_ub_cdf = 1.0\n",
    "        cut_lb_cdf = sts.norm.cdf(cut_lb, loc=mu, scale=sigma)\n",
    "    # Upper bound truncation, no lower bound truncation\n",
    "    elif (cut_lb == None) & (cut_ub != None):\n",
    "        cut_ub_cdf = sts.norm.cdf(cut_ub, loc=mu, scale=sigma)\n",
    "        cut_lb_cdf = 0.0\n",
    "    # Lower bound and upper bound truncation\n",
    "    elif (cut_lb != None) & (cut_ub != None):\n",
    "        cut_ub_cdf = sts.norm.cdf(cut_ub, loc=mu, scale=sigma)\n",
    "        cut_lb_cdf = sts.norm.cdf(cut_lb, loc=mu, scale=sigma)\n",
    "    \n",
    "    pdf_vals = (sts.norm.pdf(xvals, loc=mu, scale=sigma) /\n",
    "                (cut_ub_cdf - cut_lb_cdf))\n",
    "    \n",
    "    return pdf_vals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define function that draws N x S test score values from a truncated\n",
    "# normal distribution\n",
    "def trunc_norm_draws(unif_vals, mu, sigma, cut_lb, cut_ub):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    Draw (N x S) matrix of random draws from a truncated normal\n",
    "    distribution based on a normal distribution with mean mu and\n",
    "    standard deviation sigma and cutoffs (cut_lb, cut_ub). These draws\n",
    "    correspond to an (N x S) matrix of randomly generated draws from a\n",
    "    uniform distribution U(0,1).\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    unif_vals = (N, S) matrix, (N,) vector, or scalar in (0,1), random\n",
    "                draws from uniform U(0,1) distribution\n",
    "    mu        = scalar, mean of the nontruncated normal distribution\n",
    "                from which the truncated normal is derived\n",
    "    sigma     = scalar > 0, standard deviation of the nontruncated\n",
    "                normal distribution from which the truncated normal is\n",
    "                derived\n",
    "    cut_lb    = scalar or string, ='None' if no lower bound cutoff is\n",
    "                given, otherwise is scalar lower bound value of\n",
    "                distribution. Values below this cutoff have zero\n",
    "                probability\n",
    "    cut_ub    = scalar or string, ='None' if no upper bound cutoff is\n",
    "                given, otherwise is scalar lower bound value of\n",
    "                distribution. Values below this cutoff have zero\n",
    "                probability\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        scipy.stats.norm()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    cut_ub_cdf  = scalar in [0, 1], cdf of N(mu, sigma) at upper bound\n",
    "                  cutoff of truncated normal distribution\n",
    "    cut_lb_cdf  = scalar in [0, 1], cdf of N(mu, sigma) at lower bound\n",
    "                  cutoff of truncated normal distribution\n",
    "    unif2_vals  = (N, S) matrix, (N,) vector, or scalar in (0,1),\n",
    "                  rescaled uniform derived from original.\n",
    "    tnorm_draws = (N, S) matrix, (N,) vector, or scalar in (0,1),\n",
    "                  values drawn from truncated normal PDF with base\n",
    "                  normal distribution N(mu, sigma) and cutoffs\n",
    "                  (cut_lb, cut_ub)\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: tnorm_draws\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    # No cutoffs: truncated normal = normal\n",
    "    if (cut_lb == None) & (cut_ub == None):\n",
    "        cut_ub_cdf = 1.0\n",
    "        cut_lb_cdf = 0.0\n",
    "    # Lower bound truncation, no upper bound truncation\n",
    "    elif (cut_lb != None) & (cut_ub == None):\n",
    "        cut_ub_cdf = 1.0\n",
    "        cut_lb_cdf = sts.norm.cdf(cut_lb, loc=mu, scale=sigma)\n",
    "    # Upper bound truncation, no lower bound truncation\n",
    "    elif (cut_lb == None) & (cut_ub != None):\n",
    "        cut_ub_cdf = sts.norm.cdf(cut_ub, loc=mu, scale=sigma)\n",
    "        cut_lb_cdf = 0.0\n",
    "    # Lower bound and upper bound truncation\n",
    "    elif (cut_lb != None) & (cut_ub != None):\n",
    "        cut_ub_cdf = sts.norm.cdf(cut_ub, loc=mu, scale=sigma)\n",
    "        cut_lb_cdf = sts.norm.cdf(cut_lb, loc=mu, scale=sigma)\n",
    "    \n",
    "    unif2_vals = unif_vals * (cut_ub_cdf - cut_lb_cdf) + cut_lb_cdf\n",
    "    tnorm_draws = sts.norm.ppf(unif2_vals, loc=mu, scale=sigma)\n",
    "    \n",
    "    return tnorm_draws"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What would one simulation of 161 test scores look like from a truncated normal with mean $\\mu=300$, $\\sigma=30$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_1 = 300.0\n",
    "sig_1 = 30.0\n",
    "cut_lb_1 = 0.0\n",
    "cut_ub_1 = 450.0\n",
    "unif_vals_1 = sts.uniform.rvs(0, 1, size=161)\n",
    "draws_1 = trunc_norm_draws(unif_vals_1, mu_1, sig_1,\n",
    "                           cut_lb_1, cut_ub_1)\n",
    "print('Mean score =', draws_1.mean())\n",
    "print('Variance of scores =', draws_1.var())\n",
    "print('Standard deviation of scores =', draws_1.std())\n",
    "\n",
    "# Plot data histogram vs. simulated data histogram\n",
    "count_d, bins_d, ignored_d = \\\n",
    "    plt.hist(pts, 30, normed=True, color='b', edgecolor='black',\n",
    "             linewidth=0.8, label='Data')\n",
    "count_m, bins_m, ignored_m = \\\n",
    "    plt.hist(draws_1, 30, normed=True, color='r', edgecolor='black',\n",
    "             linewidth=0.8, label='Simulated data')\n",
    "xvals = np.linspace(0, 450, 500)\n",
    "plt.plot(xvals, trunc_norm_pdf(xvals, mu_1, sig_1, cut_lb_1, cut_ub_1),\n",
    "         linewidth=2, color='k', label='PDF')\n",
    "plt.title('Econ 381 scores: 2011-2012', fontsize=20)\n",
    "plt.xlabel('Total points')\n",
    "plt.ylabel('Percent of scores')\n",
    "plt.xlim([0, 550])  # This gives the xmin and xmax to be plotted\"\n",
    "plt.legend(loc='upper left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From that simulation, we can calculate moments from the simulated data just like we did from the actual data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def data_moments(xvals):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the two data moments for SMM\n",
    "    (mean(data), variance(data)) from both the actual data and from the\n",
    "    simulated data.\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    xvals = (N, S) matrix, (N,) vector, or scalar in (cut_lb, cut_ub),\n",
    "            test scores data, either real world or simulated. Real world\n",
    "            data will come in the form (N,). Simulated data comes in the\n",
    "            form (N,) or (N, S).\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    mean_data = scalar or (S,) vector, mean value of test scores data\n",
    "    var_data  = scalar > 0 or (S,) vector, variance of test scores data\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: mean_data, var_data\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    if xvals.ndim == 1:\n",
    "        mean_data = xvals.mean()\n",
    "        var_data = xvals.var()\n",
    "    elif xvals.ndim == 2:\n",
    "        mean_data = xvals.mean(axis=0)\n",
    "        var_data = xvals.var(axis=0)\n",
    "    \n",
    "    return mean_data, var_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_data, var_data = data_moments(pts)\n",
    "print('Data mean =', mean_data)\n",
    "print('Data variance =', var_data)\n",
    "mean_sim, var_sim = data_moments(draws_1)\n",
    "print('Sim. mean =', mean_sim)\n",
    "print('Sim. variance =', var_sim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also simulate many $(S)$ data sets of test scores, each with $N=161$ test scores. The estimate of the model moments will be the average of the simulated data moments across the simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 161\n",
    "S = 100\n",
    "mu_2 = 300.0\n",
    "sig_2 = 30.0\n",
    "cut_lb = 0.0\n",
    "cut_ub = 450.0\n",
    "unif_vals_2 = sts.uniform.rvs(0, 1, size=(N, S))\n",
    "draws_2 = trunc_norm_draws(unif_vals_2, mu_2, sig_2,\n",
    "                           cut_lb, cut_ub)\n",
    "\n",
    "mean_sim, var_sim = data_moments(draws_2)\n",
    "print(mean_sim)\n",
    "print(var_sim)\n",
    "mean_mod = mean_sim.mean()\n",
    "var_mod = var_sim.mean()\n",
    "print('Estimated model mean =', mean_mod)\n",
    "print('Estimated model variance =', var_mod)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our SMM model moments $\\hat{m}(\\tilde{scores}_i|\\mu,\\sigma)$ are an estimate of the true models moments that we got in the GMM case by integrating using the PDF of the truncated normal distribution. Our SMM moments we got by simulating the data $S$ times and taking the average of the simulated data moments across the simulations as our estimator of the model moments.\n",
    "\n",
    "Define the error vector as the vector of percent deviations of the model moments from the data moments.\n",
    "\n",
    "$$ e(\\tilde{scores}_i,scores_i|\\mu,\\sigma) \\equiv \\frac{\\hat{m}(\\tilde{scores}_i|\\mu,\\sigma) - m(scores_i)}{m(scores_i)} $$\n",
    "\n",
    "The SMM estimator for this moment vector is the following.\n",
    "\n",
    "$$ (\\hat{\\mu}_{SMM},\\hat{\\sigma}_{SMM}) = (\\mu,\\sigma):\\quad \\min_{\\mu,\\sigma} e(\\tilde{scores}_i,scores_i|\\mu,\\sigma)^T \\, W \\, e(\\tilde{scores}_i,scores_i|\\mu,\\sigma) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's define a criterion function that takes as inputs the parameters and the estimator for the weighting matrix $\\hat{W}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def err_vec(data_vals, sim_vals, mu, sigma, cut_lb, cut_ub, simple):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the vector of moment errors (in percent\n",
    "    deviation from the data moment vector) for SMM.\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    data_vals = (N,) vector, test scores data\n",
    "    sim_vals  = (N, S) matrix, S simulations of test scores data\n",
    "    mu        = scalar, mean of the nontruncated normal distribution\n",
    "                from which the truncated normal is derived\n",
    "    sigma     = scalar > 0, standard deviation of the nontruncated\n",
    "                normal distribution from which the truncated normal is\n",
    "                derived\n",
    "    cut_lb    = scalar or string, ='None' if no lower bound cutoff is\n",
    "                given, otherwise is scalar lower bound value of\n",
    "                distribution. Values below this cutoff have zero\n",
    "                probability\n",
    "    cut_ub    = scalar or string, ='None' if no upper bound cutoff is\n",
    "                given, otherwise is scalar lower bound value of\n",
    "                distribution. Values below this cutoff have zero\n",
    "                probability\n",
    "    simple    = boolean, =True if errors are simple difference, =False\n",
    "                if errors are percent deviation from data moments\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        data_moments()\n",
    "        \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    mean_data  = scalar, mean value of data\n",
    "    var_data   = scalar > 0, variance of data\n",
    "    moms_data  = (2, 1) matrix, column vector of two data moments\n",
    "    mean_model = scalar, estimated mean value from model\n",
    "    var_model  = scalar > 0, estimated variance from model\n",
    "    moms_model = (2, 1) matrix, column vector of two model moments\n",
    "    err_vec    = (2, 1) matrix, column vector of two moment error\n",
    "                 functions\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: err_vec\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    mean_data, var_data = data_moments(data_vals)\n",
    "    moms_data = np.array([[mean_data], [var_data]])\n",
    "    mean_sim, var_sim = data_moments(sim_vals)\n",
    "    mean_model = mean_sim.mean()\n",
    "    var_model = var_sim.mean()\n",
    "    moms_model = np.array([[mean_model], [var_model]])\n",
    "    if simple:\n",
    "        err_vec = moms_model - moms_data\n",
    "    else:\n",
    "        err_vec = (moms_model - moms_data) / moms_data\n",
    "    \n",
    "    return err_vec\n",
    "\n",
    "\n",
    "def criterion(params, *args):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the SMM weighted sum of squared moment errors\n",
    "    criterion function value given parameter values and an estimate of\n",
    "    the weighting matrix.\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    params    = (2,) vector, ([mu, sigma])\n",
    "    mu        = scalar, mean of the normally distributed random variable\n",
    "    sigma     = scalar > 0, standard deviation of the normally\n",
    "                distributed random variable\n",
    "    args      = length 5 tuple,\n",
    "                (xvals, unif_vals, cut_lb, cut_ub, W_hat)\n",
    "    xvals     = (N,) vector, values of the truncated normally\n",
    "                distributed random variable\n",
    "    unif_vals = (N, S) matrix, matrix of draws from U(0,1) distribution.\n",
    "                This fixes the seed of the draws for the simulations\n",
    "    cut_lb    = scalar or string, ='None' if no lower bound cutoff is\n",
    "                given, otherwise is scalar lower bound value of\n",
    "                distribution. Values below this cutoff have zero\n",
    "                probability\n",
    "    cut_ub    = scalar or string, ='None' if no upper bound cutoff is\n",
    "                given, otherwise is scalar lower bound value of\n",
    "                distribution. Values below this cutoff have zero\n",
    "                probability\n",
    "    W_hat     = (R, R) matrix, estimate of optimal weighting matrix\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        norm_pdf()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    err        = (2, 1) matrix, column vector of two moment error\n",
    "                 functions\n",
    "    crit_val   = scalar > 0, GMM criterion function value\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: crit_val\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    mu, sigma = params\n",
    "    xvals, unif_vals, cut_lb, cut_ub, W_hat = args\n",
    "    sim_vals = trunc_norm_draws(unif_vals, mu, sigma, cut_lb, cut_ub)\n",
    "    err = err_vec(xvals, sim_vals, mu, sigma, cut_lb, cut_ub,\n",
    "                  simple=False)\n",
    "    crit_val = np.dot(np.dot(err.T, W_hat), err) \n",
    "    \n",
    "    return crit_val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_test = 400\n",
    "sig_test = 70\n",
    "cut_lb = 0.0\n",
    "cut_ub = 450.0\n",
    "# unif_vals_2 = sts.uniform.rvs(0, 1, size=(N, S))\n",
    "sim_vals = trunc_norm_draws(unif_vals_2, mu_test, sig_test, cut_lb, cut_ub)\n",
    "mean_sim, var_sim = data_moments(sim_vals)\n",
    "mean_mod = mean_sim.mean()\n",
    "var_mod = var_sim.mean()\n",
    "print(mean_mod, var_mod)\n",
    "err_vec(pts, sim_vals, mu_test, sig_test, cut_lb, cut_ub, simple=False)\n",
    "crit_test = criterion(np.array([mu_test, sig_test]), pts, unif_vals_2,\n",
    "                      0.0, 450.0, np.eye(2))\n",
    "print(crit_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can perform the SMM estimation. Let's start with the identity matrix as our estimate for the optimal weighting matrix $W = I$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_init_1 = 300\n",
    "sig_init_1 = 30\n",
    "params_init_1 = np.array([mu_init_1, sig_init_1])\n",
    "W_hat1_1 = np.eye(2)\n",
    "smm_args1_1 = (pts, unif_vals_2, cut_lb, cut_ub, W_hat1_1)\n",
    "results1_1 = opt.minimize(criterion, params_init_1, args=(smm_args1_1),\n",
    "                          method='L-BFGS-B',\n",
    "                          bounds=((1e-10, None), (1e-10, None)))\n",
    "mu_SMM1_1, sig_SMM1_1 = results1_1.x\n",
    "print('mu_SMM1_1=', mu_SMM1_1, ' sig_SMM1_1=', sig_SMM1_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_data, var_data = data_moments(pts)\n",
    "print('Data mean of scores =', mean_data, ', Data variance of scores =', var_data)\n",
    "sim_vals_1 = trunc_norm_draws(unif_vals_2, mu_SMM1_1, sig_SMM1_1, cut_lb, cut_ub)\n",
    "mean_sim_1, var_sim_1 = data_moments(sim_vals_1)\n",
    "mean_model_1 = mean_sim_1.mean()\n",
    "var_model_1 = var_sim_1.mean()\n",
    "err_1 = err_vec(pts, sim_vals_1, mu_SMM1_1, sig_SMM1_1, cut_lb, cut_ub,\n",
    "                False).reshape(2,)\n",
    "\n",
    "print('Model mean 1 =', mean_model_1, ', Model variance 1 =', var_model_1)\n",
    "print('Error vector 1 =', err_1)\n",
    "print(results1_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the PDF implied by these results against the histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Plot the histogram of the data\n",
    "count, bins, ignored = plt.hist(pts, 30, normed=True,\n",
    "                                edgecolor='black', linewidth=1.2)\n",
    "plt.title('Econ 381 scores: 2011-2012', fontsize=20)\n",
    "plt.xlabel('Total points')\n",
    "plt.ylabel('Percent of scores')\n",
    "plt.xlim([0, 550])  # This gives the xmin and xmax to be plotted\"\n",
    "\n",
    "# Plot the estimated SMM PDF\n",
    "dist_pts = np.linspace(0, 450, 500)\n",
    "plt.plot(dist_pts, trunc_norm_pdf(dist_pts, mu_SMM1_1, sig_SMM1_1, 0.0, 450.0),\n",
    "         linewidth=2, color='k', label='1: $\\mu_{SMM1}$,$\\sigma_{SMM1}$')\n",
    "plt.legend(loc='upper left')\n",
    "bins"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That looks just like the ML estimate from the previous notebook. Let's see what the criterion function looks like for different values of $\\mu$ and $\\sigma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Note that this will take a few minutes because the intgr.quad() commands\n",
    "# are a little slow\n",
    "mu_vals = np.linspace(60, 700, 50)\n",
    "sig_vals = np.linspace(20, 250, 50)\n",
    "# mu_vals = np.linspace(600, 610, 50)\n",
    "# sig_vals = np.linspace(190, 196, 50)\n",
    "crit_vals = np.zeros((50, 50))\n",
    "for mu_ind in range(50):\n",
    "    for sig_ind in range(50):\n",
    "        crit_vals[mu_ind, sig_ind] = \\\n",
    "            criterion(np.array([mu_vals[mu_ind], sig_vals[sig_ind]]),\n",
    "                      pts, unif_vals_2, cut_lb, cut_ub, W_hat1_1)\n",
    "\n",
    "mu_mesh, sig_mesh = np.meshgrid(mu_vals, sig_vals)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "ax.plot_surface(sig_mesh, mu_mesh, crit_vals, rstride=8,\n",
    "                cstride=1, cmap=cmap1)\n",
    "ax.set_title('Criterion function for values of mu and sigma')\n",
    "ax.set_xlabel(r'$\\sigma$')\n",
    "ax.set_ylabel(r'$\\mu$')\n",
    "ax.set_zlabel(r'Crit. func.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to the ML problem, it looks like the criterion function is roughly equal for a specific portion increase of $\\mu$ and $\\sigma$ together. That is, with these two moments probably have a correspondence of values of $\\mu$ and $\\sigma$ that give roughly the same criterion function value. This issue has two possible solutions.\n",
    "\n",
    "1. Maybe we need the two-step variance covariance estimator to calculate a \"more\" optimal weighting matrix $W$.\n",
    "2. Maybe our two moments aren't very good moments for fitting the data.\n",
    "\n",
    "Let's first try the two-step weighting matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "err2_1 = err_vec(pts, sim_vals_1, mu_SMM1_1, sig_SMM1_1, cut_lb, cut_ub, False)\n",
    "VCV2_1 = np.dot(err2_1, err2_1.T) / pts.shape[0]\n",
    "print(VCV2_1)\n",
    "W_hat2_1 = lin.pinv(VCV2_1)  # Use the pseudo-inverse calculated by SVD because VCV2 is ill-conditioned\n",
    "print(W_hat2_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params_init2_1 = np.array([mu_SMM1_1, sig_SMM1_1])\n",
    "# W_hat3 = np.array([[1. / VCV2[0, 0], 0.], [0., 1. / VCV2[1, 1]]])\n",
    "smm_args2_1 = (pts, unif_vals_2, cut_lb, cut_ub, W_hat2_1)\n",
    "results2_1 = opt.minimize(criterion, params_init2_1, args=(smm_args2_1),\n",
    "                         method='L-BFGS-B',\n",
    "                         bounds=((1e-10, None), (1e-10, None)))\n",
    "mu_SMM2_1, sig_SMM2_1 = results2_1.x\n",
    "print('mu_SMM2_1=', mu_SMM2_1, ' sig_SMM2_1=', sig_SMM2_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using a better weighting matrix didn't improve our estimates or fit very much. This means that we did not choose good moments for fitting the data. Let's try some different moments. How about four moments to match.\n",
    "\n",
    "1. The percent of observations greater than 430 (between 430 and 450)\n",
    "2. The percent of observations between 320 and 430\n",
    "3. The percent of observations between 220 and 320\n",
    "4. The percent of observations less than 220 (between 0 and 220)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def data_moments4(xvals):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the four data moments for SMM\n",
    "    (binpct_1, binpct_2, binpct_3, binpct_4) from both the actual data\n",
    "    and from the simulated data.\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    xvals = (N, S) matrix, (N,) vector, or scalar in (cut_lb, cut_ub),\n",
    "            test scores data, either real world or simulated. Real world\n",
    "            data will come in the form (N,). Simulated data comes in the\n",
    "            form (N,) or (N, S).\n",
    "\n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION: None\n",
    "\n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    bpct_1 = scalar in [0, 1] or (S,) vector, percent of observations\n",
    "             0 <= x < 220\n",
    "    bpct_2 = scalar in [0, 1] or (S,) vector, percent of observations\n",
    "             220 <= x < 320\n",
    "    bpct_3 = scalar in [0, 1] or (S,) vector, percent of observations\n",
    "             320 <= x < 430\n",
    "    bpct_4 = scalar in [0, 1] or (S,) vector, percent of observations\n",
    "             430 <= x <= 450\n",
    "\n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "\n",
    "    RETURNS: bpct_1, bpct_2, bpct_3, bpct_4\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    if xvals.ndim == 1:\n",
    "        bpct_1 = (xvals < 220).sum() / xvals.shape[0]\n",
    "        bpct_2 = ((xvals >=220) & (xvals < 320)).sum() / xvals.shape[0]\n",
    "        bpct_3 = ((xvals >=320) & (xvals < 430)).sum() / xvals.shape[0]\n",
    "        bpct_4 = (xvals >= 430).sum() / xvals.shape[0]\n",
    "    if xvals.ndim == 2:\n",
    "        bpct_1 = (xvals < 220).sum(axis=0) / xvals.shape[0]\n",
    "        bpct_2 = (((xvals >=220) & (xvals < 320)).sum(axis=0) /\n",
    "                    xvals.shape[0])\n",
    "        bpct_3 = (((xvals >=320) & (xvals < 430)).sum(axis=0) /\n",
    "                    xvals.shape[0])\n",
    "        bpct_4 = (xvals >= 430).sum(axis=0) / xvals.shape[0]\n",
    "\n",
    "    return bpct_1, bpct_2, bpct_3, bpct_4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def err_vec4(data_vals, sim_vals, mu, sigma, cut_lb, cut_ub, simple):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the vector of moment errors (in percent\n",
    "    deviation from the data moment vector) for SMM.\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    data_vals = (N,) vector, test scores data\n",
    "    sim_vals  = (N, S) matrix, S simulations of test scores data\n",
    "    mu        = scalar, mean of the nontruncated normal distribution\n",
    "                from which the truncated normal is derived\n",
    "    sigma     = scalar > 0, standard deviation of the nontruncated\n",
    "                normal distribution from which the truncated normal is\n",
    "                derived\n",
    "    cut_lb    = scalar or string, ='None' if no lower bound cutoff is\n",
    "                given, otherwise is scalar lower bound value of\n",
    "                distribution. Values below this cutoff have zero\n",
    "                probability\n",
    "    cut_ub    = scalar or string, ='None' if no upper bound cutoff is\n",
    "                given, otherwise is scalar lower bound value of\n",
    "                distribution. Values below this cutoff have zero\n",
    "                probability\n",
    "    simple    = boolean, =True if errors are simple difference, =False\n",
    "                if errors are percent deviation from data moments\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        data_moments4()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    mean_data  = scalar, mean value of data\n",
    "    var_data   = scalar > 0, variance of data\n",
    "    moms_data  = (4, 1) matrix, column vector of two data moments\n",
    "    mean_model = scalar, mean value from model\n",
    "    var_model  = scalar > 0, variance from model\n",
    "    moms_model = (2, 1) matrix, column vector of two model moments\n",
    "    err_vec    = (2, 1) matrix, column vector of two moment error\n",
    "                 functions\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: err_vec\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    bpct_1_dat, bpct_2_dat, bpct_3_dat, bpct_4_dat = \\\n",
    "        data_moments4(data_vals)\n",
    "    moms_data = np.array([[bpct_1_dat], [bpct_2_dat], [bpct_3_dat],\n",
    "                          [bpct_4_dat]])\n",
    "    bpct_1_sim, bpct_2_sim, bpct_3_sim, bpct_4_sim = \\\n",
    "        data_moments4(sim_vals)\n",
    "    bpct_1_mod = bpct_1_sim.mean()\n",
    "    bpct_2_mod = bpct_2_sim.mean()\n",
    "    bpct_3_mod = bpct_3_sim.mean()\n",
    "    bpct_4_mod = bpct_4_sim.mean()\n",
    "    moms_model = np.array([[bpct_1_mod], [bpct_2_mod], [bpct_3_mod],\n",
    "                          [bpct_4_mod]])\n",
    "    if simple:\n",
    "        err_vec = moms_model - moms_data\n",
    "    else:\n",
    "        err_vec = (moms_model - moms_data) / moms_data\n",
    "    \n",
    "    return err_vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def criterion4(params, *args):\n",
    "    '''\n",
    "    --------------------------------------------------------------------\n",
    "    This function computes the SMM weighted sum of squared moment errors\n",
    "    criterion function value given parameter values and an estimate of\n",
    "    the weighting matrix.\n",
    "    --------------------------------------------------------------------\n",
    "    INPUTS:\n",
    "    params    = (2,) vector, ([mu, sigma])\n",
    "    mu        = scalar, mean of the normally distributed random variable\n",
    "    sigma     = scalar > 0, standard deviation of the normally\n",
    "                distributed random variable\n",
    "    args      = length 5 tuple,\n",
    "                (xvals, unif_vals, cut_lb, cut_ub, W_hat)\n",
    "    xvals     = (N,) vector, values of the truncated normally\n",
    "                distributed random variable\n",
    "    unif_vals = (N, S) matrix, matrix of draws from U(0,1) distribution.\n",
    "                This fixes the seed of the draws for the simulations\n",
    "    cut_lb    = scalar or string, ='None' if no lower bound cutoff is\n",
    "                given, otherwise is scalar lower bound value of\n",
    "                distribution. Values below this cutoff have zero\n",
    "                probability\n",
    "    cut_ub    = scalar or string, ='None' if no upper bound cutoff is\n",
    "                given, otherwise is scalar lower bound value of\n",
    "                distribution. Values below this cutoff have zero\n",
    "                probability\n",
    "    W_hat     = (R, R) matrix, estimate of optimal weighting matrix\n",
    "    \n",
    "    OTHER FUNCTIONS AND FILES CALLED BY THIS FUNCTION:\n",
    "        norm_pdf()\n",
    "    \n",
    "    OBJECTS CREATED WITHIN FUNCTION:\n",
    "    err        = (2, 1) matrix, column vector of two moment error\n",
    "                 functions\n",
    "    crit_val   = scalar > 0, GMM criterion function value\n",
    "    \n",
    "    FILES CREATED BY THIS FUNCTION: None\n",
    "    \n",
    "    RETURNS: crit_val\n",
    "    --------------------------------------------------------------------\n",
    "    '''\n",
    "    mu, sigma = params\n",
    "    xvals, unif_vals, cut_lb, cut_ub, W_hat = args\n",
    "    sim_vals = trunc_norm_draws(unif_vals, mu, sigma, cut_lb, cut_ub)\n",
    "    err = err_vec4(xvals, sim_vals, mu, sigma, cut_lb, cut_ub,\n",
    "                   simple=False)\n",
    "    \n",
    "    # These next two lines only get uncommented to diagnose a problem\n",
    "#     print(mu, sigma)\n",
    "#     print(err.reshape(4,))\n",
    "    \n",
    "    crit_val = np.dot(np.dot(err.T, W_hat), err)\n",
    "    \n",
    "    return crit_val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will execute the SMM minimization problem, but a strange issue will arise. And the issue has to do with the minimizer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_init1_2 = 300\n",
    "sig_init1_2 = 30\n",
    "params_init1_2 = np.array([mu_init1_2, sig_init1_2])\n",
    "W_hat1_2 = np.eye(4)\n",
    "# W_hat[1, 1] = 2.0\n",
    "# W_hat[2, 2] = 2.0\n",
    "smm_args1_2 = (pts, unif_vals_2, cut_lb, cut_ub, W_hat1_2)\n",
    "results1_2 = opt.minimize(criterion4, params_init1_2, args=(smm_args1_2),\n",
    "                          method='L-BFGS-B',\n",
    "                          bounds=((1e-10, None), (1e-10, None)))\n",
    "mu_SMM1_2, sig_SMM1_2 = results1_2.x\n",
    "print('mu_SMM1_2=', mu_SMM1_2, ' sig_SMM1_2', sig_SMM1_2)\n",
    "print(results1_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the optimization problem only did three function evaluations, and it decided that the parameter values that minimized the criterion function are the initial values. Something is wrong.\n",
    "\n",
    "To see what is happening in the minimizer, let's insert a line in the `criterion4()` function that prints out the values of $\\mu$ and $\\sigma$ for each function evaluation in the minimizer as well as the error vector associated with each guess of $\\mu$ and $\\sigma$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the three function evaluations are for guesses of $\\mu$ and $\\sigma$ of:\n",
    "\n",
    "* Guess 1: $\\mu$=`mu_init` and $\\sigma$=`sig_init`\n",
    "* Guess 2: $\\mu$=`mu_init + 0.00000001` and $\\sigma$=`sig_init`\n",
    "* Guess 3: $\\mu$=`mu_init` and $\\sigma$=`sig_init + 0.00000001`\n",
    "\n",
    "This is the `L-BFGS-B` method's way of computing the Jacobian or slope (gradient) matrix of the criterion function by finite difference. However, the epsilon of `0.00000001` seems to be too small. We can set this step size to be bigger by using the `minimize()` function's `options={}` argument.\n",
    "\n",
    "The `options={}` argument in the `minimize()` function is a dictionary of solver options available to each particular method. In our case, we want to look at the `options={}` arguments for the [`L-BFGS-B` method](https://docs.scipy.org/doc/scipy-0.18.1/reference/optimize.minimize-lbfgsb.html#optimize-minimize-lbfgsb) of the `scipy.minimize()` function. Looking at this documentation, we find that we can set the `eps` option to something other than its default which is `options={'eps': 1e-08}`. In our case, we want to set that epsilon value used the finite differnce estimation of the Jacobian to be something bigger. Our means and variances seem to be in the 100's, so let's see if we get a solution setting the epsilon equal to 1.0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results1_2 = opt.minimize(criterion4, params_init1_2, args=(smm_args1_2),\n",
    "                          method='L-BFGS-B',\n",
    "                          bounds=((1e-10, None), (1e-10, None)),\n",
    "                          options={'eps': 1.0})\n",
    "mu_SMM1_2, sig_SMM1_2 = results1_2.x\n",
    "print('mu_SMM1_2=', mu_SMM1_2, ' sig_SMM1_2', sig_SMM1_2)\n",
    "print(results1_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the PDF implied by these results against the histogram of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the histogram of the data\n",
    "count, bins, ignored = plt.hist(pts, 30, normed=True,\n",
    "                                edgecolor='black', linewidth=1.2)\n",
    "plt.title('Econ 381 scores: 2011-2012', fontsize=20)\n",
    "plt.xlabel('Total points')\n",
    "plt.ylabel('Percent of scores')\n",
    "plt.xlim([0, 550])  # This gives the xmin and xmax to be plotted\"\n",
    "\n",
    "# Plot the estimated SMM PDF\n",
    "dist_pts = np.linspace(cut_lb, cut_ub, 500)\n",
    "plt.plot(dist_pts, trunc_norm_pdf(dist_pts, mu_SMM1_2, sig_SMM1_2,\n",
    "         cut_lb, cut_ub), linewidth=2, color='k',\n",
    "         label='1: $\\mu_{SMM1}$,$\\sigma_{SMM1}$')\n",
    "plt.legend(loc='upper left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's print the data moments and the model moments as well as the error vector evaluated at the SMM estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bpct_1_data, bpct_2_data, bpct_3_data, bpct_4_data = data_moments4(pts)\n",
    "print(bpct_1_data, bpct_2_data, bpct_3_data, bpct_4_data)\n",
    "sim_vals1_2 = trunc_norm_draws(unif_vals_2, mu_SMM1_2, sig_SMM1_2, cut_lb, cut_ub)\n",
    "bpct_1_sim1_2, bpct_2_sim1_2, bpct_3_sim1_2, bpct_4_sim1_2 = data_moments4(sim_vals1_2)\n",
    "bpct_1_model_2 = bpct_1_sim1_2.mean()\n",
    "bpct_2_model_2 = bpct_2_sim1_2.mean()\n",
    "bpct_3_model_2 = bpct_3_sim1_2.mean()\n",
    "bpct_4_model_2 = bpct_4_sim1_2.mean()\n",
    "print(bpct_1_model_2, bpct_2_model_2, bpct_3_model_2, bpct_4_model_2)\n",
    "err1_2 = err_vec4(pts, sim_vals1_2, mu_SMM1_2, sig_SMM1_2, cut_lb, cut_ub,\n",
    "                  False)\n",
    "print('Error vector =', err1_2.reshape(4,))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how much things change if we use the two-step estimator for the optimal weighting matrix $W$ instead of the identity matrix.\n",
    "\n",
    "TWO QUESTIONs: In the weighting matrix below $\\hat{W}_{2step}$, which moment will receive the most weight? Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VCV2_2 = np.dot(err1_2, err1_2.T) / pts.shape[0]\n",
    "print(VCV2_2)\n",
    "\n",
    "# Use the pseudo-inverse calculated by SVD because VCV2 is ill-conditioned\n",
    "W_hat2_2 = lin.pinv(VCV2_2) \n",
    "print(W_hat2_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# params_init2_2 = np.array([mu_SMM1_2, sig_SMM1_2])\n",
    "params_init2_2 = np.array([400, 70])\n",
    "# W_hat[1, 1] = 2.0\n",
    "# W_hat[2, 2] = 2.0\n",
    "smm_args2_2 = (pts, unif_vals_2, cut_lb, cut_ub, W_hat2_2)\n",
    "results2_3 = opt.minimize(criterion4, params_init2_2, args=(smm_args2_2),\n",
    "                          method='SLSQP',\n",
    "                          bounds=((1e-10, None), (1e-10, None)),\n",
    "                          options={'eps': 1.0})\n",
    "mu_SMM2_3, sig_SMM2_3 = results2_3.x\n",
    "print('mu_SMM2_3=', mu_SMM2_3, ' sig_SMM2_3', sig_SMM2_3)\n",
    "print(results2_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, using the two-step estimator for the optimal weighting matrix $\\hat{W}_{2step}$ makes a difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the histogram of the data\n",
    "count, bins, ignored = plt.hist(pts, 30, normed=True,\n",
    "                                edgecolor='black', linewidth=1.2)\n",
    "plt.title('Econ 381 scores: 2011-2012', fontsize=20)\n",
    "plt.xlabel('Total points')\n",
    "plt.ylabel('Percent of scores')\n",
    "plt.xlim([0, 550])  # This gives the xmin and xmax to be plotted\"\n",
    "\n",
    "# Plot the estimated SMM PDF\n",
    "dist_pts = np.linspace(cut_lb, cut_ub, 500)\n",
    "plt.plot(dist_pts, trunc_norm_pdf(dist_pts, mu_SMM1_2, sig_SMM1_2,\n",
    "         cut_lb, cut_ub), linewidth=2, color='r',\n",
    "         label='1: $\\mu_{SMM2}$,$\\sigma_{SMM2}$')\n",
    "\n",
    "# Plot the estimated SMM PDF\n",
    "dist_pts = np.linspace(cut_lb, cut_ub, 500)\n",
    "plt.plot(dist_pts, trunc_norm_pdf(dist_pts, mu_SMM2_3, sig_SMM2_3,\n",
    "         cut_lb, cut_ub), linewidth=2, color='k',\n",
    "         label='1: $\\mu_{SMM3}$,$\\sigma_{SMM3}$')\n",
    "plt.legend(loc='upper left')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(bpct_1_data, bpct_2_data, bpct_3_data, bpct_4_data)\n",
    "sim_vals2_2 = trunc_norm_draws(unif_vals_2, mu_SMM2_2, sig_SMM2_2, cut_lb, cut_ub)\n",
    "bpct_1_sim2_2, bpct_2_sim2_2, bpct_3_sim2_2, bpct_4_sim2_2 = data_moments4(sim_vals2_2)\n",
    "bpct_1_mode2_2 = bpct_1_sim2_2.mean()\n",
    "bpct_2_mode2_2 = bpct_2_sim2_2.mean()\n",
    "bpct_3_mode2_2 = bpct_3_sim2_2.mean()\n",
    "bpct_4_mode2_2 = bpct_4_sim2_2.mean()\n",
    "print(bpct_1_mode2_2, bpct_2_mode2_2, bpct_3_mode2_2, bpct_4_mode2_2)\n",
    "err2_2 = err_vec4(pts, sim_vals1_2, mu_SMM2_2, sig_SMM2_2, cut_lb, cut_ub,\n",
    "                  False)\n",
    "print('Error vector =', err2_2.reshape(4,))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The criterion function for different values of $\\mu$ and $\\sigma$ here has a clear minimum in a certain area. But it also has some really interesting nonlinearities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Note that this will take a few minutes because the intgr.quad() commands\n",
    "# are a little slow\n",
    "mu_vals = np.linspace(370, 380, 50)\n",
    "sig_vals = np.linspace(35, 55, 50)\n",
    "# mu_vals = np.linspace(350, 370, 50)\n",
    "# sig_vals = np.linspace(85, 98, 50)\n",
    "crit_vals4 = np.zeros((50, 50))\n",
    "for mu_ind in range(50):\n",
    "    for sig_ind in range(50):\n",
    "        crit_vals4[mu_ind, sig_ind] = \\\n",
    "            criterion4(np.array([mu_vals[mu_ind], sig_vals[sig_ind]]),\n",
    "                      pts, unif_vals_2, cut_lb, cut_ub, W_hat2_2)\n",
    "\n",
    "mu_mesh, sig_mesh = np.meshgrid(mu_vals, sig_vals)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "ax.plot_surface(sig_mesh, mu_mesh, crit_vals4, rstride=8,\n",
    "                cstride=1, cmap=cmap1)\n",
    "ax.set_title('Criterion function for values of mu and sigma')\n",
    "ax.set_xlabel(r'$\\sigma$')\n",
    "ax.set_ylabel(r'$\\mu$')\n",
    "ax.set_zlabel(r'Crit. func.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2. Brock and Mirman (1972) estimation by SMM\n",
    "Give Brock and Mirman (1972) example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Identification\n",
    "An issue that we saw in the examples from the previous section is that there is some science as well as some art in choosing moments to identify the parameters in an SMM estimation as well as in GMM. Suppose the parameter vector $\\theta$ has $K$ elements, or rather, $K$ parameters to be estimated. In order to estimate $\\theta$ by GMM, you must have at least as many moments as parameters to estimate $R\\geq K$. If you have exactly as many moments as parameters to be estimated $R=K$, the model is said to be *exactly identified*. If you have more moments than parameters to be estimated $R>K$, the model is said to be *overidentified*. If you have fewer moments than parameters to be estimated $R<K$, the model is said to be *underidentified*. There are good reasons to overidentify $R>K$ the model in SMM estimation as we saw in the previous example. The main reason is that not all moments are orthogonal. That is, some moments convey roughly the same information about the data and, therefore, do not separately identify any extra parameters. So a good SMM model often is overidentified $R>K$.\n",
    "\n",
    "One last point about MM regards moment selection and verification of results. The real world has an infinite supply of potential moments that describe some part of the data. Choosing moments to estimate parameters by SMM requires understanding of the model, intuition about its connections to the real world, and artistry. A good SMM estimation will include moments that have some relation to or story about their connection to particular parameters of the model to be estimated. In addition, a good verification of a SMM estimation is to take some moment from the data that was not used in the estimation and see how well the corresponding moment from the estimated model matches that *outside moment*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Indirect Inference\n",
    "Indirect inference is a particular application of SMM with some specific characteristics. As moments to match it uses parameters of an auxiliary model that can be estimated both on the real-world data and on the simulated data. Smith (2008) gives a great summary of the topic with some examples. See also Gourieroux and Monfort (1996, ch. 4) for a textbook treatment of the topic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1. Restatement of the general SMM estimation problem\n",
    "Define a model or data generating process (DGP) as a system of equations,\n",
    "\n",
    "$$ F(x_t,z_t|\\theta)=0 $$\n",
    "\n",
    "which are functions of a vector of endogenous variables $x_t$, exogenous variables $z_t$, and parameters $\\theta$. In the general simulated method of moments (SMM), one would choose data moments $m(x_t,z_t)$ that are just statistics of the data and model moments $\\hat{m}(\\tilde{x}_t,\\tilde{z}_t|\\theta)$ that are averages of the same data moments calculated on simulated samples of the data. The SMM estimator is to choose the parameter vector $\\hat{\\theta}_{SMM}$ to minimize some distance of the model moments from the data moments.\n",
    "\n",
    "$$ \\hat{\\theta}_{SMM}=\\theta:\\quad \\min_{\\theta} ||\\hat{m}(\\tilde{x}_t,\\tilde{z}_t|\\theta) - m(x_t,z_t)|| $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2. Indirect inference estimation problem\n",
    "Indirect inference is to change the model moments from being stastics that are calculated directly from the simulated data to being statistics that are calculated indirectly from the simulated data. These indirect inference model moments are moments calculated by an auxiliary model. In practice, these auxiliary model moments are often estimated regression parameters.\n",
    "\n",
    "Let an auxiliary model be defined as $g(x_t,z_t|\\phi)=0$. The parameters of the auxiliary model will be the moments we use to identify the model parameters $\\theta$. Suppose that the model parameter vector $\\theta$ has $K$ elements. Then the auxiliary model parameter vector $\\phi$ must have $R$ elements such that $R\\geq K$. This is the typical identification restriction that the number of model moments must be at least as many as the number of model parameters being estimated.\n",
    "\n",
    "When the auxiliary model is run on real-world data $g(x_t,z_t|\\phi)=0$, the resulting values of the auxiliary model parameters are the data moments $\\hat{\\phi}(x_t,z_t)$. Note that these data moments $\\hat{\\phi}$ have a hat on them to represent that these moments are usually estimated in some way. When the auxiliary model is run on the $s$th simulation of the data given model parameters $g(\\tilde{x}_{s,t},\\tilde{z}_{s,t}|\\phi)=0$, the auxiliary model parameters are the $s$th estimate of the model moments $\\hat{\\phi}_s(\\tilde{x}_{s,t},\\tilde{z}_{s,t}|\\theta)$. The model moments are then the average of these auxiliary model parameter estimates across the simulations.\n",
    "\n",
    "$$ \\hat{\\phi}(\\tilde{x}_{t},\\tilde{z}_{t}|\\theta) = \\frac{1}{S}\\sum_{s=1}^S \\hat{\\phi}_s(\\tilde{x}_{s,t},\\tilde{z}_{s,t}|\\theta) $$\n",
    "\n",
    "The indirect inference estimation method is simply to choose a model parameter vector $\\theta$ that minimizes some distance metric between the model moments $\\hat{\\phi}(\\tilde{x}_{t},\\tilde{z}_{t}|\\theta)$ and the data moments $\\hat{\\phi}(x_t,z_t)$.\n",
    "\n",
    "$$ \\hat{\\theta}_{SMM}=\\theta:\\quad \\min_{\\theta} ||\\hat{\\phi}(\\tilde{x}_{t},\\tilde{z}_{t}|\\theta) - \\hat{\\phi}(x_t,z_t)|| $$\n",
    " \n",
    "In most examples of indirect, the data moments and model moments are some regression of endogenous variables on exogenous variables. In the univariate case, it is usually linear regression. In the multivariate case, it is usually a vector autoregression (VAR). But most examples are reduced form parameter estimation exercises. Other examples are probit, logit, and two-stage IV regressions. The key is that these statistics be computationally tractable and have convenient or accurate data availability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.3. Hypothesis testing with indirect inference\n",
    "* Wald test\n",
    "* likelihood ratio test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. References\n",
    "* Adda, Jerome and Russell Cooper, *Dynamic Economics: Quantitative Methods and Applications*, MIT Press (2003).\n",
    "* Brock, William A. and Leonard J. Mirman, \"Optimal Economic Growth and Uncertainty: The Discounted Case,\" *Journal of Economic Theory*, 4:3, pp. 479-513 (June 1972).\n",
    "* Davidson, Russell and James G. MacKinnon, *Econometric Theory and Methods*, Oxford University Press (2004).\n",
    "* Duffie, Darrell and Kenneth J. Singleton, \"Simulated Moment Estimation of Markov Models of Asset Prices\", *Econometrica*, 61:4, pp. 929-952 (July 1993).\n",
    "* Gourieroux, Christian and Alain Monfort, *Simulation-based Econometric Methods*, Oxford University Press (1996).\n",
    "* Laroque, G. and B. Salanie, \"Simulation Based Estimation Models with Lagged Latent Variables\", *Journal of Applied Econometrics*, 8:Supplement, pp. 119-133 (December 1993).\n",
    "* Lee, Bong-Soo and Beth Fisher Ingram, \"Simulation Estimation of Time Series Models\", *Journal of Econometrics*, 47:2-3, pp. 197-205 (February 1991).\n",
    "* McFadden, Daniel, \"A Method of Simulated Moments for Estimation of Discrete Response Models without Numerical Integration,\" *Econometrica*, 57:5, pp. 995-1026 (September 1989).\n",
    "* Newey, Whitney K. and Kenneth D. West, \"A Simple, Positive, Semi-definite, Heteroskedasticy and Autocorrelation Consistent Covariance Matrix,\" *Econometrica*, 55:3, pp. 703-708 (May 1987).\n",
    "* Smith, Anthony A. Jr., \"[Indirect Inference](http://www.econ.yale.edu/smith/palgrave7.pdf),\" *New Palgrave Dictionary of Economics*, 2nd edition, (2008)."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}