clc
clear all
addpath C:\dynare\4.4.3\matlab  
%Initial parameterization

%rho_Psi=0.95;
sigma=0.0072;
rho_z=0.95;
gam=0;
chi=0;
b=0;
A=0.5;
eps=0.1;
rho=1/1.4;
%Quarterly calibration
f=0.36^(1/3);
betta = 0.99; 
delta=(1.10)^(1/4)-1;
u=0.051;
s=(f*u/(1-u)-delta)/(1-delta);
theta=0.51;
 lambda = 0.5; 
B = f/(theta^(1-lambda));
 Z = 1; 
 Psi= 0.7; 
 r = (1 - betta)/betta; 
nu = 0.5;
R=1.0025;
Rbar = Psi*(1 + r)/(r*nu + Psi);
omega = 0.1;
k=0.1;
%Store variables%
save par theta u Z f;
%run dynare
dynare OBC_DSCES_392016.mod noclearall nograph
irf=[qm_e1 qc_e1 Q_e1 Qt_e1 qt_e1 a_e1 at_e1 d_e1 Omega_e1 p_e1 nbar_e1 nm_e1 nc_e1 S_e1 MC_e1 w_e1 U_e1 ...
    theta_e1 f_e1 Z_e1 V_e1 Y_e1 C_e1 wr_e1 P_e1 share_e1];
ss=oo_.steady_state;
irf_ss=log((irf+repmat(ss',40,1))./repmat(ss',40,1))% perc. deviation from steady state

a11=irf_ss(:,6);
d11=irf_ss(:,8);
p11=irf_ss(:,10);
nbar11=irf_ss(:,11);
S11=irf_ss(:,14);
w11=irf_ss(:,16);
U11=irf_ss(:,17);
Y11=irf_ss(:,22);
C11=irf_ss(:,23);
wr11=irf_ss(:,24);
P11=irf_ss(:,25);
share11=irf_ss(:,26);
%D11=D_e1;
clear irf
%%
% Calculate moments from simulated data
% Log and detrend
% Y=BK(log(Y),6,32,12);%
% U=BK(log(U),6,32,12);% cyclical component in level
% C=BK(log(C),6,32,12);
% w=BK(log(w),6,32,12);
% V=BK(log(V),6,32,12);
% d=BK(log(d),6,32,12);
% a=BK(log(a),6,32,12); %cyclical component in logs
% theta=BK(log(theta),6,32,12);
% 
% 
% data=[Y U C w V d a];
% nvar=size(data,2);
% sd=var(data)'.^(1/2);
% cor=tril(corr(data));
% cor=reshape(cor,nvar^2,1);
% cor(cor==0)=[];
% cor(cor==1)=[];
% init=zeros(3,1);
% for i=1:nvar
% acf=[init autocorr(data(:,i),2)];
% init=acf;
% end
% acf=reshape(acf,3*(nvar+1),1);
% acf(acf==0)=[];
% acf(acf==1)=[];
% acf1=acf(1:2:end);
% acf2=acf(2:2:end);
% names=char('Y', 'U', 'C', 'w', 'V', 'd', 'a')
% out=[sd, acf1, acf2]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Change to low taste for
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%variety%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
rho=1/1.1; 
 gam=14081.8138;
dynare OBC_DSCES_nt noclearall
irf=[qm_e1 qc_e1 Q_e1 Qt_e1 qt_e1 a_e1 at_e1 d_e1 Omega_e1 p_e1 nbar_e1 nm_e1 nc_e1 S_e1 MC_e1 w_e1 U_e1 ...
    theta_e1 f_e1 Z_e1 V_e1 Y_e1 C_e1 wr_e1 P_e1 share_e1];
ss=oo_.steady_state;
irf_ss=log((irf+repmat(ss',40,1))./repmat(ss',40,1))% perc. deviation from steady state
a12=irf_ss(:,6);
d12=irf_ss(:,8);
p12=irf_ss(:,10);
nbar12=irf_ss(:,11);
S12=irf_ss(:,14);
w12=irf_ss(:,16);
U12=irf_ss(:,17);
Y12=irf_ss(:,22);
C12=irf_ss(:,23);
wr12=irf_ss(:,24);
P12=irf_ss(:,25);
share12=irf_ss(:,26);
% %impulse responses
 fig=figure;
 
 
subplot(3,3,1)
plot(U11)
 a=[cellstr(num2str(get(gca,'ytick')'*100))]; 
% % Create a vector of '%' signs
%      pct = char(ones(size(a,1),1)*'%'); 
% % Append the '%' signs after the percentage values
%      new_yticks = [char(a),pct];
% % 'Reflect the changes on the plot
%      set(gca,'yticklabel',new_yticks) 
title('Unemployment');
xlabel('Periods');
hold all
plot(U12)
legend('markup=1.4','markup=1.1','Location','northeast','Fontsize',8)

subplot(3,3,2)
plot(S11)
%  a=[cellstr(num2str(get(gca,'ytick')'*100))]; 
% % Create a vector of '%' signs
%      pct = char(ones(size(a,1),1)*'%'); 
% % Append the '%' signs after the percentage values
%      new_yticks = [char(a),pct];
% % 'Reflect the changes on the plot
%      set(gca,'yticklabel',new_yticks) 
title('Number of Varieties');
xlabel('Periods');
hold all
plot(S12)

%legend('eta=3.5','eta=10','Location','northeast','Fontsize',8)

subplot(3,3,3)
plot(wr11)
%  a=[cellstr(num2str(get(gca,'ytick')'*100))]; 
% % Create a vector of '%' signs
%      pct = char(ones(size(a,1),1)*'%'); 
% % Append the '%' signs after the percentage values
%      new_yticks = [char(a),pct];
% % 'Reflect the changes on the plot
%      set(gca,'yticklabel',new_yticks) 
title('Wages');
xlabel('Periods');
hold all
plot(w12)
%legend('eta=3.5','eta=10','Location','northeast','Fontsize',8)
hold off

subplot(3,3,4)
plot(d11)
%  a=[cellstr(num2str(get(gca,'ytick')'*100))]; 
% % Create a vector of '%' signs
%      pct = char(ones(size(a,1),1)*'%'); 
% % Append the '%' signs after the percentage values
%      new_yticks = [char(a),pct];
% % 'Reflect the changes on the plot
%      set(gca,'yticklabel',new_yticks) 
title('Debt');
xlabel('Periods');
hold all
plot(d12)
hold off

subplot(3,3,5)
plot(P11)
 a=[cellstr(num2str(get(gca,'ytick')'*100))]; 
% % Create a vector of '%' signs
%      pct = char(ones(size(a,1),1)*'%'); 
% % Append the '%' signs after the percentage values
%      new_yticks = [char(a),pct];
% % 'Reflect the changes on the plot
%      set(gca,'yticklabel',new_yticks) 
title('Prices');
xlabel('Periods');
hold all
plot(p12)
hold off

subplot(3,3,6)
plot(nbar11)
 a=[cellstr(num2str(get(gca,'ytick')'*100))]; 
% % Create a vector of '%' signs
%      pct = char(ones(size(a,1),1)*'%'); 
% % Append the '%' signs after the percentage values
%      new_yticks = [char(a),pct];
% % 'Reflect the changes on the plot
%      set(gca,'yticklabel',new_yticks) 
title('Firm Size');
xlabel('Periods');
hold all
plot(nbar12)
hold off

subplot(3,3,6)
plot(a11)
 a=[cellstr(num2str(get(gca,'ytick')'*100))]; 
% % Create a vector of '%' signs
%      pct = char(ones(size(a,1),1)*'%'); 
% % Append the '%' signs after the percentage values
%      new_yticks = [char(a),pct];
% % 'Reflect the changes on the plot
%      set(gca,'yticklabel',new_yticks) 
title('Liquid Assets');
xlabel('Periods');
hold all
 plot(a12)
 hold off

subplot(3,3,7)
plot(C11)
title('Consumption');
xlabel('Periods');
hold all
plot(C12)
hold off

subplot(3,3,8)
plot(Y11)
title('Output');
xlabel('Periods')
hold all
plot(Y12)
hold off

subplot(3,3,9)
plot(share11)
title('Share of Monopolistically Competitive Sector');
xlabel('Periods')
hold all
plot(share12)
hold off
  filename='C:\Users\mariors\Documents\first field paper\summer research 2015\Work\credit unemployment variety limited commitment\prod_shock'
 print(fig,'-dpdf',filename)
 
 
 %%
 %%Exogenous debt limit
 rho=1/1.3;
 gam=13.1692;
 dbar=ss(8);
 dynare OBC_DSCES_exogenous_debt noclearall
irf=[qm_e1 qc_e1 Q_e1 Qt_e1 qt_e1 a_e1 at_e1 p_e1 nbar_e1 nm_e1 nc_e1 S_e1 MC_e1 w_e1 U_e1 ...
    theta_e1 f_e1 Z_e1 V_e1 Y_e1 C_e1 wr_e1 P_e1 share_e1];
ss=oo_.steady_state;
irf_ss=log((irf+repmat(ss',40,1))./repmat(ss',40,1))% perc. deviation from steady state
a13=irf_ss(:,6);
p13=irf_ss(:,8);
nbar13=irf_ss(:,9);
S13=irf_ss(:,12);
w13=irf_ss(:,14);
U13=irf_ss(:,15);
Y13=irf_ss(:,20);
C13=irf_ss(:,21);
wr13=irf_ss(:,22);
P13=irf_ss(:,23);
share13=irf_ss(:,24);
% %%
%  filename='C:\Users\mariors\Documents\first field paper\summer research 2015\Work\credit unemployment variety limited commitment\credit_shock'
% % print(fig,'-dpdf',filename)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%TARGETS FROM DATA%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %unemployment
% %autocorr=0.94
% % cor(U,V)=-0.89
% %std(theta)=0.38
% %s=0.034
% %std(f)=0.12
% %cor(f,thet)=0.95
% %%%%%%%%%%%%%%%%%%%%%%
% %Model
% %current problem
% %autocorr(U)=0.9022
% %std(U)=0.314
% %std(f)=0.34
% %cor(f,thet)=1
% %cor(U,V)=-0.87
% %std(theta)=0.68
% %%%%job finding rate and market tightness are 
% %too volatile relative to the data. %Capital can act
% %as a dampener.
% %std(theta)=0.69
% %need to reduce volatility and increase persistence of
% %vacancies slightly
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% m=6; %state variables: theta, S, d,a_,U,
% n=13;
%  Psi = oo_.dr.ghx;
%  omega = oo_.dr.ghu;
%  ss = oo_.dr.ys;
% 
%  T = 200; % number of periods to simulate
% 
%  es = sigma*randn(T,1);
%  Xsim = zeros(n+m,T);
% 
%  Xsim(:,1) = omega*es(1,1);
% 
%  for j = 2:T
%  Xsim(:,j) = Psi*Xsim(1:6,j-1) + omega*es(j,1);
%  end
%  % add back in steady state values so that these are expressed as log
%  % levels, not log deviatins
%  for j = 1:T
%  Xsim(:,j) = Xsim(:,j) + ss;
%  end
%  
%  fig2=figure
% subplot(3,2,1)
% plot(Xsim(12,:))
% title('Number of sellers')
% 
% subplot(3,2,2)
% plot(Xsim(6,:))
% title('Debt level')
% 
% subplot(3,2,3)
% plot(Xsim(4,:))
% title('Liquid assets')
% 
% subplot(3,2,4)
% plot(Xsim(14,:))
% title('Wages')
% 
% subplot(3,2,5)
% plot(Xsim(15,:))
% title('Unemployment')
% 
% subplot(3,2,6)
% plot(Xsim(19,:))
% title('Productivity')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Impulse Response
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Functions%%%%%%%%%%%%%%%%%%