use portlib c c TAX3UNI.F May 12 2001 c c !!! There are several programming errors that we need to fix c (in particular in the tranition laws and the calculation of the tax base) c We don't think that this will do something to the basic message but we c will have to correct these c c c When the taxfunction is linear then the solutions are linear too. c Although the program allows for higher-order approximations it exploits c the linear structure in calculating the breakup levels c c For the nonlinear taxfunctions used in the paper small non-linearities c appear as can be seen from the example input file c c c c integer ngrid,npar,i,j,k,igrid,iter1,iterbb,iter2,iter3,itergov, $itt,iturb,iphi parameter(ngrid=41,npar=3,nsimpson=11) real*8 uhi,ulo,shi,slo,pi,bbbb,phi,barg,xlab, $roexhi,roexlo,roreti,gamup,gamsw,tx,txn,gov,tt real*8 u0,ulolo,ulohi,uhihi,elo,ehi,elotemp0,ehitempup,elotemphi, $u0n,ulolon,ulohin,uhihin,elon,ehin,elotemp0n, $ehitempupn,elotemphin,elohi,elohin, $w0,wlolo,wlohi,whihi,wwo0,wwololo,wwolohi,wwohihi,wfi, $w0old,wloloold,wlohiold,whihiold,wage(6),tax(6), $g,gwolo,gwohi,out,outwo, $zswitchlo,zswitchhi,zaccept,raccept, $zaccept0,zacceptlolo,zacceptlohi,zaccepthihi, $eaccept0,eacceptlolo,eacceptlohi,eaccepthihi,eswitchlo,eswitchhi, $rswitchlo,rswitchhi, $raccept0,racceptlolo,racceptlohi,raccepthihi, $blo,bhi,bloold,bhiold,bfix, $adj1,adj2,det,coef(3,1),qx,qw,yy,zlo,zhi,zupper, $temp0,temp1,temp2,temp3,temp4,temp5,tempu,temps,accur, $ftx,fdens,fnormcum,accurcrit real*8 grid(ngrid,2) real*8 coefhi(npar),xxhi(3,3),xyhi(3,1),yhi(ngrid) real*8 coeflo(npar),xxlo(3,3),xylo(3,1),ylo(ngrid) real*8 hsw(nsimpson) real*8 t1,t2,t3,t4,t5,t6 real*8 baal character(9) today common/par/ uhi,ulo,shi,slo,pi,bbbb,phi,barg,xlab, $roexhi,roexlo,roreti,gamup,gamsw,tx common/simpson/ hsw external ftx,fdens,fnormcum pi = dacos(-1.d0) call date(today) 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 5 i = 2,nsimpson-1 if(itel.lt.0) then hsw(i) = 4./3. else hsw(i) = 2./3. endif itel = -1*itel 5 continue open(31,file='euro1uni.old', $status='unknown') open(32,file='euro1uni.new', $status='unknown') open(33,file='euro1uni.out', $status='unknown') c c a sample for euro1uni.old is attached at the end c c uniform distribution over [ulo-0.5slo,ulo+0.5*slo] and [uhi-0.5*shi,uhi+0.5*shi] c write(33,*) ' UNIFORM DISTRIBUTION ',today adj1 = 0.1d0 adj2 = 0.1d0 adj3 = 0.1d0 read(31,*) (coeflo(j),j=1,npar) read(31,*) (coefhi(j),j=1,npar) read(31,*) blo,bhi,w0,wlolo,wlohi,whihi,wfi, $wwo0,wwololo,wwolohi,wwohihi, $zswitchlo,zswitchhi, $zaccept0,zacceptlolo,zacceptlohi,zaccepthihi, $u0,ulolo,ulohi,uhihi,elo,ehi,elotemp0,ehitempup,elohi ehi=1.d0-u0-ulolo-ulohi-uhihi-elo-elotemp0-ehitempup-elohi do 448 itergov = 1,1 do 448 iphi = 1,2 do 448 iterbb = 1,3 do 448 iter3 = 1,1 do 448 iturb = 1,1 accurcrit = 0.0000001 w0 = 0. wlolo = 0. wlohi = 0. whihi = 0. wfi = 0. wwo0 = 0. wwololo = 0. wwolohi = 0. wwohihi = 0. c if(iter3.gt.1.and.u0+ulolo+ulohi+uhihi.gt.0.05) stop roexlo = 0.0105 roexhi = 0.0105 roreti = 0.005d0 bfix = 0.d0 uhi = 1.0d0 ulo = 1.0d0 shi = 3.046 slo = 3.046 gov = 0.0d0 gamsw = 0.8d0 gamup = 0.0d0 barg = 0.5d0 xlab = 0.3d0 turb = 0.0d0 tx = 1.0d0 tt = 1.0d0 itt = 0 if(itergov.eq.1) gov = 0.000 if(itergov.eq.2) gov = 0.025 if(itergov.eq.3) gov = 0.05 if(iturb.eq.1) turb = 0.0 if(iturb.eq.2) turb = 0.05 if(iturb.eq.3) turb = 0.1 if(iturb.eq.4) turb = 0.15 if(iturb.eq.5) turb = 0.2 if(iturb.eq.6) turb = 0.25 if(iturb.eq.7) turb = 0.3 if(iturb.eq.8) turb = 0.35 if(iturb.eq.9) turb = 0.4 if(iturb.eq.10) turb = 0.45 if(iturb.eq.11) turb = 1.0 if(iphi.eq.1) phi = 0.3 if(iphi.eq.2) phi = 0.5 if(iterbb.eq.1) bbbb = 0.97 if(iterbb.eq.2) bbbb = 0.98 if(iterbb.eq.3) bbbb = 0.99 c c construct the grid points for determining the functional form of g_k_(z) c do 10 i = 1,ngrid grid(i,1)=ulo - 0.5d0*slo + dble(i-1)*slo/dble(ngrid-1) grid(i,2)=uhi - 0.5d0*shi + dble(i-1)*shi/dble(ngrid-1) 10 continue do 15 j = 1,npar do 15 k = 1,npar xxlo(j,k) = 0.d0 xxhi(j,k) = 0.d0 do 15 i = 1,ngrid xxlo(j,k) = xxlo(j,k) + grid(i,1)**dble(j-1)*grid(i,1)**dble(k-1) xxhi(j,k) = xxhi(j,k) + grid(i,2)**dble(j-1)*grid(i,2)**dble(k-1) 15 continue call invert(xxlo,npar,det) call invert(xxhi,npar,det) c c c iterate to solve for the coefficients of the polynomial, the unemployment masses, c the "w"-values, and the "b"-values c do 1448 iter2 = 1,25000 accur = 0.00000 w0old = w0 wloloold = wlolo wlohiold = wlohi whihiold = whihi bloold = blo bhiold = bhi c c calculate ghi and glo as a function of z c out = tt*blo + wlolo - coeflo(1) if(out.ge.0.d0) then zacceptlolo = out/(tx+coeflo(2)) else zacceptlolo = out/(1.d0+coeflo(2)) endif out = tt*bhi + wlohi - coeflo(1) if(out.ge.0.d0) then zacceptlohi = out/(tx+coeflo(2)) else zacceptlohi = out/(1.d0+coeflo(2)) endif out = tt*bhi + (1.d0-turb)*whihi+ turb*wlohi - coefhi(1) if(out.ge.0.d0) then zaccepthihi = out/(tx+coefhi(2)) else zaccepthihi = out/(1.d0+coefhi(2)) endif zswitchlo = zacceptlolo zswitchhi = zaccepthihi out = w0 - bbbb*(tt*blo+wlolo) if(out.ge.0.d0) then temp1 = out/tx else temp1 = out endif out = w0- coeflo(1) if(out.ge.0.d0) then temp2 = out/(tx+coeflo(2)) else temp2 = out/(1.d0+coeflo(2)) endif if(temp1.lt.temp2) then zaccept0=temp1 else zaccept0=temp2 endif c c note that there are two continuation values for the guys with c no unemployment benefits who accept a job. The ones with a good shock will always stay c in a relationship. The ones with a bad shock might breakup the c next period (these guys only accept the job to become entitled c to unemployment benefits). c t1 = zaccept0 t2 = zacceptlolo t3 = zacceptlohi t4 = zaccepthihi t5 = zswitchlo t6 = zswitchhi if(zaccept0.gt.ulo+0.5*slo) then write(*,*) 'darn' stop endif if(zacceptlolo.gt.ulo+0.5*slo) then write(*,*) 'darn' stop endif if(zacceptlohi.gt.ulo+0.5*slo) then write(*,*) 'darn' stop endif if(zaccepthihi.gt.uhi+0.5*shi) then write(*,*) 'darn' stop endif if(zswitchlo.gt.ulo+0.5*slo) then write(*,*) 'darn' stop endif if(zswitchhi.gt.uhi+0.5*shi) then write(*,*) 'darn' stop endif raccept0 = fnormcum(zaccept0,ulo,slo) racceptlolo = fnormcum(zacceptlolo,ulo,slo) racceptlohi = fnormcum(zacceptlohi,ulo,slo) raccepthihi = fnormcum(zaccepthihi,uhi,shi) rswitchlo = racceptlolo rswitchhi = fnormcum(zswitchhi,uhi,shi) if(zaccept0 .lt.ulo-0.5d0*slo) zaccept0 = ulo-0.5d0*slo if(zacceptlolo.lt.ulo-0.5d0*slo) zacceptlolo = ulo-0.5d0*slo if(zacceptlohi.lt.ulo-0.5d0*slo) zacceptlohi = ulo-0.5d0*slo if(zaccepthihi.lt.uhi-0.5d0*shi) zaccepthihi = uhi-0.5d0*shi if(zswitchlo .lt.ulo-0.5d0*slo) zswitchlo = ulo-0.5d0*slo if(zswitchhi .lt.uhi-0.5d0*shi) zswitchhi = uhi-0.5d0*shi eaccept0 = 0.d0 yy = (ulo+0.5*slo-zaccept0)/dble(nsimpson-1) do 110 i = 1,nsimpson qx = dble(i-1)*yy + zaccept0 if(qx.lt.zacceptlolo) then temp0 = ftx(qx)*qx + bbbb*(tt*blo+wlolo) else temp0 = ftx(qx)*qx + g(qx,coeflo,npar) endif qw = hsw(i)*yy*temp0*fdens(slo) eaccept0 = eaccept0 + qw 110 continue c c note the two different continuation values used in 110 c eacceptlolo = 0.d0 yy = (ulo+0.5*slo-zacceptlolo)/dble(nsimpson-1) do 120 i = 1,nsimpson qx = dble(i-1)*yy + zacceptlolo temp0 = ftx(qx)*qx + g(qx,coeflo,npar) qw = hsw(i)*yy*temp0*fdens(slo) eacceptlolo = eacceptlolo + qw 120 continue eacceptlohi = 0.d0 yy = (ulo+0.5*slo-zacceptlohi)/dble(nsimpson-1) do 125 i = 1,nsimpson qx = dble(i-1)*yy + zacceptlohi temp0 = ftx(qx)*qx + g(qx,coeflo,npar) qw = hsw(i)*yy*temp0*fdens(slo) eacceptlohi = eacceptlohi + qw 125 continue eaccepthihi = 0.d0 yy = (uhi+0.5*shi-zaccepthihi)/dble(nsimpson-1) do 130 i = 1,nsimpson qx = dble(i-1)*yy + zaccepthihi temp0 = ftx(qx)*qx + g(qx,coefhi,npar) qw = hsw(i)*yy*temp0*fdens(shi) eaccepthihi = eaccepthihi + qw 130 continue eswitchlo = eacceptlolo eswitchhi = 0.d0 yy = (uhi+0.5*shi-zswitchhi)/dble(nsimpson-1) do 135 i = 1,nsimpson qx = dble(i-1)*yy + zswitchhi temp0 = ftx(qx)*qx + g(qx,coefhi,npar) qw = hsw(i)*yy*temp0*fdens(shi) eswitchhi = eswitchhi + qw 135 continue do 100 igrid = 1,ngrid zlo = grid(igrid,1) zhi = grid(igrid,2) c c calculate ghi and glo as a function of z c ylo(igrid) = bbbb*(1.d0-roexlo) $ * ( $ ( gamup) *(eswitchhi/(1.d0-rswitchhi)) $ +(1.d0-gamup)*(1.d0-gamsw)*(ftx(zlo)*zlo+g(zlo,coeflo,npar)) $ +(1.d0-gamup)*( gamsw)*eswitchlo $ ) $ + bbbb*( $ roexlo * (tt*blo+wlolo) $ +(1.d0-roexlo)*(1.d0-gamup)*gamsw*rswitchlo*(tt*blo+wlolo) $ ) c c note that we assume that an upgraded worker gets a draw above zswitchhi c given that he is still only entitled to blo this will certainly guarantee c that he doesn't break up c yhi(igrid) = bbbb*(1.d0-roexhi) $ * ( $ (1.d0-gamsw)*(ftx(zhi)*zhi+g(zhi,coefhi,npar)) $ +( gamsw)*eswitchhi $ ) $ + bbbb $ * ( $ roexhi * turb *(tt*bhi+wlohi) $ + roexhi *(1.d0-turb)*(tt*bhi+whihi) $ +(1.d0-roexhi)* turb *(tt*bhi+wlohi)*gamsw*rswitchhi $ +(1.d0-roexhi)*(1.d0-turb)*(tt*bhi+whihi)*gamsw*rswitchhi $ ) 100 continue c c calculate wlolo and whihi c out = (1.d0-turb)*(tt*bhi+whihi) +turb*(tt*bhi+wlohi) outwo = (1.d0-turb)*(tt*bhi+wwohihi)+turb*(tt*bhi+wwolohi) temp0 = xlab*bbbb*barg temp1 = temp0*(eaccept0 -(1.d0-raccept0) * w0 ) $ + bbbb*wwo0 temp2 = temp0*(eacceptlolo-(1.d0-racceptlolo)*(tt*blo+wlolo)) $ + bbbb*(tt*blo+wwololo) temp3 = temp0*(eacceptlohi-(1.d0-racceptlohi)*(tt*bhi+wlohi)) $ + bbbb*(tt*bhi+wwolohi) temp4 = temp0*(eaccepthihi-(1.d0-raccepthihi)* out ) $ + bbbb*outwo temp0 =xlab*bbbb*(1.d0-barg) temp5 $ = (u0 /(u0+ulolo+ulohi+uhihi)) * temp0 $ * (eaccept0 -(1.d0-raccept0) * w0 ) $ + (ulolo/(u0+ulolo+ulohi+uhihi)) * temp0 $ * (eacceptlolo-(1.d0-racceptlolo)*(tt*blo+wlolo)) $ + (ulohi/(u0+ulolo+ulohi+uhihi)) *temp0 $ * (eacceptlohi-(1.d0-racceptlohi)*(tt*bhi+wlohi)) $ + (uhihi/(u0+ulolo+ulohi+uhihi)) * temp0 $ * (eaccepthihi-(1.d0-raccepthihi)* out ) $ + bbbb*wfi accur=accur+dabs(temp1+temp5-w0old)+dabs(temp2+temp5-wloloold) $ + dabs(temp3+temp5-wlohiold) + dabs(temp4+temp5-whihiold) wwo0 = temp1*adj2 + (1.d0-adj2)*wwo0 wwololo = temp2*adj2 + (1.d0-adj2)*wwololo wwolohi = temp3*adj2 + (1.d0-adj2)*wwolohi wwohihi = temp4*adj2 + (1.d0-adj2)*wwohihi wfi = temp5*adj2 + (1.d0-adj2)*wfi w0 = wwo0 + wfi wlolo = wwololo + wfi wlohi = wwolohi + wfi whihi = wwohihi + wfi c c run the regressions c do 150 j = 1,npar xylo(j,1) = 0.d0 xyhi(j,1) = 0.d0 do 150 i = 1,ngrid xylo(j,1) = xylo(j,1) + grid(i,1)**dble(j-1)*ylo(i) xyhi(j,1) = xyhi(j,1) + grid(i,2)**dble(j-1)*yhi(i) 150 continue call mult(xxlo,xylo,coef,npar,npar,1) do 152 j = 1,npar accur = accur + dabs(coeflo(j)-coef(j,1)) 152 continue do 151 j = 1,npar coeflo(j) = adj1*coef(j,1)+(1.d0-adj1)*coeflo(j) 151 continue call mult(xxhi,xyhi,coef,npar,npar,1) do 153 j = 1,npar accur = accur + dabs(coefhi(j)-coef(j,1)) 153 continue do 156 j = 1,npar coefhi(j) = adj1*coef(j,1)+(1.d0-adj1)*coefhi(j) 156 continue c c update the unemployment masses and the w_k_(b) values c note that each outflow shows up as an inflow somewhere c u0n = u0 $ +roreti $ -( roreti)*u0 $ -(1.d0-roreti)*u0 *( xlab)*(1.d0-roexlo)*(1.d0-raccept0) ulolon=ulolo $ +(1.d0-roreti)*elotemp0 $ +(1.d0-roreti)*( roexlo)*elo $ +(1.d0-roreti)*(1.d0-roexlo)*elo*(1.d0-gamup)*gamsw*rswitchlo $ -( roreti)*ulolo $ -(1.d0-roreti)*ulolo*( xlab)*(1.d0-roexlo)*(1.d0-racceptlolo) $ +(1.d0-roreti)*( roexlo)*elohi $ +(1.d0-roreti)*(1.d0-roexlo)*elohi*(1.d0-gamup)*gamsw $ *(racceptlolo) ulohin=ulohi $ +(1.d0-roreti)*( roexhi)*ehi*turb $ +(1.d0-roreti)*(1.d0-roexhi)*ehi*gamsw*rswitchhi*turb $ -( roreti)*ulohi $ -(1.d0-roreti)*ulohi*( xlab)*(1.d0-roexlo)*(1.d0-racceptlohi) uhihin=uhihi $ +(1.d0-roreti)*( roexhi)*ehi*(1.d0-turb) $ +(1.d0-roreti)*(1.d0-roexhi)*ehi*gamsw*rswitchhi*(1.d0-turb) $ -( roreti)*uhihi $ -(1.d0-roreti)*uhihi*( xlab)*(1.d0-roexhi)*(1.d0-raccepthihi) elotemp0n = elotemp0 $ -( roreti)*elotemp0 $ -(1.d0-roreti)*elotemp0 $ +(1.d0-roreti)*u0*xlab*(1.d0-roexlo)*(racceptlolo-raccept0) elon=elo $ +(1.d0-roreti)*u0 *( xlab)*(1.d0-roexlo)*(1.d0-racceptlolo) $ +(1.d0-roreti)*ulolo*( xlab)*(1.d0-roexlo)*(1.d0-racceptlolo) $ -(1.d0-roreti)*(1.d0-roexlo)*elo*(1.d0-gamup)*gamsw*rswitchlo $ -(1.d0-roreti)*(1.d0-roexlo)*elo*gamup $ -( roreti)*elo $ -(1.d0-roreti)*( roexlo)*elo $ +(1.d0-roreti)*(1.d0-roexlo)*elohi*(1.d0-gamup)*gamsw $ *(1.d0-racceptlolo) elohin = elohi $ +(1.d0-roreti)*elotemphi $ -( roreti)*elohi $ -(1.d0-roreti)*(1.d0-roexlo)*elohi*(1.d0-gamup)*gamsw $ -(1.d0-roreti)*(1.d0-roexlo)*elohi*( gamup) $ -(1.d0-roreti)*( roexlo)*elohi elotemphin = elotemphi $ -( roreti)*elotemphi $ -(1.d0-roreti)*elotemphi $ +(1.d0-roreti)*ulohi*( xlab)*(1.d0-roexlo)*(1.d0-racceptlohi) ehin=ehi $ +(1.d0-roreti)*uhihi*( xlab)*(1.d0-roexhi)*(1.d0-raccepthihi) $ -( roreti)*ehi $ -(1.d0-roreti)*( roexhi)*ehi $ -(1.d0-roreti)*(1.d0-roexhi)*ehi*gamsw*rswitchhi $ +(1.d0-roreti)*ehitempup ehitempupn = ehitempup $ -( roreti)*ehitempup $ -(1.d0-roreti)*ehitempup $ +(1.d0-roreti)*(1.d0-roexlo)*elo*gamup $ +(1.d0-roreti)*(1.d0-roexlo)*elohi*( gamup) accur = accur+dabs(u0-u0n)+dabs(ulolo-ulolon)+dabs(ulohi-ulohin) $ + dabs(ulohi-ulohin)+dabs(elo-elon) +dabs(ehi-ehin) $ + dabs(elotemp0-elotemp0n)+dabs(ehitempup-ehitempupn) $ + dabs(elotemphi-elotemphin) u0 = u0n ulolo = ulolon ulohi = ulohin uhihi = uhihin elo = elon ehi = ehin elotemp0 = elotemp0n ehitempup = ehitempupn elotemphi = elotemphin elohi = elohin c c calcualate average wages for elotemp0,elotemphi,elo, and ehi. Note c despite the updating processes the wages are uniformly distributed c within [zaccept*,u*+s*/2]. Note that this is the point where c heterogeneity really becomes painful. Everything above could have c been done by including elotemp0 & elotemphi in elo. And this is c the reason why we have made assumptions to keep heterogeneity low c (like the assumption that an unemployed with high skills who gets a job offer c faces the same risk of skill-loss as a high skilled worker who receives c a new draw and the assumption that an upgrade implies getting c a shock in the region [zaccept*,u*+s*/2] c c j=1: e0 c j=2: elo c j=3: elotemphi c j=4: elohi c j=5: ehi c j=6: ehitempup do 210 j = 1,6 if(j.eq.1.and.dabs(racceptlolo-raccept0).lt.0.001) then wage(j) = 0.d0 goto 210 else zaccept = zaccept0 raccept = raccept0+1.d0-racceptlolo zupper = zacceptlolo out = w0 outwo = wwo0 endif if(j.eq.2) then zaccept = zacceptlolo raccept = racceptlolo zupper = ulo+0.5d0*slo out = tt*blo+wlolo outwo = tt*blo+wwololo endif if(j.eq.3) then zaccept = zacceptlohi raccept = racceptlohi zupper = ulo+0.5d0*slo out = tt*bhi+wlohi outwo = tt*bhi+wwolohi endif if(j.eq.4) then zaccept = zacceptlohi raccept = racceptlohi zupper = ulo+0.5d0*slo out = tt*blo+wlolo outwo = tt*blo+wwololo endif if(j.eq.5) then zaccept = zaccepthihi raccept = raccepthihi zupper = uhi+0.5*shi out = (1.d0-turb)*(tt*bhi+whihi) +turb*(tt*bhi+wlohi) outwo= (1.d0-turb)*(tt*bhi+wwohihi)+turb*(tt*bhi+wwolohi) endif if(j.eq.6) then zaccept = zaccepthihi raccept = raccepthihi zupper = uhi+0.5*shi out = tt*blo+wlolo outwo= tt*blo+wwololo endif if(j.le.4) then temps = slo do 211 k = 1,npar coef(k,1) = coeflo(k) 211 continue else temps = shi do 212 k = 1,npar coef(k,1) = coefhi(k) 212 continue endif wage(j) = 0.0d0 tax(j) = 0.0d0 yy = (zupper-zaccept)/dble(nsimpson-1) do 220 i = 1,nsimpson qx = dble(i-1)*yy + zaccept temp0 = barg*g(qx,coef,npar)-bbbb*barg*out+bbbb*outwo temp0 = barg*(tx*qx + g(qx,coef,npar) - out) $ +outwo-temp0 qw = hsw(i)*yy*fdens(temps)*temp0 wage(j) = wage(j) + qw tax(j) = tax(j) + hsw(i)*yy*fdens(temps)*ftx(qx)*qx 220 continue wage(j) = wage(j)/(1.d0-raccept) tax(j) = tax(j) /(1.d0-raccept) 210 continue blo = (elotemp0*wage(1)+elo*wage(2)+elotemphi*wage(3) $ + elohi*wage(4))*phi/ $ (elotemp0+elo+elotemphi+elohi) bhi = (ehi*wage(5)+ehitempup*wage(6))*phi/ $ (ehi+ehitempup) temp1 = ulolo*tt*blo+ulohi*tt*bhi+uhihi*tt*bhi + gov temp2 = elotemp0*tax(1)+elo*tax(2)+elotemphi*tax(3) $ + elohi*tax(4)+ehi*tax(5)+ehitempup*tax(6) txn = 1.d0 - temp1/temp2 tx = adj3*txn+(1.d0-adj3)*tx if(itt.eq.1) tt = tx blo = blo + bfix*ulo bhi = bhi + bfix*uhi accur = accur + dabs(blo-bloold) + dabs(bhi-bhiold) $ + dabs(tx-txn) blo = adj3*blo+(1.d0-adj3)*bloold bhi = adj3*bhi+(1.d0-adj3)*bhiold c write(*,'(8f10.5)') zacceptlolo,wlolo,coeflo(1),coeflo(2),blo, c $ulolo,elo,wfi c write(*,*) accur if(accur.lt.accurcrit) goto 1449 1448 continue write(*,*) 'not converged' 1449 continue if(dabs(u0+ulolo+ulohi+uhihi+elotemphi+elohi $ +elo+ehi+elotemp0+ehitempup-1.d0).gt.0.0001) then write(*,*) 'THERE IS AN ERROR IN THE PROGRAM YOU STUPID' endif if(itergov.eq.1) baal = t6 write(*,'(4f8.4,i6,4f8.4)') (u0+ulolo+ulohi+uhihi)*100.d0,roexhi, $shi,accur,iter3,t6,uhi-0.5d0*shi write(33,*) write(33,500) turb,bbbb,accur write(33,501) ulo,uhi,slo,shi write(33,502) barg,roreti*100.,roexlo*100.,roexhi*100. write(33,503) phi,gamsw,xlab,gamup write(33,*) write(33,601) blo,bhi write(33,600) u0*100,ulolo*100,ulohi*100,uhihi*100, $u0*100+ulolo*100+ulohi*100+uhihi*100 write(33,603) $elotemp0*100,elo*100,elotemphi*100,elohi*100 write(33,605)ehi*100,ehitempup*100 write(33,604) (wage(j),j=1,4) write(33,606) (wage(j),j=5,6),(1.d0-tx)*100.d0,(1.d0-tt)*100.d0, $gov write(33,*) write(33,701) raccept0*100, t1, ulo-0.5d0*slo write(33,702) racceptlolo*100,t2,ulo-0.5d0*slo write(33,703) racceptlohi*100,t3,ulo-0.5d0*slo write(33,704) raccepthihi*100,t4,uhi-0.5d0*shi write(33,705) rswitchlo*100,t5,ulo-0.5d0*slo write(33,706) rswitchhi*100,t6,uhi-0.5d0*shi write(33,*) 448 continue write(32,*) (coeflo(j),j=1,npar) write(32,*) (coefhi(j),j=1,npar) write(32,*) blo,bhi,w0,wlolo,wlohi,whihi,wfi, $wwo0,wwololo,wwolohi,wwohihi, $zswitchlo,zswitchhi, $zaccept0,zacceptlolo,zacceptlohi,zaccepthihi, $u0,ulolo,ulohi,uhihi,elo,ehi,elotemp0,ehitempup,elohi 500 format(' turb =',f7.3,' bbbb =',f7.3, $ ' accur =',f7.5) 501 format(' ulo =',f7.3,' uhi =',f7.3,' slo =',f7.3, $ ' shi =',f7.3) 502 format(' barg =',f7.3,' roreti =',f7.3,' roexlo =',f7.3, $ ' roexhi =',f7.3) 503 format(' phi =',f7.3, $ ' gamsw =',f7.3, ' xlab =',f7.3,' gamup =',f8.4) 600 format(' u0 =',f6.2,' ulolo =',f6.2,' ulohi =',f6.2, $ ' uhihi =',f6.2, ' usum =',f6.2) 601 format(' blo =',f6.3,' bhi =',f6.3) 603 format(' elotemp0 =',f6.2,' elo =',f6.2,' elotemphi =',f6.2, $ ' elohi =',f6.2) 605 format(' ehi =',f6.2,' ehitup =',f6.2) 604 format(' wlotemp0 =',f6.3,' wlo =',f6.3,' wlotemphi =',f6.3, $ ' wlohi =',f6.3) 606 format(' whi =',f6.3,' whitup =',f6.3,' txrate =',f6.3, $ ' ttrate =',f6.3,' gov =',f6.3) 701 format(' prob accept 0 =',3f12.4) 702 format(' prob accept lolo =',3f12.4) 703 format(' prob accept lohi (if turb > 0 ) =',3f12.4) 704 format(' prob accept hihi (if turb < 1 ) =',3f12.4) 705 format(' prob destruct lo =',3f12.4) 706 format(' prob destruct hi =',3f12.4) 707 format(' accuracy =',3f10.4) stop end *************************************************************************** * * polynomials * *************************************************************************** real*8 function g(z,coef,np) integer np real*8 temp,z real*8 coef(np) g = coef(1) temp = 1.d0 do 10 i = 2,np temp = temp*z g = g + coef(i)*temp 10 continue return end *************************************************************************** * * FDENS(); Returns the density of a uniform * *************************************************************************** real*8 function fdens(s) real*8 s real*8 uhi,ulo,shi,slo,pi,bbbb,phi,barg,xlab, $roexhi,roexlo,roreti,gamup,gamsw,tx common/par/ uhi,ulo,shi,slo,pi,bbbb,phi,barg,xlab, $roexhi,roexlo,roreti,gamup,gamsw,tx fdens = 1.d0/s return end *************************************************************************** * * FNORMCUM(); CDF of a uniform * *************************************************************************** real*8 function fnormcum(x,u,s) parameter(ngrid=41,npar=3,nsimpson=11) real*8 x,u,s,fdens real*8 uhi,ulo,shi,slo,pi,bbbb,phi,barg,xlab, $roexhi,roexlo,roreti,gamup,gamsw,tx real*8 hsw(nsimpson) common/par/ uhi,ulo,shi,slo,pi,bbbb,phi,barg,xlab, $roexhi,roexlo,roreti,gamup,gamsw,tx common/simpson/ hsw external fdens if(x.le.u-0.5*s) then fnormcum = 0.d0 goto 2000 endif if(x.ge.u+0.5*s) then fnormcum = 1.d0 goto 2000 endif fnormcum = (x-(u-0.5d0*s))/s 2000 continue return end *************************************************************************** * * FTX(); Returns the tax rate * *************************************************************************** real*8 function ftx(s) real*8 s real*8 uhi,ulo,shi,slo,pi,bbbb,phi,barg,xlab, $roexhi,roexlo,roreti,gamup,gamsw,tx common/par/ uhi,ulo,shi,slo,pi,bbbb,phi,barg,xlab, $roexhi,roexlo,roreti,gamup,gamsw,tx if(s.ge.0.d0) then ftx = tx else ftx = 1.d0 endif return end *************************************************************************** * SUBROUTINE MULT(E , F , G , NROWE , NROWF , NCOLF) * * * THIS SUBROUTINE CALCULATES THE PRODUCT G OF TWO MATRICES * E AND F. * *************************************************************************** IMPLICIT REAL*8(A-H,O-Z) REAL*8 E(NROWE,NROWF) , F(NROWF,NCOLF) , G(NROWE,NCOLF) DO 831 J = 1,NCOLF DO 832 I = 1,NROWE G(I,J) = 0.0 832 CONTINUE 831 CONTINUE DO 833 J = 1,NCOLF DO 834 I = 1,NROWE DO 835 K = 1,NROWF G(I,J) = E(I,K)*F(K,J) + G(I,J) 835 CONTINUE 834 CONTINUE 833 CONTINUE RETURN END *************************************************************************** * SUBROUTINE INVERT(A,N,D) * * * THIS SUBROUTINE COMPUTES THE INVERSE OF A MATRIX. * * *************************************************************************** IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(N,N) ,L(100) , M(100) COMMON // ISING D=1.D0 DO 80 K=1,N L(K)=K M(K)=K BIGA=A(K,K) DO 20 I=K,N DO 20 J=K,N IF(DABS(BIGA)-DABS(A(I,J))) 10,20,20 10 BIGA=A(I,J) L(K)=I M(K)=J 20 CONTINUE IF (DABS(BIGA).LT.(1.0E-15)) GO TO 99 J=L(K) IF(L(K)-K) 35,35,25 25 DO 30 I=1,N HOLD=-A(K,I) A(K,I)=A(J,I) 30 A(J,I)=HOLD 35 I=M(K) IF(M(K)-K) 45,45,37 37 DO 40 J=1,N HOLD=-A(J,K) A(J,K)=A(J,I) 40 A(J,I)=HOLD 45 DO 55 I=1,N IF(I-K) 50,55,50 50 A(I,K)=A(I,K)/(-A(K,K)) 55 CONTINUE DO 65 I=1,N DO 65 J=1,N IF(I-K) 57,65,57 57 IF(J-K) 60,65,60 60 A(I,J)=A(I,K)*A(K,J)+A(I,J) 65 CONTINUE DO 75 J=1,N IF(J-K) 70,75,70 70 A(K,J)=A(K,J)/A(K,K) 75 CONTINUE D=D*A(K,K) A(K,K)=1.D0/A(K,K) 80 CONTINUE K=N 100 K=K-1 IF(K) 150,150,103 103 I=L(K) IF(I-K) 120,120,105 105 DO 110 J=1,N HOLD=A(J,K) A(J,K)=-A(J,I) 110 A(J,I)=HOLD 120 J=M(K) IF(J-K) 100,100,125 125 DO 130 I=1,N HOLD=A(K,I) A(K,I)=-A(J,I) 130 A(J,I)=HOLD GO TO 100 99 CONTINUE ISING = 1 WRITE(11,991) ISING 991 FORMAT(' SINGULAR MATRIX = = > ISING = ',I4) 150 CONTINUE RETURN END *************************************************************************** * * sample for euro1uni.old * *************************************************************************** c 47.139618210339200 2.383991833942383E-001 -3.256097917741678E-004 c 47.139667194739790 2.383991833942399E-001 -3.256097917765428E-004 c 3.041951162492100E-001 3.042091926714169E-001 44.638611679559160 c 46.198278265761550 46.198348685559160 46.198391784114980 c 17.611621083910490 27.026990595648670 28.586657181851060 c 28.586727601648670 28.586770700204490 -5.144906772046096E-001 c -5.144271998388031E-001 -5.230000019073486E-001 -5.144906772046096E-001 c -5.144224470582484E-001 -5.144271998388031E-001 1.664637196250435E-002 c 3.988235245985401E-002 0.000000000000000E+000 3.801783113090964E-084 c 9.434575400633458E-001 8.841851124386825E-083 1.373551429433376E-005 c 0.000000000000000E+000 0.000000000000000E+000