File: test_multivariate_switch_univariate.py

package info (click to toggle)
statsmodels 0.12.2-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 46,076 kB
  • sloc: python: 232,741; f90: 612; sh: 389; javascript: 337; makefile: 164; asm: 156; ansic: 16; xml: 9
file content (491 lines) | stat: -rw-r--r-- 19,226 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
"""
Tests for automatic switching of the filter method from multivariate to
univariate when the forecast error covariance matrix is singular.

Author: Chad Fulton
License: Simplified-BSD

References
----------

Kim, Chang-Jin, and Charles R. Nelson. 1999.
"State-Space Models with Regime Switching:
Classical and Gibbs-Sampling Approaches with Applications".
MIT Press Books. The MIT Press.

Hamilton, James D. 1994.
Time Series Analysis.
Princeton, N.J.: Princeton University Press.
"""
import numpy as np
import pytest
import sys
import platform
import re
i386_looser_tolerances=bool(re.match('i.?86|x86',platform.uname()[4])) and np.log2(sys.maxsize)<33

from statsmodels.tsa.statespace import (
    mlemodel, sarimax, structural, varmax, dynamic_factor)
from statsmodels.tsa.statespace.tests.test_impulse_responses import TVSS
from numpy.testing import assert_allclose


def get_model(univariate, missing=None, init=None):
    if univariate:
        endog = np.array([0.5, 1.2, -0.2, 0.3, -0.1, 0.4, 1.4, 0.9])

        if missing == 'init':
            endog[0:2] = np.nan
        elif missing == 'mixed':
            endog[2:4] = np.nan
        elif missing == 'all':
            endog[:] = np.nan

        mod = mlemodel.MLEModel(endog, k_states=1, k_posdef=1)
        mod['design', 0, 0] = 1.
        mod['transition', 0, 0] = 0.5
        mod['selection', 0, 0] = 1.
        mod['state_cov', 0, 0] = 1.
        mod['state_intercept', 0, 0] = 1.
    else:
        endog = np.array([[0.5, 1.2, -0.2, 0.3, -0.1, 0.4, 1.4, 0.9],
                          [-0.2, -0.3, -0.1, 0.1, 0.01, 0.05, -0.13, -0.2]]).T

        if missing == 'init':
            endog[0:2, :] = np.nan
        elif missing == 'mixed':
            endog[2:4, 0] = np.nan
            endog[3:6, 1] = np.nan
        elif missing == 'all':
            endog[:] = np.nan

        mod = mlemodel.MLEModel(endog, k_states=3, k_posdef=2)
        mod['obs_intercept'] = np.array([0.5, 0.2])
        mod['design'] = np.array([[0.1, -0.1, 0],
                                  [0.2, 0.3, 0]])
        mod['obs_cov'] = np.array([[5, -0.2],
                                   [-0.2, 3.]])

        mod['transition', 0, 0] = 1
        mod['transition', 1:, 1:] = np.array([[0.5, -0.1],
                                              [1., 0.]])
        mod['selection', :2, :2] = np.eye(2)
        mod['state_cov'] = np.array([[1.2, 0.2],
                                     [0.2, 2.5]])
        mod['state_intercept', :2] = np.array([1., -1.])

    if init == 'diffuse':
        mod.ssm.initialize_diffuse()
    elif init == 'approximate_diffuse':
        mod.ssm.initialize_approximate_diffuse()
    elif init == 'stationary':
        mod.ssm.initialize_stationary()

    return mod


def check_filter_output(mod, periods, atol=0):
    if isinstance(mod, mlemodel.MLEModel):
        # Multivariate filter
        res_mv = mod.ssm.filter()

        # Manually perform filtering with a switch
        mod.ssm.filter()
        kfilter = mod.ssm._kalman_filter
        kfilter.seek(0, True)
        kfilter.univariate_filter[periods] = 1
        for _ in range(mod.nobs):
            next(kfilter)
        # Create the results object
        res_switch = mod.ssm.results_class(mod.ssm)
        res_switch.update_representation(mod.ssm)
        res_switch.update_filter(kfilter)

        # Univariate filter
        mod.ssm.filter_univariate = True
        res_uv = mod.ssm.filter()
    else:
        res_mv, res_switch, res_uv = mod

    # Test attributes that are the same regardless of the univariate or
    # multivariate method
    assert_allclose(res_switch.llf, res_mv.llf)
    assert_allclose(res_switch.llf, res_uv.llf)
    assert_allclose(res_switch.scale, res_mv.scale)
    assert_allclose(res_switch.scale, res_uv.scale)

    attrs = ['forecasts_error_diffuse_cov', 'predicted_state',
             'predicted_state_cov', 'predicted_diffuse_state_cov',
             'filtered_state', 'filtered_state_cov', 'llf_obs']
    for attr in attrs:
        attr_mv = getattr(res_mv, attr)
        attr_uv = getattr(res_uv, attr)
        attr_switch = getattr(res_switch, attr)
        if attr_mv is None:
            continue
        assert_allclose(attr_switch, attr_mv, atol=atol)
        assert_allclose(attr_switch, attr_uv, atol=atol)

    # Test attributes that can differ for the univariate vs multivariate method
    attrs = ['forecasts_error', 'forecasts_error_cov', 'kalman_gain']
    for attr in attrs:
        # Test all periods against the multivariate filter, except for periods
        # that were switched (it's easiest to just set those values to zero)
        actual = getattr(res_switch, attr).copy()
        desired = getattr(res_mv, attr).copy()
        actual[..., periods] = 0
        desired[..., periods] = 0
        assert_allclose(actual, desired, atol=atol)

        # Test switched periods against the univariate filter
        actual = getattr(res_switch, attr)[..., periods]
        desired = getattr(res_uv, attr)[..., periods]
        assert_allclose(actual, desired, atol=atol)


def check_smoother_output(mod, periods, atol=1e-12):
    if isinstance(mod, mlemodel.MLEModel):
        # Multivariate filter / smoother
        res_mv = mod.ssm.smooth()

        # Manually perform filtering / smoothing with a switch
        kfilter = mod.ssm._kalman_filter
        kfilter.seek(0, True)
        kfilter.univariate_filter[periods] = 1
        for _ in range(mod.nobs):
            next(kfilter)
        # Create the results object
        res_switch = mod.ssm.results_class(mod.ssm)
        res_switch.update_representation(mod.ssm)
        res_switch.update_filter(kfilter)
        mod.ssm._kalman_smoother.reset(True)
        smoother = mod.ssm._smooth()
        res_switch.update_smoother(smoother)

        # Univariate filter / smoother
        mod.ssm.filter_univariate = True
        res_uv = mod.ssm.smooth()
    else:
        res_mv, res_switch, res_uv = mod

    # Test attributes that are the same regardless of the univariate or
    # multivariate method
    attrs = ['scaled_smoothed_estimator', 'scaled_smoothed_estimator_cov',
             'smoothed_state', 'smoothed_state_cov', 'smoothed_state_autocov',
             'smoothed_state_disturbance', 'smoothed_state_disturbance_cov',
             'innovations_transition']
    for attr in attrs:
        attr_mv = getattr(res_mv, attr)
        attr_uv = getattr(res_uv, attr)
        attr_switch = getattr(res_switch, attr)
        if attr_mv is None:
            continue
        assert_allclose(attr_uv, attr_mv, atol=atol)
        assert_allclose(attr_switch, attr_mv, atol=atol)
        assert_allclose(attr_switch, attr_uv, atol=atol)

    # Test attributes that can differ for the univariate vs multivariate method
    attrs = ['smoothing_error', 'smoothed_measurement_disturbance',
             'smoothed_measurement_disturbance_cov']
    for attr in attrs:
        attr_mv = getattr(res_mv, attr)
        attr_uv = getattr(res_uv, attr)
        attr_switch = getattr(res_switch, attr)
        if attr_mv is None:
            continue

        # Test all periods against the multivariate filter, except for periods
        # that were switched (it's easiest to just set those values to zero)
        actual = attr_switch.copy()
        desired = attr_mv.copy()
        actual[..., periods] = 0
        desired[..., periods] = 0
        assert_allclose(actual, desired)


@pytest.mark.parametrize('missing', [None, 'init', 'mixed', 'all'])
def test_basic(missing):
    # Test that the multivariate filter switches to the univariate filter
    # when it runs into problems
    mod = get_model(univariate=True, missing=missing)

    # Here, because of the known initialization with P_0 = [[0]], we will also
    # have F_0 = 0.
    # Then the Kalman filter gives P_0|0 = 0, and P_1 = Q = [[1.]]
    # so that F_1 != 0, and the rest of the periods do not have a singular
    # forecast error covariance matrix.
    mod.initialize_known([0], [[0]])
    mod.ssm.filter()
    uf = np.array(mod.ssm._kalman_filter.univariate_filter)

    # As a result, we expect that in the period t=0, we had to fall back to the
    # univariate filter, while in the periods t >= 1, the multivariate filter
    # works as usual.
    # However, if the first period is missing (as in init and all), then we
    # essentially skip the forecast error and forecast error cov computation.
    # As a result, we don't need to switch to the univariate methods
    if missing in ['init', 'all']:
        assert_allclose(uf, 0)
    else:
        assert_allclose(uf[0], 1)
        assert_allclose(uf[1:], 0)


@pytest.mark.parametrize('univariate', [True, False])
@pytest.mark.parametrize('missing', [None, 'init', 'mixed', 'all'])
@pytest.mark.parametrize(
    'init', ['stationary', 'diffuse', 'approximate_diffuse'])
@pytest.mark.parametrize('periods', [np.s_[0], np.s_[4:6], np.s_[:]])
def test_filter_output(univariate, missing, init, periods):
    # Test the output when the multivariate filter switches to the univariate
    # filter
    mod = get_model(univariate, missing, init)
    check_filter_output(mod, periods, atol=1e-10 if i386_looser_tolerances else 0)


@pytest.mark.parametrize('univariate', [True, False])
@pytest.mark.parametrize('missing', [None, 'init', 'mixed', 'all'])
@pytest.mark.parametrize('init',
                         ['stationary', 'diffuse', 'approximate_diffuse'])
@pytest.mark.parametrize('periods', [np.s_[0], np.s_[4:6], np.s_[:]])
@pytest.mark.parametrize('option', [None, 'alternate_timing'])
def test_smoother_output(univariate, missing, init, periods, option):
    # Test the output when the multivariate filter switches to the univariate
    # filter

    mod = get_model(univariate, missing, init)
    if option == 'alternate_timing':
        # Can't use diffuse initialization with alternate timing
        if init == 'diffuse':
            return
        mod.ssm.timing_init_filtered = True
    atol = 1e-8 if i386_looser_tolerances else 1e-12
    # Tolerance is lower for approximate diffuse for one attribute in this case
    if missing == 'init' and init == 'approximate_diffuse':
        atol = 1e-6
    check_smoother_output(mod, periods, atol=atol)


def test_invalid_options():
    mod = get_model(univariate=True)
    mod.initialize_known([0], [[0]])

    mod.ssm.set_inversion_method(0, solve_lu=True)
    msg = ('Singular forecast error covariance matrix detected, but'
           ' multivariate filter cannot fall back to univariate'
           ' filter when the inversion method is set to anything'
           ' other than INVERT_UNIVARIATE or SOLVE_CHOLESKY.')
    with pytest.raises(NotImplementedError, match=msg):
        mod.ssm.filter()

    mod = get_model(univariate=True)
    mod.initialize_known([0], [[0]])
    mod.ssm.smooth_classical = True
    msg = ('Cannot use classical smoothing when the multivariate filter has'
           ' fallen back to univariate filtering.')
    with pytest.raises(NotImplementedError, match=msg):
        mod.ssm.smooth()

    mod = get_model(univariate=True)
    mod.initialize_known([0], [[0]])
    mod.ssm.smooth_alternative = True
    msg = ('Cannot use alternative smoothing when the multivariate filter has'
           ' fallen back to univariate filtering.')
    with pytest.raises(NotImplementedError, match=msg):
        mod.ssm.smooth()


@pytest.mark.parametrize('missing', [None, 'init', 'mixed', 'all'])
@pytest.mark.parametrize('periods', [np.s_[0], np.s_[4:6], np.s_[:]])
@pytest.mark.parametrize('use_exact_diffuse', [False, True])
def test_sarimax(missing, periods, use_exact_diffuse):
    endog = np.array([0.5, 1.2, -0.2, 0.3, -0.1, 0.4, 1.4, 0.9])
    exog = np.ones_like(endog)
    if missing == 'init':
        endog[0:2] = np.nan
    elif missing == 'mixed':
        endog[2:4] = np.nan
    elif missing == 'all':
        endog[:] = np.nan

    mod = sarimax.SARIMAX(endog, order=(1, 1, 1), trend='t',
                          seasonal_order=(1, 1, 1, 2), exog=exog,
                          use_exact_diffuse=use_exact_diffuse)
    mod.update([0.1, 0.3, 0.5, 0.2, 0.05, -0.1, 1.0])
    check_filter_output(mod, periods, atol=1e-8)
    check_smoother_output(mod, periods)


@pytest.mark.parametrize('missing', [None, 'init', 'mixed', 'all'])
@pytest.mark.parametrize('periods', [np.s_[0], np.s_[4:6], np.s_[:]])
@pytest.mark.parametrize('use_exact_diffuse', [False, True])
def test_unobserved_components(missing, periods, use_exact_diffuse):
    endog = np.array([0.5, 1.2, -0.2, 0.3, -0.1, 0.4, 1.4, 0.9])
    exog = np.ones_like(endog)
    if missing == 'init':
        endog[0:2] = np.nan
    elif missing == 'mixed':
        endog[2:4] = np.nan
    elif missing == 'all':
        endog[:] = np.nan

    mod = structural.UnobservedComponents(
        endog, 'llevel', exog=exog, seasonal=2, autoregressive=1,
        use_exact_diffuse=use_exact_diffuse)
    mod.update([1.0, 0.1, 0.3, 0.05, 0.15, 0.5])
    check_filter_output(mod, periods)
    check_smoother_output(mod, periods)


@pytest.mark.parametrize('missing', [None, 'init', 'mixed', 'all'])
@pytest.mark.parametrize('periods', [np.s_[0], np.s_[4:6], np.s_[:]])
def test_varmax(missing, periods):
    endog = np.array([[0.5, 1.2, -0.2, 0.3, -0.1, 0.4, 1.4, 0.9],
                      [-0.2, -0.3, -0.1, 0.1, 0.01, 0.05, -0.13, -0.2]]).T
    exog = np.ones_like(endog[:, 0])
    if missing == 'init':
        endog[0:2, :] = np.nan
    elif missing == 'mixed':
        endog[2:4, 0] = np.nan
        endog[3:6, 1] = np.nan
    elif missing == 'all':
        endog[:] = np.nan

    mod = varmax.VARMAX(endog, order=(1, 0), trend='t', exog=exog)
    mod.update([0.1, -0.1, 0.5, 0.1, -0.05, 0.2, 0.4, 0.25, 1.2, 0.4, 2.3])
    check_filter_output(mod, periods, atol=1e-12)
    check_smoother_output(mod, periods)


@pytest.mark.parametrize('missing', [None, 'init', 'mixed', 'all'])
@pytest.mark.parametrize('periods', [np.s_[0], np.s_[4:6], np.s_[:]])
def test_dynamic_factor(missing, periods):
    endog = np.array([[0.5, 1.2, -0.2, 0.3, -0.1, 0.4, 1.4, 0.9],
                      [-0.2, -0.3, -0.1, 0.1, 0.01, 0.05, -0.13, -0.2]]).T
    exog = np.ones_like(endog[:, 0])
    if missing == 'init':
        endog[0:2, :] = np.nan
    elif missing == 'mixed':
        endog[2:4, 0] = np.nan
        endog[3:6, 1] = np.nan
    elif missing == 'all':
        endog[:] = np.nan

    mod = dynamic_factor.DynamicFactor(endog, k_factors=1, factor_order=2,
                                       exog=exog)
    mod.update([1.0, -0.5, 0.3, -0.1, 1.2, 2.3, 0.5, 0.1])
    check_filter_output(mod, periods)
    check_smoother_output(mod, periods)


@pytest.mark.parametrize('missing', [None, 'mixed'])
def test_simulation_smoothing(missing):
    # Test that the simulation smoother works when the multivariate filter
    # switches to the univariate filter when it runs into problems
    # (see test_basic for a description of the model used here)

    # Get the model where switching will occur
    mod_switch = get_model(univariate=True, missing=missing)
    mod_switch.initialize_known([0], [[0]])
    sim_switch = mod_switch.simulation_smoother()

    # Get the model where we have specified univariate filtering (so there is
    # no need to switch)
    mod_uv = get_model(univariate=True, missing=missing)
    mod_uv.initialize_known([0], [[0]])
    mod_uv.ssm.filter_univariate = True
    sim_uv = mod_uv.simulation_smoother()

    # Test for basic simulationg of a new observed series
    np.random.seed(1234)
    simulate_switch = mod_switch.simulate([], 10)
    np.random.seed(1234)
    simulate_uv = mod_uv.simulate([], 10)
    assert_allclose(simulate_switch, simulate_uv)

    # Perform simulation smoothing
    np.random.seed(1234)
    sim_switch.simulate()
    np.random.seed(1234)
    sim_uv.simulate()

    # Make sure that switching happened in the first model but not the second
    kfilter = sim_switch._simulation_smoother.simulated_kfilter
    uf_switch = np.array(kfilter.univariate_filter, copy=True)
    assert_allclose(uf_switch[0], 1)
    assert_allclose(uf_switch[1:], 0)
    kfilter = sim_uv._simulation_smoother.simulated_kfilter.univariate_filter
    uf_uv = np.array(kfilter, copy=True)
    assert_allclose(uf_uv, 1)
    if missing == 'mixed':
        kfilter = (sim_switch._simulation_smoother
                             .secondary_simulated_kfilter.univariate_filter)
        uf_switch = np.array(kfilter, copy=True)
        assert_allclose(uf_switch[0], 1)
        assert_allclose(uf_switch[1:], 0)
        kfilter = (sim_uv._simulation_smoother
                         .secondary_simulated_kfilter.univariate_filter)
        uf_uv = np.array(kfilter, copy=True)
        assert_allclose(uf_uv, 1)

    # Test all simulation smoothing output
    attrs = ['generated_measurement_disturbance',
             'generated_state_disturbance', 'generated_obs', 'generated_state',
             'simulated_state', 'simulated_measurement_disturbance',
             'simulated_state_disturbance']
    for attr in attrs:
        assert_allclose(getattr(sim_switch, attr), getattr(sim_uv, attr))


def test_time_varying_model(reset_randomstate):
    endog = np.array([[0.5, 1.2, -0.2, 0.3, -0.1, 0.4, 1.4, 0.9],
                      [-0.2, -0.3, -0.1, 0.1, 0.01, 0.05, -0.13, -0.2]]).T

    # The basic model switches to the univariate method at observation 3,
    # because the forecast error covariance matrix will have a singular
    # component corresponding to the first endog variable
    np.random.seed(1234)
    mod_switch = TVSS(endog)
    mod_switch['design', ..., 3] = 0
    mod_switch['obs_cov', ..., 3] = 0
    mod_switch['obs_cov', 1, 1, 3] = 1.
    res_switch = mod_switch.ssm.smooth()
    kfilter = mod_switch.ssm._kalman_filter
    uf_switch = np.array(kfilter.univariate_filter, copy=True)

    # Next, this model only uses the univariate method
    np.random.seed(1234)
    mod_uv = TVSS(endog)
    mod_uv['design', ..., 3] = 0
    mod_uv['obs_cov', ..., 3] = 0
    mod_uv['obs_cov', 1, 1, 3] = 1.
    mod_uv.ssm.filter_univariate = True
    res_uv = mod_uv.ssm.smooth()
    kfilter = mod_uv.ssm._kalman_filter
    uf_uv = np.array(kfilter.univariate_filter, copy=True)

    # Finally, this model uses the multivariate method and gets around the
    # issue by setting the endog variable to NaN that would have contributed
    # to the singular part of the forecast error covariance matrix
    np.random.seed(1234)
    endog_mv = endog.copy()
    endog_mv[3, 0] = np.nan
    mod_mv = TVSS(endog_mv)
    mod_mv['design', ..., 3] = 0
    mod_mv['obs_cov', ..., 3] = 0
    mod_mv['obs_cov', 1, 1, 3] = 1.
    res_mv = mod_mv.ssm.smooth()
    kfilter = mod_mv.ssm._kalman_filter
    uf_mv = np.array(kfilter.univariate_filter, copy=True)

    # Make sure that switching happened in the switch model but not in the
    # other two models
    assert_allclose(uf_switch[:3], 0)
    assert_allclose(uf_switch[3], 1)
    assert_allclose(uf_switch[4:], 0)
    assert_allclose(uf_uv, 1)
    assert_allclose(uf_mv, 0)

    # Check filter and smoother output
    check_filter_output([res_mv, res_switch, res_uv], np.s_[3])
    check_smoother_output([res_mv, res_switch, res_uv], np.s_[3])