code: def calc_asymptotic_covariance(hessian, fisher_info_matrix):
      
      # Calculate the inverse of the hessian
      hess_inv = scipy.linalg.inv(hessian)
  
      return np.dot(hess_inv, np.dot(fisher_info_matrix, hess_inv))
docs: Parameters
      ----------
      hessian : 2D ndarray.
          It should have shape `(num_vars, num_vars)`. It is the matrix of second
          derivatives of the total loss across the dataset, with respect to each
          pair of coefficients being estimated.
      fisher_info_matrix : 2D ndarray.
          It should have a shape of `(num_vars, num_vars)`.  It is the
          approximation of the negative of the expected hessian formed by taking
          the outer product of (each observation's gradient of the loss function)
          with itself, and then summing across all observations.
  
      Returns
      -------
      huber_white_matrix : 2D ndarray.
          Will have shape `(num_vars, num_vars)`. The entries in the returned
          matrix are calculated by the following formula:
          `hess_inverse * fisher_info_matrix * hess_inverse`.