PageRenderTime 89ms CodeModel.GetById 35ms app.highlight 46ms RepoModel.GetById 2ms app.codeStats 0ms

/statsmodels/tsa/vector_ar/tests/test_var.py

http://github.com/statsmodels/statsmodels
Python | 845 lines | 806 code | 31 blank | 8 comment | 9 complexity | 20e80069c621097a3025693467ab0d9d MD5 | raw file
  1# -*- coding: utf-8 -*-
  2"""
  3Test VAR Model
  4"""
  5from statsmodels.compat.pandas import assert_index_equal
  6
  7import os
  8import sys
  9import warnings
 10from io import StringIO, BytesIO
 11
 12import numpy as np
 13import pandas as pd
 14import pytest
 15from numpy.testing import (assert_almost_equal, assert_equal,
 16                           assert_allclose)
 17
 18import statsmodels.api as sm
 19import statsmodels.tools.data as data_util
 20import statsmodels.tsa.vector_ar.util as util
 21from statsmodels.compat.python import iteritems, lrange
 22from statsmodels.tools.sm_exceptions import ValueWarning
 23from statsmodels.tsa.vector_ar.var_model import VAR, var_acf
 24
 25DECIMAL_12 = 12
 26DECIMAL_6 = 6
 27DECIMAL_5 = 5
 28DECIMAL_4 = 4
 29DECIMAL_3 = 3
 30DECIMAL_2 = 2
 31
 32
 33@pytest.fixture()
 34def bivariate_var_data(reset_randomstate):
 35    """A bivariate dataset for VAR estimation"""
 36    e = np.random.standard_normal((252, 2))
 37    y = np.zeros_like(e)
 38    y[:2] = e[:2]
 39    for i in range(2, 252):
 40        y[i] = .2 * y[i - 1] + .1 * y[i - 2] + e[i]
 41    return y
 42
 43
 44@pytest.fixture()
 45def bivariate_var_result(bivariate_var_data):
 46    """A bivariate VARResults for reuse"""
 47    mod = VAR(bivariate_var_data)
 48    return mod.fit()
 49
 50
 51class CheckVAR(object):  # FIXME: not inherited, so these tests are never run!
 52    # just so pylint will not complain
 53    res1 = None
 54    res2 = None
 55
 56    def test_params(self):
 57        assert_almost_equal(self.res1.params, self.res2.params, DECIMAL_3)
 58
 59    def test_neqs(self):
 60        assert_equal(self.res1.neqs, self.res2.neqs)
 61
 62    def test_nobs(self):
 63        assert_equal(self.res1.avobs, self.res2.nobs)
 64
 65    def test_df_eq(self):
 66        assert_equal(self.res1.df_eq, self.res2.df_eq)
 67
 68    def test_rmse(self):
 69        results = self.res1.results
 70        for i in range(len(results)):
 71            assert_almost_equal(results[i].mse_resid**.5,
 72                    eval('self.res2.rmse_'+str(i+1)), DECIMAL_6)
 73
 74    def test_rsquared(self):
 75        results = self.res1.results
 76        for i in range(len(results)):
 77            assert_almost_equal(results[i].rsquared,
 78                    eval('self.res2.rsquared_'+str(i+1)), DECIMAL_3)
 79
 80    def test_llf(self):
 81        results = self.res1.results
 82        assert_almost_equal(self.res1.llf, self.res2.llf, DECIMAL_2)
 83        for i in range(len(results)):
 84            assert_almost_equal(results[i].llf,
 85                    eval('self.res2.llf_'+str(i+1)), DECIMAL_2)
 86
 87    def test_aic(self):
 88        assert_almost_equal(self.res1.aic, self.res2.aic)
 89
 90    def test_bic(self):
 91        assert_almost_equal(self.res1.bic, self.res2.bic)
 92
 93    def test_hqic(self):
 94        assert_almost_equal(self.res1.hqic, self.res2.hqic)
 95
 96    def test_fpe(self):
 97        assert_almost_equal(self.res1.fpe, self.res2.fpe)
 98
 99    def test_detsig(self):
100        assert_almost_equal(self.res1.detomega, self.res2.detsig)
101
102    def test_bse(self):
103        assert_almost_equal(self.res1.bse, self.res2.bse, DECIMAL_4)
104
105
106def get_macrodata():
107    data = sm.datasets.macrodata.load_pandas().data[['realgdp','realcons','realinv']]
108    data = data.to_records(index=False)
109    nd = data.view((float,3), type=np.ndarray)
110    nd = np.diff(np.log(nd), axis=0)
111    return nd.ravel().view(data.dtype, type=np.ndarray)
112
113
114def generate_var():  # FIXME: make a test?
115    from rpy2.robjects import r
116    import pandas.rpy.common as prp
117    r.source('tests/var.R')
118    return prp.convert_robj(r['result'], use_pandas=False)
119
120
121def write_generate_var():  # FIXME: make a test?
122    result = generate_var()
123    np.savez('tests/results/vars_results.npz', **result)
124
125
126class RResults(object):
127    """
128    Simple interface with results generated by "vars" package in R.
129    """
130
131    def __init__(self):
132        #data = np.load(resultspath + 'vars_results.npz')
133        from .results.results_var_data import var_results
134        data = var_results.__dict__
135
136        self.names = data['coefs'].dtype.names
137        self.params = data['coefs'].view((float, len(self.names)), type=np.ndarray)
138        self.stderr = data['stderr'].view((float, len(self.names)), type=np.ndarray)
139
140        self.irf = data['irf'].item()
141        self.orth_irf = data['orthirf'].item()
142
143        self.nirfs = int(data['nirfs'][0])
144        self.nobs = int(data['obs'][0])
145        self.totobs = int(data['totobs'][0])
146
147        crit = data['crit'].item()
148        self.aic = crit['aic'][0]
149        self.sic = self.bic = crit['sic'][0]
150        self.hqic = crit['hqic'][0]
151        self.fpe = crit['fpe'][0]
152
153        self.detomega = data['detomega'][0]
154        self.loglike = data['loglike'][0]
155
156        self.nahead = int(data['nahead'][0])
157        self.ma_rep = data['phis']
158
159        self.causality = data['causality']
160
161
162_orig_stdout = None
163
164
165def setup_module():
166    global _orig_stdout
167    _orig_stdout = sys.stdout
168    sys.stdout = StringIO()
169
170
171class CheckIRF(object):
172
173    ref = None
174    res = None
175    irf = None
176    k = None
177
178    #---------------------------------------------------------------------------
179    # IRF tests
180
181    def test_irf_coefs(self):
182        self._check_irfs(self.irf.irfs, self.ref.irf)
183        self._check_irfs(self.irf.orth_irfs, self.ref.orth_irf)
184
185    def _check_irfs(self, py_irfs, r_irfs):
186        for i, name in enumerate(self.res.names):
187            ref_irfs = r_irfs[name].view((float, self.k), type=np.ndarray)
188            res_irfs = py_irfs[:, :, i]
189            assert_almost_equal(ref_irfs, res_irfs)
190
191    @pytest.mark.matplotlib
192    def test_plot_irf(self, close_figures):
193        self.irf.plot()
194        self.irf.plot(plot_stderr=False)
195
196        self.irf.plot(impulse=0, response=1)
197        self.irf.plot(impulse=0)
198        self.irf.plot(response=0)
199
200        self.irf.plot(orth=True)
201        self.irf.plot(impulse=0, response=1, orth=True)
202
203    @pytest.mark.matplotlib
204    def test_plot_cum_effects(self, close_figures):
205        self.irf.plot_cum_effects()
206        self.irf.plot_cum_effects(plot_stderr=False)
207        self.irf.plot_cum_effects(impulse=0, response=1)
208
209        self.irf.plot_cum_effects(orth=True)
210        self.irf.plot_cum_effects(impulse=0, response=1, orth=True)
211
212
213@pytest.mark.smoke
214class CheckFEVD(object):
215
216    fevd = None
217
218    #---------------------------------------------------------------------------
219    # FEVD tests
220
221    @pytest.mark.matplotlib
222    def test_fevd_plot(self, close_figures):
223        self.fevd.plot()
224
225    def test_fevd_repr(self):
226        self.fevd
227
228    def test_fevd_summary(self):
229        self.fevd.summary()
230
231    @pytest.mark.xfail(reason="FEVD.cov() is not implemented",
232                       raises=NotImplementedError, strict=True)
233    def test_fevd_cov(self):
234        # test does not crash
235        # not implemented
236        covs = self.fevd.cov()
237        raise NotImplementedError
238
239
240class TestVARResults(CheckIRF, CheckFEVD):
241
242    @classmethod
243    def setup_class(cls):
244        cls.p = 2
245
246        cls.data = get_macrodata()
247        cls.model = VAR(cls.data)
248        cls.names = cls.model.endog_names
249
250        cls.ref = RResults()
251        cls.k = len(cls.ref.names)
252        cls.res = cls.model.fit(maxlags=cls.p)
253
254        cls.irf = cls.res.irf(cls.ref.nirfs)
255        cls.nahead = cls.ref.nahead
256
257        cls.fevd = cls.res.fevd()
258
259    def test_constructor(self):
260        # make sure this works with no names
261        ndarr = self.data.view((float, 3), type=np.ndarray)
262        model = VAR(ndarr)
263        res = model.fit(self.p)
264
265    def test_names(self):
266        assert_equal(self.model.endog_names, self.ref.names)
267
268        model2 = VAR(self.data)
269        assert_equal(model2.endog_names, self.ref.names)
270
271    def test_get_eq_index(self):
272        assert type(self.res.names) is list  # noqa: E721
273
274        for i, name in enumerate(self.names):
275            idx = self.res.get_eq_index(i)
276            idx2 = self.res.get_eq_index(name)
277
278            assert_equal(idx, i)
279            assert_equal(idx, idx2)
280
281        with pytest.raises(Exception):
282            self.res.get_eq_index('foo')
283
284    @pytest.mark.smoke
285    def test_repr(self):
286        # just want this to work
287        foo = str(self.res)
288        bar = repr(self.res)
289
290    def test_params(self):
291        assert_almost_equal(self.res.params, self.ref.params, DECIMAL_3)
292
293    @pytest.mark.smoke
294    def test_cov_params(self):
295        # do nothing for now
296        self.res.cov_params
297
298    @pytest.mark.smoke
299    def test_cov_ybar(self):
300        self.res.cov_ybar()
301
302    @pytest.mark.smoke
303    def test_tstat(self):
304        self.res.tvalues
305
306    @pytest.mark.smoke
307    def test_pvalues(self):
308        self.res.pvalues
309
310    @pytest.mark.smoke
311    def test_summary(self):
312        summ = self.res.summary()
313
314    def test_detsig(self):
315        assert_almost_equal(self.res.detomega, self.ref.detomega)
316
317    def test_aic(self):
318        assert_almost_equal(self.res.aic, self.ref.aic)
319
320    def test_bic(self):
321        assert_almost_equal(self.res.bic, self.ref.bic)
322
323    def test_hqic(self):
324        assert_almost_equal(self.res.hqic, self.ref.hqic)
325
326    def test_fpe(self):
327        assert_almost_equal(self.res.fpe, self.ref.fpe)
328
329    def test_lagorder_select(self):
330        ics = ['aic', 'fpe', 'hqic', 'bic']
331
332        for ic in ics:
333            res = self.model.fit(maxlags=10, ic=ic, verbose=True)
334
335        with pytest.raises(Exception):
336            self.model.fit(ic='foo')
337
338    def test_nobs(self):
339        assert_equal(self.res.nobs, self.ref.nobs)
340
341    def test_stderr(self):
342        assert_almost_equal(self.res.stderr, self.ref.stderr, DECIMAL_4)
343
344    def test_loglike(self):
345        assert_almost_equal(self.res.llf, self.ref.loglike)
346
347    def test_ma_rep(self):
348        ma_rep = self.res.ma_rep(self.nahead)
349        assert_almost_equal(ma_rep, self.ref.ma_rep)
350
351    #--------------------------------------------------
352    # Lots of tests to make sure stuff works...need to check correctness
353
354    def test_causality(self):
355        causedby = self.ref.causality['causedby']
356
357        for i, name in enumerate(self.names):
358            variables = self.names[:i] + self.names[i + 1:]
359            result = self.res.test_causality(name, variables, kind='f')
360            assert_almost_equal(result.pvalue, causedby[i], DECIMAL_4)
361
362            rng = lrange(self.k)
363            rng.remove(i)
364            result2 = self.res.test_causality(i, rng, kind='f')
365            assert_almost_equal(result.pvalue, result2.pvalue, DECIMAL_12)
366
367            # make sure works
368            result = self.res.test_causality(name, variables, kind='wald')
369
370        # corner cases
371        _ = self.res.test_causality(self.names[0], self.names[1])
372        _ = self.res.test_causality(0, 1)
373
374        with pytest.raises(Exception):
375            self.res.test_causality(0, 1, kind='foo')
376
377    def test_causality_no_lags(self):
378        res = VAR(self.data).fit(maxlags=0)
379        with pytest.raises(RuntimeError, match="0 lags"):
380            res.test_causality(0, 1)
381
382    @pytest.mark.smoke
383    def test_select_order(self):
384        result = self.model.fit(10, ic='aic', verbose=True)
385        result = self.model.fit(10, ic='fpe', verbose=True)
386
387        # bug
388        model = VAR(self.model.endog)
389        model.select_order()
390
391    def test_is_stable(self):
392        # may not necessarily be true for other datasets
393        assert(self.res.is_stable(verbose=True))
394
395    def test_acf(self):
396        # test that it works...for now
397        acfs = self.res.acf(10)
398
399        # defaults to nlags=lag_order
400        acfs = self.res.acf()
401        assert(len(acfs) == self.p + 1)
402
403    def test_acf_2_lags(self):
404        c = np.zeros((2, 2, 2))
405        c[0] = np.array([[.2, .1], [.15, .15]])
406        c[1] = np.array([[.1, .9], [0, .1]])
407
408        acf = var_acf(c, np.eye(2), 3)
409
410        gamma = np.zeros((6, 6))
411        gamma[:2, :2] = acf[0]
412        gamma[2:4, 2:4] = acf[0]
413        gamma[4:6, 4:6] = acf[0]
414        gamma[2:4, :2] = acf[1].T
415        gamma[4:, :2] = acf[2].T
416        gamma[:2, 2:4] = acf[1]
417        gamma[:2, 4:] = acf[2]
418        recovered = np.dot(gamma[:2, 2:], np.linalg.inv(gamma[:4, :4]))
419        recovered = [recovered[:, 2 * i:2 * (i + 1)] for i in range(2)]
420        recovered = np.array(recovered)
421        assert_allclose(recovered, c, atol=1e-7)
422
423    @pytest.mark.smoke
424    def test_acorr(self):
425        acorrs = self.res.acorr(10)
426
427    @pytest.mark.smoke
428    def test_forecast(self):
429        self.res.forecast(self.res.endog[-5:], 5)
430
431    @pytest.mark.smoke
432    def test_forecast_interval(self):
433        y = self.res.endog[:-self.p:]
434        point, lower, upper = self.res.forecast_interval(y, 5)
435
436    @pytest.mark.matplotlib
437    def test_plot_sim(self, close_figures):
438        self.res.plotsim(steps=100)
439
440    @pytest.mark.matplotlib
441    def test_plot(self, close_figures):
442        self.res.plot()
443
444    @pytest.mark.matplotlib
445    def test_plot_acorr(self, close_figures):
446        self.res.plot_acorr()
447
448    @pytest.mark.matplotlib
449    def test_plot_forecast(self, close_figures):
450        self.res.plot_forecast(5)
451
452    def test_reorder(self):
453        # manually reorder
454        data = self.data.view((float,3), type=np.ndarray)
455        names = self.names
456        data2 = np.append(np.append(data[:,2,None], data[:,0,None], axis=1), data[:,1,None], axis=1)
457        names2 = []
458        names2.append(names[2])
459        names2.append(names[0])
460        names2.append(names[1])
461        res2 = VAR(data2).fit(maxlags=self.p)
462
463        #use reorder function
464        res3 = self.res.reorder(['realinv','realgdp', 'realcons'])
465
466        #check if the main results match
467        assert_almost_equal(res2.params, res3.params)
468        assert_almost_equal(res2.sigma_u, res3.sigma_u)
469        assert_almost_equal(res2.bic, res3.bic)
470        assert_almost_equal(res2.stderr, res3.stderr)
471
472    def test_pickle(self):
473        fh = BytesIO()
474        #test wrapped results load save pickle
475        del self.res.model.data.orig_endog
476        self.res.save(fh)
477        fh.seek(0,0)
478        res_unpickled = self.res.__class__.load(fh)
479        assert type(res_unpickled) is type(self.res)  # noqa: E721
480
481
482class E1_Results(object):
483    """
484    Results from L├╝tkepohl (2005) using E2 dataset
485    """
486
487    def __init__(self):
488        # Lutkepohl p. 120 results
489
490        # I asked the author about these results and there is probably rounding
491        # error in the book, so I adjusted these test results to match what is
492        # coming out of the Python (double-checked) calculations
493        self.irf_stderr = np.array([[[.125, 0.546, 0.664 ],
494                                     [0.032, 0.139, 0.169],
495                                     [0.026, 0.112, 0.136]],
496
497                                    [[0.129, 0.547, 0.663],
498                                     [0.032, 0.134, 0.163],
499                                     [0.026, 0.108, 0.131]],
500
501                                    [[0.084, .385, .479],
502                                     [.016, .079, .095],
503                                     [.016, .078, .103]]])
504
505        self.cum_irf_stderr = np.array([[[.125, 0.546, 0.664 ],
506                                         [0.032, 0.139, 0.169],
507                                         [0.026, 0.112, 0.136]],
508
509                                        [[0.149, 0.631, 0.764],
510                                         [0.044, 0.185, 0.224],
511                                         [0.033, 0.140, 0.169]],
512
513                                        [[0.099, .468, .555],
514                                         [.038, .170, .205],
515                                         [.033, .150, .185]]])
516
517        self.lr_stderr = np.array([[.134, .645, .808],
518                                   [.048, .230, .288],
519                                   [.043, .208, .260]])
520
521
522basepath = os.path.split(sm.__file__)[0]
523resultspath = basepath + '/tsa/vector_ar/tests/results/'
524
525
526def get_lutkepohl_data(name='e2'):
527    path = resultspath + '%s.dat' % name
528
529    return util.parse_lutkepohl_data(path)
530
531
532def test_lutkepohl_parse():
533    files = ['e%d' % i for i in range(1, 7)]
534
535    for f in files:
536        get_lutkepohl_data(f)
537
538
539class TestVARResultsLutkepohl(object):
540    """
541    Verify calculations using results from L├╝tkepohl's book
542    """
543
544    @classmethod
545    def setup_class(cls):
546        cls.p = 2
547        sdata, dates = get_lutkepohl_data('e1')
548
549        data = data_util.struct_to_ndarray(sdata)
550        adj_data = np.diff(np.log(data), axis=0)
551        # est = VAR(adj_data, p=2, dates=dates[1:], names=names)
552
553        cls.model = VAR(adj_data[:-16], dates=dates[1:-16], freq='BQ-MAR')
554        cls.res = cls.model.fit(maxlags=cls.p)
555        cls.irf = cls.res.irf(10)
556        cls.lut = E1_Results()
557
558    def test_approx_mse(self):
559        # 3.5.18, p. 99
560        mse2 = np.array([[25.12, .580, 1.300],
561                         [.580, 1.581, .586],
562                         [1.300, .586, 1.009]]) * 1e-4
563
564        assert_almost_equal(mse2, self.res.forecast_cov(3)[1],
565                            DECIMAL_3)
566
567    def test_irf_stderr(self):
568        irf_stderr = self.irf.stderr(orth=False)
569        for i in range(1, 1 + len(self.lut.irf_stderr)):
570            assert_almost_equal(np.round(irf_stderr[i], 3),
571                                self.lut.irf_stderr[i-1])
572
573    def test_cum_irf_stderr(self):
574        stderr = self.irf.cum_effect_stderr(orth=False)
575        for i in range(1, 1 + len(self.lut.cum_irf_stderr)):
576            assert_almost_equal(np.round(stderr[i], 3),
577                                self.lut.cum_irf_stderr[i-1])
578
579    def test_lr_effect_stderr(self):
580        stderr = self.irf.lr_effect_stderr(orth=False)
581        orth_stderr = self.irf.lr_effect_stderr(orth=True)
582        assert_almost_equal(np.round(stderr, 3), self.lut.lr_stderr)
583
584
585def test_get_trendorder():
586    results = {
587        'c' : 1,
588        'nc' : 0,
589        'ct' : 2,
590        'ctt' : 3
591    }
592
593    for t, trendorder in iteritems(results):
594        assert(util.get_trendorder(t) == trendorder)
595
596
597def test_var_constant():
598    # see 2043
599    import datetime
600    from pandas import DataFrame, DatetimeIndex
601
602    series = np.array([[2., 2.], [1, 2.], [1, 2.], [1, 2.], [1., 2.]])
603    data = DataFrame(series)
604
605    d = datetime.datetime.now()
606    delta = datetime.timedelta(days=1)
607    index = []
608    for i in range(data.shape[0]):
609        index.append(d)
610        d += delta
611
612    data.index = DatetimeIndex(index)
613
614    #with pytest.warns(ValueWarning):  #does not silence warning in test output
615    with warnings.catch_warnings():
616        warnings.simplefilter("ignore", category=ValueWarning)
617        model = VAR(data)
618    with pytest.raises(ValueError):
619        model.fit(1)
620
621
622def test_var_trend():
623    # see 2271
624    data = get_macrodata().view((float,3), type=np.ndarray)
625
626    model = sm.tsa.VAR(data)
627    results = model.fit(4) #, trend = 'c')
628    irf = results.irf(10)
629
630    data_nc = data - data.mean(0)
631    model_nc = sm.tsa.VAR(data_nc)
632    results_nc = model_nc.fit(4, trend = 'nc')
633    with pytest.raises(ValueError):
634        model.fit(4, trend='t')
635
636
637def test_irf_trend():
638    # test for irf with different trend see #1636
639    # this is a rough comparison by adding trend or subtracting mean to data
640    # to get similar AR coefficients and IRF
641    data = get_macrodata().view((float,3), type=np.ndarray)
642
643    model = sm.tsa.VAR(data)
644    results = model.fit(4) #, trend = 'c')
645    irf = results.irf(10)
646
647    data_nc = data - data.mean(0)
648    model_nc = sm.tsa.VAR(data_nc)
649    results_nc = model_nc.fit(4, trend = 'nc')
650    irf_nc = results_nc.irf(10)
651
652    assert_allclose(irf_nc.stderr()[1:4], irf.stderr()[1:4], rtol=0.01)
653
654    trend = 1e-3 * np.arange(len(data)) / (len(data) - 1)
655    # for pandas version, currently not used, if data is a pd.DataFrame
656    #data_t = pd.DataFrame(data.values + trend[:,None], index=data.index, columns=data.columns)
657    data_t = data + trend[:,None]
658
659    model_t = sm.tsa.VAR(data_t)
660    results_t = model_t.fit(4, trend = 'ct')
661    irf_t = results_t.irf(10)
662
663    assert_allclose(irf_t.stderr()[1:4], irf.stderr()[1:4], rtol=0.03)
664
665
666class TestVARExtras(object):
667
668    @classmethod
669    def setup_class(cls):
670        mdata = sm.datasets.macrodata.load_pandas().data
671        mdata = mdata[['realgdp','realcons','realinv']]
672        data = mdata.values
673        data = np.diff(np.log(data), axis=0) * 400
674        cls.res0 = sm.tsa.VAR(data).fit(maxlags=2)
675
676    def test_process(self, close_figures):
677        res0 = self.res0
678        k_ar = res0.k_ar
679        fc20 = res0.forecast(res0.endog[-k_ar:], 20)
680        mean_lr = res0.mean()
681        assert_allclose(mean_lr, fc20[-1], rtol=5e-4)
682
683        ysim = res0.simulate_var(seed=987128)
684        assert_allclose(ysim.mean(0), mean_lr, rtol=0.1)
685        # initialization does not use long run intercept, see #4542
686        assert_allclose(ysim[0], res0.intercept, rtol=1e-10)
687        assert_allclose(ysim[1], res0.intercept, rtol=1e-10)
688
689        n_sim = 900
690        ysimz = res0.simulate_var(steps=n_sim, offset=np.zeros((n_sim, 3)),
691                                  seed=987128)
692        zero3 = np.zeros(3)
693        assert_allclose(ysimz.mean(0), zero3, atol=0.4)
694        # initialization does not use long run intercept, see #4542
695        assert_allclose(ysimz[0], zero3, atol=1e-10)
696        assert_allclose(ysimz[1], zero3, atol=1e-10)
697
698        # check attributes
699        assert_equal(res0.k_trend, 1)
700        assert_equal(res0.k_exog_user, 0)
701        assert_equal(res0.k_exog, 1)
702        assert_equal(res0.k_ar, 2)
703
704        irf = res0.irf()
705
706    @pytest.mark.matplotlib
707    def test_process_plotting(self, close_figures):
708        # Partially a smoke test
709        res0 = self.res0
710        k_ar = res0.k_ar
711        fc20 = res0.forecast(res0.endog[-k_ar:], 20)
712        irf = res0.irf()
713
714        res0.plotsim()
715        res0.plot_acorr()
716
717        fig = res0.plot_forecast(20)
718        fcp = fig.axes[0].get_children()[1].get_ydata()[-20:]
719        # Note values are equal, but keep rtol buffer
720        assert_allclose(fc20[:, 0], fcp, rtol=1e-13)
721        fcp = fig.axes[1].get_children()[1].get_ydata()[-20:]
722        assert_allclose(fc20[:, 1], fcp, rtol=1e-13)
723        fcp = fig.axes[2].get_children()[1].get_ydata()[-20:]
724        assert_allclose(fc20[:, 2], fcp, rtol=1e-13)
725
726        fig_asym = irf.plot()
727        fig_mc = irf.plot(stderr_type='mc', repl=1000, seed=987128)
728
729        for k in range(3):
730            a = fig_asym.axes[1].get_children()[k].get_ydata()
731            m = fig_mc.axes[1].get_children()[k].get_ydata()
732            # use m as desired because it is larger
733            # a is for some irf much smaller than m
734            assert_allclose(a, m, atol=0.1, rtol=0.9)
735
736    def test_forecast_cov(self):
737        # forecast_cov can include parameter uncertainty if contant-only
738        res = self.res0
739
740        covfc1 = res.forecast_cov(3)
741        assert_allclose(covfc1, res.mse(3), rtol=1e-13)
742        # ignore warning, TODO: assert OutputWarning
743        with warnings.catch_warnings():
744            warnings.simplefilter("ignore")
745            covfc2 = res.forecast_cov(3, method='auto')
746        assert_allclose(covfc2, covfc1, rtol=0.05)
747        # regression test, TODO: replace with verified numbers (Stata)
748        res_covfc2 = np.array([[[ 9.45802013,  4.94142038,  37.1999646 ],
749                                [ 4.94142038,  7.09273624,   5.66215089],
750                                [37.1999646 ,  5.66215089, 259.61275869]],
751
752                               [[11.30364479,  5.72569141,  49.28744123],
753                                [ 5.72569141,  7.409761  ,  10.98164091],
754                                [49.28744123, 10.98164091, 336.4484723 ]],
755
756                               [[12.36188803,  6.44426905,  53.54588026],
757                                [ 6.44426905,  7.88850029,  13.96382545],
758                                [53.54588026, 13.96382545, 352.19564327]]])
759        assert_allclose(covfc2, res_covfc2, atol=1e-6)
760
761    def test_exog(self):
762        # check that trend and exog are equivalent for basics and varsim
763        data = self.res0.model.endog
764        res_lin_trend = sm.tsa.VAR(data).fit(maxlags=2, trend='ct')
765        ex = np.arange(len(data))
766        res_lin_trend1 = sm.tsa.VAR(data, exog=ex).fit(maxlags=2)
767        ex2 = np.arange(len(data))[:, None]**[0, 1]
768        res_lin_trend2 = sm.tsa.VAR(data, exog=ex2).fit(maxlags=2, trend='nc')
769        # TODO: intercept differs by 4e-3, others are < 1e-12
770        assert_allclose(res_lin_trend.params, res_lin_trend1.params, rtol=5e-3)
771        assert_allclose(res_lin_trend.params, res_lin_trend2.params, rtol=5e-3)
772        assert_allclose(res_lin_trend1.params, res_lin_trend2.params, rtol=1e-10)
773
774        y1 = res_lin_trend.simulate_var(seed=987128)
775        y2 = res_lin_trend1.simulate_var(seed=987128)
776        y3 = res_lin_trend2.simulate_var(seed=987128)
777        assert_allclose(y2.mean(0), y1.mean(0), rtol=1e-12)
778        assert_allclose(y3.mean(0), y1.mean(0), rtol=1e-12)
779        assert_allclose(y3.mean(0), y2.mean(0), rtol=1e-12)
780
781        h = 10
782        fc1 = res_lin_trend.forecast(res_lin_trend.endog[-2:], h)
783        exf = np.arange(len(data), len(data) + h)
784        fc2 = res_lin_trend1.forecast(res_lin_trend1.endog[-2:], h,
785                                      exog_future=exf)
786        exf2 = exf[:, None]**[0, 1]
787        fc3 = res_lin_trend2.forecast(res_lin_trend2.endog[-2:], h,
788                                      exog_future=exf2)
789        assert_allclose(fc2, fc1, rtol=1e-12)
790        assert_allclose(fc3, fc1, rtol=1e-12)
791        assert_allclose(fc3, fc2, rtol=1e-12)
792
793        fci1 = res_lin_trend.forecast_interval(res_lin_trend.endog[-2:], h)
794        exf = np.arange(len(data), len(data) + h)
795        fci2 = res_lin_trend1.forecast_interval(res_lin_trend1.endog[-2:], h,
796                                                exog_future=exf)
797        exf2 = exf[:, None]**[0, 1]
798        fci3 = res_lin_trend2.forecast_interval(res_lin_trend2.endog[-2:], h,
799                                               exog_future=exf2)
800        assert_allclose(fci2, fci1, rtol=1e-12)
801        assert_allclose(fci3, fci1, rtol=1e-12)
802        assert_allclose(fci3, fci2, rtol=1e-12)
803
804
805@pytest.mark.parametrize('attr', ['y', 'ys_lagged'])
806def test_deprecated_attributes_varresults(bivariate_var_result, attr):
807    with pytest.warns(FutureWarning):
808        getattr(bivariate_var_result, attr)
809
810
811def test_var_cov_params_pandas(bivariate_var_data):
812    df = pd.DataFrame(bivariate_var_data, columns=['x', 'y'])
813    mod = VAR(df)
814    res = mod.fit(2)
815    cov = res.cov_params()
816    assert isinstance(cov, pd.DataFrame)
817    exog_names = ('const', 'L1.x', 'L1.y', 'L2.x', 'L2.y')
818    index = pd.MultiIndex.from_product((exog_names, ('x', 'y')))
819    assert_index_equal(cov.index, cov.columns)
820    assert_index_equal(cov.index, index)
821
822
823def test_summaries_exog(reset_randomstate):
824    y = np.random.standard_normal((500, 6))
825    df = pd.DataFrame(y)
826    cols = (["endog_{0}".format(i) for i in range(2)] +
827            ["exog_{0}".format(i) for i in range(4)])
828    df.columns = cols
829    df.index = pd.date_range('1-1-1950', periods=500, freq="MS")
830    endog = df.iloc[:, :2]
831    exog = df.iloc[:, 2:]
832
833    res = VAR(endog=endog, exog=exog).fit(maxlags=0)
834    summ = res.summary().summary
835    assert 'exog_0' in summ
836    assert 'exog_1' in summ
837    assert 'exog_2' in summ
838    assert 'exog_3' in summ
839
840    res = VAR(endog=endog, exog=exog).fit(maxlags=2)
841    summ = res.summary().summary
842    assert 'exog_0' in summ
843    assert 'exog_1' in summ
844    assert 'exog_2' in summ
845    assert 'exog_3' in summ