$ontext

   L1 regression with a trend break
   
   Erwin Kalvelagen, oct 2003
   
   
$offtext

set i 'number of observations' /i1*i50/;

parameter x(i);
parameter y(i);

*
* generate test data
* x must be ordered: x(i) >= x(i-1).
*
set i1(i) /i1*i30/;
set i2(i) /i31*i50/;

x(i) = 0;
loop(i,
   x(i) = x(i-1) + uniform(0.1,10);
);
display x;

y(i1) = 200 - 2*x(i1) + normal(0,10);
y(i2) = -220 + x(i2) + normal(0,10); 

display y;


*
* coefficients to estimate
*
variables a1,b1,a2,b2,break;

*
* error terms
*
variables e1(i),e2(i);

*
* equations to fit
*
equations fit1(i),fit2(i);
fit1(i)..  y(i) =e= a1 + b1*x(i) + e1(i);
fit2(i)..  y(i) =e= a2 + b2*x(i) + e2(i);
 

*
* we use binary variables b(i) to indicate
* if x(i) is in segment 1 or segment 2.
*
binary variable b(i);

equations seg1(i),seg2(i);
scalar m1; m1 = smax(i,x(i))-smin(i,x(i));

seg1(i).. x(i) =l= break + b(i)*M1;
seg2(i).. x(i) =g= break - (1-b(i))*M1;


variable u(i) 'either e1(i) or e2(i)';
equation udef1(i),udef2(i),udef3(i),udef4(i);

*
* 10*range of y is generous for m2
*
scalar m2;
m2 = 10*[smax(i,y(i))-smin(i,y(i))];

udef1(i).. u(i) =g= e1(i) - b(i)*m2;
udef2(i).. u(i) =l= e1(i) + b(i)*m2;
udef3(i).. u(i) =g= e2(i) - (1-b(i))*m2;
udef4(i).. u(i) =l= e2(i) + (1-b(i))*m2;

*
* variable splitting
*
positive variable v(i),w(i);
equation split(i);
split(i)..  v(i) - w(i) =e= u(i);

*
* objective
*
variable l1;
equation obj;

obj.. l1 =e= sum(i, v(i) + w(i));

model m/all/;
option optcr=0;
option mip=cplex;
solve m minimizing l1 using mip;
      
parameter res(i,*);
res(i,'x') = x(i);
res(i,'b') = b.l(i);
res(i,'u') = u.l(i);
res(i,'e1') = e1.l(i);
res(i,'e2') = e2.l(i);

display break.l,res;

      

*
* produce some graphs
*
file dat0 /l0.dat/;
file dat1 /l1.dat/;
file dat2 /l2.dat/;
loop(i,
   put dat0,x(i):12:7, y(i):12:7/
   if (b.l(i)<0.5,
      put dat1,x(i):12:7, (y(i)-e1.l(i)):12:7/
   else
      put dat2,x(i):12:7, (y(i)-e2.l(i)):12:7/
   );      
);
putclose dat0;
putclose dat1;
putclose dat2;

file plt /l1.plt/;
putclose plt,
    'plot "l0.dat" title "data",'
    ' "l1.dat" title "segment 1" with lines,'
    ' "l2.dat" title "segment 2" with lines'/
    ' pause -1'/
execute 'gnuplot l1.plt';