{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# $\\ell_1$ trend filtering\n",
    "\n",
    "A derivative work by Judson Wilson, 5/28/2014.<br>\n",
    "Adapted from the CVX example of the same name, by Kwangmoo Koh, 12/10/2007\n",
    "\n",
    "Topic Reference:\n",
    "\n",
    "* S.-J. Kim, K. Koh, S. Boyd, and D. Gorinevsky, ``l_1 Trend Filtering''<br>\n",
    "http://stanford.edu/~boyd/papers/l1_trend_filter.html\n",
    "\n",
    "## Introduction\n",
    "\n",
    "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$.\n",
    "\n",
    "The $\\ell_1$ trend estimation problem can be formulated as\n",
    "    \\begin{array}{ll}\n",
    "    \\mbox{minimize}   &  (1/2)||y-x||_2^2 + \\lambda ||Dx||_1,\n",
    "    \\end{array}\n",
    "with variable $x$ , and problem data $y$ and $\\lambda$, with $\\lambda >0$.\n",
    "$D$ is the second difference matrix, with rows \n",
    "    $$\\begin{bmatrix}0 & \\cdots & 0 & -1 & 2 & -1 & 0 & \\cdots & 0 \\end{bmatrix}.$$\n",
    "CVXPY is not optimized for the $\\ell_1$ trend filtering problem.\n",
    "For large problems, use l1_tf (http://www.stanford.edu/~boyd/l1_tf/).\n",
    "\n",
    "## Formulate and solve problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     pcost       dcost       gap    pres   dres   k/t\n",
      " 0:  0.0000e+00 -1.0000e+00  1e+05  1e-01  4e-02  1e+00\n",
      " 1:  2.2350e-01  1.5374e-01  8e+03  8e-03  3e-03  7e-02\n",
      " 2:  1.9086e-01  2.3346e-01  1e+03  1e-03  4e-04  6e-02\n",
      " 3:  3.9403e-01  4.4110e-01  7e+02  7e-04  3e-04  6e-02\n",
      " 4:  3.5979e-01  4.1278e-01  3e+02  3e-04  1e-04  6e-02\n",
      " 5:  6.4154e-01  6.4522e-01  2e+01  2e-05  7e-06  4e-03\n",
      " 6:  9.0480e-01  9.0710e-01  1e+01  1e-05  4e-06  3e-03\n",
      " 7:  9.9603e-01  9.9825e-01  1e+01  1e-05  4e-06  2e-03\n",
      " 8:  1.0529e+00  1.0542e+00  6e+00  6e-06  2e-06  1e-03\n",
      " 9:  1.1994e+00  1.2004e+00  4e+00  4e-06  2e-06  1e-03\n",
      "10:  1.2689e+00  1.2693e+00  2e+00  2e-06  6e-07  4e-04\n",
      "11:  1.3728e+00  1.3729e+00  5e-01  5e-07  2e-07  1e-04\n",
      "12:  1.3802e+00  1.3803e+00  2e-01  2e-07  9e-08  6e-05\n",
      "13:  1.3965e+00  1.3965e+00  1e-01  1e-07  4e-08  3e-05\n",
      "14:  1.3998e+00  1.3998e+00  3e-02  3e-08  1e-08  8e-06\n",
      "15:  1.3999e+00  1.3999e+00  3e-02  3e-08  1e-08  7e-06\n",
      "16:  1.4011e+00  1.4011e+00  9e-03  9e-09  3e-09  2e-06\n",
      "17:  1.4013e+00  1.4013e+00  3e-03  3e-09  1e-09  8e-07\n",
      "18:  1.4014e+00  1.4014e+00  6e-04  6e-10  3e-10  2e-07\n",
      "19:  1.4014e+00  1.4014e+00  2e-04  2e-10  7e-11  4e-08\n",
      "20:  1.4014e+00  1.4017e+00  4e-05  4e-11  2e-08  1e-08\n",
      "21:  1.4014e+00  1.4015e+00  3e-06  4e-12  8e-08  7e-10\n",
      "22:  1.4014e+00  1.4013e+00  4e-08  4e-13  2e-08  9e-12\n",
      "Optimal solution found.\n",
      "Solver status: optimal\n",
      "optimal objective value: 1.4014300716775199\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import cvxpy as cp\n",
    "import scipy as scipy\n",
    "import cvxopt as cvxopt\n",
    "\n",
    "# Load time series data: S&P 500 price log.\n",
    "y = np.loadtxt(open('data/snp500.txt', 'rb'), delimiter=\",\", skiprows=1)\n",
    "n = y.size\n",
    "\n",
    "# Form second difference matrix.\n",
    "e = np.ones((1, n))\n",
    "D = scipy.sparse.spdiags(np.vstack((e, -2*e, e)), range(3), n-2, n)\n",
    "\n",
    "# Set regularization parameter.\n",
    "vlambda = 50\n",
    "\n",
    "# Solve l1 trend filtering problem.\n",
    "x = cp.Variable(shape=n)\n",
    "obj = cp.Minimize(0.5 * cp.sum_squares(y - x)\n",
    "                  + vlambda * cp.norm(D*x, 1) )\n",
    "prob = cp.Problem(obj)\n",
    "\n",
    "# ECOS and SCS solvers fail to converge before\n",
    "# the iteration limit. Use CVXOPT instead.\n",
    "prob.solve(solver=cp.CVXOPT, verbose=True)\n",
    "print('Solver status: {}'.format(prob.status))\n",
    "\n",
    "# Check for error.\n",
    "if prob.status != cp.OPTIMAL:\n",
    "    raise Exception(\"Solver did not converge!\")\n",
    "\n",
    "print(\"optimal objective value: {}\".format(obj.value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'log price')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 44836 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Show plots inline in ipython.\n",
    "%matplotlib inline\n",
    "\n",
    "# Plot properties.\n",
    "plt.rc('text', usetex=True)\n",
    "plt.rc('font', family='serif')\n",
    "font = {'weight' : 'normal',\n",
    "        'size'   : 16}\n",
    "plt.rc('font', **font)\n",
    "\n",
    "# Plot estimated trend with original signal.\n",
    "plt.figure(figsize=(6, 6))\n",
    "plt.plot(np.arange(1,n+1), y, 'k:', linewidth=1.0)\n",
    "plt.plot(np.arange(1,n+1), np.array(x.value), 'b-', linewidth=2.0)\n",
    "plt.xlabel('date')\n",
    "plt.ylabel('log price')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
