Add HP and Hamilton filter (#205)

* add hp and hamilton filter

* maybe fixing undefined type error

* does it pass the test?

* maybe this does the trick

* minor change of using CSV

* I don't want to mess my PR though...

* hopefully this works

* be more specific

* forgot to fix docs

* simplify the code and docs

* no longer need it

* I was so dumb to miss it in the pervious commit

* need this instead

* use dict instead of dataframe

src/QuantEcon.jl

@@ -137,7 +137,11 @@ export
     MVNSampler,
 
 # modeltools
-    @def_sim
+    @def_sim,
+
+# filter
+    hp_filter,
+    hamilton_filter
 
 
 include("sampler.jl")
@@ -152,6 +156,7 @@ include("markov/random_mc.jl")
 include("discrete_rv.jl")
 include("ecdf.jl")
 include("estspec.jl")
+include("filter.jl")
 include("kalman.jl")
 include("lae.jl")
 include("lqcontrol.jl")

src/filter.jl

@@ -0,0 +1,87 @@
+doc"""
+apply Hodrick-Prescott filter to `AbstractVector`.
+
+##### Arguments
+- `y::AbstractVector` : data to be detrended
+- `λ::Real` : penalty on variation in trend
+
+##### Returns
+- `y_cyclical::Vector`: cyclical component
+- `y_trend::Vector`: trend component
+"""
+function hp_filter(y::AbstractVector{T}, λ::Real) where T <: Real
+    y = Vector(y)
+    N = length(y)
+    H = spdiagm(-2 => fill(λ, N-2),
+                -1 => vcat(-2λ, fill(-4λ, N - 3), -2λ),
+                 0 => vcat(1 + λ, 1 + 5λ, fill(1 + 6λ, N-4),
+                           1 + 5λ, 1 + λ),
+                 1 => vcat(-2λ, fill(-4λ, N - 3), -2λ),
+                 2 => fill(λ, N-2))
+    y_trend = H \ y
+    y_cyclical = y - y_trend
+    return y_cyclical, y_trend
+end
+
+doc"""
+This function applies "Hamilton filter" to `AbstractVector`.
+
+http://econweb.ucsd.edu/~jhamilto/hp.pdf
+
+##### Arguments
+- `y::AbstractVector` : data to be filtered
+- `h::Integer` : Time horizon that we are likely to predict incorrectly.
+                 Original paper recommends 2 for annual data, 8 for quarterly data,
+                 24 for monthly data.
+- `p::Integer` : Number of lags in regression. Must be greater than `h`.
+Note: For seasonal data, it's desirable for `p` and `h` to be integer multiples
+      of the number of obsevations in a year.
+      e.g. For quarterly data, `h = 8` and `p = 4` are recommended.
+##### Returns
+- `y_cycle::Vector` : cyclical component
+- `y_trend::Vector` : trend component
+"""
+function hamilton_filter(y::AbstractVector, h::Integer, p::Integer)
+    y = Vector(y)
+    T = length(y)
+    y_cycle = fill(NaN, T)
+
+    # construct X matrix of lags
+    X = ones(T-p-h+1)
+    for j = 1:p
+        X = hcat(X, y[p-j+1:T-h-j+1])
+    end
+
+    # do OLS regression
+    b = (X'*X)\(X'*y[p+h:T])
+    y_cycle[p+h:T] = y[p+h:T] - X*b
+    y_trend = vcat(fill(NaN, p+h-1), X*b)
+    return y_cycle, y_trend
+end
+
+doc"""
+This function applies "Hamilton filter" to `<:AbstractVector`
+under random walk assumption.
+
+http://econweb.ucsd.edu/~jhamilto/hp.pdf
+
+##### Arguments
+- `y::AbstractVector` : data to be filtered
+- `h::Integer` : Time horizon that we are likely to predict incorrectly.
+                 Original paper recommends 2 for annual data, 8 for quarterly data,
+                 24 for monthly data.
+Note: For seasonal data, it's desirable for `h` to be an integer multiple
+      of the number of obsevations in a year.
+      e.g. For quarterly data, `h = 8` is recommended.
+##### Returns
+- `y_cycle::Vector` : cyclical component
+- `y_trend::Vector` : trend component
+"""
+function hamilton_filter(y::AbstractVector, h::Integer)
+    y = Vector(y)
+    T = length(y)
+    y_cycle = fill(NaN, T)
+    y_cycle[h+1:T] = y[h+1:T] - y[1:T-h]
+    y_trend = y - y_cycle
+    return y_cycle, y_trend
+end

test/REQUIRE

@@ -2,4 +2,4 @@ QuantEcon
 HDF5
 MAT
 JLD
-DataStructures
+DataStructures

test/runtests.jl

@@ -11,6 +11,7 @@ tests = [
         "discrete_rv",
         "ecdf",
         "estspec",
+        "filter",
         "kalman",
         "lae",
         "lqcontrol",

test/util.jl

@@ -25,4 +25,3 @@ if !(isfile(quad_data_file_name))
         warn(m)
     end
 end
-