c PROGRAM: ONESECT.F c c DESCRIPTION: Computation of equilibrium in one-sector endogenous c growth model. (See EXPLAINING CROSS-COUNTRY INCOME c DIFFERENCES by McGrattan and Schmitz for details.) c c MODEL: Consumers: c c E sum_t beta^t U(C[t]/N[t],L[t]/N[t],G1[t]/N[t])*N[t] c c subject to c c (1+tau_c[t])*C[t] + (1+tau_xk[t])*Xk[t] c + (1+tau_xh[t])*Xh[t]} c <= (1-tau_k[t])*r[t]*K[t] c +(1-tau_h[t])*w[t]*l[t]*H[t] + T[t] c c K[t+1] = (1-del_k)*K[t] + Xk[t] c H[t+1] = (1-del_h)*H[t] + Xh[t] c C[t],Xk[t],Xh[t] >= 0, 0<= L[t]/N[t] <= 1 c N[t+1] = (1+eta)* N[t] (=population) c l[t] = L[t]/N[t] (=labor input) c c Factor prices: c r[t] = F_1(K[t],l[t]*H[t],G2[t]) c w[t] = F_2(K[t],l[t]*H[t],G2[t]) c c Gov't spending: c Gj[t] = zj[t] * F(K[t],l[t]*H[t],G2[t]) c c Resource constraint: c C[t] + Xk[t] + Xh[t] + G1[t] + G2[t] c = F(K[t],l[t]*H[t],G2[t]) c c Stochastic processes: c s[t]: Markov chain with transition probabilities pi c Gj[t], tau_j[t] take on different values for each s[t]. c c c PROBLEM: find c=c(x,i), xh=xh(x,i) that solve the intertemporal c first order conditions -- once the equations have been c transformed by dividing all variables on the growth path c by H[t] and replacing H[t+1]/H[t] by gam[t]. Thus the c state variable is x[t]=K[t]/H[t] and the decision variables c are c[t]=C[t]/H[t] and xh[t]=Xht[t]/H[t]. The remaining c decision variables are obtained by way of the static first c order conditions. (See Handbook notes for details.) c c ALGORITHM: the finite element method with piecewise linear c shape functions N_a(x), ie., f(x,i) = sum N_a(x)*f^i_a, c where f^i_a's are the coefficients computed here. c A sparse-matrix solver is used to invert the Jacobian. c c FUNCTIONS: U(c,l,g) = [c*l^psi*g^zeta]^(1-sigma)/(1-sigma) c F(k,l,g) = a*k^alp_k*l^alp_h*g^(1-alp_k-alp_h) c c TEST CASE: U(c,l,g) = log(c) c F(k,l,g) = a*k^alpha*l^(1-alpha) c and del_k = del_h = 1 c ==> normalized variables all constant c c INPUT FILE: (a) onesect.inp: the following need to be specified c c Input Description Format c ------------------------------------------------------ c a parameter in F(.) real > 0 c alpk parameter in F(.) real in [0,1] c alph parameter in F(.) real in [0,1] c beta discount factor real in (0,1) c chi penalty parameter real >= 0 c delk depreciation rate real in [0,1] c delh depreciation rate real in [0,1] c eta population growth (%) real >= 0 c psi parameter in U(.) real >= 0 c sigma parameter in U(.) real > 0 c zeta parameter in U(.) real >= 0 c ievenx option for evenly 1 (=yes) or 0 c spaced x-partition c initialize option to provide 1 (=yes) or 0 c initial guesses c nbc1 number of boundary integer c conditions for c c nbc2 number of boundary integer c conditions for xh c x capital grid if ievenx=1, c 2 x 1 real vector c with x(1),x(nx) c otherwise, c nx x 1 real vector c with entire grid c ns number of policy (should agree with c states dimension statement) c policies =[tauxk,tauxh, ... ns x 7 real matrix c tauk,tauh, ... c tauc,z1,z2] c nnzp number of nonzero integer c elements in pi c i,j,pi(i,j) elements of pi nnzp triplets of c (integer,integer,real) c c initial guess for nx*ns x 1 real c consumption function vector; read c(:,1) c (if initialize=1) first, c(:,2) c second, etc. c xh initial guess for nx*ns x 1 real c investment function vector; read xh(:,1) c (if initialize=1) first, xh(:,2) c second, etc. c i,j,c points c(x(i),j) triplets of (integer, c that are to be fixed integer,real) c during computation c i,j,xh points xh(x(i),j) triplets of (integer, c that are to be fixed integer,real) c during computation c c Notes: c i. see onesect.inp for an example input file; c ii. if grid sizes change, edit this file and c change dimensions; c iii. all other inputs are set in this file; c see comment lines below. c c OUTPUT: (b) onesect.nxt: new input file with current solution c (c) onesect.dat: input file for Matlab post-processing c routine -- type onesect in Matlab. c c REQUIRED (SPARSKIT2) ilut, matvec c PROGRAMS: (PRESS) qgausl, qcksrt3 c c REFERENCES: McGrattan (1996), JEDC, 20: 19-42 c Press, et. al. (1986), NUMERICAL RECIPES c Saad (1996), ITERATIVE METHODS FOR SPARSE LINEAR SYSTEMS c c SEE ALSO: files in ~erm/work/ckm1/ffiles, ~erm/work/jones/ffiles c c Ellen McGrattan, 6-19-97 c Revised, ERM, 11-19-98 c * c c Set dimensions: c nx = length of x partition c ns = number of discrete states c na = 2*nx*ns c nq = maximum number of quadrature points per element (in x c direction); nq>=nxp(i), i=1:ne. c nm = number of columns stored for sparse matrix solver (<=50) c ldf = dimension of Jacobian vector (df) c lpi = dimension of vector for transition probabilities (pi) c parameter (nx=11,ns=800,na=17600,nq=5,nm=40,ldf=1000000,lpi=10000) implicit double precision (a-h,o-z) double precision f(na), df(ldf), wrk1(na), wrk2(na), & xa(nx), ax(nq,nx-1), wx(nq,nx-1), c0(nx,ns), c1(nx,ns), & xh0(nx,ns), xh1(nx,ns), tauxh(ns), tauxk(ns), tauk(ns), & tauh(ns), tauc(ns), z1(ns), z2(ns), pi(lpi), lab, lab0, & lei, labt, leit, lhs1, lhs2, alu(ldf), vv(na,nm+1), & sol(na) integer iwork(na), nxp(nx-1), ibc(na), ipi(ns+1), jpi(lpi), & ind(ldf), ia(ldf), ja(ldf), ip(na+1), jlu(ldf), ju(na) c nxm1 = nx-1 ne = nxm1 c c Open input and output files. c open(unit=5, file='onesect.inp') open(unit=7, file='onesect.nxt') open(unit=8, file='onesect.dat') c c Read in the parameters. c read(5,*) a read(5,*) alpk read(5,*) alph read(5,*) beta read(5,*) chi read(5,*) delk read(5,*) delh read(5,*) eta read(5,*) lab0 read(5,*) psi read(5,*) sigma read(5,*) zeta read(5,*) ievenx read(5,*) initialize read(5,*) nbc1 read(5,*) nbc2 c nbc = nbc1+nbc2 alp = alpk+alph alp1 = 1.d0-alp alpi = 1.d0/alp aa = a**alpi alphak = alpk*alpi alphah = alph*alpi alphai = alp1*alpi delk1 = 1.d0-delk delh1 = 1.d0-delh sigm1 = 1.d0-sigma zeta1 = 1.d0+zeta sigz = sigma-zeta*sigm1 betn = beta*(1.d0+eta/100.d0)**sigz c c Read in grid for capital stock. c if (ievenx.eq.1) then read(5,*) xa(1) read(5,*) xa(nx) tem = (xa(nx)-xa(1))/float(nxm1) do 10 i=2,nxm1 xa(i) = xa(i-1) + tem 10 continue else do 20 i=1,nx read(5,*) xa(i) 20 continue endif c c Read in policy variables. c read(5,*) i if (i.ne.ns) then write(*,*) 'ERROR: wrong number of policy states' stop endif do 30 i=1,ns read(5,*) tauxk(i),tauxh(i),tauk(i),tauh(i),tauc(i),z1(i),z2(i) 30 continue c c Read in transition matrix which is stored in coordinate format. c Convert to compressed sparse row format. c read(5,*) nnzp do 40 i=1,nnzp read(5,*) j,k,pi(i) ind(i) = (j-1)*ns+k 40 continue call qcksrt3(nnzp,ind,pi) ipi(1) = 1 j = 1 jpi(1) = mod(ind(1)-1,ns)+1 do 50 i=2,nnzp if ((ind(i)-1)/ns+1 .ne. j) then j = j+1 ipi(j) = i endif jpi(i) = mod(ind(i)-1,ns)+1 50 continue ipi(ns+1) = nnzp+1 c c Initialize consumption and investment functions. c if (initialize.eq.1) then do 60 j=1,ns do 60 i=1,nx read(5,*) c0(i,j) 60 continue do 65 j=1,ns do 65 i=1,nx read(5,*) xh0(i,j) 65 continue else c c If initial guess is not given, set c, xh as in test case c do 70 j=1,ns do 70 i=1,nx xh0(i,j) = beta*alpk**alpk*(1.d0-alpk)**(1.d0-alpk) c0(i,j) = (1.d0-beta)*(alpk/(1.d0-alpk))**alpk 70 continue endif c c Impose boundary conditions. c do 75 i=1,na ibc(i) = 0 75 continue do 80 k=1,nbc1 read(5,*) i,j,tem ibc(2*(i+(j-1)*nx)-1) = 1 c0(i,j) = tem 80 continue do 90 k=1,nbc2 read(5,*) i,j,tem ibc(2*(i+(j-1)*nx)) = 1 xh0(i,j) = tem 90 continue c c Set the following parameters for the algorithm: c * nxp(i), i=1,nx-1, the # of quadrature points per x interval c * crit, the stopping criterion c * maxit, the maximum number of interations. c do 100 i=1,nxm1 nxp(i) = 3 100 continue crit = 0.0000001d0 maxit = 200 c c Compute quadrature abscissas and weights. c do 110 i=1,nxm1 call qgausl (nxp(i), xa(i), xa(i+1), ax(1,i), wx(1,i)) 110 continue c c The Newton-Raphson iteration, i.e., find x such that f(x)=0 where c x are coefficients of the finite element approximation. c do 120 it=1,maxit c c Initialize f and its derivative df. c do 130 i=1,na f(i) = 0.d0 130 continue do 140 i=1,ldf df(i) = 0.d0 140 continue nnz = 0 c c Construct f(c0;xh0) = int R(x;c0,xh0) N_n(x) dx and its derivative c on each element n. c do 150 j=1,ns do 150 n=1,ne c c Element n is the interval [xa(n),xa(n+1)] c x1 = xa(n) x2 = xa(n+1) c c Compute f at all quadrature points on the element. c do 160 i=1,nxp(n) x = ax(i,n) w = wx(i,n) xi = 2.d0*(x-x1)/(x2-x1)-1.d0 bs1 = 0.5d0*(1.d0-xi) bs2 = 0.5d0*(1.d0+xi) c c Compute consumption and investment decisions using a finite c element approximation. c c = c0(n,j)*bs1 + c0(n+1,j)*bs2 xh = xh0(n,j)*bs1 + xh0(n+1,j)*bs2 c c Compute intermediate parameters. c tauxk1= 1.d0+tauxk(j) tauxh1= 1.d0+tauxh(j) tauk1 = 1.d0-tauk(j) tauh1 = 1.d0-tauh(j) tauc1 = 1.d0+tauc(j) tem1 = psi*c tem2 = tauh1*alph*aa*x**alphak*z2(j)**alphai/tauc1 c c Back out other variables from static first order conditions c and definitions. c lab = lab0 lei = 1.d0-lab do 165 k=1,10 tem = tem1 - tem2*lab**(-alphak)*lei dtem = tem2*lab**(-alphak)*(alphak*lei/lab+1.d0) lab = lab - tem/dtem lei = 1.d0-lab 165 continue ys = a*x**alpk*lab**alph g2 = (z2(j)*ys)**alpi y = ys*g2**alp1 g1 = z1(j)*y gam = delh1+xh gams = gam**sigz xk = y-g1-g2-xh-c u = (c*lei**psi*g1**zeta)**sigm1/sigm1 duc = sigm1*u/c pxk2 = chi*dmin1(xk,0.d0)**2.d0 pxh2 = chi*dmin1(xh,0.d0)**2.d0 c c Compute left-hand sides of 2 dynamic first-order conditions. c lhs1 = gams*(duc*tauxk1/tauc1-pxk2) lhs2 = gams*(duc*tauxh1/tauc1-pxh2) c c Compute total derivatives of lhs1 and lhs2 in terms of unknowns c (i.e., tl1c, tl1xh, tl2c, tl2xh). c d2uc = -sigma*duc/c ducl = -sigm1*psi*duc/lei ducg = sigm1*zeta*duc/g1 dg1l = alphah*g1/lab dl1xk = -gams*2.d0*chi*dmin1(xk,0.d0) dl1l = gams*(ducl+ducg*dg1l)*tauxk1/tauc1 dl1c = gams*d2uc*tauxk1/tauc1 dl1xh = sigz*gams/gam*(duc*tauxk1/tauc1-pxk2) dl2xk = 0.d0 dl2l = gams*(ducl+ducg*dg1l)*tauxh1/tauc1 dl2c = gams*d2uc*tauxh1/tauc1 dl2xh = sigz*gams/gam*(duc*tauxh1/tauc1-pxh2)- & gams*2.d0*chi*dmin1(xh,0.d0) dxkc = -1.d0 dxkxh = -1.d0 dxkl = alphah*(y-g1-g2)/lab tlc = -psi/(aa*alph*z2(j)**alphai*tauh1/tauc1* & (x/lab)**alphak*(alphak*lei/lab+1.d0)) txkc = dxkc+dxkl*tlc txkxh = dxkxh tl1c = dl1c + dl1l*tlc + dl1xk*txkc tl1xh = dl1xh+ dl1xk*txkxh tl2c = dl2c + dl2l*tlc tl2xh = dl2xh c c Find next period's capital stock on grid and derivatives of c next period's capital with respect to unknowns. c xt = (delk1*x+xk)/gam ii = 1 do 170 k=1,nxm1 if (xt.ge.xa(k)) ii = k 170 continue x1t = xa(ii) x2t = xa(ii+1) dxt = x2t-x1t xi = 2.d0*(xt-x1t)/dxt-1.d0 bs1t = 0.5d0*(1.d0-xi) bs2t = 0.5d0*(1.d0+xi) txtc = txkc/gam txtxh = (txkxh-xt)/gam c c Initialize variables for summation. c sum1 = 0.d0 sum2 = 0.d0 ts1c = 0.d0 ts2c = 0.d0 ts1xh = 0.d0 ts2xh = 0.d0 do 180 k=1,na wrk1(k) = 0.d0 wrk2(k) = 0.d0 180 continue c c The following loop computes the sums in the Euler equations. c do 190 l=ipi(j),ipi(j+1)-1 k = jpi(l) c c Compute next-period consumption and investment decisions c using a finite element approximation. c ct = c0(ii,k)*bs1t + c0(ii+1,k)*bs2t xht = xh0(ii,k)*bs1t + xh0(ii+1,k)*bs2t c c Compute intermediate parameters. c tauxk2 = 1.d0+tauxk(k) tauxh2 = 1.d0+tauxh(k) tauk2 = 1.d0-tauk(k) tauh2 = 1.d0-tauh(k) tauc2 = 1.d0+tauc(k) tem1 = psi*ct tem2 = tauh2*alph*aa*xt**alphak*z2(k)** & alphai/tauc2 c c Back out other variables from static first order conditions c and definitions. c labt = lab leit = 1.d0-labt do 195 m=1,15 tem = tem1 - tem2*labt**(-alphak)*leit dtem = tem2*labt**(-alphak)*(alphak*leit/ & labt+1.d0) labt = labt - tem/dtem leit = 1.d0-labt 195 continue ys = a*xt**alpk*labt**alph g2t = (z2(k)*ys)**alpi yt = ys*g2t**alp1 g1t = z1(k)*yt xkt = yt-g1t-g2t-xht-ct ut = (ct*leit**psi*g1t**zeta)**sigm1/sigm1 dutc = sigm1*ut/ct pxkt2 = chi*dmin1(xkt,0.d0)**2.d0 pxht2 = chi*dmin1(xht,0.d0)**2.d0 pxkt3 = pxkt2*dmin1(xkt,0.d0) pxht3 = pxht2*dmin1(xht,0.d0) c c Compute right-hand sides sums of 2 dynamic first-order c conditions. c sum1 = sum1 + betn*pi(l)* (dutc*(tauxk2*delk1+ & tauk2*alpk*yt/xt)/tauc2 - pxkt2*delk1) sum2 = sum2 + betn*pi(l)* (dutc*(tauxh2*delh1+ & tauh2*alph*yt)/tauc2 - pxht2*delh1 + & (sigm1*zeta1/3.d0-1.d0)*(pxht3+pxkt3)) c c Compute total derivatives of sum1 and sum2 in terms c of unknowns (i.e., ts1c, ts1xh, ts1ct, ts1xht, ts2c, c ts2xh, ts2ct, ts2xht). c d2utc = -sigma*dutc/ct dutcl = -sigm1*psi*dutc/leit dutcg = sigm1*zeta*dutc/g1t dytx = alphak*yt/xt dytl = alphah*yt/labt dg1tx = z1(k)*dytx dg1tl = z1(k)*dytl tctxt = -(c0(ii,k) -c0(ii+1,k))/dxt txhtxt = -(xh0(ii,k)-xh0(ii+1,k))/dxt tem1 = aa*alph*z2(k)**alphai*tauh2/tauc2 tem2 = tem1*(xt/labt)**alphak tem3 = tem2*(alphak*leit/labt+1.d0) dltxt = tem2*alphak*leit/xt/tem3 tltct = -psi/tem3 tltxt = dltxt+tltct*tctxt dxktxt = alphak*(yt-g1t-g2t)/xt dxktlt = alphah*(yt-g1t-g2t)/labt dxktct = -1.d0 txktct = dxktct+txktlt*tltct txktxht = -1.d0 txktxt = dxktxt+dxktct*tctxt+dxktlt*tltxt+ & txktxht*txhtxt tem = (tauxk2*delk1+tauk2*alpk*yt/xt)/tauc2 dr1xkt = -2.d0*chi*dmin1(xkt,0.d0)*delk1 dr1lt = (dutcl+dutcg*dg1tl)*tem+ & dutc*tauk2*alpk*dytl/xt/tauc2 dr1ct = d2utc*tem dr1xt = dutcg*dg1tx*tem+dutc*tauk2*alpk* & (dytx-yt/xt)/xt/tauc2 tr1xt = dr1xt+dr1xkt*txktxt+dr1lt*tltxt+ & dr1ct*tctxt tem = (tauxh2*delh1+tauh2*alph*yt)/tauc2 dr2xkt = pxkt2*(sigm1*zeta1-3.d0) dr2lt = ((dutcl+dutcg*dg1tl)*tem+ & dutc*tauh2*alph*dytl)/tauc2 dr2ct = d2utc*tem dr2xht = -2.d0*chi*dmin1(xht,0.d0)*delh1+pxht2* & (sigm1*zeta1-3.d0) dr2xt = dutcg*dg1tx*tem+dutc*tauh2*alph*dytx/tauc2 tr2xt = dr2xt+dr2xkt*txktxt+dr2lt*tltxt+ & dr2ct*tctxt+dr2xht*txhtxt ts1c = ts1c +betn*pi(l)*tr1xt*txtc ts1xh = ts1xh+betn*pi(l)*tr1xt*txtxh ts1ct = betn*pi(l)*(dr1ct +dr1xkt*txktct + & dr1lt*tltct) ts1xht = betn*pi(l)*dr1xkt*txktxht ts2c = ts2c +betn*pi(l)*tr2xt*txtc ts2xh = ts2xh+betn*pi(l)*tr2xt*txtxh ts2ct = betn*pi(l)*(dr2ct +dr2xkt*txktct + & dr2lt*tltct) ts2xht = betn*pi(l)*(dr2xht+dr2xkt*txktxht) c c Store derivatives in work vectors `wrk1' and `wrk2'. c m = 2*ii-1+2*(k-1)*nx wrk1(m) = wrk1(m) + ts1ct*bs1t wrk1(m+1) = wrk1(m+1) + ts1xht*bs1t wrk1(m+2) = wrk1(m+2) + ts1ct*bs2t wrk1(m+3) = wrk1(m+3) + ts1xht*bs2t wrk2(m) = wrk2(m) + ts2ct*bs1t wrk2(m+1) = wrk2(m+1) + ts2xht*bs1t wrk2(m+2) = wrk2(m+2) + ts2ct*bs2t wrk2(m+3) = wrk2(m+3) + ts2xht*bs2t 190 continue c c The vectors `wrk1' and `wrk2' contain derivatives of the c two residual equations. c dres1c = ts1c - tl1c dres1xh = ts1xh- tl1xh dres2c = ts2c - tl2c dres2xh = ts2xh- tl2xh m = 2*n-1+2*(j-1)*nx wrk1(m) = wrk1(m) + dres1c *bs1 wrk1(m+1) = wrk1(m+1) + dres1xh*bs1 wrk1(m+2) = wrk1(m+2) + dres1c *bs2 wrk1(m+3) = wrk1(m+3) + dres1xh*bs2 wrk2(m) = wrk2(m) + dres2c *bs1 wrk2(m+1) = wrk2(m+1) + dres2xh*bs1 wrk2(m+2) = wrk2(m+2) + dres2c *bs2 wrk2(m+3) = wrk2(m+3) + dres2xh*bs2 c c Add the residual x quadrature weights x shape functions to f(.). c res1 = sum1-lhs1 res2 = sum2-lhs2 f(m) = f(m) + res1*w*bs1 f(m+1) = f(m+1) + res2*w*bs1 f(m+2) = f(m+2) + res1*w*bs2 f(m+3) = f(m+3) + res2*w*bs2 c c If vector of derivatives (df) is too long, sort and eliminate c duplicates. c if (nnz.gt.ldf-4*na) then call qcksrt3(nnz,ind,df) k = ind(1) l = 1 do 200 m=2,nnz if (ind(m).eq.k) then df(l) = df(l)+df(m) else k = ind(m) l = l+1 ind(l) = ind(m) df(l) = df(m) endif 200 continue nnz = l endif if (nnz.gt.ldf-4*na) then write(*,*) 'WARNING: increase length of vector df' stop endif c c Add derivative of the residual x quadrature weights x shape c functions to df(.). NOTE: (i,j) indices switched. c m = 2*n-1+2*(j-1)*nx do 210 k=1,na if (dabs(wrk1(k)).gt.0.d0) then nnz = nnz+1 ind(nnz) = (m-1)*na+k df(nnz) = wrk1(k)*w*bs1 nnz = nnz+1 ind(nnz) = (m+1)*na+k df(nnz) = wrk1(k)*w*bs2 endif if (dabs(wrk2(k)).gt.0.d0) then nnz = nnz+1 ind(nnz) = m*na+k df(nnz) = wrk2(k)*w*bs1 nnz = nnz+1 ind(nnz) = (m+2)*na+k df(nnz) = wrk2(k)*w*bs2 endif 210 continue 160 continue 150 continue call qcksrt3(nnz,ind,df) k = ind(1) l = 1 do 220 m=2,nnz if (ind(m).eq.k) then df(l) = df(l)+df(m) else k = ind(m) l = l+1 ind(l) = ind(m) df(l) = df(m) endif 220 continue nnz = l c c Convert (j-1)*na+i stored in ind(.) into i and j c do 230 i=1,nnz ja(i) = mod(ind(i)-1,na)+1 ia(i) = (ind(i)-1)/na+1 230 continue c c Eliminate rows and columns of f,df corresponding to points c of c(.) and xh(.) preset by boundary conditions. c j = 0 iwork(1) = ibc(1) do 240 i=1,na-1 if (ibc(i).eq.0) then j = j+1 f(j) = f(i) endif iwork(i+1) = iwork(i)+ibc(i+1) 240 continue if (ibc(na).eq.0) then j = j+1 f(j) = f(na) endif j = 0 do 250 i=1,nnz if ((ibc(ia(i)).eq.0).and.(ibc(ja(i)).eq.0)) then j = j+1 df(j) = df(i) ja(j) = ja(i)-iwork(ja(i)) ia(j) = ia(i)-iwork(ia(i)) endif 250 continue nnz = j c c Use ia(.) to determine pointers to rows in sparse matrix. c ip(1) = 1 j = 1 do 260 i=2,nnz if (ia(i).ne.ia(i-1)) then j = j+1 ip(j) = i endif 260 continue ip(na-nbc+1) = nnz+1 c c Invert the derivative of f (i.e., df) with the SPARSKIT2 c iterative solver. c call ilu0(na-nbc,df,ja,ip,alu,jlu,ju,iwork,ierr) if (ierr.ne.0) then write(*,*) 'ERROR in ilu0, ierr = ',ierr stop endif eps = 1.0e-10 iout = 6 im = nm call pgmres(na-nbc,im,f,sol,vv,eps,maxit,iout,df,ja, & ip,alu,jlu,ju,ierr) if (ierr.eq.1) then write(*,*) 'WARNING: max. iterations in pgmres' endif c c Do Newton update. c sum1 = 0.d0 k = 0 do 270 j=1,ns do 270 i=1,nx if (ibc(2*i-1).eq.0) then k = k+1 c1(i,j) = c0(i,j) - sol(k) sum1 = sum1 + sol(k)*sol(k) endif if (ibc(2*i).eq.0) then k = k+1 xh1(i,j) = xh0(i,j) - sol(k) sum1 = sum1 + sol(k)*sol(k) endif 270 continue sum1 = dsqrt(sum1)/float(na-nbc) write(*,'(''The residual at iteration '',I3,'' is '',F11.9)') & it,sum1 c c Check to see if the solution is converged. If so, stop. c if (sum1.lt.crit) go to 999 do 280 j=1,ns do 280 i=1,nx c0(i,j) = c1(i,j) xh0(i,j) = xh1(i,j) 280 continue 120 continue 999 continue c c Write new input file with results. c write(7,'(F9.6,12X,''a > 0 '')') a write(7,'(F9.6,12X,''alpk in (0,1) '')') alpk write(7,'(F9.6,12X,''alph in (0,1) '')') alph write(7,'(F9.6,12X,''beta in (0,1) '')') beta write(7,'(E9.2,12X,''chi >= 0 '')') chi write(7,'(F9.6,12X,''delk in [0,1] '')') delk write(7,'(F9.6,12X,''delh in [0,1] '')') delh write(7,'(F9.6,12X,''eta >= 0 '')') eta write(7,'(F9.6,12X,''lab0 in [0,1] '')') lab0 write(7,'(F9.6,12X,''psi > 0 '')') psi write(7,'(F9.6,12X,''sigma > 0 '')') sigma write(7,'(F9.6,12X,''zeta > 0 '')') zeta write(7,'(8X,I1,12X,''evenly-spaced grid (1=yes,0=no)'')')ievenx write(7,'(8X,''1'',12X,''read initial decision fns (1=yes)'')') write(7,'(6X,I3,12X,''# of boundary constraints - c'')') nbc1 write(7,'(6X,I3,12X,''# of boundary constraints - xh'')') nbc2 write(7,*) write(7,'(F9.6,12X,''lower bound for capital grid'')') xa(1) if (ievenx.eq.0) then do 290 i=2,nxm1 write(7,1) xa(i) 290 continue endif write(7,'(F9.6,12X,''upper bound for capital grid'')') xa(nx) write(7,*) write(7,'(I4,17X,''# of policy states (given below)'')') ns write(7,*) do 300 i=1,ns write(7,3) tauxk(i),tauxh(i),tauk(i),tauh(i),tauc(i),z1(i),z2(i) 300 continue write(7,*) write(7,'(I5,16X, & ''# of nonzero elements of pi (given below)'')') nnzp write(7,*) i = 1 j = 2 do 310 k=1,nnzp if (k .eq. ipi(j)) then i = i + 1 j = j + 1 endif write(7,2) i,jpi(k),pi(k) 310 continue write(7,*) write(7,'(F12.9,9X,''consumption function, c(k,i)'')') c1(1,1) do 320 i=2,nx write(7,4) c1(i,1) 320 continue do 330 j=2,ns do 330 i=1,nx write(7,4) c1(i,j) 330 continue write(7,*) write(7,'(F12.9,9X,''investment function, xh(k,i)'')') xh1(1,1) do 325 i=2,nx write(7,4) xh1(i,1) 325 continue do 335 j=2,ns do 335 i=1,nx write(7,4) xh1(i,j) 335 continue write(7,*) k = 0 do 340 j=1,ns do 340 i=1,nx if (ibc(2*(i+(j-1)*nx)-1).eq.1) then if (k.eq.0) then write(7,'(I4,1X,I4,2X,F8.6,2X, & ''boundary conditions on c'')') i,j,c1(i,j) k = 1 else write(7,2) i,j,c1(i,j) endif endif 340 continue k = 0 do 345 j=1,ns do 345 i=1,nx if (ibc(2*(i+(j-1)*nx)).eq.1) then if (k.eq.0) then write(7,'(I4,1X,I4,2X,F8.6,2X, & ''boundary conditions on xh'')') i,j,xh1(i,j) k = 1 else write(7,2) i,j,xh1(i,j) endif endif 345 continue c c Write parameters and solution to data file for Matlab routine. c write(8,1) a write(8,1) alpk write(8,1) alph write(8,1) beta write(8,*) chi write(8,1) delk write(8,1) delh write(8,1) eta write(8,1) lab0 write(8,1) psi write(8,1) sigma write(8,1) zeta write(8,*) nx do 350 i=1,nx write(8,1) xa(i) 350 continue write(8,*) ns do 360 i=1,ns write(8,1) tauxk(i) 360 continue do 370 i=1,ns write(8,1) tauxh(i) 370 continue do 380 i=1,ns write(8,1) tauk(i) 380 continue do 390 i=1,ns write(8,1) tauh(i) 390 continue do 400 i=1,ns write(8,1) tauc(i) 400 continue do 410 i=1,ns write(8,5) z1(i) 410 continue do 420 i=1,ns write(8,5) z2(i) 420 continue write(8,*) nnzp i = 1 j = 2 do 430 k=1,nnzp if (k .eq. ipi(j)) then i = i + 1 j = j + 1 endif write(8,*) i 430 continue do 440 i=1,nnzp write(8,*) jpi(i) 440 continue do 450 i=1,nnzp write(8,1) pi(i) 450 continue do 460 j=1,ns do 460 i=1,nx write(8,1) c1(i,j) 460 continue do 470 j=1,ns do 470 i=1,nx write(8,1) xh1(i,j) 470 continue c c Format statements. c 1 format (F9.6) 2 format (I4,1X,I4,2X,F8.6) 3 format (5(F9.6,1X),2(1X,E9.2)) 4 format (F12.9) 5 format (E10.3) stop end