c**********************************************************************
c
c   subrotine oglqa  given zeta and agent, computes a second order
c     approximation with increments h(j), around the steady states
c     value of the utility
c
c**********************************************************************
c
      subroutine oglqa(   lgthlf,  nusho,  param,    zz,   eff,     h,
     1                    pai, r,      m,      i,     x,    xs,     s)
c

      implicit real*8(a-h,o-z)
      integer             lgthlf,  nusho,  m,i,j,jj
c
      real *8         pai( nusho,  nusho)
      real *8           r( nusho, lgthlf, 2*nusho+5, 2*nusho+5)
      real *8           h(  2*nusho+4)
      real *8       param(  8)
      real *8           x(  2*nusho+4)
      real *8          xs(  2*nusho+4)
      real *8         eff
      real *8               ubar,      u1,     u2,     z,s,zz
c
      external ret
c*************************begin****************************************
c
      ubar=ret(lgthlf,pai,xs,param,m,nusho,zz,eff,s)
      r( m, i, 1, 1) = ubar
      do 10 j= 1, 2*nusho+4
          do 8 jj= 1, 2*nusho+4
8            x(jj)=xs(jj)
          x( j) = xs( j) + h( j)
          u1 =ret( lgthlf,pai,x,param,m,nusho,zz,eff,s)
          x( j) = xs( j) - h( j)
          u2 =ret( lgthlf,pai,x,param,m,nusho,zz,eff,s)
          r( m, i, 1, j+1) = ( (u1 - u2) / ( 2. * h( j) ) ) / 2.
          r( m, i, j+1, 1) = r( m, i, 1, j+1)
          r( m, i, j+1,j+1)= ( ((u1 + u2) / 2.) - ubar ) / (h(j) * h(j))
          do 9 jj= 1, j-1
              x( j) = xs( j) + h( j)
              x(jj) = xs(jj) + h(jj)
              z=ret(lgthlf,pai,x,param,m,nusho,zz,eff,s)
              x(jj) = xs(jj) - h(jj)
              z=z-ret(lgthlf,pai,x,param,m,nusho,zz,eff,s)
              x( j) = xs( j) - h( j)
              z=z+ret(lgthlf,pai,x,param,m,nusho,zz,eff,s)
              x(jj) = xs(jj) + h(jj)
              z=(z-ret(lgthlf,pai,x,param,m,nusho,
     1                           zz,eff,s))/(8.*h(jj)*h(j))
              x(jj) = xs(jj)
              r( m, i, jj+1, j+1) = z
              r( m, i, j+1, jj+1) = z
9         continue
          x( j) = xs( j)
10    continue
      z=0.
      do 13 j= 1, 2*nusho+4
          x(j)=0.
          do 14 jj= 1, 2*nusho+4
14            x( j)= x( j) + r( m, i, jj+1, j+1) * xs( jj)
          z= z + x( j) * xs( j)
13    continue
      r( m, i, 1, 1)=r( m, i, 1, 1) + z
      z=0.
      do 15 j= 1, 2*nusho+4
15        z=z+ r( m, i, 1, j+1) * xs( j) * 2
      r( m, i, 1, 1)= r( m, i, 1, 1) - z
      do 11 j= 1, 2*nusho+4
          z=0.
          do 12 jj= 1, 2*nusho+4
12            z=z+ r( m, i, j+1, jj+1) * xs( jj)
          r( m, i, 1, j+1)= r( m, i, 1, j+1) - z
          r( m, i, j+1, 1)= r( m, i, 1, j+1)
11    continue
      return
      end
c
c*********************************************************************
c
c      function ret computes utility for variables x
c
c*********************************************************************
c
      real *8 function ret( lgthlf,pai, x, param, m, nusho,zz, eff,s)
c
      implicit real*8(a-h,o-z)
      integer       lgthlf,m,nusho,mm
c
      real *8         pai( nusho,  nusho)
      real *8   param(      8 )
      real *8       x(   2*nusho+4 )
      real *8       eff
      real *8       delta, zeta, alfa, gamma, theta, eta, one, kk
      real *8       nn, aa, ee, ll,s,zz,ss,cont
c
c*********************************************************************
c
      one   = 1.
      alfa  = param( 1)
      eta   = param( 2)
      gamma = param( 4)
      theta = param( 5)
      zeta  = zz
      grth  = param( 7)+one
      delta = param( 8)
      ss    = s
c
c***********************begin*****************************************
c
      cont = 0.
      ee   = eff
      aa   = x( 1)
      kk   = x( 2)
      nn   = x( 3)
      do 10 mm=1,nusho
10        cont=cont+x(3+mm)*x(3+nusho+mm)
      ll   = x( 2*nusho+4)
c
      ret=(one/(one-gamma)) / ( ( aa * (one-delta+(one-alfa)*zeta* nn**
     1    (alfa)/kk**alfa) +ee*(one-ll)*alfa*zeta*kk**(one-alfa)/nn**
     2    (one-alfa)-cont*grth)
     3              **theta*ll**(one-theta)-eta)**(gamma-one)
c
      return
      end
c
c*********************************************************************