c PROGRAM: EXOG.F c c DESCRIPTION: Computation of equilibrium in one-sector exogenous 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],h[t]) c c subject to c c c[t] + (1+taux[t])*x[t] = w[t]*h[t]+r[t]*k[t] +T[t] c T[t] = taux[t]*x[t] c k[t+1] = (1-delta)*k[t] + x[t] c x[t] >= 0 (penalty functions used to ensure this) c c where processes for r,w,taux assumed known c c Firms: max_k F(k,1) - w*h - r*k c c PROBLEM: find c=c(k,i) that solves the functional equation: c c U_1(c)*(1+taux(i)) = beta pi(i,j) U_1(c~) * r~ c c where c c~ = c(k~,j), c r~ = F_1(k~,h~)+(1-delta)*(1+taux(j)) c k~ = F(k,h)-c+(1-delta)*k c pi(i,j) = probability of transiting from state i to j c (see example below for derivation of pi) c k in [0,kmax] c i,j in {1,2,...,ns} c c MARKOV Let theta[t] = 1+taux[t]. The estimated process for c PROCESS: theta[t] is given by c c theta[t+1] = mu_R*(1-rho_R)+rho_R*theta[t]+sig_R*eps[t+1] c c where R is the regime that a country is in and the c probability of switching regimes depends on how c long a country has been in the regime. R = `persistent' c or `turbulent' and the probability of switching is: c c Prob(Switch from R to R' after being in R tau periods) c = a_R + b_R * (tau-1), tau = 1:20 c = a_R + b_R * 19, tau >= 20 c c The Markov chain governing the tax rate is a discrete c state description of the above continuous process with c each element being associated with a different triplet c (theta,R,tau). See createp.m for the code that maps c mu_R,rho_R,sig_R,a_R,b_R to a discrete set of tax rates c (taux) and a transition matrix (pi) for a Markov chain c governing the evolution of the tax rates. c c ALGORITHM: the finite element method with piecewise bilinear c shape functions N_a(x), ie., c(x,i) = sum N_a(x)*c^i_a, c where c^i_a's are the coefficients computed here. c c FUNCTIONS: U(c,h) = [(c*(1-h)^psi)^(1-sigma)-1]/(1-sigma) c F(k,h) = a*k^alpha*h^(1-alpha) c c TEST CASE: psi=0,delta=1,sigma=1 ==> c = (1-beta*alpha)*a*k^alpha. c c INPUT FILE: (a) exog.inp: the following need to be specified c c Input Description Format c ------------------------------------------------------ c a scale factor real > 0 c alpha capital share real in (0,1) c beta discount factor real in (0,1) c delta depreciation rate real in [0,1] c glab growth in labor force real (e.g., .01) c gtec growth in technology real (e.g., .02) c kbar scale factor for real > 0 c capital so k in [0,1] c psi utility parameter real >= 0 c sigma utility parameter real >= 0 c pen penalty parameter 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 guess for h c nbc number of boundary integer c conditions 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 npl number per line tax rates are read c in groups of npl, e.g. c if rates are listed: c .0 .1 .2 .3 c .4 .5 .6 .7 c .8 c in input file, then c npl = 4 and ns = 9 c taux investment tax grid ns x 1 real vector c (read in groups of c npl per line) c nnzp number of nonzero integer c elements of pi c i,j,pi nonzero elements of nnzp triplets of c the transition (integer,integer, c matrix, pi(i,j) real); read nonzero c elements of row 1 c first, row 2 second, c etc. 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 i,j,c points c(x(i),j) triplets of (integer, c that are to be fixed integer,real) c during computation c c Notes: c i. see pov_tra.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 (b) relinc.dat: relative incomes from Summers-Heston data. c c OUTPUT: (c) exog.nxt: new input file with solution for c(k,i) c (d) exog.dat: input file for Matlab post-processing c routine -- type exog in Matlab. c c c REQUIRED (SPARSKIT2) ilut, matvec c PROGRAMS: (PRESS) qgausl, qcksrtr, qcksrtr2, qcksrt3 c (KAHANER) uni c c REFERENCES: Lapack and BLAS, ftp to netlib2.cs.utk.edu c 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: derivs[1-8].f c c Ellen McGrattan, 4-11-98 c Revised, ERM, 11-19-98 c * c c Set dimensions: c nx = length of x partition c ns = number of discrete states for tax rate c na = nx*ns c nq = maximum number of quadrature points per element; c nq>=nxp(i), i=1:ne c nm = number of columns stored for sparse matrix inversion (max=50) c ldf = dimension of Jacobian vector (df) c lpi = dimension of vector for transition probabilities (pi) c parameter (nx=31,ns=800,na=24800,nq=5,nm=40,ldf=1600000,lpi=10000) implicit double precision (a-h,o-z) double precision f(na), df(ldf), work(na), xa(nx), ax(nq,nx-1), & wx(nq,nx-1), c0(nx,ns), c1(nx,ns), taux(ns), pi(lpi), & kbar, lhs, lei, 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='exog.inp') open(unit=7, file='exog.nxt') open(unit=8, file='exog.dat') c c Read in parameters. c read(5,*) a read(5,*) alpha read(5,*) beta read(5,*) delta read(5,*) glab read(5,*) gtec read(5,*) kbar read(5,*) psi read(5,*) sigma read(5,*) pen read(5,*) ievenx read(5,*) initialize read(5,*) nbc c alph1 = 1.d0-alpha delt1 = 1.d0-delta grate = (1.d0+glab)*(1.d0+gtec) sigm1 = 1.d0-sigma betg = beta*(1.d0+gtec)**(-sigma) c c Read in grid for capital stock (normalized to be on [0,1]). 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 investment distortions. c read(5,*) npl do 30 i=1,ns/npl read(5,*) (taux(j),j=(i-1)*npl+1,i*npl) 30 continue read(5,*) (taux(j),j=ns-mod(ns,npl)+1,ns) 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 function. c if (initialize.eq.1) then do 60 j=1,ns do 60 i=1,nx read(5,*) c0(i,j) ibc(i+(j-1)*nx) = 0 60 continue else c c If initial guess is not given, set c=(1-alpha*beta)*output. c do 70 j=1,ns do 70 i=1,nx k = i+(j-1)*nx c0(i,j) = (1.d0-alpha*betg)*a*(kbar*xa(i))**alpha/kbar ibc(k) = 0 70 continue endif c c Impose boundary conditions. c do 80 k=1,nbc read(5,*) i,j,tem ibc(i+(j-1)*nx) = 1 c0(i,j) = tem 80 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 iterations. c * gamma, the kink point used for the piecewise approximation c of h^alpha, ie. h^alpha ~= a*x, xgamma c do 100 i=1,nxm1 nxp(i) = 3 100 continue crit = 0.0000001d0 maxit = 200 gamma = 0.1d0 gama = gamma**(-alph1) gam1 = (1.d0-gamma**alpha)/(1.d0-gamma) 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) = int R(x;c0) N_n(x) dx and its derivative on c 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 using a finite element approximation. c c = c0(n,j)*bs1+c0(n+1,j)*bs2 c c Back out labor/leisure from static first order conditions c and definitions. The fixed point problem for hours c is of the form: h^alpha = constant * (1-h). Use a c piecewise linear function (i.e., y = a*x if x<=gamma and c y=b*x+(a-b)*gamma if x>gamma) to approximate h^alpha. c With this piecewise approximation, it is easy to derive c an initial guess for the fixed point of the nonlinear equation. c tem1 = alph1*a*(x*kbar)**alpha / kbar tem2 = tem1/psi/c h = tem2/(tem2+gama) if (h.gt.gamma) then h = (tem2-(gama-gam1)*gamma)/(tem2+gam1) endif lei = 1.d0-h do 165 k=1,20 tem = psi*c - tem1*h**(-alpha)*lei dtem = tem1*h**(-alpha)*(alpha*lei/h+1.d0) h = h - tem/dtem lei = 1.d0-h 165 continue c c Back out other variables from first-order conditions c and definitions. c y = a*(x*kbar)**alpha*h**alph1 / kbar u = (c*lei**psi)**sigm1/sigm1 duc = sigm1*u/c d = dmin1(y-c,0.d0) s = pen*d*d taux1 = 1.d0+taux(j) c c Compute left-hand side of dynamic first-order condition. c lhs = duc*taux1-s c c Compute total derivative of lhs with respect to unknowns (tlc). c d2uc = -sigma*duc/c duch = -sigm1*psi*duc/lei dlh = duch*taux1-2.d0*pen*d*alph1*y/h dlc = d2uc*taux1+2.d0*pen*d thc = -psi/(tem1*h**(-alpha)*(alpha*lei/h+1.d0)) tlc = dlc + dlh*thc c c Find next period's capital stock on grid. c xt = (y-c+delt1*x)/grate txtc = (alph1*y/h*thc-1.d0)/grate 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) c c Initialize variables for summation. c rhs = 0.d0 trc = 0.d0 do 180 k=1,na work(k) = 0.d0 180 continue c c The following loop computes the sum in the Euler equation: c c rhs = beta sum pi(j,k) u_1(c~) (F_1(k~,1)+(1-delta)*(1+taux)). c do 190 l=ipi(j),ipi(j+1)-1 k = jpi(l) c c Given next period's state (xt,k), compute next-period c consumption (ct) using a finite element approximation. c ct = c0(ii,k)*bs1t+c0(ii+1,k)*bs2t c c Back out labor/leisure from static first order conditions c and definitions. See note above about initial guess. c tem1 = alph1*a*(xt*kbar)**alpha / kbar tem2 = tem1/psi/ct ht = tem2/(tem2+gama) if (ht.gt.gamma) then ht = (tem2-(gama-gam1)*gamma)/(tem2+gam1) endif lei = 1.d0-ht do 195 m=1,20 tem = psi*ct - tem1*ht**(-alpha)*lei dtem = tem1*ht**(-alpha)*(alpha*lei/ht+1.d0) ht = ht - tem/dtem lei = 1.d0-ht 195 continue c c Back out other variables from first-order conditions c and definitions. c yt = a*(xt*kbar)**alpha*ht**alph1 / kbar ut = (ct*lei**psi)**sigm1/sigm1 dutct = sigm1*ut/ct d = dmin1(yt-ct,0.d0) st = pen*d*d taux1 = 1.d0+taux(k) rt = alpha*yt/xt + delt1*taux1 c c Compute right-hand side sum of dynamic first-order condition. c rhs = rhs + betg*pi(l)*(dutct*rt-delt1*st) c c Compute total derivatives of `rhs' in terms of unknown c coefficients in representation of c (i.e., trc) and c unknown coefficients in representation of ct (i.e., trct). c d2utct = -sigma*dutct/ct dutcht = -sigm1*psi*dutct/lei tctxt = -(c0(ii,k) -c0(ii+1,k))/dxt tem1 = a*alph1*(xt*kbar/ht)**alpha / kbar dhtxt = alpha*lei/xt/(alpha*lei/ht+1.d0) thtct = -psi/(tem1*(alpha*lei/ht+1.d0)) thtxt = dhtxt+thtct*tctxt drht = dutcht*rt+dutct*alpha*alph1*yt/xt/ht- & 2.d0*delt1*pen*d*alph1*yt/ht drct = d2utct*rt+2.d0*delt1*pen*d drxt =-dutct*alpha*alph1*yt/xt/xt trxt = drxt+drht*thtxt+drct*tctxt trc = trc +betg*pi(l)*trxt*txtc trct = betg*pi(l)*(drct+drht*thtct) c c Store derivatives in work vector `work'. c m = ii+(k-1)*nx work(m) = work(m) + trct*bs1t work(m+1) = work(m+1) + trct*bs2t 190 continue c c The vector `work' contains derivatives of the residual. c dres = trc-tlc m = n+(j-1)*nx work(m) = work(m) + dres*bs1 work(m+1) = work(m+1) + dres*bs2 c c Add the residual x quadrature weights x shape functions to f(.). c res = rhs-lhs f(m) = f(m) + res*w*bs1 f(m+1) = f(m+1) + res*w*bs2 c c If vector of derivatives (df) is too long, sort and eliminate c duplicates. c if (nnz.gt.ldf-2*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-2*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 = n+(j-1)*nx do 210 k=1,na if (dabs(work(k)).gt.0.d0) then nnz = nnz+1 ind(nnz) = (m-1)*na+k df(nnz) = work(k)*w*bs1 nnz = nnz+1 ind(nnz) = m*na+k df(nnz) = work(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(.) 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((j-1)*nx+i).eq.0) then k = k+1 c1(i,j) = c0(i,j) - sol(k) sum1 = sum1 + sol(k)*sol(k) else c1(i,j) = c0(i,j) 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) 280 continue 120 continue 999 continue c c Write new input file with results. c write(7,'(6X,F15.9,12X,''a > 0 '')') a write(7,'(6X,F15.9,12X,''alpha in (0,1) '')') alpha write(7,'(6X,F15.9,12X,''beta in (0,1) '')') beta write(7,'(6X,F15.9,12X,''delta in [0,1] '')') delta write(7,'(6X,F15.9,12X,''growth in labor '')') glab write(7,'(6X,F15.9,12X,''growth in technology'')') gtec write(7,'(6X,F15.9,12X,''kbar > 0 '')') kbar write(7,'(6X,F15.9,12X,''psi >= 0 '')') psi write(7,'(6X,F15.9,12X,''sigma >= 0 '')') sigma write(7,'(1X,F20.2,12X,''pen >= 0 '')') pen write(7,'(6X,I15,12X,''evenly-spaced grid (1=yes,0=no)'')') ievenx write(7,'(20X,I1,12X,''read initial guess for c (1=yes,0=no)'')')1 write(7,'(6X,I15,12X,''number of boundary constraints on c'')')nbc write(7,*) write(7,'(6X,F15.9,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,'(6X,F15.9,12X,''upper bound for capital grid'')') xa(nx) write(7,*) write(7,'(9X,I3,21X,''tax rates follow (# per line)'')') npl write(7,*) do 300 i=1,ns/npl write(7,3) (taux(j),j=(i-1)*npl+1,i*npl) 300 continue write(7,3) (taux(j),j=ns-mod(ns,npl)+1,ns) write(7,*) write(7,'(8X,I6,19X, & ''pi follows (# of nonzero elements)'')') 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,'(6X,F15.9,12X,''consumption as a function of k,i'')') & c1(1,1) do 320 i=2,nx write(7,1) c1(i,1) 320 continue do 330 j=2,ns do 330 i=1,nx write(7,1) c1(i,j) 330 continue write(7,*) k = 0 do 340 j=1,ns do 340 i=1,nx if (ibc(i+(j-1)*nx).eq.1) then if (k.eq.0) then write(7,'(I2,1X,I2,1X,F15.9,12X, & ''boundary conditions '')') i,j,c1(i,j) k = 1 else write(7,2) i,j,c1(i,j) endif endif 340 continue c c Write parameters and solution to data file for Matlab routine. c write(8,1) a write(8,1) alpha write(8,1) beta write(8,1) delta write(8,1) gamma write(8,1) glab write(8,1) gtec write(8,1) kbar write(8,1) psi write(8,1) sigma 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) taux(i) 360 continue write(8,*) nnzp i = 1 j = 2 do 400 k=1,nnzp if (k .eq. ipi(j)) then i = i + 1 j = j + 1 endif write(8,*) i 400 continue do 410 i=1,nnzp write(8,*) jpi(i) 410 continue do 420 i=1,nnzp write(8,1) pi(i) 420 continue do 430 j=1,ns do 430 i=1,nx write(8,1) c1(i,j) 430 continue c c Format statements. c 1 format (6X,F15.9) 2 format (7X,I4,1X,I4,3X,F15.9) 3 format (9X,20(1X,F11.8)) stop end