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;
