% TransitionDynamicsShocks.m
%
%   -- Called from TransitionDynamicsMaster with "CaseName", 
%      "ParamToShock" and "ParamToShockValues"
%
%  Simulate the transition dynamics of the inequality model
%  assuming entrepreneurial effort is fixed at its steady state value.
%  That is, we shut off the transition dynamics associated with entreprenuerial effort 
%  choice in the HJB equation (because solving that as well is very hard, while this is easy).
%
%  For the same reason, assume that s (research allocation) jumps to its steady state value.
%  So the only source of transition dynamics is in the evolution of the distribution of x.
%
%  Discretize: 200,000 entrepreneurs, with initial x drawn from initial SS distribution.
%   - Shock the economy (e.g. increase phiH) ==> new SS effort choice e* and new s* (==> delta* and g*)
%   - Let the discretized GBM evolve, so that the distribution of {x} evolves.
%   - Keep track of eta and y:=Y/Lbar and show how they evolve over time.
%
%  This is a rough guide to transition dynamics in the model. Goal is to illustrate two things:
%   1. Transition dynamics are fast for our estimated pbar.
%   2. Observed growth over 10 or 20 years looks different from SS growth bc of the level effect.
%
%  Graphs: show the initial steady-state for 20 years, then the shock hits in 1975.
%     How quickly does eta rise to the new steady state?
%     What does GDP look like?
%
%   (optional) SplineLength = rather than having "s" and "muH" jump to their new SS values,
%            have them rise gradually over SplineLength years.


diarychad(['TransitionDynamicsMaster_' CaseName]);

disp ' '; disp ' ';
disp '==========================================================================';
disp(['                     CaseName: ' CaseName]);
fprintf('           ParamToShockValues: %8.3f  %8.3f\n',ParamToShockValues);
fprintf('                NumGridPoints: %6.0f\n',NumGridPoints);
fprintf('                PL_percentile: %8.5f\n',PL_percentile);
fprintf('                           dt:    1 /%3.0f\n',1/dt);
disp '==========================================================================';
disp ' '; disp ' ';

% Create etaUS -- See FractalInequality.m
load WorldTopInequalityDatabase2016;
T=2015;
yrs_data=(1970:T)';
yy=yrs_data-yr0;
s3=is001(:,USA)./is01(:,USA);
etaUS=1+log10(s3(yy));


% Setup parameters
YearShock=1975
figname='TransitionDynamicsShocks';


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initial steady state
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eval([ParamToShock '=' 'ParamToShockValues(1);']);   % e.g. phiH=PhiH(1);
[Results,keep,eta_old,g_old,etaL_old]=ideamodelsolution(Omega,beta,rho,phiH,phiL,gamma,sigmaH,sigmaL,theta,lambda,zBar,Lbar,mBar,dltaBar,tau,pBar,qbar);
%tle='t eta gx100 mutldeH s beta gamma zbar phiH phiL sigmaH dlta dltabar tau PBar eH eL MBar Lbar';
s_old=Results(5);
eH=Results(16);
eL=Results(17);
muH_old=phiH*eH;
muL_old=phiL*eL;
pBar_old=pBar;

% Check Low condition
mutildeL=muL_old-.5*sigmaL^2;
mm=mutildeL/sigmaL^2;
etaL_old = 1/(-mm+sqrt(mm^2+2*dltaBar/sigmaL^2))
if etaL_old>eta_old; disp 'Wait! etaL > etaH for these parameter values!!!!'; wait; end;



% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% New steady state
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eval([ParamToShock '=' 'ParamToShockValues(2);']);   % e.g. phiH=PhiH(2);
[Results,keep,eta_new,g_new]=ideamodelsolution(Omega,beta,rho,phiH,phiL,gamma,sigmaH,sigmaL,theta,lambda,zBar,Lbar,mBar,dltaBar,tau,pBar,qbar);
%tle='t eta gx100 mutldeH s beta gamma zbar phiH phiL sigmaH dlta dltabar tau PBar eH eL MBar Lbar';
s_new=Results(5);
dlta=Results(12);
eH=Results(16);
eL=Results(17);
muH_new=phiH*eH;
muL_new=phiL*eL;
mutildeH=muH_new-.5*sigmaH^2;
mutildeL=muL_new-.5*sigmaL^2;


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Transition dynamics
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[eta,gdp,yrs,ytrend]=transitiondynamics_1shock(eta_old,mutildeH,mutildeL,sigmaH,sigmaL,pBar,qbar,dltaBar,NumGridPoints,alpha,s_old,s_new,g_old,YearShock,yrs_data,etaUS,figname,CaseName,YearsToShow,dt,PL_percentile);
chadtimer


% Save the results for later analysis
eval(['data_' ParamToShock '=table(yrs,eta,gdp,ytrend);']);


diary off
