c PROGRAM: TWOSECT.F c c DESCRIPTION: Computation of equilibrium in two-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])*q[t]*Xh[t]} c <= (1-tau_k1[t])*r1[t]*K[t]*v[t] c +(1-tau_k2[t])*r2[t]*K[t]*(1-v[t]) c +(1-tau_h1[t])*w1[t]*l[t]*H[t]*u[t] c +(1-tau_h2[t])*w2[t]*l[t]*H[t]*(1-u[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,Xk,Xh,L/N,u,v >= 0, L/N,u,v <=1 c N[t+1] = (1+eta)* N[t] (=population) c l[t] = L[t]/N[t] (=labor input) c c Factor prices: c r1[t] = F_1(K[t],l[t]*H[t],G2[t]) c w1[t] = F_2(K[t],l[t]*H[t],G2[t]) c r2[t] = G_1(K[t],l[t]*H[t],G3[t])*q[t] c w2[t] = G_2(K[t],l[t]*H[t],G3[t])*q[t] c c Gov't spending: c Gj[t] = zj[t] * F(K[t],l[t]*H[t],G2[t]) c c Resource constraints: c C[t] + Xk[t] + G1[t] + G2[t] + G3[t] c = F(K[t],l[t]*H[t],G2[t]) c Xh[t] = G(K[t],l[t]*H[t],G3[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), v=v(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 v[t]. The remaining decision c variables are obtained by way of the static first-order c 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 G(k,l,g) = b*k^thet_k*l^thet_h*g^(1-thet_k-thet_h) c c TEST CASE: U(c,l,g) = log(c) c F(k,l,g) = a*k^alpha*l^(1-alpha) c G(k,l,g) = a*k^alpha*l^(1-alpha) c and del_k= del_h= 1 c ==> normalized variables all constant c c INPUT FILES: (a) twosect.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 b parameter in G(.) real > 0 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 thetk parameter in G(.) real in [0,1] c theth parameter in G(.) real in [0,1] 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 v 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 10 real matrix c tauk1,tauk2, ... c tauh1,tauh2, ... c tauc,z1,z2,z3] 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 v initial guess for nx*ns x 1 real c capital allocation vector; read xh(:,1) c in sector 1 first, xh(:,2) c (if initialize=1) 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,v points v(x(i),j) triplets of (integer, c that are to be fixed integer,real) c during computation c c Notes: c i. see twosect.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) twosect.nxt: new input file with current solution c (c) twosect.dat: input file for Matlab post-processing c routine -- type twosect 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-23-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 c parameter (nx=11,ns=800,na=17600,nq=5,nm=40,ldf=1500000,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), & v0(nx,ns), v1(nx,ns), tauxh(ns), tauxk(ns), tauk1(ns), & tauk2(ns), tauh1(ns), tauh2(ns), tauc(ns), z1(ns), & z2(ns), z3(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='twosect.inp') open(unit=7, file='twosect.nxt') open(unit=8, file='twosect.dat') c c Read in the parameters. c read(5,*) a read(5,*) alpk read(5,*) alph read(5,*) b 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,*) thetk read(5,*) theth 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 thet = thetk+theth thet1 = 1.d0-thet theti = 1.d0/thet 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),tauk1(i),tauk2(i),tauh1(i), & tauh2(i),tauc(i),z1(i),z2(i),z3(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,*) v0(i,j) 65 continue else c c If initial guess is not given, set c, v as in test case c do 70 j=1,ns do 70 i=1,nx c0(i,j) = (1.d0-beta)*(alpk/(1.d0-alpk))**alpk v0(i,j) = 1.d0-beta*(1.d0-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 v0(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;v0) = int R(x;c0,v0) 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 capital allocation decision using c a finite element approximation. c c = c0(n,j)*bs1 + c0(n+1,j)*bs2 v = v0(n,j)*bs1 + v0(n+1,j)*bs2 c c Compute intermediate parameters. c tauxk1= 1.d0+tauxk(j) tauxh1= 1.d0+tauxh(j) tauk11= 1.d0-tauk1(j) tauk21= 1.d0-tauk2(j) tauh11= 1.d0-tauh1(j) tauh21= 1.d0-tauh2(j) tauc1 = 1.d0+tauc(j) phi = tauk11*tauh21*alpk*theth/(tauk21*tauh11* & alph*thetk) uu = 1.d0/(1.d0-phi+phi/v) tem1 = psi*c tem2 = tauh11*alph*aa*(x*v/uu)**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 u = (c*lei**psi*g1**zeta)**sigm1/sigm1 duc = sigm1*u/c ys = a*(x*v)**alpk*(lab*uu)**alph g2 = (z2(j)*ys)**alpi y = ys*g2**alp1 g1 = z1(j)*y g3 = z3(j)*y tem1 = dmax1(1.d0-v,0.d0) tem2 = dmax1(1.d0-uu,0.d0) xh = b*(x*tem1)**thetk*(lab*tem2)**theth* & g3**thet1 gam = delh1+xh gams = gam**sigz xk = y-g1-g2-g3-c q = tauk11*alpk*y*(1.d0-v) / (tauk21*thetk*xh*v) 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*q*tauxh1/tauc1-pxh2) c c Compute total derivatives of lhs1 and lhs2 in terms of unknowns c (i.e., tl1c, tl1v, tl2c, tl2v). c c NOTE: c dU = fn(c,l,g1) c y = fn(k,v,l,u) c g1 = fn(y) c xh = fn(k,v,l,u) c q = fn(k,v,l,u) c xk = fn(y,c) c u = fn(v) c l = fn(k,v,u,c) c lhs1 = fn(xk,xh,c,l,g1) c lhs2 = fn(q,xh,c,l,g1) c want: dlhs/dc, dlhs/dv c d2uc = -sigma*duc/c ducl = -sigm1*psi*duc/lei ducg = sigm1*zeta*duc/g1 dyv = alphak*y/v dyl = alphah*y/lab dyu = alphah*y/uu dg1y = z1(j) dxhv = -thetk*xh/(1.d0-v)+alphak*thet1*xh/v dxhl = (theth+alphah*thet1)*xh/lab dxhu = -theth*xh/(1.d0-uu)+alphah*thet1*xh/uu dqv = (alphak*thet-1.d0)*q/v-(1.d0-thetk)*q/(1.d0-v) dql = (alphah*thet-theth)*q/lab dqu = alphah*thet*q/uu+theth*q/(1.d0-uu) dxkc = -1.d0 dxky = 1.d0-z1(j)-z2(j)-z3(j) duv = phi/((phi+(1.d0-phi)*v)*(phi+(1.d0-phi)*v)) tem = tauh11*aa*alph*z2(j)**alphai/tauc1* & (x*v/lab/uu)**alphak ddc = psi ddv = -tem*alphak*lei/v ddl = tem*(alphak*lei/lab+1.d0) ddu = tem*alphak*lei/uu dlv = -ddv/ddl dlu = -ddu/ddl dlc = -ddc/ddl dl1xk = -gams*2.d0*chi*dmin1(xk,0.d0) dl1xh = sigz*gams/gam*(duc*tauxk1/tauc1-pxk2) dl1c = gams*d2uc*tauxk1/tauc1 dl1l = gams*ducl*tauxk1/tauc1 dl1g1 = gams*ducg*tauxk1/tauc1 dl2q = gams*duc*tauxh1/tauc1 dl2xh = sigz*gams/gam*(duc*q*tauxh1/tauc1-pxh2)- & gams*2.d0*chi*dmin1(xh,0.d0) dl2c = gams*d2uc*q*tauxh1/tauc1 dl2l = gams*ducl*q*tauxh1/tauc1 dl2g1 = gams*ducg*q*tauxh1/tauc1 tuv = duv tlc = dlc tlv = dlv+dlu*tuv tyc = dyl*tlc tyv = dyv+dyl*tlv+dyu*tuv tg1c = dg1y*tyc tg1v = dg1y*tyv txkc = dxkc+dxky*tyc txkv = dxky*tyv tqc = dql*tlc tqv = dqv+dql*tlv+dqu*tuv txhc = dxhl*tlc txhv = dxhv+dxhl*tlv+dxhu*tuv tl1c = dl1c+dl1xk*txkc+dl1xh*txhc+dl1l*tlc+dl1g1*tg1c tl1v = dl1xk*txkv+dl1xh*txhv+dl1l*tlv+dl1g1*tg1v tl2c = dl2c+dl2q*tqc+dl2xh*txhc+dl2l*tlc+dl2g1*tg1c tl2v = dl2q*tqv+dl2xh*txhv+dl2l*tlv+dl2g1*tg1v 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-xt*txhc)/gam txtv = (txkv-xt*txhv)/gam c c Initialize variables for summation. c sum1 = 0.d0 sum2 = 0.d0 ts1c = 0.d0 ts2c = 0.d0 ts1v = 0.d0 ts2v = 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 capital allocation c decisions using a finite element approximation. c ct = c0(ii,k)*bs1t + c0(ii+1,k)*bs2t vt = v0(ii,k)*bs1t + v0(ii+1,k)*bs2t c c Compute intermediate parameters. c tauxk2 = 1.d0+tauxk(k) tauxh2 = 1.d0+tauxh(k) tauk12 = 1.d0-tauk1(k) tauk22 = 1.d0-tauk2(k) tauh12 = 1.d0-tauh1(k) tauh22 = 1.d0-tauh2(k) tauc2 = 1.d0+tauc(k) phi = tauk12*tauh22*alpk*theth/(tauk22*tauh12* & alph*thetk) uut = 1.d0/(1.d0-phi+phi/vt) tem1 = psi*ct tem2 = tauh12*alph*aa*(xt*vt/uut)**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*vt)**alpk*(labt*uut)**alph g2t = (z2(k)*ys)**alpi yt = ys*g2t**alp1 g1t = z1(k)*yt g3t = z3(k)*yt xkt = yt-g1t-g2t-g3t-ct tem1 = dmax1(1.d0-vt,0.d0) tem2 = dmax1(1.d0-uut,0.d0) xht = b*(xt*tem1)**thetk*(labt*tem2)** & theth*g3t**thet1 qt = tauk12*alpk*yt*(1.d0-vt)/ & (tauk22*thetk*xht*vt) 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+ & tauk12*alpk*yt/xt/vt)/tauc2 -pxkt2*delk1) sum2 = sum2 + betn*pi(l)* (dutc*qt*(tauxh2* & delh1+tauh22*theth*xht/(1.d0-uut))/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, ts1v, ts1ct, ts1vt, ts2c, c ts2v, ts2ct, ts2vt). c c NOTE: c dU = fn(c,l,g1) c y = fn(k,v,l,u) c g1 = fn(y) c xh = fn(k,v,l,u) c q = fn(k,v,l,u) c xk = fn(y,c) c u = fn(v) c l = fn(k,v,u,c) c sum1 = fn(c,l,g1,y,k,v,xk) c sum2 = fn(c,l,g1,q,xh,u,xk) c want: dsum/dc, dsum/dv (for c,v today and tomorrow) c d2utc = -sigma*dutc/ct dutcl = -sigm1*psi*dutc/leit dutcg = sigm1*zeta*dutc/g1t dytxt = alphak*yt/xt dytvt = alphak*yt/vt dytlt = alphah*yt/labt dytut = alphah*yt/uut dg1tyt = z1(k) dxhtxt = (thetk+alphak*thet1)*xht/xt dxhtvt = -thetk*xht/(1.d0-vt)+alphak*thet1*xht/vt dxhtlt = (theth+alphah*thet1)*xht/labt dxhtut = -theth*xht/(1.d0-uut)+alphah*thet1*xht/uut dqtxt = (alphak*thet-thetk)*qt/xt dqtvt = (alphak*thet-1.d0)*qt/vt-(1.d0-thetk)*qt/ & (1.d0-vt) dqtlt = (alphah*thet-theth)*qt/labt dqtut = alphah*thet*qt/uut+theth*qt/(1.d0-uut) dxktct = -1.d0 dxktyt = 1.d0-z1(k)-z2(k)-z3(k) dutvt = phi/((phi+(1.d0-phi)*vt)* & (phi+(1.d0-phi)*vt)) tem = tauh12*aa*alph*z2(k)**alphai/tauc2* & (xt*vt/labt/uut)**alphak ddc = psi ddx = -tem*alphak*leit/xt ddv = -tem*alphak*leit/vt ddl = tem*(alphak*leit/labt+1.d0) ddu = tem*alphak*leit/uut dltxt = -ddx/ddl dltvt = -ddv/ddl dltut = -ddu/ddl dltct = -ddc/ddl tem = (tauxk2*delk1+tauk12*alpk*yt/xt/vt)/tauc2 dr1ct = d2utc*tem dr1lt = dutcl*tem dr1g1t = dutcg*tem dr1yt = dutc*tauk12*alpk/(xt*vt*tauc2) dr1xt = -dutc*tauk12*alpk*yt/(xt*xt*vt*tauc2) dr1vt = -dutc*tauk12*alpk*yt/(xt*vt*vt*tauc2) dr1xkt = -2.d0*chi*dmin1(xkt,0.d0)*delk1 tem = (tauxh2*delh1+tauh22*theth*xht/(1.d0- & uut))/tauc2 dr2ct = d2utc*qt*tem dr2lt = dutcl*qt*tem dr2g1t = dutcg*qt*tem dr2qt = dutc*tem dr2xht = dutc*qt*tauh22*theth/(1.d0-uut)/tauc2- & 2.d0*chi*dmin1(xht,0.d0)*delh1+pxht2* & (sigm1*zeta1-3.d0) dr2ut = dutc*qt*tauh22*theth*xht/(1.d0-uut)/ & (1.d0-uut)/tauc2 dr2xkt = pxkt2*(sigm1*zeta1-3.d0) tctxt = -(c0(ii,k) -c0(ii+1,k))/dxt tvtxt = -(v0(ii,k) -v0(ii+1,k))/dxt tutvt = dutvt tutxt = dutvt*tvtxt tltxt = dltxt+dltct*tctxt+dltvt*tvtxt+dltut*tutxt tltct = dltct tltvt = dltvt+dltut*tutvt tytxt = dytxt+dytvt*tvtxt+dytlt*tltxt+dytut*tutxt tytct = dytlt*tltct tytvt = dytvt+dytlt*tltvt+dytut*tutvt tg1txt = dg1tyt*tytxt tg1tct = dg1tyt*tytct tg1tvt = dg1tyt*tytvt txktxt = dxktyt*tytxt+dxktct*tctxt txktct = dxktct+dxktyt*tytct txktvt = dxktyt*tytvt tqtxt = dqtxt+dqtvt*tvtxt+dqtlt*tltxt+dqtut*tutxt tqtct = dqtlt*tltct tqtvt = dqtvt+dqtlt*tltvt+dqtut*tutvt txhtxt = dxhtxt+dxhtvt*tvtxt+dxhtlt*tltxt+dxhtut* & tutxt txhtct = dxhtlt*tltct txhtvt = dxhtvt+dxhtlt*tltvt+dxhtut*tutvt tr1xt = dr1ct*tctxt+dr1lt*tltxt+dr1g1t*tg1txt+ & dr1yt*tytxt+dr1xt+dr1vt*tvtxt+dr1xkt* & txktxt tr2xt = dr2ct*tctxt+dr2lt*tltxt+dr2g1t*tg1txt+ & dr2qt*tqtxt+dr2xht*txhtxt+dr2ut*tutxt+ & dr2xkt*txktxt ts1c = ts1c +betn*pi(l)*tr1xt*txtc ts1v = ts1v +betn*pi(l)*tr1xt*txtv ts1ct = betn*pi(l)*(dr1ct+dr1lt*tltct+dr1g1t* & tg1tct+dr1yt*tytct+dr1xkt*txktct) ts1vt = betn*pi(l)*(dr1lt*tltvt+dr1g1t* & tg1tvt+dr1yt*tytvt+dr1vt+dr1xkt*txktvt) ts2c = ts2c +betn*pi(l)*tr2xt*txtc ts2v = ts2v +betn*pi(l)*tr2xt*txtv ts2ct = betn*pi(l)*(dr2ct+dr2lt*tltct+dr2g1t* & tg1tct+dr2qt*tqtct+dr2xht*txhtct+dr2xkt* & txktct) ts2vt = betn*pi(l)*(dr2lt*tltvt+dr2g1t*tg1tvt+ & dr2qt*tqtvt+dr2xht*txhtvt+dr2ut*tutvt+ & dr2xkt*txktvt) 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) + ts1vt *bs1t wrk1(m+2) = wrk1(m+2) + ts1ct *bs2t wrk1(m+3) = wrk1(m+3) + ts1vt *bs2t wrk2(m) = wrk2(m) + ts2ct *bs1t wrk2(m+1) = wrk2(m+1) + ts2vt *bs1t wrk2(m+2) = wrk2(m+2) + ts2ct *bs2t wrk2(m+3) = wrk2(m+3) + ts2vt *bs2t 190 continue c c The vectors `wrk1' and `wrk2' contain derivatives of the c two residual equations. c dres1c = ts1c - tl1c dres1v = ts1v - tl1v dres2c = ts2c - tl2c dres2v = ts2v - tl2v m = 2*n-1+2*(j-1)*nx wrk1(m) = wrk1(m) + dres1c *bs1 wrk1(m+1) = wrk1(m+1) + dres1v *bs1 wrk1(m+2) = wrk1(m+2) + dres1c *bs2 wrk1(m+3) = wrk1(m+3) + dres1v *bs2 wrk2(m) = wrk2(m) + dres2c *bs1 wrk2(m+1) = wrk2(m+1) + dres2v *bs1 wrk2(m+2) = wrk2(m+2) + dres2c *bs2 wrk2(m+3) = wrk2(m+3) + dres2v *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 v(.) 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 v1(i,j) = v0(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) v0(i,j) = v1(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,''b > 0 '')') b 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,''thetk in (0,1) '')') thetk write(7,'(F9.6,12X,''theth in (0,1) '')') theth 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 - v'')') 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),tauk1(i),tauk2(i),tauh1(i), & tauh2(i),tauc(i),z1(i),z2(i),z3(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,''fraction of K in sector 1, v(k,i)'')')v1(1,1) do 325 i=2,nx write(7,4) v1(i,1) 325 continue do 335 j=2,ns do 335 i=1,nx write(7,4) v1(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 v'')') i,j,v1(i,j) k = 1 else write(7,2) i,j,v1(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) b 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) thetk write(8,1) theth 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) tauk1(i) 380 continue do 385 i=1,ns write(8,1) tauk2(i) 385 continue do 390 i=1,ns write(8,1) tauh1(i) 390 continue do 395 i=1,ns write(8,1) tauh2(i) 395 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 do 425 i=1,ns write(8,5) z3(i) 425 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) v1(i,j) 470 continue c c Format statements. c 1 format (F9.6) 2 format (I4,1X,I4,2X,F8.6) 3 format (2(F9.6,1X),5(F5.2,1X),3(E8.2,1X)) 4 format (F12.9) 5 format (E10.3) stop end