% Written for CVX by Kwangmoo Koh - 12/10/07
%
% The problem of estimating underlying trends in time series data arises in
% a variety of disciplines. The l1 trend filtering method produces trend
% estimates x that are piecewise linear from the time series y.
%
% The l1 trend estimation problem can be formulated as
%
%    minimize    (1/2)*||y-x||^2+lambda*||Dx||_1,
%
% with variable x , and problem data y and lambda, with lambda >0.
% D is the second difference matrix, with rows [0... -1 2 -1 ...0]
%
% CVX is not optimized for the l1 trend filtering problem.
% For large problems, use l1_tf (www.stanford.edu/~boyd/l1_tf/).

% load time series data
y = csvread('snp500.txt'); % log price of snp500
n = length(y);

% form second difference matrix
e = ones(n,1);
D = spdiags([e -2*e e], 0:2, n-2, n);

% set regularization parameter
lambda = 50;

% solve l1 trend filtering problem
cvx_begin
    variable x(n)
    minimize( 0.5*sum_square(y-x)+lambda*norm(D*x,1) )
cvx_end

% plot estimated trend with original signal
figure(1);
plot(1:n,y,'k:','LineWidth',1.0); hold on;
plot(1:n,x,'b-','LineWidth',2.0); hold off;
xlabel('date'); ylabel('log price');
 
Calling sedumi: 5998 variables, 1999 equality constraints
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 1999, order n = 3999, dim = 5999, blocks = 2
nnz(A) = 9991 + 1, nnz(ADA) = 9985, nnz(L) = 5992
Handling 1 + 1 dense columns.
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            3.92E+05 0.000
  1 :  -1.27E+01 2.07E+03 0.000 0.0053 0.9990 0.9990   1.04  1  1  1.0E+00
  2 :  -2.68E+00 4.81E+02 0.000 0.2326 0.9000 0.9000   1.00  1  1  9.9E-01
  3 :   7.05E-02 1.15E+02 0.000 0.2400 0.9000 0.9000   1.00  1  1  9.8E-01
  4 :   6.04E-01 5.10E+01 0.000 0.4419 0.9000 0.9000   1.00  1  1  9.7E-01
  5 :   7.94E-01 2.62E+01 0.000 0.5132 0.9000 0.9000   1.00  1  1  9.5E-01
  6 :   9.62E-01 1.32E+01 0.000 0.5044 0.9061 0.9000   1.00  1  1  9.1E-01
  7 :   1.10E+00 8.32E+00 0.000 0.6297 0.9000 0.9484   1.00  1  1  8.7E-01
  8 :   1.16E+00 5.01E+00 0.000 0.6024 0.9000 0.6707   1.00  1  1  8.0E-01
  9 :   1.25E+00 1.93E+00 0.000 0.3861 0.9247 0.9000   1.00  1  1  6.0E-01
 10 :   1.32E+00 1.14E+00 0.000 0.5890 0.9000 0.9290   1.00  1  1  4.5E-01
 11 :   1.36E+00 7.26E-01 0.016 0.6373 0.9000 0.9450   1.00  1  1  3.4E-01
 12 :   1.38E+00 3.01E-01 0.000 0.4153 0.9255 0.9000   1.00  1  1  1.8E-01
 13 :   1.39E+00 1.76E-01 0.000 0.5839 0.9000 0.9327   1.00  1  1  1.1E-01
 14 :   1.40E+00 1.09E-01 0.051 0.6170 0.9000 0.9442   1.00  1  1  7.1E-02
 15 :   1.40E+00 5.69E-02 0.000 0.5238 0.9088 0.9000   1.00  1  1  3.8E-02
 16 :   1.40E+00 2.57E-02 0.000 0.4508 0.9000 0.9302   1.00  1  1  1.8E-02
 17 :   1.40E+00 1.04E-02 0.000 0.4070 0.9003 0.9000   1.00  1  2  7.2E-03
 18 :   1.40E+00 4.24E-03 0.000 0.4061 0.9000 0.9194   1.00  2  2  3.0E-03
 19 :   1.40E+00 1.40E-03 0.000 0.3292 0.9000 0.9096   1.00  2  2  9.7E-04
 20 :   1.40E+00 2.78E-04 0.000 0.1994 0.9000 0.9038   1.00  2  2  1.9E-04
 21 :   1.40E+00 3.73E-05 0.000 0.1339 0.9073 0.9000   1.00  2  2  2.6E-05
 22 :   1.40E+00 1.27E-06 0.000 0.0341 0.9904 0.9900   1.00  2  2  8.9E-07
 23 :   1.40E+00 2.67E-08 0.000 0.0211 0.9902 0.9900   1.00  2  2  1.9E-08
 24 :   1.40E+00 9.49E-09 0.000 0.3550 0.9344 0.9000   1.00  2  2  6.6E-09

iter seconds digits       c*x               b*y
 24      1.6   8.2  1.4016181622e+00  1.4016181527e+00
|Ax-b| =   4.3e-14, [Ay-c]_+ =   0.0E+00, |x|=  2.1e+00, |y|=  1.5e+03

Detailed timing (sec)
   Pre          IPM          Post
2.000E-02    1.590E+00    1.000E-02    
Max-norms: ||b||=1, ||c|| = 50,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.99236.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +1.40162