c PROGRAM: JUDD.F c c DESCRIPTION: implementation of a weighted residual method with c Chebyshev polynomials for the 1-D deterministic growth c model c ECONOMIC the discrete-time stochastic growth model with c MODEL: agents choosing c(t) to maximize c c E[ sum beta^t c[t]^(1+gamma)/(1+gamma) | k_0] c t c subj. to c k[t+1] = lambda*k[t]^alpha+(1-delta)*k[t]-c[t] c c[t],k[t] >= 0 c i[t] >= 0 for sufficiently large value of chi c c PROBLEM: find c=c(x) that solves the functional equation: c c c^gamma = beta pi(i,j) c~^gamma*(alpha*lambda* c x~^(alpha-1) +1-delta) c where c c~ = c(x~), c x~ = lambda*x^alpha+(1-delta)*x-c(x) c c(0) = 0 c x in [0,xmax] c c ALGORITHM: a weighted residual method c c TEST CASE: delta=1, gamma=-1 ==> c = (1-beta*alpha)*lambda*k^alpha c c INPUT FILE: (a) judd.inp: the following need to be specified c c Input Description Format c ------------------------------------------------------ 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 gamma utility parameter real <= 0 c lambda production parameter real > 0 c kmin lower bound on k real >= 0 c kmax upper bound on k real > kmin c initialize option to provide 1 (=yes) or 0 c initial guess for a c a initial guess for n x 1 real vector c polynomial weights vector c (if initialize=1) c c OUTPUT FILES: (b) judd.nxt: new input file with solution for c(k,i) c (c) judd.out: output file with solution for a,c(k) c c REQUIRED (BLAS) dgemm,dtrsm,idamax,dger,dscal,dswap c PROGRAMS: (LAPACK) dgesv c (PRESS) qgausl c c REFERENCES: Judd (1992), JET c Lapack and BLAS, ftp to netlib2.cs.utk.edu c Press, et. al. (1986), NUMERICAL RECIPES c Taylor and Uhlig (1990), JBES c c SEE ALSO: sgm1dl.f c c Ellen McGrattan, 2-20-97 c Revised, 2-21-97, ERM c * c c Set dimensions: c n = number of polynomials in representation for consumption c m = number of quadrature points (m>n) c parameter (n=10,m=20) implicit double precision (a-h,o-z) double precision a0(n), a1(n), c(m), ct(m), x(m), xt(m), & phi(n,m), psi(n,m), f(n), df(n,n), res(m), & dres(m,n), lambda, k(m), kt(m), kmin, kmax integer iwork(n) c c Open input and output files. c open(unit=5, file='judd.inp') open(unit=7, file='judd.nxt') open(unit=8, file='judd.out') c c Read in parameters, grid, boundary conditions, and if specified, c the initial guess for consumption. c read(5,*) alpha read(5,*) beta read(5,*) delta read(5,*) gamma read(5,*) kmin read(5,*) kmax read(5,*) initialize lambda = (1.d0-beta*(1.d0-delta))/alpha/beta alph1 = alpha-1.d0 gamm1 = gamma-1.d0 delt1 = 1.d0-delta pi = 3.14159265358979d0 fm = float(m) fn = float(n) if (initialize.eq.1) then do 10 i=1,n read(5,*) a0(i) 10 continue else c c If initial guess for consumption not given, set a0(1) = c a0(2) = .5*(lambda-delta), a0(j) = 0, j>2. c a0(1) = .5d0*(lambda-delta) a0(2) = a0(1) do 20 i=3,n a0(i) = 0.d0 20 continue endif c c Calculate grid on capital stocks using roots of Chebyshev polynomial. c do 30 i=1,m fi = float(i) x(i) = dcos(pi*(fm-fi+0.5d0)/fm) k(i) = 0.5d0*(kmax-kmin)*(1.d0+x(i))+kmin 30 continue do 40 j=1,m do 40 i=1,n phi(i,j) = cos(float(i-1)*acos(x(j))) 40 continue c c Other parameters that are set inside the program c * stopping criteria c * maximum number of interations c crit = 0.0000001d0 maxit = 100 c c The Newton-Raphson iteration, i.e., find a such that g(a)=0 where c a are coefficients of the weighted residual approximation. c do 50 it=1,maxit c c Initialize f and its derivative df. c do 60 i=1,n f(i) = 0.d0 do 60 j=1,n df(i,j) = 0.d0 60 continue c c Construct f(a) = int R(x;a) phi_i(x) dx and its derivatives with c respect to a. c do 70 i=1,m c(i) = 0.d0 do 80 j=1,n c(i) = c(i) + phi(j,i)*a0(j) 80 continue kt(i) = lambda*k(i)**alpha+delt1*k(i)-c(i) xt(i) = 2.d0*(kt(i)-kmin)/(kmax-kmin)-1.d0 70 continue do 90 j=1,m do 90 i=1,n psi(i,j) = dcos(float(i-1)*dacos(xt(j))) 90 continue do 100 i=1,m sum1 = 0.d0 do 110 j=1,n sum1 = sum1 + psi(j,i)*a0(j) 110 continue ct(i) = sum1 res(i) = c(i)**gamma-beta*ct(i)**gamma*(alpha*lambda* & kt(i)**alph1+delt1) 100 continue do 120 j=1,m sum1 = 0.d0 do 130 i=1,n fim1 = float(i-1) sum1 = sum1 +fim1*dsin(fim1*dacos(xt(j)))/dsqrt(1.d0- & xt(j)*xt(j))*a0(i) 130 continue do 140 i=1,n dres(j,i) = gamma*c(j)**gamm1*phi(i,j)-beta*gamma* & ct(j)**gamm1*(alpha*lambda*kt(j)**alph1+ & delt1)*(psi(i,j)-2.d0*sum1*phi(i,j)/(kmax- & kmin))+beta*alpha*alph1*lambda*ct(i)**gamma* & kt(i)**(alph1-1.d0)*2.d0*phi(i,j)/(kmax-kmin) 140 continue 120 continue do 150 i=1,n f(i) = 0.d0 do 150 j=1,m f(i) = f(i) + phi(i,j)*res(j) 150 continue do 160 i=1,n do 160 j=1,n df(i,j) = 0.d0 do 170 l=1,m df(i,j) = df(i,j) + phi(i,l)*dres(l,j) 170 continue 160 continue c c Solve df\f and put solution in f. c call dgesv(n,1,df,n,iwork,f,n,info) c c Issue warnings if necessary. c if (info.lt.0) write(*,*) 'WARNING: illegal value for', c ' argument ',-info,' in JUDD' if (info.gt.0) write(*,*) 'WARNING: singular matrix', c ' when solving df*x=f in JUDD' c c Do Newton update. c sum1 = 0.d0 do 180 i=1,n a1(i) = a0(i) - f(i) sum1 = sum1 + f(i)*f(i) 180 continue sum1 = dsqrt(sum1)/fn c c Check to see if the solution is converged. If so, stop. c if (sum1.lt.crit) go to 999 do 190 i=1,n a0(i) = a1(i) 190 continue 50 continue 999 continue c c Write new input file with results. c write(7,'(9X,F12.6,12X,''alpha in (0,1) '')') alpha write(7,'(9X,F12.6,12X,''beta in (0,1) '')') beta write(7,'(9X,F12.6,12X,''delta in [0,1] '')') delta write(7,'(9X,F12.6,12X,''gamma <= 0 '')') gamma write(7,'(3X,F18.6,12X,''kmin >= 0 '')') kmin write(7,'(3X,F18.6,12X,''kmax > kmin '')') kmax write(7,'(20X,''1'',12X,''read initial guess for a (1=yes)'')') write(7,*) write(7,'(9X,F12.6,12X,''polynomial weights'')') a1(1) do 200 i=2,n write(7,1) a1(i) 200 continue c c Write solution to output file. c write(8,*)'RESULTS FOR JUDD' write(8,*) write(8,*) ' Parameters: ' write(8,'('' alpha = '',F7.4,'' '')') alpha write(8,'('' beta = '',F7.4,'' '')') beta write(8,'('' delta = '',F7.4,'' '')') delta write(8,'('' gamma = '',F7.4,'' '')') gamma write(8,'('' kmin = '',F7.4,'' '')') kmin write(8,'('' kmax = '',F7.4,'' '')') kmax write(8,'('' lambda = '',F7.4,'' '')') lambda write(8,*) write(8,*) ' Chebyshev Coefficients: ' do 210 i=1,min(n,9) write(8,'('' a'',I1,'' = '',F7.4,'' '')') i,a1(i) 210 continue do 220 i=10,n write(8,'('' a'',I2,'' = '',F7.4,'' '')') i,a1(i) 220 continue write(8,*) write(8,*) ' Capital and consumption: ' do 230 i=1,m write(8,2) k(i),c(i) 230 continue c c Format statements. c 1 format (9X,F12.6) 2 format (4X,F7.4,3X,F7.4) stop end