Hamilton Filtering

In this post I present a R function that implements James D. Hamilton's alternative to the Hodrick-Prescott Filter as described in his article "Why You Should Never Use the Hodrick Prescott Filter", forthcoming in the Review of Economics and Statistics. In the paper, Hamilton demonstrates how HP filtering creates artificial correlation between variables and suggests filtering using the residuals of the following regression:
$$ y_{t+h} = \beta_0 + \sum_{l=1}^p \beta_l y_{t+1-l} + \nu_{t+h}, $$
$y_{t+h}$ are the elements of $y$ that can be predicted $h$ periods ahead using its previous $l$ observations. This corresponds to the trend after fluctuations have disappeared after $h$ periods. $h$ should therefore be chosen to be 2 years where we assume macroeconomic shocks to be worn off. $l$ is recommended to include 1 year of data to remove possible seasonal components from the trend. I highly recommend reading Hamilton's blog post on the topic.
I implement this filter in the following R function based on the matlab code accompanying Hamilton's paper.

hamfilter <- function(yt, h=2*frequency(y), p=frequency(y), trend=FALSE){ # Hamilton Filter # # R code to calculate cyclical component based on h-period-ahead #forecast error from linear regression # y(t+h) = b0 + b1*y(t) + ... + bp*y(t+1-p) + e(t+h) # as recommended in # James D. Hamilton, "Why You Should Never Use the Hodrick-Prescott Filter" # Review of Economics and Statistics, forthcoming # Code based on James D. Hamilton's matlab code # inputs: # y = (T x 1) vector or time series, tth element is observation for date t # p = number of lags that should be considered in the regression of y(t+h), # default is the frequency of a time series # outputs: # yreg = (T x 1) vector or time series, tth element is cyclical component for date t # (optional) ytrend = (T x 1) vector or time series, tth element is trend component # # works with any frequency of data # works with leading and trailing missing values (not internal) # preserves class of input data (i.e. time series, vector, etc.) # require(zoo) # create final output matrices to preserve input class and length yreg <- yt yreg[] <- NaN ytrend <- yt ytrend[] <- NaN # prepare data y <- na.omit(yt) T <- length(y) lseq <- 0:-(p-1) X <- cbind(rep(1,T-p-h+1), as.matrix(lag(as.zoo(y),lseq)[((p):(T-h)),])) yh <- y[(p+h):T] # extimate and predict y(t+h) = b0 + b1*y(t) + ... + bp*y(p-1) + e(t+h) b <- solve(t(X) %*% X) %*% (t(X) %*% yh) ytrend[!is.na(yt)] <- c(rep(NaN,p+h-1),X%*%b) yreg[!is.na(yt)] <- y - c(rep(NaN,p+h-1),X%*%b) if(trend==TRUE){ return(cbind(ytrend,yreg)) } else { return(yreg) } }