// monte_svarma2SLS.prg 
// 
// This is a GAUSS program. 
// It performs the Monte Carlo simulation for the structural VARMA
// for the paper "Business Cycle Analysis and VARMA Models". 
// 
// Authors: Christian Kascha and Karel Mertens 
//

new;
library pgraph;
graphset;


#include "c:\\BCAandVARMArepfiles\\GAUSS\\varmaAUXprocs.ind";
#include "c:\\BCAandVARMArepfiles\\GAUSS\\varma2SLSprocs.ind";
#include "c:\\BCAandVARMArepfiles\\GAUSS\\LRidentification.ind";
#include "c:\\BCAandVARMArepfiles\\GAUSS\\getmodelname.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

                    // 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);

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

// SIMULATION PARAMETERS
    numb=1000;

    p=2;
    q=1;
    s=2; 

    {RR,r}=fma_restrictions(s,p,q);

    ARdumm=1; 
    nT=round(0.5*sqrt(t));

    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');

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

        {gamma_hut,vc_res,res}=varma2slsp(y,p,q,RR,r,nT,1,ARdumm);
        {AA0_e,AA_e,MM_e}=getvarma_mat(gamma_hut,RR,r,s,p,q);
        
       
 
        AA_e=inv(AA0_e)*AA_e;
        MM_e=inv(AA0_e)*MM_e;
        
        
        // check for stationarity 
        if (check_maxroot(MM_e)>1)or(check_maxroot2(AA_e)>1);
           
            {gamma_hut,vc_res,res}=varma2SLSp(y,p,q,RR,r,round(0.5*nT),1,3);
            {AA0_e,AA_e,MM_e}=getvarma_mat(gamma_hut,RR,r,s,p,q);
                 
            AA_e=inv(AA0_e)*AA_e;
            MM_e=inv(AA0_e)*MM_e;
            
            if (check_maxroot(MM_e)>1);
                MM_e=inv(MM_e);
            endif;  

            if (check_maxroot2(AA_e)>1);
                AA_e=inv(1.05)*AA_e;
            endif;     

        endif;

        cov_est=moment(res',0)/cols(res);      

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

        if dec_dum==1;
            {phiU,phiE,CC,y_ir,y_ira}=varma_BQ((0|0),AA_e,MM_e,cov_est,impuls,horizon);
        elseif dec_dumm=2;
            {phiU,phiE,CC,y_ir,y_ira}=varma_BQ((0|0),AA_e,MM_e,cov_est,impuls,horizon);
        elseif dec_dum==3; 
            {phiU,phiE,CC,y_ir,y_ira}=varma_BQ3((0|0),AA_e,MM_e,cov_est,impuls,horizon);
        endif; 
        //---------------------------------------------------

        
        CC_vec[.,i]=vec(CC[.,1]);
        PHIe_vec[.,i]=vec(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 "VARMA ";
    print "Parametrization: "$+model_name;     
    print "Decomp. Method: "$+dec_name;
    print "p:" p;
    print "q:" q;
    print "No. of Obs.: "$+ftocv(t,1,0); 
    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;




