diff --git a/shrink_covariance_matrix.py b/shrink_covariance_matrix.py
new file mode 100644
index 0000000..97ca620
--- /dev/null
+++ b/shrink_covariance_matrix.py
@@ -0,0 +1,70 @@
+import numpy as np
+import pandas
+
+
+def shrink_covariance_matrix(x, shrink=None):
+    """
+    Shrinks towards constant correlation matrix
+    If shrink is specified then this const is used for shrinkage
+
+    The notation follows Ledoit and Wolf (2004)
+    http://www.ledoit.net/honey_abstract.htm
+    This version: 06/2009
+
+    Parameters
+    ----------
+    x : N x N sample covariance matrix of stock returns
+    shrink : given shrinkage intensity factor; if none, code calculates
+
+    Returns
+    -------
+    tuple : numpy.ndarray which contains the shrunk covariance matrix
+          : float shrinkage intensity factor
+
+    """
+
+    if x is None:
+        raise ValueError('No covariance matrix defined')
+
+    if isinstance(x, pandas.DataFrame):
+        cov = x.as_matrix()
+    elif isinstance(x, np.ndarray):
+        cov = x
+    else:
+        raise ValueError('Covariance matrix passed must be numpy.ndarray')
+
+    [t, n] = np.shape(cov)
+    meanx = cov.mean(axis=0)
+    cov = cov - np.tile(meanx, (t, 1.))
+
+    sample = (1. / t) * np.dot(cov.T, cov)
+
+    var = np.diag(sample)
+    sqrtvar = np.sqrt(var)
+
+    a = np.tile(sqrtvar, (n, 1.))
+    rho = (sum(sum(sample / (a * a.T))) - n) / (n * (n - 1.))
+
+    prior = rho * (a * a.T)
+    prior[np.eye(t, n) == 1.] = var
+
+    # Frobenius-norm of matrix cov, sqrt(sum(diag(dot(cov.T, cov))))
+    c = np.linalg.norm(sample - prior, 'fro') ** 2.
+    y = cov ** 2.
+    p = np.dot((1. / t), sum(sum(np.dot(y.T, y)))) - sum(sum(sample ** 2.))
+    rdiag = np.dot((1. / t), sum(sum(y ** 2.))) - sum(var ** 2.)
+    v = np.dot((cov ** 3.).T, cov) / t - (var * sample).T
+    v[np.eye(t, n) == 1.] = 0.
+    roff = sum(sum(v * (a / a.T)))
+    r = rdiag + np.dot(rho, roff)
+
+    # compute shrinkage constant
+    if shrink:
+        shrinkage = shrink
+    else:
+        k = (p - r) / c
+        shrinkage = max(0., min(1., k / t))
+
+    sigma = np.dot(shrinkage, prior) + np.dot((1. - shrinkage), sample)
+
+    return sigma, shrinkage