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).