new;

rndseed 90210; 

#include paramete
#include setup
#include fabe


vcvm=0; corm=0; varm=0;
stdv=0; stdvry=0; stdverr=0; correrr=0;

if (iscplx(fa)==1) or (iscplx(fb)==1) or (iscplx(fc)==1) or (iscplx(fd)==1);

if (maxc(maxc(abs(imag(fa))))<=critc) and (maxc(maxc(abs(imag(fb))))<=critc)
and (maxc(maxc(abs(imag(fc))))<=critc) and (maxc(maxc(abs(imag(fd))))<=critc);
fa=real(fa); fb=real(fb); fc=real(fc); fd=real(fd);

else;

print "imaginary matrices";
end;

endif;
endif;

it=1; do while it<=iterations;

xdv=zeros(ns,horizon);

err=chol(vcms)*rndn(2,horizon);

xdv[.,1]=(zeros(3,1)|err[.,1]);

dt = 2; do while dt<=horizon;

xdv[.,dt] = fb*xdv[.,dt-1]+fd*(zeros(3,1)|err[.,dt]);

dt=dt+1; endo;

ndv=fa*xdv;
dv=xdv|ndv;
v=dv+xbar;

gdp1=zeros(horizon,1);
gdp2=zeros(horizon,1);
nxy=zeros(horizon,1);

i=1;
do while i<=horizon;
gdp1[i]=f1(v[.,i]);
gdp2[i]=f2(v[.,i]);
nxy[i]=(gdp1[i]-v[8,i]-v[11,i])/gdp1[i];
i=i+1;
endo;

@--Now HP Filter the series--@

w=(v[3:5,.])'~ln(v[6:12,.])'~ln(gdp1)~ln(gdp2)~nxy;

@
1 bonds
2 shock1     
3 shock2    
4 labor1
5 labor2    
6 cons1      
7 cons2
8 Q  
9 inv1
10 inv2
11 GDP1
12 GDP2
13 net exports
@

@whp=hpmat*w;@

whp=w;

@--Now compute business cycle statistics--@

vcvm=vcvm+(vcx(whp)/iterations);
corm=corm+(corrx(whp)/iterations);

stdverr=stdverr+(stdc(err')/iterations);
correrr=correrr+(corrx(err')/iterations);

it=it+1;
endo;


varm=diag(vcvm);
stdv=100*varm^(1/2);
stdvry=stdv/stdv[11];

@print "% stdv b, z1, z2, n1, n2, c1, c2, q, x1, x2, y1, y2"; 
stdvry';@

print "BUSINESS CYCLE STATISTICS";
format /rds 6,3;

print "% stdv b, z1, n1,  c1,  q, x1,  y1"; 
stdv[1 2 4 6 8 9 11]';

print "corr (z1,z2), (y1,y2), (c1,c2), (n1,n2), (x1,x2) (nx, y)";
corm[2,3]~corm[11,12]~corm[6,7]~corm[4,5]~corm[9,10]~corm[11,13];