May 29, 2012

Hodrick-Prescott Filter

There are loads of Hodrick-Prescott filters floating around the web for different programming languages, but wasn't able to find one very quickly for python. I wrote one up based on some MATLAB code written by Wilmer Henao.

Here are links to the filter and the example below.

import scipy as sp
from matplotlib.pyplot import *
from scipy import linalg as la
from scipy import sparse

def hp_filter(y,w):
    # make sure the inputs are the right shape
    m,n  = y.shape
    if m < n:
        y = y.T
        m = n

    a    = sp.array([w, -4*w, ((6*w+1)/2.)])
    d    = sp.tile(a, (m,1))

    d[0,1]   = -2.*w
    d[m-2,1] = -2.*w
    d[0,2]   = (1+w)/2.
    d[m-1,2] = (1+w)/2.
    d[1,2]   = (5*w+1)/2.
    d[m-2,2] = (5*w+1)/2.

    B = sparse.spdiags(d.T, [-2,-1,0], m, m)
    B = B+B.T

    # report the filtered series, s
    s =,y)
    return s

def testhp():
    # read in and assign variables
    data = sp.loadtxt('investment.txt')
    data = sp.log(data)
    y    = sp.atleast_2d(data)
    s    = hp_filter(y,1600)
    devs = y.T - s

    # plot the data, the filtered series, and abs. deviation
    plot(data, 'k')
    title("Data and HP filter")
    ylabel("log of investment")
    title("Stationary series")
    ylabel("log of investment")
    xlabel("Quarters from 1947Q1-2011Q4")

1 comment:

  1. (Tried before but my comment vanished).

    As given the problem does not exploit sparsity as you are converting to dense. It takes ages and about 5Gb of memory to solve a 10,000 point problem. Use instead a symmetric band solver such as la.solveh_banded() which runs in a fraction of a second and doesn't use huge amounts of memory.