data grow;

  do t = 1 to 70;
  y=85*exp(-4*exp(-0.05*t)) + rannor(1234);
  t2=t**2;
  t3=t**3;
  t4=t**4;
  output;
  end;

  run;

data grow2;
  set grow;
  dlog = log(y) - log(lag(y));
  if dlog le 0 then delete;
  ldlog = log(dlog);

proc print data=grow;
  var y t t2 t3 t4;

  run;

proc reg data=grow2;
    model ldlog = t t2 t3 t4;
    test t2, t3, t4;
    model ldlog = t t2 t3;
    test t2, t3;
    model ldlog = t t2;
    test t2;

    run;

proc nlin data=grow maxiter=150;
  parms alpha=80 beta=3 gamma=0.04;
  model y = alpha*exp(-beta*exp(-gamma*t));
  der.alpha=exp(-beta*exp(-gamma*t));
  der.beta=alpha*(-exp(-gamma*t))*exp(-beta*exp(-gamma*t));
  der.gamma=alpha*exp(-beta*exp(-gamma*t))*(-beta)*exp(-gamma*t)*(-t);
  output out=yfit p=predict r=resid;
  run;

/* Plot Fitted vs. Actual Values  */

proc plot data=yfit;
  plot y*t predict*t='p' /overlay;
  plot resid*t / vref=0;
  run;

/*  You get the same answer without providing proc nlin the gradient
    vector.
  proc nlin data=grow maxiter=150;
  parms alpha=80 beta=3 gamma=0.04;
  model y = alpha*exp(-beta*exp(-gamma*t));
  run;
*/

  proc iml;

/*  This segment of the program calculates the covariance matrix
    of the ML estimates, given the correlation matrix of the ML
    estimates.  */

  alpha=84.43776623;
  beta=4.14219561;
  gamma=0.05138937;

  sealpha=0.91928868466;
  sebeta=0.08852154155;
  segamma=0.001109947112;

  scale = {0.91928868466 0 0,
           0 0.08852154155 0,
           0 0 0.001109947112};

  correst = {1 -0.656687013 -0.90743212,
             -0.656687013 1 0.8864233294,
             -0.90743212 0.8864233294 1};

  covest = scale*correst*scale;

  print covest;

/*  This segment of the program calculates the standard error of
    tau for the Gompertz curve.  */

  tau = log(beta)/gamma;

  grad = {0, 4.6968, -537.9442};

  setau = sqrt(grad`*covest*grad);

  print tau setau;


/*  This segment of the program calculates the standard error of y
    at time t = 75 using the Gompertz curve.  */ 

  t = 75;

  y75 = alpha*exp(-beta*exp(-gamma*t));

  deralpha=exp(-beta*exp(-gamma*t));
  derbeta=alpha*(-exp(-gamma*t))*exp(-beta*exp(-gamma*t));
  dergamma=alpha*exp(-beta*exp(-gamma*t))*(-beta)*exp(-gamma*t)*(-t);

  print deralpha derbeta dergamma;

  g = {0.9159668, -1.638919, 509.15422};

  sey75 = sqrt(g`*covest*g);

  y75u = y75 + 1.96*sey75;
  y75l = y75 - 1.96*sey75;

  print y75 sey75 y75u y75l;

  run;