//* ******  Gali Monacelli (2016) * ******//


/*  Here there are two monetary regimes: e = 0, pi_H = 0 */


%----------------------------------------------
% Define one of two cases:
% - Inflation targeting : pi_H = 0 for all t        corresponds to CU = 0
% - Currency Union: e = 0 for all t
@#define CU = 0
%----------------------------------------------

%----------------------------------------------
% Preamble
%----------------------------------------------


var 
y           ${y}$       
c           
n
a
z
i
w
y_n         % natural output
c_n
n_n
n_gap       ${\tilde{n}}$
n_e         ${n^{e}}$
omega       % omega = w - p , real wage
omega_n
y_gap
c_gap
omega_gap
p
p_H
pi
pi_H
pi_w
s
e
s_n
s_gap
i_star
z_1_star
z_2_star
tau
r
i_ann
r_ann
//33 variables
;

varexo eps_z, eps_z_1, eps_z_2, eps_a, eps_tau;

parameters alpha beta phi 
//sigma 
upsilon rho_a rho_z rho_z_1 rho_z_2 epsilon_p epsilon_w theta_p theta_w rho_tau
// rho is a constant and thus omitted
lambda_p
lambda_w
;
%----------------------------------------------
% Parameterization: Follows GM(2016)!
%----------------------------------------------
alpha = 0.26;
beta = 0.99;
phi = 2.2;
upsilon = 0.3;
rho_a = 0.9;
rho_z = 0.9;
rho_z_1 = 0.9;
rho_z_2 = 0.9;
//sigma = 1;
epsilon_p = 3.8;
epsilon_w = 4.3;
theta_p = 0.8;
theta_w = 0.8;
rho_tau = 0.9;
lambda_p = (((1-theta_p)*(1-beta*theta_p))/theta_p)*((1-alpha)/(1-alpha + alpha*epsilon_p));
lambda_w = ((1-theta_w)*(1-beta*theta_w))/(theta_w*(1+epsilon_w*phi));




%----------------------------------------------
% Model
%----------------------------------------------

model;


@#if CU == 0
    pi_H = 0;
    @#else
      e = 0;
@#endif

// AD block
y = (1-upsilon)*c + upsilon*(2-upsilon)*s + upsilon*z_1_star;
c = (z - z_2_star) + (1-upsilon)*s;
c = c(+1) - (i - pi(+1)) + (1-rho_z)*z; // here there are different specifications... pi
e = s + p_H;
n = (1/(1-alpha))*(y-a);

// AS block
pi_H = beta*pi_H(+1) + ((lambda_p*alpha)/(1-alpha))*y_gap + lambda_p*omega_gap + lambda_p*upsilon*s_gap + lambda_p*tau;
pi_H = p_H - p_H(-1);
pi = p - p(-1);
p = p_H + upsilon*s;
pi_w = beta*pi_w(+1) + ((lambda_w*phi)/(1-alpha))*y_gap + lambda_w*c_gap - lambda_w*omega_gap;
pi_w = w - w(-1);
omega = w - p;

// gaps

c_gap = c - c_n;
y_gap = y - y_n;
omega_gap = omega - omega_n;
s_gap = s - s_n;
n_gap = n - n_e;

// natural values

n_n = (upsilon/(1+phi))*(z_1_star + z_2_star - z) - (1/(1+phi))*tau;
y_n = a + (1-alpha)*n_n;
s_n = a - (alpha + phi)*n_n - z + z_2_star - tau;
c_n = (z - z_2_star) + (1-upsilon)*s_n;
omega_n = a - alpha*n_n - tau - upsilon*s_n;
n_e = (upsilon/(1+phi))*(z_1_star + z_2_star - z);


r = i - pi_H(+1);

// exogenous processes

a = rho_a*a(-1) - eps_a;
z = rho_z*z(-1) - eps_z;                %specify with a minus!?
z_1_star = rho_z_1*z_1_star(-1) - eps_z_1;
z_2_star = rho_z_2*z_2_star(-1) - eps_z_2;
i_star =(1-rho_z_2)*z_2_star;          
tau = rho_tau*tau(-1) - eps_tau;        %specify with a minus!


// annualized values

i_ann = 4*i;
r_ann = 4*r;

//33 equations

end;


steady;
%----------------------------------------------
% Define shock variances
%----------------------------------------------

shocks;
var eps_a = 0.1;
end;


stoch_simul(periods=100000,order = 1, irf = 0, noprint)n_gap pi pi_w; % increase the periods for precision










%----------------------------------------------
% Welfare computation
%----------------------------------------------

% Get variable positions in variable list:
pi_pos = strmatch('pi',var_list_, 'exact');
pi_w_pos = strmatch('pi_w', var_list_, 'exact');
n_gap_pos = strmatch('n_gap', var_list_, 'exact');




% Loop over the interval theta_w element[0,1]:
space = 0.01:0.01:0.99;

% variances and components
    variance.pi = zeros(length(space),1);
    variance.pi_w = zeros(length(space),1);
    variance.n_gap = zeros(length(space),1);
    comp_n_gap = zeros(length(space),1);
    comp_pi = zeros(length(space),1);
    comp_pi_w = zeros(length(space),1);


for ii = 1:length(space)
    set_param_value('theta_w',space(ii));

    par.lambda_w = ((1-theta_w)*(1-beta*theta_w))/(theta_w*(1+epsilon_w*phi));
    set_param_value('lambda_w', par.lambda_w);
//    info = stoch_simul(var_list_);

stoch_simul(periods=100000,order = 1, irf = 0, noprint)n_gap pi pi_w; % increase the periods for precision
  


    variance.pi(ii) = oo_.var(pi_pos, pi_pos);
    variance.pi_w(ii) = oo_.var(pi_w_pos, pi_w_pos);
    variance.n_gap(ii) = oo_.var(n_gap_pos, n_gap_pos);




    L(ii) = ((1-upsilon)/2)*((1+phi)*variance.n_gap(ii) + (epsilon_p/(lambda_p*(1-alpha)))*variance.pi(ii) + (epsilon_w/lambda_w)*variance.pi_w(ii));


    comp_n_gap(ii) = ((1-upsilon)/2)*(1+phi)*variance.n_gap(ii);
    comp_pi(ii) = ((1-upsilon)/2)*(epsilon_p/(lambda_p*(1-alpha)))*variance.pi(ii);
    comp_pi_w(ii) = ((1-upsilon)/2)*(1+phi)*(epsilon_w/lambda_w)*variance.pi_w(ii);




end

    variance_matrix = [variance.pi , variance.pi_w, variance.n_gap];
    comp_matrix = [comp_pi, comp_pi_w, comp_n_gap];





 NN = L./L(80);
 xx=linspace(0,0.99,99);
plot(xx,NN)