#!/usr/bin/env python
# coding: utf-8

# # $\ell_1$ trend filtering
# 
# A derivative work by Judson Wilson, 5/28/2014.<br>
# Adapted from the CVX example of the same name, by Kwangmoo Koh, 12/10/2007
# 
# Topic Reference:
# 
# * S.-J. Kim, K. Koh, S. Boyd, and D. Gorinevsky, ``l_1 Trend Filtering''<br>
# http://stanford.edu/~boyd/papers/l1_trend_filter.html
# 
# ## Introduction
# 
# 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/).
# 
# ## Formulate and solve problem

# 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))


# ## Results plot

# In[2]:


import matplotlib.pyplot as plt

# Show plots inline in ipython.
get_ipython().run_line_magic('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')