subroutine qr2der(c0,param,xa,nxp,ya,nyp,f) parameter (na=189,nx=26,ny=6,nxm=5,nym=5,nep=10,nt=1000000) implicit double precision (a-g,o-z) double precision f(na), df(na,na), ft(na), dft(na,na), c work(nt), rcond, xa(nx+1), ya(ny+1), xx(nxm,nx), c wx(nxm,nx), yy(nym,ny), wy(nym,ny), x(nxm), u(nxm), c y(nym), v(nym), yyt(nep), ep(nep), c we(nep), wgt(nep), c0(na), c1(na), wrk(nep,4), c param(7) integer iwork(nt), nxp(nx), nyp(ny), nb(na), c jj(nep), kk(nep), ibc(na) c pi = 3.141592654d0 nx1 = nx+1 ny1 = ny+1 alpha = param(1) beta = param(2) delta = param(3) rho = param(4) sige = param(5) tau = param(6) gam = param(7) alph1 = alpha-1.d0 delt1 = 1.d0-delta call qgausl (nep, -2.88d0*sige, 2.88d0*sige, ep, we) tem1 = 1.d0/(2.d0*sige*sige) tem2 = 1.d0/dsqrt(pi/tem1) sum1 = 0.d0 do 50 i=1,nep wgt(i) = dexp(-ep(i)*ep(i)*tem1)*tem2*we(i) sum1 = sum1 + wgt(i) 50 continue do 60 i=1,nep wgt(i) = wgt(i)/sum1 60 continue do 80 i=1,nx call qgausl (nxp(i), xa(i), xa(i+1), x, u) do 80 j=1,nxp(i) xx(j,i) = x(j) wx(j,i) = u(j) 80 continue do 90 i=1,ny call qgausl (nyp(i), ya(i), ya(i+1), y, v) do 90 j=1,nyp(i) yy(j,i) = y(j) wy(j,i) = v(j) 90 continue ne=nx*ny crit = 0.0000001d0 maxit = 100 do 110 i=1,na f(i) = 0.d0 do 110 j=1,na df(i,j) = 0.d0 110 continue do 150 n=1,ne ix = mod(n-1,nx)+1 iy = int((n-1)/nx)+1 n1 = nxp(ix) n2 = nyp(iy) x1 = xa(ix) x2 = xa(ix+1) dx = x2-x1 y1 = ya(iy) y2 = ya(iy+1) dy = y2-y1 nd1= ix+nx1*(iy-1) nd2= nd1+1 nd4= nd1+nx1 nd3= nd4+1 do 160 i=1,n1 x(i) = xx(i,ix) u(i) = wx(i,ix) 160 continue do 170 j=1,n2 y(j) = yy(j,iy) v(j) = wy(j,iy) yp = (y2 - y(j))/dy ym = (y(j) - y1)/dy do 170 i=1,n1 xp = (x2 - x(i))/dx xm = (x(i) - x1)/dx c = c0(nd1)*xp*yp+c0(nd2)*xm*yp+ c c0(nd3)*xm*ym+c0(nd4)*xp*ym d = c-y(j)*x(i)**alpha c term in objective: -gam*(d*d*d+abs(d*d*d)) s = gam*3.d0*d*(d+dsign(d,0.d0)) dsdc = gam*6.d0*(d+dsign(d,0.d0)) xt = x(i)**alpha*y(j)-c+delt1*x(i) dftxt = alpha*xt**(alpha-1.d0) ii = 1 do 180 m=1,nx j1 = int( dmax1(xt-xa(m),0.d0)/(xt-xa(m))+.2 ) ii = ii + j1 180 continue ii = max(1,ii-1) do 190 l=1,nep yyt(l) = y(j)**rho*dexp(ep(l)) jj(l) = 1 do 200 m=1,ny j1 = int( dmax1(yyt(l)-ya(m),0.d0)/ c (yyt(l)-ya(m))+.2) jj(l) = jj(l) + j1 200 continue jj(l) = max(1,jj(l)-1) 190 continue sum1 = 0.d0 sum2 = 0.d0 sum3 = 0.d0 sum4 = 0.d0 sum5 = 0.d0 do 220 l=1,na work(l) = 0.d0 220 continue do 230 l=1,nep yt = yyt(l) x1t = xa(ii) x2t = xa(ii+1) dxt = x2t-x1t y1t = ya(jj(l)) y2t = ya(jj(l)+1) dyt = y2t-y1t xpt = (x2t-xt)/dxt xmt = (xt-x1t)/dxt ypt = (y2t-yt)/dyt ymt = (yt-y1t)/dyt nd1t = (jj(l)-1)*nx1+ii nd2t = nd1t+1 nd4t = nd1t+nx1 nd3t = nd4t+1 kk(l) = nd1t c c compute c * next period consumption: ct c * beta*du(c[t+1])/dc[t+1]: dutct c * dc[t+1]/dk[t+1]: dctxt c * dF[t+1]/dk[t+1]: dftkt c * int u'(c~)*F'(k~)f(e)de: sum1 c * d(sum1)/d(c_a): dsum c * d(sum1)/d(c_b): dsum (2nd instance) c where c_a,c_b are unknowns in c,ct respectively c ct = c0(nd1t)*xpt*ypt+c0(nd2t)*xmt*ypt+ c c0(nd3t)*xmt*ymt+c0(nd4t)*xpt*ymt dt = ct-yt*xt**alpha st = gam*3.d0*dt*(dt+dsign(dt,0.d0)) dstdc = gam*6.d0*(dt+dsign(dt,0.d0)) dutct = beta*ct**(-tau) dctxt = (-(c0(nd1t)-c0(nd2t))*ypt c +(c0(nd3t)-c0(nd4t))*ymt)/dxt dftkt = yt*dftxt+delt1 sum1 = sum1 + wgt(l)*(dutct*dftkt-beta*delt1*st) dsum = wgt(l)*dutct*(tau*dctxt*dftkt/ct c -alph1*yt*dftxt/xt)+ c wgt(l)*beta*delt1*dstdc*(dctxt-yt*dftxt) sum2 = sum2 + dsum*xp*yp sum3 = sum3 + dsum*xm*yp sum4 = sum4 + dsum*xm*ym sum5 = sum5 + dsum*xp*ym dsum = -wgt(l)*(tau*dutct*dftkt/ct+ c beta*delt1*dstdc) wrk(l,1) = dsum*xpt*ypt wrk(l,2) = dsum*xmt*ypt wrk(l,3) = dsum*xmt*ymt wrk(l,4) = dsum*xpt*ymt 230 continue tem1 = tau*c**(-tau-1.d0)+dsdc work(nd1) = sum2 + tem1*xp*yp work(nd2) = sum3 + tem1*xm*yp work(nd3) = sum4 + tem1*xm*ym work(nd4) = sum5 + tem1*xp*ym do 232 l = 1,nep nd1t = kk(l) nd2t = nd1t+1 nd4t = nd1t+nx1 nd3t = nd4t+1 work(nd1t) = work(nd1t) + wrk(l,1) work(nd2t) = work(nd2t) + wrk(l,2) work(nd3t) = work(nd3t) + wrk(l,3) work(nd4t) = work(nd4t) + wrk(l,4) 232 continue res = sum1-c**(-tau)+s tem1 = u(i)*v(j)*res f(nd1) = f(nd1) + tem1*xp*yp f(nd2) = f(nd2) + tem1*xm*yp f(nd3) = f(nd3) + tem1*xm*ym f(nd4) = f(nd4) + tem1*xp*ym 170 continue 150 continue return end