/*************************************************************/ 
/*************************************************************/
/******************* MARKOV **********************************/
/*************************************************************/
/***** PROGRAM TO CARRY OUT BOOTSTRAP ESTIMATION *************/
/***** OF CONFIDENCE INTERVALS OF AUTOREGRESSIVE *************/
/***** MARKOV PROCESSES **************************************/
/*************************************************************/
/*************************************************************/
/****** 27 SEPTEMBER 2004 ************************************/
/*************************************************************/

/*************************************************************/
/** REFERENCE: BOOTSTRAP METHODS for MARKOV PROCESSES (2003) */
/** ECONOMETRICA, VOL. 71, NO.4 ******************************/


new;
FORMAT 9, 4;
clear N, D, X, Z, HX, Q, DELX, XGRID, NGRID, NBOOT, LVL, BXCUM, A, B;
OUTPUT FILE = C:\GAUSS\MARKOV\MARKOV1.OUT ON;

@ *********************** NOTATION *********************************  @
@ N = SAMPLE SIZE
  D = NUMBER OF COLUMS OF X
  X = ORIGINAL DATA
  Z = DATA VECTOR OF X and Y
  HX = BANDWIDTH
  Q = ORDER OF MARKOV PROCESS
  DELX = WIDTH OF GRIDCELLS for BOOTSTRAP VALUES OF X - D X 1 VECTOR
  NG = NUMBER OF GRID POINTS IN EACH DIMENSION OF XGRID
  XGRID = GRID OF X POINTS for BOOTSTRAP VALUES OF X - NG^D MATRIX
  NGRID = NUMBER OF ELEMENTS IN XGRID
  NBOOT = NUMBER OF BOOTSTRAP SAMPLES
  LVL = NOMINAL LEVEL OF A TEST
@

/******************************************************************************/
/*********************SECOND ORDER KERNEL WITH SUPPORT [-1, 1] ****************/
/******************************************************************************/

@
  GIVES THE SECOND ORDER KERNEL for ARRAY V WITH SUPPORT [-1, 1] 
@

proc(1) = KERN2(V,K);

@
  if K < 1.5, COMPUTE KERNEL. OTHERWISE, COMPUTE INTEGRAL OF KERNEL
@ 

local TERMTEMP, TERM, TERM1, TERM2, TERM3, TERM4, VTERM4, INDEX, INDEX2, INDEX3, INDEX4, INDEX5,
TMP, TMP2, TMP3, TMP4, TMP5, VTM4, NV, CV, VTM4MEAN, I;

NV = ROWS(V);
CV = COLS(V);
TERM = ZEROS(NV,1);
I = 1;
if K < 1.5;
  TERM = ZEROS(NV,1);
  do until I > NV;
  TERM[I] = (15/16).*((1 - V[I,.]*V[I,.]')^2);
  TERM[I] = TERM[I].*(V[I,.]*V[I,.]' .<= 1);
  I = I + 1;
  endo;
else;
  if CV == 1;
    TERM = ZEROS(NV,1);
    TERM = (15/16)*((1/5)*V.^5 + V - (2/3)*V.^3) + 0.5;
    TERM = TERM.*(ABS(V) .<= 1) + (V .> 1);
  else;
    TMP = V[.,1];
    INDEX = 2;
    do until INDEX > CV;
    TMP = TMP.*V[.,INDEX];
    INDEX = INDEX + 1;
    endo;
    TERM1 = TMP;
    
    TERM2 = V.^4;
    TERM3 = V.^2;
    TERM2 = TERM2.*TERM1;
    TERM3 = TERM3.*TERM1;
    TMP2 = TERM2[.,1];
    TMP3 = TERM3[.,1];
    INDEX2 = 2;
    do until INDEX2 > CV;
    TMP2 = TMP2 + TERM2[.,INDEX2];
    TMP3 = TMP3 + TERM3[.,INDEX2];
    INDEX2 = INDEX2 + 1;
    endo;
    TERM2 = TMP2;
    TERM3 = TMP3; 

    VTERM4 = ZEROS(NV,CV - 1);
    INDEX3 = 1;
    do until INDEX3 > CV - 1;
    VTM4 = ZEROS(NV,CV - 1);
    INDEX4 = INDEX3 + 1;
    do until INDEX4 > CV;
    TMP4 = (V[.,INDEX3].^2).*(V[.,INDEX4].^2);
    TMP4 = TMP4.*TERM1;
    VTM4[.,INDEX3] = TMP4;
    INDEX4 = INDEX4 + 1;
    endo;
    
    VTM4 = SUMC(VTM4');
    VTERM4[.,INDEX3] = VTM4;
    INDEX3 = INDEX3 + 1;
    endo;
    TERM4 = SUMC(VTERM4');
    
    TERM = ZEROS(NV,1);
    TERM = (15/16)*(TERM1 + (1/5)*TERM2 - (2/3)*TERM3 + (2/9)*TERM4) + 0.5;
    do until I > NV;
    TERM[I] = TERM[I].*(V[I,.]*V[I,.]' .<= 1);
    I = I + 1;
    endo;
  endif;
   
endif;

retp(TERM);
endp;

/******************************************************************************/
/********* COMPUTE CONDITIONAL CDF AND ALL NECESSARY DENSITIES ****************/
/******************************************************************************/

proc(2) = FWIGGL(Z1,Z2,Z);

@
  THIS PROCEDURE COMPUTES DENSITY ESTIMATORS GWIGGL, PWIGGL and FWIGGL
  (CONDITIONAL CDF). GWIGGL IS COMPUTED AT POINT Z1. PWIGGL IS COMPUTED
  AT POINTS X[Q+1],...,X[2] and X[Q],...,X[1], WHERE X[.] IS A D-VECTOR. 
  GWIGGL, PWIGGL and FWIGGL ARE NGRID X 1.
@

@
  Z = (N - Q) X D*(Q + 1): MATRIX CONTAINING OBSERVATIONS OF THE STATE
  VARIABLE X and ITS Q LAGGED VALUES, Y. Z IS THE DATA. IMPORTANT: THE FIRST ROWS 
  OF Z SHOULD HAVE THE LOWEST TIME INDEXES, AS WELL AS THE FIRST COLUMNS. 
  THEREFORE Z = X[1:N - Q,.],...,X[Q + 1:N,.] 
  

  Z1 = Q X D VECTOR CONSISTING OF CURRENT and LAGGED VALUES OF THE 
  STATE VARIABLE. THE FIRST ROW OF Z1 SHOULD HAVE THE LOWEST TIME INDEX.

  Z2 = NGRID X D VECTOR OF POSSIBLE VALUES OF STATE VARIABLE IN THE
  NEXT PERIOD. THE TRANSITION PROBABILITY IS OF Z2 CONDTIONAL ON Z1.

@

local INDEX, NV, J, TEMP, DIF, INDEX1, DZTEMP, DYPTEMP, FWIG, PPWIG, GWIG, PWIG,
INDEX2, INDEX3, INDEX4, INDEX5, DYP, DY, DZZ, DZ;

INDEX = 1;
J = D*(INDEX - 1);
do until INDEX > Q;
if INDEX == 1;
  DZ = KERN2((Z[.,(1 + J):(J + D)] - Z1[INDEX,.])/HX,1);
  DY = DZ;
  DYP = ONES(N - Q,1);
else;
  DZ = DZ.*KERN2((Z[.,(1 + J):(J + D)] - Z1[INDEX,.])/HX,1);
  DY = DZ;
  DYP = DYP.*KERN2((Z[.,(1 + J):(J + D)] - Z1[INDEX,.])/HX,1);
endif;
INDEX = INDEX + 1;
endo;

DIF = ZEROS((N - Q)*NGRID,D);
INDEX1 = 1;
do until INDEX1 > (N - Q);
DIF[(1 + (INDEX1 - 1)*NGRID):(INDEX1*NGRID),.] = Z[INDEX1,(D*Q + 1):(D*(Q + 1))] - Z2;
INDEX1 = INDEX1 + 1;
endo;

DZTEMP = ZEROS((N - Q)*NGRID,1);
DYPTEMP = ZEROS((N - Q)*NGRID,1);
INDEX2 = 1;
do until INDEX2 > (N - Q);
INDEX3 = 1;
do until INDEX3 > NGRID;
DZTEMP[(INDEX3 + (INDEX2 - 1)*NGRID),.] = DZ[INDEX2,.];
DYPTEMP[(INDEX3 + (INDEX2 - 1)*NGRID),.] = DYP[INDEX2,.];
INDEX3 = INDEX3 + 1;
endo;
INDEX2 = INDEX2 + 1;
endo;

DZZ = DZTEMP.*KERN2(DIF/HX,2); @ CDF OF Z[.,1:Q] @
DZ = DZTEMP.*KERN2(DIF/HX,1);  @ DENSITY OF Z[.,1],...,Z[.,Q + 1] @
DYP = DYPTEMP.*KERN2(DIF/HX,1); @ DENSITY OF Z[.,2],...,Z[.,Q + 1] @

GWIG = ZEROS(NGRID,N - Q);
FWIG = ZEROS(NGRID,N - Q);
PPWIG = ZEROS(NGRID,N - Q);
INDEX4 = 1;
do until INDEX4 > NGRID; 
INDEX5 = 1;
do until INDEX5 > (N - Q);
GWIG[INDEX4,INDEX5] = DZ[(INDEX4 + (INDEX5 - 1)*NGRID),.];
FWIG[INDEX4,INDEX5] = DZZ[(INDEX4 + (INDEX5 - 1)*NGRID),.];
PPWIG[INDEX4,INDEX5] = DYP[(INDEX4 + (INDEX5 - 1)*NGRID),.];
INDEX5 = INDEX5 + 1;
endo;
INDEX4 = INDEX4 + 1;
endo;

GWIG = MEANC(GWIG');
FWIG = MEANC(FWIG');
PPWIG = MEANC(PPWIG');
PWIG = MEANC(DY);

GWIG = GWIG/(HX^(Q + 1));  @ DENSITY OF Z[.,1],...,Z[.,Q + 1] @
FWIG = FWIG/(HX^Q);  @ CDF AT Z1[Q + 1], DENSITY OF Z[.,1:Q] @
PWIG = PWIG/(HX^Q); @ DENSITY OF Z[.,1],...,Z[.,Q] @
TEMP = PWIG.*(ABS(PWIG) .> 1e-10) + 1e-10*(ABS(PWIG) .<= 1e-10);
PPWIG = PPWIG/(HX^Q);  @ DENSITY OF Z[.,2],...,Z[.,Q + 1] @
FWIG = FWIG./TEMP;

retp(PPWIG,FWIG); @ PPWIG IS THE MARGINAL DENSITY and FWIG IS THE TRANSITION CDF @
endp;

/***********************************************************************/
/************** GENERATE INITIAL BOOTSTRAP OBSERVATION *****************/
/***********************************************************************/

@
  THIS PROCEDURE GENERATES THE INITIAL STATE OF THE BOOTSTRAP MARKOV
  PROCESS BY SAMPLING Q CONSECUTIVE OBSERVATIONS FROM THE ORIGINAL DATA

  ZDAT = N X D MATRIX CONTAINING THE DATA
  Z[1,.] HAS THE LOWEST TIME INDEX
@

proc(1) = XSTART(ZDAT);

local N1, N2, ZVEC;

N1 = CEIL((N - Q + 1)*RNDU(1,1));
N2 = N1 + Q - 1;
ZVEC = ZDAT[N1:N2,.];

retp(ZVEC);
endp;

/***********************************************************************/
/************** GENERATE BOOTSTRAP SAMPLE ******************************/
/***********************************************************************/

proc(1) = BSAMPL(Z,ZDAT);

@
  THIS PROCEDURE GENERATES ALL MARKOV BOOTSTRAP OBSERVATIONS EXCEPT THE 
  FIRST. THE FIRST IS GENERATED IN THE PROCEDURE XSTART

  THE OUTPUT IS A COLUMN VECTOR OF N OBSERVATIONS OF A BOOSTRAP SCALAR X
  
  AS BEFORE, Z = (N - Q) X (D*Q + D): MATRIX CONTAINING OBSERVATIONS OF 
  THE STATE VARIABLE X and ITS Q LAGGED VALUES, Y. Z IS THE DATA. IMPORTANT: 
  THE FIRST ROWS OF Z SHOULD HAVE THE LOWEST TIME INDEXES, AS WELL AS THE 
  FIRST COLUMNS. THEREFORE Z = X[1:N - Q,.],...,X[Q + 1:N,.] 
   
  ZDAT = N X D MATRIX OF ORIGINAL DATA
  
  OUTPUT IS A SERIES OF N BOOTSTRAP OBSERVATIONS
@

local INDEX, Z1, INDEX2, PWIG, FWIG, RN, TEMP, XTRY, N0, N1, Z2, XTMP, X0;

X0 = XSTART(ZDAT);
RN = RNDU(N,1);

INDEX = Q + 1;
do until INDEX > N;
FWIG = ZEROS(NGRID,1);  
    
INDEX2 = 1;
do until INDEX2 > NGRID;
Z1 = X0;
Z2 = XGRID;
{PWIG,FWIG} = FWIGGL(Z1,Z2,Z);
INDEX2 = INDEX2 + 1;
endo;

TEMP = MININDC(ABS(FWIG - RN[INDEX]));

if INDEX == Q + 1;
  XTRY = X0|XGRID[TEMP,.];
else;
  XTRY = XTRY|XGRID[TEMP,.];
endif;

N1 = ROWS(XTRY);
N0 = N1 - Q + 1;
X0 = XTRY[N0:N1,.];
INDEX = INDEX + 1;
endo;

retp(XTRY);
endp;

/*******************************************************************************/
/******** ESTIMATION OF H(.) AND T-STATISTIC USING SAMPLE FROM POPULATION ******/
/*******************************************************************************/

@
  HERE WE ESTIMATE THE STATISTIC OF INTEREST USING THE ORIGINAL SAMPLE   
@

proc(2) = HEST(Z,HZERO);

@ HZERO IS THE VALUE OF H UNDER THE NULL HYPOTHESIS @

local MHAT, HHAT, TEMP, VHAT, SHAT, TSTAT, TMP, CY, VY;

MHAT = MEANC(X);  @ MEAN OF VALUES OF X @

@ ***** CALCULATE HHAT, WHICH IS A FUNCTION OF MHAT - WRITE THE FORMULA ****** @
HHAT = MHAT;   @ ***** WRITE FORMULA HERE ****** @
HHAT = HHAT[1];  @ ***** MODIFY INDEX TO DESIRED ESTIMATOR ***** @

@ ****** WRITE THE FORMULA TO FIND VARIANCE OF N^(1/2)*(HHAT - HZERO) ****** @
VY = MEANC(X[.,1]^2) - (MEANC(X[.,1]))^2;
CY = MEANC(Z[.,1].*Z[.,2]);
VHAT = VY + 2*((N - 1)/N)*CY; @ ***** WRITE FORMULA HERE ****** @

@ COMPUTE THE T-STATISTIC @
SHAT = SQRT(VHAT);    @ ***** MODIFY INDEX TO DESIRED ESTIMATOR ***** @
TSTAT = N^(1/2)*(HHAT - HZERO)/SHAT[1];    @ ***** MODIFY INDEX TO DESIRED ESTIMATOR ***** @

retp(HHAT,TSTAT);
endp;

/*******************************************************************************/
/*******************  ESTIMATION of H(.) USING BOOSTRAP SAMPLE *****************/
/*******************************************************************************/

proc(2) = BHEST(ZBOOT);

local BMHAT, BHHAT, BVHAT, BSHAT, XBOOT, BCY, BVY;

XBOOT = ZBOOT[.,1:D];
BMHAT = MEANC(XBOOT);  @ MEAN OF BOOTSTRAP VALUES OF X @

@ ***** CALCULATE BHHAT, WHICH IS A FUNCTION OF MHAT - WRITE THE FORMULA ****** @
BHHAT = BMHAT;   @ ***** WRITE FORMULA HERE ****** @
BHHAT = BHHAT[1];  @ ***** ADJUST INDEX TO DESIRED ESTIMATOR ****** @

@ ****** WRITE THE FORMULA TO FIND VARIANCE OF N^(1/2)*(BHHAT - HHAT) OF EACH BOOTSTRAP 
SAMPLE - THIS FORMULA SHOULD LOOK EXACTLY LIKE THE FORMULA for VHAT, EXCEPT for THE FACT
THAT NOW BOOTSTRAP SAMPLES ARE BEING USED ****** @
BVY = MEANC(X[.,1]^2) - (MEANC(X[.,1]))^2;
BCY = MEANC(Z[.,1].*Z[.,2]);
BVHAT = BVY + 2*((N - 1)/N)*BCY; @ ***** WRITE FORMULA HERE ****** @

@ COMPUTE THE STANDARD DEVIATION @
BSHAT = SQRT(BVHAT[1]);    @ ***** MODIFY INDEX TO DESIRED BOOTSTRAP ESTIMATOR ***** @

retp(BHHAT,BSHAT);
endp;

/*******************************************************************************/
/**************** BOOTSTRAP COMPUTATION OF CRITICAL VALUE **********************/
/*******************************************************************************/

proc(3) = TBOOT(HHAT,Z,ZDAT,NCRIT);

@
  AS BEFORE, Z = (N - Q) X (D*Q + D): MATRIX CONTAINING OBSERVATIONS OF 
  THE STATE VARIABLE X and ITS Q LAGGED VALUES, Y. Z IS THE DATA. IMPORTANT: 
  THE FIRST ROWS OF Z SHOULD HAVE THE LOWEST TIME INDEXES, AS WELL AS THE 
  FIRST COLUMNS. THEREFORE Z = X[1:N - Q,.],...,X[Q + 1:N,.] 
   
  ZDAT = N X D MATRIX OF ORIGINAL DATA
@

local XBOOT, INDEX, ZBOOT, BHHAT, BSHAT, IBOOT, BTCUM, BTCRIT1, BTCRIT2,
BSCUM, BHCUM;

IBOOT = 1;
BHCUM = ZEROS(NBOOT,1);  @ STORES BOOTSTRAP ESTIMATES OF H @
BTCUM = ZEROS(NBOOT,1);  @ STORES BOOTSTRAP ESTIMATES OF THE T STATISTIC @
BXCUM = ZEROS(N,D*NBOOT);  @ STORES BOOTSTRAP SAMPLES @
BSCUM = ZEROS(NBOOT,1);  @ STORES BOOTSTRAP ESTIMATES OF THE STDEV @

do until IBOOT > NBOOT;

@ GENERATE A BOOTSTRAP SAMPLE @
XBOOT = BSAMPL(Z,ZDAT);  @ THIS IS N X D @

@ GETTING THE DATA ORGANIZED IN THE APPROPRIATE WAY @
INDEX = 1;
ZBOOT = XBOOT[(Q + 1):N,.];
do until INDEX > Q;
ZBOOT = XBOOT[(Q + 1 - INDEX):(N - INDEX),.]~ZBOOT;
INDEX = INDEX + 1;
endo;

{BHHAT,BSHAT} = BHEST(ZBOOT);
BHCUM[IBOOT] = BHHAT; 
BSCUM[IBOOT] = BSHAT;
BXCUM[.,(1 + (IBOOT - 1)*D):IBOOT*D] = XBOOT;
IBOOT = IBOOT + 1;
endo;

@ CALCULATING THE BOOTSTRAP T-STSTISTIC @
BTCUM = N^(1/2)*(BHCUM - HHAT)./BSCUM;  
BTCUM = SORTC(BTCUM,1);

@ FINDING THE BOOTSTRAP CRITICAL VALUES @
BTCRIT1 = SORTC(ABS(BTCUM),1);  @ TWO - SIDED TEST @
BTCRIT1 = BTCRIT1[NCRIT];       @ TWO - SIDED TEST @
BTCRIT2 = BTCUM[NCRIT];         @ ONE - SIDED TEST @

retp(BTCUM,BTCRIT1,BTCRIT2);
endp;


/***********************************************************************/
/********************* CREATING THE GRID *******************************/
/***********************************************************************/

proc(1) = ONESTEP(A,XGRID);

local INDEX, INDEX1, RK, CK, B;

RK = ROWS(A);
CK = COLS(A);
B = ZEROS(RK*NGRID, CK + 1);
INDEX = 1;
do until INDEX > RK;
INDEX1 = 1;
do until INDEX1 > NGRID;
B[(INDEX1 + (INDEX - 1)*NGRID),.] = A[INDEX,.]~XGRID[INDEX1,(CK + 1)];
INDEX1 = INDEX1 + 1;
endo;
INDEX = INDEX + 1;
endo;

retp(B);
endp;

/***********************************************************************/
/********************* MAIN PROGRAM ************************************/
/***********************************************************************/

TSTART = DATE;

@ ****** MODIFY THE FOLLOWING TO DESIRED VALUES ******* @

NBOOT = 100;   @ ***** NUMBER OF BOOTSTRAP ITERATIONS ***** @
Q = 1;        @ ***** ORDER OF MARKOV PROCESS ***** @
NGRID = 50;   @ ***** NUMBER OF BOOTSTRAP GRID POINTS for EACH X[.,I] - TOTAL NUMBER 
              OF GRID POINTS WILL BE NGRID^D ***** @
LVL = 0.05;   @ ***** NOMINAL LEVEL OF THE TEST ***** @
HX = 3;       @ ***** BANDWIDTH ****** @
HZERO = 0;    @ ***** NULL HYPOTHSIS OF H(.) ****** @

@ ***** CHOOSE LEVEL OF TEST ***** @

TCRIT = {1.96, 1.645};  @ ***** ASYMPTOTIC CRITICAL VALUES for 0.05 LEVEL for 
  2 and 1 SIDED TESTS ****** ({2.575, 2.3264} for 0.01 LEVEL 2 and 1 SIDED) ***** @
NCRIT = CEIL(NBOOT*(1 - LVL));  @ BOOTSTRAP CRITICAL VALUE @


/***************** SETTING UP THE DATA **********************************/

@ 
  load THE DATA and SET UP ARRAYS 
  DATA SHOULD BE AN N X D MATRIX OF OBSERVATIONS OF X, I.E. DATA IS ZDAT
@

load DATASET[50,1] = C:\GAUSS\MARKOV\AR1.TXT; @ ***** SPECIFY SIZE OF MATRIX and PATH TO DATA - DATA SHOULD BE IN A .TXT FILE ***** @
X = DATASET;
N = ROWS(X);  @ NUMBER OF OBSERVATIONS @
D = COLS(X);  @ DIMENSION OF VECTOR X @

@
  STANDARDIZE THE VARIABLES
@
X = (X - MEANC(X)')./STDC(X)';

@ 
  SETUP THE VECTOR Z TO BE USED IN THE MARKOV BOOTSTRAP 
@

INDEX1 = 1;
Z = X[(Q + 1):N,.];
do until INDEX1 > Q;
Z = X[(Q + 1 - INDEX1):(N - INDEX1),.]~Z;
INDEX1 = INDEX1 + 1;
endo;

ZDAT = X;


OUTPUT ON;
"";
"SAMPLE SIZE:  ";; N;
"BOOTSTRAP ITERATIONS:  ";; NBOOT;
"NUMBER OF GRID POINTS:  ";; NGRID;
"ORDER OF MARKOV PROCESS:  ";; Q;
"";
OUTPUT OFF;

/**************** FINDING THE SAMPLE ESTIMATOR OF H(.) ********************/
{HHAT,TSTAT} = HEST(Z,HZERO);

/************** CREATING THE GRID for THE BOOTSTRAP POINTS: XGRID ********/
XMIN = MINC(X);
XMAX = MAXC(X);
XSTEP = ((XMAX) - (XMIN))/(NGRID - 1);
XGRID = ZEROS(NGRID,D);
I = 1;
do until I > D;
XGRID[.,I] = SEQA(XMIN[I],XSTEP[I],NGRID);
I = I + 1;
endo;

if D == 1;
XGRID = XGRID;
else;
A = XGRID[.,1];
II = 1;
do until II > D - 1;
{B} = ONESTEP(A,XGRID);
A = B;
II = II + 1;
endo;
XGRID = B;
endif;

NGRID = NGRID^D;

if NBOOT > 0;

/****** FINDING THE BOOTSTRAP T-STATISTICS and CRITICAL VALUES ************/
{BTCUM,BCRIT1,BCRIT2} = TBOOT(HHAT,Z,ZDAT,NCRIT);

endif;

RTIME = ETHSEC(TSTART,DATE)/6000;  @ RUNNING TIME IN MINUTES @
RTIMEH =ETHSEC(TSTART,DATE)/360000;  @ RUNNING TIME IN HOURS @
OUTPUT ON;
"";
"ASYMPTOTIC CRITICAL VALUES:   ";; TCRIT;
"BOOTSTRAP CRITICAL VALUES:    ";; BCRIT1;; BCRIT2;
"BOOTSTRAP DISTRIBUTION OF T-STATISTIC    ";;  BTCUM;
"RUNNING TIME IN MINUTES    ";; RTIME;
"RUNNING TIME IN HOURS    ";; RTIMEH;
OUTPUT OFF;