%newsas(uspop);
*-- Set some graphics options for plots;
%include goptions;
goptions hsize=6in vsize=6in;
symbol1 value=dot h=2 i=join c=black;
symbol2 value=circle i=join c=red;
axis1 label=(angle=90 rotate=0);
axis2 order=(1775 to 2025 by 25);

title 'Growth of US Population, 1790-1990';
    /*------------------------------ */
data uspop;
   retain Year 1780;      * start at 1780;
   input Pop :6.3 @@;     * just read the population data;
   year=year+10;          * increment year;
   yearsq=(year-1790)**2;
   label Pop='Population';
list;datalines;
3929 5308 7239 9638 12866 17069 23191 31443 39818 50155
62947 75994 91972 105710 122775 131669 151325 179323 203211 226542 248710
281422 . .
;
proc print;
     id year;
proc gplot;
     plot Pop * year / vaxis=axis1 haxis=axis2;
 
proc reg data=uspop;
     model Pop = year;
     plot residual. * year = '*' ;
     title2 'Linear model';
 
proc reg data=uspop;
   title2 'Quadratic model';
   id     year;
   model  Pop=year yearsq /  ss1 ss2;
   plot   residual. * year = '*';
run;
*-- Add a dummy variable for war years --;
data uspop;
   set uspop;
   waryear = (1940 <= year <= 1950);
proc reg data=uspop;
   title2 'Quadratic model + War Years';
   id     year;
   model  Pop=year yearsq waryear / ss1 ss2;
   plot   residual. * year = '*';
 
*-----------------------------------------------------------
* Try fitting transformed data: perhaps there is an
  exponential relationship,
        Y = b0 * exp (b1 * X)
  If we transform the data to logs, the model becomes linear:
      log Y = log b0  + b1 * X
----------------------------------------------------------- ;
data uspop;
   set uspop;
   logpop = LOG(Pop);          * This gives natural logs (base e);
proc gplot;
     plot logpop * year / vaxis=axis1 haxis=axis2;
     Title2 'Fitting exponential model, Pop = b0 * exp (b1*year)';
 
proc reg;
     model logpop = year;
     output out=P4   p=PREDICT2 r=RESID;
	 plot logpop * year = '*'
	 	  p.  * year = '@' / overlay;
	 plot r.  * year = '*';
data Pop4;
     set P4;
     PREDICT2 = exp(PREDICT2);  /* transform back */
proc gplot;
     plot Pop* year=1
          predict2*year=2 / overlay vaxis=axis1 haxis=axis2;
run;