@ CHOOSE STANDARD ERRORS OF SHOCKS @
SA=0.00852;

@ USING THIS WE SET UP THE LOADING MATRIX FOR THE SHOCKS @
E=ZEROS(NS+NN,NS+NN);  @ NOTE THAT CAPITAL FORMS THE FIRST ELEMENT @
E[2,2]=SA;

@ NUMBER OF EXPERIMENTS @
NEXP=100;

@ SIMULATION HORIZON @
NIR=100;

@ 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=50;

@ 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 7654321;

@ INITIALIZE MATRICES OF MOMENTS @
STD=ZEROS(ROWS(MKE)+ROWS(H),NEXP);    @ ABSOLUTE STANDARD DEVIATIONS @
STDR=ZEROS(ROWS(MKE)+ROWS(H),NEXP);   @ STD.DEV. RELATIVE TO OUTPUT @
AU=ZEROS(ROWS(MKE)+ROWS(H),NEXP);     @ AUTOCORRELATIONS @

@ CORRELATIONS (NEXT) ARE ORDERED SO THAT WE GET CORRELATIONS WITH
  DOMESTIC OUTPUT, AND THE CORRELATIONS OF N AND APN @
COR=ZEROS((ROWS(MKE)+ROWS(H)+1),NEXP); @ CORRELATIONS @

@ VECTOR OF VARIABLE NAMES @
VAR1=NAME;
VAR2=NAME|"N-APN";
COR1=1|2|3|4|5|6|7|8|5;
COR2=6|6|6|6|6|6|6|6|7;

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

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

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

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,8);
IR1=SV1.*(Z+IR);

DAT=LN(IR1[.,1:8]);

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))*(IT1-1)+1:(ROWS(MKE)+ROWS(H))*IT1];
@ ABSOLUTE STANDARD DEVIATIONS @
STD[.,STACK*(IT-1)+IT1]=100*STDC(DAT);

@ STANDARD DEVIATIONS RELATIVE TO Y @
J=1;
SY1=STD[6,STACK*(IT-1)+IT1];
STDR[.,STACK*(IT-1)+IT1]=STD[.,STACK*(IT-1)+IT1]./SY1;

@ 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;

STD1=(MEANC(STD'));
SSTD=STDC(STD');
STDR1=MEANC(STDR');
SSTDR=STDC(STDR');
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)   RELATIVE TO Y  (Sd)  AUTOCOR (Sd)";
J=1;
DO UNTIL J>ROWS(STD);
$VAR1[J,1];;STD1[J,1];;SSTD[J,1];;STDR1[J,1];;SSTDR[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[J,1];;COR1[J,1];;SCOR[J,1];
J=J+1;
ENDO;
WAIT;
OUTPUT OFF;
CLEAR ALL;