proc (2)=max(dvar,dloc,ivar,iloc,cloc,q,datfile,t,savname,qsav,&eval); local maxstp,minstp,maxls,minls,stol,dtol,gtol,rtol,ctol,hpr,ss2,fmt,fmt1,fmt2, mask,ttt,i,ss1,ls,lfe,ll,llsav,pcp,pnp,rws,swj,err,q0,ll0,df1, direc,h,hll,dll,f1,nr,pass,tt,dt,y,x,dparm,idf0,ss0,df0,lll,dm,ind; local eval:proc; "MAX.G: maximum likelihood estimation of conditional probability models By John Rust, University of Maryland. October, 2006"; /* Arguments: dvar : character mnemonic of dependent variable dloc : position of dependent variable ivar : character mnemonics of independent variables iloc : positions of independent variables cloc : position of constant term q : initial parameter estimates datfile: pathname of data file t : number of rows of data file read per pass through data set savname: pathname to store final printout of estimation results qsav : pathname to store final parameter estimates */ settings: /* display optimization algorithm parameters */ maxstp=10; minstp=.001; maxls=20; minls=0; stol=1e-2; dtol=1e-1; gtol=1e-1; rtol=.3; ss2=1.0; ctol=1e-7; hpr=0; format /rz 12,8; ""; "maximum stepsize ";; maxstp; "minimum stepsize ";; minstp; "maximum linesearch iterations ";; maxls; "minimum linesearch iterations ";; minls; "stop linesearch if change in ss < ";; stol; "stop linesearch if line derivative < ";; dtol; "switch to step size 1 when grad*direc<";; gtol; "Armijo-Goldstein linesearch tolerance ";; rtol; "initial stepsize ";; ss2; "convergence tolerance (grad*direc) ";; ctol; "print hessian matrix? (1=yes, 0=no) ";; hpr; ""; "dependent variable ";; $dvar; "independent variables"; if rows(ivar) < 7; $ivar'; else; i=1; do until i > rows(ivar)/7; $ivar[7*(i-1)+1:7*i]'; i=i+1; endo; if (rows(ivar)%7 > 0); i=(rows(ivar)-rows(ivar)%7)/7; $ivar[7*i+1:7*i+rows(ivar)%7]'; endif; endif; "Number of parameters ";; rows(ivar); ""; /* initialize variables used in program */ let fmt[1,3]="-*.*s " 8 8; let fmt1[1,3]="*.*lf" 16 9; let fmt2[1,3]="-*.*lf" 4 0; fmt=fmt2|fmt|fmt1|fmt1|fmt1|fmt1; mask=1~0~1~1~1~1; ttt=hsec; i=1; ss1=0; ls=0; lfe=0; ll=0; llsav=0; pcp=0; pnp=0; rws=0; swj=0; err=0; dm=rows(q); q0=zeros(dm,1); direc=ones(dm,1); dparm=direc; h=zeros(dm,dm); hll=h; /* open estimation data file for reading */ closeall; open f1=^datfile; nr=rowsf(f1); "reading from data file ";; $datfile; "number of rows of data ";; nr; pass=(nr-(nr%t))/t; nr=nr%t; "rows of data read/pass ";; t; "store output in file ";; $savname; "initial parameter values ";; q'; start: /* start of main program loop. Increment likelihood function evaluation counter, lfe */ lfe=lfe+1; core: /* core of program: evaluation of likelihood and derivatives */ tt=hsec; hll=zeros(dm,dm); dll=zeros(dm,1); ll=0; rws=0; if (i == -1); h=hll; endif; do until eof(f1); dt=readr(f1,t); ind=seqa(1,1,rows(dt)); ind=packr(ind~dt[.,dloc~iloc']); ind=ind[.,1]; if not(scalmiss(ind)); dt=dt[ind,.]; {ll,dll,h,hll}=eval(dt,dloc,iloc,cloc,q,ll,dll,ls,i,h,hll); endif; rws=rws+rows(dt); "\r";; rws;; endo; y=seekr(f1,1); if (i == 1); llsav=ll; endif; /* branch to linesearch code to maximize along search direction */ ""; if ls > 0; goto lsrch; endif; /* compute new search direction when linesearch terminates */ if hpr > 0; hll; endif; save hll; hll=invpd(-hll); if hpr > 0; hll; endif; direc=hll*dll; if ss2 > maxstp; "Current stepsize exceeds maxstp: reset ss2=1"; ss2=1; endif; q0=q; q=q0+ss2*direc; tt=hsec-tt; /* print summary of current iteration */ "Iteration";; i; i=i+1; if (i == 0); output file=^savname on; "Dependent variable";; $dvar;; " Data set ";; $datfile; " variable estimates std error t-stat u t-stat c"; else; " variable new parms direction gradient t-stat"; endif; if (i == 0); q0=sqrt(diag(hll*h*hll)); x=seqa(1,1,dm)~ivar~q~q0~(q./(sqrt(diag(hll))))~(q./q0); else; x=seqa(1,1,dm)~ivar~q~direc~dll~(q./(sqrt(diag(hll)))); endif; x=printfm(x,mask,fmt); print; "log likelihood ss=0";; ll/rws;; " grad*direc ";; idf0=dll'*direc; idf0; "evaluation number ";; lfe;; " cpu time ";; tt; ""; if ((idf0 < ctol) and (i == 0)); "percent naively predicted ";; pnp/rws; "percent correctly predicted";; pcp/rws; "total observations ";; rws; goto bye; elseif ((idf0 < ctol) and (i > 0)); i=-1; ls=0; ss2=1; goto start; elseif idf0 < gtol; ls=0; ss2=1; goto start; endif; /* save last stepsize and likelihood value and return to start of program */ ss0=0; if (i <= 3); ss2=1; endif; ss1=ss2; df0=idf0; lll=ll; ll0=lll; ls=1; goto start; lsrch: /* secant algorithm for linear search */ "linesearch";; ls;; "derivative time";; hsec-tt; "likelihood ss=";; ss2;; ll/rws; df1=dll'*direc; "df1,df0";; df1;; df0; /* compute new stepsize by secant method */ ss2=ss1-df1*(ss0-ss1)/(df0-df1); /* check to see if computed stepsize is OK, taking corrective action if necessary */ if ss2 < 0; "negative step";; ss2; ss2=ss1; ls=0; q=q0+ss2*direc; goto start; endif; if ss2 > maxstp; "step exceeds maximum";; ss2; ss2=ss1; ls=0; q=q0+ss2*direc; goto start; endif; if ls > maxls; "line search limit ss=";; ss2; ls=0; q=q0+ss2*direc; goto start; endif; /* Armijo upper bound step size test */ if (((ll <= ll0+rtol*idf0*ss2) or ss2 > maxstp) and (ls >= minls)); "Armijo upper limit ss=";; ss2; ls=0; q=q0+ss2*direc; goto start; endif; /* Goldstein lower bound step size test */ if (((ll > ll0+(1-rtol)*idf0*ss2) or ss2 < minstp) and (ls >= minls)); "Goldstein lower limit ss=";; ss2; ls=0; q=q0+ss2*direc; goto start; endif; /* Terminate linesearch when derviative is close to zero */ if (((abs(df1) < dtol) or (abs(ss1-ss2) < stol)) and (ls >= minls)); "small change in lf, ss=";; ss2; ls=0; q=q0+ss2*direc; goto start; endif; /* update parameters to new stepsize and return to start of program */ q=q0+ss2*direc; ls=ls+1; df0=df1; lll=ll; ss0=ss1; ss1=ss2; goto start; /* Termination of program */ bye: "Likelihood at initial parameter estimates ";; llsav/rws; "MAX converged: total execution time";; hsec-ttt; save ^qsav=q; ""; "parameters saved in file : ";; qsav; chrs(175)$+"edit ";; $savname; retp(q,hll); output off; closeall; endp;