

clc;
clear;

ss.Rlfs             =            1.033^0.25;            
par.Rre             =            1.016^0.25;
ss.Rbs              =            1.056^0.25;
ss.As               =            1;
ss.omegas           =            1;
par.epsilon         =            6;
ss.ws               =            ss.As*ss.omegas*(par.epsilon - 1)/par.epsilon;
par.zeta            =            0.05;
par.u               =            0.052;
ss.phis             =            0.14;
ss.Rs               =            (1 - ss.phis)*ss.Rlfs + ss.phis*par.Rre;
ss.tauls            =            0.05;
ss.uis              =            0.241*ss.ws;
ss.Tas              =            0;
ss.Tbs              =            0;
ss.Tss              =            0;
ss.Pis              =            1.016^0.25;
par.betas           =            ss.Pis/ss.Rs;
par.betab           =            0.98*par.betas;
ss.Lamss            =            par.betas;
ss.Lambs            =            par.betab;
ss.Qs               =            par.betas/ss.Pis;
ss.pas              =            1.1*ss.ws;
ss.Nas              =            0.35*(1 - par.u);
ss.Nns              =            1 - par.u - ss.Nas;
ss.Yns              =            ss.As*ss.omegas*ss.Nns;
ss.Ys               =            ss.Yns; 
ss.Js               =            ss.Nas;
par.sigma           =            1;
par.sigmaa          =            1;
ss.Fbs              =            0.019;
ss.Qbs              =            (1 - ss.Fbs)/ss.Rbs;
par.chi             =            0.24;
ss.bbs              =            0.26*((1 - par.u)*ss.ws*(1 - ss.tauls) + par.u*ss.uis);
par.Gamma           =            ss.Qbs*ss.bbs;
ss.lambdas          =            1 - ss.Lambs*ss.Rbs/ss.Pis;
ss.epses            =            ss.bbs/ss.Pis - ss.ws*(1 - ss.tauls) - ss.Tbs;
ss.epsus            =            ss.bbs/ss.Pis - ss.uis - ss.Tbs;
par.sigmae          =            ss.epses/norminv(ss.Fbs/(1 + 4*par.u));
par.sigmau          =            ss.epsus/norminv(5*ss.Fbs/(1 + 4*par.u));
ss.Cbs              =            (1 - par.u)*ss.ws*(1 - ss.tauls) + par.u*ss.uis + ss.Qbs*ss.bbs - ss.bbs*(1 - ss.Fbs)/ss.Pis;

par.psigma          =            0.57;
ss.cbars            =            par.psigma;
par.sigmak          =            3.2;
ss.Has              =            normcdf(log(ss.cbars)/par.sigmak + par.sigmak/2);
ss.Psis             =            normcdf(log(ss.cbars)/par.sigmak - par.sigmak/2);
% fHas   = @(x) normcdf(log(ss.cbars)/x + x/2);
% fPsis = @(x) normcdf(log(ss.cbars)/x - x/2); 
% f1 = @(x) ss.pas*ss.As  - ss.ws + par.betas*(fHas(x)*ss.cbars - fPsis(x)) - ss.cbars;
% fsolveopts = optimoptions('fsolve','OptimalityTolerance',1e-10);
% out                 =            fsolve(f1, 0, fsolveopts);
% par.sigmak          =            out;
% ss.Psis             =            fPsis(out);
% ss.Has              =            fHas(out);

ss.vs               =            (1 - ss.Has)*ss.Js;
par.psi             =            1;
ss.Yas              =            ss.As*ss.Nas;
ss.Cas              =            ss.Yas/(1 - par.chi);
ss.GDPs             =            ss.Yns + ss.pas*ss.Yas;
ss.Gs               =            0.167*ss.GDPs;
ss.Css              =            (ss.Yns - par.chi*ss.Cbs -ss.Gs - ss.Psis*ss.Js)*(1 - par.chi)^(-1);
ss.alphas           =            ss.pas*ss.Cas^par.sigmaa/ss.Css^par.sigma;
par.xi              =            ss.ws/(ss.As*(ss.Nas + ss.omegas*ss.Nns)^par.zeta);
calvo               =            0.75; 
par.peta            =            (par.epsilon-1)*calvo/((1 - calvo)*(1-par.betas*calvo));
lev                 =            1.5;
ss.res              =            ss.phis*(lev - 1)*par.chi*ss.Qbs*ss.bbs/(ss.phis + lev);
ss.lfs              =            ss.res;
ss.ds               =            ss.res/((1 - par.chi)*ss.phis*ss.Qs);
ss.Es               =            (par.chi*ss.Qbs*ss.bbs + ss.res)/lev;
ss.Phis             =            par.chi*ss.Qbs*ss.bbs/ss.Es;
par.theta           =            0.9;
Theta               =            ss.Pis^(-1)*ss.Lamss*(ss.Phis*(ss.Rbs - ss.Rlfs) + ss.Rlfs);
par.kappa           =            Theta*(1 - par.theta)/((1 - par.theta*Theta)*ss.Phis); 
ss.Omegas           =            1 - par.theta + par.theta*par.kappa*ss.Phis;
ss.mus              =            1 - ss.Omegas*ss.Pis^(-1)*ss.Lamss*ss.Rlfs/(par.kappa*ss.Phis);
ss.etas             =            ss.Lamss*ss.Pis^(-1)*ss.Omegas*(ss.Rlfs - par.Rre)/(1 - ss.mus);
par.pomega          =            ss.Es - ss.Pis^(-1)*par.theta*((1 - ss.Fbs)*par.chi*ss.bbs - (1 - par.chi)*ss.ds + par.Rre*ss.res - ss.Rlfs*ss.lfs);
ss.Ps               =            ss.Yns - ss.ws*ss.Nns + (ss.pas - ss.ws)*ss.Nas - ss.Psis*ss.Js +...
    ss.Pis^(-1)*(1 - par.theta)*((1 - ss.Fbs)*par.chi*ss.bbs - (1 - par.chi)*ss.ds + par.Rre*ss.res - ss.Rlfs*ss.lfs);
ss.bgs              =            1.55*ss.GDPs/((1 - par.chi)*ss.Qs);
ss.Tcbs             =            ss.Pis^(-1)*(ss.Rlfs*ss.lfs - par.Rre*ss.res);
par.tauk            =            ss.Ps^(-1)*(ss.Gs + (1 - par.chi)*ss.bgs/ss.Pis + par.u*ss.uis -...
    ss.tauls*ss.ws*(1 - par.u) - ss.Tcbs - (1 - par.chi)*ss.Qs*ss.bgs);
par.phiPi           =            2; 
par.phiu            =            0.25; 
par.phitau          =            0.33;
par.rholf           =            0.8;
par.slf             =            0.01/4;
par.rhophi          =            0.98;
par.sphi            =            0.01;
par.rhoA            =            0;
par.sA              =            1;
par.rhoalpha        =            0;
par.salpha          =            1;
par.rhoomega        =            0;
par.somega          =            1;
par.rhoa            =            0;
par.sa              =            1;
par.rhob            =            0;
par.sb              =            1;
par.rhoui           =            0;
par.sui             =            1;
par.rhol            =            0;
par.sl              =            1;
par.rhoG            =            0;
par.sG              =            1;

save calibrations par ss;

shock_A = zeros(1,2);
shock_alpha = zeros(1,2);
shock_omega = zeros(1,1);
shock_lf = zeros(1,1);
shock_phi = zeros(1,1);
shock_G = 0.01*ones(1,1);
shock_ui = zeros(1,1);
shock_l = zeros(1,1);
shock_b = zeros(1,1);
shock_a = zeros(1,1);

save shocks shock_A shock_alpha shock_omega shock_lf shock_phi shock_G shock_ui shock_l shock_b shock_a;
