-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathps5.f95
314 lines (282 loc) · 11.5 KB
/
ps5.f95
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
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
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
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
MODULE PS5PARA
IMPLICIT NONE
INTEGER:: I
INTEGER, PARAMETER:: N=66, RETIRE_AGE=46 ! LIFE EXPECTANCY AND RETIREMENT
REAL(KIND=8), PARAMETER:: A_MAX=15.0, A_MIN=0.1, STEP=0.1
INTEGER, PARAMETER:: NZ=2, NA=(A_MAX-A_MIN)/STEP+1
REAL(KIND=8):: ALPHA=0.36, BETA=0.97
REAL(KIND=8), PARAMETER:: THETA=0.0, GAMMA=0.42, SIGMA=2.0, DELTA=0.06, POP_GROWTH= 0.011
REAL(KIND=8), DIMENSION(2):: PROD_STATES = (/3.0, 0.5/)
REAL(KIND=8), DIMENSION(2,2):: PROD_MARKOV = TRANSPOSE(RESHAPE((/0.9261,(1-0.9261),(1-0.9811),0.9811/),(/NZ,NZ/)))
REAL(KIND=8), DIMENSION(NA):: A_GRID = (/(I*STEP, I=1,NA)/) + A_MIN - STEP
REAL(KIND=8), DIMENSION(N):: COHORT_POP
REAL(KIND=8), DIMENSION(N,NZ):: AGE_EFF=0.
REAL(KIND=8), DIMENSION(45):: EF
CONTAINS
SUBROUTINE INIT_ARRAY()
COHORT_POP = (/((1.0/((1.0+POP_GROWTH)**I)), I=1,N)/)
COHORT_POP = COHORT_POP/SUM(COHORT_POP)
OPEN(UNIT=13, FILE='ef.txt', STATUS='OLD')
DO I=1,RETIRE_AGE-1
READ(UNIT=13, FMT=*) EF(I)
ENDDO
CLOSE(UNIT=13)
DO I=1,RETIRE_AGE-1
AGE_EFF(I,1) = EF(I)*PROD_STATES(1)
AGE_EFF(I,2) = EF(I)*PROD_STATES(2)
ENDDO
END SUBROUTINE
END MODULE
MODULE PS5RES
use PS5PARA
IMPLICIT NONE
REAL(KIND=8):: WAGE, AGG_LABOR=1.0, AGG_CAP=2.0, BENEFIT, INTEREST, ERROR=100.
REAL(KIND=8), DIMENSION(NA,NZ,N):: VFUNC, LFUNC, STAT_DIST, VFUNC_OLD
INTEGER, DIMENSION(NA,NZ,N):: PFUNC=1
END MODULE
PROGRAM PS5
USE PS5PARA
USE PS5RES
IMPLICIT NONE
REAL(KIND=8):: CRIT_P=1e-3, WEL, CV
INTEGER:: SROWIDX, AGEIDX
CALL INIT_ARRAY()
PRINT*, NA
DO WHILE(ERROR>CRIT_P)
! UPDATE VARS
WAGE = (1-ALPHA)*(AGG_CAP**(ALPHA))*(AGG_LABOR**(-ALPHA))! FOC WRT LABOR
INTEREST = ALPHA*(AGG_CAP**(ALPHA-1))*(AGG_LABOR**(1-ALPHA))-DELTA! FOC WRT CAPITAL
BENEFIT = THETA*WAGE*AGG_LABOR/SUM(COHORT_POP(RETIRE_AGE:N)) ! TAX REV/RETIRED_POP
PRINT*, "=================================================="
PRINT*, "R, B,W:", INTEREST, BENEFIT, WAGE
PRINT*, "L, K:",AGG_LABOR, AGG_CAP
LFUNC = 0.
VFUNC=0.
PFUNC=1
! BACKWARD INDUCTION USING NEW VARS (BELLMAN)
CALL BACKWARD_INDUCTION()
! UPDATE STAT_DIST
CALL FIND_STAT_DIST()
! OBSERVE ERROR AND UPDATE PRICE
CALL UPDATE_PRICE()
ENDDO
CALL WRITE_ALL()
WAGE = (1-ALPHA)*(AGG_CAP**ALPHA)*(AGG_LABOR**(-ALPHA))! FOC WRT LABOR
INTEREST = ALPHA*(AGG_CAP**(ALPHA-1))*(AGG_LABOR**(1-ALPHA))-DELTA! FOC WRT CAPITAL
BENEFIT = THETA*WAGE*AGG_LABOR/SUM(COHORT_POP(RETIRE_AGE:N)) ! TAX REV/RETIRED_POP
DO SROWIDX=1,NA
DO AGEIDX=1,N
WEL=WEL+VFUNC(SROWIDX,1,AGEIDX)*STAT_DIST(SROWIDX,1,AGEIDX)+VFUNC(SROWIDX,2,AGEIDX)*STAT_DIST(SROWIDX,2,AGEIDX)
CV = CV + A_GRID(SROWIDX)*SUM(STAT_DIST(SROWIDX, :, AGEIDX))
ENDDO
ENDDO
PRINT*, "=================================================="
PRINT*, "R, B,W:", INTEREST, BENEFIT, WAGE
PRINT*, "L, K:",AGG_LABOR, AGG_CAP
PRINT*, "WELFARE", WEL
PRINT*, "WEALTH", CV
END PROGRAM PS5
! SUBROUTINES
SUBROUTINE FIND_LABOR(CHOICEIDX, SROWIDX, SCOLIDX, AGE, WORKING)
USE PS5RES
USE PS5PARA
IMPLICIT NONE
INTEGER, INTENT(IN):: CHOICEIDX, SROWIDX, SCOLIDX, AGE
REAL(KIND=8), INTENT(OUT):: WORKING
REAL(KIND=8):: LAST_CHOICE, CUR_CHOICE, OPT_LABOR, NOM
CUR_CHOICE = A_GRID(CHOICEIDX)
IF (AGE/=1) THEN
LAST_CHOICE = (1+INTEREST)*A_GRID(SROWIDX) ! DISCOUNTED
ELSE
LAST_CHOICE = 0.
ENDIF
NOM = (GAMMA*(1-THETA)*AGE_EFF(AGE,SCOLIDX)*WAGE-(1-GAMMA)*(LAST_CHOICE-CUR_CHOICE))
OPT_LABOR = NOM/((1-THETA)*WAGE*AGE_EFF(AGE,SCOLIDX))
IF (OPT_LABOR > 1.0) THEN
OPT_LABOR = 1.0
ELSE IF (OPT_LABOR < 0.) THEN
OPT_LABOR = 0.0
ENDIF
WORKING = OPT_LABOR
END SUBROUTINE
SUBROUTINE BACKWARD_INDUCTION()
USE PS5RES
USE PS5PARA
IMPLICIT NONE
INTEGER:: AGEIDX, AGE, SROWIDX
REAL(KIND=8):: CONSUM, UTIL, ERROR_VFI, CRIT_VFI=1e-3
ERROR_VFI=100.
VFUNC=0.
DO WHILE(ERROR_VFI> CRIT_VFI)
! FOR LAST AGE, ! SHOULD BE DETERMINISTIC, DO ONCE IS OK (NO CHOICE)
DO SROWIDX=1,NA ! THIS IS THE CHOICE FROM LAST PERIOD/ CURRENT ASSET
CONSUM = (1+INTEREST)*A_GRID(SROWIDX) + BENEFIT ! SPEND ALL IN LAST PERIOD
PFUNC(:,:,N) = 1 ! ACTUALLY DOES NOT MATTER
UTIL = CONSUM**((1-SIGMA)*GAMMA)/(1-SIGMA) ! NO NEXT PERIOD VALUE
VFUNC(SROWIDX,1,N) = UTIL
VFUNC(SROWIDX,2,N) = UTIL
! ASSET MUST BE ZERO FOR AGE
ENDDO
DO AGEIDX = 2, N ! FOR AGE 20-65
AGE = N- AGEIDX+1
CALL BELLMAN(AGE)
ENDDO
ERROR_VFI = MAXVAL(ABS(VFUNC-VFUNC_OLD))
VFUNC_OLD = VFUNC
ENDDO
END SUBROUTINE
SUBROUTINE BELLMAN(AGE)
USE PS5RES
USE PS5PARA
IMPLICIT NONE
INTEGER, INTENT(IN):: AGE
REAL(KIND=8):: COND_MAX_UTIL, CONSUM, UTIL, WORKING, PROD, NEXTU
INTEGER:: SROWIDX, SCOLIDX, CHOICEIDX
IF (AGE< RETIRE_AGE) THEN ! WOKRING BELLMAN
DO SROWIDX=1, NA
DO SCOLIDX=1, NZ
COND_MAX_UTIL = -1e12
PROD = AGE_EFF(AGE, SCOLIDX)
DO CHOICEIDX=1,NA ! LOOP OVER CHOICE OF ASSET PRIME
CALL FIND_LABOR(CHOICEIDX, SROWIDX, SCOLIDX, AGE, WORKING)
CONSUM = (1+INTEREST)*A_GRID(SROWIDX) + (1-THETA)*PROD*WORKING*WAGE - A_GRID(CHOICEIDX)
IF (CONSUM > 0.) THEN
NEXTU = BETA*SUM(PROD_MARKOV(SCOLIDX,:)*VFUNC(CHOICEIDX,:,AGE+1))
UTIL = (CONSUM**GAMMA*(1-WORKING)**(1-GAMMA))**(1-SIGMA)/(1-SIGMA)+NEXTU
IF (UTIL>COND_MAX_UTIL) THEN
PFUNC(SROWIDX, SCOLIDX, AGE) = CHOICEIDX
LFUNC(SROWIDX,SCOLIDX, AGE) = WORKING
COND_MAX_UTIL = UTIL
ENDIF
ENDIF
ENDDO ! END LOOP CHOICE SPACE FOR ONE STATE
IF (COND_MAX_UTIL<-1E10) THEN
PRINT*, VFUNC(PFUNC(SROWIDX,1, AGE),1,AGE+1), VFUNC(PFUNC(SROWIDX,2, AGE),1,AGE+1), AGE, SROWIDX
PRINT*, COND_MAX_UTIL, CONSUM, UTIL, (CONSUM**GAMMA*(1-WORKING)**(1.-GAMMA)), WORKING
STOP
ENDIF
VFUNC(SROWIDX, SCOLIDX, AGE) = COND_MAX_UTIL
ENDDO
ENDDO
ELSE ! BELLMAN FOR RETIRED PEOPLE
DO SROWIDX=1, NA
COND_MAX_UTIL = -1e12
DO CHOICEIDX=1, NA
CONSUM = (1+INTEREST)*A_GRID(SROWIDX) + BENEFIT - A_GRID(CHOICEIDX)
IF (CONSUM>0.) THEN
UTIL = CONSUM**((1-SIGMA)*GAMMA)/(1-SIGMA)+ BETA*(VFUNC(CHOICEIDX,1,AGE+1))
IF (UTIL>COND_MAX_UTIL) THEN
PFUNC(SROWIDX,:, AGE) = CHOICEIDX ! NO PRODUCTIVITY DIFFERENCE
COND_MAX_UTIL = UTIL
ENDIF
ENDIF
ENDDO
IF (COND_MAX_UTIL<-1E10) THEN
PRINT*, VFUNC(PFUNC(SROWIDX,1, AGE),1,AGE+1), VFUNC(PFUNC(SROWIDX,2, AGE),1,AGE+1), AGE, SROWIDX
PRINT*, COND_MAX_UTIL
STOP
ENDIF
VFUNC(SROWIDX, 1, AGE) = COND_MAX_UTIL ! PROD_STATE DO NOT AFFECT RETIRED PEOPLE
VFUNC(SROWIDX, 2, AGE) = COND_MAX_UTIL ! PROD_STATE DO NOT AFFECT RETIRED PEOPLE
ENDDO
ENDIF
END SUBROUTINE
SUBROUTINE FIND_STAT_DIST()
USE PS5PARA
USE PS5RES
IMPLICIT NONE
INTEGER:: SROWIDX, SCOLIDX, AGEIDX, LAST_AGE, LAST_CHOICE
REAL(KIND=8), DIMENSION(NA,NZ,N):: STAT_DIST_NEW
PRINT*, COUNT(STAT_DIST/=0), SUM(STAT_DIST), MAXLOC(STAT_DIST), SUM(COHORT_POP)
! EVERYONE HAVE ZERO ASSET AT STARTING AGE
STAT_DIST_NEW = 0.
STAT_DIST_NEW(1,1,1) = 0.2037 ! ALWAYS NO ASSET AT FIRST AGE,
STAT_DIST_NEW(1,2,1) = 0.7963 ! BUT 2 PROD STATES
DO AGEIDX=2, N
DO SROWIDX=1, NA ! LAST PERIOD ASSET
DO SCOLIDX=1, NZ !LAST PERIOD PRODUCTIVITY
LAST_CHOICE = PFUNC(SROWIDX, SCOLIDX, AGEIDX-1) ! TWO POSSIBLE
LAST_AGE = AGEIDX-1
STAT_DIST_NEW(LAST_CHOICE, 1, AGEIDX) = &
STAT_DIST_NEW(LAST_CHOICE, 1, AGEIDX)+ STAT_DIST_NEW(SROWIDX,SCOLIDX,LAST_AGE)*PROD_MARKOV(SCOLIDX,1)
STAT_DIST_NEW(LAST_CHOICE, 2, AGEIDX) = &
STAT_DIST_NEW(LAST_CHOICE, 2, AGEIDX) + STAT_DIST_NEW(SROWIDX,SCOLIDX,LAST_AGE)*PROD_MARKOV(SCOLIDX,2)
ENDDO
ENDDO
ENDDO
DO AGEIDX=1,N
STAT_DIST(:,:,AGEIDX) = STAT_DIST_NEW(:,:,AGEIDX)*COHORT_POP(AGEIDX)
ENDDO
END SUBROUTINE
SUBROUTINE UPDATE_PRICE()
! COMPUTE ERROR AND UPDATE PRICE
USE PS5RES
USE PS5PARA
IMPLICIT NONE
INTEGER:: SROWIDX, SCOLIDX, AGEIDX
REAL(KIND=8):: AGG_LABOR_NEW=0., AGG_CAP_NEW=0., PORTION, ERROR_LAB, ERROR_CAP
AGG_LABOR_NEW=0.
AGG_CAP_NEW=0.
DO AGEIDX=1, N ! CAL ERRORS
DO SROWIDX=1,NA
DO SCOLIDX=1,NZ
PORTION = STAT_DIST(SROWIDX, SCOLIDX, AGEIDX)
AGG_CAP_NEW = AGG_CAP_NEW+ A_GRID(SROWIDX)*PORTION
IF (AGEIDX<RETIRE_AGE) THEN
AGG_LABOR_NEW = AGG_LABOR_NEW+LFUNC(SROWIDX,SCOLIDX,AGEIDX)*PORTION*AGE_EFF(AGEIDX,SCOLIDX)
ENDIF
ENDDO
ENDDO
ENDDO
ERROR_LAB = ABS(AGG_LABOR_NEW-AGG_LABOR)
ERROR_CAP = ABS(AGG_CAP_NEW-AGG_CAP)
PRINT*, "UPDATE_PRICE ERROR", ERROR_LAB, ERROR_CAP
ERROR = MAX(ERROR_CAP, ERROR_LAB)
AGG_LABOR = 0.01*AGG_LABOR_NEW + 0.99*AGG_LABOR
AGG_CAP = 0.01*AGG_CAP_NEW + 0.99*AGG_CAP
! UPDATE INTEREST, WAGE, BENEFIT AT THE START OF NEXT PERIOD (IT IS THE SAME)
END SUBROUTINE
SUBROUTINE WRITE_ALL()
USE PS5RES
USE PS5PARA
IMPLICIT NONE
INTEGER:: SROWIDX
CHARACTER(LEN=130):: PATH="/Users/chek_choi/Downloads/fortran/"
CHARACTER(LEN=150):: FILE_NAME
FILE_NAME = TRIM(PATH)//"VFUNC"
OPEN(UNIT=1, FILE=FILE_NAME, STATUS='REPLACE') ! START WITH THE TWO VALUE FUNCTIONS
DO SROWIDX=1, NA
WRITE(UNIT=1,FMT=*) VFUNC(SROWIDX,:, 50)
ENDDO
CLOSE(UNIT=1)
FILE_NAME = TRIM(PATH)//"PFUNC"
OPEN(UNIT=2, FILE=FILE_NAME, STATUS='REPLACE') ! ALSO SAVE POLICY FUNCTIONS
DO SROWIDX=1, NA
WRITE(UNIT=2,FMT=*) A_GRID(PFUNC(SROWIDX,1, 20)), A_GRID(PFUNC(SROWIDX,2, 20))
ENDDO
CLOSE(UNIT=2)
FILE_NAME = TRIM(PATH)//"STATDIST"
OPEN(UNIT=3, FILE=FILE_NAME, STATUS='REPLACE') ! ALSO SAVE POLICY FUNCTIONS
DO SROWIDX=1, NA
WRITE(UNIT=3,FMT=*) STAT_DIST(SROWIDX, :, 2)
ENDDO
CLOSE(UNIT=3)
FILE_NAME = TRIM(PATH)//"LFUNC"
OPEN(UNIT=4, FILE=FILE_NAME, STATUS='REPLACE') ! FOR HAVING THE X-AXIS OF PLOT
DO SROWIDX=1, NA
WRITE(UNIT=4,FMT=*) LFUNC(SROWIDX, :, 50)
ENDDO
CLOSE(UNIT=4)
FILE_NAME = TRIM(PATH)//"AGRID"
OPEN(UNIT=5, FILE=FILE_NAME, STATUS='REPLACE') ! START WITH THE TWO VALUE FUNCTIONS
DO SROWIDX=1, NA
WRITE(UNIT=5,FMT=*) A_GRID(SROWIDX)
ENDDO
CLOSE(UNIT=5)
FILE_NAME = TRIM(PATH)//"COHORTPOP"
OPEN(UNIT=32, FILE=FILE_NAME, STATUS='REPLACE') ! START WITH THE TWO VALUE FUNCTIONS
DO SROWIDX=1, N
WRITE(UNIT=32,FMT=*) COHORT_POP(SROWIDX)
ENDDO
CLOSE(UNIT=32)
END SUBROUTINE