/**********************************************/
/*                BUNDLES PROGRAM             */
/**********************************************/
/*
AUTHOR: FRANCISCO J. DIAZ, PH.D. 

Department of Biostatistics
The University of Kansas Medical Center
Mail Stop 1026
3901 Rainbow Blvd.
Kansas City, KS 66160, USA 
Phone: (913) 945-7006. 
Email: fdiaz@kumc.edu

*/
/*************************************************/
/*
OBJECTIVE OF THE PROGRAM: TO PLOT BUNDLES FOR LOG-TRANSFORMED
AND NON-TRANSFORMED DATA, AS DESCRIBED IN: 

Diaz, F. J. (1999) Identifying Tail Behavior by Means of 
Residual Quantile Functions. Journal of Computational and 
Graphical Statistics 8, 493-509.

THIS PROGRAM CAN BE FREELY USED, SHARED AND COPIED, PROVIDED THAT ITS 
AUTHOR AND THE ABOVE REFERENCE ARE NOT DELETED.  
©copyright Francisco J Diaz 1999. All rights reserved. 
*/
/*************************************************/

/**************************************************/
/*
INSTRUCTIONS FOR INPUTING YOUR DATA:

THIS SAS/IML PROGRAM ASSUMES THAT YOUR DATA HAS ONLY ONE VARIABLE AND THAT THE DATA 
IS LOCATED IN A TXT OR DAT FILE SAVED IN A FOLDER NAMED "CUANTIL" IN YOUR C 
DRIVE. (THE DRIVE AND NAME OF THE FOLDER CAN BE CHANGED BY DECLARING THEM IN 
THE LIBNAME AND INFILE COMMANDS BELOW.) THE FIRST PART OF THE PROGRAM TAKES THE 
DATA FROM THE ABOVE FILE, SORT THEM FROM LEAST TO GREATEST, AND STORES THEM 
INTO A SAS DATA FILE NAMED "CUANTIL.TEMP". 

ANY OTHER SAS METHOD OF INPUTING THE DATA INTO THIS PROGRAM CAN BE USED. HOWEVER, 
THE BUNDLE IS COMPUTED BY AN IML PROCEDURE THAT ASSUMES THAT THE DATA IS IN A SAS 
DATA FILE CALLED CUANTIL.TEMP. THUS, IF YOU WANT TO USE A DIFFERENT METHOD TO INPUT 
YOUR DATA, PLEASE MAKE SURE TO INSTRUCT SAS TO CONVERT YOUR DATA INTO A SORTED, ONE-COLUMN 
SAS FILE CALLED "CUANTIL.TEMP" BEFORE PROC IML IS EXECUTED. 

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


/*YOUR DATA IS NEXT INPUT*/

OPTIONS LS=100;

DM OUTPUT 'CLEAR';
DM LOG 'CLEAR';

LIBNAME CUANTIL 'c:\CUANTIL';

DATA CUANTIL.TEMP;
INFILE 'c:\cuantil\rain.txt' delimiter=',' ;/*ENTER HERE THE NAME OF 
                                             YOUR DATA FILE*/
INPUT X1  @@;
Run;

PROC SORT DATA=CUANTIL.TEMP;
BY X1;
RUN;
/****************************************************/

/*THE BUNDLE IS NEXT COMPUTED WITH AN IML PROCEDURE*/

PROC IML ;

USE CUANTIL.TEMP;

READ ALL INTO X;
X=X[LOC(X)]; /*ZEROS AND MISSING VALUES ARE DROPPED WITH THIS INSTRUCTION*/
/*****************************************************/
X=LOG(X); /*THIS STATEMENT SHOULD BE DROPPED IF A BUNDLE FOR 
           THE NON-TRANSFORMED DATA IS NEEDED*/
/*****************************************************/
N=NROW(X);
print n;

VP={ 0.1  0.2  0.3  0.4   0.5  0.6  0.7  0.8  0.9 };


AUX=J(1,N,0);
K=1;
I=1;

DO WHILE (K<=N);
AUX[I] = SUM(X<=X[K]);
K=AUX[I]+1;
I=I+1;

END;



AUX=AUX[LOC(AUX)];

CORTE=1;

KI=AUX[1:NROW(AUX)-CORTE];
XKI=X[KI];
XKIMUNO=X[AUX[2:NROW(AUX)-CORTE+1]];



Y=XKI||KI;
CREATE CUANTIL.X FROM X [COLNAME={X1}];
APPEND FROM X;
CREATE CUANTIL.Y FROM Y [COLNAME={Y1 Y2}];
APPEND FROM Y;



DO I=1 TO NCOL(VP);
L1=VP[I]#(N-KI);
L=(INT(L1)=L1);
AUX2=X[(L1+KI)#L+(INT(L1)+1+KI)#(1-L)];



MPBS=MPBS||(AUX2-XKI); 


CUANTIL=CUANTIL//(XKI||(AUX2-XKI)||J(NROW(XKI),1,VP[I]));

END;


CREATE CUANTIL.CUANTIL FROM CUANTIL [COLNAME={XKI BS  P}];
APPEND FROM CUANTIL;


VRM=J((NROW(AUX)-CORTE),1,0);
DO I=1 TO NROW(AUX)-CORTE;
  VRM[I]=SUM((X>XKI[I])#X)/(N-KI[I])-XKI[I];
END;


P=REPEAT(-LOG(1-VP),NROW(XKI),1);

RAMILL=(MPBS/P)||VRM||XKI;

CREATE CUANTIL.RAMILL FROM RAMILL [COLNAME={X1 X2 X3 X4 X5 X6 X7 X8 X9 VRM XKI}];

APPEND FROM RAMILL;


RAMILL2=(MPBS/P)||VRM;
CREATE CUANTIL.RAMILL2 FROM RAMILL2 [COLNAME={X1 X2 X3 X4 X5 X6 X7 X8 X9 VRM }];

APPEND FROM RAMILL2;

CALL SOUND(600,1);
QUIT;
/*********************************************************/


/*THE BUNDLE IS NEXT PLOTTED*/

GOPTIONS   CBACK=WHITE COLORS=(BLACK BLACK BLACK) BORDER FBY=SWISS HBY=1.7;
OPTIONS NODATE;

SYMBOL;
FOOTNOTE;
title;

SYMBOL1 INTERPOL=JOIN VALUE=DOT HEIGHT=0.2 COLOR=black;
SYMBOL2 INTERPOL=JOIN VALUE=DOT HEIGHT=0.2 COLOR=black;
SYMBOL3 INTERPOL=JOIN VALUE=DOT HEIGHT=0.2 COLOR=black;
SYMBOL4 INTERPOL=JOIN VALUE=DOT HEIGHT=0.2 COLOR=black;
SYMBOL5 INTERPOL=JOIN VALUE=DOT HEIGHT=0.2 COLOR=black;
SYMBOL6 INTERPOL=JOIN VALUE=DOT HEIGHT=0.2 COLOR=black;
SYMBOL7 INTERPOL=JOIN VALUE=DOT HEIGHT=0.2 COLOR=black;
SYMBOL8 INTERPOL=JOIN VALUE=DOT HEIGHT=0.2 COLOR=black;
SYMBOL9 INTERPOL=JOIN VALUE=DOT HEIGHT=0.2 COLOR=black;


AXIS1 minor=none major=(number=5) label=(angle=90 f=simplex height=1.7  'Q' f=simplex h=1.02 '[n]' h=1.7 f=simplex 'p'
f=simplex  '(' f=greek
'q' f=simplex ')' f=simplex '/log' f=simplex '(1-p)' h=1.02 '-1' h=1.7 ', p=0.1, 0.2, ..., 0.9');

AXIS2 LABEL=(f=greek H=1.7  'q') MAJOR=NONE MINOR=NONE;


PROC GPLOT DATA=CUANTIL.RAMILL ;
PLOT X1*XKI X2*XKI X3*XKI X4*XKI X5*XKI X6*XKI X7*XKI X8*XKI X9*XKI  /OVERLAY
   VREF=0    LVREF=4      FRAME      VAXIS=AXIS1 HAXIS=AXIS2 ;
RUN;
quit;
/*************************************************************/