%=========================================================================
// Preamble
%=========================================================================

// Name of the endogenous variables
var  d, c, h, y, i, k, a, lambda,  tb, ca, risk, r, g, tau, dg, s, taxp, w, rk, sy, dy, dgy, primary, overall;  
 
// Name of the exogenous variable
varexo e, eg;                                    
                                             
// Name of the model parameters  
parameters  gamma, omega,                    // Preferences
            rho, sigmae,                     // Productivity shock (Persistence and volatility)
            delta,                           // Depreciation rate
            psi,                             // Elasticity of the discount factor
            psib,                            // Elasticity of portfolio adjustment cost (households)
            psid,                            // Elasticity of portfolio adjustment cost (government)  
            alpha,                           // Technology
            phi,                             // capital adjustment costs
            beta,                            // subjective dicount factor
            r_w,                             // world interst rate
            rhogg,                           // ar(1) parameter on government expenditure
            rhogd,                           // debt-parameter on government expenditure
            theta,                           // debt-parameter on taxes
            eta,                             // debt-deficit preference parameter
            h_bar,                           // Steady state level of work
            tau_bar,                         // Steady state tax-rate
            y_bar,                           // Steady state level of y
            k_bar,                           // Steady state level of k
            w_bar,                           // Steady state level of w 
            r_bar,                           // Steady state level of r   
            rk_bar,                          // Steady state level of rk  
            d_bar,                           // Steady state level of d  
            dg_bar,                          // Steady state level of dg
            gy_bar,                          // Steady state level of g/y
            dgy_bar,                         // Steady state level of gov debt to GDP
            g_bar,                           // Steady state level of g
            s_bar;                           // STeady state level of niip
            
%=========================================================================
// Calibration
%=========================================================================
        
psib        = 0.00074;                       // From Uribe (2003)
psid        = 0.00074;
gamma       = 2;


psi         = 0.000742*100;        
  
alpha       = 0.32;
phi         = 0.084;

r_w         = 0.04;

delta       = 0.1;
rho         = 0.42;

sigmae      = 0.01;
sigmaeg     = 0.01;

beta        = 1/(1+r_w);
eta         = 1;

rhogg       = 0.42;                     
rhogd       = -0.2;                 

theta       = 0.8;                   
h_bar       = 0.3;                     

tau_bar     = 0.3;                  
gy_bar      = 0.22;

omega       = 0.62;

%=========================================================================
//Assisting Steady State values
%=========================================================================

// steady state value of hours
hss         = h_bar; 

// steady state tax rate on output
tauss       = tau_bar;

// steady state capital
r_kss       = (1/beta-(1-delta)-delta*tauss)/(1-tauss);
kss         = hss/(((r_kss)/alpha)^(1/(1-alpha)));
k_bar       = kss;

// steady state of w and rk
wss         = (1-alpha)*kss^(alpha)*hss^(-alpha);
rkss        = alpha*kss^(alpha-1)*hss^(1-alpha);
w_bar       = wss; 
rk_bar      = rkss;
r_bar       = r_w;

// steady state investment
iss         = delta*kss;                                                     

// steady state output
yss         = (kss^alpha)*(hss^(1-alpha));  
y_bar       = yss;                              

// steady state consumption
css         = (omega/(1-omega))*(1-hss)*(1-tauss)*(1-alpha)*yss;

// steady state financial assets 
dss         = (1/r_w)*((1-tau_bar)*yss+tau_bar*delta*kss-css-iss);         
d_bar       = dss;

// steady state goverment consumption
gss         = gy_bar*yss;
g_bar       = gss;

// steady state trade balance
tbss        = yss-css-iss-gss;

// steady state goverment debt-stock
dgss        = (tauss*(yss-delta*kss)-gss)/r_w;
dgy_bar     = dgss/yss;
dg_bar      = dgss;

// steady state NIIP
sss         = dgss+dss;
s_bar       = sss;

// steady state of debts/GDP
dyss        =dss/yss;
dgyss       =dgss/yss;

// steady state levels of primary and overall fiscal balances
primaryss   = g_bar-tau_bar*(w_bar*h_bar+(rk_bar-delta)*k_bar);
overallss   = r_w*d_bar+g_bar-tau_bar*(w_bar*h_bar+(rk_bar-delta)*k_bar);
primary_bar = primaryss;
overall_bar = overallss;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// Description of the model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

model;

// Interest rate elastic to total foreign net position (ratio of y)
%risk             = psi*(exp(s/exp(y)-s_bar/y_bar)-1);

// Interest rate elastic to a linear combination of total foreign debt and overall balance
risk             = eta*(psi*(exp(s/exp(y)-s_bar/y_bar)-1))+(1-eta)*(psi*(exp((exp(r(-1))*d(-1)+exp(g)-exp(tau)*(exp(w)*exp(h)+(rk-delta)*exp(k)))/exp(y)-(r_w*d_bar+g_bar-tau_bar*(w_bar*h_bar+(rk_bar-delta)*k_bar))/y_bar)-1));

// Interest rate elastic to total foreign net position
%risk             = psi*(exp(s/-s_bar)-1);

// Interest rate elastic to public foreign net position
%risk             = psi*(exp(dg/exp(y)-dg_bar/y_bar)-1);

// Interest rate elastic to overall fiscal balance
%risk             = psi*(exp((exp(r(-1))*d(-1)+exp(g)-exp(tau)*(exp(w)*exp(h)+(rk-delta)*exp(k)))/exp(y)-(r_w*d_bar+g_bar-tau_bar*(w_bar*h_bar+(rk_bar-delta)*k_bar))/y_bar)-1);


// Interest rate elastic to primary fiscal balance
%risk             = psi*(exp((exp(g)-exp(tau)*(exp(w)*exp(h)+(rk-delta)*exp(k)))/exp(y)-(g_bar-tau_bar*(w_bar*h_bar+(rk_bar-delta)*k_bar))/y_bar)-1);

// Primary and overall fiscal balance
primary = exp(g)-exp(tau)*(exp(w)*exp(h)+(rk-delta)*exp(k)); 
overall = exp(r(-1))*d(-1)+exp(g)-exp(tau)*(exp(w)*exp(h)+(rk-delta)*exp(k));

// Interest rate equation
exp(r)           = r_w+risk;

// Tax rate evolution
taxp             = theta*(dg(-1)/exp(y(-1))-dgy_bar);
exp(tau)         = tau_bar+taxp;

// FOC for profit maximization (wrt k and h)
rk               = exp(a)*alpha*exp(k)^(alpha-1)*exp(h)^(1-alpha);
exp(w)           = exp(a)*(1-alpha)*exp(k)^(alpha)*exp(h)^(-alpha);

// Goverment consumption evolution 
exp(g)           = g_bar+rhogg*(exp(g(-1))-g_bar)+rhogd*(dg(-1)/exp(y(-1))-dgy_bar)+eg;   

// Goverment debt evolution
dg               = (1+exp(r(-1)))*dg(-1)+exp(g)-exp(tau)*(exp(w)*exp(h)+(rk-delta)*exp(k))+(psid/2)*(dg-dg_bar)^2;

// Budget constraint, Market clearing condition
d                = (1+exp(r(-1)))*d(-1)-(1-exp(tau))*(exp(w)*exp(h)+(rk-delta)*exp(k))+exp(c)+exp(k)-exp(k(-1))+(phi/2)*(exp(k)-exp(k(-1)))^2+(psib/2)*(d-d_bar)^2;

// Production function
exp(y)           = exp(a)*(exp(k))^(alpha)*(exp(h)^(1-alpha));
            
// Evolution of capital stock                                          
exp(k)           = exp(i)+(1-delta)*exp(k(-1))-(phi/2)*(exp(k)-exp(k(-1)))^2; 

// FOC wrt d (Euler Equation)                        
exp(lambda)*(1-psib*(d-d_bar))= beta*(1+exp(r))*exp(lambda(+1)); 

// FOC wrt c 
((exp(c)^omega)*((1-exp(h))^(1-omega)))^(-gamma)*omega*exp(c)^(omega-1)*(1-exp(h))^(1-omega)=exp(lambda);

// FOC wrt h 
((exp(c)^omega)*((1-exp(h))^(1-omega)))^(-gamma)*(1-omega)*exp(c)^(omega)*(1-exp(h))^(-omega)=exp(lambda)*(1-exp(tau))*exp(w);
                                  
// FOC wrt k                 
exp(lambda)*(1+phi*(exp(k)-exp(k(-1)))) = beta*exp(lambda(+1))*((1-exp(tau(+1)))*(rk(1)-delta)+1+phi*(exp(k(+1))-exp(k)));

// technology shock         
a      = rho*a(-1)+e; 

// niip
s      = d+dg;

// trade balance
tb     = 1-((exp(c)+exp(i)+exp(g))/exp(y));

// current account
ca     = (1/exp(y))*(-s+s(-1));         

// extra variables

sy      = s/exp(y);
dy      = d/exp(y); 
dgy     = dg/exp(y);
                          
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


initval;
risk  = 0;
taxp  = 0;
r     = log((1-beta)/beta);
tau   = log(tauss);
g     = log(gss);
dg    = dgss;
a     = 0;
e     = 0;
d     = dss;
h     = log(hss);
w     = log(wss);
rk    = rkss;
k     = log(kss);
y     = log(yss);
c     = log(css);
i     = log(iss);
s     = sss;
tb    = 1-((exp(c)+exp(i)+exp(g))/exp(y));
ca    = 0;
lambda= log(((exp(c)^omega)*((1-exp(h))^(1-omega)))^(-gamma)*omega*exp(c)^(omega-1)*(1-exp(h))^(1-omega));
sy    = sss/yss;
dy    = dyss;
dgy   = dgyss;
primary = primaryss;
overall = overallss;

end;
resid;
steady; 
check;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


shocks;

var e; 
%periods 1:2 3 4:1000;
%values 0.0 .01 0.0;
stderr 0.01;

%var eg; 
%stderr 0.01;
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M=50;
grid=linspace(0,1,M);
result=zeros(1,M);

%for j = 1:M
%eta = grid(1,j);
stoch_simul(order=1, periods=1000, irf=10)  y;
%result(1,j)=risk(1);
%end

%plot(grid,result)

stoch_simul(order=1, periods=1000, irf=10)  y c g i h k d tb ca risk sy s dg d dy dgy primary overall;