Ch16_MLE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maximum Likelihood Estimation\n",
    "\n",
    "This notebook illustrates maximum likelihood estimation and how to calculate different standard errors (from the information matrix, the gradients and the \"sandwich\" approach).\n",
    "\n",
    "The application is very simple: estimating the mean and variance of a random variable. Some of the subsequent chapters (notebooks) work with more complicated models, for instance, GARCH models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Packages and Extra Functions\n",
    "\n",
    "The notebook first implements MLE step-by-step. At the end it also presents the `MLE()` function from the (local) `FinEcmt_MLEGMM` module that wraps those calculations.\n",
    "\n",
    "For the numerical optimization we use the [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl) package and for calculating derivatives the [FiniteDiff.jl](https://github.com/JuliaDiff/FiniteDiff.jl) package. ([ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) is an alternative.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "MyModulePath = joinpath(pwd(),\"jlFiles\")\n",
    "!in(MyModulePath,LOAD_PATH) && push!(LOAD_PATH,MyModulePath);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "using LinearAlgebra, DelimitedFiles, Statistics, Optim\n",
    "using FinEcmt_OLS\n",
    "using FinEcmt_MLEGMM: MLE     #load only the MLE() function\n",
    "\n",
    "#loading and renaming some functions from FiniteDiff (to get shorter names)\n",
    "using FiniteDiff: finite_difference_hessian as hessian, finite_difference_jacobian as jacobian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "xx  = readdlm(\"Data/FFdSizePs.csv\",',',skipstart=1)\n",
    "x   = xx[:,2]                 #returns for the portfolio of the smallest firms\n",
    "xx  = nothing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Traditional Estimates\n",
    "\n",
    "of the mean $\\mu$ and the variance $\\sigma^2$.\n",
    "\n",
    "The standard errors of the mean and standard deviation are from traditional textbook formulas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mTraditional estimates and their std:\u001b[22m\u001b[39m\n",
      "\n",
      "    estimate       std\n",
      "μ      0.042     0.010\n",
      "σ²     0.840     0.013\n",
      "\n"
     ]
    }
   ],
   "source": [
    "T = length(x)\n",
    "\n",
    "(μ_trad,σ²_trad) = (mean(x),var(x,corrected=false))    #corrected=false gives 1/T formula, not 1/(T-1)\n",
    "\n",
    "std_trad = sqrt.([σ²_trad,2*σ²_trad^2]/T)   #standard errors, textbook formulas\n",
    "\n",
    "printblue(\"Traditional estimates and their std:\\n\")\n",
    "xx = [[μ_trad,σ²_trad] std_trad]\n",
    "printmat(xx;colNames=[\"estimate\",\"std\"],rowNames=[\"μ\",\"σ²\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Point Estimates from ML\n",
    "\n",
    "The next few cells define a log likelihood function and estimate the coefficients by maximizing it.\n",
    "\n",
    "## The (log) Likelihood Function for Estimating the Parameters of a N(,)\n",
    "\n",
    "### A Remark on the Code\n",
    "\n",
    "- `(μ,σ²) = par` splits up the vector `par` into two numbers (for the mean and variance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    NormalLL(par,x)\n",
    "\n",
    "Calculate log likelihood for a `N(μ,σ²)` distribution.\n",
    "\n",
    "`par = [μ,σ²]` is a vector with the parameters , `x` is a vector with data\n",
    "\"\"\"\n",
    "function NormalLL(par,x)\n",
    "    (μ,σ²) = par\n",
    "    LLt    = -(1/2)*log(2*pi) - (1/2)*log(σ²) .- (1/2)*abs2.(x.-μ)/σ²  #vector, all x[t]\n",
    "    return LLt\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Try the Likelihood Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "log likelihood value at par0: -11155.385\n"
     ]
    }
   ],
   "source": [
    "par0 = [0.0,1.0]           #initial parameter guess of [μ,σ²]\n",
    "\n",
    "LLt = NormalLL(par0,x)     #just trying the log likelihood fn\n",
    "\n",
    "printlnPs(\"log likelihood value at par0: \",sum(LLt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimize the Likelihood Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "log-likelihood at point estimate (compare with the value above): -11088.409\n",
      "\n",
      "\u001b[34m\u001b[1mParameter estimates:\u001b[22m\u001b[39m\n",
      "\n",
      "    traditional          MLE\n",
      "μ         0.042        0.042\n",
      "σ²        0.840        0.840\n",
      "\n"
     ]
    }
   ],
   "source": [
    "Sol = optimize(par->-sum(NormalLL(par,x)),par0)  #minimize -sum(LLt)\n",
    "parHat = Optim.minimizer(Sol)                    #the optimal solution\n",
    "\n",
    "printlnPs(\"log-likelihood at point estimate (compare with the value above): \",-Optim.minimum(Sol))\n",
    "\n",
    "printblue(\"\\nParameter estimates:\\n\")\n",
    "xx = [[μ_trad,σ²_trad] parHat]\n",
    "printmat(xx;colNames=[\"traditional\",\"MLE\"],rowNames=[\"μ\",\"σ²\"],width=13)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Standard Errors I: Information Matrix \n",
    "\n",
    "If the likelihood function is correctly specified, then MLE is typically asymptotically normally distributed as\n",
    "\n",
    "$\n",
    "\\sqrt{T}(\\hat{\\theta}-\\theta)  \\rightarrow^{d}N(0,V) \\: \\text{, where } \\: V=I(\\theta)^{-1}\\text{ with }\n",
    "$\n",
    "\n",
    "$\n",
    "I(\\theta) =-\\text{E}\\frac{\\partial^{2}\\ln L_t}{\\partial\\theta\\partial\\theta^{\\prime}}\n",
    "$\n",
    "\n",
    "where $I(\\theta)$ is the information matrix (*not* an identity matrix) and $\\ln L_t$  is the contribution of period $t$ to the log likelihood function.\n",
    "\n",
    "### A Remark on the Code\n",
    "\n",
    "The code below calculates numerical derivatives. It does so by noticing that $\n",
    "\\text{E}\\frac{\\partial^{2}\\ln L_t}{\\partial\\theta\\partial\\theta^{\\prime}} = \n",
    "\\frac{\\partial^{2}\\text{E}\\ln L_t}{\\partial\\theta\\partial\\theta^{\\prime}} \n",
    "$,\n",
    "so we can differentiate the mean (across data points) log likelihood value. (Clearly, we could also calculate $T$ different $2\\times 2$ matrices and then average, but that would be slower.)\n",
    "\n",
    "- `-hessian(par->mean(NormalLL(par,x)),parHat)` calculates the negative of the 2nd derivatives of the average log-likelihood function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mstandard errors:\u001b[22m\u001b[39m\n",
      "\n",
      "         traditional     MLE (InfoMat)\n",
      "μ              0.010             0.010\n",
      "σ²             0.013             0.013\n",
      "\n"
     ]
    }
   ],
   "source": [
    "Ia = -hessian(par->mean(NormalLL(par,x)),parHat)  #(-1) * 2nd derivatives of mean(LLt)\n",
    "\n",
    "Ia       = (Ia+Ia')/2         #to guarantee symmetry, fixes possible rounding errors\n",
    "vcv      = inv(Ia)/T\n",
    "std_hess = sqrt.(diag(vcv))\n",
    "\n",
    "printblue(\"standard errors:\\n\")\n",
    "xx = [std_trad std_hess]\n",
    "printmat(xx;colNames=[\"traditional\",\"MLE (InfoMat)\"],rowNames=[\"μ\",\"σ²\"],width=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Standard Errors II: Gradients and Sandwich\n",
    "\n",
    "An alternative way of calculating the information matrix is to use the outer product of the gradients:\n",
    "\n",
    "$\n",
    "J(\\theta)=\\text{E}\\left[  \\frac{\\partial\\ln L_t}{\\partial\\theta\n",
    "}\\frac{\\partial\\ln L_t}{\\partial\\theta^{\\prime}}\\right]\n",
    "$\n",
    "\n",
    "This would coincide with $I(\\theta)$ as defined above, if the model is correctly specified, but the sample estimates can differ (also when the model is correctly specified).\n",
    "\n",
    "### A Remark on the Code\n",
    "\n",
    "The code below fills row $t$ of a $T\\times 2$ matrix (called `δL`) with \n",
    "$\n",
    "\\frac{\\partial\\ln L_t}{\\partial\\theta}.\n",
    "$\n",
    "For each $t$, the outer product is a $2\\times2$ matrix, and then we average across $t$. This is done by calculating \n",
    "`J = δL'δL/T`,\n",
    "which is the same as `for t = 1:T; J = J + δL[t,:]*δL[t,:]'/T;end;`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Std from Gradients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mstandard errors:\u001b[22m\u001b[39m\n",
      "\n",
      "         traditional     MLE (InfoMat)   MLE (gradients)\n",
      "μ              0.010             0.010             0.010\n",
      "σ²             0.013             0.013             0.005\n",
      "\n"
     ]
    }
   ],
   "source": [
    "δL = jacobian(par->NormalLL(par,x),parHat)   #Tx2\n",
    "J  = δL'δL/T               #2xT * Tx2\n",
    "\n",
    "vcv       = inv(J)/T\n",
    "std_grad  = sqrt.(diag(vcv))                          #std from gradients\n",
    "\n",
    "printblue(\"standard errors:\\n\")\n",
    "xx = [std_trad std_hess std_grad]\n",
    "printmat(xx;colNames=[\"traditional\",\"MLE (InfoMat)\",\"MLE (gradients)\"],rowNames=[\"μ\",\"σ²\"],width=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Std from Sandwich\n",
    "\n",
    "We could also use the \"sandwich\" estimator\n",
    "\n",
    "$\n",
    "V=I(\\theta)^{-1}J(\\theta)I(\\theta)^{-1}.\n",
    "$\n",
    "\n",
    "When the model (likelihood function) is misspecified, then the three variance-covariance matrices may differ, and the sandwich approach is often the most robust."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mstandard errors:\u001b[22m\u001b[39m\n",
      "\n",
      "         traditional     MLE (InfoMat)   MLE (gradients)    MLE (sandwich)\n",
      "μ              0.010             0.010             0.010             0.010\n",
      "σ²             0.013             0.013             0.005             0.036\n",
      "\n"
     ]
    }
   ],
   "source": [
    "vcv       = inv(Ia) * J * inv(Ia)/T\n",
    "std_sandw = sqrt.(diag(vcv))                          #std from sandwich\n",
    "\n",
    "printblue(\"standard errors:\\n\")\n",
    "xx = [std_trad std_hess std_grad std_sandw]\n",
    "printmat(xx,colNames=[\"traditional\",\"MLE (InfoMat)\",\"MLE (gradients)\",\"MLE (sandwich)\"],rowNames=[\"μ\",\"σ²\"],width=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try this: replace the data series `x` with simulated data from a $N()$ distribution. Then, do the different standard errors get closer to each other?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A Function for MLE (extra)\n",
    "\n",
    "The next cell uses a function which combines the computations of several of the cells above. \n",
    "\n",
    "It requires a function for the log-likehood function as written above, that is, taking `(par,x)` as inputs and generating a T-vector `LLt` as output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "```\n",
       "MLE(LLtFun,par0,x,lower,upper)\n",
       "```\n",
       "\n",
       "Calculate ML point estimates of K parameters and three different types of standard errors: from the Information matrix, from the gradients and the sandwich approach.\n",
       "\n",
       "### Input\n",
       "\n",
       "  * `LLtFun::Function`:  name of log-likelihood function\n",
       "  * `par0::Vector`:      K-vector, starting guess of the parameters\n",
       "  * `x::VecOrMat`:       vector or matrix with data, unpacked (if necessary) inside LLtFun\n",
       "  * `lower::Vector`:     lower bounds on the parameters, nothing or fill(-Inf,K) if no bounds\n",
       "  * `upper::Vector`:     upper bounds on the parameters, nothing or fill(Inf,K) if no bounds\n",
       "\n",
       "### Requires\n",
       "\n",
       "  * `using FiniteDiff: finite_difference_hessian as hessian, finite_difference_jacobian as jacobian`\n",
       "\n",
       "### Notice\n",
       "\n",
       "The `LLtFun` should take `(par,x)` as inputs and generate a T-vector `LLt` as output.\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@doc2 MLE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#using CodeTracking\n",
    "#println(@code_string MLE(NormalLL,[1],[1]))    #print the source code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mpoint estimates and standard errors:\u001b[22m\u001b[39m\n",
      "\n",
      "            estimate     std (InfoMat)   std (gradients)    std (sandwich)\n",
      "μ              0.042             0.010             0.010             0.010\n",
      "σ²             0.840             0.013             0.005             0.036\n",
      "\n"
     ]
    }
   ],
   "source": [
    "(parHat,std_hess,std_grad,std_sandw) = MLE(NormalLL,par0,x)\n",
    "\n",
    "printblue(\"point estimates and standard errors:\\n\")\n",
    "xx = [parHat std_hess std_grad std_sandw]\n",
    "printmat(xx;colNames=[\"estimate\",\"std (InfoMat)\",\"std (gradients)\",\"std (sandwich)\"],rowNames=[\"μ\",\"σ²\"],width=18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 1.10.2",
   "language": "julia",
   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}