Ch04a_OLS_Diagnostics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OLS Diagnostics\n",
    "\n",
    "\n",
    "This notebook tests *(a)* the fit of a regression model; *(b)* properties of the residuals (heteroskedasticity, autocorrelation and lots more).\n",
    "\n",
    "You may also consider the [HypothesisTests.jl](https://github.com/JuliaStats/HypothesisTests.jl) package (not used here)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Packages and Extra Functions\n",
    "\n",
    "The key functions for the diagnostic tests are from the (local) `FinEcmt_OLS` module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "MyModulePath = joinpath(pwd(),\"src\")\n",
    "!in(MyModulePath,LOAD_PATH) && push!(LOAD_PATH,MyModulePath);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "using FinEcmt_OLS, DelimitedFiles, Statistics, LinearAlgebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(388,)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = readdlm(\"Data/FFmFactorsPs.csv\",',',skipstart=1)\n",
    "\n",
    "                #yearmonth, market, small minus big, high minus low\n",
    "(ym,Rme,RSMB,RHML) = (x[:,1],x[:,2]/100,x[:,3]/100,x[:,4]/100)\n",
    "x = nothing\n",
    "println(size(Rme))\n",
    "\n",
    "Y = Rme         #or copy(Rme) is independent copies are needed\n",
    "T = size(Y,1)\n",
    "X = [ones(T) RSMB RHML]\n",
    "k = size(X,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mOLS with traditional standard errors:\u001b[22m\u001b[39m\n",
      "\n",
      "         coef       std\n",
      "c       0.007     0.002\n",
      "SMB     0.217     0.073\n",
      "HML    -0.429     0.074\n",
      "\n"
     ]
    }
   ],
   "source": [
    "(b,u,_,V,R²) = OlsGM(Y,X)    #do OLS\n",
    "Stdb = sqrt.(diag(V))\n",
    "\n",
    "printblue(\"OLS with traditional standard errors:\\n\")\n",
    "xNames = [\"c\",\"SMB\",\"HML\"]\n",
    "printmat([b Stdb],colNames=[\"coef\",\"std\"],rowNames=xNames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression Diagnostics: Testing All Slope Coefficients\n",
    "\n",
    "The `OlsR2Test()` function tests all slope coefficients (or equivalently, the $R^2$) of a regression. Notice that the regression must contain an intercept for R² to be useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "```\n",
       "OlsR2Test(R²,T,df)\n",
       "```\n",
       "\n",
       "Test of all slope coefficients. Notice that the regression must contain an intercept for R² to be useful.\n",
       "\n",
       "### Input\n",
       "\n",
       "  * `R²::Number`:    R² value\n",
       "  * `T::Int`:        number of observations\n",
       "  * `df::Number`:    number of (non-constant) regressors\n",
       "\n",
       "### Output\n",
       "\n",
       "  * `RegrStat::Number`: test statistic\n",
       "  * `pval::Number`:     p-value\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@doc2 OlsR2Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "function OlsR2Test(R²,T,df)\n",
      "    RegrStat = T*R²/(1-R²)           #R\\^2[TAB]\n",
      "    pval     = ccdf(Chisq(df),RegrStat)    #same as 1-cdf()\n",
      "    return RegrStat, pval\n",
      "end\n"
     ]
    }
   ],
   "source": [
    "using CodeTracking\n",
    "println(@code_string OlsR2Test(1.0,1,25))    #print the source code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mTest of all slopes = 0:\u001b[22m\u001b[39m\n",
      "\n",
      "stat     60.165\n",
      "p-val     0.000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "df = k - 1              #number of slope coefficients\n",
    "(RegrStat,pval) = OlsR2Test(R²,T,df)\n",
    "\n",
    "printblue(\"Test of all slopes = 0:\\n\")\n",
    "printmat([RegrStat,pval],rowNames=[\"stat\",\"p-val\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression Diagnostics: Heteroskedasticity\n",
    "\n",
    "The `OlsWhitesTest()` function does White's test for heteroskedasticity. Again, the regression must have an intercept for this test to be useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "```\n",
       "OlsWhitesTest(u,x)\n",
       "```\n",
       "\n",
       "Test of heteroskedasticity. Notice that the regression must contain  an intercept for the test to be useful.\n",
       "\n",
       "### Input\n",
       "\n",
       "  * `u::Vector`:   T-vector, residuals\n",
       "  * `x::Matrix`:   Txk, regressors\n",
       "\n",
       "### Output\n",
       "\n",
       "  * `RegrStat::Number`: test statistic\n",
       "  * `pval::Number`:     p-value\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@doc2 OlsWhitesTest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#println(@code_string OlsWhitesTest([1],[1]))    #print the source code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mWhite's test (H₀: heteroskedasticity is not correlated with regressors):\u001b[22m\u001b[39m\n",
      "\n",
      "stat     77.278\n",
      "p-val     0.000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "(WhiteStat,pval) = OlsWhitesTest(u,X)\n",
    "\n",
    "printblue(\"White's test (H₀: heteroskedasticity is not correlated with regressors):\\n\")\n",
    "printmat([WhiteStat,pval],rowNames=[\"stat\",\"p-val\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression Diagnostics: Autocorrelation of the Residuals\n",
    "\n",
    "The `OlsAutoCorr()` function estimates autocorrelations, calculates the DW and Box-Pierce statistics for the input (often, the residual)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "```\n",
       "OlsAutoCorr(u,L=1)\n",
       "```\n",
       "\n",
       "Test the autocorrelation of OLS residuals\n",
       "\n",
       "### Input\n",
       "\n",
       "  * `u::Vector`:   T-vector, residuals\n",
       "  * `L::Int`:      scalar, number of lags in autocorrelation and Box-Pierce test\n",
       "\n",
       "### Output\n",
       "\n",
       "  * `AutoCorr::Matrix`:   Lx3, autocorrelation, t-stat and p-value\n",
       "  * `BoxPierce::Matrix`:  1x2, Box-Pierce statistic and p-value\n",
       "  * `DW::Number`:         DW statistic\n",
       "\n",
       "### Requires\n",
       "\n",
       "  * StatsBase, Distributions\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@doc2 OlsAutoCorr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#println(@code_string OlsAutoCorr([1],5))    #print the source code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[35m\u001b[1mTesting autocorrelation of residuals\u001b[22m\u001b[39m\n",
      "\n",
      "\u001b[34m\u001b[1mAutocorrelations (lag 1 to 3):\u001b[22m\u001b[39m\n",
      "\n",
      "lag  autocorr    t-stat     p-val\n",
      "1       0.074     1.467     0.142\n",
      "2      -0.037    -0.733     0.464\n",
      "3       0.019     0.377     0.706\n",
      "\n",
      "\n",
      "\u001b[34m\u001b[1mBoxPierce (3 lags): \u001b[22m\u001b[39m\n",
      "stat      2.831\n",
      "p-val     0.418\n",
      "\n",
      "\u001b[34m\u001b[1mDW statistic:\u001b[22m\u001b[39m\n",
      "     1.849\n"
     ]
    }
   ],
   "source": [
    "L = 3     #number of autocorrs to test\n",
    "\n",
    "(ρStats,BoxPierce,DW) = OlsAutoCorr(u,L)\n",
    "\n",
    "printmagenta(\"Testing autocorrelation of residuals\\n\")\n",
    "\n",
    "printblue(\"Autocorrelations (lag 1 to $L):\\n\")\n",
    "printmat(ρStats,colNames=[\"autocorr\",\"t-stat\",\"p-val\"],rowNames=1:L,cell00=\"lag\")\n",
    "\n",
    "printblue(\"\\nBoxPierce ($L lags): \")\n",
    "printmat(BoxPierce',rowNames=[\"stat\",\"p-val\"])\n",
    "\n",
    "printblue(\"DW statistic:\")\n",
    "printlnPs(DW)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Autocorrelation of of `X.*u`\n",
    "\n",
    "What matters most for the uncertainty about a slope coefficient is not the autocorrelation of the residual itself, but of the residual times the regressor. This is tested below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mAutocorrelations of c*u  (lag 1 to 3):\u001b[22m\u001b[39m\n",
      "lag  autocorr    t-stat     p-val\n",
      "1       0.074     1.467     0.142\n",
      "2      -0.037    -0.733     0.464\n",
      "3       0.019     0.377     0.706\n",
      "\n",
      "\u001b[34m\u001b[1mAutocorrelations of SMB*u  (lag 1 to 3):\u001b[22m\u001b[39m\n",
      "lag  autocorr    t-stat     p-val\n",
      "1       0.219     4.312     0.000\n",
      "2      -0.014    -0.268     0.789\n",
      "3       0.044     0.857     0.391\n",
      "\n",
      "\u001b[34m\u001b[1mAutocorrelations of HML*u  (lag 1 to 3):\u001b[22m\u001b[39m\n",
      "lag  autocorr    t-stat     p-val\n",
      "1       0.278     5.472     0.000\n",
      "2       0.131     2.582     0.010\n",
      "3       0.225     4.438     0.000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for i = 1:k         #iterate over different regressors\n",
    "    #local pStats\n",
    "    ρStats, = OlsAutoCorr(X[:,i].*u,L)\n",
    "    printblue(\"Autocorrelations of $(xNames[i])*u  (lag 1 to $L):\")\n",
    "    printmat(ρStats,colNames=[\"autocorr\",\"t-stat\",\"p-val\"],rowNames=1:L,cell00=\"lag\")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Measures of Fit\n",
    "\n",
    "Adjusted R², AIC, BIC "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "```\n",
       "RegressionFit(u,R²,k)\n",
       "```\n",
       "\n",
       "Calculate adjusted R², AIC and BIC from regression residuals.\n",
       "\n",
       "### Input\n",
       "\n",
       "  * `u::Vector`:      T-vector of residuals\n",
       "  * `R²::Float`:      the R² value\n",
       "  * `k::Int`:         number of regressors\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@doc2 RegressionFit\n",
    "#println(@code_string RegressionFit([1],0.0,3))    #print the source code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mMeasures of fit\u001b[22m\u001b[39m\n",
      "R²        0.134\n",
      "R²adj     0.130\n",
      "AIC      -6.285\n",
      "BIC      -6.255\n",
      "\n"
     ]
    }
   ],
   "source": [
    "(R²adj,AIC,BIC) = RegressionFit(u,R²,k)\n",
    "\n",
    "printblue(\"Measures of fit\")\n",
    "printmat([R²,R²adj,AIC,BIC];rowNames=[\"R²\",\"R²adj\",\"AIC\",\"BIC\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test of Normality\n",
    "\n",
    "of the residuals, applying the Bera-Jarque test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "```\n",
       "BeraJarqueTest(x)\n",
       "```\n",
       "\n",
       "Calculate the BJ test for each column in a matrix. Reports `(skewness,kurtosis,BJ)`.\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@doc2 BeraJarqueTest\n",
    "#println(@code_string BeraJarqueTest([1]))    #print the source code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mTest of normality\u001b[22m\u001b[39m\n",
      "                 stat   p-value\n",
      "skewness       -0.746     0.000\n",
      "kurtosis        5.583     0.000\n",
      "Bera-Jarque   143.834     0.000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "(skewness,kurtosis,BJ,pvals) = BeraJarqueTest(u)\n",
    "\n",
    "printblue(\"Test of normality\")\n",
    "xut = vcat(skewness,kurtosis,BJ)\n",
    "printmat(xut,collect(pvals);rowNames=[\"skewness\",\"kurtosis\",\"Bera-Jarque\"],colNames=[\"stat\",\"p-value\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multicollinearity\n",
    "\n",
    "by studying the correlation matrix and the variance inflation factor (VIF). A high VIF (5 to 10) might indicate issues with multicollinearity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "```\n",
       "VIF(X)\n",
       "```\n",
       "\n",
       "Calculate the variance inflation factor\n",
       "\n",
       "### Input\n",
       "\n",
       "  * `x::Matrix`:    Txk matrix with regressors\n",
       "\n",
       "### Output\n",
       "\n",
       "  * `maxVIF::Float`:     highest VIF value\n",
       "  * `allVIF::Vector`:    a k VIF values\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@doc2 VIF\n",
    "#println(@code_string VIF([1]))    #print the source code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mCorrelation matrix (checking multicollinearity)\u001b[22m\u001b[39m\n",
      "            c       SMB       HML\n",
      "c       1.000       NaN       NaN\n",
      "SMB       NaN     1.000    -0.320\n",
      "HML       NaN    -0.320     1.000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "printblue(\"Correlation matrix (checking multicollinearity)\")\n",
    "printmat(cor(X);colNames=xNames,rowNames=xNames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mVIF (checking multicollinearity)\u001b[22m\u001b[39m\n",
      "c       1.000\n",
      "SMB     1.114\n",
      "HML     1.114\n",
      "\n"
     ]
    }
   ],
   "source": [
    "(maxVIF,allVIF) = VIF(X)\n",
    "printblue(\"VIF (checking multicollinearity)\")\n",
    "printmat(allVIF;rowNames=xNames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A Convenience Function for Printing All These Tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "```\n",
       "DiagnosticsTable(X,u,R²,nlags,xNames=\"\")\n",
       "```\n",
       "\n",
       "Compute and print a number of regression diagnostic tests.\n",
       "\n",
       "### Input\n",
       "\n",
       "  * `X::Matrix`:      Txk matrix of regressors\n",
       "  * `u::Vector`:      T-vector of residuals\n",
       "  * `R²::Float`:      the R² value\n",
       "  * `nlags::Int`:     number of lags to use in autocorrelation test\n",
       "  * `xNames::Vector`: of strings, regressor names\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@doc2 DiagnosticsTable\n",
    "#println(@code_string DiagnosticsTable([1],[1],0.0,1))    #print the source code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mTest of all slopes = 0\u001b[22m\u001b[39m\n",
      "stat     60.165\n",
      "p-val     0.000\n",
      "\n",
      "\u001b[34m\u001b[1mWhite's test (H₀: heteroskedasticity is not correlated with regressors)\u001b[22m\u001b[39m\n",
      "stat     77.278\n",
      "p-val     0.000\n",
      "\n",
      "\u001b[34m\u001b[1mTesting autocorrelation of residuals (lag 1 to 3)\u001b[22m\u001b[39m\n",
      "lag  autocorr    t-stat     p-val\n",
      "1       0.074     1.467     0.142\n",
      "2      -0.037    -0.733     0.464\n",
      "3       0.019     0.377     0.706\n",
      "\n",
      "\u001b[34m\u001b[1mBoxPierce (3 lags) \u001b[22m\u001b[39m\n",
      "stat      2.831\n",
      "p-val     0.418\n",
      "\n",
      "\u001b[34m\u001b[1mDW statistic\u001b[22m\u001b[39m\n",
      "     1.849          \n",
      "\n",
      "\u001b[34m\u001b[1mAutocorrelations of c*u  (lag 1 to 3)\u001b[22m\u001b[39m\n",
      "lag  autocorr    t-stat     p-val\n",
      "1       0.074     1.467     0.142\n",
      "2      -0.037    -0.733     0.464\n",
      "3       0.019     0.377     0.706\n",
      "\n",
      "\u001b[34m\u001b[1mAutocorrelations of SMB*u  (lag 1 to 3)\u001b[22m\u001b[39m\n",
      "lag  autocorr    t-stat     p-val\n",
      "1       0.219     4.312     0.000\n",
      "2      -0.014    -0.268     0.789\n",
      "3       0.044     0.857     0.391\n",
      "\n",
      "\u001b[34m\u001b[1mAutocorrelations of HML*u  (lag 1 to 3)\u001b[22m\u001b[39m\n",
      "lag  autocorr    t-stat     p-val\n",
      "1       0.278     5.472     0.000\n",
      "2       0.131     2.582     0.010\n",
      "3       0.225     4.438     0.000\n",
      "\n",
      "\u001b[34m\u001b[1mMeasures of fit\u001b[22m\u001b[39m\n",
      "R²        0.134\n",
      "R²adj     0.130\n",
      "AIC      -6.285\n",
      "BIC      -6.255\n",
      "\n",
      "\u001b[34m\u001b[1mTest of normality\u001b[22m\u001b[39m\n",
      "                 stat   p-value\n",
      "skewness       -0.746     0.000\n",
      "kurtosis        5.583     0.000\n",
      "Bera-Jarque   143.834     0.000\n",
      "\n",
      "\u001b[34m\u001b[1mCorrelation matrix (checking multicollinearity)\u001b[22m\u001b[39m\n",
      "            c       SMB       HML\n",
      "c       1.000       NaN       NaN\n",
      "SMB       NaN     1.000    -0.320\n",
      "HML       NaN    -0.320     1.000\n",
      "\n",
      "\u001b[34m\u001b[1mVIF (checking multicollinearity)\u001b[22m\u001b[39m\n",
      "c       1.000\n",
      "SMB     1.114\n",
      "HML     1.114\n",
      "\n"
     ]
    }
   ],
   "source": [
    "DiagnosticsTable(X,u,R²,3,xNames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 1.11.1",
   "language": "julia",
   "name": "julia-1.11"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}