/* Hodrick-Prescott Filter (Matheny,vanNorden,Vigfusson) (Paul Soderlin's HPMat2 program follows this, and Lubuele's program follows that.) */ /* HPFILTER First Draft : December 21, 1989 by Ken Matheny (hodpres()) Second Draft: August 3, 1994 by Simon van Norden (hpfilter()) Third Draft : April 27, 1995 by Robert Vigfusson (hpfilter()) This program is written and submitted for public, non-commercial use. The authors make no guarantees as to its performance, and request that use of this code be acknowledged in any published research, including working papers. This program is designed to smooth the stochastic series y by use of the Hodrick-Prescott method mentioned in footnotes of the Kydland & Prescott "Time to Build" paper, and Hansen's "Indivisible Labor" paper. The method involves minimizing an expression containing the squared differences between the actual and smoothed series, and terms involving the differences between successive elements of the smoothed series. This leads to a linear filter for the original data. "w" is a smoothness parameter. The larger w, the smoother the resulting series, shp. Kydland and Prescott suggest a "w" of 1600 for quarterly data. FORMAT : {s} = hpfilter(y,w); INPUTS : y : The series to be filtered; w : The smoothness parameter (suggest 1600). OUTPUT : shp: The smoothed series. The deviations can be calculated simply as (y-shp). */ proc(1) = hpfilter(y,w); local t,a,b,c,d,m,i,j,shp; t = rows(y); a = 6*w+1; b = -4*w; c = w; d = c~b~a; d = ones(t,1)*d; d[1,1]=0; d[2,1]=0; d[1,2]=0; d[2,2] = -2*w; d[t,2] = -2*w; d[1,3] = 1+w; d[t,3] = 1+w; d[2,3]=5*w+1; d[t-1,3]=5*w+1; @ The Smoothed Series @ shp = bandsolpd(eye(t),d)*y; retp(shp); endp; /******************************** Gaston Luis Repetto grepet1@icarus.cc.uic.edu offered the following comments: For monthly data, Edward Prescott suggested (by email) to begin the analysis with 14,400 and then try to see with which value the adjusted series better mimic the behavior of the original one. "There is no such a thing as a "right" smoothing parameter. For purposes of postwar business cycles, the 1600 number was the first tried and seem to mimic well the smooth curve people typically fit throughout the data. As a year is four quarters, reducing it by four squared, that is by factor of 16, resulted in the same amount of variation in trend being permitted. This logic suggest multiplyings 1600 by 9 (three squared) to have same amount of variability in trend with monthly data as 1600 permits for quarterly data. Some good people tried to develop some formal theory, that made connections between different length time periods. The "right" method depended upon many very specific probabilistic assumptions. But this is a data representation exercise, not a measurement or estimation procedure." ***********************************/ /*----------------------------------------------------------------------- ** HPMat2 ** Author: Paul Soderlind ** An efficient HP-filter. Any comments are welcome. ** Provided without guarantees for free public, non-commercial use. ** Date: 30 May 1997 ** ** ** Purpose: HP-filtering of a several time series, assuming no missing ** values. ** ** Usage: { gy,cy } = HPMat2( y,lambda ); ** ** Input: y NxM matrix with M time series ** lambda degree of smoothness (1600 is a common value) ** ** Output: gy NxM matrix of trend component ** cy NxM matrix of cyclical component ** ** ** ** Note: (a) This procedure uses Toeplitz() to set up the K matrix, ** and the band matrix routines for inverting. This gives ** a fast code with very small storage requirements. ** ** (b) The trend component, gy, is calculated as ** cy = inv(A) * y, where A = I + lambda*K'K. The cyclical ** component is then cy = y - gy. The matrix K is as follows ** ** 1 -2 1 0 0 ... 0 0 0 ** 0 1 -2 1 0 ... 0 0 0 ** K = 0 0 1 -2 1 ... 0 0 0 ** . . . . . . . . ** . . . . . . . . ** 0 0 0 0 0 ... 1 -2 1 ** ** ** Reference: Danthine and Girardin (1989), "Business Cycles In Switzerland: ** A Comparative Study," European Economic Review 33, 31-50. ** ** ** Paul Soderlind, 30 May 1997. -----------------------------------------------------------------------*/ proc (2) = HPMat2( y,lambda ); local N,gk,K,A,bandA,CholA,gy,cy; N = rows(y); if N < 5; /*pattern of A different if N < 5*/ retp( error(0),error(0) ); endif; gk = 1|(-2)|1|zeros(2,1); /*prototype 5x5 K matrix*/ K = upmat(toeplitz(gk)); K = trimr(k,0,2); A = eye(5)+lambda*k'*k; /*prototype 5x5 A matrix*/ bandA = band(A,2); /*band from A, extending to TxT matrix*/ bandA = bandA[1:2,.] | (ones(N-4,1).*bandA[3,.]) | bandA[4:5,.]; CholA = bandchol( bandA ); /*Cholesky decomp of band representation*/ gy = bandcholsol( y,CholA ); /*solving for "trend" component*/ cy = y - gy; /*"cyclical" component*/ retp( gy,cy ); endp; /*-----------------------------------------------------------------------*/ /* ANOTHER HPFILTER: Author: Luan Lubuele Senior Economist City of New York Independent Budget Office lubuele@nwu.edu Date: 11 May 1997 Provided without guarantees for free public, non-commercial use. The HP filter code posted by M. Heiko uses the inverse method. It should do the job and is the shortest and fastest method. It would be simpler to write the code as a procedure yielding directly the trend so as to avoid holding a huge inverse in memory. There is however the issue of interdmediate swelling and loss of accuracy. Indeed the matrix to be inverted is of the size n^2 (n being the sample size), which can get big for say weekly data for a long period of time. Plus, there will be a point where Gauss will need to hold temporarely another matrix of the same size. As a result, one may run out of memory. On the other end, the use of the straight inverse function on a huge matrix may yield less accurate results. I have seen occasions when the function inv(.) failed to pick up a clear case of a matrix with incomplete rank. I am enclosing the HP-filter procedure which (1) does not need an inverse, (2) uses just storing space of the order of 5*n instead of n^2. Its shortfalls are (1) it is slow since it relies on loops, and (2) it does not allow the computation of the hptrends of many series at once. I am also enclosing the MATLAB version of the method based on the inverse, just in case there is now a way of handling sparse matrices and someone is interested in translating in Gauss. The problems discussed above are addressed by using sparse matrix methods and LU decomposition. /****************** Lubuele GAUSS CODE STARTS HERE **********************/ declare __lambda != 1600; proc (2) = hpfilter(y); /* Computes the Hodrck-Prescott filtered trend and cycles. ** Based on Prescott algorithm using the two pass kalman filter approach ** This code is a translation of Prescott,s fortran code. ** The procedure needs the global variable __lambda (default=1600). */ local s,n,ierr,t,d,v,m1,m2,v11,v12,v22,x,z,b11,b12,b22,dett,e1,e2; local istep,i1,i2,i; s=__lambda; n=rows(y); t=zeros(n,1); d=t; v=zeros(n,3); if s<= 0; s=1600; endif; if n<=3; retp(0,0); endif; istep=1; m1=y[2]; m2=y[1]; i1=3; i2=n; v11=1; v22=1; v12=0; i=i1; do while i<=i2; /* first pass */ gosub pass; endo; t[n]=m1; t[n-1]=m2; istep=-1; m1=y[n-1]; m2=y[n]; i1=n-2; i2=1; v11=1; v22=1; v12=0; i=i1; do while i>=i2; /* second backward pass */ gosub pass; endo; t[1]=m1; t[2]=m2; i=1; do while i<=n; d[i]=y[i]-t[i]; i=i+1; endo; retp(t,d); pass: x=m1; m1=2*m1-m2; m2=x; x=v11; z=v12; v11=1/s+4*(x-z)+v22; v12=2*x-z; v22=x; dett=v11*v22-v12*v12; if istep==1; v[i-1,1]=v22/dett; v[i-1,3]=v11/dett; v[i-1,2]=-v12/dett; t[i-1]=v[i-1,1]*m1+v[i-1,2]*m2; d[i-1]=v[i-1,2]*m1+v[i-1,3]*m2; else; if i>=2; b11=v11/dett; b12=-v12/dett; b22=v22/dett; e1=b11*m2+b12*m1+t[i]; e2=b12*m2+b22*m1+d[i]; b12=b12+v[i,2]; b22=b22+v[i,3]; b11=b11+v[i,1]; dett=b11*b22-b12*b12; t[i]=(-b12*e1+b11*e2)/dett; endif; endif; x=v11+1; z=(y[i]-m1)/x; m1=m1+v11*z; m2=m2+v12*z; z=v11; v11=v11-v11*v11/x; v22=v22-v12*v12/x; v12=v12-z*v12/x; i=i+istep; return; endp; /*********************** Lubuele GAUSS CODE ENDS HERE ***************************/