1import numpy as np
2
3import pytest
4from numpy.testing import assert_allclose, assert_raises
5
6from statsmodels.tsa.innovations.arma_innovations import arma_innovations
7from statsmodels.tsa.arima.datasets.brockwell_davis_2002 import dowj, lake
8from statsmodels.tsa.arima.estimators.durbin_levinson import durbin_levinson
9
10
11@pytest.mark.low_precision('Test against Example 5.1.1 in Brockwell and Davis'
12                           ' (2016)')
13def test_brockwell_davis_example_511():
14    # Note: this example is primarily tested in
15    # test_yule_walker::test_brockwell_davis_example_511.
16
17    # Difference the series
18    endog = dowj.diff().iloc[1:]
19
20    # Durbin-Levinson
21    dl, _ = durbin_levinson(endog, ar_order=2, demean=True)
22
23    assert_allclose(dl[0].params, np.var(endog))
24    assert_allclose(dl[1].params, [0.4219, 0.1479], atol=1e-4)
25    assert_allclose(dl[2].params, [0.3739, 0.1138, 0.1460], atol=1e-4)
26
27
28def check_itsmr(lake):
29    # Test against R itsmr::yw; see results/results_yw_dl.R
30    dl, _ = durbin_levinson(lake, 5)
31
32    assert_allclose(dl[0].params, np.var(lake))
33    assert_allclose(dl[1].ar_params, [0.8319112104])
34    assert_allclose(dl[2].ar_params, [1.0538248798, -0.2667516276])
35    desired = [1.0887037577, -0.4045435867, 0.1307541335]
36    assert_allclose(dl[3].ar_params, desired)
37    desired = [1.08425065810, -0.39076602696, 0.09367609911, 0.03405704644]
38    assert_allclose(dl[4].ar_params, desired)
39    desired = [1.08213598501, -0.39658257147, 0.11793957728, -0.03326633983,
40               0.06209208707]
41    assert_allclose(dl[5].ar_params, desired)
42
43    # itsmr::yw returns the innovations algorithm estimate of the variance
44    # we'll just check for p=5
45    u, v = arma_innovations(np.array(lake) - np.mean(lake),
46                            ar_params=dl[5].ar_params, sigma2=1)
47    desired_sigma2 = 0.4716322564
48    assert_allclose(np.sum(u**2 / v) / len(u), desired_sigma2)
49
50
51def test_itsmr():
52    # Note: apparently itsmr automatically demeans (there is no option to
53    # control this)
54    endog = lake.copy()
55
56    check_itsmr(endog)           # Pandas series
57    check_itsmr(endog.values)    # Numpy array
58    check_itsmr(endog.tolist())  # Python list
59
60
61def test_nonstationary_series():
62    # Test against R stats::ar.yw; see results/results_yw_dl.R
63    endog = np.arange(1, 12) * 1.0
64    res, _ = durbin_levinson(endog, 2, demean=False)
65
66    desired_ar_params = [0.92318534179, -0.06166314306]
67    assert_allclose(res[2].ar_params, desired_ar_params)
68
69
70@pytest.mark.xfail(reason='Different computation of variances')
71def test_nonstationary_series_variance():
72    # See `test_nonstationary_series`. This part of the test has been broken
73    # out as an xfail because we compute a different estimate of the variance
74    # from stats::ar.yw, but keeping the test in case we want to also implement
75    # that variance estimate in the future.
76    endog = np.arange(1, 12) * 1.0
77    res, _ = durbin_levinson(endog, 2, demean=False)
78
79    desired_sigma2 = 15.36526603
80    assert_allclose(res[2].sigma2, desired_sigma2)
81
82
83def test_invalid():
84    endog = np.arange(2) * 1.0
85    assert_raises(ValueError, durbin_levinson, endog, ar_order=2)
86    assert_raises(ValueError, durbin_levinson, endog, ar_order=-1)
87    assert_raises(ValueError, durbin_levinson, endog, ar_order=1.5)
88
89    endog = np.arange(10) * 1.0
90    assert_raises(ValueError, durbin_levinson, endog, ar_order=[1, 3])
91
92
93def test_misc():
94    # Test defaults (order = 0, demean=True)
95    endog = lake.copy()
96    res, _ = durbin_levinson(endog)
97    assert_allclose(res[0].params, np.var(endog))
98
99    # Test that integer input gives the same result as float-coerced input.
100    endog = np.array([1, 2, 5, 3, -2, 1, -3, 5, 2, 3, -1], dtype=int)
101    res_int, _ = durbin_levinson(endog, 2, demean=False)
102    res_float, _ = durbin_levinson(endog * 1.0, 2, demean=False)
103    assert_allclose(res_int[0].params, res_float[0].params)
104    assert_allclose(res_int[1].params, res_float[1].params)
105    assert_allclose(res_int[2].params, res_float[2].params)
106