# L1 Trend Filtering

## Introduction

Kim et al. (2009) propose the $$l_1$$ trend filtering method for trend estimation. The method solves an optimization problem of the form

$\begin{array}{ll} \underset{\beta}{\mbox{minimize}} & \frac{1}{2}\sum_{i=1}^m (y_i - \beta_i)^2 + \lambda ||D\beta||_1 \end{array}$

where the variable to be estimated is $$\beta$$ and we are given the problem data $$y$$ and $$\lambda$$. The matrix $$D$$ is the second-order difference matrix,

$D = \left[ \begin{matrix} 1 & -2 & 1 & & & & \\ & 1 & -2 & 1 & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & 1 & -2 & 1 & \\ & & & & 1 & -2 & 1\\ \end{matrix} \right].$

The implementation is in both C and Matlab. Hadley Wickham provides an R interface to the C code. This is on GitHub and can be installed via:

library(devtools)

## Example

We will use the example in l1tf to illustrate. The package provides the function l1tf which computes the trend estimate for a specified $$\lambda$$.

library(l1tf)
sp_data <- data.frame(x = sp500$date, y = sp500$log,
l1_50 = l1tf(sp500$log, lambda = 50), l1_100 = l1tf(sp500$log, lambda = 100))

## The CVXR version

CVXR provides all the atoms and functions necessary to formulat the problem in a few lines. For example, the $$D$$ matrix above is provided by the function diff(..., differences = 2). Notice how the formulation tracks the mathematical construct above.

## lambda = 50
y <- sp500$log lambda_1 <- 50 beta <- Variable(length(y)) objective_1 <- Minimize(0.5 * p_norm(y - beta) + lambda_1 * p_norm(diff(x = beta, differences = 2), 1)) p1 <- Problem(objective_1) betaHat_50 <- solve(p1)$getValue(beta)

## lambda = 100
lambda_2 <- 100
objective_2 <- Minimize(0.5 * p_norm(y - beta) +
lambda_2 * p_norm(diff(x = beta, differences = 2), 1))
p2 <- Problem(objective_2)
betaHat_100 <- solve(p2)$getValue(beta) NOTE Of course, CVXR is much slower since it is not optimized just for one problem. ## Comparison Plots A plot of the estimates for two values of $$\lambda$$ is shown below using both approaches. First the l1tf plot. ggplot(data = sp_data) + geom_line(mapping = aes(x = x, y = y), color = 'grey50') + labs(x = "Date", y = "SP500 log-price") + geom_line(mapping = aes(x = x, y = l1_50), color = 'red', size = 1) + geom_line(mapping = aes(x = x, y = l1_100), color = 'blue', size = 1) ## Warning in as.POSIXlt.POSIXct(x): unknown timezone 'zone/tz/2017c.1.0/ ## zoneinfo/US/Pacific' Next the corresponding CVXR plots. cvxr_data <- data.frame(x = sp500$date,
y = sp500\$log,
l1_50 = betaHat_50,
l1_100 = betaHat_100)
ggplot(data = cvxr_data) +
geom_line(mapping = aes(x = x, y = y), color = 'grey50') +
labs(x = "Date", y = "SP500 log-price") +
geom_line(mapping = aes(x = x, y = l1_50), color = 'red', size = 1) +
geom_line(mapping = aes(x = x, y = l1_100), color = 'blue', size = 1)

## Notes

The CVXR solution is not quite exactly that of l1tf: on the left it shows a larger difference for the two $$\lambda$$ values; in the middle, it is less flatter than l1tf; and on the right, it does not have as many knots as l1tf.

## Session Info

sessionInfo()
## R version 3.4.2 (2017-09-28)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] methods   stats     graphics  grDevices datasets  utils     base
##
## other attached packages:
## [1] l1tf_0.0.0.9000 CVXR_0.94-4     ggplot2_2.2.1
##
## loaded via a namespace (and not attached):
##  [1] gmp_0.5-13.1      Rcpp_0.12.13      highr_0.6
##  [4] compiler_3.4.2    plyr_1.8.4        R.methodsS3_1.7.1
##  [7] R.utils_2.6.0     tools_3.4.2       digest_0.6.12
## [10] bit_1.1-12        evaluate_0.10.1   tibble_1.3.4
## [13] gtable_0.2.0      lattice_0.20-35   rlang_0.1.2
## [16] Matrix_1.2-11     yaml_2.1.14       blogdown_0.1.7
## [19] Rmpfr_0.6-1       ECOSolveR_0.3-2   stringr_1.2.0
## [22] knitr_1.17        rprojroot_1.2     bit64_0.9-7
## [25] grid_3.4.2        R6_2.2.2          rmarkdown_1.6
## [28] bookdown_0.5      magrittr_1.5      backports_1.1.1
## [31] scales_0.5.0      htmltools_0.3.6   scs_1.1-1
## [34] colorspace_1.3-2  labeling_0.3      stringi_1.1.5
## [37] lazyeval_0.2.1    munsell_0.4.3     R.oo_1.21.0

R Markdown

## References

Kim, Seung-Jean, Kwangmoo Koh, Stephen Boyd, and Dimitry Gorinevsky. 2009. “$$l_1$$ Trend Filtering.” SIAM Review 51 (2): 339–60. doi:doi:10.1137/070690274.