# TODO: Determine which tests are valid for GLSAR, and under what conditions
# TODO: Fix issue with constant and GLS
# TODO: GLS: add options Iterative GLS, for iterative fgls if sigma is None
# TODO: GLS: default if sigma is none should be two-step GLS
# TODO: Check nesting when performing model based tests, lr, wald, lm
"""
This module implements standard regression models:

Generalized Least Squares (GLS)
Ordinary Least Squares (OLS)
Weighted Least Squares (WLS)
Generalized Least Squares with autoregressive error terms GLSAR(p)

Models are specified with an endogenous response variable and an
exogenous design matrix and are fit using their `fit` method.

Subclasses that have more complicated covariance matrices
should write over the 'whiten' method as the fit method
prewhitens the response by calling 'whiten'.

General reference for regression models:

D. C. Montgomery and E.A. Peck. "Introduction to Linear Regression
    Analysis." 2nd. Ed., Wiley, 1992.

Econometrics references for regression models:

R. Davidson and J.G. MacKinnon.  "Econometric Theory and Methods," Oxford,
    2004.

W. Green.  "Econometric Analysis," 5th ed., Pearson, 2003.
"""
from __future__ import annotations

from statsmodels.compat.pandas import Appender
from statsmodels.compat.python import lrange, lzip

from typing import TYPE_CHECKING, Literal
import warnings

import numpy as np
from scipy import optimize, stats
from scipy.linalg import cholesky, toeplitz
from scipy.linalg.lapack import dtrtri

import statsmodels.base.model as base
import statsmodels.base.wrapper as wrap
from statsmodels.emplike.elregress import _ELRegOpts

# need import in module instead of lazily to copy `__doc__`
from statsmodels.regression._prediction import PredictionResults
from statsmodels.tools.decorators import cache_readonly, cache_writable
from statsmodels.tools.sm_exceptions import (
    InvalidTestWarning,
    SingularMatrixWarning,
    ValueWarning,
)
from statsmodels.tools.tools import pinv_extended
from statsmodels.tools.validation import bool_like, float_like, string_like

from . import _prediction as pred

if TYPE_CHECKING:
    from collections.abc import Sequence

    from statsmodels.tools.typing import Float64Array

__docformat__ = "restructuredtext en"

__all__ = [
    "GLS",
    "GLSAR",
    "OLS",
    "WLS",
    "PredictionResults",
    "RegressionResultsWrapper",
]


_fit_regularized_doc = r"""
        Return a regularized fit to a linear regression model.

        Parameters
        ----------
        method : str
            Either 'elastic_net' or 'sqrt_lasso'.
        alpha : scalar or array_like
            The penalty weight.  If a scalar, the same penalty weight
            applies to all variables in the model.  If a vector, it
            must have the same length as `params`, and contains a
            penalty weight for each coefficient.
        L1_wt : scalar
            The fraction of the penalty given to the L1 penalty term.
            Must be between 0 and 1 (inclusive).  If 0, the fit is a
            ridge fit, if 1 it is a lasso fit.
        start_params : array_like
            Starting values for ``params``.
        profile_scale : bool
            If True the penalized fit is computed using the profile
            (concentrated) log-likelihood for the Gaussian model.
            Otherwise the fit uses the residual sum of squares.
        refit : bool
            If True, the model is refit using only the variables that
            have non-zero coefficients in the regularized fit.  The
            refitted model is not regularized.
        **kwargs
            Additional keyword arguments that contain information used when
            constructing a model using the formula interface.

        Returns
        -------
        statsmodels.base.elastic_net.RegularizedResults
            The regularized results.

        Notes
        -----
        The elastic net uses a combination of L1 and L2 penalties.
        The implementation closely follows the glmnet package in R.

        The function that is minimized is:

        .. math::

            0.5*RSS/n + alpha*((1-L1\_wt)*|params|_2^2/2 + L1\_wt*|params|_1)

        where RSS is the usual regression sum of squares, n is the
        sample size, and :math:`|*|_1` and :math:`|*|_2` are the L1 and L2
        norms.

        For WLS and GLS, the RSS is calculated using the whitened endog and
        exog data.

        Post-estimation results are based on the same data used to
        select variables, hence may be subject to overfitting biases.

        The elastic_net method uses the following keyword arguments:

        maxiter : int
            Maximum number of iterations
        cnvrg_tol : float
            Convergence threshold for line searches
        zero_tol : float
            Coefficients below this threshold are treated as zero.

        The square root lasso approach is a variation of the Lasso
        that is largely self-tuning (the optimal tuning parameter
        does not depend on the standard deviation of the regression
        errors).  If the errors are Gaussian, the tuning parameter
        can be taken to be

        alpha = 1.1 * np.sqrt(n) * norm.ppf(1 - 0.05 / (2 * p))

        where n is the sample size and p is the number of predictors.

        The square root lasso uses the following keyword arguments:

        zero_tol : float
            Coefficients below this threshold are treated as zero.

        The cvxopt module is required to estimate model using the square root
        lasso.

        References
        ----------
        .. [*] Friedman, Hastie, Tibshirani (2008).  Regularization paths for
           generalized linear models via coordinate descent.  Journal of
           Statistical Software 33(1), 1-22 Feb 2010.

        .. [*] A Belloni, V Chernozhukov, L Wang (2011).  Square-root Lasso:
           pivotal recovery of sparse signals via conic programming.
           Biometrika 98(4), 791-806. https://arxiv.org/pdf/1009.5689.pdf
        """


def _get_sigma(sigma, nobs):
    """
    Returns sigma (matrix, nobs by nobs) for GLS and the inverse of its
    Cholesky decomposition.  Handles dimensions and checks integrity.
    If sigma is None, returns None, None. Otherwise returns sigma,
    cholsigmainv.
    """
    if sigma is None:
        return None, None
    sigma = np.asarray(sigma).squeeze()
    if sigma.ndim == 0:
        sigma = np.repeat(sigma, nobs)
    if sigma.ndim == 1:
        if sigma.shape != (nobs,):
            raise ValueError(
                "Sigma must be a scalar, 1d of length %s or a 2d "
                "array of shape %s x %s" % (nobs, nobs, nobs)
            )
        cholsigmainv = 1 / np.sqrt(sigma)
    else:
        if sigma.shape != (nobs, nobs):
            raise ValueError(
                "Sigma must be a scalar, 1d of length %s or a 2d "
                "array of shape %s x %s" % (nobs, nobs, nobs)
            )
        cholsigmainv, info = dtrtri(
            cholesky(sigma, lower=True), lower=True, overwrite_c=True
        )
        if info > 0:
            raise np.linalg.LinAlgError(
                "Cholesky decomposition of sigma yields a singular matrix"
            )
        elif info < 0:
            raise ValueError("Invalid input to dtrtri (info = %d)" % info)
    return sigma, cholsigmainv


class RegressionModel(base.LikelihoodModel):
    """
    Base class for linear regression models. Should not be directly called.

    Intended for subclassing.
    """

    def __init__(self, endog, exog, **kwargs):
        super().__init__(endog, exog, **kwargs)
        self.pinv_wexog: Float64Array | None = None
        self._data_attr.extend(["pinv_wexog", "wendog", "wexog", "weights"])

    def initialize(self):
        """Initialize model components."""
        self.wexog = self.whiten(self.exog)
        self.wendog = self.whiten(self.endog)
        # overwrite nobs from class Model:
        self.nobs = float(self.wexog.shape[0])

        self._df_model = None
        self._df_resid = None
        self.rank = None

    @property
    def df_model(self):
        """
        The model degree of freedom.

        The dof is defined as the rank of the regressor matrix minus 1 if a
        constant is included.
        """
        if self._df_model is None:
            if self.rank is None:
                self.rank = np.linalg.matrix_rank(self.exog)
            self._df_model = float(self.rank - self.k_constant)
        return self._df_model

    @df_model.setter
    def df_model(self, value):
        self._df_model = value

    @property
    def df_resid(self):
        """
        The residual degree of freedom.

        The dof is defined as the number of observations minus the rank of
        the regressor matrix.
        """

        if self._df_resid is None:
            if self.rank is None:
                self.rank = np.linalg.matrix_rank(self.exog)
            self._df_resid = self.nobs - self.rank
        return self._df_resid

    @df_resid.setter
    def df_resid(self, value):
        self._df_resid = value

    def whiten(self, x):
        """
        Whiten method that must be overwritten by individual models.

        Parameters
        ----------
        x : array_like
            Data to be whitened.
        """
        raise NotImplementedError("Subclasses must implement.")

    def fit(
        self,
        method: Literal["pinv", "qr"] = "pinv",
        cov_type: Literal[
            "nonrobust",
            "fixed scale",
            "HC0",
            "HC1",
            "HC2",
            "HC3",
            "HAC",
            "hac-panel",
            "hac-groupsum",
            "cluster",
        ] = "nonrobust",
        cov_kwds=None,
        use_t: bool | None = None,
        **kwargs,
    ):
        """
        Full fit of the model.

        The results include an estimate of covariance matrix, (whitened)
        residuals and an estimate of scale.

        Parameters
        ----------
        method : str, optional
            Can be "pinv", "qr".  "pinv" uses the Moore-Penrose pseudoinverse
            to solve the least squares problem. "qr" uses the QR
            factorization.
        cov_type : str, optional
            See `regression.linear_model.RegressionResults` for a description
            of the available covariance estimators.
        cov_kwds : list or None, optional
            See `linear_model.RegressionResults.get_robustcov_results` for a
            description required keywords for alternative covariance
            estimators.
        use_t : bool, optional
            Flag indicating to use the Student's t distribution when computing
            p-values.  Default behavior depends on cov_type. See
            `linear_model.RegressionResults.get_robustcov_results` for
            implementation details.
        **kwargs
            Additional keyword arguments that contain information used when
            constructing a model using the formula interface.

        Returns
        -------
        RegressionResults
            The model estimation results.

        See Also
        --------
        RegressionResults
            The results container.
        RegressionResults.get_robustcov_results
            A method to change the covariance estimator used when fitting the
            model.

        Notes
        -----
        The fit method uses the pseudoinverse of the design/exogenous variables
        to solve the least squares minimization.
        """
        if method == "pinv":
            if not (
                hasattr(self, "pinv_wexog")
                and hasattr(self, "normalized_cov_params")
                and hasattr(self, "rank")
            ):

                self.pinv_wexog, singular_values = pinv_extended(self.wexog)
                self.normalized_cov_params = np.dot(
                    self.pinv_wexog, np.transpose(self.pinv_wexog)
                )

                # Cache these singular values for use later.
                self.wexog_singular_values = singular_values
                self.rank = np.linalg.matrix_rank(np.diag(singular_values))

            beta = np.dot(self.pinv_wexog, self.wendog)

        elif method == "qr":
            if not (
                hasattr(self, "exog_Q")
                and hasattr(self, "exog_R")
                and hasattr(self, "normalized_cov_params")
                and hasattr(self, "rank")
            ):
                Q, R = np.linalg.qr(self.wexog)
                self.exog_Q, self.exog_R = Q, R
                self.normalized_cov_params = np.linalg.inv(np.dot(R.T, R))

                # Cache singular values from R.
                self.wexog_singular_values = np.linalg.svd(R, 0, 0)
                self.rank = np.linalg.matrix_rank(R)
            else:
                Q, R = self.exog_Q, self.exog_R
            # Needed for some covariance estimators, see GH #8157
            self.pinv_wexog = np.linalg.pinv(self.wexog)
            # used in ANOVA
            self.effects = effects = np.dot(Q.T, self.wendog)
            beta = np.linalg.solve(R, effects)
        else:
            raise ValueError('method has to be "pinv" or "qr"')

        if self._df_model is None:
            self._df_model = float(self.rank - self.k_constant)
        if self._df_resid is None:
            self.df_resid = self.nobs - self.rank

        if isinstance(self, OLS):
            lfit = OLSResults(
                self,
                beta,
                normalized_cov_params=self.normalized_cov_params,
                cov_type=cov_type,
                cov_kwds=cov_kwds,
                use_t=use_t,
            )
        else:
            lfit = RegressionResults(
                self,
                beta,
                normalized_cov_params=self.normalized_cov_params,
                cov_type=cov_type,
                cov_kwds=cov_kwds,
                use_t=use_t,
                **kwargs,
            )
        return RegressionResultsWrapper(lfit)

    def predict(self, params, exog=None):
        """
        Return linear predicted values from a design matrix.

        Parameters
        ----------
        params : array_like
            Parameters of a linear model.
        exog : array_like, optional
            Design / exogenous data. Model exog is used if None.

        Returns
        -------
        array_like
            An array of fitted values.

        Notes
        -----
        If the model has not yet been fit, params is not optional.
        """
        # JP: this does not look correct for GLMAR
        # SS: it needs its own predict method

        if exog is None:
            exog = self.exog

        return np.dot(exog, params)

    def get_distribution(self, params, scale, exog=None, dist_class=None):
        """
        Construct a random number generator for the predictive distribution.

        Parameters
        ----------
        params : array_like
            The model parameters (regression coefficients).
        scale : scalar
            The variance parameter.
        exog : array_like
            The predictor variable matrix.
        dist_class : class
            A random number generator class.  Must take 'loc' and 'scale'
            as arguments and return a random number generator implementing
            an ``rvs`` method for simulating random values. Defaults to normal.

        Returns
        -------
        gen
            Frozen random number generator object with mean and variance
            determined by the fitted linear model.  Use the ``rvs`` method
            to generate random values.

        Notes
        -----
        Due to the behavior of ``scipy.stats.distributions objects``,
        the returned random number generator must be called with
        ``gen.rvs(n)`` where ``n`` is the number of observations in
        the data set used to fit the model.  If any other value is
        used for ``n``, misleading results will be produced.
        """
        fit = self.predict(params, exog)
        if dist_class is None:
            from scipy.stats.distributions import norm

            dist_class = norm
        gen = dist_class(loc=fit, scale=np.sqrt(scale))
        return gen


class GLS(RegressionModel):
    __doc__ = r"""
    Generalized Least Squares

    {params}
    sigma : scalar or array
        The array or scalar `sigma` is the weighting matrix of the covariance.
        The default is None for no scaling.  If `sigma` is a scalar, it is
        assumed that `sigma` is an n x n diagonal matrix with the given
        scalar, `sigma` as the value of each diagonal element.  If `sigma`
        is an n-length vector, then `sigma` is assumed to be a diagonal
        matrix with the given `sigma` on the diagonal.  This should be the
        same as WLS.
    {extra_params}

    Attributes
    ----------
    pinv_wexog : ndarray
        `pinv_wexog` is the p x n Moore-Penrose pseudoinverse of `wexog`.
    cholsimgainv : ndarray
        The transpose of the Cholesky decomposition of the pseudoinverse.
    df_model : float
        p - 1, where p is the number of regressors including the intercept.
        of freedom.
    df_resid : float
        Number of observations n less the number of parameters p.
    llf : float
        The value of the likelihood function of the fitted model.
    nobs : float
        The number of observations n.
    normalized_cov_params : ndarray
        p x p array :math:`(X^{{T}}\Sigma^{{-1}}X)^{{-1}}`
    results : RegressionResults instance
        A property that returns the RegressionResults class if fit.
    sigma : ndarray
        `sigma` is the n x n covariance structure of the error terms.
    wexog : ndarray
        Design matrix whitened by `cholsigmainv`
    wendog : ndarray
        Response variable whitened by `cholsigmainv`

    See Also
    --------
    WLS : Fit a linear model using Weighted Least Squares.
    OLS : Fit a linear model using Ordinary Least Squares.

    Notes
    -----
    If sigma is a function of the data making one of the regressors
    a constant, then the current postestimation statistics will not be correct.

    Examples
    --------
    >>> import statsmodels.api as sm
    >>> data = sm.datasets.longley.load()
    >>> data.exog = sm.add_constant(data.exog)
    >>> ols_resid = sm.OLS(data.endog, data.exog).fit().resid
    >>> res_fit = sm.OLS(ols_resid[1:], ols_resid[:-1]).fit()
    >>> rho = res_fit.params

    `rho` is a consistent estimator of the correlation of the residuals from
    an OLS fit of the longley data.  It is assumed that this is the true rho
    of the AR process data.

    >>> from scipy.linalg import toeplitz
    >>> order = toeplitz(np.arange(16))
    >>> sigma = rho**order

    `sigma` is an n x n matrix of the autocorrelation structure of the
    data.

    >>> gls_model = sm.GLS(data.endog, data.exog, sigma=sigma)
    >>> gls_results = gls_model.fit()
    >>> print(gls_results.summary())
    """.format(
        params=base._model_params_doc,
        extra_params=base._missing_param_doc + base._extra_param_doc,
    )

    def __init__(
        self, endog, exog, sigma=None, missing="none", hasconst=None, **kwargs
    ):
        if type(self) is GLS:
            self._check_kwargs(kwargs)
        # TODO: add options igls, for iterative fgls if sigma is None
        # TODO: default if sigma is none should be two-step GLS
        sigma, cholsigmainv = _get_sigma(sigma, len(endog))

        super().__init__(
            endog,
            exog,
            missing=missing,
            hasconst=hasconst,
            sigma=sigma,
            cholsigmainv=cholsigmainv,
            **kwargs,
        )

        # store attribute names for data arrays
        self._data_attr.extend(["sigma", "cholsigmainv"])

    def whiten(self, x):
        """
        GLS whiten method.

        Parameters
        ----------
        x : array_like
            Data to be whitened.

        Returns
        -------
        ndarray
            The value np.dot(cholsigmainv,X).

        See Also
        --------
        GLS : Fit a linear model using Generalized Least Squares.
        """
        x = np.asarray(x)
        if self.sigma is None or self.sigma.shape == ():
            return x
        elif self.sigma.ndim == 1:
            if x.ndim == 1:
                return x * self.cholsigmainv
            else:
                return x * self.cholsigmainv[:, None]
        else:
            return np.dot(self.cholsigmainv, x)

    def loglike(self, params):
        r"""
        Compute the value of the Gaussian log-likelihood function at params.

        Given the whitened design matrix, the log-likelihood is evaluated
        at the parameter vector `params` for the dependent variable `endog`.

        Parameters
        ----------
        params : array_like
            The model parameters.

        Returns
        -------
        float
            The value of the log-likelihood function for a GLS Model.

        Notes
        -----
        The log-likelihood function for the normal distribution is

        .. math:: -\frac{n}{2}\log\left(\left(Y-\hat{Y}\right)^{\prime}
                   \left(Y-\hat{Y}\right)\right)
                  -\frac{n}{2}\left(1+\log\left(\frac{2\pi}{n}\right)\right)
                  -\frac{1}{2}\log\left(\left|\Sigma\right|\right)

        Y and Y-hat are whitened.
        """
        # TODO: combine this with OLS/WLS loglike and add _det_sigma argument
        nobs2 = self.nobs / 2.0
        SSR = np.sum((self.wendog - np.dot(self.wexog, params)) ** 2, axis=0)
        llf = -np.log(SSR) * nobs2  # concentrated likelihood
        llf -= (1 + np.log(np.pi / nobs2)) * nobs2  # with likelihood constant
        if np.any(self.sigma):
            # FIXME: robust-enough check? unneeded if _det_sigma gets defined
            if self.sigma.ndim == 2:
                det = np.linalg.slogdet(self.sigma)
                llf -= 0.5 * det[1]
            else:
                llf -= 0.5 * np.sum(np.log(self.sigma))
            # with error covariance matrix
        return llf

    def hessian_factor(self, params, scale=None, observed=True):
        """
        Compute weights for calculating Hessian.

        Parameters
        ----------
        params : ndarray
            The parameter at which Hessian is evaluated.
        scale : None or float
            If scale is None, then the default scale will be calculated.
            Default scale is defined by `self.scaletype` and set in fit.
            If scale is not None, then it is used as a fixed scale.
        observed : bool
            If True, then the observed Hessian is returned. If false then the
            expected information matrix is returned.

        Returns
        -------
        ndarray
            A 1d weight vector used in the calculation of the Hessian.
            The hessian is obtained by `(exog.T * hessian_factor).dot(exog)`.
        """

        if self.sigma is None or self.sigma.shape == ():
            return np.ones(self.exog.shape[0])
        elif self.sigma.ndim == 1:
            return self.cholsigmainv
        else:
            return np.diag(self.cholsigmainv)

    @Appender(_fit_regularized_doc)
    def fit_regularized(
        self,
        method="elastic_net",
        alpha=0.0,
        L1_wt=1.0,
        start_params=None,
        profile_scale=False,
        refit=False,
        **kwargs,
    ):
        if not np.isscalar(alpha):
            alpha = np.asarray(alpha)
        # Need to adjust since RSS/n term in elastic net uses nominal
        # n in denominator
        if self.sigma is not None:
            if self.sigma.ndim == 2:
                var_obs = np.diag(self.sigma)
            elif self.sigma.ndim == 1:
                var_obs = self.sigma
            else:
                raise ValueError("sigma should be 1-dim or 2-dim")

            alpha = alpha * np.sum(1 / var_obs) / len(self.endog)

        rslt = OLS(self.wendog, self.wexog).fit_regularized(
            method=method,
            alpha=alpha,
            L1_wt=L1_wt,
            start_params=start_params,
            profile_scale=profile_scale,
            refit=refit,
            **kwargs,
        )

        from statsmodels.base.elastic_net import (
            RegularizedResults,
            RegularizedResultsWrapper,
        )

        rrslt = RegularizedResults(self, rslt.params)
        return RegularizedResultsWrapper(rrslt)


class WLS(RegressionModel):
    __doc__ = """
    Weighted Least Squares

    The weights are presumed to be (proportional to) the inverse of
    the variance of the observations.  That is, if the variables are
    to be transformed by 1/sqrt(W) you must supply weights = 1/W.

    {params}
    weights : array_like, optional
        A 1d array of weights.  If you supply 1/W then the variables are
        pre- multiplied by 1/sqrt(W).  If no weights are supplied the
        default value is 1 and WLS results are the same as OLS.
    {extra_params}

    Attributes
    ----------
    weights : ndarray
        The stored weights supplied as an argument.

    See Also
    --------
    GLS : Fit a linear model using Generalized Least Squares.
    OLS : Fit a linear model using Ordinary Least Squares.

    Notes
    -----
    If the weights are a function of the data, then the post estimation
    statistics such as fvalue and mse_model might not be correct, as the
    package does not yet support no-constant regression.

    Examples
    --------
    >>> import statsmodels.api as sm
    >>> Y = [1,3,4,5,2,3,4]
    >>> X = range(1,8)
    >>> X = sm.add_constant(X)
    >>> wls_model = sm.WLS(Y,X, weights=list(range(1,8)))
    >>> results = wls_model.fit()
    >>> results.params
    array([ 2.91666667,  0.0952381 ])
    >>> results.tvalues
    array([ 2.0652652 ,  0.35684428])
    >>> print(results.t_test([1, 0]))
    <T test: effect=array([ 2.91666667]), sd=array([[ 1.41224801]]),
     t=array([[ 2.0652652]]), p=array([[ 0.04690139]]), df_denom=5>
    >>> print(results.f_test([0, 1]))
    <F test: F=array([[ 0.12733784]]), p=[[ 0.73577409]], df_denom=5, df_num=1>
    """.format(
        params=base._model_params_doc,
        extra_params=base._missing_param_doc + base._extra_param_doc,
    )

    def __init__(
        self, endog, exog, weights=1.0, missing="none", hasconst=None, **kwargs
    ):
        if type(self) is WLS:
            self._check_kwargs(kwargs)
        weights = np.array(weights)
        if weights.shape == ():
            if (
                missing == "drop"
                and "missing_idx" in kwargs
                and kwargs["missing_idx"] is not None
            ):
                # patsy may have truncated endog
                weights = np.repeat(weights, len(kwargs["missing_idx"]))
            else:
                weights = np.repeat(weights, len(endog))
        # handle case that endog might be of len == 1
        if len(weights) == 1:
            weights = np.array([weights.squeeze()])
        else:
            weights = weights.squeeze()
        super().__init__(
            endog, exog, missing=missing, weights=weights, hasconst=hasconst, **kwargs
        )
        nobs = self.exog.shape[0]
        weights = self.weights
        if weights.size != nobs and weights.shape[0] != nobs:
            raise ValueError("Weights must be scalar or same length as design")

    def whiten(self, x):
        """
        Whitener for WLS model, multiplies each column by sqrt(self.weights).

        Parameters
        ----------
        x : array_like
            Data to be whitened.

        Returns
        -------
        array_like
            The whitened values sqrt(weights)*X.
        """

        x = np.asarray(x)
        if x.ndim == 1:
            return x * np.sqrt(self.weights)
        elif x.ndim == 2:
            return np.sqrt(self.weights)[:, None] * x

    def loglike(self, params):
        r"""
        Compute the value of the gaussian log-likelihood function at params.

        Given the whitened design matrix, the log-likelihood is evaluated
        at the parameter vector `params` for the dependent variable `Y`.

        Parameters
        ----------
        params : array_like
            The parameter estimates.

        Returns
        -------
        float
            The value of the log-likelihood function for a WLS Model.

        Notes
        -----
        .. math:: -\frac{n}{2}\log SSR
                  -\frac{n}{2}\left(1+\log\left(\frac{2\pi}{n}\right)\right)
                  +\frac{1}{2}\log\left(\left|W\right|\right)

        where :math:`W` is a diagonal weight matrix matrix,
        :math:`\left|W\right|` is its determinant, and
        :math:`SSR=\left(Y-\hat{Y}\right)^\prime W \left(Y-\hat{Y}\right)` is
        the sum of the squared weighted residuals.
        """
        nobs2 = self.nobs / 2.0
        SSR = np.sum((self.wendog - np.dot(self.wexog, params)) ** 2, axis=0)
        llf = -np.log(SSR) * nobs2  # concentrated likelihood
        llf -= (1 + np.log(np.pi / nobs2)) * nobs2  # with constant
        llf += 0.5 * np.sum(np.log(self.weights))
        return llf

    def hessian_factor(self, params, scale=None, observed=True):
        """
        Compute the weights for calculating the Hessian.

        Parameters
        ----------
        params : ndarray
            The parameter at which Hessian is evaluated.
        scale : None or float
            If scale is None, then the default scale will be calculated.
            Default scale is defined by `self.scaletype` and set in fit.
            If scale is not None, then it is used as a fixed scale.
        observed : bool
            If True, then the observed Hessian is returned. If false then the
            expected information matrix is returned.

        Returns
        -------
        ndarray
            A 1d weight vector used in the calculation of the Hessian.
            The hessian is obtained by `(exog.T * hessian_factor).dot(exog)`.
        """

        return self.weights

    @Appender(_fit_regularized_doc)
    def fit_regularized(
        self,
        method="elastic_net",
        alpha=0.0,
        L1_wt=1.0,
        start_params=None,
        profile_scale=False,
        refit=False,
        **kwargs,
    ):
        # Docstring attached below
        if not np.isscalar(alpha):
            alpha = np.asarray(alpha)
        # Need to adjust since RSS/n in elastic net uses nominal n in
        # denominator
        alpha = alpha * np.sum(self.weights) / len(self.weights)

        rslt = OLS(self.wendog, self.wexog).fit_regularized(
            method=method,
            alpha=alpha,
            L1_wt=L1_wt,
            start_params=start_params,
            profile_scale=profile_scale,
            refit=refit,
            **kwargs,
        )

        from statsmodels.base.elastic_net import (
            RegularizedResults,
            RegularizedResultsWrapper,
        )

        rrslt = RegularizedResults(self, rslt.params)
        return RegularizedResultsWrapper(rrslt)


class OLS(WLS):
    __doc__ = """
    Ordinary Least Squares

    {params}
    {extra_params}

    Attributes
    ----------
    weights : scalar
        Has an attribute weights = array(1.0) due to inheritance from WLS.

    See Also
    --------
    WLS : Fit a linear model using Weighted Least Squares.
    GLS : Fit a linear model using Generalized Least Squares.

    Notes
    -----
    No constant is added by the model unless you are using formulas.

    Examples
    --------
    >>> import statsmodels.api as sm
    >>> import numpy as np
    >>> duncan_prestige = sm.datasets.get_rdataset("Duncan", "carData")
    >>> Y = duncan_prestige.data['income']
    >>> X = duncan_prestige.data['education']
    >>> X = sm.add_constant(X)
    >>> model = sm.OLS(Y,X)
    >>> results = model.fit()
    >>> results.params
    const        10.603498
    education     0.594859
    dtype: float64

    >>> results.tvalues
    const        2.039813
    education    6.892802
    dtype: float64

    >>> print(results.t_test([1, 0]))
                                 Test for Constraints
    ==============================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
    ------------------------------------------------------------------------------
    c0            10.6035      5.198      2.040      0.048       0.120      21.087
    ==============================================================================

    >>> print(results.f_test(np.identity(2)))
