1# -*- coding: utf-8 -*-
2"""Test Johansen's Cointegration test against jplv, Spatial Econometrics Toolbox
3
4Created on Thu Aug 30 21:51:08 2012
5Author: Josef Perktold
6
7"""
8import os
9import warnings
10
11import numpy as np
12from numpy.testing import assert_almost_equal, assert_equal
13import pandas as pd
14import pytest
15
16from statsmodels.tools.sm_exceptions import HypothesisTestWarning
17from statsmodels.tsa.vector_ar.vecm import coint_johansen
18
19current_path = os.path.dirname(os.path.abspath(__file__))
20dta_path = os.path.join(current_path, "Matlab_results", "test_coint.csv")
21with open(dta_path, "rb") as fd:
22    dta = np.genfromtxt(fd)
23
24
25class CheckCointJoh(object):
26
27    def test_basic(self):
28        assert_equal(self.res.ind, np.arange(len(self.res.ind), dtype=int))
29        assert_equal(self.res.r0t.shape, (self.nobs_r, 8))
30
31    def test_table_trace(self):
32        table1 = np.column_stack((self.res.lr1, self.res.cvt))
33        assert_almost_equal(table1,
34                            self.res1_m.reshape(table1.shape, order='F'))
35
36    def test_table_maxeval(self):
37        table2 = np.column_stack((self.res.lr2, self.res.cvm))
38        assert_almost_equal(table2,
39                            self.res2_m.reshape(table2.shape, order='F'))
40
41    def test_normalization(self):
42        # GH 5517
43        evec = self.res.evec
44        non_zero = evec.flat != 0
45        assert evec.flat[non_zero][0] > 0
46
47
48class TestCointJoh12(CheckCointJoh):
49
50    @classmethod
51    def setup_class(cls):
52        cls.res = coint_johansen(dta, 1, 2)
53        cls.nobs_r = 173 - 1 - 2
54
55        cls.res1_m = np.array([241.985452556075,  166.4781461662553,  110.3298006342814,  70.79801574443575,  44.90887371527634,  27.22385073668511,  11.74205493173769,  3.295435325623445,           169.0618,           133.7852,           102.4674,            75.1027,            51.6492,            32.0645,            16.1619,             2.7055,           175.1584,            139.278,           107.3429,  79.34220000000001,            55.2459,            35.0116,            18.3985,             3.8415,           187.1891,           150.0778,           116.9829,            87.7748,            62.5202,            41.0815,            23.1485,             6.6349])
56        cls.res2_m = np.array([75.50730638981975,  56.14834553197396,   39.5317848898456,   25.8891420291594,  17.68502297859124,  15.48179580494741,  8.446619606114249,  3.295435325623445,            52.5858,            46.5583,            40.5244,            34.4202,            28.2398,            21.8731,            15.0006,             2.7055,            55.7302,            49.5875,            43.4183,            37.1646,            30.8151,            24.2522,            17.1481,             3.8415,            62.1741,            55.8171,            49.4095,            42.8612,             36.193,            29.2631,            21.7465,             6.6349,])
57
58        evec = np.array([
59            0.01102517075074406, -0.2185481584930077, 0.04565819524210763, -0.06556394587400775,
60            0.04711496306104131, -0.1500111976629196, 0.03775327003706507, 0.03479475877437702,
61
62            0.007517888890275335, -0.2014629352546497, 0.01526001455616041, 0.0707900418057458,
63            -0.002388919695513273, 0.04486516694838273, -0.02936314422571188, 0.009900554050392113,
64
65            0.02846074144367176, 0.02021385478834498, -0.04276914888645468, 0.1738024290422287,
66            0.07821155002012749, -0.1066523077111768, -0.3011042488399306, 0.04965189679477353,
67
68            0.07141291326159237, -0.01406702689857725, -0.07842109866080313, -0.04773566072362181,
69            -0.04768640728128824, -0.04428737926285261, 0.4143225656833862, 0.04512787132114879,
70
71            -0.06817130121837202, 0.2246249779872569, -0.009356548567565763, 0.006685350535849125,
72            -0.02040894506833539, 0.008131690308487425, -0.2503209797396666, 0.01560186979508953,
73
74            0.03327070126502506, -0.263036624535624, -0.04669882107497259, 0.0146457545413255,
75            0.01408691619062709, 0.1004753600191269, -0.02239205763487946, -0.02169291468272568,
76
77            0.08782313160608619, -0.07696508791577318, 0.008925177304198475, -0.06230900392092828,
78            -0.01548907461158638, 0.04574831652028973, -0.2972228156126774, 0.003469819004961912,
79
80            -0.001868995544352928, 0.05993345996347871, 0.01213394328069316, 0.02096614212178651,
81            -0.08624395993789938, 0.02108183181049973, -0.08470307289295617, -5.135072530480897e-005])
82        cls.evec_m = evec.reshape(cls.res.evec.shape, order='F')
83
84        cls.eig_m = np.array([0.3586376068088151, 0.2812806889719111, 0.2074818815675726,  0.141259991767926, 0.09880133062878599, 0.08704563854307619,  0.048471840356709, 0.01919823444066367])
85
86    def test_evec(self):
87        for col in range(self.evec_m.shape[1]):
88            try:
89                assert_almost_equal(self.res.evec[:, col],
90                                    self.evec_m[:, col])
91            except AssertionError:
92                assert_almost_equal(self.res.evec[:, col],
93                                    -self.evec_m[:, col])
94
95    def test_evals(self):
96        assert_almost_equal(self.res.eig, self.eig_m)
97
98
99class TestCointJoh09(CheckCointJoh):
100
101    @classmethod
102    def setup_class(cls):
103        cls.res = coint_johansen(dta, 0, 9)
104        cls.nobs_r = 173 - 1 - 9
105        #fprintf(1, '%18.16g, ', r1)
106        cls.res1_m = np.array([307.6888935095814,  205.3839229398245,  129.1330243009336,   83.3101865760208,  52.51955460357912,  30.20027050520502,  13.84158157562689, 0.4117390188204866,           153.6341,           120.3673,             91.109,            65.8202,            44.4929,            27.0669,            13.4294,             2.7055,            159.529,           125.6185,            95.7542,            69.8189,            47.8545,            29.7961,            15.4943,             3.8415,           171.0905,           135.9825,           104.9637,            77.8202,            54.6815,            35.4628,            19.9349,             6.6349])
107        #r2 = [res.lr2 res.cvm]
108        cls.res2_m = np.array([102.3049705697569,  76.25089863889085,  45.82283772491284,   30.7906319724417,  22.31928409837409,  16.35868892957814,   13.4298425568064, 0.4117390188204866,            49.2855,            43.2947,            37.2786,            31.2379,            25.1236,            18.8928,            12.2971,             2.7055,            52.3622,            46.2299,            40.0763,            33.8777,            27.5858,            21.1314,            14.2639,             3.8415,            58.6634,            52.3069,            45.8662,            39.3693,            32.7172,             25.865,              18.52,             6.6349])
109
110
111class TestCointJohMin18(CheckCointJoh):
112
113    @classmethod
114    def setup_class(cls):
115        cls.res = coint_johansen(dta, -1, 8)
116        cls.nobs_r = 173 - 1 - 8
117
118        cls.res1_m = np.array([260.6786029744658,  162.7966072512681,  105.8253545950566,  71.16133060790817,  47.68490211260372,  28.11843682526138,  13.03968537077271,   2.25398078597622,           137.9954,           106.7351,            79.5329,            56.2839,            37.0339,            21.7781,            10.4741,             2.9762,           143.6691,           111.7797,            83.9383,            60.0627,            40.1749,            24.2761,            12.3212,             4.1296,           154.7977,           121.7375,            92.7136,  67.63670000000001,            46.5716,            29.5147,             16.364,             6.9406])
119        cls.res2_m = np.array([97.88199572319769,  56.97125265621156,  34.66402398714837,  23.47642849530445,  19.56646528734234,  15.07875145448866,   10.7857045847965,   2.25398078597622,             45.893,            39.9085,            33.9271,             27.916,             21.837,            15.7175,             9.4748,             2.9762,            48.8795,            42.7679,            36.6301,            30.4428,            24.1592,            17.7961,            11.2246,             4.1296,            55.0335,            48.6606,            42.2333,            35.7359,            29.0609,            22.2519,            15.0923,             6.9406])
120
121
122class TestCointJoh25(CheckCointJoh):
123
124    @classmethod
125    def setup_class(cls):
126        with warnings.catch_warnings():
127            warnings.simplefilter("ignore", category=HypothesisTestWarning)
128            cls.res = coint_johansen(dta, 2, 5)
129        cls.nobs_r = 173 - 1 - 5
130
131        # Note: critical values not available if trend>1
132        cls.res1_m = np.array([270.1887263915158, 171.6870096307863,
133                               107.8613367358704, 70.82424032233558,
134                               44.62551818267534, 25.74352073857572,
135                               14.17882426926978, 4.288656185006764,
136                               0, 0, 0, 0, 0, 0, 0, 0, 0,
137                               0, 0, 0, 0, 0, 0, 0, 0, 0,
138                               0, 0, 0, 0, 0, 0])
139        cls.res1_m[cls.res1_m == 0] = np.nan
140        cls.res2_m = np.array([98.50171676072955, 63.82567289491584,
141                               37.03709641353485, 26.19872213966024,
142                               18.88199744409963, 11.56469646930594,
143                               9.890168084263012, 4.288656185006764,
144                               0, 0, 0, 0, 0, 0, 0, 0,
145                               0, 0, 0, 0, 0, 0, 0, 0,
146                               0, 0, 0, 0, 0, 0, 0, 0])
147        cls.res2_m[cls.res2_m == 0] = np.nan
148
149
150@pytest.mark.smoke
151def test_coint_johansen_0lag(reset_randomstate):
152    # GH 5731
153    x_diff = np.random.normal(0, 1, 1000)
154    x = pd.Series(np.cumsum(x_diff))
155    e1 = np.random.normal(0, 1, 1000)
156    y = x + 5 + e1
157    data = pd.concat([x, y], axis=1)
158    result = coint_johansen(data, det_order=-1, k_ar_diff=0)
159    assert result.eig.shape == (2,)
160