// monte_svarAI.prg 
// 
//
// The program performs the Monte Carlo simulation for the structural VAR
// 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\\varestimation.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 
                    // for all parametrizations
                    // t=500 
                    // t=1000
                    // for the KP or CKM benchmark models
   
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

p=4;                
lag_s=1;            // 1 - use p lags 
                    // 2 - use AIC crit. for lag-selection

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


// ORIGINAL DATA PROCESS AND USUAL COEFFICIENTS FOR TRUE PROCESS ... 
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

// SIMULATION PARAMETERS
    numb=1000;
    s=2;
    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]); 
       
        //--------------------------------------------------

        SS=inv(diagrv(eye(2),stdc(y')));
        y=SS*y;

        //---------------------------------------------------
        if lag_s==1;
            p_lag=p; 
        elseif lag_s==2; 
            {aic,aic_p,hq,hq_p,sc,sc_p}=ord_VAR(y,1,1,round(sqrt(t)),1);
            p_lag=aic_p; 
        endif; 

        //---------------------------------------------------
        {a_est,cov_est,res_est}=var_ols(y,1,p_lag);
        
        // check for stationarity: 
        big_A=a_est[.,2:cols(a_est)];
        if p_lag ne 1;
            ba=eye(s*(p_lag-1))~zeros(s*(p_lag-1),s);
            big_A=big_A|ba;
        endif;  
    
        if maxc(abs(eig(big_A)))>1; 
            {a_est,cov_est,res_est}=var_yw(y,p_lag);
            a_est=ones(s,1)~a_est;  
        endif;
        
        cov_est=moment(res_est',0)/cols(res_est);      
        //---------------------------------------------------

        if dec_dum==1;
            {phiU,phiE,CC,y_ir,y_ira}=varma_BQ((0|0),a_est[.,2:(p_lag*s+1)],zeros(2,2),cov_est,impuls,horizon);
        elseif dec_dum==2; 
            {phiU,phiE,CC,y_ir,y_ira}=varma_BQ2((0|0),a_est[.,2:(p_lag*s+1)],zeros(2,2),cov_est,impuls,horizon);
        elseif dec_dum==3; 
            {phiU,phiE,CC,y_ir,y_ira}=varma_BQ3((0|0),a_est[.,2:(p_lag*s+1)],zeros(2,2),cov_est,impuls,horizon);
        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 "VAR ";
    print "Parametrization: "$+model_name;     
    print "Decomp. Method: "$+dec_name;
    print "Lag Selection: "$+lags_name; 
    print "No. of Obs.: "$+ftocv(t,1,0); 
    print "True Mean StD"
    
    (CC_true*100)~CC_mean~CC_std;

    print "MSE"
    MSE;









