	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
