# 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?

• Maybe some clarification would make it easier to help you. Why do you think that the $(i,j)$ entry is equal to this sum of squared elements? I agree to your point that the output of this small function is the last sum you gave in your question. So, what exactly is your question? Oct 23, 2015 at 20:02
• @muffin1974 the final result in their function is phi.mat, which is a elementwise addition/substraction of three matrices. As far as I see, each term is for one of the terms in $((y_{it}-z_i)(y_{jt}-z_j)-s_{ij})^2$ when expanding the square. May I ask you if you downvoted the question and if so why?
– math
Oct 23, 2015 at 20:11
• I understand that the output of this function contains phi.mat and it is correct that the 3 terms are nothing but the parts evolving when expanding the square. However, although I opened the link you provided and I know the Ledoit/Wolf paper I cannot figure out what exactly you want to know, for me this question either lacks some details or is not precisely showing what you are curious about, this is why I downvoted it. However, given you can clarify your problem with this R implementation I am sure that there is some help=) Oct 23, 2015 at 20:24
• @muffin1974 I edited my question and hope the more detailed explanation is fully sufficient now. Let me know if there is still a point of ambiguity.
– math
Oct 23, 2015 at 20:35

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.

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.

• Thanks for your answer. I edited my question slightly. Its not about the square, its about the matrix multiplication itself. I dont understand how they get the correct entries by the multiplication they suggest.
– math
Oct 22, 2015 at 7:56