Add Hodrick-Prescott filter to current smoothers

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.5"
 
 [deps]
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
 Dierckx = "0.5"

README.md

@@ -7,6 +7,7 @@ The package Smoothers provides a collection of smoothing heuristics, models and
 * Locally Estimated Scatterplot Smoothing (**loess**)
 * Seasonal and Trend decomposition based on Loess (**stl**) 
 * Simple Moving Average (**sma**)
+* Hodrick-Prescott filter (**hpfilter**) for trend and cyclical components. 
 
 <img src="./docs/src/images/smoothers.png">
 
@@ -37,6 +38,9 @@ plot!(t,filter(b,a,x), label ="filter(1,[1/"*sw*",...],x)")
 # Simple Moving Average
 plot!(t, sma(x,w,true), label = "sma(x,"*sw*",true)")
 
+# Hodrick-Prescott filter
+plot!(t, hpfilter(x, 1600), label = "hpfilter(x, 1600)")
+
 ```
 <img src="./docs/src/images/smoothers_examples.png">
 

src/Smoothers.jl

@@ -2,15 +2,17 @@ module Smoothers
 
 using Dierckx
 import Base: filter
+using LinearAlgebra: diagm, I
 
-export filter, hma, loess, sma, stl
+export filter, hma, loess, sma, stl, hpfilter
 
 # Methods
 include("filter.jl")
 include("hma.jl")
 include("loess.jl")
 include("sma.jl")
 include("stl.jl")
+include("hpfilter.jl")
 
 """
 Collection of smoothing heuristics, models and smoothing related applications. The current available smoothers and applications are:
@@ -20,6 +22,7 @@ Collection of smoothing heuristics, models and smoothing related applications. T
     `loess`:  Locally Estimated Scatterplot Smoothing
     `sma`:    Simple Moving Average
     `stl`:    Seasonal and Trend decomposition based on Loess
+    `hpfilter`: Hodrick-Prescott trend decomposition. 
 
 """
 Smoothers

src/hpfilter.jl

@@ -0,0 +1,91 @@
+"""
+Package: Smoothers
+
+    hpfilter(y, lambda=1600)
+
+Hodrick-Prescott filter vector `y` using smoothing parameter `lambda`. The
+parameter `lambda` adjusts the sensibility of the trend to short-term
+fluctuations. Higher values of `lambda` give more smoothing. Hodrick and
+Prescott suggest `lambda=1600` for quarterly data. 
+- `lambda=0` gives no smoothing, so it returns `y`. 
+- `lambda`→∞ gives a linear trend.
+
+# Arguments
+- `y`: Vector of data.
+- `lambda`: Smoothing parameter.
+
+# Returns
+Vector with the trend component of `y`. 
+
+# Examples
+```julia-repl
+julia> x = 0:0.2:2pi
+0.0:0.2:6.2
+
+julia> y = sin.(x) + 0.1*randn(length(x))
+32-element Vector{Float64}:
+  0.09200627969509063
+  0.24164024275732685
+  0.49994368543454165
+  0.5122933360879012
+  0.6699190280838954
+  0.9630579533383277
+  ⋮
+ -0.7908099824097626
+ -0.5986641704672869
+ -0.6482750454579537
+ -0.2963091814085014
+  0.050926679466235095
+
+julia> ytrend = hpfilter(y)
+32-element Vector{Float64}:
+  0.8054238538000016
+  0.7828566743521803
+  0.7598436089205433
+  0.7356005112515279
+  0.7091807976393933
+  0.6794983173939206
+  ⋮
+ -0.7006212315020678
+ -0.7350501835373815
+ -0.7671301582815108
+ -0.7977981479965659
+ -0.8279168604991427
+```
+"""
+function hpfilter(y, lambda=1600)
+
+    lambda == 0 && (return y)
+    T = length(y)
+    R = eltype(y)
+
+    # Center part of the matrix
+    A_c = 6*I(T-4) + 
+        diagm(1 => repeat([-4], T-5)) + 
+        diagm(-1 => repeat([-4], T-5)) + 
+        diagm(2 => repeat([1], T-6)) + 
+        diagm(-2 => repeat([1], T-6)) 
+
+    # First and last rows
+    r1 = [1,-2, 1, zeros(R, T-3)...] 
+    r2 = [-2, 5, -4, 1, zeros(R, T-4)...]
+
+    A = zeros(R, T, T) 
+    # first two rows/cols
+    A[1, :] = r1 
+    A[:, 1] = r1' 
+    A[2, :] = r2 
+    A[:, 2] = r2' 
+
+    # center    
+    A[3:end-2, 3:end-2] = A_c 
+    # last two rows/cols
+    A[end-1, end:-1:1] = r2
+    A[end:-1:1, end-1] = r2' 
+    A[end, end:-1:1]   = r1 
+    A[end:-1:1, end]   = r1' 
+
+    # Smoothed version 
+    y_s = (I(T) + lambda*A)\y
+    y_s
+end

test/hpfilter.jl

@@ -0,0 +1,24 @@
+using Test
+using Smoothers
+
+@testset "hpfilter" begin
+
+    # Input test data
+    x = [1579.1, 1577.7, 1576.9, 1600.4, 1626.4, 1655.5, 1665.1, 1669.0, 1643.8, 1638.6]
+
+    # Computed HP filter from another source with lambda=1600
+    x_hp_trend = [1577.18535023776, 1587.47395687673, 1597.7637601718, 1608.04984805602, 1618.31426861234, 1628.53428876867, 1638.69222903504, 1648.78726349099, 1658.83507107293, 1668.86396367757]
+    # Compute HP filter with package function 
+    x_trend = hpfilter(x, 1600)
+    # Are entries approximate?
+    @test all(isapprox.(x_trend, x_hp_trend; atol=1e-8))
+
+    # Test with different data precision 
+    x_f32 = Float32.(x)
+    x_trend_f32 = hpfilter(x_f32, 1600)
+    @test eltype(x_trend_f32) == Float32
+
+    # Test with lambda=0
+    @test x == hpfilter(x, 0)
+
+end

test/runtests.jl

@@ -9,6 +9,7 @@ const tests = [
     "loess",
     "sma",
     "stl",
+    "hpfilter",
 ]
 
 printstyled("\nTest Summary List:\n", color=:underline)