/* This data is taken from Greene, Econometrics Analysis, 4th ed.,
   Appendix Table A7.2.  The nonlinear model estimated is specified
   in Example 10.3, p. 421.  */

Data YC;
 input year y c;
cards;
1950  791.8  733.2
1951  819.0  748.7
1952  844.3  771.4
1953  880.0  802.5
1954  894.0  822.7
1955  944.5  873.8
1956  989.4  899.8
1957 1012.1  919.7
1958 1028.8  932.9
1959 1067.2  979.4
1960 1091.1 1005.1
1961 1123.2 1025.2
1962 1170.2 1069.0
1963 1207.3 1108.4
1964 1291.0 1170.6
1965 1365.7 1236.4
1966 1431.3 1298.9
1967 1493.2 1337.7
1968 1551.3 1405.9
1969 1599.8 1456.7
1970 1688.1 1492.0
1971 1728.4 1538.8
1972 1797.4 1621.9
1973 1916.3 1689.6
1974 1896.6 1674.0
1975 1931.7 1711.9
1976 2001.0 1803.9
1977 2066.6 1883.8
1978 2167.4 1961.0
1979 2212.6 2004.4
1980 2214.3 2000.4
1981 2248.6 2024.2
1982 2261.5 2050.7
1983 2334.6 2145.9
1984 2468.4 2239.9
1985 2509.0 2312.6
;

/* Use OLS to get starting values for Nonlinear Least Squares.  */

proc reg data=yc;
  model c = y;

run;

/* Use the Gauss-Newton Method and LS estimates as starting point to
   obtain Nonlinear LS Estimates  */

proc nlin data=yc maxiter=150;
  parms alpha=11.1458 beta=0.89853 gamma=1.0;
  model c = alpha + beta*y**gamma;
  der.alpha=1;
  der.beta=y**gamma;
  der.gamma=beta*(y**gamma)*log(y);
  output out=ycfit p=predict r=resid;
  run;

/* Plot Fitted vs. Actual Values  */

proc plot data=ycfit;
  plot c*year predict*year='p' /overlay;
  plot resid*year / vref=0;
  run;

/* Use several starting points and the Gauss-Newton Method.  Using several
   starting points allows us to check the robustness of the NLS result.  */

proc nlin data=yc maxiter=150 best=8;
  parms alpha=6 to 16 by 2 beta=0.5 to 1.5 by 0.5 gamma=0 to 2 by 0.5;
  model c = alpha + beta*y**gamma;
  der.alpha=1;
  der.beta=y**gamma;
  der.gamma=beta*(y**gamma)*log(y);
  output out=ycfit p=predict r=resid;
  run;

/* Plot Fitted vs. Actual Values  */

proc plot data=ycfit;
  plot c*year predict*year='p' /overlay;
  plot resid*year / vref=0;
  run;

/* If we leave out the derivatives in proc nlin the method of estimation
   is by the secant (DUD) method.  See the SAS/STAT guide for further
   information on this method.  Let's see if we get roughly the same
   results using the DUD method as we got using the Gauss-Newton Method. */

proc nlin data=yc maxiter=200;
  parms alpha=11.1458 beta=0.89853 gamma=1.0;
  model c = alpha + beta*y**gamma;
  output out=ycfit p=predict r=resid;
  run;

proc plot data=ycfit;
  plot c*year predict*year='p' /overlay;
  plot resid*year / vref=0;
  run;

proc nlin data=yc maxiter=200 best=8;
  parms alpha=6 to 16 by 2 beta=0.5 to 1.5 by 0.5 gamma=0 to 2 by 0.5;
  model c = alpha + beta*y**gamma;
  output out=ycfit p=predict r=resid;
  run;

proc plot data=ycfit;
  plot c*year predict*year='p' /overlay;
  plot resid*year / vref=0;
  run;