var
r so cc p_c pi_c wc wo ytt n_to n_tn pio c_o mco yo ao po u g y a_t n_t mc_t w pi c_f c_h p_f c p_h p pi_h pi_f e s p_star q c_star pi_star r_star;

varexo
eps_r eps_e eps_ao eps_o eps_a eps_cstar eps_pistar eps_rstar eps_g eps_u;


// Declaration of parameters
//-------------------------------------

parameters
rho_ao sigma_o sigma_ao hh jj j b kapa lambda ee rhoo rhou sigma_g sigma_u rhog rho_a sigma rho_rstar phi eta alpha beta theta rho_e rho_r psi_y psi_pi psi_q gamma rho_cstar rho_pistar sigma_r sigma_a sigma_cstar sigma_pistar sigma_ner;


// Calibration of parameter values
//-------------------------------------
sigma_o=0.015;
sigma_ao=0.015;
hh=0.3;
jj = 0.03;
j = 0.3;
kapa=0.9;
alpha=0.2;
ee = 0.5;
rhoo = 0.5;
beta    = 0.99;                         // Rate of time preference (subjective discount factor)
rho_e   = 0.99;                         // Nominal exchange rate persistence
rho_cstar = 0.7;                       // World consumption persistence
rho_pistar = 0.3217;                      // World inflation persistence
rhou = 0.5;
rhog = 0.45;                  
// Setting initial values of estimated parameters using prior mean
//----------------------------------------------------------------

b = 0.5;
sigma   = 1.5;                         // Inverse of the elasticity of substitution between consumption and labor  rho
eta     = 6;                         // Elasticity of substitution between domestic tradable and imported goods theta
%gamma   = 1.5;                            // Elasticity of substitution between goods produced in different foreign countries
phi     = 1;                            // Inverse labor elasticity eta
rho_a  = 0.4;                        // Productivity persistence for domestic tradable firms
rho_ao = 0.4;
theta  = 0.39;                       // Calvo parameter for tradables firms (probability of keeping prices unchanged)
rho_r = 0.5;                          // Taylor rule interest rate persistence
psi_y = 100;                          // Taylor rule output
psi_pi = 100;                        // Taylor rule inflation
%psi_q = 2;                          // Taylor rule exchange rate
rho_rstar = 0.6;
sigma_r=0.1;                                                                                                    
sigma_a=0.0071;                                                                                                    
sigma_cstar=0.0129;
sigma_pistar=0.05;
sigma_ner=0.05;
sigma_u=0.0013;
sigma_g=0.0056;
%mu=0.5;
lambda=0.5;
// Declaration of the log-linearized model
//-----------------------------------------

model(linear);
%c(+1) - c = c_star(+1) - c_star + ((q(+1)-q)/(sigma)) - ((mu*b)/sigma);

//production function ok
ytt = (1-lambda)*y + lambda*yo;           //Eq 1 ok ρρ
y = a_t + (1-jj)*n_tn + j*yo;                             //Eq 2 ok ρρ
yo = ao + hh*n_to;                             //Eq 3 ok ρρ
n_tn = kapa*n_t;                            //Eq 25 ok ρρ
n_to = (1-kapa)*n_t;                        //Eq 25 ok ρρ
//productivity process ok ρρ
a_t = rho_a*a_t(-1) + eps_a ;               //Eq 4    ok      ρρ                             // 6. Domestic tradable firms productivity
ao = rho_ao*ao(-1) + eps_ao;                 //Eq 5 ok ρρ

//marginal cost ok ρρ
mc_t = (1-jj)*wc + jj*po - p_h -a_t;                        //Eq 8 ρρ
%mc_t = wc - p_h -a_t;                        //Eq 8 ρρ
mco = wo - p_h - ao;                         //Eq 9 ρρ
n_tn = -b*(wc - w) + n_t;                    // ok ρρ
n_to = -b*(wo - w) + n_t;                    //ok ρρ
//EE ok ρρ
c = -(1/sigma)*(phi*n_t - w + p - so) ;                                   // 10. Aggregate consumption

//government shock ok ρρ
g = rhog*g(-1) + eps_g ;

//risk sharing ok ρρ
c = c_star +  ((1-alpha)*q)/sigma + so/sigma;
%beta*b=b(-1)+y-c-((alpha)/(1-alpha))*q

//home non-oil consumption  ρρ
cc = -j*(p_c - p_h) + c_h;

//consumption of foreign goods ok ρρ
c_f = -eta*(p_f - p) + c ;  

//consumption of domestic goods ok ρρ
c_h = -eta*(p_h - p) + c ;  

//consumption of oil ok ρρ
c_o = -j*(po - p_h) + c_h ;

//oil inflation ok ok ρρ
pio = po - po(-1) ;

//Share of inflation MMMMMMAL
%pi_c = (1-alpha)*pi_h + alpha*pi_f;
pi_h = (1-lambda)*pi_c + lambda*pio;

//Phillips curve ok ρρ
pi_h = beta*pi_h(+1) + (((1 - theta)*(1 - beta*theta))/theta)*mc_t  + (((1 - theta)*(1 - beta*theta))/theta)*u;

//Markup shock ok ρρ
u=rhou*u(-1)+eps_u;

//Imported inflation ok 3.53 ρρ
pi_f = e - e(-1) + pi_star ;                                                       // 23. Foreign tradables inflation

//definition of inflation OK ρρ
p = p(-1) + pi ;

//cpi inflation MMMALLLLL
%pi = (1-lambda)*pi_c + lambda*pio;
%pi=(1-alpha)*pi_h + alpha*pi_f;
pi=(1-lambda)*pi_c + lambda*pio + alpha*(pi_f - pi_h);

//definition of home non-oil inflation ρρ
pi_c = p_c - p_c(-1);



//definition of home inflation ok ρρ
p_h = p_h(-1) + pi_h ;

//definition of foreign inflation ok ρρ
p_f = p_f(-1) + pi_f ;

//nominal exchange rate ok ρρ
e = rho_e*e(-1) + eps_e ;                                                           // 36. Nominal exchange rate AR process    

//terms of trade ok ρρ
s = e + p_star - p_h ;

//real exchange rate OK ρρ
q = (1 - alpha)*s;  

//monetary policy ok ρρ
r = rho_r*r(-1) + (1 - rho_r)*(psi_y*ytt + psi_pi*pi_h) + eps_r ;             // 39. Monetary policy under flexible exchange rate
%e - e(-1)=0;

//Aggregate demand ok ρρ
y = alpha*c_star + eta*alpha*q + (1-alpha)*c ;

//stochastic process for foreign consumption ok ρρ
c_star = rho_cstar*c_star(-1) + eps_cstar ;                                         // 44. World consumption

//stochastic process for foreign inflation ok ρρ
pi_star = rho_pistar*pi_star(-1) + eps_pistar ;                                     // 45. World inflation

//definition of foreign inflation ok ρρ
p_star = p_star(-1) + pi_star ;                                                     // 46. World price level

//stochastic process for foreign interest rate ok ρρ
r_star = rho_rstar*r_star(-1) + eps_rstar ;                                         // 47. World interest rate

//Stochastic process for oil price ok ρρ
po=rhoo*po(-1)+eps_o;

//Demand for intermediate goods for oil production
%m = (1/(1-ee))*po + (1/(ee-1))*p_f  + (1/(1-ee))*s ;

//oil production
%yo = ee*m;

//Subsidy on consumption
so + c + p = po + yo;



end;

steady(solve_algo=0);

// Calibration of Shocks
//-------------------------------------
shocks;
var eps_ao     = sigma_ao^2;
var eps_o       =sigma_o^2;
var eps_r       = sigma_r^2;                                                                                                    
var eps_a      = sigma_a^2;                                                                                                    
var eps_cstar   = sigma_cstar^2;
var eps_pistar  = sigma_pistar^2;
var eps_e       = sigma_ner^2;
var eps_u       = sigma_u^2;
var eps_g       = sigma_g^2;
end;

nobs=5500;
set_dynare_seed(5001)

// Simulation to Generate Simulated Series \\
stoch_simul(drop=0,periods=5000,order=1,irf=20,nograph) n_t;


n_t_pos=strmatch('n_t',var_list_ ,'exact');

    psi_pi_vec=[1.5 1.5 5 1.5];
    psi_y_vec=[0.125 0 0 10];

    variance.n_t=NaN(1,length(psi_pi_vec));
    L=NaN(1,length(psi_pi_vec));
    for ii=1:length(psi_pi_vec)
        set_param_value('psi_pi',psi_pi_vec(ii));
        set_param_value('psi_y',psi_y_vec(ii));
        %[info, oo_, options_] = stoch_simul(M_, options_, oo_, var_list_); %loop over stoch_simul
        %if ~info(1)
            %read out current parameter values

            par.phi=M_.params(strmatch('phi',M_.param_names,'exact'));

            variance.n_t(ii)=oo_.var(n_t_pos,n_t_pos);
            L(ii)=-0.5*(1+par.phi)*variance.n_t(ii);
        %end
      
end

   //Print result

    labels={'psi_pi';'psi_y';'Welfare loss'};
    headers={' ';' ';' ';' '};
    values=[psi_pi_vec;psi_y_vec;L];
    options_.noprint=0;
    dyntable(options_,'Exchange rate peg',headers,labels,values,size(labels,2)+2,4,10)
    options_.noprint=1;