code: def calculate_covariance_matrix(X):
      
  
      n_features = X.shape[1]    
      S = np.zeros((n_features, n_features))
      m = np.mean(X, axis=0).reshape(n_features, 1)
  
      for x in X:
          v = x.reshape(n_features, 1) - m
          S += v @ v.T
  
      return 1/(X.shape[0]-1) * S
docs: Calculates the Variance-Covariance matrix
      
      Parameters:
      -----------
      X : array-like, shape (m, n) - the data
      
      Returns:
      --------
      variance_covariance_matrix : array-like, shape(n, n)