-
Notifications
You must be signed in to change notification settings - Fork 1
/
calibration.f90
183 lines (163 loc) · 8.43 KB
/
calibration.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
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
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! THIS MODULE CONTAINS THE FUNCTIONS AND SUBROUTINES REQUIRED TO CALIBRATE THE PARAMETERS
! OF THE MODEL
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MODULE calibration
USE params
USE firmsproblem , ONLY : SOLVEPROBLEM,FIND_EQUIL_wrate
IMPLICIT NONE
REAL(rp) :: MAXERRORCALIB
CONTAINS
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! This subroutine normalize the parameter values to lie within a specific range
SUBROUTINE SETPARAMS(VECPARS,INDICATOR)
USE toolkit , ONLY : NORMALIZE
IMPLICIT NONE
INTEGER , INTENT(IN) :: INDICATOR
REAL(rp) , INTENT(INOUT) :: VECPARS(:)
INTEGER :: i ; i = 1
CALL NORMALIZE ( VECPARS(i) , fc , DBLE(0.800) , DBLE(0.000) , INDICATOR ) ; i = i + 1 ! 1
CALL NORMALIZE ( VECPARS(i) , sigma_z , DBLE(0.500) , DBLE(0.010) , INDICATOR ) ; i = i + 1 ! 2
CALL NORMALIZE ( VECPARS(i) , mu_z , DBLE(0.150) , DBLE(0.010) , INDICATOR ) ; i = i + 1 ! 3
CALL NORMALIZE ( VECPARS(i) , bigA , DBLE(4.000) , DBLE(1.000) , INDICATOR ) ; i = i + 1 ! 4
CALL NORMALIZE ( VECPARS(i) , sigma0_z , DBLE(2.000) , DBLE(0.000) , INDICATOR ) ; i = i + 1 ! 5
CALL NORMALIZE ( VECPARS(i) , rho_0 , DBLE(1.000) , DBLE(0.000) , INDICATOR ) ; i = i + 1 ! 6
CALL NORMALIZE ( VECPARS(i) , kappa_0 , DBLE(70.00) , DBLE(0.000) , INDICATOR ) ; i = i + 1 ! 7
CALL NORMALIZE ( VECPARS(i) , kappa_1 , DBLE(3.000) , DBLE(-1.00) , INDICATOR ) ; i = i + 1 ! 8
RETURN
END SUBROUTINE SETPARAMS
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! This subroutine controls the calibration process
SUBROUTINE CALIBRATE( )
USE toolkit , ONLY : SIMPLEX,LMMIN,TIMING
IMPLICIT NONE
INTEGER , PARAMETER :: NUMP = 8 , SHOW0 = 2
REAL(rp) , PARAMETER :: SHOCK0 = 0.05 , DAMP0 = 0.10
REAL(rp) , DIMENSION(NUMP) :: PARS1,PARS0,PARSS,PARSSS,SHOCKS
REAL(rp) :: TERROR,TERROR0,TERROR1
REAL(rp) :: RES(SIZE(MOMS))
REAL(rp) :: TTIME,MAXTIME
INTEGER :: ECODE,NITER,IITER,MAXITER
WRITE(*,FMT="(A)",ADVANCE="YES") ' ---------------------------------------------------------------'
WRITE(*,FMT="(A)",ADVANCE="YES") ' '
WRITE(*,FMT="(A)",ADVANCE="YES") ' CALIBRATION '
WRITE(*,FMT="(A)",ADVANCE="YES") ' '
WRITE(*,FMT="(A)",ADVANCE="YES") ' Choose calibratrion algorithm '
WRITE(*,FMT="(A)",ADVANCE="YES") ' '
WRITE(*,FMT="(A)",ADVANCE="YES") ' [ 1 ] Levenberg-Marquardt '
WRITE(*,FMT="(A)",ADVANCE="YES") ' [ 2 ] Simplex '
WRITE(*,FMT="(A)",ADVANCE="YES") ' [ 3 ] Random search '
WRITE(*,FMT="(A)",ADVANCE="YES") ' [ 4 ] Levenberg-Marquardt, shocking initial values '
WRITE(*,FMT="(A)",ADVANCE="YES") ' [ 5 ] Simplex, shocking initial values '
WRITE(*,FMT="(A)",ADVANCE="YES") ' '
WRITE(*,FMT="(A)",ADVANCE="NO" ) ' Your choice: '; READ (*,*) CALIBMOD
IF (CALIBMOD.le.2) THEN
WRITE(*,FMT="(A)",ADVANCE="NO" ) ' Max iterations: '; READ (*,*) MAXITER
ELSEIF (CALIBMOD.eq.3) THEN
WRITE(*,FMT="(A)",ADVANCE="NO" ) ' Max time (hours): '; READ (*,*) MAXTIME
ELSEIF (CALIBMOD.gt.3) THEN
WRITE(*,FMT="(A)",ADVANCE="NO" ) ' Max iterations (each param): '; READ (*,*) MAXITER
END IF
WRITE(*,FMT="(A)",ADVANCE="YES") ' '
WRITE(*,FMT="(A)",ADVANCE="YES") ' ---------------------------------------------------------------'
WRITE(*,FMT="(A)",ADVANCE="YES") ' '
! Start timing
timing0 = TIMING(1)
bigA = (DBLE(4.10)**(one-gamma))/gamma
! Set initial values
CALL SETPARAMS(PARS0,0)
! Initialize max. error
MAXERRORCALIB = 999999999.999
! Levenberg-Marquardt algorithm
IF(INT(CALIBMOD).EQ.1) THEN
CALL LMMIN(FUNC_CALIB,PARS1,RES,NITER,ECODE,PARS0,SHCK=SHOCK0,DAMP=DAMP0,IPRINT=SHOW0)
! SIMPLEX algorithm
ELSE IF (INT(CALIBMOD).EQ.2) THEN
CALL SIMPLEX(SUM_FUNC_CALIB,PARS1,TERROR,NITER,ECODE,PARS0,IPRINT=SHOW0)
! Random search
ELSE IF ( INT(CALIBMOD) .EQ. 3 ) THEN
TTIME = TIMING(3)
TERROR0 = SUM_FUNC_CALIB(PARS0) ! Initial error
TERROR1 = TERROR0 ! Best error so far
TERROR = TERROR0 ! Error current iteration
PRINT ('(I6,5(F12.6))') , NITER , TIMING(3)-TTIME , TERROR1 , TERROR , TERROR0
DO WHILE ( TIMING(3)-TTIME .LT. MAXTIME ) ; CALL RANDOM_NUMBER(SHOCKS)
PARS1(:) = PARS0(:) + PARS0(:)*SHOCK0*(SHOCKS(:)-half)
TERROR1 = SUM_FUNC_CALIB(PARS1)
PRINT ('(I6,5(F12.6))') , NITER , TIMING(3)-TTIME , TERROR1 , TERROR , TERROR0
IF (TERROR1.LE.TERROR) THEN
TERROR = TERROR1
PARS0 = PARS1
END IF
END DO
PARS1 = PARS0
! Levenberg-Marquardt or Simplex , shocking initial values
ELSE IF ( INT(CALIBMOD) .GT. 3 ) THEN ! Simplex with large shocks
IF (INT(CALIBMOD).EQ.4) CALL LMMIN(FUNC_CALIB,PARSS,RES,NITER,ECODE,PARS0,SHCK=SHOCK0,DAMP=DAMP0,IPRINT=SHOW0)
IF (INT(CALIBMOD).EQ.5) CALL SIMPLEX(SUM_FUNC_CALIB,PARSS,TERROR,NITER,ECODE,PARS0,IPRINT=SHOW0)
PARSSS = PARSS
1 TERROR = SUM_FUNC_CALIB(PARSS)
PRINT * , 0 , TERROR
DO IITER=1,NUMP
PARS0 = PARSS ; PARS0(IITER) = PARSS(IITER)*DBLE(1.50)
IF (INT(CALIBMOD).EQ.4) CALL LMMIN(FUNC_CALIB,PARS1,RES,NITER,ECODE,PARS0,SHCK=SHOCK0,DAMP=DAMP0,IPRINT=SHOW0)
IF (INT(CALIBMOD).EQ.5) CALL SIMPLEX(SUM_FUNC_CALIB,PARS1,TERROR0,NITER,ECODE,PARS0,IPRINT=SHOW0)
TERROR0 = SUM_FUNC_CALIB(PARS1)
PRINT * , IITER , TERROR0 , TERROR
IF (TERROR0.LT.TERROR) THEN
PARSSS = PARS1
TERROR = TERROR0
END IF
PARS0 = PARSS ; PARS0(IITER) = PARSS(IITER)*DBLE(0.50)
IF (INT(CALIBMOD).EQ.4) CALL LMMIN(FUNC_CALIB,PARS1,RES,NITER,ECODE,PARS0,SHCK=SHOCK0,DAMP=DAMP0,IPRINT=SHOW0)
IF (INT(CALIBMOD).EQ.5) CALL SIMPLEX(SUM_FUNC_CALIB,PARS1,TERROR0,NITER,ECODE,PARS0,IPRINT=SHOW0)
TERROR0 = SUM_FUNC_CALIB(PARS1)
PRINT * , -IITER , TERROR0 , TERROR
IF (TERROR0.LT.TERROR) THEN
PARSSS = PARS1
TERROR = TERROR0
END IF
END DO
IF (SUM(ABS(PARSSS-PARSS)).GT.TOL) THEN
PRINT * , ' ' ; PRINT * , ' ' ; PRINT * , ' Re-initialize ' ; PRINT * , ' '
PARSS = PARSSS ; GOTO 1
END IF
PARS1 = PARSSS
END IF
! Fix resulting parameter values
CALL SETPARAMS(PARS1,1)
! Solve the model with new parameters, assuming wrate=1
CALL SOLVEPROBLEM( )
theta = FINDTHETA(TOTL,TOTY-TOTLI,one)
RETURN
END SUBROUTINE CALIBRATE
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! This function returns the vector differences between the model- and data-generated
! moments. It is used when applying the Levenberg-Marquardt method
FUNCTION FUNC_CALIB(VECPARS) RESULT(ER)
IMPLICIT NONE
REAL(rp) :: VECPARS(:),SUMER
REAL(rp) , ALLOCATABLE :: ER(:)
ALLOCATE(ER(SIZE(MOMS)))
CALL SETPARAMS(VECPARS,1)
CALL SOLVEPROBLEM( )
ER(:) = MOMS(:)*WMAT(:)
SUMER = SUM(ER(:)*ER(:))
IF ( SUMER.LT.MAXERRORCALIB ) THEN
theta = FINDTHETA(TOTL,TOTY-TOTLI,one)
CALL WRITE_RESULTS(1,"calibration.txt")
MAXERRORCALIB = SUMER
END IF
RETURN
END FUNCTION FUNC_CALIB
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! This function returns the sum of squared differences between the model- and
! data-generated moments. It is used when applying the Simplex algorithm
FUNCTION SUM_FUNC_CALIB(VECPARS) RESULT(SUMER)
IMPLICIT NONE
REAL(rp) :: VECPARS(:),SUMER,ER(SIZE(MOMS))
ER = FUNC_CALIB(VECPARS) ; SUMER = SUM(ER(:)*ER(:))
RETURN
END FUNCTION SUM_FUNC_CALIB
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
END MODULE calibration