c     Parameters can be varied by directly altering them in the
c     code. Sometimes, depending on what the solution looks like,
c     the starting value for the number of firms also needs to be
c     modified in order for a solution to be reached.

      parameter (maxn = 2, maxwa = 50)
      integer i,j,k,iflag,info,lwa
      double precision g,w0,lb,p,a,b,c
      double precision l,e,wr,ej,en,de,de1,ew,ew1,wb,w0b,dudw,dudn,dndw
      double precision u,du,dpidw,dpidn
      double precision val(maxn),fval(maxn),saveval(maxn),wa(maxwa)
      double precision tol,dpmpar
      integer newloop1,newloop2

      common /params/ g,w0,lb,p,a,b,c

      external eqn

      lwa=maxwa

c      open (unit=1,file='output.txt',status='new')
c
c     initialize parameter values
c

c     set gamma
      g = .5
c     return from home production
      w0 = .25
c     endowment of labor
      lb = 1
c     net revenue product of labor
      p = 1
c     set b parameter of the Pareto distribution
      b = 1
c     set error tolerance
      tol = dsqrt(dpmpar(1))
c
c     set starting value for number of employers
c
      saveval(1)=10
c
c     set starting value for wage
c
      saveval(2)=.8*p
c
c     print header line
c
      print *, "a,b,c,p,w0,lb,n,w,l,e,ej,en,dldw,dldn,de,de1,ew,ew1
     *,dudw,dudn,dndw,du,u,euw"

c     (w0/b)^(1-g)
      w0b = (w0/b)**(1-g)

c   loop for setting Pareto distribution parameter a
      a=0.25
      newloop1 = 1
c   use old starting values
      val(1) = saveval(1)
      val(2) = saveval(2)
c     loop for setting fixed production cost
        do k = 1,20
         c = .005*k

c     call subroutine to numerically compute n and w
         call hybrd1(eqn,maxn,val,fval,tol,info,wa,lwa)

c     (w/b)^(1-g)
         wb = (val(2)/b)**(1-g)
c     cross establishment elasticity of labor supply
         ej = a*(val(1)-1)
       if (val(2).ge.w0) then
c     (w0/w)^(an)
          wr = (w0/val(2))**(a*val(1))
c     labor supply
          l = (1-g)*lb*((1-wr)/val(1)+wr/(val(1)+m))

c     establishment elasticity of labor supply
          e = a*((val(1)+m)*(val(1)-1)+m*wr)/(val(1)+m-m*wr)
c     elasticity of labor supply with respect to n
          en = (val(1)**2+2*val(1)*m*(1-wr)+m**2*(1-wr)+a*val(1)*m*
     *         (val(1)+m)*wr*log(w0/val(2)))/
     *         ((val(1)+m)*(val(1)+m*(1-wr)))
c     partial derivative of welfare wrt the wage and n
          dudw = (a*val(1)/val(2))*(((1-g)/(a*val(1)+1-g))
     *           *wb+(a*val(1)/(a*val(1)+1-g)-
     *           a*(val(1)+m)/(a*(val(1)+m)+1-g))*w0b*wr)
          dudn = (a*(1-g)/(a*(val(1)+m)+1-g)**2)*w0b*wr
     *           +(a*(1-g)/(a*val(1)+1-g)**2)*(wb-w0b*wr)
     *           +(a*a*(val(1)+m)/(a*(val(1)+m)+1-g)
     *           -a*a*val(1)/(a*val(1)+1-g))*w0b*wr*log(w0/val(2))

c     partial derivative of profits wrt n and w
          dpidn = (p-val(2))*val(2)*(val(1)**2+m*m*(1-wr)
     *            +2*val(1)*m*(1-wr)
     *            +a*val(1)*m*(val(1)+m)*wr*log(w0/val(2)))
          dpidw = val(1)*(val(1)+m)*(val(1)*val(2)+m*(val(2)*(1-wr)
     *            -a*val(1)*(p-val(2))*wr))

c     partial derivatives of labor supply with respect to the wage
c     and n
          dldw = (1-g)*lb*wr*(a*m/(val(2)*(val(1)+m)))
          dldn = (1-g)*lb*(-wr/(val(1)+m)**2-(1-wr)/val(1)**2
     *           -(a/val(1)-a/(val(1)+m))*wr*log(w0/val(2)))

c     equilibrium welfare
          u = wr*w0b*a*(val(1)+m)/(a*(val(1)+m)+1-g)
     *        +(wb-wr*w0b)*a*val(1)/(a*val(1)+1-g)
       else
c   (w/w0)^(am)
        wr=(val(2)/w0)**(a*m)

c   labor supply
        l = (1-g)*lb*wr/(val(1)+m)

c   establishment elasticity of labor supply
        e = a*(val(1)+m-1)
c   elasticity of labor supply wrt n
        en = val(1)/(val(1)+m)

c   partial derivative of welfare wrt w and n
        dudw = a*val(1)/val(2)*(1-g)*wb*wr/(a*(val(1)+m)+1-g)
          dudn = a*(1-g)*wb*wr/(a*(val(1)+m)+1-g)**2

c     partial derivative of profit wrt n and w
        dpidn = -(p-val(2))*wr/(val(1)+m)**2
        dpidw = (a*m*(p-val(2))-val(2))*wr/(val(2)*(val(1)+m))

c     partial derivative of labor supply wrt w and n
        dldw = (1-g)*lb*a*m*wr/(val(2)*(val(1)+m))
        dldn = -(1-g)*lb*wr/(val(1)+m)**2
c     equilibrium welfare
        u = a*(val(1)+m)/(a*(val(1)+m)+1-g)*wr*wb+a*m/(a*m+1-g)*
     *        (w0b-wr*wb)
       end if
c     change in total employment given a change in the minimum
c     wage
       de1 = l*val(1)/val(2)*(e-ej/en)

c     total derivative
         dndw = -dpidw/dpidn
c     general equilibrium change in welfare from a change in min
c     wage
         du = dudw+dudn*dndw
         de = (l+val(1)*dldn)*dndw+val(1)*dldw
c     employment elasticity of the minimum wage
         ew = de*val(2)/(l*val(1))
c     employment elasticity using Walsh decomposition to check
c     for consistency
         ew1 = e-ej/en
c     elasticity of welfare wrt the minimum wage
         euw = du*val(2)/u

c     print results
        if (val(1).ge.2.0) then
          print *, a,",",b,",",c,",",p,",",w0,",",lb,",",val(1)
     *    ,",",val(2),",",l,",",e,",",ej,",",en,",",dldw,",",dldn
     *    ,",",de,",",de1,",",ew,",",ew1,",",dudw,",",dudn,",",dndw
     *    ,",",du,",",u,",",euw,",",fval(1),",",fval(2),",",info
        end if
c     when at the beginning of a "c" loop, store solution for
c     use as starting values
        if (newloop1.eq.1) then
         newloop1 = 0
         saveval(1) = val(1)
         saveval(2) = val(2)
        end if
       end do
      end

c
c     subroutine for computing the value of the zero profit
c     condition and the deviation of the computed wage given
c     values for n and w
c
      subroutine eqn(nvar,x,fvec,iflag)
      integer nvar,iflag
      double precision x(nvar),fvec(nvar),l,e,wr
      double precision g,w0,lb,p,a,b,c
      common /params/ g,w0,lb,p,a,b,c

c   compute intermediate values
      if (x(2).ge.w0) then
c     (w0/w)^(an)
       wr = (w0/x(2))**(a*x(1))
c     labor supply
       l = (1-g)*lb*((1-wr)/x(1)+wr/(x(1)+m))
c     establishment level elasticity of labor supply
       e = a*((x(1)+m)*(x(1)-1)+m*wr)/(x(1)+m-m*wr)
      else
c   (w/w0)^(am)
       wr=(x(2)/w0)**(a*m)
c   labor supply
       l = (1-g)*lb*wr/(x(1)+m)
c   establishment level elasticity of labor supply
       e = a*(x(1)+m-1)
      end if
c
c     compute net profits
c
      fvec(1)=(p-x(2))*l-c
c
c     computed wage less initial value of the wage
c
      fvec(2)=p*e/(1+e)-x(2)
      end