How to use the linearmodels.compat.numpy.lstsq function in linearmodels

gls_eps = np.reshape(gls_eps, (k, gls_eps.shape[0] // k)).T
        eps = np.reshape(eps, (k, eps.shape[0] // k)).T
        cov_est = cov_est(self._wxhat, gls_eps, sigma, full_sigma, gls=True,
                          constraints=self._constraints, **cov_config)
        cov = cov_est.cov

        # Repackage results for individual equations
        individual = AttrDict()
        debiased = cov_config.get('debiased', False)
        method = 'Iterative GLS' if iter_count > 1 else 'GLS'
        for i in range(k):
            cons = int(self.has_constant.iloc[i])

            if cons:
                c = np.sqrt(self._w[i])
                ye = self._wy[i] - c @ lstsq(c, self._wy[i])[0]
            else:
                ye = self._wy[i]
            total_ss = float(ye.T @ ye)
            stats = self._common_indiv_results(i, beta, cov, gls_eps, eps,
                                               method, cov_type, cov_est, iter_count,
                                               debiased, cons, total_ss)
            key = self._eq_labels[i]
            individual[key] = stats

        # Populate results dictionary
        nobs = eps.size
        results = self._common_results(beta, cov, method, iter_count, nobs,
                                       cov_type, sigma, individual, debiased)

        # wresid is different between GLS and OLS
        wresid = []

# Covariance estimation
        cov_est = COV_EST[cov_type]
        cov_est = cov_est(self._wxhat, eps, sigma, sigma, gls=False,
                          constraints=self._constraints, **cov_config)
        cov = cov_est.cov

        individual = AttrDict()
        debiased = cov_config.get('debiased', False)
        for i in range(k):
            wy = wye = self._wy[i]
            w = self._w[i]
            cons = int(self.has_constant.iloc[i])
            if cons:
                wc = np.ones_like(wy) * np.sqrt(w)
                wye = wy - wc @ lstsq(wc, wy)[0]
            total_ss = float(wye.T @ wye)
            stats = self._common_indiv_results(i, beta, cov, eps, eps, 'OLS',
                                               cov_type, cov_est, 0, debiased, cons, total_ss)
            key = self._eq_labels[i]
            individual[key] = stats

        nobs = eps.size
        results = self._common_results(beta, cov, 'OLS', 0, nobs, cov_type,
                                       sigma, individual, debiased)
        results['wresid'] = results.resid
        results['cov_estimator'] = cov_est
        results['cov_config'] = cov_est.cov_config

        return SystemResults(results)

return True, int(loc)

    if np.any((np.ptp(x, axis=0) == 0) & ~np.all(x == 0, axis=0)):
        loc = (np.ptp(x, axis=0) == 0) & ~np.all(x == 0, axis=0)
        loc = np.argwhere(loc)
        return True, int(loc)

    n = x.shape[0]
    aug_rank = np.linalg.matrix_rank(np.c_[np.ones((n, 1)), x])
    rank = np.linalg.matrix_rank(x) if x_rank is None else x_rank

    has_const = (aug_rank == rank) and x.shape[0] > x.shape[1]
    has_const = has_const or rank < min(x.shape)
    loc = None
    if has_const:
        out = lstsq(x, np.ones((n, 1)))
        beta = out[0].ravel()
        loc = np.argmax(np.abs(beta) * x.var(0))
    return bool(has_const), loc

Notes
        -----
        The kernel covariance estimator takes the optional arguments
        ``kernel``, one of 'bartlett', 'parzen' or 'qs' (quadratic spectral)
        and ``bandwidth`` (a positive integer).
        """
        nobs, nf, nport, nrf, s1, s2, s3 = self._boundaries()
        excess_returns = not self._risk_free
        f = self.factors.ndarray
        p = self.portfolios.ndarray
        nport = p.shape[1]

        # Step 1, n regressions to get B
        fc = np.c_[np.ones((nobs, 1)), f]
        b = lstsq(fc, p)[0]  # nf+1 by np
        eps = p - fc @ b
        if excess_returns:
            betas = b[1:].T
        else:
            betas = b.T.copy()
            betas[:, 0] = 1.0

        sigma_m12 = self._sigma_m12
        lam = lstsq(sigma_m12 @ betas, sigma_m12 @ p.mean(0)[:, None])[0]
        expected = betas @ lam
        pricing_errors = p - expected.T
        # Moments
        alphas = pricing_errors.mean(0)[:, None]
        moments = self._moments(eps, betas, lam, alphas, pricing_errors)
        # Jacobian
        jacobian = self._jacobian(betas, lam, alphas)

drop_first = True
        if self.other_effects:
            oe = self._other_effect_cats.dataframe
            for c in oe:
                dummies = pd.get_dummies(oe[c], drop_first=drop_first).astype(np.float64)
                d.append(dummies.values)
                drop_first = True

        d = np.column_stack(d)
        wd = root_w * d
        if self.has_constant:
            wd -= root_w * (w.T @ d / w.sum())
            z = np.ones_like(root_w)
            d -= z * (z.T @ d / z.sum())

        x_mean = lstsq(wd, x)[0]
        y_mean = lstsq(wd, y)[0]

        # Save fitted unweighted effects to use in eps calculation
        x_effects = d @ x_mean
        y_effects = d @ y_mean

        # Purge fitted, weighted values
        x = x - wd @ x_mean
        y = y - wd @ y_mean

        ybar = root_w @ lstsq(root_w, y)[0]
        return y, x, ybar, y_effects, x_effects

df_resid = wy.shape[0] - wx.shape[1]
        cov_est, cov_config = self._choose_cov(cov_type, **cov_config)
        cov = cov_est(wy, wx, params, self.dependent.entity_ids, self.dependent.time_ids,
                      debiased=debiased, **cov_config)

        weps = wy - wx @ params
        eps = weps / root_w
        index = self.dependent.index
        fitted = pd.DataFrame(self.exog.values2d @ params, index, ['fitted_values'])
        effects = pd.DataFrame(self.dependent.values2d - np.asarray(fitted) - eps, index,
                               ['estimated_effects'])
        idiosyncratic = pd.DataFrame(eps, index, ['idiosyncratic'])
        residual_ss = float(weps.T @ weps)
        wmu = 0
        if self.has_constant:
            wmu = root_w * lstsq(root_w, wy)[0]
        wy_demeaned = wy - wmu
        total_ss = float(wy_demeaned.T @ wy_demeaned)
        r2 = 1 - residual_ss / total_ss

        res = self._postestimation(params, cov, debiased, df_resid, weps, wy, wx, root_w)
        res.update(dict(df_resid=df_resid, df_model=x.shape[1], nobs=y.shape[0],
                        residual_ss=residual_ss, total_ss=total_ss, r2=r2,
                        resids=eps, wresids=weps, index=index, sigma2_eps=sigma2_e,
                        sigma2_effects=sigma2_u, rho=rho, theta=theta_out,
                        fitted=fitted, effects=effects, idiosyncratic=idiosyncratic))

        return RandomEffectsResults(res)

loc = np.argwhere(np.all(x == 1, axis=0))
        return True, int(loc)

    if np.any((np.ptp(x, axis=0) == 0) & ~np.all(x == 0, axis=0)):
        loc = np.any((np.ptp(x, axis=0) == 0) & ~np.all(x == 0, axis=0))
        loc = np.argwhere(loc)
        return True, int(loc)

    n = x.shape[0]
    aug_rank = matrix_rank(np.c_[np.ones((n, 1)), x])
    rank = matrix_rank(x) if x_rank is None else x_rank

    has_const = bool(aug_rank == rank)
    loc = None
    if has_const:
        out = lstsq(x, np.ones((n, 1)))
        beta = out[0].ravel()
        loc = np.argmax(np.abs(beta) * x.var(0))
    return has_const, loc

def single(z: pd.DataFrame):
            exog = z.iloc[:, 1:].values
            if exog.shape[0] < exog.shape[1] or matrix_rank(exog) != exog.shape[1]:
                return pd.Series([np.nan] * len(z.columns), index=z.columns)
            dep = z.iloc[:, :1].values
            params = lstsq(exog, dep)[0]
            return pd.Series(np.r_[np.nan, params.ravel()], index=z.columns)

weights : PanelData, optional
             Weights to use in demeaning
        """
        if self.nentity > self.nobs:
            group = 'entity'
            dummy = 'time'
        else:
            group = 'time'
            dummy = 'entity'
        e = self.demean(group, weights=weights)
        d = self.dummies(dummy, drop_first=True)
        d.index = e.index
        d = PanelData(d).demean(group, weights=weights)
        d = d.values2d
        e = e.values2d
        resid = e - d @ lstsq(d, e)[0]
        resid = DataFrame(resid, index=self._frame.index,
                          columns=self._frame.columns)

        return PanelData(resid)

if self.other_effects:
            oe = self._other_effect_cats.dataframe
            for c in oe:
                dummies = pd.get_dummies(oe[c], drop_first=drop_first).astype(np.float64)
                d.append(dummies.values)
                drop_first = True

        d = np.column_stack(d)
        wd = root_w * d
        if self.has_constant:
            wd -= root_w * (w.T @ d / w.sum())
            z = np.ones_like(root_w)
            d -= z * (z.T @ d / z.sum())

        x_mean = lstsq(wd, x)[0]
        y_mean = lstsq(wd, y)[0]

        # Save fitted unweighted effects to use in eps calculation
        x_effects = d @ x_mean
        y_effects = d @ y_mean

        # Purge fitted, weighted values
        x = x - wd @ x_mean
        y = y - wd @ y_mean

        ybar = root_w @ lstsq(root_w, y)[0]
        return y, x, ybar, y_effects, x_effects