aiyagari_disc.f90 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82

Page 1

!********************************************************************* ! file name : aiyagari_disc.f90 ! programmer : makoto nakajima ([email protected]) ! date : december 8, 2006 ! description : solving aiyagari (1994) economy using discretization ! of value function. double precision program. !********************************************************************* !********************************************************************* ! module: globvar !********************************************************************* module globvar !******* prohibit implicit variable declaration ******* implicit none !******* meta parameter ******* real(8),parameter:: eps = 1.0d-10 integer,parameter:: scr = 8 integer,parameter:: out = 10

!really small number !screen output file number !output file number

!******* model parameters ******* real(8),parameter:: beta = 0.960_8 real(8),parameter:: alpha = 0.360_8 real(8),parameter:: sigma = 2.000_8 real(8),parameter:: delta = 0.100_8 real(8),parameter:: kyguess = 3.000_8

!discount factor !capital share of output !curvature of utility func !depreciation rate !guess for eqm cap/y ratio

!******* variables for approximation ******* integer,parameter:: ne=7 !number of earnings shocks real(8),dimension(1:ne):: eval !grid points of shocks real(8),dimension(1:ne,1:ne):: pe !transition matrix integer,parameter:: na = 2000 !number of grid points real(8),parameter:: apara = 2.0_8 !control distance bw grids real(8),parameter:: aub = 50.0_8 !upperbound of capital space real(8),parameter:: alb = 0.0_8 !lowerbound of capital space real(8),dimension(1:na):: aval !grids on capital space !******* endogenous objects ******* integer,parameter:: ni=ne*na real(8),dimension(1:ni):: val0,val1 real(8),dimension(1:ne,1:na):: expval integer,dimension(1:ni):: astar real(8),dimension(1:ni):: dist0,dist1

!value function !expected value !optimal decision rule !type distribution

!******* for conversion between (e,a) and integer,dimension(ni):: iafun,iefun integer,dimension(ne,na):: iifun !******* for iteration ******* integer,parameter:: maxiter0 = integer,parameter:: maxiter1 = integer,parameter:: maxiter2 = real(8),parameter:: maxerr0 = real(8),parameter:: maxerr1 = real(8),parameter:: maxerr2 = real(8),parameter:: update0 =

2000 2000 2000 1.0d-3 1.0d-5 1.0d-10 0.10_8

i *******

!max number of price iter !max number of value iter !max number of dist iter !max tolerance for price iter !max tolerance for value iter !max tolerance for dist iter !update spped for price iter

!******* prices and aggregates ******* real(8):: ky0,ky1,kl,cap,lab,ra,wa end module globvar !********************************************************************* ! program: main !********************************************************************* program main !******* use global variable ******* use globvar !******* prohibit implicit variable declaration ******* implicit none !******* local variable declaration ******* integer:: iter0,iter1,iter2,i0,e0,e1,a0 real(8):: err0,err1,err2,cash0 real(8),dimension(na):: vvec,cvec logical,dimension(na):: cfea integer,dimension(1):: adummy !******* open screen output file *******

aiyagari_disc.f90 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164

Page 2

if (scr/=6) open(scr,file='./aiyagari_scr.out') !******* prepare conversion between i and (e,a) ******* i0=0 do e0=1,ne do a0=1,na i0=i0+1 iifun(e0,a0)=i0 iefun(i0)=e0 iafun(i0)=a0 end do end do !******* set earnings shock (exogenously obtained) ******* eval(1:ne)=& (/-1.20014,-0.75736,-0.36941,0.00000,0.36941,0.75736,1.20014/) pe(1,1:ne)=(/0.2019378,0.5034201,0.2568107,0.0363963,0.0014232,& 0.0000119,0.0000000/) pe(2,1:ne)=(/0.0411231,0.3233636,0.4513294,0.1667877,0.0170059,& 0.0003894,0.0000010/) pe(3,1:ne)=(/0.0057942,0.1246568,0.4202567,0.3596454,0.0849178,& 0.0046970,0.0000321/) pe(4,1:ne)=(/0.0005483,0.0307571,0.2401232,0.4571429,0.2401232,& 0.0307571,0.0005483/) pe(5,1:ne)=(/0.0000321,0.0046970,0.0849178,0.3596454,0.4202567,& 0.1246568,0.0057942/) pe(6,1:ne)=(/0.0000010,0.0003894,0.0170059,0.1667877,0.4513294,& 0.3233636,0.0411231/) pe(7,1:ne)=(/0.0000000,0.0000119,0.0014232,0.0363963,0.2568107,& 0.5034201,0.2019378/) !******* construct array of grid points ******* aval(1)=alb aval(na)=aub do a0=2,na-1 aval(a0)=alb+(aub-alb)*(real(a0-1,8)/real(na-1,8))**apara end do !******* set initial guess for prices ******* ky0=kyguess kl=ky0**(1.0_8/(1.0_8-alpha)) ra=1.0_8-delta+alpha*kl**(alpha-1.0_8) wa=(1.0_8-alpha)*kl**alpha !******* set initial guess for value function ******* do i0=1,ni val0(i0)=max(exp(eval(iefun(i0)))*wa+aval(iafun(i0))*ra-alb,eps)& **(1.0_8-sigma)/(1.0_8-sigma) end do !******* set stupid initial guess for distribution ******* dist0=0.0_8 dist0(iifun(1,1))=1.0_8 !******* beginning of the loop for prices ******* do iter0=1,maxiter0 !******* print pout current prices ******* write(scr,"('price iter=',i5,' r =',f15.12,' w iter0,ra,wa

=',f15.12,/)") &

!******* beginning of the loop for value iteration ******* do iter1=1,maxiter1 !******* construct expected value ******* expval=0.0_8 do e0=1,ne do a0=1,na do e1=1,ne expval(e0,a0)=expval(e0,a0)& +beta*pe(e0,e1)*val0(iifun(e1,a0)) end do end do end do !******* loop for eash state ******* do i0=1,ni !******* update value ******* cash0=exp(eval(iefun(i0)))*wa+aval(iafun(i0))*ra cvec=cash0-aval

aiyagari_disc.f90 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246

Page 3

cfea=(cvec>eps) where (cfea==.true.) vvec=cvec**(1.0_8-sigma)/(1.0_8-sigma)+expval(iefun(i0),:) end where adummy=maxloc(vvec,cfea) astar(i0)=adummy(1) val1(i0)=vvec(astar(i0)) !******* end of loop for each state ******* end do !******* compute error ******* err1=maxval(abs(val0-val1)) !******* print out iteration result ******* write(scr,"(' value iter=',i5,' err=',f15.12,' tol=',f15.12,& ' conv?=',l1)") iter1,err1,maxerr1,(err1
aiyagari_disc.f90

Page 4

247 end do 248 close(out) 249 250 !******* compute new capital/output ratio ******* 251 cap=0.0_8 252 lab=0.0_8 253 do i0=1,ni 254 cap=cap+aval(iafun(i0))*dist0(i0) 255 lab=lab+exp(eval(iefun(i0)))*dist0(i0) 256 end do 257 ky1=(cap/lab)**(1.0_8-alpha) 258 259 !******* compute error ******* 260 err0=abs(ky0-ky1) 261 262 !******* print out iteration result ******* 263 write(scr,"(/,'price iter=',i5,' err=',f15.12,' tol=',f15.12,& 264 ' conv?=',l1)") iter0,err0,maxerr0,(err0