examples_homoskedastic_and_hac_noise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Code examples for Simons (2023)\n",
    "\n",
    "This notebook details the usage of the function `eigenvector_estimates_cis` applied to two examples. The first uses the built-in covariance matrix estimator based on the outer product and enforces homoskedasticity across columns of $M$. The example uses the `CovarianceMatrices` package to estimate a HAC-robust matrix.\n",
    "\n",
    "## Inference for eigenvectors with homoskedastic noise\n",
    "\n",
    "We have the model $M_{t} = M + \\varepsilon_{t}$ and take the mean $\\hat{M}$. We assume that the estimation error matrix $\\hat{M} - M$ has $\\text{i.i.d.}$ columns $ \\left[\\textbf{a}_{1} \\, \\textbf{a}_{2} \\, \\dots \\, \\textbf{a}_{p} \\right] $ with finite covariance as per $\\mathbb{E} \\textbf{a}_{i} \\textbf{a}_{j}^{\\intercal} = \\Omega \\delta_{ij}$. Therefore, $\\Psi = I_p \\otimes \\Omega$ as in the standard setup of the paper. We use "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×1×3 Array{Float64, 3}:\n",
       "[:, :, 1] =\n",
       " -1.6473512926853844\n",
       "  1.0\n",
       "\n",
       "[:, :, 2] =\n",
       " -1.9056459093950218\n",
       "  1.0\n",
       "\n",
       "[:, :, 3] =\n",
       " -1.389056675975747\n",
       "  1.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "include(\"eigenvector_estimates_cis.jl\")\n",
    "using Statistics, Random\n",
    "M = [1 -3; -1 1];\n",
    "T = 50;\n",
    "p = size(M)[1];\n",
    "q = 1;\n",
    "r=p-q;\n",
    "A = randn(MersenneTwister(1234),Float64,p,p,T) .+ M;\n",
    "eigenvector_estimates_cis(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the above with the population analogue. Note that the first entry of the third array dimension is the point estimate, followed by the lower and upper bounds according to the critical value of $1.96$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×1 Matrix{Float64}:\n",
       " -1.732050807568877"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ev,R = eigen(M,sortby = x -> (-floor(real(x), digits = 6), floor(imag(x), digits = 6)))\n",
    "R[1:r,1:q] * pinv(R[(r+1):p,1:q])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression on an intercept example with HAC covariance estimation\n",
    "\n",
    "If we confront a case where we suspect that our error sequence is not homoskedastic, we can use a customised covariance matrix. For example, we may want to use a Newey-West-type robust covariance matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×1×3 Array{Float64, 3}:\n",
       "[:, :, 1] =\n",
       " -1.6473512926853844\n",
       "  1.0\n",
       "\n",
       "[:, :, 2] =\n",
       " -3.381733547962595\n",
       "  1.0\n",
       "\n",
       "[:, :, 3] =\n",
       " 0.08703096259182597\n",
       " 1.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "using CovarianceMatrices\n",
    "T=50;\n",
    "M = [1 -3; -1 1];\n",
    "p=size(M)[1]\n",
    "A = randn(MersenneTwister(1234),Float64,2,2,T) .+ M;\n",
    "vectorised_matrices = reshape(A,(p^2,T));\n",
    "y = vectorised_matrices'\n",
    "X   = ones(T,1)\n",
    "_,K = size(X)\n",
    "b   = X\\y\n",
    "res = y - X*b\n",
    "\n",
    "kernel = QuadraticSpectralKernel{NeweyWest}()\n",
    "bw = CovarianceMatrices.optimalbandwidth(kernel, res, prewhite=false);\n",
    "Ψ   = lrvar(QuadraticSpectralKernel(bw), res, scale = T^2/(T-K))   # df adjustment is built into vcov\n",
    "eigenvector_estimates_cis(A,1,1.96,false,Ψ)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.10.0",
   "language": "julia",
   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}