% IdeaModel.m
%
%  Simulates the extended (H vs L) model with an endogenous allocation of
%  labor to research and shows comparative static results of growth and
%  inequality with respect to changes in key parameter values.
%
%  Produces the plots in Fiugre 8 "Steady State Comparative Statics"

clc
clear all;
diarychad('IdeaModel');
delete('IdeaModel.ps');
definecolors;


for i=1:8;
%for x=1:1
% Set baseline parameter values
SetParamsIdeaModel
NumPeriods=(2007-1980)+1;
oo=ones(NumPeriods,1);  


% Values we vary
if i==1; Phi=(.7:.01:1.57)'; oo=ones(length(Phi),1);  end;
if i==2; DltaBar=[(.05:.001:.12)']; oo=ones(length(DltaBar),1);   end;
if i==3; Beta=(.86:.02:2)'; oo=ones(length(Beta),1);   end;
%if i==4; MBar=(.2:.01:.7)'; oo=ones(length(MBar),1);   end;
if i==4; MBar=(.3:.005:.7)'; oo=ones(length(MBar),1);   end;
if i==5; Sigma=(.1:.005:.395)'; oo=ones(length(Sigma),1);   end;
if i==6; Tau=(0.12:.01:.75)'; oo=ones(length(Tau),1);   end;
if i==7; ZBar=(0:.01:.6)'; oo=ones(length(ZBar),1); end;
if i==8; PBar=(.82:.02:2)'; oo=ones(length(PBar),1); end;

% Generate the entire default time series for each parameter
if length(Phi)==1; Phi=Phi*oo; end;
if length(Beta)==1; Beta=Beta*oo; end;
if length(Gamma)==1; Gamma=Gamma*oo; end;
if length(ZBar)==1; ZBar=ZBar*oo; end;
if length(Sigma)==1; Sigma=Sigma*oo; end;
if length(DltaBar)==1; DltaBar=DltaBar*oo; end;
if length(MBar)==1; MBar=MBar*oo; end;
if length(Lbar)==1; Lbar=Lbar*oo; end;
if length(Tau)==1; Tau=Tau*oo; end;
if length(PBar)==1; PBar=PBar*oo; end;

SigmaH=Sigma;
SigmaL=Sigma;
PhiH=Phi;

% Makes mutildeL \approx 0 (at least at 2% growth).
% Important to do throughout the parameter space; otherwise eL<0 too easily...
PhiL=(Beta(1).*(rho+.02./log(Gamma(1))+DltaBar(1))+.5*Sigma(1).^2)./Omega(1);
if i==3; PhiL=.201; end;
if i==5; PhiL=.187; end;
PhiL=PhiL*oo;
%PhiL=(Beta.*(rho+.02./log(Gamma)+DltaBar)+.5*Sigma.^2)./Omega;

% Solve the model for all parameter combinations
[Results,keep,eta,g]=ideamodelsolution(Omega,Beta,rho,PhiH,PhiL,Gamma,SigmaH,SigmaL,theta,lambda,ZBar,Lbar,MBar,DltaBar,Tau,PBar,qbar);

xtick=[];
if i==1; x=Results(keep,9); xlab='x-technology, \phi_H'; figname='IdeaModel-phi'; end; %xtick=strmat('0.3 0.4 0.5 0.6 0.7'); end;
if i==2; x=Results(keep,13); xlab='Exogenous destruction, \delta bar'; figname='IdeaModel-dltabar';  end;
if i==3; x=Results(keep,6); xlab='Utility weight on leisure, \beta'; figname='IdeaModel-beta';  end;%xtick=strmat('0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9');
if i==4; x=Results(keep,18); xlab='Wage for failed research, m-bar'; figname='IdeaModel-mbar';  end;
if i==5; x=Results(keep,11); xlab='Stdev of shocks, \sigma'; figname='IdeaModel-sigma';  end;%xtick=strmat('0.2 0.3 0.4 0.5 0.6 0.7');
if i==6; x=Results(keep,14); xlab='Taxes, \tau'; figname='IdeaModel-tau'; end; %xtick=strmat('0 0.1 0.2 0.3 0.4 0.5'); end;
if i==7; x=Results(keep,8); xlab='Innovation blocking, z-bar'; figname='IdeaModel-zbar';  end;%xtick=strmat('0.01 0.02 0.03 0.04 0.05');
if i==8; x=Results(keep,15); xlab='Decay rate from high growth, p-bar'; figname='IdeaModel-pbar';  end; %xtick=strmat('0.01 0.02 0.03 0.04 0.05');
%x=Results(keep,14); xlab='Market size, L-bar'; figname='IdeaModel-Lbar'; 
estarL=Results(keep,17);

FS=20;   % FontSize for 2x2 plots for paper
FSb=16 ; % for chadfig2 axis label 
ShrinkSize=.93;

%FS=16
figure(1); figsetup; makefigwide
[ax,h1,h2]=plotyy(x,eta,x,g); 
darkfactor=.8; LW=4;
set(h1,'LineWidth',LW,'Color',myblue);
set(h2,'LineWidth',LW,'Color',mygreen,'LineStyle','--');
set(ax(1),'ycolor',myblue*darkfactor);
set(ax(2),'ycolor',mygreen*darkfactor);
set(ax(1),'FontSize',FS);
set(ax(2),'FontSize',FS);

%xh=xlabel(xlab,'FontWeight','bold'); %,'FontName','Helvetica','FontSize',12);
%set(xh, 'Units', 'Normalized')
%pos = get(xh, 'Position');
ytick1=strmat('0 0.25 0.50 0.75 1');
ytick2=strmat('0 1 2 3 4');
%ytick2=strmat('1.50 1.75 2.00 2.25 2.50');
yylabels(ax,'Power law inequality','Growth rate (percent)',ytick1,ytick2,xtick);
chadfig2(' ','Power law inequality',ShrinkSize,0,FSb);
chadfigyy(ax(2),'Growth rate (percent)',2.85,2,FSb);
% Uncomment this line to show negative estarL
%if any(estarL<0); plot(ax(1),x(estarL<0),eta(estarL<0),'rx'); end;
%set(xh,'String',xlab,'Position',[.5 -.1 0],'FontSize',FS,'HorizontalAlignment','center'); % restore original center position                                                                                                                
%set(xh,'String',xlab,'Position',[.5 -.1 0],'FontSize',FS,'Interpreter','latex','HorizontalAlignment','center'); % restore original center position
%set(xh,'String',xlab,'Position',pos,'FontSize',FS,'HorizontalAlignment','center'); % restore original center position
print('-depsc',figname);
%print('-dpng',figname);
print('-dpsc','-append','-bestfit','IdeaModel');


%wait;
end; % Loop over the parameters

diary off;

