program ssmlmnp c ported to the sun by vah, Feb 91 c program statement also added c c vah 18:55:22.11, Thu, Mar 7, 1991 c c 1. lower case c 2. fix '$include' statements c 3. replace suitably 'call getarg' to 'call parsecl' c 4. correct '\' in format statements c 5. comment out the '$large' statement c c-----sml.for----------------------------------------------------------- c c ===================================================== c **** multinomial multiperiod probit model **** c * estimation by smooth simulated maximum likelihood * c ===================================================== c c (c) axel boersch-supan version 2.1 1. september 1990 c c based on: - geweke, 1989 c - hajivassiliou and mcfadden, 1990 c.pl100/.rm72/ c----------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' real*4 data(maxbuf) integer*4 iseed character*8 pname(maxnpp) character*20 runtag c c the following list includes all named common blocks c common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns, * nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp common / dmode / method,model,modran,modar1,modmnp common / doffs / ns1,ns2,ns3,ns4,ns5,maps(maxalt,maxalt), * locs(maxnpp) common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp) common / dname / pname common / domeg / omega(maxnt,maxnt),domega(maxnt,maxnt,0:maxns) common / doptm / mchk0,mchk1,iter1,accu1,ialg1, * mchk2,iter2,accu2,ialg2,nwesml common / dsimu / nsimul common / dseed / iseed common / dtag / runtag common / ddata / data common / dfreq / wtsum,share(maxnt),wshare(maxnt), * amean(maxndm),wmean(maxndm) common / dxvar / x(maxnt,maxnpp),ic(maxper),wt common / ddiff / idiff2 common / dival / ival common / doutp / kanal0,ispool,ipause c c reduce stack size c common / dfunc / em1(maxnt1),em2(maxnt1),sg1(maxnt1,maxnt1), * dem1(maxnt1,maxnp),dem2(maxnt1,maxnp), * dsg1(maxnt1,maxnt1,maxnp),da(maxnp),dalf(maxnp) common / dzbar / zquer(maxnt1),covz(maxnt1,maxnt1), * dzquer(maxnt1,maxnb),dcovz(maxnt1,maxnt1,maxns) common / dopt1 / parm(maxnp),grad(maxnp),hess(maxnp,maxnp), * hinv(maxnp,maxnp) common / dopt9 / dum1(maxnp),dum2(maxnp),dum3(maxnp),dum4(maxnp), * dum5(maxnp,maxnp),dum6(maxnp,maxnp),dum7(maxnp) c c initialize size of blank common c common calfa(maxalt,maxalt),dcalfa(maxalt,maxalt,0:maxns), * cnunu(maxalt,maxalt),dcnunu(maxalt,maxalt,0:maxns), * omegax(maxnt,maxnt) c c write (*,1) maxalt,maxper,maxnt,maxnt1, * maxnb,maxns,maxnp,maxnpp, * maxbuf 1 format( *//' ============================================================' * /' **** multinomial multiperiod probit model ****' ccc * /' version 2.1 [sep 01, 1990] (c) axel boersch-supan' * /' version 2.1 - sep 01, 1990 (c) axel boersch-supan' * /' ============================================================' * /' maxalt=',i3,' maxper=',i3,' maxnt=',i3,' maxnt1=',i3 * /' maxnb =',i3,' maxns =',i3,' maxnp=',i3,' maxnpp=',i3 * /' maxbuf=',i6 * /' ============================================================' * /' *** initialization, parameter and data input ***' * /' ============================================================') call inparm(parm,np,*998) call indata(*998) c c ****************************** c * omega derivative check * c ****************************** c 100 write (*,101) 101 format( *//' ===================================================' * /' multinomial multiperiod probit model:' * /' omega derivative check' * /' ===================================================') if (ipause.lt.0) then eps=0.000001 if (mchk0.eq.0) goto 200 else write (*,*) 'enter eps [1-9: 3=0.001 etc.; 0=no check] ?' read (*,'(i1)') mchk0 if (mchk0.eq.0) goto 200 if (mchk0.gt.0) eps=0.1**mchk0 end if c call dcomeg(parm,np,eps) c c **************************************************** c * initial probit/logit estimation (clark approx) * c **************************************************** c 200 if (ipause.lt.0.and.method.eq.0) goto 300 201 write (*,202) 202 format( *//' ===============================================' * /' multinomial multiperiod probit model:' * /' initial estimation' * /' ===============================================') if (ipause.ge.0) then write (*,*) * ' enter initial estimation method [0=none, 1=logit, 2=clark] ?' read (*,'(i1)') method end if if (method.le.0.or.method.ge.4) goto 300 if (method.eq.1) write (*,203) if (method.eq.2) write (*,204) 203 format( *//' ========================' * /' scaled multinomial logit' * /' ========================') 204 format( *//' ==========================' * /' clark approximation probit' * /' ==========================') if (ispool.ne.0) write (kanal0,202) if (ispool.ne.0.and.method.eq.1) write (kanal0,203) if (ispool.ne.0.and.method.eq.2) write (kanal0,204) c c **** derivative check c ival=0 idiff2=2 if (ipause.ge.0) then write (*,*) ' derivative check of initial estimator', * ' [0=none, 1=func, 2=check] ?' read (*,'(i1)') mchk1 end if if (mchk1.ne.0) then if (ispool.ne.0) write (kanal0,'(//a)') * ' derivative check of initial estimator' call func(parm,np,fval,*998) write (*,'(a,f15.7)') ' fval=',fval if (ispool.ne.0) write (kanal0,'(a,f15.7)') ' fval=',fval if (mchk1.eq.2) then call deriv(parm,np,fval,grad,hess,idiff2,*998) write (*,'(a,30f10.4)') ' deriv=',(grad(i),i=1,np) if (ispool.ne.0) write (kanal0,'(a,30f10.4)') * ' deriv=',(grad(i),i=1,np) call fdiff(parm,np,fval,grad,*998) write (*,'(a,30f10.4)') ' fdiff=',(grad(i),i=1,np) if (ispool.ne.0) write (kanal0,'(a,30f10.4)') * ' fdiff=',(grad(i),i=1,np) end if call snooze(ipause) end if c c **** maximize likelihood c mx=1 ialg=ialg1 iter=iter1 acc =accu1 call opt(parm,grad,hess,hinv,np, * xlf,fpnorm,grdir,ier, * mx,iter,acc,ialg) call smlout(parm,hinv,np, * xlf,ier,iter,ival,ialg,fpnorm,grdir,nsimul) call snooze(ipause) k0=0 call promeg(parm,np,k0) if (ispool.ne.0) call promeg(parm,np,kanal0) call snooze(ipause) c c ************************************************* c * method of simulated likelihood estimation * c ************************************************* c 300 if (ipause.lt.0.and.nsimul.eq.0) goto 400 write (*,301) 301 format( *//' ============================================================' * /' multinomial multiperiod probit model:' * /' smooth simulated likelihood estimation' * /' ============================================================') if (ipause.ge.0) then write (*,*) ' enter number of replications [0=quit] ?' cccvaherror!!! read (*,'(i1)') nsimul read (*,'(i4)') nsimul print *,nsimul call snooze(1) end if if (nsimul.eq.0) goto 400 if (ispool.ne.0) write (kanal0,301) c method=0 c c **** derivative check c ival=0 idiff2=2 if (ipause.ge.0) then write (*,*) ' derivative check of sml estimator', * ' [0=none, 1=func, 2=check] ?' read (*,'(i1)') mchk2 end if if (mchk2.ne.0) then if (ispool.ne.0) write (kanal0,'(//a)') * ' derivative check of sml estimator' call func(parm,np,fval,*998) write (*,'(a,f15.7)') ' fval=',fval if (ispool.ne.0) write (kanal0,'(a,f15.7)') ' fval=',fval if (mchk2.eq.2) then call deriv(parm,np,fval,grad,hess,idiff2,*998) write (*,'(a,30f10.4)') ' deriv=',(grad(i),i=1,np) if (ispool.ne.0) write (kanal0,'(a,30f10.4)') * ' deriv=',(grad(i),i=1,np) call fdiff(parm,np,fval,grad,*998) write (*,'(a,30f10.4)') ' fdiff=',(grad(i),i=1,np) if (ispool.ne.0) write (kanal0,'(a,30f10.4)') * ' fdiff=',(grad(i),i=1,np) end if call snooze(ipause) end if c c **** maximize likelihood c mx=1 ialg=ialg2 iter=iter2 acc =accu2 call opt(parm,grad,hess,hinv,np, * xlf,fpnorm,grdir,ier, * mx,iter,acc,ialg) call smlout(parm,hinv,np, * xlf,ier,iter,ival,ialg,fpnorm,grdir,nsimul) call snooze(ipause) k0=0 call promeg(parm,np,k0) if (ispool.ne.0) call promeg(parm,np,kanal0) c c ***************************** c * covariance correction * c ***************************** c 400 if (ipause.lt.0.and.nwesml.eq.0) goto 900 write (*,401) 401 format( *//' ============================================================' * /' multinomial multiperiod probit model:' * /' robust (wesml) covariance computation' * /' ============================================================') if (ipause.ge.0) then write (*,*) 'robust covariances ', * '[0=none, 1=old hess, 2=new hess] ?' read (*,'(i1)') nwesml end if if (nwesml.eq.0) goto 900 c if (ispool.ne.0) write (kanal0,401) call cwesml(parm,grad,hess,hinv,np,nwesml,*900) if (nwesml.eq.2) then ier1=100+ier call smlout(parm,hess,np, * xlf,ier1,iter,ival,ialg,fpnorm,grdir,nsimul) end if ier2=200+ier call smlout(parm,hinv,np, * xlf,ier2,iter,ival,ialg,fpnorm,grdir,nsimul) call snooze(ipause) c 900 stop '*** sml finished ***' 998 stop '*** sml aborted ***' end c c c subroutine func(parm,np,xlf,*) c----------------------------------------------------------------------- c accumulates criterion function over individuals c----------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' dimension parm(np) integer*4 nn1,iseed,iseed1 logical lderiv c common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns, * nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp common / dmode / method,model,modran,modar1,modmnp common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp) common / dopt1 / parm2(maxnp),grad2(maxnp),hess2(maxnp,maxnp), * hinv2(maxnp,maxnp) common / ddiff / idiff2 common / dival / ival common / doutp / kanal0,ispool,ipause common / dseed / iseed common / dsimu / nsimul c ccccccc write (*,'(//a,i3)') '$$$ func: np,parm=',np ccccccc write (*,'(9f8.5)') (parm(i),i=1,np) c if (method.eq.0.and.nsimul.eq.0) goto 98 c kobs=0 nerr=0 lderiv=.false. call omega1(parm,np,*97) c ---accumulate likelihood over nindiv nn1=1 iseed1=iseed xlf=0.d0 do 100 iobs=1,nindiv kobs=kobs+1 call pull(nn1) call cont(xlf2,grad2,lderiv,iseed1,*95) xlf=xlf+xlf2 100 continue c c error handling, exits c 90 ival=ival+1 return 95 write (*,*) ' func: failure in obs',kobs if (ispool.ne.0) write (kanal0,*) ' func: failure in obs',kobs return 1 97 write (*,*) ' func: failure in omega' if (ispool.ne.0) write (kanal0,*) ' func: failure in omega' return 1 98 write (*,*) ' func: misspecified nsimul=0' if (ispool.ne.0) write (kanal0,*) ' func: misspecified nsimul=0' 99 return 1 end c subroutine deriv(parm,np,xlf,grad,hess,idiff2,*) c----------------------------------------------------------------------- c c evaluates function, gradient (and hessian) of the c loglikelihood-function c c idiff2=4: bhhh approximation of hessian is cumulated c [hess = -grad*grad' = - inv(cov)] c otherwise: only function and gradient is computed c [hess will not be changed] c c----------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' dimension parm(np),grad(np),grad2(maxnp),hess(np,np) integer*4 nn1,iseed,iseed1 logical lderiv c common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns, * nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp common / dmode / method,model,modran,modar1,modmnp common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp) common / dival / ival common / doutp / kanal0,ispool,ipause common / dseed / iseed c ccc xlf=0.d0 do 5 i=1,np 5 grad(i)=0.d0 if (idiff2.eq.4) then do 6 i=1,np do 6 j=1,np 6 hess(i,j)=0.d0 end if kobs=0 lderiv=.true. c ---transfer parm on parmp and omega call omega2(parm,np,*95) c ---accumulate derivatives over nindiv nn1=1 iseed1=iseed do 100 iobs=1,nindiv kobs=kobs+1 call pull(nn1) call cont(xlf2,grad2,lderiv,iseed1,*95) ccc xlf=xlf+xlf2 do 110 k=1,np 110 grad(k) = grad(k) + grad2(k) if (idiff2.eq.4) then do 150 j=1,np do 150 k=1,np 150 hess(j,k) = hess(j,k) - grad2(j)*grad2(k) end if 100 continue c c error handling, exits c 90 ival=ival+2 return 95 write (*,*) ' failure in gradient evaluation, observation ',kobs if (ispool.ne.0) write (kanal0,*) * ' failure in gradient evaluation, observation ',kobs return 1 end c c subroutine cont(xlfi,grad,lderiv,iseed1,*) c----------------------------------------------------------------------- c c likelihood contribution xlfi for observation iobs c c ** version for (nalt-1)*nper stochastic utility differences ** c ** and first derivatives in grad c c--------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' real*8 grad(maxnp) integer*4 iseed1 logical lderiv c common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns, * nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp common / dmode / method,model,modran,modar1,modmnp common / domeg / omega(maxnt,maxnt),domega(maxnt,maxnt,0:maxns) common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp) common / dxvar / x(maxnt,maxnpp),ic(maxper),wt common / dzbar / zquer(maxnt1),covz(maxnt1,maxnt1), * dzquer(maxnt1,maxnb),dcovz(maxnt1,maxnt1,maxns) c c--------------------------------------------------- c i definition der zquer, c dzquer traegt ableitungen nach freien betas c--------------------------------------------------- c it =0 it1=0 do 10 ip=1,nper ccc write (*,*) '$$$ iper,nper=',ip,nper ialt=ic(ip) ita=it+ialt do 11 i=1,nalt it=it+1 if (i.eq.ialt) goto 11 it1=it1+1 ccc write (*,*) '$$$ i,ialt=',i,ialt sum=0 do 12 k=1,nxy ccc write (*,*) '$$$ it,it1,italt=',it,it1,ita xdiff=x(it,k)-x(ita,k) sum=sum+xdiff*parmp(k) i1=loc(k) if (i1.gt.0) dzquer(it1,i1)=xdiff 12 continue zquer(it1)=sum 11 continue 10 continue c ccc write (*,'(//a)') '$$$ zquer=' ccc write (*,'(9f8.5)') (zquer(i),i=1,nt1) ccc do 13 k=1,nb ccc write (*,'(//a)') '$$$ dzquer=' ccc13 write (*,'(9f8.5)') (dzquer(i,k),i=1,nt1) c c---------------------------------------------------- c ii berechnung der covarianzmatrix von z, c dcovz traegt ableitungen nach freien sigmas c---------------------------------------------------- c do 50 ip=1,nper ialt=ic(ip) ioff1=(ip-1)*nalt ioff2=(ip-1)*nalt1 iii=ioff1+ialt do 50 jp=1,nper jalt=ic(jp) joff1=(jp-1)*nalt joff2=(jp-1)*nalt1 jjj=joff1+jalt c omegij=omega(iii,jjj) c ki=1 do 20 i=1,nalt if (i.ne.ialt) then kj=1 do 25 j=1,nalt c if (j.ne.ialt) then covz(ioff2+ki,joff2+kj)= * omega(ioff1+i,joff1+j) * -omega(ioff1+i,jjj ) * -omega(iii ,joff1+j) * +omegij c if (lderiv) then do 26 is=1,ns dcovz(ioff2+ki,joff2+kj,is)= * domega(ioff1+i,joff1+j,is) * -domega(ioff1+i,jjj ,is) * -domega(iii ,joff1+j,is) * +domega(iii ,jjj ,is) 26 continue end if c kj=kj+1 end if c 25 continue ki=ki+1 end if 20 continue c 50 continue c ccc write (*,'(//a)') '$$$ omega=' ccc do 58 i=1,nt ccc58 write (*,'(9f8.5)') (omega(i,j),j=1,nt) ccc write (*,'(//a)') '$$$ covz=' ccc do 59 i=1,nt1 ccc59 write (*,'(9f8.5)') (covz(i,j),j=1,nt1) c c------------------------------------------------------ c iii probability evaluation c------------------------------------------------------- c 900 if (method.eq.1) then call logit2(zquer,dzquer,prob,grad,lderiv,*991) else if (method.eq.2) then call clark2(zquer,dzquer,covz,dcovz,prob,grad,lderiv,*991) else call gewek2(zquer,dzquer,covz,dcovz,prob,grad,lderiv,iseed1,*991) end if c if (prob.lt.1.d-300) prob=1.d-300 cccc write (*,*) '$$$ average prob=',prob c xlfi=wt*dlog(prob) if (lderiv) then do 910 k=1,nb+ns 910 grad(k)=wt*grad(k)/prob end if c return 991 return 1 end c c c c====================================================================== c memory resident utilities c====================================================================== c c c subroutine snooze(ipause) c---------------------------------------------------------- c gives you a little break to view output c c ipause<0 skips c ipause in seconds<60. requires system clock call. c ipause=60 pauses c---------------------------------------------------------- implicit integer (i-n) integer*4 isec0,isec1,isec2 integer*2 i0,i1,i2,i3,i4 if (ipause.le.0) goto 90 if (ipause.ge.60) goto 10 ccc write (*,'(a,i2,a\)') ' ...pausing for ',ipause,' seconds' write (*,'(a,i2,a)') ' ...pausing for ',ipause,' seconds' i0=0 call gettim(i1,i2,i0,i4) isec0=i1*3600+i2*60+i0 isec2=isec0 isec=0 i3=0 do 1000 i=1,32000 call gettim(i1,i2,i3,i4) isec1=i1*3600+i2*60+i3 isec=isec1-isec0 if (isec1.ne.isec2) then ccc write (*,'(a1\)') '.' write (*,'(a1)') '.' isec2=isec1 end if if (isec.gt.ipause) goto 90 1000 continue goto 90 10 write (*,'(/'' hit enter to continue.....''/)') read (*,'(a1)') 90 write (*,'(a1)') ' ' return end c subroutine tstamp(tdate) c---------------------------------------------------------- c generates string*21 tdate with time and date c e.g.: 31 dec 1989, 12:00:00 c---------------------------------------------------------- implicit integer (i-n) integer*2 iyear,imon,iday,ihour,imin,isec,ihsec character*21 tdate character*3 month(12) character*2 s1,s2,s3,s4,s5 data month/'jan','feb','mar','apr','may','jun', * 'jul','aug','sep','oct','nov','dec'/ call getdat(iyear,imon,iday) isec=0 call gettim(ihour,imin,isec,ihsec) call int2ch(iday,s1) iyear=iyear-1900 call int2ch(iyear,s2) call int2ch(ihour,s3) call int2ch(imin,s4) call int2ch(isec,s5) tdate=s1//' '//month(imon)//' 19'//s2//', '//s3//':'//s4//':'//s5 return end c subroutine int2ch(int,ch2) c---------------------------------------------------------- c converts 2-digit integer to 2-byte character c---------------------------------------------------------- integer*2 int,i1,i2 character*1 c1,c2 character*2 ch2 i1=int/10 i2=int-i1*10 if (i1.le.9.and.i2.le.9.and.i1.ge.0.and.i2.ge.0) then c1=char(i1+48) c2=char(i2+48) ch2=c1//c2 else ch2='??' end if return end c subroutine openf(kanal,text,defaul,mode) c--------------------------------------------------------------------- c conditional and interactive opening of files c mode: 1 = old, formatted 2 = old, unformatted c 3 = new, formatted 4 = new, unformatted c--------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) character prompt*1,defaul*40,text*40,infile*40 logical exist,fmt common / doutp / kanal0,ispool,ipause c 10 write (*,11) text,defaul 11 format(//' enter filename for ',a40 * /' [default name is ',a40,']' ccc * /' >'\) * /' >') if (ipause.ge.0) then read (*,'(a40)') infile if (infile.eq.' ') infile=defaul if (infile.eq.' ') goto 10 else infile=defaul end if write (*,12) infile,text if (ispool.ne.0.and.kanal.ne.kanal0) write (kanal0,12) infile,text 12 format(/' opening ',a40 * /' for ',a40) inquire (file=infile,exist=exist) c if (mode.eq.1.or.mode.eq.3) then fmt=.true. else fmt=.false. end if c if (mode.le.2) then if (exist) then if (fmt) then open (kanal,file=infile,status='old',form='formatted') else open (kanal,file=infile,status='old',form='unformatted') end if else write (*,*) '*** error, does not exist' if (ipause.lt.0) stop '*** hlogit aborted ***' goto 10 end if else if (exist) then write (*,'('' *** warning: overwrite ? (y/n) >'')') if (ipause.ge.0) then read (*,'(a1)') prompt if (prompt.eq.'n'.or.prompt.eq.'n') goto 10 end if if (fmt) then open (kanal,file=infile,status='old',form='formatted') else open (kanal,file=infile,status='old',form='unformatted') end if else if (fmt) then open (kanal,file=infile,status='new',form='formatted') else open (kanal,file=infile,status='new',form='unformatted') end if end if end if return end c subroutine matp(a,n,nadim,k,b,nbdim,l,c,ncdim) c---------------------------------------------------------------------- c c(n,l) = a(n,k) * b(k,l) c---------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) dimension a(nadim,k),b(nbdim,l),c(ncdim,l) c do 50 i=1,n do 50 j=1,l sum=0.d0 do 51 m=1,k 51 sum=sum+a(i,m)*b(m,j) 50 c(i,j)=sum return end c real*8 function dotp(x,np,y) c---------------------------------------------------------------------- c dotp=x'y c---------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) dimension x(np),y(np) dotp=0. do 10 i=1,np 10 dotp=dotp+x(i)*y(i) return end c subroutine move(x1,np,x2) c---------------------------------------------------------------------- c vector x2=x1 c---------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) dimension x1(np),x2(np) do 10 i=1,np 10 x2(i)=x1(i) return end c subroutine movem(a1,a2,np,maxnp) c---------------------------------------------------------------------- c matrix a2=a1 c---------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) dimension a1(maxnp,np),a2(maxnp,np) do 10 i=1,np do 10 j=1,np 10 a2(i,j)=a1(i,j) return end c subroutine step(parm,np,size,direc,pnew) c---------------------------------------------------------------------- c pnew=parm+size*direc c---------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) dimension parm(np),direc(np),pnew(np) do 10 i=1,np 10 pnew(i)=parm(i)+size*direc(i) return end c subroutine gaussj(a,b,n,npmax,*) c--------------------------------------------------------------------- c solves ax=b by gauss-jordan elimination. full pivoting. c on return: a carries inv(a), b carries x c--------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) parameter (nmax=50) dimension a(npmax,n),b(n),ipiv(nmax),indxr(nmax),indxc(nmax) do 11 j=1,n ipiv(j)=0 11 continue do 22 i=1,n big=0.0 do 13 j=1,n if (ipiv(j).ne.1) then do 12 k=1,n if (ipiv(k).eq.0) then if (abs(a(j,k)).ge.big) then big=abs(a(j,k)) irow=j icol=k endif else if (ipiv(k).gt.1) then return 1 endif 12 continue endif 13 continue ipiv(icol)=ipiv(icol)+1 if (irow.ne.icol) then do 14 l=1,n dum=a(irow,l) a(irow,l)=a(icol,l) a(icol,l)=dum 14 continue dum=b(irow) b(irow)=b(icol) b(icol)=dum endif indxr(i)=irow indxc(i)=icol if (a(icol,icol).eq.0.0) return 1 pivinv=1./a(icol,icol) a(icol,icol)=1.0 do 16 l=1,n a(icol,l)=a(icol,l)*pivinv 16 continue b(icol)=b(icol)*pivinv do 21 ll=1,n if (ll.ne.icol) then dum=a(ll,icol) a(ll,icol)=0.0 do 18 l=1,n a(ll,l)=a(ll,l)-a(icol,l)*dum 18 continue b(ll)=b(ll)-b(icol)*dum endif 21 continue 22 continue do 24 l=n,1,-1 if (indxr(l).ne.indxc(l)) then do 23 k=1,n dum=a(k,indxr(l)) a(k,indxr(l))=a(k,indxc(l)) a(k,indxc(l))=dum 23 continue endif 24 continue return end cc subroutine trunc2(urand,rr,eps,prob,deps,dprob) c----------------------------------------------------------------------- c eps: right truncated normal at rr for uniform random variate urand c prob: prob(-inf < standard normal variate < rr) c c deps,dprob: first derivatives of eps and prob with respect to rr c----------------------------------------------------------------------- implicit real*8 (a-h,o-z) prob=cumphi(rr) dprob=denphi(rr) p=prob*urand dp=dprob*urand call phinv2(p,eps,deps) deps=deps*dp if (eps.le.-8.25) then cover write (*,'(a,f12.5,d25.17)') cover* ' overflow in cuminv: rr,h(rr)=',rr,prob if (eps.gt.rr) eps=rr end if return end c subroutine trunc(urand,rr,eps,prob) c----------------------------------------------------------------------- c eps: right truncated normal at rr for uniform random variate urand c prob: prob(-inf < standard normal variate < rr) c----------------------------------------------------------------------- implicit real*8 (a-h,o-z) prob=cumphi(rr) p=prob*urand call phinv(p,eps) if (eps.le.-8.25) then cover write (*,'(a,f12.5,d25.17)') cover* ' overflow in cuminv: rr,h(rr)=',rr,prob if (eps.gt.rr) eps=rr end if return end c subroutine phinv2(prob,arg,darg) c--------------------------------------------------------------------- c c inverse normal c.d.f.: prob=cumphi(cuminv(prob)) c arg =cuminv(cumphi(arg )) c c (1) formula 26.2.23 of abramovitz and stegun p.933 10th edition. c (accuracy<4.5e-4) c (2) inverse interpolation to 3rd order (a&s example 5 p.954) c c note: "darg" carries first derivative of "arg" w.r.t. "prob". c c modified from vah 24sep84/14apr85. c axel boersch-supan (c) version 25. august 1990 c c--------------------------------------------------------------------- implicit real*8 (a-h,o-z) data c0,c1,c2/2.515517d0,.802853d0,.010328d0/, * d1,d2,d3/1.432788d0,.189269d0,.001308d0/ c c crude ways of dealing with prob=0 and prob=1 c if (prob.lt.1.d-37) then arg=-8.25d0 darg=2.d15 return end if p1=1.d0-prob if (p1.lt.1.d-37) then arg=8.25d0 darg=2.d15 return end if c c flip if prob>0.5 c p=prob iflip=0 if (prob.gt.0.5d0) then p=p1 iflip=1 end if c c p<0.5 as required by inversion: c t=dsqrt(dlog(1.d0/(p**2))) dt=-1.d0/(t*p) top =c0+c1*t+c2*t**2 dtop= c1 +2*c2*t bot =1+d1*t+d2*t**2+d3*t**3 dbot= d1 +2*d2*t +3*d3*t**2 x0=t-top/bot dx0=(1-dtop/bot+(dbot/bot)*(top/bot))*dt c c up to now x0<0: c x0=-x0 dx0=-dx0 c c now x0>0 as required by inverse interpolation: c denx=denphi(x0) cumx=cumphi(x0) ddenx=-denx*x0*dx0 dcumx= denx*dx0 t=(p-cumx)/denx dt=(1-dcumx)/denx-t*ddenx/denx c ph1 = x0 + t + 0.5* x0*t**2 dph1 =dx0 +dt + 0.5*dx0*t**2 + x0*t*dt ph2 = ((1+2*x0**2)*t**3)/6.d0 dph2 = ( 4*x0*dx0*t**3)/6.d0 +((1+2*x0**2)*t**2)*dt/2.d0 arg=ph1+ph2 darg=dph1+dph2 c c flip back if prob>0.5 c if (iflip.eq.1) arg=-arg return end c subroutine phinv(prob,arg) c--------------------------------------------------------------------- c c inverse normal c.d.f.: prob=cumphi(cuminv(prob)) c arg =cuminv(cumphi(arg )) c c (1) formula 26.2.23 of abramovitz and stegun p.933 10th edition. c (accuracy<4.5e-4) c (2) inverse interpolation to 3rd order (a&s example 5 p.954) c c modified from vah 24sep84/14apr85. c axel boersch-supan (c) version 25. august 1990 c c--------------------------------------------------------------------- implicit real*8 (a-h,o-z) data c0,c1,c2/2.515517d0,.802853d0,.010328d0/, * d1,d2,d3/1.432788d0,.189269d0,.001308d0/ c c crude ways of dealing with prob=0 and prob=1 c if (prob.lt.1.d-37) then arg=-8.25d0 darg=2.d15 return end if p1=1.d0-prob if (p1.lt.1.d-37) then arg=8.25d0 darg=2.d15 return end if c c flip if prob>0.5 c p=prob iflip=0 if (prob.gt.0.5d0) then p=p1 iflip=1 end if c c p<0.5 as required by inversion: c t=dsqrt(dlog(1.d0/(p**2))) top =c0+c1*t+c2*t**2 bot =1+d1*t+d2*t**2+d3*t**3 x0=t-top/bot c c up to now x0<0: c x0=-x0 c c now x0>0 as required by inverse interpolation: c denx=denphi(x0) cumx=cumphi(x0) t=(p-cumx)/denx c ph1 = x0 + t + 0.5* x0*t**2 ph2 = ((1+2*x0**2)*t**3)/6.d0 arg=ph1+ph2 c c flip back if prob>0.5 c if (iflip.eq.1) arg=-arg return end c c c c real*8 function denphi(x) implicit real*8 (a-h,o-z) c 1/sqrt(2*pi) data constp/0.39894228040143267793d00/ t=-0.5*x*x if (t.lt.-707.0) then cover write (*,*) 'denphi(x) underflow. x=',x denphi=0 else denphi=dexp(t)*constp end if return end c real*8 function cumphi(x) implicit real*8 (a-h,o-z) data sqrt2/1.414213562d0/ if (x.ge.0) then z=x/sqrt2 cumphi=1.0-0.5*erfc(z) else z=-x/sqrt2 cumphi=0.5*erfc(z) end if return end c real*8 function erfc(x) implicit real*8 (a-h,o-z) data a0/-1.26551223d0/, * a1/ 1.00002368d0/, * a2/ 0.37409196d0/, * a3/ 0.09678418d0/, * a4/-0.18628806d0/, * a5/ 0.27886807d0/, * a6/-1.13520398d0/, * a7/ 1.48851587d0/, * a8/-0.82215223d0/, * a9/ 0.17087277d0/ z=dabs(x) t=1.0/(1.0+0.5*z) p9=a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*a9)))))))) zzz=-z*z+p9 if (zzz.lt.-707.0) then cover write (*,*) 'erfc(x) underflow. x=',x zzz=0 else zzz=t*exp(zzz) end if if (x.lt.0) zzz=2.0-zzz erfc=zzz return end c real*8 function randu(iseed) c----------------------------------------------------------------------- c generates randu pseudo uniform [0,1] random number generator. c requires a seed (iseed) which is updated and returned. c [lewis,goodman & miller,ibm systems journal,v.8#2 1969,pp.136-145] c----------------------------------------------------------------------- implicit real*8 (a-z) integer*4 iseed data xb,xm/16807.d0,2147483647.d0/ seed=iseed xseed=dmod(xb*seed,xm) randu=xseed/xm iseed=xseed return end c c subroutine chol2(a,da,c,dc,n,maxn,ns,maxns,*) c---------------------------------------------------------------------- c cholesky decomposition of symmetric nxn matrix a (dimension maxn) c results in lower triangular matrix c with a=c*c' c c *** version with first partials: ** c dc carries first partials of c w.r.t. ns deep parameters in a, c if da carries first partials of a w.r.t. deep parameters. c c modified from stoer/bulirsch c (c) axel boersch-supan version: 25. august 1990 c---------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) dimension a(maxn,maxn),c(maxn,maxn), * da(maxn,maxn,maxns),dc(maxn,maxn,maxns), * dx(10) notpos=0 do 10 i=1,n do 10 j=i,n x=a(i,j) do 101 is=1,ns 101 dx(is)=da(i,j,is) do 20 kk=1,i-1 k=i-kk x=x-a(j,k)*a(i,k) do 201 is=1,ns 201 dx(is)=dx(is)-da(j,k,is)*a(i,k)-a(j,k)*da(i,k,is) 20 continue if (i.eq.j) then if (x.le.0) goto 99 sqrtx=dsqrt(x) c(i,i)=1.0/sqrtx do 211 is=1,ns 211 dc(i,i,is)=-0.5*dx(is)/(x*sqrtx) else a(j,i)=x*c(i,i) do 221 is=1,ns 221 da(j,i,is)=dx(is)*c(i,i)+x*dc(i,i,is) end if 10 continue 30 do 31 i=1,n do 32 j=i,n if (i.eq.j) goto 32 c(j,i)=a(j,i) c(i,j)=0.0 a(j,i)=a(i,j) do 301 is=1,ns dc(j,i,is)=da(j,i,is) dc(i,j,is)=0.0 301 da(j,i,is)=da(i,j,is) 32 continue dd=c(i,i) c(i,i)=1.0/dd dd=dd**2 do 321 is=1,ns 321 dc(i,i,is)=-dc(i,i,is)/dd 31 continue if (notpos.eq.1) return 1 goto 90 99 write (*,*) 'cholesky: matrix not positive definite' notpos=1 goto 30 90 return end c subroutine chol(a,c,n,maxn,*) c---------------------------------------------------------------------- c cholesky decomposition of symmetric nxn matrix a (dimension maxn) c results in lower triangular matrix c with a=c*c' c c modified from stoer/bulirsch c (c) axel boersch-supan version: 25. august 1990 c---------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) dimension a(maxn,maxn),c(maxn,maxn) notpos=0 do 10 i=1,n do 10 j=i,n x=a(i,j) do 20 kk=1,i-1 k=i-kk x=x-a(j,k)*a(i,k) 20 continue if (i.eq.j) then if (x.le.0) goto 99 sqrtx=dsqrt(x) c(i,i)=1.0/sqrtx else a(j,i)=x*c(i,i) end if 10 continue 30 do 31 i=1,n do 32 j=i,n if (i.eq.j) goto 32 c(j,i)=a(j,i) c(i,j)=0.0 a(j,i)=a(i,j) 32 continue dd=c(i,i) c(i,i)=1.0/dd 31 continue if (notpos.eq.1) return 1 goto 90 99 write (*,*) 'cholesky: matrix not positive definite' notpos=1 goto 30 90 return end c subroutine pull(nn1) c----------------------------------------------------------------------- c pulls data for one observation from internal worksheet. c (nda items beginning with offset nn1) c c this routine in combination with sr trans is designed to be the c only access to the internal worksheet. c the nesting (pull calls trans) is necessary for the data(*) call c in order to isolate routines that must be compiled with the $large c metacommand (under ms-fortran). c c----------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' real*4 data(maxbuf) integer*4 nn1 common / ddata / data common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns, * nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp common / dmode / method,model,modran,modar1,modmnp call trans(data(nn1)) nn1=nn1+ndp return end c ccc$large subroutine trans(data) c----------------------------------------------------------------------- c transfers relevant piece of data(*) on matrix x in common dtrans. c if maxbuf is larger than one memory segment, this (and only this) c routine must be compiled with the $large metacommand (ms-fortran) c----------------------------------------------------------------------- c data set-up: x-vars (nx), choice (1), weight(nwt), y-vars (ny) c ndd=nx*nalt c ndw=nx*nalt+1+nwt c nda=nx*nalt+1+nwt+ny = number of data items per observation c ndp=nda*nper = number of data items per individual c nb=nx+nalt1*(ny+nc) = number of parameters in beta/x c----------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' c-----$large data----------------------include, if maxbuf>16000--------- real*4 data(*) common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns, * nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp common / dxvar / x(maxnt,maxnpp),ic(maxper),wt c ccccc write (*,'(a,7f9.4,5(/8f9.4))') '$$$ data=',(data(i),i=1,ndp) c mper1=0 mper2=0 do 1000 iper=1,nper c c chosen alternative and weight c ic(iper)=data(mper1+ndd+1)+.5 ialt=ic(iper) if (ialt.lt.1.or.ialt.gt.nalt) then write (*,*) 'error in cont: ialt=',data(mper1+ndd+1) write (*,*) ' iper,ioff,ndd=',iper,mper1+ndd+1,ndd end if wt=1.0 if (nwt.eq.1) wt=data(mper1+ndd+2) c c transfer observation-relevant piece of data(*) c c ---alternative-specific variables k1=0 do 10 k=1,nx do 10 i=1,nalt k1=k1+1 10 x(mper2+i,k)=data(mper1+k1) c ---agent-specific variables k1=nx do 20 k=1,ny da=data(mper1+ndw+k) do 20 k2=1,nalt1 k1=k1+1 do 22 i=1,nalt 22 x(mper2+i,k1)=0.0 x(mper2+k2,k1)=da 20 continue c ---alternative-specific constants do 30 k2=1,nalt1 k1=k1+1 do 32 i=1,nalt 32 x(mper2+i,k1)=0.0 x(mper2+k2,k1)=1.0 30 continue c mper1=mper1+nda mper2=mper2+nalt 1000 continue c cc i=0 cc do 50 i1=1,nper cc write (*,*) 'nper,iper,ic(iper) / x=',nper,i1,ic(i1) cc do 50 i2=1,nalt cc i=i+1 cc50 write (*,'(a,7f9.4/8f9.4)') '$$$ x=',(x(i,k),k=1,nxy) c return end c c$include smlfunc.for c$include smlinp.for c$include smlutil.for c$include smlopt.for c$include wf.add c--------------------------------------------------------------------- c.pl100/.rm72/ c sml-func: likelihood function c ============================= c c copyright (c) axel boersch-supan, version 10. march 1990 c c *** updated from mppmss.for: 1. july 1990 c *** updated from mslsiml.for: 27. august 1990 c *** merged from sml-siml/init.for: 1. september 1990 c c subroutines: gewek2 c clark2 c logit2 c omega1 c omega2 c c--------------------------------------------------------------------- c subroutine gewek2(zquer,dzquer,covz,dcovz, * prob,dprob,lderiv,iseed1,*) c----------------------------------------------------------------------- c c geweke algorithm to generate unbiased nt-dim choice probability c c u = n(0,covz) s.t. .mdnm/-l s u s zquer c c c axel boersch-supan (c) 25. august 1990 c c *** version with analytic first derivatives dprob c w.r.t. deep parms in covz c c----------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' real*8 dprob(maxnp), * zquer(maxnt1),covz(maxnt1,maxnt1), * dzquer(maxnt1,maxnb),dcovz(maxnt1,maxnt1,maxns) integer*4 iseed1 logical lderiv c common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns, * nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp common / dsimu / nsimul common / dfunc / cmat(maxnt1,maxnt1),e(maxnt1), * dcmat(maxnt1,maxnt1,maxns),de(maxnt1,maxnp), * dgprob(maxnp),drr(maxnp) common / doutp / kanal0,ispool,ipause c c note: w-dimension = nt1 = (nalt-1)*nper c c initialize c prob=0 np=nb+ns mxt1=maxnt1 c c *** get cholesky factor of covz (on lower triangle) c if (lderiv) then do 10 ip=1,np 10 dprob(ip)=0 mxs=maxns call chol2(covz,dcovz,cmat,dcmat,nt1,mxt1,ns,mxs,*990) else call chol(covz,cmat,nt1,mxt1,*990) end if c c loop on nsimul c c *** draw truncated normal eps and associated probability for each c alternative and period, store on e and probs c do 100 isimul=1,nsimul c j=1 dd=cmat(1,1) rr=-zquer(1)/dd urand=randu(iseed1) cccc write (*,*) 'isimul,urand=',isimul,urand c if (lderiv) then call trunc2(urand,rr,eps,zprob,deps,dzprob) e(1)=eps gprob=zprob do 110 ib=1,nb 110 drr(ib)=-dzquer(1,ib)/dd do 111 is=1,ns ip=nb+is 111 drr(ip)=-rr*dcmat(1,1,is)/dd do 115 ip=1,np de(1,ip)=deps*drr(ip) 115 dgprob(ip)=dzprob*drr(ip) else call trunc(urand,rr,eps,zprob) e(1)=eps gprob=zprob end if c do 200 j=2,nt1 dd=cmat(j,j) rr=0 do 210 i=1,j-1 210 rr=rr+cmat(j,i)*e(i) rr=(-zquer(j)-rr)/dd urand=randu(iseed1) c if (lderiv) then call trunc2(urand,rr,eps,zprob,deps,dzprob) e(j)=eps oldgpr=gprob if (zprob.lt.1.d-137) zprob=0.d0 gprob=gprob*zprob do 220 ib=1,nb drr(ib)=0 do 221 i=1,j-1 221 drr(ib)=drr(ib)+cmat(j,i)*de(i,ib) drr(ib)=(-dzquer(j,ib)-drr(ib))/dd 220 continue do 222 is=1,ns ip=nb+is drr(ip)=0 do 223 i=1,j-1 223 drr(ip)=drr(ip)+dcmat(j,i,is)*e(i) * + cmat(j,i) *de(i,ip) drr(ip)=-drr(ip)/dd-rr*dcmat(j,j,is)/dd 222 continue do 225 ip=1,np de(j,ip)=deps*drr(ip) dgprob(ip)=dgprob(ip)*zprob+oldgpr*dzprob*drr(ip) 225 continue else call trunc(urand,rr,eps,zprob) e(j)=eps if (zprob.lt.1.d-137) zprob=0.d0 gprob=gprob*zprob end if c 200 continue ccccccc write (*,*) 'isimul,prob=',isimul,gprob c c accumulate probability c prob=prob+gprob if (lderiv) then do 120 ip=1,np 120 dprob(ip)=dprob(ip)+dgprob(ip) end if c c end loop on nsimul c 100 continue c c average across nsimul c 900 prob=prob/float(nsimul) if (lderiv) then do 901 ip=1,np 901 dprob(ip)=dprob(ip)/float(nsimul) end if c return 990 write (*,*) '*** error in geweke: covz not pos.definite ***' if (ispool.ne.0) write (kanal0,*) * '*** error in geweke: covz not pos.definite ***' return 1 end c c subroutine omega1(parm,np,*) c======================================================================= c c creates covariance matrix omega (abbreviated from omega2) c c (c) axel boersch-supan, 27. august 1990 c c====================================================================== c implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' c dimension parm(np) c common calfa(maxalt,maxalt),dcalfa(maxalt,maxalt,0:maxns), * cnunu(maxalt,maxalt),dcnunu(maxalt,maxalt,0:maxns), * omegax(maxnt,maxnt) c common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns, * nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp common / dmode / method,model,modran,modar1,modmnp common / domeg / omega(maxnt,maxnt),domega(maxnt,maxnt,0:maxns) common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp) common / doffs / ns1,ns2,ns3,ns4,ns5,maps(maxalt,maxalt), * locs(maxnpp) common / doutp / kanal0,ispool,ipause c c (a) transfer updated parm vector on full parmp vector c and initialize matrices with identity/zero matrix c do 100 i=1,np 100 parmp(loc1(i))=parm(i) c do 110 i=1,nt do 111 j=1,nt omega(i,j)=0.0 111 continue omega(i,i)=1.0 110 continue c c (b) correlation matrix of unobserved utility components c (cases 2,4) c 2000 if (model.ne.2.and.model.ne.4) goto 3000 c do 210 ia=1,nalt do 210 ja=1,nalt c c fetch standard deviation and correlation of unobserved c utility components and undo transformation c ns1i=ns1+ia ns1j=ns1+ja ns2ij=ns2+maps(ia,ja) stni=parmp(ns1i) if (dabs(stni).gt.707.0) goto 9000 stni=dexp(stni) stnj=parmp(ns1j) if (dabs(stnj).gt.707.0) goto 9000 stnj=dexp(stnj) if (ia.eq.ja) then cnij=1.0 else temp=parmp(ns2ij) if (dabs(temp).gt.707.0) goto 9000 temp=dexp(temp) cnij=(temp-1)/(temp+1) end if c c compute covariance of unobserved utilities c cnij=stni*stnj*cnij cnunu(ia,ja)=cnij c 210 continue c c (c) correlation matrix of random effects c (cases 3,4,7,8) c 3000 if (model.le.2.or.model.eq.5.or.model.eq.6) goto 4000 c do 310 ia=1,nalt do 310 ja=1,nalt c c fetch standard deviation and correlation of random effects c and undo transformation c ns3i=ns3+ia ns3j=ns3+ja ns4ij=ns4+maps(ia,ja) stai=dexp(parmp(ns3i)) staj=dexp(parmp(ns3j)) if (ia.eq.ja) then caij=1.0 else temp=dexp(parmp(ns4ij)) caij=(temp-1)/(temp+1) end if c c compute covariance of unobserved utilities c caij=stai*staj*caij calfa(ia,ja)=caij c 310 continue c c (d) autoregressive process with iid nu's c (cases 5,7) c 4000 if (model.ne.5.and.model.ne.7) goto 5000 c do 400 ia=1,nalt do 400 ja=1,nalt c c fetch rho and undo transformation c ns5i=ns5+ia ns5j=ns5+ja temp=dexp(parmp(ns5i)) rhoi=(temp-1)/(temp+1) temp=dexp(parmp(ns5j)) rhoj=(temp-1)/(temp+1) c c identity for nu's c cnij=0.0 if (ia.eq.ja) cnij=1.0 c c time dependent part c do 420 it=1,nper do 420 jt=1,nper is=it-jt ita=(it-1)*nalt+ia jta=(jt-1)*nalt+ja if (jt.lt.it.and.rhoi.ne.0) then rhoit=rhoi**is omega(ita,jta)=rhoit*cnij else if (it.lt.jt.and.rhoj.ne.0) then rhojt=rhoj**(-is) omega(ita,jta)=rhojt*cnij else if (it.eq.jt) then omega(ita,jta)=cnij end if 420 continue 400 continue c c (e) autoregressive process with correlated nu's c (cases 6,8) c 5000 if (model.ne.6.and.model.ne.8) goto 8000 c do 500 ia=1,nalt do 500 ja=1,nalt c c fetch std.dev's and correlation of unobserved util's c and undo transformation c ns1i=ns1+ia ns1j=ns1+ja ns2ij=ns2+maps(ia,ja) stni=dexp(parmp(ns1i)) stnj=dexp(parmp(ns1j)) if (ia.eq.ja) then cnij=1.0 else temp=dexp(parmp(ns2ij)) cnij=(temp-1)/(temp+1) end if c c fetch rho and and undo transformations c ns5i=ns5+ia ns5j=ns5+ja temp=dexp(parmp(ns5i)) rhoi=(temp-1)/(temp+1) temp=dexp(parmp(ns5j)) rhoj=(temp-1)/(temp+1) rho1ij=1.0-rhoi*rhoj c c stationarity assumption c snij=stni*stnj/rho1ij cnij=snij*cnij c c time dependent part c do 520 it=1,nper do 520 jt=1,nper is=it-jt ita=(it-1)*nalt+ia jta=(jt-1)*nalt+ja c c 3 cases: it>jt (upper right blocks) c it=jt (diagonal blocks) c itjt (upper right blocks) c it=jt (diagonal blocks) c it' read (*,'(i1)') iquery if (iquery.ne.0) then write (*,*) 'enter nper, mod_mnp, mod_ran, mod_ar1 >' read (*,*) nper,modmnp,modran,modar1 goto 20 end if end if c 28 nbuff=nindiv*ndp if (nbuff.gt.maxbuf) then write (*,29) nindiv,ndp,maxbuf if (ispool.ne.0) write (kanal0,29) nindiv,ndp,maxbuf 29 format(' *** error: out of memory ***' * /' nindiv=',i9,' at ndp=',i4,' exceeds maxbuf=',i9) if (ipause.lt.0) return 1 write (*,*) 'enter new nindiv' read (*,*) nindiv goto 28 end if c c====================================================================== c c *** read in initialvalues **************************************** c c parmp = primary parameter vector with length npp c parm = estimated parameter vector with length np c c pname = label of variable c c loc = 0 if fixed c 1 if estimated c c note: the specification of model overrides istat. c c npp = number of primary parameters (any istat) c np = number of estimated parameters (istat=0) c c c---------------------------------------------------------------------- c c number of.... c estimated all parametrs: contents: c --------- -------------------- -------- c nb nxy=nx+nalt1*(ny+nc) beta in xbeta c ns 3*nalt+2*naltud sigma in omega c np npp stacked(beta,sigma) c --------- -------------------- -------- c parm parmp c c c---------------------------------------------------------------------- c c loc first status (1=estim, 0=fixed), then c location of i-th full parameter in estimated parm, c if parameter is estimated, 0 if fixed c e.g.: c grad(loc(ip)) is derivative associated with c full parameter component parm(ip) c c loc1 location of i-th estimated parameter in full parmp, c e.g.: c pname(loc1(ip1)) is parameter name associated with c estimated parameter parm(ip1) c c---------------------------------------------------------------------- c c parmp: c offset end length variable c ------ --- ---------------- ---------------------------------- c 1 nx nx x-variables c nx+1 nxy (ny+nc)*nalt y-interactions incl. constants c ns1+1 ns2 nalt std.dev's of unobserved utils c ns2+1 ns3 naltud correlation of unobserved utils c ns3+1 ns4 nalt std.dev's of random effects c ns4+1 ns5 naltud correlation of random effects c ns5+1 npp nalt autocorrelations c c notes: naltud = nalt*(nalt-1)/2 c c====================================================================== c 300 m13=13 call skip(kanal,m13) write (*,301) if (ispool.ne.0) write (kanal0,301) 301 format(//' initialvalues of parameters:' * /' ============================'/) do 310 k=1,npp read (kanal,311,end=997) pname(k),loc(k),parmp(k) 311 format(a8,i2,f15.9) write (*,312) pname(k),loc(k),parmp(k) if (ispool.ne.0) * write (kanal0,312) pname(k),loc(k),parmp(k) 312 format(1x,a8,1x,i2,1x,f15.9) if (k.eq.15.or.k.eq.30) call snooze(ipause) 310 continue if (npp.ne.15.and.npp.ne.30) call snooze(ipause) close(kanal) c c impose model structure. obey the following restrictions: c c nalt=4: | var1 cov1 cov2 0 | - restrictions for i,j=nalt c | var2 cov3 0 | due to differencing c | 1 0 | - restriction for var3 due to c | 0 0 0 1 | identification of beta c np=0 nb=0 ns=0 c c --- x and y var's do 320 i=1,nxy if (loc(i).ne.0) then np=np+1 nb=nb+1 loc1(np)=i loc(i)=np end if 320 continue c c ---interalternative correlations if (modmnp.eq.1) then do 330 i=1,nalt if (loc(ns1+i).ne.0) then np=np+1 ns=ns+1 loc1(np)=ns1+i loc(ns1+i)=np end if 330 continue do 331 i=1,naltud if (loc(ns2+i).ne.0) then np=np+1 ns=ns+1 loc1(np)=ns2+i loc(ns2+i)=np end if 331 continue else do 335 i=1,nalt 335 loc(ns1+i)=0 do 336 i=1,naltud 336 loc(ns2+i)=0 end if c c ---random effects if (modran.eq.1) then do 340 i=1,nalt if (loc(ns3+i).ne.0) then np=np+1 ns=ns+1 loc1(np)=ns3+i loc(ns3+i)=np end if 340 continue do 341 i=1,naltud if (loc(ns4+i).ne.0) then np=np+1 ns=ns+1 loc1(np)=ns4+i loc(ns4+i)=np end if 341 continue else do 345 i=1,nalt 345 loc(ns3+i)=0 do 346 i=1,naltud 346 loc(ns4+i)=0 end if c if (modar1.eq.1) then do 350 i=1,nalt if (loc(ns5+i).ne.0) then np=np+1 ns=ns+1 loc1(np)=ns5+i loc(ns5+i)=np end if 350 continue else do 355 i=1,nalt 355 loc(ns5+i)=0 end if c write (*,360) npp,nxy,npp-nxy,np,nb,ns if (ispool.ne.0) write (kanal0,360) npp,nxy,npp-nxy,np,nb,ns 360 format(//' number of parameters:' * /' ==========================' * /' primary parameters: ',i3 * /' of those: beta ',i3 * /' of those: sigma ',i3 * /' estimated parameters:',i3 * /' of those: beta ',i3 * /' of those: sigma ',i3) call snooze(ipause) if (np.eq.0) goto 996 if (nb.gt.maxnb) goto 9008 if (ns.gt.maxns) goto 9009 c c maps: location of correlation-matrix elements in parmp c k=0 do 380 i=1,nalt maps(i,i)=0 do 385 j=i+1,nalt k=k+1 maps(i,j)=k maps(j,i)=k 385 continue 380 continue c write (*,*) 'maps=' c do 389 i=1,nalt c write (*,'(10i3)') (maps(i,j),j=1,nalt) 389 continue c call snooze(ipause) c c locs: storage of derivatives in domega c do 390 i=1,npp 390 locs(i)=max(0,loc(i)-nb) c c reparametrize and bound away from -1,0,1, respectively c eps= 0.0001d0 elim=0.9999d0 do 401 i=1,nalt c ---- std.dev nu ppp=parmp(ns1+i) if (ppp.lt.0) goto 996 if (ppp.lt.eps) ppp=eps parmp(ns1+i)=log(ppp) c ---- std.dev alfa ppp=parmp(ns3+i) if (ppp.lt.0) goto 996 if (ppp.lt.eps) ppp=eps parmp(ns3+i)=log(ppp) c ---- rho ppp=parmp(ns5+i) if (ppp.lt.-1.0.or.ppp.gt.1.0) goto 996 if (ppp.lt.-elim) ppp=-elim if (ppp.gt. elim) ppp= elim parmp(ns5+i)=log((1+ppp)/(1-ppp)) 401 continue do 402 i=1,naltud c ---- correl nu ppp=parmp(ns2+i) if (ppp.lt.-1.0.or.ppp.gt.1.0) goto 996 if (ppp.lt.-elim) ppp=-elim if (ppp.gt. elim) ppp= elim parmp(ns2+i)=log((1+ppp)/(1-ppp)) c ---- correl alfa ppp=parmp(ns4+i) if (ppp.lt.-1.0.or.ppp.gt.1.0) goto 996 if (ppp.lt.-elim) ppp=-elim if (ppp.gt. elim) ppp= elim parmp(ns4+i)=log((1+ppp)/(1-ppp)) 402 continue c c finally: echo and transfer on actual parameter vector c 500 write (*,501) if (ispool.ne.0) write (kanal0,501) 501 format(//' parameters after transformation:' * /' ================================'/) do 510 k=1,npp write (*,312) pname(k),loc(k),parmp(k) if (ispool.ne.0) * write (kanal0,312) pname(k),loc(k),parmp(k) if (k.eq.15.or.k.eq.30) call snooze(ipause) 510 continue if (npp.ne.15.and.npp.ne.30) call snooze(ipause) do 512 i=1,np parm(i)=parmp(loc1(i)) 512 continue c c exit c 900 close (kanal) return 996 write (*,*) '*** error: initial values not admissible ***' if (ispool.ne.0) write (kanal0,*) * '*** error: initial values not admissible ***' return 1 997 write (*,*) '*** error: not enough initial values provided ***' if (ispool.ne.0) write (kanal0,*) * '*** error: not enough initial values provided ***' return 1 9003 write (*,9013) nalt,maxalt if (ispool.ne.0) write (kanal0,9013) nalt,maxalt 9013 format(' *** error: out of memory ***' * /' nalt=',i9,' exceeds max=',i9) return 1 9005 write (*,9015) nper,maxper if (ispool.ne.0) write (kanal0,9015) nper,maxper 9015 format(' *** error: out of memory ***' * /' nper=',i9,' exceeds max=',i9) return 1 9006 write (*,9016) nt,maxnt if (ispool.ne.0) write (kanal0,9016) nt,maxnt 9016 format(' *** error: out of memory ***' * /' nt=nper*nalt=',i9,' exceeds max=',i9) return 1 9007 write (*,9017) npp,maxnp if (ispool.ne.0) write (kanal0,9017) npp,maxnp 9017 format(' *** error: out of memory ***' * /' npp=',i9,' exceeds max=',i9) return 1 9008 write (*,9018) nb,maxnb if (ispool.ne.0) write (kanal0,9018) nb,maxnb 9018 format(' *** error: out of memory ***' * /' nb=',i9,' exceeds max=',i9) return 1 9009 write (*,9019) ns,maxns if (ispool.ne.0) write (kanal0,9019) ns,maxns 9019 format(' *** error: out of memory ***' * /' ns=',i9,' exceeds max=',i9) return 1 end c subroutine skip(kanal,lines) implicit integer (i-n) do 10 i=1,lines 10 read (kanal,'(a1)') return end c subroutine indata(*) c----------------------------------------------------------------------- c c initial data input for mpp-msl c - reads formatted "indata" file c - checks data for outliers (>10.000?) c - calculates weighted and unweighted means c c----------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' c$include:'sml.dim' c....................................... c....................................... real*4 data(maxbuf) character*72 fmtdat c common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns, * nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp common / ddata / data common / dfreq / wtsum,share(maxnt),wshare(maxnt), * amean(maxndm),wmean(maxndm) common / doutp / kanal0,ispool,ipause c c *** dialog (note: file have been opended in inparm) *** c kanal7=7 read (kanal7,'(a72)') fmtdat c write (*,10) nindiv,nper if (ispool.ne.0) write (kanal0,10) nindiv,nper 10 format(//' now starting to read ',i6,' individuals at nper=',i3 * /' ===================================================') nread=0 if (ipause.ge.0) then write (*,11) 11 format(' enter new number of individuals to read [0=keep] >') read (*,'(bn,i4)') nread end if if (nread.gt.nindiv) then write (*,*) * 'note: truncated to maximum memory allocation' if (ispool.ne.0) write (kanal0,*) * 'note: truncated to maximum memory allocation' nread=0 call snooze(ipause) end if if (nread.le.0) nread=nindiv c c *** offsets and lengths c ndy=1+nwt+ny c ---just to remember: c ndd=nx*nalt c nda=nx*nalt+1+nwt+ny c ---storage for means: ndm=ndd+ny c c *** initializations for means and frequencies c do 20 i=1,ndm*nper amean(i)=0.0 20 wmean(i)=0.0 do 21 i=1,nalt*nper share(i)=0.0 21 wshare(i)=0.0 c c *** loop through data c nn=0 iwt=0 nobs=0 nindiv=0 do 100 ihhd=1,nread do 101 iper=1,nper c ---check eof before all periods read ipe1=iper if (nn.gt.maxbuf) goto 993 c ---read data according to whether alt-spec var's exist if (nx.gt.0) then read (kanal7,fmtdat,err=298,end=299) * (data(nn+(j-1)*nalt+1),j=1,nx), * (data(nn+ndd+j),j=1,ndy) do 110 i=2,nalt read (kanal7,fmtdat,err=991,end=991) * (data(nn+(j-1)*nalt+i),j=1,nx) 110 continue else read (kanal7,fmtdat,err=298,end=299) * (data(nn+ndd+j),j=1,ndy) end if c ---counter nobs=nobs+1 c ---dependent variable in choice set? ic=data(nn+ndd+1) if (ic.lt.1.or.ic.gt.nalt) then write (*,121) nobs,ic if (ispool.ne.0) write (kanal0,121) nobs,ic 121 format(' invalid chosen alternative in observation',i6,i3) goto 999 end if c ---extract relevant weight wt=1.0 if (nwt.eq.1) wt=data(nn+ndd+2) if (wt.ne.1.0) iwt=1 if (wt.lt.0) then write (*,122) nobs,wt if (ispool.ne.0) write (kanal0,122) nobs,wt 122 format(' invalid weight in observation',i6,g20.10) goto 999 end if c ---accumulate sample frequencies icp=(iper-1)*nalt+ic share(icp) = share(icp) + 1.0 wshare(icp) = wshare(icp) + wt c ---check alt-specific attributes for typos and accumulate means indn=0 indm=(iper-1)*ndm do 140 k=1,nx do 140 i=1,nalt indn=indn+1 indm=indm+1 da=data(nn+indn) if (abs(da).gt.10000.0) then write (*,142) i,k,nobs,da if (ispool.ne.0) write (kanal0,142) i,k,nobs,da 142 format(' data check: x(ialt,k,iobs)=',3i5,g20.10) end if amean(indm) = amean(indm) + da wmean(indm) = wmean(indm) + da*wt 140 continue c ---check agent-specific characteristics and accumulate means indn=ndd+1+nwt do 145 k=1,ny indn=indn+1 indm=indm+1 da=data(nn+indn) if (abs(da).gt.10000.0) then write (*,146) k,nobs,da if (ispool.ne.0) write (kanal0,146) k,nobs,da 146 format(' data check: y(k,iobs)=',2i5,g20.10) end if amean(indm) = amean(indm) + da wmean(indm) = wmean(indm) + da*wt 145 continue c ---new offset nn=nn+nda 101 continue nindiv=nindiv+1 100 continue goto 300 c c *** end of file c 298 write (*,*) 'warning: err condition ended reading data' write (*,*) ' check format for data file' 299 if (nobs.eq.0) then write (*,*) ' *** error: data-file empty ***' return 1 end if close(kanal7) if (ipe1.ne.1) goto 991 c c *** sample means and frequencies c 300 wtsum=0.d0 do 320 i=1,nalt*nper 320 wtsum=wtsum+wshare(i) do 321 i=1,ndm*nper amean(i)=amean(i)/nindiv 321 wmean(i)=wmean(i)/wtsum*nper do 322 i=1,nalt*nper share(i)=share(i)/nindiv 322 wshare(i)=wshare(i)/wtsum*nper c 400 write (*,404) nindiv,nobs,wtsum if (ispool.ne.0) write (kanal0,404) nindiv,nobs,wtsum 404 format(//' unweighted number of individuals: ',i7 * /' unweighted number of observations:',i7 * /' weighted number of observations:',f9.1 * /' =========================================') c write (*,409) 'unweighted sample frequencies by period:' write (*,411) (i,i=1,nalt) do 410 iper=1,nper 410 write (*,412) iper,(share((iper-1)*nalt+i),i=1,nalt) call snooze(ipause) write (*,409) 'unweighted sample averages by period:' write (*,421) (i,i=1,ndm) do 420 iper=1,nper 420 write (*,422) iper,(amean((iper-1)*ndm+i),i=1,ndm) call snooze(ipause) c 409 format(/' ',a/' -----------------------------------------') 411 format(' ialt=' ,10i7) 412 format(' t=',i2,':',10f7.3) 421 format(' ivar=' ,10i7: ,5(/6x,10i7:)) 422 format(' t=',i2,':',10f7.3:,5(/6x,10f7.3:)) c if (iwt.eq.1) then write (*,409) 'weighted sample frequencies by period:' write (*,411) (i,i=1,nalt) do 415 iper=1,nper 415 write (*,412) iper,(wshare((iper-1)*nalt+i),i=1,nalt) call snooze(ipause) write (*,409) 'weighted sample averages by period:' write (*,421) (i,i=1,ndm) do 425 iper=1,nper 425 write (*,422) iper,(wmean((iper-1)*ndm+i),i=1,ndm) call snooze(ipause) else write (*,409) 'note: all weights are one' call snooze(ipause) end if c if (ispool.eq.0) goto 900 k=kanal0 write (k,409) 'unweighted sample frequencies by period:' write (k,411) (i,i=1,nalt) do 414 iper=1,nper 414 write (k,412) iper,(share((iper-1)*nalt+i),i=1,nalt) write (k,409) 'unweighted sample averages by period:' write (k,421) (i,i=1,ndm) do 424 iper=1,nper 424 write (k,422) iper,(amean((iper-1)*ndm+i),i=1,ndm) c if (iwt.eq.1) then write (k,409) 'weighted sample frequencies by period:' write (k,411) (i,i=1,nalt) do 419 iper=1,nper 419 write (k,412) iper,(wshare((iper-1)*nalt+i),i=1,nalt) write (k,409) 'weighted sample averages by period:' write (k,421) (i,i=1,ndm) do 429 iper=1,nper 429 write (k,422) iper,(wmean((iper-1)*ndm+i),i=1,ndm) else write (k,409) 'note: all weights are one' end if c 900 return 991 write (*,992) 992 format(/' alignment error:'/' at least one observation', * ' does not have records for all alternatives.') 999 write (*,*) '*** error in input data set ***' return 1 993 write (*,994) 994 format(/' *** fatal error: size of data set exceeds memory ***') return 1 end c----------------------------------------------------------------------- c.pl100/.rm72/ c sml-util: external utilities for method of simulated likelihood c c copyright (c) axel boersch-supan c version 26. february 1990, update from mpputil 1. july 1990 c update from mslutil 27. august 1990 c c subroutines: dcomeg derivative check of omega c fdiff numerical first derivatives c sdiff numerical second derivatives c promeg print omega c prtmat print matrix c mslout result file c cwesml wesml-robust covariance computation c c---------------------------------------------------------------------- subroutine dcomeg(parm,np,eps) c---------------------------------------------------------------------- c c checks derivatives of omega with respect to deep sigma parameters c c---------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' c$include:'sml.dim' c....................................... c....................................... c dimension parm(np) character*40 spfile,pafile c c common / domeg1/ common * calfa(maxalt,maxalt),dcalfa(maxalt,maxalt,0:maxns), * cnunu(maxalt,maxalt),dcnunu(maxalt,maxalt,0:maxns), * omegax(maxnt,maxnt) c common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns, * nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp common / domeg / omega(maxnt,maxnt),domega(maxnt,maxnt,0:maxns) common / doutp / kanal0,ispool,ipause common / dfile / spfile,pafile c c analytic derivs c write (*,*) 'derivative check at eps=',eps call omega2(parm,np,*99) write (*,'(/a)') ' omega matrix=' if (ispool.ne.0) write (kanal0,'(/a)') ' omega matrix=' do 215 i=1,nt write (*,'(12f6.3))') (omega(i,j),j=1,nt) if (ispool.ne.0) write (kanal0,'(12f6.3))') (omega(i,j),j=1,nt) do 215 j=1,nt 215 omegax(i,j)=omega(i,j) if (ipause.ge.0) pause 'hit enter for more....' c c numerical derivs c do 220 k=1,ns parm(nb+k)=parm(nb+k)+eps call omega1(parm,np,*99) devmax=-1.0 do 221 i=1,nt do 221 j=1,nt omega(i,j)=(omega(i,j)-omegax(i,j))/eps devij=dabs(omega(i,j)-domega(i,j,k)) if (devij.gt.devmax) devmax=devij 221 continue write (*,222) k,devmax if (ispool.ne.0) write (kanal0,222) k,devmax 222 format(//' analytical/numerical derivs w.r.t. s',i2.2, * ' *** max abs error=',d10.3,' ***'/' ',74('=')) do 224 i=1,nt write (*,223) ' a:',(domega(i,j,k),j=1,nt) write (*,223) ' n:',(omega(i,j),j=1,nt) if (ispool.ne.0) then write (kanal0,223) ' a:',(domega(i,j,k),j=1,nt) write (kanal0,223) ' n:',(omega(i,j),j=1,nt) end if 223 format(a3,20f8.4) 224 continue parm(nb+k)=parm(nb+k)-eps if (ipause.ge.0) write (*,227) devmax 227 format(/' ==> max absolute error = ',d15.6) if (ipause.ge.0) pause 'hit enter for more....' 220 continue return 99 write (*,*) 'error in dcomeg: inadmissible omega matrix' if (ispool.ne.0) write (kanal0,*) * 'error in dcomeg: inadmissible omega matrix' return end c c subroutine fdiff(parm,np,xlf,grad,*) c-------------------------------------------------------------------- c numerical first derivatives c-------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) dimension parm(np),grad(np) c ccccc (dummy statement, do not update xlf!!) ccccc dummy=xlf c call func(parm,np,xlf1,*98) cccc write (*,'(a,15f10.6)') 'fdiff:',(parm(j),j=1,np),xlf1 do 100 i=1,np ai=parm(i) ei=dmax1(dabs(ai)*.1d-3,.1d-5) 12 parm(i)=ai+ei call func(parm,np,xlf2,*15) cccc write (*,'(a,15f10.6)') 'fdiff:',(parm(j),j=1,np),xlf2 grad(i)=(xlf2-xlf1)/ei goto 10 15 ei=ei/10. if (ei.lt..1d-10) goto 99 goto 12 10 parm(i)=ai 100 continue return 98 write (*,*) 'error in fdiff: first func call' return 1 99 write (*,*) 'error in fdiff: can''t find suffiently small step' return 1 end c subroutine sdiff(parm,np,grad,hess,*) c-------------------------------------------------------------------- c numerical second derivatives from analytic first derivatives c-------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' c$include:'sml.dim' c....................................... c....................................... dimension parm(np),grad(np),grad1(maxnp),hess(np,np) character*21 tdate character*20 runtag common / dtag / runtag common / doutp / kanal0,ispool,ipause c idiff2=0 call deriv(parm,np,xlf,grad1,hess,idiff2,*98) call tstamp(tdate) fpnorm=dsqrt(dotp(grad1,np,grad1)) write (*,101) 0,fpnorm,tdate,runtag if (ispool.ne.0) write (kanal0,101) 0,fpnorm,tdate,runtag 101 format(' sdiff(',i2,'): fpnorm=',d12.6,'; ',a21,'; ',a20) c do 100 i=1,np ai=parm(i) ei=dmax1(dabs(ai)*.1d-3,.1d-5) 12 parm(i)=ai+ei idiff2=0 call deriv(parm,np,xlf,grad,hess,idiff2,*15) call tstamp(tdate) fpnorm=dsqrt(dotp(grad,np,grad)) write (*,101) i,fpnorm,tdate,runtag if (ispool.ne.0) write (kanal0,101) i,fpnorm,tdate,runtag do 13 j=1,np 13 hess(i,j)=(grad(j)-grad1(j))/ei goto 10 15 ei=ei/10. if (ei.lt..1d-10) goto 99 goto 12 10 parm(i)=ai 100 continue c do 20 i=1,np do 20 j=i+1,np h1=hess(i,j) h2=hess(j,i) h0=0.5*h1+0.5*h2 hess(i,j)=h0 hess(j,i)=h0 20 continue c return 98 write (*,*) 'error in sdiff: first deriv call' return 1 99 write (*,*) 'error in sdiff: can''t find suffiently small step' return 1 end c c subroutine promeg(parm,np,kanal) c---------------------------------------------------------------------- c print omega c---------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' c$include:'sml.dim' c....................................... c....................................... c dimension parm(np) character*8 fmt logical lderiv c common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns, * nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp common / dmode / method,model,modran,modar1,modmnp common / domeg / omega(maxnt,maxnt),domega(maxnt,maxnt,0:maxns) common / doutp / kanal0,ispool,ipause c lderiv=.false. call omega1(parm,np,*99) fmt='(15f5.2)' if (nt.le.12) fmt='(12f6.3)' if (nt.ge.16) fmt='(20f4.1)' if (kanal.eq.0) then write (*,'(//a)') ' omega matrix=' do 15 i=1,nt 15 write (*,fmt) (omega(i,j),j=1,nt) call snooze(ipause) else write (kanal,'(//a)') ' omega matrix=' do 25 i=1,nt 25 write (kanal,fmt) (omega(i,j),j=1,nt) end if return 99 write (*,*) 'error in promeg: inadmissible omega matrix' if (ispool.ne.0) write (kanal0,*) * 'error in promeg: inadmissible omega matrix' return end c subroutine prtmat(text,amat,n,m,ndim) c-------------------------------------------------------------------- c print matrix c-------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) character*(*) text dimension amat(ndim,m) common / doutp / kanal0,ispool,ipause c write (*,1) text,n,m,ndim if (ispool.gt.0) write (kanal0,1) text,n,m,ndim 1 format(' matrix ',a6,' [',i2.2,',',i2.2,'] [',i2.2,',*]') do 10 i=1,n write (*,'(9f8.3)') (amat(i,j),j=1,m) if (ispool.gt.0) write (kanal0,'(9f8.3)') (amat(i,j),j=1,m) 10 continue return end c c subroutine smlout(parm,hinv,np, * xlf,ier,iter,ival,ialg,fpnorm,grdir,nsimul) c-------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' c$include:'sml.dim' c....................................... c....................................... dimension parm(np),hinv(np,np) character*21 tdate character*20 text(3) character*52 text2(8) character*4 text3(4) character*8 pname(maxnpp) character*20 runtag character*40 spfile,pafile c common pvar(maxnpp),lmap(maxalt) c common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns, * nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp common / dmode / method,model,modran,modar1,modmnp common / doffs / ns1,ns2,ns3,ns4,ns5,maps(maxalt,maxalt), * locs(maxnpp) common / dname / pname common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp) common / doutp / kanal0,ispool,ipause common / dtag / runtag common / dfile / spfile,pafile c data text/'multinomial logit ', * 'clark approximation ', * 'simulated likelihood'/ c data text2/ * 'iid across periods, iid across alternatives ', * 'iid across periods, correlated across alternatives ', * 'random effects, iid across alternatives ', * 'random effects, correlated across alternatives ', * 'ar1 errors, iid across alternatives ', * 'ar1 errors, correlated across alternatives ', * 'random effects + ar1 errors, iid across alternatives', * 'random effects + ar1 errors, correlated across alts.'/ c data text3/' ???',' dfp','bfgs','bhhh'/ c kanal8=8 open(kanal8,file=pafile,status='old',access='append') if (ispool.ne.0) * open(kanal0,file=spfile,status='old',access='append') itext=3 if (method.eq.1) itext=1 if (method.eq.2) itext=2 write (*,84) model,text2(model),text(itext),text3(ialg) if (ispool.ne.0) write (kanal0,84) * model,text2(model),text(itext),text3(ialg) 84 format(//' ====================================================' * /' sml 2.1 multinomial-multiperiod probit model [' * ,i1,']' * /' ====================================================' * /' ',a52 * /' ====================================================' * /' estimation by: ',a20,' [',a4,']' * /' ====================================================' * //' variable estimate transformed std-err t-stat' * /' -------- -------- ----------- ------- ------') c c *** transfer deep parameters to parmp and variances to pvar *** c undo transformation for sigmas c npp=ns5+nalt j=0 do 171 i=1,npp pvar(i)=0.0 if (loc(i).ne.0) then j=j+1 parmp(i)=parm(j) pvar(i)=-hinv(j,j) end if 171 continue c c *** output of parameters *** c 180 do 181 i=1,npp m=loc(i) std = 0.0 tst = 0.0 if (m.gt.0) then if (pvar(i).gt.0.) then std = dsqrt(pvar(i)) tst = parmp(i)/std else std = 99.99 tst = 99.99 end if end if c c *** undo transformation *** c if (i.le.nxy) then c ---beta: [divide logit estimates by sqrt(pi/6)] pt=parmp(i) if (method.eq.1) pt=0.7797*pt else if ((i.le.ns2).or.(i.gt.ns3.and.i.le.ns4)) then c ---std.dev's nu and alfa: pt=exp(parmp(i)) else c ---rho, correl nu and alfa: pt=(exp(parmp(i))-1)/(exp(parmp(i))+1) end if c c *** write *** c 185 if (m.gt.0) write (*,186) pname(i),pt,parmp(i),std,tst if (m.eq.0) write (*,186) pname(i),pt,parmp(i) 186 format(1x,a8,f10.4,f14.4,f9.3,f8.2) if (ispool.ne.0) then if (m.gt.0) write (kanal0,186) pname(i),pt,parmp(i),std,tst if (m.eq.0) write (kanal0,186) pname(i),pt,parmp(i) end if j=loc(i) if (m.gt.0) write (kanal8,187) pname(i),j,pt,parmp(i),std,tst if (m.le.0) write (kanal8,187) pname(i),j,pt,parmp(i) 187 format(a8,i2,4f15.9) 181 continue c c *** likelihood at zero, mean square gradient c call tstamp(tdate) write (* ,192) model,modran,modar1,modmnp, * tdate,runtag,text(itext),text3(ialg), * xlf,nindiv,nper,nobs,nsimul, * iter,ival,fpnorm,grdir,ier if (ispool.ne.0) write (kanal0,192) * model,modran,modar1,modmnp, * tdate,runtag,text(itext),text3(ialg), * xlf,nindiv,nper,nobs,nsimul, * iter,ival,fpnorm,grdir,ier write (kanal8,192) model,modran,modar1,modmnp, * tdate,runtag,text(itext),text3(ialg), * xlf,nindiv,nper,nobs,nsimul, * iter,ival,fpnorm,grdir,ier 192 format(/' ===========================================' * /' sml ver.2.1 multinomial-multiperiod probit' * /' ===========================================' * /' model=',i1,15x,'(ran=',i1,' ar1=',i1,' mnp=',i1,')' * /' ===========================================' * /' ',a21,', ',a20, * /' estimation by: ',a20,' [',a4,']' * /' ===========================================' * /' function value ',f15.4 * /' number of individuals: ',i15 * /' number of periods: ',i15 * /' number of observations: ',i15 * /' number of replications: ',i15 * /' -------------------------------------------' * /' total number of iterations: ',i15 * /' total number of evaluations:',i15 * /' norm of gradient: ',d15.4 * /' gradient times direction: ',d15.4 * /' return code: ',i15 * /' ===========================================') close(kanal8) return end c c subroutine cwesml(parm,grad,hess,hinv,np,nwesml,*) c------------------------------------------------------------------- c c calculates correct covariance for wesml estimator c c covm = inv(hess) * ( bhhh + aba' ) * inv(hess) c c cwesml returns the corrected covariance matrix in hinv c c c nwesml=1: inv(hess) given in hinv c nwesml=2: recompute inv(hess) as finite differences of grad c c------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' c$include:'sml.dim' c....................................... c....................................... dimension parm(np),grad(np),hess(np,np),hinv(np,np) c common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns, * nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp common / doutp / kanal0,ispool,ipause common / ddiff / idiff2 c common scra(maxnp,maxnp) c if (nwesml.eq.1) then c ---save inverse hessian in "hess" do 10 i=1,np do 10 j=1,np 10 hess(i,j)=hinv(i,j) else if (nwesml.eq.2) then c ---compute inverse hessian in "hess" idiff2=0 call sdiff(parm,np,grad,hess,*98) ccc call prtmat('sdiff=',hess,np,np,np) call gaussj(hess,grad,np,np,*99) ccc call prtmat('inv(sdiff)=',hess,np,np,np) else ccc write (*,*) 'warning: nwesml=',nwesml return 1 end if c c calculate bhhh matrix in "hinv" c idiff2=4 call deriv(parm,np,xlf,grad,hinv,idiff2,*98) ccc call prtmat('bhhh=',hinv,np,np,np) c c pre- and post-multiply inv.hessian (hess) with bhhh (hinv) c store result in hinv c ndim=maxnp call matp(hess,np,np ,np,hinv,np,np,scra,ndim) call matp(scra,np,ndim,np,hess,np,np,hinv,np ) ccc call prtmat('wesml=',hinv,np,np,np) c return 98 write (*,*) 'error in cwesml: invalid initial values' if (ispool.ne.0) write (kanal0,*) * 'error in cwesml: invalid initial values' return 1 99 write (*,*) 'error in cwesml: hessian singular' if (ispool.ne.0) write (kanal0,*) * 'error in cwesml: hessian singular' return 1 end c---------------------------------------------------------------------- c smlopt: usual opt plus sml-specific itlog / 1. september 1990 c---------------------------------------------------------------------- subroutine opt(parm,grad,hess,hinv,np, * fval,fpnorm,grdir,ier, * mx,iter,acc,ialg) c---------------------------------------------------------------------- c.pl100/.rm72/ c driver and outer loop for nonlinear optimization routines c ========================================================= c axel boersch-supan [version feb 16, 1990] c c input: parm,np .......initial parms and dimension c mx ............1=max, 2=min c iter ..........iteration limit c acc ...........accuracy for convergence criterion c ialg ...* only 2=davidon-fletcher-powell algorithm c **sml** * and 4=newton-raphson with bhhh approximation c c output: parm,fval............ optimal parms and opt. func-value c grad ................ gradient c hess,hinv ........... hessian, inv.hessian (=-covmat) c fpnorm,grdir ........ norm of grad and grad*direc at opt c ier ................. return code c -9=hessian singular c -8=eigenvalues did not converge c -7=numerical saddlepoint c -6=cannot find improving step c -5=too many function errors c -4=*** out of memory *** c -3=function error in gradient evaluation c -2=initial values not admissible c -1=iteration limit exceeded c 0=*** unknown error *** c 1=stepsize convergence achieved c 2=gradient convergence achieved c 3=function value convergence achieved c 4=gradient*direction convergence achieved c c---------------------------------------------------------------------- c c sml-specific: c$include c common / doutp / c common / dopt9 / c c---------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' c$include:'sml.dim' c....................................... c....................................... c real*8 parm(np),grad(np),hess(np,np),hinv(np,np) character*1 prompt character*3 cmx(2) character*9 txt(6) character*21 tdate character*40 spfile,pafile character*45 text(14) c common / dopt9 / a1(maxnp),a2(maxnp),del(maxnp),fpd1(maxnp), * spd1(maxnp,maxnp),p(maxnp,maxnp),delr(maxnp) common / doutp / kanal0,ispool,ipause c common / dival / ival common / ddiff / idiff2 common / dstop / istop common / dfile / spfile,pafile c data cmx/'max','min'/ data text/ * 'hessian singular', * 'eigenvalues did not converge', * 'numerical saddlepoint', * 'cannot find improving step', * 'too many function errors', * '*** out of memory ***', * 'function error in gradient evaluation', * 'initial values not admissible', * 'iteration limit exceeded', * '*** unknown error ***', * 'stepsize convergence achieved', * 'gradient convergence achieved', * 'function value convergence achieved', * 'gradient*direction convergence achieved'/ data txt/'steepdesc','davflepow','newt-raph','newt-bhhh', * 'quadh-1up','quadh-2up'/ c c *** initializations *** c iterl=iter ispd=ialg ialg0=99 istop=1 ival=0 iteral=0 nevalf=0 do 10 i=1,np 10 hinv(i,i)=0.0 c c *** outer iteration *** c write (*,901) cmx(mx),np if (ispool.ne.0) write (kanal0,901) cmx(mx),np 901 format(/ */' =============================================================' */' general nonlinear optimization routine [version feb 16,1989]' */' ***** ',a3,'imization over',i3,' parameters *****' */' =============================================================') call func(parm,np,fval,*999) write (*,902) fval if (ispool.ne.0) write (kanal0,902) fval 902 format(' initial function value =',f20.8) 900 write (*,903) txt(ialg),iterl,acc,ialg,istop if (ispool.ne.0) write (kanal0,903) * txt(ialg),iterl,acc,ialg,istop 903 format( * ' -------------------------------------------------------------' */' default iteration parameters:',21x,'[',a9,']' */' iterations=',i3,', accuracy=',f8.6,', algorithm=',i1, * ', stopping=',i1 */' -------------------------------------------------------------' */' algorithms:' */' 1 = steepdesc = method of steepest descent' */' 2 = davflepow = davidon-fletcher-powell algorithm' */' 4 = newt-bhhh = newton-raphson with bhhh approximation' */' -------------------------------------------------------------') if (ipause.lt.0) goto 905 write (*,*) * ' hit return to accept default parameters [q=quit, c=change] >' read (*,'(a1)') prompt if (prompt.eq.'q'.or.prompt.eq.'q') goto 980 if (prompt.eq.'a'.or.prompt.eq.'a') goto 905 if (prompt.eq.' '.or.prompt.eq.'0') goto 905 write (*,904) 904 format(' enter new values for' * /' (1) iteration limit,' * /' (2) accuracy,' * /' (3) algorithm [2=dfp, 4=bhhh]:' * /' (4) stopping [1=all crit, 2=grad only]:') read (*,*) iterl,acc,ialg,istop 905 call tstamp(tdate) write (*,906) txt(ialg),iterl,acc,ialg,istop,tdate if (ispool.ne.0) write (kanal0,906) * txt(ialg),iterl,acc,ialg,istop,tdate 906 format( * ' -------------------------------------------------------------' */' actual iteration parameters:',22x,'[',a9,']' */' iterations=',i3,', accuracy=',f8.6,', algorithm=',i1, * ', stopping=',i1 */' -------------------------------------------------------------' */' time = ',a21) iterc=0 ival=0 ier=0 if (ispool.ne.0) then close(kanal0) open(kanal0,file=spfile,status='old',access='append') end if if (iterl.eq.0) goto 980 if (ialg.le.2) then idiff2=2 c ---do not initialize dfp-hessian if restart of same algorithm if (ialg.eq.ialg0) ier=99 call dfpx(parm,grad,hess,hinv,np, * fval,fpnorm,grdir, * mx,iterl,acc,ialg,iteral,ier, * a1,a2,del,fpd1,delr) else idiff2=4 call newton(parm,grad,hess,hinv,np, * fval,fpnorm,grdir, * mx,iterl,acc,ialg,iteral,ier, * a1,a2,del,fpd1,delr) end if if (ier.gt.4.or.ier.lt.-9) ier=0 910 if (ier.gt.0) then write (*,911) text(ier+10) if (ispool.ne.0) write (kanal0,911) text(ier+10) 911 format(//' optimum reached: ',a45//) else write (*,912) text(ier+10) if (ispool.ne.0) write (kanal0,912) text(ier+10) 912 format(//' failure to reach optimum: ',a45//) end if c c *** end of outer iteration *** c 950 iteral=iteral+iterc nevalf=nevalf+ival ialg0=ialg call itlog(parm,np,fval,fpnorm,grdir,iteral,'stop') write (*,951) ier 951 format(//' *** optimization terminated *** ier=',i3) if (ipause.lt.0) goto 980 write (*,*) 'hit return to quit, no to continue >' read (*,'(a1)') prompt if (prompt.eq.' '.or.prompt.eq.'0') goto 980 goto 900 999 ier=-2 goto 910 980 ival=nevalf iter=iteral return end c subroutine dfpx(parm,grad,hess,hinv,np, * fval,fpnorm,grdir, * mx,iterl,acc,ialg,iterc,ier, * dir,grad1,parm1,gam,hgam) c---------------------------------------------------------------------- c c nonlinear optimization routines 1 and 2 version feb 1990 c =========================================================== c c ialg: (1) method of steepest descent c (2) davidon-fletcher-powell algorithm c mx: maximization (1) or minimization (2) c acc: accuracy for convergence criterion c iterl iteration limit c c---------------------------------------------------------------------- implicit real*8(a-h,o-z) implicit integer (i-n) real*8 parm(np),grad(np),hess(np,np),hinv(np,np), * dir(np),parm1(np),grad1(np),gam(np),hgam(np) character*4 text c common / doutp / kanal0,ispool,ipause common / dstop / istop c data step1,stpmin,stracc,rsmall,absmal * /1.d0,1.d-10,1.0d-5,1.d-14,1.d-35/ c c start: initialize, initial function value c (unless restart with dfp) c if (ier.eq.99) then ier=0 text='more' goto 2000 end if c pm1=float(3-2*mx) ier=0 text='strt' st0=step1 grdir=0.0 delf=2*fval f1=-pm1*1.d20 f2=-pm1*1.d20 f3=-pm1*1.d20 iterl=iterc+iterl if (iterc.gt.0) goto 2000 c c initial grad and hess evaluation c ialg1=1 call deriv(parm,np,fval,grad,hess,ialg1,*9018) c c initial h-matrix for dfp and steepest descent c do 20 i=1,np do 21 j=1,np 21 hinv(i,j)=0. hinv(i,i)=-pm1 20 continue c c main iteration loop c 2000 if (iterc.ge.iterl) goto 9015 gradn=dotp(grad,np,grad) fpnorm=dsqrt(gradn) call itlog(parm,np,fval,fpnorm,grdir,iterc,text) iterc=iterc+1 c c check convergence by norm of gradient c if (fpnorm.le.acc) goto 9002 c c check convergence by function value (stationarity in three steps) c fmove=dabs(fval-f1)+dabs(f1-f2)+dabs(f2-f3) if (fmove.lt.fval*acc.and.istop.ne.2) goto 9003 f3=f2 f2=f1 f1=fval c c conjugate gradient direction c if (ialg.eq.1) goto 110 do 105 i=1,np sum=0.d0 do 106 j=1,np 106 sum=sum+hinv(i,j)*grad(j) 105 dir(i)=-sum c c watch for singularity c 109 grdir=dotp(dir,np,grad) if (grdir*pm1.gt.stracc**2) goto 120 c c steepest descent direction c 110 grdir=gradn*pm1 text='grad' do 115 i=1,np 115 dir(i)=pm1*grad(i) c c save old parameter and function values c 120 fval1=fval grdir1=grdir call move(parm,np,parm1) call move(grad,np,grad1) c c linesearch c z=1.0 if (iterc.le.np) z=dabs(st0) st0=dmax1(stpmin,dmin1(z,dabs(2.d0*delf/grdir))) call linesr(parm,np,fval,dir,pm1,st0,acc,gam,*9016) c c check convergence by parameter movement c 130 ier=1 do 131 i=1,np if (dabs(parm(i)-parm1(i)).gt.acc*dabs(parm(i))+rsmall) ier=0 131 continue if (istop.eq.2) ier=0 c c accept step c 140 delf=2.d0*(fval-fval1) 145 call deriv(parm,np,fval,grad,hess,ialg,*9018) grdir=dotp(dir,np,grad) if (ialg.eq.1) goto 2000 c c davidon-fletcher-powell c 170 z=-1.d0 call step(grad,np,z,grad1,gam) dg=st0*(grdir-grdir1) if (dg*pm1.ge.-absmal) goto 190 c c --update h-matrix c 180 do 181 i=1,np sum=0.d0 do 182 j=1,np 182 sum=sum+hinv(i,j)*gam(j) 181 hgam(i)=sum ghg=dotp(gam,np,hgam) if (ghg*pm1.ge.-absmal) ghg=.01*dg ca= 1.d0/dg cb=-1.d0/ghg c4= 1.d0/(grdir1-dg/st0**2) c3=-dg*c4/st0 call step(gam,np,c3,grad1,parm1) c3= st0/(c3*dg) do 183 i=1,np 183 dir(i)=st0*dir(i) if (dg*pm1.gt.ghg*pm1) goto 185 c c --broyden variant of dfp formula c text='dfpb' ca=1.0/(dg+ghg) cb=-ca call move(gam,np,parm1) c3=1.d0/dg c4=1.d0/grdir1 z=1.0+ghg/dg do 184 i=1,np 184 dir(i)=dir(i)*z-hgam(i) c c --original dfp formula c 185 text='dfph' do 186 i=1,np diri=dir(i)*ca hgami=cb*hgam(i) do 186 j=i,np hinv(i,j)=hinv(i,j)+diri*dir(j)+hgami*hgam(j) hinv(j,i)=hinv(i,j) 186 continue if (pm1*grdir.ge.0.5*grdir1*pm1) st0=2.0*st0 if (ier.eq.1) goto 9001 goto 2000 c c --use old h-matrix c 190 text='dfpo' grdir1=grdir st0=4.d0*st0 if (ier.eq.1) goto 9001 goto 2000 c c successful exits c 9001 ier=1 goto 9000 9002 ier=2 goto 9000 9003 ier=3 goto 9000 c c failure exits c 9015 ier=-1 goto 9000 9016 ier=-6 goto 9000 9018 ier=-3 goto 9000 9000 return end c subroutine newton(parm,grad,hess,hinv,np, * fval,fpnorm,grdir, * mx,iterl,acc,ialg,iterc,ier, * dir,grad1,parm1,gam,hgam) c---------------------------------------------------------------------- c c nonlinear optimization routines 3 and 4 version feb 1990 c =========================================================== c c ialg: (3) newton-raphson c (4) bhhh c mx: maximization (1) or minimization (2) c acc: accuracy for convergence criterion c iterl iteration limit c c---------------------------------------------------------------------- implicit real*8(a-h,o-z) implicit integer (i-n) real*8 parm(np),grad(np),hess(np,np),hinv(np,np), * dir(np),parm1(np),grad1(np),gam(np),hgam(np) character*4 text,meth(4) c common / doutp / kanal0,ispool,ipause common / dstop / istop c data meth/'grad','dfpx','newt','bhhh'/ data step1,stpmin,rsmall * /1.d0,1.d-10,1.d-14/ c c start: initialize, initial function value c hgam(1)=1.d0 pm1=float(3-2*mx) ier=0 text='strt' st0=step1 grdir=0.0 delf=2*fval f1=-pm1*1.d20 f2=-pm1*1.d20 f3=-pm1*1.d20 iterl=iterc+iterl if (iterc.gt.0) goto 2000 c c main iteration loop c 2000 if (iterc.ge.iterl) goto 9015 c c grad and hess evaluation c text=meth(ialg) call deriv(parm,np,fval,grad,hess,ialg,*9018) c c newton-raphson direction: solve dir = hess(-1)*grad c do 151 i=1,np dir(i)=grad(i) do 151 j=1,np 151 hinv(i,j)=hess(i,j) call gaussj(hinv,dir,np,np,*9019) grdir=dotp(dir,np,grad) c c iteration log c gradn=dotp(grad,np,grad) fpnorm=dsqrt(gradn) call itlog(parm,np,fval,fpnorm,grdir,iterc,text) iterc=iterc+1 c c check convergence by norm of gradient and by gradient*direction c if (fpnorm.le.acc.and.istop.ne.2) goto 9002 if (dabs(grdir).le.acc.and.istop.ne.2) goto 9004 c c check convergence by function value (stationarity in three steps) c fmove=dabs(fval-f1)+dabs(f1-f2)+dabs(f2-f3) if (fmove.lt.fval*acc.and.istop.ne.2) goto 9003 f3=f2 f2=f1 f1=fval c c save old parameter and function values c 120 fval1=fval grdir1=grdir call move(parm,np,parm1) call move(grad,np,grad1) c c linesearch c z=1.0 if (iterc.le.np) z=dabs(st0) st0=dmax1(stpmin,dmin1(z,dabs(2.d0*delf/grdir))) call linesr(parm,np,fval,dir,pm1,st0,acc,gam,*9016) c c check convergence by parameter movement c 130 ier=1 do 131 i=1,np if (dabs(parm(i)-parm1(i)).gt.acc*dabs(parm(i))+rsmall) ier=0 131 continue if (istop.eq.2) ier=0 if (ier.eq.1) goto 9001 c c next iteration c 140 delf=2.d0*(fval-fval1) goto 2000 c c successful exits c 9001 ier=1 goto 9000 9002 ier=2 goto 9000 9003 ier=3 goto 9000 9004 ier=4 goto 9000 c c failure exits c 9015 ier=-1 goto 9000 9016 ier=-6 goto 9000 9018 ier=-3 goto 9000 9019 ier=-9 goto 9000 9000 return end c subroutine linesr(parm,np,fval,dir,pm1,st0,acc,parm1,*) c---------------------------------------------------------------------- c linesearch: computes new parm and fval along dir. c---------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) dimension parm(np),dir(np),parm1(np) call brack(parm,np,fval,dir,pm1,st0,st1,f1,st2,f2,parm1,iedge) if (iedge.eq.0) then call quadin(parm,np,fval,dir,pm1,st0,st1,f1,st2,f2,acc,parm1) else call move(parm1,np,parm) end if return end c subroutine brack(parm,np,f0,dir,pm1,st0,st1,f1,st2,f2,parm1,iedge) c---------------------------------------------------------------------- c linesearch: find initial bracket [st0,st1,st2] c---------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) dimension parm(np),dir(np),parm1(np) c c compute second point: take initial stepsize in search direction c st1=st0 st0=0. 10 call step(parm,np,st1,dir,parm1) call func(parm1,np,f1,*12) goto 20 c ---failure of function evaluation 12 st1 = st1/2.d0 if (st1.lt.1.d-31) goto 99 goto 10 c c if new point is better, swap c 20 iedge = 0 if (f1*pm1.gt.f0*pm1) then f=f1 f1=f0 f0=f st=st1 st1=st0 st0=st end if c c compute third point: deepen step in improving direction c st2=2.0*st0-st1 30 call step(parm,np,st2,dir,parm1) call func(parm1,np,f2,*32) cc write (*,*) 'brack: fval=',f2 if (f2*pm1.le.f0*pm1) goto 90 goto 40 c ---failure of function evaluation 32 shold = st2 st2=(st2+st0)/2. if (shold.eq.st2) goto 95 goto 30 c c if improvement, deepen bracket even more c 40 f1=f0 f0=f2 st1=st0 st0=st2 st2=3.*st0-2.*st1 goto 30 c c iedge indicates that the linesearch lead to points where the c function was undefined (do not use inv. quadratic interpolation) 95 iedge = 1 c c normal exit: no improvement found, use inv. quadr. interpolation 90 call step(parm,np,st0,dir,parm1) return c c failure 99 iedge = 1 return end c subroutine quadin(parm,np,f0,dir,pm1,st0,st1,f1,st2,f2,acc,parm1) c---------------------------------------------------------------------- c linesearch: refine given bracket [st0,st1,st2] by inverse c quadratic interpolation c---------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) dimension parm(np),dir(np),parm1(np) data nlnsr/10/ c improv=0 ilnsr=0 dd=dotp(dir,np,dir) if (dd.eq.0) goto 101 stpac=acc/dd c c iteration (max. two improvements or nlnsr steps) c 10 if (ilnsr.gt.nlnsr) goto 100 ilnsr=ilnsr+1 c c inverse quadratic interpolation c stp=st0**2*(f1-f2)+st1**2*(f2-f0)+st2**2*(f0-f1) std=(2.*(st0*(f1-f2)+st1*(f2-f0)+st2*(f0-f1))) c c check for stepsize convergence c if (dabs(std).lt.1.d-31) goto 100 stp=stp/std if (dabs(stp-st0).le.stpac) goto 100 c c define new bracket c call step(parm,np,stp,dir,parm1) call func(parm1,np,fnew,*100) cc write (*,*) 'quadi: fval=',fnew if (fnew*pm1.gt.f0*pm1) then improv=1 if ((stp-st0)*(st1-st0).gt.0.) then st2=st0 f2=f0 st0=stp f0=fnew else st1=st0 f1=f0 st0=stp f0=fnew end if else if (improv.eq.1) goto 100 if ((stp-st0)*(st1-st0).gt.0.) then st1=stp f1=fnew else st2=stp f2=fnew end if end if goto 10 c c exit c 100 call step(parm,np,st0,dir,parm1) call move(parm1,np,parm) return 101 write (*,*) 'quadin: dir=0' call move(parm1,np,parm) return end c c c subroutine itlog(parm,np,fval,fpnorm,grdir,iter,text) c--------------------------------------------------------------------- implicit real*8 (a-h,o-z) implicit integer (i-n) ccc$include:'sml.dim' include 'ssmlmnp.inc' c$include:'sml.dim' c....................................... c....................................... dimension parm(np),parm1(10) character*4 text character*20 runtag character*21 tdate character*40 spfile,pafile common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns, * nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp common / doffs / ns1,ns2,ns3,ns4,ns5,maps(maxalt,maxalt), * locs(maxnpp) common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp) common / doutp / kanal0,ispool,ipause common / dival / ival common / dtag / runtag common / dfile / spfile,pafile c if (text.eq.'quit') return call tstamp(tdate) write (*,10) text,iter,fval,fpnorm,grdir,ival if (nb.gt.0) write (*,11) (parm(i),i=1,nb) if (ns.gt.0) then c *** undo transformation *** do 5 k=1,ns k1=nb+k i1=loc1(k1) if ((i1.le.ns2).or.(i1.gt.ns3.and.i1.le.ns4)) then c ---std.dev's nu and alfa: parm1(k)=exp(parm(k1)) else c ---rho, correl nu and alfa: parm1(k)=(exp(parm(k1))-1)/(exp(parm(k1))+1) end if 5 continue write (*,12) (parm1(k),k=1,ns) end if write (*,13) tdate,runtag if (ispool.ne.0) then write (kanal0,10) text,iter,fval,fpnorm,grdir,ival if (nb.gt.0) write (kanal0,11) (parm(i),i=1,nb) if (ns.gt.0) write (kanal0,12) (parm1(k),k=1,ns) write (kanal0,13) tdate,runtag close(kanal0) open(kanal0,file=spfile,status='old',access='append') end if 10 format(/1x,a4,'-iteration',i3,':', * ' lik=',d12.5,', gradnm=',d9.3,', grad*dir=',d9.3,i5) 11 format(' beta =',d11.4,5d12.4:/10(6x,6d12.4:/)) 12 format(' sigma=',d11.4,5d12.4:/10(6x,6d12.4:/)) 13 format(' time = ',a21,', ',a20) return end c----------------------------------------------------------------------c c system-dependent routines c c If you cannot find the correct calls for your system, c c use the generic ones, currently commented out by CDUMMY c c c c----------------------------------------------------------------------c subroutine getdat(iy,im,id) c system-dependent routine: SUN/UNIX version plus VAH library integer*2 iy,im,id integer*4 itim(5),idat(3) call timdat(itim,idat) id=idat(1) im=idat(2) iy=idat(3) cvah since timdat/unix returns year in two digits iy=iy+1900 return end subroutine gettim(ih,im,is,ihs) c system-dependent routine: SUN/UNIX version plus VAH library integer*2 ih,im,is,ihs integer*4 itim(5),idat(3) call timdat(itim,idat) ih=itim(1) im=itim(2) is=itim(3) ihs=itim(4) return end cdummy subroutine gettim(ih,im,is,ihs) cdummy c returns 0 for hour, minute, and 10th of sec, and adds 1 to cdummy c previous value of second. cdummy c ideally, should return the time as cdummy c ih=hour, im=minute, is=second, ihs=hundredth of a sec cdummy integer*2 ih,im,is,ihs cdummy integer itemp cdummy save itemp cdummy ih=0 cdummy im=0 cdummy itemp=itemp+1 cdummy is=itemp cdummy ihs=0 cdummy return cdummy end cdummy subroutine getdat(iy,im,id) cdummy integer*2 iy,im,id cdummy c ideally, should return date as cdummy c iy=year, im=month, id=day of the month cdummy iy=1991 cdummy im=1 cdummy id=1 cdummy return cdummy end cdummy subroutine parsecl(cline,istrt,iend,narg) cdummy c ideally, returns the command line when the program was cdummy c called as cline. It also parses cline and returns cdummy c number of distinct arguments, narg, as well as integer cdummy c arrays of length narg, specifying where each argument cdummy c begins (istrt) and ends (iend) cdummy character*(*) cline cdummy integer*4 istrt,iend,narg cdummy istrt=1 cdummy iend=1 cdummy narg=0 cdummy return cdummy end