/*
  MAIN PROGRAM TO IMPLEMENT THE GENERALIZED MOMENTS TEST OF THE
  PROPORTIONAL HAZARDS MODEL.  K IS THE DIMENSION OF BETA.  L IS
  THE DIMENSION OF THE TEST STATISTIC.

  REFERENCE:  J.L. HOROWITZ AND G.R. NEUMANN, A GENERALIZED MOMENTS
              SPECIFICATION TEST OF THE PROPORTIONAL HAZARDS MODEL,
              JOURNAL OF THE AMERICAN STATISTICAL ASSOCIATION, 87,
              234-240, 1992.

  ***SPECIFY THE DATA SET***
*/;
  DATASET = "DIME"; @ ***** SPECIFY PATH TO DATA ***** @
  OPEN DSET =^DATASET;
  OUTPUT FILE=LMTEST.OUT RESET;  @ ***** SPECIFY PATH TO OUTPUT FILE ***** @
  TSTART = TIME;
  "DATA SET:  ";; DATASET;
  N = ROWSF(DSET);
  "";
  "THE SIZE OF THE DATA SET IS N =  ";; N;
  "";
/*
  ***IDENTIFY COLUMNS CONTAINING THE DEPENDENT AND CENSORING VARIABLES***
*/;
  DEP = 1;
  CEN = 5;
/*
  ***SET THE NUMBER OF INDEPENDENT VARIABLES***
*/;
  K = 3;
  "THE NUMBER OF INDEPENDENT VARIABLES IS:  ";; K;
  "";
/*
  SET THE DIMENSION OF THE SCORE
*/;
  L = 1;
/*
  ***SPECIFY THE PARAMETER VALUES OBTAINED FROM PMLE***
*/;
@   LET B[2,1] = 2.703741  -0.4999706;   UNCENSORED STRIKE DATA @
   LET B[3,1] =  -0.0069  0.0565  -0.3158;  @ DIME DATA  @

  "";
  "THE PARAMETER VECTOR IS  ";; B';
  "";
  OUTPUT OFF;
  YY = READR(DSET,N);
  CLOSE(DSET);
  Y = SORTC(YY,DEP);
  "SORTING FINISHED  ";
/*
  IDENTIFY THE DEPENDENT VARIABLE
*/;
  T = Y[.,DEP];
/*
  ***IDENTIFY THE INDEPENDENT VARIABLES***
*/;
  X = Y[.,2:4];
/*
  ***IDENTIFY THE INDICATOR OF UNCENSORED OBSERVATIONS***
  ***SET DEL = ONES(N,1) IF THERE IS NO CENSORING.***
  ***THIS IS THE LAST CHANGE.***
*/;
  DEL = 1- Y[.,CEN];
  VRS = T~X~DEL;
  MNS = MEANC(VRS);
  STD = STDC(VRS);
  OUTPUT ON;
  "MEAN AND STD. DEVIATION OF T,X,DEL: ";; MNS~STD;
  OUTPUT OFF;
  CLEAR Y, YY, VRS, MNS, STD;
  E = EHAT(X,T,B,N);
  "FINISHED WITH EHAT  ";
  UHAT = (EXP(X*B)).*CUMSUMC(DEL./(N*E));
@  UHAT; @
  OUTPUT ON;
  "THE PROPORTION OF UNCENSORED OBSERVATIONS IS:  ";; MEANC(DEL);
  "HERE IS THE MEAN OF DEL - UHAT.  ALWAYS 0 IN P.H.  ";;  MEANC(DEL - UHAT);
  "HERE IS THE STANDARD DEVIATION OF UHAT - DEL:  ";; STDC(UHAT - DEL);
  "";
  OUTPUT OFF;
  S = SCR(X,T,DEL,B,UHAT,N);
  "FINISHED WITH SCR  ";
  EX = EXHAT(X,T,B,N,K);
  "FINISHED WITH EXHAT  ";
  VB = VBB(X,T,DEL,B,N,K,E,EX);
  "FINISHED WITH VBB ";
  OUTPUT ON;
  "";
  "THE VBB MATRIX FOLLOWS";
  "";
  VB;
 "";
@  VB = 0;@
  EI = EINT(T,DEL,E,N);
  OUTPUT OFF;
  "FINISHED WITH EINT  ";
  PH = PHIN(X,T,DEL,B,K,N,EX,E);
  "FINISHED WITH PHI  ";
  SSU = SU(X,T,DEL,UHAT,N,L);
  "FINISHED WITH SSU  ";
  OUTPUT ON;
  "";
  "THE SCORE FUNCTION IS:"
  "";
  SCORE = SQRT(N)*MEANC(S);
  SCORE';
  "";
/*
  COMPUTE THE VARIANCE USING THE METHOD THAT GUARANTEES IT IS POSITIVE
*/;
  PHINT = SSU'(PH'VB.*EXP(X*B))/N;
/*
  PHINT HAS DIMENSION (L X K).
*/;
  SMC = CPROC(X,T,DEL,B,K,N,L);
  CTRM = (PHINT*SMC')';
  CLEAR SMC, PHINT;
/*
  SMC RETURNS A (N X K) MATRIX THAT CONTAINS THE SUM OF THE C'S FOR
  EACH K AND N.
*/;
  ATRM = APROC(X,T,DEL,SSU,B,K,N,L);
/*
  APROC RETURNS A (N X L) MATRIX CONTAINING THE INTEGRAL OF THE A TERMS
  FOR THE VARIANCE.
*/;
  TRM =   S + CTRM + ATRM;
  MT = MEANC(TRM);
  SGM = (TRM - MT)'(TRM - MT)/N;
  "THE SIGMA MATRIX FOLLOWS";
  SGM;
  "";
  CHISQ = SCORE'INVPD(SGM)*SCORE;
  TEND=TIME;
ELAPSE=(TEND[1,1]-TSTART[1,1])*60*60 + (TEND[2,1]-TSTART[2,1])*60
  + (TEND[3,1]-TSTART[3,1]) + (TEND[1,1]-TSTART[1,1])/100;

  "";
  "THE CHI SQUARE STATISTIC IS: ";; CHISQ;
  "";
/*
  COMPUTE THE MEAN CORRECTION AND CORRECTED CHI SQUARE STATISTIC
*/;
  VP = VB*PH;
  PVP = SUMC(PH.*VP);
  CORR = -SSU'(EXP(2*X*B).*(-EI + PVP))/(N^1.5);
  CSCORE = SCORE + CORR/2;
  CORRCHI = CSCORE'INVPD(SGM)*CSCORE;
  "";
  "THE CORRECTED CHI SQUARE STATISTIC IS: ";; CORRCHI;
  "";
  "PROBABILITY OF A LARGER CHI SQUARE VALUE:  ";; CDFCHIC(CORRCHI,L);
  "";
  "TOTAL RUNNING TIME IN SECONDS:  ";;ELAPSE;
  CLEAR T,X,DEL,E,EX,VB,EI,PH,SG;
  OUTPUT OFF;
  END;

  PROC APROC(X,T,DEL,SSU,B,K,N,L);
/* THIS PROCEDURE DOES THE INTEGRAL OF THE SUM OF THE THREE A'S
*/;
  LOCAL MD,D,TMP,TRM1,TRM2,EXB,DQ,DD,R1,R2,S1,S2;
 @ COMPUTE TERM ONE@;
  EXB = EXP(X*B);
  TRM1 = SUMC(SSU.*EXB)/N;
  TRM1 = TRM1';
  MD = MEANC(EXB);
  TRM1 = DEL.*TRM1/MD;
  D = REV(CUMSUMC(REV(EXB)))/N;
  TRM1 = TRM1 - DEL.*CUMSUMC(SSU.*EXB./D)/N;
@
  NOW COMPUTE THE SECOND TERM
@
  DQ = -DEL/N;
  R1 = CUMSUMC(DQ./(D.*D));
  CLEAR DQ;
  DD = -EXB/N;
  R2 = CUMSUMC(DD./(D.*D));
  CLEAR D,DD;
  S1 = EXB.*CUMSUMC(SSU.*EXB.*R1)/N;
  S1 = S1 - DEL.*CUMSUMC(SSU.*EXB.*R2)/N;
  S2 = REV(CUMSUMC(REV(SSU.*EXB)))/N;
  S2 = S2[2:N,.]|ZEROS(1,L);
  S2 = S2.*(EXB.*R1 - DEL.*R2);
  TRM2 = S1 + S2;
  RETP(TRM1 + TRM2);
  CLEAR TRM1,TRM2,EXB,R1,R2,S1,S2;
  ENDP;

  PROC CPROC(X,T,DEL,B,K,N,L);
/*THIS PROC DOES THE SUM OF THE C'S IN THE VARIANCE.  IT RETURNS
  A (N X K) MATRIX.
*/;
  LOCAL EXB,D,DD,DDX,DQ,DX,TRM;
  EXB = EXP(X*B);
  DX = REV(CUMSUMC(REV(X.*EXB)))/N;
  TRM = DEL.*X - DEL.*DX[1,.]/MEANC(EXB);
  D = REV(CUMSUMC(REV(EXB)))/N;
  DD = -EXB/N;
  DDX = -X.*EXB/N;
  TRM = TRM - DEL.*CUMSUMC(DDX./D);
  TRM = TRM + DEL.*CUMSUMC(DX.*DD./(D.*D));
  DQ = -DEL/N;
  TRM = TRM + X.*EXB.*CUMSUMC(DQ./D);
  TRM = TRM - EXB.*CUMSUMC(DX.*DQ./(D.*D));
  RETP(TRM);
  CLEAR EXB,D,DD,DDX,DQ,DX,TRM;
  ENDP;

  PROC EHAT(X,T,B,N);
/*
  PROCEDURE TO EVALUATE EHAT(S) AS DEFINED BY TSIATIS.  RETURNS VECTOR EHAT(T)
  CONTAINING EHAT EVALUATED AT EACH T IN THE DATA SET.  DATA MUST BE SORTED
  IN ORDER OF INCREASING T.
*/;
  LOCAL PSI;
  PSI = REV(EXP(X*B));
  RETP(REV(CUMSUMC(PSI)/N));
  CLEAR PSI;
  ENDP;

  PROC EINT(T,DEL,EH,N);
/*
  PROCEDURE TO EVALUATE THE INTEGRAL OF E IN SIGMA(2,2).  RETURNS A
  (N X 1) MATRIX CONSISTING OF THE INTEGRAL UP TO EACH VALUE OF T
  IN THE DATA SET.  THE INTEGRAL UP TO THE MIN OF TWO T VALUES MUST
  BE COMPUTED IN A SEPARATE PROGRAM.  DATA MUST BE SORTED IN ORDER
  OF INCREASING T.
*/;
  RETP(-CUMSUMC(DEL./(N*EH.*EH)));
  ENDP;

  PROC EXHAT(X,T,B,N,K);
/*
  PROCEDURE TO EVALUATE EXHAT AS DEFINED BY TSIATIS.  RETURNS (N X K) VECTOR
  EXH CONTAINING EXH EVALUATED AT EACH T IN THE DATA SET.  DATA MUST BE SORTED
  IN ORDER OF INCREASING T.
*/;
  LOCAL PSI;
  PSI = REV(X.*EXP(X*B));
  RETP(REV(CUMSUMC(PSI)/N));
  CLEAR PSI;
  ENDP;

  PROC PHIN(X,T,DEL,B,K,N,EX,E);
/*
  PROCEDURE TO EVALUATE THE FUNCTION PHI IN THE COVARIANCE OF FUNCTION
  OF THE INTEGRATED HAZARD.  DATA MUST BE SORTED IN ORDER OF INCREASING
  T.  RETURNS (K X N) DIMENSIONAL MATRIX.  EACH COLUMN OF THIS MATRIX IS
  THE K-VECTOR PHI EVALUATED AT EACH POINT (X,T) IN THE DATA SET.  K IS
  THE DIMENSION OF B.
*/;
  LOCAL PH,TERM1,TERM2;
  TERM1 = X.*CUMSUMC(DEL./(N*E));
  TERM2 = CUMSUMC(DEL.*EX./(N.*E.*E));
  RETP(TERM1' - TERM2');
  CLEAR PH, TERM1,TERM2;
  ENDP;

  PROC SCR(X,T,DEL,B,UHAT,N);
/*
  PROCEDURE TO EVALUATE THE SCORE FUNCTION FOR THE LM TEST OF PROPORTIONAL
  HAZARDS.  RETURNS A (N X K) MATRIX.  THE DATA MUST BE SORTED IN
  ORDER OF INCREASING T.
*/;
  LOCAL S1;

@ SCORE BASED ON MOMENTS OF EXP(-UHAT)@;
  S1 = (1 + DEL).*EXP(-UHAT) - 1;
  RETP(S1);
  CLEAR S1;
  ENDP;

  PROC SU(X,T,DEL,UHAT,N,L);
/*
  PROCEDURE TO EVALUATE THE DERIVATIVE OF THE SCORE WITH RESPECT TO
  U.  RETURNS A (N X L) MATRIX, WHERE L IS THE DIMENSION OF THE LM
  TEST STATISTIC.
*/;
  LOCAL S1;
 @SCORE BASED ON MOMENTS OF EXP(-UHAT)@ ;
  S1 = -(1 + DEL).*EXP(-UHAT);
  RETP(S1);
  CLEAR S1;
  ENDP;

  PROC VBB(X,T,DEL,B,N,K,E,EX);
/*
  PROCEDURE TO EVALUATE THE COVARIANCE MATRIX OF B IN PROPORTIONAL HAZARDS.
  RETURNS THE (K X K) MATRIX VBB.  DATA MUST BE SORTED IN ORDER OF INCREASING
  T.
*/;
  LOCAL I, J, VB, VXX;
  LET I = 1;
  VXX = (X[N,.])'X[N,.]*EXP(X[N,.]*B)/N;
  VB = DEL[N,1]*(EX[N,.]'*EX[N,.]./(E[N,.]*E[N,.]) - VXX./E[N,.]);
  DO UNTIL I == N;
  J = N - I;
  VXX = VXX + (X[J,.])'X[J,.]*EXP(X[J,.]*B)/N;
  VB = VB + DEL[J,1]*(EX[J,.]'*EX[J,.]./(E[J,.]*E[J,.]) - VXX./E[J,.]);
  I = I + 1;
  ENDO;
  RETP(INVPD(-VB/N));
  CLEAR VXX;
  ENDP;