program bayes;

(*
   PROJECT          : Bayesian estimator for full GP1 model
   FILE-SPEC.       : bayesest.sp
   AUTHOR           : R.-D. Reiss, M. Thomas
   CREATION DATE    : Oct-9-1999
   LAST UPDATE      : Oct-9-1999
   RELATED DOCUMENTS: "A new class for Bayesian estimators in
                      Paretian excess-of-loss reinsurance" by
                      R.-D. Reiss and M. Thomas
   VERSION          : 1.0
   FUNCTIONAL DESCR.: Estimation of alpha and sigma as given in
                      formula (7) of the above document.
                      The estimator is implemented completely in 
                      the procedure "BayesFullGP1Est" and can 
                      be extracted easily
   SPECIAL FEATURES : requires Xtremes Version 3.0 Beta 4 or higher,
                      available from http://www.xtremes.de/
   INPUT            : Xtremes Univariate Data as active data set;
                      interactive input of threshold u and parameters
                      a, b, s, d of reciprocal gamma priors
   OUTPUT           : Estimated values alpha and sigma
*)

type realfunc = function (real): real;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% BayesFullGP1Est performs the estimation. The procedure is
% self-contained and can be extracted from the program easily.
%
% Parameters are:
%
%  x:            vector with normalized excesses
%  s, d:         parameters of reciprocal gamma prior for alpha
%  a, b:         parameters of reciprocal gamma prior for sigma
%  alphaest, 
%  sigmaest      return estimated values
%
% To use another prior distribution, change the definition of p.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

procedure BayesFullGP1Est (x: realvector; s, d, a, b: real;
                           var alphaest, sigmaest: real);

    const sigma0 = 0.00001;	% ranges for numerical integration
          sigma1 = 30.0;

    var k: integer;
        help: real;
    
    function p (x: real): real;
        begin
            return gammadensity (a, b / x) * b / sqr (x)
        end;

    % Adaptive Gaussian integration taken from Klugman, S.A. (1992):	
    % Bayesian Statistics in Actuarial Science. Kluwer, Boston.        

    function GaussIntegrate (f: realfunc; a, b: real): real;

        var x, w: realvector;
            p: real;

        function integrate (a, b: real): real;
            begin
                return sum ((b - a) * 0.5 * w * f ((b + a) * 0.5 + (b - a) * 0.5 * x))
            end;

        function improve (a, b, prev, err: real; level: integer): real;
            var d, v1, v2: real;
            begin
                v1 := integrate (a, 0.5 * (a + b));
                v2 := integrate (0.5 * (a + b), b);
                d := prev - v1 - v2;
                if (abs (d) < 1e6 * err) or (level = 7) then 
                    return v1 + v2
                else 
                    return improve (a, 0.5 * (a + b), v1, 0.5 * err, level + 1) +
                           improve (0.5 * (a + b), b, v2, 0.5 * err, level + 1)
            end;

        begin
            x := combine (-0.973906528517172, -0.865063366688985,
                          -0.679409568299024, -0.433395394129247,
                          -0.148874338981631,  0.148874338981631,
                           0.433395394129247,  0.679409568299024,
                           0.865063366688985,  0.973906528517172);
            w := combine ( 0.066671344308688,  0.149451349050581,
                           0.219086362515982,  0.269266719309996,
                           0.295524334714753,  0.295524334714753,        
                           0.269266719309996,  0.219086362515982,
                           0.149451349050581,  0.066671344308688);
            p := integrate (a, b);
            return improve (a, b, p, p / 1e10, 1)
        end;

    function dprime (sigma: real): real;
        begin
            return d + sum (log (1 + x / sigma))
        end;

    function ptilde (sigma: real): real;
        var help: real;
        begin
            help := -k * log (sigma) - (s+k) * log (dprime (sigma)) -sum (log (1.0 + x / sigma));
            return exp (help) * p (sigma)
        end;

    function g (sigma: real): real;
        begin
            return sigma * ptilde (sigma)
        end;

    function h (sigma: real): real;
	begin
	    return (s + k) * ptilde (sigma) / dprime (sigma)
	end;

    begin
        k := size (x);
        help := GaussIntegrate (ptilde, sigma0, sigma1);
	alphaest := GaussIntegrate (h, sigma0, sigma1) / help;
        sigmaest := GaussIntegrate (g, sigma0, sigma1) / help;
    end;
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Interactive part: Reads active data set of Xtremes and asks
% for threshold and parameters of priors.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

procedure InteractiveEstimation;

    var alphaest, sigmaest: real;
        x: realvector;
        u, s, d, a, b: real;

    procedure AskParameter;
        var v: vector of real;
        begin
            v := DialogBox ('Parameters', 's|d|u|a|b', combine (1.0, 1.0, 1.0, 1.0, 1.0));
            s := v [1]; d := v [2];
            u := v [3];
            a := v [4]; b := v [5];
        end;

    begin
        AskParameter;
        x := (ColumnData (1) - u) / u;
        BayesFullGP1Est (x, s, d, a, b, alphaest, sigmaest);
	writeln (s, ' ', d, ' ', a, ' ', b);
        MessageBox ('alpha* = ' + str (alphaest) + ', sigma* = ' + str (sigmaest));
    end;

begin
    InteractiveEstimation
end.