source('codes/00_import_libraries.R')


#-----------------------------------------------------------------------#
# shanken correction: with intercept in the cross section ---------------
#-----------------------------------------------------------------------#

fn_shanken_correction = function(factor, asset, error_test = F) {
  
  # this function perfoms the two-pass test with Shanken (1992) correction
  # discussed in Cochrane (2005)
  # first, run a series of time series regressions of asset returns onto factors
  # second, run one cross-sectional regression of mean asset returns onto 
  # factor betas
  
  # INPUT: 
  #     + factor: a TxK matrix of factor returns
  #     + asset: a TxN matrix of asset returns
  #     + error_test: test of whether mean pricing error is zero
  # OUTPUT:
  #     + a list containing: (1) risk premium estimates, (2) cross-sectional R2,
  #                          (3) test of mean pricing is zero (error_test = T)
  
  factor_names = colnames(factor)
  
  factor   = as.matrix(factor)
  asset    = as.matrix(asset)
  mean_ret = apply(asset, 2, mean)
  
  n  = ncol(asset)
  tt = nrow(asset)
  k  = ncol(factor)
  
  
  ### first pass: time-series
  beta  = apply(asset, 2, function(x) coef(lm(x ~ factor))[-1])
  error = apply(asset, 2, function(x) resid(lm(x ~ factor)))
  
  if (k > 1) {
    beta = t(beta)
  }
  
  
  ### second pass: cross-section
  mod  = lm(mean_ret ~ beta)
  lam  = coef(mod)
  r2   = summary(mod)$adj.r.squared*100
  mpe  = mean(abs(resid(mod)))*100
  
  x = cbind(1, beta)
  a = solve(t(x) %*% x) %*% t(x)
  sig_f = var(factor)
  sigf_bar = rbind(
    rep(0, k+1),
    cbind( rep(0, k),  sig_f )
  )
  sig   = cov(error)
  
  
  ### variance of lambda
  v_lam = ( ( (1 + (t(lam[-1]) %*% solve(sig_f) %*% lam[-1])[1,1] ) * 
                a %*% sig %*% t(a) + sigf_bar ) / tt )
  se_lam = sqrt(diag(v_lam))
  t_lam  = lam / se_lam
  p_lam  = pnorm(abs(t_lam), lower.tail = F)*2
  
  
  ### output
  coefs = tibble(
    term      = c('(Intercept)', factor_names),
    estimate  = lam*100,
    std.error = se_lam,
    statistic = t_lam,
    p.value   = p_lam
  )
  
  
  out = list(
    coefs = coefs,
    r2    = r2,
    mpe   = mpe,
    n     = n,
    t     = tt
  )
  
  
  ### chi2 test of pricing errors
  if (error_test == T) {
    
    ### variance of pricing error
    m_x   = diag(n) - x %*% solve(t(x) %*% x) %*% t(x)
    pe    = m_x %*% mean_ret # pricing errors
    v_err = ( ( 1 + (t(lam[-1]) %*% solve(sig_f) %*% lam[-1])[1,1] ) * 
                m_x %*% sig %*% m_x ) / tt
    
    chi2 = t(pe) %*% pinv(v_err) %*% pe
    p_chi2 = pchisq(chi2[1, 1], df = n-k, lower.tail = F)
    
    out = list(
      coefs  = coefs,
      r2     = r2,
      p_chi2 = p_chi2,
      mpe    = mpe,
      n     = n,
      t     = tt
    )
  }
  
  return(out)
}