MultipleRegression.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "MultipleRegression.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyMdxMSzd3yuNf50TfwKnH/f",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/pragneshrana/Algo_Trial/blob/master/MultipleRegression.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bJgvWBcLLjY2",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "outputId": "d99eca0b-c60a-408c-e290-ed6b4ca2e718"
      },
      "source": [
        "'''\n",
        "Code to extract the OLS summary of different parameter.\n",
        "TO save in excel file\n",
        "RUN model and append the result\n",
        "'''\n",
        "import statsmodels.api as sm\n",
        "import numpy as np\n",
        "\n",
        "Y = [1,3,4,5,2,3,4]\n",
        "X = range(1,8)\n",
        "X = sm.add_constant(X)\n",
        "model = sm.OLS(Y,X)\n",
        "results = model.fit()\n",
        "print('Summary:',results.summary())\n",
        "# print('params:',results.params)\n",
        "# print('tvalues:',results.tvalues)\n",
        "# print('ttest:',results.t_test([1, 0]))\n",
        "# print('ftest:',results.f_test(np.identity(2)))\n",
        "print('Conition Number:',results.condition_number)\n",
        "print('Eigen values:',results.eigenvals)\n",
        "print('Rsquared:',results.rsquared)\n",
        "print('F-statistic:', results.fvalue)\n",
        "print('Prob (F-statistic):', results.f_pvalue)\n",
        "print('Log-Likelihood:', None)\n",
        "print('AIC:', results.aic)\n",
        "print('BIC:', results.bic)\n",
        "print('SSR:',results.ssr)\n",
        "print('Confidence intervasl:',results.conf_int())\n",
        "print('Number of observations :',results.nobs)\n",
        "print('Predicted Values:',results.fittedvalues)\n",
        "print('Residual:',results.resid)\n",
        "print('Scale:',results.scale)\n",
        "print('Centered_tss:',results.centered_tss)\n",
        "print('ESS:',results.ess)\n",
        "print('Rsquared_adj:',results.rsquared_adj)\n",
        "print('MSE Model:',results.mse_model)\n",
        "print('MSE_residuals:',results.mse_resid)\n",
        "print('MSE Total:',results.mse_total)\n",
        "print(' Residuals Pearson:',results.resid_pearson)\n",
        "print('Residual Transformed:',)\n",
        "from statsmodels.stats.stattools import (jarque_bera, omni_normtest, durbin_watson)\n",
        "wresid_result = results.wresid\n",
        "jb, jbpv, skew, kurtosis = jarque_bera(wresid_result)\n",
        "omni, omnipv = omni_normtest(wresid_result)\n",
        "dw = durbin_watson(wresid_result)\n",
        "print('Omnibus:',omni)\n",
        "print('Prob(Omnibus):',omnipv)\n",
        "print('Skewnes:',skew)\n",
        "print('Kurtosis:',kurtosis)\n",
        "print('Jarque-Bera:',jb)\n",
        "print('JB Probability:',jbpv)\n",
        "print('Durbin-Watson:',dw)\n",
        "\n",
        "\n"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Summary:                             OLS Regression Results                            \n",
            "==============================================================================\n",
            "Dep. Variable:                      y   R-squared:                       0.161\n",
            "Model:                            OLS   Adj. R-squared:                 -0.007\n",
            "Method:                 Least Squares   F-statistic:                    0.9608\n",
            "Date:                Tue, 28 Jan 2020   Prob (F-statistic):              0.372\n",
            "Time:                        14:25:56   Log-Likelihood:                -10.854\n",
            "No. Observations:                   7   AIC:                             25.71\n",
            "Df Residuals:                       5   BIC:                             25.60\n",
            "Df Model:                           1                                         \n",
            "Covariance Type:            nonrobust                                         \n",
            "==============================================================================\n",
            "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
            "------------------------------------------------------------------------------\n",
            "const          2.1429      1.141      1.879      0.119      -0.789       5.075\n",
            "x1             0.2500      0.255      0.980      0.372      -0.406       0.906\n",
            "==============================================================================\n",
            "Omnibus:                          nan   Durbin-Watson:                   1.743\n",
            "Prob(Omnibus):                    nan   Jarque-Bera (JB):                0.482\n",
            "Skew:                           0.206   Prob(JB):                        0.786\n",
            "Kurtosis:                       1.782   Cond. No.                         10.4\n",
            "==============================================================================\n",
            "\n",
            "Warnings:\n",
            "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
            "Conition Number: 10.40388203202207\n",
            "Eigen values: [145.65434845   1.34565155]\n",
            "Rsquared: 0.16118421052631593\n",
            "F-statistic: 0.9607843137254913\n",
            "Prob (F-statistic): 0.37200486130017757\n",
            "Log-Likelihood: None\n",
            "AIC: 25.70718166636078\n",
            "BIC: 25.599001964471405\n",
            "SSR: 9.107142857142856\n",
            "Confidence intervasl: [[-0.78920715  5.07492144]\n",
            " [-0.40562951  0.90562951]]\n",
            "Number of observations : 7.0\n",
            "Predicted Values: [2.39285714 2.64285714 2.89285714 3.14285714 3.39285714 3.64285714\n",
            " 3.89285714]\n",
            "Residual: [-1.39285714  0.35714286  1.10714286  1.85714286 -1.39285714 -0.64285714\n",
            "  0.10714286]\n",
            "Scale: 1.8214285714285712\n",
            "Centered_tss: 10.857142857142858\n",
            "ESS: 1.7500000000000018\n",
            "Rsquared_adj: -0.0065789473684207955\n",
            "MSE Model: 1.7500000000000018\n",
            "MSE_residuals: 1.8214285714285712\n",
            "MSE Total: 1.8095238095238095\n",
            " Residuals Pearson: [-1.03204944  0.26462806  0.82034699  1.37606592 -1.03204944 -0.47633051\n",
            "  0.07938842]\n",
            "Residual Transformed:\n",
            "Omnibus: nan\n",
            "Prob(Omnibus): nan\n",
            "Skewnes: 0.20591728159799946\n",
            "Kurtosis: 1.7818800461361015\n",
            "Jarque-Bera: 0.48224864608789914\n",
            "JB Probability: 0.7857439342381092\n",
            "Durbin-Watson: 1.7431372549019608\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/statsmodels/stats/stattools.py:71: ValueWarning: omni_normtest is not valid with less than 8 observations; 7 samples were given.\n",
            "  \"samples were given.\" % int(n), ValueWarning)\n",
            "/usr/local/lib/python3.6/dist-packages/statsmodels/stats/stattools.py:71: ValueWarning: omni_normtest is not valid with less than 8 observations; 7 samples were given.\n",
            "  \"samples were given.\" % int(n), ValueWarning)\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ej10OxgKMSe9",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "h8m-X_TPln4D",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 386
        },
        "outputId": "f9d91f90-089b-4a8b-f776-1f0264b66ec3"
      },
      "source": [
        "##stat Library \n",
        "# TODO: Determine which tests are valid for GLSAR, and under what conditions\n",
        "# TODO: Fix issue with constant and GLS\n",
        "# TODO: GLS: add options Iterative GLS, for iterative fgls if sigma is None\n",
        "# TODO: GLS: default if sigma is none should be two-step GLS\n",
        "# TODO: Check nesting when performing model based tests, lr, wald, lm\n",
        "\"\"\"\n",
        "This module implements standard regression models:\n",
        "\n",
        "Generalized Least Squares (GLS)\n",
        "Ordinary Least Squares (OLS)\n",
        "Weighted Least Squares (WLS)\n",
        "Generalized Least Squares with autoregressive error terms GLSAR(p)\n",
        "\n",
        "Models are specified with an endogenous response variable and an\n",
        "exogenous design matrix and are fit using their `fit` method.\n",
        "\n",
        "Subclasses that have more complicated covariance matrices\n",
        "should write over the 'whiten' method as the fit method\n",
        "prewhitens the response by calling 'whiten'.\n",
        "\n",
        "General reference for regression models:\n",
        "\n",
        "D. C. Montgomery and E.A. Peck. \"Introduction to Linear Regression\n",
        "    Analysis.\" 2nd. Ed., Wiley, 1992.\n",
        "\n",
        "Econometrics references for regression models:\n",
        "\n",
        "R. Davidson and J.G. MacKinnon.  \"Econometric Theory and Methods,\" Oxford,\n",
        "    2004.\n",
        "\n",
        "W. Green.  \"Econometric Analysis,\" 5th ed., Pearson, 2003.\n",
        "\"\"\"\n",
        "\n",
        "\n",
        "from statsmodels.compat.python import lrange, lzip\n",
        "from statsmodels.compat.pandas import Appender\n",
        "\n",
        "import numpy as np\n",
        "from scipy.linalg import toeplitz\n",
        "from scipy import stats\n",
        "from scipy import optimize\n",
        "\n",
        "from statsmodels.tools.tools import pinv_extended\n",
        "from statsmodels.tools.decorators import (cache_readonly,\n",
        "                                          cache_writable)\n",
        "import statsmodels.base.model as base\n",
        "import statsmodels.base.wrapper as wrap\n",
        "from statsmodels.emplike.elregress import _ELRegOpts\n",
        "import warnings\n",
        "from statsmodels.tools.sm_exceptions import InvalidTestWarning\n",
        "\n",
        "# need import in module instead of lazily to copy `__doc__`\n",
        "from statsmodels.regression._prediction import PredictionResults\n",
        "from . import _prediction as pred\n",
        "\n",
        "__docformat__ = 'restructuredtext en'\n",
        "\n",
        "__all__ = ['GLS', 'WLS', 'OLS', 'GLSAR', 'PredictionResults',\n",
        "           'RegressionResultsWrapper']\n",
        "\n",
        "\n",
        "_fit_regularized_doc =\\\n",
        "        r\"\"\"\n",
        "        Return a regularized fit to a linear regression model.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        method : str\n",
        "            Either 'elastic_net' or 'sqrt_lasso'.\n",
        "        alpha : scalar or array_like\n",
        "            The penalty weight.  If a scalar, the same penalty weight\n",
        "            applies to all variables in the model.  If a vector, it\n",
        "            must have the same length as `params`, and contains a\n",
        "            penalty weight for each coefficient.\n",
        "        L1_wt : scalar\n",
        "            The fraction of the penalty given to the L1 penalty term.\n",
        "            Must be between 0 and 1 (inclusive).  If 0, the fit is a\n",
        "            ridge fit, if 1 it is a lasso fit.\n",
        "        start_params : array_like\n",
        "            Starting values for ``params``.\n",
        "        profile_scale : bool\n",
        "            If True the penalized fit is computed using the profile\n",
        "            (concentrated) log-likelihood for the Gaussian model.\n",
        "            Otherwise the fit uses the residual sum of squares.\n",
        "        refit : bool\n",
        "            If True, the model is refit using only the variables that\n",
        "            have non-zero coefficients in the regularized fit.  The\n",
        "            refitted model is not regularized.\n",
        "        **kwargs\n",
        "            Additional keyword arguments that contain information used when\n",
        "            constructing a model using the formula interface.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        statsmodels.base.elastic_net.RegularizedResults\n",
        "            The regularized results.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        The elastic net uses a combination of L1 and L2 penalties.\n",
        "        The implementation closely follows the glmnet package in R.\n",
        "\n",
        "        The function that is minimized is:\n",
        "\n",
        "        .. math::\n",
        "\n",
        "            0.5*RSS/n + alpha*((1-L1\\_wt)*|params|_2^2/2 + L1\\_wt*|params|_1)\n",
        "\n",
        "        where RSS is the usual regression sum of squares, n is the\n",
        "        sample size, and :math:`|*|_1` and :math:`|*|_2` are the L1 and L2\n",
        "        norms.\n",
        "\n",
        "        For WLS and GLS, the RSS is calculated using the whitened endog and\n",
        "        exog data.\n",
        "\n",
        "        Post-estimation results are based on the same data used to\n",
        "        select variables, hence may be subject to overfitting biases.\n",
        "\n",
        "        The elastic_net method uses the following keyword arguments:\n",
        "\n",
        "        maxiter : int\n",
        "            Maximum number of iterations\n",
        "        cnvrg_tol : float\n",
        "            Convergence threshold for line searches\n",
        "        zero_tol : float\n",
        "            Coefficients below this threshold are treated as zero.\n",
        "\n",
        "        The square root lasso approach is a variation of the Lasso\n",
        "        that is largely self-tuning (the optimal tuning parameter\n",
        "        does not depend on the standard deviation of the regression\n",
        "        errors).  If the errors are Gaussian, the tuning parameter\n",
        "        can be taken to be\n",
        "\n",
        "        alpha = 1.1 * np.sqrt(n) * norm.ppf(1 - 0.05 / (2 * p))\n",
        "\n",
        "        where n is the sample size and p is the number of predictors.\n",
        "\n",
        "        The square root lasso uses the following keyword arguments:\n",
        "\n",
        "        zero_tol : float\n",
        "            Coefficients below this threshold are treated as zero.\n",
        "\n",
        "        The cvxopt module is required to estimate model using the square root\n",
        "        lasso.\n",
        "\n",
        "        References\n",
        "        ----------\n",
        "        .. [*] Friedman, Hastie, Tibshirani (2008).  Regularization paths for\n",
        "           generalized linear models via coordinate descent.  Journal of\n",
        "           Statistical Software 33(1), 1-22 Feb 2010.\n",
        "\n",
        "        .. [*] A Belloni, V Chernozhukov, L Wang (2011).  Square-root Lasso:\n",
        "           pivotal recovery of sparse signals via conic programming.\n",
        "           Biometrika 98(4), 791-806. https://arxiv.org/pdf/1009.5689.pdf\n",
        "        \"\"\"\n",
        "\n",
        "\n",
        "def _get_sigma(sigma, nobs):\n",
        "    \"\"\"\n",
        "    Returns sigma (matrix, nobs by nobs) for GLS and the inverse of its\n",
        "    Cholesky decomposition.  Handles dimensions and checks integrity.\n",
        "    If sigma is None, returns None, None. Otherwise returns sigma,\n",
        "    cholsigmainv.\n",
        "    \"\"\"\n",
        "    if sigma is None:\n",
        "        return None, None\n",
        "    sigma = np.asarray(sigma).squeeze()\n",
        "    if sigma.ndim == 0:\n",
        "        sigma = np.repeat(sigma, nobs)\n",
        "    if sigma.ndim == 1:\n",
        "        if sigma.shape != (nobs,):\n",
        "            raise ValueError(\"Sigma must be a scalar, 1d of length %s or a 2d \"\n",
        "                             \"array of shape %s x %s\" % (nobs, nobs, nobs))\n",
        "        cholsigmainv = 1/np.sqrt(sigma)\n",
        "    else:\n",
        "        if sigma.shape != (nobs, nobs):\n",
        "            raise ValueError(\"Sigma must be a scalar, 1d of length %s or a 2d \"\n",
        "                             \"array of shape %s x %s\" % (nobs, nobs, nobs))\n",
        "        cholsigmainv = np.linalg.cholesky(np.linalg.inv(sigma)).T\n",
        "    return sigma, cholsigmainv\n",
        "\n",
        "\n",
        "class RegressionModel(base.LikelihoodModel):\n",
        "    \"\"\"\n",
        "    Base class for linear regression models. Should not be directly called.\n",
        "\n",
        "    Intended for subclassing.\n",
        "    \"\"\"\n",
        "    def __init__(self, endog, exog, **kwargs):\n",
        "        super(RegressionModel, self).__init__(endog, exog, **kwargs)\n",
        "        self._data_attr.extend(['pinv_wexog', 'wendog', 'wexog', 'weights'])\n",
        "\n",
        "    def initialize(self):\n",
        "        \"\"\"Initialize model components.\"\"\"\n",
        "        self.wexog = self.whiten(self.exog)\n",
        "        self.wendog = self.whiten(self.endog)\n",
        "        # overwrite nobs from class Model:\n",
        "        self.nobs = float(self.wexog.shape[0])\n",
        "\n",
        "        self._df_model = None\n",
        "        self._df_resid = None\n",
        "        self.rank = None\n",
        "\n",
        "    @property\n",
        "    def df_model(self):\n",
        "        \"\"\"\n",
        "        The model degree of freedom.\n",
        "\n",
        "        The dof is defined as the rank of the regressor matrix minus 1 if a\n",
        "        constant is included.\n",
        "        \"\"\"\n",
        "        if self._df_model is None:\n",
        "            if self.rank is None:\n",
        "                self.rank = np.linalg.matrix_rank(self.exog)\n",
        "            self._df_model = float(self.rank - self.k_constant)\n",
        "        return self._df_model\n",
        "\n",
        "    @df_model.setter\n",
        "    def df_model(self, value):\n",
        "        self._df_model = value\n",
        "\n",
        "    @property\n",
        "    def df_resid(self):\n",
        "        \"\"\"\n",
        "        The residual degree of freedom.\n",
        "\n",
        "        The dof is defined as the number of observations minus the rank of\n",
        "        the regressor matrix.\n",
        "        \"\"\"\n",
        "\n",
        "        if self._df_resid is None:\n",
        "            if self.rank is None:\n",
        "                self.rank = np.linalg.matrix_rank(self.exog)\n",
        "            self._df_resid = self.nobs - self.rank\n",
        "        return self._df_resid\n",
        "\n",
        "    @df_resid.setter\n",
        "    def df_resid(self, value):\n",
        "        self._df_resid = value\n",
        "\n",
        "    def whiten(self, x):\n",
        "        \"\"\"\n",
        "        Whiten method that must be overwritten by individual models.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        x : array_like\n",
        "            Data to be whitened.\n",
        "        \"\"\"\n",
        "        raise NotImplementedError(\"Subclasses must implement.\")\n",
        "\n",
        "    def fit(self, method=\"pinv\", cov_type='nonrobust', cov_kwds=None,\n",
        "            use_t=None, **kwargs):\n",
        "        \"\"\"\n",
        "        Full fit of the model.\n",
        "\n",
        "        The results include an estimate of covariance matrix, (whitened)\n",
        "        residuals and an estimate of scale.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        method : str, optional\n",
        "            Can be \"pinv\", \"qr\".  \"pinv\" uses the Moore-Penrose pseudoinverse\n",
        "            to solve the least squares problem. \"qr\" uses the QR\n",
        "            factorization.\n",
        "        cov_type : str, optional\n",
        "            See `regression.linear_model.RegressionResults` for a description\n",
        "            of the available covariance estimators.\n",
        "        cov_kwds : list or None, optional\n",
        "            See `linear_model.RegressionResults.get_robustcov_results` for a\n",
        "            description required keywords for alternative covariance\n",
        "            estimators.\n",
        "        use_t : bool, optional\n",
        "            Flag indicating to use the Student's t distribution when computing\n",
        "            p-values.  Default behavior depends on cov_type. See\n",
        "            `linear_model.RegressionResults.get_robustcov_results` for\n",
        "            implementation details.\n",
        "        **kwargs\n",
        "            Additional keyword arguments that contain information used when\n",
        "            constructing a model using the formula interface.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        RegressionResults\n",
        "            The model estimation results.\n",
        "\n",
        "        See Also\n",
        "        --------\n",
        "        RegressionResults\n",
        "            The results container.\n",
        "        RegressionResults.get_robustcov_results\n",
        "            A method to change the covariance estimator used when fitting the\n",
        "            model.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        The fit method uses the pseudoinverse of the design/exogenous variables\n",
        "        to solve the least squares minimization.\n",
        "        \"\"\"\n",
        "        if method == \"pinv\":\n",
        "            if not (hasattr(self, 'pinv_wexog') and\n",
        "                    hasattr(self, 'normalized_cov_params') and\n",
        "                    hasattr(self, 'rank')):\n",
        "\n",
        "                self.pinv_wexog, singular_values = pinv_extended(self.wexog)\n",
        "                self.normalized_cov_params = np.dot(\n",
        "                    self.pinv_wexog, np.transpose(self.pinv_wexog))\n",
        "\n",
        "                # Cache these singular values for use later.\n",
        "                self.wexog_singular_values = singular_values\n",
        "                self.rank = np.linalg.matrix_rank(np.diag(singular_values))\n",
        "\n",
        "            beta = np.dot(self.pinv_wexog, self.wendog)\n",
        "\n",
        "        elif method == \"qr\":\n",
        "            if not (hasattr(self, 'exog_Q') and\n",
        "                    hasattr(self, 'exog_R') and\n",
        "                    hasattr(self, 'normalized_cov_params') and\n",
        "                    hasattr(self, 'rank')):\n",
        "                Q, R = np.linalg.qr(self.wexog)\n",
        "                self.exog_Q, self.exog_R = Q, R\n",
        "                self.normalized_cov_params = np.linalg.inv(np.dot(R.T, R))\n",
        "\n",
        "                # Cache singular values from R.\n",
        "                self.wexog_singular_values = np.linalg.svd(R, 0, 0)\n",
        "                self.rank = np.linalg.matrix_rank(R)\n",
        "            else:\n",
        "                Q, R = self.exog_Q, self.exog_R\n",
        "\n",
        "            # used in ANOVA\n",
        "            self.effects = effects = np.dot(Q.T, self.wendog)\n",
        "            beta = np.linalg.solve(R, effects)\n",
        "        else:\n",
        "            raise ValueError('method has to be \"pinv\" or \"qr\"')\n",
        "\n",
        "        if self._df_model is None:\n",
        "            self._df_model = float(self.rank - self.k_constant)\n",
        "        if self._df_resid is None:\n",
        "            self.df_resid = self.nobs - self.rank\n",
        "\n",
        "        if isinstance(self, OLS):\n",
        "            lfit = OLSResults(\n",
        "                self, beta,\n",
        "                normalized_cov_params=self.normalized_cov_params,\n",
        "                cov_type=cov_type, cov_kwds=cov_kwds, use_t=use_t)\n",
        "        else:\n",
        "            lfit = RegressionResults(\n",
        "                self, beta,\n",
        "                normalized_cov_params=self.normalized_cov_params,\n",
        "                cov_type=cov_type, cov_kwds=cov_kwds, use_t=use_t,\n",
        "                **kwargs)\n",
        "        return RegressionResultsWrapper(lfit)\n",
        "\n",
        "    def predict(self, params, exog=None):\n",
        "        \"\"\"\n",
        "        Return linear predicted values from a design matrix.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        params : array_like\n",
        "            Parameters of a linear model.\n",
        "        exog : array_like, optional\n",
        "            Design / exogenous data. Model exog is used if None.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        array_like\n",
        "            An array of fitted values.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        If the model has not yet been fit, params is not optional.\n",
        "        \"\"\"\n",
        "        # JP: this does not look correct for GLMAR\n",
        "        # SS: it needs its own predict method\n",
        "\n",
        "        if exog is None:\n",
        "            exog = self.exog\n",
        "\n",
        "        return np.dot(exog, params)\n",
        "\n",
        "    def get_distribution(self, params, scale, exog=None, dist_class=None):\n",
        "        \"\"\"\n",
        "        Construct a random number generator for the predictive distribution.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        params : array_like\n",
        "            The model parameters (regression coefficients).\n",
        "        scale : scalar\n",
        "            The variance parameter.\n",
        "        exog : array_like\n",
        "            The predictor variable matrix.\n",
        "        dist_class : class\n",
        "            A random number generator class.  Must take 'loc' and 'scale'\n",
        "            as arguments and return a random number generator implementing\n",
        "            an ``rvs`` method for simulating random values. Defaults to normal.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        gen\n",
        "            Frozen random number generator object with mean and variance\n",
        "            determined by the fitted linear model.  Use the ``rvs`` method\n",
        "            to generate random values.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        Due to the behavior of ``scipy.stats.distributions objects``,\n",
        "        the returned random number generator must be called with\n",
        "        ``gen.rvs(n)`` where ``n`` is the number of observations in\n",
        "        the data set used to fit the model.  If any other value is\n",
        "        used for ``n``, misleading results will be produced.\n",
        "        \"\"\"\n",
        "        fit = self.predict(params, exog)\n",
        "        if dist_class is None:\n",
        "            from scipy.stats.distributions import norm\n",
        "            dist_class = norm\n",
        "        gen = dist_class(loc=fit, scale=np.sqrt(scale))\n",
        "        return gen\n",
        "\n",
        "\n",
        "class GLS(RegressionModel):\n",
        "    __doc__ = r\"\"\"\n",
        "    Generalized Least Squares\n",
        "\n",
        "    %(params)s\n",
        "    sigma : scalar or array\n",
        "        The array or scalar `sigma` is the weighting matrix of the covariance.\n",
        "        The default is None for no scaling.  If `sigma` is a scalar, it is\n",
        "        assumed that `sigma` is an n x n diagonal matrix with the given\n",
        "        scalar, `sigma` as the value of each diagonal element.  If `sigma`\n",
        "        is an n-length vector, then `sigma` is assumed to be a diagonal\n",
        "        matrix with the given `sigma` on the diagonal.  This should be the\n",
        "        same as WLS.\n",
        "    %(extra_params)s\n",
        "\n",
        "    Attributes\n",
        "    ----------\n",
        "    pinv_wexog : ndarray\n",
        "        `pinv_wexog` is the p x n Moore-Penrose pseudoinverse of `wexog`.\n",
        "    cholsimgainv : ndarray\n",
        "        The transpose of the Cholesky decomposition of the pseudoinverse.\n",
        "    df_model : float\n",
        "        p - 1, where p is the number of regressors including the intercept.\n",
        "        of freedom.\n",
        "    df_resid : float\n",
        "        Number of observations n less the number of parameters p.\n",
        "    llf : float\n",
        "        The value of the likelihood function of the fitted model.\n",
        "    nobs : float\n",
        "        The number of observations n.\n",
        "    normalized_cov_params : ndarray\n",
        "        p x p array :math:`(X^{T}\\Sigma^{-1}X)^{-1}`\n",
        "    results : RegressionResults instance\n",
        "        A property that returns the RegressionResults class if fit.\n",
        "    sigma : ndarray\n",
        "        `sigma` is the n x n covariance structure of the error terms.\n",
        "    wexog : ndarray\n",
        "        Design matrix whitened by `cholsigmainv`\n",
        "    wendog : ndarray\n",
        "        Response variable whitened by `cholsigmainv`\n",
        "\n",
        "    See Also\n",
        "    --------\n",
        "    WLS : Fit a linear model using Weighted Least Squares.\n",
        "    OLS : Fit a linear model using Ordinary Least Squares.\n",
        "\n",
        "    Notes\n",
        "    -----\n",
        "    If sigma is a function of the data making one of the regressors\n",
        "    a constant, then the current postestimation statistics will not be correct.\n",
        "\n",
        "    Examples\n",
        "    --------\n",
        "    >>> import statsmodels.api as sm\n",
        "    >>> data = sm.datasets.longley.load(as_pandas=False)\n",
        "    >>> data.exog = sm.add_constant(data.exog)\n",
        "    >>> ols_resid = sm.OLS(data.endog, data.exog).fit().resid\n",
        "    >>> res_fit = sm.OLS(ols_resid[1:], ols_resid[:-1]).fit()\n",
        "    >>> rho = res_fit.params\n",
        "\n",
        "    `rho` is a consistent estimator of the correlation of the residuals from\n",
        "    an OLS fit of the longley data.  It is assumed that this is the true rho\n",
        "    of the AR process data.\n",
        "\n",
        "    >>> from scipy.linalg import toeplitz\n",
        "    >>> order = toeplitz(np.arange(16))\n",
        "    >>> sigma = rho**order\n",
        "\n",
        "    `sigma` is an n x n matrix of the autocorrelation structure of the\n",
        "    data.\n",
        "\n",
        "    >>> gls_model = sm.GLS(data.endog, data.exog, sigma=sigma)\n",
        "    >>> gls_results = gls_model.fit()\n",
        "    >>> print(gls_results.summary())\n",
        "    \"\"\" % {'params': base._model_params_doc,\n",
        "           'extra_params': base._missing_param_doc + base._extra_param_doc}\n",
        "\n",
        "    def __init__(self, endog, exog, sigma=None, missing='none', hasconst=None,\n",
        "                 **kwargs):\n",
        "        # TODO: add options igls, for iterative fgls if sigma is None\n",
        "        # TODO: default if sigma is none should be two-step GLS\n",
        "        sigma, cholsigmainv = _get_sigma(sigma, len(endog))\n",
        "\n",
        "        super(GLS, self).__init__(endog, exog, missing=missing,\n",
        "                                  hasconst=hasconst, sigma=sigma,\n",
        "                                  cholsigmainv=cholsigmainv, **kwargs)\n",
        "\n",
        "        # store attribute names for data arrays\n",
        "        self._data_attr.extend(['sigma', 'cholsigmainv'])\n",
        "\n",
        "    def whiten(self, x):\n",
        "        \"\"\"\n",
        "        GLS whiten method.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        x : array_like\n",
        "            Data to be whitened.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        ndarray\n",
        "            The value np.dot(cholsigmainv,X).\n",
        "\n",
        "        See Also\n",
        "        --------\n",
        "        GLS : Fit a linear model using Generalized Least Squares.\n",
        "        \"\"\"\n",
        "        x = np.asarray(x)\n",
        "        if self.sigma is None or self.sigma.shape == ():\n",
        "            return x\n",
        "        elif self.sigma.ndim == 1:\n",
        "            if x.ndim == 1:\n",
        "                return x * self.cholsigmainv\n",
        "            else:\n",
        "                return x * self.cholsigmainv[:, None]\n",
        "        else:\n",
        "            return np.dot(self.cholsigmainv, x)\n",
        "\n",
        "\n",
        "    def loglike(self, params):\n",
        "        r\"\"\"\n",
        "        Compute the value of the Gaussian log-likelihood function at params.\n",
        "\n",
        "        Given the whitened design matrix, the log-likelihood is evaluated\n",
        "        at the parameter vector `params` for the dependent variable `endog`.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        params : array_like\n",
        "            The model parameters.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        float\n",
        "            The value of the log-likelihood function for a GLS Model.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        The log-likelihood function for the normal distribution is\n",
        "\n",
        "        .. math:: -\\frac{n}{2}\\log\\left(\\left(Y-\\hat{Y}\\right)^{\\prime}\n",
        "                   \\left(Y-\\hat{Y}\\right)\\right)\n",
        "                  -\\frac{n}{2}\\left(1+\\log\\left(\\frac{2\\pi}{n}\\right)\\right)\n",
        "                  -\\frac{1}{2}\\log\\left(\\left|\\Sigma\\right|\\right)\n",
        "\n",
        "        Y and Y-hat are whitened.\n",
        "        \"\"\"\n",
        "        # TODO: combine this with OLS/WLS loglike and add _det_sigma argument\n",
        "        nobs2 = self.nobs / 2.0\n",
        "        SSR = np.sum((self.wendog - np.dot(self.wexog, params))**2, axis=0)\n",
        "        llf = -np.log(SSR) * nobs2      # concentrated likelihood\n",
        "        llf -= (1+np.log(np.pi/nobs2))*nobs2  # with likelihood constant\n",
        "        if np.any(self.sigma):\n",
        "            # FIXME: robust-enough check? unneeded if _det_sigma gets defined\n",
        "            if self.sigma.ndim == 2:\n",
        "                det = np.linalg.slogdet(self.sigma)\n",
        "                llf -= .5*det[1]\n",
        "            else:\n",
        "                llf -= 0.5*np.sum(np.log(self.sigma))\n",
        "            # with error covariance matrix\n",
        "        return llf\n",
        "\n",
        "\n",
        "    def hessian_factor(self, params, scale=None, observed=True):\n",
        "        \"\"\"\n",
        "        Compute weights for calculating Hessian.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        params : ndarray\n",
        "            The parameter at which Hessian is evaluated.\n",
        "        scale : None or float\n",
        "            If scale is None, then the default scale will be calculated.\n",
        "            Default scale is defined by `self.scaletype` and set in fit.\n",
        "            If scale is not None, then it is used as a fixed scale.\n",
        "        observed : bool\n",
        "            If True, then the observed Hessian is returned. If false then the\n",
        "            expected information matrix is returned.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        ndarray\n",
        "            A 1d weight vector used in the calculation of the Hessian.\n",
        "            The hessian is obtained by `(exog.T * hessian_factor).dot(exog)`.\n",
        "        \"\"\"\n",
        "\n",
        "        if self.sigma is None or self.sigma.shape == ():\n",
        "            return np.ones(self.exog.shape[0])\n",
        "        elif self.sigma.ndim == 1:\n",
        "            return self.cholsigmainv\n",
        "        else:\n",
        "            return np.diag(self.cholsigmainv)\n",
        "\n",
        "    @Appender(_fit_regularized_doc)\n",
        "    def fit_regularized(self, method=\"elastic_net\", alpha=0.,\n",
        "                        L1_wt=1., start_params=None, profile_scale=False,\n",
        "                        refit=False, **kwargs):\n",
        "        # Need to adjust since RSS/n term in elastic net uses nominal\n",
        "        # n in denominator\n",
        "        if self.sigma is not None:\n",
        "            alpha = alpha * np.sum(1 / np.diag(self.sigma)) / len(self.endog)\n",
        "\n",
        "        rslt = OLS(self.wendog, self.wexog).fit_regularized(\n",
        "            method=method, alpha=alpha,\n",
        "            L1_wt=L1_wt,\n",
        "            start_params=start_params,\n",
        "            profile_scale=profile_scale,\n",
        "            refit=refit, **kwargs)\n",
        "\n",
        "        from statsmodels.base.elastic_net import (\n",
        "            RegularizedResults, RegularizedResultsWrapper)\n",
        "        rrslt = RegularizedResults(self, rslt.params)\n",
        "        return RegularizedResultsWrapper(rrslt)\n",
        "\n",
        "\n",
        "class WLS(RegressionModel):\n",
        "    __doc__ = \"\"\"\n",
        "    Weighted Least Squares\n",
        "\n",
        "    The weights are presumed to be (proportional to) the inverse of\n",
        "    the variance of the observations.  That is, if the variables are\n",
        "    to be transformed by 1/sqrt(W) you must supply weights = 1/W.\n",
        "\n",
        "    %(params)s\n",
        "    weights : array_like, optional\n",
        "        A 1d array of weights.  If you supply 1/W then the variables are\n",
        "        pre- multiplied by 1/sqrt(W).  If no weights are supplied the\n",
        "        default value is 1 and WLS results are the same as OLS.\n",
        "    %(extra_params)s\n",
        "\n",
        "    Attributes\n",
        "    ----------\n",
        "    weights : ndarray\n",
        "        The stored weights supplied as an argument.\n",
        "\n",
        "    See Also\n",
        "    --------\n",
        "    GLS : Fit a linear model using Generalized Least Squares.\n",
        "    OLS : Fit a linear model using Ordinary Least Squares.\n",
        "\n",
        "    Notes\n",
        "    -----\n",
        "    If the weights are a function of the data, then the post estimation\n",
        "    statistics such as fvalue and mse_model might not be correct, as the\n",
        "    package does not yet support no-constant regression.\n",
        "\n",
        "    Examples\n",
        "    --------\n",
        "    >>> import statsmodels.api as sm\n",
        "    >>> Y = [1,3,4,5,2,3,4]\n",
        "    >>> X = range(1,8)\n",
        "    >>> X = sm.add_constant(X)\n",
        "    >>> wls_model = sm.WLS(Y,X, weights=list(range(1,8)))\n",
        "    >>> results = wls_model.fit()\n",
        "    >>> results.params\n",
        "    array([ 2.91666667,  0.0952381 ])\n",
        "    >>> results.tvalues\n",
        "    array([ 2.0652652 ,  0.35684428])\n",
        "    >>> print(results.t_test([1, 0]))\n",
        "    <T test: effect=array([ 2.91666667]), sd=array([[ 1.41224801]]), t=array([[ 2.0652652]]), p=array([[ 0.04690139]]), df_denom=5>\n",
        "    >>> print(results.f_test([0, 1]))\n",
        "    <F test: F=array([[ 0.12733784]]), p=[[ 0.73577409]], df_denom=5, df_num=1>\n",
        "    \"\"\" % {'params': base._model_params_doc,\n",
        "           'extra_params': base._missing_param_doc + base._extra_param_doc}\n",
        "\n",
        "    def __init__(self, endog, exog, weights=1., missing='none', hasconst=None,\n",
        "                 **kwargs):\n",
        "        weights = np.array(weights)\n",
        "        if weights.shape == ():\n",
        "            if (missing == 'drop' and 'missing_idx' in kwargs and\n",
        "                    kwargs['missing_idx'] is not None):\n",
        "                # patsy may have truncated endog\n",
        "                weights = np.repeat(weights, len(kwargs['missing_idx']))\n",
        "            else:\n",
        "                weights = np.repeat(weights, len(endog))\n",
        "        # handle case that endog might be of len == 1\n",
        "        if len(weights) == 1:\n",
        "            weights = np.array([weights.squeeze()])\n",
        "        else:\n",
        "            weights = weights.squeeze()\n",
        "        super(WLS, self).__init__(endog, exog, missing=missing,\n",
        "                                  weights=weights, hasconst=hasconst, **kwargs)\n",
        "        nobs = self.exog.shape[0]\n",
        "        weights = self.weights\n",
        "        # Experimental normalization of weights\n",
        "        weights = weights / np.sum(weights) * nobs\n",
        "        if weights.size != nobs and weights.shape[0] != nobs:\n",
        "            raise ValueError('Weights must be scalar or same length as design')\n",
        "\n",
        "    def whiten(self, x):\n",
        "        \"\"\"\n",
        "        Whitener for WLS model, multiplies each column by sqrt(self.weights).\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        x : array_like\n",
        "            Data to be whitened.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        array_like\n",
        "            The whitened values sqrt(weights)*X.\n",
        "        \"\"\"\n",
        "\n",
        "        x = np.asarray(x)\n",
        "        if x.ndim == 1:\n",
        "            return x * np.sqrt(self.weights)\n",
        "        elif x.ndim == 2:\n",
        "            return np.sqrt(self.weights)[:, None] * x\n",
        "\n",
        "    def loglike(self, params):\n",
        "        r\"\"\"\n",
        "        Compute the value of the gaussian log-likelihood function at params.\n",
        "\n",
        "        Given the whitened design matrix, the log-likelihood is evaluated\n",
        "        at the parameter vector `params` for the dependent variable `Y`.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        params : array_like\n",
        "            The parameter estimates.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        float\n",
        "            The value of the log-likelihood function for a WLS Model.\n",
        "\n",
        "        Notes\n",
        "        --------\n",
        "        .. math:: -\\frac{n}{2}\\log SSR\n",
        "                  -\\frac{n}{2}\\left(1+\\log\\left(\\frac{2\\pi}{n}\\right)\\right)\n",
        "                  -\\frac{1}{2}\\log\\left(\\left|W\\right|\\right)\n",
        "\n",
        "        where :math:`W` is a diagonal weight matrix matrix and\n",
        "        :math:`SSR=\\left(Y-\\hat{Y}\\right)^\\prime W \\left(Y-\\hat{Y}\\right)` is\n",
        "        the sum of the squared weighted residuals.\n",
        "        \"\"\"\n",
        "        nobs2 = self.nobs / 2.0\n",
        "        SSR = np.sum((self.wendog - np.dot(self.wexog, params))**2, axis=0)\n",
        "        llf = -np.log(SSR) * nobs2      # concentrated likelihood\n",
        "        llf -= (1+np.log(np.pi/nobs2))*nobs2  # with constant\n",
        "        llf += 0.5 * np.sum(np.log(self.weights))\n",
        "        return llf\n",
        "\n",
        "    def hessian_factor(self, params, scale=None, observed=True):\n",
        "        \"\"\"\n",
        "        Compute the weights for calculating the Hessian.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        params : ndarray\n",
        "            The parameter at which Hessian is evaluated.\n",
        "        scale : None or float\n",
        "            If scale is None, then the default scale will be calculated.\n",
        "            Default scale is defined by `self.scaletype` and set in fit.\n",
        "            If scale is not None, then it is used as a fixed scale.\n",
        "        observed : bool\n",
        "            If True, then the observed Hessian is returned. If false then the\n",
        "            expected information matrix is returned.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        ndarray\n",
        "            A 1d weight vector used in the calculation of the Hessian.\n",
        "            The hessian is obtained by `(exog.T * hessian_factor).dot(exog)`.\n",
        "        \"\"\"\n",
        "\n",
        "        return self.weights\n",
        "\n",
        "    @Appender(_fit_regularized_doc)\n",
        "    def fit_regularized(self, method=\"elastic_net\", alpha=0.,\n",
        "                        L1_wt=1., start_params=None, profile_scale=False,\n",
        "                        refit=False, **kwargs):\n",
        "        # Docstring attached below\n",
        "\n",
        "        # Need to adjust since RSS/n in elastic net uses nominal n in\n",
        "        # denominator\n",
        "        alpha = alpha * np.sum(self.weights) / len(self.weights)\n",
        "\n",
        "        rslt = OLS(self.wendog, self.wexog).fit_regularized(\n",
        "            method=method, alpha=alpha,\n",
        "            L1_wt=L1_wt,\n",
        "            start_params=start_params,\n",
        "            profile_scale=profile_scale,\n",
        "            refit=refit, **kwargs)\n",
        "\n",
        "        from statsmodels.base.elastic_net import (\n",
        "            RegularizedResults, RegularizedResultsWrapper)\n",
        "        rrslt = RegularizedResults(self, rslt.params)\n",
        "        return RegularizedResultsWrapper(rrslt)\n",
        "\n",
        "\n",
        "\n",
        "class OLS(WLS):\n",
        "    __doc__ = \"\"\"\n",
        "    Ordinary Least Squares\n",
        "\n",
        "    %(params)s\n",
        "    %(extra_params)s\n",
        "\n",
        "    Attributes\n",
        "    ----------\n",
        "    weights : scalar\n",
        "        Has an attribute weights = array(1.0) due to inheritance from WLS.\n",
        "\n",
        "    See Also\n",
        "    --------\n",
        "    WLS : Fit a linear model using Weighted Least Squares.\n",
        "    GLS : Fit a linear model using Generalized Least Squares.\n",
        "\n",
        "    Notes\n",
        "    -----\n",
        "    No constant is added by the model unless you are using formulas.\n",
        "\n",
        "    Examples\n",
        "    --------\n",
        "    >>> import statsmodels.api as sm\n",
        "    >>> Y = [1,3,4,5,2,3,4]\n",
        "    >>> X = range(1,8)\n",
        "    >>> X = sm.add_constant(X)\n",
        "    >>> model = sm.OLS(Y,X)\n",
        "    >>> results = model.fit()\n",
        "    >>> results.params\n",
        "    array([ 2.14285714,  0.25      ])\n",
        "\n",
        "    >>> results.tvalues\n",
        "    array([ 1.87867287,  0.98019606])\n",
        "\n",
        "    >>> print(results.t_test([1, 0]))\n",
        "    <T test: effect=array([ 2.14285714]), sd=array([[ 1.14062282]]), t=array([[ 1.87867287]]), p=array([[ 0.05953974]]), df_denom=5>\n",
        "    >>> print(results.f_test(np.identity(2)))\n",
        "    <F test: F=array([[ 19.46078431]]), p=[[ 0.00437251]], df_denom=5, df_num=2>\n",
        "    \"\"\" % {'params': base._model_params_doc,\n",
        "           'extra_params': base._missing_param_doc + base._extra_param_doc}\n",
        "\n",
        "    # TODO: change example to use datasets.  This was the point of datasets!\n",
        "    def __init__(self, endog, exog=None, missing='none', hasconst=None,\n",
        "                 **kwargs):\n",
        "        super(OLS, self).__init__(endog, exog, missing=missing,\n",
        "                                  hasconst=hasconst, **kwargs)\n",
        "        if \"weights\" in self._init_keys:\n",
        "            self._init_keys.remove(\"weights\")\n",
        "\n",
        "    def loglike(self, params, scale=None):\n",
        "        \"\"\"\n",
        "        The likelihood function for the OLS model.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        params : array_like\n",
        "            The coefficients with which to estimate the log-likelihood.\n",
        "        scale : float or None\n",
        "            If None, return the profile (concentrated) log likelihood\n",
        "            (profiled over the scale parameter), else return the\n",
        "            log-likelihood using the given scale value.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        float\n",
        "            The likelihood function evaluated at params.\n",
        "        \"\"\"\n",
        "        nobs2 = self.nobs / 2.0\n",
        "        nobs = float(self.nobs)\n",
        "        resid = self.endog - np.dot(self.exog, params)\n",
        "        if hasattr(self, 'offset'):\n",
        "            resid -= self.offset\n",
        "        ssr = np.sum(resid**2)\n",
        "        if scale is None:\n",
        "            # profile log likelihood\n",
        "            llf = -nobs2*np.log(2*np.pi) - nobs2*np.log(ssr / nobs) - nobs2\n",
        "        else:\n",
        "            # log-likelihood\n",
        "            llf = -nobs2 * np.log(2 * np.pi * scale) - ssr / (2*scale)\n",
        "        return llf\n",
        "\n",
        "\n",
        "    def whiten(self, x):\n",
        "        \"\"\"\n",
        "        OLS model whitener does nothing.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        x : array_like\n",
        "            Data to be whitened.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        array_like\n",
        "            The input array unmodified.\n",
        "\n",
        "        See Also\n",
        "        --------\n",
        "        OLS : Fit a linear model using Ordinary Least Squares.\n",
        "        \"\"\"\n",
        "        return x\n",
        "\n",
        "\n",
        "    def score(self, params, scale=None):\n",
        "        \"\"\"\n",
        "        Evaluate the score function at a given point.\n",
        "\n",
        "        The score corresponds to the profile (concentrated)\n",
        "        log-likelihood in which the scale parameter has been profiled\n",
        "        out.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        params : array_like\n",
        "            The parameter vector at which the score function is\n",
        "            computed.\n",
        "        scale : float or None\n",
        "            If None, return the profile (concentrated) log likelihood\n",
        "            (profiled over the scale parameter), else return the\n",
        "            log-likelihood using the given scale value.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        ndarray\n",
        "            The score vector.\n",
        "        \"\"\"\n",
        "\n",
        "        if not hasattr(self, \"_wexog_xprod\"):\n",
        "            self._setup_score_hess()\n",
        "\n",
        "        xtxb = np.dot(self._wexog_xprod, params)\n",
        "        sdr = -self._wexog_x_wendog + xtxb\n",
        "\n",
        "        if scale is None:\n",
        "            ssr = self._wendog_xprod - 2 * np.dot(self._wexog_x_wendog.T,\n",
        "                                                  params)\n",
        "            ssr += np.dot(params, xtxb)\n",
        "            return -self.nobs * sdr / ssr\n",
        "        else:\n",
        "            return -sdr / scale\n",
        "\n",
        "\n",
        "    def _setup_score_hess(self):\n",
        "        y = self.wendog\n",
        "        if hasattr(self, 'offset'):\n",
        "            y = y - self.offset\n",
        "        self._wendog_xprod = np.sum(y * y)\n",
        "        self._wexog_xprod = np.dot(self.wexog.T, self.wexog)\n",
        "        self._wexog_x_wendog = np.dot(self.wexog.T, y)\n",
        "\n",
        "    def hessian(self, params, scale=None):\n",
        "        \"\"\"\n",
        "        Evaluate the Hessian function at a given point.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        params : array_like\n",
        "            The parameter vector at which the Hessian is computed.\n",
        "        scale : float or None\n",
        "            If None, return the profile (concentrated) log likelihood\n",
        "            (profiled over the scale parameter), else return the\n",
        "            log-likelihood using the given scale value.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        ndarray\n",
        "            The Hessian matrix.\n",
        "        \"\"\"\n",
        "\n",
        "        if not hasattr(self, \"_wexog_xprod\"):\n",
        "            self._setup_score_hess()\n",
        "\n",
        "        xtxb = np.dot(self._wexog_xprod, params)\n",
        "\n",
        "        if scale is None:\n",
        "            ssr = self._wendog_xprod - 2 * np.dot(self._wexog_x_wendog.T,\n",
        "                                                  params)\n",
        "            ssr += np.dot(params, xtxb)\n",
        "            ssrp = -2*self._wexog_x_wendog + 2*xtxb\n",
        "            hm = self._wexog_xprod / ssr - np.outer(ssrp, ssrp) / ssr**2\n",
        "            return -self.nobs * hm / 2\n",
        "        else:\n",
        "            return -self._wexog_xprod / scale\n",
        "\n",
        "\n",
        "    def hessian_factor(self, params, scale=None, observed=True):\n",
        "        \"\"\"\n",
        "        Calculate the weights for the Hessian.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        params : ndarray\n",
        "            The parameter at which Hessian is evaluated.\n",
        "        scale : None or float\n",
        "            If scale is None, then the default scale will be calculated.\n",
        "            Default scale is defined by `self.scaletype` and set in fit.\n",
        "            If scale is not None, then it is used as a fixed scale.\n",
        "        observed : bool\n",
        "            If True, then the observed Hessian is returned. If false then the\n",
        "            expected information matrix is returned.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        ndarray\n",
        "            A 1d weight vector used in the calculation of the Hessian.\n",
        "            The hessian is obtained by `(exog.T * hessian_factor).dot(exog)`.\n",
        "        \"\"\"\n",
        "\n",
        "        return np.ones(self.exog.shape[0])\n",
        "\n",
        "\n",
        "    @Appender(_fit_regularized_doc)\n",
        "    def fit_regularized(self, method=\"elastic_net\", alpha=0.,\n",
        "                        L1_wt=1., start_params=None, profile_scale=False,\n",
        "                        refit=False, **kwargs):\n",
        "\n",
        "        # In the future we could add support for other penalties, e.g. SCAD.\n",
        "        if method not in (\"elastic_net\", \"sqrt_lasso\"):\n",
        "            msg = \"Unknown method '%s' for fit_regularized\" % method\n",
        "            raise ValueError(msg)\n",
        "\n",
        "        # Set default parameters.\n",
        "        defaults = {\"maxiter\":  50, \"cnvrg_tol\": 1e-10,\n",
        "                    \"zero_tol\": 1e-8}\n",
        "        defaults.update(kwargs)\n",
        "\n",
        "        if method == \"sqrt_lasso\":\n",
        "            from statsmodels.base.elastic_net import (\n",
        "                RegularizedResults, RegularizedResultsWrapper\n",
        "            )\n",
        "            params = self._sqrt_lasso(alpha, refit, defaults[\"zero_tol\"])\n",
        "            results = RegularizedResults(self, params)\n",
        "            return RegularizedResultsWrapper(results)\n",
        "\n",
        "        from statsmodels.base.elastic_net import fit_elasticnet\n",
        "\n",
        "        if L1_wt == 0:\n",
        "            return self._fit_ridge(alpha)\n",
        "\n",
        "        # If a scale parameter is passed in, the non-profile\n",
        "        # likelihood (residual sum of squares divided by -2) is used,\n",
        "        # otherwise the profile likelihood is used.\n",
        "        if profile_scale:\n",
        "            loglike_kwds = {}\n",
        "            score_kwds = {}\n",
        "            hess_kwds = {}\n",
        "        else:\n",
        "            loglike_kwds = {\"scale\": 1}\n",
        "            score_kwds = {\"scale\": 1}\n",
        "            hess_kwds = {\"scale\": 1}\n",
        "\n",
        "        return fit_elasticnet(self, method=method,\n",
        "                              alpha=alpha,\n",
        "                              L1_wt=L1_wt,\n",
        "                              start_params=start_params,\n",
        "                              loglike_kwds=loglike_kwds,\n",
        "                              score_kwds=score_kwds,\n",
        "                              hess_kwds=hess_kwds,\n",
        "                              refit=refit,\n",
        "                              check_step=False,\n",
        "                              **defaults)\n",
        "\n",
        "\n",
        "    def _sqrt_lasso(self, alpha, refit, zero_tol):\n",
        "\n",
        "        try:\n",
        "            import cvxopt\n",
        "        except ImportError:\n",
        "            msg = 'sqrt_lasso fitting requires the cvxopt module'\n",
        "            raise ValueError(msg)\n",
        "\n",
        "        n = len(self.endog)\n",
        "        p = self.exog.shape[1]\n",
        "\n",
        "        h0 = cvxopt.matrix(0., (2*p+1, 1))\n",
        "        h1 = cvxopt.matrix(0., (n+1, 1))\n",
        "        h1[1:, 0] = cvxopt.matrix(self.endog, (n, 1))\n",
        "\n",
        "        G0 = cvxopt.spmatrix([], [], [], (2*p+1, 2*p+1))\n",
        "        for i in range(1, 2*p+1):\n",
        "            G0[i, i] = -1\n",
        "        G1 = cvxopt.matrix(0., (n+1, 2*p+1))\n",
        "        G1[0, 0] = -1\n",
        "        G1[1:, 1:p+1] = self.exog\n",
        "        G1[1:, p+1:] = -self.exog\n",
        "\n",
        "        c = cvxopt.matrix(alpha / n, (2*p + 1, 1))\n",
        "        c[0] = 1 / np.sqrt(n)\n",
        "\n",
        "        from cvxopt import solvers\n",
        "        solvers.options[\"show_progress\"] = False\n",
        "\n",
        "        rslt = solvers.socp(c, Gl=G0, hl=h0, Gq=[G1], hq=[h1])\n",
        "        x = np.asarray(rslt['x']).flat\n",
        "        bp = x[1:p+1]\n",
        "        bn = x[p+1:]\n",
        "        params = bp - bn\n",
        "\n",
        "        if not refit:\n",
        "            return params\n",
        "\n",
        "        ii = np.flatnonzero(np.abs(params) > zero_tol)\n",
        "        rfr = OLS(self.endog, self.exog[:, ii]).fit()\n",
        "        params *= 0\n",
        "        params[ii] = rfr.params\n",
        "\n",
        "        return params\n",
        "\n",
        "    def _fit_ridge(self, alpha):\n",
        "        \"\"\"\n",
        "        Fit a linear model using ridge regression.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        alpha : scalar or array_like\n",
        "            The penalty weight.  If a scalar, the same penalty weight\n",
        "            applies to all variables in the model.  If a vector, it\n",
        "            must have the same length as `params`, and contains a\n",
        "            penalty weight for each coefficient.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        Equivalent to fit_regularized with L1_wt = 0 (but implemented\n",
        "        more efficiently).\n",
        "        \"\"\"\n",
        "\n",
        "        u, s, vt = np.linalg.svd(self.exog, 0)\n",
        "        v = vt.T\n",
        "        q = np.dot(u.T, self.endog) * s\n",
        "        s2 = s * s\n",
        "        if np.isscalar(alpha):\n",
        "            sd = s2 + alpha * self.nobs\n",
        "            params = q / sd\n",
        "            params = np.dot(v, params)\n",
        "        else:\n",
        "            vtav = self.nobs * np.dot(vt, alpha[:, None] * v)\n",
        "            d = np.diag(vtav) + s2\n",
        "            np.fill_diagonal(vtav, d)\n",
        "            r = np.linalg.solve(vtav, q)\n",
        "            params = np.dot(v, r)\n",
        "\n",
        "        from statsmodels.base.elastic_net import RegularizedResults\n",
        "        return RegularizedResults(self, params)\n",
        "\n",
        "\n",
        "\n",
        "class GLSAR(GLS):\n",
        "    __doc__ = \"\"\"\n",
        "    Generalized Least Squares with AR covariance structure\n",
        "\n",
        "    %(params)s\n",
        "    rho : int\n",
        "        The order of the autoregressive covariance.\n",
        "    %(extra_params)s\n",
        "\n",
        "    Notes\n",
        "    -----\n",
        "    GLSAR is considered to be experimental.\n",
        "    The linear autoregressive process of order p--AR(p)--is defined as:\n",
        "    TODO\n",
        "\n",
        "    Examples\n",
        "    --------\n",
        "    >>> import statsmodels.api as sm\n",
        "    >>> X = range(1,8)\n",
        "    >>> X = sm.add_constant(X)\n",
        "    >>> Y = [1,3,4,5,8,10,9]\n",
        "    >>> model = sm.GLSAR(Y, X, rho=2)\n",
        "    >>> for i in range(6):\n",
        "    ...     results = model.fit()\n",
        "    ...     print(\"AR coefficients: {0}\".format(model.rho))\n",
        "    ...     rho, sigma = sm.regression.yule_walker(results.resid,\n",
        "    ...                                            order=model.order)\n",
        "    ...     model = sm.GLSAR(Y, X, rho)\n",
        "    ...\n",
        "    AR coefficients: [ 0.  0.]\n",
        "    AR coefficients: [-0.52571491 -0.84496178]\n",
        "    AR coefficients: [-0.6104153  -0.86656458]\n",
        "    AR coefficients: [-0.60439494 -0.857867  ]\n",
        "    AR coefficients: [-0.6048218  -0.85846157]\n",
        "    AR coefficients: [-0.60479146 -0.85841922]\n",
        "    >>> results.params\n",
        "    array([-0.66661205,  1.60850853])\n",
        "    >>> results.tvalues\n",
        "    array([ -2.10304127,  21.8047269 ])\n",
        "    >>> print(results.t_test([1, 0]))\n",
        "    <T test: effect=array([-0.66661205]), sd=array([[ 0.31697526]]), t=array([[-2.10304127]]), p=array([[ 0.06309969]]), df_denom=3>\n",
        "    >>> print(results.f_test(np.identity(2)))\n",
        "    <F test: F=array([[ 1815.23061844]]), p=[[ 0.00002372]], df_denom=3, df_num=2>\n",
        "\n",
        "    Or, equivalently\n",
        "\n",
        "    >>> model2 = sm.GLSAR(Y, X, rho=2)\n",
        "    >>> res = model2.iterative_fit(maxiter=6)\n",
        "    >>> model2.rho\n",
        "    array([-0.60479146, -0.85841922])\n",
        "    \"\"\" % {'params': base._model_params_doc,\n",
        "           'extra_params': base._missing_param_doc + base._extra_param_doc}\n",
        "    # TODO: Complete docstring\n",
        "\n",
        "    def __init__(self, endog, exog=None, rho=1, missing='none', hasconst=None,\n",
        "                 **kwargs):\n",
        "        # this looks strange, interpreting rho as order if it is int\n",
        "        if isinstance(rho, np.int):\n",
        "            self.order = rho\n",
        "            self.rho = np.zeros(self.order, np.float64)\n",
        "        else:\n",
        "            self.rho = np.squeeze(np.asarray(rho))\n",
        "            if len(self.rho.shape) not in [0, 1]:\n",
        "                raise ValueError(\"AR parameters must be a scalar or a vector\")\n",
        "            if self.rho.shape == ():\n",
        "                self.rho.shape = (1,)\n",
        "            self.order = self.rho.shape[0]\n",
        "        if exog is None:\n",
        "            # JP this looks wrong, should be a regression on constant\n",
        "            # results for rho estimate now identical to yule-walker on y\n",
        "            # super(AR, self).__init__(endog, add_constant(endog))\n",
        "            super(GLSAR, self).__init__(endog, np.ones((endog.shape[0], 1)),\n",
        "                                        missing=missing, hasconst=None,\n",
        "                                        **kwargs)\n",
        "        else:\n",
        "            super(GLSAR, self).__init__(endog, exog, missing=missing,\n",
        "                                        **kwargs)\n",
        "\n",
        "    def iterative_fit(self, maxiter=3, rtol=1e-4, **kwargs):\n",
        "        \"\"\"\n",
        "        Perform an iterative two-stage procedure to estimate a GLS model.\n",
        "\n",
        "        The model is assumed to have AR(p) errors, AR(p) parameters and\n",
        "        regression coefficients are estimated iteratively.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        maxiter : int, optional\n",
        "            The number of iterations.\n",
        "        rtol : float, optional\n",
        "            Relative tolerance between estimated coefficients to stop the\n",
        "            estimation.  Stops if max(abs(last - current) / abs(last)) < rtol.\n",
        "        **kwargs\n",
        "            Additional keyword arguments passed to `fit`.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        RegressionResults\n",
        "            The results computed using an iterative fit.\n",
        "        \"\"\"\n",
        "        # TODO: update this after going through example.\n",
        "        converged = False\n",
        "        i = -1  # need to initialize for maxiter < 1 (skip loop)\n",
        "        history = {'params': [], 'rho': [self.rho]}\n",
        "        for i in range(maxiter - 1):\n",
        "            if hasattr(self, 'pinv_wexog'):\n",
        "                del self.pinv_wexog\n",
        "            self.initialize()\n",
        "            results = self.fit()\n",
        "            history['params'].append(results.params)\n",
        "            if i == 0:\n",
        "                last = results.params\n",
        "            else:\n",
        "                diff = np.max(np.abs(last - results.params) / np.abs(last))\n",
        "                if diff < rtol:\n",
        "                    converged = True\n",
        "                    break\n",
        "                last = results.params\n",
        "            self.rho, _ = yule_walker(results.resid,\n",
        "                                      order=self.order, df=None)\n",
        "            history['rho'].append(self.rho)\n",
        "\n",
        "        # why not another call to self.initialize\n",
        "        # Use kwarg to insert history\n",
        "        if not converged and maxiter > 0:\n",
        "            # maxiter <= 0 just does OLS\n",
        "            if hasattr(self, 'pinv_wexog'):\n",
        "                del self.pinv_wexog\n",
        "            self.initialize()\n",
        "\n",
        "        # if converged then this is a duplicate fit, because we did not\n",
        "        # update rho\n",
        "        results = self.fit(history=history, **kwargs)\n",
        "        results.iter = i + 1\n",
        "        # add last fit to history, not if duplicate fit\n",
        "        if not converged:\n",
        "            results.history['params'].append(results.params)\n",
        "            results.iter += 1\n",
        "\n",
        "        results.converged = converged\n",
        "\n",
        "        return results\n",
        "\n",
        "\n",
        "    def whiten(self, x):\n",
        "        \"\"\"\n",
        "        Whiten a series of columns according to an AR(p) covariance structure.\n",
        "\n",
        "        Whitening using this method drops the initial p observations.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        x : array_like\n",
        "            The data to be whitened.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        ndarray\n",
        "            The whitened data.\n",
        "        \"\"\"\n",
        "        # TODO: notation for AR process\n",
        "        x = np.asarray(x, np.float64)\n",
        "        _x = x.copy()\n",
        "\n",
        "        # the following loops over the first axis,  works for 1d and nd\n",
        "        for i in range(self.order):\n",
        "            _x[(i + 1):] = _x[(i + 1):] - self.rho[i] * x[0:-(i + 1)]\n",
        "        return _x[self.order:]\n",
        "\n",
        "\n",
        "\n",
        "def yule_walker(x, order=1, method=\"unbiased\", df=None, inv=False,\n",
        "                demean=True):\n",
        "    \"\"\"\n",
        "    Estimate AR(p) parameters from a sequence using the Yule-Walker equations.\n",
        "\n",
        "    Unbiased or maximum-likelihood estimator (mle)\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    x : array_like\n",
        "        A 1d array.\n",
        "    order : int, optional\n",
        "        The order of the autoregressive process.  Default is 1.\n",
        "    method : str, optional\n",
        "       Method can be 'unbiased' or 'mle' and this determines\n",
        "       denominator in estimate of autocorrelation function (ACF) at\n",
        "       lag k. If 'mle', the denominator is n=X.shape[0], if 'unbiased'\n",
        "       the denominator is n-k.  The default is unbiased.\n",
        "    df : int, optional\n",
        "       Specifies the degrees of freedom. If `df` is supplied, then it\n",
        "       is assumed the X has `df` degrees of freedom rather than `n`.\n",
        "       Default is None.\n",
        "    inv : bool\n",
        "        If inv is True the inverse of R is also returned.  Default is\n",
        "        False.\n",
        "    demean : bool\n",
        "        True, the mean is subtracted from `X` before estimation.\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    rho : ndarray\n",
        "        AR(p) coefficients computed using the Yule-Walker method.\n",
        "    sigma : float\n",
        "        The estimate of the residual standard deviation.\n",
        "\n",
        "    See Also\n",
        "    --------\n",
        "    burg : Burg's AR estimator.\n",
        "\n",
        "    Notes\n",
        "    -----\n",
        "    See https://en.wikipedia.org/wiki/Autoregressive_moving_average_model for\n",
        "    further details.\n",
        "\n",
        "    Examples\n",
        "    --------\n",
        "    >>> import statsmodels.api as sm\n",
        "    >>> from statsmodels.datasets.sunspots import load\n",
        "    >>> data = load(as_pandas=False)\n",
        "    >>> rho, sigma = sm.regression.yule_walker(data.endog, order=4,\n",
        "    ...                                        method=\"mle\")\n",
        "\n",
        "    >>> rho\n",
        "    array([ 1.28310031, -0.45240924, -0.20770299,  0.04794365])\n",
        "    >>> sigma\n",
        "    16.808022730464351\n",
        "    \"\"\"\n",
        "    # TODO: define R better, look back at notes and technical notes on YW.\n",
        "    # First link here is useful\n",
        "    # http://www-stat.wharton.upenn.edu/~steele/Courses/956/ResourceDetails/YuleWalkerAndMore.htm\n",
        "    method = str(method).lower()\n",
        "    if method not in [\"unbiased\", \"mle\"]:\n",
        "        raise ValueError(\"ACF estimation method must be 'unbiased' or 'MLE'\")\n",
        "    x = np.array(x, dtype=np.float64)\n",
        "    if demean:\n",
        "        x -= x.mean()\n",
        "    n = df or x.shape[0]\n",
        "\n",
        "    # this handles df_resid ie., n - p\n",
        "    adj_needed = method == \"unbiased\"\n",
        "\n",
        "    if x.ndim > 1 and x.shape[1] != 1:\n",
        "        raise ValueError(\"expecting a vector to estimate AR parameters\")\n",
        "    r = np.zeros(order+1, np.float64)\n",
        "    r[0] = (x ** 2).sum() / n\n",
        "    for k in range(1, order+1):\n",
        "        r[k] = (x[0:-k] * x[k:]).sum() / (n - k * adj_needed)\n",
        "    R = toeplitz(r[:-1])\n",
        "\n",
        "    rho = np.linalg.solve(R, r[1:])\n",
        "    sigmasq = r[0] - (r[1:]*rho).sum()\n",
        "    if inv:\n",
        "        return rho, np.sqrt(sigmasq), np.linalg.inv(R)\n",
        "    else:\n",
        "        return rho, np.sqrt(sigmasq)\n",
        "\n",
        "\n",
        "def burg(endog, order=1, demean=True):\n",
        "    \"\"\"\n",
        "    Compute Burg's AP(p) parameter estimator.\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    endog : array_like\n",
        "        The endogenous variable.\n",
        "    order : int, optional\n",
        "        Order of the AR.  Default is 1.\n",
        "    demean : bool, optional\n",
        "        Flag indicating to subtract the mean from endog before estimation.\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    rho : ndarray\n",
        "        The AR(p) coefficients computed using Burg's algorithm.\n",
        "    sigma2 : float\n",
        "        The estimate of the residual variance.\n",
        "\n",
        "    See Also\n",
        "    --------\n",
        "    yule_walker : Estimate AR parameters using the Yule-Walker method.\n",
        "\n",
        "    Notes\n",
        "    -----\n",
        "    AR model estimated includes a constant that is estimated using the sample\n",
        "    mean (see [1]_). This value is not reported.\n",
        "\n",
        "    References\n",
        "    ----------\n",
        "    .. [1] Brockwell, P.J. and Davis, R.A., 2016. Introduction to time series\n",
        "        and forecasting. Springer.\n",
        "\n",
        "    Examples\n",
        "    --------\n",
        "    >>> import statsmodels.api as sm\n",
        "    >>> from statsmodels.datasets.sunspots import load\n",
        "    >>> data = load(as_pandas=True)\n",
        "    >>> rho, sigma2 = sm.regression.linear_model.burg(data.endog, order=4)\n",
        "\n",
        "    >>> rho\n",
        "    array([ 1.30934186, -0.48086633, -0.20185982,  0.05501941])\n",
        "    >>> sigma2\n",
        "    271.2467306963966\n",
        "    \"\"\"\n",
        "    # Avoid circular imports\n",
        "    from statsmodels.tsa.stattools import levinson_durbin_pacf, pacf_burg\n",
        "\n",
        "    endog = np.squeeze(np.asarray(endog))\n",
        "    if endog.ndim != 1:\n",
        "        raise ValueError('endog must be 1-d or squeezable to 1-d.')\n",
        "    order = int(order)\n",
        "    if order < 1:\n",
        "        raise ValueError('order must be an integer larger than 1')\n",
        "    if demean:\n",
        "        endog = endog - endog.mean()\n",
        "    pacf, sigma = pacf_burg(endog, order, demean=demean)\n",
        "    ar, _ = levinson_durbin_pacf(pacf)\n",
        "    return ar, sigma[-1]\n",
        "\n",
        "\n",
        "class RegressionResults(base.LikelihoodModelResults):\n",
        "    r\"\"\"\n",
        "    This class summarizes the fit of a linear regression model.\n",
        "\n",
        "    It handles the output of contrasts, estimates of covariance, etc.\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    model : RegressionModel\n",
        "        The regression model instance.\n",
        "    params : ndarray\n",
        "        The estimated parameters.\n",
        "    normalized_cov_params : ndarray\n",
        "        The normalized covariance parameters.\n",
        "    scale : float\n",
        "        The estimated scale of the residuals.\n",
        "    cov_type : str\n",
        "        The covariance estimator used in the results.\n",
        "    cov_kwds : dict\n",
        "        Additional keywords used in the covariance specification.\n",
        "    use_t : bool\n",
        "        Flag indicating to use the Student's t in inference.\n",
        "    **kwargs\n",
        "        Additional keyword arguments used to initialize the results.\n",
        "\n",
        "    Attributes\n",
        "    ----------\n",
        "    pinv_wexog\n",
        "        See model class docstring for implementation details.\n",
        "    cov_type\n",
        "        Parameter covariance estimator used for standard errors and t-stats.\n",
        "    df_model\n",
        "        Model degrees of freedom. The number of regressors `p`. Does not\n",
        "        include the constant if one is present.\n",
        "    df_resid\n",
        "        Residual degrees of freedom. `n - p - 1`, if a constant is present.\n",
        "        `n - p` if a constant is not included.\n",
        "    het_scale\n",
        "        adjusted squared residuals for heteroscedasticity robust standard\n",
        "        errors. Is only available after `HC#_se` or `cov_HC#` is called.\n",
        "        See HC#_se for more information.\n",
        "    history\n",
        "        Estimation history for iterative estimators.\n",
        "    model\n",
        "        A pointer to the model instance that called fit() or results.\n",
        "    params\n",
        "        The linear coefficients that minimize the least squares\n",
        "        criterion.  This is usually called Beta for the classical\n",
        "        linear model.\n",
        "    \"\"\"\n",
        "\n",
        "    _cache = {}  # needs to be a class attribute for scale setter?\n",
        "\n",
        "    def __init__(self, model, params, normalized_cov_params=None, scale=1.,\n",
        "                 cov_type='nonrobust', cov_kwds=None, use_t=None, **kwargs):\n",
        "        super(RegressionResults, self).__init__(\n",
        "            model, params, normalized_cov_params, scale)\n",
        "\n",
        "        self._cache = {}\n",
        "        if hasattr(model, 'wexog_singular_values'):\n",
        "            self._wexog_singular_values = model.wexog_singular_values\n",
        "        else:\n",
        "            self._wexog_singular_values = None\n",
        "\n",
        "        self.df_model = model.df_model\n",
        "        self.df_resid = model.df_resid\n",
        "\n",
        "        if cov_type == 'nonrobust':\n",
        "            self.cov_type = 'nonrobust'\n",
        "            self.cov_kwds = {\n",
        "                'description': 'Standard Errors assume that the ' +\n",
        "                'covariance matrix of the errors is correctly ' +\n",
        "                'specified.'}\n",
        "            if use_t is None:\n",
        "                use_t = True  # TODO: class default\n",
        "            self.use_t = use_t\n",
        "        else:\n",
        "            if cov_kwds is None:\n",
        "                cov_kwds = {}\n",
        "            if 'use_t' in cov_kwds:\n",
        "                # TODO: we want to get rid of 'use_t' in cov_kwds\n",
        "                use_t_2 = cov_kwds.pop('use_t')\n",
        "                if use_t is None:\n",
        "                    use_t = use_t_2\n",
        "                # TODO: warn or not?\n",
        "            self.get_robustcov_results(cov_type=cov_type, use_self=True,\n",
        "                                       use_t=use_t, **cov_kwds)\n",
        "        for key in kwargs:\n",
        "            setattr(self, key, kwargs[key])\n",
        "\n",
        "    def conf_int(self, alpha=.05, cols=None):\n",
        "        \"\"\"\n",
        "        Compute the confidence interval of the fitted parameters.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        alpha : float, optional\n",
        "            The `alpha` level for the confidence interval. The default\n",
        "            `alpha` = .05 returns a 95% confidence interval.\n",
        "        cols : array_like, optional\n",
        "            Columns to included in returned confidence intervals.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        array_like\n",
        "            The confidence intervals.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        The confidence interval is based on Student's t-distribution.\n",
        "        \"\"\"\n",
        "        # keep method for docstring for now\n",
        "        ci = super(RegressionResults, self).conf_int(alpha=alpha, cols=cols)\n",
        "        return ci\n",
        "\n",
        "\n",
        "    @cache_readonly\n",
        "    def nobs(self):\n",
        "        \"\"\"Number of observations n.\"\"\"\n",
        "        return float(self.model.wexog.shape[0])\n",
        "\n",
        "    @cache_readonly\n",
        "    def fittedvalues(self):\n",
        "        \"\"\"The predicted values for the original (unwhitened) design.\"\"\"\n",
        "        return self.model.predict(self.params, self.model.exog)\n",
        "\n",
        "    @cache_readonly\n",
        "    def wresid(self):\n",
        "        \"\"\"\n",
        "        The residuals of the transformed/whitened regressand and regressor(s).\n",
        "        \"\"\"\n",
        "        return self.model.wendog - self.model.predict(\n",
        "            self.params, self.model.wexog)\n",
        "\n",
        "    @cache_readonly\n",
        "    def resid(self):\n",
        "        \"\"\"The residuals of the model.\"\"\"\n",
        "        return self.model.endog - self.model.predict(\n",
        "            self.params, self.model.exog)\n",
        "\n",
        "    # TODO: fix writable example\n",
        "    @cache_writable()\n",
        "    def scale(self):\n",
        "        \"\"\"\n",
        "        A scale factor for the covariance matrix.\n",
        "\n",
        "        The Default value is ssr/(n-p).  Note that the square root of `scale`\n",
        "        is often called the standard error of the regression.\n",
        "        \"\"\"\n",
        "        wresid = self.wresid\n",
        "        return np.dot(wresid, wresid) / self.df_resid\n",
        "\n",
        "\n",
        "    @cache_readonly\n",
        "    def ssr(self):\n",
        "        \"\"\"Sum of squared (whitened) residuals.\"\"\"\n",
        "        wresid = self.wresid\n",
        "        return np.dot(wresid, wresid)\n",
        "\n",
        "    @cache_readonly\n",
        "    def centered_tss(self):\n",
        "        \"\"\"The total (weighted) sum of squares centered about the mean.\"\"\"\n",
        "        model = self.model\n",
        "        weights = getattr(model, 'weights', None)\n",
        "        sigma = getattr(model, 'sigma', None)\n",
        "        if weights is not None:\n",
        "            mean = np.average(model.endog, weights=weights)\n",
        "            return np.sum(weights * (model.endog - mean)**2)\n",
        "        elif sigma is not None:\n",
        "            # Exactly matches WLS when sigma is diagonal\n",
        "            iota = np.ones_like(model.endog)\n",
        "            iota = model.whiten(iota)\n",
        "            mean = model.wendog.dot(iota) / iota.dot(iota)\n",
        "            err = model.endog - mean\n",
        "            err = model.whiten(err)\n",
        "            return np.sum(err**2)\n",
        "        else:\n",
        "            centered_endog = model.wendog - model.wendog.mean()\n",
        "            return np.dot(centered_endog, centered_endog)\n",
        "\n",
        "    @cache_readonly\n",
        "    def uncentered_tss(self):\n",
        "        \"\"\"\n",
        "        Uncentered sum of squares.\n",
        "\n",
        "        The sum of the squared values of the (whitened) endogenous response\n",
        "        variable.\n",
        "        \"\"\"\n",
        "        wendog = self.model.wendog\n",
        "        return np.dot(wendog, wendog)\n",
        "\n",
        "    @cache_readonly\n",
        "    def ess(self):\n",
        "        \"\"\"\n",
        "        The explained sum of squares.\n",
        "\n",
        "        If a constant is present, the centered total sum of squares minus the\n",
        "        sum of squared residuals. If there is no constant, the uncentered total\n",
        "        sum of squares is used.\n",
        "        \"\"\"\n",
        "\n",
        "        if self.k_constant:\n",
        "            return self.centered_tss - self.ssr\n",
        "        else:\n",
        "            return self.uncentered_tss - self.ssr\n",
        "\n",
        "    @cache_readonly\n",
        "    def rsquared(self):\n",
        "        \"\"\"\n",
        "        R-squared of the model.\n",
        "\n",
        "        This is defined here as 1 - `ssr`/`centered_tss` if the constant is\n",
        "        included in the model and 1 - `ssr`/`uncentered_tss` if the constant is\n",
        "        omitted.\n",
        "        \"\"\"\n",
        "        if self.k_constant:\n",
        "            return 1 - self.ssr/self.centered_tss\n",
        "        else:\n",
        "            return 1 - self.ssr/self.uncentered_tss\n",
        "\n",
        "    @cache_readonly\n",
        "    def rsquared_adj(self):\n",
        "        \"\"\"\n",
        "        Adjusted R-squared.\n",
        "\n",
        "        This is defined here as 1 - (`nobs`-1)/`df_resid` * (1-`rsquared`)\n",
        "        if a constant is included and 1 - `nobs`/`df_resid` * (1-`rsquared`) if\n",
        "        no constant is included.\n",
        "        \"\"\"\n",
        "        return 1 - (np.divide(self.nobs - self.k_constant, self.df_resid)\n",
        "                    * (1 - self.rsquared))\n",
        "\n",
        "    @cache_readonly\n",
        "    def mse_model(self):\n",
        "        \"\"\"\n",
        "        Mean squared error the model.\n",
        "\n",
        "        The explained sum of squares divided by the model degrees of freedom.\n",
        "        \"\"\"\n",
        "        if np.all(self.df_model == 0.0):\n",
        "            return np.full_like(self.ess, np.nan)\n",
        "        return self.ess/self.df_model\n",
        "\n",
        "    @cache_readonly\n",
        "    def mse_resid(self):\n",
        "        \"\"\"\n",
        "        Mean squared error of the residuals.\n",
        "\n",
        "        The sum of squared residuals divided by the residual degrees of\n",
        "        freedom.\n",
        "        \"\"\"\n",
        "        if np.all(self.df_resid == 0.0):\n",
        "            return np.full_like(self.ssr, np.nan)\n",
        "        return self.ssr/self.df_resid\n",
        "\n",
        "    @cache_readonly\n",
        "    def mse_total(self):\n",
        "        \"\"\"\n",
        "        Total mean squared error.\n",
        "\n",
        "        The uncentered total sum of squares divided by the number of\n",
        "        observations.\n",
        "        \"\"\"\n",
        "        if np.all(self.df_resid + self.df_model == 0.0):\n",
        "            return np.full_like(self.centered_tss, np.nan)\n",
        "        if self.k_constant:\n",
        "            return self.centered_tss / (self.df_resid + self.df_model)\n",
        "        else:\n",
        "            return self.uncentered_tss / (self.df_resid + self.df_model)\n",
        "\n",
        "    @cache_readonly\n",
        "    def fvalue(self):\n",
        "        \"\"\"\n",
        "        F-statistic of the fully specified model.\n",
        "\n",
        "        Calculated as the mean squared error of the model divided by the mean\n",
        "        squared error of the residuals if the nonrobust covariance is used.\n",
        "        Otherwise computed using a Wald-like quadratic form that tests whether\n",
        "        all coefficients (excluding the constant) are zero.\n",
        "        \"\"\"\n",
        "        if hasattr(self, 'cov_type') and self.cov_type != 'nonrobust':\n",
        "            # with heteroscedasticity or correlation robustness\n",
        "            k_params = self.normalized_cov_params.shape[0]\n",
        "            mat = np.eye(k_params)\n",
        "            const_idx = self.model.data.const_idx\n",
        "            # TODO: What if model includes implicit constant, e.g. all\n",
        "            #       dummies but no constant regressor?\n",
        "            # TODO: Restats as LM test by projecting orthogonalizing\n",
        "            #       to constant?\n",
        "            if self.model.data.k_constant == 1:\n",
        "                # if constant is implicit, return nan see #2444\n",
        "                if const_idx is None:\n",
        "                    return np.nan\n",
        "\n",
        "                idx = lrange(k_params)\n",
        "                idx.pop(const_idx)\n",
        "                mat = mat[idx]  # remove constant\n",
        "                if mat.size == 0:  # see  #3642\n",
        "                    return np.nan\n",
        "            ft = self.f_test(mat)\n",
        "            # using backdoor to set another attribute that we already have\n",
        "            self._cache['f_pvalue'] = ft.pvalue\n",
        "            return ft.fvalue\n",
        "        else:\n",
        "            # for standard homoscedastic case\n",
        "            return self.mse_model/self.mse_resid\n",
        "\n",
        "    @cache_readonly\n",
        "    def f_pvalue(self):\n",
        "        \"\"\"The p-value of the F-statistic.\"\"\"\n",
        "        # Special case for df_model 0\n",
        "        if self.df_model == 0:\n",
        "            return np.full_like(self.fvalue, np.nan)\n",
        "        return stats.f.sf(self.fvalue, self.df_model, self.df_resid)\n",
        "\n",
        "    @cache_readonly\n",
        "    def bse(self):\n",
        "        \"\"\"The standard errors of the parameter estimates.\"\"\"\n",
        "        return np.sqrt(np.diag(self.cov_params()))\n",
        "\n",
        "    @cache_readonly\n",
        "    def aic(self):\n",
        "        r\"\"\"\n",
        "        Akaike's information criteria.\n",
        "\n",
        "        For a model with a constant :math:`-2llf + 2(df\\_model + 1)`. For a\n",
        "        model without a constant :math:`-2llf + 2(df\\_model)`.\n",
        "        \"\"\"\n",
        "        return -2 * self.llf + 2 * (self.df_model + self.k_constant)\n",
        "\n",
        "    @cache_readonly\n",
        "    def bic(self):\n",
        "        r\"\"\"\n",
        "        Bayes' information criteria.\n",
        "\n",
        "        For a model with a constant :math:`-2llf + \\log(n)(df\\_model+1)`.\n",
        "        For a model without a constant :math:`-2llf + \\log(n)(df\\_model)`.\n",
        "        \"\"\"\n",
        "        return (-2 * self.llf + np.log(self.nobs) * (self.df_model +\n",
        "                                                     self.k_constant))\n",
        "\n",
        "    @cache_readonly\n",
        "    def eigenvals(self):\n",
        "        \"\"\"\n",
        "        Return eigenvalues sorted in decreasing order.\n",
        "        \"\"\"\n",
        "        if self._wexog_singular_values is not None:\n",
        "            eigvals = self._wexog_singular_values ** 2\n",
        "        else:\n",
        "            eigvals = np.linalg.linalg.eigvalsh(np.dot(self.model.wexog.T,\n",
        "                                                       self.model.wexog))\n",
        "        return np.sort(eigvals)[::-1]\n",
        "\n",
        "    @cache_readonly\n",
        "    def condition_number(self):\n",
        "        \"\"\"\n",
        "        Return condition number of exogenous matrix.\n",
        "\n",
        "        Calculated as ratio of largest to smallest eigenvalue.\n",
        "        \"\"\"\n",
        "        eigvals = self.eigenvals\n",
        "        return np.sqrt(eigvals[0]/eigvals[-1])\n",
        "\n",
        "    # TODO: make these properties reset bse\n",
        "    def _HCCM(self, scale):\n",
        "        H = np.dot(self.model.pinv_wexog,\n",
        "                   scale[:, None] * self.model.pinv_wexog.T)\n",
        "        return H\n",
        "\n",
        "    @cache_readonly\n",
        "    def cov_HC0(self):\n",
        "        \"\"\"\n",
        "        Heteroscedasticity robust covariance matrix. See HC0_se.\n",
        "        \"\"\"\n",
        "        self.het_scale = self.wresid**2\n",
        "        cov_HC0 = self._HCCM(self.het_scale)\n",
        "        return cov_HC0\n",
        "\n",
        "    @cache_readonly\n",
        "    def cov_HC1(self):\n",
        "        \"\"\"\n",
        "        Heteroscedasticity robust covariance matrix. See HC1_se.\n",
        "        \"\"\"\n",
        "        self.het_scale = self.nobs/(self.df_resid)*(self.wresid**2)\n",
        "        cov_HC1 = self._HCCM(self.het_scale)\n",
        "        return cov_HC1\n",
        "\n",
        "    @cache_readonly\n",
        "    def cov_HC2(self):\n",
        "        \"\"\"\n",
        "        Heteroscedasticity robust covariance matrix. See HC2_se.\n",
        "        \"\"\"\n",
        "        # probably could be optimized\n",
        "        wexog = self.model.wexog\n",
        "        h = np.diag(wexog @ self.normalized_cov_params @ wexog.T)\n",
        "        self.het_scale = self.wresid**2/(1-h)\n",
        "        cov_HC2 = self._HCCM(self.het_scale)\n",
        "        return cov_HC2\n",
        "\n",
        "    @cache_readonly\n",
        "    def cov_HC3(self):\n",
        "        \"\"\"\n",
        "        Heteroscedasticity robust covariance matrix. See HC3_se.\n",
        "        \"\"\"\n",
        "        wexog = self.model.wexog\n",
        "        h = np.diag(wexog @ self.normalized_cov_params @ wexog.T)\n",
        "        self.het_scale = (self.wresid / (1 - h))**2\n",
        "        cov_HC3 = self._HCCM(self.het_scale)\n",
        "        return cov_HC3\n",
        "\n",
        "    @cache_readonly\n",
        "    def HC0_se(self):\n",
        "        \"\"\"\n",
        "        White's (1980) heteroskedasticity robust standard errors.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        Defined as sqrt(diag(X.T X)^(-1)X.T diag(e_i^(2)) X(X.T X)^(-1)\n",
        "        where e_i = resid[i].\n",
        "\n",
        "        When HC0_se or cov_HC0 is called the RegressionResults instance will\n",
        "        then have another attribute `het_scale`, which is in this case is just\n",
        "        resid**2.\n",
        "        \"\"\"\n",
        "        return np.sqrt(np.diag(self.cov_HC0))\n",
        "\n",
        "    @cache_readonly\n",
        "    def HC1_se(self):\n",
        "        \"\"\"\n",
        "        MacKinnon and White's (1985) heteroskedasticity robust standard errors.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        Defined as sqrt(diag(n/(n-p)*HC_0).\n",
        "\n",
        "        When HC1_se or cov_HC1 is called the RegressionResults instance will\n",
        "        then have another attribute `het_scale`, which is in this case is\n",
        "        n/(n-p)*resid**2.\n",
        "        \"\"\"\n",
        "        return np.sqrt(np.diag(self.cov_HC1))\n",
        "\n",
        "    @cache_readonly\n",
        "    def HC2_se(self):\n",
        "        \"\"\"\n",
        "        MacKinnon and White's (1985) heteroskedasticity robust standard errors.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        Defined as (X.T X)^(-1)X.T diag(e_i^(2)/(1-h_ii)) X(X.T X)^(-1)\n",
        "        where h_ii = x_i(X.T X)^(-1)x_i.T\n",
        "\n",
        "        When HC2_se or cov_HC2 is called the RegressionResults instance will\n",
        "        then have another attribute `het_scale`, which is in this case is\n",
        "        resid^(2)/(1-h_ii).\n",
        "        \"\"\"\n",
        "        return np.sqrt(np.diag(self.cov_HC2))\n",
        "\n",
        "    @cache_readonly\n",
        "    def HC3_se(self):\n",
        "        \"\"\"\n",
        "        MacKinnon and White's (1985) heteroskedasticity robust standard errors.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        Defined as (X.T X)^(-1)X.T diag(e_i^(2)/(1-h_ii)^(2)) X(X.T X)^(-1)\n",
        "        where h_ii = x_i(X.T X)^(-1)x_i.T.\n",
        "\n",
        "        When HC3_se or cov_HC3 is called the RegressionResults instance will\n",
        "        then have another attribute `het_scale`, which is in this case is\n",
        "        resid^(2)/(1-h_ii)^(2).\n",
        "        \"\"\"\n",
        "        return np.sqrt(np.diag(self.cov_HC3))\n",
        "\n",
        "    @cache_readonly\n",
        "    def resid_pearson(self):\n",
        "        \"\"\"\n",
        "        Residuals, normalized to have unit variance.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        array_like\n",
        "            The array `wresid` normalized by the sqrt of the scale to have\n",
        "            unit variance.\n",
        "        \"\"\"\n",
        "\n",
        "        if not hasattr(self, 'resid'):\n",
        "            raise ValueError('Method requires residuals.')\n",
        "        eps = np.finfo(self.wresid.dtype).eps\n",
        "        if np.sqrt(self.scale) < 10 * eps * self.model.endog.mean():\n",
        "            # do not divide if scale is zero close to numerical precision\n",
        "            from warnings import warn\n",
        "            warn(\"All residuals are 0, cannot compute normed residuals.\",\n",
        "                 RuntimeWarning)\n",
        "            return self.wresid\n",
        "        else:\n",
        "            return self.wresid / np.sqrt(self.scale)\n",
        "\n",
        "    def _is_nested(self, restricted):\n",
        "        \"\"\"\n",
        "        Parameters\n",
        "        ----------\n",
        "        restricted : Result instance\n",
        "            The restricted model is assumed to be nested in the current\n",
        "            model. The result instance of the restricted model is required to\n",
        "            have two attributes, residual sum of squares, `ssr`, residual\n",
        "            degrees of freedom, `df_resid`.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        nested : bool\n",
        "            True if nested, otherwise false\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        A most nests another model if the regressors in the smaller\n",
        "        model are spanned by the regressors in the larger model and\n",
        "        the regressand is identical.\n",
        "        \"\"\"\n",
        "\n",
        "        if self.model.nobs != restricted.model.nobs:\n",
        "            return False\n",
        "\n",
        "        full_rank = self.model.rank\n",
        "        restricted_rank = restricted.model.rank\n",
        "        if full_rank <= restricted_rank:\n",
        "            return False\n",
        "\n",
        "        restricted_exog = restricted.model.wexog\n",
        "        full_wresid = self.wresid\n",
        "\n",
        "        scores = restricted_exog * full_wresid[:, None]\n",
        "        score_l2 = np.sqrt(np.mean(scores.mean(0) ** 2))\n",
        "        # TODO: Could be improved, and may fail depending on scale of\n",
        "        # regressors\n",
        "        return np.allclose(score_l2, 0)\n",
        "\n",
        "    def compare_lm_test(self, restricted, demean=True, use_lr=False):\n",
        "        \"\"\"\n",
        "        Use Lagrange Multiplier test to test a set of linear restrictions.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        restricted : Result instance\n",
        "            The restricted model is assumed to be nested in the\n",
        "            current model. The result instance of the restricted model\n",
        "            is required to have two attributes, residual sum of\n",
        "            squares, `ssr`, residual degrees of freedom, `df_resid`.\n",
        "        demean : bool\n",
        "            Flag indicating whether the demean the scores based on the\n",
        "            residuals from the restricted model.  If True, the covariance of\n",
        "            the scores are used and the LM test is identical to the large\n",
        "            sample version of the LR test.\n",
        "        use_lr : bool\n",
        "            A flag indicating whether to estimate the covariance of the model\n",
        "            scores using the unrestricted model. Setting the to True improves\n",
        "            the power of the test.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        lm_value : float\n",
        "            The test statistic which has a chi2 distributed.\n",
        "        p_value : float\n",
        "            The p-value of the test statistic.\n",
        "        df_diff : int\n",
        "            The degrees of freedom of the restriction, i.e. difference in df\n",
        "            between models.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        The LM test examines whether the scores from the restricted model are\n",
        "        0. If the null is true, and the restrictions are valid, then the\n",
        "        parameters of the restricted model should be close to the minimum of\n",
        "        the sum of squared errors, and so the scores should be close to zero,\n",
        "        on average.\n",
        "        \"\"\"\n",
        "        import statsmodels.stats.sandwich_covariance as sw\n",
        "        from numpy.linalg import inv\n",
        "\n",
        "        if not self._is_nested(restricted):\n",
        "            raise ValueError(\"Restricted model is not nested by full model.\")\n",
        "\n",
        "        wresid = restricted.wresid\n",
        "        wexog = self.model.wexog\n",
        "        scores = wexog * wresid[:, None]\n",
        "\n",
        "        n = self.nobs\n",
        "        df_full = self.df_resid\n",
        "        df_restr = restricted.df_resid\n",
        "        df_diff = (df_restr - df_full)\n",
        "\n",
        "        s = scores.mean(axis=0)\n",
        "        if use_lr:\n",
        "            scores = wexog * self.wresid[:, None]\n",
        "            demean = False\n",
        "\n",
        "        if demean:\n",
        "            scores = scores - scores.mean(0)[None, :]\n",
        "        # Form matters here.  If homoskedastics can be sigma^2 (X'X)^-1\n",
        "        # If Heteroskedastic then the form below is fine\n",
        "        # If HAC then need to use HAC\n",
        "        # If Cluster, should use cluster\n",
        "\n",
        "        cov_type = getattr(self, 'cov_type', 'nonrobust')\n",
        "        if cov_type == 'nonrobust':\n",
        "            sigma2 = np.mean(wresid**2)\n",
        "            xpx = np.dot(wexog.T, wexog) / n\n",
        "            s_inv = inv(sigma2 * xpx)\n",
        "        elif cov_type in ('HC0', 'HC1', 'HC2', 'HC3'):\n",
        "            s_inv = inv(np.dot(scores.T, scores) / n)\n",
        "        elif cov_type == 'HAC':\n",
        "            maxlags = self.cov_kwds['maxlags']\n",
        "            s_inv = inv(sw.S_hac_simple(scores, maxlags) / n)\n",
        "        elif cov_type == 'cluster':\n",
        "            # cluster robust standard errors\n",
        "            groups = self.cov_kwds['groups']\n",
        "            # TODO: Might need demean option in S_crosssection by group?\n",
        "            s_inv = inv(sw.S_crosssection(scores, groups))\n",
        "        else:\n",
        "            raise ValueError('Only nonrobust, HC, HAC and cluster are ' +\n",
        "                             'currently connected')\n",
        "\n",
        "        lm_value = n * (s @ s_inv @ s.T)\n",
        "        p_value = stats.chi2.sf(lm_value, df_diff)\n",
        "        return lm_value, p_value, df_diff\n",
        "\n",
        "\n",
        "    def compare_f_test(self, restricted):\n",
        "        \"\"\"\n",
        "        Use F test to test whether restricted model is correct.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        restricted : Result instance\n",
        "            The restricted model is assumed to be nested in the\n",
        "            current model. The result instance of the restricted model\n",
        "            is required to have two attributes, residual sum of\n",
        "            squares, `ssr`, residual degrees of freedom, `df_resid`.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        f_value : float\n",
        "            The test statistic which has an F distribution.\n",
        "        p_value : float\n",
        "            The p-value of the test statistic.\n",
        "        df_diff : int\n",
        "            The degrees of freedom of the restriction, i.e. difference in\n",
        "            df between models.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        See mailing list discussion October 17,\n",
        "\n",
        "        This test compares the residual sum of squares of the two\n",
        "        models.  This is not a valid test, if there is unspecified\n",
        "        heteroscedasticity or correlation. This method will issue a\n",
        "        warning if this is detected but still return the results under\n",
        "        the assumption of homoscedasticity and no autocorrelation\n",
        "        (sphericity).\n",
        "        \"\"\"\n",
        "\n",
        "        has_robust1 = getattr(self, 'cov_type', 'nonrobust') != 'nonrobust'\n",
        "        has_robust2 = (getattr(restricted, 'cov_type', 'nonrobust') !=\n",
        "                       'nonrobust')\n",
        "\n",
        "        if has_robust1 or has_robust2:\n",
        "            warnings.warn('F test for comparison is likely invalid with ' +\n",
        "                          'robust covariance, proceeding anyway',\n",
        "                          InvalidTestWarning)\n",
        "\n",
        "        ssr_full = self.ssr\n",
        "        ssr_restr = restricted.ssr\n",
        "        df_full = self.df_resid\n",
        "        df_restr = restricted.df_resid\n",
        "\n",
        "        df_diff = (df_restr - df_full)\n",
        "        f_value = (ssr_restr - ssr_full) / df_diff / ssr_full * df_full\n",
        "        p_value = stats.f.sf(f_value, df_diff, df_full)\n",
        "        return f_value, p_value, df_diff\n",
        "\n",
        "\n",
        "    def compare_lr_test(self, restricted, large_sample=False):\n",
        "        \"\"\"\n",
        "        Likelihood ratio test to test whether restricted model is correct.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        restricted : Result instance\n",
        "            The restricted model is assumed to be nested in the current model.\n",
        "            The result instance of the restricted model is required to have two\n",
        "            attributes, residual sum of squares, `ssr`, residual degrees of\n",
        "            freedom, `df_resid`.\n",
        "\n",
        "        large_sample : bool\n",
        "            Flag indicating whether to use a heteroskedasticity robust version\n",
        "            of the LR test, which is a modified LM test.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        lr_stat : float\n",
        "            The likelihood ratio which is chisquare distributed with df_diff\n",
        "            degrees of freedom.\n",
        "        p_value : float\n",
        "            The p-value of the test statistic.\n",
        "        df_diff : int\n",
        "            The degrees of freedom of the restriction, i.e. difference in df\n",
        "            between models.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        The exact likelihood ratio is valid for homoskedastic data,\n",
        "        and is defined as\n",
        "\n",
        "        .. math:: D=-2\\\\log\\\\left(\\\\frac{\\\\mathcal{L}_{null}}\n",
        "           {\\\\mathcal{L}_{alternative}}\\\\right)\n",
        "\n",
        "        where :math:`\\\\mathcal{L}` is the likelihood of the\n",
        "        model. With :math:`D` distributed as chisquare with df equal\n",
        "        to difference in number of parameters or equivalently\n",
        "        difference in residual degrees of freedom.\n",
        "\n",
        "        The large sample version of the likelihood ratio is defined as\n",
        "\n",
        "        .. math:: D=n s^{\\\\prime}S^{-1}s\n",
        "\n",
        "        where :math:`s=n^{-1}\\\\sum_{i=1}^{n} s_{i}`\n",
        "\n",
        "        .. math:: s_{i} = x_{i,alternative} \\\\epsilon_{i,null}\n",
        "\n",
        "        is the average score of the model evaluated using the\n",
        "        residuals from null model and the regressors from the\n",
        "        alternative model and :math:`S` is the covariance of the\n",
        "        scores, :math:`s_{i}`.  The covariance of the scores is\n",
        "        estimated using the same estimator as in the alternative\n",
        "        model.\n",
        "\n",
        "        This test compares the loglikelihood of the two models.  This\n",
        "        may not be a valid test, if there is unspecified\n",
        "        heteroscedasticity or correlation. This method will issue a\n",
        "        warning if this is detected but still return the results\n",
        "        without taking unspecified heteroscedasticity or correlation\n",
        "        into account.\n",
        "\n",
        "        This test compares the loglikelihood of the two models.  This\n",
        "        may not be a valid test, if there is unspecified\n",
        "        heteroscedasticity or correlation. This method will issue a\n",
        "        warning if this is detected but still return the results\n",
        "        without taking unspecified heteroscedasticity or correlation\n",
        "        into account.\n",
        "\n",
        "        is the average score of the model evaluated using the\n",
        "        residuals from null model and the regressors from the\n",
        "        alternative model and :math:`S` is the covariance of the\n",
        "        scores, :math:`s_{i}`.  The covariance of the scores is\n",
        "        estimated using the same estimator as in the alternative\n",
        "        model.\n",
        "        \"\"\"\n",
        "        # TODO: put into separate function, needs tests\n",
        "\n",
        "        # See mailing list discussion October 17,\n",
        "\n",
        "        if large_sample:\n",
        "            return self.compare_lm_test(restricted, use_lr=True)\n",
        "\n",
        "        has_robust1 = (getattr(self, 'cov_type', 'nonrobust') != 'nonrobust')\n",
        "        has_robust2 = (\n",
        "            getattr(restricted, 'cov_type', 'nonrobust') != 'nonrobust')\n",
        "\n",
        "        if has_robust1 or has_robust2:\n",
        "            warnings.warn('Likelihood Ratio test is likely invalid with ' +\n",
        "                          'robust covariance, proceeding anyway',\n",
        "                          InvalidTestWarning)\n",
        "\n",
        "        llf_full = self.llf\n",
        "        llf_restr = restricted.llf\n",
        "        df_full = self.df_resid\n",
        "        df_restr = restricted.df_resid\n",
        "\n",
        "        lrdf = (df_restr - df_full)\n",
        "        lrstat = -2*(llf_restr - llf_full)\n",
        "        lr_pvalue = stats.chi2.sf(lrstat, lrdf)\n",
        "\n",
        "        return lrstat, lr_pvalue, lrdf\n",
        "\n",
        "\n",
        "    def get_robustcov_results(self, cov_type='HC1', use_t=None, **kwargs):\n",
        "        \"\"\"\n",
        "        Create new results instance with robust covariance as default.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        cov_type : str\n",
        "            The type of robust sandwich estimator to use. See Notes below.\n",
        "        use_t : bool\n",
        "            If true, then the t distribution is used for inference.\n",
        "            If false, then the normal distribution is used.\n",
        "            If `use_t` is None, then an appropriate default is used, which is\n",
        "            `True` if the cov_type is nonrobust, and `False` in all other\n",
        "            cases.\n",
        "        **kwargs\n",
        "            Required or optional arguments for robust covariance calculation.\n",
        "            See Notes below.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        RegressionResults\n",
        "            This method creates a new results instance with the\n",
        "            requested robust covariance as the default covariance of\n",
        "            the parameters.  Inferential statistics like p-values and\n",
        "            hypothesis tests will be based on this covariance matrix.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        The following covariance types and required or optional arguments are\n",
        "        currently available:\n",
        "\n",
        "        - 'fixed scale' and optional keyword argument 'scale' which uses\n",
        "            a predefined scale estimate with default equal to one.\n",
        "        - 'HC0', 'HC1', 'HC2', 'HC3' and no keyword arguments:\n",
        "            heteroscedasticity robust covariance\n",
        "        - 'HAC' and keywords\n",
        "\n",
        "            - `maxlag` integer (required) : number of lags to use\n",
        "            - `kernel` callable or str (optional) : kernel\n",
        "                  currently available kernels are ['bartlett', 'uniform'],\n",
        "                  default is Bartlett\n",
        "            - `use_correction` bool (optional) : If true, use small sample\n",
        "                  correction\n",
        "\n",
        "        - 'cluster' and required keyword `groups`, integer group indicator\n",
        "\n",
        "            - `groups` array_like, integer (required) :\n",
        "                  index of clusters or groups\n",
        "            - `use_correction` bool (optional) :\n",
        "                  If True the sandwich covariance is calculated with a small\n",
        "                  sample correction.\n",
        "                  If False the sandwich covariance is calculated without\n",
        "                  small sample correction.\n",
        "            - `df_correction` bool (optional)\n",
        "                  If True (default), then the degrees of freedom for the\n",
        "                  inferential statistics and hypothesis tests, such as\n",
        "                  pvalues, f_pvalue, conf_int, and t_test and f_test, are\n",
        "                  based on the number of groups minus one instead of the\n",
        "                  total number of observations minus the number of explanatory\n",
        "                  variables. `df_resid` of the results instance is adjusted.\n",
        "                  If False, then `df_resid` of the results instance is not\n",
        "                  adjusted.\n",
        "\n",
        "        - 'hac-groupsum' Driscoll and Kraay, heteroscedasticity and\n",
        "            autocorrelation robust standard errors in panel data\n",
        "            keywords\n",
        "\n",
        "            - `time` array_like (required) : index of time periods\n",
        "            - `maxlag` integer (required) : number of lags to use\n",
        "            - `kernel` callable or str (optional). The available kernels\n",
        "              are ['bartlett', 'uniform']. The default is Bartlett.\n",
        "            - `use_correction` False or string in ['hac', 'cluster'] (optional).\n",
        "              If False the the sandwich covariance is calculated without small\n",
        "              sample correction. If `use_correction = 'cluster'` (default),\n",
        "              then the same small sample correction as in the case of\n",
        "              `covtype='cluster'` is used.\n",
        "            - `df_correction` bool (optional) The adjustment to df_resid, see\n",
        "              cov_type 'cluster' above\n",
        "              # TODO: we need more options here\n",
        "\n",
        "        - 'hac-panel' heteroscedasticity and autocorrelation robust standard\n",
        "            errors in panel data.\n",
        "            The data needs to be sorted in this case, the time series\n",
        "            for each panel unit or cluster need to be stacked. The\n",
        "            membership to a timeseries of an individual or group can\n",
        "            be either specified by group indicators or by increasing\n",
        "            time periods.\n",
        "\n",
        "            keywords\n",
        "\n",
        "            - either `groups` or `time` : array_like (required)\n",
        "              `groups` : indicator for groups\n",
        "              `time` : index of time periods\n",
        "            - `maxlag` integer (required) : number of lags to use\n",
        "            - `kernel` callable or str (optional)\n",
        "                  currently available kernels are ['bartlett', 'uniform'],\n",
        "                  default is Bartlett\n",
        "            - `use_correction` False or string in ['hac', 'cluster'] (optional)\n",
        "                  If False the sandwich covariance is calculated without\n",
        "                  small sample correction.\n",
        "            - `df_correction` bool (optional)\n",
        "                  adjustment to df_resid, see cov_type 'cluster' above\n",
        "                  # TODO: we need more options here\n",
        "\n",
        "        Reminder:\n",
        "        `use_correction` in \"hac-groupsum\" and \"hac-panel\" is not bool,\n",
        "        needs to be in [False, 'hac', 'cluster']\n",
        "\n",
        "        TODO: Currently there is no check for extra or misspelled keywords,\n",
        "        except in the case of cov_type `HCx`\n",
        "        \"\"\"\n",
        "        import statsmodels.stats.sandwich_covariance as sw\n",
        "        from statsmodels.base.covtype import normalize_cov_type, descriptions\n",
        "\n",
        "        cov_type = normalize_cov_type(cov_type)\n",
        "\n",
        "        if 'kernel' in kwargs:\n",
        "            kwargs['weights_func'] = kwargs.pop('kernel')\n",
        "        if 'weights_func' in kwargs and not callable(kwargs['weights_func']):\n",
        "            kwargs['weights_func'] = sw.kernel_dict[kwargs['weights_func']]\n",
        "\n",
        "        # TODO: make separate function that returns a robust cov plus info\n",
        "        use_self = kwargs.pop('use_self', False)\n",
        "        if use_self:\n",
        "            res = self\n",
        "        else:\n",
        "            res = self.__class__(\n",
        "                self.model, self.params,\n",
        "                normalized_cov_params=self.normalized_cov_params,\n",
        "                scale=self.scale)\n",
        "\n",
        "        res.cov_type = cov_type\n",
        "        # use_t might already be defined by the class, and already set\n",
        "        if use_t is None:\n",
        "            use_t = self.use_t\n",
        "        res.cov_kwds = {'use_t': use_t}  # store for information\n",
        "        res.use_t = use_t\n",
        "\n",
        "        adjust_df = False\n",
        "        if cov_type in ['cluster', 'hac-panel', 'hac-groupsum']:\n",
        "            df_correction = kwargs.get('df_correction', None)\n",
        "            # TODO: check also use_correction, do I need all combinations?\n",
        "            if df_correction is not False:  # i.e. in [None, True]:\n",
        "                # user did not explicitely set it to False\n",
        "                adjust_df = True\n",
        "\n",
        "        res.cov_kwds['adjust_df'] = adjust_df\n",
        "\n",
        "        # verify and set kwargs, and calculate cov\n",
        "        # TODO: this should be outsourced in a function so we can reuse it in\n",
        "        #       other models\n",
        "        # TODO: make it DRYer   repeated code for checking kwargs\n",
        "        if cov_type in ['fixed scale', 'fixed_scale']:\n",
        "            res.cov_kwds['description'] = descriptions['fixed_scale']\n",
        "\n",
        "            res.cov_kwds['scale'] = scale = kwargs.get('scale', 1.)\n",
        "            res.cov_params_default = scale * res.normalized_cov_params\n",
        "        elif cov_type.upper() in ('HC0', 'HC1', 'HC2', 'HC3'):\n",
        "            if kwargs:\n",
        "                raise ValueError('heteroscedasticity robust covariance '\n",
        "                                 'does not use keywords')\n",
        "            res.cov_kwds['description'] = descriptions[cov_type.upper()]\n",
        "            res.cov_params_default = getattr(self, 'cov_' + cov_type.upper())\n",
        "        elif cov_type.lower() == 'hac':\n",
        "            # TODO: check if required, default in cov_hac_simple\n",
        "            maxlags = kwargs['maxlags']\n",
        "            res.cov_kwds['maxlags'] = maxlags\n",
        "            weights_func = kwargs.get('weights_func', sw.weights_bartlett)\n",
        "            res.cov_kwds['weights_func'] = weights_func\n",
        "            use_correction = kwargs.get('use_correction', False)\n",
        "            res.cov_kwds['use_correction'] = use_correction\n",
        "            res.cov_kwds['description'] = descriptions['HAC'].format(\n",
        "                maxlags=maxlags,\n",
        "                correction=['without', 'with'][use_correction])\n",
        "\n",
        "            res.cov_params_default = sw.cov_hac_simple(\n",
        "                self, nlags=maxlags, weights_func=weights_func,\n",
        "                use_correction=use_correction)\n",
        "        elif cov_type.lower() == 'cluster':\n",
        "            # cluster robust standard errors, one- or two-way\n",
        "            groups = kwargs['groups']\n",
        "            if not hasattr(groups, 'shape'):\n",
        "                groups = np.asarray(groups).T\n",
        "\n",
        "            if groups.ndim >= 2:\n",
        "                groups = groups.squeeze()\n",
        "\n",
        "            res.cov_kwds['groups'] = groups\n",
        "            use_correction = kwargs.get('use_correction', True)\n",
        "            res.cov_kwds['use_correction'] = use_correction\n",
        "            if groups.ndim == 1:\n",
        "                if adjust_df:\n",
        "                    # need to find number of groups\n",
        "                    # duplicate work\n",
        "                    self.n_groups = n_groups = len(np.unique(groups))\n",
        "                res.cov_params_default = sw.cov_cluster(\n",
        "                    self, groups, use_correction=use_correction)\n",
        "\n",
        "            elif groups.ndim == 2:\n",
        "                if hasattr(groups, 'values'):\n",
        "                    groups = groups.values\n",
        "\n",
        "                if adjust_df:\n",
        "                    # need to find number of groups\n",
        "                    # duplicate work\n",
        "                    n_groups0 = len(np.unique(groups[:, 0]))\n",
        "                    n_groups1 = len(np.unique(groups[:, 1]))\n",
        "                    self.n_groups = (n_groups0, n_groups1)\n",
        "                    n_groups = min(n_groups0, n_groups1)  # use for adjust_df\n",
        "\n",
        "                # Note: sw.cov_cluster_2groups has 3 returns\n",
        "                res.cov_params_default = sw.cov_cluster_2groups(\n",
        "                    self, groups, use_correction=use_correction)[0]\n",
        "            else:\n",
        "                raise ValueError('only two groups are supported')\n",
        "            res.cov_kwds['description'] = descriptions['cluster']\n",
        "\n",
        "        elif cov_type.lower() == 'hac-panel':\n",
        "            # cluster robust standard errors\n",
        "            res.cov_kwds['time'] = time = kwargs.get('time', None)\n",
        "            res.cov_kwds['groups'] = groups = kwargs.get('groups', None)\n",
        "            # TODO: nlags is currently required\n",
        "            # nlags = kwargs.get('nlags', True)\n",
        "            # res.cov_kwds['nlags'] = nlags\n",
        "            # TODO: `nlags` or `maxlags`\n",
        "            res.cov_kwds['maxlags'] = maxlags = kwargs['maxlags']\n",
        "            use_correction = kwargs.get('use_correction', 'hac')\n",
        "            res.cov_kwds['use_correction'] = use_correction\n",
        "            weights_func = kwargs.get('weights_func', sw.weights_bartlett)\n",
        "            res.cov_kwds['weights_func'] = weights_func\n",
        "            if groups is not None:\n",
        "                groups = np.asarray(groups)\n",
        "                tt = (np.nonzero(groups[:-1] != groups[1:])[0] + 1).tolist()\n",
        "                nobs_ = len(groups)\n",
        "            elif time is not None:\n",
        "                time = np.asarray(time)\n",
        "                # TODO: clumsy time index in cov_nw_panel\n",
        "                tt = (np.nonzero(time[1:] < time[:-1])[0] + 1).tolist()\n",
        "                nobs_ = len(time)\n",
        "            else:\n",
        "                raise ValueError('either time or groups needs to be given')\n",
        "            groupidx = lzip([0] + tt, tt + [nobs_])\n",
        "            self.n_groups = n_groups = len(groupidx)\n",
        "            res.cov_params_default = sw.cov_nw_panel(self, maxlags, groupidx,\n",
        "                                                     weights_func=weights_func,\n",
        "                                                     use_correction=use_correction)\n",
        "            res.cov_kwds['description'] = descriptions['HAC-Panel']\n",
        "\n",
        "        elif cov_type.lower() == 'hac-groupsum':\n",
        "            # Driscoll-Kraay standard errors\n",
        "            res.cov_kwds['time'] = time = kwargs['time']\n",
        "            # TODO: nlags is currently required\n",
        "            # nlags = kwargs.get('nlags', True)\n",
        "            # res.cov_kwds['nlags'] = nlags\n",
        "            # TODO: `nlags` or `maxlags`\n",
        "            res.cov_kwds['maxlags'] = maxlags = kwargs['maxlags']\n",
        "            use_correction = kwargs.get('use_correction', 'cluster')\n",
        "            res.cov_kwds['use_correction'] = use_correction\n",
        "            weights_func = kwargs.get('weights_func', sw.weights_bartlett)\n",
        "            res.cov_kwds['weights_func'] = weights_func\n",
        "            if adjust_df:\n",
        "                # need to find number of groups\n",
        "                tt = (np.nonzero(time[1:] < time[:-1])[0] + 1)\n",
        "                self.n_groups = n_groups = len(tt) + 1\n",
        "            res.cov_params_default = sw.cov_nw_groupsum(\n",
        "                self, maxlags, time, weights_func=weights_func,\n",
        "                use_correction=use_correction)\n",
        "            res.cov_kwds['description'] = descriptions['HAC-Groupsum']\n",
        "        else:\n",
        "            raise ValueError('cov_type not recognized. See docstring for ' +\n",
        "                             'available options and spelling')\n",
        "\n",
        "        if adjust_df:\n",
        "            # Note: df_resid is used for scale and others, add new attribute\n",
        "            res.df_resid_inference = n_groups - 1\n",
        "\n",
        "        return res\n",
        "\n",
        "\n",
        "    @Appender(pred.get_prediction.__doc__)\n",
        "    def get_prediction(self, exog=None, transform=True, weights=None,\n",
        "                       row_labels=None, **kwargs):\n",
        "\n",
        "        return pred.get_prediction(\n",
        "            self, exog=exog, transform=transform, weights=weights,\n",
        "            row_labels=row_labels, **kwargs)\n",
        "\n",
        "\n",
        "    def summary(self, yname=None, xname=None, title=None, alpha=.05):\n",
        "        \"\"\"\n",
        "        Summarize the Regression Results.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        yname : str, optional\n",
        "            Name of endogenous (response) variable. The Default is `y`.\n",
        "        xname : list[str], optional\n",
        "            Names for the exogenous variables. Default is `var_##` for ## in\n",
        "            the number of regressors. Must match the number of parameters\n",
        "            in the model.\n",
        "        title : str, optional\n",
        "            Title for the top table. If not None, then this replaces the\n",
        "            default title.\n",
        "        alpha : float\n",
        "            The significance level for the confidence intervals.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        Summary\n",
        "            Instance holding the summary tables and text, which can be printed\n",
        "            or converted to various output formats.\n",
        "\n",
        "        See Also\n",
        "        --------\n",
        "        statsmodels.iolib.summary.Summary : A class that holds summary results.\n",
        "        \"\"\"\n",
        "        from statsmodels.stats.stattools import (\n",
        "            jarque_bera, omni_normtest, durbin_watson)\n",
        "\n",
        "        jb, jbpv, skew, kurtosis = jarque_bera(self.wresid)\n",
        "        omni, omnipv = omni_normtest(self.wresid)\n",
        "\n",
        "        eigvals = self.eigenvals\n",
        "        condno = self.condition_number\n",
        "\n",
        "        # TODO: Avoid adding attributes in non-__init__\n",
        "        self.diagn = dict(jb=jb, jbpv=jbpv, skew=skew, kurtosis=kurtosis,\n",
        "                          omni=omni, omnipv=omnipv, condno=condno,\n",
        "                          mineigval=eigvals[-1])\n",
        "\n",
        "        # TODO not used yet\n",
        "        # diagn_left_header = ['Models stats']\n",
        "        # diagn_right_header = ['Residual stats']\n",
        "\n",
        "        # TODO: requiring list/iterable is a bit annoying\n",
        "        #   need more control over formatting\n",
        "        # TODO: default do not work if it's not identically spelled\n",
        "\n",
        "        top_left = [('Dep. Variable:', None),\n",
        "                    ('Model:', None),\n",
        "                    ('Method:', ['Least Squares']),\n",
        "                    ('Date:', None),\n",
        "                    ('Time:', None),\n",
        "                    ('No. Observations:', None),\n",
        "                    ('Df Residuals:', None),\n",
        "                    ('Df Model:', None),\n",
        "                    ]\n",
        "\n",
        "        if hasattr(self, 'cov_type'):\n",
        "            top_left.append(('Covariance Type:', [self.cov_type]))\n",
        "\n",
        "        rsquared_type = '' if self.k_constant else ' (uncentered)'\n",
        "        top_right = [('R-squared' + rsquared_type + ':',\n",
        "                      [\"%#8.3f\" % self.rsquared]),\n",
        "                     ('Adj. R-squared' + rsquared_type + ':',\n",
        "                      [\"%#8.3f\" % self.rsquared_adj]),\n",
        "                     ('F-statistic:', [\"%#8.4g\" % self.fvalue]),\n",
        "                     ('Prob (F-statistic):', [\"%#6.3g\" % self.f_pvalue]),\n",
        "                     ('Log-Likelihood:', None),\n",
        "                     ('AIC:', [\"%#8.4g\" % self.aic]),\n",
        "                     ('BIC:', [\"%#8.4g\" % self.bic])\n",
        "                     ]\n",
        "\n",
        "        diagn_left = [('Omnibus:', [\"%#6.3f\" % omni]),\n",
        "                      ('Prob(Omnibus):', [\"%#6.3f\" % omnipv]),\n",
        "                      ('Skew:', [\"%#6.3f\" % skew]),\n",
        "                      ('Kurtosis:', [\"%#6.3f\" % kurtosis])\n",
        "                      ]\n",
        "\n",
        "        diagn_right = [('Durbin-Watson:',\n",
        "                        [\"%#8.3f\" % durbin_watson(self.wresid)]\n",
        "                        ),\n",
        "                       ('Jarque-Bera (JB):', [\"%#8.3f\" % jb]),\n",
        "                       ('Prob(JB):', [\"%#8.3g\" % jbpv]),\n",
        "                       ('Cond. No.', [\"%#8.3g\" % condno])\n",
        "                       ]\n",
        "\n",
        "        if title is None:\n",
        "            title = self.model.__class__.__name__ + ' ' + \"Regression Results\"\n",
        "\n",
        "        # create summary table instance\n",
        "        from statsmodels.iolib.summary import Summary\n",
        "        smry = Summary()\n",
        "        smry.add_table_2cols(self, gleft=top_left, gright=top_right,\n",
        "                             yname=yname, xname=xname, title=title)\n",
        "        smry.add_table_params(self, yname=yname, xname=xname, alpha=alpha,\n",
        "                              use_t=self.use_t)\n",
        "\n",
        "        smry.add_table_2cols(self, gleft=diagn_left, gright=diagn_right,\n",
        "                             yname=yname, xname=xname,\n",
        "                             title=\"\")\n",
        "\n",
        "        # add warnings/notes, added to text format only\n",
        "        etext = []\n",
        "        if hasattr(self, 'cov_type'):\n",
        "            etext.append(self.cov_kwds['description'])\n",
        "        if self.model.exog.shape[0] < self.model.exog.shape[1]:\n",
        "            wstr = \"The input rank is higher than the number of observations.\"\n",
        "            etext.append(wstr)\n",
        "        if eigvals[-1] < 1e-10:\n",
        "            wstr = \"The smallest eigenvalue is %6.3g. This might indicate \"\n",
        "            wstr += \"that there are\\n\"\n",
        "            wstr += \"strong multicollinearity problems or that the design \"\n",
        "            wstr += \"matrix is singular.\"\n",
        "            wstr = wstr % eigvals[-1]\n",
        "            etext.append(wstr)\n",
        "        elif condno > 1000:  # TODO: what is recommended?\n",
        "            wstr = \"The condition number is large, %6.3g. This might \"\n",
        "            wstr += \"indicate that there are\\n\"\n",
        "            wstr += \"strong multicollinearity or other numerical \"\n",
        "            wstr += \"problems.\"\n",
        "            wstr = wstr % condno\n",
        "            etext.append(wstr)\n",
        "\n",
        "        if etext:\n",
        "            etext = [\"[{0}] {1}\".format(i + 1, text)\n",
        "                     for i, text in enumerate(etext)]\n",
        "            etext.insert(0, \"Warnings:\")\n",
        "            smry.add_extra_txt(etext)\n",
        "\n",
        "        return smry\n",
        "\n",
        "\n",
        "    def summary2(self, yname=None, xname=None, title=None, alpha=.05,\n",
        "                 float_format=\"%.4f\"):\n",
        "        \"\"\"\n",
        "        Experimental summary function to summarize the regression results.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        yname : str\n",
        "            The name of the dependent variable (optional).\n",
        "        xname : list[str], optional\n",
        "            Names for the exogenous variables. Default is `var_##` for ## in\n",
        "            the number of regressors. Must match the number of parameters\n",
        "            in the model.\n",
        "        title : str, optional\n",
        "            Title for the top table. If not None, then this replaces the\n",
        "            default title.\n",
        "        alpha : float\n",
        "            The significance level for the confidence intervals.\n",
        "        float_format : str\n",
        "            The format for floats in parameters summary.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        Summary\n",
        "            Instance holding the summary tables and text, which can be printed\n",
        "            or converted to various output formats.\n",
        "\n",
        "        See Also\n",
        "        --------\n",
        "        statsmodels.iolib.summary2.Summary\n",
        "            A class that holds summary results.\n",
        "        \"\"\"\n",
        "        # Diagnostics\n",
        "        from statsmodels.stats.stattools import (jarque_bera,\n",
        "                                                 omni_normtest,\n",
        "                                                 durbin_watson)\n",
        "\n",
        "        from collections import OrderedDict\n",
        "        jb, jbpv, skew, kurtosis = jarque_bera(self.wresid)\n",
        "        omni, omnipv = omni_normtest(self.wresid)\n",
        "        dw = durbin_watson(self.wresid)\n",
        "        eigvals = self.eigenvals\n",
        "        condno = self.condition_number\n",
        "        eigvals = np.sort(eigvals)  # in increasing order\n",
        "        diagnostic = OrderedDict([\n",
        "            ('Omnibus:',  \"%.3f\" % omni),\n",
        "            ('Prob(Omnibus):', \"%.3f\" % omnipv),\n",
        "            ('Skew:', \"%.3f\" % skew),\n",
        "            ('Kurtosis:', \"%.3f\" % kurtosis),\n",
        "            ('Durbin-Watson:', \"%.3f\" % dw),\n",
        "            ('Jarque-Bera (JB):', \"%.3f\" % jb),\n",
        "            ('Prob(JB):', \"%.3f\" % jbpv),\n",
        "            ('Condition No.:', \"%.0f\" % condno)\n",
        "            ])\n",
        "\n",
        "        # Summary\n",
        "        from statsmodels.iolib import summary2\n",
        "        smry = summary2.Summary()\n",
        "        smry.add_base(results=self, alpha=alpha, float_format=float_format,\n",
        "                      xname=xname, yname=yname, title=title)\n",
        "        smry.add_dict(diagnostic)\n",
        "\n",
        "        # Warnings\n",
        "        if eigvals[-1] < 1e-10:\n",
        "            warn = \"The smallest eigenvalue is %6.3g. This might indicate that\\\n",
        "            there are strong multicollinearity problems or that the design\\\n",
        "            matrix is singular.\" % eigvals[-1]\n",
        "            smry.add_text(warn)\n",
        "        if condno > 1000:\n",
        "            warn = \"* The condition number is large (%.g). This might indicate \\\n",
        "            strong multicollinearity or other numerical problems.\" % condno\n",
        "            smry.add_text(warn)\n",
        "\n",
        "        return smry\n",
        "\n",
        "\n",
        "\n",
        "class OLSResults(RegressionResults):\n",
        "    \"\"\"\n",
        "    Results class for for an OLS model.\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    model : RegressionModel\n",
        "        The regression model instance.\n",
        "    params : ndarray\n",
        "        The estimated parameters.\n",
        "    normalized_cov_params : ndarray\n",
        "        The normalized covariance parameters.\n",
        "    scale : float\n",
        "        The estimated scale of the residuals.\n",
        "    cov_type : str\n",
        "        The covariance estimator used in the results.\n",
        "    cov_kwds : dict\n",
        "        Additional keywords used in the covariance specification.\n",
        "    use_t : bool\n",
        "        Flag indicating to use the Student's t in inference.\n",
        "    **kwargs\n",
        "        Additional keyword arguments used to initialize the results.\n",
        "\n",
        "    See Also\n",
        "    --------\n",
        "    RegressionResults\n",
        "        Results store for WLS and GLW models.\n",
        "\n",
        "    Notes\n",
        "    -----\n",
        "    Most of the methods and attributes are inherited from RegressionResults.\n",
        "    The special methods that are only available for OLS are:\n",
        "\n",
        "    - get_influence\n",
        "    - outlier_test\n",
        "    - el_test\n",
        "    - conf_int_el\n",
        "    \"\"\"\n",
        "\n",
        "    def get_influence(self):\n",
        "        \"\"\"\n",
        "        Calculate influence and outlier measures.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        OLSInfluence\n",
        "            The instance containing methods to calculate the main influence and\n",
        "            outlier measures for the OLS regression.\n",
        "\n",
        "        See Also\n",
        "        --------\n",
        "        statsmodels.stats.outliers_influence.OLSInfluence\n",
        "            A class that exposes methods to examine observation influence.\n",
        "        \"\"\"\n",
        "        from statsmodels.stats.outliers_influence import OLSInfluence\n",
        "        return OLSInfluence(self)\n",
        "\n",
        "\n",
        "    def outlier_test(self, method='bonf', alpha=.05, labels=None,\n",
        "                     order=False, cutoff=None):\n",
        "        \"\"\"\n",
        "        Test observations for outliers according to method.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        method : str\n",
        "            The method to use in the outlier test.  Must be one of:\n",
        "\n",
        "            - `bonferroni` : one-step correction\n",
        "            - `sidak` : one-step correction\n",
        "            - `holm-sidak` :\n",
        "            - `holm` :\n",
        "            - `simes-hochberg` :\n",
        "            - `hommel` :\n",
        "            - `fdr_bh` : Benjamini/Hochberg\n",
        "            - `fdr_by` : Benjamini/Yekutieli\n",
        "\n",
        "            See `statsmodels.stats.multitest.multipletests` for details.\n",
        "        alpha : float\n",
        "            The familywise error rate (FWER).\n",
        "        labels : None or array_like\n",
        "            If `labels` is not None, then it will be used as index to the\n",
        "            returned pandas DataFrame. See also Returns below.\n",
        "        order : bool\n",
        "            Whether or not to order the results by the absolute value of the\n",
        "            studentized residuals. If labels are provided they will also be\n",
        "            sorted.\n",
        "        cutoff : None or float in [0, 1]\n",
        "            If cutoff is not None, then the return only includes observations\n",
        "            with multiple testing corrected p-values strictly below the cutoff.\n",
        "            The returned array or dataframe can be empty if t.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        array_like\n",
        "            Returns either an ndarray or a DataFrame if labels is not None.\n",
        "            Will attempt to get labels from model_results if available. The\n",
        "            columns are the Studentized residuals, the unadjusted p-value,\n",
        "            and the corrected p-value according to method.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        The unadjusted p-value is stats.t.sf(abs(resid), df) where\n",
        "        df = df_resid - 1.\n",
        "        \"\"\"\n",
        "        from statsmodels.stats.outliers_influence import outlier_test\n",
        "        return outlier_test(self, method, alpha, labels=labels,\n",
        "                            order=order, cutoff=cutoff)\n",
        "\n",
        "\n",
        "    def el_test(self, b0_vals, param_nums, return_weights=0, ret_params=0,\n",
        "                method='nm', stochastic_exog=1):\n",
        "        \"\"\"\n",
        "        Test single or joint hypotheses using Empirical Likelihood.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        b0_vals : 1darray\n",
        "            The hypothesized value of the parameter to be tested.\n",
        "        param_nums : 1darray\n",
        "            The parameter number to be tested.\n",
        "        return_weights : bool\n",
        "            If true, returns the weights that optimize the likelihood\n",
        "            ratio at b0_vals. The default is False.\n",
        "        ret_params : bool\n",
        "            If true, returns the parameter vector that maximizes the likelihood\n",
        "            ratio at b0_vals.  Also returns the weights.  The default is False.\n",
        "        method : str\n",
        "            Can either be 'nm' for Nelder-Mead or 'powell' for Powell.  The\n",
        "            optimization method that optimizes over nuisance parameters.\n",
        "            The default is 'nm'.\n",
        "        stochastic_exog : bool\n",
        "            When True, the exogenous variables are assumed to be stochastic.\n",
        "            When the regressors are nonstochastic, moment conditions are\n",
        "            placed on the exogenous variables.  Confidence intervals for\n",
        "            stochastic regressors are at least as large as non-stochastic\n",
        "            regressors. The default is True.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        tuple\n",
        "            The p-value and -2 times the log-likelihood ratio for the\n",
        "            hypothesized values.\n",
        "\n",
        "        Examples\n",
        "        --------\n",
        "        >>> import statsmodels.api as sm\n",
        "        >>> data = sm.datasets.stackloss.load(as_pandas=False)\n",
        "        >>> endog = data.endog\n",
        "        >>> exog = sm.add_constant(data.exog)\n",
        "        >>> model = sm.OLS(endog, exog)\n",
        "        >>> fitted = model.fit()\n",
        "        >>> fitted.params\n",
        "        >>> array([-39.91967442,   0.7156402 ,   1.29528612,  -0.15212252])\n",
        "        >>> fitted.rsquared\n",
        "        >>> 0.91357690446068196\n",
        "        >>> # Test that the slope on the first variable is 0\n",
        "        >>> fitted.el_test([0], [1])\n",
        "        >>> (27.248146353888796, 1.7894660442330235e-07)\n",
        "        \"\"\"\n",
        "        params = np.copy(self.params)\n",
        "        opt_fun_inst = _ELRegOpts()  # to store weights\n",
        "        if len(param_nums) == len(params):\n",
        "            llr = opt_fun_inst._opt_nuis_regress(\n",
        "                [],\n",
        "                param_nums=param_nums,\n",
        "                endog=self.model.endog,\n",
        "                exog=self.model.exog,\n",
        "                nobs=self.model.nobs,\n",
        "                nvar=self.model.exog.shape[1],\n",
        "                params=params,\n",
        "                b0_vals=b0_vals,\n",
        "                stochastic_exog=stochastic_exog)\n",
        "            pval = 1 - stats.chi2.cdf(llr, len(param_nums))\n",
        "            if return_weights:\n",
        "                return llr, pval, opt_fun_inst.new_weights\n",
        "            else:\n",
        "                return llr, pval\n",
        "        x0 = np.delete(params, param_nums)\n",
        "        args = (param_nums, self.model.endog, self.model.exog,\n",
        "                self.model.nobs, self.model.exog.shape[1], params,\n",
        "                b0_vals, stochastic_exog)\n",
        "        if method == 'nm':\n",
        "            llr = optimize.fmin(opt_fun_inst._opt_nuis_regress, x0,\n",
        "                                maxfun=10000, maxiter=10000, full_output=1,\n",
        "                                disp=0, args=args)[1]\n",
        "        if method == 'powell':\n",
        "            llr = optimize.fmin_powell(opt_fun_inst._opt_nuis_regress, x0,\n",
        "                                       full_output=1, disp=0,\n",
        "                                       args=args)[1]\n",
        "\n",
        "        pval = 1 - stats.chi2.cdf(llr, len(param_nums))\n",
        "        if ret_params:\n",
        "            return llr, pval, opt_fun_inst.new_weights, opt_fun_inst.new_params\n",
        "        elif return_weights:\n",
        "            return llr, pval, opt_fun_inst.new_weights\n",
        "        else:\n",
        "            return llr, pval\n",
        "\n",
        "\n",
        "    def conf_int_el(self, param_num, sig=.05, upper_bound=None,\n",
        "                    lower_bound=None, method='nm', stochastic_exog=True):\n",
        "        \"\"\"\n",
        "        Compute the confidence interval using Empirical Likelihood.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        param_num : float\n",
        "            The parameter for which the confidence interval is desired.\n",
        "        sig : float\n",
        "            The significance level.  Default is 0.05.\n",
        "        upper_bound : float\n",
        "            The maximum value the upper limit can be.  Default is the\n",
        "            99.9% confidence value under OLS assumptions.\n",
        "        lower_bound : float\n",
        "            The minimum value the lower limit can be.  Default is the 99.9%\n",
        "            confidence value under OLS assumptions.\n",
        "        method : str\n",
        "            Can either be 'nm' for Nelder-Mead or 'powell' for Powell.  The\n",
        "            optimization method that optimizes over nuisance parameters.\n",
        "            The default is 'nm'.\n",
        "        stochastic_exog : bool\n",
        "            When True, the exogenous variables are assumed to be stochastic.\n",
        "            When the regressors are nonstochastic, moment conditions are\n",
        "            placed on the exogenous variables.  Confidence intervals for\n",
        "            stochastic regressors are at least as large as non-stochastic\n",
        "            regressors.  The default is True.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        lowerl : float\n",
        "            The lower bound of the confidence interval.\n",
        "        upperl : float\n",
        "            The upper bound of the confidence interval.\n",
        "\n",
        "        See Also\n",
        "        --------\n",
        "        el_test : Test parameters using Empirical Likelihood.\n",
        "\n",
        "        Notes\n",
        "        -----\n",
        "        This function uses brentq to find the value of beta where\n",
        "        test_beta([beta], param_num)[1] is equal to the critical value.\n",
        "\n",
        "        The function returns the results of each iteration of brentq at each\n",
        "        value of beta.\n",
        "\n",
        "        The current function value of the last printed optimization should be\n",
        "        the critical value at the desired significance level. For alpha=.05,\n",
        "        the value is 3.841459.\n",
        "\n",
        "        To ensure optimization terminated successfully, it is suggested to do\n",
        "        el_test([lower_limit], [param_num]).\n",
        "\n",
        "        If the optimization does not terminate successfully, consider switching\n",
        "        optimization algorithms.\n",
        "\n",
        "        If optimization is still not successful, try changing the values of\n",
        "        start_int_params.  If the current function value repeatedly jumps\n",
        "        from a number between 0 and the critical value and a very large number\n",
        "        (>50), the starting parameters of the interior minimization need\n",
        "        to be changed.\n",
        "        \"\"\"\n",
        "        r0 = stats.chi2.ppf(1 - sig, 1)\n",
        "        if upper_bound is None:\n",
        "            upper_bound = self.conf_int(.01)[param_num][1]\n",
        "        if lower_bound is None:\n",
        "            lower_bound = self.conf_int(.01)[param_num][0]\n",
        "\n",
        "        def f(b0):\n",
        "            return self.el_test(np.array([b0]), np.array([param_num]),\n",
        "                                method=method,\n",
        "                                stochastic_exog=stochastic_exog)[0] - r0\n",
        "\n",
        "        lowerl = optimize.brenth(f, lower_bound,\n",
        "                                 self.params[param_num])\n",
        "        upperl = optimize.brenth(f, self.params[param_num],\n",
        "                                 upper_bound)\n",
        "        #  ^ Seems to be faster than brentq in most cases\n",
        "        return (lowerl, upperl)\n",
        "\n",
        "\n",
        "\n",
        "class RegressionResultsWrapper(wrap.ResultsWrapper):\n",
        "\n",
        "    _attrs = {\n",
        "        'chisq': 'columns',\n",
        "        'sresid': 'rows',\n",
        "        'weights': 'rows',\n",
        "        'wresid': 'rows',\n",
        "        'bcov_unscaled': 'cov',\n",
        "        'bcov_scaled': 'cov',\n",
        "        'HC0_se': 'columns',\n",
        "        'HC1_se': 'columns',\n",
        "        'HC2_se': 'columns',\n",
        "        'HC3_se': 'columns',\n",
        "        'norm_resid': 'rows',\n",
        "    }\n",
        "\n",
        "    _wrap_attrs = wrap.union_dicts(base.LikelihoodResultsWrapper._attrs,\n",
        "                                   _attrs)\n",
        "\n",
        "    _methods = {}\n",
        "\n",
        "    _wrap_methods = wrap.union_dicts(\n",
        "                        base.LikelihoodResultsWrapper._wrap_methods,\n",
        "                        _methods)\n",
        "\n",
        "\n",
        "wrap.populate_wrapper(RegressionResultsWrapper,\n",
        "                      RegressionResults)"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "error",
          "ename": "ImportError",
          "evalue": "ignored",
          "traceback": [
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[0;31mImportError\u001b[0m                               Traceback (most recent call last)",
            "\u001b[0;32m<ipython-input-11-abfb13fb12a2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     29\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     30\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mstatsmodels\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpython\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mlrange\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlzip\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 31\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mstatsmodels\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpandas\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mAppender\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     32\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     33\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mnumpy\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;31mImportError\u001b[0m: cannot import name 'Appender'",
            "",
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0;32m\nNOTE: If your import is failing due to a missing package, you can\nmanually install dependencies using either !pip or !apt.\n\nTo view examples of installing some common dependencies, click the\n\"Open Examples\" button below.\n\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Pwenv_wcRiPg",
        "colab_type": "text"
      },
      "source": [
        "LOG-LIKELIHOOD ESTIMATORS: <br>\n",
        "\n",
        "$-\\frac{n}{2}\\log\\left(\\left(Y-\\hat{Y}\\right)^{\\prime}              \\left(Y-\\hat{Y}\\right)\\right)\n",
        "                  -\\frac{n}{2}\\left(1+\\log\\left(\\frac{2\\pi}{n}\\right)\\right)\n",
        "                  -\\frac{1}{2}\\log\\left(\\left|\\Sigma\\right|\\right)$\n"
      ]
    }
  ]
}