{
 "cells": [
  {
   "outputs": [],
   "cell_type": "markdown",
   "source": [
    "# Maximum Likelihood Estimation: The Normal Linear Model"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "markdown",
   "source": [
    "The following tutorial will introduce maximum likelihood estimation\n",
    "in Julia for the normal linear model.\n",
    "\n",
    "The normal linear model (sometimes referred to as the OLS model) is\n",
    "the workhorse of regression modeling and is utilized across a number\n",
    "of diverse fields.  In this tutorial, we will utilize simulated data\n",
    "to demonstrate how Julia can be used to recover the parameters of\n",
    "interest.\n",
    "\n",
    "The first order of business is to use the `Optim` package\n",
    "and also include the `NLSolversBase` routine:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Optim, NLSolversBase\n",
    "srand(0);                            # Fix random seed generator for reproducibility"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "outputs": [],
   "cell_type": "markdown",
   "source": [
    "The first item that needs to be addressed is the data generating process or DGP.\n",
    "The following code will produce data from a nomral linear model:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "n = 500                             # Number of observations\n",
    "nvar = 2                            # Number of variables\n",
    "β = ones(nvar) * 3.0                # True coefficients\n",
    "x = [ones(n) randn(n, nvar - 1)]    # X matrix of explanatory variables plus constant\n",
    "ε = randn(n) * 0.5                  # Error variance\n",
    "y = x * β + ε;                      # Generate Data"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "outputs": [],
   "cell_type": "markdown",
   "source": [
    "In the above example, we have 500 observations, 2 explanatory\n",
    "variables plus an intercept, an error variance equal to 0.5,\n",
    "coefficients equal to 3.0, and all of these are subject to change by\n",
    "the user. Since we know the true value of these parameters, we\n",
    "should obtain these values when we maximize the likelihood function.\n",
    "\n",
    "The next step in our tutorial is to define a Julia function for the\n",
    "likelihood function. The following function defines the likelihood\n",
    "function for the normal linear model:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Log_Likelihood (generic function with 1 method)"
     },
     "metadata": {},
     "execution_count": 3
    }
   ],
   "cell_type": "code",
   "source": [
    "function Log_Likelihood(X, Y, β, log_σ)\n",
    "    σ = exp(log_σ)\n",
    "    llike = -n/2*log(2π) - n/2* log(σ^2) - (sum((Y - X * β).^2) / (2σ^2))\n",
    "    llike = -llike\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "outputs": [],
   "cell_type": "markdown",
   "source": [
    "The log likelihood function accepts 4 inputs: the matrix of\n",
    "explanatory variables (X), the dependent variable (Y), the β's, and\n",
    "the error varicance. Note that we exponentiate the error variance in\n",
    "the second line of the code because the error variance cannot be\n",
    "negative and we want to avoid this situation when maximizing the\n",
    "likelihood.\n",
    "\n",
    "The next step in our tutorial is to optimize our function. We first\n",
    "use the `TwiceDifferentiable` command in order to obtain the Hessian\n",
    "matrix later on, which will be used to help form t-statistics:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Results of Optimization Algorithm\n * Algorithm: Newton's Method\n * Starting Point: [1.0,1.0,1.0]\n * Minimizer: [3.002788633849947,2.964549617572727, ...]\n * Minimum: 3.851229e+02\n * Iterations: 7\n * Convergence: true\n   * |x - x'| ≤ 0.0e+00: false \n     |x - x'| = 2.62e-08 \n   * |f(x) - f(x')| ≤ 0.0e+00 |f(x)|: false\n     |f(x) - f(x')| = 1.33e-15 |f(x)|\n   * |g(x)| ≤ 1.0e-08: true \n     |g(x)| = 3.41e-13 \n   * Stopped by an increasing objective: false\n   * Reached Maximum Number of Iterations: false\n * Objective Calls: 31\n * Gradient Calls: 31\n * Hessian Calls: 7"
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "cell_type": "code",
   "source": [
    "func = TwiceDifferentiable(vars -> Log_Likelihood(x, y, vars[1:nvar], vars[nvar + 1]),\n",
    "                           ones(nvar+1); autodiff=:forward);\n",
    "# The above statment accepts 4 inputs: the x matrix, the dependent\n",
    "# variable y, and a vector of β's and the error variance.  The\n",
    "# `vars[1:nvar]` is how we pass the vector of β's and the `vars[nvar +\n",
    "# 1]` is how we pass the error variance. You can think of this as a\n",
    "# vector of parameters with the first 2 being β's and the last one is\n",
    "# the error variance.\n",
    "#\n",
    "# The `ones(nvar+1)` are the starting values for the parameters and\n",
    "# the `autodiff=:forward` command performs forward mode automatic\n",
    "# differentiation.\n",
    "#\n",
    "# The actual optimization of the likelihood function is accomplished\n",
    "# with the following command:\n",
    "opt = optimize(func, ones(nvar+1))"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "outputs": [],
   "cell_type": "markdown",
   "source": [
    "The first input to the command is the function we wish to optimize\n",
    "and the second input are the starting values.\n",
    "\n",
    "After a brief period of time, you should see output of the\n",
    "optimization routine, with the parameter  estimates being very close\n",
    "to our simulated values.\n",
    "\n",
    "The optimization routine stores several quantities and we can obtain\n",
    "the maximim likelihood estimates with the following command:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "3-element Array{Float64,1}:\n  3.00279 \n  2.96455 \n -0.648693"
     },
     "metadata": {},
     "execution_count": 5
    }
   ],
   "cell_type": "code",
   "source": [
    "parameters = Optim.minimizer(opt)"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "outputs": [],
   "cell_type": "markdown",
   "source": [
    "!!! Note\n",
    "    Fieldnames for all of the quantities can be obtained with the following command:\n",
    "    fieldnames(opt)\n",
    "\n",
    "Since we paramaterized our likelihood to use the exponentiated\n",
    "value, we need to exponentiate it to get back to our original log\n",
    "scale:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "0.5227286513837306"
     },
     "metadata": {},
     "execution_count": 6
    }
   ],
   "cell_type": "code",
   "source": [
    "parameters[nvar+1] = exp(parameters[nvar+1])"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "outputs": [],
   "cell_type": "markdown",
   "source": [
    "In order to obtain the correct Hessian matrix, we have to \"push\" the\n",
    "actual parameter values that maximizes the likelihood function since\n",
    "the `TwiceDifferentiable` command uses the next to last values to\n",
    "calculate the Hessian:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "3×3 Array{Float64,2}:\n 175.766        -12.0877        5.67464e-14\n -12.0877       182.437         5.88344e-15\n   5.67464e-14    5.88344e-15  96.0542     "
     },
     "metadata": {},
     "execution_count": 7
    }
   ],
   "cell_type": "code",
   "source": [
    "numerical_hessian = hessian!(func,parameters)"
   ],
   "metadata": {},
   "execution_count": 7
  },
  {
   "outputs": [],
   "cell_type": "markdown",
   "source": [
    "We can now invert our Hessian matrix to obtain the variance-covariance matrix:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "3×3 Array{Float64,2}:\n  0.00571544    0.000378687  -3.39973e-18\n  0.000378687   0.00550643   -5.60994e-19\n -3.39973e-18  -5.60994e-19   0.0104108  "
     },
     "metadata": {},
     "execution_count": 8
    }
   ],
   "cell_type": "code",
   "source": [
    "var_cov_matrix = inv(numerical_hessian)"
   ],
   "metadata": {},
   "execution_count": 8
  },
  {
   "outputs": [],
   "cell_type": "markdown",
   "source": [
    "In this example, we are only interested in the statistical\n",
    "significance of the coefficient estimates so we obtain those with\n",
    "the following command:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "2-element Array{Float64,1}:\n 3.00279\n 2.96455"
     },
     "metadata": {},
     "execution_count": 9
    }
   ],
   "cell_type": "code",
   "source": [
    "β = parameters[1:nvar]"
   ],
   "metadata": {},
   "execution_count": 9
  },
  {
   "outputs": [],
   "cell_type": "markdown",
   "source": [
    "We now need to obtain those elements of the variance-covariance\n",
    "matrix needed to obtain our t-statistics, and we can do this with\n",
    "the following commands:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "2-element Array{Float64,1}:\n 0.00571544\n 0.00550643"
     },
     "metadata": {},
     "execution_count": 10
    }
   ],
   "cell_type": "code",
   "source": [
    "temp = diag(var_cov_matrix)\n",
    "temp1 = temp[1:nvar]"
   ],
   "metadata": {},
   "execution_count": 10
  },
  {
   "outputs": [],
   "cell_type": "markdown",
   "source": [
    "The t-statistics are formed by dividing element-by-element the\n",
    "coefficients by their standard errors, or the square root of the\n",
    "diagonal elements of the variance-covariance matrix:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "2-element Array{Float64,1}:\n 39.7191\n 39.9506"
     },
     "metadata": {},
     "execution_count": 11
    }
   ],
   "cell_type": "code",
   "source": [
    "t_stats = β./sqrt.(temp1)"
   ],
   "metadata": {},
   "execution_count": 11
  },
  {
   "outputs": [],
   "cell_type": "markdown",
   "source": [
    "From here, one may examine other statistics of interest using the\n",
    "output from the optimization routine."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "markdown",
   "source": [
    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.4"
  },
  "kernelspec": {
   "name": "julia-0.6",
   "display_name": "Julia 0.6.4",
   "language": "julia"
  }
 },
 "nbformat": 4
}