Week6_Empirical_Exercise_Matlab.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\textbf{Empirical Exercise-Week 6}$  Page 402, Problem 8.1\n",
    "Heteroskedasticity and Labor Market Discrimination "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "clear all;\n",
    "Tab1 = readtable('cps5.csv');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ans =\n",
      "\n",
      "  1x23 table\n",
      "\n",
      "    age    asian    black    divorced    educ    exper    faminc    female    hrswork    insure    married    mcaid    mcare    metro    midwest    nchild    northeast    single    south    union    wage     west    white\n",
      "    ___    _____    _____    ________    ____    _____    ______    ______    _______    ______    _______    _____    _____    _____    _______    ______    _________    ______    _____    _____    _____    ____    _____\n",
      "\n",
      "    45       0        0         0         13      26      39200       1         38         1          1         0        0        0         0         0           1          0         0        1      14.38     0        1  \n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "Tab1(1,:)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "wage = table2array(Tab1(:,21)); educ = table2array(Tab1(:,5));\n",
    "exper = table2array(Tab1(:,6)); female = table2array(Tab1(:,8));\n",
    "aframer = table2array(Tab1(:,3)); metro = table2array(Tab1(:,14));\n",
    "south = table2array(Tab1(:,19)); midwest = table2array(Tab1(:,15));\n",
    "west = table2array(Tab1(:,22));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "y = log(wage);\n",
    "x = [ones(length(y),1) educ exper exper.^2 female aframer metro south midwest west];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "BINT =\n",
      "\n",
      "    1.1384    1.2643\n",
      "    0.0978    0.1047\n",
      "    0.0271    0.0322\n",
      "   -0.0005   -0.0004\n",
      "   -0.1842   -0.1468\n",
      "   -0.1447   -0.0783\n",
      "    0.0949    0.1431\n",
      "   -0.0723   -0.0192\n",
      "   -0.0916   -0.0363\n",
      "   -0.0348    0.0216\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "[B BINT R STATS] = regress(y,x);\n",
    "BINT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You see that the coefficient for female = 1 has a depressive effect on wages.  Given equal education, experience, location, women earn less in the United States.\n",
    "Lets do a Goldfeld-Quandt test on equalty of variances for the residuals of wages for men and women in the US workforce. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "ehat = R;\n",
    "index1 = female == 1;\n",
    "ehat_women = ehat(index1,:);\n",
    "index2 = female == 0;\n",
    "ehat_men = ehat(index2,:);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ans =\n",
      "\n",
      "    0.2113    0.2215\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "Var_women = var(ehat_women);  Var_men = var(ehat_men);\n",
    "[Var_women Var_men]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "F_stat =\n",
      "\n",
      "    1.1829\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "dgf_women = length(ehat_women)-10; dgf_men = length(ehat_men)-10;\n",
    "F_stat =  (Var_women/dgf_women)/(Var_men/dgf_men)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ans =\n",
      "\n",
      "   2.3886e-09\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "Fprob = fcdf(F_stat, dgf_women, dgf_men);\n",
    "1 - Fprob"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can reject the hypothesis of equaltiy of variances between the men and women wage-equation residuals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "New we do the N $R^2$ test for Heteroskeasticity.  We regress the squared residuals on three variables, on Metro, Female, Black and their cross terms. If the null is true, no heteroskedascity, the term N $R^2$ has a $\\chi^{2}$ distribution with K*(K+3)/2 degress of freedom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "yy = ehat .^2;\n",
    "x1 = ones(length(yy),1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "xfa = female .* aframer;\n",
    "xfm = female .* metro;\n",
    "xma = metro .* aframer;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "xx = [x1 female aframer metro xfa xfm xma];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "[B2 BINT2 R2 RINT2 STATS2] = regress(yy,xx);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Chistat =\n",
      "\n",
      "   32.2982\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "Chistat =  STATS2(1)* length(ehat);\n",
    "dgf2 = 7*10/2;\n",
    "Chistat\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ans =\n",
      "\n",
      "    0.5992\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "Chip = chi2cdf(Chistat, dgf2);\n",
    "1 - Chip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we caonnot reject heteroskedasticty with these variables in this test.  Let repeat it with all of the explanatory variables, as your book suggests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ans =\n",
      "\n",
      "   4.7186e-04\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "[B4 BINT4 R4 RINT4 STATS4] = regress(yy,x);\n",
    "dgf4 =  10 * 13/2;\n",
    "Chistat =  length(yy)*STATS4(1);\n",
    "Chip = chi2cdf(Chistat, dgf4);\n",
    "1-Chip\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now reject homoskecasticty when we use all of the regressors.\n",
    "White test is similar, with quadratic expansion on the variables. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "xxx= [ones(length(yy),1) educ  educ.^2 exper exper.^2 educ.*exper ...\n",
    "female aframer metro south midwest west female.*aframer metro.*aframer ...\n",
    "metro.*south metro.*female];\n",
    "[B5 BINT5 R5 RINT5 STATS5] = regress(yy,xxx);\n",
    "dgf5 = 15;\n",
    "White = length(yy) * STATS5(1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ans =\n",
      "\n",
      "     0\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "White_test = chi2cdf(White,dgf5);\n",
    "1 - White_test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can reject homoskedasticity with the White test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " HAC Heteroscedasticity and autocorrelation consistent covariance estimators\n",
      " \n",
      "  Syntax:\n",
      " \n",
      "    [EstCoeffCov,se,coeff] = hac(X,y)\n",
      "    [EstCoeffCov,se,coeff] = hac(Tbl)\n",
      "    [EstCoeffCov,se,coeff] = hac(Mdl)\n",
      "    [EstCoeffCov,se,coeff] = hac(...,param,val,...)\n",
      " \n",
      "  Description:\n",
      " \n",
      "    HAC computes robust covariance estimates for ordinary least squares\n",
      "    (OLS) coefficient estimates of multiple linear regression models\n",
      " \n",
      "    y = X*coeff + e\n",
      " \n",
      "    under general forms of heteroscedasticity and autocorrelation in the\n",
      "    innovations process e. Estimates are of the form\n",
      "                \n",
      "    PhiHat = X'*OmegaHat*X\n",
      "    EstCoeffCov = c*inv(X'*X)*PhiHat*inv(X'*X)\n",
      " \n",
      "    where OmegaHat is an estimate of the innovations covariance and c is an\n",
      "    optional small-sample correction.\n",
      " \n",
      "  Input Arguments:  \n",
      " \n",
      "    X - numObs-by-numPreds matrix of predictor data for a multiple linear\n",
      "        regression model.\n",
      " \n",
      "    y - numObs-by-1 vector of response data for a multiple linear\n",
      "        regression model.\n",
      " \n",
      "    Tbl - numObs-by-numPreds+1 tabular array of data for a multiple linear\n",
      "        regression model, with predictor data X in the first numPreds\n",
      "        columns and response data y in the last column.\n",
      " \n",
      "    Mdl - Linear model of the form Mdl = LinearModel.fit(X,y).\n",
      " \n",
      "    Observations with missing (NaN) values in the predictors or the\n",
      "    response are removed from the data.\n",
      " \n",
      "  Optional Input Parameter Name/Value Pairs:\n",
      " \n",
      "    NAME        VALUE\n",
      " \n",
      "    'varNames'  String vector or cell vector of character vectors, of\n",
      "                length numPreds, to be used as variable names in displays\n",
      "                of the results. Names should include the intercept term\n",
      "                (for example, 'Const') and higher-order terms (for example,\n",
      "                'x1^2' or 'x1:x2'), if present. The default for matrix X is\n",
      "                {'x1','x2',...}. The default for tabular array Tbl is\n",
      "                Tbl.Properties.VariableNames. The default for linear models\n",
      "                Mdl is Mdl.CoefficientNames.\n",
      " \n",
      "    'intercept' Logical value indicating whether or not to add an intercept\n",
      "                when fitting the model. The default is true. When the input\n",
      "                is already a fitted model, Mdl, 'intercept' is ignored.\n",
      " \n",
      "    'type' \t\tType of estimator. Values are 'HC' and 'HAC'. 'HC' computes\n",
      "                heteroscedasticity-consistent estimates, as described in\n",
      "                [9], [7], and [3]. 'HAC' computes heteroscedasticity-and-\n",
      "                autocorrelation-consistent estimates, as described in [8],\n",
      "                [5], [2], and [1]. The default is 'HAC'.\n",
      " \n",
      "    'weights' \tWeights used to compute PhiHat = X'*OmegaHat*X. Values are\n",
      "                a numerical vector w of length numObs or a string or\n",
      "                character vector indicating a method for computing w from\n",
      "                the data.\n",
      " \n",
      "                When 'type' is 'HC', OmegaHat = diag(w). Elements w(i)\n",
      "                estimate the innovations variance at each observation i.\n",
      "                Data-driven w are computed from model residuals u(i), their\n",
      "                leverages h(i), and their number of degrees of freedom dfe.\n",
      "                Values are:\n",
      " \n",
      "                o 'CLM' w(i) = (1/dfe)*sum(u(i)^2). This is the Classical\n",
      "                        Linear Model estimator, used in the absence of\n",
      "                        heteroscedasticity.\n",
      " \n",
      "                o 'HC0' w(i) = u(i)^2. This is the White estimator in [9].\n",
      "                        This is the default when 'type' is 'HC'.\n",
      " \n",
      "                o 'HC1' w(i) = (numObs/dfe)*u(i)^2. This is the HC1\n",
      "                        estimator in [7].\n",
      " \n",
      "                o 'HC2' w(i) = (u(i)^2)/(1-h(i)). This is the HC2\n",
      "                        estimator in [7].\n",
      " \n",
      "                o 'HC3' w(i) = (u(i)^2)/((1-h(i))^2). This is the HC3\n",
      "                        estimator in [7].\n",
      " \n",
      "                o 'HC4' w(i) = (u(i)^2)/((1-h(i))^d(i)), where d(i) is\n",
      "                        min(4,h(i)/mean(h)). This is the estimator in [3].\n",
      " \n",
      "                o w     Vector of user-specified weights, of length numObs.\n",
      "                        If data contain a missing value at observation i,\n",
      "                        w(i) is removed during estimation.\n",
      " \n",
      "                When 'type' is 'HAC', component products forming PhiHat,\n",
      "                X(i,:)'*(u(i)*u(j))*X(j,:), are weighted by w(l), where \n",
      "                l = abs(i-j) is the lag length between observations.\n",
      "                Elements w(l) represent the strength of autocorrelations at\n",
      "                each lag length. Data-driven w(l) = k(l/b) are computed by\n",
      "                a kernel density estimator k, with b specified by the\n",
      "                'bandwidth' parameter. Values are:\n",
      " \n",
      "                o 'TR' Truncated kernel:\n",
      " \n",
      "                              | 1  for abs(x) <= 1\n",
      "                       k(x) = |\n",
      "                              | 0  otherwise\n",
      " \n",
      "                       This produces the White estimator in [10, p. 152].\n",
      " \n",
      "                o 'BT' Bartlett kernel:\n",
      " \n",
      "                              | 1-abs(x)  for abs(x) <= 1\n",
      "                       k(x) = |\n",
      "                              | 0         otherwise\n",
      " \n",
      "                       This produces the Newey-West estimator in [8]. This\n",
      "                       is the default when 'type' is 'HAC'.\n",
      " \n",
      "                o 'PZ' Parzen kernel:\n",
      " \n",
      "                              | 1-6*x^2+6*abs(x)^3  for 0 <= abs(x) <= 1/2\n",
      "                       k(x) = | 2*(1-abs(x))^3      for 1/2 <= abs(x) <= 1\n",
      "                              | 0                   otherwise\n",
      " \n",
      "                       This produces the Gallant estimator in [5, p. 533].\n",
      " \n",
      "                o 'TH' Tukey-Hanning kernel:\n",
      " \n",
      "                              | (1+cos(pi*x))/2  for abs(x) <= 1\n",
      "                       k(x) = |\n",
      "                              | 0                otherwise\n",
      " \n",
      "                       This estimator was introduced by Andrews in [1].\n",
      " \n",
      "                o 'QS' Quadratic spectral kernel:\n",
      " \n",
      "                       k(x) = [25/(12*pi^2*x^2)]*...\n",
      "                              [sin(6*pi*x/5)/(6*pi*x/5)-cos(6*pi*x/5)]\n",
      " \n",
      "                       This estimator was introduced by Andrews in [1].\n",
      " \n",
      "                o w    Vector of user-specified weights, of length numObs.\n",
      "                       If data contain a missing value at observation i,\n",
      "                       w(i) is removed during estimation.\n",
      " \n",
      "    'bandwidth' Scalar bandwidth parameter b used by the kernel estimator\n",
      "                in 'weights' when 'type' is 'HAC'. If 'type' is 'HC',\n",
      "                'bandwidth' is ignored. The default is data-driven, as in\n",
      "                [1]. Values of 'AR1' or 'ARMA11' specify the model used for\n",
      "                data-driven bandwidth selection. The default model is\n",
      "                'AR1'. The method used to estimate the 'AR1' model can be\n",
      "                specified by values of 'AR1OLS', to use OLS, or 'AR1MLE',\n",
      "                to use maximum likelihood. The default method is 'AR1MLE'.\n",
      " \n",
      "    'smallT'    Logical flag indicating whether or not to apply the small-\n",
      "                sample correction c = T/dfe, where T is the effective\n",
      "                sample size and dfe is the number of degrees of freedom of\n",
      "                the model residuals. Values are true and false. The default\n",
      "                is false when 'type' is 'HC' and true when 'type' is 'HAC'.\n",
      " \n",
      "    'whiten'    Nonnegative integer specifying the lag order of a vector\n",
      "                autoregressive model used as a prewhitening filter when\n",
      "                'type' is 'HAC', as in [2]. If 'type' is 'HC', 'whiten' is\n",
      "                ignored. The default is 0, which bypasses the filter.\n",
      " \n",
      "    'display'   String or character vector to control the display of\n",
      "                results to the command window. The display shows outputs in\n",
      "                tabular form. Values are 'cov' to display only EstCoeffCov,\n",
      "                'full' to display all of coeff, se, and EstCoeffCov, or\n",
      "                'off' to turn off the display. The default is 'cov'.\n",
      "  \n",
      "  Output Arguments:\n",
      " \n",
      "    EstCoeffCov - numPreds-by-numPreds matrix of coefficient covariance\n",
      "                  estimates.\n",
      " \n",
      "    se - numPreds-by-1 vector of coefficient standard error estimates. This\n",
      "         is sqrt(diag(EstCoeffCov)).\n",
      " \n",
      "    coeff - numPreds-by-1 vector of OLS coefficient estimates.\n",
      " \n",
      "  Notes:\n",
      " \n",
      "    o Estimates formed by HAC are often called \"sandwich estimates,\" with\n",
      "      X'*OmegaHat*X the \"meat\" and inv(X'*X) the \"bread.\"\n",
      " \n",
      "    o For kernels with unit-interval support, the bandwidth parameter b is\n",
      "      often called the \"lag-truncation parameter,\" since w(l) = k(l/b) = 0\n",
      "      for lags l > b.\n",
      " \n",
      "    o HAC with default settings computes a Newey-West estimator using a\n",
      "      Bartlett kernel with a data-driven bandwidth. To obtain the standard\n",
      "      Newey-West estimate described in [8], set 'bandwidth' to maxLag+1,\n",
      "      with maxLag = floor(4*(T/100)^(2/9)).\n",
      " \n",
      "    o The original White HC estimator, HC0, is justified asymptotically.\n",
      "      HC1, HC2, HC3, and HC4 are meant to improve small-sample performance.\n",
      "      [6] and [3] recommend HC3 and HC4, respectively, in the presence of\n",
      "      influential observations.\n",
      " \n",
      "    o HAC estimators formed with the truncated kernel are not guaranteed to\n",
      "      be positive semidefinite in finite samples. [8] proposes the Bartlett\n",
      "      kernel as a remedy, but the resulting estimator is suboptimal in\n",
      "      terms of its rate of consistency. The quadratic spectral kernel\n",
      "      achieves an optimal rate of consistency.\n",
      " \n",
      "    o The default estimation method for HAC bandwidth selection, 'AR1MLE',\n",
      "      is generally more accurate, but slower, than the 'AR1' alternative,\n",
      "      'AR1OLS'. The 'ARMA11' model must be estimated by maximum likelihood.\n",
      "      Bandwidth-selection models may exhibit sensitivity to the relative\n",
      "      scale of the predictors in X.\n",
      " \n",
      "    o [2] recommends prewhitening for HAC estimators to reduce bias. The\n",
      "      procedure tends to increase estimator variance and mean-squared\n",
      "      error, but can improve confidence interval coverage probabilities and\n",
      "      reduce the over-rejection of t statistics.\n",
      " \n",
      "  Example:\n",
      " \n",
      "    % Compute OLS coefficients and Newey-West standard errors for a\n",
      "    % regression of nominal GNP (GNPN) on the consumer price index (CPI),\n",
      "    % real wages (WR), and the money stock (MS):\n",
      " \n",
      "    load Data_NelsonPlosser\n",
      "    Tbl = DataTable(:,[8,10,11,2]);\n",
      "    T = sum(~any(ismissing(Tbl),2));\n",
      "    maxLag = floor(4*(T/100)^(2/9));\n",
      "    EstCoeffCov = hac(Tbl,'bandwidth',maxLag+1,'display','full');\n",
      " \n",
      "  References:\n",
      "  \n",
      "    [1] Andrews, D. W. K. \"Heteroskedasticity and Autocorrelation\n",
      "        Consistent Covariance Matrix Estimation.\" Econometrica. v. 59,\n",
      "        1991, pp. 817-858.\n",
      " \n",
      "    [2] Andrews, D. W. K., and J. C. Monohan. \"An Improved\n",
      "        Heteroskedasticity and Autocorrelation Consistent Covariance Matrix\n",
      "        Estimator.\" Econometrica. v. 60, 1992, pp. 953-966.\n",
      " \n",
      "    [3] Cribari-Neto, F. \"Asymptotic Inference Under Heteroskedasticity of\n",
      "        Unknown Form.\" Computational Statistics & Data Analysis. v. 45,\n",
      "        2004, pp. 215-233.\n",
      " \n",
      "    [4] den Haan, W. J., and A. Levin. \"A Practitioner's Guide to Robust\n",
      "        Covariance Matrix Estimation.\" In Handbook of Statistics. Edited by\n",
      "        G. S. Maddala and C. R. Rao. Amsterdam: Elsevier, 1997.\n",
      " \n",
      "    [5] Gallant, A. R. Nonlinear Statistical Models. Hoboken, NJ: John\n",
      "        Wiley & Sons, Inc., 1987.\n",
      " \n",
      "    [6] Long, J. S., and L. H. Ervin. \"Using Heteroscedasticity-Consistent\n",
      "        Standard Errors in the Linear Regression Model.\" The American\n",
      "        Statistician. v. 54, 2000, pp. 217-224.\n",
      " \n",
      "    [7] MacKinnon, J. G., and H. White. \"Some Heteroskedasticity-Consistent\n",
      "        Covariance Matrix Estimators with Improved Finite Sample\n",
      "        Properties.\" Journal of Econometrics. v. 29, 1985, pp. 305-325.\n",
      " \n",
      "    [8] Newey, W. K., and K. D. West. \"A Simple, Positive-Definite,\n",
      "        Heteroskedasticity and Autocorrelation Consistent Covariance\n",
      "        Matrix.\" Econometrica. v. 55, 1987, pp. 703-708.\n",
      " \n",
      "    [9] White, H. \"A Heteroskedasticity-Consistent Covariance Matrix and a\n",
      "        Direct Test for Heteroskedasticity.\" Econometrica. v. 48, 1980, pp.\n",
      "        817-838.\n",
      " \n",
      "   [10] White, H. Asymptotic Theory for Econometricians. New York: Academic\n",
      "        Press, 1984.\n",
      " \n",
      "  See also FITLM, LSCOV.\n",
      "\n",
      "    Documentation for hac\n",
      "       doc hac\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help hac  % Use the heteroskeastic estimators from the econometrics toolbox;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Estimator type: HC\n",
      "Estimation method: HC0\n",
      "Effective sample size: 9799\n",
      "Small sample correction: off\n",
      "\n",
      "Coefficient Covariances:\n",
      "\n",
      "       |  Const      x1       x2       x3       x4       x5       x6       x7       x8       x9   \n",
      "--------------------------------------------------------------------------------------------------\n",
      " Const |  0.0011  -0.0001  -0.0000   0.0000  -0.0000  -0.0000  -0.0001  -0.0001  -0.0001  -0.0001 \n",
      " x1    | -0.0001   0.0000  -0.0000   0.0000  -0.0000   0.0000  -0.0000   0.0000   0.0000   0.0000 \n",
      " x2    | -0.0000  -0.0000   0.0000  -0.0000   0.0000   0.0000  -0.0000  -0.0000  -0.0000  -0.0000 \n",
      " x3    |  0.0000   0.0000  -0.0000   0.0000  -0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 \n",
      " x4    | -0.0000  -0.0000   0.0000  -0.0000   0.0001  -0.0000   0.0000  -0.0000  -0.0000   0.0000 \n",
      " x5    | -0.0000   0.0000   0.0000   0.0000  -0.0000   0.0003  -0.0000  -0.0000  -0.0000   0.0000 \n",
      " x6    | -0.0001  -0.0000  -0.0000   0.0000   0.0000  -0.0000   0.0001  -0.0000   0.0000  -0.0000 \n",
      " x7    | -0.0001   0.0000  -0.0000   0.0000  -0.0000  -0.0000  -0.0000   0.0002   0.0001   0.0001 \n",
      " x8    | -0.0001   0.0000  -0.0000   0.0000  -0.0000  -0.0000   0.0000   0.0001   0.0002   0.0001 \n",
      " x9    | -0.0001   0.0000  -0.0000   0.0000   0.0000   0.0000  -0.0000   0.0001   0.0001   0.0002 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "  [EstCoeffCov,se,coeff] = hac(x(:,2:end), y,'type','hc');\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ans =\n",
      "\n",
      "    1.2014    0.0328   36.6527\n",
      "    0.1012    0.0019   53.1431\n",
      "    0.0296    0.0013   22.5391\n",
      "   -0.0004    0.0000  -16.1615\n",
      "   -0.1655    0.0095  -17.4517\n",
      "   -0.1115    0.0161   -6.9333\n",
      "    0.1190    0.0116   10.2814\n",
      "   -0.0458    0.0139   -3.2931\n",
      "   -0.0639    0.0137   -4.6615\n",
      "   -0.0066    0.0145   -0.4529\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "[coeff se coeff./se] %  coefficients std errors, tstats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that being female or African American have negative effects on wages."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we turn to FGLS estimators, Feasible Generalized Least Squares Estimators\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " FGLS Feasible Generalized Least Squares\n",
      " \n",
      "  Syntax:\n",
      " \n",
      "    [coeff,se,EstCoeffCov] = fgls(X,y)\n",
      "    [coeff,se,EstCoeffCov] = fgls(Tbl)\n",
      "    [coeff,se,EstCoeffCov] = fgls(...,param,val,...)\n",
      "    [coeff,se,EstCoeffCov,iterPlots] = fgls(...)\n",
      "    [coeff,se,EstCoeffCov,iterPlots] = fgls(ax,...)\n",
      " \n",
      "  Description:\n",
      " \n",
      "    FGLS computes generalized least squares (GLS) estimates of coefficients\n",
      "    in multiple linear regression models\n",
      " \n",
      "    y = X*coeff + e\n",
      " \n",
      "    by first estimating the covariance of the innovations process e.\n",
      "    Coefficient estimates are of the form\n",
      "                \n",
      "    coeff = inv(X'*inv(OmegaHat)*X)*X'*inv(OmegaHat)*y\n",
      " \n",
      "    where OmegaHat is an innovations covariance estimate using a specified\n",
      "    model. The process can be iterated.\n",
      " \n",
      "  Input Arguments:\n",
      " \n",
      "    X - numObs-by-numPreds matrix of predictor data for a multiple linear                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         \n",
      "        regression model.\n",
      " \n",
      "    y - numObs-by-1 vector of response data for a multiple linear\n",
      "        regression model.\n",
      " \n",
      "    Tbl - numObs-by-numPreds+1 tabular array of data for a multiple linear\n",
      "        regression model, with predictor data X in the first numPreds\n",
      "        columns and response data y in the last column.\n",
      " \n",
      "    ax - Vector of axes objects in which to plot, of length equal to the\n",
      "        number of plots specified by the 'plot' parameter. If unspecified,\n",
      "        FGLS creates separate figures for each plot.\n",
      " \n",
      "    Observations with missing (NaN) values in the predictors or the\n",
      "    response are removed from the data.\n",
      " \n",
      "  Optional Input Parameter Name/Value Pairs:\n",
      " \n",
      "    NAME        VALUE\n",
      " \n",
      "    'varNames'\tCell vector of variable name strings of length numPreds to\n",
      "                be used in displays of the results. Names should include\n",
      "                the intercept term (for example, 'Const'), if present. The\n",
      "                default for matrix X is {'x1','x2',...}. The default for\n",
      "                tabular array Tbl is Tbl.Properties.VariableNames.\n",
      " \n",
      "    'intercept' Logical value indicating whether or not to add an intercept\n",
      "                when fitting the model. The default is true. The number of\n",
      "                model coefficients, numCoeffs, is numPreds plus\n",
      "                'intercept'.\n",
      " \n",
      "    'innovMdl'  Model used to estimate the innovations covariance,\n",
      "                OmegaHat.\n",
      " \n",
      "                Diagonal models (of heteroscedasticity) are of the form\n",
      "                OmegaHat = diag(w), where w is a vector of length numObs.\n",
      "                Elements w(i) estimate the innovations variance at each\n",
      "                observation i. Data-driven w are computed from model\n",
      "                residuals u(i), their leverages h(i), and their number of\n",
      "                degrees of freedom dfe. Values are:\n",
      " \n",
      "                o 'CLM' w(i) = (1/dfe)*sum(u(i)^2). This is the Classical\n",
      "                        Linear Model estimator, used in the absence of\n",
      "                        heteroscedasticity.\n",
      " \n",
      "                o 'HC0' w(i) = u(i)^2. This is the White estimator in [5].\n",
      " \n",
      "                o 'HC1' w(i) = (numObs/dfe)*u(i)^2. This is the HC1\n",
      "                        estimator in [4].\n",
      " \n",
      "                o 'HC2' w(i) = (u(i)^2)/(1-h(i)). This is the HC2\n",
      "                        estimator in [4].\n",
      " \n",
      "                o 'HC3' w(i) = (u(i)^2)/((1-h(i))^2). This is the HC3\n",
      "                        estimator in [4].\n",
      " \n",
      "                o 'HC4' w(i) = (u(i)^2)/((1-h(i))^d(i)), where d(i) is\n",
      "                        min(4,h(i)/mean(h)). This is the estimator in [1].\n",
      " \n",
      "                Full covariance models (of heteroscedasticity and\n",
      "                autocorrelation) must be specified with a feasible number\n",
      "                of parameters for estimation. Values are:\n",
      " \n",
      "                o 'AR' models innovations as an AR(p) process, with number\n",
      "                       of lags p specified by the 'arLags' parameter. This\n",
      "                       is the default.\n",
      " \n",
      "    'arLags'    When 'innovMdl' is 'AR', 'arLags' is the number of\n",
      "                autoregressive lags in the innovations model. The default\n",
      "                is 1. For diagonal models, 'arLags' is ignored.\n",
      " \n",
      "    'InnovCov0' Initial innovations covariance. This replaces the data-\n",
      "                driven estimate of OmegaHat in the first iteration of GLS.\n",
      "                For diagonal models, 'InnovCov0' is a length numObs vector\n",
      "                of innovations variances. For full models, 'InnovCov0' is a\n",
      "                numObs-by-numObs matrix of innovations covariances. The\n",
      "                default is empty, to indicate that OmegaHat is to be\n",
      "                estimated in the first iteration.\n",
      " \n",
      "    'numIter'   Number of iterations of the FGLS algorithm. Each iteration\n",
      "                estimates the innovations covariance OmegaHat from the\n",
      "                residual series according to 'innovMdl', then computes GLS\n",
      "                estimates of the model coefficients. The default is 1.\n",
      " \n",
      "    'resCond'   Logical value indicating whether or not to scale residuals\n",
      "                in each iteration, to improve the conditioning of OmegaHat.\n",
      "                The default is false.\n",
      " \n",
      "    'display'   String to control the display of results to the command\n",
      "                window. The display shows estimation results in tabular\n",
      "                form. Values are 'iter' to display each iteration, 'final'\n",
      "                to display only the final estimates, or 'off' to turn off\n",
      "                the display. The default is 'off'.\n",
      " \n",
      "    'plot'      String or string vector to control iteration plots. Values\n",
      "                are 'coeff' to plot coefficient estimates, 'se' to plot\n",
      "                coefficient standard error estimates, 'mse' to plot\n",
      "                residual mean-squared errors, 'all' to produce all plots,\n",
      "                or 'off' to turn off plotting. The default is 'off'.\n",
      " \n",
      "  Output Arguments:\n",
      " \n",
      " \tcoeff - numCoeffs-by-1 vector of FGLS coefficient estimates.\n",
      " \n",
      " \tse - numCoeffs-by-1 vector of coefficient standard error estimates.\n",
      " \n",
      "    EstCoeffCov - numCoeffs-by-numCoeffs matrix of coefficient covariance\n",
      "        estimates. se is sqrt(diag(EstCoeffCov)).\n",
      " \n",
      "    iterPlots - Structure array of handles to plotted graphics objects.\n",
      " \n",
      "  Notes:\n",
      " \n",
      "    o Standard GLS estimates are obtained with FGLS by setting InnovCov0\n",
      "      to a known covariance matrix, and numIter = 1. For weighted least\n",
      "      squares (WLS), set InnovCov0 to a vector of inverse weights (for\n",
      "      example, variance estimates).\n",
      " \n",
      "    o In the presence of nonspherical innovations, GLS produces efficient\n",
      "      estimates relative to OLS, and consistent coefficient covariances,\n",
      "      conditional on the innovations covariance. The degree to which these\n",
      "      properties are maintained by FGLS depends on the accuracy of both the\n",
      "      model and estimation of the innovations covariance.\n",
      " \n",
      "    o Asymptotic distributions of FGLS estimators are unchanged by repeated\n",
      "      iteration. However, iteration may change finite sample distributions.\n",
      " \n",
      "    o FGLS estimates have the form\n",
      " \n",
      "      coeff = inv(X'*inv(OmegaHat)*X)*X'*inv(OmegaHat)*y.\n",
      " \n",
      "      However, FGLS uses methods that are faster and more stable, and are\n",
      "      applicable to rank-deficient cases.\n",
      " \n",
      "    o In repeated iteration, scale differences in the residuals may produce\n",
      "      badly-conditioned OmegaHat, and numerical instability, in specific\n",
      "      models. Setting 'resCond' to true improves conditioning.\n",
      " \n",
      "    o General ARMA models of the innovations covariance can be converted to\n",
      "      the default 'AR' form using a lag operator representation and the\n",
      "      division operator in the LagOp class. You can also obtain the\n",
      "      coefficients of the AR representation of an ARMA model using ARMA2AR.\n",
      " \n",
      "    o Traditional FGLS methods, such as the Cochrane-Orcutt procedure, use\n",
      "      low order 'AR' models. These methods, however, estimate parameters in\n",
      "      the innovations covariance matrix using OLS, where FGLS uses maximum\n",
      "      likelihood estimation (MLE) [2].\n",
      " \n",
      "  Example:\n",
      " \n",
      "    % Compare OLS and FGLS estimates in a model with AR(2) errors:\n",
      " \n",
      "    numObs = 500;\n",
      "    beta = [1 2 3]';\n",
      "    phi =[0.2 0.1];\n",
      "    sigma2 = 0.5;\n",
      "    e = simulate(arima('Constant',0,'AR',phi,'Variance',sigma2),numObs);\n",
      "    X = randn(numObs,2);\n",
      "    y = [ones(numObs,1),X]*beta + e;\n",
      "  \n",
      "    fgls(X,y,'innovMdl','AR','arLags',2,'display','final');\n",
      " \n",
      "  References:\n",
      " \n",
      "    [1] Cribari-Neto, F. \"Asymptotic Inference Under Heteroskedasticity of\n",
      "        Unknown Form.\" Computational Statistics & Data Analysis. v. 45,\n",
      "        2004, pp. 215-233.\n",
      " \n",
      "    [2] Hamilton, J. D. Time Series Analysis. Princeton, NJ: Princeton\n",
      "        University Press, 1994.\n",
      " \n",
      "    [3] Judge, G. G., W. E. Griffiths, R. C. Hill, H. Lutkepohl, and T. C.\n",
      "        Lee. The Theory and Practice of Econometrics. New York, NY: John\n",
      "        Wiley & Sons, Inc., 1985.\n",
      " \n",
      "    [4] MacKinnon, J. G., and H. White. \"Some Heteroskedasticity-Consistent\n",
      "        Covariance Matrix Estimators with Improved Finite Sample\n",
      "        Properties.\" Journal of Econometrics. v. 29, 1985, pp. 305-325.\n",
      " \n",
      "    [5] White, H. \"A Heteroskedasticity-Consistent Covariance Matrix and a\n",
      "        Direct Test for Heteroskedasticity.\" Econometrica. v. 48, 1980, pp.\n",
      "        817-838.\n",
      " \n",
      "  See also LSCOV, HAC\n",
      "\n",
      "    Documentation for fgls\n",
      "       doc fgls\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help fgls\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    " [coeff1,se1,EstCoeffCov1] = fgls(x(:,2:end),y);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ans =\n",
      "\n",
      "    1.2235    0.0330   37.1014\n",
      "    0.0996    0.0018   55.7920\n",
      "    0.0296    0.0013   22.6168\n",
      "   -0.0004    0.0000  -16.7470\n",
      "   -0.1648    0.0093  -17.6437\n",
      "   -0.1136    0.0174   -6.5376\n",
      "    0.1193    0.0132    9.0308\n",
      "   -0.0461    0.0146   -3.1535\n",
      "   -0.0644    0.0152   -4.2279\n",
      "   -0.0072    0.0156   -0.4661\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "[coeff1 se1 coeff1./se1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the FGLS are similar to the White estimates.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "ans =\n",
      "\n",
      "    1.2014    1.2235\n",
      "    0.1012    0.0996\n",
      "    0.0296    0.0296\n",
      "   -0.0004   -0.0004\n",
      "   -0.1655   -0.1648\n",
      "   -0.1115   -0.1136\n",
      "    0.1190    0.1193\n",
      "   -0.0458   -0.0461\n",
      "   -0.0639   -0.0644\n",
      "   -0.0066   -0.0072\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "[coeff coeff1]  % Slight differences"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The White estimators are more restrictive than the FGLS estimators.  However it is always good to check for robustness with respect to estimaiton procedure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Matlab",
   "language": "matlab",
   "name": "matlab"
  },
  "language_info": {
   "codemirror_mode": "octave",
   "file_extension": ".m",
   "help_links": [
    {
     "text": "MetaKernel Magics",
     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "matlab",
   "version": "0.16.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}