%let mini=1;			*Creates the macro variable mini equal to 1.;

data one;







infile 'c:\data.dta';			*Inputs the data set data.dta where Hd=dominant height,;







input plot age Hd Ai Hdi;		* Ai=base age, and Hdi=dominant height at the base age.;

proc sort; 







by plot age;

data two;  







set one; 







by plot age;







if last.plot and last.age then do;







call symput('maxi', trim(left(plot)));	*Creates the macro variable maxi equal to the;







end;				* total number of plots.;

data three;  







set two; 







by plot age;







array p(&maxi) p1-p&maxi;	*This data step creates an array of dummy;







array S(&maxi) S1-S&maxi;	*variables (p1, p2, …, pn) as shown in;







retain S1-S&maxi;		*Table 2 and creates a horizontal vector;







p(plot)=1;			*of observed site index values (S1, S2, …, Sn),;







do i=1 to &maxi;			*as shown in Table 3, at the specified base age.;







if p(i)=. then p(i)=0;  







end;  







S(plot)=Hdi;

proc sort; 







by descending plot;

data start;







set three; 







by descending plot;		*This data step creates the temporary SAS data set start; 







keep S1-S&maxi alpha beta;	*that contains the starting values for each parameter.;







if _n_ =1;







alpha=.8 /*starting value for alpha*/;  







beta=1.4 /*starting value for beta*/;

%macro fits;			*This is the SAS macro fits that contains;







proc model data=three;		* the model fitting procedure.;







exogenous age Ai %do i=&mini %to &maxi %by 1; 







p&i %end;;







endogenous Hd;  







parms S1-S&maxi alpha beta;

Hd=(%do i=&mini %to %eval(&maxi-1) %by 1; 







p&i*S&i+ %end; 







p&maxi*S&maxi)*((1-exp(-alpha*age))/(1-exp(-alpha*Ai)))**beta;

fit Hd / converge=.0000001 estdata=start;

run;







%mend fits;

%fits;

run;







quit;

The necessary SAS code for the dummy variable fitting approach is given above.   

%macro GLOBALR(EQ,T,BASE);                      	/* macro to carry out the global fit */







  libname PARMS 'D:\PMRC\DOCS\PARMS';           /* SAS library to output parameter estimtes */







  proc datasets;







    delete WORK04;







  run;







  %if (&EQ=1) %then %do;







    proc model data=WORK02 method=MARQUARDT maxiter=100 noprint;







      parms b=0.08 c=1.4;







      TEM1=1-exp(-b*AGE);







      TEM2=1-exp(-b*&BASE);







      TEM=TEM1/TEM2;







      HD=S*TEM**c;







      fit HD / outest=WORK04;







      title "CHAPMAN-RICHARDS - Global Fit";







    run;







    data PARMS.G_&T;                            /* output parameter estimates to SAS dataset */







      set WORK04(obs=1);                        /* files referenced by iteration number */







      keep b c;







    run;







  %end;







%mend GLOBALR;

%macro BYPLOT(EQ,STRT_FAC,T,BASE);              /* macro to carry out the site specific */







  libname PARMS 'D:\PMRC\DOCS\PARMS';           /* parameter estimation */







  proc datasets;    







delete WORK06;







  run;







  %let STRT_FAC=%upcase(&STRT_FAC);             /* STRT_FAC is the stratification factor */







  %if (&EQ=1) %then %do;







    proc model data=WORK05 method=MARQUARDT maxiter=100 noprint;







      %if (&STRT_FAC^=) %then %do;







        by &STRT_FAC;







      %end;







      parms SPARM=45;







      TEM1=1-exp(-b*AGE); 







	TEM2=1-exp(-b*&BASE); 







	TEM=TEM1/TEM2;  







	HD=SPARM*TEM**c;







      fit HD / outest=WORK06;







      title "CHAPMAN-RICHARDS - By PLOT Fit";







    run;







    data PARMS.P_&T;                            		/* output estimates to SAS dataset */







      %if (&STRT_FAC^=) %then %do;              	/* files referenced by iteration number */







        set WORK06;  







        by &STRT_FAC;







        if first.&STRT_FAC then output;







        else delete;







      %end;







      %else %do;







        set WORK06(obs=1);







      %end;







      keep &STRT_FAC SPARM;







    run;







    %end;







%mend BYPLOT;

%macro FITIT(BTOLER,CTOLER,BASE);               /* Main fitting macro */







  libname PARMS 'D:\PMRC\DOCS\PARMS';

  %GLOBALR(1,0,&BASE);                          /* call global regression macro */

  data PARMS1;







  set PARMS.G_0;







  run;

  %do i=1 %to 10;

    %let j=%eval(&i-1);                         /* replace global parameters on input */







    data WORK05;                                /* dataset with new global parameters */







      set WORK02;                               /* from the fitting macro */







      if _N_=1 then set PARMS1;







      %if &i gt 1 %then %do;







        BCHK=abs(Pb-b);                         /* check parameter estimates for convergence */







        CCHK=abs(Pc-c);







        %if BCHK lt &BTOLER and CCHK lt &CTOLER %then %do;







          put 'ENDING EXECUTION';







          stop;







        %end;







      %end;







    run;

    %BYPLOT(1,SFACT,&i,&BASE);           /* call by plot regression to estimate */







                                         /* the site specific parameters */







    data WORK07;







    set PARMS.P_&i;







    run;

    proc sort;







      by SFACT;







    run;

    data WORK08;                         /* replace site index values with estimates */







      merge WORK02 WORK07;







        by SFACT;







      S=SPARM;







      drop SPARM;







    run;

    data WORK02;







      set WORK08;







    run;

    %GLOBALR(1,&i,&BASE);







    data PARM0;







      set PARMS.G_&j;







      Pb=b;







       Pc=c;







      drop b c;







    run;







    data PARM1;







      set PARMS.G_&i;







    run;







    data PARMS1;







      merge PARM0 PARM1;







    run;







  %end;

%mend FITIT;

/***************************/







/*                         */







/* data WORK02;            */







/* WORK02 dataset contains */







/* PLOT ID, AGE, Hd, S10   */







/*                         */







/***************************/







%FITIT(0.00001,0.00001,10);







quit;