* * drw3net.f = drw3.f without the specification of the subdirectories * * drw3.f, a rewritten version of vac5.f with additional * subroutines to calculate the steady state, to check for * accuracy, and to calculate summary statistics * * vac5.f, like vac4.f but with exogenous separation before * and after production takes place. * * vac4.f, like vac3.f but with new matching function * * vac3.f, like vac1.f but with variable elasticity in the * matching function * * vac2.f, program with exogenous separation after production * has taken place * * vac1.f, with endogenous vacancies and exogenous seperation * * mp1.f, like d1.f but positive value of being in the matching pool * * d1.f, like depr1.f but only parameterizing the difference * of being matched and being in the pool * * depr1.f, complete basis * * depr2.f, capital choice is made before idiosyncratic shock * is observed * implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap dimension root(3000),roothx(1000),roothw(1000),rootlx(1000), $ rootlw(1000) external fdens,equil,equilcap c c npar is number of parameters for the policy functions c open(33,file='drwcoef.dat',status='unknown') open(34,file='drwcoef.new',status='unknown') open(41,file='drw.dat',status='unknown') open(36,file='cheb0.dat',status='unknown') open(39,file='herm.dat',status='unknown') open(40,file='herm.pro',status='unknown') open(55,file='for055.dat',status='unknown') open(99,file='for099.dat',status='unknown') c c read chebyshev and hermite nodes c read(36,*) (root(j),j=1,870) read(39,*) (roothx(j),j=1,780) read(40,*) (roothw(j),j=1,780) read(33,*) ((aold(j,k),j=1,npar),k=1,4) read(41,*) i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson, $ alpha,bta,depr,gam,ro,sig,xvalue,pmatch,eps, $ bl1,bh1,bl3,bh3,fstan,acc,xacc,itotal,i2bound,bvalue,share, $ prob,eta,auto2,abba,auto1,istst,it1,it2,it3,imom,icreate, $ impulse,icap pi = 3.141592654 c c c construct Simpson weights (note that the number of Simpson c nodes has to be odd) c c hsw(1) = 1./3. hsw(nsimpson) = 1./3. itel = -1 do 905 i = 2,nsimpson-1 if(itel.lt.0) then hsw(i) = 4./3. else hsw(i) = 2./3. endif itel = -1*itel 905 continue write(*,'(10f7.3)') (( aold(i,kk),i=1,npar),kk=1,4) if(nhermite.gt.30) then write(*,*) 'not enough nodes in herm.dat and herm.pro' stop endif do 910 i = 1,nhermite istart = (nhermite-5)*30 hhx(i) = roothx(istart+i) hhw(i) = roothw(istart+i) 910 continue c c numerically calculate the integral that will be useful c for calculating the interest rate c asum = 0. yy = (8.*fstan)/dble(nsimpson-1.) do 1911 i = 1,nsimpson qx = dble(i-1)*yy - 4.*fstan qw = hsw(i)*yy*fdens(qx) asum = asum + exp(qx)*qw 1911 continue c c calculate the steady state c capstst = (1.-bta*(1.-depr))/(bta*alpha*asum) capstst = capstst**(1./(alpha-1.)) write(*,*) 'capstst (incorporates asum)', capstst c c Construct Gridpoints c c stbar(.,.) contains the boundaries c stbar(1,1) = bl1 stbar(1,2) = bh1 stbar(2,1) = -i2bound*(sig**2./(1.-ro**2.))**0.5 stbar(2,2) = -stbar(2,1) stbar(3,1) = bl3 stbar(3,2) = bh3 c c below we first find the appropriate roots, and then c find the corresponding unscaled numbers c do 10 is1 = 1,is1max istart = (is1max-2)*30 st(is1,1) = root(istart+is1) 10 continue do 20 is1 = 1,is1max bb = st(is1,1) stu(is1,1) = scinv(bb,1) 20 continue do 11 is2 = 1,is2max istart = (is2max-2)*30 st(is2,2) = root(istart+is2) 11 continue do 21 is2 = 1,is2max bb = st(is2,2) stu(is2,2) = scinv(bb,2) 21 continue do 12 is3 = 1,is3max istart = (is3max-2)*30 st(is3,3) = root(istart+is3) 12 continue do 22 is3 = 1,is3max bb = st(is3,3) stu(is3,3) = scinv(bb,3) 22 continue c c Construct the chebyshev polynomials at the grid points c (Since these have to be calculated many times, it is c worthwhile to save them c do 40 is3 = 1,is3max do 40 is2 = 1,is2max do 40 is1 = 1,is1max jt = ispoint(is1,is2,is3) ss1 = st(is1,1) ss2 = st(is2,2) ss3 = st(is3,3) do 41 j1 = 1,i1max+1 do 41 j2 = 1,i2max+1 do 41 j3 = 1,i3max+1 i1 = j1-1 i2 = j2-1 i3 = j3-1 j = ipoint(i1,i2,i3) cheb(jt,j) = hh(ss1,ss2,ss3,i1,i2,i3) 41 continue 40 continue write(*,*) ipoint(i1max,i2max,i3max),npar jmax = is1max*is2max*is3max DO 110 J = 1,NPAR XX(J,J) = 0.0 110 CONTINUE DO 130 J = 1,NPAR DO 130 I = 1,JMAX XX(J,J) = XX(J,J) + CHEB(I,J)*CHEB(I,J) 130 CONTINUE c c c note that the orthogonality property of the chebyshev polynomial c is used to calculate the inverse of the x'x matrix c c DO 135 J = 1,NPAR XX(J,J) = 1./XX(J,J) 135 CONTINUE c c c Iterating on Consumer Problem c c do 300 iter = 1,itotal call sbta1() accur = 0. do 310 kk = 1,4 do 310 i = 1,npar if(abs(aold(i,kk)).gt.0.001) then accur = accur + abs( (aold(i,kk)-anew(i,kk))/aold(i,kk) ) else accur = accur + abs( (aold(i,kk)-anew(i,kk)) ) endif 310 continue write(*,*) accur,iter if(accur.lt.acc.or.iter.eq.itotal) goto 305 do 311 kk = 1,4 do 311 i = 1,npar aold(i,kk) = eps*anew(i,kk) + (1.-eps)*aold(i,kk) 311 continue 300 continue 305 continue do 600 kk = 1,4 write(34,*) (aold(j,kk),j=1,npar) 600 continue if(istst.eq.1) call sbta2() if(impulse.eq.1) call sbta3() if(it1.eq.1) call sacc1() if(it2.eq.1) call sacc2() if(it3.eq.1) call sacc3() if(imom .eq.1) call sbta4() stop end *************************************************************************** * * Sbta1(); Subroutine to calculate at each grid point numerically * the four conditional expectations and to project the calculated * values on the Chebyshev polynomial terms * *************************************************************************** subroutine sbta1() implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap dimension y(20000,4),xy(4000) external fdens,equil,equilcap do 10 is1 = 1,is1max do 10 is2 = 1,is2max do 10 is3 = 1,is3max jt = ispoint(is1,is2,is3) s1old = st(is1,1) s2old = st(is2,2) s3old = st(is3,3) pemployed = stu(is3,3) cap = stu(is1,1) the = exp(stu(is2,2)) gg = polgrid(is1,is2,is3,2) ww = polgrid(is1,is2,is3,3) wf = polgrid(is1,is2,is3,4) exp1 = polgrid(is1,is2,is3,1) call poltrunc() call polprod() expold = exp1 theold = the cap = capnew pemployed = pempnew call polexp() y(jt,1) = log(expcons) y(jt,2) = log(expmatch) y(jt,3) = log(expwork) y(jt,4) = log(expfirm) 10 continue JMAX = IS1MAX*IS2MAX*IS3MAX do 600 kk = 1,4 DO 110 J = 1,NPAR XY(J) = 0.0 110 CONTINUE DO 140 J = 1,NPAR DO 140 I = 1,JMAX XY(J) = XY(J) + CHEB(I,J)*Y(I,kk) 140 CONTINUE DO 145 J = 1,NPAR ANEW(J,KK) = XX(J,J)*XY(J) 145 CONTINUE 600 continue end *************************************************************************** * * Sbta2(); Subroutine to calculate the steady state values * *************************************************************************** subroutine sbta2() implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap external fdens,equil,equilcap pemployed = stu(2,3) cap = stu(2,1) theold = exp(stu(2,2)) idum = -200597 do 10 itime = 1,10000 shock = 0. c write(*,*) cap,pemployed the = ro*log(theold)+shock*sig s2old = sca(the,2) the = exp(the) s1old = sca(cap,1) s3old = sca(pemployed,3) gg = pol(s1old,s2old,s3old,2) ww = pol(s1old,s2old,s3old,3) wf = pol(s1old,s2old,s3old,4) exp1 = pol(s1old,s2old,s3old,1) call poltrunc() call polprod() accz = abs(pempnew-pemployed)+abs(truncpr-old2) $ + abs(pmatch-old3)+abs(vacpr-old4) old2 = truncpr old3 = pmatch old4 = vacpr if(itime.eq.1000.or.accz.le.0.0001) then write(*,701) unemp/pemployed write(*,702) vac/pemployed write(*,703) truncpr write(*,704) pmatch write(*,705) vacpr write(*,706) bvalue write(*,707) xvalue write(*,708) prob write(*,709) eta write(*,710) fstan write(55,701) unemp/pemployed write(55,702) vac/pemployed write(55,703) truncpr write(55,704) pmatch write(55,705) vacpr write(55,706) bvalue write(55,707) xvalue write(55,708) prob write(55,709) eta write(55,710) fstan goto 11 endif * * calculate next period's fraction of people employed * pemployed = pempnew cap = capnew theold = the 10 continue 11 continue 701 format('un/pemployed: ',f14.6) 702 format('vac/pemployed: ',f14.6) 703 format('truncpr: ',f14.6) 704 format('pmatch: ',f14.6) 705 format('pvac: ',f14.6) 706 format('bvalue: ',f14.6) 707 format('xvalue: ',f14.6) 708 format('prob: ',f14.6) 709 format('eta: ',f14.6) 710 format('fstan: ',f14.6) end *************************************************************************** * * Sbta3(); Subroutine to calculate impulse response function * *************************************************************************** subroutine sbta3() implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap external fdens,equil,equilcap open(31,file='impu.dat', $status='unknown') pemployed = stu(is3max,3) cap = stu(is1max,1) theold = exp(stu(is2max,2)) idum = -200597 do 10 itime = 1,1000 shock = 0. if(itime.eq.800) shock = -3. the = ro*log(theold)+shock*sig s2old = sca(the,2) the = exp(the) s1old = sca(cap,1) s3old = sca(pemployed,3) gg = pol(s1old,s2old,s3old,2) ww = pol(s1old,s2old,s3old,3) wf = pol(s1old,s2old,s3old,4) exp1 = pol(s1old,s2old,s3old,1) call poltrunc() call polprod() rotot = 1.-auto1*auto2*truncpr roex = truncpr*(1.-auto1*auto2) repost = roex yx2 = pempnew yx3 = (rotot - repost*vacpr) $ * 2.*pemployed/(pemployed+pempnew) yx4 = vacpr*((vac/pemployed)-repost) $ * 2.*pemployed/(pemployed+pempnew) rentadj = dble(icap)*truncpr*auto1 + (1.-dble(icap))*1. efrental = rentadj*rental if(itime.ge.781) then write(31,'(i6,17f12.7)') itime,the,prod,yx2,yx3,yx4, $ pmatch,vacpr,vac,cap, $ efrental,rental,truncpr,cons1,bta*consold/cons1, $ gg,ww,wf endif * calculate next period's fraction of people employed * pemployed = pempnew cap = capnew theold = the consold = cons1 10 continue 11 continue end *************************************************************************** * * SACC1(); Subroutine to check for accuracy for a range of * different values for the capital stock * *************************************************************************** subroutine sacc1() implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap external fdens,equil,equilcap open(91,file $ ='for091.dat',status='unknown') open(92,file $ ='for092.dat',status='unknown') open(93,file $ ='for093.dat',status='unknown') open(94,file $ ='for094.dat',status='unknown') open(95,file $ ='for095.dat',status='unknown') open(96,file $ ='for096.dat',status='unknown') open(97,file $ ='for097.dat',status='unknown') open(98,file $ ='for098.dat',status='unknown') do 10 is2 = 1,is2max do 10 is3 = 1,is3max write(*,*) is2,is3 do 10 is1 = 1,51 cap = stbar(1,1)+dble(is1-1)*(stbar(1,2)-stbar(1,1))/51. s1old = sca(cap,1) s2old = st(is2,2) s3old = st(is3,3) pemployed = stu(is3,3) the = exp(stu(is2,2)) exp1 = pol(s1old,s2old,s3old,1) gg = pol(s1old,s2old,s3old,2) ww = pol(s1old,s2old,s3old,3) wf = pol(s1old,s2old,s3old,4) pol1 = exp1 pol2 = gg pol3 = ww pol4 = wf call poltrunc() call polprod() expold = exp1 theold = the cap = capnew pemployed = pempnew call polexp() write(91,'(f14.6)') pol1 write(92,'(f14.6)') pol2 write(93,'(f14.6)') pol3 write(94,'(f14.6)') pol4 write(95,'(f14.6)') expcons write(96,'(f14.6)') expmatch write(97,'(f14.6)') expwork write(98,'(f14.6)') expfirm 10 continue close(unit=91) close(unit=92) close(unit=93) close(unit=94) close(unit=95) close(unit=96) close(unit=97) close(unit=98) end *************************************************************************** * * SACC2(); Subroutine to check for accuracy for a range of * different values for the productivity shock * *************************************************************************** subroutine sacc2() implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap external fdens,equil,equilcap open(81,file $ ='for081.dat',status='unknown') open(82,file $ ='for082.dat',status='unknown') open(83,file $ ='for083.dat',status='unknown') open(84,file $ ='for084.dat',status='unknown') open(85,file $ ='for085.dat',status='unknown') open(86,file $ ='for086.dat',status='unknown') open(87,file $ ='for087.dat',status='unknown') open(88,file $ ='for088.dat',status='unknown') do 10 is1 = 1,is1max do 10 is3 = 1,is3max write(*,*) is1,is3 do 10 is2 = 1,51 the = stbar(2,1)+dble(is2-1)*(stbar(2,2)-stbar(2,1))/51. s2old = sca(the,2) s1old = st(is1,1) s3old = st(is3,3) pemployed = stu(is3,3) cap = stu(is1,1) the = exp(the) exp1 = pol(s1old,s2old,s3old,1) gg = pol(s1old,s2old,s3old,2) ww = pol(s1old,s2old,s3old,3) wf = pol(s1old,s2old,s3old,4) pol1 = exp1 pol2 = gg pol3 = ww pol4 = wf call poltrunc() call polprod() expold = exp1 theold = the cap = capnew pemployed = pempnew call polexp() write(81,'(f14.6)') pol1 write(82,'(f14.6)') pol2 write(83,'(f14.6)') pol3 write(84,'(f14.6)') pol4 write(85,'(f14.6)') expcons write(86,'(f14.6)') expmatch write(87,'(f14.6)') expwork write(88,'(f14.6)') expfirm 10 continue close(unit=81) close(unit=82) close(unit=83) close(unit=84) close(unit=85) close(unit=86) close(unit=87) close(unit=88) end *************************************************************************** * * SACC3(); Subroutine to check for accuracy for a range of * different values for the fraction of employed matches * *************************************************************************** subroutine sacc3() implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap external fdens,equil,equilcap open(71,file $ ='for071.dat',status='unknown') open(72,file $ ='for072.dat',status='unknown') open(73,file $ ='for073.dat',status='unknown') open(74,file $ ='for074.dat',status='unknown') open(75,file $ ='for075.dat',status='unknown') open(76,file $ ='for076.dat',status='unknown') open(77,file $ ='for077.dat',status='unknown') open(78,file $ ='for078.dat',status='unknown') do 10 is2 = 1,is2max do 10 is1 = 1,is1max write(*,*) is2,is1 do 10 is3 = 1,51 pemployed = stbar(3,1)+dble(is3-1)*(stbar(3,2)-stbar(3,1))/51. s3old = sca(pemployed,3) s2old = st(is2,2) s1old = st(is1,1) cap = stu(is1,1) the = exp(stu(is2,2)) exp1 = pol(s1old,s2old,s3old,1) gg = pol(s1old,s2old,s3old,2) ww = pol(s1old,s2old,s3old,3) wf = pol(s1old,s2old,s3old,4) pol1 = exp1 pol2 = gg pol3 = ww pol4 = wf call poltrunc() call polprod() expold = exp1 theold = the cap = capnew pemployed = pempnew call polexp() write(71,'(f14.6)') pol1 write(72,'(f14.6)') pol2 write(73,'(f14.6)') pol3 write(74,'(f14.6)') pol4 write(75,'(f14.6)') expcons write(76,'(f14.6)') expmatch write(77,'(f14.6)') expwork write(78,'(f14.6)') expfirm 10 continue close(unit=71) close(unit=72) close(unit=73) close(unit=74) close(unit=75) close(unit=76) close(unit=77) close(unit=78) end *************************************************************************** * * Sbta4(); Subroutine to calculate the summary statistics * *************************************************************************** subroutine sbta4() implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap common/stat/ y(500,10),cor1(2,9),cor2(2,9),ave(10),stan(10), $ cor3(2,9),autocor(10,10),nobs,idump dimension cor1ave(2,9),cor2ave(2,9),cor3ave(1,9),ratave(10), $ autoave(10,10),hpy(150),hpt(150),hpd(150),hpv(150,3) dimension cor1se(2,9),cor2se(2,9),cor3se(1,9),ratse(10), $ autose(10,10) external fdens,equil,equilcap idum = -200597 montetot = 100 s = 1600 nobs = 67 idump = 200 nw = 3 iopt = 0 do 211 jj = 1,8 ratave(jj) = 0. ratse(jj) = 0. 211 continue do 212 jj = 1,2 do 212 kk = 1,9 cor1ave(jj,kk) = 0. cor2ave(jj,kk) = 0. cor3ave(1,kk) = 0. cor1se(jj,kk) = 0. cor2se(jj,kk) = 0. cor3se(jj,kk) = 0. 212 continue do 213 jj = 1,8 do 213 kk = 1,4 autoave(jj,kk) = 0. autose(jj,kk) = 0. 213 continue do 115 monte = 1,montetot write(*,*) monte pemployed = stu(2,3) cap = stu(2,1) theold = exp(stu(2,2)) do 10 itime = 1,idump+nobs shock = gasdev(idum) the = ro*log(theold)+shock*sig s2old = sca(the,2) the = exp(the) s1old = sca(cap,1) s3old = sca(pemployed,3) gg = pol(s1old,s2old,s3old,2) ww = pol(s1old,s2old,s3old,3) wf = pol(s1old,s2old,s3old,4) exp1 = pol(s1old,s2old,s3old,1) call poltrunc() call polprod() rotot = 1.-auto1*auto2*truncpr roex = truncpr*(1.-auto1*auto2) repost = roex y(itime,1) = prod y(itime,2) = pempnew y(itime,3) = (rotot - repost*vacpr) $ * 2.*pemployed/(pemployed+pempnew) y(itime,4) = vacpr*((vac/pemployed)-repost) $ * 2.*pemployed/(pemployed+pempnew) y(itime,5) = cons1 y(itime,6) = capnew-(1.-depr)*cap y(itime,7) = prod/cap**alpha y(itime,8) = the if(itime.gt.idump) then do 444 ikk = 1,8 if(y(itime,ikk).le.0.) then write(*,*) itime,ikk write(99,*) itime,ikk y(itime,ikk) = 0.0001 endif y(itime,ikk) = log(y(itime,ikk)) 444 continue endif pemployed = pempnew cap = capnew theold = the 10 continue do 900 j = 1,8 do 910 i = 1,nobs hpy(i) = y(i+idump,j) 910 continue call hpfilt(hpy,hpt,hpd,hpv,nobs,nw,s,iopt) do 920 i = 1,nobs y(i+idump,j) = hpd(i) 920 continue 900 continue call avesub() call stansub() call corsub() call autosub() ratave(1) = ratave(1) + stan(1) ratave(2) = ratave(2) + stan(2)/stan(1) ratave(3) = ratave(3) + stan(3)/stan(2) ratave(4) = ratave(4) + stan(4)/stan(2) ratave(8) = ratave(8) + stan(1)/stan(8) ratse(1) = ratse(1) + stan(1)**2. ratse(2) = ratse(2) + (stan(2)/stan(1))**2. ratse(3) = ratse(3) + (stan(3)/stan(2))**2. ratse(4) = ratse(4) + (stan(4)/stan(2))**2. ratse(8) = ratse(8) + (stan(1)/stan(8))**2. do 11 jj = 5,7 ratave(jj) = ratave(jj) + stan(jj)/stan(1) ratse(jj) = ratse(jj) + (stan(jj)/stan(1))**2. 11 continue do 12 kk = 1,9 cor3ave(1,kk) = cor3ave(1,kk) + cor3(1,kk) cor3se(1,kk) = cor3se(1,kk) + cor3(1,kk)**2. do 12 jj = 1,2 cor1ave(jj,kk) = cor1ave(jj,kk) + cor1(jj,kk) cor2ave(jj,kk) = cor2ave(jj,kk) + cor2(jj,kk) cor1se(jj,kk) = cor1se(jj,kk) + cor1(jj,kk)**2. cor2se(jj,kk) = cor2se(jj,kk) + cor2(jj,kk)**2. 12 continue do 13 jj = 1,8 do 13 kk = 1,4 autoave(jj,kk) = autoave(jj,kk) + autocor(jj,kk) autose(jj,kk) = autose(jj,kk) + autocor(jj,kk)**2. 13 continue 115 continue do 111 jj = 1,8 ratave(jj) = ratave(jj)/dble(montetot) ratse(jj) = ratse(jj)/dble(montetot) ratse(jj) = (ratse(jj)-ratave(jj)**2.)**0.5 111 continue do 112 kk = 1,9 cor3ave(1,kk) = cor3ave(1,kk)/dble(montetot) cor3se(1,kk) = cor3se(1,kk)/dble(montetot) cor3se(1,kk) = (cor3se(1,kk)-cor3ave(1,kk)**2.)**0.5 do 112 jj = 1,2 cor1ave(jj,kk) = cor1ave(jj,kk)/dble(montetot) cor2ave(jj,kk) = cor2ave(jj,kk)/dble(montetot) cor1se(jj,kk) = cor1se(jj,kk)/dble(montetot) cor2se(jj,kk) = cor2se(jj,kk)/dble(montetot) cor1se(jj,kk) = (cor1se(jj,kk)-cor1ave(jj,kk)**2.)**0.5 cor2se(jj,kk) = (cor2se(jj,kk)-cor2ave(jj,kk)**2.)**0.5 112 continue do 113 jj = 1,8 do 113 kk = 1,4 autoave(jj,kk) = autoave(jj,kk)/dble(montetot) autose(jj,kk) = autose(jj,kk)/dble(montetot) autose(jj,kk) = (autose(jj,kk)-autoave(jj,kk)**2.)**0.5 113 continue write(*,701) ratave(1) write(*,702) ratave(2) write(*,703) ratave(3) write(*,704) ratave(4) write(*,705) ratave(5) write(*,706) ratave(6) write(*,707) ratave(7) write(*,708) ratave(8) write(*,801) (autoave(1,kk),kk=1,4) write(*,802) (autoave(2,kk),kk=1,4) write(*,803) (autoave(3,kk),kk=1,4) write(*,804) (autoave(4,kk),kk=1,4) write(*,805) (autoave(5,kk),kk=1,4) write(*,806) (autoave(6,kk),kk=1,4) write(*,807) (autoave(7,kk),kk=1,4) write(*,808) (autoave(8,kk),kk=1,4) write(*,721) (cor1ave(1,kk),kk=1,9) write(*,722) (cor1ave(2,kk),kk=1,9) write(*,723) (cor2ave(1,kk),kk=1,9) write(*,724) (cor2ave(2,kk),kk=1,9) write(*,725) (cor3ave(1,kk),kk=1,9) write(55,701) ratave(1) write(55,702) ratave(2) write(55,703) ratave(3) write(55,704) ratave(4) write(55,705) ratave(5) write(55,706) ratave(6) write(55,707) ratave(7) write(55,708) ratave(8) write(55,801) (autoave(1,kk),kk=1,4) write(55,802) (autoave(2,kk),kk=1,4) write(55,803) (autoave(3,kk),kk=1,4) write(55,804) (autoave(4,kk),kk=1,4) write(55,805) (autoave(5,kk),kk=1,4) write(55,806) (autoave(6,kk),kk=1,4) write(55,807) (autoave(7,kk),kk=1,4) write(55,808) (autoave(8,kk),kk=1,4) write(55,721) (cor1ave(1,kk),kk=1,9) write(55,722) (cor1ave(2,kk),kk=1,9) write(55,723) (cor2ave(1,kk),kk=1,9) write(55,724) (cor2ave(2,kk),kk=1,9) write(55,725) (cor3ave(1,kk),kk=1,9) write(*,701) ratse(1) write(*,702) ratse(2) write(*,703) ratse(3) write(*,704) ratse(4) write(*,705) ratse(5) write(*,706) ratse(6) write(*,707) ratse(7) write(*,708) ratse(8) write(*,801) (autose(1,kk),kk=1,4) write(*,802) (autose(2,kk),kk=1,4) write(*,803) (autose(3,kk),kk=1,4) write(*,804) (autose(4,kk),kk=1,4) write(*,805) (autose(5,kk),kk=1,4) write(*,806) (autose(6,kk),kk=1,4) write(*,807) (autose(7,kk),kk=1,4) write(*,808) (autose(8,kk),kk=1,4) write(*,721) (cor1se(1,kk),kk=1,9) write(*,722) (cor1se(2,kk),kk=1,9) write(*,723) (cor2se(1,kk),kk=1,9) write(*,724) (cor2se(2,kk),kk=1,9) write(*,725) (cor3se(1,kk),kk=1,9) write(55,701) ratse(1) write(55,702) ratse(2) write(55,703) ratse(3) write(55,704) ratse(4) write(55,705) ratse(5) write(55,706) ratse(6) write(55,707) ratse(7) write(55,708) ratse(8) write(55,801) (autose(1,kk),kk=1,4) write(55,802) (autose(2,kk),kk=1,4) write(55,803) (autose(3,kk),kk=1,4) write(55,804) (autose(4,kk),kk=1,4) write(55,805) (autose(5,kk),kk=1,4) write(55,806) (autose(6,kk),kk=1,4) write(55,807) (autose(7,kk),kk=1,4) write(55,808) (autose(8,kk),kk=1,4) write(55,721) (cor1se(1,kk),kk=1,9) write(55,722) (cor1se(2,kk),kk=1,9) write(55,723) (cor2se(1,kk),kk=1,9) write(55,724) (cor2se(2,kk),kk=1,9) write(55,725) (cor3se(1,kk),kk=1,9) 701 format(' output stan. dev. ',2f14.4) 702 format(' employment/output ',2f14.4) 703 format(' destruction/employment',2f14.4) 704 format(' creation/employment',2f14.4) 705 format(' consumption/output ',2f14.4) 706 format(' investment/output ',2f14.4) 707 format(' solow/output ',2f14.4) 708 format('output/productivity ',2f14.4) 721 format('destruction & output ',9f6.2) 722 format(' creation & output ',9f6.2) 723 format('destruction & employment ',9f6.2) 724 format(' creation & employment ',9f6.2) 725 format(' creation & destruction',9f6.2) 801 format(' output',4f7.3) 802 format(' employement',4f7.3) 803 format(' destruction',4f7.3) 804 format(' creation',4f7.3) 805 format(' consumption',4f7.3) 806 format(' investment',4f7.3) 807 format(' solow',4f7.3) 808 format('productivity',4f7.3) end *************************************************************************** * * AVESUB(); subroutine to calculate the average * *************************************************************************** subroutine avesub() implicit real*8(a-h,o-z),integer(i-n) common/stat/ y(500,10),cor1(2,9),cor2(2,9),ave(10),stan(10), $ cor3(2,9),autocor(10,10),nobs,idump do 10 j = 1,8 ave(j) = 0. 10 continue do 11 j = 1,8 do 11 i = idump+1,idump+nobs ave(j) = ave(j) + y(i,j) 11 continue do 12 j = 1,8 ave(j) = ave(j)/dble(nobs) 12 continue end *************************************************************************** * * STANSUB(); subroutine to calculate the standard deviation * *************************************************************************** subroutine stansub() implicit real*8(a-h,o-z),integer(i-n) common/stat/ y(500,10),cor1(2,9),cor2(2,9),ave(10),stan(10), $ cor3(2,9),autocor(10,10),nobs,idump do 10 j = 1,8 stan(j) = 0. 10 continue do 11 j = 1,8 do 11 i = idump+1,idump+nobs stan(j) = stan(j) + (y(i,j)-ave(j))**2. 11 continue do 12 j = 1,8 stan(j) = (stan(j)/dble(nobs))**.5 12 continue end *************************************************************************** * * CORSUB(); subroutine to calculate the correlation coefficient * *************************************************************************** subroutine corsub() implicit real*8(a-h,o-z),integer(i-n) common/stat/ y(500,10),cor1(2,9),cor2(2,9),ave(10),stan(10), $ cor3(2,9),autocor(10,10),nobs,idump dimension cov1(2,9),cov2(2,9),cov3(1,9) do 10 k = 1,9 cov3(1,k) = 0.0 do 10 j = 1,2 cov1(j,k) = 0.0 cov2(j,k) = 0.0 10 continue do 15 k = 1,9 if(k.le.5) then ibegin = idump+1+5-k iend = idump+nobs else ibegin = idump+1 iend = idump+nobs+5-k endif do 11 i = ibegin,iend cov3(1,k) = cov3(1,k) $ + (y(i-5+k,4)-ave(4))*(y(i,3)-ave(3)) do 11 j = 1,2 cov1(j,k) = cov1(j,k) $ + (y(i-5+k,j+2)-ave(j+2))*(y(i,1)-ave(1)) cov2(j,k) = cov2(j,k) $ + (y(i-5+k,j+2)-ave(j+2))*(y(i,2)-ave(2)) 11 continue cor3(1,k) = cov3(1,k)/dble(nobs-abs(k-5)) cor3(1,k) = cor3(1,k)/(stan(3)*stan(4)) do 12 j = 1,2 cor1(j,k) = cov1(j,k)/dble(nobs-abs(k-5)) cor2(j,k) = cov2(j,k)/dble(nobs-abs(k-5)) cor1(j,k) = cor1(j,k)/(stan(j+2)*stan(1)) cor2(j,k) = cor2(j,k)/(stan(j+2)*stan(2)) 12 continue 15 continue end *************************************************************************** * * AUTOSUB(); subroutine to calculate autocorrelation coefficients * *************************************************************************** subroutine autosub() implicit real*8(a-h,o-z),integer(i-n) common/stat/ y(500,10),cor1(2,9),cor2(2,9),ave(10),stan(10), $ cor3(2,9),autocor(10,10),nobs,idump dimension autocov(10,10) do 10 k = 1,4 do 10 j = 1,8 autocov(j,k) = 0. 10 continue do 11 k = 1,4 do 11 j = 1,8 do 12 i = idump+1+k,idump+nobs autocov(j,k) = autocov(j,k) $ + (y(i,j)-ave(j))*(y(i-k,j)-ave(j)) 12 continue autocov(j,k) = autocov(j,k)/dble(nobs-k) autocor(j,k) = autocov(j,k)/(stan(j)**2.) 11 continue end *************************************************************************** * * SCA(); subroutine to scale the state variable such that they are * between minus one and one * *************************************************************************** real*8 function sca(x,i) implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap sca = ( 2.*x - stbar(i,2) - stbar(i,1) ) $ / (stbar(i,2)-stbar(i,1)) end *************************************************************************** * * SCINV(); Inverse of SCA() * *************************************************************************** real*8 function scinv(x,i) implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap scinv = ( stbar(i,2) + stbar(i,1) $ + x*(stbar(i,2) - stbar(i,1)) ) / 2. end * * * Define Chebyshev Polynomials * * real*8 function h(x,i) implicit real*8(a-h,o-z),integer(i-n) dimension a(200) a(1) = 1. a(2) = x if(i.le.1) then h = a(i+1) goto 99 endif if(i.gt.1) then do 10 j = 2,i a(j+1) = 2.*x*a(j) - a(j-1) 10 continue h = a(i+1) endif 99 continue end real*8 function hh(s1,s2,s3,i1,i2,i3) implicit real*8(a-h,o-z),integer(i-n) hh = h(s1,i1)*h(s2,i2)*h(s3,i3) end *************************************************************************** * * IPOINT(); Subroutine to allocate a number to a polynomial term * *************************************************************************** integer function ipoint(i1,i2,i3) implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap ipoint = 1 + i2 $ + i1*(i2max+1) $ + i3*(i1max+1)*(i2max+1) c c i* goes from 0 to i*max c c with this convention increasing i1max adds parameters at the c end of the parameter vector c end *************************************************************************** * * IPOINT(); Subroutine to allocate a number to a gridpoint * *************************************************************************** integer function ispoint(is1,is2,is3) implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap ispoint = is1 $ + (is2-1)*is1max $ + (is3-1)*is1max*is2max c c is* goes from 1 to is*max c end *************************************************************************** * * POL(); Function to calculate the parameterized expectation for * arbitrary values of the _scaled_ state variables * *************************************************************************** real*8 function pol(s1,s2,s3,kk) implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap pol = 0. do 41 j1 = 1,i1max+1 do 41 j2 = 1,i2max+1 do 41 j3 = 1,i3max+1 i1 = j1-1 i2 = j2-1 i3 = j3-1 j = ipoint(i1,i2,i3) pol = pol + aold(j,kk)*hh(s1,s2,s3,i1,i2,i3) 41 continue pol = exp(pol) end *************************************************************************** * * POLGRID(); Same as POL() but does it at the grid points for which * we have stored all the Chebyshev terms * *************************************************************************** real*8 function polgrid(is1,is2,is3,kk) implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap jt = ispoint(is1,is2,is3) polgrid = 0. do 20 jj = 1,npar polgrid = polgrid + aold(jj,kk)*cheb(jt,jj) 20 continue polgrid = exp(polgrid) end *************************************************************************** * * RTBIS(); Nonlinear Equation Solver from Numerical Recipes * *************************************************************************** REAL*8 FUNCTION rtbis(func,x1,x2,xacc) IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N) EXTERNAL func PARAMETER (JMAX=400) 60 continue fmid=func(x2) f=func(x1) if(f*fmid.ge.0.) then write(*,*) 'root must be bracketed in rtbis' stop endif if(f.lt.0.)then rtbis=x1 dx=x2-x1 else rtbis=x2 dx=x1-x2 endif do 11 j=1,JMAX dx=dx*.5 xmid=rtbis+dx fmid=func(xmid) if(fmid.le.xacc) rtbis=xmid c if(ibaal.eq.1) write(*,'(3f14.6)') fmid,xmid,dx c if(abs(dx).lt.xacc .or. fmid.eq.0.) return IF(ABS(FMID).LT.XACC)RETURN c c c this is a slight adjustment of the numerical recipes subroutine c due to the numerical errors caused by the numerical integration c procedure there is the possibility that the function is not c continuous and that there is no zero. The n.r. subroutine c would converge anyway. This algorithm would hit the bound c of maximum number of iterations c c 11 continue pause 'too many bisections in rtbis' END C (C) Copr. 1986-92 Numerical Recipes Software D.Q=. *************************************************************************** * * FDENS(); Returns the density of a normal * *************************************************************************** real*8 function fdens(x) implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap fdens = exp( -0.5*x**2./fstan**2.)/(fstan*(2.*pi)**0.5) return end *************************************************************************** * * FNORMCUM(); CDF of a normal * *************************************************************************** real*8 function fnormcum(x) implicit real*8(a-h,o-z),integer(i-n) parameter(npar=27) common/c1/ stu(1000,10),stbar(10,2),st(1000,10), $ cheb(4000,64),xx(npar,npar),hhx(1000),hhw(1000),hsw(1000), $ aold(npar,4),anew(npar,4), $ alpha,bta,depr,gam,ro,sig,pi,asum,xvalue,fstan,share, $ prob,eta,auto1,auto2,abba,bvalue,xacc,acc, $ cap,the,pemployed,capnew,theold,pempnew,prod,cons1,unemp, $ rentemp,rental,pmatch,vacpr,vac,trunc,trunclog,truncpr, $ exp1,expold,expcons,expmatch,expwork,expfirm,gg,ww,wf common/i1/ i1max,i2max,i3max,is1max,is2max,is3max, $ nhermite,nsimpson,icreate,icap external fdens if(x.le.0.) then ff = 0. yy = (x+4.*fstan)/dble(nsimpson-1) do 1912 i = 1,nsimpson qx = dble(i-1)*yy - 4.*fstan qw = hsw(i)*yy*fdens(qx) ff = ff + qw 1912 continue fnormcum = ff else ff = 0. yy = (4.*fstan-x)/dble(nsimpson-1) do 1913 i = 1,nsimpson qx = dble(i-1)*yy + x qw = hsw(i)*yy*fdens(qx) ff = ff + qw 1913 continue fnormcum = 1.-ff endif return end ************************************************************************* * * THE FOLLOWING TWO NUMERICAL RECIPES SUBROUTINES GENERATE RANDOM NUMBERS * ************************************************************************* REAL*8 FUNCTION GASDEV(IDUM) c returns a normally distributed deviate with zero mean and unit c variance, using RAN1(IDUM) as the source of uniform deviates IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N) save iset,gset DATA ISET/0/ IF (ISET.EQ.0) THEN 1 V1=2.*RAN1(IDUM)-1. V2=2.*RAN1(IDUM)-1. R=V1**2.+V2**2. IF(R.GE.1..OR.R.EQ.0.)GO TO 1 FAC=SQRT(-2.*LOG(R)/R) GSET=V1*FAC GASDEV=V2*FAC ISET=1 ELSE GASDEV=GSET ISET=0 ENDIF RETURN END REAL*8 FUNCTION RAN1(IDUM) c c returns a uniform random deviae between 0.0 and 1.0. c set IDUM to any negative value to initialize or reinitialize c the sequence IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N) DIMENSION R(97) PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=1./M1) PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=1./M2) PARAMETER (M3=243000,IA3=4561,IC3=51349) save iff,ix1,ix2,ix3,r DATA IFF /0/ IF(IDUM.LT.0.OR.IFF.EQ.0) THEN IFF = 1 IX1 = MOD(IC1-IDUM,M1) IX1 = MOD(IA1*IX1+IC1,M1) IX2 = MOD(IX1,M2) IX1 = MOD(IA1*IX1+IC1,M1) IX3 = MOD(IX1,M3) DO 11 J=1,97 IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) R(J)=(DBLE(IX1)*DBLE(IX2)*RM2)*RM1 11 CONTINUE IDUM = 1 ENDIF IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) IX3=MOD(IA3*IX3+IC3,M3) J=1+(97*IX3)/M3 IF(J.GT.97.OR.J.LT.1)PAUSE RAN1=R(J) R(J)=(DBLE(IX1)+DBLE(IX2)*RM2)*RM1 RETURN END *************************************************************************** * * HPFILT(); * Kalman smoothing routine for HP filter written by E Prescott. * y=data series, d=deviations from trend, t=trend, n=no. obs, * s=smoothing parameter (eg, 1600 for std HP). * Array v is scratch area and must have dimension at least 3n. * If IOPT=1 and n and s are the same as for the previous call, * the numbers in v are not recomputed. This reduces execution * time by about a third. Note that if this option is exercised, * v cannot be used for other purposes between calls. * *************************************************************************** subroutine hpfilt(y,t,d,v,n,nw,s,iopt) implicit real*8(a-h,o-z),integer(i-n) real*8 y(n),t(n),v(n,nw),d(n) double precision m1,m2,v11,v12,v22,x,z,b11,b12,b22,det,e1,e2 data ss,nn/0.0,0/ c c compute sequences of covariance matrix for f[x(t),x(t-1) | y(