27.0 OPTCONTROL Command The OPTCONTROL command branches to the Gregory Chow optimal control routines. If targets for endogenous and control variables are supplied, the OPTCONTROL command calculates the optimal levels of the control variables. The model is: y(t) = A(t) y(t-1) + C(t) x(t) + b(t) + u(t) (1) where y(t) is a vector of p (=NY + NX) state varibles, x(t) is a vector of NX control variables and A(t), C(t) and b(t) are known matrices. The optimal control problem is to minimize the loss function W = SUM((y(t) - A(t)'(K(t)(y(t) - A(t)). The original program control language has been maintained and is discussed in Chow (1981 Chapter 5). Prior to running this option the user must prepare a SUBROUTINE MODELS. B34S will dynamically link with this routine. For further details on the limitations and advantages of the B34S dynamic link facility, see the NONLIN help document. A sample job, part of which is shown below, is contained in BEC4346.#B34S.TEST.LIB(CHOW2). Prior to running the OPTCONTROL command, it is suggested that this job be run. Complete references are: The currect PC, and UNIX versions of B34S do not support dynamic linking. The OPTCONTROL command can be used if B34S is relinked. B34S can be relinked if the user has access to the B34S libraries and the Lahey FORTRAN compiler. - Chow, Gregory, "Econometric Analysis by Control Methods", Wiley 1981, especially chapter 5. - Chow, Gregory, "Analysis and Control of Dynamic Economic Systems," Wiley, 1975. The OPTCONTROL command also will run using the B34S command language. This option is discussed also. Form of the OPTCONTROL command using Chow (1981) input form. B34SEXEC OPTCONTROL options parameters$ PGMCARDS$ optimal control program cards here B34SRETURN$ B34SEEND$ OPTCONTROL sentence parameters. IJUMP Calls the subroutine MODELS which must be linked into B34S. The IJUMP option is the only way that the OPTCONTROL command can be run on PC, and UNIX operating systems. INPUT=n1 Sets input unit. Default = 5. If INPUT set NE 5, different input control card set can be accessed. LIST=n2 Sets output unit. Default = 6. INOUT2=n3 Sets unit to store solution. Default = 11. Used if NROUND=2 on &OPT card. INOUT3=n4 Sets unit to store partial derivatives. Default=12. Used if NROUND=2. INOUT4=n5 Sets unit to store big and little g (G(t), g(t)). Default=13 INOUT5=n6 Sets unit to store A(t) + C(t)G(t) and V(t) for each period. Default = 14. script=name execution commands. Usually this is not needed to be changed model =mname Sets model name. Default = models. Note: INOUT2 - INOUT5 use DP unformatted files. Coding user SUBROUTINE MODELS. Detailed instructions for coding user SUBROUTINE MODELS is contained in Chow (1981) Chapter 5. In that reference two routines are described. In the B34S implementation of the Chow program, there is only ONE user SUBROUTINE. User routine MODELS has arguments SUBROUTINE MODELS(II,NY,NX,NW,D,Y,YL,XL,X,Z,IBAD) II If II = 0 routine branches to NS structural equations If II NE 0 routine branches to NID identities of the lagged variables. NY = # of structural equations (NS) + # of identities. NX = # of control variables. NW = # of uncontrolled exogenous variables. D = array containing new values of the endogenous variables. Y = vector of endogenous variables. YL = vector of lagged endogenous variables. XL = vector of lagged control variables. X = vector of control variables. Z = vector of exogenous variables. If there are no exogenous variables, Z and NW have to be specified. IBAD= 0 for normal return. = 1 for abnormal return. The user sets IBAD=1 if some endogenous variable (i. e. D(k) ) gets outside the acceptable range. This problem will occure if the trial solution is too poor to continue with the problem. Usage note: All real variables in SUBROUTINE MODELS are REAL*8. A detailed discussion of this routine is contain below. *************************************** Sample job. /$ /$ Uses I/O and User Fortran /$ b34sexec options open('test.f') disp=unknown unit(77); b34srun; b34sexec options rewind(77); b34srun; b34sexec options copyf(4,77); pgmcards; c c Load Interface files which hard wire models name c include 'c:\b34slm\optrun.f' subroutine models(ii,ny,nx,nw,d,y,yl,xl,x,z,ibad) c c built in test problem c implicit real*8(a-h,o-z) dimension y(ny),x(nx),yl(ny),xl(nx),d(ny),z(nw) c c this combines both routines in one ::::::::::::::::: c ibad=0 c if set ne 0 will mean non standard return c if(ii.gt.0)go to 999 c c c imsl i/o debug code - note that d and y can be c in same location c c call dwrrrn('d at 1 ',1,ny,d, 1,0) c call dwrrrn('y at 1 ',1,ny,y, 1,0) c call dwrrrn('yl at 1 ',1,ny,yl,1,0) c call dwrrrn('xl at 1 ',1,nx,xl,1,0) c call dwrrrn('x at 1 ',1,nx,x, 1,0) c call dwrrrn('z at 1 ',1,nw,z, 1,0) c d(22) = .0512d0*(y(9)+x(3)) d(21) = .248d0*(y(14)-y(20)-y(3)) + .2695d0*yl(15)*(yl(14)- 1 yl(20)-yl(3))/y(15) 2 + .4497d0*y(4)-5.7416d0 d(20) = .4497d0*y(4)+2.7085d0 d(19) = .1549d0*y(6)+.131*x(1)-6.9076d0 d(18) = .0924d0*y(13)-1.3607d0 d(17) = yl(17) + y(3) d(16) = yl(16) + y(2) - y(5) d(15) = 1.062d0*y(8)*y(7)/(y(6)+x(1)) d(14) = y(13)-y(18)-y(5)-y(6)-x(1)-y(9)-x(3) d(13) = y(1) + y(2) + x(2) d(12) = .26d0*y(6)-2.55d0-.26d0*(y(15)-yl(15)) + .61d0*yl(12) d(11) = .14d0*(y(6) + x(1) - y(19)+y(14)-y(21)-y(3)+y(9) +x(3) 1 - y(22))+54.083939d0 d(10) = 1.39d0*y(15)+32.d0 d(9) = .054d0*(y(6)+x(1)-y(19)+y(14)-y(21)-y(3) ) 1 + .012d0*z(1)*y(10)/y(15) d(8) = yl(8) +4.11d0-.74d0*(z(3)-y(7)-z(4)-z(5)) + .52d0*( yl(15) 1 -yl(23))+ .54d0*z(6) d(7) = x(4)-(z(4)+z(5))/1.062d0 +(26.08d0+y(13)-x(1)-.08d0*y(16) 1 -.08d0*yl(16) -2.05*z(6) )/ 2.304537d0 d(6) = -1.4d+00 + .24d0*(y(13)-x(1))+ .24d0*(yl(13) -xl(1)) 1 + .29d0*z(6) d(5) = 7.25d0+.05d0*(y(16)+yl(16)) +.044d0*(y(13)-x(1)) d(4) = -7.6d0 +.68*y(14) d(3) = -3.53d0+.72d0*(y(4)-y(20))-.028d0*yl(17) d(2) = -16.71d0+.78d0*(yl(14)-yl(21)+yl(9)+xl(3)-yl(22) 1 + yl(5)) - .073d0*yl(16) +.14d0*yl(12) d(1) = -22.26 + .55*(y(6)+x(1)-y(19)) + .41d0*(y(14)-y(21)-y(3)) 1 + .34d0*(y(9)+x(3)-y(22))+.26d0*yl(1) +.072d0*yl(11) +.26d0* 2 z(2) go to 99 c 999 continue c c statements for modeli are next :::::::::::::::::::::::: c d(23)=yl(15) c 99 continue c c call dwrrrn('d at 2 ',1,ny,d, 1,0) c call dwrrrn('y at 2 ',1,ny,y, 1,0) c call dwrrrn('yl at 2 ',1,ny,yl,1,0) c call dwrrrn('xl at 2 ',1,nx,xl,1,0) c call dwrrrn('x at 2 ',1,nx,x, 1,0) c call dwrrrn('z at 2 ',1,nw,z, 1,0) c return end b34sreturn; b34srun; b34sexec options close(77); b34srun; b34sexec options dodos('lf95 test.f') dounix('f77 test.f'); b34srun; /$ ijump runs the internal test setup -- not used usually B34SEXEC OPTCONTROL script='test > testbug' model='models' /$ /$ ijump $ OPTIONS NY=23,NX=4,NS=22,NW=6,NPD=5 OSUP=3,DAMP=1.0,ITERL2=150 EPS1=0.0001,EPS2=0.0001,UGZ=T NROUND=1,ITERL1=5,GAMMA=6 $ IDENTVEC=15 COMMENT('SETS LAG OF VAR 15')$ VARBLVEC=(1 2 3 4 5 6 7 8 9 10 0 0 13 14 15 16 0 18 19 20 21 22 23 0 25 0 27 0 0 30 31 0 33 34 35 36 37 38 39 0 0 42 43 44 45 46 0 48 0 50 51 52 53) COMMENT=('SETS PARTIALS')$ V = (.9718 8.6789 29.4632 .49 .4886 1.0661 1.0126 23.6293 .9120 .0003 1.6669 246.9927 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.) COMMENT=('VAR-COV')$ KCAP=(0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0) COMMENT=('Y7, Y15, Y26 Y27')$ Z0 =(111.4 24.3 1.9 16.51 19.35 78.65 56.0 326.2 7.3 303.0 95.2 38.10 172.0 35.17 202.4 41.5 .19 14.51 8.63 10.14 13.72 .38 197.5 15.12 0.0 .1187 9.393) COMMANT=('INITIAL TARGET VALUES')$ ZRATE=(1.0 1.0 1.0 1.0 1.0 1.0 1.02 1.035 1.0 1.0 1.0 1.0 1.05 1.05 1.01 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.03 ) COMMENT=('TARGET GROWTH RATE')$ Y0 =(111.4 24.3 1.9 16.51 19.35 78.65 56.0 326.2 7.3 303.0 95.2 38.1 172.0 35.17 202.4 41.5 .19 14.51 8.63 10.14 13.72 .38 197.5 15.12 0.0 .1187 9.393) COMMENT=('INITIAL Y-X VECTOR')$ W =(171.859985 159.6 67.63 6.53 4.25 24.00 171.859985 162.114990 68.629990 6.649999 4.222 25.000000 171.859985 164.703979 69.593994 6.767999 4.203 26.000000 171.859985 167.292984 70.557983 6.885999 4.187 27.000000 171.859985 169.969672 71.545795 7.003061 4.173 28.000000) COMMENT=('EXOGENOUS VARIABLES FOR ALL PERIODS')$ X =(15.700000 33.500000 0.117000 9.709999 16.170996 34.504991 0.120510 10.001296 16.656111 35.540109 0.124125 10.301326 17.155779 36.606281 0.127849 10.610357 17.670435 37.704432 0.131684 10.928657) COMMENT=('POLICY GUESS')$ Y =(111.4 24.3 1.9 16.51 19.35 78.65 56.0 326.2 7.3 303.0 95.2 38.1 172.0 35.17 202.4 41.5 .19 14.51 8.63 10.14 13.72 .38 197.5) COMMENT=('INITIAL Y')$ GAMVAR=(7 13 14 15 26 27) COMMENT=('ROW COLS OF GAMMA TO PRINT')$ B34SRETURN$ B34SEEND$ Listing of optrun.f Interface File C Last change: HHS 10 May 2002 8:46 pm c Driver for OPTCONTROL -- Use include to place at c top of FORTRAN of users model subroutine c include 'c:\b34slm\optrun.f' c implicit real*8(a-h,o-z) parameter(maxny=500,maxnx=500,maxnw=500) logical xx dimension d(maxny),y(maxny),yl(maxny), * xl(maxnx),x(maxnx),z(maxnw) c inquire(file='models',exist=xx) c if(.not.xx)then write(6,*)'Model file not found' stop endif c open(file='models',unit=77,err=991) rewind(unit=77) read(unit=77,fmt='(6i10)',end=992)ny,nx,nw,ii,ibad,isame c if(ny.gt.maxny.or.ny.gt.maxnx.or.nw.gt.maxnw)then write(6,*)'Limits in optrun exceeded' write(6,*)'ny = ',ny,' maxny = ',maxny write(6,*)'nx = ',nx,' maxnx = ',maxnx write(6,*)'nw = ',nw,' maxnw = ',maxnw stop endif c read(unit=77,fmt='(3e25.16)',end=993,err=993)(d(i),i=1,ny), * ( y(i),i=1,ny), * (yl(i),i=1,ny), * (xl(i),i=1,nx), * ( x(i),i=1,nx), * ( z(i),i=1,nw) c c = note that y ius updated as is d c if(isame.eq.0)then call models(ii,ny,nx,nw,d,y,yL,xL,x,z,ibad) else call dcopy(ny,y,1,d,1) call models(ii,ny,nx,nw,d,d,yL,xL,x,z,ibad) endif c rewind(unit=77) write(unit=77,fmt='(6i10)', err=993)ny,nx,nw,ii,ibad,isame if(isame.eq.0)then write(unit=77,fmt='(3e25.16)',err=993)( d(i),i=1,ny), * ( y(i),i=1,ny), *(yl(i),i=1,ny),(xl(i),i=1,nx),(x(i),i=1,nx),(z(i),i=1,nw) else write(unit=77,fmt='(3e25.16)',err=993)( d(i),i=1,ny), * ( d(i),i=1,ny), *(yl(i),i=1,ny),(xl(i),i=1,nx),(x(i),i=1,nx),(z(i),i=1,nw) call dcopy(ny,d,1,y,1) endif c rewind(unit=77) c go to 999 991 ibad=1 go to 998 992 ibad=2 go to 998 993 ibad=3 998 write(6,*)'OPTRUN.F Error # ',ibad go to 999 999 close(unit=77) stop end ********************************************************************* Coding of User Supplied SUBROUTINE MODELS and discussion of B34S input options Problem of second and higher order lags Let the model be written as a system of simultaneous structural equations. y(t) = f(y(t),y(t-1),x(t),x(t-1),w(t)) + u(t), (3) where y(t) = vector of ny endogenous variables x(t) = vector of nx control variables w(t) = vector of nw exogenous variables not subject to control NPD = number of time periods in the planning horizon u(t) = vector of random error terms Because the program does not accept lagged endogenous variables of order higher than the first, the user should eliminate the variables with lags of two or more periods by introducing identities. For example, let there be 50 endogenous variables in the model to begin with. The number of simultaneous equations NS is 50. If the model consists of y6(t-2), an identity y51(t) = y6(t-1) should be introduced. This identity permits the user to write y6(t-2) as y51(t-1) and get rid of the second-order lag. If y6(t-3) is also present, another identity y52(t) = y51(t-1) can be used, permitting the user to write y6(t-3) as y51(t-1). Let 40 additional identities of this kind be required in our example to eliminate all endogenous variables with lags of two or more periods. There will then be 90 endogenous variables in the model, to be included in the vector y(t). If the model consists of control variables lagged two or more periods, more identities and endogenous variables will be required. Let there be NX = 4 control variables in our example, x1(t),...,x4(t). To eliminate the variable x1(t-2), introduce the identity y91(t(= x1(t-1) and write x1(t-2) as y91(t-1). Similarly, to eliminate x1(t-3) introduce the identity y92(t) = y91(t-3) and write x1(t-3) as y92(t-1). If 10 identities of this type are required, there will be all together 100 variables in the vector y(t). NY, the number of endogenous variables in y(t), will be set equal to 100. The computer program will automatically make up a vector consisting of the 104 elements of y(t) and x(t). This augmented vector will serve as the argument in the welfare function, its last four elements in our example being included to serve the possible need to penalize variations in the instruments or control variables. The expanded model takes the form y1(t) = f(1)(y(t),y(t-1),x(t-1),w(t)) + u1(t) . . . NS structural . . . equations . . . y50(t)= f(ns)(y(t),y(t-1),x(t),x(t-1),w(t)) + u50(t) y51(t) = y6(t-1) . . NID identities . . of lagged values . . y(100,t) = y99(t-1) y(101,t) = x1(t) . . NX identities of . . control variables . . y(104,t) = x(4,t) (4) There are NS simultaneous structural equations, NID identities of lagged values, and NX identities for the control variables; altogether the augmented y vector has P elements, with P = NY + NX, and NY = NS + NID. u(t) is a random vector with mean zero and covariance matrix V, to be supplied by the user. The program calculates a linear approximation of the model explaining the augmented y vector shown in (2) y(t) = A(t)y(t-1) + C1x(t) + b(t), (5) from which it derives a feedback control equation x(t) = G(t)y(t-1) + g(t) (6) so as to minimize the expectation of the welfare loss NPD E0W = E0 SUM (y(t) - z(t))'K(t)(y(t) - z(t)), (7) t=1 where NPD is the number of periods, zt is a vector of targets, and K(t) = K*EXKCAP(t) is the weighting matrix. NPD, z(t), KCAP, and K are to be specified by the user. The model in the form of equation 21.4 is coded in two sections of the user SUBROUTINE MODELS. If II = 0, the routine will branch to the N simultaneous structural equations. If II NE 0 the routine will branch to the NID identities of the lagged values. The last NX identities of the control variables are not coded by the user, but are "remembered" by the program. SUBROUTINE MODELS must be coded in double precision. For an example sample program as start of this section. include 'c:\b34slm\optrun.f' subroutine models(ii,ny,nx,nw,d,y,yl,xl,x,z,ibad) implicit real*8(a-h,0-z) dimension y(ny),x(nx),yl(ny),xl(nx),d(ny),z(nw) ibad=0 if(ii.gt.0)go to 999 c structural equations code here . . . go to 99 999 continue c identitiy code here . . 99 return end The include file places the inmterface file at the start of the fortran source. While a dynamic link such as a dll or unix dynamic library is faster, the I/O implementation here is 100% portable across operating systems. Discussion of arguments II = 0 for structural equations branch of models. = 1 for identities branch of models. NY = # of structual equations (NS) + # of identities NX = # of control variables NW = # of uncontrollable exogenous variables D = the array that contains the new values of the endogenous variables resulting from the calculations. Warning: as D is built, y is changing in many calls. For example d(1)=.3d+8 * y(6) d(2)=.3d+2 * y(1) in general is not the same as d(2)=.3d+2 * y(1) d(1)=.3d+8 * y(6) Y = the vector of endogenous variables YL = the vector of lagged endogenous variables. X = the vector of control variables. XL = the vector of lagged control variables. Z = the vector of exogenous variables; if there are no exogenous variables in the model, Z and NW should nevertheless be coded as dummy parameters. IBAD = 0 for normal return = 1 for abnormal return. The USER sets IBAD = 1, if for example, some endogenous variable (i.e., D(k)) gets outside the acceptable range. The meaning would be that the trial solution is too poor to continue with the problem. ********************************************************************* General Discussion of input into OPTCONTROL using B34S input Form of the command. B34SEXEC OPTCONTROL options parameters$ OPTIONS options parameters $ IDENTVEC=( ) $ VARBLVEC=( ) $ V =( ) $ KCAP =( ) $ Z0 =( ) $ ZRATE =( ) $ Y0 =( ) $ W =( ) $ X =( ) $ Y =( ) $ GAMVAR =( ) $ B34SEEND$ Required sentences: OPTIONS, IDENTVEC, V, KCAP, Y0, X Except for the OPTIONS sentence, all other sentences have the optional COMMENT parameter. For futher discussion of the setup, see Chow (1981) from which this section has been adapted. Parameters on the OPTIONS sentence. NS = n1 number of structural equations in the model; default = 1. NY = n2 number of endogenous variables of the vector y(t) (5); that is, NY = NS + the number of identities of lagged values; default = 1. NX = n3 number of control variables in the vector xt in (5); default = 1. NW = n4 number of uncontrollable exogenous variables in the vector w(t) in (5); default = 0. NPD= n5 number of time periods for the plan; default = 1. OSUP = 0 for full print-out of A(t),C(t),b(t),G(t),H(t),g(t),h(t) at every time period. [These symbols are defined in Chow (1975) and in Chow (1981) Chapter 1. = 1 for print-out of G(t),H(t),g(t),h(t) only at every time period. = 2 for print-out of G(1),H(1),g(1),h(1) at time period 1 only; all other time periods have no print out. = 3 same as in OSUP = 2 except that the program will also omit printing H1 and h1. This is the default setting GAUSS= 0 will omit the printing of the Gauss-Siedel solution for y(t) used for the first linearization of the model. This is the default setting. = 1 will print the above Gauss-Siedel solution. = 2 will print and save the Gauss-Siedel solution and the program will terminate at that point. GAMMA= 0 will omit the calculation of the covariance matrices gamma(t) (default value). = 1 when the optimal solution y(t) converges, the program will caluculate the covariance matrices t of the variables y(t). The program will also compute t when the iteration limit ITERLI is reached. The program will not print the entire gama(t) but will print only those rows of gama(t) that correspond to the variables targetted in the weighting matrix K. = n where 2 LE n LE (NS + NX). As in GAMMA = 1 the program will calculate gamma(t) when the optimal solution y(t) converges or iteration limit is reached; the print-out gamma(t) is controlled by the input vector GAMVAR, which is supplied by the user; n indicates the number of entries in GAMVAR. PLOT = 0 will plot the values of the means of y(t) against target values z(t) only for those variables that have nonzero diagonal weights in the K matrix; the means of the control variables x(t) are also plotted (default value). = 1 will plot the means of all the variables of the optimal solution y(t) as well as the control variables x(t) against their targets z(t). UGZ = T if the targets z(t) are to grow at constant percentage rates; in this case the user will supply the initial z0, the Z0 vector, and the rate of growth, the ZRATE vector, and the program will compute zt = Z0*(ZRATE)t. = F if the targets z(t) for all periods are to be supplied by the user in the input matrix Z (default value). OFDIAV= T if the covariance matrix V of the random vector u(t) in (3) have nonzero off-diagonal elements; in this case the user must supply the lower triangle of the symmetric V. = F if the off-diagonal elements of V are zero; only the diagonal elements of V are required for input (default value). OFDIAG= T if the K matrix in the loss function have nonzero off- diagonal elements; in this case the user must supply the lower traingle of the symmetric K matrix. = F if the off-diagonal elements of K are zero; only the diagonal elements of K are required for input (default value). EXKCAP = discount factor for modifying the K matrix for each time period t according to the formula K(t )= K*(EXKCAP)**t, where K is supplied by the user; default = 1.0. NROUND = 1 if the tentative path for yt required in the first iteration is to be computed by the Gauss-Siedel method; this marks the first job run in the iterative process to calculate the optimal path yt which might require more than one job run (default value). = 2 (or 3,4, etc.) marks the second (or third, fourth, etc.) job run in the iterative process to find the optimal path yt; the solution for yt and xt from the previous job will be read from an output device named INOUT2, to serve as the tentative path for the first iteration of the current run. All other input data remain the same as in the previous job. DAMP =r1 A factor used to dampen the changes in successive iterations of the Gauss-Siedel solution. y(t)(k+1) = y(t)(k) + DAMP*(y'(t)(k+1) - y(t)(k)) where y'(t)(k+1) = solution obtained at iteration k +1, and y(t)(k+1) = solution actually used to compute y(t) in the (k + 1)th iteration. DAMP may be made small if solution tends to oscillate from iteration to iteration; default = 1.0. DAMPX = a factor used to dampen the changes in successive iterations of the control variables xt. x(t)(k+1) = x(t)(k) + DAMPX*(x'(t)(k+1) - x(t)(k)) where x'(t)(k+1) = solution obtained at iteration k+1 and x(t)(k+1) = solution actually used to compute the optimal path y(t) for the (k + 1)th iteration; default = 1.0. ITERL1 = the maximum number of times the model will be linearized for calculating optimal control path; default = 1. The method starts with a guess for x1,...,xnpd, which, by the use of the Gauss-Siedel method, if NROUND = 1, implies a solution for y1,...,ynpd. Around this tentative path the model is linearized and an optimal control path is obtained for the linearized model. This optimal path becomes the initial guess of x1,...,xnpd in the next linearization. ITERL1 refers to the maximum number of times the model will be linearized. ITERL2 = the maximum number of iterations allowed for solving the model by Gauss-Siedel; default = 10. EPS1 = the convergence criterion for optimal policy; default = 0.001. | (k) (k-1) | | y(i) - y(i) | when |_____________________| LE EPS1, | (k-1) | | y(i) | for each i = 1,...,NS that corresponds to a nonzero diagonal element of the K matrix, k being the iteration count, the program terminates. The solution yit will be plotted against their targets zit. If the convergence criterion is not met another iteration will be performed. EPS2 = convergence criterion for the solution of the non- linear model by the Gauss-Siedel method; default = 0.001, Gauss-Siedel solves the system of equations (1) until | (k+1) (k) | | y(i) - y(i) | |_________________ _| LE EPS2, | (k) | | y(i) | for i = 1,...,NY,k being the iteration count. FDFRAC = parameter for computing step sizes in evaluating the first derivatives; the step size for the variable yi(t) is dyi = max(|FDFRAC.yi(t)|,FDMIN); default = 0.001. FDMIN = minimum step size allowed in the above formula; default = 0.001. Input Matrices are read with named sentences. The COMMENT parameter, which allows passing of a 40 character comment, is supported on all of the below listed sentences All data matrices except V and KCAP are to be read in row-wise., IDENTVEC, VARBLVEC, and GAMVAR are the only vectors whose elements are integers; all other vectors have floating-point numbers as elements. The dimension P is defined as P = NY + NX,. NID = NY - NS. IDENTVEC, V KCAP and X are required sentences. SENTENCE Dimension Description IDENTVEC (NID) A vector in which elements are the variables that appear in the identities of lagged values in the model; its construction will be described in the next section. VARBLVEC (NS + P + NX) A vector that identifies those variables with respect to which the first partial derivatives are computed by the program; its con- struction will be described in the next section. Note: If a derivative is known to be zero, it can be specified on VARBLVEC and time can be saved. V (NS) if OFDIAV = F Variance-covariance matrix of ut [NS(NS + 1)/2] in (3); if OFDIAV = F, only the if OFDIAV = T diagonal of V is entered in the input; if OFDIAV = T, the lower triangle of V is entered by col- umns, starting with the left-most column of ns elements each suc- ceeding column has one less element than the previous one. KCAP (P) if OFDIAG = F The weighting matrix K in the wel- [P(P + 1)/2] fare function; if OFDIAG = F, only if OFDIAG = T the diagonal of K is entered as input; if OFDIAG = T, the lower triangle is entered by columns, as described for V above. Z0 (P) Initial values of targets for yt and xt; required only if UGZ = T. ZRATE (P) Growth rate for targets which the program will compute based on Z0 and ZRATE; required only if UGZ = T Z (NPD,P) Targets to be read in row-wise for all periods; required only if UGZ = Y0 (P) Initial values for the vectors y(t) and x(t) in period t = 0. W (NPD,NW) Exogenous variables for all periods, to be read in row-wise; required only if NW > 0. X (NPD,NX) A tentative policy for all periods, to be read in row-wise. This policy will be used to compute the tentative path y(t) by Gauss-Siedel. Y (NY) A trial solution of the y vector for the first period, to be used in the first iteration of the Gauss-Siedel solution for period 1; this vector may take the same value as Y0. GAMVAR (N) A vector whose elements are the row/column numbers of t to be printed; required only if GAMMA = n, where 2 LE N LE (NS + NX) and is the dimension of GAMVAR. Construction of IDENTVEC and VARBLVEC The use of IDENTVEC and VARBLVEC is discussed in a Chow (1981). Before constructing IDENTVEC, the user must first eliminate all higher- order lags in the nonlinear model by adding identities of lagged values to the model, as discussed in the section discussing problems of higher order lags. In this example the model has 50 simultaneous structural equations, 50 identities of lagged values, and four control variables; that is, NS = 50, NID = 50, NY = NS + NID = 100, and NX = 4, and the identities of lagged values are: y51(t) = y6(t-1) 6 y52(t) = y51(t-1) 51 y53(t) = y52(t-1) 52 . . . . . . y91(t) = x1(t-1) 101 y92(t) = y91(t-1) 91 y93(t) = x3(t-1) 103 . . . . . . y100(t) =y99(t-1) 99 The numbers listed on the right are the variable numbers that will appear in IDENTVEC: for yk(t-1) in an identity, the number listed is simply k; for xk(t-1), the number listed is NY + k. The IDENTVEC for this model will look like this: IDENTVEC=(6 51 52,................,101 91 103..99)$ The dimension of IDENTVEC is NID, the number of identities of lagged values. The construction of VARBLVEC is a bit more involved, but the user should not have too much trouble if the steps listed are followed: 1. List the integers from 1 to NS + P + NS. In our example from the numbers would be from 1 to 158. 2. Above this list of numbers write the following variables in order, with one variable corresponding to one number: y1(t),...,yNS(t) y1(t-1),...yNY(t-1),x1(t-1),...,xNX(t-1) x1(t),...,xNX(t) In our example, we would have: y1(t) y2(t) ... y50(t) 1 2 ... 50 y1(t-1) y2(t-1) ... y100(t-1), x1(t-1),...,x4(t-1) 51 52 ... 150 151 ,..., 154 x1(t) ... x4(t) 155 ... 158 3. Look through the right-hand side of each structural equation of the model. Whenever one of the variables listed in step 2 appears, circle the number underneath the variable in the list. In our example, we might obtain 1 2 3 4 ... 50 51 52 53 ... . . . 155 156 157 158 Any number in the list not underlined means the variable above it appears nowhere in the right-hand side of the structural equations; any variable of the list in step 2 that appears somewhere in the right-hand side of the structural equations should have its corresponding number circled. In our example, y3(t), y50(t),y1(t-1),x2(t) are some of the variables not in the structural model. 4. VARBLVEC is simply the list of integers from step 2, with the underlined numbers replaced by 0; thus VARBLVEC has NS + P + NX entries. In our example, VARBLVEC would look like this: VARBLVEC=(1 2 0 4 ... ......... 0 52 53 .....155 0 157 158)$ Note: It is vitally important that IDENTVEC and VARBLVEC be con- structed correctly because the program depends on them to calculate the optimal solution yt. If the econometric model is not very large or the saving of computer space and time is not essential, the user may simply fill in the vector VARBLVEC with all nonzero elements. If the VARBLVEC sentence is not supplied, it will default to all non zero elements. Sample input file showing b34s commands /$ /$ uses new control language /$ b34sexec optcontrol mod=ch13mod $ options ny=23,nx=4,ns=22,nw=6,npd=5 osup=3,damp=1.0,iterl2=150 eps1=0.0001,eps2=0.0001,ugz=t nround=1,iterl1=5,gamma=6 $ identvec=15 comment('sets lag of var 15')$ varblvec=(1 2 3 4 5 6 7 8 9 10 0 0 13 14 15 16 0 18 19 20 21 22 23 0 25 0 27 0 0 30 31 0 33 34 35 36 37 38 39 0 0 42 43 44 45 46 0 48 0 50 51 52 53) comment=('sets partials')$ v = (.9718 8.6789 29.4632 .49 .4886 1.0661 1.0126 23.6293 .9120 .0003 1.6669 246.9927 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.) comment=('var-cov')$ kcap=(0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0) comment=('y7, y15, y26 y27')$ z0 =(111.4 24.3 1.9 16.51 19.35 78.65 56.0 326.2 7.3 303.0 95.2 38.10 172.0 35.17 202.4 41.5 .19 14.51 8.63 10.14 13.72 .38 197.5 15.12 0.0 .1187 9.393) commant=('initial target values')$ zrate=(1.0 1.0 1.0 1.0 1.0 1.0 1.02 1.035 1.0 1.0 1.0 1.0 1.05 1.05 1.01 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.03 ) comment=('target growth rate')$ y0 =(111.4 24.3 1.9 16.51 19.35 78.65 56.0 326.2 7.3 303.0 95.2 38.1 172.0 35.17 202.4 41.5 .19 14.51 8.63 10.14 13.72 .38 197.5 15.12 0.0 .1187 9.393) comment=('initial y-x vector')$ w =(171.859985 159.6 67.63 6.53 4.25 24.00 171.859985 162.114990 68.629990 6.649999 4.222 25.000000 171.859985 164.703979 69.593994 6.767999 4.203 26.000000 171.859985 167.292984 70.557983 6.885999 4.187 27.000000 171.859985 169.969672 71.545795 7.003061 4.173 28.000000) comment=('exogenous variables for all periods')$ x =(15.700000 33.500000 0.117000 9.709999 16.170996 34.504991 0.120510 10.001296 16.656111 35.540109 0.124125 10.301326 17.155779 36.606281 0.127849 10.610357 17.670435 37.704432 0.131684 10.928657) comment=('policy guess')$ y =(111.4 24.3 1.9 16.51 19.35 78.65 56.0 326.2 7.3 303.0 95.2 38.1 172.0 35.17 202.4 41.5 .19 14.51 8.63 10.14 13.72 .38 197.5 15.12) comment=('initial y')$ gamvar=(7 13 14 15 26 27) comment=('row cols of gamma to print')$ b34seend$ OUTPUT AND ERROR MESSAGES The program prints the following on unit LIST for each job: 1. A list of all options, including the default values if they have not been overridden. 2. All the input matrices in the order they are read in. 3. If GAUSS = 1 or 2, the values of the Gauss-Siedel solution for yt used for the first linearization of the model. 4. The A,b,C,G,g,h,H matrices, starting with the last period. Their frequency of print-out from period to period depends on the value of the OSUP parameter. 5. The values of the optimal solution yt. 6. Total welfare cost and its deterministic and stochastic components. 7. If the convergence criterion EPS1 is met before the number of iterations becomes greater than ITERL1, a plot of the values of the optimal solution y(t) against the target values z(t). Whether all of y(t) variables or only those that have nonzero weights in the diagonal of the K matrix are plotted depends on the value of the PLOT parameter. The control variables x(t) are always plotted against their targets. 8. TABL1 and TABL3. The only useful information these two tables give to the user is the number of entries in TABL3, which should be used as the value for p*. 9. The t matrices if GAMMA > 0. If GAMMA = 1, the rows of t correpsonding to those elements targetted in K are printed; if GAMMA > 2, the input vector GAMVAR controls the printing of gamma(t) The following is a list of messages that may be printed: "PROGRAM TERMINATING BECAUSE GAUSS-SIEDEL FAILED TO CONVERGE FOR VARIABLE I, PERIOD T, IN K ITERATIONS" -In calculating the Gauss-Siedel solution for yt, if any variable yi,t does not converge after ITERL2 iterations, this message appears and the progam terminates with a condition code of 29. "PROGRAM TERMINATING BECAUSE OPTIMAL SOLUTION FAILED TO CONVERGE IN K ITERATIONS, ROUND NO.N"--If the optimal solution yt does not converge after ITERL1 iterations, this message appears and the program terminates with a condition code of 29. In addition, the optimal solution yt and xt are saved on the output device assigned by INOUT2 of the main program. "STORAGE ALLOCATION TO THIS POINT = XXXX WORDS OUT OF YYYY (ZZK UNUSED)"--This message appears throughout the printed output, the user may reduce the value of NDIM in the main program accordingly. "***INSUFFICIENT MEMORY AREA TO ALLOCATE WORK ARRAYS, XXXX WORDS ALLOCATED TO STORAGE BY MAIN PROGRAM, YYYY WORDS REQUIRED SO FAR, INCLUDING THE FOLLOWING ZZZZ WORDS CURRENTLY REQUESTED AT THE POINT INDICATED BELOW"--A trace table follows this message showing the subroutine that led to it. The program will continue and allocates more space when required but will eventually terminate after messages like the one following are printed. "STORAGE ALLOCATION TO THIS POINT = XXXX WORDS OUT OF YYYY (ZZK REQUIRED)"--When several messages of this type appear before the program terminates, the user should set B34S work away to 2 times the greatest value of XXXX in these messages for the next job (see C.17). "LOOKING FOR 'AAAA' DATA, FOUND THE FOLLOWING: ..." "LOOKING FOR 'AAAA' DATA, FOUND END-OF-FILE" "END-OF-FILE WHILE READING DATA FOR 'AAAA'"--These messages appear when the user fails to supply enough data for the matrix named AAAA or when the entire matrix has been omitted from the input deck. the program terminates with a condition code of 130, 131, or 132. Final Comments on SS Option The SS option is a very powerful program which the developer of B34S was fortunate to obtain from Professor Gregory Chow. The SS option code has been heavily modified by H. H. Stokes, who is responsible for the accuracy of the program. Users interested in utilizing this option should consult Chow (1975) and Chow (1981).