program judd90 ! ! DESCRIPTION: implementation of a weighted residual method with ! Chebyshev polynomials for the 1-D deterministic growth ! model (for FORTRAN 90) ! ECONOMIC the discrete-time stochastic growth model with ! MODEL: agents choosing c(t) to maximize ! ! E[ sum beta^t c[t]^(1+gamma)/(1+gamma) | k_0] ! t ! subj. to ! k[t+1] = lambda*k[t]^alpha+(1-delta)*k[t]-c[t] ! c[t],k[t] >= 0 ! i[t] >= 0 for sufficiently large value of chi ! ! PROBLEM: find c=c(x) that solves the functional equation: ! ! c^gamma = beta pi(i,j) c~^gamma*(alpha*lambda* ! x~^(alpha-1) +1-delta) ! where ! c~ = c(x~), ! x~ = lambda*x^alpha+(1-delta)*x-c(x) ! c(0) = 0 ! x in [0,xmax] ! ! ALGORITHM: a weighted residual method ! ! TEST CASE: delta=1, gamma=-1 ==> c = (1-beta*alpha)*lambda*k^alpha ! ! INPUT FILE: (a) judd.inp: the following need to be specified ! ! Input Description Format ! ------------------------------------------------------ ! alpha capital share real in (0,1) ! beta discount factor real in (0,1) ! delta depreciation rate real in [0,1] ! gamma utility parameter real <= 0 ! lambda production parameter real > 0 ! kmin lower bound on k real >= 0 ! kmax upper bound on k real > kmin ! initialize option to provide 1 (=yes) or 0 ! initial guess for a ! a initial guess for n x 1 real vector ! polynomial weights vector ! (if initialize=1) ! ! OUTPUT FILES: (b) judd.nxt: new input file with solution for c(k,i) ! (c) judd.out: output file with solution for a,c(k) ! ! REQUIRED (BLAS) dgemm,dtrsm,idamax,dger,dscal,dswap ! PROGRAMS: (LAPACK) dgesv ! (PRESS) qgausl ! ! REFERENCES: Judd (1992), JET ! Lapack and BLAS, ftp to netlib2.cs.utk.edu ! Press, et. al. (1986), NUMERICAL RECIPES ! Taylor and Uhlig (1990), JBES ! ! SEE ALSO: judd.f, sgm1dl.f ! ! Ellen McGrattan, 2-20-97 ! Revised, 2-22-97, ERM ! ! ! Set dimensions: ! n = number of polynomials in representation for consumption ! m = number of quadrature points (m>n) ! ! implicit none integer, parameter :: n = 10, m = 20, maxit = 100 real, parameter :: pi = 3.14159265358979 real, parameter :: crit = 0.0000001 real, dimension(n,1) :: a0 = 0.0, a1, f, sn real, dimension(m,1) :: c, ct, x, xt, k, kt, res, sm real, dimension(1,n) :: en = 1.0 real, dimension(1,m) :: em = 1.0 real, dimension(m,n) :: dres, dcta real, dimension(n,m) :: phi, psi, dpsi real, dimension(n,n) :: df real :: alpha, alph1, beta, delta, delt1 real :: gamma, gamm1, kmin, kmax, lambda real :: fi, fim1, fm, fn, sum1 integer, dimension(n) :: iwork integer :: i, j, l, it, info, initialize ! ! Open input and output files. ! open(unit=5, file='judd.inp') open(unit=7, file='judd.nxt') open(unit=8, file='judd.out') ! ! Read in parameters, grid, boundary conditions, and if specified, ! the initial guess for consumption. ! read(5,*) alpha read(5,*) beta read(5,*) delta read(5,*) gamma read(5,*) kmin read(5,*) kmax read(5,*) initialize lambda = (1.0-beta*(1.0-delta))/alpha/beta alph1 = alpha-1.0 gamm1 = gamma-1.0 delt1 = 1.0-delta fm = float(m) fn = float(n) sn(:,1) = (/ (i,i=0,n-1) /) sm(:,1) = (/ (i,i=1,m) /) if (initialize.eq.1) then read(5,*) a0 else ! ! If initial guess for consumption not given, set a0(1) = ! .9*(lambda-delta), a0(2) = .1*(lambda-delta), a0(j) = 0, j>2. ! a0(1,1) = 0.9*(lambda-delta) a0(2,1) = 0.1*(lambda-delta) endif ! ! Calculate grid on capital stocks using roots of Chebyshev polynomial. ! x = cos(pi*(fm-sm+0.5)/fm) k = 0.5*(kmax-kmin)*(1.0+x)+kmin phi = cos( matmul(sn, acos(transpose(x))) ) ! ! The Newton-Raphson iteration, i.e., find a such that g(a)=0 where ! a are coefficients of the weighted residual approximation. ! newton: do ! ! Initialize f and its derivative df. ! f = 0.0 df = 0.0 ! ! Construct f(a) = int R(x,a) phi_i(x) dx and its derivatives with ! respect to a. ! c = matmul(transpose(phi),a0) kt = lambda*k**alpha+delt1*k-c xt = 2.0*(kt-kmin)/(kmax-kmin)-1.0 psi = cos( matmul(sn, acos(transpose(xt))) ) ct = matmul(transpose(psi),a0) res = c**gamma-beta*ct**gamma*(alpha*lambda*kt**alph1+delt1) dpsi= matmul(sn,em) * sin(matmul(sn, acos(transpose(xt))))/ & matmul(transpose(en), sqrt(em-transpose(xt*xt))) dcta= transpose(psi)-2.0*(matmul(matmul(transpose(dpsi),a0),en)) & * transpose(phi)/(kmax-kmin) dres= gamma*matmul(c**gamm1,en) * transpose(phi) - beta*gamma* & matmul(ct**gamm1*(alpha*lambda*kt**alph1+delt1),en) * & dcta + 2.0*beta*alpha*alph1*lambda*matmul(ct**gamm1* & kt**(alph1-1.0),en) * transpose(phi)/(kmax-kmin) f = matmul(phi,res) df = matmul(phi,dres) ! ! Solve df\f and put solution in f. ! call sgesv(n,1,df,n,iwork,f,n,info) ! ! Issue warnings if necessary. ! if (info.lt.0) write(*,*) 'WARNING: illegal value for', & ' argument ',-info,' in JUDD' if (info.gt.0) write(*,*) 'WARNING: singular matrix', & ' when solving df*x=f in JUDD' ! ! Do Newton update. ! sum1 = sum(f*f) a1 = a0-f sum1 = sqrt(sum1)/fn ! ! Check to see if the solution is converged. If so, stop. ! if (sum1.lt.crit) exit newton a0 = a1 end do newton ! ! Write new input file with results. ! 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,1) do i=2,n write(7,1) a1(i,1) end do ! ! Write solution to output file. ! 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 i=1,min(n,9) write(8,'('' a'',I1,'' = '',F7.4,'' '')') i,a1(i,1) end do do i=10,n write(8,'('' a'',I2,'' = '',F7.4,'' '')') i,a1(i,1) end do write(8,*) write(8,*) ' Capital and consumption: ' do i=1,m write(8,2) k(i,1),c(i,1) end do ! ! Format statements. ! 1 format (9X,F12.6) 2 format (4X,F7.4,3X,F7.4) end program judd90