A derivative work by Judson Wilson, 5/28/2014.
Adapted from the CVX example of the same name, by Kwangmoo Koh, 12/10/2007
Topic Reference:
The problem of estimating underlying trends in time series data arises in a variety of disciplines. The $\ell_1$ trend filtering method produces trend estimates $x$ that are piecewise linear from the time series $y$.
The $\ell_1$ trend estimation problem can be formulated as \begin{array}{ll} \mbox{minimize} & (1/2)||y-x||_2^2 + \lambda ||Dx||_1, \end{array} with variable $x$ , and problem data $y$ and $\lambda$, with $\lambda >0$. $D$ is the second difference matrix, with rows $$\begin{bmatrix}0 & \cdots & 0 & -1 & 2 & -1 & 0 & \cdots & 0 \end{bmatrix}.$$ CVXPY is not optimized for the $\ell_1$ trend filtering problem. For large problems, use l1_tf (http://www.stanford.edu/~boyd/l1_tf/).
In [1]:
import numpy as np
import cvxpy as cp
import scipy as scipy
import cvxopt as cvxopt
# Load time series data: S&P 500 price log.
y = np.loadtxt(open('data/snp500.txt', 'rb'), delimiter=",", skiprows=1)
n = y.size
# Form second difference matrix.
e = np.ones((1, n))
D = scipy.sparse.spdiags(np.vstack((e, -2*e, e)), range(3), n-2, n)
# Set regularization parameter.
vlambda = 50
# Solve l1 trend filtering problem.
x = cp.Variable(shape=n)
obj = cp.Minimize(0.5 * cp.sum_squares(y - x)
+ vlambda * cp.norm(D*x, 1) )
prob = cp.Problem(obj)
# ECOS and SCS solvers fail to converge before
# the iteration limit. Use CVXOPT instead.
prob.solve(solver=cp.CVXOPT, verbose=True)
print('Solver status: {}'.format(prob.status))
# Check for error.
if prob.status != cp.OPTIMAL:
raise Exception("Solver did not converge!")
print("optimal objective value: {}".format(obj.value))
In [2]:
import matplotlib.pyplot as plt
# Show plots inline in ipython.
%matplotlib inline
# Plot properties.
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
font = {'weight' : 'normal',
'size' : 16}
plt.rc('font', **font)
# Plot estimated trend with original signal.
plt.figure(figsize=(6, 6))
plt.plot(np.arange(1,n+1), y, 'k:', linewidth=1.0)
plt.plot(np.arange(1,n+1), np.array(x.value), 'b-', linewidth=2.0)
plt.xlabel('date')
plt.ylabel('log price')
Out[2]: