I have calculated weights of selected assets in a market-neutral portfolio (presumably with min variance) using PCA and simple data covariance matrix.
The question is :
It is obvious that Cov Matrix contains estimation error so i need to shrink it, the best choice probably would be to implement Ledoit-Wolf estimator but seems that i do not completely understand the formula for Shrinkage Intensity here, could somebody provide a link displaying how to use it on finite samples?
As an alternative i also found this explanation why covariance matrix is often inverted before solving eigenproblem : http://scikit-learn.org/stable/modules/covariance.html#sparse-inverse-covariance
As far as i understood - inverse of covariance matrix it is so called precision matrix (or sparse estimator) which helps to eliminate spurious relationships caused by estimation error of covariance matrix and thus i would like to ask - can simply inverted covariance matrix be a replacement for more complex estimator like Ledoit-Wolf?
P. S. in general i need a good estimator for covariance matrix but in case of Ledoit-Wolf would be good to see how it is calculated on a real data to understand how exactly it works and how to implement it e.g. using C++.
Update : as for the question about general algorithm of the code below, as far as i understand Ledoit-Wolf estimator is kind of mix of the smart averaging and weighting coefficients and the main idea is to decrease spurious variances by recalculating them with averaged values where each values has its own intensity.
The more basic algorithm (just a suggestion) would be to set growing weighting coefficients from recent return to older ones (sorry, do not remember exact book where it was described).
Unfortunately, i tried to apply this estimator on Forex market and could not find significant differences between original covariance matrix and the one with applied estimator. This is part of another question but i decided to use existing library for PCA instead of improving my own so i switched to alglib.
Edit. Ok, i implemented mentioned estimator and to make sure whether it is correct or not i am trying to do the same thing in R. At the moment i have following implementation.
struct SSeries // supposed to be an array of structures of different assets
{
string mName;
double mMean;
double mPoints[]; // returns of specific asset
};
// series - array of returns, size N x M, N - number of assets, M - returns
// covs - covariance matrix, N X N
// estimators - buffer for calculated estimation matrix
// order = N
// depth = M
void GetEstimator(SSeries& series[], SSeries& covs[], SSeries& estimators[], const int order, const int depth)
{
double
R = 0, // average return (first 2 lines with formulas in Appendix A)
A = 0, // value of "П" in the formula below (Pi)
B = 0, // value of "p"
C = 0, // value of "v" (gamma)
cross = 0;
// Shrink Target
for (int k = 0; k < order; k++)
{
double summary = 0;
for (int n = k + 1; n < order; n++)
{
double divisor = covs[k].mPoints[k] * covs[n].mPoints[n] > 0 ? 1 : -1; // just to make sure that MathSqrt will apply to positive number (square root)
summary += covs[k].mPoints[n] / MathSqrt(covs[k].mPoints[k] * covs[n].mPoints[n] * divisor) * divisor;
}
R += summary;
}
R = R * 2 / ((order - 1) * order); // average return
// Shrink Intensity - weighting coefficient for items in covariance matrix k = (П - p) / v
for (int k = 0; k < order; k++)
{
double summaryA = 0;
for (int n = 0; n < order; n++)
{
double Akn = 0; // equivalent of П[i][j]
for (int i = 0; i < depth; i++)
{
double X = series[k].mPoints[i] - series[k].mMean; // substract mean from first series
double Y = series[n].mPoints[i] - series[n].mMean; // substract mean from second series
Akn += (X * Y - covs[k].mPoints[n]) * (X * Y - covs[k].mPoints[n]);
}
double divisor = covs[k].mPoints[k] * covs[n].mPoints[n] > 0 ? 1 : -1; // make value positiove before taking its square root
estimators[k].mPoints[n] = k == n ? covs[k].mPoints[n] : R * MathSqrt(covs[k].mPoints[k] * covs[n].mPoints[n] * divisor) * divisor; // for diagonal elements use variance instead of covariance
summaryA += Akn / depth; // calculation of "П" in the formula above
cross += k == n ? summaryA : 0; // calculate only variances (diagonal)
}
A += summaryA; // "П" - sum of asymptotic variances of the entries of the sample covariance matrix scaled by
}
for (int k = 0; k < order; k++)
{
double summaryB = 0, summaryC = 0; // calculation of "p" and "v" using existing cov matrix and recently calculated average return R
for (int n = 0; n < order; n++)
{
double aCovKk = 0, aCovNn = 0;
for (int i = 0; i < depth; i++)
{
double X = series[k].mPoints[i] - series[k].mMean;
double Y = series[n].mPoints[i] - series[n].mMean;
aCovKk += (X * X - covs[k].mPoints[k]) * (X * Y - covs[k].mPoints[n]);
aCovNn += (Y * Y - covs[n].mPoints[n]) * (X * Y - covs[k].mPoints[n]);
}
double divisorA = covs[n].mPoints[n] / covs[k].mPoints[k] > 0 ? 1 : -1;
double divisorB = covs[k].mPoints[k] / covs[n].mPoints[n] > 0 ? 1 : -1;
aCovKk /= aCovKk;
aCovNn /= aCovNn;
summaryB += k == n ? 0 : (R / 2) * (aCovKk * MathSqrt(covs[n].mPoints[n] / covs[k].mPoints[k] * divisorA) * divisorA + aCovNn * MathSqrt(covs[k].mPoints[k] / covs[n].mPoints[n] * divisorB) * divisorB);
summaryC += (estimators[k].mPoints[n] - covs[k].mPoints[n]) * (estimators[k].mPoints[n] - covs[k].mPoints[n]);
}
B += summaryB;
C += summaryC;
}
B += cross; // "p" supposed to be a sum of asymptotic variances and covariances
double intensity = MathMax(0, MathMin(1, ((A - B) / C) / depth)); // formula for weighting coefficient
for (int k = 0; k < order; k++)
{
for (int n = 0; n < order; n++)
{
estimators[k].mPoints[n] = estimators[k].mPoints[n] * intensity + (1 - intensity) * covs[k].mPoints[n]; // apply calculated intensity - formula (2) in the chapter "3.3 Shrinkage Constant"
}
}
}
