%macro bearsim(fallharv=0,sprharv=0,doubharv=0.1,probbadb=0.1,
       spropfem=33, fpropfem=45);
* fallharv       /* Fall Harvest Rate % */;
* sprharv        /* Spring Harvest Rate % */;
* doubharv=0.1   /* Probability of double harvest */;
* probbadb=0.1   /* Probability of bad berry year */;
* spropfem=33    /* Selection of females in spring harvest -- 33% */;
* fpropfem=45    /* Selection of females in fall harvest -- 45% */;
%let maxage=20;      /* Maximum age */
%let yearsim=30;     /* Years in simulation */
%let startyr=1980;   /* Starting year of simulation */
%let lastyear=%eval(&startyr+&yearsim);
%let num_sim=100;    /* Reps to simulate */
title 'Colorado Black Bear Stochastic Simulation Model, Fig. 1 Data';
title2 "Summary Statistics for &yearsim Year simulation";
title3
"Maxage = &maxage, Years in sim = &yearsim, Start = &startyr, End = &lastyear";
title4
"Bad Berry Years with Pr = &probbadb , Double Harvest with Pr = &doubharv";
title5
"Spring Harvest Rate &sprharv %,  Fall Harvest Rate &fallharv %";
title6 "Harvest Ratio % Females for Spring &spropfem  Fall &fpropfem ";
data bearsim;

   /* House keeping details */
   seed=0;
   drop S_admal S_adfem S_4ymal S_4yfem S_3ymal S_3yfem
        S_2ymal S_2yfem S_1ywom S_1ywmom S_0ywom S_0ywmom
        rate2yr rate3yr rate4yr rate5yr rate6yr ratemoms
        mom1cubs mom2cubs mom3cubs sexratio;
   drop i j seed deaths age lived nomoved breed breed2;
   drop sprhrvm sprhrvf fallhrvm fallhrvf;
   array subadlt{0:1} cubs_wom year_wom;
   array fem{2:&maxage,0:6}
         fem2_0-fem2_6   fem3_0-fem3_6   fem4_0-fem4_6
         fem5_0-fem5_6   fem6_0-fem6_6   fem7_0-fem7_6
         fem8_0-fem8_6   fem9_0-fem9_6   fem10_0-fem10_6
         fem11_0-fem11_6 fem12_0-fem12_6 fem13_0-fem13_6
         fem14_0-fem14_6 fem15_0-fem15_6 fem16_0-fem16_6
         fem17_0-fem17_6 fem18_0-fem18_6 fem19_0-fem19_6
         fem20_0-fem20_6 ;
   array males{2:&maxage} males2-males&maxage;

   keep year sprharv fallharv spropfem fpropfem
        cubs yearling adult_m adult_f pop r
        msprhrv fsprhrv mfallhrv ffallhrv
        cubs_wom year_wom;

   /* Probability of a bad berry year */   probbadb=&probbadb;
   /* Probability of double harvest */     doubharv=&doubharv;

   /* History of berry years */
   array bryhist{0:5};
   do i=0 to 5;
      bryhist{i}=0;
      if ranuni(seed)<probbadb then bryhist{i}=1;
   end;

   /* Set survival rates */
   array S_female{2:5} S_2yfem S_3yfem S_4yfem S_adfem;
   S_adfem=0.96;           /* Adult Females 5+         */
   S_4yfem=0.96;           /* 4 year olds Females      */
   S_3yfem=0.94;           /* 3 year olds Females      */
   S_2yfem=0.94;           /* 2 year olds Females      */
   array S_males {2:5} S_2ymal S_3ymal S_4ymal S_admal;
   S_admal=0.70;           /* Adult Males 5+           */
   S_4ymal=0.70;           /* 4 year olds Males        */
   S_3ymal=0.76;           /* 3 year olds Males        */
   S_2ymal=0.76;           /* 2 year olds Males        */
   S_1ywom=0.94;           /* yearlings without mother */
   S_1ywmom=0.94;          /* yearlings with mother    */
              /* Good Berry Years */
   S_0ywomG=0.56;          /* Cubs without mother      */
   S_0ywmG=0.56;           /* Cubs with mother         */
              /*Bad Berry Years */
   S_0ywomB=0.33;          /* Cubs without mother      */
   S_0ywmB=0.33;           /* Cubs with mother         */


   /* Set reproduction parameters */
   array reprate{2:7} rate2yr rate3yr rate4yr rate5yr rate6yr ratemoms;
         /* Good Berry Years */
   Grate2yr=0.001;    /* Proportion of barren 2 year olds having cubs */
   Grate3yr=0.22;     /* Proportion of barren 3 year olds having cubs */
                      /* 2 of 9, page 30 */
   Grate4yr=0.29;     /* Proportion of barren 4 year olds having cubs */
                      /* 2 of 7, page 30 */
   Grate5yr=0.60;     /* Proportion of barren 5 year olds having cubs */
                      /* 3 of 5, page 30 */
   Grate6yr=0.50;     /* Proportion of barren 6+ year olds having cubs */
                      /* 1 of 2, page 30 */
   Gratemom=0.80;     /* Proportion of adults with yearlings having cubs */
                      /* 0.8 gives mean litter interval of 2.2 years */
         /* Bad Berry Years */
   Brate2yr=0.001;    /* Proportion of barren 2 year olds having cubs */
   Brate3yr=0.001;    /* Proportion of barren 3 year olds having cubs */
   Brate4yr=0.001;    /* Proportion of barren 4 year olds having cubs */
   Brate5yr=0.001;    /* Proportion of barren 5 year olds having cubs */
   Brate6yr=0.20;     /* Proportion of barren 6+ year olds having cubs */
   Bratemom=0.20;     /* Proportion of adults with yearlings having cubs */
         /* Year following a bad year */
   Frate6yr=0.90;     /* Proportion of barren 6+ year olds having cubs */
   Fratemom=0.90;     /* Proportion of adults with yearlings having cubs */

   /* Set proportions of 1, 2, and 3 cubs per mother */
   mom1cubs=0.20;    /* Proportion of mothers with 1 cub */
   mom2cubs=0.60;    /* Proportion of mothers with 2 cubs */
   mom3cubs=1-mom1cubs-mom2cubs;  /* Rest have 3 cubs */

   /* Set sex ratio -- Proportion of males in yearling population */
   sexratio=0.50;

   /* Harvest rates -- % females */
   spropfem=&spropfem/100;
   fpropfem=&fpropfem/100;
   put 'Selection Constants for Females in Harvest ' spropfem= fpropfem=;

   do isim=1 to &num_sim;

      /* Set initial conditions -- Cubs have just been born */
CUBS_WOM=0; YEAR_WOM=82;
FEM2_0=532; FEM2_1=0; FEM2_2=0; FEM2_3=0; FEM2_4=0; FEM2_5=0; FEM2_6=0;
FEM3_0=517; FEM3_1=0; FEM3_2=0; FEM3_3=0; FEM3_4=0; FEM3_5=0; FEM3_6=0;
FEM4_0=334; FEM4_1=15; FEM4_2=45; FEM4_3=15; FEM4_4=0; FEM4_5=0; FEM4_6=0;
FEM5_0=261; FEM5_1=16; FEM5_2=48; FEM5_3=16; FEM5_4=32; FEM5_5=18; FEM5_6=2;
FEM6_0=126; FEM6_1=32; FEM6_2=96; FEM6_3=32; FEM6_4=35; FEM6_5=19; FEM6_6=2;
FEM7_0=86; FEM7_1=23; FEM7_2=71; FEM7_3=23; FEM7_4=62; FEM7_5=33; FEM7_6=4;
FEM8_0=73; FEM8_1=24; FEM8_2=71; FEM8_3=24; FEM8_4=50; FEM8_5=27; FEM8_6=3;
FEM9_0=67; FEM9_1=20; FEM9_2=62; FEM9_3=20; FEM9_4=54; FEM9_5=29; FEM9_6=4;
FEM10_0=53; FEM10_1=18; FEM10_2=54; FEM10_3=18; FEM10_4=42; FEM10_5=23; FEM10_6=3;
FEM11_0=50; FEM11_1=16; FEM11_2=49; FEM11_3=16; FEM11_4=43; FEM11_5=24; FEM11_6=3;
FEM12_0=42; FEM12_1=13; FEM12_2=39; FEM12_3=13; FEM12_4=35; FEM12_5=19; FEM12_6=2;
FEM13_0=39; FEM13_1=13; FEM13_2=40; FEM13_3=13; FEM13_4=32; FEM13_5=18; FEM13_6=2;
FEM14_0=36; FEM14_1=11; FEM14_2=34; FEM14_3=11; FEM14_4=29; FEM14_5=16; FEM14_6=2;
FEM15_0=31; FEM15_1=10; FEM15_2=32; FEM15_3=11; FEM15_4=26; FEM15_5=14; FEM15_6=2;
FEM16_0=28; FEM16_1=9; FEM16_2=26; FEM16_3=9; FEM16_4=22; FEM16_5=12; FEM16_6=1;
FEM17_0=25; FEM17_1=8; FEM17_2=25; FEM17_3=8; FEM17_4=21; FEM17_5=11; FEM17_6=1;
FEM18_0=23; FEM18_1=7; FEM18_2=22; FEM18_3=7; FEM18_4=18; FEM18_5=10; FEM18_6=1;
FEM19_0=19; FEM19_1=7; FEM19_2=20; FEM19_3=7; FEM19_4=16; FEM19_5=9; FEM19_6=1;
FEM20_0=18; FEM20_1=6; FEM20_2=18; FEM20_3=6; FEM20_4=15; FEM20_5=8; FEM20_6=1;
MALES2=530; MALES3=418; MALES4=267; MALES5=186; MALES6=118; MALES7=76; MALES8=50;
MALES9=34; MALES10=21; MALES11=14; MALES12=9; MALES13=6; MALES14=4; MALES15=3;
MALES16=2; MALES17=1; MALES18=1; MALES19=0; MALES20=0;

      /* Cumulate some summary statistics on initial population */
      /* and output, but only for the first simulation.         */
      if isim=1 then do;
         year=&startyr;
         link summary;
         output;
      end;

      do year=&startyr+1 to &lastyear;

         /* Set harvest parameters */
         if ranuni(seed)<doubharv then do;
            sprharv=(&sprharv/100)*2;    /* Overall Spring harvest */
            fallharv=(&fallharv/100)*2;  /* Overall Fall harvest */
         end;
         else do;
            sprharv=(&sprharv/100);      /* Overall Spring harvest */
            fallharv=(&fallharv/100);    /* Overall Fall harvest */
         end;

         /* Increment history of berry failures */
         do i=4 to 0 by -1;
            bryhist{i+1}=bryhist{i};
         end;
         if ranuni(seed)<probbadb then bryhist{0}=1;
         else bryhist{0}=0;

         /* Initialize parameters depending on Good or Bad Berry Year */
         /* Current year */
         if bryhist{0}>0.5 then do;
            S_0ywmom=S_0ywmB; S_0ywom=S_0ywomB;
            rate2yr=Brate2yr; rate3yr=Brate3yr; rate4yr=Brate4yr;
            rate5yr=Brate5yr; rate6yr=Brate6yr; ratemoms=Bratemom;
         end;
         else do;
            S_0ywmom=S_0ywmG; S_0ywom=S_0ywomG;
            rate2yr=Grate2yr; rate3yr=Grate3yr; rate4yr=Grate4yr;
            rate5yr=Grate5yr; rate6yr=Grate6yr; ratemoms=Gratemom;
            if bryhist{1}>0.5 then do;  /* Was a bad year last year */
               rate6yr=Frate6yr; ratemoms=Fratemom;
            end;
                 /* Now put in residual effect of previous years */
            if bryhist{2}>0.5 then do;  /* Was a bad year 2 years ago */
               rate2yr=Brate2yr; end;
            if bryhist{3}>0.5 then do;  /* Was a bad year 3 years ago */
               rate3yr=Brate3yr; end;
            if bryhist{4}>0.5 then do;  /* Was a bad year 4 years ago */
               rate4yr=Brate4yr; end;
            if bryhist{5}>0.5 then do;  /* Was a bad year 5 years ago */
               rate5yr=Brate5yr; end;
         end;

         /* Spring harvest */
         msprhrv=0;
         fsprhrv=0;
         if sprharv>0 then do;
            adult_m=0;
            adult_f=0;
            do age=2 to &maxage;
               adult_m=adult_m+males{age};
               do i=0 to 6;
                  adult_f=adult_f+fem{age,i};
               end;
            end;
            harvest=sprharv*(adult_m+adult_f);
            /* Never harvest more than 90% of the population */
            sprhrvm=min(0.9,(1-spropfem)*harvest/adult_m);
            sprhrvf=min(0.9,spropfem*harvest/adult_f);
            do age=2 to &maxage;
               /* Males */
               if males{age}>0 & sprhrvm>0 then do;
                  deaths=ranbin(seed,males{age},sprhrvm);
                  males{age}=males{age}-deaths;
                  msprhrv=msprhrv+deaths;
               end;
               /* Females */
               do i=0 to 6;
                  if fem{age,i} > 0  & sprhrvf>0 then do;
                     deaths=ranbin(seed,fem{age,i},sprhrvf);
                     fem{age,i}=fem{age,i}-deaths;
                     /* Cubs are all assumed to die */
                     /* If  yearlings are present,
                        move these to no mother category */
                     if i>=4 then subadlt{1}=subadlt{1}+deaths*(i-3);
                     fsprhrv=fsprhrv+deaths;
                  end;
               end;
            end;
         end;

         /* Fall harvest */
         mfallhrv=0;
         ffallhrv=0;
         if fallharv>0 then do;
            adult_m=0;
            adult_f=0;
            do age=2 to &maxage;
               adult_m=adult_m+males{age};
               do i=0 to 6;
                  adult_f=adult_f+fem{age,i};
               end;
            end;
            harvest=fallharv*(adult_m+adult_f);  /* Expected total harvest */
            /* Never harvest more than 90% of the population */
            fallhrvm=min(0.9,(1-fpropfem)*harvest/adult_m);
            fallhrvf=min(0.9,fpropfem*harvest/adult_f);
            do age=2 to &maxage;
               /* Adult Males */
               if males{age}>0 & fallhrvm>0 then do;
                  deaths=ranbin(seed,males{age},fallhrvm);
                  males{age}=males{age}-deaths;
                  mfallhrv=mfallhrv+deaths;
               end;
               /* Females */
               do i=0 to 6;
                  if fem{age,i} > 0 & fallhrvm>0 then do;
                     deaths=ranbin(seed,fem{age,i},fallhrvf);
                     fem{age,i}=fem{age,i}-deaths;
                     /* If cubs or yearlings are present,
                        move these to no mother category */
                     if i>=1 & i<=3 then subadlt{0}=subadlt{0}+deaths*i;
                     else if i>=4 then subadlt{1}=subadlt{1}+deaths*(i-3);
                     ffallhrv=ffallhrv+deaths;
                  end;
               end;
            end;
         end;

         /* Survival from natural processes */
         do age=2 to &maxage;
            /* Process adult survival */
            /* Males */
            if males{age}>0 then
            males{age}=ranbin(seed,males{age},S_males{min(age,5)});
            /* Females */
            do i=0 to 6;
               if fem{age,i} > 0 then do;
                  deaths=ranbin(seed,fem{age,i},1-S_female{min(age,5)});
                  fem{age,i}=fem{age,i}-deaths;
                  /* If cubs or yearlings are present,
                     move these to no mother category */
                  if i>=1 & i<=3 then subadlt{0}=subadlt{0}+deaths*i;
                  else if i>=4 then subadlt{1}=subadlt{1}+deaths*(i-3);
               end;
            end;
            /* Process cub survival for cubs with mothers */
            do i=1 to 3;
               nomoved=0;
               do j=1 to fem{age,i};
                  lived=ranbin(seed,i,S_0ywmom);
                  if lived < i then do;
                     fem{age,lived}=fem{age,lived}+1;
                     nomoved=nomoved+1;
                  end;
               end;
               fem{age,i}=fem{age,i}-nomoved;
            end;

            /* Process yearling survival for yearlings with mothers */
            do i=4 to 6;
               nomoved=0;
               do j=1 to fem{age,i};
                  lived=ranbin(seed,i-3,S_1ywmom);
                  if lived < i-3 then do;
                     if lived > 0 then fem{age,lived+3}=fem{age,lived+3}+1;
                     else fem{age,0}=fem{age,0}+1;
                     nomoved=nomoved+1;
                  end;
               end;
               fem{age,i}=fem{age,i}-nomoved;
            end;
         end; /* All adult mortality processed, including offspring */

         /* Process subadlt mortality for animals without mothers */
         if subadlt{1} > 0 then
            subadlt{1}=ranbin(seed,subadlt{1},S_1ywom); /* Yearlings w/o mom */
         if subadlt{0} > 0 then
            subadlt{0}=ranbin(seed,subadlt{0},S_0ywom); /* Cubs w/0 mom */

         /* Increment ages and reproduce */
         /* Maxage year old females die: cubs and yearlings loose mother */
         subadlt{0}=subadlt{0}+fem{&maxage,1}+fem{&maxage,2}*2+fem{&maxage,3}*3;
         subadlt{1}=subadlt{1}+fem{&maxage,4}+fem{&maxage,5}*2+fem{&maxage,6}*3;
         do age=&maxage-1 to 2 by -1;
            /* Clear next age category to 0 */
            males{age+1}=0;
            do i=0 to 6;
               fem{age+1,i}=0;
            end;

            /* Increment male ages */
            males{age+1}=males{age};

            /* Now do females */
            /* Wean yearlings */
            subadlt{1}=subadlt{1}+fem{age,4}+fem{age,5}*2+fem{age,6}*3;
            /* Age cubs up one */
            do i=1 to 3;
               fem{age+1,i+3}=fem{age,i};
            end;

            /* Reproduction for barren females */
            if fem{age,0} > 0 then do;
               breed=ranbin(seed,fem{age,0},reprate{min(age,6)});
               fem{age+1,0}=fem{age+1,0}+fem{age,0}-breed; /* Non-breeders */
               if breed>0 then do;
                  breed2=ranbin(seed,breed,mom1cubs);          /* 1 cub */
                  fem{age+1,1}=fem{age+1,1}+breed2;
                  breed=breed-breed2;
                  if breed>0 then do;
                     breed2=ranbin(seed,breed,
                        mom2cubs/(mom2cubs+mom3cubs));         /* 2 cubs */
                     fem{age+1,2}=fem{age+1,2}+breed2;
                     fem{age+1,3}=fem{age+1,3}+(breed-breed2); /* 3 cubs */
                  end;
               end;
            end;

            /* Reproduction for females previously with yearlings */
            if fem{age,4}+fem{age,5}+fem{age,6} > 0 then do;
               breed=ranbin(seed,fem{age,4}+fem{age,5}+fem{age,6},ratemoms);
                                                        /* Non-breeders */
               fem{age+1,0}=fem{age+1,0}+fem{age,4}+fem{age,5}+fem{age,6}-breed;
               if breed>0 then do;
                  breed2=ranbin(seed,breed,mom1cubs);          /* 1 cub */
                  fem{age+1,1}=fem{age+1,1}+breed2;
                  breed=breed-breed2;
                  if breed>0 then do;
                     breed2=ranbin(seed,breed,
                        mom2cubs/(mom2cubs+mom3cubs));         /* 2 cubs */
                     fem{age+1,2}=fem{age+1,2}+breed2;
                     fem{age+1,3}=fem{age+1,3}+(breed-breed2); /* 3 cubs */
                  end;
               end;
            end;
         end;

         /* Split yearlings into males and females and age one year */
         if subadlt{1}>0 then
         males{2}=ranbin(seed,subadlt{1},sexratio);
         fem{2,0}=subadlt{1}-males{2};

         /* Age cubs */
         subadlt{1}=subadlt{0};

         /* No motherless cubs on day cubs are born */
         subadlt{0}=0;

         /* Cumulate some summary statistics and output */
         if year=&lastyear then do;
            link summary;
            sprharv=(&sprharv/100);      /* Overall Spring harvest */
            fallharv=(&fallharv/100);    /* Overall Fall harvest */
            output;
         end;
      end;
   end;
   return;
   summary:  /* Define some summary parameters in a subroutine */
      cubs=0;
      yearling=subadlt{1};
      adult_m=0;
      adult_f=0;
      do age=2 to &maxage;
         cubs=cubs+fem{age,1}+fem{age,2}*2+fem{age,3}*3;
         yearling=yearling+fem{age,4}+fem{age,5}*2+fem{age,6}*3;
         adult_m=adult_m+males{age};
         do i=0 to 6;
            adult_f=adult_f+fem{age,i};
         end;
      end;
      pop=cubs+yearling+adult_m+adult_f;
      r=((pop/10000)**(1/&yearsim)-1)*100;
   return;
run;
proc sort;
   by fallharv sprharv spropfem fpropfem;
proc means data=bearsim n mean std stderr min max;
   by fallharv sprharv spropfem fpropfem;
   where year=&lastyear;
   var cubs yearling adult_m adult_f pop r
       msprhrv fsprhrv mfallhrv ffallhrv;
   output out=means mean=cubs yearling adult_m adult_f pop r
       msprhrv fsprhrv mfallhrv ffallhrv
       stderr=se_cubs se_yearl se_ad_m se_ad_f se_pop se_r
       se_mshrv se_fshrv se_mfhrv se_ffhrv;
proc append base=library.sprgharv data=means;
run;
%mend bearsim;

options pagesize=55 linesize=78;
libname library '.';

%bearsim(spropfem=25,sprharv=10);
%bearsim(spropfem=30,sprharv=10);
%bearsim(spropfem=35,sprharv=10);
%bearsim(spropfem=40,sprharv=10);
%bearsim(spropfem=45,sprharv=10);
%bearsim(spropfem=25,sprharv=12.5);
%bearsim(spropfem=30,sprharv=12.5);
%bearsim(spropfem=35,sprharv=12.5);
%bearsim(spropfem=40,sprharv=12.5);
%bearsim(spropfem=45,sprharv=12.5);
%bearsim(spropfem=25,sprharv=15);
%bearsim(spropfem=30,sprharv=15);
%bearsim(spropfem=35,sprharv=15);
%bearsim(spropfem=40,sprharv=15);
%bearsim(spropfem=45,sprharv=15);
%bearsim(spropfem=25,sprharv=17.5);
%bearsim(spropfem=30,sprharv=17.5);
%bearsim(spropfem=35,sprharv=17.5);
%bearsim(spropfem=40,sprharv=17.5);
%bearsim(spropfem=45,sprharv=17.5);
%bearsim(spropfem=25,sprharv=20);
%bearsim(spropfem=30,sprharv=20);
%bearsim(spropfem=35,sprharv=20);
%bearsim(spropfem=40,sprharv=20);
%bearsim(spropfem=45,sprharv=20);

endsas;

%bearsim(fpropfem=30,fallharv=10);
%bearsim(fpropfem=35,fallharv=10);
%bearsim(fpropfem=40,fallharv=10);
%bearsim(fpropfem=45,fallharv=10);
%bearsim(fpropfem=50,fallharv=10);
%bearsim(fpropfem=30,fallharv=12.5);
%bearsim(fpropfem=35,fallharv=12.5);
%bearsim(fpropfem=40,fallharv=12.5);
%bearsim(fpropfem=45,fallharv=12.5);
%bearsim(fpropfem=50,fallharv=12.5);
%bearsim(fpropfem=30,fallharv=15);
%bearsim(fpropfem=35,fallharv=15);
%bearsim(fpropfem=40,fallharv=15);
%bearsim(fpropfem=45,fallharv=15);
%bearsim(fpropfem=50,fallharv=15);
%bearsim(fpropfem=30,fallharv=17.5);
%bearsim(fpropfem=35,fallharv=17.5);
%bearsim(fpropfem=40,fallharv=17.5);
%bearsim(fpropfem=45,fallharv=17.5);
%bearsim(fpropfem=50,fallharv=17.5);
%bearsim(fpropfem=30,fallharv=20);
%bearsim(fpropfem=35,fallharv=20);
%bearsim(fpropfem=40,fallharv=20);
%bearsim(fpropfem=45,fallharv=20);
%bearsim(fpropfem=50,fallharv=20);

/*
bearsim(fallharv=0,sprharv=0,doubharv=0.1,probbadb=0.1,spropfem=33,fpropfem=45)
*/
run;