// monte_sss.prg 
// 
// This is a GAUSS program. 
// It performs the Monte Carlo simulation for the structural state 
// space model estimated by CCA for the paper "Business Cycle Analysis 
// and VARMA Models". 
// 
// Authors: Christian Kascha and Karel Mertens 
//




new;
library pgraph;
graphset;

#include "c:\\BCAandVARMArepfiles\\GAUSS\\LRidentification.ind";
#include "c:\\BCAandVARMArepfiles\\GAUSS\\getmodelname.ind";
#include "c:\\BCAandVARMArepfiles\\GAUSS\\CCAprocs.ind";

//------------------------------------------------------------
// SETTINGS:

t=180;              // Sample length
                    // t=180 as in CKM or CEV 

   
model_num=8;        // specifies the model: 
                    // 1 KP benchmark model  
                    // 2 KP w indivisible labor 
                    // 3 KP w Frisch E. 0.63
                    // 4 CKM benchmark model
                    // 5 CKM w indivisble labor
                    // 6 CKM w Frisch E. 0.63
                    // 7 CKM w sigma_l/2
                    // 8 CKM w sigma_l/3

lag_s=1;            // 1 - use f0,n0 lags 
                    // 2 - use IC crit
                    // for f,p selection

f0=8;              // used lags if lag_s==1;
n0=8;


s=2;                 // used state dimension

                    // Decomposition of variance-covariance matrix: 
dec_dum=1;          // 1 - use Cholesky decomposition
                    // 2 - solve by "hand" 
                    // 3 - solve and invert matrices by "hand"
                                


//------------------------------------------------------------
// GET SERIES and TRUE PARAMETERS

model_name=getmodelname(model_num);
dec_name=getdecname(dec_dum);
lags_name=getlagsname(lag_s,f0); 

filename1="c:\\BCAandVARMArepfiles\\output\\"$+"series_"$+model_name$+ftocv(t,1,0)$+".dat";
filename2="c:\\BCAandVARMArepfiles\\output\\"$+"coeff_"$+model_name$+".dat";

load series_mat[2000,t]= ^filename1;
load param[2,7]= ^filename2;

const=param[.,1];
ar_mat=param[.,2:3];
ma_mat=param[.,4:5];
cov_mat=param[.,6:7];
cov_mat=cov_mat*(cov_mat');


// GET TRUE IMPULSES
impuls=100|0;
horizon=10;
{phiU_true,phiE_true,CC_true,true_ir,true_ira}=varma_BQ(const,ar_mat,ma_mat,cov_mat,impuls,horizon);

CC_true=vec(CC_true[.,1]);
PhiE_true=PhiE_true[.,1]; 



/*-----------------------------------------------------------------------------------------------------------------------------------------------*/
// SIMULATION

    numb=1000;

    CC_vec=zeros(2,numb);
    PHIe_vec=zeros(2,numb); 

        


    for i (1,numb,1);
        
        y=series_mat[((i-1)*2+1):(i*2),.];
        y=y-meanc(y');
        SS=inv(diagrv(eye(2),stdc(y')));
        y=SS*y;

        //---------------------------------------------------
        if lag_s==1;
            nu=n0;
            fu=f0; 
        elseif lag_s==2; 
            nu=0;
            fu=0;  
        endif; 
        
        
        {A,B,C,res,x_ss,f,p,n_hat}=ss_CCAm(y,0,0,nu,fu,s);
        

        // check for stationarity etc. ...
        if (maxc(abs(eig(A)))>1) or (maxc(abs(eig(A-B*C)))>1);
             {A,B,C,res,x_ss,f,p,n_hat}=ss_CCArYW(y,0,0,nu+6,fu+6,s);
        endif; 

        if (maxc(abs(eig(A-B*C)))>1)or(maxc(abs(eig(A)))>1);
            1; 
        endif; 

        covmat=(res*(res'))/cols(res); 

        if dec_dum==1;
            {phiU,phiE,CC,y_ir,y_ira}=ss_BQ(A,C,B,covmat,horizon,impuls,eye(2));
        elseif dec_dum==2; 
            {phiU,phiE,CC,y_ir,y_ira}=ss_BQ2(A,C,B,covmat,horizon,impuls,eye(2));
        elseif dec_dum==3; 
            {phiU,phiE,CC,y_ir,y_ira}=ss_BQ3(A,C,B,covmat,horizon,impuls,eye(2));
        endif; 

        CC_vec[.,i]=vec(inv(SS)*CC[.,1]);
        PHIe_vec[.,i]=vec(inv(SS)*phiE[.,1]);


    endfor;







/*-----------------------------------------------------------------------------------------------------------------------------------------------*/
//   SOME STATISTICS:   
    format 3,3;
    MSE=CC_vec-ones(1,numb).*.CC_true;
    MSE=MSE.^2;
    MSE=meanc(MSE');

    CC_mean=meanc(CC_vec')*100;
    CC_std=stdc(CC_vec')*100;



    PHIe_vectr=trimm(PHIe_vec',40);
    PHIe_vectr=PHIe_vectr'; 

    MSEphi=PHIe_vectr-ones(1,cols(PHIe_vectr)).*.PHie_true;
    MSEphi=MSEphi.^2;
    MSEphi=meanc(MSEphi');

    PHIe_mean=meanc(PHIe_vectr')*100;
    PHIe_std=stdc(PHIe_vectr')*100;



    print "CCA ";
    print "Parametrization: "$+model_name;     
    print "Decomp. Method: "$+dec_name;
    print "Lag Selection: "$+lags_name; 
    print "No. of Obs.: "$+ftocv(t,1,0); 
    print "f:" f;
    print "p:" p;
    print "Impact --------------------------"; 
    print "True Mean StD"
    
    (CC_true*100)~CC_mean~CC_std;

    print "MSE"
    MSE;

    print "Long Run Effect -----------------"; 
    print "True Mean StD"
    
    (PHIe_true*100)~PHIe_mean~PHIe_std;

    print "MSE"
    MSEphi;
 



//------------------------------------------------------------------


proc(1)=trimm(x,k);
// computes columnwise trimmed x, where
// the lower k % and the upper k % of the column
// are thrown away   
local trimmc,xr,xtr,xi,kh,T,upperbound,lowerbound;  
local i,xd; 
    
    xtr={}; 
    kh=round(k/2); 
    T=rows(x); 
    upperbound=T-kh; 
    lowerbound=kh;     

    for i (1,cols(x),1); 

         xi=x[.,i]; 
         xr=rankindx(xi,1); 
         xd=(xr.>upperbound)+(xr.<lowerbound);            
         xi=delif(xi,xd);           
         xtr=xtr~xi;   

    endfor;

    retp(xtr); 
endp;












