% Distribution of Multiplier

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

mean_alphag=.1882 ;
var_alphag=0.00310445897051282;
mean_omega=0.1823 ;
var_omega=0.00172161400805937;
mean_gammap= 0.2464;
var_gammap=0.00212035378565345;
mean_gammaw=0.2810 ;
var_gammaw=0.00289445876459664;
mean_etapi=0.5004 ;
var_etapi=0.00396639646149869;
mean_etay=0.0649 ;
var_etay=0.000134845723028720;
mean_etarg=0.8297;
var_etarg=0.0160662588797905;
mean_iotap=0.9399 ;
var_iotap=0.0333207458436709;
mean_chi= 0.8836;
var_chi = 0.00349229803740436 ; 
mean_zeta= 1.4919;
var_zeta = 0.0896967502484217 ;
mean_phi= 2.2369;
var_phi= 0.159571089810899;
mean_rhocg= 0.7717;
var_rhocg = 0.00205151863562418; 
mean_rhoig= 0.4976;
var_rhoig = 0.00508752454576657 ; 

% Gamma Distribution

aa_iotap = mean_iotap^2/var_iotap;
bb_iotap = var_iotap/mean_iotap;
aa_zeta = mean_zeta^2/var_zeta;
bb_zeta = var_zeta/mean_zeta;
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_gammap = ((1-mean_gammap)*mean_gammap^2/var_gammap) - mean_gammap;
bb_gammap = (((1-mean_gammap)*mean_gammap/var_gammap)-1)*(1-mean_gammap);
aa_gammaw = ((1-mean_gammaw)*mean_gammaw^2/var_gammaw) - mean_gammaw;
bb_gammaw = (((1-mean_gammaw)*mean_gammaw/var_gammaw)-1)*(1-mean_gammaw);
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_zeta = makedist('Gamma','a', aa_zeta, 'b', bb_zeta);
pd_iotap = makedist('Gamma','a', aa_iotap, 'b', bb_iotap);
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_gammap = makedist('Beta','a', aa_gammap, 'b', bb_gammap);
pd_gammaw = makedist('Beta','a', aa_gammaw, 'b', bb_gammaw);
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_gammap= random(pd_gammap,20000,1); 
rand_gammaw= random(pd_gammaw,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_alphag= random(pd_alphag,20000,1);   
rand_zeta= random(pd_zeta,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); 


for i=1:20000

dynare keynes_multiplier.mod;

% Calculation of Impact Multiplier

gdp_cg_mul(i)= (gdp_eps_cg(1)*exp(oo_.steady_state(1)))/(.01*exp(oo_.steady_state(4)));
cp_cg_mul(i)= (cp_eps_cg(1)*exp(oo_.steady_state(2)))/(.01*exp(oo_.steady_state(4)));
cpr_cg_mul(i)= (cpr_eps_cg(1)*exp(oo_.steady_state(6)))/(.01*exp(oo_.steady_state(4)));
cpnr_cg_mul(i)= (cpnr_eps_cg(1)*exp(oo_.steady_state(7)))/(.01*exp(oo_.steady_state(4)));
ip_cg_mul(i)= (ip_eps_cg(1)*exp(oo_.steady_state(3)))/(.01*exp(oo_.steady_state(4)));
gdp_ig_mul(i)= (gdp_eps_ig(1)*exp(oo_.steady_state(1)))/(.01*exp(oo_.steady_state(5)));
cp_ig_mul(i)= (cp_eps_ig(1)*exp(oo_.steady_state(2)))/(.01*exp(oo_.steady_state(5)));
cpr_ig_mul(i)= (cpr_eps_ig(1)*exp(oo_.steady_state(6)))/(.01*exp(oo_.steady_state(5)));
cpnr_ig_mul(i)= (cpnr_eps_ig(1)*exp(oo_.steady_state(7)))/(.01*exp(oo_.steady_state(5)));
ip_ig_mul(i)= (ip_eps_ig(1)*exp(oo_.steady_state(3)))/(.01*exp(oo_.steady_state(5)));
impact_cg_mul(i,:) = [gdp_cg_mul(i) cp_cg_mul(i) cpr_cg_mul(i) cpnr_cg_mul(i) ip_cg_mul(i)];
impact_ig_mul(i,:) =  [gdp_ig_mul(i) cp_ig_mul(i) cpr_ig_mul(i) cpnr_ig_mul(i) ip_ig_mul(i)];


% Calculation of PV of Five Year Multiplier


rds = rd_eps_cg*exp(oo_.steady_state(8)); % Product of interst ratesover time
rds1(1)=1;
for ii=2:50
rds1(ii) = rds(ii)*rds1(ii-1);
end

eefec_dis (i,:)=   rds1; % Effective Discounting


% Present Values of Macro-Aggregates Over Five Years

for kj=1:5
    
PVY_CG5(kj)= (gdp_eps_cg(kj)*exp(oo_.steady_state(1)))/ rds1(kj);
PVCP_CG5(kj) = (cp_eps_cg(kj)*exp(oo_.steady_state(2)))/ rds1(kj);
PVCPR_CG5(kj) = (cpr_eps_cg(kj)*exp(oo_.steady_state(6)))/ rds1(kj);
PVCPNR_CG5(kj) = (cpnr_eps_cg(kj)*exp(oo_.steady_state(7)))/ rds1(kj);
PVIP_CG5(kj) = (ip_eps_cg(kj)*exp(oo_.steady_state(3)))/ rds1(kj);
PVCG_CG5(kj) = (cg_eps_cg(kj)*exp(oo_.steady_state(4)))/ rds1(kj);
PVGDP_IG5(kj)= (gdp_eps_ig(kj)*exp(oo_.steady_state(1)))/ rds1(kj); 
PVCP_IG5(kj)= (cp_eps_ig(kj)*exp(oo_.steady_state(2)))/ rds1(kj); 
PVCPR_IG5(kj)= (cpr_eps_ig(kj)*exp(oo_.steady_state(6)))/ rds1(kj); 
PVCPNR_IG5(kj)= (cpnr_eps_ig(kj)*exp(oo_.steady_state(7)))/ rds1(kj); 
PVINVP_IG5(kj)= (ip_eps_ig(kj)*exp(oo_.steady_state(3)))/ rds1(kj); 
PVIG_IG5(kj) = (ig_eps_ig(kj)*exp(oo_.steady_state(5)))/ rds1(kj);
end

% Present Value Multipliers (5 Years)

pv_gdp_cg_mul5(i) = sum(PVY_CG5)/sum(PVCG_CG5);
pv_cp_cg_mul5(i) = sum(PVCP_CG5)/sum(PVCG_CG5);
pv_cpr_cg_mul5(i) = sum(PVCPR_CG5)/sum(PVCG_CG5);
pv_cpnr_cg_mul5(i) = sum(PVCPNR_CG5)/sum(PVCG_CG5);
pv_ip_cg_mul5(i) = sum(PVIP_CG5)/sum(PVCG_CG5);
pv_gdp_ig_mul5(i) = sum(PVGDP_IG5)/sum(PVIG_IG5);
pv_cp_ig_mul5(i) = sum(PVCP_IG5)/sum(PVIG_IG5);
pv_cpr_ig_mul5(i) = sum(PVCPR_IG5)/sum(PVIG_IG5);
pv_cpnr_ig_mul5(i) = sum(PVCPNR_IG5)/sum(PVIG_IG5);
pv_ip_ig_mul5(i) = sum(PVINVP_IG5)/sum(PVIG_IG5);


pv_cg_mul5(i,:) = [pv_gdp_cg_mul5(i) pv_cp_cg_mul5(i) pv_cpr_cg_mul5(i) pv_cpnr_cg_mul5(i) pv_ip_cg_mul5(i)];
pv_ig_mul5(i,:) =  [pv_gdp_ig_mul5(i) pv_cp_ig_mul5(i) pv_cpr_ig_mul5(i) pv_cpnr_ig_mul5(i) pv_ip_ig_mul5(i)];


% Present Values of Macro-Aggregates Over Fifteen Years

for kj1=1:15
    
PVY_CG15(kj1)= (gdp_eps_cg(kj1)*exp(oo_.steady_state(1)))/ rds1(kj1);
PVCP_CG15(kj1) = (cp_eps_cg(kj1)*exp(oo_.steady_state(2)))/ rds1(kj1);
PVCPR_CG15(kj1) = (cpr_eps_cg(kj1)*exp(oo_.steady_state(6)))/ rds1(kj1);
PVCPNR_CG15(kj1) = (cpnr_eps_cg(kj1)*exp(oo_.steady_state(7)))/ rds1(kj1);
PVIP_CG15(kj1) = (ip_eps_cg(kj1)*exp(oo_.steady_state(3)))/ rds1(kj1);
PVCG_CG15(kj1) = (cg_eps_cg(kj1)*exp(oo_.steady_state(4)))/ rds1(kj1);
PVY_IG15(kj1)= (gdp_eps_ig(kj1)*exp(oo_.steady_state(1)))/ rds1(kj1); 
PVCP_IG15(kj1)= (cp_eps_ig(kj1)*exp(oo_.steady_state(2)))/ rds1(kj1); 
PVCPR_IG15(kj1)= (cpr_eps_ig(kj1)*exp(oo_.steady_state(6)))/ rds1(kj1); 
PVCPNR_IG15(kj1)= (cpnr_eps_ig(kj1)*exp(oo_.steady_state(7)))/ rds1(kj1); 
PVIP_IG15(kj1)= (ip_eps_ig(kj1)*exp(oo_.steady_state(3)))/ rds1(kj1); 
PVIG_IG15(kj1) = (ig_eps_ig(kj1)*exp(oo_.steady_state(5)))/ rds1(kj1);

end

% Present Value Multipliers (15 Years)

pv_gdp_cg_mul15(i) = sum(PVY_CG15)/sum(PVCG_CG15);
pv_cp_cg_mul15(i) = sum(PVCP_CG15)/sum(PVCG_CG15);
pv_cpr_cg_mul15(i) = sum(PVCPR_CG15)/sum(PVCG_CG15);
pv_cpnr_cg_mul15(i) = sum(PVCPNR_CG15)/sum(PVCG_CG15);
pv_ip_cg_mul15(i) = sum(PVIP_CG15)/sum(PVCG_CG15);
pv_gdp_ig_mul15(i) = sum(PVY_IG15)/sum(PVIG_IG15);
pv_cp_ig_mul15(i) = sum(PVCP_IG15)/sum(PVIG_IG15);
pv_cpr_ig_mul15(i) = sum(PVCPR_IG15)/sum(PVIG_IG15);
pv_cpnr_ig_mul15(i) = sum(PVCPNR_IG15)/sum(PVIG_IG15);
pv_ip_ig_mul15(i) = sum(PVIP_IG15)/sum(PVIG_IG15);


pv_cg_mul15(i,:) = [pv_gdp_cg_mul15(i) pv_cp_cg_mul15(i) pv_cpr_cg_mul15(i) pv_cpnr_cg_mul15(i) pv_ip_cg_mul15(i)];
pv_ig_mull5(i,:) =  [pv_gdp_ig_mul15(i) pv_cp_ig_mul15(i) pv_cpr_ig_mul15(i) pv_cpnr_ig_mul5(i) pv_ip_ig_mul15(i)];





end

% Calculation of moments of fiscal multipliers

mean_impact_cg_mul = mean(impact_cg_mul, 1);
impact_cg_mul90 = prctile(impact_cg_mul, 90, 1);
impact_cg_mul10 = prctile(impact_cg_mul, 10,1);
max_impact_cg_mul = max(impact_cg_mul);
min_impact_cg_mul = min(impact_cg_mul);

mean_impact_ig_mul = mean(impact_ig_mul, 1);
impact_ig_mul90 = prctile(impact_ig_mul, 90, 1);
impact_ig_mul10 = prctile(impact_ig_mul, 90, 1);
max_impact_ig_mul = max(impact_ig_mul);
min_impact_ig_mul = min(impact_ig_mul);

mean_pv_cg_mul5 = mean(pv_cg_mul5, 1);
pv_cg_mul590 = prctile(pv_cg_mul15, 90, 1);
pv_cg_mul510 = prctile(pv_cg_mul15, 10,1);
max_pv_cg_mul5 = max(pv_cg_mul15);
min_pv_cg_mul5 = min(pv_cg_mul15);

mean_pv_ig_mul5 = mean(pv_ig_mul5, 1);
pv_ig_mul590 = prctile(pv_ig_mul5, 90, 1);
pv_ig_mul510 = prctile(pv_ig_mul5, 90, 1);
max_pv_ig_mul5 = max(pv_ig_mul5);
min_pv_ig_mul5 = min(pv_ig_mul5);

mean_pv_cg_mul15 = mean(pv_cg_mul15, 1);
pv_cg_mul1590 = prctile(pv_cg_mul15, 90, 1);
pv_cg_mul1510 = prctile(pv_cg_mul15, 10,1);
max_pv_cg_mul15 = max(pv_cg_mul15);
min_pv_cg_mul15 = min(pv_cg_mul15);

mean_pv_invg_mul15 = mean(pv_ig_mull5, 1);
pv_ig_mul1590 = prctile(pv_ig_mull5 , 90, 1);
pv_ig_mul1510 = prctile(pv_ig_mull5, 10, 1);
max_pv_ig_mul15 = max(pv_ig_mull5);
min_pv_ig_mul15 = min(pv_ig_mull5);

% Impact Multiplier

impact_cg_mult = [mean_impact_cg_mul; impact_cg_mul10; impact_cg_mul90; min_impact_cg_mul; max_impact_cg_mul ]; 
impact_invg_mult = [mean_impact_ig_mul; impact_ig_mul10; impact_ig_mul90; min_impact_ig_mul; max_impact_ig_mul ]; 

% PV 5 Years Multiplier

pv_cg_mult5 = [mean_pv_cg_mul5; pv_cg_mul510; pv_cg_mul590; min_pv_cg_mul5; max_pv_cg_mul5 ]; 
pv_invg_mult5 = [mean_pv_ig_mul5; pv_ig_mul510; pv_ig_mul590; min_pv_ig_mul5; max_pv_ig_mul5 ]; 

% PV 15 Years Multiplier

pv_cg_mult15 = [mean_pv_cg_mul15; pv_cg_mul1510; pv_cg_mul1590; min_pv_cg_mul15 ; max_pv_cg_mul15 ]; 
pv_invg_mult15 = [mean_pv_invg_mul15; pv_ig_mul1510; pv_ig_mul1590; min_pv_ig_mul15; max_pv_ig_mul15 ]; 



save calculated_multiplier1.mat