code: def jacobian(func, x, nf, nobs, *args, **kwargs):
      
      nargs = len(x)  # degrees of freedom
      f = lambda x_: func(x_, *args, **kwargs)
      j = np.zeros((nargs, nf, nobs))  # matrix of zeros
      for n in xrange(nargs):
          j[n] = partial_derivative(f, x, n, nargs)
          # better to transpose J once than transpose partial derivative each time
          # j[:,:,n] = df.T
      return j.T
docs: Estimate Jacobian matrices :math:`\\frac{\\partial f_i}{\\partial x_{j,k}}`
      where :math:`k` are independent observations of :math:`x`.
  
      The independent variable, :math:`x`, must be a numpy array with exactly 2
      dimensions. The first dimension is the number of independent arguments,
      and the second dimensions is the number of observations.
  
      The function must return a Numpy array with exactly 2 dimensions. The first
      is the number of returns and the second dimension corresponds to the number
      of observations. If the input argument is 2-D then the output should also
      be 2-D
  
      Constant arguments can be passed as additional positional arguments or
      keyword arguments. If any constant argument increases the number of
      observations of the return value, tile the input arguments to match.
  
      Use :func:`numpy.atleast_2d` or :func:`numpy.reshape` to get the
      correct dimensions for scalars.
  
      :param func: function
      :param x: independent variables grouped by observation
      :param nf: number of return in output (1st dimension)
      :param nobs: number of observations in output (2nd dimension)
      :return: Jacobian matrices for each observation