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? – muffin1974 Oct 23 '15 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 '15 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=) – muffin1974 Oct 23 '15 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 '15 at 20:35

e   means <- t(returns) %*% ones / T

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.