proc RegTree depvar start end
*
*   RegTree builds a CART (classification and regression tree) for a dependent variable using
*   a set of explanatory variables.
*
*   Syntax:
*
*   @regtree( options ) depvar start end
*   # explanatory variables in regression format
*
*   Parameters:
*
*   depvar      dependent variable
*
*   start end   range to estimate [maximum allowed by variables involved]
*
*   Options:
*
*   trials=maximum number of divisions [20]
*   maxrsq=maximum R**2 [.98]
*
*   Both options are designed to limit the division process to avoid overfitting.
*
*   Output:
*
*   Each division will produce two lines of output - one anchoring the branch for the lower
*   values, the other, the higher. (This is somewhat crude, but showing the information in a
*   "tree" format would require fairly complex graphics.) Subsequent divisions of a branch
*   will be indented one space farther in. Each branch is labelled with the regressor on
*   which the split is based and the criterion, the number of observations and the mean of
*   the dependent variable for that branch.
*
*   Revision schedule:
*     June, 2002, Written by Tom Doan, Estima
*
type series depvar
type integer start end
*
option integer trials    20
option real    maxrsq    .98
*
local equation rtree
local vect[series] v
local vect[int]    vnum
local vect[int]    treepos
local rect[int]    regtable
local vect[int]    reglist
local rect         treeinfo
*
local integer      i j nreg startl endl test time biggest besttest
local integer      ngroups
local real         basess currss maxss minss ss ss1 ss2 sstest break bestbreak
local series       smpl groups
local integer      nobs1 nobs2 bestn1 bestn2
local real         mean1 mean2 bestm1 bestm2
local vect         posvalues positions
local vect[string] reglabels
local vect[labels] comparison
local real         level
local real         x count
local real         base width
local int          ntrees
*
*  Get explanatory variables
*
enter(varying) reglist
*
*  Figure out the usable range
*
inquire(regressorlist) startl>>start endl>>end
# reglist depvar
*
*  Define an equation so we can use the %eqn functions
*
equation rtree depvar
# reglist
*
compute nreg      = %eqnsize(rtree)
compute regtable  = %eqntable(rtree)
compute reglabels = %eqnreglabels(rtree)
*
dim v(nreg) vnum(nreg)
do i=1,nreg
   set v(i) startl endl = regtable(1,i){regtable(2,i)}
end do i
*
dim treeinfo(2*(trials+1),8)
dim treepos(trials+1)
*
set groups startl endl = 1
compute ngroups=1
*
do j=1,trials
*
*     Figure out which group has the highest within sum of squares. That will be the
*     group which we divide.
*
   compute maxss=0.0
   do i=1,ngroups
      set smpl startl endl = groups==i
      stats(noprint,smpl=smpl) depvar startl endl
      compute ss=%variance*(%nobs-1)
      if ss>maxss
         compute maxss=ss,biggest=i
   end do i
   if ngroups==1
       compute basess=maxss,currss=maxss
*
*
*
   compute minss=-1.0
*
*     Loop over regressors, and, within each regressor, over the data points. Test
*     the effect of dividing the group based upon that value of the regressor. Reject
*     any split which leaves zero observations in either partition. Locate the one which
*     gives the lowest sum of within sum of squares for the two partitions.
*
   do test=1,nreg
      do time=startl,endl
*
*           Skip data points which aren't in the group being divided
*
         if groups(time)<>biggest
            next
*
         compute break=v(test)(time)
         compute count=0.0
         set smpl startl endl = x=(groups==biggest.and.v(test)<break),count=count+x,x
         if count==0
            next
         stats(noprint,smpl=smpl) depvar startl endl
         compute ss1=%if(%nobs>1,%variance*(%nobs-1),0.0),nobs1=%nobs,mean1=%mean
         compute count=0.0
         set smpl startl endl = x=(groups==biggest.and.v(test)>=break),count=count+x,x
         if count==0
            next
         stats(noprint,smpl=smpl) depvar startl endl
         compute ss2=%if(%nobs>1,%variance*(%nobs-1),0.0),nobs2=%nobs,mean2=%mean
         compute sstest=ss1+ss2
         if minss<0.or.sstest<minss {
             compute minss=sstest,besttest=test,bestbreak=break
             compute bestn1=nobs1,bestn2=nobs2,bestm1=mean1,bestm2=mean2
         }
      end do time
   end do test
*
*    Split the "biggest" group. The values below the break point keep the old tag, while
*    the ones above get assigned the newer value
*
   compute ngroups=ngroups+1
   set groups startl endl = %if(groups<>biggest.or.v(besttest)<bestbreak,groups,ngroups)
*
*    Figure out the positioning for the partitions, and save the observation counts and
*    group means.
*
   if ngroups==2 {
      compute ntrees=0
      compute base=.50,width=.50,level=0
   }
   else {
      compute base =treeinfo(treepos(biggest),1),width=treeinfo(treepos(biggest),2)
      compute level=treeinfo(treepos(biggest),3)
   }
   compute treepos(biggest)=ntrees+1,treepos(ngroups)=ntrees+2
   compute treeinfo(ntrees+1,1)=base+.5*width
   compute treeinfo(ntrees+2,1)=base-.5*width
   compute treeinfo(ntrees+1,2)=.5*width
   compute treeinfo(ntrees+2,2)=.5*width
   compute treeinfo(ntrees+1,3)=level+1
   compute treeinfo(ntrees+2,3)=level+1
   compute treeinfo(ntrees+1,4)=besttest
   compute treeinfo(ntrees+2,4)=besttest
   compute treeinfo(ntrees+1,5)=bestbreak
   compute treeinfo(ntrees+2,5)=bestbreak
   compute treeinfo(ntrees+1,6)=1
   compute treeinfo(ntrees+2,6)=2
   compute treeinfo(ntrees+1,7)=bestn1
   compute treeinfo(ntrees+2,7)=bestn2
   compute treeinfo(ntrees+1,8)=bestm1
   compute treeinfo(ntrees+2,8)=bestm2
   compute ntrees=ntrees+2
   compute currss=currss+minss-maxss
   if currss/basess<(1-maxrsq)
      break
end do j
*
compute comparison=||"<",">="||
dim posvalues(ntrees)
ewise posvalues(i)=treeinfo(i,1)
compute positions=%ranks(posvalues)
do i=1,ntrees
   do j=1,ntrees
      if positions(j)==i {
         display @fix(treeinfo(j,3)) reglabels(fix(treeinfo(j,4))) comparison(fix(treeinfo(j,6))) treeinfo(j,5) $
            'Obs=' fix(treeinfo(j,7)) 'Mean=' treeinfo(j,8)
         break
      }
   end do j
end do i
end regtree

