1C ------------------------------------------------------------- 2C CHECKDESIGN cheks design.dll or design.m 3C Developed by A.Rossi, C.Planas and G.Fiorentini 4C 5C State-space format: y(t) = c(t)z(t) + H(t)x(t) + G(t)u(t) 6C x(t) = a(t) + F(t)x(t-1) + R(t)u(t) 7C 8C y(t) (ny x 1) ny = # of endogenous series 9C z(t) (nz x 1) nz = # of exogenous series 10C x(t) (nx x 1) nx = # of continous states 11C u(t) (nu x 1) nu = # of shocks 12C c(t) (ny x nz x ns1) ns1 = # of states for c(t) 13C H(t) (ny x nx x ns2) ns2 = # of states for H(t) 14C G(t) (ny x nu x ns3) ns3 = # of states for G(t) 15C a(t) (nx x ns4) ns4 = # of states for a(t) 16C F(t) (nx x nx x ns5) ns5 = # of states for F(t) 17C R(t) (nx x nu x ns6) ns6 = # of states for R(t) 18C 19C Copyright (C) 2010-2014 European Commission 20C 21C This file is part of Program DMM 22C 23C DMM is free software developed at the Joint Research Centre of the 24C European Commission: you can redistribute it and/or modify it under 25C the terms of the GNU General Public License as published by 26C the Free Software Foundation, either version 3 of the License, or 27C (at your option) any later version. 28C 29C DMM is distributed in the hope that it will be useful, 30C but WITHOUT ANY WARRANTY; without even the implied warranty of 31C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 32C GNU General Public License for more details. 33C 34C You should have received a copy of the GNU General Public License 35C along with DMM. If not, see <http://www.gnu.org/licenses/>. 36C ------------------------------------------------------------- 37 SUBROUTINE CHECKDESIGN(ny,nz,nx,nu,ns,nt,d,theta,pdll,PATH,NMLNAME) 38 39 USE dfwin 40 INTERFACE 41 SUBROUTINE DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R) 42 INTEGER ny,nz,nx,nu,ns(6),nt 43 DOUBLE PRECISION theta(nt) 44 DOUBLE PRECISION c(ny,max(1,nz),ns(1)),H(ny,nx,ns(2)), 45 1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6)) 46 END SUBROUTINE 47 END INTERFACE 48 POINTER (pdll,fittizia) ! ASSOCIATE pointer pdll alla DLL ad una varibile fittizia 49 POINTER (pdesign,DESIGN) 50 51C INPUT 52 INTEGER ny,nz,nx,nu,ns(6),nt,d(2) 53 DOUBLE PRECISION theta(nt) 54 CHARACTER*200 NMLNAME,PATH,FILEOUT 55 56C LOCALS 57 INTEGER I,J,maxnz,IFAIL,ESTABLE 58 DOUBLE PRECISION WORK(4*nx),WR(nx),WI(nx),VR(1), 59 1 VI(1),W(nx) !WRY(ny),WORK1(64*ny) 60 DOUBLE PRECISION,ALLOCATABLE::c(:,:,:),H(:,:,:), 61 1 G(:,:,:),a(:,:),F(:,:,:),R(:,:,:) !,HRG(:,:),HRGRH(:,:) 62 CHARACTER*3 CJ 63C EXTERNAL SUBROUTINES 64 EXTERNAL DGEEV 65 66 ALLOCATE(c(ny,max(nz,1),ns(1)),H(ny,nx,ns(2)), 67 1 G(ny,nu,ns(3)),a(nx,ns(4)),F(nx,nx,ns(5)),R(nx,nu,ns(6))) !,HRG(ny,nu),HRGRH(ny,ny)) 68 69 pdesign = getprocaddress(pdll, "design_"C) 70 CALL DESIGN(ny,nz,nx,nu,ns,nt,theta,c,H,G,a,F,R) 71 72 maxnz = max(1,nz) 73 FILEOUT = TRIM(PATH)//TRIM(NMLNAME)//'.CHK' 74 OPEN(11,FILE = FILEOUT, ACCESS='SEQUENTIAL') 75 76C write ny,nx,etc 77 WRITE(11,1000) nt,ny,nx,nu,nz 781000 FORMAT(' nt = ',I3,'; ny = ',I3,'; nx = ',I3,'; nu = ',I3, 79 # '; nz = ',I3,';') 80 81C write theta 82 WRITE(11,*) 'theta(1:nt) = [' 83 WRITE(11,*) ' ' 84 WRITE(11,'(<nt>(F20.10))') theta(1:nt) 85 WRITE(11,*) ']' 86 87C write c(ny,max(1,nz),ns(1)) 88 DO J = 1,ns(1) 89 WRITE(11,*) ' ' 90 IF (J.LE.9) THEN 91 WRITE(CJ,'(I1)') J 92 ELSEIF ((J.GE.10).AND.(J.LE.99)) THEN 93 WRITE(CJ,'(I2)') J 94 ELSEIF ((J.GE.100).AND.(J.LE.999)) THEN 95 WRITE(CJ,'(I3)') J 96 ENDIF 97 WRITE(11,*) 'c(1:ny,1:nz,'//TRIM(CJ)// ') = [' 98 WRITE(11,*) ' ' 99 WRITE(11,'(<maxnz>(F20.10))') (c(I,1:maxnz,J),I=1,ny) 100 WRITE(11,*) ']' 101 ENDDO 102 103C write H(ny,nx,ns(2)) 104 DO J = 1,ns(2) 105 WRITE(11,*) ' ' 106 IF (J.LE.9) THEN 107 WRITE(CJ,'(I1)') J 108 ELSEIF ((J.GE.10).AND.(J.LE.99)) THEN 109 WRITE(CJ,'(I2)') J 110 ELSEIF ((J.GE.100).AND.(J.LE.999)) THEN 111 WRITE(CJ,'(I3)') J 112 ENDIF 113 WRITE(11,*) 'H(1:ny,1:nx,'//TRIM(CJ)// ') = [' 114 WRITE(11,*) ' ' 115 WRITE(11,'(<nx>(F20.10))') (H(I,1:nx,J),I=1,ny) 116 WRITE(11,*) ']' 117 ENDDO 118 119C write G(ny,nu,ns(3)) 120 DO J = 1,ns(3) 121 WRITE(11,*) ' ' 122 IF (J.LE.9) THEN 123 WRITE(CJ,'(I1)') J 124 ELSEIF ((J.GE.10).AND.(J.LE.99)) THEN 125 WRITE(CJ,'(I2)') J 126 ELSEIF ((J.GE.100).AND.(J.LE.999)) THEN 127 WRITE(CJ,'(I3)') J 128 ENDIF 129 WRITE(11,*) 'G(1:ny,1:nu,'//TRIM(CJ)// ') = [' 130 WRITE(11,*) ' ' 131 WRITE(11,'(<nu>(F20.10))') (G(I,1:nu,J),I=1,ny) 132 WRITE(11,*) ']' 133 ENDDO 134 135C write a(nx,ns(4)) 136 DO J = 1,ns(4) 137 WRITE(11,*) ' ' 138 IF (J.LE.9) THEN 139 WRITE(CJ,'(I1)') J 140 ELSEIF ((J.GE.10).AND.(J.LE.99)) THEN 141 WRITE(CJ,'(I2)') J 142 ELSEIF ((J.GE.100).AND.(J.LE.999)) THEN 143 WRITE(CJ,'(I3)') J 144 ENDIF 145 WRITE(11,*) 'a(1:nx,'//TRIM(CJ)// ') = [' 146 WRITE(11,*) ' ' 147 WRITE(11,'(<1>(F20.10))') (a(I,J),I=1,nx) 148 WRITE(11,*) ']' 149 ENDDO 150 151C write F(nx,nx,ns(5)) 152 DO J = 1,ns(5) 153 WRITE(11,*) ' ' 154 IF (J.LE.9) THEN 155 WRITE(CJ,'(I1)') J 156 ELSEIF ((J.GE.10).AND.(J.LE.99)) THEN 157 WRITE(CJ,'(I2)') J 158 ELSEIF ((J.GE.100).AND.(J.LE.999)) THEN 159 WRITE(CJ,'(I3)') J 160 ENDIF 161 WRITE(11,*) 'F(1:nx,1:nx,'//TRIM(CJ)// ') = [' 162 WRITE(11,*) ' ' 163 WRITE(11,'(<nx>(F20.10))') (F(I,1:nx,J),I=1,nx) 164 WRITE(11,*) ']' 165 ENDDO 166 167C write R(nx,nu,ns(6)) 168 DO J = 1,ns(6) 169 WRITE(11,*) ' ' 170 IF (J.LE.9) THEN 171 WRITE(CJ,'(I1)') J 172 ELSEIF ((J.GE.10).AND.(J.LE.99)) THEN 173 WRITE(CJ,'(I2)') J 174 ELSEIF ((J.GE.100).AND.(J.LE.999)) THEN 175 WRITE(CJ,'(I3)') J 176 ENDIF 177 WRITE(11,*) 'R(1:nx,1:nu,'//TRIM(CJ)// ') = [' 178 WRITE(11,*) ' ' 179 WRITE(11,'(<nu>(F20.10))') (R(I,1:nu,J),I=1,nx) 180 WRITE(11,*) ']' 181 ENDDO 182 183C Check unstable eigenvalues of F 184 DO J = 1,ns(5) 185 IF (d(2).GT.0) THEN 186 IFAIL=-1 187C CALL F02EBF('N',d(2),F(1:d(2),1:d(2),J),d(2), 188C 1 WR(1:d(2)),WI(1:d(2)),VR,1,VI,1,WORK,4*nx,IFAIL) 189 CALL DGEEV('N','N',d(2),F(1:d(2),1:d(2),J),d(2), 190 1 WR(1:d(2)),WI(1:d(2)),VR,1,VI,1,WORK,4*nx,IFAIL) 191 192 ESTABLE = 0 193 DO I = 1,d(2) 194 W(I) = WR(I)**2+WI(I)**2 195 ESTABLE = ESTABLE + ABS(W(I).GE.1.D0) 196 ENDDO 197 IF (ESTABLE.NE.d(2)) THEN 198 WRITE(11,*) ' ' 199 WRITE(11,*) 'WARNING: the number of unstable eigenvalues for ' 200 WRITE(11,*) 'F(1:nx,1:nx,'//TRIM(CJ)// 'is not equal to d(2) ' 201 WRITE(11,*) 'or non-stationary states are not placed in the' 202 WRITE(11,*) 'first d(2) positions.' 203 ENDIF 204 ENDIF 205 206C Check stable eigenvalues of F 207 IF (nx-d(2).GT.0) THEN 208 IFAIL=-1 209c CALL F02EBF('N',nx-d(2),F(d(2)+1:nx,d(2)+1:nx,J), 210c 1 nx-d(2),WR,WI,VR,1,VI,1,WORK,4*nx,IFAIL) 211 CALL DGEEV('N','N',nx-d(2),F(d(2)+1:nx,d(2)+1:nx,J), 212 # nx-d(2),WR,WI,VR,1,VI,1,WORK,4*nx,IFAIL) 213 214 ESTABLE = 0 215 DO I = 1,nx-d(2) 216 W(I) = WR(I)**2+WI(I)**2 217 ESTABLE = ESTABLE + ABS(W(I).LT.1.D0) 218 ENDDO 219 IF (ESTABLE.NE.(nx-d(2))) THEN 220 WRITE(11,*) ' ' 221 WRITE(11,*) 'WARNING: the number of stable eigenvalues for ' 222 WRITE(11,*) 'F(1:nx,1:nx,'//TRIM(CJ)//'is not equal to nx-d(2)' 223 WRITE(11,*) 'or non-stationary states are not placed in the ' 224 WRITE(11,*) 'first d(2) positions.' 225 ENDIF 226 ENDIF 227 ENDDO 228 229 CLOSE(11) 230 DEALLOCATE(c,H,G,a,F,R) 231 232 RETURN 233 PAUSE 234 END 235 236C Check rank{(HR+G)*(HR+G)'} this check is wrong!!! 237c DO J = 1,ns(2) !H 238c DO JJ = 1,ns(3) !G 239c DO JJJ = 1,ns(6) !R 240 241c DO I =1,ny 242c DO K =1,nu 243c HRG(I,K) = SUM(H(I,1:nx,J)*R(1:nx,K,JJJ))+G(I,K,JJ) 244c ENDDO 245c ENDDO 246 247c DO I =1,ny 248c DO K =1,ny 249c HRGRH(I,K) = SUM(HRG(I,1:nu)*HRG(K,1:nu)) 250c ENDDO 251c ENDDO 252 253c IFAIL = -1 254c CALL F02FAF('N','U',ny,HRGRH,ny,WRY(1:ny),WORK1,64*ny,IFAIL) 255c SRANK = 0 256c DO 10 I=1,ny 257c10 IF (WRY(I).GT.1.D-12) SRANK=SRANK+1 258 259c IF (SRANK.LT.ny) THEN 260c WRITE(11,*) ' ' 261c WRITE(11,*) 'WARNING: the rank of the system computed looking ' 262c WRITE(11,*) 'at rank{(HR+G)*transpose(HR+G)} is less than ny ' 263c ENDIF 264 265c ENDDO 266c ENDDO 267c ENDDO 268 269