Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/statsmodels/tsa/filters/bk_filter.py : 21%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
2import numpy as np
3from scipy.signal import fftconvolve
5from statsmodels.tools.validation import array_like, PandasWrapper
8def bkfilter(x, low=6, high=32, K=12):
9 """
10 Filter a time series using the Baxter-King bandpass filter.
12 Parameters
13 ----------
14 x : array_like
15 A 1 or 2d ndarray. If 2d, variables are assumed to be in columns.
16 low : float
17 Minimum period for oscillations, ie., Baxter and King suggest that
18 the Burns-Mitchell U.S. business cycle has 6 for quarterly data and
19 1.5 for annual data.
20 high : float
21 Maximum period for oscillations BK suggest that the U.S.
22 business cycle has 32 for quarterly data and 8 for annual data.
23 K : int
24 Lead-lag length of the filter. Baxter and King propose a truncation
25 length of 12 for quarterly data and 3 for annual data.
27 Returns
28 -------
29 ndarray
30 The cyclical component of x.
32 See Also
33 --------
34 statsmodels.tsa.filters.cf_filter.cffilter
35 The Christiano Fitzgerald asymmetric, random walk filter.
36 statsmodels.tsa.filters.bk_filter.hpfilter
37 Hodrick-Prescott filter.
38 statsmodels.tsa.seasonal.seasonal_decompose
39 Decompose a time series using moving averages.
40 statsmodels.tsa.seasonal.STL
41 Season-Trend decomposition using LOESS.
43 Notes
44 -----
45 Returns a centered weighted moving average of the original series. Where
46 the weights a[j] are computed ::
48 a[j] = b[j] + theta, for j = 0, +/-1, +/-2, ... +/- K
49 b[0] = (omega_2 - omega_1)/pi
50 b[j] = 1/(pi*j)(sin(omega_2*j)-sin(omega_1*j), for j = +/-1, +/-2,...
52 and theta is a normalizing constant ::
54 theta = -sum(b)/(2K+1)
56 References
57 ----------
58 Baxter, M. and R. G. King. "Measuring Business Cycles: Approximate
59 Band-Pass Filters for Economic Time Series." *Review of Economics and
60 Statistics*, 1999, 81(4), 575-593.
62 Examples
63 --------
64 >>> import statsmodels.api as sm
65 >>> import pandas as pd
66 >>> dta = sm.datasets.macrodata.load_pandas().data
67 >>> index = pd.DatetimeIndex(start='1959Q1', end='2009Q4', freq='Q')
68 >>> dta.set_index(index, inplace=True)
70 >>> cycles = sm.tsa.filters.bkfilter(dta[['realinv']], 6, 24, 12)
72 >>> import matplotlib.pyplot as plt
73 >>> fig, ax = plt.subplots()
74 >>> cycles.plot(ax=ax, style=['r--', 'b-'])
75 >>> plt.show()
77 .. plot:: plots/bkf_plot.py
78 """
79 # TODO: change the docstring to ..math::?
80 # TODO: allow windowing functions to correct for Gibb's Phenomenon?
81 # adjust bweights (symmetrically) by below before demeaning
82 # Lancosz Sigma Factors np.sinc(2*j/(2.*K+1))
83 pw = PandasWrapper(x)
84 x = array_like(x, 'x', maxdim=2)
85 omega_1 = 2. * np.pi / high # convert from freq. to periodicity
86 omega_2 = 2. * np.pi / low
87 bweights = np.zeros(2 * K + 1)
88 bweights[K] = (omega_2 - omega_1) / np.pi # weight at zero freq.
89 j = np.arange(1, int(K) + 1)
90 weights = 1 / (np.pi * j) * (np.sin(omega_2 * j) - np.sin(omega_1 * j))
91 bweights[K + j] = weights # j is an idx
92 bweights[:K] = weights[::-1] # make symmetric weights
93 bweights -= bweights.mean() # make sure weights sum to zero
94 if x.ndim == 2:
95 bweights = bweights[:, None]
96 x = fftconvolve(x, bweights, mode='valid')
97 # get a centered moving avg/convolution
99 return pw.wrap(x, append='cycle', trim_start=K, trim_end=K)