c PROGRAM: TAXCHNGE.F c c DESCRIPTION: Program that accompanies THE OPTIMUM QUANTITY OF c DEBT by Aiyagari and McGrattan. The optimum c quantity of debt is computed for a general equilibrium c model with borrowing constraints and precautionary c saving. The model an augmented version of Aiyagari's c model (QJE 1994) -- augmented to permit growth c and to include government debt, taxes, and government c consumption. This version has a different tax rate on c labor and capital, TFP as an input, and total government c consumption, transfers, and debt as inputs (rather than c ratios with respect to GDP). c c PROBLEM: Step 1. find a~=A(a,e), that solves the functional c equation: c c R(a,e) = beta~ sum_j pi(i,j){(1+rbar)u'(c~,l~)+s}- c (1+grate)*u'(c,l) c = 0 c c where c c = w*e*(1-l) +(1+rbar)a-T-A(a,e) c c~ = w*e~*(1-l~)+(1+rbar)A(a,e)-T-A(a~,e~) c T = g+(rbar-grate)*b-taxes c l = ith value of Markov chain c l~ = jth value of Markov chain c pi = the transition matrix for l c w,r = after tax wage and interest rates c s = 3*zeta*d*(d+|d|), d = -a c c Step 2. find r such that E[a(r)] = K(r)/Y(r)+b where c E[a(r)] is the mean of the distribution of assets, c K(r) is the capital stock, Y(r) is output, and c K and r satisfy c c r = alpha*(K/L)^(alpha-1)-delta, rbar=r*(1-tauk) c c ALGORITHM: the finite element method with piecewise shape c functions N(x) to compute A(a,l) and a bisection c method to find the equilibrium interest rate (r). c c MODEL: the discrete-time Lucas asset-pricing model with c an individual agent choosing c[t],a[t+1] to maximize c c E sum_t beta^t u(c[t],l[t]) c c subject to c c c[t] + a[t+1] = (1+rbar)*a[t] + w*e[t]*(1-l[t])-T[t] c a[t]>= 0 c T[t] = G[t]+(1+rbar)*B[t]-B[t+1]- c tauk*{f(K[t],L[t])-w*L[t]-delta*K[t]} c taul*{w*L[t]} c e[t]: Markov chain w/ transition matrix pi c c Assuming that c[t], G[t], B[t], T[t], etc. are growing c at same rate (e.g., grate), we can divide all by GNP c or Y[0]*(1+grate)^t for each date t: c c E sum_t (beta~)^t {u(c~[t],l[t])-zeta*(d[t]^3+|d[t]^3|)} c |--> penalty function c subject to c c c~[t] + (1+grate)*a~[t+1] c = (1+rbar)*a~[t]+w*e[t]*(1-l[t])-tax c beta~ = beta*(1+grate)^(eta*(1-mu)) c d[t] = -a~[t] c tax = g+(rbar-grate)*b-tauk*(1-w*L/Y-delta*K/Y) c -taul*w*L/Y, Y=f(K,L) c c and in equilibrium c E a~[t] = K/f(K,L) + b c rbar = (1-tauk)*(f_1(K,L)-delta) c w = (1-taul)*f_2(K,L)/f(K,L) c c FUNCTIONS: u(c,l) = [c^eta*l^(1-eta)]^(1-mu)/(1-mu) c f(K,L) = K^alpha*L^(1-alpha) c c TEST CASE: Stokey and Lucas (1989), ex. 5.17. (See example1.tex.) c c INPUTS: (a) optimal.inp: c alpha c beta c delta c eta c g c grate c mu c transfer c lstax c k0 c ncases c edu(i) i=1,ns c pi(i,j) i=1,ns j=1,ns c xf(i) i=1,nf c repeat | b,bndl,bndu,nbc c ncases+1 | xa(i),a0(i,j) i=1,na j=1,ns c times | c c where c alpha,beta,delta,eta,g,grate,mu,transfer are c parameters defined above c lstax: 1 if lump-sum tax or 0 if proportional c k0: the penalty function for enforcing asset c bound uses zeta=10^k0 c ncases: number of debt levels being compared c to benchmark c edu: idiosyncratic productivity levels c pi: transition probabilities for edu c xf's: grid pts. for the F[t]'s c b: debt/GDP ratio c bndl,bndu: upper and lower bounds for r c nbc: number of boundary conditions c xa's: grid pts. for the a~[t]'s c a0's: initial guesses for asset holdings c indices: i,j pairs such that a0(i,j) stays fixed c na: number of points in the grid for a~[t] c nf: number of points in the grid for F[t] c ns: number of values that edu takes c c NOTE: if the user changes the grid size, he c must edit this program and change na, c nf,ns,nas=na*ns,and nfs=nf*ns with nf>na c c OUTPUTS: (b) optimal.out: general results c (c) optimal.nxt: new input file c c SUBROUTINES: (a) SPARSKIT2: ilut,matvec c (b) other: qgausl,qcksrt3 c c Ellen McGrattan, 9-5-93 c Revised, 10-29-97 ERM c parameter (na=61,nf=101,ns=7,nas=427,nfs=707,nxm=3,nt=500000, & ldh=500000, nm=60) implicit double precision (a-h,o-z) double precision h(nfs), dh(ldh), work(nt), xa(na), xf(nf), & absa(nxm,na-1), wgta(nxm,na-1), & absf(nxm,nf-1), wgtf(nxm,nf-1), a0(na,ns), a1(na,ns), & edu(ns), pi(ns,ns), mu, mu1, lei0, lei, leit, lhs, & ai(nf,ns), f0(nf,ns), f1(nf,ns), atem(nf,ns), & alu(ldh), vv(nfs,nm+1), sol(nfs), chk(100,6) integer iwork(nt), nxp(nf-1), ibc(nas), & ind(ldh), ia(ldh), ja(ldh), ip(nfs+1), jlu(ldh), ju(nfs) c c Stop if the number of grid points for A exceeds that of F. c if (nf.lt.na) then write(*,*) 'WARNING: nf must be greater than or equal to na.' stop endif c c Open input and output files c open(unit=5, file='taxchnge.inp') open(unit=7, file='taxchnge.out') open(unit=8, file='taxchnge.nxt') open(unit=9, file='taxchnge.tmp') c c Read in parameters for benchmark, initial guess for A(.,.) and grid c read(5,*) alpha read(5,*) beta read(5,*) delta read(5,*) eta read(5,*) g read(5,*) grate read(5,*) mu read(5,*) transfer read(5,*) tauk read(5,*) tfp read(5,*) lstax read(5,*) k0 read(5,*) ncases read(5,*) (edu(i),i=1,ns) do 10 i=1,ns read(5,*) (pi(i,j),j=1,ns) 10 continue read(5,*) (xf(i),i=1,nf) c c General results written to unit 7 (assumes b is varied across cases). c write(7,*) 'RESULTS FOR OPTIMAL.F' write(7,*) write(7,*) ' Economic parameters that remain fixed' write(7,'('' alpha = '',F12.6,'' '')') alpha write(7,'('' beta = '',F12.6,'' '')') beta write(7,'('' delta = '',F12.6,'' '')') delta write(7,'('' eta = '',F12.6,'' '')') eta write(7,'('' g = '',F12.6,'' '')') g write(7,'('' grate = '',F12.6,'' '')') grate write(7,'('' mu = '',F12.6,'' '')') mu write(7,'('' transfer = '',F12.6,'' '')') transfer write(7,'('' tauk = '',F12.6,'' '')') tauk write(7,'('' tfp = '',F12.6,'' '')') tfp write(7,*) write(7,*) ' Markov chain for education ' write(7,*) ' values that education take on: ' write(7,3) (edu(i),i=1,ns) write(7,*) ' transition probabilities: ' do 20 i=1,ns write(7,3) (pi(i,j),j=1,ns) 20 continue write(7,*) write(7,8) ' Debt Welfare E[a~] Int. Rate', & ' Tax Rate Value Fn. Tot. Hours Output' write(*,8) ' Debt Welfare E[a~] Int. Rate', & ' Tax Rate Value Fn. Tot. Hours Output' write(9,8) ' Debt Welfare E[a~] Int. Rate', & ' Tax Rate Value Fn. Tot. Hours Output' c c Write parameters and solutions to the next input file. c write(8,'('' '',F12.8,'' alpha '')') alpha write(8,'('' '',F12.8,'' beta '')') beta write(8,'('' '',F12.8,'' delta '')') delta write(8,'('' '',F12.8,'' eta '')') eta write(8,'('' '',F12.8,'' g '')') g write(8,'('' '',F12.8,'' grate '')') grate write(8,'('' '',F12.8,'' mu '')') mu write(8,'('' '',F12.8,'' transfer'')') transfer write(8,'('' '',F12.8,'' tauk '')') tauk write(8,'('' '',F12.8,'' tfp '')') tfp write(8,'('' '',I12 ,'' lump sum tax?'')')lstax write(8,'('' '',I12 ,'' k0 '')') k0 write(8,'('' '',I12 ,'' ncases '')') ncases write(8,*) write(8,1) (edu(i),i=1,ns) do 30 i=1,ns write(8,1) (pi(i,j),j=1,ns) 30 continue write(8,*) i1 = mod(nf,7) j1 = nf/7 write(8,*) do 35 i=1,j1 write(8,4) (xf(j),j=(i-1)*7+1,i*7) 35 continue write(8,4) (xf(j),j=nf-i1+1,nf) c c Other parameters that are set inside the program c * nxp(i), i=1,nx are the # of quadrature points per interval c * stopping criteria (crit*) c * maximum number of interations for various loops (maxit*) c do 40 i=1,nf-1 nxp(i) = 3 40 continue crit1 = 0.0000001d0 crit2 = 0.000000001d0 crit3 = 0.00000000001d0 maxit1 = 20 maxit2 = 50 maxit3 = 10 c c Initialize the cumulutive distribution function for asset holdings (f0). c do 50 i=1,ns work(i) = 1.d0/float(ns) 50 continue do 60 j=1,500 do 70 i=1,ns sum1 = 0.d0 do 80 k=1,ns sum1 = sum1 + pi(k,i)*work(k) 80 continue work(i+ns) = sum1 70 continue do 90 i=1,ns work(i) = work(i+ns) work(i+ns) = 0.d0 90 continue 60 continue do 100 i=1,nf do 100 j=1,ns f0(i,j) = work(j) 100 continue c c Set intermediate parameters. c alph1 = alpha-1.d0 zeta = float(10**k0)-1.d0 eta1 = 1.d0-eta grat1 = 1.d0+grate mu1 = 1.d0-mu betg = beta*grat1**(eta*mu1) bet1 = 1.d0-betg taul = 0.d0 tau1 = 1.d0 tax = -transfer c c Iterate over the alternative cases (i.e. different levels of debt/GDP). c do 110 it=1,ncases+1 read(5,*) b,bndl,bndu,nbc c c Initialize or reinitialize certain variables. c k = 0 do 120 j=1,ns do 120 i=1,na k = k+1 ibc(k) = 0 a1(i,j) = 0.d0 120 continue taul = 0.d0 tau1 = 1.d0 tax = -transfer bndu0 = bndu bndl0 = bndl lei0 = eta c c Continue reading. c do 130 i=1,na read(5,*) xa(i),(a0(i,j),j=1,ns) 130 continue read(5,*) (iwork(i),i=1,2*nbc) do 140 i=1,nbc j = iwork(2*i-1) k = iwork(2*i) ibc(j+(k-1)*na) = 1 140 continue c c Set abscissas and weights for quadrature on elements. c do 150 i=1,na-1 call qgausl (nxp(i), xa(i), xa(i+1), absa(1,i), wgta(1,i)) 150 continue do 160 i=1,nf-1 call qgausl (nxp(i), xf(i), xf(i+1), absf(1,i), wgtf(1,i)) 160 continue c c Iterate over the aggregate labor supply until equilibrium found. c r = 0.5*(bndu+bndl) if (lstax.eq.0) then taul = (g-tax+r*b-grate*b-tauk*r*(b+alpha/(r+delta)))/ & (-alph1) tau1 = 1.d0-taul endif aggli = 1.d0/(1.d0+eta1*(1.d0-g-alpha*(grate+delta)/ & (r+delta))/(-eta*tau1*alph1)) del = 0.0001*aggli do 170 it1=1,maxit1 do 180 it0=0,1 bndu = bndu0 bndl = bndl0 aggl = aggli+float(it0)*del c c Iterate over the after-tax interest rate until equilibrium found. c do 190 it2=1,maxit2 nx = na-1 c c Set or reset the interest rate and all parameters that c depend on it. c r = 0.5*(bndu+bndl) c c Adjust aggregates by output, which depends on r,aggl c y = (tfp*(alpha/(r+delta))**alpha)**(-1/alph1)*aggl aggk = alpha/(r+delta)*y c c Set taxes and prices c if (lstax.eq.1) then tax = g+(r-grate)*b else taul = (g-tax+r*b-grate*b-tauk*r*(b+aggk))/(-alph1*y) tau1 = 1.d0-taul endif rbar = r*(1.d0-tauk) w = -tau1*alph1*y/aggl r1 = 1.d0+rbar amn = dmax1(0.d0,(-w*edu(1)+tax)/(rbar-grate)) c c The Newton-Raphson iteration, i.e., find x such that c h(x)=0 where x are coefficients of the finite element c approximation. c do 200 it3=1,maxit3 do 210 i=1,nfs h(i) = 0.d0 210 continue do 220 i=1,ldh dh(i) = 0.d0 220 continue nnz = 0 c c Construct h(a0) = int R(a,l;a0) dx and its derivative. c do 230 j=1,ns do 240 n=1,nx x1 = xa(n) x2 = xa(n+1) c c Compute h at all quadrature points on the element. c do 250 i=1,nxp(n) x = absa(i,n) xwgt = wgta(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 asset holdings using finite element approximation. c a = a0(n,j)*bs1+a0(n+1,j)*bs2 c c Back out other variables from first-order conditions c and definitions. c lei = lei0 do 260 k=1,15 c = r1*x+w*edu(j)*(1.d0-lei)-tax-grat1*a u = (c**eta*lei**eta1)**mu1 d = dmin1(1.d0-lei,0.d0) f = eta1*u/lei-eta*u*w*edu(j)/c-zeta*d*d fj = -2.d0*eta1*eta*mu1*u/lei/c*w*edu(j)+ & eta1*(eta1*mu1-1.d0)*u/lei/lei+eta* & (eta*mu1-1.d0)*u*w*w*edu(j)*edu(j)/c/c+ & 2.d0*zeta*d lei = lei-f/fj 260 continue c c Compute left-hand side of dynamic first-order condition. c lhs = eta*mu1*u/c*grat1/mu1 c c Compute total derivative of lhs with respect to c unknowns (tla). c dfa = -grat1*(eta1*eta*mu1*u/lei/c-eta*(eta* & mu1-1.d0)*u*w*edu(j)/c/c) dla = -dfa/fj dca = -w*edu(j)*dla-grat1 dua = eta*u/c*dca+eta1*u/lei*dla tla = (dua-u/mu1/c*dca)*eta*mu1*grat1/c c c Find next period's asset holdings on grid. c xt = a ii = 1 do 270 k=1,nx if (xt.ge.xa(k)) ii = k 270 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) c c Compute sum_j pi(i,j){(1+r)u'(c~)+s} and its derivatives. c rhs = 0.d0 tra = 0.d0 do 280 k=1,nas work(k) = 0.d0 280 continue do 290 k=1,ns at = a0(ii,k)*bs1t+a0(ii+1,k)*bs2t leit = lei0 do 300 l=1,15 ct = r1*xt+w*edu(k)*(1.d0-leit)-tax- & grat1*at ut = (ct**eta*leit**eta1)**mu1 d = dmin1(1.d0-leit,0.d0) f = eta1*ut/leit-eta*ut*w*edu(k)/ct- & zeta*d*d fj = -2.d0*eta1*eta*mu1*ut/leit/ct*w* & edu(k)+eta1*(eta1*mu1-1.d0)*ut/ & leit/leit+eta*(eta*mu1-1.d0)*ut* & w*w*edu(k)*edu(k)/ct/ct+2.d0*zeta*d leit = leit-f/fj 300 continue dftct = -eta1*eta*mu1*ut/leit/ct+eta* & (eta*mu1-1.d0)*ut*w*edu(k)/ct/ct dltat = -grat1*dftct/fj dctat = -w*edu(k)*dltat-grat1 datxt = (-a0(ii,k)+a0(ii+1,k))/dxt dltxt = -dftct*(grat1*datxt-r1)/fj dctxt = -w*edu(k)*dltxt+r1-grat1*datxt dutat = betg*ut*(eta*dctat/ct+eta1*dltat/ & leit) dutxt = betg*ut*(eta*dctxt/ct+eta1*dltxt/ & leit) ut = betg*ut/mu1 d = dmin1(xt-amn,0.d0) st = zeta*d*d dstat = 2.d0*zeta*d rhs = rhs + pi(j,k)*(eta*mu1*r1*ut/ct+ & betg*st) tra = tra + pi(j,k)*(eta*mu1*r1*(dutxt- & ut/ct*dctxt)/ct+betg*dstat) trat = pi(j,k)*(eta*mu1*r1*(dutat-ut/ct* & dctat)/ct) m = ii+(k-1)*na work(m) = work(m) + trat*bs1t work(m+1) = work(m+1) + trat*bs2t 290 continue c c Residual = beta~ sum_j pi(i,j){(1+r)u'(c~,l~)+s} c -(1+grate)*u'(c,l). c res = rhs-lhs c c Add the derivatives of the residual to the work vector. c dres = tra-tla m = n+(j-1)*na work(m) = work(m) + dres*bs1 work(m+1) = work(m+1) + dres*bs2 c c Add the residual x quadrature weights x shape c functions to h(.). c res = rhs-lhs h(m) = h(m) + res*xwgt*bs1 h(m+1) = h(m+1) + res*xwgt*bs2 c c If vector of derivatives (dh) is too long, sort c and eliminate duplicates. c if (nnz.gt.ldh-2*nas) then call qcksrt3(nnz,ind,dh) k = ind(1) l = 1 do 310 m=2,nnz if (ind(m).eq.k) then dh(l) = dh(l)+dh(m) else k = ind(m) l = l+1 ind(l) = ind(m) dh(l) = dh(m) endif 310 continue nnz = l endif if (nnz.gt.ldh-2*nas) then write(*,*) 'WARNING: increase length of dh' write(9,*) 'WARNING: increase length of dh' stop endif c c Add derivative of the residual x quadrature c weights x shape functions to dh(.). NOTE: (i,j) c indices switched. c m = n+(j-1)*na do 320 k=1,nas if (dabs(work(k)).gt.0.d0) then nnz = nnz+1 ind(nnz) = (m-1)*nas+k dh(nnz) = work(k)*xwgt*bs1 nnz = nnz+1 ind(nnz) = m*nas+k dh(nnz) = work(k)*xwgt*bs2 endif 320 continue 250 continue 240 continue 230 continue call qcksrt3(nnz,ind,dh) k = ind(1) l = 1 do 330 m=2,nnz if (ind(m).eq.k) then dh(l) = dh(l)+dh(m) else k = ind(m) l = l+1 ind(l) = ind(m) dh(l) = dh(m) endif 330 continue nnz = l c c Convert (j-1)*nas+i stored in ind(.) into i and j c do 340 i=1,nnz ja(i) = mod(ind(i)-1,nas)+1 ia(i) = (ind(i)-1)/nas+1 340 continue c c Eliminate rows and columns of h, dh corresponding to points c of a(.) preset by boundary conditions. c j = 0 iwork(1) = ibc(1) do 350 i=1,nas-1 if (ibc(i).eq.0) then j = j+1 h(j) = h(i) endif iwork(i+1) = iwork(i)+ibc(i+1) 350 continue if (ibc(nas).eq.0) then j = j+1 h(j) = h(nas) endif j = 0 do 360 i=1,nnz if ((ibc(ia(i)).eq.0).and.(ibc(ja(i)).eq.0)) then j = j+1 dh(j) = dh(i) ja(j) = ja(i)-iwork(ja(i)) ia(j) = ia(i)-iwork(ia(i)) endif 360 continue nnz = j c c Use ia(.) to determine pointers to rows in sparse matrix. c ip(1) = 1 j = 1 do 370 i=2,nnz if (ia(i).ne.ia(i-1)) then j = j+1 ip(j) = i endif 370 continue ip(nas-nbc+1) = nnz+1 c c Invert the derivative of h (i.e., dh) with the SPARSKIT2 c iterative solver. c call ilu0(nas-nbc,dh,ja,ip,alu,jlu,ju,iwork,ierr) if (ierr.ne.0) then write(*,*) 'ERROR in ilu0, ierr = ',ierr write(9,*) 'ERROR in ilu0, ierr = ',ierr stop endif eps = 1.0e-10 iout = 6 im = nm call pgmres(nas-nbc,im,h,sol,vv,eps,maxit2,iout,dh,ja, & ip,alu,jlu,ju,ierr) if (ierr.eq.1) then write(*,*) 'WARNING: max. iterations in pgmres' write(9,*) 'WARNING: max. iterations in pgmres' endif c c Do Newton update. c sum1 = 0.d0 k = 0 do 380 j=1,ns do 380 i=1,na if (ibc((j-1)*na+i).eq.0) then k = k+1 a1(i,j) = a0(i,j) - sol(k) sum1 = sum1 + sol(k)*sol(k) else a1(i,j) = a0(i,j) endif 380 continue sum1 = dsqrt(sum1)/float(nas-nbc) write(9,'(''The residual at iteration '',I3,'' is '', & F11.9)') it3,sum1 c c Check to see if the solution is converged. If so, stop. c if (sum1.lt.crit3) go to 996 if (it3.eq.maxit3) & write(9,'(''WARNING: iteration limit of '',I2, & '' reached in inner loop.'')') it3 do 390 j=1,ns do 390 i=1,na a0(i,j) = a1(i,j) 390 continue 200 continue 996 continue c c For the current interest rate, we have optimal decision c functions. The next major computational task is to compute c the invariant distribution for assets. c nx = nf-1 c c Construct a = A^{-1}(a',l) from solution just computed, a'=A(a,l). c do 400 j=1,ns do 410 i=1,na-1 if (a1(i,j).gt.a1(i+1,j)) then write(*,*) 'ERROR: functions not increasing' write(*,*) 'ERROR: functions written to unit 55' write(9,*) 'ERROR: functions not increasing' write(9,*) 'ERROR: functions written to unit 55' write(55,*) b write(55,*) do 420 l=1,ns do 420 k=1,na write(55,*) a1(k,l) 420 continue write(55,*) go to 110 endif 410 continue do 430 i=1,nf ii = 1 do 440 k=1,na-1 if (xf(i).ge.xa(k)) ii = k 440 continue x1t = xa(ii) x2t = xa(ii+1) dxt = x2t-x1t xi = 2.d0*(xf(i)-x1t)/dxt-1.d0 bs1t = 0.5d0*(1.d0-xi) bs2t = 0.5d0*(1.d0+xi) atem(i,j) = a1(ii,j)*bs1t+a1(ii+1,j)*bs2t 430 continue do 450 i=1,nf ii = 1 do 460 k=1,nx if (xf(i).ge.atem(k,j)) ii = k 460 continue a1t = atem(ii,j) a2t = atem(ii+1,j) dxt = a2t-a1t xi = 2.d0*(xf(i)-a1t)/dxt-1.d0 bs1t = 0.5d0*(1.d0-xi) bs2t = 0.5d0*(1.d0+xi) ai(i,j)= xf(ii)*bs1t+xf(ii+1)*bs2t 450 continue 400 continue c c The Newton-Raphson iteration. c do 470 it3=1,maxit3 do 480 i=1,nfs h(i) = 0.d0 480 continue do 490 i=1,ldh dh(i) = 0.d0 490 continue nnz = 0 c c Construct h(f0) = int Q(a,l;f0) dx and its derivative, where c c Q(a',j;F) = sum_i pi(i,j) F(A^-1(a',i),i)chi(a'>A(0,i))-F(a',j). c do 500 j=1,ns do 510 n=1,nx x1 = xf(n) x2 = xf(n+1) dx = x2-x1 do 520 i=1,nxp(n) x = absf(i,n) xwgt = wgtf(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 cdf using a finite element approximation. c f = f0(n,j)*bs1+f0(n+1,j)*bs2 c c Compute sum_i pi(i,j) F(A^{-1}(a',i),i)chi(a'>A(0,i)). c rhs = 0.d0 do 530 k=1,nfs work(k) = 0.d0 530 continue do 540 k=1,ns on = 0.d0 if (x.ge.atem(1,k)) on = 1.d0 xt = ai(n,k)*bs1+ai(n+1,k)*bs2 ii = 1 do 550 l=1,nx if (xt.ge.xf(l)) ii = l 550 continue x1t = xf(ii) x2t = xf(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) ft = f0(ii,k)*bs1t+f0(ii+1,k)*bs2t rhs = rhs + pi(k,j)*ft*on trft = pi(k,j)*on m = ii+(k-1)*nf work(m) = work(m) + trft*bs1t work(m+1) = work(m+1) + trft*bs2t 540 continue c c The vector `work' contains derivatives of the residual. c m = n+(j-1)*nf work(m) = work(m) - bs1 work(m+1) = work(m+1) - bs2 c c Add the residual x quadrature weights x shape c functions to h(.), where c c Res = sum_i pi(i,j)F(A^-1(a',i),i)chi(a'>A(0,i))-F(a',j) c res = rhs-f h(m) = h(m) + res*xwgt*bs1 h(m+1) = h(m+1) + res*xwgt*bs2 c c If vector of derivatives (dh) is too long, sort c and eliminate duplicates. c if (nnz.gt.ldh-2*nfs) then call qcksrt3(nnz,ind,dh) k = ind(1) l = 1 do 555 m=2,nnz if (ind(m).eq.k) then dh(l) = dh(l)+dh(m) else k = ind(m) l = l+1 ind(l) = ind(m) dh(l) = dh(m) endif 555 continue nnz = l endif if (nnz.gt.ldh-2*nfs) then write(*,*) 'WARNING: increase length of dh' write(9,*) 'WARNING: increase length of dh' stop endif c c Add derivative of the residual x quadrature c weights x shape functions to dh(.). NOTE: (i,j) c indices switched. c m = n+(j-1)*nf do 560 k=1,nfs if (dabs(work(k)).gt.0.d0) then nnz = nnz+1 ind(nnz) = (m-1)*nfs+k dh(nnz) = work(k)*xwgt*bs1 nnz = nnz+1 ind(nnz) = m*nfs+k dh(nnz) = work(k)*xwgt*bs2 endif 560 continue 520 continue 510 continue 500 continue call qcksrt3(nnz,ind,dh) k = ind(1) l = 1 do 570 m=2,nnz if (ind(m).eq.k) then dh(l) = dh(l)+dh(m) else k = ind(m) l = l+1 ind(l) = ind(m) dh(l) = dh(m) endif 570 continue nnz = l c c Convert (j-1)*nfs+i stored in ind(.) into i and j c do 580 i=1,nnz ja(i) = mod(ind(i)-1,nfs)+1 ia(i) = (ind(i)-1)/nfs+1 580 continue c c Fix value of f0(nf,ns) to probability of being in state ns c by eliminating the last equation in the system of equations. c icnt = 0 do 590 i=nnz,1,-1 if ((ia(i).eq.nfs).or.(ja(i).eq.nfs)) then do 600 j=i+1,nnz-icnt ia(j-1) = ia(j) ja(j-1) = ja(j) dh(j-1) = dh(j) 600 continue icnt = icnt+1 endif 590 continue nnz = nnz-icnt c c Use ia(.) to determine pointers to rows in sparse matrix. c ip(1) = 1 j = 1 do 610 i=2,nnz if (ia(i).ne.ia(i-1)) then j = j+1 ip(j) = i endif 610 continue ip(nfs) = nnz+1 c c Invert the derivative of h (i.e., dh) with the SPARSKIT2 c iterative solver. c call ilu0(nfs-1,dh,ja,ip,alu,jlu,ju,iwork,ierr) if (ierr.ne.0) then write(*,*) 'ERROR in ilu0, ierr = ',ierr write(9,*) 'ERROR in ilu0, ierr = ',ierr stop endif eps = 1.0e-10 iout = 6 im = nm call pgmres(nfs-1,im,h,sol,vv,eps,maxit2,iout,dh,ja, & ip,alu,jlu,ju,ierr) if (ierr.eq.1) then write(*,*) 'WARNING: max. iterations in pgmres' write(9,*) 'WARNING: max. iterations in pgmres' endif c c Do Newton update. c sum1 = 0.d0 k = 0 do 620 j=1,ns-1 do 620 i=1,nf k = k + 1 f1(i,j) = f0(i,j) - sol(k) sum1 = sum1 + sol(k)*sol(k) 620 continue do 630 i=1,nf-1 k = k + 1 f1(i,ns) = f0(i,ns) - sol(k) sum1 = sum1 + sol(k)*sol(k) 630 continue sum1 = dsqrt(sum1)/float(nfs) f1(nf,ns) = f0(nf,ns) c c If distance between f0 and f1 sufficiently small, stop. c if (sum1.lt.crit3) go to 997 do 640 j=1,ns do 640 i=1,nf f0(i,j) = f1(i,j) 640 continue 470 continue 997 continue c c Set any negative elements of f1 equal to 0 and force c monotonicity. c do 650 j=1,ns do 650 i=1,nf f1(i,j) = dmax1(0.d0,f1(i,j)) 650 continue do 660 j=1,ns do 670 i=3,nf-1 if (f1(i,j).lt.f1(i-1,j)) then tem = (f1(i-1,j)-f1(i-2,j))/(xf(i-1)-xf(i-2)) if (tem.gt.2.0) then do 680 k=i+1,nf if (f1(k,j).ge.f1(k-1,j)) then do 690 l=i-1,k-1 f1(l,j) = f1(k-1,j) 690 continue go to 670 endif 680 continue else f1(i,j) = f1(i-1,j) endif endif 670 continue 660 continue c c Compute mean asset holdings (sum_i int x dF(x,i)). c amean = 0.d0 do 700 j=1,ns do 700 i=1,nx amean = amean+0.5*(f1(i+1,j)-f1(i,j))*(xf(i+1)+xf(i)) 700 continue if ((aggk+b).gt.amean) then bndl = r else bndu = r endif sum1 = dabs(r-0.5*(bndu+bndl)) tem = 0.00001 if ((dabs(r-bndu0).lt.tem).or.(dabs(r-bndl0).lt.tem)) & then write(*,*) 'WARNING: after-tax int. rate near bound' write(9,*) 'WARNING: after-tax int. rate near bound' write(*,*) r,bndl0,bndu0 go to 110 endif write(9,'(''The interest rate and residual at iteration'', & I4,'' is '',F11.9,'' and '',F11.9)') it2,r,sum1 if (sum1.lt.crit2) go to 999 190 continue 999 continue c c Update estimate of aggregate labor supply. c hrs = 0.d0 do 710 j=1,ns do 710 i=1,nx lei = lei0 leit = lei0 do 720 m=1,15 c = r1*xf(i)+w*edu(j)*(1.d0-lei)-tax-grat1* & atem(i,j) ct = r1*xf(i+1)+w*edu(j)*(1.d0-leit)-tax-grat1* & atem(i+1,j) u = (c**eta*lei**eta1)**mu1 ut = (ct**eta*leit**eta1)**mu1 d = dmin1(1.d0-lei,0.d0) dt = dmin1(1.d0-leit,0.d0) f = eta1*u/lei-eta*u*w*edu(j)/c-zeta*d*d ft = eta1*ut/leit-eta*ut*w*edu(j)/ct-zeta*dt*dt fj = -2.d0*eta1*eta*mu1*u/lei/c*w*edu(j)+eta1*(eta1* & mu1-1.d0)*u/lei/lei+eta*(eta*mu1-1.d0)*u*w*w* & edu(j)*edu(j)/c/c+2.d0*zeta*d fjt = -2.d0*eta1*eta*mu1*ut/leit/ct*w*edu(j)+eta1* & (eta1*mu1-1.d0)*ut/leit/leit+eta*(eta*mu1- & 1.d0)*ut*w*w*edu(j)*edu(j)/ct/ct+2.d0*zeta*dt lei = lei-f/fj leit = leit-ft/fjt 720 continue hrs = hrs+0.5*(f1(i+1,j)-f1(i,j))*edu(j)*(2.d0- & lei-leit) if (i.eq.1) hrs = hrs + f1(1,j)*edu(j)*(1.d0-lei) 710 continue if (it0.eq.0) then fl = aggl-hrs dfl = fl else dfl = (aggl-hrs-dfl)/del aggl = aggli-fl/dfl if (dabs(aggli-aggl).lt.crit1) go to 998 aggli= aggl endif 180 continue write(9,'(''The labor supply at iteration'', & I4,'' is '',F11.9)') it1,aggli 170 continue 998 continue c c We now have the equilibrium interest rate, optimal decision c functions given this interest rate, and the invariant distribution c for the asset holdings. The next major computational task is to c compute the value function, V(x,i), which satisfies Bellman's c equation, and social welfare, sum_i int V(x,i) dF(x,i). The c short cut to this calculation is to compute: c c sum_i int U_i(a) dF(a,i) / (1-beta). c value = 0.d0 hrs = 0.d0 con = 0.d0 alei = 0.d0 do 730 j=1,ns do 730 i=1,nx lei = lei0 leit = lei0 do 740 m=1,15 c = r1*xf(i)+w*edu(j)*(1.d0-lei)-tax-grat1*atem(i,j) ct = r1*xf(i+1)+w*edu(j)*(1.d0-leit)-tax-grat1* & atem(i+1,j) u = (c**eta*lei**eta1)**mu1 ut = (ct**eta*leit**eta1)**mu1 d = dmin1(1.d0-lei,0.d0) dt = dmin1(1.d0-leit,0.d0) f = eta1*u/lei-eta*u*w*edu(j)/c-zeta*d*d ft = eta1*ut/leit-eta*ut*w*edu(j)/ct-zeta*dt*dt fj = -2.d0*eta1*eta*mu1*u/lei/c*w*edu(j)+eta1*(eta1*mu1- & 1.d0)*u/lei/lei+eta*(eta*mu1-1.d0)*u*w*w*edu(j)* & edu(j)/c/c+2.d0*zeta*d fjt = -2.d0*eta1*eta*mu1*ut/leit/ct*w*edu(j)+eta1*(eta1* & mu1-1.d0)*ut/leit/leit+eta*(eta*mu1-1.d0)*ut*w*w* & edu(j)*edu(j)/ct/ct+2.d0*zeta*dt lei = lei-f/fj leit = leit-ft/fjt 740 continue value = value+0.5*(f1(i+1,j)-f1(i,j))*(u+ut)/bet1/mu1 hrs = hrs+0.5*(f1(i+1,j)-f1(i,j))*edu(j)*(2.d0-lei-leit) con = con+0.5*(f1(i+1,j)-f1(i,j))*(c+ct) alei = alei+0.5*(f1(i+1,j)-f1(i,j))*(lei+leit) if (i.eq.1) then value = value + f1(1,j)*u/bet1/mu1 hrs = hrs + f1(1,j)*edu(j)*(1.d0-lei) con = con + f1(1,j)*c alei = alei + f1(1,j)*lei endif 730 continue chk(it,1) = b chk(it,2) = con chk(it,3) = y-g-(grate+delta)*aggk chk(it,4) = hrs chk(it,5) = aggl chk(it,6) = alei if (it.eq.1) then r0 = r aggl0 = aggl val0 = value write(7,7) b,0.0,amean,r,taul,value,aggl,y write(*,7) b,0.0,amean,r,taul,value,aggl,y write(9,7) b,0.0,amean,r,taul,value,aggl,y else c c Compute welfare gain in terms of consumption (x), e.g., c c (1+x)^(eta(1-mu)) = (Y/Y0)^(eta(1-mu)) * (V/V0), c c where Y=output, V=value, and 0 indicates values for the baseline case. c gain = ((r+delta)/(r0+delta))**(alpha/alph1)*aggl/aggl0* & (value/val0)**(1.d0/(eta*mu1))-1.d0 c c Print results of each case. c write(7,7) b,100*gain,amean,r,taul,value,aggl,y write(*,7) b,100*gain,amean,r,taul,value,aggl,y write(9,7) b,100*gain,amean,r,taul,value,aggl,y endif write(8,*) write(8,*) write(8,2) b,r-0.002,r+0.002,nbc write(8,*) do 750 i=1,na write(8,5) xa(i),(a1(i,j),j=1,ns) 750 continue write(8,*) do 760 i=1,na ii = 0 do 770 j=1,ns if (ibc(i+(j-1)*na).eq.1) then ii = ii+2 iwork(ii-1) = i iwork(ii) = j endif 770 continue if (ii.gt.0) write(8,6) (iwork(j),j=1,ii) 760 continue 110 continue c c Check to see that the following are approximately true: c c E[c] = 1-g-alpha*(g+delta)/(r+delta) c E[e*l] = 1-L c write(7,*) write(7,*) ' Debt Mean of Consumption', & ' Mean of Hours Mean of' write(7,*) ' Using F wo/ F', & ' Using F wo/ F Leisure' do 790 i=1,ncases+1 write(7,4) (chk(i,j),j=1,6) 790 continue c c Format statements. c 1 format (7(1X, F10.6)) 2 format (F8.4,2(1X,F10.6),1X,I3) 3 format (4X, 7(1X, F7.4)) 4 format (7(1X, F10.5)) 5 format (F9.6,7(1X,F9.5)) 6 format (8(1X,I2)) 7 format (5(1X, F10.5),1X,F12.4,1X,F9.5,1X,F9.5) 8 format (a43,a43) stop end