Hide keyboard shortcuts

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

1 

2import numpy as np 

3from scipy.signal import fftconvolve 

4 

5from statsmodels.tools.validation import array_like, PandasWrapper 

6 

7 

8def bkfilter(x, low=6, high=32, K=12): 

9 """ 

10 Filter a time series using the Baxter-King bandpass filter. 

11 

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. 

26 

27 Returns 

28 ------- 

29 ndarray 

30 The cyclical component of x. 

31 

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. 

42 

43 Notes 

44 ----- 

45 Returns a centered weighted moving average of the original series. Where 

46 the weights a[j] are computed :: 

47 

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,... 

51 

52 and theta is a normalizing constant :: 

53 

54 theta = -sum(b)/(2K+1) 

55 

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. 

61 

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) 

69 

70 >>> cycles = sm.tsa.filters.bkfilter(dta[['realinv']], 6, 24, 12) 

71 

72 >>> import matplotlib.pyplot as plt 

73 >>> fig, ax = plt.subplots() 

74 >>> cycles.plot(ax=ax, style=['r--', 'b-']) 

75 >>> plt.show() 

76 

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 

98 

99 return pw.wrap(x, append='cycle', trim_start=K, trim_end=K)