program bayes;

(*
   PROJECT          : Bayesian estimator for full GP1 model
   FILE-SPEC.       : bayessim.sp
   AUTHOR           : R.-D. Reiss, M. Thomas
   CREATION DATE    : Oct-9-1999
   LAST UPDATE      : Dec-22-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.: Simulation of alpha* as given in 
                      formula (7) of the above document.
		      See section 4 of the paper for a description
		      of the parameters.
   SPECIAL FEATURES : requires Xtremes Version 3.0 Beta 4 or higher,
                      available from http://www.xtremes.de/
   INPUT            : interactive input of Monte Carlo sample size (k),
		      parameters of Pareto distribution (alpha, sigma),
		      hyperparameters (s, d, a and b),
                      number of Monte Carlo simulations (N).
   OUTPUT           : Creates active data set with simulated values
		      of the shape parameter.
*)

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: Asks for parameters of simulation and creates
% data set with estimated shape parameters called bayessim.dat.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

procedure InteractiveSimulation;

    var alpha, sigma: real;
        s, d, a, b: real;
        N, k: integer;

    procedure AskParameter;
        var v: vector of real;
        begin
            v := DialogBox ('Parameters', 'k|alpha|sigma|s|d|a|b|N', combine (20.0, 4.0, 2.0, 4.0, 1.0, 3.0, 4.0, 4000));
	    k := round (v [1]);
	    alpha := v [2]; sigma := v [3];
            s := v [4]; d := v [5];
            a := v [6]; b := v [7];
	    N := round (v [8]);
        end;

    procedure DoSimulation;
 	var i: integer;
	    x, y: vector of real;
	    alphaest, sigmaest: real;
	begin
	    for i := 1 to N do begin
		x := (exp (log (1.0 - Random (k)) * (-1.0 / alpha)) - 1) * sigma;
		BayesFullGP1Est (x, s, d, a, b, alphaest, sigmaest);
		y := combine (y, alphaest)
	    end;
	    createunivariate (y, 'bayessim.dat', '')
	end;
	
    begin
        AskParameter;
	DoSimulation
    end;

begin
    InteractiveSimulation
end.