@ STANDARD PROGRAMME FOR COMPUTATION OF POLICY FUNCTIONS
  USING MATRIX RICATTI ITERATIONS @

@ ------------------------------------------------------------ @
@ THIS IS AN EXAMPLE FOR THE SIMPLE STOCHASTIC GROWTH MODEL    @
@ ------------------------------------------------------------ @

TX=HSEC;

@ ------------------------------------------------------------ @
@                INITIALIZING: SPECIFICATION OF PARAMETERS
                 AND STEADY-STATE RELATIONSHIPS                @
@ ------------------------------------------------------------ @

TT=1;                   @ TIME ENDOWMENT @
NBAR=0.3;               @ SHARE OF TIME USED ON THE MARKET @
RR=0.04/4;              @ REAL INTEREST RATE @
SN=0.64;                @ LABOR SHARE OF INCOME @
SK=1-SN;                @ CAPITAL SHARE OF INCOME @
GAMMAX=1.00;           @ GROWTH RATE @
BSTAR=GAMMAX/(1+RR);    @ DISCOUNT FACTOR @
SIGMA=2;                @ RELATIVE RISK AVERSION PARAMETER @
DELTA=0.025;            @ DEPRECIATION RATE @
KYRATIO=BSTAR*(1-SN)/(GAMMAX-BSTAR*(1-DELTA));
                        @ CAPITAL-OUTPUT RATIO @
SG=0.25;                @ GOVERNMENT SPENDING SHARE OF OUTPUT @
IKRATIO=GAMMAX-(1-DELTA);@ INVESTMENT-CAPITAL RATIO @
SI=IKRATIO*KYRATIO;     @ INVESTMENT SHARE OF OUTPUT @
SC=1-SI-SG;             @ CONSUMPTION SHARE OF OUTPUT @
MU=0.0;                 @ GOV'T CONSUMPTION EFFECT ON PRIVATE MARGINAL
                          UTILITY @
THETA=(1/SN)*(NBAR/(1-NBAR))*(SC+MU*SG)
       /(1+(1/SN)*(NBAR/(1-NBAR))*(SC+MU*SG));
                        @ SHARE PARAMETER IN UTILITY FUNCTION @

@ ------------------------------------------------------------ @
@                       NORMALIZATIONS                         @
@ ------------------------------------------------------------ @

ZSS=1;

@ ------------------------------------------------------------ @
@                     STEADY STATE VALUES                      @
@ ------------------------------------------------------------ @

YSS=ZSS*((KYRATIO)^(SK/SN))*NBAR;
CSS=SC*YSS;
ISS=SI*YSS;
GSS=SG*YSS;
KSS=YSS*KYRATIO;
UTILITY=(1/(1-SIGMA))*((CSS+MU*GSS)^(THETA*(1-SIGMA)))
     *((TT-NBAR)^((1-THETA)*(1-SIGMA)));

@ ------------------------------------------------------------ @
@              SPECIFIKATION OF STOCHASTIC SHOCKS              @
@ ------------------------------------------------------------ @

SIGMAG=0.016*GSS;
RG=0.95;
SIGMAZ=0.012*ZSS;
RZ=0.98;

@ ------------------------------------------------------------ @
@                ORDERING OF VARIABLES                         @
@ ------------------------------------------------------------ @

@ 1  2  3  4  5
  K  Z  G  N  I @

VARN="K|Z|N|I";
VARNS="KONST"|"K"|"Z"|"G";
VARNC="N"|"I";

Z=KSS|ZSS|GSS|NBAR|ISS;   @ VECTOR WITH STEADY-STATE VALUES @
 
 

"2: APPROXIMATING THE OBJECTIVE FUNCTION";

@ ------------------------------------------------------------ @
@    THIS PROCEDURE IS USED TO DEFINE THE OBJECTIVE FUNCTION
     AFTER THIS IT IS USED TO COMPUTE THE JACOBIAN AND THE
     HESSIAN EVALUATED AT THE STEADY-STATE                     @
@ ------------------------------------------------------------ @

PROC GRADF(Z);

LOCAL Y, X1, X2, X3, X4, X5;
X1= Z[1,1]; X2= Z[2,1]; X3 = Z[3,1]; X4 = Z[4,1]; X5 = Z[5,1];

Y =(1/(1-SIGMA))*((X2*(X1^SK)*(X4^SN)-(1-MU)*X3-X5)^(THETA*(1-SIGMA)))*
   ((TT-X4)^((1-THETA)*(1-SIGMA)));

RETP(Y);
ENDP;

H = GRADP(&GRADF,Z);   @ THE GRADIENT OF THE FUNCTION DEFINED
                         IN GRADF IS COMPUTED AND EVALUATED AT
                         THE STEADY-STATE @

J = 0.5*HESSP(&GRADF,Z);   @ EQUIVALENTLY FOR THE HESSIAN @

@ ------------------------------------------------------------ @
@   IF YOU WISH TO CHECK YOUR RESULT, CHECK IF THE
    GRADIENT wrt. N EQUALS ZERO (WHY?)                         @
@ ------------------------------------------------------------ @
 

@ ------------------------------------------------------------

   NOW WE START THE APPROXIMATION OF THE MODEL

#  U IS APPROXIMATED BY: Z = X'QX + u'Ru + 2u'WX

   THIS LEAVES US WITH THE MAXIMIZATION PROBLEM:

#  Max sum(BSTAR^t)*Z(t) st. X(t+1) = AX(t) + Bu(t) + Ge(t)

   THE SOLUTION IS GIVEN BY:

#  u = -FX  WHERE  F = (INVERS(R + BSTAR*B'PB))(W + BSTAR*B'PA)

   WHERE P IS FOUND BY ITERATING ON:

#  P(i+1) = Q + BSTAR*A'P(i)A
            - (INV(BSTAR*A'P(i)B + W'))(R + BSTAR*B'P(i)B)(W + BSTAR*B'P(i)A)

   FOR AN INITIAL VALUE OF P(0) GIVEN BY A BEGATIVE DEFINIT MATRIX
   THIS SHOULD CONVERGE TO A UNIQUE LIMIT IN A FINITE NUMBER OF
   STEPS, SEE BERTSEKAS, 1976. P(0) SHOULD BE CHOSEN AS
   -EPSILON*EYE WHERE EPSILON IS A SMALL NUMBER AND EYE IS AN
   IDENTITY MATRIX. A GOOD GUESS CAN BE USED AS WELL.

   WE HAVE THAT:

   X(t) = { CONSTANT  K  Z  G }
   u(t) = { N  I }
                                                          4   2
                                                          q | w 4
   q, r, and w can be found from the Hessian: j = 0.5*H = w'| r 2

------------------------------------------------------------------- @

"3: DEFINITION OF MATRICES FOR MATRIX RICATTI ITERATIONER";

A = zeros(4,4); B = zeros(4,2); G=ZEROS(4,2);
Q = zeros(4,4); R = zeros(2,2); W = zeros(4,2);
P = zeros(4,4);

@ ------------------------------------------------------------ @
@                   DEFINITION OF Q, R AND W                   @
@ ------------------------------------------------------------ @
 

HH = 0.5*H-Z'*J;
R = J[4:5,4:5];
Q[1,1] = UTILITY - H*Z + Z'*J*Z;
Q[1,2:4] = HH[1,1:3]; Q[2:4,1] = Q[1,2:4]';
Q[2:4,2:4] = J[1:3,1:3];
W[2:4,1:2] = J[1:3,4:5];
W[1,1:2] = HH[1,4:5];

@ ------------------------------------------------------------ @
@                    DEFINITION OF A-MATRIX                    @
@ ------------------------------------------------------------ @
A[1,1] = 1;
A[2,2] = (1-DELTA)/GAMMAX;
A[3,1]=(1-RZ)*ZSS; A[3,3] = RZ;
A[4,1]=(1-RG)*GSS; A[4,4] = RG;

@ ------------------------------------------------------------ @
@                    DEFINITION OF B-MATRIX                    @
@ ------------------------------------------------------------ @
B[2,2] = 1/GAMMAX;

@ ------------------------------------------------------------ @
@                    DEFINITION OF G-MATRIX                    @
@ ------------------------------------------------------------ @
g[3,1] = 1; g[4,2] = 1;

@ ------------------------------------------------------------ @
@          LET'S ROCK AND ROLL: GET THE INTERATIONS STARTED    @
@ ------------------------------------------------------------ @

"4. START RICATTI";
DD = 1; I = 0;     @ DD IS THE SUP NORM USED FOR THE CONVERGENCE CRIT. @
LOADM PR;
PP=-0.000001*EYE(4); @ CHOICE OF P(0) @
PP=PR;             @ THIS INDICATES THAT I HAVE A GOOD GUESS @

DO UNTIL DD<0.00000000000000001;
locate 5,5;
P = Q + BSTAR*A'*PP*A - (BSTAR*A'*PP*B + W)*(INV(R + BSTAR*B'*PP*B))
                    *(W'+BSTAR*B'*PP*A);
PPP = P - PP;

D = MAXC(ABS(PPP));
DD = MAXC(D);
PP = P;
I = I+1;
I;;DD;
ENDO;

@ ------------------------------------------------------------ @
@                      CONVERGENCE IS OBTAINED                 @
@ ------------------------------------------------------------ @

pr=p;
save pr;

"   SUCCESS: IT HAS CONVERGED! NUMBER OF ITERATIONS ";;I-1;

@ ------------------------------------------------------------ @
@   THE UNIQUE SOLUTION FOR P HAS BEEN FOUND. THUS, WE ALSO
    KNOW F. THE CLOSED LOOP SOLUTION IS GIVEN BY:

  X(t+1) = (A-BF)X(t) + Ge(t)

    THE SYSTEM IS STABLE UNLESS THE GREATEST EIGENVALUE OF
    (A-BF) IS LESS THAN 1 (BSTAR^(-1) TO BE MORE PRECISE)
    THIS IS CHECKED BELOW                                      @
@ ------------------------------------------------------------ @
 

"5. DERIVATION OF OPTIMAL POLICY FUNCTIONS ";

FF1=(INV(R + BSTAR*B'*P*B));
FF2=(BSTAR*B'*P*A + W');
FF=FF1*FF2;

O = A - B*FF;
"6. CHECK FOR STABILITY, EGENVALUES SHOULD BE LESS THAN /BSTAR ";
"         EGENVALUES:";
{X1,X2,X3,X4}=EIGRG2(O);
X1;
" THAT'S IT - COMPUTATION TIME (SECONDS) ";;(HSEC-TX)/100;

SAVE FF; SAVE O;
 

OUTPUT FILE=STUD1.OUT RESET;
OUTPUT ON;
"---------------------------------------------------------";
"          THE SIMPLEST RBC MODEL                         ";
"---------------------------------------------------------";
"";
"----------- STEADY STATE RELATIONS AND PARAMETERS -------";
"";
"SHARE OF TIME USED ON THE MARKET          ";;NBAR;
"REAL INTEREST RATE                        ";;RR;
"LABOR SHARE OF INCOME                     ";;SN;
"CAPITAL SHARE OF INCOME                   ";;Sk;
"GROSS GROWTH RATE                         ";;GAMMAX;
"DISCOUNT FACTOR                           ";;BSTAR;
"RELATIVE RISK AVERSION                    ";;SIGMA;
"DEPRECIATION RATE                         ";;DELTA;
"CAPITAL-OUTPUT RATIO                      ";;KYRATIO;
"INVESTMENT-CAPITAL RATIO                  ";;IKRATIO;
"GOVERNMENT SHARE OF OUTPUT                ";;SG;
"INVESTMENT SHARE OF OUTPUT                ";;SI;
"PRIVATE CONSUMPTION SHARE OF OUTPUT       ";;SC;
"GOV'T CONS. EFFECT ON PRIV. MARG. UTILITY ";;MU;
"THE PREFERENCE PARAMETER THETA            ";;THETA;
"STANDARD DEVIATION G-SHOCKS               ";;SIGMAG;
"PERSISTENCE OF G-SHOCKS                   ";;RG;
"STANDARD DEV. OF Z-SHOCKS                 ";;SIGMAZ;
"PERSISTENCE OF  Z-SHOCKS                     ";;RZ;
"---------------------------------------------------------";
"DECISION RULES";
"1.STATE-VARIABLES (ENDOGENOUS AND EXOGENOUS)";
I=1;
DO UNTIL I>ROWS(VARNS);
$VARNS[I,1];;"=";
J=1;
DO UNTIL J>ROWS(VARNS);
O[I,J];;"*";;$VARNS[J,1];
J=J+1;
ENDO;
I=I+1;
ENDO;
"2.CONTROL-VARIABLES ";
I=1;
DO UNTIL I>ROWS(VARNC);
$VARNC[I,1];;"=";
J=1;
DO UNTIL J>ROWS(VARNS);
-FF[I,J];;"*";;$VARNS[J,1];
J=J+1;
ENDO;
I=I+1;
ENDO;
 

OUTPUT OFF;