I am interested in implementing Hodrick–Prescott filter (HP filter) on my time series data in Mathematica.
I know that there is already a version of two-sided original filter for Mathematica (see Demo 5161), however wanted to know if anyone has implemented the one-sided version of the filter and would like to share?
So far I found only an implementation for Matlab.
All a one-sided filter is, is the last data point of a two-sided filter, for each subsample of the time series that begins with its first period, in other words, for each sample
$x_1....x_k$, for $k = l...T$, where $l$ is a small number of which the HP filter can still be computed, and $T$ is the total length of the time series.
The code due to Ludsteck and Schlicht available on MathLibrary still works in recent versions if you are setting your smoothing parameter explicitly. You can then write a small function to brute-force a one-sided filter like this. Here I am assuming $l$ above is 3.
OneSidedHPFilter[data:{__?NumericQ}, a_?Positive] :=
Table[Last[HPFilter[Take[data, i], a]], {i, 3, Length[data]}]
There is also an old, old implementation I wrote in the late 1990s available on my web site, but I no longer recommend it as it is about three times as slow as the Schlicht version. The version on the web site also has an obsolete reference to the old Statistics`DataSmoothing package, but it works fine if you delete that reference from the package definition.
As far as I can tell, the only problem with the Ludsteck-Schlicht code is that it has references to FilterOptions that needs to be changed to FilterRules or an OptionsPattern/OptionValue construct. See this post by Leonid Shifrin for more information.
I have fixed the Ludsteck–Schlicht code to make it version 9 compatible. It is available at http://www.verbeia.com/mathematica/mma/HPFilter.m
I added the one-sided case based on the implementation above, including a version that works out the optimal smoothing parameter for the whole sample, and uses that for the subsamples.
testdata = Accumulate[RandomVariate[NormalDistribution[0, 0.5], 200]];
Get["HPFilter`HPFilter`"]
hptest = HPFilter[testdata, 1600];
ListLinePlot[{testdata, hptest}]
If you get an error message with the optimal-smoothing version, it still works.
optimalhp = HPFilter[testdata];
(* Eigenvalues::arh: Because finding 200 out of the 200 eigenvalues and/or eigenvectors is
likely to be faster with dense matrix methods, the sparse input matrix will be converted.
If fewer eigenvalues and/or eigenvectors would be sufficient, consider restricting this
number using the second argument to Eigenvalues. >> *)
ListLinePlot[{testdata, First@optimalhp}]
The following shows the issue with the one-sided version: it inherently has a phase shift in it because it only uses lagging information. Note that if you start with an initial sample of 3, you have to drop the first two points in the original data to get it to match up.
oneside = OneSidedHPFilter[testdata, 1600];
ListLinePlot[{Drop[testdata, 2], oneside}]
