@ 1. CHOOSE THE MOMENTS OF THE INNOVATION PROCESSES
  YOU HAVE TO:
  A) CHOOSE THE STANDARD ERRORS
  B) CHOOSE THE COVARIANCE STRUCTURE, IF YOU HAVE NON-ZEROS
     OUTSIDE THE DIAGONAL, YOU HAVE TO MAKE A CHOLESKI
     DECOMPOSITION @

@ A) CHOOSE STANDARD ERRORS OF SHOCKS @

@ B) CHOOSE CORRELATIONS @
 

@ USING THIS WE SET UP THE LOADING MATRIX FOR THE SHOCKS @
E=ZEROS(6,6);             @ NOTE THAT CAPITAL FORMS THE TWO FIRST ELEMENTS @
/* FILL IN THE MATRIX */
 

@ NUMBER OF EXPERIMENTS @
NEXP=;

@ SIMULATION HORIZON @
NIR=;

@ CHOOSE THE STACKING OF EXPERIMENTS - HERE THERE IS A TRADE OFF:
  THE MORE EXPERIMENTS THAT ARE STACKED THE FASTER BUT ONE
  RISKS RUNNING OUT OF MEMORY - REMEMBER THAT NUMBER OF EXPERIMENTS
  MUST BE AN INTEGER TIMES STACKING @
STACK=;

@ MATRIX OF STEADY STATE VALUES @
SV1=ONES(NIR,ROWS(MKE)+ROWS(H)).*SSVAL';

@ INCLUDE HPFILTER @
UU=0;
#INCLUDE HP1;

@ FINALLY CHOOSE THE SEED FOR THE RANDOM NUMBER GENERATOR @
RNDSEED;

@ INITIALIZE MATRICES OF MOMENTS @
STD=ZEROS(ROWS(MKE)+ROWS(H),NEXP);    @ ABSOLUTE STANDARD DEVIATIONS @
AU=ZEROS(ROWS(MKE)+ROWS(H),NEXP);     @ AUTOCORRELATIONS @

@ CORRELATIONS ARE ORDERED IN ANY WAY THAT YOU WANT @
COR=ZEROS(,NEXP);  @ CORRELATIONS @

@ VECTOR OF VARIABLE NAMES @
VAR1 = { };                             @ FOR STD,AU @
VAR2 = { };                             @ FOR COR @

@ CORRELATIONS ARE COMPUTED BETWEEN VARIABLE I OF COR1 AND VARIABLE
  I OF COR2 @
COR1=;
COR2=;

LOCATE 1,5;"SIMULATING ROUND NO:";

IT=1;
DO UNTIL IT>NEXP/STACK;
DAT1=ZEROS(NIR,ROWS(MKE)+ROWS(H)+1);
IT1=1;
DO UNTIL IT1>STACK;

LOCATE 1,28;STACK*(IT-1)+IT1;
 

@ NOW WE COMPUTE THE EQUILIBRIUM SEQUENCE USING A DO-LOOP.
  THIS CAN BE DONE MORE EFFICIENTLY USING RECSERAR BUT AT
  THE RISK OF MAKING MISTAKES - TRY IT OUT IF YOU WISH @
IR=ZEROS(NIR,ROWS(MKE)+ROWS(H));
SS=E*RNDN(NIR+1,ROWS(E))';
S=SS[.,1];
I=1;
DO UNTIL I>NIR;
IR[I,1:ROWS(MKE)]=S';
IR[I,ROWS(MKE)+1:COLS(IR)]=(H*S)';
S=MKE*S+SS[.,I+1];
I=I+1;
ENDO;
 

Z=(ONES(NIR,ROWS(MKE)+ROWS(H)));
IR1=SV1.*(Z+IR);

DAT=LN(IR1);

DAT1=DAT1~DAT;

IT1=IT1+1;
ENDO;

DAT=DAT1[.,ROWS(MKE)+ROWS(H)+1:COLS(DAT1)];
UU=ZEROS(ROWS(DAT),COLS(DAT));
CALL HP1(DAT);

CLEAR DAT1;
CLEAR DAT;
IT1=1;
DO UNTIL IT1>STACK;
DAT=UU[.,(ROWS(MKE)+ROWS(H)+1)*(IT1-1)+1:(ROWS(MKE)+ROWS(H)+1)*IT1];

@ ABSOLUTE STANDARD DEVIATIONS @
STD[.,STACK*(IT-1)+IT1]=100*STDC(DAT);
 

@ AUTOCORRELATIONS @
X=DAT[2:NIR,.]; X1=DAT[1:NIR-1,.];
X2=X'*X1;
X3=SQRT(DIAG(X'*X))*(SQRT(DIAG(X1'*X1)))';
X4=X2./X3;
AU[.,STACK*(IT-1)+IT1]=DIAG(X4);

@ CORRELATION MATRIX @
X=DAT'*DAT;
X1=SQRT(DIAG(X));
X2=X1*X1';
X3=X./X2;
J=1;
DO UNTIL J>ROWS(COR);
COR[J,STACK*(IT-1)+IT1]=X3[COR1[J,1],COR2[J,1]];
J=J+1;
ENDO;

IT1=IT1+1;
ENDO;
CLEAR DAT;

IT=IT+1;
ENDO;

@ COMPUTATIONS OF THE MEAN MOMENTS OF THE SIMULATIONS
  NORMALIZE THE STDC-MEANSURES BY SQRT(NEXP) IF YOU WANT
  STANDARD ERRORS @

STD1=(MEANC(STD'));
SSTD=STDC(STD');
AU1=MEANC(AU');
SAU=STDC(AU');
COR1=MEANC(COR');
SCOR=STDC(COR');
 
 

LOCATE 2,5; " TOTAL TIME";;(HSEC-T)/100;
WAIT;

@ PRINTING RESULTS OUT @
OUTPUT FILE=KPR.OUT;
OUTPUT ON;
"DATE";;(DATE)';
"NO. OF SIMULATIONS";;NEXP;
"";
" STANDARD DEVIATIONS AND AUTOCORRELATIONS";
"VARIABLE   ABSOLUTE  (Sd)   AUTOCOR (Sd)";
J=1;
DO UNTIL J>ROWS(STD);
$VAR1[1,J];;STD1[J,1];;SSTD[J,1];;AU1[J,1];;SAU[J,1];
J=J+1;
ENDO;
WAIT;
" CONTEMPORANOUS CORRELATIONS (WITH Y UNLESS OTHERWISE STATED)";
"VARIABLE   COR  (Sd) ";
J=1;
DO UNTIL J>ROWS(COR1);
$VAR2[1,J];;COR1[J,1];;SCOR[J,1];
J=J+1;
ENDO;
WAIT;
OUTPUT OFF;
CLEAR ALL;