c			  C:\DRW2\SIMPLE\IMP2.f
c
c
c               Program to calculate:
c
c			  (a) impulse response function for one-time shock in the amount
c				  of exogenousbreakups
c
c               This program: 
c
c               - uses simple two point distribution {0,hu}
c               - the fraction that receives hu, gamh, depends on H
c               - to calculate the impulse response functions, 
c				it is assumed that H never drops so much that 
c                 agents who receive hu would like to quit
c                 this means that the endogenous breakup probability
c                 is equal to 1-gamh
c
c               September 17 1998
c

		use msimsl
		implicit real*8(a-h,o-z),integer(i-n)
		common/simpson/hsw(1001),nsimpson
		common/par/xvalue,alpha,prob,depr,bbbb,xi,rox,gam,eta,
     $          hstst,vcost,ststg,ststw,ststge,ubar
		common/var/qn(1000),qh(1000),qg(1000),qw(1000),
     $          qge(1000),qlabx(1000),qv(1000),qubar(1000),
     $          qy(1000),qa(1000),qnn(1000),qhu(1000),qgam(1000)                
		common/iiii/iend,ifunc

		dimension xguess(3),xout(3)

		external fncimp,fnccaph1,fgam,fgaminv


		open(99,file='for099.dat',
     $    status='unknown')


			alpha  = 0.3
			prob   = 0.5
			depr   = 0.0
			bbbb   = 0.96
			xi	   = 0.25


c
c			x = 3, ifunc = 2
c
c 			gam    = 0.069777
c			eta    = 32.5
c
c			gam    = 0.059246
c			eta    = 25.
c
c			gam    = 0.053831
c			eta    = 20.
c
c			gam    = 0.065874
c			eta    = 30.
c
c			gam    = 0.042246
c			eta    = 5.
c			iii    = 11
c			rox    = 0.999
c
c			x = 3, ifunc = 3
c			gam    = 0.323145
c			eta    = 0.25
c			iii    = 22
c			rox    = 0.99
c
c
c			gam     = 0.400925
c			eta     = 0.2
c			iii     = 15
c			rox     = 0.99
c			
c			when x = 3 then hstst  = 74.6980935014833
c
c
c			x = 1 values, ifunc = 2
c
c			gam    = 0.075705
c			eta    = 40.
c			gam    = 0.086873
c			eta    = 45.
c			gam    = 0.080809
c			eta    = 42.5
c			hmean  = 78.89381


			ifunc   = 3
			gam     = 0.400925
			eta     = 0.2
			
			hstst   = 74.6980935014833
			xvalue  = 3.
			zero    = 0.

			ststgam =  fgam(hstst)
			ubar    =  1-ststgam
			ststn   =  prob/(prob+ubar-prob*ubar)
			ststv   =  (1.-ststn+ubar*ststn)*(prob/xi)**2.
			ststge  =  ststgam*bbbb*xvalue
			vcost   =  xi*ststge*((1.-ststn+ubar*ststn)/ststv)**0.5             

			ststm   =  xvalue-ststge
			ststhu  =  hstst/ststgam
			ststy   =  ststgam*((1.-depr)*ststhu)**alpha

			cw = -prob*ststge - (1.-prob)*bbbb*depr*hstst
			cw =  cw /(1.-(1.-prob)*bbbb)
			sw =  prob/(1.-(1.-prob)*bbbb)
			ststg = bbbb*ststy-bbbb*depr*hstst
     $        + (1.-ststgam)*bbbb*cw
			ststg = ststg/(1.-bbbb*ststgam-(1.-ststgam)*bbbb*sw)
			ststw = cw + sw*ststg

			write(*,'(4f14.6)') ststn,ststv,vcost,ststge
			write(*,'(5f14.6)') ststg,ststw,ststhu,ststy,ststgam

			pause							  

			iii = 15
c			iii = 11

			rox   = 0.99
			qn(1) = ststn


c
c
c			solve for "fixed H" impulse response functions
c
c


		do 190 ii = 1,iii+1
		qa(ii) = 1.
190		continue


			do 200 ik = 1,iii+1
			qh(ik)     = hstst
			qgam(ik)   = fgam(hstst)
			qubar(ik)  = 1.-qgam(ik)
			qhu(ik)    = hstst/qgam(ik)
			qy(ik)     = qa(ik)*qgam(ik)*((1.-depr)*qhu(ik))**alpha
200			continue

			qnn(1) = qn(1)*rox

			do 210 ik = 1,iii
			qge(ik)    = bbbb*qgam(ik+1)*xvalue
			qlabx(ik)  = xi*xi*qge(ik)/vcost
			qn(ik+1)    = qgam(ik)*qnn(ik) 
     $                   + qlabx(ik)*(1.-qgam(ik)*qnn(ik))
			qnn(ik+1)  = qn(ik+1)
210			continue

			
			do 220 i = 1,iii
			write(*,*) i
			write(99,'(i6,8f14.8)') i,qa(i),qubar(i),qn(i),qh(i),
     $		qy(i),qy(i)*qnn(i)/(qy(iii+1)*qnn(iii+1))
			write(*,'(i6,8f14.8)') i,qa(i),qubar(i),qn(i),qh(i),
     $		qy(i),qy(i)*qnn(i)/(qy(iii+1)*qnn(iii+1))
220			continue
		
			pause





c
c
c               solve for impulse response functions
c
c
		npar  = 3
		itmax = 1000            
		xguess(1) = ststg
		xguess(2) = ststw
		xguess(3) = ststge	 
     	

		do 5555 iend = 4,iii

		write(*,*) iend
		pause


		qn(1) = ststn

c
c		rox only hits qn(1) in period 1 before qh(1) is chosen
c


		errrel = 0.0000001		

		call dneqnf(fncimp,errrel,npar,itmax,xguess,xout,fnorm)    
		
		xguess(1) = xout(1)
		xguess(2) = xout(2)
		xguess(3) = xout(3)
5555		continue

c
c               actually print out impulse response functions 
c

		do 400 i = 1,iii
		write(99,'(i6,8f14.8)') i,qa(i),qubar(i),qn(i),qh(i),
     $    qy(i),qy(i)*qnn(i)/(qy(iii+1)*qnn(iii+1))
		write(*,'(i6,8f14.8)') i,qa(i),qubar(i),qn(i),qh(i),
     $    qy(i),qy(i)*qnn(i)/(qy(iii+1)*qnn(iii+1))
400       continue
		pause

		do 401 i = 1,iii
		write(*,'(i6,7f10.6)') i,qg(i),qw(i),qge(i),qlabx(i),qv(i)
		write(99,'(i6,7f14.8)') i,qg(i),qw(i),qge(i),qlabx(i),qv(i)
401       continue


		stop
		end


		subroutine fncimp(xin,fout,npar)
		implicit real*8(a-h,o-z),integer(i-n)
		common/simpson/hsw(1001),nsimpson
		common/par/xvalue,alpha,prob,depr,bbbb,xi,rox,gam,eta,
     $          hstst,vcost,ststg,ststw,ststge,ubar
		common/var/qn(1000),qh(1000),qg(1000),qw(1000),
     $          qge(1000),qlabx(1000),qv(1000),qubar(1000),
     $          qy(1000),qa(1000),qnn(1000),qhu(1000),qgam(1000)                
		common/iiii/iend,ifunc


		dimension fout(npar),xin(npar)

		external fnccaph1,fgam,fgaminv

				
		qg(1)  = xin(1)
		qw(1)  = xin(2)
		qge(1) = xin(3)


c		write(*,'(i5,3f20.15)') iend,xin(1),xin(2),xin(3)
			
c
c
c		solve for capital H
c
c
		ain = 0.80*hstst
		bin = 1.05*hstst

		errabs = 0.0
		errrel = 0.00000000001
		maxfn  = 1000


		call dzbren(fnccaph1,errabs,errrel,ain,bin,maxfn)

		qh(1) = bin			
		qubar(1) = ubar
		qgam(1)  = 1.-qubar(1)
		qhu(1)   = qh(1)/qgam(1)

		
		qy(1)	= qa(1)*qgam(1)*((1.-depr)*qhu(1))**alpha
						
		
c
c		define number of relationships after exogenous breakups
c		
	
		qnn(1) = rox*qn(1)

c		write(*,'(i6,5f14.6)') 1,qg(1),qw(1),qge(1),qhu(1),qh(1)				
		
		do 200 i = 1,iend
c		write(*,*) iend,i		
		unemp      = 1.-qnn(i)+qnn(i)*qubar(i)
		qv(i)      = unemp*(xi*qge(i)/vcost)**2.
		qlabx(i)   = xi*(qv(i)/unemp)**0.5
		qn(i+1)    = (1.-qubar(i))*qnn(i) + xi*unemp**0.5*qv(i)**0.5
		qubar(i+1) = (bbbb*xvalue-qge(i))/(bbbb*xvalue)		

		
		qgam(i+1)  = 1.-qubar(i+1)
		temp       = qgam(i+1)
		qh(i+1)    = fgaminv(temp)
		
		qhu(i+1)   = qh(i+1)/qgam(i+1)

		qw(i+1)    = qw(i) - qlabx(i)*(qg(i)-qge(i))
     $               + (1.-qlabx(i))*bbbb*depr*qh(i+1)
		qw(i+1)    = qw(i+1)/((1.-qlabx(i))*bbbb)

		qy(i+1)    = qa(i+1)*qgam(i+1)*((1.-depr)*qhu(i+1))**alpha

		qg(i+1)    = qg(i)-bbbb*qy(i+1)+bbbb*depr*qh(i+1)
     $               - qubar(i+1)*bbbb*qw(i+1)		
		qg(i+1)    = qg(i+1)/(bbbb*(1.-qubar(i+1)))

		qnn(i+1) = qn(i+1)


c
c		solve for qge(i+1) using the intertemporal condition
c
c
		
		rrrr     = (1./bbbb)-1
		qge(i+1) = qnn(i+1)*qy(i+1)-depr*qh(i+1)
     $             -  rrrr*qh(i+1)
		qge(i+1) = qge(i+1)/(qnn(i+1)*qgam(i+1))
		qge(i+1) = xvalue-qge(i+1)
c		write(*,'(i6,5f14.6)') i+1,qg(i+1),qw(i+1),qge(i+1),
c     $    qhu(i+1),qh(i+1)				
200		continue
c		pause
		fout(1) = qg(iend+1)  - ststg
		fout(2) = qw(iend+1)  - ststw
		fout(3) = qge(iend+1) - ststge

	
c		write(*,'(6f12.6)') xin(1),xin(2),xin(3),
c     $    qg(iend+1),qw(iend+1),qge(iend+1)
		write(*,'(6f12.6)') qg(iend+1),qw(iend+1),qge(iend+1),
     $    fout(1),fout(2),fout(3)


c		pause	
		end

		real*8 function fnccaph1(hhhh)
		implicit real*8(a-h,o-z),integer(i-n)
		common/simpson/hsw(1001),nsimpson
		common/par/xvalue,alpha,prob,depr,bbbb,xi,rox,gam,eta,
     $          hstst,vcost,ststg,ststw,ststge,ubar
		common/var/qn(1000),qh(1000),qg(1000),qw(1000),
     $          qge(1000),qlabx(1000),qv(1000),qubar(1000),
     $          qy(1000),qa(1000),qnn(1000),qhu(1000),qgam(1000)                
		common/iiii/iend,ifunc


			
		hmean    = hhhh
		gamtemp =  fgam(hmean)
		ubar    =  1-gamtemp
		hutemp  =  hmean/gamtemp
		ytemp   =  qa(1)*((1.-depr)*hutemp)**alpha

				
c
c		calculate return
c

		xm  = xvalue - qge(1)
	
		xint  = rox*qn(1)*gamtemp*(ytemp-xm)
     $          - depr*hmean
		xint  = xint/hmean

		fnccaph1 = (1./bbbb) - 1. - xint
		
c		write(*,'(i6,2f14.6)') 1,hhhh,fnccaph1
		
		return
		end


		real*8 function fgam(hhhh)
		implicit real*8(a-h,o-z),integer(i-n)
		common/simpson/hsw(1001),nsimpson
		common/par/xvalue,alpha,prob,depr,bbbb,xi,rox,gam,eta,
     $          hstst,vcost,ststg,ststw,ststge,ubar
		common/var/qn(1000),qh(1000),qg(1000),qw(1000),
     $          qge(1000),qlabx(1000),qv(1000),qubar(1000),
     $          qy(1000),qa(1000),qnn(1000),qhu(1000),qgam(1000)                
		common/iiii/iend,ifunc


		if(ifunc.eq.1) fgam = (gam*hhhh)**eta/(1.+(gam*hhhh)**eta)
		if(ifunc.eq.2) fgam = 1./(1.+exp(-gam*(hhhh-eta)))
		if(ifunc.eq.3) fgam = gam*hhhh**eta
		
		return
		end


		real*8 function fgaminv(gggg)
		implicit real*8(a-h,o-z),integer(i-n)
		common/simpson/hsw(1001),nsimpson
		common/par/xvalue,alpha,prob,depr,bbbb,xi,rox,gam,eta,
     $          hstst,vcost,ststg,ststw,ststge,ubar
		common/var/qn(1000),qh(1000),qg(1000),qw(1000),
     $          qge(1000),qlabx(1000),qv(1000),qubar(1000),
     $          qy(1000),qa(1000),qnn(1000),qhu(1000),qgam(1000)                
		common/iiii/iend,ifunc


		if(ifunc.eq.1) then
		fgaminv = (gggg/(1.-gggg))**(1./eta)
		fgaminv = fgaminv/gam
		endif

		if(ifunc.eq.2) then
		fgaminv = log(((1./gggg)-1.))
		fgaminv = fgaminv/(-gam)
		fgaminv = fgaminv+eta
		endif

		if(ifunc.eq.3) then
		fgaminv = (gggg/gam)**(1./eta)
		endif

		return
		end