:math:`\ell_1` trend filtering
==============================

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:

-  S.-J. Kim, K. Koh, S. Boyd, and D. Gorinevsky, \`\`l\_1 Trend
   Filtering'' 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 :math:`\ell_1` trend filtering method
produces trend estimates :math:`x` that are piecewise linear from the
time series :math:`y`.

The :math:`\ell_1` trend estimation problem can be formulated as

.. math::

   \begin{array}{ll}
       \mbox{minimize}   &  (1/2)||y-x||_2^2 + \lambda ||Dx||_1,
       \end{array}

with variable :math:`x` , and problem data :math:`y` and
:math:`\lambda`, with :math:`\lambda >0`. :math:`D` is the second
difference matrix, with rows

.. math:: \begin{bmatrix}0 & \cdots & 0 & -1 & 2 & -1 & 0 & \cdots & 0 \end{bmatrix}.

CVXPY is not optimized for the :math:`\ell_1` trend filtering problem.
For large problems, use l1\_tf (https://www.stanford.edu/~boyd/l1\_tf/).

Formulate and solve problem
---------------------------

.. code:: python

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


.. parsed-literal::

         pcost       dcost       gap    pres   dres   k/t
     0:  0.0000e+00 -1.0000e+00  1e+05  1e-01  4e-02  1e+00
     1:  2.2350e-01  1.5374e-01  8e+03  8e-03  3e-03  7e-02
     2:  1.9086e-01  2.3346e-01  1e+03  1e-03  4e-04  6e-02
     3:  3.9403e-01  4.4110e-01  7e+02  7e-04  3e-04  6e-02
     4:  3.5979e-01  4.1278e-01  3e+02  3e-04  1e-04  6e-02
     5:  6.4154e-01  6.4522e-01  2e+01  2e-05  7e-06  4e-03
     6:  9.0480e-01  9.0710e-01  1e+01  1e-05  4e-06  3e-03
     7:  9.9603e-01  9.9825e-01  1e+01  1e-05  4e-06  2e-03
     8:  1.0529e+00  1.0542e+00  6e+00  6e-06  2e-06  1e-03
     9:  1.1994e+00  1.2004e+00  4e+00  4e-06  2e-06  1e-03
    10:  1.2689e+00  1.2693e+00  2e+00  2e-06  6e-07  4e-04
    11:  1.3728e+00  1.3729e+00  5e-01  5e-07  2e-07  1e-04
    12:  1.3802e+00  1.3803e+00  2e-01  2e-07  9e-08  6e-05
    13:  1.3965e+00  1.3965e+00  1e-01  1e-07  4e-08  3e-05
    14:  1.3998e+00  1.3998e+00  3e-02  3e-08  1e-08  8e-06
    15:  1.3999e+00  1.3999e+00  3e-02  3e-08  1e-08  7e-06
    16:  1.4011e+00  1.4011e+00  9e-03  9e-09  3e-09  2e-06
    17:  1.4013e+00  1.4013e+00  3e-03  3e-09  1e-09  8e-07
    18:  1.4014e+00  1.4014e+00  6e-04  6e-10  3e-10  2e-07
    19:  1.4014e+00  1.4014e+00  2e-04  2e-10  7e-11  4e-08
    20:  1.4014e+00  1.4017e+00  4e-05  4e-11  2e-08  1e-08
    21:  1.4014e+00  1.4015e+00  3e-06  4e-12  8e-08  7e-10
    22:  1.4014e+00  1.4013e+00  4e-08  4e-13  2e-08  9e-12
    Optimal solution found.
    Solver status: optimal
    optimal objective value: 1.4014300716775199


Results plot
------------

.. code:: python

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




.. parsed-literal::

    Text(0, 0.5, 'log price')




.. image:: l1_trend_filter_files/l1_trend_filter_3_1.png