c
c       xx5.f, January 1, 1995  nper = mutiple of 10
c
c       info about the variance of bj(j) in state space
c
C ****************************************************************

	implicit real*8(A-H,O-Z),integer(i-n)
	common/series/ y(100,1),bj(100),
     $  q(60),ya(60),ww(100),
     $  cj(100),polc(2,2,3,25,15,1),polb(2,2,3,25,15,1),
     $  polq(2,2,3,25,15,1),polbold(2,2,3,25,15,1),
     $  polcold(2,2,3,25,15,1),polqold(2,2,3,25,15,1),
     $  prede(100),s1(2),s2(20),i3(100),s4(30),s5(20),i6(10),
     $  s3(100),s6(100)
	common/par/ bound,nper,pagg,ye,yu,del,gam,
     $  m1,m2,m3,m4,m5,m6,kk,m1t,m2t,m3t,m4t,m5t,m6t
	common/ser/ jmax(10),yy(100,10)
	dimension ex(2,2),bjn(100),accmax(3),prob(2,2,2,3),
     $  intacc(3,6),vmean(100),var(100),varn(100)


		external func

	 open(33,file='todd.par',status='old')
c        open(34,file='todd.in',status='old')
	 open(35,file='todd.new',status='unknown')
	 read(33,*) bound,acc,del,nper,idum,numiter,gam,znu
	 if(nper.gt.100) write(*,*) 'nper too big, adjust dimensions '
	 if(nper.gt.100) stop
c
c        example for todd.par
c

c        -60. 0.0001 0.96 10 -333921 400 -1. .99 30.  


c
c
c               Note this program is not opening the file with the
c               initial guesses
c
c

		fl   = 0.9
		fu   = 1.1 
		facc = 0.00001
		ye   = 100.
		yu   = 50.
		pb0  = 0.5
		pb1  = 0.4
		pb2  = 0.1
		pg0   = 0.65
		pg1   = 0.3
		pg2   = 0.05 
		pagg = 0.9375

c
c
c
c               prob(i,j,k,l) is the probability that for next period
c
c               m1 = j, m2 = k and m3 = l given that m1 = i today
c
c
c

		prob(1,1,1,1) = pagg*pg0*0.
		prob(1,1,2,1) = pagg*pg0*1.
		prob(1,1,1,2) = pagg*pg1*0.1
		prob(1,1,2,2) = pagg*pg1*0.9
		prob(1,1,1,3) = pagg*pg2*0.2
		prob(1,1,2,3) = pagg*pg2*0.8

		prob(1,2,1,1) = (1.-pagg)*pb0*0.
		prob(1,2,2,1) = (1.-pagg)*pb0*1.
		prob(1,2,1,2) = (1.-pagg)*pb1*0.1
		prob(1,2,2,2) = (1.-pagg)*pb1*0.9
		prob(1,2,1,3) = (1.-pagg)*pb2*0.2
		prob(1,2,2,3) = (1.-pagg)*pb2*0.8

		prob(2,1,1,1) = (1.-pagg)*pg0*0.
		prob(2,1,2,1) = (1.-pagg)*pg0*1.
		prob(2,1,1,2) = (1.-pagg)*pg1*0.1
		prob(2,1,2,2) = (1.-pagg)*pg1*0.9
		prob(2,1,1,3) = (1.-pagg)*pg2*0.2
		prob(2,1,2,3) = (1.-pagg)*pg2*0.8

		prob(2,2,1,1) = pagg*pb0*0.
		prob(2,2,2,1) = pagg*pb0*1.
		prob(2,2,1,2) = pagg*pb1*0.1
		prob(2,2,2,2) = pagg*pb1*0.9
		prob(2,2,1,3) = pagg*pb2*0.2
		prob(2,2,2,3) = pagg*pb2*0.8

c
c
c               setting up grid
c
c               s1: aggregate state
c               s2: individual income status 
c               i3: number unemployed agents
c               s4: individual bond level
c               s5: "pseudo" sd of bondholdings excluding agent 1

 
		m1t = 2
		m2t = 2 
		m3t = 3   
		m4t =25 
		m5t =11 
 
		s1(1) =  10
		s1(2) = -10

		s2(1) = yu
		s2(2) = ye

		i3(1) =    0 
		i3(2) =    nper/10
		i3(3) =    nper/5

		do 22 jj = 1,m4t
		s4(jj) = bound+2.*abs(bound)*dble(jj-1)/dble(m4t-1)
		write(*,*) jj,s4(jj)
22              continue

		do 23 jj = 1,m5t
		s5(jj) = dble(jj)*5.
23              continue

	

c
c               initializing the stuff for the cons and bond interpolation
c

		jmax(1) = m4t
		jmax(2) = m5t
		do 206 ii = 1,m4t
		yy(ii,1) = s4(ii)
206             continue
		do 207 ii = 1,m5t
		yy(ii,2) = s5(ii)
207             continue


c
c
c               initializing the policy functions
c
c
		do 440 m1 = 1,m1t
		do 440 m2 = 1,m2t
		do 440 m3 = 1,m3t
		do 440 m4 = 1,m4t
		do 440 m5 = 1,m5t
c                read(34,*) ij1,ij2,ij3,ij4,ij5,
c     $          polcold(m1,m2,m3,m4,m5,1), 
c     $          polbold(m1,m2,m3,m4,m5,1), 
c     $          polqold(m1,m2,m3,m4,m5,1) 
	       polcold(m1,m2,m3,m4,m5,1) = s2(m2)
	       polbold(m1,m2,m3,m4,m5,1) = 0.0001 
	       polqold(m1,m2,m3,m4,m5,1) = del
     
440             continue

c
c               this program does not read in an initial guess
c               but has as the initial guess the "never borrow"
c               policy functions
c
c
c
c               start iterating
c
c

		do 520 kk = 2,numiter
		write(*,*) kk
		mcount = 0
		mtot = m1t*m2t*m3t*m4t*m5t

		do 510 m1 = 1,m1t
		do 510 m2 = 1,m2t
		do 510 m3 = 1,m3t
		do 510 m4 = 1,m4t
		do 510 m5 = 1,m5t

		mcount = mcount + 1
		
		
		if(m2.eq.1.and.m3.eq.1) then
		polc(m1,m2,m3,m4,m5,1) = 999. 
		polb(m1,m2,m3,m4,m5,1) = 999. 
		polq(m1,m2,m3,m4,m5,1) = 999.
		endif 
		if(m2.eq.1.and.m3.eq.1) goto 510 

		bj(1) = s4(m4)
c
c               calculating the amount of debt for the other people
c

		ztot  =  dble(nper-1)
		zhalf =  dble(nper)/2.
		zhalf1 = dble(nper)/2. - 1
		root2 = -(ztot*s5(m5)**2.
     $                /(zhalf**2./zhalf1+zhalf) )**0.5
		root1 = -zhalf*root2/zhalf1
		xyz = -s4(m4)/dble(nper-1)
		bond2 = root1+xyz
		bond1 = root2+xyz


		jik = 10
		do 5100 j = 2,nper
		if(jik.gt.0) then
		bj(j) = bond1
		else
		bj(j) = bond2
		endif
		jik = -jik
5100            continue
		if(m2.eq.1) then
		iun = i3(m3) - 1
		else
		iun = i3(m3)
		endif

c
c               note that iun.ge.0
c

		do 5130 j = 2,nper-iun
		y(j,1) = ye
		ww(j)  = y(j,1)+bj(j)
5130            continue

		if(iun.eq.0) goto 5140

		do 5131 j = nper-iun+1,nper
		y(j,1) = yu
		ww(j)  = y(j,1)+bj(j)
5131            continue
		
5140            continue

		y(1,1) = s2(m2)
		ww(1) = s2(m2)+bj(1)

c
c               to calculate next period's consumption you need to know
c               the realizations to the stochastic shocks (aggregate and
c               individual) and today's savings level.  For today's
c               savings level we choose the level from the last iteration 
c               (=polbold).  Note that this is not necessary, but prevents 
c               the need to use a non-linear equation solver here.

		do 390 j = 1,nper
		vmean(j) = -bj(j)/dble(nper-1)
390             continue
		do 395 j = 1,nper
		var(j)   = 0.0
		do 395 jk = 1,nper
		if(j.ne.jk) var(j) = var(j) + (bj(jk)-vmean(j))**2.
395             continue
		do 397 j = 1,nper
		var(j) = (var(j)/dble(nper-1))**0.5
397             continue



		bjn(1) = pol(m1,m2,m3,bj(1),s5(m5),2)

		do 400 j = 2,nper

		if(y(j,1).eq.yu) then
		m2j = 1
		else
		m2j = 2
		endif   

		s5j = var(j)

		bjn(j) = pol(m1,m2j,m3,bj(j),s5j,2)

400             continue

c
c               calculate some aggregate statistics for next period's
c               bond holdings
c
c


		do 1390 j = 1,nper
		vmean(j) = -bjn(j)/dble(nper-1)
1390            continue
		do 1395 j = 1,nper
		varn(j)   = 0.0
		do 1395 jk = 1,nper
		if(j.ne.jk) varn(j) = varn(j) + (bjn(jk)-vmean(j))**2.
1395            continue
		do 1397 j = 1,nper
		varn(j) = (varn(j)/dble(nper-1))**0.5
1397            continue


c
c               calculating conditional expectations
c

		do 5 j = 1,nper

		s5jn = varn(j)
c
		prede(j) = 0.0

		do 20 k1 = 1,m1t
		do 20 k2 = 1,m2t
		do 20 k3 = 1,m3t
		aaa = (pol(k1,k2,k3,bjn(j),s5jn,1))**gam
		prede(j) = prede(j)+prob(m1,k1,k2,k3)*aaa
20              continue

		prede(j) = del*prede(j)
5               continue


		polq(m1,m2,m3,m4,m5,1) = rtsec(func,fl,fu,facc)
		polc(m1,m2,m3,m4,m5,1) 
     $          = (prede(1)/polq(m1,m2,m3,m4,m5,1))**(1./gam)
		polb(m1,m2,m3,m4,m5,1) 
     $          = (  ww(1) - polc(m1,m2,m3,m4,m5,1) ) 
     $          / polq(m1,m2,m3,m4,m5,1)
		if(polb(m1,m2,m3,m4,m5,1).lt.bound) then
		polb(m1,m2,m3,m4,m5,1) = bound
		polc(m1,m2,m3,m4,m5,1) = ww(1) 
     $          - bound*polq(m1,m2,m3,m4,m5,1)
		endif


510             continue

		do 560 mm = 1,5
		do 560 nn = 1,3
		intacc(nn,mm) = 0
560             continue

		do 570 nn = 1,3
		accmax(nn) = 0.
570             continue
	
		do 550 m1 = 1,m1t
		do 550 m2 = 1,m2t
		do 550 m3 = 1,m3t
		do 550 m4 = 1,m4t
		do 550 m5 = 1,m5t

		do 550 nn =  1,3

		if(nn.eq.1) then
		cnew = polq(m1,m2,m3,m4,m5,1)
		cold = polqold(m1,m2,m3,m4,m5,1)
		endif
		if(nn.eq.2) then
		cnew = polb(m1,m2,m3,m4,m5,1)
		cold = polbold(m1,m2,m3,m4,m5,1)
		endif
		if(nn.eq.3) then
		cnew = polc(m1,m2,m3,m4,m5,1)
		cold = polcold(m1,m2,m3,m4,m5,1)
		endif

		if(cold.ne.0.) dev = abs(log(abs(cnew/cold)))

		if(nn.eq.2) then
		if(abs(cnew).le.1.) dev = (cnew-cold) /10.
		endif

		if(dev.gt.accmax(nn)) then
		accmax(nn) = dev
		intacc(nn,1) = m1
		intacc(nn,2) = m2
		intacc(nn,3) = m3
		intacc(nn,4) = m4
		intacc(nn,5) = m5
		endif
550             continue

		write(*,*) (accmax(nn),nn=1,3)
		write(*,*) (intacc(1,mm),mm=1,1)
		write(*,*) (intacc(2,mm),mm=1,1)
		write(*,*) (intacc(3,mm),mm=1,1)

		write(*,*) accmax(1)+accmax(2)+accmax(3),acc
		if(accmax(1)+accmax(2)+accmax(3).lt.acc) goto 525


		do 580 m1 = 1,m1t
		do 580 m2 = 1,m2t
		do 580 m3 = 1,m3t
		do 580 m4 = 1,m4t
		do 580 m5 = 1,m5t
	polcold(m1,m2,m3,m4,m5,1)=znu*polc(m1,m2,m3,m4,m5,1)
     #  +(1.-znu)*polcold(m1,m2,m3,m4,m5,1)
	polbold(m1,m2,m3,m4,m5,1)=znu*polb(m1,m2,m3,m4,m5,1)
     #  +(1.-znu)*polbold(m1,m2,m3,m4,m5,1)
	polqold(m1,m2,m3,m4,m5,1)=znu*polq(m1,m2,m3,m4,m5,1)
     #  +(1.-znu)*polqold(m1,m2,m3,m4,m5,1)
580             continue


520             continue


525                     continue


		do 530 m1 = 1,m1t
		do 530 m2 = 1,m2t
		do 530 m3 = 1,m3t
		do 530 m4 = 1,m4t
		do 530 m5 = 1,m5t
		write(35,'(5i4,3f16.8)') m1,m2,m3,m4,m5,
     $          polc(m1,m2,m3,m4,m5,1), 
     $          polb(m1,m2,m3,m4,m5,1), 
     $          polq(m1,m2,m3,m4,m5,1) 
530             continue



		stop
		end





		real*8 function func(qq)
		implicit real*8(a-h,o-z),integer(i-n)
		common/series/ y(100,1),bj(100),
     $          q(60),ya(60),ww(100),
     $          cj(100),polc(2,2,3,25,15,1),polb(2,2,3,25,15,1),
     $          polq(2,2,3,25,15,1),polbold(2,2,3,25,15,1),
     $          polcold(2,2,3,25,15,1),polqold(2,2,3,25,15,1),
     $          prede(100),s1(2),s2(20),i3(100),s4(30),s5(20),i6(10),
     $          s3(100),s6(100)
		common/par/ bound,nper,pagg,ye,yu,del,gam,
     $          m1,m2,m3,m4,m5,m6,kk,m1t,m2t,m3t,m4t,m5t,m6t


		dimension save(100)


		do 220 j = 1,nper


		cj(j) = (prede(j)/qq)**(1./gam)    
		save(j) = (  ww(j) - cj(j) ) / qq
		if(save(j).lt.bound) then
		save(j) =  bound
		endif
220             continue

c
c
c       check aggregate demand conditions
c
c

		btot = 0.0
		do 235 j = 1,nper
		btot = btot + save(j)
235             continue

		func  =  btot

		end




		real*8 function pol(m1p,m2p,m3p,s4p,s5p,kkk) 
		implicit real*8(a-h,o-z),integer(i-n)
		common/series/ y(100,1),bj(100),
     $          q(60),ya(60),ww(100),
     $          cj(100),polc(2,2,3,25,15,1),polb(2,2,3,25,15,1),
     $          polq(2,2,3,25,15,1),polbold(2,2,3,25,15,1),
     $          polcold(2,2,3,25,15,1),polqold(2,2,3,25,15,1),
     $          prede(100),s1(2),s2(20),i3(100),s4(30),s5(20),i6(10),
     $          s3(100),s6(100)
		common/par/ bound,nper,pagg,ye,yu,del,gam,
     $          m1,m2,m3,m4,m5,m6,kk,m1t,m2t,m3t,m4t,m5t,m6t
		common/ser/ jmax(10),yy(100,10)
		dimension z3(2,2,2),z2(2,2),z1(2),
     $          xx(4),t(4),ll(4,2),z4(2,2,2,2)



		call locate(s4,m4t,s4p,j1)
		call locate(s5,m5t,s5p,j2)

		ll(1,1) = j1
		ll(2,1) = j2
		xx(1) = s4p
		xx(2) = s5p

		do 10 ii = 1,2

		if(ll(ii,1).eq.0) then
		ll(ii,1) = 1
		endif 

		if(ll(ii,1).eq.jmax(ii)) then
		ll(ii,1) = ll(ii,1) -1
		endif

		ll(ii,2) = ll(ii,1) + 1

		t(ii) = (xx(ii)-yy(ll(ii,1),ii)) 
     $                / (yy(ll(ii,2),ii)-yy(ll(ii,1),ii))

		
10              continue
	

		do 20 l1 = 1,2
		do 20 l2 = 1,2
		if(kkk.eq.1) z4(l1,l2,1,1) 
     $          = polcold(m1p,m2p,m3p,ll(1,l1),ll(2,l2),1  )
		if(kkk.eq.2) z4(l1,l2,1,1) 
     $          = polbold(m1p,m2p,m3p,ll(1,l1),ll(2,l2),1  )
20              continue

		do 30 l2 = 1,2
		z3(l2,1,1) = (1.-t(1))*z4(1,l2,1,1)
     $                      +    t(1) *z4(2,l2,1,1)
30              continue

		pol =    (1.-t(2))*z3(1,1,1)
     $              +        t(2) *z3(2,1,1)

		if(kkk.eq.2) then
		if(pol.lt.bound) pol = bound
		endif
		
		end


		subroutine locate(xx,n,x,j)
		implicit real*8(a-h,o-z),integer(i-n)
		dimension xx(n)
		jl = 0
		ju = n+1
10              if(ju-jl.gt.1) then
		jm=(ju+jl)/2
			if((xx(n).gt.xx(1)).eqv.(x.gt.xx(jm)))then
			jl=jm
			else
			ju=jm
			endif
		goto 10
		endif
		j=jl
		return
		end



		real*8 function perm(j,k,nper)
		implicit real*8(a-h,o-z),integer(i-n)


		if(j.eq.0) then
		perm = 1
		goto 100
		endif


		s2 = 1
		s3 = 1


		if(j.ge.k) then
		nhigh = j
		nlow  = k
		else
		nhigh = k
		nlow  = j
		endif

		do 10 i = nhigh+1,nper
		s3 = s3*dble(i)
10              continue
		do 20 i = 1,nlow
		s2 = s2*dble(i)
20              continue
		perm = s3/s2


100             continue
		

		end

		





      REAL*8 FUNCTION RTBIS(FUNC,X1,X2,XACC)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
      PARAMETER (JMAX=100)
      EXTERNAL FUNC
      FMID=FUNC(X2)
      F=FUNC(X1)
      IF(F*FMID.GE.0.) PAUSE 'Root must be bracketed for bisection.'
      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.0.)RTBIS=XMID
	IF(ABS(DX).LT.XACC .OR. FMID.EQ.0.) RETURN
11    CONTINUE
      PAUSE 'too many bisections'
      END
      REAL*8 FUNCTION RTSEC(FUNC,X1,X2,XACC)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
      PARAMETER (MAXIT=30)
      EXTERNAL FUNC
      FL=FUNC(X1)
      F=FUNC(X2)
      IF(ABS(FL).LT.ABS(F))THEN
	RTSEC=X1
	XL=X2
	SWAP=FL
	FL=F
	F=SWAP
      ELSE
	XL=X1
	RTSEC=X2
      ENDIF
      DO 11 J=1,MAXIT
	DX=(XL-RTSEC)*F/(F-FL)
	XL=RTSEC
	FL=F
	RTSEC=RTSEC+DX
	F=FUNC(RTSEC)
	IF(ABS(DX).LT.XACC.OR.F.EQ.0.)RETURN
11    CONTINUE
      PAUSE 'RTSEC exceed maximum iterations'
      END