Estimating the Probability Distribution of Insurance Payments

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

                    SAS Sample Library

        Name: hcdmex01.sas
 Description: Example program from SAS/ETS High Performance
              Procedures User's Guide, The HPCDM Procedure
       Title: Estimating the Probability Distribution of Insurance Payments
     Product: SAS High Performance Econometrics Software
        Keys: Compound Distribution Modeling,
              Aggregate Loss Modeling
        PROC: HPCDM
       Notes:

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


/* Simulate data for losses incurred by several policyholders
   in some period of time */
data losses(keep=policyholderId age gender carType annualMiles
                 education carSafety income noloss lossamount);
   call streaminit(12345);
   array cx{5} age gender carType annualMiles education;
   array cbeta{6} _TEMPORARY_ (1 -0.75 1 0.6 -1 -0.25);

   array sx{3} carType carSafety income;
   array sbeta{4} _TEMPORARY_ (3.5 1.5 -0.8 0.6);

   alpha = 1/3; theta = 1/alpha;
   Sigma = 1;
   do policyholderId=1 to 5000;
      /* simulate policyholder and vehicle attributes */
      age = MAX(int(rand('NORMAL', 35, 15)),16)/50;

      if (rand('UNIFORM') < 0.5) then gender = 1; * female;
      else gender = 2; * male;

      if (rand('UNIFORM') < 0.7) then carType = 1; * sedan;
      else carType = 2; * SUV;

      annualMiles = MAX(1000, int(rand('NORMAL', 12000, 5000)))/5000;

      educationLevel = rand('UNIFORM');
      if (educationLevel < 0.5) then education = 1; *high school graduate;
      else if (educationLevel < 0.85) then education = 2; *college graduate;
      else education = 3; *advanced degree;

      carSafety = rand('UNIFORM'); /* scaled to be between 0 & 1 */

      income = MAX(15000,int(rand('NORMAL', education*30000, 50000)))/100000;

      /* simulate number of losses incurred by this policyholder */
      cxbeta = cbeta(1);
      do i=1 to dim(cx);
         cxbeta = cxbeta + cx(i) * cbeta(i+1);
      end;
      Mu = exp(cxbeta);
      p = theta/(Mu+theta);
      numloss = rand('NEGB',p,theta);

      /* simulate severity of each of the losses */
      if (numloss > 0) then do;
         noloss = 0;
         do iloss=1 to numloss;
            Mu = sbeta(1);
            do i=1 to dim(sx);
               Mu = Mu + sx(i) * sbeta(i+1);
            end;
            lossamount = exp(Mu) * rand('LOGNORMAL')**Sigma;
            output;
         end;
      end;
      else do;
         noloss = 1;
         lossamount = .;
         output;
      end;
   end;
run;

/* Aggregate number of losses for each policyholder */
data losscounts(keep=age gender carType annualMiles education numloss);
   set losses;
   by policyholderId;
   retain numloss 0;
   if (noloss ne 1) then
      numloss = numloss + 1;
   if (last.policyholderId) then do;
      output;
      numloss = 0;
   end;
run;

/* Fit count regression model for the number of losses */
proc countreg data=losscounts;
   model numloss = age gender carType annualMiles education / dist=negbin;
   store work.countregmodel;
run;

/* Fit severity scale regression model for the loss severity */
proc severity data=losses plots=none outest=work.sevregest;
   loss lossamount;
   scalemodel carType carSafety income;
   dist logn;
run;

/* Generate the scenario data set for multiple policyholders */
data groupOfPolicies(keep=policyholderId age gender carType annualMiles
                          education carSafety income
                          lowDeductible highDeductible limit annualLimit);
   call streaminit(67897);

   do policyholderId=1 to 5;
      age = MAX(int(rand('NORMAL', 35, 15)),16)/50;

      if (rand('UNIFORM') < 0.5) then gender = 1; * female;
      else gender = 2; * male;

      if (rand('UNIFORM') < 0.7) then carType = 1; * sedan;
      else carType = 2; * SUV;

      annualMiles = MAX(1000, int(rand('NORMAL', 12000, 5000)))/5000;

      educationLevel = rand('UNIFORM');
      if (educationLevel < 0.5) then education = 1; *high school graduate;
      else if (educationLevel < 0.85) then education = 2; *college graduate;
      else education = 3; *advanced degree;

      carSafety = rand('UNIFORM'); /* scaled to be between 0 & 1 */

      income = MAX(15000,int(rand('NORMAL', education*30000, 50000)))/100000;

      lowDeductible = 100*(1+floor(rand('UNIFORM')*5));
      highDeductible = lowDeductible + 500*(1+floor(rand('UNIFORM')*2));
      limit = 2500*(1+floor(rand('UNIFORM')*3));
      annualLimit = 10000*(1+floor(rand('UNIFORM')*2));

      output;
   end;
run;

/* Simulate the aggregate loss distribution and aggregate adjusted
   loss distribution for the scenario with multiple policyholders */
proc hpcdm data=groupOfPolicies nreplicates=10000 seed=13579 print=all
           countstore=work.countregmodel severityest=work.sevregest
           plots=(edf pdf) nperturbedSamples=50
           adjustedseverity=amountPaid;
   severitymodel logn;

   if (_sev_ <= lowDeductible) then
      amountPaid = 0;
   else do;
      if (_sev_ <= highDeductible) then
         amountPaid = highDeductible *
            (_sev_-lowDeductible)/(highDeductible-lowDeductible);
      else
         amountPaid = MIN(_sev_, limit); /* imposes per-loss payment limit */
   end;
run;


/* Simulate the aggregate loss distribution and aggregate adjusted
   loss distribution for the modified set of policy provisions */
proc hpcdm data=groupOfPolicies nreplicates=10000 seed=13579 print=all
           countstore=work.countregmodel severityest=work.sevregest
           plots=none nperturbedSamples=50
           adjustedseverity=amountPaid;
   severitymodel logn;

   if (_sev_ <= lowDeductible) then
      amountPaid = 0;
   else do;
      if (_sev_ <= highDeductible) then
         amountPaid = highDeductible *
            (_sev_-lowDeductible)/(highDeductible-lowDeductible);
      else
         amountPaid = MIN(_sev_, limit); /* imposes per-loss payment limit */

      /* impose policyholder's annual limit */
      amountPaid = MIN(amountPaid, MAX(0,annualLimit - _cumadjsevforobs_));

      /* impose group's annual limit */
      amountPaid = MIN(amountPaid, MAX(0,15000 - _cumadjsev_));
   end;
run;