{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# R\n",
    "\n",
    "using the [RCall.jl](https://juliainterop.github.io/RCall.jl/stable/) package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CovNWFn"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Printf, DelimitedFiles, LinearAlgebra, Statistics\n",
    "\n",
    "include(\"jlFiles/printmat.jl\")\n",
    "include(\"jlFiles/OlsNW.jl\")         #functions for OLS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample size:    (388,)\n"
     ]
    }
   ],
   "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",
    "\n",
    "printlnPs(\"Sample size:\",size(Rme))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Do OLS (in Julia)\n",
    "\n",
    "use the function sin the file OlsNW.jl to do OLS. Report point estimates and standard errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mOLS Results (assuming iid residuals):\u001b[22m\u001b[39m\n",
      "\n",
      "            b   std_iid\n",
      "c       0.007     0.002\n",
      "SMB     0.217     0.073\n",
      "HML    -0.429     0.074\n",
      "\n"
     ]
    }
   ],
   "source": [
    "Y = Rme\n",
    "T = size(Y,1)\n",
    "X = [ones(T) RSMB RHML]\n",
    "\n",
    "(b,u,Yhat,V,R2) = OlsGMFn(Y,X)\n",
    "std_iid = sqrt.(diag(V))\n",
    "\n",
    "printblue(\"OLS Results (assuming iid residuals):\\n\")\n",
    "xNames = [\"c\",\"SMB\",\"HML\"]\n",
    "printmat([b std_iid],colNames=[\"b\",\"std_iid\"],rowNames=xNames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting Started with RCall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "using RCall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RObject{VecSxp}\n",
      "\n",
      "Call:\n",
      "lm(formula = Y ~ X + 0)\n",
      "\n",
      "Residuals:\n",
      "     Min       1Q   Median       3Q      Max \n",
      "-0.20224 -0.02477  0.00335  0.02663  0.11840 \n",
      "\n",
      "Coefficients:\n",
      "    Estimate Std. Error t value Pr(>|t|)    \n",
      "X1  0.006983   0.002205   3.167  0.00166 ** \n",
      "X2  0.216968   0.073565   2.949  0.00338 ** \n",
      "X3 -0.429088   0.073710  -5.821 1.23e-08 ***\n",
      "---\n",
      "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
      "\n",
      "Residual standard error: 0.04295 on 385 degrees of freedom\n",
      "Multiple R-squared:  0.1488,\tAdjusted R-squared:  0.1422 \n",
      "F-statistic: 22.44 on 3 and 385 DF,  p-value: 2.07e-13\n",
      "\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "RObject{VecSxp}\n",
       "\n",
       "Call:\n",
       "lm(formula = Y ~ X + 0)\n",
       "\n",
       "Coefficients:\n",
       "       X1         X2         X3  \n",
       " 0.006983   0.216968  -0.429088  \n",
       "\n"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@rput X Y\n",
    "\n",
    "ResultsR = reval(\"summary(mod <- lm(Y ~ X+0))\")       #print summary of regression\n",
    "println(ResultsR)\n",
    "\n",
    "resultsR = reval(\"mod <- lm(Y ~ X+0)\")                #get all output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[:coefficients, :residuals, :effects, :rank, Symbol(\"fitted.values\"), :assign, :qr, Symbol(\"df.residual\"), :xlevels, :call, :terms, :model]\n"
     ]
    }
   ],
   "source": [
    "println(names(resultsR))           #print all keys (field names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Task 1\n",
    "\n",
    "Print the Julia and Python estimates (of the coefficients) in a table so we can compare directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Comparing the estimates in Julia and R\n",
      "     0.007     0.007\n",
      "     0.217     0.217\n",
      "    -0.429    -0.429\n",
      "\n"
     ]
    }
   ],
   "source": [
    "b_R = rcopy(resultsR[:coefficients])     #the numerical results are now a Julia vector\n",
    "\n",
    "println(\"Comparing the estimates in Julia and R\")\n",
    "printmat([b b_R])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Task 2\n",
    "\n",
    "Print the smallest and largest values of the difference between the residuals according to Julia and those according to Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-1.942890293094024e-16, 6.453171330633722e-16)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "resid_R = rcopy(resultsR[:residuals]) \n",
    "\n",
    "printmat(extrema(resid_R - u))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OLS (in Julia) with Robust Standard Errors\n",
    "\n",
    "Use standard errors that are robust to heteroskedastcity and autocorrelation (2 lags)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[34m\u001b[1mOLS Results (robust std):\u001b[22m\u001b[39m\n",
      "\n",
      "            b    std_nw\n",
      "c       0.007     0.002\n",
      "SMB     0.217     0.129\n",
      "HML    -0.429     0.118\n",
      "\n"
     ]
    }
   ],
   "source": [
    "(b,u,Yhat,V,R2) = OlsNWFn(Y,X,2)\n",
    "std_nw = sqrt.(diag(V))\n",
    "\n",
    "printblue(\"OLS Results (robust std):\\n\")\n",
    "xNames = [\"c\",\"SMB\",\"HML\"]\n",
    "printmat([b std_nw],colNames=[\"b\",\"std_nw\"],rowNames=xNames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Task 3 \n",
    "\n",
    "Now redo the R estimation with the same sort of robust standard errors. Hint: the `NeweyWest` in the `sandwich`package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "@rlibrary sandwich"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            bstd_nw (R)\n",
      "c       0.007     0.002\n",
      "SMB     0.217     0.129\n",
      "HML    -0.429     0.118\n",
      "\n"
     ]
    }
   ],
   "source": [
    "reval(\"mod <- lm(Y ~ X+0)\")\n",
    "\n",
    "VCV_nwR = reval(NeweyWest(resultsR,lag=2,prewhite=0))\n",
    "std_nwR = sqrt.(diag(rcopy(VCV_nwR)))\n",
    "\n",
    "printmat([b std_nwR],colNames=[\"b\",\"std_nw (R)\"],rowNames=xNames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.7.0",
   "language": "julia",
   "name": "julia-1.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}