c     Solves the Bewley-Aiyagari model
c	Adds entrepreneur sector and occupational choice
c     3-state case
      implicit real*8 (a-h,o-z)
      parameter (nk2pts=5000,nkpts=150,nread=0,nread2=0,nread3=0,nx=20)
	common/cfs/f0(nk2pts,3,2),f1(nk2pts,3,2)
	common/cv2/y2(nkpts,3),optk12(nkpts),optk02(nkpts),optk22(nkpts),
     &ye2(nkpts,3),optke12(nkpts),optke02(nkpts),optke22(nkpts),
     &ioccup12(nkpts),ioccupe12(nkpts),ioccup02(nkpts),ioccupe02(nkpts),
     &ioccup22(nkpts),ioccupe22(nkpts)
      common/cprice/r
	common/cwages/ylo,ymi,yhi
	dimension xkpts(nkpts),x(nx),y(nx)
	external zeta2

      one = 1.0D+00
	zero=0.0D+00

	ylo=0.57D+00
	ymi=0.93D+00
	yhi=1.51D+00

      if (nread2 .eq. 1) then
         open(2,file='entrep.f1',form='formatted')
         do 66 i = 1,nk2pts
            read(2,547) dum,f1(i,1,1),f1(i,2,1), f1(i,3,1),
     &				  f1(i,1,2),f1(i,2,2),f1(i,3,2)
 547        format(7f18.10)
 66      continue
         close(2)
      else

         do 49 i = 1,nk2pts
            do 48 j = 1,3
			do 47 k=1,2
               f1(i,j,k) = 0.0D+00
 47			continue
 48         continue
 49      continue
         f1(275,1,1) = one/6.0D+00
         f1(275,2,1) = one/6.0D+00
	   f1(275,3,1) = one/6.0D+00
	   f1(275,1,2) = one/6.0D+00
         f1(275,2,2) = one/6.0D+00
	   f1(275,3,2) = one/6.0D+00
      endif   

      if (nread .eq. 1) then
         open(2,file='entrep.out',form='formatted')
         do 98 i = 1,nkpts
            read(2,700) dum,y2(i,1),y2(i,2),y2(i,3),
     &	  ye2(i,1),ye2(i,2),ye2(i,3),
     &  optk02(i),optk12(i),optk22(i),optke02(i),optke12(i),optke22(i)
 700        format(13f18.10)
 98      continue
         close(2)
      else
	   open(2,file='xkpts150.in',form='formatted')
         do 97 i = 1,nkpts
	      read(2,390) xkpts(i)
 390        format(f18.10)
            do 96 j = 1,3
               y2(i,j) = dlog(xkpts(i))
			 ye2(i,j) = dlog(xkpts(i))
 96         continue
		  optk22(i) = xkpts(i)
            optk12(i) = xkpts(i)
            optk02(i) = xkpts(i)
		  optke22(i) = xkpts(i)
		  optke12(i) = xkpts(i)
            optke02(i) = xkpts(i)
 97      continue
         close(2)
      endif
	

	   x=0.033946D+00

	   y=zeta2(x)

	   open(2,file='entrep.out3',form='formatted')
	   do 361 i = 1,nx
c	   write(2,490) y(i)
 361	   continue
	   write(2,490) y
 490	   format(f18.10)
	   close(2)

      end

      double precision function zeta2(x)
      implicit real*8 (a-h,o-z)
      common/cprice/r	
c	common/cparam/beta,z1
c	common/cres/xkgdpc,aratio
	common/cwages/ylo,ymi,yhi
      parameter (nkpts=150,delta=0.06D+00,alphac=0.366D+00,
     &        hfix=0.3271D+00,xkbor=0.0D+00,phi=0.5D+00,acoef=1.0D+00,
     &           mxiter=400,nk2pts=5000,mxpol=200,npol=0,beta=0.9D+00)
      common/cv2/y2(nkpts,3),optk12(nkpts),optk02(nkpts),optk22(nkpts),
     &ye2(nkpts,3),optke12(nkpts),optke02(nkpts),optke22(nkpts),
     &ioccup12(nkpts),ioccupe12(nkpts),ioccup02(nkpts),ioccupe02(nkpts),
     &ioccup22(nkpts),ioccupe22(nkpts)
      common/cvs/v0(nkpts),v0dp(nkpts),v1(nkpts),
     &           v1dp(nkpts),v2(nkpts),v2dp(nkpts),
     &           ve0(nkpts),ve0dp(nkpts),ve1(nkpts),
     &           ve1dp(nkpts),ve2(nkpts),ve2dp(nkpts)
      common/cfs/f0(nk2pts,3,2),f1(nk2pts,3,2)
      common/cprob/piw(3,3),pie(3,3),qw(3),qe(3)
      dimension y(nkpts,3),ye(nkpts,3),
     &          optk1(nkpts),xkpts(nkpts),optkg1(nk2pts),
     &          vgrid(nk2pts,3),vegrid(nk2pts,3),
     &          optk0(nkpts),xk2pts(nk2pts),optkg0(nk2pts),
     &optk2(nkpts),optkg2(nk2pts),optke2(nkpts),optkge2(nk2pts),
     &optke1(nkpts),optkge1(nk2pts),optke0(nkpts),optkge0(nk2pts),
     &xp(3,nkpts),xprofit(3),xxp2(3,nk2pts),xprofit2(3),
     &xk(3,nkpts),xl(3,nkpts),xy(3,nkpts),xxk2(3,nk2pts),
     &xxl2(3,nk2pts),xxy2(3,nk2pts),
     &ioccup1(nkpts),ioccupe1(nkpts),ioccup0(nkpts),ioccupe0(nkpts),
     &ioccup2(nkpts),ioccupe2(nkpts),
     &ioccupg1(nk2pts),ioccupge1(nk2pts),
     &ioccupg0(nk2pts),ioccupge0(nk2pts),
     &ioccupg2(nk2pts),ioccupge2(nk2pts)	

      r = x

      one = 1.0D+00
      xkcut1 = 0.2D+00
      xkcut2 = 0.4D+00
	xkcut3 = 0.6D+00

c	increasing pie(1,1) can greatly increase the gini coeff		
	pie(1,1)=0.62D+00
	pie(1,2)=0.2D+00
	pie(1,3)=one-pie(1,1)-pie(1,2)

	pie(2,1)=0.1D+00
	pie(2,2)=0.7D+00
	pie(2,3)=one-pie(2,1)-pie(2,2)

	pie(3,1)=0.0D+00
	pie(3,2)=0.1D+00
	pie(3,3)=one-pie(3,1)-pie(3,2)

	piw(1,1)=0.75D+00
	piw(1,2)=0.24D+00
	piw(1,3)=0.01D+00

	piw(2,1)=0.19D+00
	piw(2,2)=0.62D+00
	piw(2,3)=0.19D+00
	
	piw(3,1)=0.01D+00
	piw(3,2)=0.24D+00
	piw(3,3)=0.75D+00


	qw(1)=0.31D+00
	qw(2)=0.38D+00
	qw(3)=0.31D+00

c	qw(1)=0.6D+00
c	qw(2)=0.3D+00
c	qw(3)=0.1D+00

	qe(1)=0.6D+00
	qe(2)=0.3D+00
	qe(3)=0.1D+00

      open(2,file='xkpts150.in',form='formatted')
      do 999 i = 1,nkpts
         read(2,570) xkpts(i)
 570     format(f18.10)
 999  continue

      close(2)
      xklow = xkbor
      xkhigh = 1250.0D+00
      xk2pts(1) = xkpts(1)
      xkinc = (xkpts(nkpts)-xkpts(1))/dble(real(nk2pts-1))
      do 393 i = 2,nk2pts
         xk2pts(i) = xk2pts(i-1)+xkinc
 393  continue

      w = (one-alphac)*(((r+delta)/alphac)**(alphac/(alphac-one)))
	w3=yhi*w
	w2=ymi*w
	w1=ylo*w

	vmxdif = 100.0D+00
	vmxdife = 100.0D+00

      do 30 kkk = 1,mxiter

	   if (vmxdif .lt. 1.0D-00 .and. vmxdife .lt. 1.0D+00)then
	      ncalop = 1
	   else
	      ncalop = 1
	   endif                 

         do 12 i = 1,nkpts
            do 11 j = 1,3
               y(i,j) = y2(i,j)
			 ye(i,j) = ye2(i,j)
 11         continue
            optk1(i) = optk12(i)
            optk0(i) = optk02(i)
		  optk2(i) = optk22(i)
		  optke1(i) = optke12(i)
            optke0(i) = optke02(i)
		  optke2(i) = optke22(i)
 12      continue

         call interp(xkpts,y,v0,v1,v2,v0dp,v1dp,v2dp)
	   call interp(xkpts,ye,ve0,ve1,ve2,ve0dp,ve1dp,ve2dp)

	   call prof1(r,w,xk,xl,xy,xp)
C	  
	  open(22,file='profit1.out',form='formatted')	
		do 81 i = 1,nkpts
          write(22,499) xkpts(i),xk(1,i),xk(2,i),xk(3,i),xl(1,i),
     &	xl(2,i),xl(3,i),xy(1,i),xy(2,i),
     &	xy(3,i),xp(1,i),xp(2,i),xp(3,i)      
 499      format(13f18.10)
 81       continue
		close(22)
	 
         do 13 i = 1,nkpts
            xkcur = xkpts(i)

		do 133 j=1,3
		  xprofit(j)=xp(j,i)
 133	    continue

c	worker and have low productivity
            neps = 0
	      nocc=1
            if (xkcur .le. xkcut1 .or. ncalop .eq. 1) then
               call calopt(xkcur,neps,nocc,xkp0,fval0,niter,
     &             1.0D-10,r,w,xkhigh,xklow,grad0,xkpts,xprofit,iocc0)
            else
               call altopt(xkcur,neps,nocc,xkp0,fval0,niter,
     &             1.0D-10,optk0(i),grad0,xkpts,r,w,xprofit,iocc0)
            endif
            optk02(i) = xkp0
            y2(i,1) = fval0
		  ioccup02(i)=iocc0

c	worker and have mid productivity
            neps = 1
		  nocc=1
            if (xkcur .le. xkcut1 .or. ncalop .eq. 1) then
               call calopt(xkcur,neps,nocc,xkp1,fval1,niter,
     &           1.0D-10,r,w,xkhigh,xklow,grad1,xkpts,xprofit,iocc1)
            else
               call altopt(xkcur,neps,nocc,xkp1,fval1,niter,
     &           1.0D-10,optk1(i),grad1,xkpts,r,w,xprofit,iocc1)
            endif
            optk12(i) = xkp1
            y2(i,2) = fval1
		  ioccup12(i)=iocc1

c	worker and have high productivity
            neps = 2
		  nocc=1
            if (xkcur .le. xkcut3 .or. ncalop .eq. 1) then
               call calopt(xkcur,neps,nocc,xkp2,fval2,niter,
     &           1.0D-10,r,w,xkhigh,xklow,grad2,xkpts,xprofit,iocc2)
            else
               call altopt(xkcur,neps,nocc,xkp2,fval2,niter,
     &           1.0D-10,optk2(i),grad2,xkpts,r,w,xprofit,iocc2)
            endif
            optk22(i) = xkp2
            y2(i,3) = fval2
		  ioccup22(i)=iocc2
c----------------------------------------------------------------------
c	entrepreneur and have bad shock
            neps = 0
	      nocc=0
		  if (xkcur .le. xkcut1 .or. ncalop .eq. 1) then
               call calopt(xkcur,neps,nocc,xkp3,fval3,niter,
     &            1.0D-10,r,w,xkhigh,xklow,grade3,xkpts,xprofit,iocc3)
            else
               call altopt(xkcur,neps,nocc,xkp3,fval3,niter,
     &            1.0D-10,optke0(i),grade3,xkpts,r,w,xprofit,iocc3)
            endif
            optke02(i) = xkp3
            ye2(i,1) = fval3
		  ioccupe02(i)=iocc3

c	entrepreneur and have mid shock
            neps = 1
		  nocc=0
            if (xkcur .le. xkcut2 .or. ncalop .eq. 1) then
               call calopt(xkcur,neps,nocc,xkp4,fval4,niter,
     &           1.0D-10,r,w,xkhigh,xklow,grade4,xkpts,xprofit,iocc4)
            else
               call altopt(xkcur,neps,nocc,xkp4,fval4,niter,
     &           1.0D-10,optke1(i),grade4,xkpts,r,w,xprofit,iocc4)
            endif
            optke12(i) = xkp4
            ye2(i,2) = fval4
		  ioccupe12(i)=iocc4

c	entrepreneur and have high shock
		  neps = 2
		  nocc=0
            if (xkcur .le. xkcut3 .or. ncalop .eq. 1) then
               call calopt(xkcur,neps,nocc,xkp5,fval5,niter,
     &           1.0D-10,r,w,xkhigh,xklow,grade5,xkpts,xprofit,iocc5)
            else
               call altopt(xkcur,neps,nocc,xkp5,fval5,niter,
     &           1.0D-10,optke2(i),grade5,xkpts,r,w,xprofit,iocc5)
            endif
            optke22(i) = xkp5
            ye2(i,3) = fval5
		  ioccupe22(i)=iocc5

 13      continue  

	   if (npol .eq. 0) goto 2
         do 43 ii = 1,mxpol

	      call interp(xkpts,y2,v0,v1,v2,v0dp,v1dp,v2dp)
		  call interp(xkpts,ye2,ve0,ve1,ve2,ve0dp,ve1dp,ve2dp)

	      do 42 i = 1,nkpts

	         xkcur = xkpts(i)

			do 136 j=1,3
			 xprofit(j)=xp(j,i)
 136			continue

c	worker and low productivity
		     con1 = (r+one)*xkcur+w1*hfix-optk02(i)			
	         call splint(xkpts,v0,v0dp,optk02(i),yval1,yp,
     &                     ydp,khi,klo,0,nkpts,1,0,0)
			 call splint(xkpts,v1,v1dp,optk02(i),yval2,yp,
     &                     ydp,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,v2,v2dp,optk02(i),yval3,yp,
     &                     ydp,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,ve0,ve0dp,optk02(i),yvale1,ype,
     &                     ydpe,khi,klo,0,nkpts,1,0,0)
			 call splint(xkpts,ve1,ve1dp,optk02(i),yvale2,ype,
     &                     ydpe,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,ve2,ve2dp,optk02(i),yvale3,ype,
     &                     ydpe,khi,klo,1,nkpts,1,0,0)
			wvalue=piw(1,1)*yval1+piw(1,2)*yval2+piw(1,3)*yval3			
			evalue=qe(1)*yvale1+qe(2)*yvale2+qe(3)*yvale3-phi  
	
			if (wvalue .gt. evalue)then
			y2(i,1) = dlog(con1)+beta*wvalue
			else
			y2(i,1) = dlog(con1)+beta*evalue
			endif

c	worker and mid productivity
		     con2 = (r+one)*xkcur+w2*hfix-optk12(i)				
	         call splint(xkpts,v0,v0dp,optk12(i),yval1,yp,
     &                     ydp,khi,klo,0,nkpts,1,0,0)
			 call splint(xkpts,v1,v1dp,optk12(i),yval2,yp,
     &                     ydp,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,v2,v2dp,optk12(i),yval3,yp,
     &                     ydp,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,ve0,ve0dp,optk12(i),yvale1,ype,
     &                     ydpe,khi,klo,0,nkpts,1,0,0)
			 call splint(xkpts,ve1,ve1dp,optk12(i),yvale2,ype,
     &                     ydpe,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,ve2,ve2dp,optk12(i),yvale3,ype,
     &                     ydpe,khi,klo,1,nkpts,1,0,0)
			wvalue=piw(2,1)*yval1+piw(2,2)*yval2+piw(2,3)*yval3			
			evalue=qe(1)*yvale1+qe(2)*yvale2+qe(3)*yvale3-phi 	
			if (wvalue .gt. evalue)then
			y2(i,2) = dlog(con2)+beta*wvalue
			else
			y2(i,2) = dlog(con2)+beta*evalue
			endif

c	worker and high productivity
		     con3 = (r+one)*xkcur+w3*hfix-optk22(i)				
	         call splint(xkpts,v0,v0dp,optk22(i),yval1,yp,
     &                     ydp,khi,klo,0,nkpts,1,0,0)
			 call splint(xkpts,v1,v1dp,optk22(i),yval2,yp,
     &                     ydp,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,v2,v2dp,optk22(i),yval3,yp,
     &                     ydp,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,ve0,ve0dp,optk22(i),yvale1,ype,
     &                     ydpe,khi,klo,0,nkpts,1,0,0)
			 call splint(xkpts,ve1,ve1dp,optk22(i),yvale2,ype,
     &                     ydpe,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,ve2,ve2dp,optk22(i),yvale3,ype,
     &                     ydpe,khi,klo,1,nkpts,1,0,0)
			wvalue=piw(3,1)*yval1+piw(3,2)*yval2+piw(3,3)*yval3			
			evalue=qe(1)*yvale1+qe(2)*yvale2+qe(3)*yvale3-phi
			if (wvalue .gt. evalue)then
			y2(i,3) = dlog(con3)+beta*wvalue
			else
			y2(i,3) = dlog(con3)+beta*evalue
			endif
c---------------------------------------------------------------------------	        
c	entrepreneur and low productivity
	cone1 = (r+one)*xkcur+xprofit(1)-optke02(i)						
	         call splint(xkpts,v0,v0dp,optke02(i),yval1,yp,
     &                     ydp,khi,klo,0,nkpts,1,0,0)
			 call splint(xkpts,v1,v1dp,optke02(i),yval2,yp,
     &                     ydp,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,v2,v2dp,optke02(i),yval3,yp,
     &                     ydp,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,ve0,ve0dp,optke02(i),yvale1,ype,
     &                     ydpe,khi,klo,0,nkpts,1,0,0)
			 call splint(xkpts,ve1,ve1dp,optke02(i),yvale2,ype,
     &                     ydpe,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,ve2,ve2dp,optke02(i),yvale3,ype,
     &                     ydpe,khi,klo,1,nkpts,1,0,0)
		wvalue=qw(1)*yval1+qw(2)*yval2+qw(3)*yval3
		evalue=pie(1,1)*yvale1+pie(1,2)*yvale2+pie(1,3)*yvale3-phi	
			 if (wvalue .gt. evalue)then
			 ye2(i,1) = dlog(cone1)+beta*wvalue
			 else
			 ye2(i,1) = dlog(cone1)+beta*evalue
			 endif

c	entrepreneur and mid productivity
	cone2 = (r+one)*xkcur+xprofit(2)-optke12(i)
			 call splint(xkpts,v0,v0dp,optke12(i),yval1,yp,
     &                     ydp,khi,klo,0,nkpts,1,0,0)
			 call splint(xkpts,v1,v1dp,optke12(i),yval2,yp,
     &                     ydp,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,v2,v2dp,optke12(i),yval3,yp,
     &                     ydp,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,ve0,ve0dp,optke12(i),yvale1,ype,
     &                     ydpe,khi,klo,0,nkpts,1,0,0)
			 call splint(xkpts,ve1,ve1dp,optke12(i),yvale2,ype,
     &                     ydpe,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,ve2,ve2dp,optke12(i),yvale3,ype,
     &                     ydpe,khi,klo,1,nkpts,1,0,0)
		wvalue=qw(1)*yval1+qw(2)*yval2+qw(3)*yval3
		evalue=pie(2,1)*yvale1+pie(2,2)*yvale2+pie(2,3)*yvale3-phi 
			 if (wvalue .gt. evalue)then
			 ye2(i,2) = dlog(cone2)+beta*wvalue
			 else
			 ye2(i,2) = dlog(cone2)+beta*evalue
			 endif
c	entrepreneur and high productivity
	cone3 = (r+one)*xkcur+xprofit(3)-optke22(i)
			 call splint(xkpts,v0,v0dp,optke22(i),yval1,yp,
     &                     ydp,khi,klo,0,nkpts,1,0,0)
			 call splint(xkpts,v1,v1dp,optke22(i),yval2,yp,
     &                     ydp,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,v2,v2dp,optke22(i),yval3,yp,
     &                     ydp,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,ve0,ve0dp,optke22(i),yvale1,ype,
     &                     ydpe,khi,klo,0,nkpts,1,0,0)
			 call splint(xkpts,ve1,ve1dp,optke22(i),yvale2,ype,
     &                     ydpe,khi,klo,1,nkpts,1,0,0)
			 call splint(xkpts,ve2,ve2dp,optke22(i),yvale3,ype,
     &                     ydpe,khi,klo,1,nkpts,1,0,0)
		wvalue=qw(1)*yval1+qw(2)*yval2+qw(3)*yval3
		evalue=pie(3,1)*yvale1+pie(3,2)*yvale2+pie(3,3)*yvale3-phi 
			 if (wvalue .gt. evalue)then
			 ye2(i,3) = dlog(cone3)+beta*wvalue
			 else
			 ye2(i,3) = dlog(cone3)+beta*evalue
			 endif

 42         continue
 43      continue

 2       vmxdif = 0.0D+00
	   vmxdife= 0.0D+00
	   cmxdif= 0.0D+00
         dmxdif = 0.0D+00
         emxdif = 0.0D+00
	   cmxdife = 0.0D+00
	   dmxdife = 0.0D+00
         emxdife = 0.0D+00

         do 15 i = 1,nkpts
            do 14 j = 1,3
               vdiff = dabs(y2(i,j)-y(i,j))
               if (vdiff .gt. vmxdif) vmxdif = vdiff
			 vdiffe = dabs(ye2(i,j)-ye(i,j))
               if (vdiffe .gt. vmxdife) vmxdife = vdiffe
 14         continue
		  cdiff = dabs(optk22(i)-optk2(i))
            if (cdiff .gt. cmxdif) cmxdif = cdiff
            ddiff = dabs(optk12(i)-optk1(i))
            if (ddiff .gt. dmxdif) dmxdif = ddiff
            ediff = dabs(optk02(i)-optk0(i))
            if (ediff .gt. emxdif) emxdif = ediff

		  cdiffe = dabs(optke22(i)-optke2(i))
            if (cdiffe .gt. cmxdife) cmxdife = cdiffe
		  ddiffe = dabs(optke12(i)-optke1(i))
            if (ddiffe .gt. dmxdife) dmxdife = ddiffe
            ediffe = dabs(optke02(i)-optke0(i))
            if (ediffe .gt. emxdife) emxdife = ediffe
 15      continue

         if (vmxdif .lt. 1.0D-06 .and. dmxdif .lt. 1.0D-08
     &       .and. emxdif .lt. 1.0D-08 .and. vmxdife .lt. 1.0D-06      
     &	   .and. dmxdife .lt. 1.0D-07
     &       .and. emxdife .lt. 1.0D-08 .and. 
     &		cmxdif .lt. 1.0D-08 .and. cmxdife .lt. 1.0D-08
     &		.and. kkk .gt. 2) goto 8

         if (multpl(kkk,1) .eq. 1) then
      write(6,888) kkk,mxiter,vmxdif,vmxdife,cmxdif,dmxdif,emxdif,
     &cmxdife,dmxdife,emxdife
 888        format(2i5,8f15.10)
         endif
 30   continue

 8    write(6,888) kkk,mxiter,vmxdif,vmxdife,cmxdif,dmxdif,emxdif,
     &cmxdife,dmxdife,emxdife

      open(3,file='entrep.out',form='formatted')
         do 928 i = 1,nkpts
            write(3,700)xkpts(i),y2(i,1),y2(i,2),y2(i,3),
     &	  ye2(i,1),ye2(i,2),ye2(i,3),
     &  optk02(i),optk12(i),optk22(i),optke02(i),optke12(i),optke22(i)
 700        format(13f18.10)
 928      continue
      close(3)   

	open(4,file='occupation.out',form='formatted')
      do 47 i = 1,nkpts
         write(4,454)xkpts(i),ioccup22(i),ioccup12(i),ioccup02(i),
     &ioccupe22(i),ioccupe12(i),ioccupe02(i)
 454     format(f18.10,6i15)
 47   continue
      close(4)   

	call prof2(r,w,xxk2,xxl2,xxy2,xxp2)

	open(21,file='profit2.out',form='formatted')	
		do 83 i = 1,nk2pts
          write(21,449) xk2pts(i),xxk2(1,i),xxk2(2,i),xxk2(3,i),
     &	xxl2(1,i),xxl2(2,i),xxl2(3,i),xxy2(1,i),xxy2(2,i),xxy2(3,i),
     &	xxp2(1,i),xxp2(2,i),xxp2(3,i)     
 449      format(13f18.10)
 83       continue
	close(21)

      do 77 i = 1,nk2pts         

		do 134 j=1,3
		  xprofit2(j)=xxp2(j,i)
 134	    continue

c	   ncalop=1

         xkcur = xk2pts(i)
c----------------------------------- worker with high productivity
	   neps = 2
	   nocc=1
         if (xkcur .le. xkcut3 .or. ncalop .eq. 1) then
            call calopt(xkcur,neps,nocc,xkp2,fval2,niter,
     &        1.0D-10,r,w,xkhigh,xklow,grad2,xkpts,xprofit2,iocc2)
         else
	      if (i .eq. 1) then
	         sk = xkcur
	      else
	         sk = optkg2(i-1)
	      endif
            call altopt(xkcur,neps,nocc,xkp2,fval2,niter,
     &        1.0D-10,sk,grad2,xkpts,r,w,xprofit2,iocc2)
         endif
         optkg2(i) = xkp2
         vgrid(i,3) = fval2
	   ioccupg2(i)=iocc2

         neps = 1
	   nocc=1
         if (xkcur .le. xkcut1 .or. ncalop .eq. 1) then
            call calopt(xkcur,neps,nocc,xkp1,fval1,niter,
     &        1.0D-10,r,w,xkhigh,xklow,grad1,xkpts,xprofit2,iocc1)
         else
	      if (i .eq. 1) then
	         sk = xkcur
	      else
	         sk = optkg1(i-1)
	      endif
            call altopt(xkcur,neps,nocc,xkp1,fval1,niter,
     &        1.0D-10,sk,grad1,xkpts,r,w,xprofit2,iocc1)
         endif
         optkg1(i) = xkp1
         vgrid(i,2) = fval1
	   ioccupg1(i)=iocc1

         neps = 0
	   nocc=1
         if (xkcur .le. xkcut2 .or. ncalop .eq. 1) then
            call calopt(xkcur,neps,nocc,xkp0,fval0,niter,
     &        1.0D-10,r,w,xkhigh,xklow,grad0,xkpts,xprofit2,iocc0)
         else
	      if (i .eq. 1) then
	         sk = xkcur
	      else
	         sk = optkg0(i-1)
	      endif
            call altopt(xkcur,neps,nocc,xkp0,fval0,niter,
     &        1.0D-10,sk,grad0,xkpts,r,w,xprofit2,iocc0)
         endif
         optkg0(i) = xkp0
         vgrid(i,1) = fval0
	   ioccupg0(i)=iocc0

c	--------------------------------------------------entrepreneur and has good shock
	   neps = 2
	   nocc=0
         if (xkcur .le. xkcut1 .or. ncalop .eq. 1) then
            call calopt(xkcur,neps,nocc,xkp3,fval3,niter,
     &         1.0D-10,r,w,xkhigh,xklow,grade3,xkpts,xprofit2,iocc3)
         else
	      if (i .eq. 1) then
	         sk = xkcur
	      else
	         sk = optkge2(i-1)
	      endif
            call altopt(xkcur,neps,nocc,xkp3,fval3,niter,
     &         1.0D-10,sk,grade3,xkpts,r,w,xprofit2,iocc3)
         endif
         optkge2(i) = xkp3
         vegrid(i,3) = fval3
	   ioccupge2(i)=iocc3	

c	entrepreneur and has mid shock
	   neps = 1
	   nocc=0
         if (xkcur .le. xkcut2 .or. ncalop .eq. 1) then
            call calopt(xkcur,neps,nocc,xkp4,fval4,niter,
     &         1.0D-10,r,w,xkhigh,xklow,grade4,xkpts,xprofit2,iocc4)
         else
	      if (i .eq. 1) then
	         sk = xkcur
	      else
	         sk = optkge1(i-1)
	      endif
            call altopt(xkcur,neps,nocc,xkp4,fval4,niter,
     &         1.0D-10,sk,grade4,xkpts,r,w,xprofit2,iocc4)
         endif
         optkge1(i) = xkp4
         vegrid(i,2) = fval4
	   ioccupge1(i)=iocc4	

         neps = 0
	   nocc=0
         if (xkcur .le. xkcut1 .or. ncalop .eq. 1) then
            call calopt(xkcur,neps,nocc,xkp5,fval5,niter,
     &        1.0D-10,r,w,xkhigh,xklow,grade5,xkpts,xprofit2,iocc5)
         else
	      if (i .eq. 1) then
	         sk = xkcur
	      else
	         sk = optkge0(i-1)
	      endif
            call altopt(xkcur,neps,nocc,xkp5,fval5,niter,
     &        1.0D-10,sk,grade5,xkpts,r,w,xprofit2,iocc5)
         endif
         optkge0(i) = xkp5
         vegrid(i,1) = fval5
	   ioccupge0(i)=iocc5

 77   continue

c      open(5,file='vgrid.sav',form='formatted')

c      do 71 i = 1,nk2pts

c         write(5,603) xk2pts(i),(vgrid(i,j),j=1,2),(vegrid(i,j),j=1,2)

c 603     format(5f18.10)

c 71   continue

c      close(5)

      open(9,file='entrepgrd.out',form='formatted')
      do 99 i = 1,nk2pts
         write(9,400)xk2pts(i),optkg2(i),optkg1(i),optkg0(i),
     &optkge2(i),optkge1(i),optkge0(i)
 400     format(7f18.10)
 99   continue
      close(9)

	open(7,file='occupgrd.out',form='formatted')
      do 95 i = 1,nk2pts
         write(7,402)xk2pts(i),ioccupg2(i),ioccupg1(i),ioccupg0(i),
     &ioccupge2(i),ioccupge1(i),ioccupge0(i)
 402     format(f18.10,7i15)
 95   continue
      close(7)

      write(6,7984)

 7984 format('computing invariant distribution')

      call calcf(optkg2,optkg1,optkg0,optkge2,optkge1,optkge0,
     &ioccupg2,ioccupg1,ioccupg0,ioccupge2,ioccupge1,ioccupge0,
     &xk2pts,xxk2,xxl2,xxy2,xxp2,xkmean,xkcmean,xlcmean,xkemean,ymean)

      rnew = acoef*alphac*(xkcmean**(alphac-one))*(xlcmean
     &  **(one-alphac))-delta

      wnew = acoef*(one-alphac)*(xkcmean**alphac)*
     &       (xlcmean**(-alphac))
	
	gdp=acoef*(xkcmean**alphac)*(xlcmean**(one-alphac))+ymean

	gdpc=acoef*(xkcmean**alphac)*(xlcmean**(one-alphac))

	xkgdp=xkmean/gdp

	xkgdpc=xkcmean/gdpc

	xkgdpe=xkemean/ymean

	aratio=xkemean/xkcmean

      write(6,509) xkmean,gdp,xkgdp,rnew,r,
     &			 xkemean,xkcmean,aratio,xlcmean,hfix,
     &			 ymean,gdpc,xkgdpe,xkgdpc
 509  format(5f18.10)

      open(16,file='entrep.f1',form='formatted')
      do 71 i = 1,nk2pts
         write(16,479) xk2pts(i),f1(i,1,1),f1(i,2,1),f1(i,3,1),
     &	f1(i,1,2),f1(i,2,2),f1(i,3,2)
 479     format(7f18.10)
 71      continue
      close(16)

	open(6,file='mean.out',form='formatted')
	   write(6,478)xkmean,gdp,xkgdp,rnew,r,
     &			 xkmean-xkcmean,xkcmean,aratio,xlcmean,hfix
 478     format(5f18.10)
      close(6)

      zeta2 = rnew-r

      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C

C   subroutines and functions

C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine splint(xa,ya,y2a,x,y,yp,ydp,khi,klo,kvals,n,
     &                  ndoy,ndoyp,ndoydp)
      implicit real*8 (a-h,o-z)
      dimension xa(n),ya(n),y2a(n)

      one = 1.0D+00
      if (kvals .eq. 0) then
         klo = 1
         khi = n
 1       if ((khi-klo) .gt. 1) then
            k = (khi+klo)/2
            if (xa(k) .gt. x) then
               khi = k
            else
               klo = k
            endif
            goto 1
         endif
      endif

      h = xa(khi)-xa(klo)
      a = (xa(khi)-x)/h
      b = (x-xa(klo))/h
      asq = a*a
      bsq = b*b

      if (ndoy .eq. 1)
     &   y = a*ya(klo)+b*ya(khi)+
     &       ((asq*a-a)*y2a(klo)+
     &       (bsq*b-b)*y2a(khi))*(h*h)/6.0D+00
      if (ndoyp .eq. 1)
     &   yp = (ya(khi)-ya(klo))/h-
     &        (3.0D+00*asq-one)/6.0D+00*h*y2a(klo)+
     &        (3.0D+00*bsq-one)/6.0D+00*h*y2a(khi)
      if (ndoydp .eq. 1)
     &   ydp = a*y2a(klo)+b*y2a(khi)      
      return
      end

      subroutine interp(x1,y,y0,y1,y2,y0dp,y1dp,y2dp)
      implicit real*8 (a-h,o-z)
      parameter (nkpts=150,n=nkpts)
      dimension x1(n),y(n,3),a(n),b(n),c(n),r(n),u(n),xx(n),
     &          yy(n),y1(n),y0(n),y2(n),y1dp(n),y0dp(n),y2dp(n)
      do 72 i = 1,n
         xx(i) = x1(i)
 72   continue
      a(1) = 0.0D+00
      a(n) = -1.0D+00
      b(1) = 1.0D+00
      b(n) = 1.05D+00
      c(1) = -1.05D+00
      c(n) = 0.0D+00
      r(1) = 0.0D+00
      r(n) = 0.0D+00
      do 81 k = 1,3
         do 30 i = 1,n
            yy(i) = y(i,k)
            if (k .eq. 1) y0(i) = yy(i)
            if (k .eq. 2) y1(i) = yy(i)
		  if (k .eq. 3) y2(i) = yy(i)
 30      continue

         do 21 i = 2,n-1
            a(i) = (xx(i)-xx(i-1))/6.0D+00
            b(i) = (xx(i+1)-xx(i-1))/3.0D+00
            c(i) = (xx(i+1)-xx(i))/6.0D+00
            r(i) = (yy(i+1)-yy(i))/(xx(i+1)-xx(i))-
     &             (yy(i)-yy(i-1))/(xx(i)-xx(i-1))
 21      continue
         call tridag(a,b,c,r,u,n)
         do 22 i = 1,n
            if (k .eq. 1) y0dp(i) = u(i)
            if (k .eq. 2) y1dp(i) = u(i)
		  if (k .eq. 3) y2dp(i) = u(i)
 22      continue
 81   continue

      return
      end

      subroutine tridag(a,b,c,r,u,n)
      implicit real*8 (a-h,o-z)
      parameter (toler=1.0D-12,nmax=5000)
      dimension gam(nmax),a(n),b(n),c(n),r(n),u(n)
      bet = b(1)
	u(1) = r(1)/bet
      do 11 j = 2,n
         gam(j) = c(j-1)/bet
         bet = b(j)-a(j)*gam(j)
         if (dabs(bet) .lt. toler) then
            write(6,100)
 100        format('failure in tridag')
            pause
         endif
         u(j) = (r(j)-a(j)*u(j-1))/bet
 11   continue

      do 12 j = n-1,1,-1
         u(j) = u(j)-gam(j+1)*u(j+1)
 12   continue
      return
      end

      integer function multpl(i,m)
      implicit real*8 (a-h,o-z)
      parameter (toler=1.0D-12)
      remain = dble(real(i))/dble(real(m))-i/m
      if (remain .le. toler) then
         multpl = 1
      else
         multpl = 0
      endif

      return
      end

      subroutine calopt(xk,neps,nocc,ynew,value,niter,toler,
     &                  r,w,ykhigh,yklow,grad,xkpts,xp,iocc)
c	xp is profit of entrepreneur produced by subroutine 'profit1m'
      implicit real*8 (a-h,o-z)
      parameter (mxiter=200,nkpts=150,hfix=0.3271D+00,delta=0.06D+00,
     &phi=0.5D+00,beta=0.9D+00)
      common/cprob/piw(3,3),pie(3,3),qw(3),qe(3)
	common/cwages/ylo,ymi,yhi
c	common/cparam/beta,z1
c	common/cres/xkgdpc,aratio
      common/cvs/v0(nkpts),v0dp(nkpts),
     &           v1(nkpts),v1dp(nkpts),ve0(nkpts),ve0dp(nkpts),
     &           ve1(nkpts),ve1dp(nkpts),
     &           v2(nkpts),v2dp(nkpts),
     &           ve2(nkpts),ve2dp(nkpts)
      dimension xkpts(nkpts),xp(3)
      one = 1.0D+00
      niter = 0

      if (nocc .eq. 1 .and. neps .eq. 2) then
         y = r*xk+w*hfix*yhi
      elseif (nocc .eq. 1 .and. neps .eq. 1) then
	   y = r*xk+w*hfix*ymi
      elseif (nocc .eq. 1 .and. neps .eq. 0) then
         y = r*xk+w*hfix*ylo
	elseif (nocc .eq. 0 .and. neps .eq. 0) then
	   y=xp(1)+r*xk
	elseif (nocc .eq. 0 .and. neps .eq. 1) then
	   y=xp(2)+r*xk
	else
	   y=xp(3)+r*xk
      endif

      if (nocc .eq. 1 .and. neps .eq. 0) then
         piw1 = piw(1,1)
	   piw2= piw(1,2)
      elseif (nocc .eq. 1 .and. neps .eq. 1)then
         piw1 = piw(2,1)
	   piw2 = piw(2,2)
	elseif (nocc .eq. 1 .and. neps .eq. 2)then
         piw1 = piw(3,1)
	   piw2 = piw(3,2)
	elseif (nocc .eq. 0 .and. neps .eq. 0)then
         pie1 = pie(1,1)
	   pie2 = pie(1,2)
	elseif (nocc .eq. 0 .and. neps .eq. 1)then
         pie1 = pie(2,1)
	   pie2 = pie(2,2)
	else
         pie1 = pie(3,1)
	   pie2 = pie(3,2)
      endif


      xklow = yklow
      xkhigh = ykhigh

      do 10 i = 1,mxiter
         xkcur = (xkhigh+xklow)/2.0D+00
         x = xkcur-(one)*xk
         if (x .ge. y) then
            xkhigh = xkcur
            goto 10
         endif

c	adds comparision of value functions because of occupational choice
         c = y-x
         upc = -one/c

         	
		call splint(xkpts,v0,v0dp,xkcur,yval1,yp,ydp,khi,klo,
     &            0,nkpts,1,0,0)
		call splint(xkpts,v1,v1dp,xkcur,yval2,yp,ydp,khi,klo,
     &            1,nkpts,1,0,0)
		call splint(xkpts,v2,v2dp,xkcur,yval3,yp,ydp,khi,klo,
     &            1,nkpts,1,0,0)	

		call splint(xkpts,ve0,ve0dp,xkcur,yvale1,ype,ydpe,khi,klo,
     &            0,nkpts,1,0,0) 
		call splint(xkpts,ve1,ve1dp,xkcur,yvale2,ype,ydpe,khi,klo,
     &            1,nkpts,1,0,0)
		call splint(xkpts,ve2,ve2dp,xkcur,yvale3,ype,ydpe,khi,klo,
     &            1,nkpts,1,0,0)
		   
		if (nocc .eq. 1) then	 
		wvalue=piw1*yval1+piw2*yval2+(one-piw1-piw2)*yval3
		evalue=qe(1)*yvale1+qe(2)*yvale2+(one-qe(1)-qe(2))*yvale3-phi	

		if (wvalue .gt. evalue)then
         call splint(xkpts,v0,v0dp,xkcur,yval,yp1,ydp,khi,klo,
     &               0,nkpts,0,1,0)
	   call splint(xkpts,v1,v1dp,xkcur,yval,yp2,ydp,khi,klo,
     &               1,nkpts,0,1,0)   
	   call splint(xkpts,v2,v2dp,xkcur,yval,yp3,ydp,khi,klo,
     &               1,nkpts,0,1,0)

	    grad = upc+beta*(piw1*yp1+piw2*yp2+(one-piw1-piw2)*yp3)
	   else
         call splint(xkpts,ve0,ve0dp,xkcur,yvale,ype1,ydpe,khi,klo,
     &               0,nkpts,0,1,0)
	   call splint(xkpts,ve1,ve1dp,xkcur,yvale,ype2,ydpe,khi,klo,
     &               1,nkpts,0,1,0)
	   call splint(xkpts,ve2,ve2dp,xkcur,yvale,ype3,ydpe,khi,klo,
     &               1,nkpts,0,1,0)
         grad = upc+beta*(qe(1)*ype1+qe(2)*ype2+(one-qe(1)-qe(2))*ype3)
	   endif

		else	 
		wvalue=qw(1)*yval1+qw(2)*yval2+(one-qw(1)-qw(2))*yval3
		evalue=pie1*yvale1+pie2*yvale2+(one-pie1-pie2)*yvale3-phi

		if (wvalue .gt. evalue)then
		call splint(xkpts,v0,v0dp,xkcur,yval,yp1,ydp,khi,klo,
     &               0,nkpts,0,1,0)
		call splint(xkpts,v1,v1dp,xkcur,yval,yp2,ydp,khi,klo,
     &               1,nkpts,0,1,0)   
	    call splint(xkpts,v2,v2dp,xkcur,yval,yp3,ydp,khi,klo,
     &               1,nkpts,0,1,0)
	    grad = upc+beta*(qw(1)*yp1+qw(2)*yp2+(one-qw(1)-qw(2))*yp3)
		else
		call splint(xkpts,ve0,ve0dp,xkcur,yvale,ype1,ydpe,khi,klo,
     &               0,nkpts,0,1,0)
		call splint(xkpts,ve1,ve1dp,xkcur,yvale,ype2,ydpe,khi,klo,
     &               1,nkpts,0,1,0)
	    call splint(xkpts,ve2,ve2dp,xkcur,yvale,ype3,ydpe,khi,klo,
     &               1,nkpts,0,1,0)
		grad = upc+beta*(pie1*ype1+pie2*ype2+(one-pie1-pie2)*ype3)
		endif

		endif		  
	 	   
	
         if (grad .ge. 0.0D+00) then
            xklow = xkcur
         else
            xkhigh = xkcur
         endif

         if ((xkhigh-xklow) .le. toler) then
            niter = i
            goto 9
         endif

 10   continue

 9		xkcur = (xkhigh+xklow)/2.0D+00
		x = xkcur-(one)*xk
		c = y-x

		if (c .le. 0.0D+00) then
		 write(6,787) c
 787		 format('failure in calopt: negative c',f15.10)
		 pause
		endif

		utilc = dlog(c)  
	   call splint(xkpts,v0,v0dp,xkcur,yval1,yp,ydp,khi,klo,
     &            0,nkpts,1,0,0)
	   call splint(xkpts,v1,v1dp,xkcur,yval2,yp,ydp,khi,klo,
     &            1,nkpts,1,0,0)
	   call splint(xkpts,v2,v2dp,xkcur,yval3,yp,ydp,khi,klo,
     &            1,nkpts,1,0,0)
	   call splint(xkpts,ve0,ve0dp,xkcur,yvale1,ype,ydpe,khi,klo,
     &            0,nkpts,1,0,0)
	   call splint(xkpts,ve1,ve1dp,xkcur,yvale2,ype,ydpe,khi,klo,
     &            1,nkpts,1,0,0)     	
	   call splint(xkpts,ve2,ve2dp,xkcur,yvale3,ype,ydpe,khi,klo,
     &            1,nkpts,1,0,0)

		if (nocc .eq. 1) then
		wvalue=piw1*yval1+piw2*yval2+(one-piw1-piw2)*yval3
		evalue=qe(1)*yvale1+qe(2)*yvale2+(one-qe(1)-qe(2))*yvale3-phi  
		else
		wvalue=qw(1)*yval1+qw(2)*yval2+(one-qw(1)-qw(2))*yval3
		evalue=pie1*yvale1+pie2*yvale2+(one-pie1-pie2)*yvale3-phi  
		endif
	
		if (wvalue .gt. evalue)then
		value = utilc+beta*wvalue
		iocc=1
		else
		value = utilc+beta*evalue
		iocc=0
		endif

		ynew = xkcur

      return
      end

      subroutine altopt(xk,neps,nocc,ynew,value,niter,toler,
     &                  startk,grad,xkpts,r,w,xp,iocc)
      implicit real*8 (a-h,o-z)
      parameter (mxiter=200,hfix=0.3271D+00,delta=0.06D+00,nkpts=150,
     &phi=0.5D+00,beta=0.9D+00)
      common/cprob/piw(3,3),pie(3,3),qw(3),qe(3)
	common/cwages/ylo,ymi,yhi
c	common/cparam/beta,z1
c	common/cres/xkgdpc,aratio
      common/cvs/v0(nkpts),v0dp(nkpts),v2(nkpts),v2dp(nkpts),
     &           v1(nkpts),v1dp(nkpts),ve0(nkpts),ve0dp(nkpts),
     &           ve1(nkpts),ve1dp(nkpts),ve2(nkpts),ve2dp(nkpts)
      dimension xkpts(nkpts),xp(3)

      one = 1.0D+00
      niter = 0

	if (nocc .eq. 1 .and. neps .eq. 2) then
         y = r*xk+w*hfix*yhi
      elseif (nocc .eq. 1 .and. neps .eq. 1) then
	   y = r*xk+w*hfix*ymi
      elseif (nocc .eq. 1 .and. neps .eq. 0) then
         y = r*xk+w*hfix*ylo
	elseif (nocc .eq. 0 .and. neps .eq. 0) then
	   y=xp(1)+r*xk
	elseif (nocc .eq. 0 .and. neps .eq. 1) then
	   y=xp(2)+r*xk
	else
	   y=xp(3)+r*xk
      endif

      if (nocc .eq. 1 .and. neps .eq. 0) then
         piw1 = piw(1,1)
	   piw2= piw(1,2)
      elseif (nocc .eq. 1 .and. neps .eq. 1)then
         piw1 = piw(2,1)
	   piw2 = piw(2,2)
	elseif (nocc .eq. 1 .and. neps .eq. 2)then
         piw1 = piw(3,1)
	   piw2 = piw(3,2)
	elseif (nocc .eq. 0 .and. neps .eq. 0)then
         pie1 = pie(1,1)
	   pie2 = pie(1,2)
	elseif (nocc .eq. 0 .and. neps .eq. 1)then
         pie1 = pie(2,1)
	   pie2 = pie(2,2)
	else
         pie1 = pie(3,1)
	   pie2 = pie(3,2)
      endif

      xk2p = startk

      do 10 i = 1,mxiter
         xk2 = xk2p
         x = xk2-(one)*xk
         c = y-x

         upc = -one/c
         udpc = -one/(c*c)
		call splint(xkpts,v0,v0dp,xk2,yval1,yp,ydp,khi,klo,
     &            0,nkpts,1,0,0)
		call splint(xkpts,v1,v1dp,xk2,yval2,yp,ydp,khi,klo,
     &            1,nkpts,1,0,0)
		call splint(xkpts,v2,v2dp,xk2,yval3,yp,ydp,khi,klo,
     &            1,nkpts,1,0,0)
		call splint(xkpts,ve0,ve0dp,xk2,yvale1,ype,ydpe,khi,klo,
     &            0,nkpts,1,0,0)
		call splint(xkpts,ve1,ve1dp,xk2,yvale2,ype,ydpe,khi,klo,
     &            1,nkpts,1,0,0)       
		call splint(xkpts,ve2,ve2dp,xk2,yvale3,ype,ydpe,khi,klo,
     &            1,nkpts,1,0,0) 
		

		if (nocc .eq. 1) then	 
		wvalue=piw1*yval1+piw2*yval2+(one-piw1-piw2)*yval3
		evalue=qe(1)*yvale1+qe(2)*yvale2+(one-qe(1)-qe(2))*yvale2-phi	 
	 	   
	   if (wvalue .gt. evalue)then
         call splint(xkpts,v0,v0dp,xk2,yval,yp1,ydp1,khi,klo,
     &               0,nkpts,0,1,1)   
	   call splint(xkpts,v1,v1dp,xk2,yval,yp2,ydp2,khi,klo,
     &               1,nkpts,0,1,1)
	   call splint(xkpts,v2,v2dp,xk2,yval,yp3,ydp3,khi,klo,
     &               1,nkpts,0,1,1)
	    grad = upc+beta*(piw1*yp1+piw2*yp2+(one-piw1-piw2)*yp3)
		hess = udpc+beta*(piw1*ydp1+piw2*ydp2+(one-piw1-piw2)*ydp3)
	   else
         call splint(xkpts,ve0,ve0dp,xk2,yvale,ype1,ydpe1,khi,klo,
     &               0,nkpts,0,1,1)
	   call splint(xkpts,ve1,ve1dp,xk2,yvale,ype2,ydpe2,khi,klo,
     &               1,nkpts,0,1,1)
	   call splint(xkpts,ve2,ve2dp,xk2,yvale,ype3,ydpe3,khi,klo,
     &               1,nkpts,0,1,1)
         grad = upc+beta*(qe(1)*ype1+qe(2)*ype2+(one-qe(1)-qe(2))*ype3)
	   hess = udpc+beta*(qe(1)*ydpe1+qe(2)*ydpe2
     &			+(one-qe(1)-qe(2))*ydpe3)
	   endif
		
		else	 
		wvalue=qw(1)*yval1+qw(2)*yval2+(one-qw(1)-qw(2))*yval3
		evalue=pie1*yvale1+pie2*yvale2+(one-pie1-pie2)*yvale3-phi
		if (wvalue .gt. evalue)then
		call splint(xkpts,v0,v0dp,xk2,yval,yp1,ydp1,khi,klo,
     &               0,nkpts,0,1,1)   
		call splint(xkpts,v1,v1dp,xk2,yval,yp2,ydp2,khi,klo,
     &               1,nkpts,0,1,1)
	    call splint(xkpts,v2,v2dp,xk2,yval,yp3,ydp3,khi,klo,
     &               1,nkpts,0,1,1)
	    grad = upc+beta*(qw(1)*yp1+qw(2)*yp2+(one-qw(1)-qw(2))*yp3)
		hess = udpc+beta*(qw(1)*ydp1+qw(2)*ydp2+(one-qw(1)-qw(2))*ydp3)
		else
		call splint(xkpts,ve0,ve0dp,xk2,yvale,ype1,ydpe1,khi,klo,
     &               0,nkpts,0,1,1)
		call splint(xkpts,ve1,ve1dp,xk2,yvale,ype2,ydpe2,khi,klo,
     &               1,nkpts,0,1,1)
	    call splint(xkpts,ve2,ve2dp,xk2,yvale,ype3,ydpe3,khi,klo,
     &               1,nkpts,0,1,1)
		grad = upc+beta*(pie1*ype1+pie2*ype2+(one-pie1-pie2)*ype3)
		hess = udpc+beta*(pie1*ydpe1+pie2*ydpe2+(one-pie1-pie2)*ydpe3)
		endif

		endif

         xk2p = xk2-grad/hess

         if (dabs(xk2p-xk2) .le. toler) then
            niter = i
            goto 9
         endif

 10   continue

 9    if (c .gt. 0.0D+00) then
         utilc = dlog(c)
	   call splint(xkpts,v0,v0dp,xk2,yval1,yp,ydp,khi,klo,
     &            0,nkpts,1,0,0)	
	   call splint(xkpts,v1,v1dp,xk2,yval2,yp,ydp,khi,klo,
     &            1,nkpts,1,0,0)
	   call splint(xkpts,v2,v2dp,xk2,yval3,yp,ydp,khi,klo,
     &            1,nkpts,1,0,0)
	   call splint(xkpts,ve0,ve0dp,xk2,yvale1,ype,ydpe,khi,klo,
     &            0,nkpts,1,0,0)
	   call splint(xkpts,ve1,ve1dp,xk2,yvale2,ype,ydpe,khi,klo,
     &            1,nkpts,1,0,0)
	   call splint(xkpts,ve2,ve2dp,xk2,yvale3,ype,ydpe,khi,klo,
     &            1,nkpts,1,0,0)

		if (nocc .eq. 1) then
		wvalue=piw1*yval1+piw2*yval2+(one-piw1-piw2)*yval3
		evalue=qe(1)*yvale1+qe(2)*yvale2+(one-qe(1)-qe(2))*yvale3-phi  
		else
		wvalue=qw(1)*yval1+qw(2)*yval2+(one-qw(1)-qw(2))*yval3
		evalue=pie1*yvale1+pie2*yvale2+(one-pie1-pie2)*yvale3-phi  
		endif
	
		if (wvalue .gt. evalue)then
		value = utilc+beta*wvalue
		iocc=0
		else
		value = utilc+beta*evalue
		iocc=1
		endif
         ynew = xk2         
      else
         write(6,123) xk,neps,c
 123     format('failure in altopt: negative c',f8.3,i3,f12.8)
         pause
      endif

      return
      end

      subroutine hunt(xx,n,x,jlo)
      implicit real*8 (a-h,o-z)
      dimension xx(n)
      logical ascnd

      ascnd = xx(n) .ge. xx(1)
      if (jlo .le. 0 .or. jlo .gt. n) then
         jlo = 0
         jhi = n+1
         goto 3
      endif

      inc = 1     

      if (x .ge. xx(jlo) .eqv. ascnd) then
 1       jhi = jlo+inc
         if (jhi .gt. n) then
            jhi = n+1
         elseif (x .ge. xx(jhi) .eqv. ascnd) then
            jlo = jhi
            inc = inc+inc
            goto 1
         endif
      else
         jhi = jlo
 2       jlo = jhi-inc
         if (jlo .lt. 1) then
            jlo = 0
         elseif (x .lt. xx(jlo) .eqv. ascnd) then
            jhi = jlo
            inc = inc+inc
            goto 2
         endif
      endif
 3    if (jhi-jlo .eq. 1) then
         if (x .eq. xx(n)) jlo = n-1
         if (x .eq. xx(1)) jlo = 1
         return
      endif

      jm = (jhi+jlo)/2

      if (x .ge. xx(jm) .eqv. ascnd) then
         jlo = jm
      else
         jhi = jm
      endif
      goto 3
      end

      subroutine calcf(optkg2,optkg1,optkg0,optkge2,optkge1,optkge0,
     &ioccupg2,ioccupg1,ioccupg0,ioccupge2,ioccupge1,ioccupge0,
     &xk2pts,xxk,xxl,xxy,xxp,xkmean,xkcmean,xlcmean,xkemean,ymean)
      implicit real*8 (a-h,o-z)
      parameter (nk2pts=5000,mxiter=10000000,
     &hfix=0.3271D+00)
      common/cfs/f0(nk2pts,3,2),f1(nk2pts,3,2)
	common/cwages/ylo,ymi,yhi
      common/cprob/piw(3,3),pie(3,3),qw(3),qe(3)
      dimension xk2pts(nk2pts),optkg1(nk2pts),optkg2(nk2pts),
     &optkg0(nk2pts),prob0(nk2pts),prob1(nk2pts),prob2(nk2pts),
     &jlos0(nk2pts),jlos1(nk2pts),jlos2(nk2pts),
     &wghts1(nk2pts),wghts0(nk2pts),wghts3(nk2pts),
     &optkge2(nk2pts),optkge1(nk2pts),optkge0(nk2pts),
     &jlose2(nk2pts),jlose1(nk2pts),jlose0(nk2pts),
     &probe0(nk2pts),probe1(nk2pts),probe2(nk2pts),
     &wghtse2(nk2pts),wghtse1(nk2pts),wghtse0(nk2pts),
     &xxk(3,nk2pts),xxl(3,nk2pts),xxp(3,nk2pts),xxy(3,nk2pts),
     &jlos(nk2pts,3,2),wghts(nk2pts,3,2),
     &ioccupg1(nk2pts),ioccupge1(nk2pts),
     &ioccupg0(nk2pts),ioccupge0(nk2pts),
     &ioccupg2(nk2pts),ioccupge2(nk2pts),
     &ioccupg(3,nk2pts),ioccupge(3,nk2pts),
     &    xpt(nk2pts),ypt(nk2pts),ydp(nk2pts),
     &    aa(nk2pts),bb(nk2pts),cc(nk2pts),rr(nk2pts)

      one = 1.0D+00
	xkinc = xk2pts(2)-xk2pts(1)
	
      do 81 ii = 1,nk2pts
         prob0(ii) = piw(1,1)
         prob1(ii) = piw(2,1)
	   prob2(ii) = piw(3,1)
	   probe0(ii) = pie(1,1)
         probe1(ii) = pie(2,1)
	   probe2(ii) = pie(3,1)

	ioccupg(1,ii)=ioccupg0(ii)
	ioccupg(2,ii)=ioccupg1(ii)
	ioccupg(3,ii)=ioccupg2(ii)
	ioccupge(1,ii)=ioccupge0(ii)
	ioccupge(2,ii)=ioccupge1(ii)
	ioccupge(3,ii)=ioccupge2(ii)

	   call hunt(xk2pts,nk2pts,optkg2(ii),jlo2)
         call hunt(xk2pts,nk2pts,optkg1(ii),jlo1)
         call hunt(xk2pts,nk2pts,optkg0(ii),jlo0)
	   call hunt(xk2pts,nk2pts,optkge2(ii),jloe2)
	   call hunt(xk2pts,nk2pts,optkge1(ii),jloe1)
         call hunt(xk2pts,nk2pts,optkge0(ii),jloe0)
		
	   if (jlo2 .lt. 1) then
            jlo2 = 1
            wght2 = 1.0D+00
         elseif (jlo2 .ge. nk2pts) then
            jlo2 = nk2pts-1
            wght2 = 0.0D+00
         else   
            wght2 = one-(optkg2(ii)-xk2pts(jlo2))/xkinc
         endif
         wghts(ii,3,1) = wght2
	   jlos(ii,3,1) = jlo2

         if (jlo1 .lt. 1) then
            jlo1 = 1
            wght1 = 1.0D+00
         elseif (jlo1 .ge. nk2pts) then
            jlo1 = nk2pts-1
            wght1 = 0.0D+00
         else   
            wght1 = one-(optkg1(ii)-xk2pts(jlo1))/xkinc
         endif
         wghts(ii,2,1) = wght1
	   jlos(ii,2,1) = jlo1

         if (jlo0 .lt. 1) then
            jlo0 = 1
            wght0 = 1.0D+00
         elseif (jlo0 .ge. nk2pts) then
            jlo0 = nk2pts-1
            wght0 = 0.0D+00
         else
            wght0 = one-(optkg0(ii)-xk2pts(jlo0))/xkinc
c     &               (xk2pts(jlo0+1)-xk2pts(jlo0))
         endif
	   wghts(ii,1,1) = wght0
	   jlos(ii,1,1) = jlo0	  

	   if (jloe2 .lt. 1) then
            jloe2 = 1
            wghte2 = 1.0D+00
         elseif (jloe2 .ge. nk2pts) then
            jloe2 = nk2pts-1
            wghte2 = 0.0D+00
         else   
            wghte2 = one-(optkge2(ii)-xk2pts(jloe2))/xkinc
         endif
         wghts(ii,3,2) = wghte2
	   jlos(ii,3,2) = jloe2

	   if (jloe1 .lt. 1) then
            jloe1 = 1
            wghte1 = 1.0D+00
         elseif (jloe1 .ge. nk2pts) then
            jloe1 = nk2pts-1
            wghte1 = 0.0D+00
         else   
            wghte1 = one-(optkge1(ii)-xk2pts(jloe1))/xkinc
         endif
         wghts(ii,2,2) = wghte1
	   jlos(ii,2,2) = jloe1

         if (jloe0 .lt. 1) then
            jloe0 = 1
            wghte0 = 1.0D+00
         elseif (jloe0 .ge. nk2pts) then
            jloe0 = nk2pts-1
            wghte0 = 0.0D+00
         else
            wghte0 = one-(optkge0(ii)-xk2pts(jloe0))/xkinc
         endif
	   wghts(ii,1,2) = wghte0
	   jlos(ii,1,2) = jloe0
 81   continue
	
      do 70 k = 1,mxiter

         do 12 i = 1,nk2pts
            do 11 j = 1,3
			do 10 kk=1,2
               f0(i,j,kk) = f1(i,j,kk)
               f1(i,j,kk) = 0.0D+00
 10			continue
 11         continue
 12      continue

		do 62 iz=1,3

		do 68 ii = 1,nk2pts

		if(ioccupge(iz,ii) .eq. 0)then

		do 768 i2=1,3
		f1(jlos(ii,iz,2),i2,2)=pie(iz,i2)*f0(ii,iz,2)*wghts(ii,i2,2)
     &							+f1(jlos(ii,iz,2),i2,2)
	    f1(jlos(ii,iz,2)+1,i2,2)=pie(iz,i2)*f0(ii,iz,2)*(one-
     &						wghts(ii,i2,2))+f1(jlos(ii,iz,2)+1,i2,2)
 768		continue

		else

		do 767 i2=1,3
		f1(jlos(ii,iz,2),i2,1)=qw(i2)*f0(ii,iz,2)*wghts(ii,i2,1)
     &							+f1(jlos(ii,iz,2),i2,1)
	    f1(jlos(ii,iz,2)+1,i2,1)=qw(i2)*f0(ii,iz,2)*(one-
     &						wghts(ii,i2,1))+f1(jlos(ii,iz,2)+1,i2,1)
 767		continue

		endif

 68		continue

		do 66 ii = 1,nk2pts
		if(ioccupg(iz,ii) .eq. 0)then

		do 758 i2=1,3
		f1(jlos(ii,iz,1),i2,2)=qe(i2)*f0(ii,iz,1)*wghts(ii,i2,2)
     &							+f1(jlos(ii,iz,1),i2,2)
	    f1(jlos(ii,iz,1)+1,i2,2)=qe(i2)*f0(ii,iz,1)*(one-
     &						wghts(ii,i2,2))+f1(jlos(ii,iz,1)+1,i2,2)
 758		continue

		else

		do 753 i2=1,3
		f1(jlos(ii,iz,1),i2,1)=piw(iz,i2)*f0(ii,iz,1)*wghts(ii,i2,1)
     &							+f1(jlos(ii,iz,1),i2,1)
	    f1(jlos(ii,iz,1)+1,i2,1)=piw(iz,i2)*f0(ii,iz,1)*(one-
     &						wghts(ii,i2,1))+f1(jlos(ii,iz,1)+1,i2,1)
 753		continue
		endif

 66		continue

 62		continue

         fmxdif = 0.0D+00
         fsum = 0.0D+00
         do 98 i = 1,nk2pts
            do 97 j = 1,3
			do 96 kk=1,2
               fsum = fsum+f1(i,j,kk)
               fdiff = dabs(f1(i,j,kk)-f0(i,j,kk))
               if (fdiff .gt. fmxdif) fmxdif = fdiff
 96			continue
 97         continue
 98      continue

         if (dabs(fsum-one) .gt. 1.0D-06) then
            write(6,543) k,fsum
 543        format('does not integrate to 1: ',i10,f18.10)
            read *
         endif

         if (multpl(k,10000) .eq. 1) then
            write(6,703) k,mxiter,fmxdif
 703        format(2i10,f18.10)
         endif

         if (fmxdif .le. 1.0D-08) goto 1
 70   continue

 1    xkmean = 0.0D+00
      do 32 i = 1,nk2pts
      xkmean=xkmean+xk2pts(i)*(f1(i,1,1)+f1(i,2,1)+f1(i,3,1)+
     &			f1(i,1,2)+f1(i,2,2)+f1(i,3,2))
 32   continue
	
	xkemean = 0.0D+00
      do 35 i = 1,nk2pts
      xkemean=xkemean+xxk(1,i)*f1(i,1,2)+xxk(2,i)*f1(i,2,2)+
     &				xxk(3,i)*f1(i,3,2)
 35   continue	
	xkcmean=xkmean-xkemean

	ymean = 0.0D+00
      do 33 i = 1,nk2pts
      ymean=ymean+xxy(1,i)*f1(i,1,2)+xxy(2,i)*f1(i,2,2)+
     &				xxy(3,i)*f1(i,3,2)
 33   continue	


	aggh = 0.0D+00
      do 37 i = 1,nk2pts
      aggh=aggh+hfix*(ylo*f1(i,1,1)+ymi*f1(i,2,1)+yhi*f1(i,3,1))
 37   continue

	agghr = 0.0D+00
      do 30 i = 1,nk2pts
      agghr=agghr+hfix*(f1(i,1,1)+f1(i,2,1)+f1(i,3,1))
 30   continue

	xfirm = 0.0D+00
      do 48 i = 1,nk2pts
      xfirm=xfirm+f1(i,1,2)+f1(i,2,2)+f1(i,3,2)
 48   continue

	xworker = 0.0D+00
      do 49 i = 1,nk2pts
      xworker=xworker+f1(i,1,1)+f1(i,2,1)+f1(i,3,1)
 49   continue

	xlemean = 0.0D+00
      do 34 i = 1,nk2pts
      xlemean=xlemean+xxl(1,i)*f1(i,1,2)+xxl(2,i)*f1(i,2,2)+
     &				xxl(3,i)*f1(i,3,2)
 34   continue
 	
	xlcmean=aggh-xlemean

	exitf=0.0D+00
	entryf=0.0D+00
	do 44 i=1,nk2pts
		do 45 j=1,3
		if(ioccupge(j,i) .eq. 1)then
		exitf=exitf+f1(i,j,2)
		endif
		if(ioccupg(j,i) .eq. 0)then
		entryf=entryf+f1(i,j,1)
		endif
45		continue
44	continue
	
c	compute the gini coefficient	
	sumx = 0.0D+00
      sumy = 0.0D+00
      xpt(1) = 0.0D+00
      ypt(1) = 0.0D+00
      nhpts = 1
      do 89 j = 2,nk2pts
         sumx=sumx+f1(j-1,1,1)+f1(j-1,2,1)+f1(j-1,3,1)+
     &		f1(j-1,1,2)+f1(j-1,2,2)+f1(j-1,3,2)
         sumy=sumy+(f1(j-1,1,1)+f1(j-1,2,1)+f1(j-1,3,1)+
     &		f1(j-1,1,2)+f1(j-1,2,2)+f1(j-1,3,2))*xk2pts(j-1)/xkmean
         if (f1(j-1,1,1)+f1(j-1,2,1)+f1(j-1,3,1)+
     &		f1(j-1,1,2)+f1(j-1,2,2)+f1(j-1,3,2) .gt. 1.0D-10) then
            nhpts = nhpts+1
            xpt(nhpts) = sumx
            ypt(nhpts) = sumy
         endif
 89   continue

      call calcsd(xpt,ypt,ydp,aa,bb,cc,rr,nhpts,-one,-one)
      call integr(xpt,ypt,ydp,nhpts,xint)
      gini = one-2.0D+00*xint
	
	write(6,399) xkmean,aggh,agghr,gini,
     &			 exitf,entryf,xfirm,xworker
      open(2,file='entrep1.mom',form='formatted')
	write(2,399) xkmean,aggh,agghr,gini,
     &			 exitf,entryf,xfirm,xworker
 399  format(4f18.10)
      close(2)

      return
      end

	
	subroutine integr(xa,ya,y2a,n,xint)
      implicit real*8 (a-h,o-z)
      dimension xa(n),ya(n),y2a(n)
      one = 1.0D+00
      one24 = one/24.0D+00
      sum = 0.0D+00
      do 10 i = 1,n-1
         diff = xa(i+1)-xa(i)
         sum = sum+diff*( (ya(i)+ya(i+1))/2.0D+00 -
     &         one24*diff*diff*(y2a(i)+y2a(i+1)))
 10   continue

      xint = sum

      return
      end

      subroutine calcsd(x,y,ydp,a,b,c,r,n,an,c1)
      implicit real*8 (a-h,o-z)
      dimension x(n),y(n),ydp(n),a(n),b(n),c(n),r(n)
      a(1) = 0.0D+00
      a(n) = an
      b(1) = 1.0D+00
      b(n) = 1.0D+00
      c(1) = c1
      c(n) = 0.0D+00
      r(1) = 0.0D+00
      r(n) = 0.0D+00

      do 21 i = 2,n-1
         a(i) = (x(i)-x(i-1))/6.0D+00
         b(i) = (x(i+1)-x(i-1))/3.0D+00
         c(i) = (x(i+1)-x(i))/6.0D+00
         r(i) = (y(i+1)-y(i))/(x(i+1)-x(i))-
     &          (y(i)-y(i-1))/(x(i)-x(i-1))
 21   continue

      call tridag(a,b,c,r,ydp,n)

      return
      end

c------------------------- compute the demands and profit in the entrepreneur sector
	subroutine prof1(r,w,xk,xl,xy,xp)
	implicit real*8 (a-h,o-z)
c	common/cparam/beta,z1
c	common/cres/xkgdp,aratio
      parameter (nkpts=150,
     &tol=1.0D-08,nc=3,ns1=3,ns2=nkpts,depreciation=0.06D+00,
     &alpha=0.366D+00,xbeta=0.88D+00-alpha,theta=0.6D+00)
c     ,z1=1.0D+00)
      dimension xk(ns1,ns2),xl(ns1,ns2),xy(ns1,ns2),xp(ns1,ns2),
     &	a(nkpts),z(ns1),delta(ns1),hw(ns1),he(ns1)
      one=1.0D+00
	zero=0.0D+00

	ri=r+depreciation
	
c	grid
	open(2,file='xkpts150.in',form='formatted')
         do 97 i = 1,nkpts
	      read(2,390) a(i)
 390        format(f18.10)
 97	   continue 
	close(2)

c	z(1)=2.45D+00
	z(1)=1.0D+00
	z(2)=10.0D+00**(one-alpha-xbeta)*z(1)
	z(3)=100.0D+00**(one-alpha-xbeta)*z(1)	
	
c	profit function
		do 114 iz=1,ns1
c		  xkinit=(z(iz)*((alpha/ri)**(one-beta))*((beta/w)**beta))**
c     &		(one/(one-alpha-beta))
c
c		  pro=z(iz)*((beta*ri/(alpha*w))**beta)*(xkinit**
c     &		(alpha+beta))-((one+(beta/alpha))*ri*xkinit)
		xkinit=(z(iz)*((alpha/ri)**(one-xbeta))*((xbeta/w)**xbeta))
     &				**(one/(one-alpha-xbeta))  
	pro=z(iz)*((xbeta*ri/(alpha*w))**xbeta)*(xkinit**(alpha+xbeta))
     &			-((one+(xbeta/alpha))*ri*xkinit)
	      
c	search for the constrained solution
		do 18 ia=1,ns2
		if ((theta*pro/(one+ri))+a(ia) .lt. xkinit) then
			diff=1.0D+00
			it=0
			bound=xkinit
			procons=pro
		
		do while (diff .gt. 1.0D-8 .or. it .gt. 1000)
		it=it+1
		if (it .eq. 1000) then
		write(6,543) it
 543      format('too many iterations: ',i10,f18.10)
		endif

		bound=(theta*procons/(one+ri))+a(ia)
		xlprov = (w/(xbeta*z(iz)*bound**alpha))**(one/(xbeta-one))
c		prov=z(iz)*((beta*ri/(alpha*w))**beta)*(bound**(alpha+beta))
c     &		-((one+(beta/alpha))*ri*bound)
		prov  = z(iz)*bound**alpha*xlprov**xbeta-w*xlprov-ri*bound
		diff=dabs(prov-procons)
		procons=prov
		end do

		xk(iz,ia)=bound
		xl(iz,ia)=xlprov
		xp(iz,ia)=procons
	xy(iz,ia)=z(iz)*bound**alpha*xlprov**xbeta
		
		else

		xk(iz,ia)=xkinit
		xl(iz,ia)=ri*xkinit/w
		xp(iz,ia)=pro
		xy(iz,ia)=z(iz)*((xbeta*ri/(alpha*w))**xbeta)*
     &(xkinit**(alpha+xbeta))
		endif

 18		continue

 114	continue				
 
 
c	open(22,file='profit1.out',form='formatted')	
c		do 81 i = 1,nkpts
c          write(22,499) a(i),xk(1,i),xk(2,i),xl(1,i),xl(2,i),
c     &	xy(1,i),xy(2,i),xp(1,i),xp(2,i)     
c 499      format(9f18.10)
c 81       continue
	close(22)
		
	return
	end

	subroutine prof2(r,w,xk,xl,xy,xp)
	implicit real*8 (a-h,o-z)
c	common/cparam/beta,z1
c	common/cres/xkgdp,aratio
      parameter (nkpts=150,nk2pts=5000,
     &tol=1.0D-08,nc=3,ns1=3,ns2=nk2pts,depreciation=0.06D+00,
     &alpha=0.366D+00,xbeta=0.88D+00-alpha,theta=0.6D+00)
c	,z1=1.0D+00)
      dimension ge(ns1,ns2,nc),gl(ns1,ns2,nc),
     &	ve(ns1,ns2),vl(ns1,ns2),vep(ns1,ns2),vlp(ns1,ns2),
     &	xmu(ns1,ns1,ns2),xmuprov(ns1,ns1,ns2),
     &	xmu1(ns2),xmu2(ns2),xmu3(ns2),xmu4(ns2),
     &	xk(ns1,ns2),xl(ns1,ns2),xy(ns1,ns2),xp(ns1,ns2),
     &a(ns2),z(ns1),delta(ns1),xkpts(nkpts)
      one=1.0D+00
	zero=0.0D+00
	ri=r+depreciation

c	z(1)=2.45D+00
	z(1)=1.0D+00
	z(2)=10.0D+00**(one-alpha-xbeta)*z(1)
	z(3)=100.0D+00**(one-alpha-xbeta)*z(1)

c	grid
	open(2,file='xkpts150.in',form='formatted')
      do 999 i = 1,nkpts
         read(2,570) xkpts(i)
 570     format(f18.10)
 999  continue
      close(2)

	xkbor=0.0D+00
      xklow = xkbor
      xkhigh = 1250.0D+00

      a(1) = xkpts(1)
      xkinc = (xkpts(nkpts)-xkpts(1))/dble(real(nk2pts-1))
      do 393 i = 2,nk2pts
         a(i) = a(i-1)+xkinc
 393  continue
	
c	profit function
	do 114 iz=1,ns1
c		  xkinit=(z(iz)*((alpha/ri)**(one-beta))*((beta/w)**beta))**
c     &		(one/(one-alpha-beta))
c
c		  pro=z(iz)*((beta*ri/(alpha*w))**beta)*(xkinit**
c     &		(alpha+beta))-((one+(beta/alpha))*ri*xkinit)
		xkinit=(z(iz)*((alpha/ri)**(one-xbeta))*((xbeta/w)**xbeta))
     &				**(one/(one-alpha-xbeta))  
	pro=z(iz)*((xbeta*ri/(alpha*w))**xbeta)*(xkinit**(alpha+xbeta))
     &			-((one+(xbeta/alpha))*ri*xkinit)
	      
c	search for the constrained solution
		do 18 ia=1,ns2
		if ((theta*pro/(one+ri))+a(ia) .lt. xkinit) then
			diff=1.0D+00
			it=0
			bound=xkinit
			procons=pro
		
		do while (diff .gt. 1.0D-8 .or. it .gt. 1000)
		it=it+1
		if (it .eq. 1000) then
		write(6,543) it
 543      format('too many iterations: ',i10,f18.10)
		endif

		bound=(theta*procons/(one+ri))+a(ia)
		xlprov = (w/(xbeta*z(iz)*bound**alpha))**(one/(xbeta-one))
c		prov=z(iz)*((beta*ri/(alpha*w))**beta)*(bound**(alpha+beta))
c     &		-((one+(beta/alpha))*ri*bound)
		prov  = z(iz)*bound**alpha*xlprov**xbeta-w*xlprov-ri*bound
		diff=dabs(prov-procons)
		procons=prov
		end do

		xk(iz,ia)=bound
		xl(iz,ia)=xlprov
		xp(iz,ia)=procons
	xy(iz,ia)=z(iz)*bound**alpha*xlprov**xbeta
		
		else

		xk(iz,ia)=xkinit
		xl(iz,ia)=ri*xkinit/w
		xp(iz,ia)=pro
	xy(iz,ia)=z(iz)*((xbeta*ri/(alpha*w))**xbeta)
     &*(xkinit**(alpha+xbeta))

c		xp(iz,ia)=pro
c	xy(iz,ia)=z(iz)*((beta*ri/(alpha*w))**beta)*(xkinit**(alpha+beta))
c		xl(iz,ia)=(beta*ri/(alpha*w))*xk(iz,ia)
		endif

 18		continue

 114	continue				
 			
 
c	open(21,file='profit2.out',form='formatted')	
c		do 83 i = 1,nk2pts
c          write(21,499) a(i),xk(1,i),xk(2,i),xl(1,i),xl(2,i),
c     &	xy(1,i),xy(2,i),xp(1,i),xp(2,i)     
c 499      format(9f18.10)
c 83       continue
c	close(21)

	return
	end





