from collections import OrderedDict
import contextlib
import warnings
import numpy as np
import pandas as pd
from scipy.stats import norm
from statsmodels.base.data import PandasData
from statsmodels.tools.decorators import cache_readonly
from statsmodels.tools.eval_measures import aic, aicc, bic, hqic
from statsmodels.tools.sm_exceptions import PrecisionWarning
from statsmodels.tools.numdiff import (
_get_epsilon,
approx_fprime,
approx_fprime_cs,
approx_hess_cs,
)
from statsmodels.tools.tools import pinv_extended
import statsmodels.tsa.base.tsa_model as tsbase
from statsmodels.tsa.statespace.tools import _safe_cond
class StateSpaceMLEModel(tsbase.TimeSeriesModel):
"""
This is a temporary base model from ETS, here I just copy everything I need
from statespace.mlemodel.MLEModel
"""
def __init__(
self, endog, exog=None, dates=None, freq=None, missing="none", **kwargs
):
# TODO: this was changed from the original, requires some work when
# using this as base class for state space and exponential smoothing
super().__init__(
endog=endog, exog=exog, dates=dates, freq=freq, missing=missing
)
# Store kwargs to recreate model
self._init_kwargs = kwargs
# Prepared the endog array: C-ordered, shape=(nobs x k_endog)
self.endog, self.exog = self.prepare_data(self.data)
self.use_pandas = isinstance(self.data, PandasData)
# Dimensions
self.nobs = self.endog.shape[0]
# Setup holder for fixed parameters
self._has_fixed_params = False
self._fixed_params = None
self._params_index = None
self._fixed_params_index = None
self._free_params_index = None
@staticmethod
def prepare_data(data):
raise NotImplementedError
def clone(self, endog, exog=None, **kwargs):
raise NotImplementedError
def _validate_can_fix_params(self, param_names):
for param_name in param_names:
if param_name not in self.param_names:
raise ValueError(
'Invalid parameter name passed: "%s".' % param_name
)
@property
def k_params(self):
return len(self.param_names)
@contextlib.contextmanager
def fix_params(self, params):
"""
Fix parameters to specific values (context manager)
Parameters
----------
params : dict
Dictionary describing the fixed parameter values, of the form
`param_name: fixed_value`. See the `param_names` property for valid
parameter names.
Examples
--------
>>> mod = sm.tsa.SARIMAX(endog, order=(1, 0, 1))
>>> with mod.fix_params({'ar.L1': 0.5}):
res = mod.fit()
"""
# Initialization (this is done here rather than in the constructor
# because param_names may not be available at that point)
if self._fixed_params is None:
self._fixed_params = {}
self._params_index = OrderedDict(
zip(self.param_names, np.arange(self.k_params))
)
# Cache the current fixed parameters
cache_fixed_params = self._fixed_params.copy()
cache_has_fixed_params = self._has_fixed_params
cache_fixed_params_index = self._fixed_params_index
cache_free_params_index = self._free_params_index
# Validate parameter names and values
all_fixed_param_names = (
set(params.keys()) | set(self._fixed_params.keys())
)
self._validate_can_fix_params(all_fixed_param_names)
# Set the new fixed parameters, keeping the order as given by
# param_names
self._fixed_params.update(params)
self._fixed_params = OrderedDict(
[
(name, self._fixed_params[name])
for name in self.param_names
if name in self._fixed_params
]
)
# Update associated values
self._has_fixed_params = True
self._fixed_params_index = [
self._params_index[key] for key in self._fixed_params.keys()
]
self._free_params_index = list(
set(np.arange(self.k_params)).difference(self._fixed_params_index)
)
try:
yield
finally:
# Reset the fixed parameters
self._has_fixed_params = cache_has_fixed_params
self._fixed_params = cache_fixed_params
self._fixed_params_index = cache_fixed_params_index
self._free_params_index = cache_free_params_index
def fit_constrained(self, constraints, start_params=None, **fit_kwds):
"""
Fit the model with some parameters subject to equality constraints.
Parameters
----------
constraints : dict
Dictionary of constraints, of the form `param_name: fixed_value`.
See the `param_names` property for valid parameter names.
start_params : array_like, optional
Initial guess of the solution for the loglikelihood maximization.
If None, the default is given by Model.start_params.
**fit_kwds : keyword arguments
fit_kwds are used in the optimization of the remaining parameters.
Returns
-------
results : Results instance
Examples
--------
>>> mod = sm.tsa.SARIMAX(endog, order=(1, 0, 1))
>>> res = mod.fit_constrained({'ar.L1': 0.5})
"""
with self.fix_params(constraints):
res = self.fit(start_params, **fit_kwds)
return res
@property
def start_params(self):
"""
(array) Starting parameters for maximum likelihood estimation.
"""
if hasattr(self, "_start_params"):
return self._start_params
else:
raise NotImplementedError
@property
def param_names(self):
"""
(list of str) List of human readable parameter names (for parameters
actually included in the model).
"""
if hasattr(self, "_param_names"):
return self._param_names
else:
try:
names = ["param.%d" % i for i in range(len(self.start_params))]
except NotImplementedError:
names = []
return names
@classmethod
def from_formula(
cls, formula, data, subset=None, drop_cols=None, *args, **kwargs
):
"""
Not implemented for state space models
"""
raise NotImplementedError
def _wrap_data(self, data, start_idx, end_idx, names=None):
# TODO: check if this is reasonable for statespace
# squeezing data: data may be:
# - m x n: m dates, n simulations -> squeeze does nothing
# - m x 1: m dates, 1 simulation -> squeeze removes last dimension
# - 1 x n: don't squeeze, already fine
# - 1 x 1: squeeze only second axis
if data.ndim > 1 and data.shape[1] == 1:
data = np.squeeze(data, axis=1)
if self.use_pandas:
if data.shape[0]:
_, _, _, index = self._get_prediction_index(start_idx, end_idx)
else:
index = None
if data.ndim < 2:
data = pd.Series(data, index=index, name=names)
else:
data = pd.DataFrame(data, index=index, columns=names)
return data
def _wrap_results(
self,
params,
result,
return_raw,
cov_type=None,
cov_kwds=None,
results_class=None,
wrapper_class=None,
):
if not return_raw:
# Wrap in a results object
result_kwargs = {}
if cov_type is not None:
result_kwargs["cov_type"] = cov_type
if cov_kwds is not None:
result_kwargs["cov_kwds"] = cov_kwds
if results_class is None:
results_class = self._res_classes["fit"][0]
if wrapper_class is None:
wrapper_class = self._res_classes["fit"][1]
res = results_class(self, params, result, **result_kwargs)
result = wrapper_class(res)
return result
def _score_complex_step(self, params, **kwargs):
# the default epsilon can be too small
# inversion_method = INVERT_UNIVARIATE | SOLVE_LU
epsilon = _get_epsilon(params, 2., None, len(params))
kwargs['transformed'] = True
kwargs['complex_step'] = True
return approx_fprime_cs(params, self.loglike, epsilon=epsilon,
kwargs=kwargs)
def _score_finite_difference(self, params, approx_centered=False,
**kwargs):
kwargs['transformed'] = True
return approx_fprime(params, self.loglike, kwargs=kwargs,
centered=approx_centered)
def _hessian_finite_difference(self, params, approx_centered=False,
**kwargs):
params = np.array(params, ndmin=1)
warnings.warn('Calculation of the Hessian using finite differences'
' is usually subject to substantial approximation'
' errors.', PrecisionWarning)
if not approx_centered:
epsilon = _get_epsilon(params, 3, None, len(params))
else:
epsilon = _get_epsilon(params, 4, None, len(params)) / 2
hessian = approx_fprime(params, self._score_finite_difference,
epsilon=epsilon, kwargs=kwargs,
centered=approx_centered)
# TODO: changed this to nobs_effective, has to be changed when merging
# with statespace mlemodel
return hessian / (self.nobs_effective)
def _hessian_complex_step(self, params, **kwargs):
"""
Hessian matrix computed by second-order complex-step differentiation
on the `loglike` function.
"""
# the default epsilon can be too small
epsilon = _get_epsilon(params, 3., None, len(params))
kwargs['transformed'] = True
kwargs['complex_step'] = True
hessian = approx_hess_cs(
params, self.loglike, epsilon=epsilon, kwargs=kwargs)
# TODO: changed this to nobs_effective, has to be changed when merging
# with statespace mlemodel
return hessian / (self.nobs_effective)
class StateSpaceMLEResults(tsbase.TimeSeriesModelResults):
r"""
Class to hold results from fitting a state space model.
Parameters
----------
model : MLEModel instance
The fitted model instance
params : ndarray
Fitted parameters
Attributes
----------
model : Model instance
A reference to the model that was fit.
nobs : float
The number of observations used to fit the model.
params : ndarray
The parameters of the model.
"""
def __init__(self, model, params, scale=1.0):
self.data = model.data
self.endog = model.data.orig_endog
super().__init__(model, params, None, scale=scale)
# Save the fixed parameters
self._has_fixed_params = self.model._has_fixed_params
self._fixed_params_index = self.model._fixed_params_index
self._free_params_index = self.model._free_params_index
# TODO: seems like maybe self.fixed_params should be the dictionary
# itself, not just the keys?
if self._has_fixed_params:
self._fixed_params = self.model._fixed_params.copy()
self.fixed_params = list(self._fixed_params.keys())
else:
self._fixed_params = None
self.fixed_params = []
self.param_names = [
"%s (fixed)" % name if name in self.fixed_params else name
for name in (self.data.param_names or [])
]
# Dimensions
self.nobs = self.model.nobs
self.k_params = self.model.k_params
self._rank = None
@cache_readonly
def nobs_effective(self):
raise NotImplementedError
@cache_readonly
def df_resid(self):
return self.nobs_effective - self.df_model
@cache_readonly
def aic(self):
"""
(float) Akaike Information Criterion
"""
return aic(self.llf, self.nobs_effective, self.df_model)
@cache_readonly
def aicc(self):
"""
(float) Akaike Information Criterion with small sample correction
"""
return aicc(self.llf, self.nobs_effective, self.df_model)
@cache_readonly
def bic(self):
"""
(float) Bayes Information Criterion
"""
return bic(self.llf, self.nobs_effective, self.df_model)
@cache_readonly
def fittedvalues(self):
# TODO
raise NotImplementedError
@cache_readonly
def hqic(self):
"""
(float) Hannan-Quinn Information Criterion
"""
# return (-2 * self.llf +
# 2 * np.log(np.log(self.nobs_effective)) * self.df_model)
return hqic(self.llf, self.nobs_effective, self.df_model)
@cache_readonly
def llf(self):
"""
(float) The value of the log-likelihood function evaluated at `params`.
"""
raise NotImplementedError
@cache_readonly
def mae(self):
"""
(float) Mean absolute error
"""
return np.mean(np.abs(self.resid))
@cache_readonly
def mse(self):
"""
(float) Mean squared error
"""
return self.sse / self.nobs
@cache_readonly
def pvalues(self):
"""
(array) The p-values associated with the z-statistics of the
coefficients. Note that the coefficients are assumed to have a Normal
distribution.
"""
pvalues = np.zeros_like(self.zvalues) * np.nan
mask = np.ones_like(pvalues, dtype=bool)
mask[self._free_params_index] = True
mask &= ~np.isnan(self.zvalues)
pvalues[mask] = norm.sf(np.abs(self.zvalues[mask])) * 2
return pvalues
@cache_readonly
def resid(self):
raise NotImplementedError
@cache_readonly
def sse(self):
"""
(float) Sum of squared errors
"""
return np.sum(self.resid ** 2)
@cache_readonly
def zvalues(self):
"""
(array) The z-statistics for the coefficients.
"""
return self.params / self.bse
def _get_prediction_start_index(self, anchor):
"""Returns a valid numeric start index for predictions/simulations"""
if anchor is None or anchor == "start":
iloc = 0
elif anchor == "end":
iloc = self.nobs
else:
iloc, _, _ = self.model._get_index_loc(anchor)
if isinstance(iloc, slice):
iloc = iloc.start
iloc += 1 # anchor is one before start of prediction/simulation
if iloc < 0:
iloc = self.nobs + iloc
if iloc > self.nobs:
raise ValueError("Cannot anchor simulation outside of the sample.")
return iloc
def _cov_params_approx(
self, approx_complex_step=True, approx_centered=False
):
evaluated_hessian = self.nobs_effective * self.model.hessian(
params=self.params,
transformed=True,
includes_fixed=True,
method="approx",
approx_complex_step=approx_complex_step,
approx_centered=approx_centered,
)
# TODO: Case with "not approx_complex_step" is not hit in
# tests as of 2017-05-19
if len(self.fixed_params) > 0:
mask = np.ix_(self._free_params_index, self._free_params_index)
if len(self.fixed_params) < self.k_params:
(tmp, singular_values) = pinv_extended(evaluated_hessian[mask])
else:
tmp, singular_values = np.nan, [np.nan]
neg_cov = np.zeros_like(evaluated_hessian) * np.nan
neg_cov[mask] = tmp
else:
(neg_cov, singular_values) = pinv_extended(evaluated_hessian)
self.model.update(self.params, transformed=True, includes_fixed=True)
if self._rank is None:
self._rank = np.linalg.matrix_rank(np.diag(singular_values))
return -neg_cov
@cache_readonly
def cov_params_approx(self):
"""
(array) The variance / covariance matrix. Computed using the numerical
Hessian approximated by complex step or finite differences methods.
"""
return self._cov_params_approx(
self._cov_approx_complex_step, self._cov_approx_centered
)
def test_serial_correlation(self, method, lags=None):
"""
Ljung-Box test for no serial correlation of standardized residuals
Null hypothesis is no serial correlation.
Parameters
----------
method : {'ljungbox','boxpierece', None}
The statistical test for serial correlation. If None, an attempt is
made to select an appropriate test.
lags : None, int or array_like
If lags is an integer then this is taken to be the largest lag
that is included, the test result is reported for all smaller lag
length.
If lags is a list or array, then all lags are included up to the
largest lag in the list, however only the tests for the lags in the
list are reported.
If lags is None, then the default maxlag is min(10, nobs//5) for
non-seasonal time series and min (2*m, nobs//5) for seasonal time
series.
Returns
-------
output : ndarray
An array with `(test_statistic, pvalue)` for each endogenous
variable and each lag. The array is then sized
`(k_endog, 2, lags)`. If the method is called as
`ljungbox = res.test_serial_correlation()`, then `ljungbox[i]`
holds the results of the Ljung-Box test (as would be returned by
`statsmodels.stats.diagnostic.acorr_ljungbox`) for the `i` th
endogenous variable.
See Also
--------
statsmodels.stats.diagnostic.acorr_ljungbox
Ljung-Box test for serial correlation.
Notes
-----
For statespace models: let `d` = max(loglikelihood_burn, nobs_diffuse);
this test is calculated ignoring the first `d` residuals.
Output is nan for any endogenous variable which has missing values.
"""
if method is None:
method = 'ljungbox'
if self.standardized_forecasts_error is None:
raise ValueError('Cannot compute test statistic when standardized'
' forecast errors have not been computed.')
if method == 'ljungbox' or method == 'boxpierce':
from statsmodels.stats.diagnostic import acorr_ljungbox
if hasattr(self, "loglikelihood_burn"):
d = np.maximum(self.loglikelihood_burn, self.nobs_diffuse)
# This differs from self.nobs_effective because here we want to
# exclude exact diffuse periods, whereas self.nobs_effective
# only excludes explicitly burned (usually approximate diffuse)
# periods.
nobs_effective = self.nobs - d
else:
nobs_effective = self.nobs_effective
output = []
# Default lags for acorr_ljungbox is 40, but may not always have
# that many observations
if lags is None:
seasonal_periods = getattr(self.model, "seasonal_periods", 0)
if seasonal_periods:
lags = min(2 * seasonal_periods, nobs_effective // 5)
else:
lags = min(10, nobs_effective // 5)
for i in range(self.model.k_endog):
if hasattr(self, "filter_results"):
x = self.filter_results.standardized_forecasts_error[i][d:]
else:
x = self.standardized_forecasts_error
results = acorr_ljungbox(
x, lags=lags, boxpierce=(method == 'boxpierce'),
return_df=False)
if method == 'ljungbox':
output.append(results[0:2])
else:
output.append(results[2:])
output = np.c_[output]
else:
raise NotImplementedError('Invalid serial correlation test'
' method.')
return output
def test_heteroskedasticity(self, method, alternative='two-sided',
use_f=True):
r"""
Test for heteroskedasticity of standardized residuals
Tests whether the sum-of-squares in the first third of the sample is
significantly different than the sum-of-squares in the last third
of the sample. Analogous to a Goldfeld-Quandt test. The null hypothesis
is of no heteroskedasticity.
Parameters
----------
method : {'breakvar', None}
The statistical test for heteroskedasticity. Must be 'breakvar'
for test of a break in the variance. If None, an attempt is
made to select an appropriate test.
alternative : str, 'increasing', 'decreasing' or 'two-sided'
This specifies the alternative for the p-value calculation. Default
is two-sided.
use_f : bool, optional
Whether or not to compare against the asymptotic distribution
(chi-squared) or the approximate small-sample distribution (F).
Default is True (i.e. default is to compare against an F
distribution).
Returns
-------
output : ndarray
An array with `(test_statistic, pvalue)` for each endogenous
variable. The array is then sized `(k_endog, 2)`. If the method is
called as `het = res.test_heteroskedasticity()`, then `het[0]` is
an array of size 2 corresponding to the first endogenous variable,
where `het[0][0]` is the test statistic, and `het[0][1]` is the
p-value.
Notes
-----
The null hypothesis is of no heteroskedasticity. That means different
things depending on which alternative is selected:
- Increasing: Null hypothesis is that the variance is not increasing
throughout the sample; that the sum-of-squares in the later
subsample is *not* greater than the sum-of-squares in the earlier
subsample.
- Decreasing: Null hypothesis is that the variance is not decreasing
throughout the sample; that the sum-of-squares in the earlier
subsample is *not* greater than the sum-of-squares in the later
subsample.
- Two-sided: Null hypothesis is that the variance is not changing
throughout the sample. Both that the sum-of-squares in the earlier
subsample is not greater than the sum-of-squares in the later
subsample *and* that the sum-of-squares in the later subsample is
not greater than the sum-of-squares in the earlier subsample.
For :math:`h = [T/3]`, the test statistic is:
.. math::
H(h) = \sum_{t=T-h+1}^T \tilde v_t^2
\Bigg / \sum_{t=d+1}^{d+1+h} \tilde v_t^2
where :math:`d` = max(loglikelihood_burn, nobs_diffuse)` (usually
corresponding to diffuse initialization under either the approximate
or exact approach).
This statistic can be tested against an :math:`F(h,h)` distribution.
Alternatively, :math:`h H(h)` is asymptotically distributed according
to :math:`\chi_h^2`; this second test can be applied by passing
`asymptotic=True` as an argument.
See section 5.4 of [1]_ for the above formula and discussion, as well
as additional details.
TODO
- Allow specification of :math:`h`
References
----------
.. [1] Harvey, Andrew C. 1990. *Forecasting, Structural Time Series*
*Models and the Kalman Filter.* Cambridge University Press.
"""
if method is None:
method = 'breakvar'
if self.standardized_forecasts_error is None:
raise ValueError('Cannot compute test statistic when standardized'
' forecast errors have not been computed.')
if method == 'breakvar':
# Store some values
if hasattr(self, "filter_results"):
squared_resid = (
self.filter_results.standardized_forecasts_error**2
)
d = np.maximum(self.loglikelihood_burn, self.nobs_diffuse)
# This differs from self.nobs_effective because here we want to
# exclude exact diffuse periods, whereas self.nobs_effective
# only excludes explicitly burned (usually approximate diffuse)
# periods.
nobs_effective = self.nobs - d
else:
squared_resid = self.standardized_forecasts_error**2
if squared_resid.ndim == 1:
squared_resid = np.asarray(squared_resid)
squared_resid = squared_resid[np.newaxis, :]
nobs_effective = self.nobs_effective
d = 0
squared_resid = np.asarray(squared_resid)
test_statistics = []
p_values = []
for i in range(self.model.k_endog):
h = int(np.round(nobs_effective / 3))
numer_resid = squared_resid[i, -h:]
numer_resid = numer_resid[~np.isnan(numer_resid)]
numer_dof = len(numer_resid)
denom_resid = squared_resid[i, d:d + h]
denom_resid = denom_resid[~np.isnan(denom_resid)]
denom_dof = len(denom_resid)
if numer_dof < 2:
warnings.warn('Early subset of data for variable %d'
' has too few non-missing observations to'
' calculate test statistic.' % i)
numer_resid = np.nan
if denom_dof < 2:
warnings.warn('Later subset of data for variable %d'
' has too few non-missing observations to'
' calculate test statistic.' % i)
denom_resid = np.nan
test_statistic = np.sum(numer_resid) / np.sum(denom_resid)
# Setup functions to calculate the p-values
if use_f:
from scipy.stats import f
pval_lower = lambda test_statistics: f.cdf( # noqa:E731
test_statistics, numer_dof, denom_dof)
pval_upper = lambda test_statistics: f.sf( # noqa:E731
test_statistics, numer_dof, denom_dof)
else:
from scipy.stats import chi2
pval_lower = lambda test_statistics: chi2.cdf( # noqa:E731
numer_dof * test_statistics, denom_dof)
pval_upper = lambda test_statistics: chi2.sf( # noqa:E731
numer_dof * test_statistics, denom_dof)
# Calculate the one- or two-sided p-values
alternative = alternative.lower()
if alternative in ['i', 'inc', 'increasing']:
p_value = pval_upper(test_statistic)
elif alternative in ['d', 'dec', 'decreasing']:
test_statistic = 1. / test_statistic
p_value = pval_upper(test_statistic)
elif alternative in ['2', '2-sided', 'two-sided']:
p_value = 2 * np.minimum(
pval_lower(test_statistic),
pval_upper(test_statistic)
)
else:
raise ValueError('Invalid alternative.')
test_statistics.append(test_statistic)
p_values.append(p_value)
output = np.c_[test_statistics, p_values]
else:
raise NotImplementedError('Invalid heteroskedasticity test'
' method.')
return output
def test_normality(self, method):
"""
Test for normality of standardized residuals.
Null hypothesis is normality.
Parameters
----------
method : {'jarquebera', None}
The statistical test for normality. Must be 'jarquebera' for
Jarque-Bera normality test. If None, an attempt is made to select
an appropriate test.
See Also
--------
statsmodels.stats.stattools.jarque_bera
The Jarque-Bera test of normality.
Notes
-----
For statespace models: let `d` = max(loglikelihood_burn, nobs_diffuse);
this test is calculated ignoring the first `d` residuals.
In the case of missing data, the maintained hypothesis is that the
data are missing completely at random. This test is then run on the
standardized residuals excluding those corresponding to missing
observations.
"""
if method is None:
method = 'jarquebera'
if self.standardized_forecasts_error is None:
raise ValueError('Cannot compute test statistic when standardized'
' forecast errors have not been computed.')
if method == 'jarquebera':
from statsmodels.stats.stattools import jarque_bera
if hasattr(self, "loglikelihood_burn"):
d = np.maximum(self.loglikelihood_burn, self.nobs_diffuse)
else:
d = 0
output = []
for i in range(self.model.k_endog):
if hasattr(self, "fiter_results"):
resid = self.filter_results.standardized_forecasts_error[
i, d:
]
else:
resid = self.standardized_forecasts_error
mask = ~np.isnan(resid)
output.append(jarque_bera(resid[mask]))
else:
raise NotImplementedError('Invalid normality test method.')
return np.array(output)
def summary(
self,
alpha=0.05,
start=None,
title=None,
model_name=None,
display_params=True,
):
"""
Summarize the Model
Parameters
----------
alpha : float, optional
Significance level for the confidence intervals. Default is 0.05.
start : int, optional
Integer of the start observation. Default is 0.
model_name : str
The name of the model used. Default is to use model class name.
Returns
-------
summary : Summary instance
This holds the summary table and text, which can be printed or
converted to various output formats.
See Also
--------
statsmodels.iolib.summary.Summary
"""
from statsmodels.iolib.summary import Summary
# Model specification results
model = self.model
if title is None:
title = "Statespace Model Results"
if start is None:
start = 0
if self.model._index_dates:
ix = self.model._index
d = ix[start]
sample = ["%02d-%02d-%02d" % (d.month, d.day, d.year)]
d = ix[-1]
sample += ["- " + "%02d-%02d-%02d" % (d.month, d.day, d.year)]
else:
sample = [str(start), " - " + str(self.nobs)]
# Standardize the model name as a list of str
if model_name is None:
model_name = model.__class__.__name__
# Diagnostic tests results
try:
het = self.test_heteroskedasticity(method="breakvar")
except Exception: # FIXME: catch something specific
het = np.array([[np.nan] * 2])
try:
lb = self.test_serial_correlation(method="ljungbox")
except Exception: # FIXME: catch something specific
lb = np.array([[np.nan] * 2]).reshape(1, 2, 1)
try:
jb = self.test_normality(method="jarquebera")
except Exception: # FIXME: catch something specific
jb = np.array([[np.nan] * 4])
# Create the tables
if not isinstance(model_name, list):
model_name = [model_name]
top_left = [("Dep. Variable:", None)]
top_left.append(("Model:", [model_name[0]]))
for i in range(1, len(model_name)):
top_left.append(("", ["+ " + model_name[i]]))
top_left += [
("Date:", None),
("Time:", None),
("Sample:", [sample[0]]),
("", [sample[1]]),
]
top_right = [
("No. Observations:", [self.nobs]),
("Log Likelihood", ["%#5.3f" % self.llf]),
]
if hasattr(self, "rsquared"):
top_right.append(("R-squared:", ["%#8.3f" % self.rsquared]))
top_right += [
("AIC", ["%#5.3f" % self.aic]),
("BIC", ["%#5.3f" % self.bic]),
("HQIC", ["%#5.3f" % self.hqic]),
]
if hasattr(self, "filter_results"):
if (
self.filter_results is not None
and self.filter_results.filter_concentrated
):
top_right.append(("Scale", ["%#5.3f" % self.scale]))
else:
top_right.append(("Scale", ["%#5.3f" % self.scale]))
if hasattr(self, "cov_type"):
top_left.append(("Covariance Type:", [self.cov_type]))
format_str = lambda array: [ # noqa:E731
", ".join(["{0:.2f}".format(i) for i in array])
]
diagn_left = [
("Ljung-Box (Q):", format_str(lb[:, 0, -1])),
("Prob(Q):", format_str(lb[:, 1, -1])),
("Heteroskedasticity (H):", format_str(het[:, 0])),
("Prob(H) (two-sided):", format_str(het[:, 1])),
]
diagn_right = [
("Jarque-Bera (JB):", format_str(jb[:, 0])),
("Prob(JB):", format_str(jb[:, 1])),
("Skew:", format_str(jb[:, 2])),
("Kurtosis:", format_str(jb[:, 3])),
]
summary = Summary()
summary.add_table_2cols(
self, gleft=top_left, gright=top_right, title=title
)
if len(self.params) > 0 and display_params:
summary.add_table_params(
self, alpha=alpha, xname=self.param_names, use_t=False
)
summary.add_table_2cols(
self, gleft=diagn_left, gright=diagn_right, title=""
)
# Add warnings/notes, added to text format only
etext = []
if hasattr(self, "cov_type") and "description" in self.cov_kwds:
etext.append(self.cov_kwds["description"])
if self._rank < (len(self.params) - len(self.fixed_params)):
cov_params = self.cov_params()
if len(self.fixed_params) > 0:
mask = np.ix_(self._free_params_index, self._free_params_index)
cov_params = cov_params[mask]
etext.append(
"Covariance matrix is singular or near-singular,"
" with condition number %6.3g. Standard errors may be"
" unstable." % _safe_cond(cov_params)
)
if etext:
etext = [
"[{0}] {1}".format(i + 1, text) for i, text in enumerate(etext)
]
etext.insert(0, "Warnings:")
summary.add_extra_txt(etext)
return summary