@ =================================================================== @
@                    PROGRAM BY MORTEN RAVN                           @
@ =================================================================== @
@ =================================================================== @
@                         A TWO COUNTRY MODEL                         @
@    OF PRODUCT DIFFERENTIATION - USED BY BACKUS, KEHOE AND KYDLAND
     AER 1994, THE MODEL IS SOLVED USING RICATTI EQUATION ITERATIONS
     AND IS ANALYZED IN LOGS IN ORDER TO IMPROVE PRECISION            @
@ =================================================================== @
OUTPUT OFF;
"1. PROGRAMME STARTS";
TX=HSEC;
DCRIT=0.00000000000001;           @ CONVERGENCE CRITERION @

@ =================================================================== @
@                       ECONOMIC PARAMETER VALUES                     @
@ =================================================================== @

@ FACTOR SHARES (SK,SN) @

SN=0.64; SK=1-SN;

@ GROWTH RATE OF TECHNICAL CHANGE AND DEPRECIATION RATE @
GAMMAX=1.000; DELTAK=0.025;

@ DISCOUNT FACTOR AND REAL INTEREST RATE @
RR=0.04/4;
BSTAR=GAMMAX/(1+RR);

@ STEADY STATE SHARE OF TIME TO MARKET ACTIVITIES @
NSS=0.300;

@ GOVERNMENT POLICY PARAMETERS @
SG=0.2;
 

@ THE CAPITAL-OUTPUT RATIO @

KYRATIO=(BSTAR*SK)/(GAMMAX-BSTAR*(1-DELTAK));

@ STEADY STATE INVESTMENT SHARE @
SI=KYRATIO*(GAMMAX-(1-DELTAK));

@ IMPLIED CONSUMPTION SHARE @
SC=1-SI-SG;

@ GOVERNMENT EFFECT ON PRIVATE UTILITY @
MUG=0.0;

@ ELASTICITY OF SUBSTITUTION IN ARMINGTON-AGGREGATOR @
ELA=1.5;
ALFA=1/ELA;

@ IMPORT SHARE @
MS=0.15;

@ SHARES IN ARMINGTON-AGGREGATOR @
OMEGA1=(1-MS)^ALFA;
OMEGA2=MS^ALFA;

@ ==================================================================== @
@                        STEADY STATE CALCULATIONS                     @
@ ==================================================================== @
YSS=1;
ZSS=(YSS^SN)/((KYRATIO^SK)*NSS^SN);
A1SS=YSS*(1-MS);
B1SS=YSS*MS;
A2SS=B1SS;
B2SS=A1SS;
CSS=SC*YSS;
SS=1-SC-SG;                  @ SAVINGS SHARE @
SSS=SS*YSS;
XSS=SI*YSS;
ISS=XSS;
GSS=SG*YSS;
KSS=KYRATIO*YSS;
APSS=YSS/NSS;

@ VECTOR WITH STEADY STATE VALUES@
SSVAL=NSS|NSS|A1SS|B2SS|ZSS|ZSS|GSS|GSS|KSS|KSS|YSS|YSS|
      A2SS|B1SS|XSS|XSS|yss|yss|CSS|CSS|1|1|A2SS|
      B1SS|B1SS|A2SS|APSS|APSS|1|1|1|1;

@ VECTORS WITH NAMES OF VARIABLES @
VARus="N"|"N*"|"A1"|"B2"|"K'"|"K*'";            @ CONTROL VARIABLES - K'
                                                  DENOTES K NEXT PERIOD @
VARss="1"|"Z"|"Z*"|"G"|"G*"|"K"|"K*";           @ STATE VARIABLES @
varux="y"|"y*"|"a2"|"b1"|"i"|"i*"|"v"|"v*"|"c"|"c*"|"ph"|"pf"|
       "eh"|"ef"|"imh"|"imf"|"aph"|"apf"|"q1"|"q2"|"qf1"|"qf2";
                                                @ EXTRA CONTROL VARIABLES THAT
                                                  DO NOT ENTER PROBLEM @

@ PARAMETERS OF PREFERENCES
  SIGMA IS THE ABSOULTE VALUE OF THE ELASTICITY OF MARGINAL UTILITY
  OF CONSUMPTION @
SIGMA=2;

@ THETA WE DETERMINE SO THAT IT IS CONSISTENT WITH THE STEADY
  STATE NUMBER OF HOURS WORKED @
THETA=(((SC+MUG*SG)*NSS)/(SN*(1-NSS)))
      /(1+((SC+MUG*SG)*NSS)/(SN*(1-NSS)));

@ PARAMETERS OF THE STOCHASTIC PROCESSES FOR TECHNOLOGY AND GOV'T
  SPENDING @
RZZ=0.906; RZAZ=0.088;
@ RZZ IS THE PERSISTENCE PARAMETER, RZAZ IS THE CROSS-COUNTRY SPILLOVER @
Rg=0.95;

UTILITY=1/(1-SIGMA)*(CSS^((1-SIGMA)*THETA))*((1-NSS)^((1-SIGMA)*(1-THETA)));

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

@ 1  2  3   4  5   6  7   8  9   10  11  12  13
  1  Z  Z*  G  G*  K  K*  N  N*  A1  B2  K'  K*'@
@ THE FIRST VARIABLE HERE IS JUST A CONSTANT @

@ the steady-state values of the above variables @
Z=ln(ZSS|ZSS|GSS|GSS|KSS|KSS|
  NSS|NSS|A1SS|B2SS|KSS|KSS);

NC=6; NS=6;

"2: APPROXIMATING THE OBJECTIVE FUNCTION";

@ ------------------------------------------------------------ @
@  IN THIS PROCEDURE WE DEFINE THE OBJECTIVE FUNCTION          @
@ ------------------------------------------------------------ @
PROC GRADF(Z);

LOCAL Y, X1, X2, X3, X4, X5, x6, x7, x8, x9, X10, X11, X12;
X1 = Z[1,1]; X2 = Z[2,1]; X3 = Z[3,1]; X4 = Z[4,1]; X5 = Z[5,1];
X6 = Z[6,1]; X7 = Z[7,1]; X8 = Z[8,1]; X9 = Z[9,1]; X10= Z[10,1];
X11=Z[11,1]; X12=Z[12,1];

Y = .5/(1-SIGMA)*(
     (((OMEGA1*(exp(X9)^(1-ALFA))+
       OMEGA2*((exp(X2)*(exp(X8)^SN)
       *(exp(X6)^SK)-exp(X10))^(1-ALFA)))^(1/(1-ALFA))
       -GAMMAX*exp(X11)+(1-DELTAK)*exp(X5)-exp(X3))^(THETA*(1-SIGMA)))
     *((1-exp(X7))^((1-THETA)*(1-SIGMA)))
    +(((OMEGA1*(exp(X10)^(1-ALFA))+
       OMEGA2*((exp(X1)*(exp(X7)^SN)
      *(exp(X5)^SK)-exp(X9))^(1-ALFA)))^(1/(1-ALFA))
       -GAMMAX*exp(X12)+(1-DELTAK)*exp(X6)-exp(X4))^(THETA*(1-SIGMA)))
     *((1-exp(X8))^((1-THETA)*(1-SIGMA))) );
 

RETP(Y);
ENDP;
H = GRADP(&GRADF,Z);   @ THIS COMPUTES THE GRADIENTS OF THE ABOVE FUNCTION @

J = 0.5*HESSP(&GRADF,Z);   @ COMPUTATION OF THE HESSIAN MATRIX @
 

"3: DEFINITION OF MATRICES FOR THE RICATTI ITERATIONS";

NSC=NS+1;

A = zeros(NSC,NSC); B = zeros(NSC,NC);
Q11 = zeros(NSC,NSC); Q12=ZEROS(NSC,NC);
Q21 = ZEROS(NC,NSC); Q22=ZEROS(NC,NC);

P = zeros(NSC,NSC);

@ ------------------------------------------------------------ @
@                   DEFINITION OF Q                            @
@ ------------------------------------------------------------ @

HH = 0.5*H-Z'*J;

Q11[1,1] = UTILITY - H*Z + Z'*J*Z;
Q11[1,2:NSC] = HH[1,1:NS];  Q11[2:NSC,1]=Q11[1,2:NSC]';
Q11[2:NSC,2:NSC] = J[1:NS,1:NS];
Q12[1,1:NC] = HH[1,NS+1:NS+NC];  Q21[1:NC,1]=Q12[1,1:NC]';
Q12[2:NSC,1:NC]=J[1:NS,NS+1:NS+NC];
Q21[1:NC,2:NSC]=Q12[2:NSC,1:NC]';
Q22=J[NS+1:NS+NC,NS+1:NS+NC];

qtilde=(q11~q12)|(q21~q22);

@ ------------------------------------------------------------ @
@                    DEFINITION OF THE A-MATRIX                @
@ ------------------------------------------------------------ @
A[1,1] = 1;
A[2,1] = (1-RZZ-RZAZ)*(ln(ZSS)); A[2,2]=RZZ; A[2,3]=RZAZ;
A[3,1] = (1-RZZ-RZAZ)*(ln(ZSS)); A[3,2]=RZAZ; A[3,3]=RZZ;
A[4,1] = (1-rg)*(ln(GSS));  a[4,4]=rg;
A[5,1] = (1-rg)*(ln(GSS));  a[5,5]=rg;

@ ------------------------------------------------------------ @
@                    DEFINITION OF B-MATRIX                    @
@ ------------------------------------------------------------ @
B[6,5] = 1;
B[7,6] = 1;

@ ------------------------------------------------------------ @
@               DEFINITION OF DXU AND EXU-MATRISSER
                THESE DO NOT ENTER THE PROBLEM BUT ARE
                USED TO FIND THE POLICY RULES FOR THE
                EXTRA CONTROLS                                 @
@ ------------------------------------------------------------ @
DXU=ZEROS(rows(varux),rows(varus)); EXU=ZEROS(rows(varux),rows(varss));
 

@varux="y"|"y*"|"a2"|"b1"|"i"|"i*"|"v"|"v*"|"c"|"c*"|"ph"|"pf"|
       "eh"|"ef"|"imh"|"imf"|"aph"|"apf"|"q1"|"q2"|"qf1"|"qf2";
VARus="N"|"N*"|"A1"|"B2"|"K'"|"K*'";@
DXU[1,1]=sn;
DXU[2,2]=sn;
dxu[3,.]=yss/a2ss*dxu[1,.]; dxu[3,3]=dxu[3,3]-a1ss/a2ss;
dxu[4,.]=yss/a2ss*dxu[2,.]; dxu[4,4]=dxu[4,4]-a1ss/a2ss;
DXU[5,5]=KSS/ISS;
DXU[6,6]=KSS/ISS;
dxu[7,.]=dxu[4,.]*omega2*(b1ss/yss)^(1-alfa);
dxu[7,3]=dxu[7,3]+omega1*(a1ss/yss)^(1-alfa);
dxu[8,.]=dxu[3,.]*omega2*(b1ss/yss)^(1-alfa);
dxu[8,4]=dxu[8,4]+omega1*(a1ss/yss)^(1-alfa);
dxu[9,.]=dxu[7,.]*yss/css-dxu[5,.]*iss/css;
dxu[10,.]=dxu[8,.]*yss/css-dxu[6,.]*iss/css;
dxu[11,.]=dxu[4,.]*(-alfa); dxu[11,3]=dxu[11,3]+alfa;
dxu[12,.]=dxu[3,.]*(-alfa); dxu[12,4]=dxu[12,4]+alfa;
dxu[13,.]=dxu[3,.]; dxu[14,.]=dxu[4,.];
dxu[15,.]=dxu[11,.]+dxu[4,.];
dxu[16,.]=dxu[12,.]+dxu[3,.];
dxu[17,.]=dxu[1,.]; dxu[17,1]=dxu[17,1]-1;
dxu[18,.]=dxu[2,.]; dxu[18,2]=dxu[18,2]-1;
dxu[19,.]=dxu[11,.]*(1-omega2^(1/alfa)/(omega1^(1/alfa)+omega2^(1/alfa)))
         +dxu[12,.]*(omega2^(1/alfa))/(omega1^(1/alfa)+omega2^(1/alfa));
dxu[20,.]=dxu[12,.]*(1-omega2^(1/alfa)/(omega1^(1/alfa)+omega2^(1/alfa)))
         +dxu[11,.]*(omega2^(1/alfa))/(omega1^(1/alfa)+omega2^(1/alfa));
dxu[21,.]=dxu[11,.]*(1-ms)+dxu[12,.]*ms;
dxu[22,.]=dxu[12,.]*(1-ms)+dxu[11,.]*ms;
 

@VARss="1"|"Z"|"Z*"|"G"|"G*"|"K"|"K*";@
exu[1,2]=1; exu[1,6]=sk;
exu[2,3]=1; exu[2,7]=sk;
exu[3,.]=yss/a2ss*exu[1,.];
exu[3,1]=exu[3,1]+ln(a2ss)-yss/a2ss*ln(yss)+a1ss/a2ss*ln(a1ss);
exu[4,.]=yss/a2ss*exu[2,.];
exu[4,1]=exu[4,1]+ln(a2ss)-yss/a2ss*ln(yss)+a1ss/a2ss*ln(a1ss);
exu[5,1]=ln(iss)-deltak*kss/iss*ln(kss); exu[5,6]=-kss/iss*(1-deltak);
exu[6,1]=ln(iss)-deltak*kss/iss*ln(kss); exu[6,7]=-kss/iss*(1-deltak);
exu[7,.]=exu[4,.]*omega2*(b1ss/yss)^(1-alfa);
exu[7,1]=exu[7,1]+ln(yss)-ln(a1ss)*omega1*(a1ss/yss)^(1-alfa)
         -ln(b1ss)*omega2*(b1ss/yss)^(1-alfa);
exu[8,.]=exu[3,.]*omega2*(b1ss/yss)^(1-alfa);
exu[8,1]=exu[8,1]+ln(yss)-ln(a1ss)*omega1*(a1ss/yss)^(1-alfa)
         -ln(b1ss)*omega2*(b1ss/yss)^(1-alfa);
exu[9,.]=exu[7,.]*yss/css-exu[5,.]*iss/css;
exu[9,1]=exu[9,1]+ln(css)-yss/css*ln(yss)+iss/css*ln(iss)+gss/css*ln(gss);
exu[9,4]=exu[9,4]-gss/css;
exu[10,.]=exu[8,.]*yss/css-exu[6,.]*iss/css;
exu[10,1]=exu[10,1]+ln(css)-yss/css*ln(yss)+iss/css*ln(iss)+gss/css*ln(gss);
exu[10,5]=exu[10,5]-gss/css;
exu[11,.]=exu[4,.]*(-alfa); exu[11,1]=exu[11,1]+ln(omega2/omega1);
exu[12,.]=exu[3,.]*(-alfa); exu[12,1]=exu[12,1]+ln(omega2/omega1);
exu[13,.]=exu[3,.]; exu[14,.]=exu[4,.];
exu[15,.]=exu[11,.]+exu[4,.];
exu[16,.]=exu[12,.]+exu[3,.];
exu[17,.]=exu[1,.]; exu[18,.]=exu[2,.];
exu[19,.]=exu[11,.]*(1-omega2^(1/alfa)/(omega1^(1/alfa)+omega2^(1/alfa)))
         +exu[12,.]*(omega2^(1/alfa))/(omega1^(1/alfa)+omega2^(1/alfa));
exu[20,.]=exu[12,.]*(1-omega2^(1/alfa)/(omega1^(1/alfa)+omega2^(1/alfa)))
         +exu[11,.]*(omega2^(1/alfa))/(omega1^(1/alfa)+omega2^(1/alfa));
exu[21,.]=exu[11,.]*(1-ms)+exu[12,.]*ms;
exu[22,.]=exu[12,.]*(1-ms)+exu[11,.]*ms;

@ ------------------------------------------------------------ @
@          STARTING THE MATRIX-RICATTI ITERATIONERNE           @
@ ------------------------------------------------------------ @

"4. START RICATTI";
DD = 1; I = 0;
PP=-0.0001*EYE(NSC);     @ HERE YOU CHOOSE THE INITIAL VALUE OF P(0) @
@LOADM POLD; PP=POLD;@

DO until DD<DCRIT;
 locate 5,5;
  ATILDE=Q21+BSTAR*B'*PP*A;
  BTILDE=Q22+BSTAR*B'*PP*B;

  p=Q11+BSTAR*A'*PP*A-(ATILDE')*(INV(BTILDE))*ATILDE;
  ppp=p-pp;
  D = MAXC(ABS(PPP));
  DD = MAXC(D);
  PP = P;
 I = I+1;
 I;;DD;
ENDO;

POLD=PP; SAVE POLD;

@ ------------------------------------------------------------ @
@                      CONVERGENCE HAS BEEN OBTAINED           @
@ ------------------------------------------------------------ @

"   RICATTI ITERATIONS ARE DONE. NO. OF ITERATIONS: ";;I-1;
  ATILDE=Q21+BSTAR*B'*P*A;
  BTILDE=Q22+BSTAR*B'*P*B;
 

"5. NOW WE COMPUTE THE POLICY FUNCTIONS";

@ (a) for control variables @
AU=-(INV(BTILDE))*ATILDE;

@ (B) for STATE variables @
AS=A+B*AU;

@ (B) for implied variables @
AXU=DXU*AU+EXU;
 

"6. CHECK FOR STABILITY, EIGENVALUES SHOULD BE LESS THAN 1/BSTAR";
"         EIGENVALUES:";
{X1,X2,X3,X4}=EIGRG2(AS);
X1;
" WE ARE DONE - COMPUTATION TIME: ";;(HSEC-TX)/100;
WAIT;

FORMAT /LDN 12,7;
OUTPUT FILE=STUD1.OUT RESET;
OUTPUT ON;

"---------------------------------------------------------";
"                    THE BKK MODEL                        ";
"---------------------------------------------------------";
"----------- STEADY STATE RELATIONER OG PARAMETRE --------";
"";
"TIME SPENT AT MARKET ACTIVITIES (SHARE)   ";;NSS;
"REAL INTEREST RATE                        ";;(1-bstar)/bstar;
"LABOUR SHARE                              ";;SN;
"CAPITAL SHARE                             ";;SK;
"DISCOUNT FACTOR                           ";;BSTAR;
"DEPRECIATION RATE                         ";;DELTAK;
"CAPITAL-OUTPUT RATE                       ";;KYRATIO;
"INVESTMENT-CAPITAL RATE                   ";;SI/KYRATIO;
"INVESTMENT SHARE OF OUTPUT                ";;SI;
"CONSUMPTION SHARE OF OUTPUT               ";;SC;
"GOVERNMENT SHARE OF OUTPUT                ";;SG;
"GOVERNMENT EFFECT ON PRIVATE UTILITY      ";;MUG;
"ELASTICITY OF SUBSTITUTION                ";;ELA;
"IMPORT SHARE                              ";;MS;
"SHARE PARAMETER, OMEGA1                   ";;OMEGA1;
"SHARE PARAMETER, OMEGA2                   ";;OMEGA2;
"PERSISTENCE OF TECHN. SHOCKS              ";;RZZ;
"SPILLOVERS OF TECHN. SHOCKS               ";;RZAZ;
"PERSISTENCE OF GOV. SPENDING SHOCKS       ";;RG;
"---------------------------------------------------------";
"DECISION RULES (ALL VARIABLES IN LOGS)";
"1.AGGREGATE STATE-VARIABLES ";
I=1;
DO UNTIL I>ROWS(VARSS);
$VARSS[I,1];;"=";
J=1;
DO UNTIL J>ROWS(VARSS);
AS[I,J];;"*";;$VARSS[J,1];
J=J+1;
ENDO;
I=I+1;
ENDO;
"2.CONTROL VARIABLES ";
I=1;
DO UNTIL I>ROWS(VARUS);
$VARUS[I,1];;"=";
J=1;
DO UNTIL J>ROWS(VARSS);
AU[I,J];;"*";;$VARSS[J,1];
J=J+1;
ENDO;
I=I+1;
ENDO;
 
 

"3.EXTRA KONTROL-VARIABLE ";
I=1;
DO UNTIL I>ROWS(VARUX);
$VARUX[I,1];;"=";
J=1;
DO UNTIL J>ROWS(VARSS);
AXU[I,J];;"*";;$VARSS[J,1];
J=J+1;
ENDO;
I=I+1;
ENDO;
OUTPUT OFF;
 

"test of steady state reproduction ";
aslr=eye(nsc);
i=1; do until i>2000;
aslr=as*aslr;
i=i+1;
endo;
zz=1|z[1:ns,1];
zzz=aslr*zz;
"steady state of state variables";
zz';
"implied steady state";
zzz';
"steady state of control variables";
z[ns+1:ns+nc,1]';
"implied steady state";
zzzz=(au*aslr*zz);
"steady state of extra control variables";
ln(ssval[rows(ssval)-rows(varux)+1:rows(ssval),1]');
"implied steady state";
zzzz=(axu*aslr*zz);
zzzz';