Repository: statsmodels/statsmodels
Pull request: BUG/REF: fixing a few VAR problems
WIP 
trying to get a few fixes into 0.9

plot now return fig, rename test_whiteness methods closes #4537
(old test_whiteness is now plot_acorr)

I'm still thinking about #4535

Problem statement:
BUG: VAR test_whiteness is a plot function but uses sample_acorr
The results. plot produces by old VAR.test_whiteness doesn't make sense.
In 7 years nobody complained. !?

i.e. I'm estimating a VAR with three lags, but lag 1 and lag 2 are significant in the autocorrelation plot.

` acorrs = self.sample_acorr(nlags)`

Hints from issue discussion:
Actually, the doc string says sample (y) autocorrelation.
So, it's just that the method name is misleading and indicates something different (whiteness of endog while I would have expected residual whiteness as diagnostic check on the model).
Also, there is a separate plot_sample_acorr method

Commit message: BUG/DOC fix exog in VARProcess.simulate_var, docstrings

Patch to review:
diff --git a/statsmodels/tsa/vector_ar/irf.py b/statsmodels/tsa/vector_ar/irf.py
index f9e0e1c5293..267b2243f19 100644
--- a/statsmodels/tsa/vector_ar/irf.py
+++ b/statsmodels/tsa/vector_ar/irf.py
@@ -163,10 +163,12 @@ def plot(self, orth=False, impulse=None, response=None,
                                            seed=seed,
                                            component=component)
 
-        plotting.irf_grid_plot(irfs, stderr, impulse, response,
-                               self.model.names, title, signif=signif,
-                               subplot_params=subplot_params,
-                               plot_params=plot_params, stderr_type=stderr_type)
+        fig = plotting.irf_grid_plot(irfs, stderr, impulse, response,
+                                     self.model.names, title, signif=signif,
+                                     subplot_params=subplot_params,
+                                     plot_params=plot_params,
+                                     stderr_type=stderr_type)
+        return fig
 
     def plot_cum_effects(self, orth=False, impulse=None, response=None,
                          signif=0.05, plot_params=None,
@@ -222,10 +224,14 @@ def plot_cum_effects(self, orth=False, impulse=None, response=None,
         if not plot_stderr:
             stderr = None
 
-        plotting.irf_grid_plot(cum_effects, stderr, impulse, response,
-                               self.model.names, title, signif=signif,
-                               hlines=lr_effects, subplot_params=subplot_params,
-                               plot_params=plot_params, stderr_type=stderr_type)
+        fig = plotting.irf_grid_plot(cum_effects, stderr, impulse, response,
+                                   self.model.names, title, signif=signif,
+                                   hlines=lr_effects,
+                                   subplot_params=subplot_params,
+                                   plot_params=plot_params,
+                                   stderr_type=stderr_type)
+        return fig
+
 
 class IRAnalysis(BaseIRAnalysis):
     """
diff --git a/statsmodels/tsa/vector_ar/plotting.py b/statsmodels/tsa/vector_ar/plotting.py
index b75b4b5256a..422d010571a 100644
--- a/statsmodels/tsa/vector_ar/plotting.py
+++ b/statsmodels/tsa/vector_ar/plotting.py
@@ -33,12 +33,12 @@ def plot_mts(Y, names=None, index=None):
     k = Y.shape[1]
     rows, cols = k, 1
 
-    plt.figure(figsize=(10, 10))
+    fig = plt.figure(figsize=(10, 10))
 
     for j in range(k):
         ts = Y[:, j]
 
-        ax = plt.subplot(rows, cols, j+1)
+        ax = fig.add_subplot(rows, cols, j+1)
         if index is not None:
             ax.plot(index, ts)
         else:
@@ -47,6 +47,8 @@ def plot_mts(Y, names=None, index=None):
         if names is not None:
             ax.set_title(names[j])
 
+    return fig
+
 def plot_var_forc(prior, forc, err_upper, err_lower,
                   index=None, names=None, plot_stderr=True,
                   legend_options=None):
@@ -234,6 +236,8 @@ def irf_grid_plot(values, stderr, impcol, rescol, names, title,
         sz = subplot_params.get('fontsize', 12)
         ax.set_title(subtitle_temp % (names[j], names[i]), fontsize=sz)
 
+    return fig
+
 
 def _get_irf_plot_config(names, impcol, rescol):
     nrows = ncols = k = len(names)
diff --git a/statsmodels/tsa/vector_ar/svar_model.py b/statsmodels/tsa/vector_ar/svar_model.py
index 868e6806c56..4ac5de8f271 100644
--- a/statsmodels/tsa/vector_ar/svar_model.py
+++ b/statsmodels/tsa/vector_ar/svar_model.py
@@ -586,6 +586,7 @@ def __init__(self, endog, endog_lagged, params, sigma_u, lag_order,
         else:
             trendorder = None
         self.k_trend = k_trend
+        self.k_exog = k_trend  # now required by VARProcess
         self.trendorder = trendorder
 
         self.exog_names = util.make_lag_names(names, lag_order, k_trend)
diff --git a/statsmodels/tsa/vector_ar/tests/test_var.py b/statsmodels/tsa/vector_ar/tests/test_var.py
index 6b877e01441..4a87326a50d 100644
--- a/statsmodels/tsa/vector_ar/tests/test_var.py
+++ b/statsmodels/tsa/vector_ar/tests/test_var.py
@@ -160,6 +160,8 @@ def teardown_module():
 try:
     import matplotlib
     have_matplotlib = True
+    import matplotlib.pyplot as plt
+    plt.switch_backend('Agg')
 except ImportError:
     pass
 
@@ -639,6 +641,83 @@ def test_irf_trend():
     assert_allclose(irf_t.stderr()[1:4], irf.stderr()[1:4], rtol=0.03)
 
 
-if __name__ == '__main__':
-    import pytest
-    pytest.main([__file__, '-vvs', '-x', '--pdb'])
+class TestVARExtras(object):
+
+    @classmethod
+    def setup_class(cls):
+        mdata = sm.datasets.macrodata.load_pandas().data
+        mdata = mdata[['realgdp','realcons','realinv']]
+        data = mdata.values
+        data = np.diff(np.log(data), axis=0) * 400
+        cls.res0 = sm.tsa.VAR(data).fit(maxlags=2)
+
+    def test_process(self):
+        res0 = self.res0
+        k_ar = res0.k_ar
+        fc20 = res0.forecast(res0.endog[-k_ar:], 20)
+        mean_lr = res0.mean()
+        assert_allclose(mean_lr, fc20[-1], rtol=5e-4)
+
+        ysim = res0.simulate_var(seed=987128)
+        assert_allclose(ysim.mean(0), mean_lr, rtol=0.1)
+        # initialization does not use long run intercept, see #4542
+        assert_allclose(ysim[0], res0.intercept, rtol=1e-10)
+        assert_allclose(ysim[1], res0.intercept, rtol=1e-10)
+
+        n_sim = 900
+        ysimz = res0.simulate_var(steps=n_sim, offset=np.zeros((n_sim, 3)),
+                                  seed=987128)
+        zero3 = np.zeros(3)
+        assert_allclose(ysimz.mean(0), zero3, atol=0.4)
+        # initialization does not use long run intercept, see #4542
+        assert_allclose(ysimz[0], zero3, atol=1e-10)
+        assert_allclose(ysimz[1], zero3, atol=1e-10)
+
+        # check attributes
+        assert_equal(res0.k_trend, 1)
+        assert_equal(res0.k_exog_user, 0)
+        assert_equal(res0.k_exog, 1)
+        assert_equal(res0.k_ar, 2)
+
+        irf = res0.irf()
+
+        # SMOKE test
+        if have_matplotlib:
+            import matplotlib.pyplot as plt
+            fig = res0.plotsim()
+            plt.close(fig)
+            fig = res0.plot_acorr()
+            plt.close(fig)
+            fig = res0.plot_forecast(10)
+            plt.close(fig)
+            plt.close('all')
+            fig_asym = irf.plot()
+            fig_mc = irf.plot(stderr_type='mc')
+
+            for k in range(3):
+                a = fig_asym.axes[1].get_children()[k].get_ydata()
+                m = fig_mc.axes[1].get_children()[k].get_ydata()
+                assert_allclose(m, a, atol=0.03, rtol=0.01)
+            plt.close(fig_asym)
+            plt.close(fig_mc)
+            plt.close('all')
+
+    def test_exog(self):
+        # check that trend and exog are equivalent for basics and varsim
+        data = self.res0.model.endog
+        res_lin_trend = sm.tsa.VAR(data).fit(maxlags=2, trend='ct')
+        ex = np.arange(len(data))
+        res_lin_trend1 = sm.tsa.VAR(data, exog=ex).fit(maxlags=2)
+        ex2 = np.arange(len(data))[:, None]**[0, 1]
+        res_lin_trend2 = sm.tsa.VAR(data, exog=ex2).fit(maxlags=2, trend='nc')
+        # TODO: intercept differs by 4e-3, others are < 1e-12
+        assert_allclose(res_lin_trend.params, res_lin_trend1.params, rtol=5e-3)
+        assert_allclose(res_lin_trend.params, res_lin_trend2.params, rtol=5e-3)
+        assert_allclose(res_lin_trend1.params, res_lin_trend2.params, rtol=1e-10)
+
+        y1 = res_lin_trend.simulate_var(seed=987128)
+        y2 = res_lin_trend1.simulate_var(seed=987128)
+        y3 = res_lin_trend2.simulate_var(seed=987128)
+        assert_allclose(y2.mean(0), y1.mean(0), rtol=1e-12)
+        assert_allclose(y3.mean(0), y1.mean(0), rtol=1e-12)
+        assert_allclose(y3.mean(0), y2.mean(0), rtol=1e-12)
diff --git a/statsmodels/tsa/vector_ar/tests/test_var_jmulti.py b/statsmodels/tsa/vector_ar/tests/test_var_jmulti.py
index 07b794775c6..e16b02070eb 100644
--- a/statsmodels/tsa/vector_ar/tests/test_var_jmulti.py
+++ b/statsmodels/tsa/vector_ar/tests/test_var_jmulti.py
@@ -499,7 +499,7 @@ def test_whiteness():
                 print("\n" + dt_s_tup_to_string(dt) + ": ", end="")
             lags = results_ref[ds][dt]["whiteness"]["tested order"]
 
-            obtained = results_sm[ds][dt].test_whiteness_new(nlags=lags)
+            obtained = results_sm[ds][dt].test_whiteness(nlags=lags)
             # test statistic
             err_msg = build_err_msg(ds, dt, "WHITENESS OF RESIDUALS - "
                                             "TEST STATISTIC")
@@ -511,7 +511,7 @@ def test_whiteness():
             desired = results_ref[ds][dt]["whiteness"]["p-value"]
             assert_allclose(obtained.pvalue, desired, rtol, atol, False, err_msg)
 
-            obtained = results_sm[ds][dt].test_whiteness_new(nlags=lags,
+            obtained = results_sm[ds][dt].test_whiteness(nlags=lags,
                                                              adjusted=True)
             # test statistic (adjusted Portmanteau test)
             err_msg = build_err_msg(ds, dt, "WHITENESS OF RESIDUALS - "
diff --git a/statsmodels/tsa/vector_ar/tests/test_vecm.py b/statsmodels/tsa/vector_ar/tests/test_vecm.py
index bc27ef313c9..b2950c88a81 100644
--- a/statsmodels/tsa/vector_ar/tests/test_vecm.py
+++ b/statsmodels/tsa/vector_ar/tests/test_vecm.py
@@ -714,6 +714,7 @@ def test_var_to_vecm():
             sigma_u = results_sm[ds][dt].sigma_u
             coefs = results_sm[ds][dt].var_rep
             intercept = np.zeros(len(sigma_u))
+            # Note: _params_info k_trend, k_exog, ... is inferred with defaults
             var = VARProcess(coefs, intercept, sigma_u)
             vecm_results = var.to_vecm()
             obtained_pi = vecm_results["Pi"]
diff --git a/statsmodels/tsa/vector_ar/util.py b/statsmodels/tsa/vector_ar/util.py
index daf38299de5..5eb42305751 100644
--- a/statsmodels/tsa/vector_ar/util.py
+++ b/statsmodels/tsa/vector_ar/util.py
@@ -191,16 +191,54 @@ def acf_to_acorr(acf):
 
 def varsim(coefs, intercept, sig_u, steps=100, initvalues=None, seed=None):
     """
-    Simulate simple VAR(p) process with known coefficients, intercept, white
-    noise covariance, etc.
+    Simulate VAR(p) process, given coefficients and assuming Gaussian noise
+
+    Parameters
+    ----------
+    coefs : ndarray
+        Coefficients for the VAR lags of endog.
+    intercept : None or ndarray 1-D (neqs,) or (steps, neqs)
+        This can be either the intercept for each equation or an offset.
+        If None, then the VAR process has a zero intercept.
+        If intercept is 1-D, then the same (endog specific) intercept is added
+        to all observations.
+        If intercept is 2-D, then it is treated as an offset and is added as
+        an observation specific intercept to the autoregression. In this case,
+        the intercept/offset should have same number of rows as steps, and the
+        same number of columns as endogenous variables (neqs).
+    sig_u : ndarray
+        Covariance matrix of the residuals or innovations.
+        If sig_u is None, then an identity matrix is used.
+    steps : None or int
+        number of observations to simulate, this includes the initial
+        observations to start the autoregressive process.
+        If offset is not None, then exog of the model are used if they were
+        provided in the model
+    seed : None or integer
+        If seed is not None, then it will be used with for the random
+        variables generated by numpy.random.
+
+    Returns
+    -------
+    endog_simulated : nd_array
+        Endog of the simulated VAR process
+
     """
     rs = np.random.RandomState(seed=seed)
     rmvnorm = rs.multivariate_normal
     p, k, k = coefs.shape
+    if sig_u is None:
+        sig_u = np.eye(k)
     ugen = rmvnorm(np.zeros(len(sig_u)), sig_u, steps)
     result = np.zeros((steps, k))
     if intercept is not None:
-        result[p:] = intercept + ugen[p:]
+        # intercept can be 2-D like an offset variable
+        if np.ndim(intercept) > 1:
+            if not len(intercept) == len(ugen):
+                raise ValueError('2-D intercept needs to have length `steps`')
+        # add intercept/offset also to intial values
+        result += intercept
+        result[p:] += ugen[p:]
     else:
         result[p:] = ugen[p:]
 
diff --git a/statsmodels/tsa/vector_ar/var_model.py b/statsmodels/tsa/vector_ar/var_model.py
index 72eb09d2154..247c9bed849 100644
--- a/statsmodels/tsa/vector_ar/var_model.py
+++ b/statsmodels/tsa/vector_ar/var_model.py
@@ -759,22 +759,52 @@ class VARProcess(object):
     Parameters
     ----------
     coefs : ndarray (p x k x k)
-    intercept : ndarray (length k)
+        coefficients for lags of endog, part or params reshaped
+    coefs_exog : ndarray
+        parameters for trend and user provided exog
     sigma_u : ndarray (k x k)
+        residual covariance
     names : sequence (length k)
+    _params_info : dict
+        internal dict to provide information about the composition of `params`,
+        specifically `k_trend` (trend order) and `k_exog_user` (the number of
+        exog variables provided by the user).
+        If it is None, then coefs_exog are assumed to be for the intercept and
+        trend.
 
     Returns
     -------
     **Attributes**:
     """
-    def __init__(self, coefs, exog, sigma_u, names=None):
+    def __init__(self, coefs, coefs_exog, sigma_u, names=None, _params_info=None):
         self.k_ar = len(coefs)
         self.neqs = coefs.shape[1]
         self.coefs = coefs
-        self.exog = exog
+        self.coefs_exog = coefs_exog
+        # Note reshaping 1-D coefs_exog to 2_D makes unit tests fail
         self.sigma_u = sigma_u
         self.names = names
 
+        if _params_info is None:
+            _params_info = {}
+        self.k_exog_user = _params_info.get('k_exog_user', 0)
+        if self.coefs_exog is not None:
+            k_ex = self.coefs_exog.shape[0] if self.coefs_exog.ndim != 1 else 1
+            k_c = k_ex - self.k_exog_user
+        else:
+            k_c = 0
+        self.k_trend = _params_info.get('k_trend', k_c)
+        # TODO: we need to distinguish exog including trend and exog_user
+        self.k_exog = self.k_trend + self.k_exog_user
+
+        if self.k_trend > 0:
+            if coefs_exog.ndim == 2:
+                self.intercept = coefs_exog[:, 0]
+            else:
+                self.intercept = coefs_exog
+        else:
+            self.intercept = np.zeros(self.neqs)
+
     def get_eq_index(self, name):
         "Return integer position of requested equation name"
         return util.get_index(self.names, name)
@@ -802,22 +832,95 @@ def is_stable(self, verbose=False):
         """
         return is_stable(self.coefs, verbose=verbose)
 
-    def plotsim(self, steps=1000):
+    def simulate_var(self, steps=None, offset=None, seed=None):
+        """
+        simulate the VAR(p) process for the desired number of steps
+
+        Parameters
+        ----------
+        steps : None or int
+            number of observations to simulate, this includes the initial
+            observations to start the autoregressive process.
+            If offset is not None, then exog of the model are used if they were
+            provided in the model
+        offset : None or ndarray (steps, neqs)
+            If not None, then offset is added as an observation specific
+            intercept to the autoregression. If it is None and either trend
+            (including intercept) or exog were used in the VAR model, then
+            the linear predictor of those components will be used as offset.
+            This should have the same number of rows as steps, and the same
+            number of columns as endogenous variables (neqs).
+        seed : None or integer
+            If seed is not None, then it will be used with for the random
+            variables generated by numpy.random.
+
+        Returns
+        -------
+        endog_simulated : nd_array
+            Endog of the simulated VAR process
+
+        """
+        steps_ = None
+        if offset is None:
+            if self.k_exog_user > 0 or self.k_trend > 1:
+                # if more than intercept
+                # endog_lagged contains all regressors, trend, exog_user
+                # and lagged endog, trimmed initial observations
+                offset = self.endog_lagged[:,:self.k_exog].dot(
+                                                     self.coefs_exog.T)
+                steps_ = self.endog_lagged.shape[0]
+            else:
+                offset = self.intercept
+        else:
+            steps_ = offset.shape[0]
+
+        # default, but over written if exog or offset are used
+        if steps is None:
+            if steps_ is None:
+                steps = 1000
+            else:
+                steps = steps_
+        else:
+            if steps_ is not None and steps != steps_:
+                raise ValueError('if exog or offset are used, then steps must'
+                                 'be equal to their length or None')
+
+        y = util.varsim(self.coefs, offset, self.sigma_u, steps=steps, seed=seed)
+        return y
+
+    def plotsim(self, steps=None, offset=None, seed=None):
         """
         Plot a simulation from the VAR(p) process for the desired number of
         steps
         """
-        Y = util.varsim(self.coefs, self.exog, self.sigma_u, steps=steps)
-        plotting.plot_mts(Y)
+        y = self.simulate_var(steps=steps, offset=offset, seed=seed)
+        return plotting.plot_mts(y)
+
+    def intercept_longrun(self):
+        r"""Long run intercept of stable VAR process
+
+        Lütkepohl eq. 2.1.23
+
+        .. math:: \mu = (I - A_1 - \dots - A_p)^{-1} \alpha
+
+        where \alpha is the intercept (parameter of the constant)
+        """
+        return np.linalg.solve(self._char_mat, self.intercept)
 
     def mean(self):
-        r"""Mean of stable process
+        r"""Long run intercept of stable VAR process
+
+        Warning: trend and exog except for intercept are ignored for this.
+        This might change in future versions.
 
         Lütkepohl eq. 2.1.23
 
         .. math:: \mu = (I - A_1 - \dots - A_p)^{-1} \alpha
+
+        where \alpha is the intercept (parameter of the constant)
         """
-        return np.linalg.solve(self._char_mat, self.exog)
+        return self.intercept_longrun()
+
 
     def ma_rep(self, maxn=10):
         r"""Compute MA(:math:`\infty`) coefficient matrices
@@ -889,7 +992,9 @@ def acorr(self, nlags=None):
 
     def plot_acorr(self, nlags=10, linewidth=8):
         "Plot theoretical autocorrelation function"
-        plotting.plot_full_acorr(self.acorr(nlags=nlags), linewidth=linewidth)
+        fig = plotting.plot_full_acorr(self.acorr(nlags=nlags),
+                                       linewidth=linewidth)
+        return fig
 
     def forecast(self, y, steps, exog_future=None):
         """Produce linear minimum MSE forecasts for desired number of steps
@@ -1101,6 +1206,7 @@ def __init__(self, endog, endog_lagged, params, sigma_u, lag_order,
         k_trend = util.get_trendorder(trend)
         self.exog_names = util.make_lag_names(names, lag_order, k_trend, exog)
         self.params = params
+        self.exog = exog
         # print(params.shape)
         # print(params.round(3))
 
@@ -1109,28 +1215,39 @@ def __init__(self, endog, endog_lagged, params, sigma_u, lag_order,
         # Each matrix needs to be transposed
         endog_start = k_trend
         if exog is not None:
-            endog_start += exog.shape[1]
+            k_exog_user = exog.shape[1]
+            endog_start += k_exog_user
+        else:
+            k_exog_user = 0
         reshaped = self.params[endog_start:]
         reshaped = reshaped.reshape((lag_order, neqs, neqs))
         # Need to transpose each coefficient matrix
         coefs = reshaped.swapaxes(1, 2).copy()
 
         self.coefs_exog = params[:endog_start].T
-        self.k_trend = self.coefs_exog.shape[1]
-
-        # print(coefs.round(3))
-        super(VARResults, self).__init__(coefs, exog, sigma_u, names=names)
+        #self.k_trend = self.coefs_exog.shape[1]
+        self.k_exog = self.coefs_exog.shape[1]
+        self.k_exog_user = k_exog_user
+
+        # maybe change to params class, distinguish exog_all versus exog_user
+        # see issue #4535
+        _params_info = {'k_trend': k_trend,
+                        'k_exog_user': k_exog_user,
+                        'k_ar': lag_order}
+        super(VARResults, self).__init__(coefs, self.coefs_exog, sigma_u,
+                                         names=names,
+                                         _params_info=_params_info)
 
     def plot(self):
         """Plot input time series
         """
-        plotting.plot_mts(self.y, names=self.names, index=self.dates)
+        return plotting.plot_mts(self.y, names=self.names, index=self.dates)
 
     @property
     def df_model(self):
         """Number of estimated parameters, including the intercept / trends
         """
-        return self.neqs * self.k_ar + self.k_trend
+        return self.neqs * self.k_ar + self.k_exog
 
     @property
     def df_resid(self):
@@ -1158,8 +1275,9 @@ def sample_acorr(self, nlags=1):
 
     def plot_sample_acorr(self, nlags=10, linewidth=8):
         "Plot theoretical autocorrelation function"
-        plotting.plot_full_acorr(self.sample_acorr(nlags=nlags),
-                                 linewidth=linewidth)
+        fig = plotting.plot_full_acorr(self.sample_acorr(nlags=nlags),
+                                       linewidth=linewidth)
+        return fig
 
     def resid_acov(self, nlags=1):
         """
@@ -1246,7 +1364,7 @@ def _cov_alpha(self):
         Estimated covariance matrix of model coefficients w/o exog
         """
         # drop exog
-        return self.cov_params[self.k_trend*self.neqs:, self.k_trend*self.neqs:]
+        return self.cov_params[self.k_exog*self.neqs:, self.k_exog*self.neqs:]
 
     @cache_readonly
     def _cov_sigma(self):
@@ -1275,12 +1393,12 @@ def stderr(self):
 
     @cache_readonly
     def stderr_endog_lagged(self):
-        start = self.k_trend
+        start = self.k_exog
         return self.stderr[start:]
 
     @cache_readonly
     def stderr_dt(self):
-        end = self.k_trend
+        end = self.k_exog
         return self.stderr[:end]
 
     @cache_readonly
@@ -1292,12 +1410,12 @@ def tvalues(self):
 
     @cache_readonly
     def tvalues_endog_lagged(self):
-        start = self.k_trend
+        start = self.k_exog
         return self.tvalues[start:]
 
     @cache_readonly
     def tvalues_dt(self):
-        end = self.k_trend
+        end = self.k_exog
         return self.tvalues[:end]
 
     @cache_readonly
@@ -1309,12 +1427,12 @@ def pvalues(self):
 
     @cache_readonly
     def pvalues_endog_lagged(self):
-        start = self.k_trend
+        start = self.k_exog
         return self.pvalues[start:]
 
     @cache_readonly
     def pvalues_dt(self):
-        end = self.k_trend
+        end = self.k_exog
         return self.pvalues[:end]
 
     # todo: --------------------------------------------------------------------
@@ -1325,8 +1443,9 @@ def plot_forecast(self, steps, alpha=0.05, plot_stderr=True):
         """
         mid, lower, upper = self.forecast_interval(self.y[-self.k_ar:], steps,
                                                    alpha=alpha)
-        plotting.plot_var_forc(self.y, mid, lower, upper, names=self.names,
-                               plot_stderr=plot_stderr)
+        fig = plotting.plot_var_forc(self.y, mid, lower, upper,
+                                     names=self.names, plot_stderr=plot_stderr)
+        return fig
 
     # Forecast error covariance functions
 
@@ -1509,12 +1628,12 @@ def bpow(i):
 
     def _bmat_forc_cov(self):
         # B as defined on p. 96 of Lut
-        upper = np.zeros((self.k_trend, self.df_model))
-        upper[:, :self.k_trend] = np.eye(self.k_trend)
+        upper = np.zeros((self.k_exog, self.df_model))
+        upper[:, :self.k_exog] = np.eye(self.k_exog)
 
         lower_dim = self.neqs * (self.k_ar - 1)
         I = np.eye(lower_dim)
-        lower = np.column_stack((np.zeros((lower_dim, self.k_trend)), I,
+        lower = np.column_stack((np.zeros((lower_dim, self.k_exog)), I,
                                  np.zeros((lower_dim, self.neqs))))
 
         return np.vstack((upper, self.params.T, lower))
@@ -1656,7 +1775,7 @@ def test_causality(self, caused, causing=None, kind='f', signif=0.05):
 
         # number of restrictions
         num_restr = len(causing) * len(caused) * p
-        num_det_terms = self.k_trend
+        num_det_terms = self.k_exog
 
         # Make restriction matrix
         C = np.zeros((num_restr, k * num_det_terms + k**2 * p), dtype=float)
@@ -1804,9 +1923,9 @@ def test_inst_causality(self, causing, signif=0.05):
                                     crit_value, pvalue, df, signif,
                                     test="inst", method="wald")
 
-    def test_whiteness_new(self, nlags=10, signif=0.05, adjusted=False):
+    def test_whiteness(self, nlags=10, signif=0.05, adjusted=False):
         """
-        Residual whiteness tests using Portmanteau
+        Residual whiteness tests using Portmanteau test
 
         Parameters
         ----------
@@ -1846,35 +1965,41 @@ def test_whiteness_new(self, nlags=10, signif=0.05, adjusted=False):
         return WhitenessTestResults(statistic, crit_value, pvalue, df, signif,
                                     nlags, adjusted)
 
-    # TODO: This is not a formal test.  Deprecate in favor of `test_whiteness_new`
-    def test_whiteness(self, nlags=10, plot=True, linewidth=8):
+    def plot_acorr(self, nlags=10, resid=True, linewidth=8):
         """
-        Test white noise assumption. Sample (Y) autocorrelations are compared
-        with the standard :math:`2 / \sqrt{T}` bounds.
+        Plot autocorrelation of sample (endog) or residuals
+
+        Sample (Y) or Residual autocorrelations are plotted together with the
+        standard :math:`2 / \sqrt{T}` bounds.
 
         Parameters
         ----------
-        plot : boolean, default True
-            Plot autocorrelations with 2 / sqrt(T) bounds
-        """
-        acorrs = self.sample_acorr(nlags)
-        bound = 2 / np.sqrt(self.nobs)
+        nlags : int
+            number of lags to display (excluding 0)
+        resid: boolean
+            If True, then the autocorrelation of the residuals is plotted
+            If False, then the autocorrelation of endog is plotted.
+        linewidth : int
+            width of vertical bars
 
-        # TODO: this probably needs some UI work
+        Returns
+        -------
+        fig : matplotlib figure instance
 
-        if (np.abs(acorrs) > bound).any():
-            print('FAIL: Some autocorrelations exceed %.4f bound. '
-                   'See plot' % bound)
+        """
+        if resid:
+            acorrs = self.resid_acorr(nlags)
         else:
-            print('PASS: No autocorrelations exceed %.4f bound' % bound)
+            acorrs = self.sample_acorr(nlags)
+
+        bound = 2 / np.sqrt(self.nobs)
 
-        if plot:
-            fig = plotting.plot_full_acorr(acorrs[1:],
-                                           xlabel=np.arange(1, nlags+1),
-                                           err_bound=bound,
-                                           linewidth=linewidth)
-            fig.suptitle(r"ACF plots with $2 / \sqrt{T}$ bounds "
-                         "for testing whiteness assumption")
+        fig = plotting.plot_full_acorr(acorrs[1:],
+                                       xlabel=np.arange(1, nlags+1),
+                                       err_bound=bound,
+                                       linewidth=linewidth)
+        fig.suptitle(r"ACF plots for residuals with $2 / \sqrt{T}$ bounds ")
+        return fig
 
     def test_normality(self, signif=0.05):
         """
@@ -1914,7 +2039,7 @@ def info_criteria(self):
         nobs = self.nobs
         neqs = self.neqs
         lag_order = self.k_ar
-        free_params = lag_order * neqs ** 2 + neqs * self.k_trend
+        free_params = lag_order * neqs ** 2 + neqs * self.k_exog
 
         ld = logdet_symm(self.sigma_u_mle)
 
@@ -2077,6 +2202,7 @@ def plot(self, periods=None, figsize=(10, 10), **plot_kwds):
         handles, labels = ax.get_legend_handles_labels()
         fig.legend(handles, labels, loc='upper right')
         plotting.adjust_subplots(right=0.85)
+        return fig
 
 # -------------------------------------------------------------------------------
 


Code review comments:
diff_hunk: @@ -586,6 +586,7 @@ def __init__(self, endog, endog_lagged, params, sigma_u, lag_order,
         else:
             trendorder = None
         self.k_trend = k_trend
+        self.k_exog = k_trend  # now required by VARProcess
path: statsmodels/tsa/vector_ar/svar_model.py
text: @user1:
The "now" will be ambiguous in a few years.  Maybe "as of 0.9"?

@author:
"now" in the repository can always be discovered in the blame view to check the commit that added or changed a line, as long as the history is not made noisy.

@user1:
Writing that comment already took more effort than adding the suggested clarification.  The casual reader shouldn't have to check the blame when there's a self-documenting alternative.

diff_hunk: @@ -1101,6 +1206,7 @@ def __init__(self, endog, endog_lagged, params, sigma_u, lag_order,
         k_trend = util.get_trendorder(trend)
         self.exog_names = util.make_lag_names(names, lag_order, k_trend, exog)
         self.params = params
+        self.exog = exog
         # print(params.shape)
         # print(params.round(3))
path: statsmodels/tsa/vector_ar/var_model.py
text: @user1:
delete