Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1# TODO: Determine which tests are valid for GLSAR, and under what conditions 

2# TODO: Fix issue with constant and GLS 

3# TODO: GLS: add options Iterative GLS, for iterative fgls if sigma is None 

4# TODO: GLS: default if sigma is none should be two-step GLS 

5# TODO: Check nesting when performing model based tests, lr, wald, lm 

6""" 

7This module implements standard regression models: 

8 

9Generalized Least Squares (GLS) 

10Ordinary Least Squares (OLS) 

11Weighted Least Squares (WLS) 

12Generalized Least Squares with autoregressive error terms GLSAR(p) 

13 

14Models are specified with an endogenous response variable and an 

15exogenous design matrix and are fit using their `fit` method. 

16 

17Subclasses that have more complicated covariance matrices 

18should write over the 'whiten' method as the fit method 

19prewhitens the response by calling 'whiten'. 

20 

21General reference for regression models: 

22 

23D. C. Montgomery and E.A. Peck. "Introduction to Linear Regression 

24 Analysis." 2nd. Ed., Wiley, 1992. 

25 

26Econometrics references for regression models: 

27 

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

29 2004. 

30 

31W. Green. "Econometric Analysis," 5th ed., Pearson, 2003. 

32""" 

33 

34 

35from statsmodels.compat.python import lrange, lzip 

36from statsmodels.compat.pandas import Appender 

37 

38import numpy as np 

39from scipy.linalg import toeplitz 

40from scipy import stats 

41from scipy import optimize 

42 

43from statsmodels.tools.tools import pinv_extended 

44from statsmodels.tools.decorators import (cache_readonly, 

45 cache_writable) 

46import statsmodels.base.model as base 

47import statsmodels.base.wrapper as wrap 

48from statsmodels.emplike.elregress import _ELRegOpts 

49import warnings 

50from statsmodels.tools.sm_exceptions import InvalidTestWarning 

51 

52# need import in module instead of lazily to copy `__doc__` 

53from statsmodels.regression._prediction import PredictionResults 

54from . import _prediction as pred 

55 

56__docformat__ = 'restructuredtext en' 

57 

58__all__ = ['GLS', 'WLS', 'OLS', 'GLSAR', 'PredictionResults', 

59 'RegressionResultsWrapper'] 

60 

61 

62_fit_regularized_doc =\ 

63 r""" 

64 Return a regularized fit to a linear regression model. 

65 

66 Parameters 

67 ---------- 

68 method : str 

69 Either 'elastic_net' or 'sqrt_lasso'. 

70 alpha : scalar or array_like 

71 The penalty weight. If a scalar, the same penalty weight 

72 applies to all variables in the model. If a vector, it 

73 must have the same length as `params`, and contains a 

74 penalty weight for each coefficient. 

75 L1_wt : scalar 

76 The fraction of the penalty given to the L1 penalty term. 

77 Must be between 0 and 1 (inclusive). If 0, the fit is a 

78 ridge fit, if 1 it is a lasso fit. 

79 start_params : array_like 

80 Starting values for ``params``. 

81 profile_scale : bool 

82 If True the penalized fit is computed using the profile 

83 (concentrated) log-likelihood for the Gaussian model. 

84 Otherwise the fit uses the residual sum of squares. 

85 refit : bool 

86 If True, the model is refit using only the variables that 

87 have non-zero coefficients in the regularized fit. The 

88 refitted model is not regularized. 

89 **kwargs 

90 Additional keyword arguments that contain information used when 

91 constructing a model using the formula interface. 

92 

93 Returns 

94 ------- 

95 statsmodels.base.elastic_net.RegularizedResults 

96 The regularized results. 

97 

98 Notes 

99 ----- 

100 The elastic net uses a combination of L1 and L2 penalties. 

101 The implementation closely follows the glmnet package in R. 

102 

103 The function that is minimized is: 

104 

105 .. math:: 

106 

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

108 

109 where RSS is the usual regression sum of squares, n is the 

110 sample size, and :math:`|*|_1` and :math:`|*|_2` are the L1 and L2 

111 norms. 

112 

113 For WLS and GLS, the RSS is calculated using the whitened endog and 

114 exog data. 

115 

116 Post-estimation results are based on the same data used to 

117 select variables, hence may be subject to overfitting biases. 

118 

119 The elastic_net method uses the following keyword arguments: 

120 

121 maxiter : int 

122 Maximum number of iterations 

123 cnvrg_tol : float 

124 Convergence threshold for line searches 

125 zero_tol : float 

126 Coefficients below this threshold are treated as zero. 

127 

128 The square root lasso approach is a variation of the Lasso 

129 that is largely self-tuning (the optimal tuning parameter 

130 does not depend on the standard deviation of the regression 

131 errors). If the errors are Gaussian, the tuning parameter 

132 can be taken to be 

133 

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

135 

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

137 

138 The square root lasso uses the following keyword arguments: 

139 

140 zero_tol : float 

141 Coefficients below this threshold are treated as zero. 

142 

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

144 lasso. 

145 

146 References 

147 ---------- 

148 .. [*] Friedman, Hastie, Tibshirani (2008). Regularization paths for 

149 generalized linear models via coordinate descent. Journal of 

150 Statistical Software 33(1), 1-22 Feb 2010. 

151 

152 .. [*] A Belloni, V Chernozhukov, L Wang (2011). Square-root Lasso: 

153 pivotal recovery of sparse signals via conic programming. 

154 Biometrika 98(4), 791-806. https://arxiv.org/pdf/1009.5689.pdf 

155 """ 

156 

157 

158def _get_sigma(sigma, nobs): 

159 """ 

160 Returns sigma (matrix, nobs by nobs) for GLS and the inverse of its 

161 Cholesky decomposition. Handles dimensions and checks integrity. 

162 If sigma is None, returns None, None. Otherwise returns sigma, 

163 cholsigmainv. 

164 """ 

165 if sigma is None: 

166 return None, None 

167 sigma = np.asarray(sigma).squeeze() 

168 if sigma.ndim == 0: 

169 sigma = np.repeat(sigma, nobs) 

170 if sigma.ndim == 1: 

171 if sigma.shape != (nobs,): 

172 raise ValueError("Sigma must be a scalar, 1d of length %s or a 2d " 

173 "array of shape %s x %s" % (nobs, nobs, nobs)) 

174 cholsigmainv = 1/np.sqrt(sigma) 

175 else: 

176 if sigma.shape != (nobs, nobs): 

177 raise ValueError("Sigma must be a scalar, 1d of length %s or a 2d " 

178 "array of shape %s x %s" % (nobs, nobs, nobs)) 

179 cholsigmainv = np.linalg.cholesky(np.linalg.inv(sigma)).T 

180 return sigma, cholsigmainv 

181 

182 

183class RegressionModel(base.LikelihoodModel): 

184 """ 

185 Base class for linear regression models. Should not be directly called. 

186 

187 Intended for subclassing. 

188 """ 

189 def __init__(self, endog, exog, **kwargs): 

190 super(RegressionModel, self).__init__(endog, exog, **kwargs) 

191 self._data_attr.extend(['pinv_wexog', 'wendog', 'wexog', 'weights']) 

192 

193 def initialize(self): 

194 """Initialize model components.""" 

195 self.wexog = self.whiten(self.exog) 

196 self.wendog = self.whiten(self.endog) 

197 # overwrite nobs from class Model: 

198 self.nobs = float(self.wexog.shape[0]) 

199 

200 self._df_model = None 

201 self._df_resid = None 

202 self.rank = None 

203 

204 @property 

205 def df_model(self): 

206 """ 

207 The model degree of freedom. 

208 

209 The dof is defined as the rank of the regressor matrix minus 1 if a 

210 constant is included. 

211 """ 

212 if self._df_model is None: 

213 if self.rank is None: 

214 self.rank = np.linalg.matrix_rank(self.exog) 

215 self._df_model = float(self.rank - self.k_constant) 

216 return self._df_model 

217 

218 @df_model.setter 

219 def df_model(self, value): 

220 self._df_model = value 

221 

222 @property 

223 def df_resid(self): 

224 """ 

225 The residual degree of freedom. 

226 

227 The dof is defined as the number of observations minus the rank of 

228 the regressor matrix. 

229 """ 

230 

231 if self._df_resid is None: 

232 if self.rank is None: 

233 self.rank = np.linalg.matrix_rank(self.exog) 

234 self._df_resid = self.nobs - self.rank 

235 return self._df_resid 

236 

237 @df_resid.setter 

238 def df_resid(self, value): 

239 self._df_resid = value 

240 

241 def whiten(self, x): 

242 """ 

243 Whiten method that must be overwritten by individual models. 

244 

245 Parameters 

246 ---------- 

247 x : array_like 

248 Data to be whitened. 

249 """ 

250 raise NotImplementedError("Subclasses must implement.") 

251 

252 def fit(self, method="pinv", cov_type='nonrobust', cov_kwds=None, 

253 use_t=None, **kwargs): 

254 """ 

255 Full fit of the model. 

256 

257 The results include an estimate of covariance matrix, (whitened) 

258 residuals and an estimate of scale. 

259 

260 Parameters 

261 ---------- 

262 method : str, optional 

263 Can be "pinv", "qr". "pinv" uses the Moore-Penrose pseudoinverse 

264 to solve the least squares problem. "qr" uses the QR 

265 factorization. 

266 cov_type : str, optional 

267 See `regression.linear_model.RegressionResults` for a description 

268 of the available covariance estimators. 

269 cov_kwds : list or None, optional 

270 See `linear_model.RegressionResults.get_robustcov_results` for a 

271 description required keywords for alternative covariance 

272 estimators. 

273 use_t : bool, optional 

274 Flag indicating to use the Student's t distribution when computing 

275 p-values. Default behavior depends on cov_type. See 

276 `linear_model.RegressionResults.get_robustcov_results` for 

277 implementation details. 

278 **kwargs 

279 Additional keyword arguments that contain information used when 

280 constructing a model using the formula interface. 

281 

282 Returns 

283 ------- 

284 RegressionResults 

285 The model estimation results. 

286 

287 See Also 

288 -------- 

289 RegressionResults 

290 The results container. 

291 RegressionResults.get_robustcov_results 

292 A method to change the covariance estimator used when fitting the 

293 model. 

294 

295 Notes 

296 ----- 

297 The fit method uses the pseudoinverse of the design/exogenous variables 

298 to solve the least squares minimization. 

299 """ 

300 if method == "pinv": 

301 if not (hasattr(self, 'pinv_wexog') and 

302 hasattr(self, 'normalized_cov_params') and 

303 hasattr(self, 'rank')): 

304 

305 self.pinv_wexog, singular_values = pinv_extended(self.wexog) 

306 self.normalized_cov_params = np.dot( 

307 self.pinv_wexog, np.transpose(self.pinv_wexog)) 

308 

309 # Cache these singular values for use later. 

310 self.wexog_singular_values = singular_values 

311 self.rank = np.linalg.matrix_rank(np.diag(singular_values)) 

312 

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

314 

315 elif method == "qr": 

316 if not (hasattr(self, 'exog_Q') and 

317 hasattr(self, 'exog_R') and 

318 hasattr(self, 'normalized_cov_params') and 

319 hasattr(self, 'rank')): 

320 Q, R = np.linalg.qr(self.wexog) 

321 self.exog_Q, self.exog_R = Q, R 

322 self.normalized_cov_params = np.linalg.inv(np.dot(R.T, R)) 

323 

324 # Cache singular values from R. 

325 self.wexog_singular_values = np.linalg.svd(R, 0, 0) 

326 self.rank = np.linalg.matrix_rank(R) 

327 else: 

328 Q, R = self.exog_Q, self.exog_R 

329 

330 # used in ANOVA 

331 self.effects = effects = np.dot(Q.T, self.wendog) 

332 beta = np.linalg.solve(R, effects) 

333 else: 

334 raise ValueError('method has to be "pinv" or "qr"') 

335 

336 if self._df_model is None: 

337 self._df_model = float(self.rank - self.k_constant) 

338 if self._df_resid is None: 

339 self.df_resid = self.nobs - self.rank 

340 

341 if isinstance(self, OLS): 

342 lfit = OLSResults( 

343 self, beta, 

344 normalized_cov_params=self.normalized_cov_params, 

345 cov_type=cov_type, cov_kwds=cov_kwds, use_t=use_t) 

346 else: 

347 lfit = RegressionResults( 

348 self, beta, 

349 normalized_cov_params=self.normalized_cov_params, 

350 cov_type=cov_type, cov_kwds=cov_kwds, use_t=use_t, 

351 **kwargs) 

352 return RegressionResultsWrapper(lfit) 

353 

354 def predict(self, params, exog=None): 

355 """ 

356 Return linear predicted values from a design matrix. 

357 

358 Parameters 

359 ---------- 

360 params : array_like 

361 Parameters of a linear model. 

362 exog : array_like, optional 

363 Design / exogenous data. Model exog is used if None. 

364 

365 Returns 

366 ------- 

367 array_like 

368 An array of fitted values. 

369 

370 Notes 

371 ----- 

372 If the model has not yet been fit, params is not optional. 

373 """ 

374 # JP: this does not look correct for GLMAR 

375 # SS: it needs its own predict method 

376 

377 if exog is None: 

378 exog = self.exog 

379 

380 return np.dot(exog, params) 

381 

382 def get_distribution(self, params, scale, exog=None, dist_class=None): 

383 """ 

384 Construct a random number generator for the predictive distribution. 

385 

386 Parameters 

387 ---------- 

388 params : array_like 

389 The model parameters (regression coefficients). 

390 scale : scalar 

391 The variance parameter. 

392 exog : array_like 

393 The predictor variable matrix. 

394 dist_class : class 

395 A random number generator class. Must take 'loc' and 'scale' 

396 as arguments and return a random number generator implementing 

397 an ``rvs`` method for simulating random values. Defaults to normal. 

398 

399 Returns 

400 ------- 

401 gen 

402 Frozen random number generator object with mean and variance 

403 determined by the fitted linear model. Use the ``rvs`` method 

404 to generate random values. 

405 

406 Notes 

407 ----- 

408 Due to the behavior of ``scipy.stats.distributions objects``, 

409 the returned random number generator must be called with 

410 ``gen.rvs(n)`` where ``n`` is the number of observations in 

411 the data set used to fit the model. If any other value is 

412 used for ``n``, misleading results will be produced. 

413 """ 

414 fit = self.predict(params, exog) 

415 if dist_class is None: 

416 from scipy.stats.distributions import norm 

417 dist_class = norm 

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

419 return gen 

420 

421 

422class GLS(RegressionModel): 

423 __doc__ = r""" 

424 Generalized Least Squares 

425 

426 %(params)s 

427 sigma : scalar or array 

428 The array or scalar `sigma` is the weighting matrix of the covariance. 

429 The default is None for no scaling. If `sigma` is a scalar, it is 

430 assumed that `sigma` is an n x n diagonal matrix with the given 

431 scalar, `sigma` as the value of each diagonal element. If `sigma` 

432 is an n-length vector, then `sigma` is assumed to be a diagonal 

433 matrix with the given `sigma` on the diagonal. This should be the 

434 same as WLS. 

435 %(extra_params)s 

436 

437 Attributes 

438 ---------- 

439 pinv_wexog : ndarray 

440 `pinv_wexog` is the p x n Moore-Penrose pseudoinverse of `wexog`. 

441 cholsimgainv : ndarray 

442 The transpose of the Cholesky decomposition of the pseudoinverse. 

443 df_model : float 

444 p - 1, where p is the number of regressors including the intercept. 

445 of freedom. 

446 df_resid : float 

447 Number of observations n less the number of parameters p. 

448 llf : float 

449 The value of the likelihood function of the fitted model. 

450 nobs : float 

451 The number of observations n. 

452 normalized_cov_params : ndarray 

453 p x p array :math:`(X^{T}\Sigma^{-1}X)^{-1}` 

454 results : RegressionResults instance 

455 A property that returns the RegressionResults class if fit. 

456 sigma : ndarray 

457 `sigma` is the n x n covariance structure of the error terms. 

458 wexog : ndarray 

459 Design matrix whitened by `cholsigmainv` 

460 wendog : ndarray 

461 Response variable whitened by `cholsigmainv` 

462 

463 See Also 

464 -------- 

465 WLS : Fit a linear model using Weighted Least Squares. 

466 OLS : Fit a linear model using Ordinary Least Squares. 

467 

468 Notes 

469 ----- 

470 If sigma is a function of the data making one of the regressors 

471 a constant, then the current postestimation statistics will not be correct. 

472 

473 Examples 

474 -------- 

475 >>> import statsmodels.api as sm 

476 >>> data = sm.datasets.longley.load(as_pandas=False) 

477 >>> data.exog = sm.add_constant(data.exog) 

478 >>> ols_resid = sm.OLS(data.endog, data.exog).fit().resid 

479 >>> res_fit = sm.OLS(ols_resid[1:], ols_resid[:-1]).fit() 

480 >>> rho = res_fit.params 

481 

482 `rho` is a consistent estimator of the correlation of the residuals from 

483 an OLS fit of the longley data. It is assumed that this is the true rho 

484 of the AR process data. 

485 

486 >>> from scipy.linalg import toeplitz 

487 >>> order = toeplitz(np.arange(16)) 

488 >>> sigma = rho**order 

489 

490 `sigma` is an n x n matrix of the autocorrelation structure of the 

491 data. 

492 

493 >>> gls_model = sm.GLS(data.endog, data.exog, sigma=sigma) 

494 >>> gls_results = gls_model.fit() 

495 >>> print(gls_results.summary()) 

496 """ % {'params': base._model_params_doc, 

497 'extra_params': base._missing_param_doc + base._extra_param_doc} 

498 

499 def __init__(self, endog, exog, sigma=None, missing='none', hasconst=None, 

500 **kwargs): 

501 # TODO: add options igls, for iterative fgls if sigma is None 

502 # TODO: default if sigma is none should be two-step GLS 

503 sigma, cholsigmainv = _get_sigma(sigma, len(endog)) 

504 

505 super(GLS, self).__init__(endog, exog, missing=missing, 

506 hasconst=hasconst, sigma=sigma, 

507 cholsigmainv=cholsigmainv, **kwargs) 

508 

509 # store attribute names for data arrays 

510 self._data_attr.extend(['sigma', 'cholsigmainv']) 

511 

512 def whiten(self, x): 

513 """ 

514 GLS whiten method. 

515 

516 Parameters 

517 ---------- 

518 x : array_like 

519 Data to be whitened. 

520 

521 Returns 

522 ------- 

523 ndarray 

524 The value np.dot(cholsigmainv,X). 

525 

526 See Also 

527 -------- 

528 GLS : Fit a linear model using Generalized Least Squares. 

529 """ 

530 x = np.asarray(x) 

531 if self.sigma is None or self.sigma.shape == (): 

532 return x 

533 elif self.sigma.ndim == 1: 

534 if x.ndim == 1: 

535 return x * self.cholsigmainv 

536 else: 

537 return x * self.cholsigmainv[:, None] 

538 else: 

539 return np.dot(self.cholsigmainv, x) 

540 

541 def loglike(self, params): 

542 r""" 

543 Compute the value of the Gaussian log-likelihood function at params. 

544 

545 Given the whitened design matrix, the log-likelihood is evaluated 

546 at the parameter vector `params` for the dependent variable `endog`. 

547 

548 Parameters 

549 ---------- 

550 params : array_like 

551 The model parameters. 

552 

553 Returns 

554 ------- 

555 float 

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

557 

558 Notes 

559 ----- 

560 The log-likelihood function for the normal distribution is 

561 

562 .. math:: -\frac{n}{2}\log\left(\left(Y-\hat{Y}\right)^{\prime} 

563 \left(Y-\hat{Y}\right)\right) 

564 -\frac{n}{2}\left(1+\log\left(\frac{2\pi}{n}\right)\right) 

565 -\frac{1}{2}\log\left(\left|\Sigma\right|\right) 

566 

567 Y and Y-hat are whitened. 

568 """ 

569 # TODO: combine this with OLS/WLS loglike and add _det_sigma argument 

570 nobs2 = self.nobs / 2.0 

571 SSR = np.sum((self.wendog - np.dot(self.wexog, params))**2, axis=0) 

572 llf = -np.log(SSR) * nobs2 # concentrated likelihood 

573 llf -= (1+np.log(np.pi/nobs2))*nobs2 # with likelihood constant 

574 if np.any(self.sigma): 

575 # FIXME: robust-enough check? unneeded if _det_sigma gets defined 

576 if self.sigma.ndim == 2: 

577 det = np.linalg.slogdet(self.sigma) 

578 llf -= .5*det[1] 

579 else: 

580 llf -= 0.5*np.sum(np.log(self.sigma)) 

581 # with error covariance matrix 

582 return llf 

583 

584 def hessian_factor(self, params, scale=None, observed=True): 

585 """ 

586 Compute weights for calculating Hessian. 

587 

588 Parameters 

589 ---------- 

590 params : ndarray 

591 The parameter at which Hessian is evaluated. 

592 scale : None or float 

593 If scale is None, then the default scale will be calculated. 

594 Default scale is defined by `self.scaletype` and set in fit. 

595 If scale is not None, then it is used as a fixed scale. 

596 observed : bool 

597 If True, then the observed Hessian is returned. If false then the 

598 expected information matrix is returned. 

599 

600 Returns 

601 ------- 

602 ndarray 

603 A 1d weight vector used in the calculation of the Hessian. 

604 The hessian is obtained by `(exog.T * hessian_factor).dot(exog)`. 

605 """ 

606 

607 if self.sigma is None or self.sigma.shape == (): 

608 return np.ones(self.exog.shape[0]) 

609 elif self.sigma.ndim == 1: 

610 return self.cholsigmainv 

611 else: 

612 return np.diag(self.cholsigmainv) 

613 

614 @Appender(_fit_regularized_doc) 

615 def fit_regularized(self, method="elastic_net", alpha=0., 

616 L1_wt=1., start_params=None, profile_scale=False, 

617 refit=False, **kwargs): 

618 # Need to adjust since RSS/n term in elastic net uses nominal 

619 # n in denominator 

620 if self.sigma is not None: 

621 alpha = alpha * np.sum(1 / np.diag(self.sigma)) / len(self.endog) 

622 

623 rslt = OLS(self.wendog, self.wexog).fit_regularized( 

624 method=method, alpha=alpha, 

625 L1_wt=L1_wt, 

626 start_params=start_params, 

627 profile_scale=profile_scale, 

628 refit=refit, **kwargs) 

629 

630 from statsmodels.base.elastic_net import ( 

631 RegularizedResults, RegularizedResultsWrapper) 

632 rrslt = RegularizedResults(self, rslt.params) 

633 return RegularizedResultsWrapper(rrslt) 

634 

635 

636class WLS(RegressionModel): 

637 __doc__ = """ 

638 Weighted Least Squares 

639 

640 The weights are presumed to be (proportional to) the inverse of 

641 the variance of the observations. That is, if the variables are 

642 to be transformed by 1/sqrt(W) you must supply weights = 1/W. 

643 

644 %(params)s 

645 weights : array_like, optional 

646 A 1d array of weights. If you supply 1/W then the variables are 

647 pre- multiplied by 1/sqrt(W). If no weights are supplied the 

648 default value is 1 and WLS results are the same as OLS. 

649 %(extra_params)s 

650 

651 Attributes 

652 ---------- 

653 weights : ndarray 

654 The stored weights supplied as an argument. 

655 

656 See Also 

657 -------- 

658 GLS : Fit a linear model using Generalized Least Squares. 

659 OLS : Fit a linear model using Ordinary Least Squares. 

660 

661 Notes 

662 ----- 

663 If the weights are a function of the data, then the post estimation 

664 statistics such as fvalue and mse_model might not be correct, as the 

665 package does not yet support no-constant regression. 

666 

667 Examples 

668 -------- 

669 >>> import statsmodels.api as sm 

670 >>> Y = [1,3,4,5,2,3,4] 

671 >>> X = range(1,8) 

672 >>> X = sm.add_constant(X) 

673 >>> wls_model = sm.WLS(Y,X, weights=list(range(1,8))) 

674 >>> results = wls_model.fit() 

675 >>> results.params 

676 array([ 2.91666667, 0.0952381 ]) 

677 >>> results.tvalues 

678 array([ 2.0652652 , 0.35684428]) 

679 >>> print(results.t_test([1, 0])) 

680 <T test: effect=array([ 2.91666667]), sd=array([[ 1.41224801]]), t=array([[ 2.0652652]]), p=array([[ 0.04690139]]), df_denom=5> 

681 >>> print(results.f_test([0, 1])) 

682 <F test: F=array([[ 0.12733784]]), p=[[ 0.73577409]], df_denom=5, df_num=1> 

683 """ % {'params': base._model_params_doc, 

684 'extra_params': base._missing_param_doc + base._extra_param_doc} 

685 

686 def __init__(self, endog, exog, weights=1., missing='none', hasconst=None, 

687 **kwargs): 

688 weights = np.array(weights) 

689 if weights.shape == (): 

690 if (missing == 'drop' and 'missing_idx' in kwargs and 

691 kwargs['missing_idx'] is not None): 

692 # patsy may have truncated endog 

693 weights = np.repeat(weights, len(kwargs['missing_idx'])) 

694 else: 

695 weights = np.repeat(weights, len(endog)) 

696 # handle case that endog might be of len == 1 

697 if len(weights) == 1: 

698 weights = np.array([weights.squeeze()]) 

699 else: 

700 weights = weights.squeeze() 

701 super(WLS, self).__init__(endog, exog, missing=missing, 

702 weights=weights, hasconst=hasconst, **kwargs) 

703 nobs = self.exog.shape[0] 

704 weights = self.weights 

705 # Experimental normalization of weights 

706 weights = weights / np.sum(weights) * nobs 

707 if weights.size != nobs and weights.shape[0] != nobs: 

708 raise ValueError('Weights must be scalar or same length as design') 

709 

710 def whiten(self, x): 

711 """ 

712 Whitener for WLS model, multiplies each column by sqrt(self.weights). 

713 

714 Parameters 

715 ---------- 

716 x : array_like 

717 Data to be whitened. 

718 

719 Returns 

720 ------- 

721 array_like 

722 The whitened values sqrt(weights)*X. 

723 """ 

724 

725 x = np.asarray(x) 

726 if x.ndim == 1: 

727 return x * np.sqrt(self.weights) 

728 elif x.ndim == 2: 

729 return np.sqrt(self.weights)[:, None] * x 

730 

731 def loglike(self, params): 

732 r""" 

733 Compute the value of the gaussian log-likelihood function at params. 

734 

735 Given the whitened design matrix, the log-likelihood is evaluated 

736 at the parameter vector `params` for the dependent variable `Y`. 

737 

738 Parameters 

739 ---------- 

740 params : array_like 

741 The parameter estimates. 

742 

743 Returns 

744 ------- 

745 float 

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

747 

748 Notes 

749 -------- 

750 .. math:: -\frac{n}{2}\log SSR 

751 -\frac{n}{2}\left(1+\log\left(\frac{2\pi}{n}\right)\right) 

752 -\frac{1}{2}\log\left(\left|W\right|\right) 

753 

754 where :math:`W` is a diagonal weight matrix matrix and 

755 :math:`SSR=\left(Y-\hat{Y}\right)^\prime W \left(Y-\hat{Y}\right)` is 

756 the sum of the squared weighted residuals. 

757 """ 

758 nobs2 = self.nobs / 2.0 

759 SSR = np.sum((self.wendog - np.dot(self.wexog, params))**2, axis=0) 

760 llf = -np.log(SSR) * nobs2 # concentrated likelihood 

761 llf -= (1+np.log(np.pi/nobs2))*nobs2 # with constant 

762 llf += 0.5 * np.sum(np.log(self.weights)) 

763 return llf 

764 

765 def hessian_factor(self, params, scale=None, observed=True): 

766 """ 

767 Compute the weights for calculating the Hessian. 

768 

769 Parameters 

770 ---------- 

771 params : ndarray 

772 The parameter at which Hessian is evaluated. 

773 scale : None or float 

774 If scale is None, then the default scale will be calculated. 

775 Default scale is defined by `self.scaletype` and set in fit. 

776 If scale is not None, then it is used as a fixed scale. 

777 observed : bool 

778 If True, then the observed Hessian is returned. If false then the 

779 expected information matrix is returned. 

780 

781 Returns 

782 ------- 

783 ndarray 

784 A 1d weight vector used in the calculation of the Hessian. 

785 The hessian is obtained by `(exog.T * hessian_factor).dot(exog)`. 

786 """ 

787 

788 return self.weights 

789 

790 @Appender(_fit_regularized_doc) 

791 def fit_regularized(self, method="elastic_net", alpha=0., 

792 L1_wt=1., start_params=None, profile_scale=False, 

793 refit=False, **kwargs): 

794 # Docstring attached below 

795 

796 # Need to adjust since RSS/n in elastic net uses nominal n in 

797 # denominator 

798 alpha = alpha * np.sum(self.weights) / len(self.weights) 

799 

800 rslt = OLS(self.wendog, self.wexog).fit_regularized( 

801 method=method, alpha=alpha, 

802 L1_wt=L1_wt, 

803 start_params=start_params, 

804 profile_scale=profile_scale, 

805 refit=refit, **kwargs) 

806 

807 from statsmodels.base.elastic_net import ( 

808 RegularizedResults, RegularizedResultsWrapper) 

809 rrslt = RegularizedResults(self, rslt.params) 

810 return RegularizedResultsWrapper(rrslt) 

811 

812 

813class OLS(WLS): 

814 __doc__ = """ 

815 Ordinary Least Squares 

816 

817 %(params)s 

818 %(extra_params)s 

819 

820 Attributes 

821 ---------- 

822 weights : scalar 

823 Has an attribute weights = array(1.0) due to inheritance from WLS. 

824 

825 See Also 

826 -------- 

827 WLS : Fit a linear model using Weighted Least Squares. 

828 GLS : Fit a linear model using Generalized Least Squares. 

829 

830 Notes 

831 ----- 

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

833 

834 Examples 

835 -------- 

836 >>> import statsmodels.api as sm 

837 >>> Y = [1,3,4,5,2,3,4] 

838 >>> X = range(1,8) 

839 >>> X = sm.add_constant(X) 

840 >>> model = sm.OLS(Y,X) 

841 >>> results = model.fit() 

842 >>> results.params 

843 array([ 2.14285714, 0.25 ]) 

844 

845 >>> results.tvalues 

846 array([ 1.87867287, 0.98019606]) 

847 

848 >>> print(results.t_test([1, 0])) 

849 <T test: effect=array([ 2.14285714]), sd=array([[ 1.14062282]]), t=array([[ 1.87867287]]), p=array([[ 0.05953974]]), df_denom=5> 

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

851 <F test: F=array([[ 19.46078431]]), p=[[ 0.00437251]], df_denom=5, df_num=2> 

852 """ % {'params': base._model_params_doc, 

853 'extra_params': base._missing_param_doc + base._extra_param_doc} 

854 

855 # TODO: change example to use datasets. This was the point of datasets! 

856 def __init__(self, endog, exog=None, missing='none', hasconst=None, 

857 **kwargs): 

858 super(OLS, self).__init__(endog, exog, missing=missing, 

859 hasconst=hasconst, **kwargs) 

860 if "weights" in self._init_keys: 

861 self._init_keys.remove("weights") 

862 

863 def loglike(self, params, scale=None): 

864 """ 

865 The likelihood function for the OLS model. 

866 

867 Parameters 

868 ---------- 

869 params : array_like 

870 The coefficients with which to estimate the log-likelihood. 

871 scale : float or None 

872 If None, return the profile (concentrated) log likelihood 

873 (profiled over the scale parameter), else return the 

874 log-likelihood using the given scale value. 

875 

876 Returns 

877 ------- 

878 float 

879 The likelihood function evaluated at params. 

880 """ 

881 nobs2 = self.nobs / 2.0 

882 nobs = float(self.nobs) 

883 resid = self.endog - np.dot(self.exog, params) 

884 if hasattr(self, 'offset'): 

885 resid -= self.offset 

886 ssr = np.sum(resid**2) 

887 if scale is None: 

888 # profile log likelihood 

889 llf = -nobs2*np.log(2*np.pi) - nobs2*np.log(ssr / nobs) - nobs2 

890 else: 

891 # log-likelihood 

892 llf = -nobs2 * np.log(2 * np.pi * scale) - ssr / (2*scale) 

893 return llf 

894 

895 def whiten(self, x): 

896 """ 

897 OLS model whitener does nothing. 

898 

899 Parameters 

900 ---------- 

901 x : array_like 

902 Data to be whitened. 

903 

904 Returns 

905 ------- 

906 array_like 

907 The input array unmodified. 

908 

909 See Also 

910 -------- 

911 OLS : Fit a linear model using Ordinary Least Squares. 

912 """ 

913 return x 

914 

915 def score(self, params, scale=None): 

916 """ 

917 Evaluate the score function at a given point. 

918 

919 The score corresponds to the profile (concentrated) 

920 log-likelihood in which the scale parameter has been profiled 

921 out. 

922 

923 Parameters 

924 ---------- 

925 params : array_like 

926 The parameter vector at which the score function is 

927 computed. 

928 scale : float or None 

929 If None, return the profile (concentrated) log likelihood 

930 (profiled over the scale parameter), else return the 

931 log-likelihood using the given scale value. 

932 

933 Returns 

934 ------- 

935 ndarray 

936 The score vector. 

937 """ 

938 

939 if not hasattr(self, "_wexog_xprod"): 

940 self._setup_score_hess() 

941 

942 xtxb = np.dot(self._wexog_xprod, params) 

943 sdr = -self._wexog_x_wendog + xtxb 

944 

945 if scale is None: 

946 ssr = self._wendog_xprod - 2 * np.dot(self._wexog_x_wendog.T, 

947 params) 

948 ssr += np.dot(params, xtxb) 

949 return -self.nobs * sdr / ssr 

950 else: 

951 return -sdr / scale 

952 

953 def _setup_score_hess(self): 

954 y = self.wendog 

955 if hasattr(self, 'offset'): 

956 y = y - self.offset 

957 self._wendog_xprod = np.sum(y * y) 

958 self._wexog_xprod = np.dot(self.wexog.T, self.wexog) 

959 self._wexog_x_wendog = np.dot(self.wexog.T, y) 

960 

961 def hessian(self, params, scale=None): 

962 """ 

963 Evaluate the Hessian function at a given point. 

964 

965 Parameters 

966 ---------- 

967 params : array_like 

968 The parameter vector at which the Hessian is computed. 

969 scale : float or None 

970 If None, return the profile (concentrated) log likelihood 

971 (profiled over the scale parameter), else return the 

972 log-likelihood using the given scale value. 

973 

974 Returns 

975 ------- 

976 ndarray 

977 The Hessian matrix. 

978 """ 

979 

980 if not hasattr(self, "_wexog_xprod"): 

981 self._setup_score_hess() 

982 

983 xtxb = np.dot(self._wexog_xprod, params) 

984 

985 if scale is None: 

986 ssr = self._wendog_xprod - 2 * np.dot(self._wexog_x_wendog.T, 

987 params) 

988 ssr += np.dot(params, xtxb) 

989 ssrp = -2*self._wexog_x_wendog + 2*xtxb 

990 hm = self._wexog_xprod / ssr - np.outer(ssrp, ssrp) / ssr**2 

991 return -self.nobs * hm / 2 

992 else: 

993 return -self._wexog_xprod / scale 

994 

995 def hessian_factor(self, params, scale=None, observed=True): 

996 """ 

997 Calculate the weights for the Hessian. 

998 

999 Parameters 

1000 ---------- 

1001 params : ndarray 

1002 The parameter at which Hessian is evaluated. 

1003 scale : None or float 

1004 If scale is None, then the default scale will be calculated. 

1005 Default scale is defined by `self.scaletype` and set in fit. 

1006 If scale is not None, then it is used as a fixed scale. 

1007 observed : bool 

1008 If True, then the observed Hessian is returned. If false then the 

1009 expected information matrix is returned. 

1010 

1011 Returns 

1012 ------- 

1013 ndarray 

1014 A 1d weight vector used in the calculation of the Hessian. 

1015 The hessian is obtained by `(exog.T * hessian_factor).dot(exog)`. 

1016 """ 

1017 

1018 return np.ones(self.exog.shape[0]) 

1019 

1020 @Appender(_fit_regularized_doc) 

1021 def fit_regularized(self, method="elastic_net", alpha=0., 

1022 L1_wt=1., start_params=None, profile_scale=False, 

1023 refit=False, **kwargs): 

1024 

1025 # In the future we could add support for other penalties, e.g. SCAD. 

1026 if method not in ("elastic_net", "sqrt_lasso"): 

1027 msg = "Unknown method '%s' for fit_regularized" % method 

1028 raise ValueError(msg) 

1029 

1030 # Set default parameters. 

1031 defaults = {"maxiter": 50, "cnvrg_tol": 1e-10, 

1032 "zero_tol": 1e-8} 

1033 defaults.update(kwargs) 

1034 

1035 if method == "sqrt_lasso": 

1036 from statsmodels.base.elastic_net import ( 

1037 RegularizedResults, RegularizedResultsWrapper 

1038 ) 

1039 params = self._sqrt_lasso(alpha, refit, defaults["zero_tol"]) 

1040 results = RegularizedResults(self, params) 

1041 return RegularizedResultsWrapper(results) 

1042 

1043 from statsmodels.base.elastic_net import fit_elasticnet 

1044 

1045 if L1_wt == 0: 

1046 return self._fit_ridge(alpha) 

1047 

1048 # If a scale parameter is passed in, the non-profile 

1049 # likelihood (residual sum of squares divided by -2) is used, 

1050 # otherwise the profile likelihood is used. 

1051 if profile_scale: 

1052 loglike_kwds = {} 

1053 score_kwds = {} 

1054 hess_kwds = {} 

1055 else: 

1056 loglike_kwds = {"scale": 1} 

1057 score_kwds = {"scale": 1} 

1058 hess_kwds = {"scale": 1} 

1059 

1060 return fit_elasticnet(self, method=method, 

1061 alpha=alpha, 

1062 L1_wt=L1_wt, 

1063 start_params=start_params, 

1064 loglike_kwds=loglike_kwds, 

1065 score_kwds=score_kwds, 

1066 hess_kwds=hess_kwds, 

1067 refit=refit, 

1068 check_step=False, 

1069 **defaults) 

1070 

1071 def _sqrt_lasso(self, alpha, refit, zero_tol): 

1072 

1073 try: 

1074 import cvxopt 

1075 except ImportError: 

1076 msg = 'sqrt_lasso fitting requires the cvxopt module' 

1077 raise ValueError(msg) 

1078 

1079 n = len(self.endog) 

1080 p = self.exog.shape[1] 

1081 

1082 h0 = cvxopt.matrix(0., (2*p+1, 1)) 

1083 h1 = cvxopt.matrix(0., (n+1, 1)) 

1084 h1[1:, 0] = cvxopt.matrix(self.endog, (n, 1)) 

1085 

1086 G0 = cvxopt.spmatrix([], [], [], (2*p+1, 2*p+1)) 

1087 for i in range(1, 2*p+1): 

1088 G0[i, i] = -1 

1089 G1 = cvxopt.matrix(0., (n+1, 2*p+1)) 

1090 G1[0, 0] = -1 

1091 G1[1:, 1:p+1] = self.exog 

1092 G1[1:, p+1:] = -self.exog 

1093 

1094 c = cvxopt.matrix(alpha / n, (2*p + 1, 1)) 

1095 c[0] = 1 / np.sqrt(n) 

1096 

1097 from cvxopt import solvers 

1098 solvers.options["show_progress"] = False 

1099 

1100 rslt = solvers.socp(c, Gl=G0, hl=h0, Gq=[G1], hq=[h1]) 

1101 x = np.asarray(rslt['x']).flat 

1102 bp = x[1:p+1] 

1103 bn = x[p+1:] 

1104 params = bp - bn 

1105 

1106 if not refit: 

1107 return params 

1108 

1109 ii = np.flatnonzero(np.abs(params) > zero_tol) 

1110 rfr = OLS(self.endog, self.exog[:, ii]).fit() 

1111 params *= 0 

1112 params[ii] = rfr.params 

1113 

1114 return params 

1115 

1116 def _fit_ridge(self, alpha): 

1117 """ 

1118 Fit a linear model using ridge regression. 

1119 

1120 Parameters 

1121 ---------- 

1122 alpha : scalar or array_like 

1123 The penalty weight. If a scalar, the same penalty weight 

1124 applies to all variables in the model. If a vector, it 

1125 must have the same length as `params`, and contains a 

1126 penalty weight for each coefficient. 

1127 

1128 Notes 

1129 ----- 

1130 Equivalent to fit_regularized with L1_wt = 0 (but implemented 

1131 more efficiently). 

1132 """ 

1133 

1134 u, s, vt = np.linalg.svd(self.exog, 0) 

1135 v = vt.T 

1136 q = np.dot(u.T, self.endog) * s 

1137 s2 = s * s 

1138 if np.isscalar(alpha): 

1139 sd = s2 + alpha * self.nobs 

1140 params = q / sd 

1141 params = np.dot(v, params) 

1142 else: 

1143 vtav = self.nobs * np.dot(vt, alpha[:, None] * v) 

1144 d = np.diag(vtav) + s2 

1145 np.fill_diagonal(vtav, d) 

1146 r = np.linalg.solve(vtav, q) 

1147 params = np.dot(v, r) 

1148 

1149 from statsmodels.base.elastic_net import RegularizedResults 

1150 return RegularizedResults(self, params) 

1151 

1152 

1153class GLSAR(GLS): 

1154 __doc__ = """ 

1155 Generalized Least Squares with AR covariance structure 

1156 

1157 %(params)s 

1158 rho : int 

1159 The order of the autoregressive covariance. 

1160 %(extra_params)s 

1161 

1162 Notes 

1163 ----- 

1164 GLSAR is considered to be experimental. 

1165 The linear autoregressive process of order p--AR(p)--is defined as: 

1166 TODO 

1167 

1168 Examples 

1169 -------- 

1170 >>> import statsmodels.api as sm 

1171 >>> X = range(1,8) 

1172 >>> X = sm.add_constant(X) 

1173 >>> Y = [1,3,4,5,8,10,9] 

1174 >>> model = sm.GLSAR(Y, X, rho=2) 

1175 >>> for i in range(6): 

1176 ... results = model.fit() 

1177 ... print("AR coefficients: {0}".format(model.rho)) 

1178 ... rho, sigma = sm.regression.yule_walker(results.resid, 

1179 ... order=model.order) 

1180 ... model = sm.GLSAR(Y, X, rho) 

1181 ... 

1182 AR coefficients: [ 0. 0.] 

1183 AR coefficients: [-0.52571491 -0.84496178] 

1184 AR coefficients: [-0.6104153 -0.86656458] 

1185 AR coefficients: [-0.60439494 -0.857867 ] 

1186 AR coefficients: [-0.6048218 -0.85846157] 

1187 AR coefficients: [-0.60479146 -0.85841922] 

1188 >>> results.params 

1189 array([-0.66661205, 1.60850853]) 

1190 >>> results.tvalues 

1191 array([ -2.10304127, 21.8047269 ]) 

1192 >>> print(results.t_test([1, 0])) 

1193 <T test: effect=array([-0.66661205]), sd=array([[ 0.31697526]]), t=array([[-2.10304127]]), p=array([[ 0.06309969]]), df_denom=3> 

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

1195 <F test: F=array([[ 1815.23061844]]), p=[[ 0.00002372]], df_denom=3, df_num=2> 

1196 

1197 Or, equivalently 

1198 

1199 >>> model2 = sm.GLSAR(Y, X, rho=2) 

1200 >>> res = model2.iterative_fit(maxiter=6) 

1201 >>> model2.rho 

1202 array([-0.60479146, -0.85841922]) 

1203 """ % {'params': base._model_params_doc, 

1204 'extra_params': base._missing_param_doc + base._extra_param_doc} 

1205 # TODO: Complete docstring 

1206 

1207 def __init__(self, endog, exog=None, rho=1, missing='none', hasconst=None, 

1208 **kwargs): 

1209 # this looks strange, interpreting rho as order if it is int 

1210 if isinstance(rho, np.int): 

1211 self.order = rho 

1212 self.rho = np.zeros(self.order, np.float64) 

1213 else: 

1214 self.rho = np.squeeze(np.asarray(rho)) 

1215 if len(self.rho.shape) not in [0, 1]: 

1216 raise ValueError("AR parameters must be a scalar or a vector") 

1217 if self.rho.shape == (): 

1218 self.rho.shape = (1,) 

1219 self.order = self.rho.shape[0] 

1220 if exog is None: 

1221 # JP this looks wrong, should be a regression on constant 

1222 # results for rho estimate now identical to yule-walker on y 

1223 # super(AR, self).__init__(endog, add_constant(endog)) 

1224 super(GLSAR, self).__init__(endog, np.ones((endog.shape[0], 1)), 

1225 missing=missing, hasconst=None, 

1226 **kwargs) 

1227 else: 

1228 super(GLSAR, self).__init__(endog, exog, missing=missing, 

1229 **kwargs) 

1230 

1231 def iterative_fit(self, maxiter=3, rtol=1e-4, **kwargs): 

1232 """ 

1233 Perform an iterative two-stage procedure to estimate a GLS model. 

1234 

1235 The model is assumed to have AR(p) errors, AR(p) parameters and 

1236 regression coefficients are estimated iteratively. 

1237 

1238 Parameters 

1239 ---------- 

1240 maxiter : int, optional 

1241 The number of iterations. 

1242 rtol : float, optional 

1243 Relative tolerance between estimated coefficients to stop the 

1244 estimation. Stops if max(abs(last - current) / abs(last)) < rtol. 

1245 **kwargs 

1246 Additional keyword arguments passed to `fit`. 

1247 

1248 Returns 

1249 ------- 

1250 RegressionResults 

1251 The results computed using an iterative fit. 

1252 """ 

1253 # TODO: update this after going through example. 

1254 converged = False 

1255 i = -1 # need to initialize for maxiter < 1 (skip loop) 

1256 history = {'params': [], 'rho': [self.rho]} 

1257 for i in range(maxiter - 1): 

1258 if hasattr(self, 'pinv_wexog'): 

1259 del self.pinv_wexog 

1260 self.initialize() 

1261 results = self.fit() 

1262 history['params'].append(results.params) 

1263 if i == 0: 

1264 last = results.params 

1265 else: 

1266 diff = np.max(np.abs(last - results.params) / np.abs(last)) 

1267 if diff < rtol: 

1268 converged = True 

1269 break 

1270 last = results.params 

1271 self.rho, _ = yule_walker(results.resid, 

1272 order=self.order, df=None) 

1273 history['rho'].append(self.rho) 

1274 

1275 # why not another call to self.initialize 

1276 # Use kwarg to insert history 

1277 if not converged and maxiter > 0: 

1278 # maxiter <= 0 just does OLS 

1279 if hasattr(self, 'pinv_wexog'): 

1280 del self.pinv_wexog 

1281 self.initialize() 

1282 

1283 # if converged then this is a duplicate fit, because we did not 

1284 # update rho 

1285 results = self.fit(history=history, **kwargs) 

1286 results.iter = i + 1 

1287 # add last fit to history, not if duplicate fit 

1288 if not converged: 

1289 results.history['params'].append(results.params) 

1290 results.iter += 1 

1291 

1292 results.converged = converged 

1293 

1294 return results 

1295 

1296 def whiten(self, x): 

1297 """ 

1298 Whiten a series of columns according to an AR(p) covariance structure. 

1299 

1300 Whitening using this method drops the initial p observations. 

1301 

1302 Parameters 

1303 ---------- 

1304 x : array_like 

1305 The data to be whitened. 

1306 

1307 Returns 

1308 ------- 

1309 ndarray 

1310 The whitened data. 

1311 """ 

1312 # TODO: notation for AR process 

1313 x = np.asarray(x, np.float64) 

1314 _x = x.copy() 

1315 

1316 # the following loops over the first axis, works for 1d and nd 

1317 for i in range(self.order): 

1318 _x[(i + 1):] = _x[(i + 1):] - self.rho[i] * x[0:-(i + 1)] 

1319 return _x[self.order:] 

1320 

1321 

1322def yule_walker(x, order=1, method="unbiased", df=None, inv=False, 

1323 demean=True): 

1324 """ 

1325 Estimate AR(p) parameters from a sequence using the Yule-Walker equations. 

1326 

1327 Unbiased or maximum-likelihood estimator (mle) 

1328 

1329 Parameters 

1330 ---------- 

1331 x : array_like 

1332 A 1d array. 

1333 order : int, optional 

1334 The order of the autoregressive process. Default is 1. 

1335 method : str, optional 

1336 Method can be 'unbiased' or 'mle' and this determines 

1337 denominator in estimate of autocorrelation function (ACF) at 

1338 lag k. If 'mle', the denominator is n=X.shape[0], if 'unbiased' 

1339 the denominator is n-k. The default is unbiased. 

1340 df : int, optional 

1341 Specifies the degrees of freedom. If `df` is supplied, then it 

1342 is assumed the X has `df` degrees of freedom rather than `n`. 

1343 Default is None. 

1344 inv : bool 

1345 If inv is True the inverse of R is also returned. Default is 

1346 False. 

1347 demean : bool 

1348 True, the mean is subtracted from `X` before estimation. 

1349 

1350 Returns 

1351 ------- 

1352 rho : ndarray 

1353 AR(p) coefficients computed using the Yule-Walker method. 

1354 sigma : float 

1355 The estimate of the residual standard deviation. 

1356 

1357 See Also 

1358 -------- 

1359 burg : Burg's AR estimator. 

1360 

1361 Notes 

1362 ----- 

1363 See https://en.wikipedia.org/wiki/Autoregressive_moving_average_model for 

1364 further details. 

1365 

1366 Examples 

1367 -------- 

1368 >>> import statsmodels.api as sm 

1369 >>> from statsmodels.datasets.sunspots import load 

1370 >>> data = load(as_pandas=False) 

1371 >>> rho, sigma = sm.regression.yule_walker(data.endog, order=4, 

1372 ... method="mle") 

1373 

1374 >>> rho 

1375 array([ 1.28310031, -0.45240924, -0.20770299, 0.04794365]) 

1376 >>> sigma 

1377 16.808022730464351 

1378 """ 

1379 # TODO: define R better, look back at notes and technical notes on YW. 

1380 # First link here is useful 

1381 # http://www-stat.wharton.upenn.edu/~steele/Courses/956/ResourceDetails/YuleWalkerAndMore.htm 

1382 method = str(method).lower() 

1383 if method not in ["unbiased", "mle"]: 

1384 raise ValueError("ACF estimation method must be 'unbiased' or 'MLE'") 

1385 x = np.array(x, dtype=np.float64) 

1386 if demean: 

1387 x -= x.mean() 

1388 n = df or x.shape[0] 

1389 

1390 # this handles df_resid ie., n - p 

1391 adj_needed = method == "unbiased" 

1392 

1393 if x.ndim > 1 and x.shape[1] != 1: 

1394 raise ValueError("expecting a vector to estimate AR parameters") 

1395 r = np.zeros(order+1, np.float64) 

1396 r[0] = (x ** 2).sum() / n 

1397 for k in range(1, order+1): 

1398 r[k] = (x[0:-k] * x[k:]).sum() / (n - k * adj_needed) 

1399 R = toeplitz(r[:-1]) 

1400 

1401 rho = np.linalg.solve(R, r[1:]) 

1402 sigmasq = r[0] - (r[1:]*rho).sum() 

1403 if inv: 

1404 return rho, np.sqrt(sigmasq), np.linalg.inv(R) 

1405 else: 

1406 return rho, np.sqrt(sigmasq) 

1407 

1408 

1409def burg(endog, order=1, demean=True): 

1410 """ 

1411 Compute Burg's AP(p) parameter estimator. 

1412 

1413 Parameters 

1414 ---------- 

1415 endog : array_like 

1416 The endogenous variable. 

1417 order : int, optional 

1418 Order of the AR. Default is 1. 

1419 demean : bool, optional 

1420 Flag indicating to subtract the mean from endog before estimation. 

1421 

1422 Returns 

1423 ------- 

1424 rho : ndarray 

1425 The AR(p) coefficients computed using Burg's algorithm. 

1426 sigma2 : float 

1427 The estimate of the residual variance. 

1428 

1429 See Also 

1430 -------- 

1431 yule_walker : Estimate AR parameters using the Yule-Walker method. 

1432 

1433 Notes 

1434 ----- 

1435 AR model estimated includes a constant that is estimated using the sample 

1436 mean (see [1]_). This value is not reported. 

1437 

1438 References 

1439 ---------- 

1440 .. [1] Brockwell, P.J. and Davis, R.A., 2016. Introduction to time series 

1441 and forecasting. Springer. 

1442 

1443 Examples 

1444 -------- 

1445 >>> import statsmodels.api as sm 

1446 >>> from statsmodels.datasets.sunspots import load 

1447 >>> data = load(as_pandas=True) 

1448 >>> rho, sigma2 = sm.regression.linear_model.burg(data.endog, order=4) 

1449 

1450 >>> rho 

1451 array([ 1.30934186, -0.48086633, -0.20185982, 0.05501941]) 

1452 >>> sigma2 

1453 271.2467306963966 

1454 """ 

1455 # Avoid circular imports 

1456 from statsmodels.tsa.stattools import levinson_durbin_pacf, pacf_burg 

1457 

1458 endog = np.squeeze(np.asarray(endog)) 

1459 if endog.ndim != 1: 

1460 raise ValueError('endog must be 1-d or squeezable to 1-d.') 

1461 order = int(order) 

1462 if order < 1: 

1463 raise ValueError('order must be an integer larger than 1') 

1464 if demean: 

1465 endog = endog - endog.mean() 

1466 pacf, sigma = pacf_burg(endog, order, demean=demean) 

1467 ar, _ = levinson_durbin_pacf(pacf) 

1468 return ar, sigma[-1] 

1469 

1470 

1471class RegressionResults(base.LikelihoodModelResults): 

1472 r""" 

1473 This class summarizes the fit of a linear regression model. 

1474 

1475 It handles the output of contrasts, estimates of covariance, etc. 

1476 

1477 Parameters 

1478 ---------- 

1479 model : RegressionModel 

1480 The regression model instance. 

1481 params : ndarray 

1482 The estimated parameters. 

1483 normalized_cov_params : ndarray 

1484 The normalized covariance parameters. 

1485 scale : float 

1486 The estimated scale of the residuals. 

1487 cov_type : str 

1488 The covariance estimator used in the results. 

1489 cov_kwds : dict 

1490 Additional keywords used in the covariance specification. 

1491 use_t : bool 

1492 Flag indicating to use the Student's t in inference. 

1493 **kwargs 

1494 Additional keyword arguments used to initialize the results. 

1495 

1496 Attributes 

1497 ---------- 

1498 pinv_wexog 

1499 See model class docstring for implementation details. 

1500 cov_type 

1501 Parameter covariance estimator used for standard errors and t-stats. 

1502 df_model 

1503 Model degrees of freedom. The number of regressors `p`. Does not 

1504 include the constant if one is present. 

1505 df_resid 

1506 Residual degrees of freedom. `n - p - 1`, if a constant is present. 

1507 `n - p` if a constant is not included. 

1508 het_scale 

1509 adjusted squared residuals for heteroscedasticity robust standard 

1510 errors. Is only available after `HC#_se` or `cov_HC#` is called. 

1511 See HC#_se for more information. 

1512 history 

1513 Estimation history for iterative estimators. 

1514 model 

1515 A pointer to the model instance that called fit() or results. 

1516 params 

1517 The linear coefficients that minimize the least squares 

1518 criterion. This is usually called Beta for the classical 

1519 linear model. 

1520 """ 

1521 

1522 _cache = {} # needs to be a class attribute for scale setter? 

1523 

1524 def __init__(self, model, params, normalized_cov_params=None, scale=1., 

1525 cov_type='nonrobust', cov_kwds=None, use_t=None, **kwargs): 

1526 super(RegressionResults, self).__init__( 

1527 model, params, normalized_cov_params, scale) 

1528 

1529 self._cache = {} 

1530 if hasattr(model, 'wexog_singular_values'): 

1531 self._wexog_singular_values = model.wexog_singular_values 

1532 else: 

1533 self._wexog_singular_values = None 

1534 

1535 self.df_model = model.df_model 

1536 self.df_resid = model.df_resid 

1537 

1538 if cov_type == 'nonrobust': 

1539 self.cov_type = 'nonrobust' 

1540 self.cov_kwds = { 

1541 'description': 'Standard Errors assume that the ' + 

1542 'covariance matrix of the errors is correctly ' + 

1543 'specified.'} 

1544 if use_t is None: 

1545 use_t = True # TODO: class default 

1546 self.use_t = use_t 

1547 else: 

1548 if cov_kwds is None: 

1549 cov_kwds = {} 

1550 if 'use_t' in cov_kwds: 

1551 # TODO: we want to get rid of 'use_t' in cov_kwds 

1552 use_t_2 = cov_kwds.pop('use_t') 

1553 if use_t is None: 

1554 use_t = use_t_2 

1555 # TODO: warn or not? 

1556 self.get_robustcov_results(cov_type=cov_type, use_self=True, 

1557 use_t=use_t, **cov_kwds) 

1558 for key in kwargs: 

1559 setattr(self, key, kwargs[key]) 

1560 

1561 def conf_int(self, alpha=.05, cols=None): 

1562 """ 

1563 Compute the confidence interval of the fitted parameters. 

1564 

1565 Parameters 

1566 ---------- 

1567 alpha : float, optional 

1568 The `alpha` level for the confidence interval. The default 

1569 `alpha` = .05 returns a 95% confidence interval. 

1570 cols : array_like, optional 

1571 Columns to included in returned confidence intervals. 

1572 

1573 Returns 

1574 ------- 

1575 array_like 

1576 The confidence intervals. 

1577 

1578 Notes 

1579 ----- 

1580 The confidence interval is based on Student's t-distribution. 

1581 """ 

1582 # keep method for docstring for now 

1583 ci = super(RegressionResults, self).conf_int(alpha=alpha, cols=cols) 

1584 return ci 

1585 

1586 @cache_readonly 

1587 def nobs(self): 

1588 """Number of observations n.""" 

1589 return float(self.model.wexog.shape[0]) 

1590 

1591 @cache_readonly 

1592 def fittedvalues(self): 

1593 """The predicted values for the original (unwhitened) design.""" 

1594 return self.model.predict(self.params, self.model.exog) 

1595 

1596 @cache_readonly 

1597 def wresid(self): 

1598 """ 

1599 The residuals of the transformed/whitened regressand and regressor(s). 

1600 """ 

1601 return self.model.wendog - self.model.predict( 

1602 self.params, self.model.wexog) 

1603 

1604 @cache_readonly 

1605 def resid(self): 

1606 """The residuals of the model.""" 

1607 return self.model.endog - self.model.predict( 

1608 self.params, self.model.exog) 

1609 

1610 # TODO: fix writable example 

1611 @cache_writable() 

1612 def scale(self): 

1613 """ 

1614 A scale factor for the covariance matrix. 

1615 

1616 The Default value is ssr/(n-p). Note that the square root of `scale` 

1617 is often called the standard error of the regression. 

1618 """ 

1619 wresid = self.wresid 

1620 return np.dot(wresid, wresid) / self.df_resid 

1621 

1622 @cache_readonly 

1623 def ssr(self): 

1624 """Sum of squared (whitened) residuals.""" 

1625 wresid = self.wresid 

1626 return np.dot(wresid, wresid) 

1627 

1628 @cache_readonly 

1629 def centered_tss(self): 

1630 """The total (weighted) sum of squares centered about the mean.""" 

1631 model = self.model 

1632 weights = getattr(model, 'weights', None) 

1633 sigma = getattr(model, 'sigma', None) 

1634 if weights is not None: 

1635 mean = np.average(model.endog, weights=weights) 

1636 return np.sum(weights * (model.endog - mean)**2) 

1637 elif sigma is not None: 

1638 # Exactly matches WLS when sigma is diagonal 

1639 iota = np.ones_like(model.endog) 

1640 iota = model.whiten(iota) 

1641 mean = model.wendog.dot(iota) / iota.dot(iota) 

1642 err = model.endog - mean 

1643 err = model.whiten(err) 

1644 return np.sum(err**2) 

1645 else: 

1646 centered_endog = model.wendog - model.wendog.mean() 

1647 return np.dot(centered_endog, centered_endog) 

1648 

1649 @cache_readonly 

1650 def uncentered_tss(self): 

1651 """ 

1652 Uncentered sum of squares. 

1653 

1654 The sum of the squared values of the (whitened) endogenous response 

1655 variable. 

1656 """ 

1657 wendog = self.model.wendog 

1658 return np.dot(wendog, wendog) 

1659 

1660 @cache_readonly 

1661 def ess(self): 

1662 """ 

1663 The explained sum of squares. 

1664 

1665 If a constant is present, the centered total sum of squares minus the 

1666 sum of squared residuals. If there is no constant, the uncentered total 

1667 sum of squares is used. 

1668 """ 

1669 

1670 if self.k_constant: 

1671 return self.centered_tss - self.ssr 

1672 else: 

1673 return self.uncentered_tss - self.ssr 

1674 

1675 @cache_readonly 

1676 def rsquared(self): 

1677 """ 

1678 R-squared of the model. 

1679 

1680 This is defined here as 1 - `ssr`/`centered_tss` if the constant is 

1681 included in the model and 1 - `ssr`/`uncentered_tss` if the constant is 

1682 omitted. 

1683 """ 

1684 if self.k_constant: 

1685 return 1 - self.ssr/self.centered_tss 

1686 else: 

1687 return 1 - self.ssr/self.uncentered_tss 

1688 

1689 @cache_readonly 

1690 def rsquared_adj(self): 

1691 """ 

1692 Adjusted R-squared. 

1693 

1694 This is defined here as 1 - (`nobs`-1)/`df_resid` * (1-`rsquared`) 

1695 if a constant is included and 1 - `nobs`/`df_resid` * (1-`rsquared`) if 

1696 no constant is included. 

1697 """ 

1698 return 1 - (np.divide(self.nobs - self.k_constant, self.df_resid) 

1699 * (1 - self.rsquared)) 

1700 

1701 @cache_readonly 

1702 def mse_model(self): 

1703 """ 

1704 Mean squared error the model. 

1705 

1706 The explained sum of squares divided by the model degrees of freedom. 

1707 """ 

1708 if np.all(self.df_model == 0.0): 

1709 return np.full_like(self.ess, np.nan) 

1710 return self.ess/self.df_model 

1711 

1712 @cache_readonly 

1713 def mse_resid(self): 

1714 """ 

1715 Mean squared error of the residuals. 

1716 

1717 The sum of squared residuals divided by the residual degrees of 

1718 freedom. 

1719 """ 

1720 if np.all(self.df_resid == 0.0): 

1721 return np.full_like(self.ssr, np.nan) 

1722 return self.ssr/self.df_resid 

1723 

1724 @cache_readonly 

1725 def mse_total(self): 

1726 """ 

1727 Total mean squared error. 

1728 

1729 The uncentered total sum of squares divided by the number of 

1730 observations. 

1731 """ 

1732 if np.all(self.df_resid + self.df_model == 0.0): 

1733 return np.full_like(self.centered_tss, np.nan) 

1734 if self.k_constant: 

1735 return self.centered_tss / (self.df_resid + self.df_model) 

1736 else: 

1737 return self.uncentered_tss / (self.df_resid + self.df_model) 

1738 

1739 @cache_readonly 

1740 def fvalue(self): 

1741 """ 

1742 F-statistic of the fully specified model. 

1743 

1744 Calculated as the mean squared error of the model divided by the mean 

1745 squared error of the residuals if the nonrobust covariance is used. 

1746 Otherwise computed using a Wald-like quadratic form that tests whether 

1747 all coefficients (excluding the constant) are zero. 

1748 """ 

1749 if hasattr(self, 'cov_type') and self.cov_type != 'nonrobust': 

1750 # with heteroscedasticity or correlation robustness 

1751 k_params = self.normalized_cov_params.shape[0] 

1752 mat = np.eye(k_params) 

1753 const_idx = self.model.data.const_idx 

1754 # TODO: What if model includes implicit constant, e.g. all 

1755 # dummies but no constant regressor? 

1756 # TODO: Restats as LM test by projecting orthogonalizing 

1757 # to constant? 

1758 if self.model.data.k_constant == 1: 

1759 # if constant is implicit, return nan see #2444 

1760 if const_idx is None: 

1761 return np.nan 

1762 

1763 idx = lrange(k_params) 

1764 idx.pop(const_idx) 

1765 mat = mat[idx] # remove constant 

1766 if mat.size == 0: # see #3642 

1767 return np.nan 

1768 ft = self.f_test(mat) 

1769 # using backdoor to set another attribute that we already have 

1770 self._cache['f_pvalue'] = ft.pvalue 

1771 return ft.fvalue 

1772 else: 

1773 # for standard homoscedastic case 

1774 return self.mse_model/self.mse_resid 

1775 

1776 @cache_readonly 

1777 def f_pvalue(self): 

1778 """The p-value of the F-statistic.""" 

1779 # Special case for df_model 0 

1780 if self.df_model == 0: 

1781 return np.full_like(self.fvalue, np.nan) 

1782 return stats.f.sf(self.fvalue, self.df_model, self.df_resid) 

1783 

1784 @cache_readonly 

1785 def bse(self): 

1786 """The standard errors of the parameter estimates.""" 

1787 return np.sqrt(np.diag(self.cov_params())) 

1788 

1789 @cache_readonly 

1790 def aic(self): 

1791 r""" 

1792 Akaike's information criteria. 

1793 

1794 For a model with a constant :math:`-2llf + 2(df\_model + 1)`. For a 

1795 model without a constant :math:`-2llf + 2(df\_model)`. 

1796 """ 

1797 return -2 * self.llf + 2 * (self.df_model + self.k_constant) 

1798 

1799 @cache_readonly 

1800 def bic(self): 

1801 r""" 

1802 Bayes' information criteria. 

1803 

1804 For a model with a constant :math:`-2llf + \log(n)(df\_model+1)`. 

1805 For a model without a constant :math:`-2llf + \log(n)(df\_model)`. 

1806 """ 

1807 return (-2 * self.llf + np.log(self.nobs) * (self.df_model + 

1808 self.k_constant)) 

1809 

1810 @cache_readonly 

1811 def eigenvals(self): 

1812 """ 

1813 Return eigenvalues sorted in decreasing order. 

1814 """ 

1815 if self._wexog_singular_values is not None: 

1816 eigvals = self._wexog_singular_values ** 2 

1817 else: 

1818 eigvals = np.linalg.linalg.eigvalsh(np.dot(self.model.wexog.T, 

1819 self.model.wexog)) 

1820 return np.sort(eigvals)[::-1] 

1821 

1822 @cache_readonly 

1823 def condition_number(self): 

1824 """ 

1825 Return condition number of exogenous matrix. 

1826 

1827 Calculated as ratio of largest to smallest eigenvalue. 

1828 """ 

1829 eigvals = self.eigenvals 

1830 return np.sqrt(eigvals[0]/eigvals[-1]) 

1831 

1832 # TODO: make these properties reset bse 

1833 def _HCCM(self, scale): 

1834 H = np.dot(self.model.pinv_wexog, 

1835 scale[:, None] * self.model.pinv_wexog.T) 

1836 return H 

1837 

1838 @cache_readonly 

1839 def cov_HC0(self): 

1840 """ 

1841 Heteroscedasticity robust covariance matrix. See HC0_se. 

1842 """ 

1843 self.het_scale = self.wresid**2 

1844 cov_HC0 = self._HCCM(self.het_scale) 

1845 return cov_HC0 

1846 

1847 @cache_readonly 

1848 def cov_HC1(self): 

1849 """ 

1850 Heteroscedasticity robust covariance matrix. See HC1_se. 

1851 """ 

1852 self.het_scale = self.nobs/(self.df_resid)*(self.wresid**2) 

1853 cov_HC1 = self._HCCM(self.het_scale) 

1854 return cov_HC1 

1855 

1856 @cache_readonly 

1857 def cov_HC2(self): 

1858 """ 

1859 Heteroscedasticity robust covariance matrix. See HC2_se. 

1860 """ 

1861 # probably could be optimized 

1862 wexog = self.model.wexog 

1863 h = np.diag(wexog @ self.normalized_cov_params @ wexog.T) 

1864 self.het_scale = self.wresid**2/(1-h) 

1865 cov_HC2 = self._HCCM(self.het_scale) 

1866 return cov_HC2 

1867 

1868 @cache_readonly 

1869 def cov_HC3(self): 

1870 """ 

1871 Heteroscedasticity robust covariance matrix. See HC3_se. 

1872 """ 

1873 wexog = self.model.wexog 

1874 h = np.diag(wexog @ self.normalized_cov_params @ wexog.T) 

1875 self.het_scale = (self.wresid / (1 - h))**2 

1876 cov_HC3 = self._HCCM(self.het_scale) 

1877 return cov_HC3 

1878 

1879 @cache_readonly 

1880 def HC0_se(self): 

1881 """ 

1882 White's (1980) heteroskedasticity robust standard errors. 

1883 

1884 Notes 

1885 ----- 

1886 Defined as sqrt(diag(X.T X)^(-1)X.T diag(e_i^(2)) X(X.T X)^(-1) 

1887 where e_i = resid[i]. 

1888 

1889 When HC0_se or cov_HC0 is called the RegressionResults instance will 

1890 then have another attribute `het_scale`, which is in this case is just 

1891 resid**2. 

1892 """ 

1893 return np.sqrt(np.diag(self.cov_HC0)) 

1894 

1895 @cache_readonly 

1896 def HC1_se(self): 

1897 """ 

1898 MacKinnon and White's (1985) heteroskedasticity robust standard errors. 

1899 

1900 Notes 

1901 ----- 

1902 Defined as sqrt(diag(n/(n-p)*HC_0). 

1903 

1904 When HC1_se or cov_HC1 is called the RegressionResults instance will 

1905 then have another attribute `het_scale`, which is in this case is 

1906 n/(n-p)*resid**2. 

1907 """ 

1908 return np.sqrt(np.diag(self.cov_HC1)) 

1909 

1910 @cache_readonly 

1911 def HC2_se(self): 

1912 """ 

1913 MacKinnon and White's (1985) heteroskedasticity robust standard errors. 

1914 

1915 Notes 

1916 ----- 

1917 Defined as (X.T X)^(-1)X.T diag(e_i^(2)/(1-h_ii)) X(X.T X)^(-1) 

1918 where h_ii = x_i(X.T X)^(-1)x_i.T 

1919 

1920 When HC2_se or cov_HC2 is called the RegressionResults instance will 

1921 then have another attribute `het_scale`, which is in this case is 

1922 resid^(2)/(1-h_ii). 

1923 """ 

1924 return np.sqrt(np.diag(self.cov_HC2)) 

1925 

1926 @cache_readonly 

1927 def HC3_se(self): 

1928 """ 

1929 MacKinnon and White's (1985) heteroskedasticity robust standard errors. 

1930 

1931 Notes 

1932 ----- 

1933 Defined as (X.T X)^(-1)X.T diag(e_i^(2)/(1-h_ii)^(2)) X(X.T X)^(-1) 

1934 where h_ii = x_i(X.T X)^(-1)x_i.T. 

1935 

1936 When HC3_se or cov_HC3 is called the RegressionResults instance will 

1937 then have another attribute `het_scale`, which is in this case is 

1938 resid^(2)/(1-h_ii)^(2). 

1939 """ 

1940 return np.sqrt(np.diag(self.cov_HC3)) 

1941 

1942 @cache_readonly 

1943 def resid_pearson(self): 

1944 """ 

1945 Residuals, normalized to have unit variance. 

1946 

1947 Returns 

1948 ------- 

1949 array_like 

1950 The array `wresid` normalized by the sqrt of the scale to have 

1951 unit variance. 

1952 """ 

1953 

1954 if not hasattr(self, 'resid'): 

1955 raise ValueError('Method requires residuals.') 

1956 eps = np.finfo(self.wresid.dtype).eps 

1957 if np.sqrt(self.scale) < 10 * eps * self.model.endog.mean(): 

1958 # do not divide if scale is zero close to numerical precision 

1959 from warnings import warn 

1960 warn("All residuals are 0, cannot compute normed residuals.", 

1961 RuntimeWarning) 

1962 return self.wresid 

1963 else: 

1964 return self.wresid / np.sqrt(self.scale) 

1965 

1966 def _is_nested(self, restricted): 

1967 """ 

1968 Parameters 

1969 ---------- 

1970 restricted : Result instance 

1971 The restricted model is assumed to be nested in the current 

1972 model. The result instance of the restricted model is required to 

1973 have two attributes, residual sum of squares, `ssr`, residual 

1974 degrees of freedom, `df_resid`. 

1975 

1976 Returns 

1977 ------- 

1978 nested : bool 

1979 True if nested, otherwise false 

1980 

1981 Notes 

1982 ----- 

1983 A most nests another model if the regressors in the smaller 

1984 model are spanned by the regressors in the larger model and 

1985 the regressand is identical. 

1986 """ 

1987 

1988 if self.model.nobs != restricted.model.nobs: 

1989 return False 

1990 

1991 full_rank = self.model.rank 

1992 restricted_rank = restricted.model.rank 

1993 if full_rank <= restricted_rank: 

1994 return False 

1995 

1996 restricted_exog = restricted.model.wexog 

1997 full_wresid = self.wresid 

1998 

1999 scores = restricted_exog * full_wresid[:, None] 

2000 score_l2 = np.sqrt(np.mean(scores.mean(0) ** 2)) 

2001 # TODO: Could be improved, and may fail depending on scale of 

2002 # regressors 

2003 return np.allclose(score_l2, 0) 

2004 

2005 def compare_lm_test(self, restricted, demean=True, use_lr=False): 

2006 """ 

2007 Use Lagrange Multiplier test to test a set of linear restrictions. 

2008 

2009 Parameters 

2010 ---------- 

2011 restricted : Result instance 

2012 The restricted model is assumed to be nested in the 

2013 current model. The result instance of the restricted model 

2014 is required to have two attributes, residual sum of 

2015 squares, `ssr`, residual degrees of freedom, `df_resid`. 

2016 demean : bool 

2017 Flag indicating whether the demean the scores based on the 

2018 residuals from the restricted model. If True, the covariance of 

2019 the scores are used and the LM test is identical to the large 

2020 sample version of the LR test. 

2021 use_lr : bool 

2022 A flag indicating whether to estimate the covariance of the model 

2023 scores using the unrestricted model. Setting the to True improves 

2024 the power of the test. 

2025 

2026 Returns 

2027 ------- 

2028 lm_value : float 

2029 The test statistic which has a chi2 distributed. 

2030 p_value : float 

2031 The p-value of the test statistic. 

2032 df_diff : int 

2033 The degrees of freedom of the restriction, i.e. difference in df 

2034 between models. 

2035 

2036 Notes 

2037 ----- 

2038 The LM test examines whether the scores from the restricted model are 

2039 0. If the null is true, and the restrictions are valid, then the 

2040 parameters of the restricted model should be close to the minimum of 

2041 the sum of squared errors, and so the scores should be close to zero, 

2042 on average. 

2043 """ 

2044 import statsmodels.stats.sandwich_covariance as sw 

2045 from numpy.linalg import inv 

2046 

2047 if not self._is_nested(restricted): 

2048 raise ValueError("Restricted model is not nested by full model.") 

2049 

2050 wresid = restricted.wresid 

2051 wexog = self.model.wexog 

2052 scores = wexog * wresid[:, None] 

2053 

2054 n = self.nobs 

2055 df_full = self.df_resid 

2056 df_restr = restricted.df_resid 

2057 df_diff = (df_restr - df_full) 

2058 

2059 s = scores.mean(axis=0) 

2060 if use_lr: 

2061 scores = wexog * self.wresid[:, None] 

2062 demean = False 

2063 

2064 if demean: 

2065 scores = scores - scores.mean(0)[None, :] 

2066 # Form matters here. If homoskedastics can be sigma^2 (X'X)^-1 

2067 # If Heteroskedastic then the form below is fine 

2068 # If HAC then need to use HAC 

2069 # If Cluster, should use cluster 

2070 

2071 cov_type = getattr(self, 'cov_type', 'nonrobust') 

2072 if cov_type == 'nonrobust': 

2073 sigma2 = np.mean(wresid**2) 

2074 xpx = np.dot(wexog.T, wexog) / n 

2075 s_inv = inv(sigma2 * xpx) 

2076 elif cov_type in ('HC0', 'HC1', 'HC2', 'HC3'): 

2077 s_inv = inv(np.dot(scores.T, scores) / n) 

2078 elif cov_type == 'HAC': 

2079 maxlags = self.cov_kwds['maxlags'] 

2080 s_inv = inv(sw.S_hac_simple(scores, maxlags) / n) 

2081 elif cov_type == 'cluster': 

2082 # cluster robust standard errors 

2083 groups = self.cov_kwds['groups'] 

2084 # TODO: Might need demean option in S_crosssection by group? 

2085 s_inv = inv(sw.S_crosssection(scores, groups)) 

2086 else: 

2087 raise ValueError('Only nonrobust, HC, HAC and cluster are ' + 

2088 'currently connected') 

2089 

2090 lm_value = n * (s @ s_inv @ s.T) 

2091 p_value = stats.chi2.sf(lm_value, df_diff) 

2092 return lm_value, p_value, df_diff 

2093 

2094 def compare_f_test(self, restricted): 

2095 """ 

2096 Use F test to test whether restricted model is correct. 

2097 

2098 Parameters 

2099 ---------- 

2100 restricted : Result instance 

2101 The restricted model is assumed to be nested in the 

2102 current model. The result instance of the restricted model 

2103 is required to have two attributes, residual sum of 

2104 squares, `ssr`, residual degrees of freedom, `df_resid`. 

2105 

2106 Returns 

2107 ------- 

2108 f_value : float 

2109 The test statistic which has an F distribution. 

2110 p_value : float 

2111 The p-value of the test statistic. 

2112 df_diff : int 

2113 The degrees of freedom of the restriction, i.e. difference in 

2114 df between models. 

2115 

2116 Notes 

2117 ----- 

2118 See mailing list discussion October 17, 

2119 

2120 This test compares the residual sum of squares of the two 

2121 models. This is not a valid test, if there is unspecified 

2122 heteroscedasticity or correlation. This method will issue a 

2123 warning if this is detected but still return the results under 

2124 the assumption of homoscedasticity and no autocorrelation 

2125 (sphericity). 

2126 """ 

2127 

2128 has_robust1 = getattr(self, 'cov_type', 'nonrobust') != 'nonrobust' 

2129 has_robust2 = (getattr(restricted, 'cov_type', 'nonrobust') != 

2130 'nonrobust') 

2131 

2132 if has_robust1 or has_robust2: 

2133 warnings.warn('F test for comparison is likely invalid with ' + 

2134 'robust covariance, proceeding anyway', 

2135 InvalidTestWarning) 

2136 

2137 ssr_full = self.ssr 

2138 ssr_restr = restricted.ssr 

2139 df_full = self.df_resid 

2140 df_restr = restricted.df_resid 

2141 

2142 df_diff = (df_restr - df_full) 

2143 f_value = (ssr_restr - ssr_full) / df_diff / ssr_full * df_full 

2144 p_value = stats.f.sf(f_value, df_diff, df_full) 

2145 return f_value, p_value, df_diff 

2146 

2147 def compare_lr_test(self, restricted, large_sample=False): 

2148 """ 

2149 Likelihood ratio test to test whether restricted model is correct. 

2150 

2151 Parameters 

2152 ---------- 

2153 restricted : Result instance 

2154 The restricted model is assumed to be nested in the current model. 

2155 The result instance of the restricted model is required to have two 

2156 attributes, residual sum of squares, `ssr`, residual degrees of 

2157 freedom, `df_resid`. 

2158 

2159 large_sample : bool 

2160 Flag indicating whether to use a heteroskedasticity robust version 

2161 of the LR test, which is a modified LM test. 

2162 

2163 Returns 

2164 ------- 

2165 lr_stat : float 

2166 The likelihood ratio which is chisquare distributed with df_diff 

2167 degrees of freedom. 

2168 p_value : float 

2169 The p-value of the test statistic. 

2170 df_diff : int 

2171 The degrees of freedom of the restriction, i.e. difference in df 

2172 between models. 

2173 

2174 Notes 

2175 ----- 

2176 The exact likelihood ratio is valid for homoskedastic data, 

2177 and is defined as 

2178 

2179 .. math:: D=-2\\log\\left(\\frac{\\mathcal{L}_{null}} 

2180 {\\mathcal{L}_{alternative}}\\right) 

2181 

2182 where :math:`\\mathcal{L}` is the likelihood of the 

2183 model. With :math:`D` distributed as chisquare with df equal 

2184 to difference in number of parameters or equivalently 

2185 difference in residual degrees of freedom. 

2186 

2187 The large sample version of the likelihood ratio is defined as 

2188 

2189 .. math:: D=n s^{\\prime}S^{-1}s 

2190 

2191 where :math:`s=n^{-1}\\sum_{i=1}^{n} s_{i}` 

2192 

2193 .. math:: s_{i} = x_{i,alternative} \\epsilon_{i,null} 

2194 

2195 is the average score of the model evaluated using the 

2196 residuals from null model and the regressors from the 

2197 alternative model and :math:`S` is the covariance of the 

2198 scores, :math:`s_{i}`. The covariance of the scores is 

2199 estimated using the same estimator as in the alternative 

2200 model. 

2201 

2202 This test compares the loglikelihood of the two models. This 

2203 may not be a valid test, if there is unspecified 

2204 heteroscedasticity or correlation. This method will issue a 

2205 warning if this is detected but still return the results 

2206 without taking unspecified heteroscedasticity or correlation 

2207 into account. 

2208 

2209 This test compares the loglikelihood of the two models. This 

2210 may not be a valid test, if there is unspecified 

2211 heteroscedasticity or correlation. This method will issue a 

2212 warning if this is detected but still return the results 

2213 without taking unspecified heteroscedasticity or correlation 

2214 into account. 

2215 

2216 is the average score of the model evaluated using the 

2217 residuals from null model and the regressors from the 

2218 alternative model and :math:`S` is the covariance of the 

2219 scores, :math:`s_{i}`. The covariance of the scores is 

2220 estimated using the same estimator as in the alternative 

2221 model. 

2222 """ 

2223 # TODO: put into separate function, needs tests 

2224 

2225 # See mailing list discussion October 17, 

2226 

2227 if large_sample: 

2228 return self.compare_lm_test(restricted, use_lr=True) 

2229 

2230 has_robust1 = (getattr(self, 'cov_type', 'nonrobust') != 'nonrobust') 

2231 has_robust2 = ( 

2232 getattr(restricted, 'cov_type', 'nonrobust') != 'nonrobust') 

2233 

2234 if has_robust1 or has_robust2: 

2235 warnings.warn('Likelihood Ratio test is likely invalid with ' + 

2236 'robust covariance, proceeding anyway', 

2237 InvalidTestWarning) 

2238 

2239 llf_full = self.llf 

2240 llf_restr = restricted.llf 

2241 df_full = self.df_resid 

2242 df_restr = restricted.df_resid 

2243 

2244 lrdf = (df_restr - df_full) 

2245 lrstat = -2*(llf_restr - llf_full) 

2246 lr_pvalue = stats.chi2.sf(lrstat, lrdf) 

2247 

2248 return lrstat, lr_pvalue, lrdf 

2249 

2250 def get_robustcov_results(self, cov_type='HC1', use_t=None, **kwargs): 

2251 """ 

2252 Create new results instance with robust covariance as default. 

2253 

2254 Parameters 

2255 ---------- 

2256 cov_type : str 

2257 The type of robust sandwich estimator to use. See Notes below. 

2258 use_t : bool 

2259 If true, then the t distribution is used for inference. 

2260 If false, then the normal distribution is used. 

2261 If `use_t` is None, then an appropriate default is used, which is 

2262 `True` if the cov_type is nonrobust, and `False` in all other 

2263 cases. 

2264 **kwargs 

2265 Required or optional arguments for robust covariance calculation. 

2266 See Notes below. 

2267 

2268 Returns 

2269 ------- 

2270 RegressionResults 

2271 This method creates a new results instance with the 

2272 requested robust covariance as the default covariance of 

2273 the parameters. Inferential statistics like p-values and 

2274 hypothesis tests will be based on this covariance matrix. 

2275 

2276 Notes 

2277 ----- 

2278 The following covariance types and required or optional arguments are 

2279 currently available: 

2280 

2281 - 'fixed scale' and optional keyword argument 'scale' which uses 

2282 a predefined scale estimate with default equal to one. 

2283 - 'HC0', 'HC1', 'HC2', 'HC3' and no keyword arguments: 

2284 heteroscedasticity robust covariance 

2285 - 'HAC' and keywords 

2286 

2287 - `maxlag` integer (required) : number of lags to use 

2288 - `kernel` callable or str (optional) : kernel 

2289 currently available kernels are ['bartlett', 'uniform'], 

2290 default is Bartlett 

2291 - `use_correction` bool (optional) : If true, use small sample 

2292 correction 

2293 

2294 - 'cluster' and required keyword `groups`, integer group indicator 

2295 

2296 - `groups` array_like, integer (required) : 

2297 index of clusters or groups 

2298 - `use_correction` bool (optional) : 

2299 If True the sandwich covariance is calculated with a small 

2300 sample correction. 

2301 If False the sandwich covariance is calculated without 

2302 small sample correction. 

2303 - `df_correction` bool (optional) 

2304 If True (default), then the degrees of freedom for the 

2305 inferential statistics and hypothesis tests, such as 

2306 pvalues, f_pvalue, conf_int, and t_test and f_test, are 

2307 based on the number of groups minus one instead of the 

2308 total number of observations minus the number of explanatory 

2309 variables. `df_resid` of the results instance is adjusted. 

2310 If False, then `df_resid` of the results instance is not 

2311 adjusted. 

2312 

2313 - 'hac-groupsum' Driscoll and Kraay, heteroscedasticity and 

2314 autocorrelation robust standard errors in panel data 

2315 keywords 

2316 

2317 - `time` array_like (required) : index of time periods 

2318 - `maxlag` integer (required) : number of lags to use 

2319 - `kernel` callable or str (optional). The available kernels 

2320 are ['bartlett', 'uniform']. The default is Bartlett. 

2321 - `use_correction` False or string in ['hac', 'cluster'] (optional). 

2322 If False the the sandwich covariance is calculated without small 

2323 sample correction. If `use_correction = 'cluster'` (default), 

2324 then the same small sample correction as in the case of 

2325 `covtype='cluster'` is used. 

2326 - `df_correction` bool (optional) The adjustment to df_resid, see 

2327 cov_type 'cluster' above 

2328 # TODO: we need more options here 

2329 

2330 - 'hac-panel' heteroscedasticity and autocorrelation robust standard 

2331 errors in panel data. 

2332 The data needs to be sorted in this case, the time series 

2333 for each panel unit or cluster need to be stacked. The 

2334 membership to a timeseries of an individual or group can 

2335 be either specified by group indicators or by increasing 

2336 time periods. 

2337 

2338 keywords 

2339 

2340 - either `groups` or `time` : array_like (required) 

2341 `groups` : indicator for groups 

2342 `time` : index of time periods 

2343 - `maxlag` integer (required) : number of lags to use 

2344 - `kernel` callable or str (optional) 

2345 currently available kernels are ['bartlett', 'uniform'], 

2346 default is Bartlett 

2347 - `use_correction` False or string in ['hac', 'cluster'] (optional) 

2348 If False the sandwich covariance is calculated without 

2349 small sample correction. 

2350 - `df_correction` bool (optional) 

2351 adjustment to df_resid, see cov_type 'cluster' above 

2352 # TODO: we need more options here 

2353 

2354 Reminder: 

2355 `use_correction` in "hac-groupsum" and "hac-panel" is not bool, 

2356 needs to be in [False, 'hac', 'cluster'] 

2357 

2358 TODO: Currently there is no check for extra or misspelled keywords, 

2359 except in the case of cov_type `HCx` 

2360 """ 

2361 import statsmodels.stats.sandwich_covariance as sw 

2362 from statsmodels.base.covtype import normalize_cov_type, descriptions 

2363 

2364 cov_type = normalize_cov_type(cov_type) 

2365 

2366 if 'kernel' in kwargs: 

2367 kwargs['weights_func'] = kwargs.pop('kernel') 

2368 if 'weights_func' in kwargs and not callable(kwargs['weights_func']): 

2369 kwargs['weights_func'] = sw.kernel_dict[kwargs['weights_func']] 

2370 

2371 # TODO: make separate function that returns a robust cov plus info 

2372 use_self = kwargs.pop('use_self', False) 

2373 if use_self: 

2374 res = self 

2375 else: 

2376 res = self.__class__( 

2377 self.model, self.params, 

2378 normalized_cov_params=self.normalized_cov_params, 

2379 scale=self.scale) 

2380 

2381 res.cov_type = cov_type 

2382 # use_t might already be defined by the class, and already set 

2383 if use_t is None: 

2384 use_t = self.use_t 

2385 res.cov_kwds = {'use_t': use_t} # store for information 

2386 res.use_t = use_t 

2387 

2388 adjust_df = False 

2389 if cov_type in ['cluster', 'hac-panel', 'hac-groupsum']: 

2390 df_correction = kwargs.get('df_correction', None) 

2391 # TODO: check also use_correction, do I need all combinations? 

2392 if df_correction is not False: # i.e. in [None, True]: 

2393 # user did not explicitely set it to False 

2394 adjust_df = True 

2395 

2396 res.cov_kwds['adjust_df'] = adjust_df 

2397 

2398 # verify and set kwargs, and calculate cov 

2399 # TODO: this should be outsourced in a function so we can reuse it in 

2400 # other models 

2401 # TODO: make it DRYer repeated code for checking kwargs 

2402 if cov_type in ['fixed scale', 'fixed_scale']: 

2403 res.cov_kwds['description'] = descriptions['fixed_scale'] 

2404 

2405 res.cov_kwds['scale'] = scale = kwargs.get('scale', 1.) 

2406 res.cov_params_default = scale * res.normalized_cov_params 

2407 elif cov_type.upper() in ('HC0', 'HC1', 'HC2', 'HC3'): 

2408 if kwargs: 

2409 raise ValueError('heteroscedasticity robust covariance ' 

2410 'does not use keywords') 

2411 res.cov_kwds['description'] = descriptions[cov_type.upper()] 

2412 res.cov_params_default = getattr(self, 'cov_' + cov_type.upper()) 

2413 elif cov_type.lower() == 'hac': 

2414 # TODO: check if required, default in cov_hac_simple 

2415 maxlags = kwargs['maxlags'] 

2416 res.cov_kwds['maxlags'] = maxlags 

2417 weights_func = kwargs.get('weights_func', sw.weights_bartlett) 

2418 res.cov_kwds['weights_func'] = weights_func 

2419 use_correction = kwargs.get('use_correction', False) 

2420 res.cov_kwds['use_correction'] = use_correction 

2421 res.cov_kwds['description'] = descriptions['HAC'].format( 

2422 maxlags=maxlags, 

2423 correction=['without', 'with'][use_correction]) 

2424 

2425 res.cov_params_default = sw.cov_hac_simple( 

2426 self, nlags=maxlags, weights_func=weights_func, 

2427 use_correction=use_correction) 

2428 elif cov_type.lower() == 'cluster': 

2429 # cluster robust standard errors, one- or two-way 

2430 groups = kwargs['groups'] 

2431 if not hasattr(groups, 'shape'): 

2432 groups = np.asarray(groups).T 

2433 

2434 if groups.ndim >= 2: 

2435 groups = groups.squeeze() 

2436 

2437 res.cov_kwds['groups'] = groups 

2438 use_correction = kwargs.get('use_correction', True) 

2439 res.cov_kwds['use_correction'] = use_correction 

2440 if groups.ndim == 1: 

2441 if adjust_df: 

2442 # need to find number of groups 

2443 # duplicate work 

2444 self.n_groups = n_groups = len(np.unique(groups)) 

2445 res.cov_params_default = sw.cov_cluster( 

2446 self, groups, use_correction=use_correction) 

2447 

2448 elif groups.ndim == 2: 

2449 if hasattr(groups, 'values'): 

2450 groups = groups.values 

2451 

2452 if adjust_df: 

2453 # need to find number of groups 

2454 # duplicate work 

2455 n_groups0 = len(np.unique(groups[:, 0])) 

2456 n_groups1 = len(np.unique(groups[:, 1])) 

2457 self.n_groups = (n_groups0, n_groups1) 

2458 n_groups = min(n_groups0, n_groups1) # use for adjust_df 

2459 

2460 # Note: sw.cov_cluster_2groups has 3 returns 

2461 res.cov_params_default = sw.cov_cluster_2groups( 

2462 self, groups, use_correction=use_correction)[0] 

2463 else: 

2464 raise ValueError('only two groups are supported') 

2465 res.cov_kwds['description'] = descriptions['cluster'] 

2466 

2467 elif cov_type.lower() == 'hac-panel': 

2468 # cluster robust standard errors 

2469 res.cov_kwds['time'] = time = kwargs.get('time', None) 

2470 res.cov_kwds['groups'] = groups = kwargs.get('groups', None) 

2471 # TODO: nlags is currently required 

2472 # nlags = kwargs.get('nlags', True) 

2473 # res.cov_kwds['nlags'] = nlags 

2474 # TODO: `nlags` or `maxlags` 

2475 res.cov_kwds['maxlags'] = maxlags = kwargs['maxlags'] 

2476 use_correction = kwargs.get('use_correction', 'hac') 

2477 res.cov_kwds['use_correction'] = use_correction 

2478 weights_func = kwargs.get('weights_func', sw.weights_bartlett) 

2479 res.cov_kwds['weights_func'] = weights_func 

2480 if groups is not None: 

2481 groups = np.asarray(groups) 

2482 tt = (np.nonzero(groups[:-1] != groups[1:])[0] + 1).tolist() 

2483 nobs_ = len(groups) 

2484 elif time is not None: 

2485 time = np.asarray(time) 

2486 # TODO: clumsy time index in cov_nw_panel 

2487 tt = (np.nonzero(time[1:] < time[:-1])[0] + 1).tolist() 

2488 nobs_ = len(time) 

2489 else: 

2490 raise ValueError('either time or groups needs to be given') 

2491 groupidx = lzip([0] + tt, tt + [nobs_]) 

2492 self.n_groups = n_groups = len(groupidx) 

2493 res.cov_params_default = sw.cov_nw_panel(self, maxlags, groupidx, 

2494 weights_func=weights_func, 

2495 use_correction=use_correction) 

2496 res.cov_kwds['description'] = descriptions['HAC-Panel'] 

2497 

2498 elif cov_type.lower() == 'hac-groupsum': 

2499 # Driscoll-Kraay standard errors 

2500 res.cov_kwds['time'] = time = kwargs['time'] 

2501 # TODO: nlags is currently required 

2502 # nlags = kwargs.get('nlags', True) 

2503 # res.cov_kwds['nlags'] = nlags 

2504 # TODO: `nlags` or `maxlags` 

2505 res.cov_kwds['maxlags'] = maxlags = kwargs['maxlags'] 

2506 use_correction = kwargs.get('use_correction', 'cluster') 

2507 res.cov_kwds['use_correction'] = use_correction 

2508 weights_func = kwargs.get('weights_func', sw.weights_bartlett) 

2509 res.cov_kwds['weights_func'] = weights_func 

2510 if adjust_df: 

2511 # need to find number of groups 

2512 tt = (np.nonzero(time[1:] < time[:-1])[0] + 1) 

2513 self.n_groups = n_groups = len(tt) + 1 

2514 res.cov_params_default = sw.cov_nw_groupsum( 

2515 self, maxlags, time, weights_func=weights_func, 

2516 use_correction=use_correction) 

2517 res.cov_kwds['description'] = descriptions['HAC-Groupsum'] 

2518 else: 

2519 raise ValueError('cov_type not recognized. See docstring for ' + 

2520 'available options and spelling') 

2521 

2522 if adjust_df: 

2523 # Note: df_resid is used for scale and others, add new attribute 

2524 res.df_resid_inference = n_groups - 1 

2525 

2526 return res 

2527 

2528 @Appender(pred.get_prediction.__doc__) 

2529 def get_prediction(self, exog=None, transform=True, weights=None, 

2530 row_labels=None, **kwargs): 

2531 

2532 return pred.get_prediction( 

2533 self, exog=exog, transform=transform, weights=weights, 

2534 row_labels=row_labels, **kwargs) 

2535 

2536 def summary(self, yname=None, xname=None, title=None, alpha=.05): 

2537 """ 

2538 Summarize the Regression Results. 

2539 

2540 Parameters 

2541 ---------- 

2542 yname : str, optional 

2543 Name of endogenous (response) variable. The Default is `y`. 

2544 xname : list[str], optional 

2545 Names for the exogenous variables. Default is `var_##` for ## in 

2546 the number of regressors. Must match the number of parameters 

2547 in the model. 

2548 title : str, optional 

2549 Title for the top table. If not None, then this replaces the 

2550 default title. 

2551 alpha : float 

2552 The significance level for the confidence intervals. 

2553 

2554 Returns 

2555 ------- 

2556 Summary 

2557 Instance holding the summary tables and text, which can be printed 

2558 or converted to various output formats. 

2559 

2560 See Also 

2561 -------- 

2562 statsmodels.iolib.summary.Summary : A class that holds summary results. 

2563 """ 

2564 from statsmodels.stats.stattools import ( 

2565 jarque_bera, omni_normtest, durbin_watson) 

2566 

2567 jb, jbpv, skew, kurtosis = jarque_bera(self.wresid) 

2568 omni, omnipv = omni_normtest(self.wresid) 

2569 

2570 eigvals = self.eigenvals 

2571 condno = self.condition_number 

2572 

2573 # TODO: Avoid adding attributes in non-__init__ 

2574 self.diagn = dict(jb=jb, jbpv=jbpv, skew=skew, kurtosis=kurtosis, 

2575 omni=omni, omnipv=omnipv, condno=condno, 

2576 mineigval=eigvals[-1]) 

2577 

2578 # TODO not used yet 

2579 # diagn_left_header = ['Models stats'] 

2580 # diagn_right_header = ['Residual stats'] 

2581 

2582 # TODO: requiring list/iterable is a bit annoying 

2583 # need more control over formatting 

2584 # TODO: default do not work if it's not identically spelled 

2585 

2586 top_left = [('Dep. Variable:', None), 

2587 ('Model:', None), 

2588 ('Method:', ['Least Squares']), 

2589 ('Date:', None), 

2590 ('Time:', None), 

2591 ('No. Observations:', None), 

2592 ('Df Residuals:', None), 

2593 ('Df Model:', None), 

2594 ] 

2595 

2596 if hasattr(self, 'cov_type'): 

2597 top_left.append(('Covariance Type:', [self.cov_type])) 

2598 

2599 rsquared_type = '' if self.k_constant else ' (uncentered)' 

2600 top_right = [('R-squared' + rsquared_type + ':', 

2601 ["%#8.3f" % self.rsquared]), 

2602 ('Adj. R-squared' + rsquared_type + ':', 

2603 ["%#8.3f" % self.rsquared_adj]), 

2604 ('F-statistic:', ["%#8.4g" % self.fvalue]), 

2605 ('Prob (F-statistic):', ["%#6.3g" % self.f_pvalue]), 

2606 ('Log-Likelihood:', None), 

2607 ('AIC:', ["%#8.4g" % self.aic]), 

2608 ('BIC:', ["%#8.4g" % self.bic]) 

2609 ] 

2610 

2611 diagn_left = [('Omnibus:', ["%#6.3f" % omni]), 

2612 ('Prob(Omnibus):', ["%#6.3f" % omnipv]), 

2613 ('Skew:', ["%#6.3f" % skew]), 

2614 ('Kurtosis:', ["%#6.3f" % kurtosis]) 

2615 ] 

2616 

2617 diagn_right = [('Durbin-Watson:', 

2618 ["%#8.3f" % durbin_watson(self.wresid)] 

2619 ), 

2620 ('Jarque-Bera (JB):', ["%#8.3f" % jb]), 

2621 ('Prob(JB):', ["%#8.3g" % jbpv]), 

2622 ('Cond. No.', ["%#8.3g" % condno]) 

2623 ] 

2624 

2625 if title is None: 

2626 title = self.model.__class__.__name__ + ' ' + "Regression Results" 

2627 

2628 # create summary table instance 

2629 from statsmodels.iolib.summary import Summary 

2630 smry = Summary() 

2631 smry.add_table_2cols(self, gleft=top_left, gright=top_right, 

2632 yname=yname, xname=xname, title=title) 

2633 smry.add_table_params(self, yname=yname, xname=xname, alpha=alpha, 

2634 use_t=self.use_t) 

2635 

2636 smry.add_table_2cols(self, gleft=diagn_left, gright=diagn_right, 

2637 yname=yname, xname=xname, 

2638 title="") 

2639 

2640 # add warnings/notes, added to text format only 

2641 etext = [] 

2642 if hasattr(self, 'cov_type'): 

2643 etext.append(self.cov_kwds['description']) 

2644 if self.model.exog.shape[0] < self.model.exog.shape[1]: 

2645 wstr = "The input rank is higher than the number of observations." 

2646 etext.append(wstr) 

2647 if eigvals[-1] < 1e-10: 

2648 wstr = "The smallest eigenvalue is %6.3g. This might indicate " 

2649 wstr += "that there are\n" 

2650 wstr += "strong multicollinearity problems or that the design " 

2651 wstr += "matrix is singular." 

2652 wstr = wstr % eigvals[-1] 

2653 etext.append(wstr) 

2654 elif condno > 1000: # TODO: what is recommended? 

2655 wstr = "The condition number is large, %6.3g. This might " 

2656 wstr += "indicate that there are\n" 

2657 wstr += "strong multicollinearity or other numerical " 

2658 wstr += "problems." 

2659 wstr = wstr % condno 

2660 etext.append(wstr) 

2661 

2662 if etext: 

2663 etext = ["[{0}] {1}".format(i + 1, text) 

2664 for i, text in enumerate(etext)] 

2665 etext.insert(0, "Warnings:") 

2666 smry.add_extra_txt(etext) 

2667 

2668 return smry 

2669 

2670 def summary2(self, yname=None, xname=None, title=None, alpha=.05, 

2671 float_format="%.4f"): 

2672 """ 

2673 Experimental summary function to summarize the regression results. 

2674 

2675 Parameters 

2676 ---------- 

2677 yname : str 

2678 The name of the dependent variable (optional). 

2679 xname : list[str], optional 

2680 Names for the exogenous variables. Default is `var_##` for ## in 

2681 the number of regressors. Must match the number of parameters 

2682 in the model. 

2683 title : str, optional 

2684 Title for the top table. If not None, then this replaces the 

2685 default title. 

2686 alpha : float 

2687 The significance level for the confidence intervals. 

2688 float_format : str 

2689 The format for floats in parameters summary. 

2690 

2691 Returns 

2692 ------- 

2693 Summary 

2694 Instance holding the summary tables and text, which can be printed 

2695 or converted to various output formats. 

2696 

2697 See Also 

2698 -------- 

2699 statsmodels.iolib.summary2.Summary 

2700 A class that holds summary results. 

2701 """ 

2702 # Diagnostics 

2703 from statsmodels.stats.stattools import (jarque_bera, 

2704 omni_normtest, 

2705 durbin_watson) 

2706 

2707 from collections import OrderedDict 

2708 jb, jbpv, skew, kurtosis = jarque_bera(self.wresid) 

2709 omni, omnipv = omni_normtest(self.wresid) 

2710 dw = durbin_watson(self.wresid) 

2711 eigvals = self.eigenvals 

2712 condno = self.condition_number 

2713 eigvals = np.sort(eigvals) # in increasing order 

2714 diagnostic = OrderedDict([ 

2715 ('Omnibus:', "%.3f" % omni), 

2716 ('Prob(Omnibus):', "%.3f" % omnipv), 

2717 ('Skew:', "%.3f" % skew), 

2718 ('Kurtosis:', "%.3f" % kurtosis), 

2719 ('Durbin-Watson:', "%.3f" % dw), 

2720 ('Jarque-Bera (JB):', "%.3f" % jb), 

2721 ('Prob(JB):', "%.3f" % jbpv), 

2722 ('Condition No.:', "%.0f" % condno) 

2723 ]) 

2724 

2725 # Summary 

2726 from statsmodels.iolib import summary2 

2727 smry = summary2.Summary() 

2728 smry.add_base(results=self, alpha=alpha, float_format=float_format, 

2729 xname=xname, yname=yname, title=title) 

2730 smry.add_dict(diagnostic) 

2731 

2732 # Warnings 

2733 if eigvals[-1] < 1e-10: 

2734 warn = "The smallest eigenvalue is %6.3g. This might indicate that\ 

2735 there are strong multicollinearity problems or that the design\ 

2736 matrix is singular." % eigvals[-1] 

2737 smry.add_text(warn) 

2738 if condno > 1000: 

2739 warn = "* The condition number is large (%.g). This might indicate \ 

2740 strong multicollinearity or other numerical problems." % condno 

2741 smry.add_text(warn) 

2742 

2743 return smry 

2744 

2745 

2746class OLSResults(RegressionResults): 

2747 """ 

2748 Results class for for an OLS model. 

2749 

2750 Parameters 

2751 ---------- 

2752 model : RegressionModel 

2753 The regression model instance. 

2754 params : ndarray 

2755 The estimated parameters. 

2756 normalized_cov_params : ndarray 

2757 The normalized covariance parameters. 

2758 scale : float 

2759 The estimated scale of the residuals. 

2760 cov_type : str 

2761 The covariance estimator used in the results. 

2762 cov_kwds : dict 

2763 Additional keywords used in the covariance specification. 

2764 use_t : bool 

2765 Flag indicating to use the Student's t in inference. 

2766 **kwargs 

2767 Additional keyword arguments used to initialize the results. 

2768 

2769 See Also 

2770 -------- 

2771 RegressionResults 

2772 Results store for WLS and GLW models. 

2773 

2774 Notes 

2775 ----- 

2776 Most of the methods and attributes are inherited from RegressionResults. 

2777 The special methods that are only available for OLS are: 

2778 

2779 - get_influence 

2780 - outlier_test 

2781 - el_test 

2782 - conf_int_el 

2783 """ 

2784 

2785 def get_influence(self): 

2786 """ 

2787 Calculate influence and outlier measures. 

2788 

2789 Returns 

2790 ------- 

2791 OLSInfluence 

2792 The instance containing methods to calculate the main influence and 

2793 outlier measures for the OLS regression. 

2794 

2795 See Also 

2796 -------- 

2797 statsmodels.stats.outliers_influence.OLSInfluence 

2798 A class that exposes methods to examine observation influence. 

2799 """ 

2800 from statsmodels.stats.outliers_influence import OLSInfluence 

2801 return OLSInfluence(self) 

2802 

2803 def outlier_test(self, method='bonf', alpha=.05, labels=None, 

2804 order=False, cutoff=None): 

2805 """ 

2806 Test observations for outliers according to method. 

2807 

2808 Parameters 

2809 ---------- 

2810 method : str 

2811 The method to use in the outlier test. Must be one of: 

2812 

2813 - `bonferroni` : one-step correction 

2814 - `sidak` : one-step correction 

2815 - `holm-sidak` : 

2816 - `holm` : 

2817 - `simes-hochberg` : 

2818 - `hommel` : 

2819 - `fdr_bh` : Benjamini/Hochberg 

2820 - `fdr_by` : Benjamini/Yekutieli 

2821 

2822 See `statsmodels.stats.multitest.multipletests` for details. 

2823 alpha : float 

2824 The familywise error rate (FWER). 

2825 labels : None or array_like 

2826 If `labels` is not None, then it will be used as index to the 

2827 returned pandas DataFrame. See also Returns below. 

2828 order : bool 

2829 Whether or not to order the results by the absolute value of the 

2830 studentized residuals. If labels are provided they will also be 

2831 sorted. 

2832 cutoff : None or float in [0, 1] 

2833 If cutoff is not None, then the return only includes observations 

2834 with multiple testing corrected p-values strictly below the cutoff. 

2835 The returned array or dataframe can be empty if t. 

2836 

2837 Returns 

2838 ------- 

2839 array_like 

2840 Returns either an ndarray or a DataFrame if labels is not None. 

2841 Will attempt to get labels from model_results if available. The 

2842 columns are the Studentized residuals, the unadjusted p-value, 

2843 and the corrected p-value according to method. 

2844 

2845 Notes 

2846 ----- 

2847 The unadjusted p-value is stats.t.sf(abs(resid), df) where 

2848 df = df_resid - 1. 

2849 """ 

2850 from statsmodels.stats.outliers_influence import outlier_test 

2851 return outlier_test(self, method, alpha, labels=labels, 

2852 order=order, cutoff=cutoff) 

2853 

2854 def el_test(self, b0_vals, param_nums, return_weights=0, ret_params=0, 

2855 method='nm', stochastic_exog=1): 

2856 """ 

2857 Test single or joint hypotheses using Empirical Likelihood. 

2858 

2859 Parameters 

2860 ---------- 

2861 b0_vals : 1darray 

2862 The hypothesized value of the parameter to be tested. 

2863 param_nums : 1darray 

2864 The parameter number to be tested. 

2865 return_weights : bool 

2866 If true, returns the weights that optimize the likelihood 

2867 ratio at b0_vals. The default is False. 

2868 ret_params : bool 

2869 If true, returns the parameter vector that maximizes the likelihood 

2870 ratio at b0_vals. Also returns the weights. The default is False. 

2871 method : str 

2872 Can either be 'nm' for Nelder-Mead or 'powell' for Powell. The 

2873 optimization method that optimizes over nuisance parameters. 

2874 The default is 'nm'. 

2875 stochastic_exog : bool 

2876 When True, the exogenous variables are assumed to be stochastic. 

2877 When the regressors are nonstochastic, moment conditions are 

2878 placed on the exogenous variables. Confidence intervals for 

2879 stochastic regressors are at least as large as non-stochastic 

2880 regressors. The default is True. 

2881 

2882 Returns 

2883 ------- 

2884 tuple 

2885 The p-value and -2 times the log-likelihood ratio for the 

2886 hypothesized values. 

2887 

2888 Examples 

2889 -------- 

2890 >>> import statsmodels.api as sm 

2891 >>> data = sm.datasets.stackloss.load(as_pandas=False) 

2892 >>> endog = data.endog 

2893 >>> exog = sm.add_constant(data.exog) 

2894 >>> model = sm.OLS(endog, exog) 

2895 >>> fitted = model.fit() 

2896 >>> fitted.params 

2897 >>> array([-39.91967442, 0.7156402 , 1.29528612, -0.15212252]) 

2898 >>> fitted.rsquared 

2899 >>> 0.91357690446068196 

2900 >>> # Test that the slope on the first variable is 0 

2901 >>> fitted.el_test([0], [1]) 

2902 >>> (27.248146353888796, 1.7894660442330235e-07) 

2903 """ 

2904 params = np.copy(self.params) 

2905 opt_fun_inst = _ELRegOpts() # to store weights 

2906 if len(param_nums) == len(params): 

2907 llr = opt_fun_inst._opt_nuis_regress( 

2908 [], 

2909 param_nums=param_nums, 

2910 endog=self.model.endog, 

2911 exog=self.model.exog, 

2912 nobs=self.model.nobs, 

2913 nvar=self.model.exog.shape[1], 

2914 params=params, 

2915 b0_vals=b0_vals, 

2916 stochastic_exog=stochastic_exog) 

2917 pval = 1 - stats.chi2.cdf(llr, len(param_nums)) 

2918 if return_weights: 

2919 return llr, pval, opt_fun_inst.new_weights 

2920 else: 

2921 return llr, pval 

2922 x0 = np.delete(params, param_nums) 

2923 args = (param_nums, self.model.endog, self.model.exog, 

2924 self.model.nobs, self.model.exog.shape[1], params, 

2925 b0_vals, stochastic_exog) 

2926 if method == 'nm': 

2927 llr = optimize.fmin(opt_fun_inst._opt_nuis_regress, x0, 

2928 maxfun=10000, maxiter=10000, full_output=1, 

2929 disp=0, args=args)[1] 

2930 if method == 'powell': 

2931 llr = optimize.fmin_powell(opt_fun_inst._opt_nuis_regress, x0, 

2932 full_output=1, disp=0, 

2933 args=args)[1] 

2934 

2935 pval = 1 - stats.chi2.cdf(llr, len(param_nums)) 

2936 if ret_params: 

2937 return llr, pval, opt_fun_inst.new_weights, opt_fun_inst.new_params 

2938 elif return_weights: 

2939 return llr, pval, opt_fun_inst.new_weights 

2940 else: 

2941 return llr, pval 

2942 

2943 def conf_int_el(self, param_num, sig=.05, upper_bound=None, 

2944 lower_bound=None, method='nm', stochastic_exog=True): 

2945 """ 

2946 Compute the confidence interval using Empirical Likelihood. 

2947 

2948 Parameters 

2949 ---------- 

2950 param_num : float 

2951 The parameter for which the confidence interval is desired. 

2952 sig : float 

2953 The significance level. Default is 0.05. 

2954 upper_bound : float 

2955 The maximum value the upper limit can be. Default is the 

2956 99.9% confidence value under OLS assumptions. 

2957 lower_bound : float 

2958 The minimum value the lower limit can be. Default is the 99.9% 

2959 confidence value under OLS assumptions. 

2960 method : str 

2961 Can either be 'nm' for Nelder-Mead or 'powell' for Powell. The 

2962 optimization method that optimizes over nuisance parameters. 

2963 The default is 'nm'. 

2964 stochastic_exog : bool 

2965 When True, the exogenous variables are assumed to be stochastic. 

2966 When the regressors are nonstochastic, moment conditions are 

2967 placed on the exogenous variables. Confidence intervals for 

2968 stochastic regressors are at least as large as non-stochastic 

2969 regressors. The default is True. 

2970 

2971 Returns 

2972 ------- 

2973 lowerl : float 

2974 The lower bound of the confidence interval. 

2975 upperl : float 

2976 The upper bound of the confidence interval. 

2977 

2978 See Also 

2979 -------- 

2980 el_test : Test parameters using Empirical Likelihood. 

2981 

2982 Notes 

2983 ----- 

2984 This function uses brentq to find the value of beta where 

2985 test_beta([beta], param_num)[1] is equal to the critical value. 

2986 

2987 The function returns the results of each iteration of brentq at each 

2988 value of beta. 

2989 

2990 The current function value of the last printed optimization should be 

2991 the critical value at the desired significance level. For alpha=.05, 

2992 the value is 3.841459. 

2993 

2994 To ensure optimization terminated successfully, it is suggested to do 

2995 el_test([lower_limit], [param_num]). 

2996 

2997 If the optimization does not terminate successfully, consider switching 

2998 optimization algorithms. 

2999 

3000 If optimization is still not successful, try changing the values of 

3001 start_int_params. If the current function value repeatedly jumps 

3002 from a number between 0 and the critical value and a very large number 

3003 (>50), the starting parameters of the interior minimization need 

3004 to be changed. 

3005 """ 

3006 r0 = stats.chi2.ppf(1 - sig, 1) 

3007 if upper_bound is None: 

3008 upper_bound = self.conf_int(.01)[param_num][1] 

3009 if lower_bound is None: 

3010 lower_bound = self.conf_int(.01)[param_num][0] 

3011 

3012 def f(b0): 

3013 return self.el_test(np.array([b0]), np.array([param_num]), 

3014 method=method, 

3015 stochastic_exog=stochastic_exog)[0] - r0 

3016 

3017 lowerl = optimize.brenth(f, lower_bound, 

3018 self.params[param_num]) 

3019 upperl = optimize.brenth(f, self.params[param_num], 

3020 upper_bound) 

3021 # ^ Seems to be faster than brentq in most cases 

3022 return (lowerl, upperl) 

3023 

3024 

3025class RegressionResultsWrapper(wrap.ResultsWrapper): 

3026 

3027 _attrs = { 

3028 'chisq': 'columns', 

3029 'sresid': 'rows', 

3030 'weights': 'rows', 

3031 'wresid': 'rows', 

3032 'bcov_unscaled': 'cov', 

3033 'bcov_scaled': 'cov', 

3034 'HC0_se': 'columns', 

3035 'HC1_se': 'columns', 

3036 'HC2_se': 'columns', 

3037 'HC3_se': 'columns', 

3038 'norm_resid': 'rows', 

3039 } 

3040 

3041 _wrap_attrs = wrap.union_dicts(base.LikelihoodResultsWrapper._attrs, 

3042 _attrs) 

3043 

3044 _methods = {} 

3045 

3046 _wrap_methods = wrap.union_dicts( 

3047 base.LikelihoodResultsWrapper._wrap_methods, 

3048 _methods) 

3049 

3050 

3051wrap.populate_wrapper(RegressionResultsWrapper, 

3052 RegressionResults)