Estimation Details Example for PROC MODEL

/*--------------------------------------------------------------

                    SAS Sample Library

        Name: moded.sas
 Description: Example program from SAS/ETS User's Guide,
              The MODEL Procedure
       Title: Estimation Details Example for PROC MODEL
     Product: SAS/ETS Software
        Keys: nonlinear simultaneous equation models
        PROC: MODEL
       Notes:

--------------------------------------------------------------*/

ods graphics on;

/*- Examples from the Estimation Details Section of the MODEL
    Procedure -*/

data t;  /* Sum of two normals  */
   format date monyy.;
   do t = 0 to 9.9 by 0.1;
      date = intnx( 'month', '1jun90'd, (t*10)-1 );
      y =  0.1 * (rannor(123)-10) +
            .5 * (rannor(123)+10);
      output;
   end;
run;

ods select Model.Liklhood.ResidSummary
           Model.Liklhood.ParameterEstimates;

proc model data=t time=t itprint;
   dependent y;
   parm a 5;

   y = a;
   obj = resid.y * resid.y;
   errormodel y ~ general( obj )
   cdf=(empirical=(tails=( normal percent=10)));

   fit y / outsn=s out=r;
   id  date;

   solve y / data=t(where=(date='1aug98'd))
             residdata=r sdata=s
             random=200 seed=6789 out=monte ;
run;


/*--- Generate the pdf ---*/

proc kde data=monte;
   univar y / plots=density;
run;

ods graphics off;

/*--- Troubleshooting Convergence Problems ---*/

data test;
  do i = 1 to 20;
     x = 5 * ranuni(1234);
     y = 10 + 2 * sqrt(x) + .5 * rannor(1234);
     output;
     end;
run;

proc model data=test;
   y = a + b * x ** c;
   label a = "Intercept"
         b = "Coefficient of Transformed X"
         c = "Power Transformation Parameter";
   fit y;
run;

proc model data=test;
   y = a + b * x ** c;
   label a = "Intercept"
         b = "Coefficient of Transformed X"
         c = "Power Transformation Parameter";
   fit y start=(c=5);
run;

proc model data=test;
   parms a b c;
   y = a + b * x ** c;
   label a = "Intercept"
         b = "Coefficient of Transformed X"
         c = "Power Transformation Parameter";
   fit y start=(c=1) / startiter itprint;
run;

proc model data=test;
   parms a b c;
   y = a + b * x ** c;
   label a = "Intercept"
         b = "Coefficient of Transformed X"
         c = "Power Transformation Parameter";
   fit y start=(c=1 .7 .5 .3 0) / startiter itprint;
run;

/*--- Linear Dependencies ---*/

data test3;
   do t=1 to 10;
      x1 = sqrt(t) ;
      y1 = -2 * x1 * x1  ;
      y1 = y1 + 5 * rannor(1);
      output;
   end;
run;

proc model data=test3;
   exogenous x1;
   parms b1 a1 c1 ;
   y1 = a1 * x1 + b1 * x1 * x1  + c1 * x1 * x1 *x1;
   fit y1 / collin;
run;

/*--- Iteration History ---*/

data test2;
   do t=1 to 50;
      x1 = sqrt(t) ;
      x2 = rannor(10) * 10;
      y1 = -2 * x2 * x2 - 50 / x2 - x1 * x1;
      y2 = 2* y1 + 2 * x2 * x2 + 50 / x2 + 25 * rannor(1);
      y1 = y1 + 25 * rannor(1);
      output;
   end;
run;

proc model data=test2;
   y1 = a1 * x2 * x2 - exp( d1*x1);
   y2 = a2 * x1 * x1 + b2 * exp( d2*x2);
   fit y1 y2 / itall XPX I ;
run;


/*--- Computer Resource Requirements ---*/

proc model data=test2 details;
   exogenous x1 x2;
   parms b1 100 a1 a2 b2 2.5 c2 55;
   y1 = a1 * y2 + b1 * x1 * x1;
   y2 = a2 * y1 + b2 * x2 * x2 + c2 / x2;
   fit y1 y2 / n3sls memoryuse;
   inst b1 b2 c2 x1 ;
run;

/*--- Testing for Normality ---*/

proc model data=test2;
   y1 = a1 * x2 * x2 - exp( d1*x1);
   y2 = a2 * x1 * x1 + b2 * exp( d2*x2);
   fit y1 y2 / normal ;
run;

/*--- Heteroscedasticity ---*/

/* Data from U.S. Dept. of Commerce (1979, p. 157) */

data schools;
   length state $ 3;
   input state $ exp in;
   income = 1.0e-4 * in;
   income2 = income*income;
   x0 = 1;
datalines;
ala 275 6247
ark 275 6183
con 531 8914
fla 316 7505
ida 304 6813
iow 431 7873
la  316 6640
mas 427 8063
mis 259 5736
neb 294 7391
nj  423 8818
nc 335 6607
okl 320 6951
ri 342 7526
ten 268 6489
vrt 353 6541
wva 320 6456
ala 821 10851
cal 387 8850
del 424 8604
ga 265  6700
ill 437 8745
kan 355 8001
man 327 6333
mic 466 8442
mo  274 7342
nev 359 9032
nmx 388 6505
ndk 311 7478
org 397 7839
sc  315 6242
tex 315 7697
va 356 7624
arz 339 7374
col 452 8001
dc  428 10022
haw 403 8380
ind 345 7696
ky 260 6615
md 427 8306
min 477 7847
mon 433 7051
nh 279 7277
ny 447 8267
oho 322 7812
pa 412 7733
Sdk 321 6841
uta 417 6622
was 415 8450
wyo 500 9096
;

proc model data=schools;
   parms const inc inc2;

   exp = const + inc * income + inc2 * income * income;
   incsq = income * income;

   fit exp / white breusch=(1 income incsq);
run;

data test;
   do t=1 to 25;
      y = 250 * (exp( -0.2 * t ) - exp( -0.8 * t )) +
          sqrt( 9 / t ) * rannor(1);
      output;
   end;
run;

proc model data=test;
   parms b1 0.1 b2 0.9;
   y = 250 * ( exp( -b1 * t ) - exp( -b2 * t ) );
   fit y;
run;

proc model data=test;
   parms b1 0.1 b2 0.9;
   y = 250 * ( exp( -b1 * t ) - exp( -b2 * t ) );
   fit y;
   weight t;
run;

proc model data=test;
   parms b1 0.1 b2 0.9;
   y = 250 * ( exp( -b1 * t ) - exp( -b2 * t ) );
   fit y / gmm;
   instruments b1 b2;
run;

/*--- Godfrey Output Example ---*/
data godfrey;
   e1 = 0; e2 = 0; e3 = 0;
   do i=1 to 100;
      u = rannor(123213);
      e = .5*e1 + .35*e2 + u; /* Create AR(2) error */
      y = 3 + e;
      e3 = e2;
      e2 = e1;
      e1 = e;
      output;
   end;
run;
Title 'Godfrey Test Output';
proc model data=godfrey;
   parms a b;
   y = a  + b * lag(y);
   fit y / godfrey=3;
run;
title;

/*---  Ordinary Differential Equations  ---*/
/*--- Introduction to ODEs ---*/


data t;
   time=0; output;
   time=1; output;
   time=2; output;
run;

proc model data=t ;
   dependent y 0 z 1;
   parm b -2 c -4;
      /* Solve  y''=b y' + c y --------------*/

   dert.y = z;
   dert.z = b * dert.y + c * y;

   solve y z / dynamic solveprint;
run;

/*--- Estimation of Differential Equations ---*/

data fish;
   input day conc;
datalines;
0.0    0.0
1.0    0.15
2.0    0.2
3.0    0.26
4.0    0.32
6.0    0.33
;

proc model data=fish;
   parm ku ke;

   dert.conc = ku - ke * conc;

   fit conc / time=day;
run;

proc model data=fish;
   parm ku .3 ke .3;

   dert.conc = ku - ke * conc;

   fit conc / time = day dynamic;
run;

proc model data=fish;
   parm ku .3 ke .3 conc0 0;

   dert.conc = ku - ke * conc;

   fit conc initial=(conc = conc0) / time = day dynamic;
run;

proc model data=fish;
   parm ku .3  ke .3;

   conc = (ku/ ke)*( 1 -exp(-ke * day));

   fit conc;
run;

data t2;
   y=0; time=0; output;
   y=2; time=1; output;
   y=3; time=2; output;
run;

proc model data=t2;
   independent x 0;
   dependent y;
   parm a 5;
   dert.y = a * x;
   solve x / out=goaldata;
run;

proc print data=goaldata;
run;

/*---  Hausman Specification Test  ---*/

data one;
  retain date '1mar79'd;
  date = intnx('month',date,1);
  format date monyy.;
  c = 1;
  intercep = 1;
  input y1 y2 y3 z1 x1 z3 z4 by;
  datalines;
  2.4994  0.4788   2.7918 2.5722   2.1757  0.7777  1.0735 1
  1.5564  3.4489   2.3558 4.6991   3.0314  2.5722  2.1757 1
  4.9640  4.6048   4.3040 4.8526   4.7281  4.6991  3.0314 1
  5.1814  5.6219   3.9692 5.4906   6.4729  4.8526  4.7281 1
  7.8424  4.8692   4.4436 4.6224   4.8871  5.4906  6.4729 1
  7.3664  5.8843   7.5060 6.9653   6.5458  4.6224  4.8871 1
  7.2268  7.2121   7.4833 7.4776   7.4505  6.9653  6.5458 1
  7.9623  9.5777   8.2232 0.6157   9.4022  7.4776  7.4505 1
  8.5682  9.6402  12.0452 0.2958   8.4917 10.6157  9.4022 1
  9.3103 11.6103  10.5198 8.5683  11.3056 10.2958  8.4917 1
 12.1934 10.9687  10.8176 1.8121  13.4194  8.5683 11.3056 1
 12.7798 13.3972  14.3034 3.8031  12.1540 11.8121 13.4194 1
 14.2420 15.0411  13.2027 5.5366  13.4508 13.8031 12.1540 1
 14.9286 14.5269  14.9174 6.0554  15.3293 15.5366 13.4508 1
 16.6976 14.2443  15.0952 6.7391  16.0680 16.0554 15.3293 1
 16.5999 16.3579  17.5859 6.9499  16.4259 16.7391 16.0680 1
 16.2460 16.9044  17.6270 8.6441  19.0637 16.9499 16.4259 1
 19.6290 19.2028  17.7218 9.3554  18.7674 18.6441 19.0637 1
 19.9118 21.0103  19.6083 1.8921  21.4469 19.3554 18.7674 1
 21.5054 20.6841  20.4052 0.5382  21.6456 21.8921 21.4469 1
 24.2957 23.0539  21.3546 1.2352  23.2216 20.5382 21.6456 1
 21.3397 21.3554  22.1597 2.0794  22.2892 21.2352 23.2216 1
 24.7737 24.6336  25.6975 4.3533  24.3338 22.0794 22.2892 1
 24.7168 26.2163  25.9412 5.0443  25.6809 24.3533 24.3338 1
 26.6521 25.8959  27.3117 6.3538  24.9775 25.0443 25.6809 1
 26.5674 26.2813  26.3400 5.4808  27.2362 26.3538 24.9775 1
 27.9294 28.1228  28.3723 9.0160  29.2323 25.4808 27.2362 1
 29.8447 28.9544  28.5548 0.0021  28.8931 29.0160 29.2323 1
 29.9070 30.3443  29.3153 8.8005  28.7509 30.0021 28.8931 1
 29.1622 28.0716  29.0494 1.7125  30.2890 28.8005 28.7509 1
 32.5227 32.3140  33.3892 9.9792  30.9055 31.7125 30.2890 1
 31.7739 31.6445  30.0925 2.5175  34.5828 29.9792 30.9055 1
 33.4504 34.5492  33.5850 4.0684  31.8198 32.5175 34.5828 1
 33.8820 34.8228  34.6487 5.9318  34.1898 34.0684 31.8198 1
 35.8297 35.0169  36.6375 5.5199  34.9203 35.9318 34.1898 1
 38.0351 38.1078  37.2694 7.7646  37.0174 35.5199 34.9203 1
 38.3661 37.6245  38.9162 8.9295  37.5227 37.7646 37.0174 1
 38.5417 40.6043  37.2441 0.0436  38.2555 38.9295 37.5227 1
 40.5407 39.3718  41.3623 9.9658  40.1357 40.0436 38.2555 1
 40.7421 41.1084  40.0556 0.4766  40.2052 39.9658 40.1357 1
  2.4994  0.4788   2.7918 2.5722   2.1757  0.7777  1.0735 2
  1.5564  3.4489   2.3558 4.6991   3.0314  2.5722  2.1757 2
  4.9640  4.6048   4.3040 4.8526   4.7281  4.6991  3.0314 2
  5.1814  5.6219   3.9692 5.4906   6.4729  4.8526  4.7281 2
  7.8424  4.8692   4.4436 4.6224   4.8871  5.4906  6.4729 2
  7.3664  5.8843   7.5060 6.9653   6.5458  4.6224  4.8871 2
  7.2268  7.2121   7.4833 7.4776   7.4505  6.9653  6.5458 2
  7.9623  9.5777   8.2232 0.6157   9.4022  7.4776  7.4505 2
  8.5682  9.6402  12.0452 0.2958   8.4917 10.6157  9.4022 2
  9.3103 11.6103  10.5198 8.5683  11.3056 10.2958  8.4917 2
 12.1934 10.9687  10.8176 1.8121  13.4194  8.5683 11.3056 2
 12.7798 13.3972  14.3034 3.8031  12.1540 11.8121 13.4194 2
 14.2420 15.0411  13.2027 5.5366  13.4508 13.8031 12.1540 2
 14.9286 14.5269  14.9174 6.0554  15.3293 15.5366 13.4508 2
 16.6976 14.2443  15.0952 6.7391  16.0680 16.0554 15.3293 2
 16.5999 16.3579  17.5859 6.9499  16.4259 16.7391 16.0680 2
 16.2460 16.9044  17.6270 8.6441  19.0637 16.9499 16.4259 2
 19.6290 19.2028  17.7218 9.3554  18.7674 18.6441 19.0637 2
 19.9118 21.0103  19.6083 1.8921  21.4469 19.3554 18.7674 2
 21.5054 20.6841  20.4052 0.5382  21.6456 21.8921 21.4469 2
 24.2957 23.0539  21.3546 1.2352  23.2216 20.5382 21.6456 2
 21.3397 21.3554  22.1597 2.0794  22.2892 21.2352 23.2216 2
 24.7737 24.6336  25.6975 4.3533  24.3338 22.0794 22.2892 2
 24.7168 26.2163  25.9412 5.0443  25.6809 24.3533 24.3338 2
 26.6521 25.8959  27.3117 6.3538  24.9775 25.0443 25.6809 2
 26.5674 26.2813  26.3400 5.4808  27.2362 26.3538 24.9775 2
 27.9294 28.1228  28.3723 9.0160  29.2323 25.4808 27.2362 2
 29.8447 28.9544  28.5548 0.0021  28.8931 29.0160 29.2323 2
 29.9070 30.3443  29.3153 8.8005  28.7509 30.0021 28.8931 2
 29.1622 28.0716  29.0494 1.7125  30.2890 28.8005 28.7509 2
 32.5227 32.3140  33.3892 9.9792  30.9055 31.7125 30.2890 2
 31.7739 31.6445  30.0925 2.5175  34.5828 29.9792 30.9055 2
 33.4504 34.5492  33.5850 4.0684  31.8198 32.5175 34.5828 2
 33.8820 34.8228  34.6487 5.9318  34.1898 34.0684 31.8198 2
 35.8297 35.0169  36.6375 5.5199  34.9203 35.9318 34.1898 2
 38.0351 38.1078  37.2694 7.7646  37.0174 35.5199 34.9203 2
 38.3661 37.6245  38.9162 8.9295  37.5227 37.7646 37.0174 2
 38.5417 40.6043  37.2441 0.0436  38.2555 38.9295 37.5227 2
 40.5407 39.3718  41.3623 9.9658  40.1357 40.0436 38.2555 2
 40.7421 41.1084  40.0556 0.4766  40.2052 39.9658 40.1357 2
;

proc model data=one out=fiml2;
   endogenous y1 y2;

   y1 = py2 * y2 + px1 * x1 + interc;
   y2 = py1* y1 + pz1 * z1 + d2;

   fit y1 y2 / ols 2sls hausman;
   instruments x1 z1;
run;

 /*--- Chow Tests ---*/

data exp;
   x=0;
   do time=1 to 100;
      if time=50 then x=1;
      y = 35 * exp( 0.01 * time ) + rannor( 123 ) + x * 5;
      output;
   end;
run;

proc model data=exp;
   parm zo 35 b;
      dert.z = b * z;
      y=z;
   fit y init=(z=zo) / chow =(40 50 60) pchow=90;
run;

 /*--- Profile Likelihood Confidence Intervals ---*/

data exp;
   do time = 1 to 20;
      y = 35 * exp( 0.01 * time ) + 5*rannor( 123 );
   output;
   end;
run;

proc model data=exp;
   parm zo 35 b;
      dert.z = b * z;
      y=z;
   fit y init=(z=zo) / prl=both;
   test zo = 40.475437 ,/ lr;
run;

data a;
   input inv cons;
   y = cons + inv;
   datalines;
   1  10
   2  15
   3  29
   6  51
   7  66
   ;

proc model data=a;
   endogenous cons;
   parameters a 0.9 b 0.1;
   cons = a * y + b;
   fit / ols fiml;
quit;

proc model data=a;
   endogenous cons;
   parameters a 0.9 b 0.1;
   y = cons + inv;
   cons = a * y + b;
   fit / ols fiml;
quit;

 /*---  Choice of Instruments  ---*/

proc model data=test2;
   exogenous x1 x2;
   parms b1 a1 a2 b2 2.5 c2 55;
   y1 = a1 * y2 + b1 * exp(x1);
   y2 = a2 * y1 + b2 * x2 * x2 + c2 / x2;
   fit y1 y2 / n2sls;
   inst b1 b2 c2 x1 ;
run;

 /*---  Autoregressive Moving Average Error Processes  ---*/

 /*--- univariate autoregression ---*/

data in;
  do t = 1 to 20;
     x = 5 * ranuni(1234);
     x1 = x;
     x2 = x*x;
     y = 10 + 2 * x;
     output;
     end;
run;

proc model data=in;
   parms a b p1 p2;
   y = a + b * x1 + p1 * zlag1(y - (a + b * x1)) +
         p2 * zlag2(y - (a + b * x1));
   fit y;
run;

proc model data=in;
   parms a b;
   y = a + b * x1;
   %ar( y, 2 );
   fit y;
run;

proc model data=in;
   parms a b c;
   y = a + b * x1 + c * x2;
   %ar( y, 2 )
   fit y / list;
run;

proc model data=in;
   parms a b c;
   y = a + b * x1 + c * x2;
   %ar( y, 13, , 1 12 13 )
   fit y / list;
run;

proc model data=in;
   y = a + b * x1 + c * x2;
   %ar( y, 2, m=uls )
   fit y;
run;

proc model data=in;
   y = a + b * x1 + c * x2;
   %ar( y, 2, m=cls2 )
   fit y;
run;

proc model data=in;
   parms a b c;
   y = a + b * x1 + c * x2;
   %ar( y, 5, type=v )
   fit y / list;
run;

 /*--- use IML module to simulate a MA process  ---*/
proc iml;
   phi = { 1 .2 };
   theta = { 1 .3 0 .5 };
   y = armasim( phi, theta, 0, .1, 200, 32565 );
   create madat2 from y[colname='y'];
   append from y;
quit;

title 'Maximum Likelihood ARMA(1, (1 3))';
proc model data=madat2;
   y=0;
   %ar( y, 1, ,  M=ml )
   %ma( y, 3, , 1 3,  M=ml )  /* %MA always after %AR */
   fit y;
run;
title;

proc model data=in;
   parms a b;
   y = a + b * x1;
   %ar( y, 2 );
   fit y;
run;

proc model data=in;
   parms a b;
   eq.y = y - (a + b * x1);
   %eqar( y, 2, eq.y );
   fit y;
run;

 /*---  Distributed Lag Models and the %PDL Macro  ---*/

data in;
  do t = 1 to 20;
     x = 5 * ranuni(1234);
     x1 = x;
     x2 = x*x;
     y = 10 + 2 * x;
     output;
     end;
run;

proc model data=in list;
   parms int pz;
   %pdl(xpdl,5,2);
   y = int + pz * z + %pdl(xpdl,x);
   %ar(y,2,M=ULS);
   id i;
fit y / out=model1 outresid converge=1e-6;
run;

/* Generate test data */
data gmm2;
    do t=1 to 50;
       x1 = sqrt(t) ;
       x2 = rannor(10) * 10;
       y1 = -.002 * x2 * x2 - .05 / x2 - 0.001 * x1 * x1;
       y2 = 0.002* y1 + 2 * x2 * x2 + 50 / x2 + 5 * rannor(1);
       y1 = y1 + 5 * rannor(1);
       z1 = 1; z2 = x1 * x1; z3 = x2 * x2; z4 = 1.0/x2;
       output;
    end;
run;

proc model data=gmm2 ;
   exogenous x1 x2;
   parms a1 a2 b1 2.5 b2 c2 55 d1;
   inst b1 b2 c2 x1 x2;
   y1 = a1 * y2 + b1 * x1 * x1 + d1;
   y2 = a2 * y1 + b2 * x2 * x2 + c2 / x2 + d1;

   fit y1 y2 / 3sls gmm kernel=(qs,1,0.2) outest=gmmest;

   fit y1 y2 / fiml type=gmm estdata=gmmest;
run;

proc print data=gmmest;
run;


 /*--- SDATA= input data set ---*/

proc model data=gmm2 ;
   exogenous x1 x2;
   parms a1 a2 b1 2.5 b2 c2 55 d1;
   inst b1 b2 c2 x1 x2;
   y1 = a1 * y2 + b1 * x1 * x1 + d1;
   y2 = a2 * y1 + b2 * x2 * x2 + c2 / x2 + d1;

   fit y1 y2 / 3sls gmm kernel=(qs,1,0.2)
                outest=gmmest outs=gmms;
run;

proc print data=gmms;
run;


 /*--- VDATA= input data set ---*/

proc model data=gmm2;
   exogenous x1 x2;
   parms a1 a2 b2 b1 2.5 c2 55 d1;
   inst b1 b2 c2 x1 x2;
   y1 = a1 * y2 + b1 * x1 * x1 + d1;
   y2 = a2 * y1 + b2 * x2 * x2 + c2 / x2 + d1;

   fit y1 y2 / gmm outv=gmmv;
run;

proc print data=gmmv(obs=15);
run;

proc model data=gmm2;
   exogenous x1 x2;
   parms a1 a2 b2 b1 2.5 c2 55 d1;
   inst b1 b2 c2 x1 x2;
   y1 = a1 * y2 + b1 * x1 * x1 + d1;
   y2 = a2 * y1 + b2 * x2 * x2 + c2 / x2 + d1;

   fit y1 y2 / 3sls gmm out=resid outall ;
run;

proc print data=resid(obs=20);
run;