// 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". 
//
// The program estimates the structural coefficients using the
// alternative information set: (D log Y_t, log L_t)
// 
// 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=1;        // 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]);



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


    numb=1000;

    CC_vec=zeros(2,numb);

        


    for i (1,numb,1);
        
        y=series_mat[((i-1)*2+1):(i*2),.];
        // transform series -------------------------------- 
        y2=y[2,.]; 
        y1=y[1,2:t]+(y2[.,2:t]-y2[.,1:t-1]);
        y=y1|(y2[.,2:t]); 
       
        //---------------------------------------------------
        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; 

        iSSCC=vec(inv(SS)*CC[.,1]);        
        // CONVERSION TO OLD VALUES ... 
        CC_veci=(issCC[1]-issCC[2])|(issCC[2]);


        CC_vec[.,i]=CC_veci;


    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;

    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 "True Mean StD"

    (CC_true*100)~CC_mean~CC_std;
    
    print "MSE"
    MSE;

 
