/* SAS example for piecewise regression */
/* We use the raw materials data set from class */


DATA rawmater;
  INPUT y x1;
  LABEL  y = 'Cost'
        x1 = 'Shipment Size';
  IF x1 > 500 then x2 = 1; 
    ELSE x2 = 0;
  x2star = (x1 - 500)*x2;
cards;
  2.57  650
  4.40  340
  4.52  400
  1.39  800
  4.75  300
  3.55  570
  2.49  720
  4.27  480
  4.42  425
  4.10  460
  2.21  675
  2.90  600
;
run;

/* Original scatter plot of Y against x1: */

PROC GPLOT DATA=rawmater;
PLOT y*x1;
run;

/* Adding a new data point with an                    */
/* X value at the changepoint (for plotting purposes) */

DATA onemorepoint;
INPUT y x1;
  IF x1 > 500 then x2 = 1; 
    ELSE x2 = 0;
  x2star = (x1 - 500)*x2;
cards;
. 500
;
run;
DATA rawmater;
SET rawmater onemorepoint;
run;

/* Fitting the piecewise regression function with a */
/* changepoint at X = 500:                          */

PROC REG DATA = rawmater;
  MODEL y = x1 x2star;
  /* Recall x2star = (x1 - 500)*x2 here!! */
  OUTPUT out=new p=yhat;
run;

/* Sorting the data by x1   */
/* to prepare for plotting: */

PROC SORT DATA = new;
  BY x1;
run;

/* Plotting the piecewise regression function: */
 
goptions reset=all; 

symbol1 i=join c=black;
symbol2 c=black v=circle;
axis1 label=(angle=90 'Unit Cost');
PROC GPLOT data = new;
  PLOT yhat*x1 y*x1/ overlay vaxis=axis1;
run;


/***************************************************/


/* Allowing for a discontinuity in the regression function: */

DATA rawmater;
  INPUT y x1;
  LABEL  y = 'Cost'
        x1 = 'Lot Size';
  IF x1 > 500 then x2 = 1; 
    ELSE x2 = 0;
  IF x1 > 500 then x3 = 1;   /* Including a second indicator */
    ELSE x3 = 0;
  x2star = (x1 - 500)*x2;
cards;
  2.57  650
  4.40  340
  4.52  400
  1.39  800
  4.75  300
  3.55  570
  2.49  720
  4.27  480
  4.42  425
  4.10  460
  2.21  675
  2.90  600
;
run;

/* Adding 3 new data points with X values near */
/* the changepoint (for plotting purposes)     */

DATA onemorepoint;
INPUT y x1;
  IF x1 > 500 then x2 = 1; 
    ELSE x2 = 0;
  IF x1 > 500 then x3 = 1;   /* Including a second indicator */
    ELSE x3 = 0;
  x2star = (x1 - 500)*x2;
cards;
. 499.99
. 500
. 500.01
;
run;
DATA rawmater;
SET rawmater onemorepoint;
run;

/* Fitting a discontinuous piecewise regression function with a */
/* changepoint at X = 500:                                      */

PROC REG DATA = rawmater;
  MODEL y = x1 x2star x3;
  /* Recall x2star = (x1 - 500)*x2 here!! */
  OUTPUT out=new2 p=yhat;
run;
/* Sorting the data by x1   */
/* to prepare for plotting: */
PROC SORT DATA = new2;
  BY x1;
run;

/* Plotting the discontinuous piecewise regression function: */

goptions reset=all; 

symbol1 i=join c=black;
symbol2 c=black v=circle;
axis1 label=(angle=90 'Unit Cost');
PROC GPLOT data = new2;
  PLOT yhat*x1 y*x1/ overlay vaxis=axis1;
run;