Hodrick-Prescott (HP) Filter in R
Previously, the HP filter was addressed using Python. This post presents an R version.
The Hodrick-Prescott filter models a time series (\(y_t\)) as \[\begin{align} y_t = \tau_t + c_t \end{align}\] where \(\tau_t\) and \(c_t\) are trend and cyclical components respectively.
To decompose in this way, we solve the following minimization problem. \[\begin{align} \min_{\tau} \left( \sum_{t=1}^{T} (y_t - \tau_t)^2 + \lambda \sum_{t=2}^{T-1} (\Delta \tau_{t+1} - \Delta \tau_t)^2 \right) \end{align}\] \(\lambda\) is typically recommended as 100, 1600, and 14400 for annual, quarterly, and monthly frequency, respectively.
R code
The hpfilter R library provides two convenient functions for the HP filter: hp1() for the one-sided filter and hp2() for the two-sided filter. They are easy to use, as demonstrated in the code below.
It's interesting that plotting results requires more lines of code compared to the HP filter code (just two lines), which is the main topic of this post.
#========================================================# # Quantitative Financial Econometrics & Derivatives # ML/DL using R, Python, Tensorflow by Sang-Heon Lee # # https://shleeai.blogspot.com #--------------------------------------------------------- # Hodrick-Prescott (HP) filter - 1 sided & 2 sided #========================================================# graphics.off(); rm(list = ls()) library(quantmod) library(hpfilter) #------------------------------------------ # Download monthly AAPL stock data #------------------------------------------ aapl <- getSymbols("AAPL", from = "2011-01-01", to = "2022-12-31", auto.assign = FALSE, periodicity = "monthly") df <- as.data.frame(Ad(aapl)) date <- index(aapl) #as.data.frame(index(aapl)) y <- df; colnames(df) <- "AAPL" #------------------------------------------ # Apply the HP filter (1-sided, 2-sided) #------------------------------------------ yt1 = hp1(y, lambda = 14400); yc1 = y - yt1 yt2 = hp2(y, lambda = 14400); yc2 = y - yt2 #------------------------------------------ # Plot the three resulting series #------------------------------------------ x11(width=7, height=6); plot(y$AAPL, type="l", col=1, lty=1, lwd=3, ylab="", xlab="", xaxt = "n", ylim=c(min(yc1, yt1)*1.2,max(yc1, yt1)*1.1)) lines(yt1$AAPL, lwd=5, col=rgb(1,0,0, alpha=0.5)) lines(yt2$AAPL, lwd=5, col=rgb(0,0,1, alpha=0.5)) # Plotting polygons with 50% transparency polygon(c(1, seq(yc1$AAPL), length(yc1$AAPL)), c(0, yc1$AAPL, 0), col = rgb(1, 0, 0, alpha = 0.3)) polygon(c(1, seq(yc2$AAPL), length(yc2$AAPL)), c(0, yc2$AAPL, 0), col = rgb(0, 0, 1, alpha = 0.3)) # Add x-axis with selected dates sel_dates <- date[seq(18, length(date), by = 24)] axis(1, at = which(date %in% sel_dates), labels = substr(sel_dates,1,4)) # Legend legend("topleft", cex=0.8, lty = 1, col = c(1,2,4,2,4), lwd=c(3,5,5,3,3), c("Data", "Trend (1-sided)", "Trend (2-sided)", "Cycle (1-sided)", "Cycle (2-sided)")) | cs |
The outputs of the two HP filters are the same as those in the Python cases.
No comments:
Post a Comment