SEVERITY Utility Functions

/*--------------------------------------------------------------

                    SAS Sample Library

        Name: svrtutil.sas
 Description: Program defining utility functions for use with
              SEVERITY.
       Title: SEVERITY Utility Functions
     Product: SAS/ETS Software
        Keys: fitting continuous distributions
        PROC: SEVERITY
       Notes: If you run this sample program without any modification, then
              the Sasuser.Svrtdist library contains functions and subroutines
              that are identical to those in the Sashelp.Svrtdist library.

--------------------------------------------------------------*/

proc fcmp outlib=sasuser.svrtdist.utility;
     /*======================================================================
      *
      * Arguments:
      * dim      - Dimension of the input array 'x'.
      * x        - Input array that contains the unique values in the sample.
      * nx       - Input array that contains the number of times each value in
      *            array 'x' appears in the sample.
      * nRaw     - Number of raw moments requested. Upon return, 'raw' would
      *            contain first nRaw moments.
      * raw      - An array that would contain requested raw moments upon
      *            return.
      *            You need not allocate any space prior to this call.
      *======================================================================*/
     subroutine SVRTUTIL_RAWMOMENTS(dim, x[*], nx[*], nRaw, raw[*]);
         outargs raw;
         if (dim < 1 or nRaw < 1) then return;

         do j=1 to nRaw;
            raw[j] = 0;
         end;
         n = 0;
         do i=1 to dim;
             if (not(missing(x[i]) or nx[i] <= 0)) then do;
                 t1 = x[i];
                 do j=1 to nRaw;
                    raw[j] = raw[j] + nx[i] * t1;
                    t1 = t1 * x[i];
                 end;
                 n = n + nx[i];
             end;
         end;

         if (n <= 0) then do;
             do j=1 to nRaw;
                raw[j] = .;
             end;
         end;
         else do;
             do j=1 to nRaw;
                raw[j] = raw[j] / n;
             end;
         end;
     endsub;

     /*======================================================================
      * Computes and returns the requested percentile of the given sample
      * given an estimate of its empirical distribution function (EDF).
      * This is an empirical estimate of the percentile.
      *
      * Arguments:
      * p        - Requested percentile expressed as a fraction in the range
      *            [0,1]. Function returns 100*p th empirical percentile.
      * dim      - Dimension of the input array 'x'.
      * x        - Input array that contains the values in the sample. This
      *            array must be sorted in ascending order. It is assumed that
      *            x contains distinct values.
      * F        - Input array that contains an estimate of the EDF. F[i] is
      *            expected to contain an empirical estimate of Pr(X <= x[i]).
      * Ftype    - Type of estimator that was used to compute F
      *            If 1, standard estimator
      *            If 2, a product-limit estimator (Kaplan-Meier)
      *            If 3, Turnbull's estimator
      *            Types 1 & 2 implies individual data, whereas Type 3 implies
      *            grouped (interval-censored) data.
      *            Further, type 1 implies no censoring or truncation
      *            was used and type 2 implies one form of censoring
      *            (left or right, but not both) was used with optional
      *            truncation.
      *            For type 1, smoothed estimate is returned as described in
      *            Klugman et al. (1998), p. 35.
      *            For type 2, x and F are assumed to be provided at individual
      *            data points such that EDF is constant and is equal to F[i]
      *            for the semi-closed interval [x[i],x[i+1]).
      *            For grouped data, x and F are for the disjoint intervals
      *            such that F is constant between intervals and the increase
      *            in F from the left endpoint of the interval to the right is
      *            arbitrary. This function assumes that F is a linear function
      *            within an interval.
      *            Note that for some intervals, left- and right-
      *            endpoints might be identical to indicate a point-interval,
      *            which allows for a mix of individual and grouped data.
      *======================================================================*/
     function SVRTUTIL_PERCENTILE(p, dim, x[*], F[*], Ftype);
         eps = 2.220446E-16; /* constant('MACEPS') */
         if (dim < 1 or MISSING(p) or p < 0 or p > 1+eps) then return (.);

         if (Ftype=1) then do;
            /* ungrouped data (no censoring/truncation) */
            dimp1 = dim+1;
            if (p >= 1/dimp1-eps) then do;
                if (p <= dim/dimp1+eps) then do;
                    /* return smoothed estimate as per "Loss Models" by
                       Klugman et. al. (1998) pg. 35 */
                    h = dimp1*p;
                    g = floor(h);
                    h = h - g;
                    if (g >= dim) then
                        return (x[g]);
                    return (x[g]*(1-h)+x[g+1]*h);
                end;
                else do;
                    return (x[dim]);
                end;
            end;
            else do;
                return (max(0,(x[1]-eps)/2)); /* assume F is 0 in [0,x[1]) */
            end;
         end;
         else do;
             if (Ftype=2) then do;
                 /* ungrouped data (truncation and one form of censoring).
                    Kaplan-Meier estimate */
                 if (p <= F[1]+eps) then do;
                     if (p < F[1]-eps) then /* p < F[1] */
                         return (max(0,(x[1]-eps)/2)); /* F is 0 in [0,x[1]) */
                     else do; /* p = F[1] */
                         if (dim > 1) then
                            return ((x[1]+x[2])/2);
                         else
                            return (x[1]);
                     end;
                 end;

                 if (p > F[dim]+eps) then
                     return (x[dim]); /* assume F is flat after last point */

                 ileft = 2;
                 iright = dim;
                 i = int((dim+2)/2);
                 do while (i >= 2 and i <= dim);
                     if (p < F[i]-eps) then do; /* p < F[i] */
                         if (p > F[i-1]+eps) then /* F[i-1] < y < F[i] */
                             return (x[i]); /* p is along the step; pctl=x[i] */
                         if (p >= F[i-1]-eps) then /* p = F[i-1] */
                             return ((x[i-1]+x[i])/2);
                         /* p < F[i-1]; move search to left half of
                            current interval */
                         iright = i-1;
                     end;
                     else do;
                         if (p <= F[i]+eps) then do; /* p = F[i] */
                             if (i < dim) then
                                 return ((x[i]+x[i+1])/2);
                             else
                                 return (x[i]);
                         end;
                         /* p > F[i]; move search to right half of
                            current interval */
                         ileft = i+1;
                     end;
                     i = int((ileft+iright)/2);
                 end;
             end;
             else do;
                 /* grouped data (Turnbull's estimator) */
                 if (p <= F[1]+eps) then do;
                     if (p < F[1]-eps) then /* p < F[1] */
                         return (max(0,(x[1]-eps)/2)); /* F is 0 in [0,x[1]) */
                     else /* p = F[1] */
                        return (x[1]);
                 end;

                 if (p > F[dim]+eps) then
                     return (x[dim]); /* assume F is flat after last point */

                 ileft = 1;
                 iright = dim/2;
                 i = 2*int((ileft+iright)/2)-1;
                 do while (i >= 1 and i <= dim-1);
                    if (p < F[i]-eps) then do;
                       /* p < F[i]; move search to left half of
                          current interval */
                       iright = int((i+1)/2) - 1;
                    end;
                    else do;
                        /* p >= F[i] */
                        if (p <= F[i]+eps) then do; /* p = F[i] */
                            if (F[i+1] <= F[i]+eps) then do;
                                return ((x[i]+x[i+1])/2);
                            end;
                            else do;
                                if (i > 1) then
                                    return ((x[i-1]+x[i])/2);
                                else
                                    return (x[1]/2);
                            end;
                        end;
                        /* p > F[i] */
                        if (p < F[i+1]-eps) then do; /* F[i] < p < F[i+1] */
                            return (x[i]+(p-F[i])*(x[i+1]-x[i])/(F[i+1]-F[i]));
                        end;
                        /* p >= F[i+1] */
                        if (p <= F[i+1]+eps) then do; /* p = F[i+1] */
                            if (i < dim-1) then
                                return ((x[i+1]+x[i+2])/2);
                            else
                                return (x[i+1]);
                        end;
                        /* p > F[i+1]; move search to right half of
                           current interval */
                        ileft = int((i+1)/2) + 1;
                    end;
                    i = 2*int((ileft+iright)/2)-1;
                 end;
             end;
         end;
     endsub;

     /* Utility subroutine used by SVRTUTIL_SORT.
      * Given two adjacent partitions, with entries numbered a..(b-1) and
      * b..(c-1), it swaps the smaller partition with the outer entries of
      * the larger. a, b, c should be 0-based indices. */
     subroutine svrtutil_zs_vecswp(x[*], a, b, c);
       outargs x;
       __a = a; __b = b; __c = c;

       __n = MIN((__b - __a), (__c - __b));
       __c = __c - __n;
       do while(__n > 0);
          __i = __a+1; __j = __c+1;
          __t = x[__i]; x[__i] = x[__j]; x[__j] = __t; /* zs_swp */
          __n = __n - 1;
          __a = __i;
          __c = __j;
       end;
     endsub;

     /*======================================================================
      * Sorts the given input sample of numeric values in ascending or
      * descending order.
      *
      * Arguments:
      * dim          - Dimension of the input array 'x'. This array will be
      *                sorted in-place.
      * x            - Input array that contains the values in the sample.
      * isDescending - Desired order. If 0, then the array is sorted in
      *                ascending order. If non-zero, then the array is sorted
      *                in descending order.
      *======================================================================*/
     subroutine SVRTUTIL_SORT(dim, x[*], isDescending);
       outargs x;

       array _st_base[64] /nosymbols;
       array _st_num[64] /nosymbols;
       array _st_mth[64] /nosymbols;
       array _med[13] /nosymbols;

       if (isDescending) then do;
          /* change signs so ascending sort code can be used */
          do _n=1 to dim;
             x[_n] = -x[_n];
          end;
       end;

       _n = dim;
       _a = 0;
       _st_n = 0;
       _mth = 1; /* 1: sort precheck, 2: swap-limited insertion sort */

       do while(1);
          bigloop_begin:;
          /* process section [_a, _a+_n) */

          /* insertion sort on small (or non-random) arrays */
          if (_n < 8 or _mth ne 0) then do;
             if (_n < 8) then do;
                _s = 8*8;
             end;
             else do;
                if (_mth = 1) then
                   _s = 1;
                else
                   _s = _n;
             end;
             do _m = _a+1 to _a+_n-1;
                do _l = _m to _a+1 by -1;
                   /* zs_cmp */
                   if (not(x[_l] > x[_l+1])) then do; /* equal or smaller */
                      go to loopend1;
                   end;
                   else do; /* greater */
                      _j = _l+1;
                      __t = x[_l]; x[_l] = x[_j]; x[_j] = __t; /* zs_swp */
                      _s = _s - 1;
                      if (_s = 0) then go to _punt_isort;
                   end;
                end;
                loopend1:;
             end;

             if (_st_n = 0) then
                go to bigloop_done;

             _st_n = _st_n - 1;
             _a = _st_base[_st_n+1];
             _n = _st_num[_st_n+1];
             _mth = _st_mth[_st_n+1];
             go to bigloop_begin;
          end;

          _punt_isort:;
          /* choose partition value */
          _l = _a; _m = _a + int(_n/2); _h = _a + _n-1;
          _med[1] = _l; _med[2] = _m; _med[3] = _h;
          _v = 2;
          if (_n > 40) then do;
             /* pseudo-median of 9 */
             _s = int(_n/8);
             _med[1] = _l;       _med[2] = _l+_s;   _med[3] = _l+2*_s;
             _med[4] = _m - _s;   _med[5] = _m;     _med[6] = _m + _s;
             _med[7] = _h - 2*_s; _med[8] = _h - _s; _med[9] = _h;
             _v = 8;
          end;

          do _s = 0 to _v-1 by 3;
             _l = _med[_s+1]; _m = _med[_s+2]; _h = _med[_s+3];

             _v = _v + 1;

             /* set _med[_v+1] to the index of the median of values at
                indices _l, _m, _h */
             /* zs_cmp */
             if (x[_l+1] = x[_m+1]) then do;
                _med[_v+1] = _m; go to loopend2;
             end;
             else do;
                if (x[_l+1] > x[_m+1]) then do;
                    __t = _l; _l = _m; _m = __t;
                end;
                /* zs_cmp */
                if (not(x[_m+1] > x[_h+1])) then do;
                    /* equal or smaller */
                    _med[_v+1] = _m; go to loopend2;
                end;
                else do;
                    /* zs_cmp */
                    if (not(x[_l+1] > x[_h+1])) then do;
                        /* equal or smaller */
                        _med[_v+1] = _h; go to loopend2;
                    end;
                    else do;
                        _med[_v+1] = _l;
                    end;
                end;
             end;
          end;
          loopend2:;

          _m = _med[_v+1];

          _orig_a = _a; _v = _a;
          _i = _v+1; _j = _m+1;
          __t = x[_i]; x[_i] = x[_j]; x[_j] = __t; /* zs_swp */
          _c = _a + _n-1; _d = _c;
          _a = _a + 1; _b = _a; /* skip over partition value */
          _mth = 2;

          do while(1);
             do while (_b <= _c);
                /* zs_cmp */
                if (x[_b+1] > x[_v+1]) then do;
                    go to zs_b_high;
                end;
                else do;
                    if (x[_b+1] = x[_v+1]) then do;
                        /* equal */
                        _i = _a+1; _j = _b+1;
                        __t = x[_i]; x[_i] = x[_j]; x[_j] = __t; /* zs_swp */
                        _a = _a + 1;
                    end;
                    /* zs_b_low: */
                    _b = _b + 1;
                end;
             end;
             go to loopend3;

             zs_b_high:;
             do while (_b <= _c);
                /* zs_cmp */
                if (x[_c+1] < x[_v+1]) then do;
                    go to zs_c_low;
                end;
                else do;
                    if (x[_c+1] = x[_v+1]) then do;
                        /* equal */
                        _i = _c+1; _j = _d+1;
                        __t = x[_i]; x[_i] = x[_j]; x[_j] = __t; /* zs_swp */
                        _d = _d - 1;
                    end;
                    /* zs_c_high: */
                    _c = _c - 1;
                end;
             end;
             go to loopend3;

             zs_c_low:;
             _i = _b+1; _j = _c+1;
             __t = x[_i]; x[_i] = x[_j]; x[_j] = __t; /* zs_swp */
             _b = _b + 1;
             _c = _c - 1;
             _mth = 0;
          end;
          loopend3:;

          /* swap equal items back to middle */
          call svrtutil_zs_vecswp(x, _orig_a, _a, _b);
          call svrtutil_zs_vecswp(x, _c+1, _d+1, _orig_a + _n);

          /* push larger partition onto stack */
          _d = _d - _c; _c = _orig_a+_n-_d;
          _b = _b - _a; _a = _orig_a; _n = _b;
          if (_d < _b) then do;
             _a = _c;
             _n = _d;
             _c = _orig_a;
             _d = _b;
          end;
          _st_base[_st_n+1] = _c;
          _st_num[_st_n+1] = _d;
          _st_mth[_st_n+1] = _mth;
          _st_n = _st_n + 1;
       end;

       bigloop_done:;
       if (isDescending) then do;
          /* change signs back to original signs */
          do _n=1 to dim;
             x[_n] = -x[_n];
          end;
       end;
     endsub;

     /* Utility function used by SVRTUTIL_HILLCUTOFF function, to implement
      * one step of the Hill's algorithm. */
     function svrtutil_hillutil_minMSEindex(n,x[*],b,m,samp[*],mse[*],Qopt);
        outargs samp, mse, Qopt;

        do k=1 to m-1;
            mse[k] = 0;
        end;
        do ib=1 to b;
            /* form a bootstrap sample (all elements of x are assumed to
               non-missing and positive) */
            do is = 1 to m;
                ix = floor(rand('uniform')*n)+1;
                samp[is] = x[ix];
            end;
            /* sort the sample */
            call svrtutil_sort(m, samp, 0);
            /* compute mse for each k */
            do k=1 to m-1;
                cvGamma = 0; /* control variate for gamma estimate */
                gamma = 0; /* Hill's tail index estimate for given k */
                _t1 = log(samp[m-k]);
                do i=1 to k;
                    _t2 = log(samp[m-i+1]) - _t1;
                    cvGamma = cvGamma + _t2**2;
                    gamma = gamma + _t2;
                end;
                cvGamma = cvGamma/k;
                gamma = gamma/k;
                _mse = (cvGamma - 2*gamma**2)**2;

                mse[k] = mse[k] + _mse;
            end;
        end;
        /* find k that minimizes MSE averaged over all bootstrap samples */
        _kopt = 1;
        Qopt = mse[1]/b;
        do k=2 to m-1;
            _Q = mse[k]/b;
            if(_Q < Qopt) then do;
                _kopt = k;
                Qopt = _Q;
            end;
        end;
        return(_kopt);
     endsub;

     /*======================================================================
      * Returns the tail cutoff point using the Hill algorithm described
      * in the following paper:
      * Danielsson, J., L. De Haan, L. Peng, and C. G. de Vries. 2001.
      * "Using a Bootstrap Method to Choose the Sample Fraction in Tail Index
      * Estimation." Journal of Multivariate Analysis 76:226-248.
      *
      * Arguments:
      * n        - Number of severity values in the array x. Must be > 5.
      * x        - Array of severity values.
      * b        - Number of bootstrap samples taken to compute the MSE.
      *            If < 10, default value of 50 is used.
      * s        - Maximum number of values to try for searching optimal sample
      *            size (m1) to be used for step 1 of the algorithm. If < 2,
      *            default of 10 will be used. Equally spaced values in the
      *            range [n**0.75,n) are tried.
      * status   - Upon return, this argument contains 0 or 1; 0 indicates
      *            that the algorithm computed the cutoff point successfully.
      *            1 indicates that the algorithm could not compute the cutoff
      *            point successfully, in which case, 5th largest value of the
      *            input sample is returned as the approximate cutoff point.
      *======================================================================*/
     function SVRTUTIL_HILLCUTOFF(n, x[*], b, s, status);
        outargs status;

        status = .;

        _n = int(n);
        if (_n < 6) then
            /* argument error */
            return (.A);

        /* form an internal copy of input array that contains only the
           non-missing, positive values */
        array _x[1] /nosymbols;
        call dynamic_array(_x, _n);
        nvalid = 0;
        do ix=1 to _n;
            if (not(missing(x[ix])) and x[ix] > 0) then do;
                nvalid = nvalid + 1;
                _x[nvalid] = x[ix];
            end;
        end;
        if (nvalid < m1) then
            /* insufficient number of valid values. cannot proceed. */
            return (.I);

        /*------- validate & set b value --------*/
        _b = int(b);
        if (_b < 10) then
            _b = 50; /* default number of bootstrap samples */

        /* following m1low setting ensures that m2 >= sqrt(nvalid);
           alternate is to use int(sqrt(nvalid))+1, but that
           causes m2 to start from close to 1. */
        m1low = MAX(2, floor(nvalid**0.75));
        m1high = nvalid - 1;
        if (s < 2) then
            m1steps = 10;
        else
            m1steps = s;
        m1incr = int((m1high-m1low)/m1steps)-1;

        minQratio = .;
        kopt = .;
        array samp[1] /nosymbols;
        array mse[1] /nosymbols;

        /*---- iterate over m1 (sample size for the first step) ----*/
        do m1 = m1low to m1high by m1incr;
            /*------- Step 1 (samples of size m1) --------*/
            call dynamic_array(samp, m1);
            call dynamic_array(mse, m1);

            k1 = svrtutil_hillutil_minMSEindex(nvalid, _x, _b, m1, samp,
                                               mse, Q1);
            if (k1 > 1) then do;
                /*------- Step 2 (samples of size m2) --------*/
                m2 = int(m1*m1/nvalid);
                k2 = svrtutil_hillutil_minMSEindex(nvalid, _x, _b, m2, samp,
                                                   mse, Q2);

                if (Q2 = 0) then
                    Qratio = .;
                else
                    Qratio = Q1**2/Q2;

                /*------- Compute cutoff index --------*/
                _lk1_by_lm1 = log(k1)/log(m1);
                k = (k1**2/k2) * (2/_lk1_by_lm1 - 1)**(-2 + 2*_lk1_by_lm1);

                if (not(missing(Qratio))) then do;
                    /*----- Check if Qratio is minimal for current m1 -----*/
                    if (missing(minQratio) or Qratio < minQratio) then do;
                        kopt = k;
                        minQratio = Qratio;
                    end;
                end;
            end;
        end; /*----- end loop over m1 -----*/

        /*------- return element at rank kopt+1 as cutoff point -------*/
        call svrtutil_sort(nvalid, _x, 0);
        if (not(missing(kopt)) and kopt < nvalid) then do;
            status = 0;
            return(_x[int(kopt+1)]);
        end;
        else do;
            /* the algorithm failed to compute valid cutoff index;
               return fifth largest value in the sample. */
put 'NOTE: SVRTUTIL_HILLCUTOFF failed to compute a valid cutoff point.' @@;
put 'Returning the 5th largest value in the sample as the cutoff point.';

            status = 1;
            return (_x[nvalid-4]);

            /* alternate is to return 95th percentile, as follows:
            array _xuniq[1] /nosymbols;
            array _F[1] /nosymbols;
            call dynamic_array(_xuniq, nvalid);
            call dynamic_array(_F, nvalid);
            j = 0;
            eps = 2.220446E-16; * constant('MACEPS') *;
            do i=1 to nvalid;
                k = i+1;
                do while(k <= nvalid);
                    if (abs(_x[i]-_x[k]) <= eps) then
                        k = k+1;
                    else
                        go to _eqloopend;
                end;
                _eqloopend:;
                j = j+1;
                _xuniq[j] = _x[i];
                _F[j] = (k-1)/nvalid;
            end;

            return(SVRTUTIL_PERCENTILE(0.95, j, _xuniq, _F));
            */
        end;
    endsub;

     /*======================================================================
      * Computes and returns the empirical distribution function estimate for
      * the specified value of the response variable by using the EDF for
      * the input data sample. percentile of the given sample
      * given an estimate of its empirical distribution function (EDF).
      * This is an empirical estimate of the percentile.
      *
      * Arguments:
      * y        - Response variable value for which you need the EDF.
      * dim      - Dimension of the input array 'x'.
      * x        - Input array that contains the values in the sample.
      *            Must be sorted in ascending order. It is assumed that x
      *            contains distinct values.
      * F        - Input array that contains an estimate of the EDF for the
      *            sample in 'x'. F[i] is expected to contain an empirical
      *            estimate of Pr(X <= x[i]).
      * Ftype    - Type of estimator that was used to compute F
      *            If 1, standard estimator
      *            If 2, a product-limit estimator (Kaplan-Meier)
      *            If 3, Turnbull's estimator
      *            Types 1 & 2 implies individual data, whereas Type 3 implies
      *            grouped (interval-censored) data.
      *            This function treats types 1 and 2 identically.
      *            For types 1 & 2, x and F are provided at individual data
      *            points such that EDF is constant and equal to F[i] for the
      *            semi-closed interval [x[i],x[i+1]).
      *            For grouped data, x and F are for the disjoint intervals
      *            such that F is constant between intervals and the increase
      *            in F from the left endpoint of the interval to the right is
      *            arbitrary. This function assumes that F is a linear function
      *            within an interval.
      *            Note that for some intervals, left- and right-
      *            endpoints might be identical to indicate a point-interval,
      *            which allows for a mix of individual and grouped data.
      *======================================================================*/
    function SVRTUTIL_EDF(y, dim, x[*], F[*], Ftype);
        if (missing(y)) then
            return(.);

        eps = 2.220446E-16; /* constant('MACEPS') */
        if (Ftype = 1 or Ftype = 2) then do;
            /* individual data */
            if (y <= x[1]+eps) then do; /* y <= x[1] */
                if (y < x[1]-eps) then /* y < x[1] */
                    return (0); /* assume F is 0 in [0,x[1]) */
                else /* y = x[1] */
                    return (F[1]);
            end;

            if (y > x[dim]+eps) then /* y > x[dim] */
                return (F[dim]); /* assume that F is flat after last point */

            /* at this point, dim > 1 and y is in (x[1], x[dim]] interval */
            edf = .;
            ileft = 2;
            iright = dim;
            i = int((dim+2)/2);
            do while (i >= 2 and i <= dim);
                if (y < x[i]-eps) then do;
                    if (y >= x[i-1]-eps) then do; /* x[i-1] <= y < x[i] */
                        edf = F[i-1]; /* F is F[i-1] in [x[i-1],x[i]) */
                        go to doneEDF;
                    end;
                    /* y < x[i-1]; move search to left half of
                       current interval */
                    iright = i-1;
                end;
                else do;
                    if (y <= x[i]+eps) then do; /* y = x[i] */
                        edf = F[i]; /* F jumps to F[i] at x[i] */
                        go to doneEDF;
                    end;
                    /* y > x[i]; move search to right half of
                       current interval */
                    ileft = i+1;
                end;
                i = int((ileft+iright)/2);
            end;
        end;
        else do;
            /* grouped data */
            if (y <= x[1]+eps) then do;
                /* y <= x[1] */
                if (y < x[1]-eps) then /* y < x[1] */
                    return (0); /* assume F is 0 in [0,x[1]) */
                else do;
                    /* y = x[1] */
                    if (1 < dim) then do;
                        if (x[2] <= x[1]+eps) then
                            edf = F[2];
                        else
                            edf = F[1];
                    end;
                    else
                        edf = F[1];
                    go to doneEDF;
                end;
            end;

            /* at this point, dim > 1 and y is in (x[1], x[dim]] interval */
            if (y > x[dim]+eps) then /* y > x[dim] */
                return (F[dim]); /* assume that F is flat after last point */

            edf = .;
            ileft = 2;
            iright = dim;
            i = int((dim+2)/2);

            do while (i >= 2 and i <= dim);
                if (y < x[i]-eps) then do;
                    if (y > x[i-1]+eps) then do; /* x[i-1] < y < x[i] */
                        edf = F[i-1] + (y-x[i-1])*(F[i]-F[i-1])/(x[i]-x[i-1]);
                        go to doneEDF;
                    end;
                    if (y >= x[i-1]-eps) then do; /* y = x[i-1] */
                        if (x[i] <= x[i-1]+eps) then
                            edf = F[i];
                        else
                            edf = F[i-1];
                        go to doneEDF;
                    end;
                    /* y < x[i-1]; move search to left half of
                       current interval */
                    iright = i-1;
                end;
                else do;
                    if (y <= x[i]+eps) then do; /* y = x[i] */
                        if (i < dim) then do;
                            if (x[i+1] <= x[i]+eps) then
                                edf = F[i+1];
                            else
                                edf = F[i];
                        end;
                        else
                            edf = F[i];
                        go to doneEDF;
                    end;
                    /* y > x[i]; move search to right half of
                       current interval */
                    ileft = i+1;
                end;
                i = int((ileft+iright)/2);
            end;
        end;
        doneEDF:
            return(edf);
    endsub;

    function SVRTUTIL_EMPLIMMOMENT(k, u, dim, x[*], F[*], Ftype);
        if (missing(k) or missing(u)) then
            return (.);
        if (k < 1) then
            return (.);

        eps = 2.220446E-16; /* constant('MACEPS') */

        if (u < x[1]-eps) then do;
            elm = u**k;
            go to doneELMOM;
        end;

        kp1 = k+1;
        kokp1 = k/(k+1);

        elm = x[1]**k;
        if (Ftype = 1 or Ftype = 2) then do;
            /* individual data */
            do i = 1 to dim-1;
                if (u <= x[i]+eps) then do;
                    go to doneELMOM;
                end;

                if (u < x[i+1]-eps) then do;
                    elm = elm + (1-F[i])*(u**k - x[i]**k);
                    go to doneELMOM;
                end;

                /* u >= x[i+1] */
                elm = elm + (1-F[i])*(x[i+1]**k - x[i]**k);
            end;
        end;
        else do;
            /* grouped data */
            do i = 1 to dim-1;
                if (u <= x[i]+eps) then do;
                    go to doneELMOM;
                end;

                if (x[i+1] > x[i]+eps) then do;
                    si = (F[i+1]-F[i])/(x[i+1]-x[i]);
                    if (si > eps) then do;
                        if (u < x[i+1]-eps) then do;
                            elm = elm + (1-F[i]+si*x[i])*(u**k - x[i]**k)
                                      - si*kokp1*(u**kp1 - x[i]**kp1);
                            go to doneELMOM;
                        end;

                        /* u >= x[i+1] */
                        elm = elm + (1-F[i]+si*x[i])*(x[i+1]**k - x[i]**k)
                                  - si*kokp1*(x[i+1]**kp1 - x[i]**kp1);
                    end;
                    else do;
                        /* si=0 */
                        if (u < x[i+1]-eps) then do;
                            elm = elm + (1-F[i])*(u**k - x[i]**k);
                            go to doneELMOM;
                        end;

                        /* u >= x[i+1] */
                        elm = elm + (1-F[i])*(x[i+1]**k - x[i]**k);
                    end;
                end;
            end;
        end;

        if (u > x[dim]+eps) then do;
            elm = elm + (1-F[dim])*(u**k - x[dim]**k);
        end;

        doneELMOM:
            return(elm);
    endsub;
quit;