c       
c       a3agg.f, july 18 1995
c
c		like any3agg.f, but this program does not store 
c		the income levels so the number of agents can
c		be a lot higher
c
c
C ****************************************************************

	implicit real*8(A-H,O-Z)
	common y(10000),bj(10000),w1(101000),y1(101000),b1(101000),
     $  wdist(101000,2),agg(101000),cj(10000),c1(101000),
     $  q(101000),ya(101000),bet(15),sbet(15),
     $  umar(101000),temp(101000),bold(10000),bound,nper,np,num,i

	dimension prob(2,2)

		external func
		
		fl = 0.9
		fu = 1.1 
		facc = 0.00001
		ye   = 0.6228  
		yu   = 0.3772  


		open(31,file='any3a.in',status='old')
		open(32,file='any3a.new',status='unknown')
		open(33,file='any3a.par',status='old')
c
c
c		examples for input files:
c
c
c 
c  -.25  0.0001  0.96 .9   9  10000 10000  -333921  30   200   1   0.74228  -1.
c
c
c
c   5.13728858133224 -.669820846152612 .9549851025174561 -.865868672314076
c    .4312968801727971
c   -.331606048235028 .1509652888760695 -2.55757347431404
c   -.645164367738549 
c
c
c
	read(33,*) bound,acc,del,eps,np,nobs,nper,idum,numiter,
     $             iterq,num,prob(1,1),gam
     

		aggh  = 1.025
		aggl  = 0.975
		pagg  = 0.80

		nhalf = nper/2

c
c		nper has to be even
c

		ikeep = idum 

		read(31,*) (bet(j),j=1,np)

		do 140 j = 1,np
		sbet(j) = bet(j)
140		continue

		agg(1) = aggh

		do 141 i = 2,nobs

		xxx = ran1(idum)

		if(xxx.lt.pagg) then
		agg(i) = agg(i-1)
		else
		agg(i) = 2. -  agg(i-1)
		endif
141     continue

		do 300 iter = 1,numiter

		idum = -18795 

c
c
c               generating the series
c
c

		do 205 j = 1,nper
		bj(j) = 0.0
205     continue


		c1(1) = 0.5

		do 100 j = 1,nhalf
		y(j) = ye
100     continue

		do 210 i = 2,nobs
		
c
c               generate exogenous process
c


		do 110 j = 1,nhalf

		xxx = ran1(idum)                
c		if(i.eq.20.and.j.eq.2) write(*,*) xxx	
		if( y(j).gt.0.5 ) then

			if(xxx.lt.prob(1,1)) y(j) = ye
			if(xxx.ge.prob(1,1)) y(j) = yu

		else

			if(xxx.lt.prob(1,1)) y(j) = yu
			if(xxx.ge.prob(1,1)) y(j) = ye

		endif
110             continue


		do 120 j = nhalf+1,nper        
		
		y(j)  = 1. - y(j-nhalf)

120     continue


c
c
c       calculating information about the wealth distribution
c
c

		y1(i) = y(1)

		wdist(i,1) = 0.0
		wdist(i,2) = 0.0
		do 145 j = 1,nper
		wdist(i,1) = wdist(i,1) + bj(j)**2.
		wdist(i,2) = wdist(i,2) + bj(j)*y(j)
145     continue

		wdist(i,1) = wdist(i,1)/dble(nper)
		wdist(i,2) = wdist(i,2)/dble(nper)

		b1(i) =    bj(1)
		w1(i) =    (agg(i)*y(1) + bj(1) )  

		xxx  = bet(7)*wdist(i,1) + bet(8)*wdist(i,2) +
     $                 bet(9)*agg(i)

		xxx  = bet(1)*exp( xxx)

		do 214 j = 1,nper
		 
		temp(j)  = bet(2)*y(j)  +  bet(3)*bj(j)*y(j)  +
     $                     bet(4)*bj(j)   +  bet(5)*bj(j)**2.     +
     $                     bet(6)*bj(j)**3. 
     
		temp(j)  = xxx*exp(temp(j))
		bold(j)  = bj(j)
214     continue

c
c

		q(i) = rtsec(func,fl,fu,facc)
		do 220 j = 1,nper
		cj(j) = (temp(j)/q(i))**(-1.)    
		bj(j) = (  agg(i)*y(j) + bold(j) - cj(j) ) / q(i)
		if(bj(j).lt.bound) then
		bj(j) =  bound
		cj(j) = agg(i)*y(j) + bold(j) - bound*q(i)
		endif
220             continue
		
		umar(i) = del*cj(1)**gam*c1(i-1)**(-gam-1.) 
		c1(i) = cj(1)
210             continue


c
c               calculating sbeta
c
c
c               write(*,*)  '  start sbeta  '

		call sbeta(nobs)

c
c
c               updating beta
c
c



               write(*,500) (bet(i), i=1,5 )
               write(*,500) (bet(i), i=6,10)
		
		diff = 0.
		do 310 i=1,np
		diff = diff + abs( ( bet(i)-sbet(i) ) / bet(i) )
310     bet(i) = bet(i)*(1-eps) + sbet(i)*eps

               write(*,500) (sbet(i), i=1,5 )
               write(*,500) (sbet(i), i=6,10)
		
		write(*,'(a7,i6,f10.5)')  'diff= ' ,iter,diff
		if ( diff.lt.acc )
     $          go to 400 



300     continue

	write(*,'(a)')   '  iterations on beta did not converge'
	write(*,'(a)')   '  -----------------------------------'


400     continue


		write(*,500) (bet(i), i=1,5 )
		write(*,500) (bet(i), i=6,10)
		write(32,*) (bet(i), i=1,5 )
		write(32,*) (bet(i), i=6,10)

c
c
c       CALCULATING SOME SAMPLE STATISTICS
c
c

		cmean  = 0.0
		cstan  = 0.0
		rmean  = 0.0
		rstan  = 0.0
		ymean  = 0.0
		ystan  = 0.0
		bmean  = 0.0
		bstan  = 0.0
		sm     = 0.0
		ss     = 0.0

		do 720 i = 101,nobs
		wdist(i,1) = sqrt(wdist(i,1))
		cmean  = cmean  + log(c1(i)/c1(i-1))
		rmean  = rmean  - log(q(i))
		ymean  = ymean  + log(y1(i)/y1(i-1))
		bmean  = bmean  + b1(i)
		sm     = sm     + wdist(i,1) 
		write(99,*) i,agg(i),q(i)
720             continue

		cmean  = cmean/dble(nobs-100)
		rmean  = rmean/dble(nobs-100)
		ymean  = ymean/dble(nobs-100)
		bmean  = bmean/dble(nobs-100)
		sm     = sm   /dble(nobs-100)
	

		do 730  i = 101,nobs
		cstan  = cstan  + (log(c1(i)/c1(i-1)) - cmean)**2.
		rstan  = rstan  + (-log(q(i))-rmean)**2.
		ystan  = ystan  + (log(y1(i)/y1(i-1))-ymean)**2.
		bstan  = bstan  + (b1(i) - bmean)**2.
		ss     = ss     + (wdist(i,1)-sm)**2.
730             continue

		cstan  = sqrt(cstan/dble(nobs-100)  )
		rstan  = sqrt(rstan/dble(nobs-100)  )
		ystan  = sqrt(ystan/dble(nobs-100)  )
		bstan  = sqrt(bstan/dble(nobs-100)  )
		ss     = sqrt(ss   /dble(nobs-100)  )

		write(*,*)
		write(*,*)
		write(*,510) del,-log(del)*100.,bound
		write(*,*)
		write(*,516) np,nobs,nper,-gam
		write(*,*)        
		write(*,515) ikeep,iterq
		write(*,*)
		write(*,*)
		write(*,511) cmean,cstan
		write(*,512) ymean,ystan
		write(*,514) rmean,rstan
		write(*,517) bmean,bstan
		write(*,518) sm,ss
		write(*,*)
		write(*,*)


448             continue


500     format(' bet  :',6f13.6)
502     format(' not solved bond price iteration at obs.: ',i6)
510     format('   del = ',f6.3,'   discount rate (%) = ',f6.3,
     $  '   bound = ',f7.2)
511     format(' agent 1s dlog consumption:     ',2f10.6)
517     format(' agent 1s bond holdings:        ',2f10.6)
512     format(' agent 1s dlog income:          ',2f10.6)
518     format(' sd cross-sectional dist:       ',2f10.6)
514     format(' interest rate:                 ',2f10.6)
516     format('   np = ',i3,'     nobs = ',i8,'   nper =',i5,
     $  '   RRA = ',f7.2)  
515     format('   idum = ',i8,'    iterq = ',i5)


		stop
		end

C
C
C       
	SUBROUTINE INVERT(A,N,D)
C
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 A(N,N) , L(20) , M(20)
      common /ss/ ISING
      D=1.
      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(ABS(BIGA)-ABS(A(I,J))) 10,20,20
   10 BIGA=A(I,J)
      L(K)=I
      M(K)=J
   20 CONTINUE
      IF (ABS(BIGA).LT.(1.0E-7)) 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./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
  150 CONTINUE
      RETURN
      END
C    ************************************
      SUBROUTINE MULT(E , F , G , NROWE , NROWF , NCOLF)
C
C
C     THIS SUBROUTINE CALCULATES THE PRODUCT G OF TWO MATRICES
C     E AND F.
C
C
      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

C
C   *****************************************
C
C       THE FOLLOWING TWO SUBROUTINES GENERATE RANDOM NUMBERS
C
C   *****************************************
C
	REAL*8 FUNCTION GASDEV(IDUM)
c
c       retruns 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


C  *************************************************************
C
C     Calculation of S(beta)
C
C  *************************************************************

	Subroutine sbeta(time)
	implicit real*8(A-H,O-Z),integer(i-n)
	integer t,time
	common y(10000),bj(10000),w1(101000),y1(101000),b1(101000),
     $  wdist(101000,2),agg(101000),cj(10000),c1(101000),
     $  q(101000),ya(101000),bet(15),sbet(15),
     $  umar(101000),temp(101000),bold(10000),bound,nper,np,num,i

	real*8  dep, ind(9), coef(9), xx(9,9), xy(9)

c
c     a NLLS regression is done by running an OLS regression on the derivati-
c     ves of pred. 'ind' contains these derivatives, and 'dep'
c     has the dependent  variables in the corresponding OLS regression.
c
c
c     num  is the number of iterations.  very often you can speed up
c     convergence by setting num = 1.  The inaccuracy this creates
c     disappears in subsequent iterations on S(beta)
c

      Do 20 iter=1,num
	
	
cwrite(*,*) (sbet(kk),kk=1,np)


	 do 10 i= 1,np
	    xy(i) = 0.
	    do 10 j=i,np
  10           xx(i,j) = 0.

	 Do 30 t= 150, time-1
	    ind(1) = exp(
     $                   sbet(2)*y1(t)        +
     $                   sbet(3)*b1(t)*y1(t)  +
     $                   sbet(4)*b1(t)         +
     $                   sbet(5)*b1(t)**2.     +
     $                   sbet(6)*b1(t)**3.     +
     $                                   sbet(7)*wdist(t,1)    +
     $                                   sbet(8)*wdist(t,2)    +
     $                                   sbet(9)*agg(t)        ) 


	    prov   = ind(1) * sbet(1)
	    ind(2) = prov * y1(t)
	    ind(3) = prov * b1(t)*y1(t)
	    ind(4) = prov * b1(t) 
	    ind(5) = prov * b1(t)**2.
	    ind(6) = prov * b1(t)**3.
		ind(7) = prov * wdist(t,1)
		ind(8) = prov * wdist(t,2)    
		ind(9) = prov * agg(t)

		dep = umar(t+1)
		
		do 33 j = 2,np
		dep = dep + ind(j)*sbet(j)
33              continue

	 do 35 i= 1,np
	    xy(i) = xy(i) + ind(i)*dep
	    do 35 j=i,np
  35          xx(i,j) = xx(i,j) + ind(i) * ind(j)
  30     continue

	 do 40 i=2,np
	    do 40 j=1,i-1
  40           xx(i,j) = xx(j,i)

	 call invert(xx,np,d)
	 call mult(xx,xy,coef,np,np,1)

	 diff = 0.
	 do 50 i=1,np
	    diff = diff + abs( (sbet(i)-coef(i))/sbet(i) )
	  sbet(i) = coef(i)
	
50      continue

	 if (diff.lt.(.000001)) go to 100
   20    continue
      write(*,'(a)') ' NLLS not solved'
      write(*,'(a)') ' ---------------'
  100 continue
      end




      SUBROUTINE SORT(N,RA)
      implicit real*8(a-h,o-z)
      DIMENSION RA(N)
      L=N/2+1
      IR=N
10    CONTINUE
	IF(L.GT.1)THEN
	  L=L-1
	  RRA=RA(L)
	ELSE
	  RRA=RA(IR)
	  RA(IR)=RA(1)
	  IR=IR-1
	  IF(IR.EQ.1)THEN
	    RA(1)=RRA
	    RETURN
	  ENDIF
	ENDIF
	I=L
	J=L+L
20      IF(J.LE.IR)THEN
	  IF(J.LT.IR)THEN
	    IF(RA(J).LT.RA(J+1))J=J+1
	  ENDIF
	  IF(RRA.LT.RA(J))THEN
	    RA(I)=RA(J)
	    I=J
	    J=J+J
	  ELSE
	    J=IR+1
	  ENDIF
	GO TO 20
	ENDIF
	RA(I)=RRA
      GO TO 10
      END


	real*8 function func(qq)
	implicit real*8(a-h,o-z),integer(i-n)
	common y(10000),bj(10000),w1(101000),y1(101000),b1(101000),
     $  wdist(101000,2),agg(101000),cj(10000),c1(101000),
     $  q(101000),ya(101000),bet(15),sbet(15),
     $  umar(101000),temp(101000),bold(10000),bound,nper,np,num,i


	do 220 j = 1,nper
	cj(j) = qq/temp(j)    
	bj(j) = (  agg(i)*y(j) + bold(j) - cj(j) ) / qq
	if(bj(j).lt.bound) then
	bj(j) =  bound
	endif
220     continue

c
c
c       check aggregate demand conditions
c
c

	btot = 0.0
	do 235 j = 1,nper
	btot = btot + bj(j)
235     continue
	
	func  =  btot
	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


		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