add Baxter King filter

src/QuantEcon.jl

@@ -147,7 +147,8 @@ export
 
 # filter
     hp_filter,
-    hamilton_filter
+    hamilton_filter,
+    bk_filter
 
 
 include("sampler.jl")

src/filter.jl

@@ -85,3 +85,34 @@ function hamilton_filter(y::AbstractVector, h::Integer)
     y_trend = y - y_cycle
     return y_cycle, y_trend
 end
+
+@doc doc"""
+This function applies "Baxter King filter" to `y <:AbstractVector`.
+
+Reference: https://www.mitpressjournals.org/doi/abs/10.1162/003465399558454?casa_token=pfn7A97wi0QAAAAA:78AlfSOyiz__a6_snyH7jwydLu8uGZV-U3ZVY5Vo3dTi0r7da5uuVtUyUR87-uZBAuarQAYMiywJ76o
+
+##### Arguments
+- `y::AbstractVector` : data to be filtered
+- `wu::Real` : upper cutoff frequencies
+- `wl::Real` : lower cutoff frequencies
+- `K::Integer` : number of leads and lags of the moving average
+
+##### Returns
+- `y_cycle::Vector` : cyclical component
+- `y_trend::Vector` : trend component
+"""
+function bk_filter(y::AbstractVector, wu::Real, wl::Real, K::Integer)
+    T = length(y)
+    w1 = 2pi/wu
+    w2 = 2pi/wl
+    b = vcat((w2-w1)/pi, [(sin(w2*i)-sin(w1*i))/(i*pi) for i=1:K])
+    theta = (b[1] + 2*sum(b[2:end]))/(2K+1)
+    B = b .- theta
+    y_cycle = Vector{Union{Missing, Float64}}(undef, T)
+    for t = K+1:T-K
+        y_cycle[t] = y[t]*B[1] + 
+            dot(y[t-1:-1:t-K], B[2:end]) + dot(y[t+1:t+K], B[2:end])
+    end
+    y_trend = y - y_cycle
+    return y_cycle, y_trend
+end

test/test_filter.jl

@@ -20,4 +20,11 @@
         @test isapprox(data["ham_c"], data["ham_c_mat"], nans=true, rtol=1e-7, atol=1e-7)
         @test isapprox(data["ham_rw_c"], data["ham_rw_c_mat"], nans=true)
     end
+
+    @testset "test baxter king filter" begin
+        bk_cyc, _ = bk_filter([1,5,4,7,3,2,4,8], 24, 6, 2)
+        cyc_py = [missing, missing, 0.5401001384504593, 0.34785164507573785,
+                  -0.05672021307788492, -0.887951783526198, missing, missing]
+        @test all((ismissing.(bk_cyc) .& ismissing.(cyc_py)).| isapprox.(bk_cyc, cyc_py))
+    end
 end