% dynare nonlinear_general_ti_calvo.mod
close all

maxit_=100;

ii = 1;

for ti = [1.0^0.25 1.02^0.25 1.04^0.25]


var w r C x psi phi N s i pi vf ;
parameters ALFAI ALFAPAI ALFAY sigma newphi gamma alpha beta epselon teta chi mu rho rhoA 
auxCss auxCss1 vfss vfss1 chi sigma dn alpha beta epselon gamma newphi gss gss1 rss teta 
phiss xss sss Css wss Nss psiss rss iss phiss1 xss1 sss1 Css1 wss1 Nss1 psiss1 iss1 Ass Ass1 ;
varexo a rgm intssv cssv;   // rgm is (gss) pi_bar

sigma=1;        % elasticity of marginal utility of consumption
newphi=1;       % elasticity of labour supply
gamma=0;        % real wage rigidity parameter
alpha=0;        % decreasing returns to scale parameter Yit = At*Nit^(1-alpha)
beta=0.99;      % subjective discount factor
epselon=10;     % Dixit-Stiglitz elasticity of substitution among goods
teta=0.75;      % Calvo parameter
chi=0;          % degree of price indexation
mu=0;           % degree of backward looking indexation 
rho=0;          % persistence of the inflation target shock
rhoA=0;         % persistence of technology shock

% TAYLOR RULE CALIBRATION
ALFAI=0;        % degree of interest rate inertia in Taylor rule
ALFAPAI=1.5;    % coefficient on inflation gap
ALFAY=0.5;      % coefficient on output gap



%Calibration of dn
dn=3^(newphi+sigma+alpha*(1-sigma))*(1-alpha)*(epselon-1)/epselon;

% initial ss inflation
gss0=ti;   // pi_bar is set in trend inflation


% final ss inflation
%gss1=(ti^4-0.04)^0.25;
gss1=1^0.25;


% DEFINE: pi = inflation rate; x = ratio between optimal resetted price and the price level = pi*/p ( this is small p_star in paper)

// let consider the case when inflation target (pibar=0)
Ass1=1;                                                   
rss=1/beta-1;                                                                                                                        // equa (1ss or iss) as pi_bar=1
iss1=((1+rss)*gss1)-1;
xss1=((1-teta*gss1^((epselon-1)*(1-chi)))/(1-teta))^(1/(1-epselon));                                                                 // equa (3ss or piss) 
sss1=((1-teta)*(xss1)^(-epselon/(1-alpha)))/(1-teta*gss1^(epselon*(1-chi)/(1-alpha)));                                               // equa (SSS or 8ss)
auxCss1=(1-alpha)*(epselon-1)*(1-teta*beta*gss1^(epselon*(1-chi)/(1-alpha)))/(epselon*dn*(1-teta*beta*gss1^((epselon-1)*(1-chi))));  // See apendix page 10 that define this variable
Css1=((xss1^((1-alpha+epselon*alpha)/(1-alpha)))*auxCss1*Ass1^((newphi+1)/(1-alpha))*sss1^(-newphi))^((1-alpha)/(newphi+sigma+alpha*(1-sigma))); // equation yss
Nss1=sss1*(Css1/Ass1)^(1/(1-alpha));                                                                           // equation (nss or 7ss)
wss1=dn*Nss1^newphi*Css1^sigma;                                                                                // equation (wss or 2ss)
phiss1=(Css1^(1-sigma))/(1-teta*beta*gss1^((epselon-1)*(1-chi)));                                              // equation (phiss or 6ss)
psiss1=(Ass1^(-1/(1-alpha))*Css1^((1/(1-alpha))-sigma)*wss1)/(1-teta*beta*gss1^(epselon*(1-chi)/(1-alpha)));   // equation (psiss or 5ss)
vfss1=(1-beta)^(-1)*(log(Css1)-dn*Nss1^(1+newphi)/(1+newphi));                                                 // the value function (utility function) as sigma = 1
%vfss1=(1-beta)^(-1)*(Css1^(1-sigma)/(1-sigma)-dn*Nss1^(1+newphi)/(1+newphi));                                 // the value function in the case sigma > 1

// let consider a general case as gss (pi_bar) take values : 0, 2, 4
Ass=1;
rss=1/beta-1;
iss=((1+rss)*gss0)-1;
xss=((1-teta*gss0^((epselon-1)*(1-chi)))/(1-teta))^(1/(1-epselon));
sss=((1-teta)*(xss)^(-epselon/(1-alpha)))/(1-teta*gss0^(epselon*(1-chi)/(1-alpha)));
auxCss=(1-alpha)*(epselon-1)*(1-teta*beta*gss0^(epselon*(1-chi)/(1-alpha)))/(epselon*dn*(1-teta*beta*gss0^((epselon-1)*(1-chi))));
Css=((xss^((1-alpha+epselon*alpha)/(1-alpha)))*auxCss*Ass^((newphi+1)/(1-alpha))*sss^(-newphi))^((1-alpha)/(newphi+sigma+alpha*(1-sigma)));
Nss=sss*(Css/Ass)^(1/(1-alpha));
wss=dn*Nss^newphi*Css^sigma;
phiss=(Css^(1-sigma))/(1-teta*beta*gss0^((epselon-1)*(1-chi))); 
psiss=(Ass^(-1/(1-alpha))*Css^((1/(1-alpha))-sigma)*wss)/(1-teta*beta*gss0^(epselon*(1-chi)/(1-alpha)));
vfss=(1-beta)^(-1)*(log(Css)-dn*Nss^(1+newphi)/(1+newphi));
%vfss=(1-beta)^(-1)*(Css^(1-sigma)/(1-sigma)-dn*Nss^(1+newphi)/(1+newphi));



// the non-linear system
model;
C(+1)^sigma=beta*(1+r)*C^sigma;
w=w(-1)^gamma*(dn*(N^newphi)*C^sigma)^(1-gamma);
1=(teta*rgm^((1-epselon)*chi*(1-mu))*pi(-1)^((1-epselon)*mu*chi)*pi^(epselon-1)+(1-teta)*x^(1-epselon));
x^(1+(epselon*alpha)/(1-alpha))*phi=epselon/((epselon-1)*(1-alpha))*psi;
psi=a^(-1/(1-alpha))*C^((1/(1-alpha))-sigma)*w+teta*beta*(rgm^(-epselon*(1-mu)*chi/(1-alpha)))*pi^(-mu*chi*epselon/(1-alpha))*(pi(+1))^(epselon/(1-alpha))*psi(+1);
phi=C^(1-sigma)+teta*beta*rgm^((1-epselon)*(1-mu)*chi)*pi^(chi*mu*(1-epselon))*(pi(+1))^(epselon-1)*phi(+1);
N=s*(C/a)^(1/(1-alpha));
s=(1-teta)*x^(-epselon/(1-alpha))+teta*rgm^(chi*(1-mu)*epselon/(alpha-1))*pi(-1)^(-epselon*mu*chi/(1-alpha))*pi^(epselon/(1-alpha))*s(-1);
i/intssv=(i(-1)/intssv)^ALFAI * ((pi/rgm)^ALFAPAI *(C/cssv)^ALFAY)^(1-ALFAI);
i=(1+r)*pi(+1);
vf=(log(C)-dn*N^(1+newphi)/(1+newphi))+beta*vf(+1);
end;

%if sigma different from one then use
%vf=(C^(1-sigma)/(1-sigma)-dn*N^(1+newphi)/(1+newphi))+beta*vf(+1);



initval;
rgm=gss0;
intssv=iss+1;
cssv=Css;
a=Ass;
i=iss+1;
pi=gss0;
r=rss;
C=Css;
phi=phiss;
x=xss;
s=sss;
w=wss;
N=Nss;
psi=psiss;
vf=vfss;
end;

%steady;
check;

endval;
rgm=gss1;
intssv=1+iss1;
cssv=Css1;
a=Ass1;
i=1+iss1;
pi=gss1;
r=rss;
C=Css1;
phi=phiss1;
x=xss1;
s=sss1;
w=wss1;
N=Nss1;
psi=psiss1;
vf=vfss1;
end;

%steady;
check;



simul (periods=100);

annrint=100*(((1+r).^4)-1);                 // annual inerest
annnomrint=100*(((1+i).^4)-1);              // annual nominal interest
annpi=100*(pi.^4-1);                        // annual inflation    

output_path(:,ii)  = 100*(C/Css1-ones(size(C,1),1));
s_path(:,ii)       = s;
w_path(:,ii)       = 100*(w/wss1-ones(size(C,1),1));
n_path(:,ii)       = 100*(N/Nss1-ones(size(C,1),1));
c_path(:,ii)       = 100*(C/Css1-ones(size(C,1),1));
annpi_path(:,ii)   = annpi;
annrint_path(:,ii) = annrint;
annnomrint_path(:,ii)=annnomrint;

ii = ii +1;

end


%==========================================================================
% PLOTS
%==========================================================================

perio=16;
tt=0:1:(perio-1);

figure(1)
set(0,'DefaultAxesFontName','Times','DefaultAxesFontSize',11);

subplot(2,2,1)
plot(tt,output_path(1:perio,1), 'k-',tt,output_path(1:perio,2), 'k-.',tt,output_path(1:perio,3), 'k-+','LineWidth',2);
grid on
hold on
xlabel('Quarters');
title('OUTPUT PATH (% dev. from new SS)')
legend('\pi= 4%','\pi = 6%','\pi = 8%')

subplot(2,2,2)
plot (tt,w_path(1:perio,1),'k-',tt,w_path(1:perio,2),'k-.',tt,w_path(1:perio,3),'k-+','LineWidth',2);
hold on
xlabel('Quarters');
legend('\pi= 4%','\pi = 6%','\pi = 8%')
grid on 
title('REAL WAGE (% dev. from new SS)')

subplot(2,2,3)
plot(tt,annpi_path(1:perio,1), 'k-',tt,annpi_path(1:perio,2), 'k-.',tt,annpi_path(1:perio,3), 'k-+','LineWidth',2);
grid on
hold on
xlabel('Quarters');
title('INFLATION IN % (annualized)')
legend('\pi= 4%','\pi = 6%','\pi = 8%')

subplot(2,2,4)
plot (tt,annnomrint_path(1:perio,1),'k-',tt,annnomrint_path(1:perio,2),'k-.',tt,annnomrint_path(1:perio,3),'k-+','LineWidth',2);
hold on
xlabel('Quarters');
legend('\pi= 4%','\pi = 6%','\pi = 8%')
grid on 
title('NOM. INTEREST RATE IN % (annualized)');  

return

figure(2)
set(0,'DefaultAxesFontName','Times','DefaultAxesFontSize',11);

subplot(2,2,1)
plot(tt,mc_path(1:perio,1), 'k-',tt,mc_path(1:perio,2), 'k-.',tt,mc_path(1:perio,3), 'k-+','LineWidth',2);
grid on
hold on
xlabel('Quarters');
title('REAL MARGINAL COSTS (% dev. from new SS)')
legend('\pi= 4%','\pi = 6%','\pi = 8%')

subplot(2,2,2)
%plot(tt,mc_path(1:perio,1), 'k-',tt,mc_path(1:perio,2), 'k-.',tt,mc_path(1:perio,3), 'k-+','LineWidth',2);
plot(tt,s_path(1:perio,1), 'k-',tt,s_path(1:perio,2), 'k-.',tt,s_path(1:perio,3), 'k-+','LineWidth',2);
grid on
hold on
xlabel('Quarters');
%title('REAL MARGINAL COSTS (% dev. from new SS)')
title('PRICE DISPERSION')
legend('\pi= 4%','\pi = 6%','\pi = 8%')

subplot(2,2,3)
plot (tt,c_path(1:perio,1),'k-',tt,c_path(1:perio,2),'k-.',tt,c_path(1:perio,3),'k-+','LineWidth',2);
hold on
xlabel('Quarters');
legend('\pi= 4%','\pi = 6%','\pi = 8%')
grid on 
title('CONSUMPTION (% dev. from new SS)')


subplot(2,2,4)
plot (tt,n_path(1:perio,1),'k-',tt,n_path(1:perio,2),'k-.',tt,n_path(1:perio,3),'k-+','LineWidth',2);
hold on
xlabel('Quarters');
legend('\pi= 4%','\pi = 6%','\pi = 8%')
grid on 
title('HOURS (% dev. from new SS)')