% Distribution of Multiplier

% Calculation of shape (aa) and scale (bb) parameters

mean_alphag=.19127 ;
var_alphag=0.00309297743171596;
mean_omega=0.1816 ;
var_omega=0.00173573471500182;
mean_etapi=0.5010 ;
var_etapi=0.00406014799220082;
mean_etay=0.0650 ;
var_etay=0.000134845723028720;
mean_etarg=0.8273;
var_etarg=0.0157985444332002;
mean_iotap=0.9438 ;
var_iotap=0.0334740151412075;
mean_iotaw=0.8875 ;
var_iotaw=0.0334695730748423;
mean_chi= 0.8826;
var_chi = 0.00340413831891982 ;
mean_phi= 2.2280;
var_phi= 0.148941604366296;
mean_rhocg= 0.7713;
var_rhocg = 0.00223042663088622;
mean_rhoig= 0.5005;
var_rhoig = 0.00500729188967140 ;

% Gamma Distribution

aa_iotap = mean_iotap^2/var_iotap;
bb_iotap = var_iotap/mean_iotap;
aa_iotaw = mean_iotaw^2/var_iotaw;
bb_iotaw = var_iotaw/mean_iotaw;
aa_phi = mean_phi^2/var_phi;
bb_phi = var_phi/mean_phi;

% Beta Distribution

aa_chi = ((1-mean_chi)*mean_chi^2/var_chi) - mean_chi;
bb_chi = (((1-mean_chi)*mean_chi/var_chi)-1)*(1-mean_chi);
aa_rhocg = ((1-mean_rhocg)*mean_rhocg^2/var_rhocg) - mean_rhocg;
bb_rhocg = (((1-mean_rhocg)*mean_rhocg/var_rhocg)-1)*(1-mean_rhocg);
aa_rhoig = ((1-mean_rhoig)*mean_rhoig^2/var_rhoig) - mean_rhoig;
bb_rhoig = (((1-mean_rhoig)*mean_rhoig/var_rhoig)-1)*(1-mean_rhoig);
aa_alphag = ((1-mean_alphag)*mean_alphag^2/var_alphag) - mean_alphag;
bb_alphag = (((1-mean_alphag)*mean_alphag/var_alphag)-1)*(1-mean_alphag);
aa_omega = ((1-mean_omega)*mean_omega^2/var_omega) - mean_omega;
bb_omega = (((1-mean_omega)*mean_omega/var_omega)-1)*(1-mean_omega);
aa_etapi = ((1-mean_etapi)*mean_etapi^2/var_etapi) - mean_etapi;
bb_etapi = (((1-mean_etapi)*mean_etapi/var_etapi)-1)*(1-mean_etapi);
aa_etay = ((1-mean_etay)*mean_etay^2/var_etay) - mean_etay;
bb_etay = (((1-mean_etay)*mean_etay/var_etay)-1)*(1-mean_etay);
aa_etarg = ((1-mean_etarg)*mean_etarg^2/var_etarg) - mean_etarg;
bb_etarg = (((1-mean_etarg)*mean_etarg/var_etarg)-1)*(1-mean_etarg);

% Generation of Random Draws


pd_iotap = makedist('Gamma','a', aa_iotap, 'b', bb_iotap);
pd_iotaw = makedist('Gamma','a', aa_iotaw, 'b', bb_iotaw);
pd_phi = makedist('Gamma','a', aa_phi, 'b', bb_phi);
pd_chi = makedist('Beta','a', aa_chi, 'b', bb_chi);
pd_rhocg = makedist('Beta','a', aa_rhocg, 'b', bb_rhocg);
pd_rhoig = makedist('Beta','a', aa_rhoig, 'b', bb_rhoig);
pd_alphag = makedist('Beta','a', aa_alphag, 'b', bb_alphag);
pd_omega = makedist('Beta','a', aa_omega, 'b', bb_omega);
pd_etapi = makedist('Beta','a', aa_etapi, 'b', bb_etapi);
pd_etay = makedist('Beta','a', aa_etay, 'b', bb_etay);
pd_etarg = makedist('Beta','a', aa_etarg, 'b', bb_etarg);

rng('default') % For reproducibility
%s = rng;

rand_omega= random(pd_omega,20000,1);
rand_etapi= random(pd_etapi,20000,1);
rand_etay= random(pd_etay,20000,1);
rand_etarg= random(pd_etarg,20000,1);
rand_iotap= random(pd_iotap,20000,1);
rand_iotaw= random(pd_iotaw,20000,1);
rand_alphag= random(pd_alphag,20000,1);
rand_chi= random(pd_chi,20000,1);
rand_phi= random(pd_phi,20000,1);
rand_rhocg= random(pd_rhocg,20000,1);
rand_rhoig= random(pd_rhoig,20000,1);

gdp_cg_mul=NaN(1,500);
cp_cg_mul=NaN(1,500) ;
cpr_cg_mul= NaN(1,500);
cpnr_cg_mul=NaN(1,500) ;
ip_cg_mul=NaN(1,500);
gdp_steady=NaN(1,500);

first_time = 1;

for i=1:500
    
    if first_time
        dynare Keynes_multiplier.mod nostrict noclearall;
        
        first_time = 0;
    else
        set_param_value('aalpha_g', rand_alphag(i));
        set_param_value('oomega', rand_omega(i));
        set_param_value('etapi', rand_etapi(i));
        set_param_value('etay', rand_etay(i));
        set_param_value('etarg', rand_etarg(i));
        set_param_value('iotap', rand_iotap(i));
        set_param_value('iotaw', rand_iotaw(i));
        set_param_value('cchi', rand_chi(i));
        set_param_value('pphi', rand_phi(i));
        set_param_value('rho_cg', rand_rhocg(i));
        set_param_value('rho_ig', rand_rhoig(i));
        
        info = stoch_simul(var_list_);
        
        if info(1)
            disp('Failed Stoch_Simul')
        else
            gdp_cg_mul(i)= (oo_.irfs.gdp_eps_cg(1)*exp(oo_.steady_state(1)))/(.01*exp(oo_.dr.ys(4)));
            cp_cg_mul(i)= (oo_.irfs.cp_eps_cg(1)*exp(oo_.steady_state(2)))/(.01*exp(oo_.dr.ys(4)));
            cpr_cg_mul(i)= (oo_.irfs.cpr_eps_cg(1)*exp(oo_.steady_state(6)))/(.01*exp(oo_.dr.ys(4)));
            cpnr_cg_mul(i)= (oo_.irfs.cpnr_eps_cg(1)*exp(oo_.steady_state(7)))/(.01*exp(oo_.dr.ys(4)));
            ip_cg_mul(i)= (oo_.irfs.ip_eps_cg(1)*exp(oo_.steady_state(3)))/(.01*exp(oo_.dr.ys(4)));
            gdp_steady(i)=exp(oo_.dr.ys(1));
            
        end
        
    end
    
    
    
end


