Efficiently Computing Ledoit–Wolf Covariance Shrinkage Terms with Matrix Products

# Implementation of Ledoit Wolf shrinkage estimator within R package tawny


# Implementation of Ledoit Wolf shrinkage estimator within R package tawny












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 by user12348 (score 3, accepted)

https://quant.stackexchange.com/a/21432

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 by Drew (score -1)

https://quant.stackexchange.com/a/21357

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.