Model Definition for Generalized Pareto Distribution

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

                    SAS Sample Library

        Name: svrtgpd.sas
 Description: Example Program from SAS/ETS User's Guide,
              The SEVERITY Procedure
       Title: Model Definition for Generalized Pareto Distribution
     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.
              Further, if you run this sample without any modification to the
              names of functions and subroutines, then PROC FCMP writes
              warnings of the following nature to the SAS log:
              WARNING: Function <name> was defined in a previous package.
                       Function <name> as defined in the current program will
                       be used as default when the package is not specified.
              You can ignore such warnings; they appear because the functions
              defined in this sample are already defined in the input library
              Sashelp.Svrtdist.

--------------------------------------------------------------*/
proc fcmp library=sashelp.svrtdist outlib=sasuser.svrtdist.models;
    function GPD_DESCRIPTION() $32;
        length model $32;
        model = "Generalized Pareto Distribution";
        return(model);
    endsub;

    function GPD_PARMCOUNT();
        return(2);
    endsub;

    function GPD_LOGPDF(x, Theta, Xi);
        /* Theta : Scale */
        /* Xi :    Shape */
        if (x < 0) then do;
            return (-180.218266); /* constant('LOGSMALL') */
        end;
        else do;
            if (Xi > 2.220446E-16) then /* constant('MACEPS') */
                logpdf = (-1-(1/Xi)) * log1px(Xi*(x/Theta)) - log(Theta);
            else
                logpdf = -x/Theta - log(Theta); /* for Xi ~ 0, exponential */
            return (logpdf);
        end;
    endsub;
    function GPD_PDF(x, Theta, Xi);
        /* Theta : Scale */
        /* Xi :    Shape */
        if (x < 0) then do;
            return (0);
        end;
        else do;
            logpdf = GPD_LOGPDF(x, Theta, Xi);
            if (missing(logpdf)) then
                return (.);
            if (logpdf < 174.673089) then /* constant('LOGBIG') */
                return ( exp(logpdf) );
            else
                return ( 7.237005E75 ); /* constant('BIG') */
        end;
    endsub;

    function GPD_LOGCDF(x, Theta, Xi);
        /* Theta : Scale */
        /* Xi :    Shape */
        if (x < 0) then do;
            return (-180.218266); /* constant('LOGSMALL') */
        end;
        else do;
            if (Xi > 2.220446E-16) then /* constant('MACEPS') */
                c = (1 + Xi*(x/Theta))**(-1/Xi);
            else
                c = exp(-x/Theta); /* for Xi ~ 0, exponential */
            if (c < 1) then
                return (log1px(-c));
            else
                return (-180.218266); /* constant('LOGSMALL') */
        end;
    endsub;
    function GPD_CDF(x, Theta, Xi);
        /* Theta : Scale */
        /* Xi :    Shape */
        if (x < 0) then do;
            return (0);
        end;
        else do;
            if (Xi > 2.220446E-16) then /* constant('MACEPS') */
                return (1 - (1 + Xi*(x/Theta))**(-1/Xi));
            else
                return (1 - exp(-x/Theta)); /* for Xi ~ 0, exponential */
        end;
    endsub;

    function GPD_LOGSDF(x, Theta, Xi);
        /* Theta : Scale */
        /* Xi :    Shape */
        if (x < 0) then do;
            return (0);
        end;
        else do;
            if (Xi > 2.220446E-16) then /* constant('MACEPS') */
                return ((-1/Xi)*log1px(Xi*(x/Theta)));
            else
                return (-x/Theta); /* for Xi ~ 0, exponential */
        end;
    endsub;
    function GPD_SDF(x, Theta, Xi);
        /* Theta : Scale */
        /* Xi :    Shape */
        if (x < 0) then do;
            return (1);
        end;
        else do;
            if (Xi > 2.220446E-16) then /* constant('MACEPS') */
                return ((1 + Xi*(x/Theta))**(-1/Xi));
            else
                return (exp(-x/Theta)); /* for Xi ~ 0, exponential */
        end;
    endsub;

    subroutine GPD_PARMINIT(dim, x[*], nx[*], F[*], Ftype, Theta, Xi);
        outargs Theta, Xi;
        array m[2] / nosymbols;

        /* Use Method of Moments */
        call svrtutil_rawmoments(dim, x, nx, 2, m);

        if (missing(m[1])) then do;
            Theta = .;
            Xi    = .;
        end;
        else do;
            t1 = m[2] - m[1]**2;
            t2 = 2*t1 - m[2];

            eps = 2.220446E-16; /* constant('MACEPS') */
            if (t1 < eps or t2 < eps) then do;
                Theta = m[1]/2;
                Xi = 1/2;
            end;
            else do;
                Theta = (m[1]*m[2])/(2*t1);
                Xi    = t2/(2*t1);
            end;
        end;
    endsub;

    function GPD_QUANTILE(p, Theta, Xi);
        /* Theta : Scale */
        /* Xi :    Shape */
        return (((1-p)**(-Xi)-1)*Theta/Xi);
    endsub;

    function GPD_MEAN(x, Theta, Xi);
        if not(Xi < 1) then
            return (.); /* first moment does not exist */
        return (Theta/(1 - Xi));
    endsub;
quit;