Implementation of Ledoit Wolf shrinkage estimator within R package tawny

Question

I want to implement the shrinkage intensity given by Ledoit and Wolf, see here page 13. They define $y_{it}$ with $1\le i\le N$ and $1\le t\le t$ be the return on stock $i$ at time $t$. Moreover, $z_i:=\bar{y}_i:=\frac{1}{T}\sum_{t=1}^Ty_{it}$ the mean estimator of the $i$-th stock. As they explain the optimal shrinkage intensity estimator is given by (for simplicity I drop the hat notation)

$$\kappa = \frac{\pi-\rho}{\gamma}$$

I have a question about the implementation of $\pi$. They define

$$\pi_{ij} = \frac{1}{T}\sum_{t=1}^T((y_{it}-z_i)(y_{jt}-z_j)-s_{ij})^2$$

I wanted to implement this in a efficient way. The package 'tawny' has a preimplemented function of this. I checked the source code, which can be found here in the file shrinkage.R. They use (I copy):

# Sum of the asymptotic variances
# returns : T x N (zoo) - Matrix of asset returns
# sample : N x N - Sample covariance matrix
# Used internally.
# S <- cov.sample(ys)
# ys.p <- shrinkage.p(ys, S)
shrinkage.p <- function(returns, sample)
{
  T <- nrow(returns)
  N <- ncol(returns)
  ones <- rep(1,T)
  means <- t(returns) %*% ones / T
  z <- returns - matrix(rep(t(means), T), ncol=N, byrow=TRUE)

  term.1 <- t(z^2) %*% z^2
  term.2 <- 2 * sample * (t(z) %*% z)
  term.3 <- sample^2
  phi.mat <- (term.1 - term.2 + term.3) / T

  phi <- list()
  phi$sum <- sum(phi.mat)
      phi$diags <- diag(phi.mat)
  phi
}

The input is the transposed $y$ as they define it to be a $T\times N$ matrix. Their final output I'm interested in is the following:

phi.mat <- (term.1 - term.2 + term.3) / T

If their code works properly (which for sure it does) its true that $(i,j)$ entry of phi.mat is equal $\pi_{i,j}$. However, this is the point I cant figure out why this is true. Since in R addition / subtraction of matrices is elementwise it should be true that term.1 is equal to

$$\sum_{t=1}^T(y_{it}-z_i)^2(y_{jt}-z_j)^2$$

and term.2 equal to

$$2\sum_{t=1}^T(y_{it}-z_i)(y_{jt}-z_j)s_{ij}$$

and last term.3 equal to

$$\sum_{t=1}^T s_{ij}^2$$

I understand term.3. However, term.1 and term.2 are not that clear to me. I guess understanding one of them will help me understand the other. So lets focus on term.1:

Question Why is it true that $(i,j)$ entry of their code term.1 is equal to $\sum_{t=1}^T(y_{it}-z_i)^2(y_{jt}-z_j)^2$? Isn't the matrix product, which they suggest, equal to
$$\sum_{t=1}^T(y_{it}-z_i)(y_{tj}-z_j)$$ and the matrix of returns is clearly not symmetric?

Answer 1

The question you asked can be explained by these two lines of the code

e   means <- t(returns) %*% ones / T
z <- returns - matrix(rep(t(means), T), ncol=N, byrow=TRUE)
term.1 <- t(z^2) %*% z^2 e

Here returns is TxN which gives you matrix ${y_{nt}}$ where n has i and j elements ; means is TxN of matrix ${z_i}$ , same mean for each asset for the time series. Therefore, line 2 of code above gives you z variable which is matrix ${y_{nt}-z_i}$. When you take a transpose and multiply you get a sum of cross multiplication elements. Elements were squared when you do z^2. therefore you are taking a transpose of matrix of squared elements like $({y_{nt}-z_i})^2$ and multiplying by matrix $({y_{nt}-z_i})^2$. That is how you get cross multiplication of the squared elements for term.1. Likewise for term 2, the code uses t(z) and not t(z^2) which explains the elements are not squared. I hope this explains it.

Answer 2

I believe the "^2" operation is performed element by element in this instance. The matrix square in R is ... um... different. Memory fails me on the exact notation.