|
|
OPTIONS NOCENTER;
TITLE 'Nonlinear least-squares estimation procedure';
data traject;set thing.panel80;
if (partyid1 ge 5) then vote1=0;else vote1=1;
vote2=vote-1;if (vote2 gt 1) then delete;
if ((vote1=1) and (vote2=0)) then poltype=1;else poltype=2;
*if (partyid1 ge 5) then context=prepcont;
*if (partyid1 le 3) then context=pdemcont;
*if (partyid1=4) then delete;
if (rfriends=1) then neighbor=-3;else neighbor=.;
if (mfriends=1) then neighbor=0;
if (dfriends=1) then neighbor=3;
RAEsF=1-((prepcont**2)+(pdemcont**2));
convRAE=((prepcont**2)+(pdemcont**2));
context1=(RAEsF-0.40)/(0.81-0.40);
context2=((convRAE-0.19)/(0.59-0.19));
**** selecting only mixed friends people;
*if mfriends=1;
**** context choice;
context=pdemcont;
**** selecting only carter to reagan switchers;
*if poltype=1;
****rescales data;
array b reafeel1 reafeel3 carfeel1 carfeel3 dempart1 dempart3;
do over b;b=b/100;end;
PIDvote=1-((partyid1-1)/6);
***deletes problem data***;
prob=1;
array a vote1 vote2 reafeel1 reafeel3 carfeel1 carfeel3
partyid1 dempart1 dempart3;
do over a;if (a=.) then prob=.;end;
if prob=1;
**************************;
*** variable transforms;
changeC=carfeel3-carfeel1;
changeR=reafeel3-reafeel1;
changeD=dempart3-dempart1;
**************************;
**** renaming some variables for IML use;
l=carfeel1;
lnext=carfeel3;
r=reafeel1;rnext=reafeel3;
c=dempart1;cnext=dempart3;
n=0;nnext=0;
Proc Standard OUT=traject MEAN=0 STD=1;
VAR context;
Proc Freq;
Tables vote1*vote2;
Proc Means Data=traject;
*Proc Contents;
Proc Reg;
model changeC=changeR changeD partyid1;
model changeC=changeR changeD;
Proc Reg;
model carfeel3=carfeel1 changeR changeD partyid1;
model carfeel3=carfeel1 changeR changeD;
PROC IML;
RESET;
USE traject;
READ ALL;
SHOW NAMES;
*quit;
eligave=(1);
eligtot=sum(eligave);
pop=1;
weight=eligave/eligtot;
*totpop=eligtot;
totpop=NROW(r);
RSQDEV=pop#((rnext-r)##2);
RSSDEV=RSQDEV(|+,|);
LSQDEV=pop#((lnext-l)##2);
LSSDEV=LSQDEV(|+,|);
CSQDEV=pop#((cnext-c)##2);
CSSDEV=CSQDEV(|+,|);
NSQDEV=pop#((NNEXT-N)##2);
NSSDEV=NSQDEV(|+,|);
***To remember the initial conditions;
R0=R;C0=C;L0=L;N0=N;
*quit;
* In this program, there are nine state variables;
* They are r, rnext, c, cnext, l, lnext, n, nnext;
* Also, population and conditioning variables may be needed;
*******************************************************;
*** THIS COUNTS THE NUMBER OF CASES ****;
CASES=NROW(R);
PRINT CASES rssdev lssdev cssdev;
****************************************;
zing=0;print zing;
x1=uniform(zing);print x1;
*STOP;
*******************************************************;
START;
*This section sets conditioning variables;
DO C44 = 4 TO 4 ;
IF C44 = 1 THEN DO;
D1 = PROT;
END;
IF C44 = 2 THEN DO;
D1 = CATH;
END;
IF C44 = 3 THEN DO;
D1 = URBAN;
END;
IF C44 = 4 THEN DO;
D1 = context;
END;
*****************************************************;
****INITIAL PARAMETER GUESSES************************;
*****************************************************;
mag1=0.50;
pp=J(30,1,0);
Do m88=1 to 30;
pp[m88]=UNIFORM(x1);
* pp[m88]=(2#(UNIFORM(x1)))-1;
pp[m88]=mag1#pp[m88];
End;
*** special cases for select parameters;
pp[5]=6;
pp[6]=0.5;
pp[7]=1;
pp[8]=1;
pp[9]=0.05;
pp[11]=1;
pp[12]=0;
*** put unconditioned parameter estimates here;
* pp[1] = 0.26724;
* pp[2] = -0.38401;
* pp[3] = 0.19066;
* pp[4] = -0.38721;
* pp[5] = 5.98186;
* pp[6] = 0.83441;
* pp[7] = 1.03085;
* pp[8] = 0.88157;
* pp[9] = 0.41774;
* pp[10]= -0.13497;
* pp[11]= 1.24597;
*print pp;
*stop;
*****************************************************;
* Some necessary and useful numbers;
Y99 = 1;GCOUNT=0;
I = 0.00001 ;
h = 0.1 ;
z1=J(30,1, 0.001);
p=J(30,1,0);
a=J(30,1,0);
parmfit=J(30,1,0);
csp=J(30,1,0);
sim=J(30,1,0);
chisq=J(30,1,0);
* pnum is the total number of parameters in the model;
pnum=12;
* parmeq identifies which equation each parameter is in;
parmeq={1,1,2,2,3,3,3,3,3,3,3,3,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0};
TESTLIST=0;
MWISHO = 0;
ITCOUNT = 0 ;
TE = 0 ;
TE2=0;
TIMEUP=0;
E1 = 0 ;
DP1 = 0 ;
* Make typedata=1 for proportions and 2 for individuals;
typedata=2;
D2 = 1 ;
D3 = 1 ;
**********************************************************;
***** IF D3=1 YOU GET THE UNCONDITIONAL MODEL *******;
***** IF D3=0 YOU GET THE CONDITIONED MODEL *******;
IF D3 = 1 THEN DO;
D1 = 1 ;
D2 = 0 ;
p=pp;
Do m88=1 to 30;
pp[m88]=0;
END;END;
IF D3=0 THEN DO;
*Use the following assignments if you want random conditioned
initial guesses;
goto norand2;
MAG2=0.30;
Do m88=1 to 30;
* p[m88]=UNIFORM(x1);
p[m88]=(2#(UNIFORM(x1)))-1;
p[m88]=mag2#p[m88];
End;
norand2:
END;
**********************************************************;
GOTO ESTIMATE;
**********************************************************;
**** THE FOLLOWING SUBROUTINE IS A FOURTH-ORDER
**** RUNGE-KUTTA ALGORITHM *****;
**********************************************************;
MODELFIT:
IF DP1=1 THEN DO;
PRINT 'MODELFIT BEGINNING';
END;
R=R0;C=C0;L=L0;N=N0;
DO U1=1 to 10;
time=u1;
if (u1=1) then do;timeone=1;end;
m2=R;m1=C;m3=L;
LINK EQS;
RK1=h#DRDT;CK1=h#DCDT;LK1=h#DLDT;
m2=R+(.5#RK1);m1=C+(.5#CK1);m3=L+(.5#LK1);
LINK EQS;
RK2=h#DRDT;CK2=h#DCDT;LK2=h#DLDT;
m2=R+(.5#RK2);m1=C+(.5#CK2);m3=L+(.5#LK2);
LINK EQS;
RK3=h#DRDT;CK3=h#DCDT;LK3=h#DLDT;
m2=R + RK3;m1=C + CK3;m3=L + LK3;
LINK EQS;
RK4=h#DRDT;CK4=h#DCDT;LK4=h#DLDT;
RNEW=R+((1/6)#(RK1+(2#RK2)+(2#RK3)+RK4));
CNEW=C+((1/6)#(CK1+(2#CK2)+(2#CK3)+CK4));
LNEW=L+((1/6)#(LK1+(2#LK2)+(2#LK3)+LK4));
NNEW=N;
R=RNEW;C=CNEW;L=LNEW;
N=NNEW;
IF E1=1 THEN link printraj;
END;
if (typedata=1) then link propdata;
if (typedata=2) then link indidata;
link compile1;
link rsq;
IF DP1=1 THEN DO;
PRINT 'MODELFIT ENDING';
END;
RETURN;
******************************************;
***** HERE IS WHERE THE ACTUAL MODEL **;
***** EQUATIONS ARE PLACED. ****;
EQS:
* THE VARIABLES ARE M1(C), M2(R), and M3(L);
link special1;
drdt = (pp[1]+(p[1]#D1))+((pp[2]+(p[2]#D1))#m2);
dcdt = (pp[3]+(p[3]#D1))+((pp[4]+(p[4]#D1))#m1);
dldt = (pp[9]+(p[9]#D1))#((((pp[5]+(p[5]#D1))
#((pp[10]+(p[10]#D1))+w2)
#(1-(pp[11]+(p[11]#D1))#w2)#(m3-(pp[6]+(p[6]#D1))))
-((pp[7]+(p[7]#D1))
#((m3-(pp[6]+(p[6]#D1)))##3))
+((pp[8]+(p[8]#D1))#(w1+(pp[12]+(p[12]#D1))))) - m3);
RETURN;
special1:
w1=(1+m1-m2)/2;
w2=(m2+m1)/2;
return;
******************************************;
* Prints the variable trajectories;
printraj:
PRINT 'TRAJECTORY FOR PREDICTED R, C, AND L';
R1=eligave#R;C1=eligave#C;L1=eligave#L;
TRAJECT=R1||(C1||L1);
TRAJ=TRAJECT(|+,|);
TRAJMEAN=TRAJ/TOTPOP;
PRINT TRAJMEAN;
Return;
*******************************************;
*** Computes predicted values for use in chi square routine;
propdata:
R2=eligave#R;C2=eligave#C;L2=eligave#L;N2=eligave#N;
PRED=R2||(C2||(L2||N2));
PREDIC=PRED(|+,|);
PREDICT=PREDIC/TOTPOP;
IF DP1=1 THEN DO;
PRINT 'PREDICTED RIGHTIST, CENTRIST, LEFTIST, AND NONVOTER MEANS';
PRINT PREDICT;
END;
Return;
indidata:
number=nrow(R);
Rpredvot=J(number,1,0);
Lpredvot=J(number,1,0);
count=1;
do until (count>number);
if (L[count] <= R[count]) then Rpredvot[count]=1;
else Rpredvot[count]=0;
if (L[count] > R[count]) then Lpredvot[count]=1;
else Lpredvot[count]=0;
count=count+1;
end;
Rvotes=sum(Rpredvot);
Lvotes=sum(Lpredvot);
PREDIC=(Rvotes||Lvotes);
PREDICT=(Rvotes||Lvotes);
IF DP1=1 THEN DO;
PRINT 'PREDICTED NUMBER OF RVOTES AND LVOTES';
PRINT PREDICT;
END;
Return;
******************************************;
****** THE FOLLOWING ARE THE SUBROUTINES
****** FOR THE R-SQUARES AND THE SAVING
****** OF THE BEST PARAMETER ESTIMATES;
******************************************;
compile1:
if timeone=1 then do;
resr=0#pop;resc=0#pop;resl=0#pop;resn=0#pop; timeone=0;end;
resr=(pop#((rnext-rnew)##2));
resc=(pop#((cnext-cnew)##2));
resl=(pop#((lnext-lnew)##2));
resn=(pop#((nnext-nnew)##2));
return;
RSQ:
RESIDR=sum(resr);
RESIDL=sum(resl);
RESIDC=sum(resc);
RESIDN=sum(resn);
RRSQUARE=1-(RESIDR/RSSDEV);
LRSQUARE=1-(RESIDL/LSSDEV);
CRSQUARE=1-(RESIDC/CSSDEV);
NRSQUARE=0;*NRSQUARE=1-(RESIDN/NSSDEV);
*RSQUARE=(RRSQUARE+LRSQUARE+CRSQUARE+NRSQUARE)/4;
RSQUARE=(RRSQUARE+CRSQUARE+LRSQUARE)/3;
*RSQUARE=(RRSQUARE+CRSQUARE)/2;
IF DP1=1 THEN DO;
*PRINT 'SQUARED, SUMMED, AND WEIGHTED RESIDUALS';
*PRINT 'RESIDR RESIDC RESIDL RESIDN';
*SQRESID=RESIDR||(RESIDC||(RESIDL||RESIDN));
*PRINT SQRESID;
PRINT 'RRSQUARE CRSQUARE LRSQUARE NRSQUARE RSQUARE';
RSQFITS=RRSQUARE||(CRSQUARE||(LRSQUARE||(NRSQUARE||RSQUARE)));
PRINT RSQFITS;
END;
RETURN;
**************************************************************;
BESTPAR:
bestp=p;
BESTRSQ = RSQUARE;
BESTRRSQ = RRSQUARE;
BESTCRSQ = CRSQUARE;
BESTLRSQ = LRSQUARE;
PARMS = p;
IF DP1 = 1 THEN DO;
PRINT , PARMS;
END;
RETURN;
********************************************************;
SURFACE:
* DP1 = 0 ;
PRINT 'SURFACE' 'BEGINNING';
PRINT 'SURFACE' 'BEGINNING';
Do m88=1 to pnum;
p[m88]=p[m88] - I;
link modelfit;
link putfit;
p[m88]=p[m88] + (2#I);
link modelfit;
If (parmeq[m88]=1) then do;parmfit[m88]=(rrsquare - fit2);end;
If (parmeq[m88]=2) then do;parmfit[m88]=(crsquare - fit3);end;
If (parmeq[m88]=3) then do;parmfit[m88]=(lrsquare - fit4);end;
p[m88]=p[m88]-I;
End;
PARTIALS = PARMFIT / ( 2 # I );
TEST = SSQ(PARTIALS);
TESTLIST=(TESTLIST//TEST);
PRINT test, p, PARTIALS;
LINK MODELFIT;
PRINT 'RRSQUARE' 'CRSQUARE' 'LRSQUARE' 'RSQUARE';
RSQFITS=RRSQUARE||(CRSQUARE||(LRSQUARE||RSQUARE));
PRINT , RSQFITS;
PRINT 'SURFACE' 'ENDING';
DP1 = 0 ;
RETURN;
PUTFIT:
FIT1 = RSQUARE;
FIT2 = RRSQUARE;
FIT3 = CRSQUARE;
FIT4 = LRSQUARE;
RETURN;
PUTPRED:
PREDIC1 = PREDIC;
PREDS1 = RSQUARE;
RETURN;
CSFITS:
if (typedata=1) then do;
PREDIC1=PREDIC1+1;PREDIC=PREDIC+1;end;
CHSQALL1 = ((PREDIC - PREDIC1)##2) / abs(PREDIC1);
CHSQALL = SUM(CHSQALL1);
SIMALL = PREDS1 - RSQUARE;
PRINT CHSQALL SIMALL;
RETURN;
ESTIMATE:
DP1 = 1;TE = 1;E1 = 0;
print p;
LINK MODELFIT;TE = 0;E1 = 0;
* STOP;
* GOTO CHISQUAR;
* GOTO TATA;
LINK BESTPAR;
DP1 = 0 ;
DO Y99 = 1 TO 30;
ZCOUNT = 0 ;
ITCOUNT = 0 ;
LINK SURFACE;
* DP1 = 1;
Z2 = Z1;
* print test, z2;
PRINT 'ESTIMATE' 'BEGINNING';
E1 = 0 ;
LINK MODELFIT;
E1 = 0 ;
BEGIN:
FIT1 = RSQUARE;
NEWPARM = p + ( PARTIALS # Z2 );
p = NEWPARM;
LINK MODELFIT;
IF RSQUARE > FIT1 THEN DO;
LINK BESTPAR;
END;
ITCOUNT = ITCOUNT + 1 ;
IF RSQUARE > FIT1 THEN GOTO BEGIN;
IF ZCOUNT<2 THEN DO; PRINT ITCOUNT;END;
ITCOUNT = 0 ;
ZCOUNT = ZCOUNT + 1 ;
* PRINT , ZCOUNT;
Z2 = Z2 / 10 ;
IF ZCOUNT > 4 THEN GOTO JUMP1;
p=bestp;
GOTO BEGIN;
JUMP1:
p=bestp;
END;
PRINT 'ESTIMATE' 'ENDING';
TE = 1 ;LINK MODELFIT;TE = 0;
GOTO CHISQUAR;
*********************************************************;
CHISQUAR:
TE2=1;
do m88=1 to pnum;
csp[m88] = p[m88];LINK MODELFIT;LINK PUTPRED;
p[m88] = 0 ;LINK MODELFIT;LINK CSFITS;
CHISQ[m88] = CHSQALL;SIM[m88]=SIMALL;p[m88] = csp[m88];
end;
**** NOW THE OVERALL MODEL CHI SQUARE;
csp = p;LINK MODELFIT;LINK PUTPRED;
do m88=1 to 11;
p[m88] = 0 ;end;
LINK MODELFIT;LINK CSFITS;
CSMODEL = CHSQALL;SIMMODEL=SIMALL;p = csp;
*********************************************************;
**** PREPARATION FOR OUTPUT ****************************;
*********************************************************;
DP1 = 1 ;
LINK MODELFIT;
DP1 = 0 ;
ALLFIT1 = ( RRSQUARE || (CRSQUARE||(LRSQUARE||RSQUARE)));
SYSPAR = bestp;
SYSCHI = chisq;
SYSSIM = sim;
ALLFITS = SHAPE( ALLFIT1 ,0 , 1 );
SYSCHISQ = SHAPE( SYSCHI ,0 , 1 );
SYSSIMS = SHAPE( SYSSIM ,0 , 1 );
SYSPARMS = SHAPE( SYSPAR ,0 , 1 );
SYSEST = (SYSPARMS||(SYSCHISQ||SYSSIMS));
PRINT , C44 , SYSPARMS , SYSCHISQ , SYSEST , TESTLIST;
PRINT CSMODEL, SIMMODEL;
DP1=1;
MWISHO = 1;
LINK MODELFIT;
DP1=0;
STUFF={PARMS CHISQ SIMON};
ROW = {'p1' 'p2' 'p3' 'p4' 'p5' 'p6' 'p7' 'p8' 'p9' 'p10'
'p11' 'p12' 'p13' 'p14' 'p15' 'p16' 'p17' 'p18' 'p19' 'p20'
'p21' 'p22' 'p23' 'p24' 'p25' 'p26' 'p27' 'p28' 'p29' 'p30'};
CREATE BETAS FROM SYSEST (|COLNAME=STUFF ROWNAME=ROW|);
APPEND FROM SYSEST (|ROWNAME=ROW|);
CLOSE BETAS;
FIT={'FITS'};
ROWW = {'eq1' 'eq2' 'eq3' 'ave'};
CREATE FITS FROM ALLFITS (|COLNAME=FIT ROWNAME=ROWW|);
APPEND FROM ALLFITS (|ROWNAME=ROWW|);
CLOSE FITS;
END;
FINISH;RUN;
QUIT;
PROC PRINT DATA=BETAS;VAR ROW PARMS CHISQ SIMON;
PROC PRINT DATA=FITS;VAR ROWW FITS;
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
