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 = sp.dot(la.inv(B.todense()),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 subplot(211) plot(data, 'k') plot(s,'r-') title("Data and HP filter") ylabel("log of investment") subplot(212) plot(devs,'k') title("Stationary series") ylabel("log of investment") xlabel("Quarters from 1947Q1-2011Q4") show()

(Tried before but my comment vanished).

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